A New Thermal Damage-Controlled Protocol for Thermal Ablation Modeled with Modified Porous Media-Based Bioheat Equation with Variable Porosity

Thermal ablation of tumors is a minimally invasive technique more and more employed in cancer treatments. The main shortcomings of this technique are, on the one hand, the risk of an incomplete ablation, and on the other hand, the destruction of the surrounding healthy tissue. In this work, thermal ablation of a spherical hepatocellular carcinoma tumor (HCC) surrounded by healthy tissue is modeled. A modified porous media-based bioheat model is employed, including porosity variability from tumor core to healthy tissue, following experimental in vivo measures. Moreover, three different protocols are investigated: a constant heating protocol, a pulsating protocol, and a new developed damage-controlled protocol. The proposed damage-controlled protocol changes the heating source from constant to pulsating according to the thermal damage probability on the tumor rim. The equations are numerically solved by means of the commercial software COMSOL Multiphysics, and the outcomes show that the new proposed protocol is able to achieve the complete ablation in less time than the completely pulsating protocol, and to reach tissue temperature on the tumor rim 10 °C smaller than the constant protocol. These results are relevant to develop and improve different patient-based and automated protocols which can be embedded in medical devices’ software or in mobile applications, supporting medical staff with innovative technical solutions.


Introduction
Nowadays, thermal ablation of tumors is considered an effective technique for fighting cancer. The interest in this technique is growing because of the advantages it involves compared to surgery, chemotherapy, and radiotherapy: it is a so called minimally invasive treatment which yields minor sides effects, shorter hospital stays, and costs [1][2][3]. This technique consists of applying a focused heat source in the tumor region by means of various energy types such as microwaves (MW), radiofrequencies (RF), ultrasounds (HIFU), and laser. During the application, tissue temperature increases until the complete necrosis of the tumoral tissue occurs. The main risks of using this technique are, on the one hand, to achieve an incomplete tumor ablation, thus a consequent tumor recurrence, on the other hand, the destruction of the surrounding healthy tissue if too high temperatures are reached [4,5]. Moreover, the treatment is not patient-based personalized, so it does not consider differences in gender, organ type, tumor dimensions, or position. Finally, software included in medical devices such as generators are not able to predict accurately the power and duration needed to achieve the final ablation zones, probably because they are based on ex vivo measures and too simple bioheat models.
Thus, an accurate engineering modeling of heat transfer and the development of patient-based medical protocols in thermal ablation treatment is essential to predict thermal fields obtained and the consequent ablation zones. Throughout the years, different bioheat models have been implemented [6,7] since Pennes' equation [8]. This first bioheat equation is still widely used for its simplicity, but it is less accurate than other models due to its main shortcomings. Indeed, it does not consider the blood flow direction and it assumes the thermal equilibrium between tissue and venous blood. Furthermore, arterial blood temperature is assumed to be 37 • C and invariable during the treatment. Other bioheat models have been developed, especially the very effective models based on the porous media theory, implemented by Xuan and Roetzel [9] on the basis of the porous media theory [10]. The comparison between these kinds of bioheat models showed that the best tradeoff between accuracy and simplicity is achieved by the porous media-based, especially when the local thermal non-equilibrium (LTNE) equations are employed [11][12][13]. The LTNE model, indeed, considers the biological medium made up of two different phases existing together, namely the tissue phase (the extravascular fraction made up of cells and interstitial spaces) and the blood phase (the blood flowing in the vessels throughout the tissue). The two phases are considered at different temperatures. One of the most relevant parameters of this model is the porosity, namely the fraction of blood in the total volume of tissue. The porosity, indeed, varies among different organs and can significantly affect the results, especially in high vascularized tissues such as liver and kidney [14][15][16]. Thus, LTNE equations have been recently modified to include the porosity variation in the tissue [13,17]. Furthermore, as concerns the different protocols employed in thermal ablation, various pulsating energy sources have been implemented [18,19] and proposed in radiofrequency ablation (RFA) [20][21][22][23], or in microwave ablation (MWA) [24,25].
In this work, three distinct protocols are investigated in a thermal ablation treatment of a hepatocellular carcinoma tumor (HCC): a constant heating source protocol, a pulsating protocol, and a new developed damage-controlled protocol. The developed damagecontrolled protocol uses the thermal damage probability monitoring on the tumor rim to control the heating source, which changes from constant to pulsating when the first point on tumor rim is considered completely ablated. The employed numerical model is a modified porous media-based local thermal non-equilibrium bioheat model which considers the porosity variability between the core and rim of tumor and the surrounding healthy tissue. The variation is modeled using a quadratic function based on experimental in vivo measures on the hepatic blood volume in HCC implanted in rabbits [16]. The aim is to develop new accurate and patient-based bioheat models and protocols in thermal ablation treatments, improving the current procedures and medical devices. In the future, indeed, automated personalized protocols can be obtained, and the complete tumor ablation can be achieved according to specific patients' and tumors' characteristics.

Mathematical Model
The biological medium here investigated consists of a spherical hepatocellular carcinoma tumor (HCC) having 1.20 cm radius (r 1 ) bounded by an external 5.72 cm radius sphere of healthy liver tissue, as displayed in Figure 1. The tumor's dimension is chosen as in [13] according to the mean tumors' size of the 32 HCC tumors treated with microwave thermal ablation according to the study by Amabile et al. [26] led on 30 patients. The symmetry of the computational domain shown in Figure 1b justifies a 2D axisymmetric analysis, which yields the best tradeoff between simplicity and accuracy of the model compared to real cases. So, on the one hand, the computational effort is minor than a 3D analysis [27], and on the other hand, particular scenarios can be implemented such as differences in tumor shapes, the presence of large veins near the tumor, and so on. The developed bioheat model is a porous media-based [9,10], so in the entire biological domain, two different phases coexist: the tissue phase (the extravascular fraction made up of cells and interstitial spaces) and the blood phase (the blood flowing in the vessels throughout Processes 2022, 10, 236 3 of 13 the tissue). The model is modified as in our previous works [13,28] in order to take into account the water content vaporization in both the tissue phase and blood phase separately. In addition, the hepatic blood volume fraction, namely the porosity, is considered to vary between the core and rim of tumor and the surrounding healthy tissue, as experimentally measured in [16].
Processes 2022, 10, x FOR PEER REVIEW 3 of 14 such as differences in tumor shapes, the presence of large veins near the tumor, and so on. The developed bioheat model is a porous media-based [9,10], so in the entire biological domain, two different phases coexist: the tissue phase (the extravascular fraction made up of cells and interstitial spaces) and the blood phase (the blood flowing in the vessels throughout the tissue). The model is modified as in our previous works [13,28] in order to take into account the water content vaporization in both the tissue phase and blood phase separately. In addition, the hepatic blood volume fraction, namely the porosity, is considered to vary between the core and rim of tumor and the surrounding healthy tissue, as experimentally measured in [16]. Moreover, the LTNE assumptions are employed, so the two phases are considered at different temperature values, and two distinct energy equations are used as follows.
For the tissue phase: for the blood phase: where εvar is the variable porosity; ρ the density; c the specific heat; Tt and Tb are tissue and blood temperatures, respectively; u the blood velocity vector; k the thermal conductivity; t is the time; hc is the interfacial heat transfer coefficient between tissue phase and blood phase, assumed to be constant and equal to 170 W m −2 K −1 as in [29]; and a is the volumetric transfer area between tissue and blood, defined as: where d is the diameter of the considered blood vessels. Furthermore, the volume average technique is used [10] in order to consider volume averaged properties, so the symbol <> refers to the average volume of a generic variable and it is neglected from this point onwards.
As concerns the blood flow, the velocity of blood into the vessel is supposed uniform in the four main domain directions (i.e., r, −r, z, −z) to obtain a sort of realistic in vivo isotropic small vessels network. Terminal arteries are assumed for blood vessels, so their diameter and the consequent blood velocity are d = 3 × 10 −5 m and 4 × 10 −3 m s −1 , respectively [30]. This is a reasonable choice considering that blood exchange is assumed Moreover, the LTNE assumptions are employed, so the two phases are considered at different temperature values, and two distinct energy equations are used as follows.
For the tissue phase: for the blood phase: where ε var is the variable porosity; ρ the density; c the specific heat; T t and T b are tissue and blood temperatures, respectively; u the blood velocity vector; k the thermal conductivity; t is the time; h c is the interfacial heat transfer coefficient between tissue phase and blood phase, assumed to be constant and equal to 170 W m −2 K −1 as in [29]; and a is the volumetric transfer area between tissue and blood, defined as: where d is the diameter of the considered blood vessels. Furthermore, the volume average technique is used [10] in order to consider volume averaged properties, so the symbol <> refers to the average volume of a generic variable and it is neglected from this point onwards.
As concerns the blood flow, the velocity of blood into the vessel is supposed uniform in the four main domain directions (i.e., r, −r, z, −z) to obtain a sort of realistic in vivo isotropic small vessels network. Terminal arteries are assumed for blood vessels, so their diameter and the consequent blood velocity are d = 3 × 10 −5 m and 4 × 10 −3 m s −1 , respectively [30]. This is a reasonable choice considering that blood exchange is assumed to occur mainly downstream of terminal arteries and before the arterioles, according to Chen and Holmes [31]. In order to take into account the effect of thermal damage on blood velocity, namely the blood coagulation and the consequent blood flow interruption, the coefficient β is inserted in Equation (2) related to the blood phase. This coefficient, indeed, can be 0 or 1 depending on if the tissue is completely damaged or not, respectively. In order to evaluate the thermal damage level of the tissue, the Arrhenius' model is adopted, considering an exponential relationship between tissue exposure temperature, time, and the kinetic parameters that consider the effects of thermal degradation of proteins [32] as follows: where A is the frequency factor, ∆E the activation energy for the irreversible damage reaction, and R is the universal gas constant. In this work, A = 7.39 × 10 39 s −1 and ∆E = 2.577 × 10 5 J mol −1 for healthy liver tissue, while for the HCC tumor A = 3.247 × 10 43 s −1 and ∆E = 2.814 × 10 5 J mol −1 [33]. Moreover, thermal damage function Ω can be related to the cell death probability as: In this work, in order to achieve accurate values of thermal damage, the 99% cell death probably is considered, corresponding to Ω = 4.6. Thus, the coefficient β becomes zero as Ω = 4.6, so no blood flows in the coagulated tissue and blood velocity becomes zero.
In addition, water content vaporization in the two phases is included separately, employing the enthalpy method as in [34], thus the term (ρc) in Equations (1) and (2) becomes: where ρ l and c l are density and specific heat evaluated for the liquid phase, ρ g and c g are density and specific heat for the gas phase, h fg is the product of water latent heat of vaporization and water density at 100 • C, and C w represents the water percentage in healthy tissue (68%), HCC tumor tissue (81%) or blood (79%), respectively, as in literature [35,36]. As concerns ∆T b,t , temperature difference is supposed 1 • C [34]. The external power density value Q ext is chosen to model a MWA treatment, and it represents the mean total power dissipation density achieved during a 600 s microwave thermal ablation treatment, as it will described in the following section about validation of the model. It results from the electrical field vector E distribution obtained in our previous study conducted for a microwave ablation case for the same volumes of tumor and healthy tissue [13]: where σ e is the temperature dependent effective electric conductivity at 2.45 GHz evaluated as in [13] and |E| is the Euclidian norm of E. The resulting mean external power density values for tumor domain and healthy tissue domain are 2.49 × 10 6 (W m −3 ) and 3.19 × 10 4 (W m −3 ), respectively. The porosity is the most relevant parameter of the model because the results are significantly influenced by its value [10,12,13,17], especially in high vascularized organs, such as liver. The higher the porosity, indeed, the higher the convective term related to mass blood flow and the consequent heat sink effect, which causes the tissue temperature decrease [12]. In order to simulate accurately the porosity variation inside the domain, a porosity function is implemented according to experimental in vivo measures on implanted HCC tumors in rabbits [16]. Inside the tumor, indeed, the blood volume fraction varies from 0.07 at the tumor core level to 0.23 at the rim, while the maximum value of 0.31 is observed in the healthy tissue. Thus, a quadratic function is supposed to represent well the porosity variation inside the tumor, as follows: where the coefficients a 1 , b 1 , and c 1 are obtained imposing the following conditions: thus, substituting the porosity values ε min , ε int , and ε max the coefficients a 1 , b 1 , and c 1 become: where r 1 = 1.20 cm is the radius of the tumor. In Figure 2, the quadratic porosity function is shown in both r and z directions of the domain. As concerns thermal properties for tumor tissue, healthy tissue, and blood, Table 1 resumes the values chosen in this work [5,33]. As regards boundary and initial conditions, both tissue and blood temperatures are considered constant on the external sphere boundary and equal to 37 • C. This choice is justified by the domain's dimension, large enough to ignore boundary effects from surroundings. Moreover, the axial symmetry of the problem allows to fix the adiabatic condition on the r = 0 symmetry axis and the initial temperature is 37 • C for both tissue and blood.
thus, substituting the porosity values εmin, εint, and εmax the coefficien where r1 = 1.20 cm is the radius of the tumor. In Figure 2, the quad is shown in both r and z directions of the domain. As concerns tumor tissue, healthy tissue, and blood, Table 1 resumes the valu [5,33]. As regards boundary and initial conditions, both tissue and considered constant on the external sphere boundary and equal justified by the domain's dimension, large enough to ignore b surroundings. Moreover, the axial symmetry of the problem allo condition on the r = 0 symmetry axis and the initial temperature and blood.   Table 1. Thermal properties employed for the different materials [5,33].

Numerical Solution and Validation
The problem is numerically solved with the commercial software COMSOL Multiphysics (COMSOL, Burlington, MA, USA) based on the finite element method. A triangular mesh of 26,842 elements is chosen, after a grid convergence test made on temperature profiles at tumor rim (r = 1.2 cm; z = 0.0 cm), obtained using distinct grids with increasing element numbers, as displayed in Figure 3. The solver employed is the direct PARDISO and second-order Lagrangian elements are used for their discretization. The absolute tolerance is 0.0001 for transient solver, and the intermediate backward differentiation formulas (BDF) method is employed for the time step, with initial and maximum steps of 0.001 s and 1 s, respectively. As made for the grid, the time step choice is verified on temperature evolution at tumor rim as made for the grid test, as shown in Figure 4.

Numerical Solution and Validation
The problem is numerically solved with the commercial softwa Multiphysics (COMSOL, Burlington, MA, USA) based on the finite eleme triangular mesh of 26,842 elements is chosen, after a grid convergence temperature profiles at tumor rim (r = 1.2 cm; z = 0.0 cm), obtained using with increasing element numbers, as displayed in Figure 3. The solver em direct PARDISO and second-order Lagrangian elements are used for their d The absolute tolerance is 0.0001 for transient solver, and the intermedi differentiation formulas (BDF) method is employed for the time step, wi maximum steps of 0.001 s and 1 s, respectively. As made for the grid, the tim is verified on temperature evolution at tumor rim as made for the grid test Figure 4.

Numerical Solution and Validation
The problem is numerically solved with the commercial softwar Multiphysics (COMSOL, Burlington, MA, USA) based on the finite elemen triangular mesh of 26,842 elements is chosen, after a grid convergence te temperature profiles at tumor rim (r = 1.2 cm; z = 0.0 cm), obtained using d with increasing element numbers, as displayed in Figure 3. The solver emp direct PARDISO and second-order Lagrangian elements are used for their d The absolute tolerance is 0.0001 for transient solver, and the intermedia differentiation formulas (BDF) method is employed for the time step, wit maximum steps of 0.001 s and 1 s, respectively. As made for the grid, the tim is verified on temperature evolution at tumor rim as made for the grid test, Figure 4.

Validation of the Model
The validation of the developed bioheat model is made on the mean transverse ablation diameters obtained in vivo by Amabile et al. [26] after 300 s and 600 s of microwave thermal ablation treatment on 32 HCC tumors of 2.4 cm mean diameter using the same 14 G HS Amica-Gen Probe modeled in our previous numerical study [13]. Moreover, the same protocol (constant 60 W power applied for 600 s) is used to simulate the electromagnetic field from which the mean total power dissipation density is evaluated and included in the present model as the external power density source, as described previously by Equation (7). So, the comparison is between the results obtained by the numerical models and the experimental in vivo outcomes considering the same tumor dimensions, microwave antenna, power level, and duration of treatment. In Table 2, the transverse ablation zone's diameters achieved are displayed and compared with experimental measures performed by Amabile et al. [26], and in Figure 5, the thermal damage probability is shown and compared to the final ablation radius experimentally obtained. Both Table 2 and Figure 5 confirm that the achieved outcomes with the proposed modified LTNE bioheat model with variable porosity agree with the experimental measures. The differences in terms of ablation diameters, indeed, is only −9.1% after 300 s and +2.7% after 600 s.

Validation of the Model
The validation of the developed bioheat model is made on the mean transverse ablation diameters obtained in vivo by Amabile et al. [26] after 300 s and 600 s of microwave thermal ablation treatment on 32 HCC tumors of 2.4 cm mean diameter using the same 14 G HS Amica-Gen Probe modeled in our previous numerical study [13]. Moreover, the same protocol (constant 60 W power applied for 600 s) is used to simulate the electromagnetic field from which the mean total power dissipation density is evaluated and included in the present model as the external power density source, as described previously by Equation (7). So, the comparison is between the results obtained by the numerical models and the experimental in vivo outcomes considering the same tumor dimensions, microwave antenna, power level, and duration of treatment. In Table  2, the transverse ablation zone's diameters achieved are displayed and compared with experimental measures performed by Amabile et al. [26], and in Figure 5, the thermal damage probability is shown and compared to the final ablation radius experimentally obtained. Both Table 2 and Figure 5 confirm that the achieved outcomes with the proposed modified LTNE bioheat model with variable porosity agree with the experimental measures. The differences in terms of ablation diameters, indeed, is only −9.1% after 300 s and +2.7% after 600 s.  [26] 3.3 ± 0.5 3.7 ± 0.3

Description of the Different Applied Protocols
In this work, three different thermal ablation protocols are investigated in order to achieve the complete ablation of the tumor plus a 0.5 cm margin to exclude the recurrence of the tumor in the periablational zone, as described in [5]. The first protocol consists of the already described constant generation applied during the entire treatment, the second one uses a pulsating power density source for all the duration of thermal ablation, while the last protocol is a mixed constant-pulsating based on the thermal damage function. As concerns the pulsating external heating source, a cosine function is employed since larger target tissue ablation can be achieved by means of this type of wave compared with the half-square wave of the commercial generators [37]. In addition, the cosine function

Description of the Different Applied Protocols
In this work, three different thermal ablation protocols are investigated in order to achieve the complete ablation of the tumor plus a 0.5 cm margin to exclude the recurrence of the tumor in the periablational zone, as described in [5]. The first protocol consists of the already described constant generation applied during the entire treatment, the second one uses a pulsating power density source for all the duration of thermal ablation, while the last protocol is a mixed constant-pulsating based on the thermal damage function. As concerns the pulsating external heating source, a cosine function is employed since larger target tissue ablation can be achieved by means of this type of wave compared with the half-square wave of the commercial generators [37]. In addition, the cosine function allows to obtain cos(t = 0) = 1, so as to avoid the "0" value at the simulation start. So, the pulsating power density is expressed as follows: where ω is the pulsation and t the time. The selected pulsation value is ω = 0.876 s −1 , which corresponds to a period of 7.17 s and avoids too wide oscillations in terms of temperature according to our previous work [19]. The third developed protocol, instead, consists in keeping the external heating source constant until the maximum value of thermal damage function Ω becomes 4.6 on any point of the tumor rim (corresponding to the 99% cell death probability). As this condition is reached, the pulsating function is applied to the heat source to achieve the complete ablation, as described by the following equation:

Results and Discussion
The outcomes are presented in terms of ablation volume evolution during the treatment in order to obtain the different treatment duration and tissue temperature achieved during the three different protocols in two different critical points of the domain. The specific objective in this work is to achieve the complete ablation plus a 0.5 cm margin to exclude a future recurrence of the tumor in the periablational zone, as described in [5]. Thus, the objective final ablation radius is r = 1.7 cm, which corresponds to the spherical volume V obj = 20.6 cm 3 . First of all, the three different protocols led to three distinct treatment durations as shown in the following Table 3: Table 3. Different treatment durations for the three evaluated protocols.

Protocol t fin [s]
Constant heat source 425 Damage-controlled heat source 600 Pulsating heat source 715 From the table, it is clear that the constant heat source applied during the entire application yields a shorter treatment compared to the pulsating protocol (+290 s), and the mixed constant-pulsating protocol with the damage-controlled heat source (+175 s). This result is confirmed in Figure 6, where the percentage of ablation volume on the total objective volume is displayed for the different duration of treatments, resulting in a linear evolution after an initial logarithmic increase for all the cases. From the figure, the greater slope for the constant protocol than the other two can be noticed, which confirms the shorter duration of this type of protocol.
In Figures 7 and 8, tissue temperature evolution at the specific point of the tumor rim (r = 1.2 cm and z = 0.0 cm), and at the target point to obtain the security plus 0.5 cm ablation margin (r = 1.7 cm and z = 0.0 cm) are displayed for the three different protocols. From Figure 7, the highest temperature is achieved by means of the constant protocol at the end of treatment (86 • C), while both the pulsating and the damage-controlled protocols yield equal final temperatures of 76 • C, which are 10 • C lower than constant protocol's temperature. This temperature difference is relevant to determine the necessary treatment duration to achieve the complete ablation target. Processes 2022, 10, x FOR PEER REVIEW 9 of 14 Figure 6. Percentage of ablated volume on the total objective volume during the distinct protocols.
In Figures 7 and 8, tissue temperature evolution at the specific point of the tumor rim (r = 1.2 cm and z = 0.0 cm), and at the target point to obtain the security plus 0.5 cm ablation margin (r = 1.7 cm and z = 0.0 cm) are displayed for the three different protocols. From Figure 7, the highest temperature is achieved by means of the constant protocol at the end of treatment (86 °C), while both the pulsating and the damage-controlled protocols yield equal final temperatures of 76 °C, which are 10 °C lower than constant protocol's temperature. This temperature difference is relevant to determine the necessary treatment duration to achieve the complete ablation target. In Figure 8, indeed, temperature evolution at target ablation point shows similar values, the highest final temperature of the constant protocol (52 °C) differs only 2 °C from the other two protocols (50 °C), but they are obtained with the different durations already shown in Table 3. Thus, the damage-controlled developed protocol is able to achieve the complete ablation in a smaller time than the completely pulsating protocol, reducing at the same time tissue temperature on the tumor rim compared to the constant protocol. This type of protocol allows to control the energy source on the basis of thermal damage level, avoiding on the one hand, too high temperature peaks, which can lead to an undesired quick destruction of the surrounding healthy tissue, and on the other hand, too  In Figures 7 and 8, tissue temperature evolution at the specific point of the tumor rim (r = 1.2 cm and z = 0.0 cm), and at the target point to obtain the security plus 0.5 cm ablation margin (r = 1.7 cm and z = 0.0 cm) are displayed for the three different protocols. From Figure 7, the highest temperature is achieved by means of the constant protocol at the end of treatment (86 °C), while both the pulsating and the damage-controlled protocols yield equal final temperatures of 76 °C, which are 10 °C lower than constant protocol's temperature. This temperature difference is relevant to determine the necessary treatment duration to achieve the complete ablation target. In Figure 8, indeed, temperature evolution at target ablation point shows similar values, the highest final temperature of the constant protocol (52 °C) differs only 2 °C from the other two protocols (50 °C), but they are obtained with the different durations already shown in Table 3. Thus, the damage-controlled developed protocol is able to achieve the complete ablation in a smaller time than the completely pulsating protocol, reducing at the same time tissue temperature on the tumor rim compared to the constant protocol. This type of protocol allows to control the energy source on the basis of thermal damage level, avoiding on the one hand, too high temperature peaks, which can lead to an undesired quick destruction of the surrounding healthy tissue, and on the other hand, too  In order to have a clearer overview of all the results achieved, Figure 9 shows temperature fields obtained at the end of each protocol. From the figure, damagecontrolled and pulsating protocols results in lower temperatures in tumor domain compared to the constant one. The goal of obtaining a complete ablation of the target In Figure 8, indeed, temperature evolution at target ablation point shows similar values, the highest final temperature of the constant protocol (52 • C) differs only 2 • C from the other two protocols (50 • C), but they are obtained with the different durations already shown in Table 3. Thus, the damage-controlled developed protocol is able to achieve the complete ablation in a smaller time than the completely pulsating protocol, reducing at the same time tissue temperature on the tumor rim compared to the constant protocol. This type of protocol allows to control the energy source on the basis of thermal damage level, avoiding on the one hand, too high temperature peaks, which can lead to an undesired quick destruction of the surrounding healthy tissue, and on the other hand, too long duration, and consequently extended pain for the patient. Moreover, using the developed porous media-based numerical model, different protocols can be adapted to the different tissues' features, tumors' dimensions, and position. Information on blood velocity, the blood vessels' diameters, and arrangements can be achieved in clinical practice before or even during the treatment.
In order to have a clearer overview of all the results achieved, Figure 9 shows temperature fields obtained at the end of each protocol. From the figure, damage-controlled and pulsating protocols results in lower temperatures in tumor domain compared to the constant one. The goal of obtaining a complete ablation of the target volume, instead, is reached quicker than the pulsating protocol when the damage-controlled one is used. In order to have a clearer overview of all the results achieved, Figure 9 sho temperature fields obtained at the end of each protocol. From the figure, dama controlled and pulsating protocols results in lower temperatures in tumor dom compared to the constant one. The goal of obtaining a complete ablation of the tar volume, instead, is reached quicker than the pulsating protocol when the dama controlled one is used.

Conclusions
In this work, thermal ablation of a hepatocellular carcinoma tumor is modeled investigate three different protocols which differ for the external heating sou application: a constant heating protocol, a pulsating protocol, and a new develo damage-controlled protocol. The goal is to obtain a complete thermal ablation of

Conclusions
In this work, thermal ablation of a hepatocellular carcinoma tumor is modeled to investigate three different protocols which differ for the external heating source application: a constant heating protocol, a pulsating protocol, and a new developed damage-controlled protocol. The goal is to obtain a complete thermal ablation of the tumor region of 1.2 cm, plus a 0.5 cm margin to exclude the recurrence of the tumor in the periablational zone. The entire domain consists of an internal sphere with tumor properties bounded by an external sphere representing the healthy tissue. The investigation is made employing a new developed porous media-based local thermal non-equilibrium bioheat model modified in order to include the variability of the blood volume fraction, namely the porosity, inside the domain. The porosity variation between the tumor core, tumor rim, and healthy tissue is based on experimental in vivo measures on the hepatic blood volume in HCC implanted in rabbits. The equations are numerically solved by means of the commercial software COMSOL Multiphysics, and the outcomes show different final treatment durations and temperatures. The constant heating protocol, indeed, results in smaller application time but higher temperature than the other two protocols, so the risk of healthy tissue ablation becomes higher. The developed damage-controlled protocol, instead, uses the thermal damage probability monitoring on the tumor rim to control the heating source changing from constant to pulsating. This protocol is able to achieve the complete ablation in less time than the completely pulsating protocol, and at the same time to reach tissue temperature on the tumor rim 10 • C smaller than the constant protocol, as happens with the totally pulsating protocol. These results are relevant to model more and more accurately thermal ablation treatments, in order to achieve the complete tumor ablation and avoid the risk of its recurrence. Moreover, different damage-controlled protocols can be developed and improved in the future, implementing procedures based on real data from imaging techniques. Patient-based and automated protocols can be developed and employed in software embedded in medical devices or in mobile applications.