Mathematical modeling of vaporization during laser-induced thermotherapy in liver tissue

Laser-induced thermotherapy (LITT) is a minimally invasive method causing tumor destruction due to heat ablation and coagulative effects. Computer simulations can play an important role to assist physicians with the planning and monitoring of the treatment. Our recent study with ex-vivo porcine livers has shown that the vaporization of the water in the tissue must be taken into account when modeling LITT. We extend the model used for simulating LITT to account for vaporization using two different approaches. Results obtained with these new models are then compared with the measurements from the original study.


Introduction
Thermal ablation methods briefly generate cytotoxic temperatures in tumorous tissue in order to destroy it. These minimally invasive methods are used for treating cancer, e.g., in lung, liver, or prostate, when surgical resection is either not possible or too dangerous for the patient. All of these methods utilize the fact that tumorous tissue is more susceptible to heat than healthy tissue to destroy as little healthy tissue as possible. Among the most common thermal ablation methods are LITT, radio-frequency ablation, and microwave ablation.
The principle of LITT [1] is based on the local supply of energy via an optical fiber, located in a water-cooled applicator. This applicator is placed directly into the tumorous tissue. The LITT treatment can take place under MRI control because the laser applicator is sourced by an optical fiber and does not include any metal parts. Therefore the patient is not exposed to radiation, in contrast to other treatments that can only be carried out under CT control.
For the therapy planning of LITT, accurate numerical simulations are needed to guide the practitioner in deciding when to stop the treatment. Mathematical models for this have been proposed, e.g., in [2,3]. The liver consists of about 80% water which vaporizes if the temperatures during the treatment become sufficiently large. The vaporization of this water is currently not included in these models but our study in [4,5] suggests that this effect is relevant for an accurate simulation. In this study the ex-vivo experiments with a larger power of 34 W show a good agreement between measured and simulated temperature until the temperature reaches approximately 100 • C. Then, the measured temperature stagnates while the simulated one rises further (cf. [4], Fig. 3). We presume that this happens due to phase change of water which was not included in the model we used.
In this paper we use the measurements from [4] and compare two models for the vaporization. One of them is the effective specific heat (ESH) model introduced in [6] which modifies the specific heat coefficient to account for the phase change. The other one is the enthalpy model which uses an additional state equation to model the phase transition. We compare the models to experimental data with ex-vivo porcine livers from [4].
Of course the presence of vapor makes the situation far more complicated. The vapor expands, pressure builds up and the vapor has its own dynamics within the tissue. Once the vapor reaches a cooler region it may condensate again. There are many approaches studying this in detail [7,8]. The drawback of course is that such a detailed approach makes the model far more complex, at the costs of computational time, and introduces new tissue dependent parameters, which may not be easily available. Therefore, in this study we use and extremly simplified approach to model the vapor which was proposed by [6] and does not include any physically motivated transport mechanism for the vapor. One purpose of this study is to investigate if this simplified approach may be sufficient for modeling LITT or if more advanced models are necessary to account for the vapor (see also Remark 2). This paper is structured as follows. Our existing mathematical model for simulating LITT including heat and radiative transfer is described in Sect. 2. This model is based on the work of [2] and we have also used it in [4]. In Sect. 3 we modify and extend this model in such a way that it also covers the effect of vaporization during the treatment. Therefore, we consider both the ESH model of [6] as well as an enthalpy model for vaporization. Afterwards, we present the details of the numerical solution of our models in Sect. 4. Finally, the models are validated with measurement data obtained from experiments made with ex-vivo porcine liver tissue (cf. [4]) in Sect. 5.

Mathematical model
We denote by Ω ⊂ R 3 the geometry of the liver and by Γ = ∂Ω its boundary. The latter consists of the radiating surface of the adjacent applicator Γ rad , the cooled surface of the applicator Γ cool , and the ambient surface of the liver Γ amb (see Fig. 1). The mathematical model is described by a system of partial differential equations (PDEs) for the heat transfer inside the liver, the radiative transfer from the applicator into the liver tissue, and a model for tissue damage (cf. [2][3][4]).

Heat transfer
The heat transfer in the liver tissue is modeled by the well-known bio-heat equation (cf. [9]) where T = T(x, t) denotes the temperature of the tissue, depending on the position x ∈ Ω and the time t ∈ (0, τ ). Here, the end time of the simulation is denoted by τ > 0. Further, Figure 1 Sketch of the geometry including the water-cooled applicator with radiating laser fiber C p is the specific heat capacity, ρ the density of the tissue, and κ the thermal conductivity. The perfusion rate due to blood flow is denoted by ξ b and the blood temperature by T b .
Note that in the current ex-vivo study the perfusion rate ξ b is set to zero. Finally, Q rad is the energy source term due to the irradiation of the laser fiber and the initial tissue temperature distribution is given by T init . For the heat transfer between the tissue and its surroundings, given by the ambient surface and the applicator, the following Robin type boundary conditions are used Here, n is the outer unit normal vector on Γ . Additionally, α cool and α amb are the heat transfer coefficients for the water-cooled part of the applicator and the surroundings of the liver, respectively. The temperature of the cooling water is denoted by T cool and T amb is the ambient temperature.
Remark 1 Please note that the temperature T cool of the water coolant is assumed to be known and constant in this study. This is of course a simplification because the cooling water is heated up on its way through the applicator. However, measurements of the cooling temperature before and after the applicator in [4,Fig. 2] show that the temperature of the coolant does not increase by more than 5 • C. Therefore, setting T cool to the measured inlet coolant temperature should approximate the problem. Of course it is also possible to model the flow through the applicator in detail as done in [2].
We come back to this bio-heat equation in Sect. 3, where we modify it such that it also covers the effect of vaporization of water in the tissue. The radiative source term Q rad is defined in the next section by (5).

Radiative transfer
The irradiation of laser light is modeled by the radiative transfer equation where the radiative intensity I = I(s, x) depends on a direction s ∈ S 2 on the (unit) sphere and the position x ∈ Ω, and μ a and μ s are the absorption and scattering coefficients, re-spectively. In particular, as that radiative transfer happens significantly faster than temperature transfer, we neglect the time-dependence and use this quasi-stationary model. The scattering phase function P(s · s ) is given by the Henyey-Greenstein term which reads (cf. [10]) Here, g ∈ [-1, 1] is the so-called anisotropy factor that describes backward scattering for g = -1, isotropic scattering in case g = 0 and forward scattering for g = 1.
Due to the high dimensionality of the radiative transfer equation (2), we use the so-called P 1 -approximation to model the radiative energy, the details of which can be found, e.g., in [11]. Introducing the ansatz where q(x) = 1 4π S 2 I(s, x)s ds is radiative flux vector, one obtains the much simpler threedimensional diffusion equation where ϕ = ϕ(x) is the radiative energy and the diffusion coefficient D is given by To derive the boundary conditions we use Marshak's procedure as described in, e.g., [11]. We obtain Robin type boundary conditions where q app is the laser power entering the tissue and A Γ rad the surface area of the radiating part of the applicator. The former can be written as whereq is the configured laser power and the factor (1β q ) models the absorption of energy by the coolant (cf. [4]). Moreover, the parameter b in (4) is given as b = 0.5 on Γ amb and b = 0 on Γ cool . From the numerical point of view the system given by (3) and (4) is much easier to solve than the original system given by (2). Finally, the radiative energy is used to define the source term for the bio-heat equation in the following way

Tissue damage and its influence on optical parameters
The optical parameters μ a , μ s and g are very sensitive to changes of tissue's state. In particular, once the coagulation of cells starts, their optical parameters change and, as a result, the radiation cannot enter the tissue as deeply as before. Therefore, we model the damage of the tissue as in, e.g., [2,3] with the help of the Arrhenius law, which is given by with so-called frequency factor A, activation energy E a , and universal gas constant R. This describes the change of optical parameters due to coagulation in the following way where the subscripts n and c indicate properties of native and coagulated tissue, respectively (cf. [2]).

Mathematical modeling of vaporization
Vaporization of water inside organic materials plays an important role in many different fields, e.g., in medicine or the food industry. To model the temperature distribution in such materials correctly, it is important to take the vaporization into account as a significant amount of energy is necessary for the phase transition from water to vapor. The basic principle is the following (see, e.g., [12]). If energy in the form of heat is added to water (under constant pressure), the water's temperature increases as long as it is below the vaporization temperature, i.e., below 100 • C. However, as soon as the water reaches this temperature, the temperature does not increase further, although heat is still added to the water. At this point, the water starts to boil and eventually vaporizes after a sufficient amount of energy was added to it. Finally, the temperature of the emerging water vapor increases again after all water has been vaporized. This happens due to the fact that the energy added to the water at its boiling point is used to change its phase and not to increase its temperature, until all water is vaporized. In the following, we discuss two vaporization models. First, we take a look at the effective specific heat (ESH) model introduced in [6] which uses a varying specific heat capacity to model the phase change. In this model the phase transition is spread over a reasonably small interval around 100 • C. This simplification makes it possible to model the phase transition using a single PDE. Second, we propose an enthalpy model with an additional state equation for the enthalpy. For this model, the transition happens at a single temperature, namely at 100 • C.

The effective specific heat (ESH) model
The ESH model introduced in [6] considers the following modified bio-heat equation with the same initial and boundary conditions as (1). Here, Q vap is a source term that models the vaporization of water and Q cond is the source term for the condensation (see Sect. 3.2). In [6] this has the following form where λ denotes the latent heat of water and W is the tissue water density. Using the chain rule we see that Substituting this into (8) and (7) gives the following modified heat equation where the effective specific heat capacity C p is given by Since ∂W ∂T < 0 for vaporization (the water content decreases with temperature), we have that C p ≥ C p .
Based on experiments that measured water content of bovine liver as a function of temperature in [13] the following function is used to describe the tissue water density (cf. [6,13]) where S(T) is the cubic C 1 spline that interpolates between the two exponential functions, (approximately) given by The function W and its derivative are depicted in Fig. 2. In particular, we get that the effective specific heat is very large in an area around 100 • C. Therefore, it holds that which models the vaporization of the tissue water.

Simple condensation model for ESH model
In [6] it is discussed that, in addition to the vaporization of water, one also needs to consider the effect of condensation of the water vapor in order to obtain an accurate model. There, it was assumed that the water vapor diffuses into a region of lesser temperatures where it condensates and releases its latent heat obtained through the vaporization. The authors of [6] describe their model for this in the following way. They say that they first calculate the total amount of water that was vaporized in the last time step. From this, the amount of latent heat generated is computed. Finally, this is added uniformly to the tissue region whose temperature is between 60 • C and 80 • C. We have implemented this simple condensation model in the following way. We compute the total amount of latent heat which is currently consumed through the vaporization of water bȳ where [Q vap ] = W. Additionally, we define the condensation region as for given temperature boundaries T -< T + < 100 • C. Uniformly distributingQ vap over the condensation region then yields the condensation source term In particular, this implies that our model is energy conserving. This is of course a very rough condensation model because there is no real transport mechanism for the vapor involved at all. Any vapor will instantaneously condensate in another region with lower temperature. This simple model shows promising results but there is also room for improvement as discussed in Sect. 5.4.
Remark 2 Clearly this approach for dealing with the vapor is a severe simplification of what actually happens: the expanding vapor builds up pressure and moves through the tissue, thus adding a fluid dynamical component to the problem. The vapor transport in the simple model is purely artificial and not motivated by physics. The trade-off is between a very simple model and a more accurate one which is also far more complex. One goal of this study is to investigate if the simple model may be sufficient for LITT or if more advanced models for the vapor are needed [7,8].

Enthalpy model
In the this section, we present the details of the second model for vaporization, which is based on an enthalpy formulation. It consists of two coupled equations, one for the temperature of the tissue and one for its enthalpy. For the temperature, we have the following, modified bio-heat equation where λ vap = 0.8λ is the proportion of the enthalpy of vaporization corresponding to the tissue's water content of 80%. Further, the (volumetric) enthalpy of the water H, [H] = J m -3 , is modeled by the following ODE Equation (9) has the same initial and boundary conditions as (1), and the initial condition of the enthalpy is given by H = 0 in Ω, i.e., no vaporization had happened before the treatment. The term Q cond describes a heat source due to the condensation of water vapor in regions with temperatures below 100 • C, similar to the one of the ESH model (cf. Sect. 3.2). Observe that the modified bio-heat equation (9) coincides with the classical bioheat equation (1) and we also have H = 0, i.e., no vaporization is happening, as long as we have that T < 100 • C everywhere. This changes as soon as T = 100 • C at some pointx ∈ Ω. Then, we see that the bio-heat equation (9) gives ∂T ∂t (x) = 0 and, therefore, T(x) = 100 • C in case 0 ≤ H(x) < ρλ vap , i.e., the temperature at a point does not change until the enthalpy exceeds the enthalpy of vaporization ρλ vap . In the meantime, the energy that would usually lead to an increase in temperature now only increases the enthalpy, which models the phase change of the water in the tissue. Finally, as soon as the enthalpy reaches the enthalpy of vaporization, all water is vaporized and the bio-heat equation is valid again.

Simple condensation model for enthalpy model
Similar to Sect. 3.2 the simple condensation model suggested in [6] is used. In contrast to the ESH model, the total amount of latent heat can be computed from the change of enthalpy in the following waȳ Again, the condensation region is defined by where vol(Ω cond ) denotes the volume of Ω cond . With this, we get that the temperature increase due to condensation corresponds to the energy used to change the phase of the water, uniformly distributed over Ω cond . Finally, note that the numerical discretization of this model is described in Sect. 4.2.

Numerical methods
In this section, we detail the numerical methods used to discretize and solve the governing equations.

Numerical solution of the governing PDEs
The mathematical model for radiative heat transfer and the models for vaporization described above were used to simulate the behavior of ex-vivo porcine liver tissue during LITT. The computational geometry was generated using Open Cascade (Open Cascade SAS, Guyancourt, France) and the mesh was created with the help of GMSH, version 2.11.0 (cf. [14]). The governing equations were solved with the finite element method in Python, version 2.7, using the package FEniCS, version 2017.2 (cf. [15,16]). For the numerical solution of the PDEs, we first (semi-)discretize the bio-heat equation in time using the implicit Euler method. Then, we use piecewise linear Lagrange elements for the spatial discretization of the temperature and radiative energy. The resulting sequence of linear systems was then solved with the help of PETSc (cf. [17]), where we used the conjugate gradient method with a relative tolerance of 1 × 10 -10 . Afterwards, the damage function is computed using a right-hand Riemann sum to discretize the time integral of (6).

Discretization of the enthalpy model
In the following we describe our discretization of the enthalpy model. In particular, to compute the temperature distribution from time t to t + t we proceed as follows. We first solve (3) to obtain the radiative energy at t + t. With this, we compute the temperature distribution at t + t from (1). Subsequently, we iterate over the nodes of the finite element mesh and check, whether the temperature exceeds 100 • C. At these nodes, the temperature is set to 100 • C and from the excess temperature we compute the corresponding increase in enthalpy. If the enthalpy surpasses the limit of ρλ vap , we return this surplus in the form of heat to the corresponding nodes. After doing so, we integrate the (local) changes in enthalpy over Ω to compute the total change of enthalpy H. Therefore, we can now compute the source termQ vap of (11) as follows which is then used as the source term for the next time step, simulating the release of enthalpy by the condensation of the water vapor. Then, the new tissue damage is computed from (6) and the procedure is continued until we reach the end time τ .

Results and discussion
We use the experiments from the study of [4] to test the vaporization models. In this study LITT was applied to ex-vivo porcine livers and the resulting temperature was measured with a probe. The experiment was repeated nine times with different laser powers and different flow rates for the applicator cooling system. For the study in [4], the authors used the mathematical model introduced in Sect. 2 which was derived from the one presented in [2]. However, the model did not take into account the vaporization of water in the tissue. While the general agreement between experiment and simulation was good, there were notably two outliers, namely the cases P34F47 and P34F70, for which the highest laser power was used. For these cases, the simulated probe temperature would rise to well above 100 • C, while the measured probe temperature would reach a plateau below 100 • C. Therefore, in [4] the authors suspected that the missing vaporization model was the reason for this discrepancy. Now, we test this hypothesis by repeating the simulations with the previously introduced modified models that now include vaporization and condensation effects.

Experimental setting
For the validation of our models, we use the measurements from the experiments made in [4]. For these, livers were obtained from pigs which had been slaughtered approximately 6 hours prior to the experiment. The temperature of the samples was room temperature at the beginning of the experiments. A laser applicator from Somatex® Medical Technologies (Teltow, Germany) was placed into the middle of the liver sample. An optical fiber from the same company with a diffuser part of 3 cm at its tip was inserted into the applicator for delivering the laser energy from a Nd:YAG laser device (MY30; Martin Medizintechnik, Tuttlingen, Germany; wavelength 1064 nm) to the tissue. The applicator was equipped with a cooling water circulation system to protect the optical fiber and prevent the burning of tissue in close proximity to the applicator. A temperature probe was introduced into the porcine liver and placed close to the applicator in order to generate temperature measurements for validating the models of LITT. The setup for the nine test cases is shown in Table 1. The laser was applied with different powers, namely 22, 28, and 34 W, and different flow rates of the applicator cooling system. However, it is assumed that the effect of the coolant flow rate is negligible (cf. [4]). Furthermore, the position of the temperature probe is characterized by its radial distance d r to the applicator axis as well as its distance d z from the applicator tip, where the applicator itself is contained in the half space with z ≥ 0. We now simulate this experiment again using the vaporization models introduced previously, and compare the results with the measurement data as well as with the results obtained by the original model which does not consider vaporization. The optical, thermal, and damage parameters used for the simulation are listed in Table 2. They are taken from [4] and the references therein (cf. [18][19][20][21]). For the condensation region Ω cond we have chosen the points where the temperature was between T -= 60 • C and T + = 80 • C, as proposed in [6].

Figure 3
Comparison of the models for the case P34F47

The case P34F47
Let us start the investigation of the vaporization models with the test case P34F47 of [4], where a laser power of 34 W was used. The results for this case are shown in Fig. 3, where the measurement from the temperature probe, the results for the model of [4] and the results for both vaporization models of Sect. 3 are shown. For this specific case, the probe temperature simulated without a vaporization model rises well above 100 • C while the measured temperature reaches a plateau below 100 • C (see Fig. 3). In contrast, both vaporization models do not overestimate the temperature to the end of the treatment and predict the occurring plateau correctly. This is further visualized in Fig. 4, where the difference of the models from the measurement over the entire treatment is depicted. These results show that all models are reasonably close to the measured temperature until up to about 80 • C. After that point the model without vaporization overestimates the temperature significantly. The models that include vaporization give considerably better results since their predicted temperatures match the measured ones more closely throughout the whole treatment.

All nine test cases
After investigating the vaporization models in the context of the previous test case, where the original model without vaporization of [2,4] failed to predict the correct temperatures, we now investigate the other test cases from the study of [4]. The corresponding results are shown in Fig. 5, where the measured and simulated temperature at the probe is shown, and Fig. 6, which visualizes the difference of the simulated temperatures to the measurement. In general, the vaporization models are good in estimating the final temperature of the experiment. Especially for the cases P34F47 and P34F70, which could not be simulated accurately in [4], the vaporization model performs significantly better and does not overestimate the temperature to the end of the treatment. However, during the middle of the experiment the vaporization models tend to overestimate the temperatures. This can be seen, e.g., for the cases P22F70 and P28F70 (cf. Fig. 5). We suspect that the simple condensation model is responsible for this discrepancy as we explain in Sect. 5.4. Altogether, the ESH and the enthalpy model both show comparable but slightly different temperature curves. Especially the overestimation of the temperature during the middle of the experiment is usually higher for the enthalpy model. To compensate one could think about adjusting parameters, like the exact amount of water in the liver tissue. However, a first step should be to improve the simple condensation model.

Discussion of the simple condensation model
The consideration behind the simple condensation model described in Sects. 3.2 and 3.4 is solely to preserve the conservation of energy. Therefore, all the water which was vaporized at a certain time is assumed to instantaneously condensate in the condensation region This consideration is strictly global and does not involve any form of transport mechanism for the vapor. Hence, it is possible that vapor which was generated in one region can instantaneously condensate in another region. Through this mechanism temperature can be shifted from one region to another without any delay. This effect is possibly the reason for the overestimated temperatures during the middle of the experiment. This can be seen, e.g., for the case P28F70 (cf. Fig. 6), where all simulated temperatures are the same until about 400 s into the experiment. At that point the simulated temperatures rise much faster for the models that include vaporization than for the one without it. We suspect that at this point of the experiment, vaporization occurs in tissue close to the applicator. Due to the instantaneous transport mechanism of the simplified condensation model heat is then added to regions further away from the applicator, where the applicator is placed. This results in the non-physical temperature increase that can be seen in this case. Additionally the simple condensation model is also rather sensitive with respect to the choice of the condensation region as can be seen in Fig. 7, where the temperature at the probe for the case P34F47 is shown for the condensation region given by T -= 60 • C and T + = 80 • C in Fig. 7(a) as well as for T -= 70 • C and T + = 90 • C in Fig. 7(b).
To resolve this issue, the transport of vapor within the tissue must be taken into account. This could for example be done by adding an additional diffusion equation similar to the bio-heat equation to the state system. Therefore, an effective diffusion coefficient for the vapor must be known or estimated from measurements. Alternatively, a more complex solution would be to model the tissue as porous medium and to use a pressure based formulation for the vapor transport.

Conclusion
LITT is a minimally-invasive method in the field of interventional oncology used for treating liver cancer. Mathematical modeling and computer simulation are important features for treatment planning and simulating the necrosis of the tissue. The numerical simulation is in good agreement with temperature measurements for ex-vivo porcine liver. In particular, the incorporation of vaporization of water in liver tissue improves the simulation. Still a refinement of the simple and artificial model for the vapor might be necessary. Due to its global nature, this model allows for an undelayed flow of temperature from a hot region to a colder one. This is probably the reason for the overestimated temperatures during the middle of the experiments. An additional physically motivated transport meachanism for the vapor might be necessary. In order to use simulations for the monitoring and assistance during the treatment of humans it is important to model the blood perfusion, because blood vessels have a significant cooling effect. One approach can be to identify the blood perfusion rate from MR thermometry during the beginning of the treatment and use this information to correctly simulate the remaining treatment (cf. [22]).