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

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.


Introduction
Irreversible electroporation (IRE) is a relatively new ablation modality for the treatment of unresectable cancers [1]. Contrary to contemporary ablation techniques such as radiofrequency ablation (RFA) and microwave ablation (MWA), IRE is proposed as a nonthermal ablation technique [2]. 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 [3]. 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 [4].
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 [5]. 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 [6]. Therefore, the application of IRE is widely investigated for several indications, such as renal and prostate tumors [7,8], pancreatic cancer [9][10][11][12], and perihilar cholangiocarcinoma [10].
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 [13][14][15][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 [3,[17][18][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 [20]. 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 mildhyperthermic 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?

Methods
This systematic review was performed according to the Preferred Reporting Items for Systematic Reviews and Meta (PRISMA)-analyses guidelines [21,22]. 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:

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.

Eligibility criteria
The focus of the studies is about IRE with rectangular pulses and a pulse length t P ! 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) t P <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.

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.

Definitions of mild-hyperthermic and thermally ablative effects
Assuming the tissues to have an initial physiological temperature (T init ) of 37 C and a thermally ablative threshold (T th ) of 50 C (DT ! 13 C), the following assumptions were made: The thermal elevation of 3 DT [ C] < 13 (40 T [ C] < 50) results in mild-hyperthermic effects.
The thermal elevation of DT ! 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 DT [ C] < 13 for the mild-hyperthermic effects and DT ! 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 T th . This because during IRE the duration of exposure to temperatures around 50 C could be in the order of minutes [15,16], which is sufficiently long that thermal damage may indeed occur [14,23] and this is thus effectively a clinically relevant threshold value associated with a high likelihood of occurrence of thermal ablation.

Analysis of research questions
To answer the first two research questions, we extracted the maximum simulated temperature increase (DT max [ 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 DT max was not reported. In case of parametric study, in which DT max was produced from several pulse parameters, the DT max 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 S 3DT13,R [m 2 ], S DT13,R [m 2 ], and S E-IRE(th),R [m 2 ] were extracted, where S 3DT13,R is the total surface area of the simulated temperature increase within the range 3 DT [ C] < 13 that represents the region with mild-hyperthermic temperature increase, S DT13,R is the total surface area of the simulated temperature increase with DT ! 13 C that represents the thermally ablated region, and S E-IRE(th),R is the total surface area of the simulated electric field with E ! E IRE(th) that represents the IRE-TR, E [VÁm À1 ] is the magnitude of the electric-field strength and E IRE(th) [VÁm À1 ] is the electric-field threshold of IRE; the minimum E required to ablate the target tissue. Even though E IRE(th) was only used as a minimum required E to extract S E-IRE(th),R for this analysis, it must be noticed that E IRE(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 E IRE(th) values correspond to the E IRE(th) from the study from which S E-IRE(th),R was extracted, the pulse parameters on which the values of DT max , S 3DT13,R and S DT13,R depend, are different for each study, in case of E IRE(th) range, e.g. [E IRE(th),1 , E IRE(th),2 ], the average of this range was used to calculate S E-IRE(th),R . Specifically, . if E IRE(th) was not reported in a certain study, then the S E-IRE(th),R was not calculated and, accordingly, included in the analysis.
The extraction of the surface areas was done excluding surface areas of the electrodes (S elect [m 2 ]). 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 DT-contours showing the conditions mentioned above, then we manually selected S E [m 2 ], S DT [m 2 ] and S elect using the 'Freehand selections' tool; S E is the selected surface area of the electric field strength and S DT 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 R 3DT13 [%] and R DT13 [%] were calculated, where R 3DT13 is the ratio between the region with mild-hyperthermic temperature increase and IRE-TR, and R DT13 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.

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 DT max . For the third research question, the medians were calculated of R 3DT13 , and R DT13 . 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 S E-IRE(th) , R 3DT13 and R DT13 , and the energy density (u [JÁm À3 ]; see Subsection 3.3.), DT max , R 3DT13 , and R DT13 .

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.

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 [25,30,50,64], two-dimensional [3,24, [47,75], and QuickField (Terra Analysis, Denmark) [54]. 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 Figure 2. Here, V P [V] is the electric potential (the voltage) of the pulses, t P [s] is the duration of a single pulse, n P [-] is the pulse number, N P [-] is the total number of pulses, s P [s] is the duration between two pulses, and f P [Hz] is the pulse frequency, which is calculated as 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 ere's law, into the steady-state form 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 Equation (2) and transform it into the Laplace's equation with r 2 [m À2 ] as Laplace operator.

Electric-field model and threshold. For both
Equations 2 and 3, the calculated electric-potential distribution is used to determine the electric-field distribution in the volume of the medium using ] as the electric field [26,43,68,72,74,79,80]. Assuming the volume of the medium to be an arbitrary tissue, the electric-field distribution yields IRE when a minimum electric-field strength E IRE(th) is obtained, providing the corresponding duration, number, shape, and the frequency of the pulses, and the tissue type. In Supplementary [79]. See Figure 3. 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 Figure 2. Overview of IRE pulses, with V P [V] as the electric potential (the voltage) of the pulses, t P [s] as the duration of a single pulse, n P as the pulse number, N P as the total number of pulses, s P [s] as the duration between two pulses, and f P [Hz] as the pulse frequency. practice, since the simulated electric field can provide an upper limit to the temperature distribution [3].

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 when IRE is used [3,24,[26][27][28][29]33,34,36,40,41,44,48,53,55,58,59,66,76]. 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.

Estimation of time durations of mild-hyperthermic
and thermally ablative temperatures. Time durations of mild-hyperthermic (Dt MH ) and thermally ablative (Dt TA ) 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 (N P Á (t P þ s 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 Figure 5(A)). 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  , which is divided in normothermia (from t init to t MH ), moderate hyperthermia (MH; from t MH to t TA ), and thermal ablation (TA; from t TA to t end ).
where DT [ C] is the temperature increase during the time period Dt [s] and t P Á f P is a dimensionless factor that averages the Joule Heating term over the entire intra-and interpulse duration (see Supplementary Appendix 3). For illustration, Figure 5(B) shows the temperature profile of Equation 7 as a function of time over the period t endt init (or N P Á (t P þ s P )) with normothermic (from t init to t MH ), mildhyperthermic (from t MH to t TA ) and thermally ablative zones (from t TA to t end ). According to this figure, the ratios can be calculated by and where RDt MH and RDt TA are dimensionless ratios between mild-hyperthermic and thermally ablative time durations, and the total treatment time. Duo to the linearity of Figure  5(B), After combining Equations 8 and 9 with 10, the time durations can be replaced by the temperature differences as follows and Since the required data to calculate the total treatment time for each study are available, the calculation of RDt MH and RDt TA allows us to estimate the time durations of mild hyperthermia and thermal ablation using Equations 8 and 9.

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 [40]. 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 [40,53,58]. In the case of relatively short exposure, damage is relatively limited up to 50 C-60 C, at which level the damage dramatically increases [40,58]. For instance, 73.4 C is considered the threshold for instantaneous thermal damage in liver tissue that results in tissue whitening [53]. 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 [15,16], with a significant likelihood of resulting in thermal damage [14,23]. Multiple studies quantified the amount of thermal damage using the Arrhenius equation [3,26,28,30,33,34,36,40,41,53,55,58,59,80]. 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 [18,[82][83][84], uses Maxwell-Boltzmann statistics to describe the conversion of biological substances at temperature T K (t) from a native state to a thermally damaged state at a rate constant where X(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, U a [JÁmol À1 ] is the activation energy required for the molecules to denature, _ R [JÁmol À1 ÁK À1 ] is the ideal gas constant ( _ R ¼ 8.314 [JÁmol À1 ÁK À1 ]), and T K (t) [K] is the timedependent absolute temperature [80]. Then, the thermal damage is calculated by integrating the Arrhenius Equation over a certain period of time, resulting in where t init [s] is the initial time point of the heating process [3,26,28,30,33,34,36,40,41,53,55,58,59,80]. 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 T K . This estimation is calculated by hypothetically expressing both the intact and damaged substances in where N(t init ) 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 [33,55,59]. Then, the probability of thermal damage (the estimation of the damaged fraction), Y [%], in the volume of the tissue constituent is calculated [33,53,55,59,80] using The severeness of injuries due to the heating process can be represented by a threshold X th . For instance, it has been shown that the threshold for burn injuries in blood-perfused skin is X th ¼ 0.53 (Y ¼ 41%) while it increases to X th ¼ 10,000 (Y ¼ 100%) for third-degree burns [40,58]. In order to describe the thermal damage processes, such as cell death, microvascular blood flow stasis, or protein coagulation, suitable A, and U a values must be chosen for each process, which are tissue specific and experimentally determined [53].
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 t 43 [min]) [25,26,29,31,58,76]. CEM43 C is a concept, introduced by Sapareto and Dewey [19], 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 where T dt 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 T dt 43 C and R ¼ 0.5 when T dt >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 [85,86]. 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 [87]. 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 [88]. More information about both models can be found in [19,88,89].

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 [33,34,36,41], brain [31,40,57,66,69,73], breast [28,29,35], cervix [78], eye [44], in vitro cell suspension [30,43,47,60,61,75], kidney [45,63,65], liver [3,24,[37][38][39]46,48,53,54,62,71,72,74,76,77,79,80], pancreas [52,70], prostate [28,32,49,50,56,64,67,68], rectal wall [51], skin [25,58], small intestine [59], and subcutaneous tissue [42]. 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. r, q, c p , k, Q m , w b , A, and U a . Specifically, this was the case for q, c p , Q m , w b , A, and U a in all of the included studies. However, some studies assumed isotropic nonlinear tissue properties, e.g. r(E) [42,45,53,54,56,57,59,63,66,68-71, 77,79], r(T) [3,40,43,45,60], r(E,T) [31,37,40,45,62,73,76,80], r(E, t) [46,78], and k(T) [3]. In particular, most of the studies fitted the electrical conductivity in the form  (22) where r init [SÁm À1 ] is the initial electrical conductivity of the target medium before application of IRE, r max [SÁm À1 ] is the electrical conductivity of the maximally permeabilized medium, and a 1 and a 4 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 [90,91] 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 [63,92]. In addition, the inclusion of the dynamic behavior can reduce the variability of E IRE(th) as was demonstrated by [63] in canine kidneys (575 Â 10 2 ± 67 VÁm À1 ; versus 501 Â 10 2 ± 121Á10 2 VÁm À1 obtained from the model with statistic electrical conductivity) resulting in a proper prediction of E IRE(th) . An additional effect of the dynamically increasing electrical conductivity is the corresponding increase of the electric-field strength as was demonstrated by [63,68]. A larger electrical conductivity and electric-field strength result in a larger temperature increase (see Equations 5 and 6). 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.

Synthesis of data, meta-analysis, and metaregression
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 with u [JÁm À3 ] as the energy density, and D [m] as the distance between electrodes (see Figure 6). Important: Please note that the u of an included study is not plotted if one of the extracted parameters in Equation 23 is missing.

Organ-dependent electrical and thermal properties
The calculations of DT max , R 3DT13 , R DT13 , and u are based on organ-dependent electrical and thermal characteristics, such as E IRE(th) , r, and other properties provided in Equation 5. For these characteristics, a range of selected values was summarized in Table 1 and plotted in Figure 7 both for cancerous and healthy tissues. Please note that k init [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

Thermal effects as function of IRE settings
In Figure 8(A), we evaluated DT max , and u obtained from the reported pulse parameters in the included studies, while in  Please note that the brackets ([min, max]) indicate the range between and including the minimum and maximal selected values. Also note that the rare instances of values that exceeded 1 SÁm À1 for the initial electrical conductivity (r init [SÁm À1 ]) were 1.5 SÁm À1 for aqueous and vitreous in the eye [44], and 2.31 SÁm À1 for cancerous breast tissue [28]. In case of blood perfusion (w b [kgÁm À3 Ás À1 ]), the value that exceeded 20 kgÁm À3 Ás À1 was 212 kgÁm À3 Ás À1 for both healthy and cancerous pancreas tissues [70]. If w b could not be extracted nor calculated due to the absence of blood density (q b [kgÁm À3 ]), then w b was calculated by assuming q b ¼ 1060 kgÁm À3 . a 2.31 SÁm À1 was reported as the maximum value for cancerous breast tissue [28]. b 1.5 SÁm À1 was reported as the maximum value for aqueous and vitreous in the eye [44]. c 212 kgÁm À3 Ás À1 was reported as the maximum value for both healthy and cancerous pancreas tissues [70].  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 (r init [SÁm À1 ]) before the application of IRE (the r init values that exceeded 1 SÁm À1 were 1.5 SÁm À1 for aqueous and vitreous in the eye [44], and 2.31 SÁm À1 for cancerous breast tissue [28]), (B) mass density (q [kgÁm À3 ]), (C) specific heat capacity (c p [JÁkg À1 Á C À1 ]), (D) initial thermal conductivity before the start of heating process (k init [WÁm À1 Á C À1 ]), (E) blood perfusion (w b [kgÁm À3 Ás À1 ]; the w b value that exceeded 20 kgÁm À3 Ás À1 was 212 kgÁm À3 Ás À1 for both healthy and cancerous pancreas tissues [70]; if w b could not be extracted nor calculated due to the absence of blood density (q b [kgÁm À3 ]), then w b was calculated by assuming q b ¼ 1060 kgÁm À3 ), and (F) Electric-field threshold of irreversible electroporation (E IRE(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.  Assumption: In case of multiple target tissues, e.g. the heterogeneous small intestine [59], the tissue in the center of the volume was assumed to be the target.

Thermal effects as function of irreversible permeabilization effect
Furthermore, we presented the ratios between the mildhyperthermic and thermally ablative effects, and IRE-TR, R 3DT13 and R DT13 , as a function of S E-IRE(th) in Figure 9.

Thermal effects as function of energy density
The data obtained from Figures 8(C) and 9 were used in Figure 10 to show the relationships between u and DT max , and u, R 3DT13 and R DT13 . These figures indicate that an increase in the u results in an increase of DT max and R 3DT13 , while R DT13 remains almost constant. Figure 11 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 DT max . This was done for the studies included in meta-analysis with experimentally applied pulses. Please note that Dt MH and Dt TA were not calculated if they were provided by the included studies. Also, the data of Dt MH and Dt TA after the IRE treatment were included in Figure 11 only if they were provided by the included studies.

Estimation of time durations of thermal effects
To repeat, Dt MH and Dt TA represent the time durations of 3 DT max [ C] < 13 (mild hyperthermia) and DT max [ C] ! 13 (thermal ablation). The condition 'DT max [ C] ! 13' at the spatial location with maximal simulated temperature increase (see Figure 5) indicates the presence of mild hyperthermia for a duration of Dt MH prior to thermal ablation. Therefore, if the data from the included studies met the condition "DT max [ C] ! 13", then both values of Dt MH and Dt TA are plotted. Figure 11 shows that during the IRE treatment, the time duration for mild hyperthermia was MEDIAN(

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 E IRE(th) of the tissues and, consequently, enhances the IPE as was shown ex vivo for 43 C in [20]. 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 [19]. 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 E IRE(th) values in this review corresponds to the E IRE(th) from the study from which S E-IRE(th),R was extracted, and second, the pulse parameters on which the values of DT max , S 3DT13,R and S DT13,R 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 R 3DT13 and R DT13 . 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 Figure 8(A) shows that the pulses produce MEDIAN(DT max ) ¼ 9.25 C with 0.27 DT max [ C] 135, but the exact predicted maximum temperature depends on the pulse parameters and other DT max [ C] 135, which was represented by the linearly scaled red-yellow color bar with the range 0 DT max [ C] 40. The DT max values that exceeded 40 C were 50 C at (750 V, 9 Â 10 À3 s) [34], 55 C at (1500 V, 1 Â 10 À2 s) [64], 85 C at (2000 V, 3 Â 10 À4 s) [55], and 135 C at (2000 V, 3 Â 10 À4 s) [55]. Furthermore, the range of u was 1.19 Â 10 5 u [JÁm À3 ] 11.42 Â 10 8 , which was represented by the linearly scaled blue-magenta color bar with the range 0 u [JÁm À3 ] 10 Â 10 8 . The u values that exceeded 10 Â 10 8 JÁm À3 were two times the value 11.42 Â 10 8 JÁm À3 at (2000 V, 9 Â 10 À3 s) [44]. (B) Here, we only plotted the data that was obtained from the pulse parameters which were applied in the experiments. The range of DT max was 1 DT max [ C] 39, and the range of the energy density was 2.92 Â 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 DT max . 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 DT max . For instance, DT max , 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 Figure  7(A-E). These figures show differences between the values, even within the same tissue category.
The majority of the studies in Figure 8(A) 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 DT max values [34,39,44,55]. Therefore, in Figure 8(B,C) we only included DT max of the experimentally applied pulses with MEDIAN(DT max ) ¼ 9.49 C and the range 1 DT max [ 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 [13][14][15]. Please note that the duration of these effects depends on the selected experiment. For example, according to Figure  Nevertheless, DT max does not specify the size of these effects. Accordingly, in Figure 9 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 DT max , which implies that the increase of IRE-TR raises the DT max , and therefore, enlarges thermally ablative effect.

Energy density
Concerning the energy density of the experimentally applied pulses in Figure 10(A) and considering Equation 23 one would expect to find a positive linear relationship between the energy density and DT max . Indeed, by using the linear regression analysis we do observe that an increase in the energy density results in an increase in DT max . 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 Figure 10(B), which shows the relationship between energy density and the ratios R 3DT13 and R DT13 .

Selection of model parameters
According to Figure 7(A-E), 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 Figure 7(B-D), one would expect the perfusion term w b to have a great influence on DT max . 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, Figure 7(F) shows variation in the selection or calculation of E IRE(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 E IRE(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 [94] 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 [94], 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.

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 DT max . 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 DT max 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. [95]. 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 DT max . By Figure 11. The estimated time durations for both mild hyperthermia and thermal ablation as function DT max . The term 'Extracted' means that Dt MH and Dt TA values were directly provided by the included studies, while the term 'Calculated' signifies that Dt MH and Dt TA 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 Dt represents the duration of mildhyperthermic effect (3 C DT max < 13 C) and thermal ablative effect (DT max ! 13 C) of an included study. The DT max only represents the maximal temperature increase obtained from the included study. Also, note that the value of Dt MH was 0 s for DT max < 3 C and the value of Dt TA was 0 s for DT max < 13 C.
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 [33,41].

Inclusion of vascular structures
The majority of the included studies did not incorporate vascular structures with blood vessels in the models [3,24,26,[39][40][41][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 [54,62,71,96]. 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 [54]. 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 [97]. 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 [97][98][99][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 [64]. 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.

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 r(T), in which r 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 [50]. However, the increase of r stops after a certain temperature due to the vaporization process of water in the medium. Therefore, a certain threshold should be included in r(T) and r(E, T) to avoid r to further increase after a certain temperature. This temperature is typically 100 C, which was applied by several RFA studies [101,102]. Despite these uncertainties in the modeling studies, the data obtained in this review still provides useful insight in the possible thermal elevation by IRE.

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 'nonthermal' and start focusing on researching the mild-hyperthermic contribution to further optimize IRE-treatment.

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.