3,724
Views
40
CrossRef citations to date
0
Altmetric
Review

Mathematical modeling of the thermal effects of irreversible electroporation for in vitro, in vivo, and clinical use: a systematic review

ORCID Icon, ORCID Icon, , ORCID Icon, ORCID Icon, , ORCID Icon, ORCID Icon, ORCID Icon & ORCID Icon show all
Pages 486-505 | Received 24 Sep 2019, Accepted 01 Apr 2020, Published online: 18 May 2020

Abstract

Introduction

Irreversible electroporation (IRE) is a relatively new ablation method for the treatment of unresectable cancers. Although the main mechanism of IRE is electric permeabilization of cell membranes, the question is to what extent thermal effects of IRE contribute to tissue ablation.

Aim

This systematic review reviews the mathematical models used to numerically simulate the heat-generating effects of IRE, and uses the obtained data to assess the degree of mild-hyperthermic (temperatures between 40 °C and 50 °C) and thermally ablative (TA) effects (temperatures exceeding 50 °C) caused by IRE within the IRE-treated region (IRE-TR).

Methods

A systematic search was performed in medical and technical databases for original studies reporting on numerical simulations of IRE. Data on used equations, study design, tissue models, maximum temperature increase, and surface areas of IRE-TR, mild-hyperthermic, and ablative temperatures were extracted.

Results

Several identified models, including Laplace equation for calculation of electric field distribution, Pennes Bioheat Equation for heat transfer, and Arrhenius model for thermal damage, were applied on various electrode and tissue models. Median duration of combined mild-hyperthermic and TA effects is 20% of the treatment time. Based on the included studies, mild-hyperthermic temperatures occurred in 30% and temperatures ≥50 °C in 5% of the IRE-TR.

Conclusions

Simulation results in this review show that significant mild-hyperthermic effects occur in a large part of the IRE-TR, and direct thermal ablation in comparatively small regions. Future studies should aim to optimize clinical IRE protocols, maintaining a maximum irreversible permeabilized region with minimal TA effects.

1. Introduction

Irreversible electroporation (IRE) is a relatively new ablation modality for the treatment of unresectable cancers [Citation1]. Contrary to contemporary ablation techniques such as radiofrequency ablation (RFA) and microwave ablation (MWA), IRE is proposed as a nonthermal ablation technique [Citation2]. Instead, IRE uses ultra-short, high-voltage pulses between electrodes to induce an altered membrane potential, which leads to the irreversible permeabilization of the membranes of tumor cells through the creation of nanopores [Citation3]. These nanopores allow for increased transmembrane transport, which causes a disturbance in the cell’s homeostasis and consequent cell death through accidental and regulated cell death [Citation4].

Since the proposed mechanism of IRE is nonthermal, the technique is supposedly less prone to a heat-sink effect, which is known to hamper the efficacy of RFA and MWA when performing ablations in the presence of large blood vessels [Citation5]. Another possible advantage of IRE over thermal modalities includes the preservation of vital structures (e.g. blood vessels, bile ducts) traversing the ablation zones, since fibrous and collagen structures are typically not affected by IRE [Citation6]. Therefore, the application of IRE is widely investigated for several indications, such as renal and prostate tumors [Citation7,Citation8], pancreatic cancer [Citation9–12], and perihilar cholangiocarcinoma [Citation10].

However, contrary to the current assumption that IRE is a nonthermal ablation modality, several studies do report on the presence of thermal ablative (TA) temperatures (T ≥ 50 °C) within the IRE ablation zones during treatment, with histological signs of thermal damage [Citation13–16]. Moreover, it is presently unknown whether the reported hyperthermic effects of IRE, caused by slightly less elevated (mild-hyperthermic (MH)) temperatures 40 ≤ T [°C] < 50 [Citation3,Citation17–19], have a beneficial impact on the irreversible permeabilization effect (IPE; we chose IPE instead of IRE effect to distinguish between only the permeabilization effect, and the permeabilization and the thermal effects jointly produced by IRE), or could be seen as an unwanted side-effect. Potentially, temperature increase during IRE ablation might be needed to enhance the IPE and therefore contribute to a successful ablation [Citation20]. On the other hand, IRE could cause thermal damage to vital structures if applied in the vicinity of thermally vulnerable structures. This leaves the true mechanism of work of IRE to be contentious.

Investigating the degree of mild-hyperthermic and thermally ablative effects in the IRE-treated region (IRE-TR) could provide more insight into the IRE process. Determination of these ratios in vivo or in clinical applications is challenging, since a limited amount of IRE settings and parameters can be explored. Simulations using mathematical models could therefore be very useful. In such simulations, the IRE settings and model parameters can be varied extensively, resulting in improvement of our understanding of the IRE mechanism.

Therefore, the aim of this systematic review is to summarize and review the mathematical models used to numerically simulate the heat-generating effect of IRE, and to assess the degree of mild-hyperthermic and thermally ablative effects of IRE in relation to the IRE-TR, as predicted by the simulation results provided by the studies that were included in this review. Accordingly, we investigated whether the analytical/numerical model calculations reported for IRE in the literature could be used to answer the following research questions:

  • Does IRE produce a thermal elevation that results in mild-hyperthermic effects?

  • Does IRE produce a thermal elevation that results in thermally ablative effects?

  • What are the relative sizes expressed in percentages of the regions with the mild-hyperthermic and thermally ablative effects in the IRE-TR?

2. Methods

This systematic review was performed according to the Preferred Reporting Items for Systematic Reviews and Meta (PRISMA)-analyses guidelines [Citation21,Citation22]. In this article, the main data and equations were extracted, analyzed, and reviewed. Additional description of data or equations was added to Supplementary Appendices 1–6:

2.1. Study selection

The medical and technical databases and search engines (Cochrane Library, Embase, IEEE Xplore Digital Library, PubMed, ScienceDirect, Scopus, and Web of Science) were searched for journal articles from inception up to and including 8 February 2019 using search terms related to ‘electroporation’, ‘mathematical models’, and ‘temperature’. A full description of the search strategy can be found in Supplementary Appendix 1. After removing the duplicates and excluding studies written in any other language than the English language, we searched for the terms ‘irreversible electroporation’ OR ‘ablation’ using EndNote X9, restricting the search terms to ‘Any Field’. Subsequently, the titles, the abstracts, and the full-text were screened for eligibility based on predefined inclusion and exclusion criteria that can be found in Sub section 2.2. Any disagreement was resolved by discussion between the review authors (PA) and (EvV). The flow diagram of the study selection is shown in .

Figure 1. Flow diagram of the study selection.

Figure 1. Flow diagram of the study selection.

2.2. Eligibility criteria

The focus of the studies is about IRE with rectangular pulses and a pulse length tP ≥ 50 × 10−6 s. They must include analytical and/or numerical analysis of IRE and its heat-generating effect, performed in mammalian subject at macroscopic scale (i.e. tissue, scaffolds, cell suspension, phantoms including mammalian cell suspension). The inclusion criterion in the meta-analysis was the inclusion of a calculated temperature increase with the corresponding pulse parameters, and electrical and thermal properties of the target tissue. Studies that do not satisfy the criterion of the meta-analysis, or studies with the focus on the combination of IRE with other electroporation treatments (e.g. reversible electroporation, high-frequency IRE) or thermal ablations were included in qualitative synthesis.

Studies were excluded, if: (1) they exclude analytical and/or numerical analysis of the heat-generating effect of IRE, and the electrical properties of the target subject are not reported, (2) they applied IRE to nonmammalian subject, (3) they only include analytical and/or numerical analysis at microscopic scale (i.e. single cell, multiple cells, planar lipid membrane), (4) tP <50 × 10−6 s or the IRE pulses are nonrectangular, (5) their focus related to general description of electroporation, pulses with low electric field strength, reversible electroporation (including agent/drug/gene delivery) or thermal ablation, (6) they are a review, (7) their focus is about design of pulse generator.

2.3. Data collection

Parameters and equations, which were used for the analytical and/or numerical analysis of IRE and its heat-generating effect, were extracted from studies included in qualitative synthesis. Data of the included studies were extracted by one author (PA) and checked by a second author (EvV). Any disagreement was resolved by discussion between the review authors. Abbreviations, parameters, quantities, and study characteristics extracted from the included studies were explained in Supplementary Appendix 2. The extracted equations were reviewed and summarized in Section 3 and extensively described in Supplementary Appendix 3. After explaining the analysis of certain study characteristics in Supplementary Appendix 5 (see Subsection 2.4), we summarized the extracted parameters and study characteristics in four tables that are shown in Supplementary Appendix 6.

2.4. Methods of analysis

2.4.1. Definitions of mild-hyperthermic and thermally ablative effects

Assuming the tissues to have an initial physiological temperature (Tinit) of 37 °C and a thermally ablative threshold (Tth) of 50 °C (ΔT ≥ 13 °C), the following assumptions were made:

  • The thermal elevation of 3 ≤ ΔT [°C] < 13 (40 ≤ T [°C] < 50) results in mild-hyperthermic effects.

  • The thermal elevation of ΔT ≥ 13 °C (T ≥ 50 °C) results in thermally ablative effects.

  • Accordingly, in the IRE-TR the relative sizes were calculated based on the regions with the thermal elevations of 3 ≤ ΔT [°C] < 13 for the mild-hyperthermic effects and ΔT ≥ 13 °C for the thermally ablative effects.

Please note that even though the value 50 °C may not instantaneously result in thermal damage, it was still chosen as a Tth. This because during IRE the duration of exposure to temperatures around 50 °C could be in the order of minutes [Citation15,Citation16], which is sufficiently long that thermal damage may indeed occur [Citation14,Citation23] and this is thus effectively a clinically relevant threshold value associated with a high likelihood of occurrence of thermal ablation.

2.4.2. Analysis of research questions

To answer the first two research questions, we extracted the maximum simulated temperature increase (ΔTmax [°C]) with the corresponding pulse parameters from studies included in meta-analysis. If the data extracted from the pulse parameters and/or the maximum simulated temperature was incomplete, the ΔTmax was not reported. In case of parametric study, in which ΔTmax was produced from several pulse parameters, the ΔTmax at the extreme ends of the figures were included.

For the third research question, if a study included 2D plots of electric field, and temperature distributions at the same position of the same model, then S3ΔT13,Σ [m2], SΔT13,Σ [m2], and SE-IRE(th),Σ [m2] were extracted, where S3ΔT13,Σ is the total surface area of the simulated temperature increase within the range 3 ≤ ΔT [°C] < 13 that represents the region with mild-hyperthermic temperature increase, SΔT13,Σ is the total surface area of the simulated temperature increase with ΔT ≥ 13 °C that represents the thermally ablated region, and SE-IRE(th),Σ is the total surface area of the simulated electric field with E ≥ EIRE(th) that represents the IRE-TR, E [V⋅m−1] is the magnitude of the electric-field strength and EIRE(th) [V⋅m−1] is the electric-field threshold of IRE; the minimum E required to ablate the target tissue. Even though EIRE(th) was only used as a minimum required E to extract SE-IRE(th),Σ for this analysis, it must be noticed that EIRE(th) in reality also depends on pulse parameters (pulse voltage, pulse shape, pulse length, pulse number, and pulse frequency), temperature, and electrical conductivity of the target. It is also important to note that:

  • the EIRE(th) values correspond to the EIRE(th) from the study from which SE-IRE(th),Σ was extracted,

  • the pulse parameters on which the values of ΔTmax, S3ΔT13,Σ and SΔT13,Σ depend, are different for each study,

  • in case of EIRE(th) range, e.g. [EIRE(th),1, EIRE(th),2], the average of this range was used to calculate SE-IRE(th),Σ. Specifically, EIRE(th) = ½ (EIRE(th),1 + EIRE(th),2).

  • if EIRE(th) was not reported in a certain study, then the SE-IRE(th),Σ was not calculated and, accordingly, included in the analysis.

The extraction of the surface areas was done excluding surface areas of the electrodes (Select [m2]). Using the software package ImageJ (National Institutes of Health, Bethesda, MD), we firstly imported the desired figure in ImageJ and assigned the pixels to distance ratio of the image using the known distance between the electrodes. Next, if the imported figure includes E- and ΔT-contours showing the conditions mentioned above, then we manually selected SE [m2], SΔT [m2] and Select using the ‘Freehand selections’ tool; SE is the selected surface area of the electric field strength and SΔT is the selected surface area of the temperature increase. In case of no contours, we used the ‘Threshold Color’ tool to set the threshold of the surface color to the minimum value mentioned above. Circular or rectangular electrodes were manually selected using ‘Oval selections’ tool, and ‘Rectangle’ tool. After the extraction, the ratios R3ΔT13 [%] and RΔT13 [%] were calculated, where R3ΔT13 is the ratio between the region with mild-hyperthermic temperature increase and IRE-TR, and RΔT13 is the ratio between the thermally ablated region and IRE-TR. The accurate description of the extracted surface areas and calculated ratios are shown in Supplementary Appendix 5.

2.5. Statistical analysis

Statistical analysis was applied to the data of the studies included in the meta-analysis to answer the three research questions. Specifically, we answered the first two research questions by calculating the median, the minimum, and the maximum of ΔTmax. For the third research question, the medians were calculated of R3ΔT13, and RΔT13. The answers of all research questions were based on studies with experimentally applied pulse parameters. As an additional step, the linear regression was applied to model the relationship between SE-IRE(th), R3ΔT13 and RΔT13, and the energy density (u [J⋅m−3]; see Subsection 3.3.), ΔTmax, R3ΔT13, and RΔT13.

2.6. Additional assumptions, definitions, methods of analysis, and notes

To preserve the continuous course, all remaining assumptions, definitions, methods of analysis, and notes were defined in Section 3 when the need arose.

3. Results

3.1. Study selection

After removal of duplicates and screening 3815 articles by title, abstract, and full-text, we included 58 articles [Citation3,Citation24–80] in the qualitative synthesis, and 26 articles [Citation3,Citation24–28,Citation31,Citation33,Citation34, Citation36,Citation37,Citation39–44,Citation48,Citation49,Citation55,Citation59,Citation60,Citation64,Citation66,Citation73,Citation80]) in the quantitative synthesis (meta-analysis and meta-regression). The full overview of the study selection is shown in . The parameters extracted from the included studies are shown in Supplementary Appendix 6, with the definitions explained in Supplementary Appendix 2.

3.2. Review of mathematical description of IRE

In this section, we summarized and reviewed the mathematical equations used to describe the thermal behavior of IRE in multiple models with various types of model geometries, electrodes, pulse parameters, tissue properties, and boundary conditions. The extended description of the mathematical equations can be found in Supplementary Appendix 3. In addition, more details about the applied model parameters were included in Supplementary Appendix 6. In the included studies, the subjects were modeled in one-dimensional [Citation25,Citation30,Citation50,Citation64], two-dimensional [Citation3,Citation24,Citation26–28,Citation32–34,Citation36–39,Citation41,Citation51–54,Citation59,Citation65,Citation66,Citation72,Citation74,Citation75,Citation78], and/or three-dimensional space [Citation26,Citation29,Citation31,Citation35,Citation40,Citation42–49,Citation55–58,Citation60–63,Citation67–71,Citation73,Citation77,Citation79,Citation80]. The calculations were performed analytically, or numerically using Finite Element Method or Finite Difference Method using the software packages e.g. CFdesign (Blue Ridge Numerics, Inc., Charlottesville, VA) [Citation55], COMSOL Multiphysics (COMSOL AB, Stockholm, Sweden) [Citation26,Citation28,Citation31,Citation32,Citation34–46,Citation48,Citation49,Citation51–53,Citation56–63, Citation65–67,Citation70,Citation71,Citation73,Citation74,Citation76–78], FEMLAB (COMSOL AB, Stockholm, Sweden) [Citation3,Citation24,Citation27], IRENA [Citation79,Citation81], Marc/Mentat (MSC Software Corp., Santa Ana, CA) [Citation47,Citation75], and QuickField (Terra Analysis, Denmark) [Citation54]. Generally, numerical calculations were performed on three-dimensional subjects using the software package COMSOL Multiphysics, which uses Finite Element Method. Furthermore, an example of possible pulse sequence is shown in . Here, VP [V] is the electric potential (the voltage) of the pulses, tP [s] is the duration of a single pulse, nP [–] is the pulse number, NP [–] is the total number of pulses, τP [s] is the duration between two pulses, and fP [Hz] is the pulse frequency, which is calculated as (1) fP=1tP+τP(1)

Figure 2. Overview of IRE pulses, with VP [V] as the electric potential (the voltage) of the pulses, tP [s] as the duration of a single pulse, nP as the pulse number, NP as the total number of pulses, τP [s] as the duration between two pulses, and fP [Hz] as the pulse frequency.

Figure 2. Overview of IRE pulses, with VP [V] as the electric potential (the voltage) of the pulses, tP [s] as the duration of a single pulse, nP as the pulse number, NP as the total number of pulses, τP [s] as the duration between two pulses, and fP [Hz] as the pulse frequency.

For convenience, overview of definitions, symbols, and mathematical notations were summarized in Supplementary Appendix 2.

3.2.1. Electric-field distribution

3.2.1.1. Electric-potential models

The IRE can be mathematically described at macroscopic scale by calculating the electric-field [Citation3,Citation24,Citation26–29,Citation31,Citation33,Citation34,Citation36–48,Citation50–59,Citation61–63,Citation65,Citation66, Citation68–70,Citation72–75,Citation78–80], and temperature distribution [Citation3,Citation24, Citation26–31,Citation33,Citation34,Citation36,Citation37,Citation39–44,Citation48,Citation53,Citation55,Citation58,Citation59,Citation61,Citation62,Citation64,Citation66,Citation70,Citation73,Citation76,Citation80] in a volume of a certain medium (e.g. tissue). This description starts with the calculation of the electric-potential distribution, which is obtained by simplifying the combination of Maxwell’s equations: Faraday’s law of induction and Ampère’s law, into the steady-state form (2) ·σΦ=0(2) where ∇ ⋅ [m−1] is the divergence, σ [S⋅m−1] is the electrical conductivity, ∇ [m−1] is the gradient and Φ [V] is the scalar electric potential [Citation3,Citation24,Citation26,Citation27,Citation29,Citation31,Citation33,Citation34,Citation36–48,Citation50–59,Citation62,Citation63, Citation65,Citation66,Citation68–70,Citation75,Citation78–80]. Here, the vector quantities are expressed in bold and italic, while the scalar quantities have remained unchanged. Furthermore, the electrical conductivity of a medium can be assumed to be homogeneous, heterogeneous, isotropic, anisotropic, linear, and/or nonlinear (see Supplementary Appendix 2 Table A2.2 for the definitions of the medium properties and Supplementary Appendix 3 for extended description of medium properties). This electrical model can be used when several media are modeled with various electrical characteristics, such as cancerous tissue in a healthy tissue model, and tissue that includes blood vessels and nerves. Also, it is possible to choose a volume with homogenous, isotropic, and linear medium properties for simplicity. In this case, we can remove the electrical-conductivity term from EquationEquation (2) and transform it into the Laplace’s equation (3) 2Φ=0(3) with ∇2 [m−2] as Laplace operator.

3.2.1.2. Electric-field model and threshold

For both EquationEquations 2 and Equation3, the calculated electric-potential distribution is used to determine the electric-field distribution in the volume of the medium using (4) E=Φ(4) with E [V⋅m−1] as the electric field [Citation26,Citation43,Citation68,Citation72,Citation74,Citation79,Citation80]. Assuming the volume of the medium to be an arbitrary tissue, the electric-field distribution yields IRE when a minimum electric-field strength EIRE(th) is obtained, providing the corresponding duration, number, shape, and the frequency of the pulses, and the tissue type. In Supplementary Appendix 6, a range of EIRE(th) and its corresponding pulse parameters are shown for several tissue types. Additional details about the mathematical description of the electric-field component of IRE are shown in Supplementary Appendix 3.

3.2.1.3. Boundary conditions

For the calculation of the electric potential, a boundary condition (BC) needs to be applied to the boundaries between the electrodes (metallic part) and the media (BEM), and the boundaries at the outer surface (BOS) of the model. These boundaries could have a constant electric potential (Dirichlet BC; always applied by the included studies to BEM) [Citation3,Citation24,Citation26–29,Citation31–48,Citation50–55,Citation57–61, Citation64–66,Citation68–76,Citation78–80], zero electric current (Neumann BC; generally applied to BOS) [Citation3,Citation24,Citation26,Citation28,Citation29,Citation31–34,Citation36–45,Citation48,Citation53, Citation57–61,Citation66,Citation69,Citation70,Citation73–75,Citation78–80], or limited electric current to reduce the influence of the boundaries on the domain by virtually increasing the size of the domain while maintaining the number of nodes (Robin BC) [Citation79]. See .

Figure 3. Examples of applied boundary conditions for the calculations of the electric potential distribution with arrow representing electric-current flow through the outer surface of the model (IBOS [A]). The red and the blue circles represent BEM with the electrodes given in the color gray, and black lines represent the BOS of the model. (A) shows Dirichlet boundary conditions at BEM and BOS, and (B) shows Dirichlet boundary conditions at BEM, and Neumann boundary conditions at BOS if IBOS = 0 A or Robin boundary conditions if IBOS is limited.

Figure 3. Examples of applied boundary conditions for the calculations of the electric potential distribution with arrow representing electric-current flow through the outer surface of the model (IBOS [A]). The red and the blue circles represent BEM with the electrodes given in the color gray, and black lines represent the BOS of the model. (A) shows Dirichlet boundary conditions at BEM and BOS, and (B) shows Dirichlet boundary conditions at BEM, and Neumann boundary conditions at BOS if IBOS = 0 A or Robin boundary conditions if IBOS is limited.

The preference of the boundary conditions depends on the type of the study. In case of active electrodes, obviously Dirichlet BC should be applied, while Neumann BC could be applied to the nonmetallic parts of the electrodes. For BOS, it is recommended to apply the Neumann BC for the clinical practice, since the simulated electric field can provide an upper limit to the temperature distribution [Citation3].

3.2.2. Temperature distribution

3.2.2.1. Bioheat transfer models

Several heat transfer models were used for the calculation of thermal distributions (see Supplementary Appendix 3 for an extended summary). Nevertheless, the most utilized model is the Pennes bioheat equation (5) ρcpTt=k·TρbωbcbTTart+Qm+J·E(5) where ρ [kg⋅m−3] is the mass density, cp [J⋅kg−1⋅°C−1] is the specific heat capacity, T [°C] is the temperature, t [s] is the time, k [W⋅m−1⋅°C−1] is the thermal conductivity, cb [J⋅kg−1⋅°C−1] is the specific heat capacity of the blood, ρb [kg⋅m−3] is the blood density, ωb [s−1] is the blood perfusion rate (blood perfusion wb [kg⋅m−3⋅s−1] = ρb ⋅ ωb), Tart [°C] is the arterial blood temperature, Qm [W⋅m−3] is the metabolic heat generation, and J · E [W⋅m−3] is the heat generation rate per unit volume (the Joule Heating) [Citation3,Citation24,Citation26–29,Citation31,Citation33,Citation34,Citation36, Citation37,Citation39,Citation41–44,Citation48,Citation53,Citation55,Citation58,Citation61,Citation62,Citation66,Citation70,Citation73,Citation76,Citation80]. Over the intrapulses of IRE, the Joule Heating term can be simplified into (6) J·E=σE2(6) when IRE is used [Citation3,Citation24,Citation26–29,Citation33,Citation34,Citation36,Citation40,Citation41,Citation44,Citation48,Citation53,Citation55,Citation58,Citation59,Citation66,Citation76].

Even though the majority of the studies applied the Pennes bioheat model, the choice of a suitable heat transfer model depends on the modeling conditions. Yet, for the clinical practice the Pennes bioheat model is recommended due to its ability to calculate the temperature distribution in a domain with heterogeneous and nonlinear electrical and thermal characteristics, considering the metabolic heat generation and heat sink effect.

3.2.2.2. Boundary conditions

To calculate the thermal distribution, boundary conditions need to be applied to the BEM and BOS of the model. These boundaries could have a constant temperature (Dirichlet BC) [Citation28,Citation33,Citation34,Citation36,Citation41,Citation53,Citation55], a limited heat flow (Robin BC) [Citation25,Citation29,Citation30,Citation37,Citation40,Citation43–45,Citation59], or no heat flow at all (Neumann BC) [Citation3,Citation24,Citation26,Citation27,Citation29,Citation31,Citation39,Citation40,Citation43,Citation58, Citation61,Citation73,Citation80] (see ). As for the heat transfer models, the proper choice of the boundary conditions depends on the modeling conditions. Nevertheless, for the clinical practice the Neumann BC is recommended to both BEM and BOS to produce an upper limit to the temperature distribution [Citation3].

Figure 4. Examples of applied boundary conditions for the calculations of the temperature distribution. The brown circles represent BEM with the electrodes given in the color gray, and black lines represent the BOS of the model. Arrows represent heat flux vector that is normal to the BEM (qBEM [W⋅m−2]) and the BOS (qBOS [W⋅m−2]). At both BEM and BOS (A) shows Dirichlet boundary conditions, and (B) shows Neumann boundary conditions if qBEM = 0 W⋅m−2 and qBOS = 0 W⋅m−2, or Robin boundary conditions if qBEM and qBOS are limited.

Figure 4. Examples of applied boundary conditions for the calculations of the temperature distribution. The brown circles represent BEM with the electrodes given in the color gray, and black lines represent the BOS of the model. Arrows represent heat flux vector that is normal to the BEM (qBEM [W⋅m−2]) and the BOS (qBOS [W⋅m−2]). At both BEM and BOS (A) shows Dirichlet boundary conditions, and (B) shows Neumann boundary conditions if qBEM = 0 W⋅m−2 and qBOS = 0 W⋅m−2, or Robin boundary conditions if qBEM and qBOS are limited.
3.2.2.3. Estimation of time durations of mild-hyperthermic and thermally ablative temperatures

Time durations of mild-hyperthermic (ΔtMH) and thermally ablative (ΔtTA) temperatures as a function of time were generally absent in the included studies. Nevertheless, these durations are relevant and the values can be estimated by determining the ratios between the time durations of these effects and the total treatment time (NP ⋅ (tP + τP)) [s]. Even though different model properties with various electrode shapes were utilized by the included studies, for simplicity the ratios were calculated by assuming a homogeneous, isotropic, and linear model within two infinite plate electrodes to obtain a homogeneous electric-field strength (see ). To further simplify the assumptions, the heat conduction, convection, and the metabolic heat generation were neglected within the plate electrodes, allowing the bioheat transfer equation to be simplified into (7) ΔT=J·E·tP·fPρcpΔt(7) where ΔT [°C] is the temperature increase during the time period Δt [s] and tP · fP is a dimensionless factor that averages the Joule Heating term over the entire intra- and interpulse duration (see Supplementary Appendix 3). For illustration, shows the temperature profile of EquationEquation 7 as a function of time over the period tend – tinit (or NP · (tP + τP)) with normothermic (from tinit to tMH), mild-hyperthermic (from tMH to tTA) and thermally ablative zones (from tTA to tend). According to this figure, the ratios can be calculated by (8) RΔtMH=ΔtMHNP(tP+τP)(8) and (9) RΔtTA=ΔtTANP(tP+τP)(9) where RΔtMH and RΔtTA are dimensionless ratios between mild-hyperthermic and thermally ablative time durations, and the total treatment time. Duo to the linearity of , (10) ΔtNP(tP+τP)=ΔTΔTmax (10)

Figure 5. (A) Assumed model that includes infinite electrode plates to calculate the temperature increase during IRE. The electrodes are presented as gray rectangles with infinite length and IRE-TR as a yellow rectangle. The black triangle represents a possible location in the model with maximum temperature increase. (B) An example of a temperature increase profile based on EquationEquation 7, which is divided in normothermia (from tinit to tMH), moderate hyperthermia (MH; from tMH to tTA), and thermal ablation (TA; from tTA to tend).

Figure 5. (A) Assumed model that includes infinite electrode plates to calculate the temperature increase during IRE. The electrodes are presented as gray rectangles with infinite length and IRE-TR as a yellow rectangle. The black triangle represents a possible location in the model with maximum temperature increase. (B) An example of a temperature increase profile based on EquationEquation 7(7) ΔT=J·E·tP·fPρcpΔt(7) , which is divided in normothermia (from tinit to tMH), moderate hyperthermia (MH; from tMH to tTA), and thermal ablation (TA; from tTA to tend).

After combining EquationEquations 8 and Equation9 with Equation10, the time durations can be replaced by the temperature differences as follows (11) RΔtMH=ΔTMHΔTmax (11) with (12) ΔTMHΔT310 3ΔT<13ΔT13(12) and (13) RΔtTA=ΔTTAΔTmax(13) with (14) ΔTTAΔT13  ΔT13(14)

Since the required data to calculate the total treatment time for each study are available, the calculation of RΔtMH and RΔtTA allows us to estimate the time durations of mild hyperthermia and thermal ablation using EquationEquations 8 and Equation9.

3.2.3. Thermal damage

Exposing biological cells or tissues to temperatures higher than their physiological temperatures for extended periods of time can result in thermal damage. In other words, irreversible changes in tissue structure or function that includes cell death, microvascular blood flow stasis, and/or protein coagulation processes [Citation40]. The thermal damage starts at the mild-hyperthermic temperature 42 °C–43 °C with exposure period between several seconds to hours, depending on the temperature level [Citation40,Citation53,Citation58]. In the case of relatively short exposure, damage is relatively limited up to 50 °C–60 °C, at which level the damage dramatically increases [Citation40,Citation58]. For instance, 73.4 °C is considered the threshold for instantaneous thermal damage in liver tissue that results in tissue whitening [Citation53]. Nevertheless, 50 °C could be generally assumed as a typical and relevant threshold one should attempt to avoid in the clinic as when temperatures around 50 °C occur during IRE the duration of exposure to these elevated temperatures could be in order of minutes [Citation15,Citation16], with a significant likelihood of resulting in thermal damage [Citation14,Citation23]. Multiple studies quantified the amount of thermal damage using the Arrhenius equation [Citation3,Citation26,Citation28,Citation30,Citation33,Citation34,Citation36, Citation40,Citation41,Citation53,Citation55,Citation58,Citation59,Citation80]. This equation is a first-order reaction rate constant that describes the rate of a chemical reaction. The Arrhenius equation, firstly applied by Henriques and Moritz in 1947 for skin burn analysis in pigs [Citation18,Citation82–84], uses Maxwell–Boltzmann statistics to describe the conversion of biological substances at temperature TK(t) from a native state to a thermally damaged state at a rate constant (15) Ωtt=A expUaṘTKt(15) where Ω(t) is the dimensionless time dependent accumulated thermal damage, A [s−1] is the pre-exponential factor (or collision frequency), a measure for molecular collision frequency, Ua [J⋅mol−1] is the activation energy required for the molecules to denature, Ṙ [J⋅mol−1⋅K−1] is the ideal gas constant (Ṙ = 8.314 [J⋅mol−1⋅K−1]), and TK(t) [K] is the time-dependent absolute temperature [Citation80]. Then, the thermal damage is calculated by integrating the Arrhenius Equation over a certain period of time, resulting in (16) Ωt=tinittA expUaṘTKtdt(16) where tinit [s] is the initial time point of the heating process [Citation3,Citation26,Citation28,Citation30,Citation33,Citation34,Citation36,Citation40,Citation41,Citation53,Citation55,Citation58,Citation59,Citation80]. The power of this analysis is to provide an estimation of the damaged fraction of a certain volume of a tissue constituent that is exposed to temperature TK. This estimation is calculated by hypothetically expressing both the intact and damaged substances in (17) Ωt=lnNtinitNt(17) where N(tinit) and N(t) represent the amount of the substances (expressed in number or concentration) in the tissue volume just before the temperature exposure and at time t [Citation33,Citation55,Citation59]. Then, the probability of thermal damage (the estimation of the damaged fraction), Υ [%], in the volume of the tissue constituent is calculated [Citation33,Citation53,Citation55,Citation59,Citation80] using (18) Υ=1001exp(Ωt(18)

The severeness of injuries due to the heating process can be represented by a threshold Ωth. For instance, it has been shown that the threshold for burn injuries in blood-perfused skin is Ωth = 0.53 (Υ = 41%) while it increases to Ωth = 10,000 (Υ = 100%) for third-degree burns [Citation40,Citation58]. In order to describe the thermal damage processes, such as cell death, microvascular blood flow stasis, or protein coagulation, suitable A, and Ua values must be chosen for each process, which are tissue specific and experimentally determined [Citation53]. Examples of such values can be found in Supplementary Appendix 6.

Another approach to describe the thermal damage is to calculate the cumulative equivalent minutes at 43 °C (CEM43 °C [min]; also defined as the thermal isoeffect dose t43 [min]) [Citation25,Citation26,Citation29,Citation31,Citation58,Citation76]. CEM43°C is a concept, introduced by Sapareto and Dewey [Citation19], which scales the exposure time of a time-varying temperature to an equivalent exposure time at the reference temperature 43 °C. Specifically, during the equivalent exposure time the tissue is assumed to be constantly held at 43 °C. CEM43°C can be described as (19) CEM43°C=tinittR43Tdtdt(19) where Tdt is the average temperature over the time differential dt, R is the factor to compensate for a 1 °C temperature change with R = 0.25 when Tdt ≤ 43 °C and R = 0.5 when Tdt >43°C. In this case, the severeness of injuries due to the heating process can be represented by a threshold CEM43°C(th). Examples of such values can be found in Supplementary Appendix 6 and [Citation85,Citation86].

The use of a suitable thermal damage model depends on both the temperature range to which the target had been exposed to, and the desired expression of the thermal damage outcome. E.g., CEM43°C can provide a reasonable prediction of cell death for the temperature range 40 °C–47 °C [Citation87]. However, it does not provide the probability of thermal damage as in the Arrhenius equation. On the other hand, the Arrhenius equation may overestimate the cell death at the start of the mild-hyperthermic temperature range [Citation88]. More information about both models can be found in [Citation19,Citation88,Citation89].

3.2.4. Tissue properties

The electric-field and temperature distributions of IRE have been calculated for multiple tissue types. For example, IRE has been applied to artery [Citation33,Citation34,Citation36,Citation41], brain [Citation31,Citation40,Citation57,Citation66, Citation69,Citation73], breast [Citation28,Citation29,Citation35], cervix [Citation78], eye [Citation44], in vitro cell suspension [Citation30,Citation43,Citation47,Citation60,Citation61,Citation75], kidney [Citation45,Citation63,Citation65], liver [Citation3,Citation24,Citation37–39,Citation46,Citation48,Citation53,Citation54,Citation62,Citation71,Citation72,Citation74,Citation76,Citation77,Citation79,Citation80], pancreas [Citation52,Citation70], prostate [Citation28,Citation32,Citation49,Citation50,Citation56,Citation64,Citation67,Citation68], rectal wall [Citation51], skin [Citation25,Citation58], small intestine [Citation59], and subcutaneous tissue [Citation42]. Although some studies used the same tissue type, often different tissue properties were applied, which were classified as homogeneous or heterogeneous, isotropic or anisotropic, and linear or nonlinear. The definitions of these terms can be found in Supplementary Appendix 2 Table A2.2. For instance, the tissue properties were generally assumed to be isotropic linear, e.g. σ, ρ, cp, k, Qm, wb, A, and Ua. Specifically, this was the case for ρ, cp, Qm, wb, A, and Ua in all of the included studies. However, some studies assumed isotropic nonlinear tissue properties, e.g. σ(E) [Citation42,Citation45,Citation53,Citation54,Citation56,Citation57,Citation59,Citation63,Citation66,Citation68–71, Citation77,Citation79], σ(T) [Citation3,Citation40,Citation43,Citation45,Citation60], σ(E,T) [Citation31,Citation37,Citation40,Citation45,Citation62,Citation73,Citation76,Citation80], σ(E, t) [Citation46,Citation78], and k(T) [Citation3]. In particular, most of the studies fitted the electrical conductivity in the form (20) σE=σinit+σmaxσinit11+a1expEa2a3(20) with (21) a2=EIREth+EREth2(21) and (22) a3=EIREthEREtha4(22) where σinit [S⋅m−1] is the initial electrical conductivity of the target medium before application of IRE, σmax [S⋅m−1] is the electrical conductivity of the maximally permeabilized medium, and a1 and a4 are the sigmoid function parameters. Additional detail about nonlinear electrical conductivities is provided in Supplementary Appendix 3.

The majority of the included studies applied a static (constant) electrical conductivity to calculate the electric field distribution. However, as shown by [Citation90,Citation91] for example, the electrical conductivity has a dynamic behavior due to its increase during the electroporation. By including a dynamic electrical conductivity, the numerical models can be better correlated to the experimental data, such as the electric current, as shown by [Citation63,Citation92]. In addition, the inclusion of the dynamic behavior can reduce the variability of EIRE(th) as was demonstrated by [Citation63] in canine kidneys (575 × 102 ± 67 V⋅m−1; versus 501 × 102 ± 121⋅102 V⋅m−1 obtained from the model with statistic electrical conductivity) resulting in a proper prediction of EIRE(th). An additional effect of the dynamically increasing electrical conductivity is the corresponding increase of the electric-field strength as was demonstrated by [Citation63,Citation68]. A larger electrical conductivity and electric-field strength result in a larger temperature increase (see EquationEquations 5 and Equation6). To determine the thermal effects and to properly estimate the possible lesion size of the treatment, it is important to include the dynamic behavior of the electrical conductivity in the treatment planning.

3.3. Synthesis of data, meta-analysis, and meta-regression

This section summarizes the extracted data from the included studies, focusing on the maximum simulated temperature increase, and the energy density that was calculated as (23) u=σinitVPD2tPNP(23) with u [J⋅m−3] as the energy density, and D [m] as the distance between electrodes (see ). Important: Please note that the u of an included study is not plotted if one of the extracted parameters in EquationEquation 23 is missing.

Figure 6. Examples of electrode configurations applied in studies to calculate the electric-field distribution, with the colors gray and brown as metallic and insulation parts, respectively. These configurations represent: (A) spherical electrodes, (B) needle and surface electrodes, (C) circular plate electrode enclosed by a ring electrode, (D, E) plate electrodes, (F) endovascular electrodes; contralateral electrodes have the same electric potential, (G) needle electrodes, (H) bipolar needle electrode; contralateral electrodes have different electric potential. In each configuration, the distances d [m] and D [m] were defined for the calculation of the energy density. Several studies assumed the needle and bipolar electrodes to be cylinders in 3D models. Other studies assumed the needle and spherical electrodes to be circles in 2D models.

Figure 6. Examples of electrode configurations applied in studies to calculate the electric-field distribution, with the colors gray and brown as metallic and insulation parts, respectively. These configurations represent: (A) spherical electrodes, (B) needle and surface electrodes, (C) circular plate electrode enclosed by a ring electrode, (D, E) plate electrodes, (F) endovascular electrodes; contralateral electrodes have the same electric potential, (G) needle electrodes, (H) bipolar needle electrode; contralateral electrodes have different electric potential. In each configuration, the distances d [m] and D [m] were defined for the calculation of the energy density. Several studies assumed the needle and bipolar electrodes to be cylinders in 3D models. Other studies assumed the needle and spherical electrodes to be circles in 2D models.

3.3.1. Organ-dependent electrical and thermal properties

The calculations of ΔTmax, R3ΔT13, RΔT13, and u are based on organ-dependent electrical and thermal characteristics, such as EIRE(th), σ, and other properties provided in EquationEquation 5. For these characteristics, a range of selected values was summarized in and plotted in both for cancerous and healthy tissues. Please note that kinit [W⋅m−1⋅°C−1] represents the initial thermal conductivity before the start of the heating process. To repeat, the electrical and thermal characteristics extracted from the included studies are provided in Supplementary Appendix 6 Table A6.2. Additional electrical and thermal characteristics can be found in the IT’IS database for thermal and electromagnetic parameters of biological tissues [Citation93].

Figure 7. Parameters of different tissue types, as selected by the included studies (Supplementary Appendix 6 – Table A6.2). For (A)–(E), red circles represent cancerous tissue and green diamonds represent healthy tissue. In all figures, brain tissue represents gray and white matter, skin tissue represents subcutaneous tissue, fat layer, dermis and stratum epidermis, eye tissue represents aqueous, cornea, lens, retina, sclera and vitreous, colon tissue represents both colon and rectal wall, blood vessel tissue represents artery, blood and blood vessel, and finally, muscle tissue represents muscles with longitudinal and transversal properties. The selected tissue properties in this figure show (A) the initial electrical conductivity (σinit [S⋅m−1]) before the application of IRE (the σinit values that exceeded 1 S⋅m−1 were 1.5 S⋅m−1 for aqueous and vitreous in the eye [Citation44], and 2.31 S⋅m−1 for cancerous breast tissue [Citation28]), (B) mass density (ρ [kg⋅m−3]), (C) specific heat capacity (cp [J⋅kg−1⋅°C−1]), (D) initial thermal conductivity before the start of heating process (kinit [W⋅m−1⋅°C−1]), (E) blood perfusion (wb [kg⋅m−3⋅s−1]; the wb value that exceeded 20 kg⋅m−3⋅s−1 was 212 kg⋅m−3⋅s−1 for both healthy and cancerous pancreas tissues [Citation70]; if wb could not be extracted nor calculated due to the absence of blood density (ρb [kg⋅m−3]), then wb was calculated by assuming ρb = 1060 kg⋅m−3), and (F) Electric-field threshold of irreversible electroporation (EIRE(th) [V⋅m−1]) as function of pulse number. Please note that the y-axis in (F) can be read in [V⋅cm−1], and that empty fields are present in (B)–(F) due to the lack of data.

Table 1. Summary of organ-dependent electrical and thermal characteristics for cancerous and healthy tissues (see Supplementary Appendix 6 – Table A6.2).

3.3.2. Thermal effects as function of IRE settings

In , we evaluated ΔTmax, and u obtained from the reported pulse parameters in the included studies, while in we only considered the more realistic subset of parameters that were actually applied in experiments. Furthermore, show that the pulses produce MEDIAN(ΔTmax) = 9.25 °C with the range 0.27 ≤ ΔTmax [°C] ≤ 135, and MEDIAN(u) ≈ 0.37 × 108 J⋅m−3 with the range 1.19 × 105 ≤ u [J⋅m−3] ≤ 11.42 × 108. In addition, shows that experimentally applied pulses produce MEDIAN(ΔTmax) = 9.49 °C with the range 1 ≤ ΔTmax [°C] ≤ 39, and MEDIAN(u) = 0.81 × 108 J⋅m−3 with the range 2.92 × 106 ≤ u [J⋅m−3] ≤ 8.22 × 108. In a similar way, we presented ΔTmax, and u in as in , with the exception of having ‘Total treatment time (NP ⋅ fP−1 [s])’ instead of ‘Total pulsing time’ with a MEDIAN(Total treatment time) = 80 s with the range of 12.5 ≤ Total treatment time [s] ≤ 540 (0.21 ≤ Total treatment time [min] ≤ 9).

Figure 8. (A) The study characteristics ‘Tissue type’, ‘ΔTmax’, ‘σinit’, ‘Electrode type’, and ‘u’ were plotted as function of ‘Pulse voltage (VP [V])’ and ‘Total pulsing time (tP ⋅ NP [s])’ in logarithmic scale. For a proper clarification of the data points, an extended version can be found in Supplementary Appendix 4 Figure A4.1. The range of ΔTmax was 0.27 ≤ ΔTmax [°C] ≤ 135, which was represented by the linearly scaled red-yellow color bar with the range 0 ≤ ΔTmax [°C] ≤ 40. The ΔTmax values that exceeded 40 °C were 50 °C at (750 V, 9 × 10−3 s) [Citation34], 55 °C at (1500 V, 1 × 10−2 s) [Citation64], 85 °C at (2000 V, 3 × 10−4 s) [Citation55], and 135 °C at (2000 V, 3 × 10−4 s) [Citation55]. Furthermore, the range of u was 1.19 × 105 ≤ u [J⋅m−3] ≤ 11.42 × 108, which was represented by the linearly scaled blue-magenta color bar with the range 0 ≤ u [J⋅m−3] ≤ 10 × 108. The u values that exceeded 10 × 108 J⋅m−3 were two times the value 11.42 × 108 J⋅m−3 at (2000 V, 9 × 10−3 s) [Citation44]. (B) Here, we only plotted the data that was obtained from the pulse parameters which were applied in the experiments. The range of ΔTmax was 1 ≤ ΔTmax [°C] ≤ 39, and the range of the energy density was 2.92 × 106 ≤ u [J⋅m−3] ≤ 8.22 × 108. (C) is similar to (B), with the exception of having ‘Total treatment time (NP ⋅ fP−1 [s])’ in logarithmic scale instead of ‘Total pulsing time’. Please note the change in Y-axis in (C) in comparison to (B), and that the electrical conductivity is the conductivity before the treatment.

Figure 8. (A) The study characteristics ‘Tissue type’, ‘ΔTmax’, ‘σinit’, ‘Electrode type’, and ‘u’ were plotted as function of ‘Pulse voltage (VP [V])’ and ‘Total pulsing time (tP ⋅ NP [s])’ in logarithmic scale. For a proper clarification of the data points, an extended version can be found in Supplementary Appendix 4 Figure A4.1. The range of ΔTmax was 0.27 ≤ ΔTmax [°C] ≤ 135, which was represented by the linearly scaled red-yellow color bar with the range 0 ≤ ΔTmax [°C] ≤ 40. The ΔTmax values that exceeded 40 °C were 50 °C at (750 V, 9 × 10−3 s) [Citation34], 55 °C at (1500 V, 1 × 10−2 s) [Citation64], 85 °C at (2000 V, 3 × 10−4 s) [Citation55], and 135 °C at (2000 V, 3 × 10−4 s) [Citation55]. Furthermore, the range of u was 1.19 × 105 ≤ u [J⋅m−3] ≤ 11.42 × 108, which was represented by the linearly scaled blue-magenta color bar with the range 0 ≤ u [J⋅m−3] ≤ 10 × 108. The u values that exceeded 10 × 108 J⋅m−3 were two times the value 11.42 × 108 J⋅m−3 at (2000 V, 9 × 10−3 s) [Citation44]. (B) Here, we only plotted the data that was obtained from the pulse parameters which were applied in the experiments. The range of ΔTmax was 1 ≤ ΔTmax [°C] ≤ 39, and the range of the energy density was 2.92 × 106 ≤ u [J⋅m−3] ≤ 8.22 × 108. (C) is similar to (B), with the exception of having ‘Total treatment time (NP ⋅ fP−1 [s])’ in logarithmic scale instead of ‘Total pulsing time’. Please note the change in Y-axis in (C) in comparison to (B), and that the electrical conductivity is the conductivity before the treatment.

Assumption: In case of multiple target tissues, e.g. the heterogeneous small intestine [Citation59], the tissue in the center of the volume was assumed to be the target.

3.3.3. Thermal effects as function of irreversible permeabilization effect

Furthermore, we presented the ratios between the mild-hyperthermic and thermally ablative effects, and IRE-TR, R3ΔT13 and RΔT13, as a function of SE-IRE(th) in . Again, here we included ΔTmax, and u of the experimentally applied pulse parameters. Despite the small number of data points, this figure shows that an increase in the IRE-TR results in a constant mild-hyperthermic effects with MEDIAN(R3ΔT13) = 31.75% and the range 0 ≤ R3ΔT13 [%] ≤ 88.67. In addition, it also results in the increase of the thermally ablative effects with MEDIAN(RΔT13) = 4.57% with the range of 0 ≤ RΔT13 [%] ≤ 41.3.

Figure 9. Linear regression plots of R3ΔT13 and RΔT13 data points as function of SE-IRE(th). For a proper clarification of the data points, an extended version can be found in Supplementary Appendix 4 Figure A4.2. The ΔTmax and the u ranges were 1.9 ≤ ΔTmax [°C] ≤ 39, and 0.29×107 ≤ u [J⋅m−3] ≤ 8.44 × 107. Please note that the scale of the color bar of the energy density in this figure is different than the scales in . In addition, the symbols at (4.83 × 10−5 m2, 0.32%) and (4.83 × 10−5 m2, 9.19%) have a gray background since the ΔTmax was not reported.

Figure 9. Linear regression plots of R3ΔT13 and RΔT13 data points as function of SE-IRE(th). For a proper clarification of the data points, an extended version can be found in Supplementary Appendix 4 Figure A4.2. The ΔTmax and the u ranges were 1.9 ≤ ΔTmax [°C] ≤ 39, and 0.29×107 ≤ u [J⋅m−3] ≤ 8.44 × 107. Please note that the scale of the color bar of the energy density in this figure is different than the scales in Figure 8. In addition, the symbols at (4.83 × 10−5 m2, 0.32%) and (4.83 × 10−5 m2, 9.19%) have a gray background since the ΔTmax was not reported.

3.3.4. Thermal effects as function of energy density

The data obtained from and Citation9 were used in to show the relationships between u and ΔTmax, and u, R3ΔT13 and RΔT13. These figures indicate that an increase in the u results in an increase of ΔTmax and R3ΔT13, while RΔT13 remains almost constant.

Figure 10. (A) Linear regression plot of ΔTmax data points as function of the energy density. (B) Linear regression plots of R3ΔT13 and RΔT13 data points as function of the energy density. Please note the change in X-axis in (B) in comparison to (A).

Figure 10. (A) Linear regression plot of ΔTmax data points as function of the energy density. (B) Linear regression plots of R3ΔT13 and RΔT13 data points as function of the energy density. Please note the change in X-axis in (B) in comparison to (A).

3.3.5. Estimation of time durations of thermal effects

illustrates the estimated time durations of both mild-hyperthermic and thermally ablative temperatures during IRE at the spatial point with the maximally simulated temperature increase ΔTmax. This was done for the studies included in meta-analysis with experimentally applied pulses. Please note that ΔtMH and ΔtTA were not calculated if they were provided by the included studies. Also, the data of ΔtMH and ΔtTA after the IRE treatment were included in only if they were provided by the included studies. To repeat, ΔtMH and ΔtTA represent the time durations of 3 ≤ ΔTmax [°C] < 13 (mild hyperthermia) and ΔTmax [°C] ≥ 13 (thermal ablation). The condition ‘ΔTmax [°C] ≥ 13’ at the spatial location with maximal simulated temperature increase (see ) indicates the presence of mild hyperthermia for a duration of ΔtMH prior to thermal ablation. Therefore, if the data from the included studies met the condition “ΔTmax [°C] ≥ 13”, then both values of ΔtMH and ΔtTA are plotted. shows that during the IRE treatment, the time duration for mild hyperthermia was MEDIAN(ΔtMH) = 15 s with 0 ≤ ΔtMH [s] ≤ 100, while the time duration for thermal ablation was MEDIAN(ΔtTA) = 0 s with 0 ≤ ΔtTA [s] ≤ 133.3. The median values include both calculated and extracted values. Furthermore, the studies from which ΔtMH and ΔtTA were extracted during and after IRE had medians of MEDIAN(ΔtMH) = 0 s with 0 ≤ ΔtMH [s] ≤ 180.7 and MEDIAN(ΔtTA) = 0 s 0 ≤ ΔtTA [s] ≤ 44.9.

Figure 11. The estimated time durations for both mild hyperthermia and thermal ablation as function ΔTmax. The term ‘Extracted’ means that ΔtMH and ΔtTA values were directly provided by the included studies, while the term ‘Calculated’ signifies that ΔtMH and ΔtTA were calculated using the assumptions stated in the Subsection 3.2.2.3. Estimation of time durations of mild-hyperthermic and thermally ablative temperatures. Please note that Δt represents the duration of mild-hyperthermic effect (3 °C ≤ ΔTmax < 13 °C) and thermal ablative effect (ΔTmax ≥ 13 °C) of an included study. The ΔTmax only represents the maximal temperature increase obtained from the included study. Also, note that the value of ΔtMH was 0 s for ΔTmax < 3 °C and the value of ΔtTA was 0 s for ΔTmax < 13 °C.

Figure 11. The estimated time durations for both mild hyperthermia and thermal ablation as function ΔTmax. The term ‘Extracted’ means that ΔtMH and ΔtTA values were directly provided by the included studies, while the term ‘Calculated’ signifies that ΔtMH and ΔtTA were calculated using the assumptions stated in the Subsection 3.2.2.3. Estimation of time durations of mild-hyperthermic and thermally ablative temperatures. Please note that Δt represents the duration of mild-hyperthermic effect (3 °C ≤ ΔTmax < 13 °C) and thermal ablative effect (ΔTmax ≥ 13 °C) of an included study. The ΔTmax only represents the maximal temperature increase obtained from the included study. Also, note that the value of ΔtMH was 0 s for ΔTmax < 3 °C and the value of ΔtTA was 0 s for ΔTmax < 13 °C.

4. Discussion

4.1. Research questions

According to the predictions of the numerical models, the relative sizes of both locally mild-hyperthermic and thermally ablative effects in the IRE-TR, with respect to the IRE-TR, have approximate medians of 30 and 5%, respectively. Accordingly, IRE produces sufficient heat to thermally ablate small parts of the IRE-TR, and to raise the temperature within the range of 40 ≤ T [°C] < 50 in a third of the IRE-TR. This temperature range probably reduces the EIRE(th) of the tissues and, consequently, enhances the IPE as was shown ex vivo for 43 °C in [Citation20]. Mouse ectopic pancreatic cancer treated using an integrated IRE system with controllable laser heating at 42 °C was compared with standard IRE without additional heating. A complete tumor regression was observed in more than 50% of the tumor-bearing mice, while in the standard IRE group no complete regression was observed [Citation19]. Results of this review predict that the IPE including the mild-hyperthermic effect dominates the IRE-TR with an approximate median of 95% (30% by combined contributions of mild-hyperthermic and IPE effects, and 65% by IPE effect only), implying that the IRE is based on the IPE with sufficient heat as adjuvant to enhance this effect. Again, it is important to notice that: first, the EIRE(th) values in this review corresponds to the EIRE(th) from the study from which SE-IRE(th),Σ was extracted, and second, the pulse parameters on which the values of ΔTmax, S3ΔT13,Σ and SΔT13,Σ depend, are different for each study.

The calculations of these medians were based on treated regions simulated by various types of setups. Specifically, the choice of the modeling parameters such as electrode type, electrical and thermal properties, and the type of BC can result in a different outcome of R3ΔT13 and RΔT13. Nevertheless, most of the studies used for the calculation of these medians attempt to mimic the real world by utilizing dynamic electrical conductivity and perfusion. In addition, a Neumann boundary condition was mostly used for the electrical and thermal outer boundary conditions to estimate the upper limit of the induced temperatures. Therefore, these medians are likely to represent the thermal effects on the safe side.

The maximum temperature increase in ) shows that the pulses produce MEDIAN(ΔTmax) = 9.25 °C with 0.27 ≤ ΔTmax [°C] ≤ 135, but the exact predicted maximum temperature depends on the pulse parameters and other factors assumed in the model. One would expect that an increase in the pulse voltage and the total pulsing time would result in an increase of ΔTmax. However, these figures do not show a clear correlation between pulse voltage, pulsing time and temperature increase, since the temperature increase also depends on the pulse frequency (higher pulse frequency results in larger temperature increase) and on the model configuration. An adjustment of these properties results in different ΔTmax. For instance, ΔTmax, including the ratios of the mild-hyperthermic and thermally ablative regions with respect to IRE-TR, increases as a result of the increase of the pulse voltage, the pulse duration, the pulse frequency, the electrical and the thermal conductivity of the tissue and the number of electrode pairs, and by assuming Neumann boundary conditions at the electrode-medium interface and the outer surface of the model. Besides this, other factors yielding a higher predicted maximum temperature are a decrease in the distance between the electrodes, a decrease in the specific heat capacity and density, and exclusion of the blood perfusion term. Example of variation in both electrical and thermal properties is shown in . These figures show differences between the values, even within the same tissue category.

The majority of the studies in ) did not actually apply the simulated pulses in experiments. In fact, some authors performed primarily a parametric study in which pulse, setup or tissue characteristics were varied, resulting in occasionally very large ΔTmax values [Citation34,Citation39,Citation44,Citation55]. Therefore, in ) we only included ΔTmax of the experimentally applied pulses with MEDIAN(ΔTmax) = 9.49 °C and the range 1 ≤ ΔTmax [°C] ≤ 39 for a median duration of 80 s (1.33 min) with (0.2 ≤ Total treatment time [min] ≤ 9), implying the presence of both mild-hyperthermic, and thermally ablative effects in the IRE-TR, which accords the findings in [Citation13–15]. Please note that the duration of these effects depends on the selected experiment. For example, according to the duration of mild hyperthermia during IRE has a median of 15 s with the range 0 ≤ Δt [s] ≤ 100, while the duration of thermal ablation has a median of 0 s with the range 0 ≤ Δt [s] ≤ 133.3. These durations were calculated at a single point where ΔTmax was obtained (only during the temperature increase). Accordingly, the durations of both effects can be longer when considering the decay of the temperature post-IRE treatment.

Nevertheless, ΔTmax does not specify the size of these effects. Accordingly, in we plotted the ratios of the mild-hyperthermic and thermally ablative effects with respect to IRE-TR as a function of the IRE-TR. As expected, the increase of the IRE-TR results in the increase of the ΔTmax, which implies that the increase of IRE-TR raises the ΔTmax, and therefore, enlarges thermally ablative effect.

4.2. Energy density

Concerning the energy density of the experimentally applied pulses in and considering EquationEquation 23 one would expect to find a positive linear relationship between the energy density and ΔTmax. Indeed, by using the linear regression analysis we do observe that an increase in the energy density results in an increase in ΔTmax. However, the data points do not show a clear correlation, because of the many different system configurations used as mentioned above. This also applies to , which shows the relationship between energy density and the ratios R3ΔT13 and RΔT13.

4.3. Selection of model parameters

According to , the selection of tissue properties varies between studies, even within the same tissue category. One reason for these variations is the inclusion of subcategories within the main categories. For example, ‘Brain’ can represent gray and white matter, ‘Blood vessel’ represents artery, blood and blood vessel, and ‘Muscle’ can represent longitudinal and transversal muscles in breast, prostate, and small intestine. Other reasons that can explain the variations are: (1) the anisotropicity and the heterogeneity of a tissue sample, (2) the location of the selected sample in an organ, (3) the animal species, and (4) the applied method to determine a property of a tissue type.

Due to the similarities of the values of the thermal properties in , one would expect the perfusion term wb to have a great influence on ΔTmax. Nevertheless, no correlation was found between these two parameters due to variation in the selection of other model properties, such as the boundary conditions and electrode type.

Finally, shows variation in the selection or calculation of EIRE(th) between studies, even within the same tissue categories. Reasons that explain these variations are: (1) the application of different pulse parameters, (2) the choice of static or dynamic electrical conductivity, (3) the choice between the modeling of homogeneous or heterogeneous tissue type, (4) the applied electrode configuration; in comparison to plate electrodes, we expect larger mild-hyperthermic region in case of needle electrodes due to higher electric field density, with as a result a larger treated region and lower EIRE(th).

Studies including the perfusion term in the bioheat transfer models applied constant blood perfusion. However, the blood perfusion is dynamic and depends on the magnitude and the duration of the thermal exposure, the tissue type, the size of the cancerous tissue, and whether a small or large animal model, or human model was used, as was shown in a review paper performed by Rossmann and Haemmerich [Citation94] that summarized on the temperature and tissue-dependent thermal properties, including the blood perfusion, of several biological tissues at mild-hyperthermic and ablative temperatures. Based on figure 6 of [Citation94], and since the medians of both simulation time and temperature increase of the included studies approximated 1.3 min (80 s) and 10 °C (47 °C relative to the assumed body temperature), we believe that temperature-dependent perfusion should be considered at T ≥ 43 °C for short total treatment time and at T ≥ 41 °C for long total treatment time. The latter would be preferred in treatment planning due to the presence of several needle pairs (more than 3 needle pairs in general; each pair with an approximation of 90 pulses and 10 test pulses) that may allow the extension of the total treatment time for more than 10 min. Nevertheless, the data about IRE-targeted organs are still limited. Therefore, further research should be performed to determine the influence of the magnitude and duration of thermal exposure.

4.4. Estimation of time durations of thermal effects

According to the extracted data and our simplified analysis, the mild-hyperthermic effect has a median time duration of 15 s, while thermally ablative effect has a median time duration of 0 s (ranges from 0 to ∼130 s) during the IRE treatment at the point with maximal temperature increase. Assuming a total treatment time of 80 s, this means that the combined medians of the time duration of the thermal effects would occur for approximately 20% of the treatment time (mild-hyperthermic effect for ∼20% and thermally ablative effect for 0%). However, these estimations relate to the duration during the IRE treatment itself and both effects will also last post IRE treatment. As a result, larger medians could be estimated for the total time durations. It is important to notice that the durations were calculated at the spatial point with ΔTmax. Since mild-hyperthermic temperatures spatially dominates the TA temperatures, the overall time duration of the mild-hyperthermic temperature, especially around the zones with T ≥ 50 °C (e.g. zones neighboring the electrodes), is longer than the time duration of TA effects. Furthermore, the studies included in the analysis used two electrodes in general. In the clinic, the number of electrode pairs is generally higher and could increase above 6 (more than 4 electrodes could be used), possibly resulting in further increase of ΔTmax and the mild-hyperthermic and ablative time durations. In particular, this applies for mild-hyperthermic effect, since an adequate order of the electrode pairs could avoid the extension of TA effects as was shown by O’Brien et al. [Citation95]. O’Brien et al. showed in an ex vivo study using a perfused porcine liver model that cycled pulsing schemes increase the treatment zone size and maintain a low tissue temperature compared to conventional pulsing schemes.

Our analysis can include estimation errors, since (1) the temperature profile has a temporal exponential profile if our simplified assumptions were not taken into consideration, and (2) only some of the simplified assumptions were assumed by the included studies. However, by comparing the extracted time durations with the calculated ones, a deviation of ≤22% was found. Finally, one must note that the included studies in which multiple electrode pairs were used did not consider the delays between the pulse sequences. In the clinic, delays are applied between the pulse sequences of different electrode pairs, both to recharge the pulse generator and to reduce the produced ΔTmax. By taking these delays into consideration, one would expect a reduction in thermal effects. Nevertheless, only a slight reduction of the combined median time duration of the thermal effects is expected since only two out of the 12 included studies applied IRE with 4 electrodes which involve activating four different electrode pairs during IRE [Citation33,Citation41].

4.5. Inclusion of vascular structures

The majority of the included studies did not incorporate vascular structures with blood vessels in the models [Citation3,Citation24,Citation26, Citation39–42]. Due to the high electrical conductivity of blood, the inclusion of vascular structures in or in close proximity of the treatment volume can result in distortion of the distribution of the electric field strength inside the treatment volume. Consequently, neighboring areas with relatively high or low electric field strength can emerge, resulting in undesirable overtreatment or undertreatment of those areas with possible clinical complications or treatment failure as a consequence. Several experimental and in silico studies showed that the distortion of the distribution of the electric field strength resulted in reduction of the electric field strength in areas inside the treated zone adjacent to the vascular structure; the ‘electric field sink’ effect [Citation54,Citation62,Citation71,Citation96]. For example, Goldberg et al. experimentally showed undertreated zones adjacent to vascular structures in the liver of rats after application of IRE between two electrode plates [Citation54]. We therefore recommend the inclusion of vascular structures in treatment planning to avoid possible undesired poor quality treatments and to accurately predict the local electric field and temperature distribution.

Furthermore, large blood vessels can yield a substantial cooling effect in the ablated region, similar as for RFA and MWA thermal ablation, which can result in a local temperature decrease of more than 20 °C [Citation97]. Several experimental and in silico studies showed a much lower temperature increase in the vicinity of a large blood vessel during RFA or MWA, dependent on the orientation of the vessel and its distance to the applicator [Citation97–100]. This demonstrates the relevance of accounting for large blood vessels also in IRE models, since a local cooling effect may significantly affect the produced thermal increase in the close vicinity of large vessels, which may result locally in only an irreversible permeabilization effect. Van Gemert et al. demonstrated in a 1D analytical IRE model that a large blood vessel at ∼5 mm distance from a needle can produce a local cooling effect up to a few degrees [Citation64]. A simplified model not including large vessels would be sufficient when for example evaluating maximum temperature levels achievable during IRE. However, when accurate predictions of the local temperature distribution are required, the effect of nearby large blood vessels should be accounted for in the models.

4.6. Challenges of simulations

The results presented in this review are based on calculations extracted from several studies with various simulated setups. Specifically, the extracted data may include errors, due to e.g. inaccurate modeling of the media (homogeneous instead of heterogeneous), inaccurate modeling of the simulated setup, inaccurate choice of boundary conditions, inaccurate assignment/modeling of electrical and thermal properties, and inaccurate choice of pulse characteristics (e.g. not modeling the delays within a single pulse sequence and between separate pulse sequences of different electrode pairs will result in overestimations of temperature). Another example is the modeling of σ(T), in which σ was modeled as a function of temperature (see Supplementary Appendix 4 for more details). Specifically, the electrical conductivity can increase by 1–3% per 1 °C [Citation50]. However, the increase of σ stops after a certain temperature due to the vaporization process of water in the medium. Therefore, a certain threshold should be included in σ(T) and σ(E, T) to avoid σ to further increase after a certain temperature. This temperature is typically 100 °C, which was applied by several RFA studies [Citation101,Citation102]. Despite these uncertainties in the modeling studies, the data obtained in this review still provides useful insight in the possible thermal elevation by IRE.

4.7. Take-home message

This review demonstrates the likelihood of direct thermal ablation in comparatively small regions and mild-hyperthermic effects in a large part of the IRE-TR, predicting the significant presence of hyperthermia in IRE. Even though the thermally ablated regions were relatively small, care must be taken for large temperature increases, especially in case of adjacent critical structures. This review focused on studies that performed numerical calculations, which are very useful to obtain insight in the general impact of various pulse parameters on the effect of IRE. Numerical calculations typically do not model nanopore effects and assume a simplified anatomy, not accounting for variations in tissue properties. Nevertheless, results clearly indicate that IRE can have significant thermal effects. Therefore, further in vivo research should investigate thermal effects of IRE to determine the relative size of the irreversible permeabilized region in which mild-hyperthermic and thermally ablative temperatures are achieved and significantly contribute to tumor ablation. These insights will help to optimize clinical protocols of IRE.

The expression ‘non-thermal IRE’ is generally used, (1) because IRE is a technique based on the creation of nanopores in the cell membrane, and (2) to distinguish IRE from thermal ablation techniques such as RFA and MWA. However, due to the significant presence of thermal effects (including 40 ≤ T [°C] < 50), and since literature data demonstrate that thermal effects clearly do contribute to the tumor ablation, one might consider avoiding the expression ‘non-thermal’ and start focusing on researching the mild-hyperthermic contribution to further optimize IRE-treatment.

Supplemental material

Supplemental Material

Download PDF (1.6 MB)

Supplemental Material

Download PDF (2.6 MB)

Disclosure statement

The author(s) declared the following potential conflicts of interest with respect to the research, authorship, and/or publication of this article: K. P. van Lienden, T. M. van Gulik and M.R. Meijerink 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.

Additional information

Funding

This research was financially supported by the Dutch Cancer Society [grant UvA 2014-7244].

References

  • Rubinsky B, Onik G, Mikus P. Irreversible electroporation: a new ablation modality – clinical implications. Technol Cancer Res Treat. 2007;6(1):37–48.
  • Maor E, Ivorra A, Rubinsky B. Non thermal irreversible electroporation: novel technology for vascular smooth muscle cells ablation. PloS One. 2009;4(3):e4757.
  • Davalos RV, Mir LM, Rubinsky B. Tissue ablation with irreversible electroporation. Ann Biomed Eng. 2005;33(2):223–231.
  • Mercadal B, Beitel-White N, Aycock KN, et al. Dynamics of cell death after conventional IRE and H-FIRE treatments. Ann Biomed Eng. 2020;8:1451–1462.
  • Pillai K, et al. Heat sink effect on tumor ablation characteristics as observed in monopolar radiofrequency, bipolar radiofrequency, and microwave, using ex vivo calf liver model. Medicine. 2015;94(9):e580.
  • Vogel JA, van Veldhuisen E, Agnass P, et al. Time-dependent impact of irreversible electroporation on pancreas, liver, blood vessels and nerves: a systematic review of experimental studies. PloS One. 2016;11(11):e0166987.
  • Scheltema MJV, van den Bos W, de Bruin DM, et al. Focal vs extended ablation in localized prostate cancer with irreversible electroporation; a multi-center randomized controlled trial. BMC Cancer. 2016;16(1):299.
  • Buijs M, van Lienden KP, Wagstaff PG, et al. Irreversible electroporation for the ablation of renal cell carcinoma: a prospective, human, in vivo study protocol (IDEAL Phase 2b). JMIR Res Protoc. 2017;6(2):e21.
  • Martin RCG, Kwon D, Chalikonda S, et al. Treatment of 200 locally advanced (stage III) pancreatic adenocarcinoma patients with irreversible electroporation safety and efficacy. Ann Surg. 2015;262(3):486–494.
  • Coelen RJS, Vogel JA, Vroomen LGPH, et al. Ablation with irreversible electroporation in patients with advanced perihilar cholangiocarcinoma (ALPACA): a multicentre phase I/II feasibility study protocol. BMJ Open. 2017;7(9):e015810.
  • Narayanan G, Hosein PJ, Beulaygue IC, et al. Percutaneous image-guided irreversible electroporation for the treatment of unresectable, locally advanced pancreatic adenocarcinoma. J Vasc Interv Radiol. 2017;28(3):342–348.
  • Vogel JA, Rombouts SJ, de Rooij T, et al. Induction chemotherapy followed by resection or irreversible electroporation in locally advanced pancreatic cancer (IMPALA): a prospective cohort study. Ann Surg Oncol. 2017;24(9):2734–2743.
  • van den Bos W, Scheffer HJ, Vogel JA, et al. Thermal energy during irreversible electroporation and the influence of different ablation parameters. J Vasc Interv Radiol. 2016;27(3):433–443.
  • Faroja M, Ahmed M, Appelbaum L, et al. Irreversible electroporation ablation: is all the damage nonthermal? Radiology. 2013;266(2):462–470.
  • Wagstaff PGK, et al. Irreversible electroporation of the porcine kidney: temperature development and distribution. Urol Oncol Semin Orig Investig. 2015;33(4):168.e1–7.
  • Agnass P, et al. Thermodynamic profiling during irreversible electroporation in porcine liver and pancreas: a case study series. J Clin Transl Res. 2020;5(3).
  • Bakker A, van der Zee J, van Tienhoven G, et al. Temperature and thermal dose during radiotherapy and hyperthermia for recurrent breast cancer are related to clinical outcome and thermal toxicity: a systematic review. Int J Hyperthermia. 2019;36(1):1024–1039.
  • Moritz AR, Henriques FC. Studies of thermal injury: II. The relative importance of time and surface temperature in the causation of cutaneous burns. Am J Pathol. 1947;23(5):695–720.
  • Sapareto SA, Dewey WC. Thermal dose determination in cancer-therapy. Int J Radiat Oncol Biol Phys. 1984;10(6):787–800.
  • Edelblute CM, Hornef J, Burcus NI, et al. Controllable moderate heating enhances the therapeutic efficacy of irreversible electroporation for pancreatic cancer. Sci Rep. 2017;7(1):11767.
  • Hutton B, Salanti G, Caldwell DM, et al. The PRISMA extension statement for reporting of systematic reviews incorporating network meta-analyses of health care interventions: checklist and explanations. Ann Intern Med. 2015;162(11):777–784.
  • Moher A, Liberati J, Tetzlaff DG. Altman D, The PRISMA Group. Preferred reporting items for systematic reviews and meta-analyses: the PRISMA statement. PLoS Med. 2009;6(7):e1000097.
  • Dewhirst MW, Viglianti BL, Lora-Michiels M, et al. Basic principles of thermal dosimetry and thermal thresholds for tissue damage from hyperthermia. Int J Hyperthermia. 2003;19(3):267–294.
  • Edd JF, Horowitz L, Davalos RV, et al. In vivo results of a new focal tissue ablation technique: irreversible electroporation. IEEE Trans Biomed Eng. 2006;53(7):1409–1415.
  • Al-Sakere B, André F, Bernat C, et al. Tumor ablation with irreversible electroporation. PloS One. 2007;2(11):e1135.
  • Edd JF, Davalos RV. Mathematical modeling of irreversible electroporation for treatment planning. Technol Cancer Res Treat. 2007;6(4):275–286.
  • Davalos RV, Rubinsky B. Temperature considerations during irreversible electroporation. Int J Heat Mass Transf. 2008;51(23-24):5617–5622.
  • Daniels C, Rubinsky B. Electrical field and temperature model of nonthermal irreversible electroporation in heterogeneous tissues. J Biomech Eng Trans ASME. 2009;131(7):071006.
  • Neal RE, Davalos RV. The feasibility of irreversible electroporation for the treatment of breast cancer and other heterogeneous systems. Ann Biomed Eng. 2009;37(12):2615–2625.
  • Shafiee H, Garcia PA, Davalos RV. A preliminary study to delineate irreversible electroporation from thermal damage using the Arrhenius equation. J Biomech Eng Trans ASME. 2009;131(7):074509.
  • Garcia PA, Rossmeisl JH, Neal RE, et al. Intracranial nonthermal irreversible electroporation: in vivo analysis. J Membrane Biol. 2010;236(1):127–136.
  • Golberg A, Rubinsky B. A statistical model for multidimensional irreversible electroporation cell death in tissue. BioMed Eng Online. 2010;9(1):13.
  • Maor E, Ivorra A, Mitchell JJ, et al. Vascular smooth muscle cells ablation with endovascular nonthermal irreversible electroporation. J Vasc Interv Radiol. 2010;21(11):1708–1715.
  • Maor E, Rubinsky B. Endovascular nonthermal irreversible electroporation: a finite element analysis. J Biomech Eng Trans ASME. 2010;132(3):031008.
  • Neal RE, Singh R, Hatcher HC, et al. Treatment of breast cancer through the application of irreversible electroporation using a novel minimally invasive single needle electrode. Breast Cancer Res Treat. 2010;123(1):295–301.
  • Phillips M, Maor E, Rubinsky B. Nonthermal irreversible electroporation for tissue decellularization. J Biomech Eng Trans ASME. 2010;132(9):091003.
  • Sano MB, Neal RE, Garcia PA, et al. Towards the creation of decellularized organ constructs using irreversible electroporation and active mechanical perfusion. BioMed Eng Online. 2010;9(1):83.
  • Zhang Y, Guo Y, Ragin AB, et al. MR imaging to assess immediate response to irreversible electroporation for targeted ablation of liver tissues: preclinical feasibility studies in a rodent model. Radiology. 2010;256(2):424–432.
  • Adeyanju OO, Al-Angari HM, Sahakian AV. The improvement of irreversible electroporation therapy using saline-irrigated electrodes: a theoretical study. Technol Cancer Res Treat. 2011;10(4):347–360.
  • Garcia PA, Rossmeisl JH, Neal RE, et al. A parametric study delineating irreversible electroporation from thermal damage based on a minimally invasive intracranial procedure. BioMed Eng Online. 2011;10(1):34.
  • Phillips M, Maor E, Rubinsky B. Principles of tissue engineering with nonthermal irreversible electroporation. J Heat Transf Trans ASME. 2011;133(1):011004.
  • Županič A, Miklavčič D. Tissue heating during tumor ablation with irreversible electroporation. Electrotechnical Rev. 2011;78(1-2):42–47.
  • Arena CB, Szot CS, Garcia PA, et al. A three-dimensional in vitro tumor platform for modeling therapeutic irreversible electroporation. Biophys J. 2012;103(9):2033–2042.
  • Mandel Y, Rubinsky B. Treatment of uveal melanoma by nonthermal irreversible electroporation: electrical and bioheat finite element model of the human eye. J Heat Transf Trans ASME. 2012;134(11):111101.
  • Neal RE, Garcia PA, Robertson JL, et al. Experimental characterization and numerical modeling of tissue electrical conductivity during pulsed electric fields for irreversible electroporation treatment planning. IEEE Trans Biomed Eng. 2012;59(4):1076–1085.
  • Sahakian AV, Al-Angari HM, Adeyanju OO. Electrode activation sequencing employing conductivity changes in irreversible electroporation tissue ablation. IEEE Trans Biomed Eng. 2012;59(3):604–607.
  • Kurata K, Ueno R, Matsushita M, et al. Experimental and analytical studies on contact irreversible electroporation for superficial tumor treatment. JBSE. 2013;8(4):306–318.
  • Mandel Y, Malki G, Adawi E, et al. Hemorrhage control of liver injury by short electrical pulses. PloS One. 2013;8(8):e49852.
  • Neal RE, Smith RL, Kavnoudias H, et al. The effects of metallic implants on electroporation therapies: feasibility of irreversible electroporation for brachytherapy salvage. Cardiovasc Intervent Radiol. 2013;36(6):1638–1645.
  • Qin ZP, Jiang J, Long G, et al. Irreversible electroporation: an in vivo study with dorsal skin fold chamber. Ann Biomed Eng. 2013;41(3):619–629.
  • Srimathveeravalli G, Wimmer T, Monette S, et al. Evaluation of an endorectal electrode for performing focused irreversible electroporation ablations in the swine rectum. J Vasc Interv Radiol. 2013;24(8):1249–1256.
  • Wimmer T, Srimathveeravalli G, Gutta N, et al. Comparison of simulation-based treatment planning with imaging and pathology outcomes for percutaneous CT-guided irreversible electroporation of the porcine pancreas: a pilot study. J Vasc Interv Radiol. 2013;24(11):1709–1718.
  • Garcia PA, Davalos RV, Miklavcic D. A numerical investigation of the electric and thermal cell kill distributions in electroporation-based therapies in tissue. PloS One. 2014;9(8):e103083.
  • Golberg A, Bruinsma BG, Uygun BE, et al. Tissue heterogeneity in structure and conductivity contribute to cell survival during irreversible electroporation ablation by electric field sinks. Sci Rep. 2015;5(1):8485.
  • Kurata K, Nomura S, Takamatsu H. Three-dimensional analysis of irreversible electroporation: estimation of thermal and non-thermal damage. Int J Heat Mass Transf. 2014;72:66–74.
  • Neal RE, Millar JL, Kavnoudias H, et al. In vivo characterization and numerical simulation of prostate properties for non-thermal irreversible electroporation ablation. Prostate. 2014;74(5):458–468.
  • Neal RE, Rossmeisl JH, D’Alfonso V, et al. In vitro and numerical support for combinatorial irreversible electroporation and electrochemotherapy glioma treatment. Ann Biomed Eng. 2014;42(3):475–487.
  • Nickfarjam A, Firoozabadi SMP. Parametric study of irreversible electroporation with different needle electrodes: electrical and thermal analysis. Int J Hyperthermia. 2014;30(5):335–347.
  • Phillips, M. The effect of small intestine heterogeneity on irreversible electroporation treatment planning. J Biomech Eng Trans ASME. 2014;136(9):091009.
  • Dermol J, Miklavcic D. Mathematical models describing Chinese hamster ovary cell death due to electroporation in vitro. J Membrane Biol. 2015;248(5):865–881.
  • Ivey JW, Latouche EL, Sano MB, et al. Targeted cellular ablation based on the morphology of malignant cells. Sci Rep. 2015;5(1):17157.
  • Kos B, Voigt P, Miklavcic D, et al. Careful treatment planning enables safe ablation of liver tumors adjacent to major blood vessels by percutaneous irreversible electroporation (IRE). Radiol Oncol. 2015;49(3):234–241.
  • Neal RE, Garcia PA, Kavnoudias H, et al. In vivo irreversible electroporation kidney ablation: experimentally correlated numerical models. IEEE Trans Biomed Eng. 2015;62(2):561–569.
  • van Gemert MJC, Wagstaff PGK, de Bruin DM, et al. Irreversible electroporation: just another form of thermal therapy? Prostate. 2015;75(3):332–335.
  • Wimmer T, Srimathveeravalli G, Gutta N, et al. Planning irreversible electroporation in the porcine kidney: are numerical simulations reliable for predicting empiric ablation outcomes? Cardiovasc Intervent Radiol. 2015;38(1):182–190.
  • Sharabi S, Kos B, Last D, et al. A statistical model describing combined irreversible electroporation and electroporation-induced blood-brain barrier disruption. Radiol Oncol. 2016;50(1):28–38.
  • Srimathveeravalli G, Cornelis F, Mashni J, et al. Comparison of ablation defect on MR imaging with computer simulation estimated treatment zone following irreversible electroporation of patient prostate. SpringerPlus. 2016;5(1):219.
  • Campelo S, et al. An evaluation of irreversible electroporation thresholds in human prostate cancer and potential correlations to physiological measurements. APL Bioeng. 2017;1(1):016101.
  • Garcia PA, Kos B, Rossmeisl JH, et al. Predictive therapeutic planning for irreversible electroporation treatment of spontaneous malignant glioma. Med Phys. 2017;44(9):4968–4980.
  • Latouche EL, Sano MB, Lorenzo MF, et al. Irreversible electroporation for the ablation of pancreatic malignancies: a patient-specific methodology. J Surg Oncol. 2017;115(6):711–717.
  • Qasrawi R, Silve L, Burdio F, et al. Anatomically realistic simulations of liver ablation by irreversible electroporation: impact of blood vessels on ablation volumes and undertreatment. Technol Cancer Res Treat. 2017;16(6):783–792.
  • Sung CK, Kim HB, Jung JH, et al. Histological and mathematical analysis of the irreversibly electroporated liver tissue. Technol Cancer Res Treat. 2017;16(4):488–496.
  • Wasson EM, Ivey JW, Verbridge SS, et al. The feasibility of enhancing susceptibility of glioblastoma cells to IRE using a calcium adjuvant. Ann Biomed Eng. 2017;45(11):2535–2547.
  • Yao CG, Dong SL, Zhao YJ, et al. Bipolar microsecond pulses and insulated needle electrodes for reducing muscle contractions during irreversible electroporation. IEEE Trans Biomed Eng. 2017;64(12):2924–2937.
  • Yoshimatsu S, Yoshida M, Kurata K, et al. Development of contact irreversible electroporation using a comb-shaped miniature electrode. J Therm Sci Technol. 2017;12(2):JTST0023–JTST0023.
  • Castellvi Q, Mercadal B, Moll X, et al. Avoiding neuromuscular stimulation in liver irreversible electroporation using radiofrequency electric fields. Phys Med Biol. 2018;63(3):035027.
  • Ritter A, Bruners P, Isfort P, et al. Electroporation of the liver: more than 2 concurrently active, curved electrodes allow new concepts for irreversible electroporation and electrochemotherapy. Technol Cancer Res Treat. 2018;17:153303381880999.
  • Yang YJ, Moser MAJ, Zhang E, et al. Development of a statistical model for cervical cancer cell death with irreversible electroporation in vitro. PloS One. 2018;13(4):e0195561.
  • Gallinato O, Denis de Senneville B, Seror O, et al. Numerical workflow of irreversible electroporation for deep-seated tumor. Phys Med Biol. 2019;64(5):055016.
  • O’Brien TJ, Bonakdar M, Bhonsle S, et al. Effects of internal electrode cooling on irreversible electroporation using a perfused organ model. Int J Hyperthermia. 2018;35(1):44–55.
  • Gallinato OPC. 2019. IRENA: a finite volume method based software for the numerical assessment of clinical IRE (version 1.0) APP IDDN.FR.001.080021.000.S.P.2019.000.31230 [Online]. Available from: https://team.inria.fr/monc/software/.
  • Henriques FC, Moritz AR. Studies of thermal injury: I. The conduction of heat to and through skin and the temperatures attained therein. A theoretical and an experimental investigation. Am J Pathol. 1947;23(4):530–549.
  • Moritz AR. Studies of thermal injury: III. The pathology and pathogenesis of cutaneous burns. An experimental study. Am J Pathol. 1947;23(6):915–941.
  • Henriques FC Jr. Studies of thermal injury: V. The predictability and the significance of thermally induced rate processes leading to irreversible epidermal injury. Arch Pathol. 1947;43(5):489–502.
  • Yarmolenko PS, Moon EJ, Landon C, et al. Thresholds for thermal damage to normal tissues: an update. Int J Hyperthermia. 2011;27(4):320–343.
  • van Rhoon GC, Samaras T, Yarmolenko PS, et al. CEM43A degrees C thermal dose thresholds: a potential guide for magnetic resonance radiofrequency exposure levels? Eur Radiol. 2013;23(8):2215–2227.
  • Mouratidis PXE, Rivens I, Civale J, et al. Relationship between thermal dose and cell death for “rapid” ablative and “slow” hyperthermic heating. Int J Hyperthermia. 2019;36(1):229–243.
  • Pearce J. Relationship between Arrhenius models of thermal damage and the CEM 43 thermal dose (SPIE BiOS). SPIE. 2009;7181:718104.
  • Pearce JA. Comparative analysis of mathematical models of cell death and thermal damage processes. Int J Hyperthermia. 2013;29(4):262–280.
  • Ivorra A, Rubinsky B. In vivo electrical impedance measurements during and after electroporation of rat liver. Bioelectrochemistry. 2007;70(2):287–295.
  • Pavlin M, Kandušer M, Reberšek M, et al. Effect of cell electroporation on the conductivity of a cell suspension. Biophys J. 2005;88(6):4378–4390.
  • Corovic S, Lackovic I, Sustaric P, et al. Modeling of electric field distribution in tissues during electroporation. BioMed Eng Online. 2013;12(1):16.
  • Hasgall PA, et al. 2018. IT’IS Database for thermal and electromagnetic parameters of biological tissues. Available from: https://itis.swiss/virtual-population/tissue-properties/database/database-summary/.
  • Rossmanna C, Haemmerich D. Review of temperature dependence of thermal properties, dielectric properties, and perfusion of biological tissues at hyperthermic and ablation temperatures. Crit Rev Biomed Eng. 2014;42(6):467–492.
  • J. O’Brien T, et al. Cycled pulsing to mitigate thermal damage for multi-electrode irreversible electroporation therapy. Int J Hyperthermia. 2019;36(1):953–963.
  • Marcan M, Kos B, Miklavcic D. Effect of blood vessel segmentation on the outcome of electroporation-based treatments of liver tumors. PloS One. 2015;10(5):e0125591.
  • Ringe KI, Lutat C, Rieder C, et al. Experimental evaluation of the heat sink effect in hepatic microwave ablation. PloS One. 2015;10(7):e0134301.
  • Zorbas G, Samaras T. A study of the sink effect by blood vessels in radiofrequency ablation. Comput Biol Med. 2015;57:182–186.
  • Jiang K, Chen J, Liu Y, et al. Heat-irrigate effect’ of radiofrequency ablation on relevant regional hepatocyte in living swine liver-initial study on pathology. Cell Biochem Biophys. 2015;72(1):37–41.
  • Chen RD, Lu F, Wu F, et al. An analytical solution for temperature distributions in hepatic radiofrequency ablation incorporating the heat-sink effect of large vessels. Phys Med Biol. 2018;63(23):235026.
  • Gonzalez-Suarez A, Berjano E. Comparative analysis of different methods of modeling the thermal effect of circulating blood flow during RF cardiac ablation. IEEE Trans Biomed Eng. 2016;63(2):250–259.
  • Trujillo M, Bon J, Berjano E. Computational modelling of internally cooled wet (ICW) electrodes for radiofrequency ablation: impact of rehydration, thermal convection and electrical conductivity. Int J Hyperthermia. 2017;33(6):624–634.