Comparison of enthalpy method and water fraction method to mathematically model water vaporization during RF ablation

8 Journal of Advances in Biomedical Engineering and Technology, 2014, 1, 8-14 © 2014 Cosmos Scholars Publishing House Comparison of Enthalpy Method and Water Fraction Method to Mathematically Model Water Vaporization During RF Ablation Yuliia Fatieieva1, Pedro Almendárez2, Ricardo Romero-Méndez2, Enrique Berjano3 and Macarena Trujillo4,* 1Zaporizhzhye National University, Zaporizhzhye, Ukraine 2Facultad de Ingeniería, Universidad Autónoma de San Luis Potosí, SLP, Mexico 3Biomedical Synergy, Electronic Engineering Department, Universitat Politècnica de València, Spain 4Instituto Universitario de Matemática Pura y Aplicada, Universitat Politècnica de València, Spain Abstract: During high-temperature energy-based therapies such as radiofrequency ablation (RFA) the target tissue reaches temperatures around 100ºC, which causes tissue dehydration by water vaporization. In order to be as realistic as possible, RFA theoretical models should include the formulation of these phenomena. There are currently two fixed mesh methods of modeling the electrical and thermal effects produced by water vaporization: the enthalpy method and the water fraction method. Our objective was to compare both methods, especially to assess the thermal and electrical performance in terms of electrical impedance progress during heating, distributions of temperature, and temperature progress at some specific locations. The results showed the performance of both methods to be qualitatively analogous, with similar impedance progress, temperature distributions and temperature progress. They were hence equally able to mimic the thermal and electrical performance in a pulsed protocol, i.e. during the period without applying RF power. The main difference between the methods was the time at which impedance started to rise. All these findings suggest that the two methods offer equivalent results in RFA modeling. However, since the enthalpy method has one less problem to be solved (dynamic volume fraction of liquid water in the tissue) it is less complex, has a lower computational cost and therefore seems to be more suitable for modeling RFA with dry or internally cooled electrodes, i.e. those in which there is no interstitial saline infusion. However, the water fraction method would be more appropriate in the case of RFA with externally irrigated electrodes. Keywords: Enthalpy method, Mathematical modeling, Radiofrequency ablation, Tissue dehydration, Water vaporization. INTRODUCTION High-temperature thermal therapies are clinical procedures, which use some type of energy (laser, ultrasound, microwave or radiofrequency -RF-) to selectively destroy biological tissue by raising its temperature over 50ºC, or even to 100ºC. In the case of RF energy, this term could also include typical electrosurgical procedures such as cutting and coagulation, along with radiofrequency ablation (RFA), used to treat cancer in different organs [1-6]. RF power (100-200 W) is usually delivered to the tissue by creating a flow of electrical current (500 kHz) between a small active electrode in contact with the target tissue and a larger dispersive electrode on the patient’s back or thigh. Selective heating occurs around the active electrode in the target tissue, where the highest current density is located. When tissue is exposed to temperatures over 50ºC for several minutes, or 60ºC for a few seconds, the cells begin to die due to protein *Address correspondence to this author at the: Departamento de Matemática Aplicada, Escuela Técnica Superior de Arquitectura, Universidad Politécnica de Valencia Camino de Vera, 46022 Valencia Spain, Tel: +34963877007; Fax: +34963877669 ; E-mail: [email protected] denaturation and thermal lesions are created. In RFA, when the tissue temperature achieves a value 100ºC, there is a dehydration phenomenon due to water vaporization. Since the electrical current is completely dependent on tissue electrical conductivity (  ), and as  depends on the tissue water content, water vaporization during RFA is expected to significantly affect electrical performance, and in fact vaporization is known to influence the efficiency of RF power supply. Theoretical modeling and computational simulation have become powerful tools to study quickly and cheaply the electrical and thermal performance of RFbased ablative therapies [7]. As the worth of theoretical models increases with their realism [8], it is crucial to address the mathematical modeling of the vaporization phenomenon in high-temperature thermal therapies. Once the theoretical model has been set (geometry, governing equations, and initial and boundary conditions), numerical methods are usually needed to solve the problem due to the non-linearity of the equations and the presence of complex geometries. Both equations and tissue characteristics have to be suitably modified to model water vaporization. For this, Water Vaporization Modeling During RF Ablation Journal of Advances in Biomedical Engineering and Technology, 2014, Vol. 1, No.1 9 there are two methodological approximations: 1) fixed mesh and 2) moving mesh methods [9-10]. The main difference between them is that in the former the mesh remains constant in every step of the resolution, while in moving mesh the mesh varies according to the material phase change. The enthalpy method is a fixed mesh method usually employed to model water vaporization during heating [11-13]. Zhu et al. [14] recently proposed a new fixed mesh method that considers local boiling phenomena (vaporization) and water diffusion. Briefly, they considered a new variable ml, which is the mass of water by tissue volume, and which modifies the expressions describing the electrical and thermal conductivities. As a consequence, the electrical and thermal conductivities of the tissue depend on both temperature and water content. Introducing a new variable complicates the problem to be numerically treated since authors add the water content problem to the usual electrical-thermal coupled problem used in RFA. In fact, in this new method the equations describing the electrical field, heat transfer and water vaporization were solved by Comsol Multiphysics (COMSOL, Burlington MA, USA) while the thermal and electrical conductivity equations were solved by Matlab (Mathworks, Natick, MA, USA) and both results were then processed in Comsol. In contrast, the enthalpy method does not involve the introduction of new variables and hence is easier to solve. The addition of an additional problem in Comsol Multyphysics in the case of water fraction method implies an extra computational time. The question that arises is how the results of the enthalpy and water fraction methods will differ. Our objective was thus to compare the enthalpy method and the water fraction method proposed by Zhu et al. [14], especially to assess the thermal and electrical performance of both methods, identify the advantages and drawbacks of each, and to consider their applicability to the modeling of water vaporization in RFA. METHODS Description of the Model We used the same model geometry as that proposed by Zhu et al. [14] (Figure 1), which consists of a cylindrical active electrode in contact with the hepatic tissue surface. The tissue was also considered to be a cylindrical volume with larger diameter and height than the electrode. The dispersive electrode was modeled as an electrical condition on boundaries at a distance from the active electrode. The problem presented axial symmetry and hence a twodimensional analysis was possible. Figure 1: Geometry of the model used (out of scale). The cylindrical active electrode is 10 mm in diameter and is 10 mm high. It is placed on the surface of a fragment of hepatic tissue also of cylindrical volume, 40 mm in diameter and 40 mm high. We used a finer mesh in the interface between the electrode and tissue where a great temperature gradient was expected. The mesh size was automatically created by Comsol Multiphysics. We checked the suitable size of the mesh by a convergence test using the maximum temperature as a control parameter. The time step was also calculated for suitability by means of a convergence test, and a value of 0.1 s was obtained. Governing Equations The model is based on a coupled electric-thermal problem. The governing equation for the thermal problem was the bioheat equation:  cT t =(kT)+  bcb  b(TbT)+Qm+q (1) where  is density, c specific heat, T temperature, t time, k thermal conductivity,  b blood density, cb blood specific heat,  b blood perfusion coefficient, and Tb blood temperature. The value and dimensions of the constants used in the model are available in Table 1 [14]. It was assumed that once tissue temperature exceeds 50ºC there is no blood perfusion and thus  b = 0. The value of the blood perfusion coefficient was 0.008 s-1 up to 50ºC. The metabolic heat (Qm) was defined as a piecewise function: 10 Journal of Advances in Biomedical Engineering and Technology, 2014, Vol. 1, No.1 Fatieieva et al. Table 1: Material Parameters [14] Material  (S/m) k (W/mK)  (kg/m3) c (J/kgK) Electrode 7.4  106 15 8  103 480 1080 3455 Liver 0.336* 0.419 370** 2156** *Referenced at 37ºC; ** at gas phase. Qm=Qm0(1 +0.1(TT0)) T50º C 0T>50º C    (2) Qm0 being the metabolic rate in physiological state (33,8 kW/m3, at 37ºC) and T0=37ºC [14]. The heat source q was solved from the electrical problem with q=  E2 , where  is the electrical conductivity and E is an electric field, which is obtained from E=V where V is the voltage. This voltage was obtained from the Laplace equation   V=0 , which was the governing equation of the electrical problem. At RF frequencies (500 kHz) and over the distance of interest, the biological medium can be considered almost totally resistive and a quasi-static approach is therefore possible to solve the electrical problem [15]. This means that the electrical variables (V and E) are DC values in the simulations which correspond with the root-mean-square values of the RF signal. Boundary and Initial Conditions The boundary conditions for the thermal problem were: natural convection over all the outer tissue surfaces (thermal transfer coefficient of 25 W/m2K and bulk temperature of 20ºC, which would mimic a RF heating conducted at room temperature), excluding the contact surface between electrode and tissue. Natural convection (convective heat transfer coefficient of 0.1 W/m2K and bulk temperature of 20ºC) was also used on all outer electrode surfaces, which would mimic a thermal isolation condition. The initial temperature of electrode and tissue was 37ºC (body temperature). The boundary conditions for the electrical problem were: electrical insulation over all tissue surfaces, excluding the lower tissue surface, in which the condition was V=0 V (dispersive electrode) and V=35 V on all electrode surfaces. The model was solved numerically by the Finite Element Method using COMSOL Multiphysics software. Enthalpy Method to Model Water Vaporization In our study we only modeled the water vaporization using the enthalpy method [11]. In brief, its formulation assumes a change in the governing equation of the thermal problem (Eq. 1) as follows: (  h) t =(kT)+  b c b  b (T b T)+Q m +q (3) h being the enthalpy. For biological tissues the first term in Eq. (3) is computed as (  h) t =T t  lcl0T99º C H  m99 <T100º C  gcgT>100º C      (4) where the subscript l refers to the liquid tissue phase and g to the gas tissue phase (see Table 1), H is the product of the water latent heat and the percentage of water content in hepatic tissue (68%) and  m is the water density at 99ºC (958 kg/m3). In addition to the enthalpy method we used a piecewise function [4] for electrical conductivity. It is well-known that  increases with temperature up to 100ºC. Around 99-100ºC the water is lost by vaporization and  drops abruptly. We considered a  linear increase of 1.6%/ºC for T < 99ºC, which was the method used in [14]:  (T)=  0 (1 +0.016(T37)) 0 <T99°C 0.669312 99 <T100°C 0.669312 0.133849(T100) 100 <T105°C 0.66931210 4 T>105°C        (5) For thermal conductivity we considered a linear increase of 0.3%/ºC for T < 100ºC, which was the same assumption as that made in [14] and a constant value for T > 100ºC. Outcome Variables The thermal performance of both methods was compared using temperature distributions at different heating times. Temperature variations as used by Zhu et al. [14] were considered at three different locations on the central axis: 1 and 10 mm away from the electrode surface and on the surface itself (electrodetissue interface). The electrical performance of both methods was compared by impedance evolution. Rolloff time, i.e. when impedance suddenly rises, is a crucial phenomenon which denotes when the electrode is entirely surrounded by dehydrated tissue [16]. Finally, one of the issues assessed by Zhu et al. Water Vaporization Modeling During RF Ablation Journal of Advances in Biomedical Engineering and Technology, 2014, Vol. 1, No.1 11 [14] when they proposed the water fraction method was the potential of the model to mimic the rehydration process in a zone of the tissue previously dehydrated by water vaporization. For this they conducted a simulation in which RF power was cut off after impedance passed 500  and restarted 5 minutes later and allowed to exceed 500 . To compare both methods, we also reproduced this protocol with the enthalpy method. RESULTS We analyzed the same outcome variables as those reported by Zhu et al. [14]. Figure 2A shows the temperature distribution at t = 600 s obtained from the enthalpy method. For purposes of comparison, Figure 2B reproduces the temperature distribution at t = 600 s obtained by the water fraction method. Both distributions were found to be very similar, although the enthalpy method showed a higher temperature: note the hot points at the electrode edges, which are more pronounced than those in the water fraction method. The maximum tissue temperature at t = 600 s obtained with the enthalpy method was 103ºC, while Zhu et al. [14] reported a maximum value of 100ºC. Figure 2: Temperature distribution (scale in ºC) for t = 600 s obtained with the enthalpy method (A) and the water fraction method (B). (Image B is reproduced from Zhu et al. [14]). Figure 3 shows the temperature evolution obtained with both methods at three different points: 0, 1 and 10 mm away from the electrode contact surface, in the tissue on the central axis. In general, the temperature evolution was similar in both methods, although the temperature rose faster at points near the electrode in the enthalpy method (they reached the maximum temperature 25-30 s before). In contrast, the temperature at 10 mm from the electrode rose more slowly in the enthalpy method (around 50 s delay). In neither method did the maximum temperature exceed the vaporization point (100ºC). Figure 3: Temperature evolution at three points on the central axis of the tissue cylinder at 0, 1 and 10 mm from the contact surface of the electrode, computed by the enthalpy method (dashed lines) and by the water fraction method (solid lines). (The plots of the water fraction method are reproduced from Zhu et al. [14]). Figure 4 shows the impedance evolution obtained from both methods. The performance can be seen to be very similar up to the time when impedance starts to rise, which occurs sooner in the water fraction method. Roll-off time, defined as the time when impedance is 30  greater than the initial value, was at 600 s for the water fraction method and 705 s for the enthalpy method. Figure 4 also shows that once impedance started to rise (380 s in the water fraction method and 590 s in the enthalpy method) the rate of increase was much higher in the enthalpy method. Figure 4: Impedance evolution computed by the enthalpy method (dashed lines) and water fraction method (solid lines). (The plot of the water fraction method is reproduced from Zhu et al. [14]). We simulated the protocol used by Zhu et al. [14], in which RF power was cut off after impedance passed 500  and restarted 5 minutes later and allowed to exceed 500  once more. Figure 5 shows impedance evolution for both methods. In the water fraction method, the first phase of RF heating ended at 656 s and the second at 1300 s, after the 5 min cooling period [14]. This performance was repeated in the enthalpy method but for times of 751 s and 1456 s. 12 Journal of Advances in Biomedical Engineering and Technology, 2014, Vol. 1, No.1 Fatieieva et al. Figure 5: Impedance evolution computed by the enthalpy method (dashed lines) and water fraction method (solid lines). Point 1 is when the first heating phase stops due to impedance passing 500 ; point 2 represents the restart of the second RF phase (after 5 minutes of natural cooling of tissue with RF power switched off); and point 3 is the end of the second phase. (The plot of the water fraction method is reproduced from Zhu et al. [14]). Overall, impedance evolution was comparable in both methods; during the first phase impedance gradually dropped until the abrupt rise and then passed 500 , before the RF was switched off. During the 5minute period when RF power was not applied, impedance first fell rapidly and then gradually rose until RF power was again applied. During this second phase, impedance also fell and then abruptly passed 500  sooner than in the first phase. In these simulations it could be seen once more that the rate of increase was much higher in the enthalpy method (just before points 1 and 3 in Figure 5). In addition, the drop was faster in the enthalpy method once RF power had been switched off (just after point 1 in Figure 5). Zhu et al. [14] observed that just after the second phase of the water fraction method (t = 1300 s) the region where temperature over 50°C was 28.5% less than in the first heating phase. This was also observed in the enthalpy method, as shown in Figure 6, in which the region over 50°C was 20% less than in the first phase. These percentage differences correspond with the area calculated on the plane of each plot of Figure 6. DISCUSSION AND CONCLUSIONS Accurate modeling of energy-based thermal therapies should mimic the electrical and thermal phenomena caused by the high temperatures reached in them. Basically, temperatures around 100ºC produce water vaporization and hence tissue dehydration. This local water loss involves a drop in local electrical conductivity and hence a rise in total electrical impedance, which under certain circumstances can cause roll-off and the subsequent cutoff of RF power. The idea behind most methods proposed to model water vaporization is to include the latent heat associated with the phase transition from liquid to gas in the governing equation. In addition to this, most methods somehow take into account the impact of gas Figure 6: Temperature distributions at the end of the first (A) and second (B) heating phase in the enthalpy method and in the water fraction method (C and D). (The plot of the water fraction method is reproduced from Zhu et al. [14]). formation on the electrical and thermal conductivities. In the case of the enthalpy method used in this study, this is done by modifying the value of electrical conductivity once the temperature reaches 100ºC, specifically by reducing its value drastically by two or four orders of magnitude from 100 to 105ºC (Figure 7A). Conversely, in the water fraction method, electrical conductivity depends not only on temperature (in a similar way to the enthalpy method for temperatures below 100ºC), but also on the mass of liquid water per unit volume of tissue, and both factors are assumed to be independent of each other [14]. From the mathematical point of view, the main difference between the two methods is that the water fraction method adds a third problem (dynamic volume fraction of liquid water in the tissue) to electrical and thermal problems. This leads to the conclusion that with the water fraction method there is no two-way relation between electrical conductivity and temperature (as with the enthalpy method), but once tissue temperature reaches 100ºC Water Vaporization Modeling During RF Ablation Journal of Advances in Biomedical Engineering and Technology, 2014, Vol. 1, No.1 13 and consequently water vaporization occurs, electrical conductivity will take values determined not only by the local temperature but also by the local water content (see dashed lines in Figure 7B). Figure 7: Conceptual explanation of the effect of temperature (T) on electrical conductivity () in the enthalpy method (A) and water fraction method (B). In the enthalpy method there is a two-way relationship between  and T, described with a piecewise function as in Equation (5). Conversely, in the water fraction method, electrical conductivity depends not only on temperature (in a similar way to the enthalpy method for temperatures below 100ºC), but also on the mass of liquid water per unit volume of tissue, which impedes predicting the local value of  from the local temperature. For temperatures much lower than 100ºC the tissue can be expected to be completely rehydrated, so that the relation between  and T will be similar to that in the enthalpy method (conceptually represented as dashed lines converging on the solid line). The results of our comparative study showed a similar performance between both methods in terms of temperature distributions and impedance evolution. The maximum temperatures reached with both methods did not exceed the vaporization point (100ºC), which was due to the own modeling methods. This limit value was not due to neither the appearance of the roll-off (the maximum temperatures occurred before the starting point of the roll-off), nor the limited applied RF power (which was enough to reach 100ºC at many points). The main difference between the methods is the time at which impedance starts to rise. In the enthalpy method the increase in temperature happened 100 s later than in the water fraction method (Figure 4). The initial rise in impedance was similar in both methods and is related to the similar rise in temperature at points near the electrode (Figure 3). The enthalpy method was also able to mimic the same thermal and electrical performance as that obtained with the water fraction method during a pulsed protocol of applying RF power. Both the impedance evolution (Figure 5) and temperature distributions (Figure 6) were qualitatively similar in both methods. The most noticeable difference was that the rate of change in electrical impedance was much higher in the enthalpy method, both when rising (just before roll-off) and falling (just after roll-off). All these findings suggest that both methods could offer equivalent results in RFA modeling. From the viewpoint of the mathematical formulation, the enthalpy method has the advantage of not needing to add an extra problem (calculating the dynamic volume fraction of liquid water in the tissue), which means a lower computational cost. This is especially important when the RFA models quantify the thermal lesion from the computation of the thermal damage index, which is an additional function. The main advantage of the water fraction method is that it allows water content evolution to be modeled and thus computed at all points in the tissue, which could be very useful in modeling RFA with externally irrigated (wet) electrodes which infuse saline into the tissue. The interstitial infusion of hypertonic saline can modify the background electrical properties [17] and also allows desiccated tissue to be rehydrated once water vaporization occurs. This rehydration could be mathematically implemented by varying the dynamic volume fraction of liquid water, and the mathematical relation between electrical conductivity and temperature modeled by the water fraction method could show a performance with a type of hysteresis as schematized in Figure 7B (due to the simultaneous dependence on the water content), which has been observed in some experimental results [18]. We observed small differences in the absolute temperature values obtained from both methods (Figure 3). This could be due to the different meshing used to build the finite element models and is especially important in the case of circular electrodes in contact with the tissue surface, as it is known that the electrical current density is analytically infinite at the edges [19], so that its value in this zone in a finite element model will be highly influenced by mesh size. From a clinical application point of view, the results suggest that both the enthalpy method and the water fraction method could be employed to mimic water vaporization in some RFA models. Since the enthalpy method eliminates the problem of the dynamic volume fraction of liquid water in the tissue, it is less complex and has a lower computational cost and therefore seems to be a more suitable method of modeling RFA with dry or internally cooled electrodes, i.e. those without interstitial saline infusion. However, the water 14 Journal of Advances in Biomedical Engineering and Technology, 2014, Vol. 1, No.1 Fatieieva et al. fraction method seems to be more appropriate for RFA with externally irrigated electrodes. FINANCIAL SUPPORT This work received financial support from the Spanish “Plan Nacional de I+D+I del Ministerio de Ciencia e Innovación” Grant No. TEC2011-27133-C02-01. DECLARATION OF INTEREST None REFERENCES [1] Benoist S, Nordlinger B. Radiofrequency ablation in liver tumors. Ann Oncol 2004; 15:313-7. [2] McAchran SE, Lesani OA, Resnick MI. Radiofrequency ablation of renal tumors: past, present and future. Urology 2005; 66:15-22. [3] Di Staso M, Zugaro L, Gravina GL, Bonfili P, Marampon F, Di Nicola L, et al. A feasibility study of percutaneous radiofrequency ablation followed by radiotherapy in the management of painful osteolytic bone metastases. Eur Radiol 2011; 21:2004-10. [4] Sharma R, Wagner JL, Hwang RF. Ablative therapies of the breast. Surg Oncol Clin N Am 2011; 20:317-39. [5] Savoie PH, Lopez L, Simonin O, Loubat M, Bladou F, Serment G, et al. Two-years follow-up of radiofrequency thermotherapy for urination disorders due to benign prostatic hyperplasia. Prog Urol 2009; 19:501-6. [6] Akeboshi M, Yamakado K, Nakatsuka A, Hataji O, Taguchi O, Takao M, et al. Percutaneous radiofrequency ablation of lung neoplasms: initial therapeutic response. J Vasc Interv Radiol 2004; 15:463-70. [7] Berjano EJ. Theoretical modeling of epicardial radiofrequency ablation: State-of-the-art and challenges for the future. Biomed Eng Online 2006; 5:24. [8] Tungjitkusolmun S, Woo EJ, Cao H, Tsai JZ, Vorperian VR, Webster JG. Thermal-electrical finite element modeling for radio frequency cardiac ablation: effects of changes in myocardial properties. Med Biol Eng Comput 2000; 38:562-8. [9] Baines MJ, Hubbard ME, Jimack PK, Mahmood R. A movingmesh finite element method and its application to the numerical solution of phase-change problems. Commun Comp Physics 2009; 6:595-624. [10] Patz T, Preusser T. Composite finite elements for a phase change model. SIAM J Sci Comput 2012; 31:672-91. [11] Abraham JP, Sparrow EM. A thermal-ablation bioheat model including liquid-to-vapor phase change, pressureand necrosis-dependent perfusion, and moisture-dependent properties. Int J Heat Mass Transfer 2007; 50:2537-44. [12] Byeongman J, Aksan A. Prediction of the extent of thermal damage in the cornea during conductive thermokeratoplasty. J Therm Biol 2010; 35:167-74. [13] Trujillo M, Berjano E. Review of the mathematical functions used to model the temperature dependence of electrical and thermal conductivities of biological tissue in radiofrequency ablation. Int J Hyperthermia 2013; 29:590-7. [14] Zhu Q, Shen Y, Zhang A, Xu LX. Numerical study of the influence of water evaporation on radiofrequency ablation. Biomed Eng Online 2013, 12:127. [15] Doss JD. Calculation of electric fields in conductive media. Med Phys 1982; 9:566-73. [16] Trujillo M1, Alba J, Berjano E. Relationship between roll-off occurrence and spatial distribution of dehydrated tissue during RF ablation with cooled electrodes. Int J Hyperthermia 2012; 28:62-8. [17] Solazzo SA, Liu Z, Lobo SM, Ahmed M, Hines-Peralta AU, Lenkinski RE, et al. Radiofrequency ablation: importance of background tissue electrical conductivity--an agar phantom and computer modeling study. Radiology 2005; 236:495-502. [18] Esrick MA, McRae DA. The effect of hyperthermia-induced tissue conductivity changes on electrical impedance temperature mapping. Phys Med Biol 1994; 39:133-44. [19] Wiley JD, Webster JG. Analysis and control of the current distribution under circular dispersive electrodes. IEEE Trans Biomed Eng 1982; 29:381-5. Received on 16-07-2014 Accepted on 31-07-2014 Published on 17-10-2014 © 2014 Fatieieva et al.; Licensee Cosmos Scholars Publishing House. This is an open access article licensed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/), which permits unrestricted, non-commercial use, distribution and reproduction in any medium, provided the work is properly cited.