Computational Modeling of Thermal Ablation Zones in the Liver: A Systematic Review

Simple Summary Thermal ablation is an established treatment for primary and secondary liver tumors. As ablation treatment planning is a fast-emerging field, accurate and patient-specific ablation zone simulation may contribute to higher efficacy of thermal ablation. Computational modeling could facilitate these simulations. This systematic review aims to identify, evaluate, and summarize the findings of the literature on existing computational models for thermal liver ablation planning and compare their accuracy. The literature shows a wide variety of computational modeling and validation methods. Additional research, with a focus on shape-based outcome metrics, is warranted to determine which model demonstrates superior accuracy and suitability for clinical practice. More insight into parameter personalization is required to enable patient-specific ablation planning. Abstract Purpose: This systematic review aims to identify, evaluate, and summarize the findings of the literature on existing computational models for radiofrequency and microwave thermal liver ablation planning and compare their accuracy. Methods: A systematic literature search was performed in the MEDLINE and Web of Science databases. Characteristics of the computational model and validation method of the included articles were retrieved. Results: The literature search identified 780 articles, of which 35 were included. A total of 19 articles focused on simulating radiofrequency ablation (RFA) zones, and 16 focused on microwave ablation (MWA) zones. Out of the 16 articles simulating MWA, only 2 used in vivo experiments to validate their simulations. Out of the 19 articles simulating RFA, 10 articles used in vivo validation. Dice similarity coefficients describing the overlap between in vivo experiments and simulated RFA zones varied between 0.418 and 0.728, with mean surface deviations varying between 1.1 mm and 8.67 mm. Conclusion: Computational models to simulate ablation zones of MWA and RFA show considerable heterogeneity in model type and validation methods. It is currently unknown which model is most accurate and best suitable for use in clinical practice.


Introduction
Percutaneous thermal ablation is an established minimally invasive treatment for primary and secondary liver tumors [1,2].Radiofrequency ablation (RFA) and microwave ablation (MWA) are currently the most widely applied thermal ablation techniques to treat liver malignancies.Both techniques aim to induce tissue heating of at least 55-60 • C to necrotize the tumor along with an ablative margin of normal liver parenchyma around the tumor of at least 5 mm [3,4].RFA applies a rapidly alternating current that excites the ions in the liver tissue, causing frictional heating.In MWA, electromagnetic waves cause polar Cancers 2023, 15, 5684.https://doi.org/10.3390/cancers15235684https://www.mdpi.com/journal/cancersmolecules, predominantly water, to realign with the oscillating field, which generates heat through kinetic energy [5].
The obtained ablative margin in thermal ablation is correlated with local recurrence rates.In a study by Laimer et al. conducted on patients with hepatocellular carcinoma, each millimeter increase in the minimal ablative margin resulted in a 30% risk reduction for local recurrence [4].No recurrences occurred when an ablative margin of >5 mm was obtained, but this was only achieved in 37.5% of the ablations.These results are in accordance with several other studies investigating the correlation between ablative margin and local recurrences in primary and secondary liver tumors [3,4,[6][7][8][9].In these studies, the percentage of the intended ablative margin of >5 mm varied between 2.7% and 51.4%.These low rates indicate a discrepancy between the predicted and created ablation zone, either due to shape and size deviation of the created ablation zone or inaccurate needle positioning.
Currently, the ablation zone size is predicted according to the manufacturer's specifications.With each system, a chart is provided containing 2D ellipse predictions for several settings (such as ablation time and power).In recent years, treatment planning tools have become available that use these predictions to help the operator choose the probe trajectory and ablation setting [10].Yet, these predictions are mostly based on preclinical animal experiments [11].Various factors, such as vascular proximity, tissue perfusion, tumor location, and underlying fibrosis or cirrhosis, can lead to differences in the shape and volume of the ablation zone compared to predictions.Computational modeling and ex vivo experiments have demonstrated that the aforementioned tumor and liver characteristics affect the heat conductivity and thus the dimensions of the actual ablation zone [12][13][14][15].To increase the rate of adequate ablations with sufficient margins and thereby reduce the risk of local recurrence, patient-specific therapy planning is essential.
The computation of ablation necrosis volumes generally consists of two steps.First, increase in temperature during thermal ablation can be simulated by computational models based on electromagnetic and bioheat equations [16,17].Cell-death models then convert these temperatures to a volume of necrosis.A wide variety of complexity is observed in these models [16][17][18].Whereas some computational models rely only on a bioheat equation with fixed parameters, other models incorporate more complexities, such as patient-specific anatomy and temperature-dependent parameters.
Computational models are currently mainly used for in silico testing and optimization of devices as well as to study the effect of different parameters and techniques on ablation treatment [19][20][21].As treatment planning is a fast-emerging field, accurate and patientspecific ablation zone simulation may contribute to higher efficacy of thermal ablation [10].Computational modeling potentially opens up an opportunity for patient-specific ablation zone simulation.This systematic review aims to identify, evaluate, and summarize the findings of the literature on existing computational models for radiofrequency and microwave thermal liver ablation planning and compare their accuracy.

Search Strategy
Studies were identified by searching the electronic databases MEDLINE and Web of Science on 17 April 2023.The search queries were based on synonyms of the keywords "Thermal ablation", "Liver neoplasm", and "Computational modeling".The complete search strategies used can be found in Appendix A. The study has not been registered in PROSPERO.

Study Selection
After duplicate removal, abstracts were screened, followed by full-text assessment.Articles were found eligible if (i) ablation zone simulation was performed for either (ii) percutaneous RFA or MWA in (iii) liver tissue and if (iv) the model was quantitatively validated using ex vivo or in vivo experiments with (v) either ablation zone dimensions, ablation zone volumes, or volume-based metrics reported as outcome measure(s).Furthermore, articles should focus on simulating ablation zones for treatment planning and not for device or protocol optimization.Reviews, systematic reviews, letters to the editor, and articles written in other languages than English were excluded.Two authors (G.C.M.v.E.and P.H.) independently assessed the articles according to these criteria.In case of disagreement, consensus was reached by discussion.

Data Extraction
Included articles were categorized based on the thermal ablation technique used, either RFA or MWA.For each included computational model, the applied biological heat transfer model and the cell death model were extracted.Furthermore, it was noted whether perfusion, blood vessels, water vaporization, temperature-dependent thermal parameters, and/or an image-based anatomical model were incorporated.Data extraction for the validation of computational models included study type (in vivo, ex vivo, or clinical studies), number of ablations, ablation settings, ground truth comparison, outcome metrics, and validation results.The relative volume deviation (RVD) or relative diameter deviation (RDD) between the simulated and experimentally obtained ablation zones was calculated based on reported ablation zone volumes or diameters, if this outcome was not reported on already, to ensure homogenous size-based outcome metrics.

Study Selection
The search strategy identified 849 articles after removal of duplicates.A total of 673 articles were excluded after abstract screening.A total of 176 articles were full-text assessed, resulting in the exclusion of 141 articles.A total of 35 articles met the inclusion criteria and were included in this systematic review .Figure 1  Out of the 35 articles included, 16 focused on simulating MWA and 19 on RFA.Tables 1 and 2 contain details about the computational model used in these articles.Figure 2 gives a schematic overview of a general model structure.Frequently used equations in the models can be found in Appendix B and frequently used outcome metrics in Appendix C. Out of the 35 articles included, 16 focused on simulating MWA and 19 on RFA.Tables 1 and 2 contain details about the computational model used in these articles.Figure 2 gives a schematic overview of a general model structure.Frequently used equations in the models can be found in Appendix B and frequently used outcome metrics in Appendix C. Out of the 35 articles included, 16 focused on simulating MWA and 19 on RFA.Tables 1 and 2 contain details about the computational model used in these articles.Figure 2 gives a schematic overview of a general model structure.Frequently used equations in the models can be found in Appendix B and frequently used outcome metrics in Appendix C.      Cancers 2023, 15, 5684 9 of 28

MWA Data Analysis
Sixteen articles presented computational models for MWA.Thirteen articles used Pennes' bioheat equation as bioheat model [26,29,30,[35][36][37][38]44,46,52,53,55,56], while one article used the transient heat transfer equation [33], one the local thermal non-equilibrium equation [48], and one used a proprietary heat transfer model [31].As a cell death model, the Arrhenius thermal damage model was used in four articles, the three-state cell death equation was used in two articles, and an isothermal contour (54 • C or 60 • C) was used in eight articles.Two articles used both the Arrhenius thermal damage model and an isothermal contour of 52 • C and 54 • C, respectively [30,53].
Seven articles included perfusion in their model, two blood vessels, thirteen water vaporization, and thirteen temperature-dependent tissue parameters.Two models used CT-based anatomy models.Gao

MWA Ex Vivo Validation
Table 3 gives an overview of the ex vivo validation of the included computational models for MWA simulation.Fourteen articles used ex vivo experiments in animals or phantoms [26,[29][30][31]33,[35][36][37][38]44,46,52,53,55].Two articles used the Dice similarity coefficient (DSC) to express their results and found similar scores between 0.74 and 0.82 [31,33].One article used the Jaccard similarity index and found results of 0.866 and 0.934.However, these results might be biased since the electrical and thermal conductivities were reconstructed after the experiments to best fit the model [29].Sing et al. used the experiments of Wu et al. to validate their simulated ablation zone [44,55]

MWA Data Analysis
Sixteen articles presented computational models for MWA.Thirteen articles used Pennes' bioheat equation as bioheat model [26,29,30,[35][36][37][38]44,46,52,53,55,56], while one article used the transient heat transfer equation [33], one the local thermal non-equilibrium equation [48], and one used a proprietary heat transfer model [31].As a cell death model, the Arrhenius thermal damage model was used in four articles, the three-state cell death equation was used in two articles, and an isothermal contour (54 °C or 60 °C) was used in eight articles.Two articles used both the Arrhenius thermal damage model and an isothermal contour of 52 °C and 54 °C, respectively [30,53].
Seven articles included perfusion in their model, two blood vessels, thirteen water vaporization, and thirteen temperature-dependent tissue parameters.Two models used CT-based anatomy models.Gao

MWA In Vivo Validation
Table 4 presents an overview of the two articles that used in vivo validation in patients [48,56].Tucci et al. modeled four different blood vessels and compared them to the in vivo experiments of Amabile et al. [48,58].They concluded that their model, including terminal arteries, showed a good agreement with the ablation zones achieved in the clinical study.Zhai et al. performed a study on nine patients [56].Ablation simulation had an RVD of ±7.0% compared to clinically obtained ablation volumes.However, the study has a limited sample size without uniform image analysis.

RFA Data Analysis
Nineteen articles presented computational models for RFA.Thirteen articles used Pennes' bioheat equation as a bioheat model [25,27,28,32,34,[39][40][41][42]45,[49][50][51], while one article used the heat transfer equation [54], one the split volume bioheat equation [43], and one article compared three different bioheat-models, i.e., Pennes' bioheat equation, the local thermal equilibrium equation, and the local thermal non-equilibrium equation [47].Three articles by Audigier et al. used a combination of Pennes' bioheat equation with the Wulff-Klinger model [22][23][24].As a cell death model, the Arrhenius thermal damage model was used in eight articles, the three-state cell death equation in seven, and an isothermal contour in three articles.Subramanian et al. used their own thermal damage formula [45].Thirteen articles included perfusion in their model, twelve articles included blood vessels, six included water vaporization, and eight included temperature-dependent tissue parameters.Nine articles created a CT-based anatomical model for their simulation.Eight of them segmented the liver, tumor, and blood vessels, while Ooi et al.only derived a liver contour from the CT scan [42].Next to an anatomical model, Moche et al. used dynamic CT measurements to derive perfusion values [41].

RFA Ex Vivo Validation
Table 5 contains the ex vivo validation of the RFA ablation zone models [27,28,32,34,42,45,49,51,54].All experiments were performed on animal livers, either bovine or porcine, or phantoms.The experimental ablation zones were measured after tissue sample sections along the probe axis.The ablation settings differed in all experiments.A visual overview of the longitudinal and transverse RDD is given in Figure 4. Figure 5 contains a combined overview of the MWA and RFA ex vivo experiments.

RFA In Vivo Validation
An overview of the ten articles on RFA ablation zone simulation using in vivo validation is given in Table 6.These consist of four in vivo animal experiments [24,39,43,47], five retrospective clinical studies [22,23,25,40,50], and one prospective clinical study [41].The prospective study of Moche et al. found a DSC of 0.62 ± 0.14 with a surface deviation of 3.4 ± 1.7 mm [41].They concluded that the real-time simulation of RFA-induced tissue necrosis in the liver was fast (3.5 ± 1.9 min) and accurate enough for therapy planning.Mariappan et al. used the same computational model in their retrospective study and found similar results in 23 ablations, with a slightly lower surface deviation (2.50 mm versus 3.4 mm) [40].The simulation accuracy increased by using patient-specific CT-based perfusion values.The results obtained by Audigier et al. indicated lower accuracy of their simulation model compared to the previously mentioned studies [22][23][24].They found a DSC of 0.44 and a surface deviation of 5.3 ± 3.6 mm [24].This difference in results might be explained by the reconstruction of the ablation probe location.Audigier et al. did not reconstruct the clinically used probe location but assumed the center of the tumor as the

RFA In Vivo Validation
An overview of the ten articles on RFA ablation zone simulation using in vivo validation is given in Table 6.These consist of four in vivo animal experiments [24,39,43,47], five retrospective clinical studies [22,23,25,40,50], and one prospective clinical study [41].The prospective study of Moche et al. found a DSC of 0.62 ± 0.14 with a surface deviation of 3.4 ± 1.7 mm [41].They concluded that the real-time simulation of RFA-induced tissue necrosis in the liver was fast (3.5 ± 1.9 min) and accurate enough for therapy planning.Mariappan et al. used the same computational model in their retrospective study and found similar results in 23 ablations, with a slightly lower surface deviation (2.50 mm versus 3.4 mm) [40].The simulation accuracy increased by using patient-specific CT-based perfusion values.The results obtained by Audigier et al. indicated lower accuracy of their simulation model compared to the previously mentioned studies [22][23][24].They found a DSC of 0.44 and a surface deviation of 5.3 ± 3.6 mm [24].This difference in results might be explained by the reconstruction of the ablation probe location.Audigier et al. did not reconstruct the clinically used probe location but assumed the center of the tumor as the

RFA In Vivo Validation
An overview of the ten articles on RFA ablation zone simulation using in vivo validation is given in Table 6.These consist of four in vivo animal experiments [24,39,43,47], five retrospective clinical studies [22,23,25,40,50], and one prospective clinical study [41].The prospective study of Moche et al. found a DSC of 0.62 ± 0.14 with a surface deviation of 3.4 ± 1.7 mm [41].They concluded that the real-time simulation of RFA-induced tissue necrosis in the liver was fast (3.5 ± 1.9 min) and accurate enough for therapy planning.Mariappan et al. used the same computational model in their retrospective study and found similar results in 23 ablations, with a slightly lower surface deviation (2.50 mm versus 3.4 mm) [40].Two articles compared their computational model to more simplistic models [25,39].Hoffer et al. demonstrated an improved accuracy for computational models compared to the manufacturers' chart, with mean surface deviations of 1.1 mm and 2.5 mm, respectively [39].Audigier et al. also found that patient-specific models resulted in higher DSC and surface deviations compared to simplistic ellipse models [25].

Discussion
A wide variety of computational modeling and validation methods was found in this systematic review, making a full-fledged comparison between different models for ablation simulation complex.The studies included in this review used different bioheat equations, parameters, validation methods, and outcome measures or metrics, making it hard to draw conclusions on which computational model performs best and could possibly be implemented in clinical practice.
Considerable differences were identified in the RFA and MWA models, which can be explained by the different heating mechanisms of the two techniques.Thirteen of the sixteen (81%) MWA models included the effect of water vaporization and temperaturedependent parameters, compared to six of the nineteen RFA models (32%).The temperaturedependency of dielectric parameters (electrical conductivity and permittivity) and thermal conductivity is two-folded.First, tissue parameters change as a result of protein denaturation at temperatures above 60 • C [60-62].In addition, water vaporization has an isolating effect.This leads to a decrease in conductivity and permittivity above 100 • C in RFA as well as in MWA [60-62].However, RFA power is usually temperature-or impedance-controlled to avoid vaporization and gas formation, whereas tissue temperatures in MWA frequently exceed 100 • C. Therefore, temperature-dependent parameters, as well as water vaporization, have a greater potential influence on MWA.This could explain the discrepancy in incorporating these parameters in MWA models and RFA models.Moreover, blood vessels are included more frequently in RFA models compared to MWA models (12 out of 19 (63%) and 2 out of 16 (12.5%),respectively).Since the "heat-sink effect" caused by blood vessels near the ablation region is considered to have a greater impact on RFA, this finding seems logical [63,64].When designing an ablation simulation model, these differences in the heating mechanism should be taken into consideration.
The various models applied different bioheat equations, which is the basis for the simulation and therefore affects the simulation accuracy.The majority of the models are purely based on Pennes' bioheat equation due to its simplicity and feasibility.Nevertheless, this equation has an important limitation.The equation only considers microvascular perfusion, assuming a constant blood temperature of 37 • C without flow directionality.However, the blood temperature of vessels within and surrounding the ablation zone will increase during ablation [42].Other bioheat equations that incorporate changes in liver perfusion overcome this limitation while maintaining or even increasing model performances [47,48].In the liver, which is a well-perfused organ with variations in perfusion due to underlying liver diseases, these more advanced bioheat equations are likely to be more accurate.
In addition to the bioheat equation, there are numerous other model design choices affecting the simulation accuracy.This systematic review does not present all possible characteristics of the models.For example, it does not incorporate tissue contraction or two-compartment models.The latter models the tumor and liver as separate compartments, each with its own tissue-specific parameters.This limitation emphasizes the complexity of ablation zone prediction and the number of parameters potentially affecting the ablation zone.Theoretically, the highest accuracy could be achieved by the inclusion of all model characteristics and parameters, but there is a trade-off between accuracy on one hand and complexity and required computational resources on the other.With the recent advancement in artificial intelligence and machine learning, the development of more complex and accurate models becomes feasible.
The chosen input parameter values required for the computational model also contribute to the accuracy achieved.Some input parameters, e.g., thermal conductivity, electrical conductivity, and density, are dependent on tissue composition.They differ between tumor and liver parenchyma as well as per patient as a result of underlying liver disease, i.e., cirrhosis, hepatic steatosis, and previous treatments like systemic treatment or transarterial radioembolization or chemoembolization.
Several studies concluded that thermal conductivity significantly decreases with an increasing fat content [29,[65][66][67].According to the computational models of Deshazer et al. and Servin et al., a greater fat content in the liver leads to larger ablation zone volumes in MWA [66,67].Due to the lower conductivity of surrounding liver parenchyma, heat loss to peripheral parenchyma is limited since low thermally conducting tissue retains high temperatures.This leads to larger ablation zones in MWA.RFA relies more on indirect heating, i.e., heat conductivity, compared to direct heating in MWA.Therefore, smaller ablation zones are observed in livers with increasing fat content.Diminished conductivity of the parenchyma results in higher temperatures at the site of the tumor edge, while the liver parenchyma surrounding the tumor has a lower temperate increase [65].Hypothetically, this "oven-effect" would lead to higher ablation temperatures within the tumor while increasing the risk of narrow ablative margins in RFA.In vivo experiments are needed to test these hypotheses.
In contrast to the widely studied effect of hepatic steatosis in thermal liver ablation, limited research has been conducted on the relation between ablation volumes and cirrhosis or fibrosis.This might be related to the complexity of parenchymal changes due to cirrhosis and fibrosis.Deshazer et al. simulated the difference in perfusion between normal liver tissue and cirrhotic liver tissue and found larger ablation zones in cirrhotic livers compared to normal liver tissue [67].Further studies are needed to investigate heat propagation and conductivity in cirrhotic and fibrotic liver tissue to be able to include these as patient-specific input parameters in liver ablation simulation.
The simulation accuracy is determined by model validation.The included studies varied in validation methods using in or ex vivo experiments with different outcome measures or metrics.Although the bottleneck in current practice is that the manufacturer's prediction is mostly based on ex vivo experiments, 23 out of 38 included models are still validated with ex vivo experiments [11].Moreover, in most studies, the axis length and total volume of the simulated and actual ablation are compared without taking into account the shape and relative position.These variations and limitations in the model validation make it difficult to compare the different simulation models and to identify model parameters that have the largest impact on model accuracy.
To decide on the optimal balance between complexity, accuracy, and computational time, a comparative clinical study should be conducted.This could be a retrospective or prospective study in which the clinically obtained ablation zone is compared to different simulations.This validation necessitates a probe position scan to simulate the ablation at the corresponding probe position, as well as a post-ablation CT scan to segment the clinically obtained ablation zone.To compare the shape and position of the simulated and clinically obtained ablation zone, image registration of the post-ablation and probe position CT scan is required.Variations in liver shape and position should be minimized by identical patient positioning and breathing phase in both scans, which can be controlled with techniques such as high-frequency jet ventilation [68].Nevertheless, image registration remains a challenge due to ablation-induced deformations of the liver.Therefore, registration might introduce inaccuracies to the measurements.After image registration, volumetric outcome metrics, such as the DSC and surface deviations, could be used to compare the overlap of the simulated and clinically obtained ablation zones.These outcome metrics are more expressive in clinical practice compared to size-based metrics since the ablation zone should cover the tumor, including an ablative margin of at least 5 mm.Therefore, the simulation should not only be an accurate prediction of size and volume but also of shape and position relative to the ablation probe.In addition, the maximum surface deviation should be computed since it is indicative of the boundary discrepancies as well as the negative maximum surface deviation and the false positive value (Appendix C, Figure A2c,g), which might be a good measure for the technical success of the ablation.The different models could be compared to determine which combination of included parameters and models results in the most accurate prediction.
The clinical application of ablation zone simulation requires reliable models.Three of the included articles originate from the ClinicIMPPACT project, which aims to bring a planning and simulation tool for RFA into the clinical practice [40,41,50,69].In a clinical study, Moche et al. evaluated this application prospectively and concluded that the model performs sufficiently for clinical implementation with a mean surface deviation of 3.4 mm [41].Nevertheless, it is generally recommended to aim for an ablative margin of at least 5 mm, but the exact relationship between ablative margins and local recurrence is still unknown [70].Therefore, it is questionable if a mean surface deviation of 3.4 mm provides sufficient accuracy.Furthermore, to conclude on the accuracy needed for clinical implementation of computational models, the reliability of the current clinical practice, i.e., the manufacturer's prediction, should be known.Hence, the match between the manufacturer's prediction and the clinically obtained ablation zone should be quantified.Of interest is a study by Hoffer et al.They compared their computational model with the manufacturer's chart and performed in vivo ablations in six swine and found their computational model to have higher accuracy [39].However, clinical studies with a large patient population are still lacking.
Percutaneous thermal ablation techniques have evolved over the most recent decades with new innovations.One advancement is the use of treatment planning software combined with probe navigation using either stereotactic thermal ablation or robotic probe positioning to increase thermal ablation efficacy.However, this planning software still relies on the manufacturer's predicted ablation zones, lacking patient-specific parameters.The full potential of this innovation can only be exploited when reliable patient-specific ablation zone predictions are available [10].The use of image-based two-compartment models may be an important step towards a more accurate, patient-specific model.The comparative study of Audigier et al. has indicated that patient-specific modeling results in more accurate predictions compared to a spherical model [25].Moreover, Mariappan et al. concluded that the simulation yields better accuracy when personalized perfusion values are given as input in the simulation model [40].These values can be extracted from CT images, which also enables two-compartment modeling that incorporates the liver vasculature.These results emphasize the potential of computational modeling for patient-specific ablation planning.However, more insight into parameter personalization is required to enable patient-specific ablation planning and implement this into clinical practice.

Conclusions
Computational models simulating ablation zones in MWA and RFA show considerable heterogeneity in model type and validation methods.It is currently unknown which model is most accurate and best suitable for use in clinical practice.However, several studies have demonstrated a good correlation between simulated ablation zones and in vivo ablations.More research on patient-specific parameters is needed to develop more accurate models that can be used for individualized treatment planning.
LTE and LTNE are models that more extensively model the effect of blood vessels and perfusion in ablation.The LTNE assumes two states, the tissue state and blood state, and therefore consists of two equations.The LTE assumes that tissue and blood have the same temperature and therefore combines the two functions.The subscript "b" represents the characteristics of blood, while subscript "t" represents the characteristics of tissue.
The external generated heat (Q_Extern) Equations used to determine the electromagnetic energy deposition.
Helmholtz harmonic wave equation.

Blood flow models Darcy's law
Used to determine the flow velocity in capillaries and blood vessels for more extensive perfusion modeling and vascular-mediated cooling.

Navier-Stokes equation
Used to determine the flow velocity in capillaries and blood vessels for more extensive perfusion modeling and vascular-mediated cooling.

Newton formula
Describes the heat transfer between blood vessels and tissue.Under constant blood flow conditions, the Newton formula is described as ℎ =     Figure A1c represents the advancement, the distance between the probe tip and the anterior boundary of the ablation zone.The black line represents the advancement of the created ablation zone, which is larger compared to the advancement of the simulated ablation zone (blue line).

Shape-based outcome metrics
Cancers 2023, 15, x FOR PEER REVIEW 26 of 30 The RVD represents the difference in volume of the simulated and experimentally obtained ablation zone.A positive RVD means the simulated ablation zone is larger compared to the created ablation zone in the experiments and vice versa.
Figure A1c represents the advancement, the distance between the probe tip and the anterior boundary of the ablation zone.The black line represents the advancement of the created ablation zone, which is larger compared to the advancement of the simulated ablation zone (blue line).
Shape-based outcome metrics Figure A2a-c represents shape-based boundary outcome metrics.Figure A2a presents a schematic illustration of the surface deviation.The surface deviation is the average Hausdorff distance between the simulated and experimentally created ablation zone surfaces.It is a frequently used outcome measure for boundary error.The maximum surface deviation is given in Figure A2b, while the arrow in Figure A2c is the maximum negative surface deviation.
Figure A2d-g represent shape-based overlap outcome metrics.Figure A2d illustrates the computation of the DSC: two times the volume of overlap (shaded blue area) of the simulated and created ablation zone divided by the sum of both volumes.The sensitivity, Figure A2e, is the volume of overlap divided by the experimentally created ablation zone, and the positive predictive value (PPV), Figure A2f, is the volume of overlap divided by the simulated ablation zone.Figure A2g illustrates the false positive value.It represents the undertreated ablation volume (shaded blue areas), i.e., volume which was simulated as ablated, but in the experiments was not, relative to the total simulated volume.These outcome metrics are all volumetric measurements.Figure A2d-g represent shape-based overlap outcome metrics.Figure A2d illustrates the computation of the DSC: two times the volume of overlap (shaded blue area) of the simulated and created ablation zone divided by the sum of both volumes.The sensitivity, Figure A2e, is the volume of overlap divided by the experimentally created ablation zone, and the positive predictive value (PPV), Figure A2f, is the volume of overlap divided by the simulated ablation zone.Figure A2g illustrates the false positive value.It represents the undertreated ablation volume (shaded blue areas), i.e., volume which was simulated as ablated, but in the experiments was not, relative to the total simulated volume.These outcome metrics are all volumetric measurements.

30 Figure 1 .
Figure 1.Preferred Reporting Items for Systematic Reviews and Meta-Analysis (PRISMA) flow diagram describing the study-selection process.

Figure 1 .
Figure 1.Preferred Reporting Items for Systematic Reviews and Meta-Analysis (PRISMA) flow diagram describing the study-selection process.

Figure 1 .
Figure 1.Preferred Reporting Items for Systematic Reviews and Meta-Analysis (PRISMA) flow diagram describing the study-selection process.

Figure 2 .
Figure 2. Schematic overview of the basic structure of a computational model for simulating thermal liver ablation.The optional dependencies are shown by the dotted lines.Figure 2. Schematic overview of the basic structure of a computational model for simulating thermal liver ablation.The optional dependencies are shown by the dotted lines.

Figure 2 .
Figure 2. Schematic overview of the basic structure of a computational model for simulating thermal liver ablation.The optional dependencies are shown by the dotted lines.Figure 2. Schematic overview of the basic structure of a computational model for simulating thermal liver ablation.The optional dependencies are shown by the dotted lines.
. The main differences between the two models were the use of the three-state cell death model and the incorporation of tissue shrinkage within the model of Sing et al.The latter could explain why simulations by Sing et al. resulted in a smaller longitudinal diameter compared to Wu et al.: 26.24 mm (RDD: −13.4%) and 29.7 mm (RDD = −2.0%),respectively.However, the experiments of Sing et al. resulted in a greater overestimation of the transverse diameter (RDD: 5.2% versus 4.7%).Figure 3 visually gives an overview of the models validated with the longitudinal and transverse RDD.Cancers 2023, 15, x FOR PEER REVIEW 9 of 30 . The main differences between the two models were the use of the three-state cell death model and the incorporation of tissue shrinkage within the model of Sing et al.The latter could explain why simulations by Sing et al. resulted in a smaller longitudinal diameter compared to Wu et al.: 26.24 mm (RDD: −13.4%) and 29.7 mm (RDD = −2.0%),respectively.However, the experiments of Sing et al. resulted in a greater overestimation of the transverse diameter (RDD: 5.2% versus 4.7%).Figure 3 visually gives an overview of the models validated with the longitudinal and transverse RDD.

Figure 3 .
Figure 3. Scatterplot of the relative error in longitudinal and transverse diameters of the modeled MWA zones compared to ex vivo validation[26,30,35,37,38,44,46,52,53,55].In case of an experimental diameter of 30 mm, a relative error of 0.1 means the simulated diameter was 33 mm.

Figure 3 .
Figure 3. Scatterplot of the relative error in longitudinal and transverse diameters of the modeled MWA zones compared to ex vivo validation[26,30,35,37,38,44,46,52,53,55].In case of an experimental diameter of 30 mm, a relative error of 0.1 means the simulated diameter was 33 mm.

*
Relative differences are results of the computational model compared to the experiments.Advancement = the distance from the antenna tip to the boundary of the ablated zone, CT = Computed Tomography, DSC = dice similarity coefficient, MRT = Magnetic Resonance Thermometry, RDD = relative diameter deviation, RVD = relative volume deviation, SAR = specific absorption rate.

*
Relative differences are results of the computational model compared to the experiments.CT = Computed Tomography, DSC = dice similarity coefficient, RDD = relative diameter deviation, RVD = relative volume deviation.

Figure 4 .
Figure 4. Scatterplot of the relative error in longitudinal and transverse diameters of the modeled RFA zones compared to ex vivo validation[27,28,32,51].In case of an experimental diameter of 30 mm, a relative error of 0.1 means the simulated diameter was 33 mm.

Figure 5 .
Figure 5. Scatterplot of the relative error in longitudinal and transverse diameters of the modeled MWA and RFA zones compared to ex vivo validation.In case of an experimental diameter of 30 mm, a relative error of 0.1 means the simulated diameter was 33 mm.

Figure 4 .
Figure 4. Scatterplot of the relative error in longitudinal and transverse diameters of the modeled RFA zones compared to ex vivo validation[27,28,32,51].In case of an experimental diameter of 30 mm, a relative error of 0.1 means the simulated diameter was 33 mm.

Figure 4 .
Figure 4. Scatterplot of the relative error in longitudinal and transverse diameters of the modeled RFA zones compared to ex vivo validation [27,28,32,51].In case of an experimental diameter of 30 mm, a relative error of 0.1 means the simulated diameter was 33 mm.

Figure 5 .
Figure 5. Scatterplot of the relative error in longitudinal and transverse diameters of the modeled MWA and RFA zones compared to ex vivo validation.In case of an experimental diameter of 30 mm, a relative error of 0.1 means the simulated diameter was 33 mm.

Figure 5 .
Figure 5. Scatterplot of the relative error in longitudinal and transverse diameters of the modeled MWA and RFA zones compared to ex vivo validation.In case of an experimental diameter of 30 mm, a relative error of 0.1 means the simulated diameter was 33 mm.
The simulation accuracy increased by using patient-specific CT-based perfusion values.The results obtained by Audigier et al. indicated lower accuracy of their simulation model compared to the previously mentioned studies [22-24].They found a DSC of 0.44 and a surface deviation of 5.3 ± 3.6 mm [24].This difference in results might be explained by the reconstruction of the ablation probe location.Audigier et al. did not reconstruct the clinically used probe location but assumed the center of the tumor as the probe location in their simulation, which introduces an inaccuracy in the measurements.On the other hand, Moche et al. and Mariappan et al. reconstructed the probe location as used in clinical practice by using image registration [40,41].
± 0.84 mm * Relative differences are results of the computational model compared to the experiments.CT = Computed Tomography, DSC = dice similarity coefficient, PPV = positive predictive value, RDD = relative diameter deviation, RVD = relative volume deviation.
Navier-Stokes equationUsed to determine the flow velocity in capillaries and blood vessels for more extensive perfusion modeling and vascular-mediated cooling.ρb δu b δt + ρ b u b •∇u b = −∇P + µ∇ 2 u bNewton formula Describes the heat transfer between blood vessels and tissue.Under constant blood flow conditions, the Newton formula is described ash b = Nu D k b D Table A3.Parameters appearing in blood flow models with their unit and description.Heat transfer coefficient Nu D Local Nusselt number k b W/mK Blood vessel thermal conductivity D M Blood vessel diameter represents these diameters.A positive RDD means the simulated diameters are larger compared to the diameters obtained in the experiments.The RVD represents the difference in volume of the simulated and experimentally obtained ablation zone.A positive RVD means the simulated ablation zone is larger compared to the created ablation zone in the experiments and vice versa.

Figure A2 .
Figure A2.Schematic illustration of the shape-based outcome metrics; (a-c) Representation of surface-based outcome metrics; (a) average surface deviation, (b) maximum surface deviation, (c) maximum negative surface deviation.(d-g) Representation of overlap-based outcome metrics; (d) dice similarity coefficient (DSC), (e) sensitivity, (f) positive predictive value (PPV), and (g) negative predictive value.Note: visualization is in 2D, while all outcome metrics were based on 3D measurements.

Figure A2 .
Figure A2.Schematic illustration of the shape-based outcome metrics; (a-c) Representation of surfacebased outcome metrics; (a) average surface deviation, (b) maximum surface deviation, (c) maximum negative surface deviation.(d-g) Representation of overlap-based outcome metrics; (d) dice similarity coefficient (DSC), (e) sensitivity, (f) positive predictive value (PPV), and (g) negative predictive value.Note: visualization is in 2D, while all outcome metrics were based on 3D measurements.

Figure
Figure A2a-c represents shape-based boundary outcome metrics.FigureA2apresents a schematic illustration of the surface deviation.The surface deviation is the average Hausdorff distance between the simulated and experimentally created ablation zone surfaces.It is a frequently used outcome measure for boundary error.The maximum surface deviation is given in FigureA2b, while the arrow in FigureA2cis the maximum negative surface deviation.FigureA2d-g represent shape-based overlap outcome metrics.FigureA2dillustrates the computation of the DSC: two times the volume of overlap (shaded blue area) of the simulated and created ablation zone divided by the sum of both volumes.The sensitivity, FigureA2e, is the volume of overlap divided by the experimentally created ablation zone, and the positive predictive value (PPV), FigureA2f, is the volume of overlap divided by the simulated ablation zone.FigureA2gillustrates the false positive value.It represents the undertreated ablation volume (shaded blue areas), i.e., volume which was simulated as ablated, but in the experiments was not, relative to the total simulated volume.These outcome metrics are all volumetric measurements.

Table 1 .
Characteristics of the used computational model for microwave ablation zone modeling from the included articles (x = incorporated in the model).

Table 2 .
Characteristics of the used computational model for radio-frequency ablation zone modeling from the included articles (x = incorporated in the model).
et al. used CT data to extract tumor geometry to model tumor coverage, while Zhai et al. created a complete CT-based 3D model for simulating the ablation [36,56].
et al. used CT data to extract tumor geometry to model tumor coverage, while Zhai et al. created a complete CT-based 3D model for simulating the ablation [36,56].

Table 3 .
Ex vivo validation of computational models modeling microwave ablation.

Table 4 .
In vivo validation of computational models modeling microwave ablation.

Table 5 .
Ex vivo validation of computational models modeling radiofrequency ablation.

Table 6 .
In vivo validation of computational models modeling radiofrequency ablation.

Table A1 .
Parameters appearing in the tissue bioheat models with their unit and description.In case parameters might be temperature-dependent, this is denoted in the last column.

Table A2 .
Parameters appearing in external heat equations with their unit and description.In case parameters might be temperature-dependent, this is denoted in the last column.

Table A2 .
Parameters appearing in external heat equations with their unit and description.In case parameters might be temperature-dependent, this is denoted in the last column.

Table A3 .
Parameters appearing in blood flow models with their unit and description.