Numerical study of the influence of water evaporation on radiofrequency ablation

Background Radiofrequency ablation is a promising minimal invasive treatment for tumor. However, water loss due to evaporation has been a major issue blocking further RF energy transmission and correspondently eliminating the therapeutic outcome of the treatment. Method A 2D symmetric cylindrical mathematical model coupling the transport of the electrical current, heat, and the evaporation process in the tissue, has been developed to simulate the treatment process and investigate the influence of the excessive evaporation of the water on the treatment. Results Our results show that the largest specific absorption rate (Q SAR ) occurs at the edge of the circular surface of the electrode. When excessive evaporation takes place, the water dehydration rate in this region is the highest, and after a certain time, the dehydrated tissue blocks the electrical energy transmission in the radial direction. It is found that there is an interval as long as 65 s between the beginning of the evaporation and the increase of the tissue impedance. The model is further used to investigate whether purposely terminating the treatment for a while allowing diffusion of the liquid water into the evaporated region would help. Results show it has no obvious improvement enlarging the treatment volume. Treatment with the cooled-tip electrode is also studied. It is found that the cooling conditions of the inside agent greatly affect the water loss pattern. When the convection coefficient of the cooling agent increases, excessive evaporation will start from near the central axis of the tissue cylinder instead of the edge of the electrode, and the coagulation volume obviously enlarges before a sudden increase of the impedance. It is also found that a higher convection coefficient will extend the treatment time. Though the sudden increase of the tissue impedance could be delayed by a larger convection coefficient; the rate of the impedance increase is also more dramatic compared to the case with smaller convection coefficient. Conclusion The mathematical model simulates the water evaporation and diffusion during radiofrequency ablation and may be used for better clinical design of radiofrequency equipment and treatment protocol planning.


Background
Radiofrequency ablation (RFA) is one of the promising therapeutic procedures of primary and secondary malignant tumor due to its minimal invasiveness, less side-effect and immunology stimulation as compared with other treatment modalities [1,2]. The high frequency alternating electrical current forces agitation of the ions in the tissue, resulting in frictional heat of the tissue and correspondent abrupt increase of temperature. Previous studies show that, when the local temperature reaches the boiling temperature of water 100°C, excessive evaporation takes place and increases the tissue impedance [3]. The roll-off has been a major obstacle for its application [4], especially for tumors larger than 3 cm [5][6][7][8]. Studies have been focused on improving therapeutic effect of RFA, either through better design of the device (for example, shape/geometry of electrodes, cooled-tip electrode) [9,10], or optimizing the RF treatment protocols [11,12].
Multiple probes, wet electrode and cooled-tip electrode are techniques proposed to enlarge the lesion size [13][14][15][16]. These designs have been proven to be able to improve the therapeutic outcome, however a better design and planning of the treatment protocols is still needed for clinical applications of RF. When multiple probes create enlarged coagulation volume, there may exist survived tumor cells in the intermediate region of the multiple lesions [15,16]. Design of the wet electrode enables perfusion of the hypertonic saline solution into the tumor region with uncontrollable and irregular enlarged treatment area [17,18]. Design of the cooled-tip electrodes allows circulation of the cooling agents inside the electrode, which helps lower the temperature on the contact surface, thus slows down local water loss rate, and eventually increases the coagulation volume up to 50% [19]. However there also exists regions with evaporation which blocks the RF energy transportation [20,21]. The existence of the dehydration of the tissue not only leads to the final formation of charring and decreases the electrical conductivity of the tissue, it also adsorbs large quantity of energy. Understanding how the evaporation process influences the transportation of the electrical currents [9] and thus adsorption of the energy will be the key issue to be able to further optimize the treatment protocols.
Mathematical modeling has been an efficient and low cost modality in both design and treatment planning [22]. Models have been developed to simulate the temperature distributions of tissue during the treatment [4,[23][24][25][26]. The complexity of the tumor composition, the influence of heterogeneous electrical conductivity, thermal conductivity, blood perfusion rate, and the specific distribution of blood vessels on the results have been studied [27][28][29][30][31][32][33][34][35].
However, most of the models have skipped the period when the excessive evaporation taking place during their simulation. Ai et al. [36] have incorporated the latent heat of evaporation into the specific heat capacity of the tissue in the bioheat transfer equation but diffusion of the water has been neglected. Yang et al. have proposed an empirical equation considering the influence of the evaporation process [37]. The dependence of the liver water content has been measured and a function describing this dependence established based on the experimental findings. The evaporation heat is further determined from the decreasing rate of liquid water fraction, and added as the heat sink in the heat transfer equation [38]. Both of these two models have ignored the influence of water loss on the electrical field and heat transfer as well. Patz et al. [39,40] has modeled the process of gas bubbles emerging, transportation, and merging with each other during the RF treatment. In their model, as soon as the tissue reached the boiling temperature, evaporation takes place, and an interface is immediately formed separating the gas region and the liquid region, which is inconsistent with Haemmerich's experiment findings [38].
In this study, a numerical model is developed investigating transportation of the electrical currents and heat inside the tissue during RF heating. Local boiling phenomena (evaporation) and diffusion of the water are both considered. These processes are coupled together. It is expected to help accurately predict deposition of RF energy before or after carbonization, and thus for effective prevention as needed. Also numerically studied is simulation of the treatment using different treatment protocols with the cooled-tip electrode.

Theoretical modeling
The electrodes of RFA used in this model are similar to those as described previously [41]. The geometry of the tissue and the electrode are illustrated in Figure 1. The treating surface of the electrode is a 10 mm diameter circle. The grounding pad is placed at the bottom of tissue. The tissue is assumed to be a symmetric cylinder with a diameter of 40 mm and a height of 40 mm, which is estimated from a typical Wistar mouse's height and the width. And the tissue is assumed to be homogenous. Treating face of the electrode closely contacts the upper surface of tissue. The tissue and the electrode are assumed to be axially symmetric cylinders, thus only changes in the axial and radial directions are considered. So it is simplified to a 2-D model, and cross-section of electrode and tissue is studied.
The RF frequency used is 460 kHz. As the wavelength of the electromagnetic field is much longer than the electrode, the quasi-static approximation of the electrical field in the tissue is applied [26,42], where V is the electrical potential, and σ t is the electrical conductivity of tissue. Neumann boundary condition is used to describe the insulated boundary of tissue and the electrode: The absorption rate of RF energy in tissue, Q SAR (W/m 3 ) is, Pennes' equation is used to describe the heat transfer inside the tissue [43], where ρ t , c t , k t stand for the density, heat capacity and thermal conductivity of the tissue; T b , w b , ρ b are the arterial blood temperature, blood perfusion rate, the density of blood; Q meta is the metabolic rate of tissue, and Q SAR is the specific absorption rate of tissue. A natural convection boundary condition is used for the outer surface of the tissue cylinder which includes the upper surface excluding the contact surface between the electrode and tissue, the side surface, and the bottom surface of the modeled tissue cylinder, as shown in Figure 1, where h s is the natural convection coefficient of air, and a value of 25 W/m 2 K is used. The blood flow stops when the local temperature is higher than 50°C [44], And the metabolic rate of tissue is lineally dependent on temperature [45], where Q meta0 is the metabolic rate in physiological state (T 0 = 37°C). With temperature rising, ions and polar molecules in the tissue oscillate faster, which helps heat transfer and electromagnet field propagate. The physical phenomenon reflects in the changes of the thermal and electrical conductivity [46]. Their dependence on temperature has been studied by Bhattacharya [47] and Pop [48] experimentally. And it can be described by a linear function with an independent variable-temperature, where T 0 is the reference temperature which is the initial temperature of tissue, usually the body temperature; σ 0 and k 0 are the electrical conductivity and thermal conductivity at T 0 , separately. As heating continues, temperature of the surface contacting with the electrode will reach the boiling point of water. According to Webster's experimental study [37], the excessive evaporation inside the tissue takes place between 100°C and 105°C. Therefore, the boiling temperature of pure water under 1 atm, 100°C, is used as the criterion for the occurrence of the evaporation inside the tissue. Thus according to energy conservation, where Q eva is the energy absorbed by evaporation per unit volume of tissue.
And the dynamic volume fraction of liquid water in tissue is, where m l is the mass of liquid water per unit volume of tissue, h fg is the latent heat of water, Q eva is the heat absorbed per unit volume of tissue by evaporation, and D l is the water diffusion coefficient, which exponentially increases with temperature according to Mill's [49] and Krynicki's study [49,50], The liquid water is considered as incompressible and the tissue volume assumed to be constant during the treatment. After evaporation, the space previously taken by the liquid water in tissue is filled with steam. The volume is occupied by a liquid-gas mixture. The liquid and gas fractions of water could be determined by: where ϕ g and ϕ w are the fractions of water in gas and liquid phase, m l0 is the initial mass of liquid water per unit volume of tissue. The modified equation describing the electrical conductivity of the two-phase (liquid and gas) system is used [51]: where A is a constant that primarily depends on the shape and the orientation of dispersed particles, for randomly distributed bubbles in tissue, A is 1.5; B is a constant which represents the relative conductivity of the two components, and the value is 2/3 for water and water vapor; ϕ m is the maximum packing fraction of particles, and it is 0.637 for randomly packed sphere gas bubbles. Thermal conductivity of the mixture is a linear function of water concentration [52]: The electrical and thermal conductivities of the tissue depend both on temperature and the water content, the influence of both factors are assumed to be independent of each other, thus Eq. (15) is proposed to describe its dependence on both temperature and water content: The equations describing the electrical field, heat transfer and water evaporation Eq.

Results and discussion
The parameters used in the model are listed in Table 1 The temperature distribution at t = 600 s is illustrated in Figure 2. It is found that the simulated highest temperature reaches 119.79°C if evaporation is not considered, while the highest temperature inside tissue is about 100.24°C when the evaporation process considered. If using 50°C as the critical temperature assessing cell viability, it can be found from the results that without considering the evaporation process, the 2-D lesion area will be 12.41% larger. Obviously, the model without considering evaporation overestimates the temperature profile in tissue and correspondently the tumor lesion size. Constant of mixture ϕ m 0.637 [51] Electrical conductivity σ 0 0.336 S/m [42,54] Blood perfusion (volume rate of blood flow per unit volume of tissue) w b 0.008 s -1 [42] Original water fraction in the tissue (mass of water per unit volume of tissue) m l0 778 kg/m 3 [38] Electrode voltage V

V
The transient temperature curves at the monitoring points are illustrated in Figure 3. The points are at the center line of the tissue cylinder with 0 mm, 1 mm, and 10 mm away from the contact surface. As seen from Figure 3, the temperatures at the first and the second monitor points which are closer to the contact surface (0 mm, 1 mm) keep rising before the temperature reaches a turning point, where the temperature starts to keep at a constant value 100°C, with the adsorbed RF energy being exhausted by evaporation. The temperature at the third point which is 10 mm away from the contact surface keeps rising during the whole RF heating process, and no evaporation takes place at this point. To find out how the evaporation process influences the transportation of RF energy, the electrical field E and Q SAR of tissue are obtained, and results are illustrated in Figure 4. The field of the specific heat adsorption rate Q SAR and the electrical field intensity E are established at the instant the RF treatment begins (Figure 4a). The largest electrical field intensity locates at the edge of the electrode contact surface, and the E-lines form semi-circle shapes with their centers at the edge of the electrode. After being heated for 600 s, though the locations with the largest value of Q SAR and the electrical field intensity are still at the edge of the electrode, the shapes of the E-lines are no longer semi circles (Figure 4b). Comparing with the distributions at the time point t = 0+, the electric field in the tissue after being heated for 10 minutes is more concentrated in the contacted area. And as shown in Figure 3, at this time point, excessive evaporation has taken place, it is the evaporation process that has reshaped the distribution of the Q SAR and electrical intensity E.
Also shown in Figure 4, the region with the largest Q SAR value is the place where the liquid water loses the most. With decrease of the water content in this region, the local electrical conductivity decreases, the water loss region forms a barrier for transportation of the electric currents in the radial direction. The trapped electrical field propagates more in the axial direction near the center line of the cylinder. This explains the unique shape of the E lines shown in Figure 4b.
Tissue impedance is one of the most frequently used clinical criteria to decide when to adjust RF power or to stop the treatment [55]. The tissue impedance during the treatment is also calculated with results shown in Figure 5. It slowly decreases at the beginning of the treatment, and after reaches the lowest point (the turning point), it increases quickly. According to the simulated water content results, the water evaporation emerges at about 334 s after the treatment starts, while the reflection point of the impedance increasing from its former decreasing trend starts at 399 s. There is a time interval of about 65 s between the occurrence of water evaporation and the inflection point of the impedance. This may due to the overwhelming influence of temperature on the tissue's electrical conductivity when the water content is not significantly decreased at the beginning of the heating.
According to the above simulation results, without adjusting the treatment protocol, there will form regions with serious water loss as barriers for RF energy delivery. However there is an interval between the occurrence of the evaporation and the beginning of the impedance increase. If at a certain time after the evaporation starts, terminating the heating allowing some time for the liquid water in the adjacent region to rehydrate the region may help enlarge the coagulation depth of the treatment. To find out its effectiveness, a heating process as illustrated in Figure 6a is studied. The first RF heating ends at 656 s as the impedance exceeds 500 ohm. Then a 5 min of natural cooling of tissue is allowed before restarting the RF heating. The second heating lasts for 344 s until the impedance exceeds 500 ohm again. During the natural cooling process, it can be seen that the tissue impedance first drops dramatically, and then increases slowly. The change of the impedance is clearly related to the water content inside the tissue. The temperature distribution at these time points (t = 656 s, t = 956 s, t = 1300 s) are given in Figure 6b-d. Unexpected, after the second heating (t = 1300 s), the heated region whose temperature is greater than 50°C is 28.5% less than that resulted from the first heating. This result suggests that by just waiting for the water to flow back to the dehydrated region is not an effective way improving the treatment outcome. During the intermittent thawing process when water flows back (about 64.3% rehydrated after 5 min, results not shown), the local temperature also drops ( Figure 6d). The second heating only increases the temperature of the same region to the previous point and fails to further propagate before the gas gap forms again and significantly increases the impedance. Continuously hydrating the water loss region of the tissue without decreasing the local temperature shall be more effective, adding hot sterile saline is an option if the process is well controlled.
To avoid local carbonization due to dehydration, electrodes with cooling agent flowing inside have been used [56,57]. Circulation of the cooling agent dissipates heat from the electrode wall and thus decreases the temperature in the adjacent region, which helps reduce local water evaporation. It has been proved effective in improving thermal ablation outcome of large tumor. To help optimize the treatment protocols, outcome of three different cooling strategies are predicted using the proposed model. The convective coefficient between the electrode wall and the circulating fluid h e is assumed to be 25 W/m 2 K, 100 W/m 2 and infinity. The temperature of the cooling agent is 20°C.
The temperature distribution and water content in the tissue during the treatment with these conditions are calculated and illustrated in Figure 7. It is found that during the treatment with the cooled-tip electrodes, there are certain locations inside the tissue whose temperature would reach the boiling point for all the cooling conditions studied, and under the first two conditions (with limited convection coefficient, there would form a complete enclosed gas gap in the tissue (shown in Figure 7a and b). This completely blocks further energy transmission of RF. It takes 205 s and 476.5 s for the complete gas gap to form for the first two cooling conditions, respectively. While for the infinite convective coefficient, that is, the temperature of the electrode wall is maintained at a constant value, thermal equilibrium achieves at 1400 s. Circulation of the cooling agents inside the electrode not only extends the therapy time but also the therapy area. The fraction of the tissue area with temperature higher than 50°C divided by the total tissue area in the 2D-section are 22.0%, 41.5% and 33.9% respectively under the three conditions. It is the largest for the process with moderate convection coefficient (100 W/m 2 K) of the cooling agent.
Another finding is that the flowing condition of the cooling fluids also alters the pattern of water loss. For the flow with a lower convection coefficient, the water loss first occurs at the edge of contact surface. While for the medium and extremely high convection coefficients, the water loss starts from the center of tissue at about 2.5 mm below the contact surface. The results suggest that the treatment region can be shaped by controlling the flowing conditions inside the electrode.
Further investigation of the transient temperature distribution during the treatment (Figure 8) has found that for some time after evaporation takes place, there is a decrease of the local temperature. The monitoring point in Figure 8a is the center of the contact surface. At the time when the temperature at this point starts to decrease, the gas vapor may have formed a barrier for the electrical field and also decreases the thermal conduction, while at the same time, the cooling agent continuously circulates, and it lowers the local temperature of tissue.
From Figure 8b, it can be found that impedance decreases slightly before the vapor point, and increases suddenly while gas forms a gap. Comparing Figure 8a   One interesting phenomenon observed is that, at the end of treatment, the impedance during the treatment with the medium convection coefficient increases more rapidly than that under a lower convection coefficient. As seen from Figure 8b, it costs 10.6 s for the impedance increasing from 115 ohm to 500 ohm under a larger convective coefficient, while it costs 25.6 s under the other condition. The result suggests that when cooled-tip electrode is applied, monitoring the rising of the impedance is necessary.
Besides, though the larger the convection coefficient, the longer the treatment time is, with further increase of the convection coefficient, a lot more heat is dissipated by the convection. When the cooling process exceeds the heating effect of the RF, the ablated size would decrease. There would be an optimal flow conditions inside the probe for a certain shaped ablation requirement. The model presented in this paper is expected to be useful in the treatment planning.
In devices developed in [41,58,59], low temperature nitrogen is used as cooling agent inside the treatment probe. As the nitrogen has a small thermal capacity, the temperature of the flow also changes, by coupling with the flow and heat transfer of the fluid, the model could be used for its better planning.
However, in the model, free diffusion of water is used, and the effect of the tissue structure on the diffusion would result in slower liquid water diffusion, and correspondently an earlier forming of a complete gas region. The pressure difference in tissue is also neglected. With further understanding of gas and water transport during the evaporation, the model will be more accurate to predict temperature above 100°C and dynamical impedance variation. With the patient information of tumor geometry and vessel distribution by imaging technique, it is possible to predict temperate distribution by RFA, and to optimize the placement of electrode.

Conclusions
In this study, a model incorporating evaporation, water diffusion and its influence on thermal and electrical conductivity has been developed to simulate the RF treatment. The distribution of the temperature, electrical field, Q SAR and water content in the tissue under different treatment conditions have been obtained numerically. The results show that the evaporation will not only maintain the highest temperature of tissue at 100°C, but also influence the distribution of the electrical field and the Q SAR . A time lag between the emerging of the evaporation and increase of the tissue impedance is found, which is consistent with others' experiment findings. It is also predicted that terminating the heating for a while to allow time for liquid water to flow back to the dehydrated region during the treatment would not enlarge the coagulation zone. The numerical results also confirm the advantage of the cooled-tip electrode. Increase in the convective coefficient of the circulating cooling agent inside the electrode would extend the treatment time, reduce the occurrence of the local evaporation and reshape the treated region, but not necessarily enlarge the size of the ablated region. Water loss emerges at different region of tissue under different convection condition. And the tissue impedance would increase more sharply before charring occurs under the medium convection condition. Thus, it is necessary to monitor the increase rate of impedance during the treatment.