Abstract
Radiofrequency (RF) ablation with internally cooled needle-like electrodes is widely used in medical techniques such as tumor ablation. The device consists of a metallic electrode with an internal liquid cooling system that cools the electrode surface. Theoretical modeling is a rapid and inexpensive way of studying different aspects of the RF ablation process by the bioheat equation, and the analytical approach provides an exact solution to the thermal problem. Our aim was to solve analytically the RF ablation transient time problem with a needle-like internally cooled cylindrical electrode while considering the blood perfusion term. The results showed that the maximal tissue temperature is reached ≈3 mm from the electrode, which confirms previous experimental findings. We also observed that the temperature distributions were similar for three coolant temperature values (5, 15 and 25). The differences were only notable in temperature very close to the probe. Finally, considering the 50 line as a thermal lesion mark, we found that lesion diameter was around 2 cm, which is exactly that observed experimentally in perfused hepatic tissue and slightly smaller than that observed in nonperfused (ex vivo) hepatic tissue.
1. Introduction
Radiofrequency (RF) ablation with internally cooled needle-like electrodes is widely used for medical techniques such as tumor ablation [1], treatment of autonomously functioning thyroid nodules [2], and cardiac ablation to cure arrhythmias [3]. The device consists of an internally liquid cooled metallic electrode that cools the electrode surface (see Figure 1(a)). Theoretical models provide a fast and inexpensive way of studying different aspects of the RF ablation process with the bioheat equation [4]. Unlike numerical solutions such as those based on the finite element method, the analytical approach provides an exact solution to the thermal problem in a simplified scenario. In the case of a needle-like electrode, the model is based on an ideal conductor of infinite length totally immersed in homogeneous tissue (see Figure 1(b))). As far as we know, only Haemmerich et al. [5] have developed an analytical model of RF ablation with cooled-needle cylindrical electrodes, but this model had two important limitations: (1) it only considered the steady-state solution and (2) did not include the blood perfusion term, which is crucial in RF ablation of well-perfused organs such as the liver. A similar study including the blood perfusion term in the steady-state case was developed by Yue et al. in [6]. The aim of our study was thus to solve analytically the transient time problem of RF ablation with a needle-like internally cooled cylindrical electrode in well-perfused tissue, that is, considering the blood perfusion term, which is essential in the bioheat equation. As far as we are aware, this is the first full analytical solution obtained for this problem.
(a)
(b)
2. Governing Equations
Temperature distribution in the tissue during RF ablation is mathematically obtained by solving the Pennes bioheat equation [4]: where , and are the density, specific heat, and thermal conductivity of the tissue; respectively, , and are the density, specific heat, and perfusion coefficient of the blood, is the blood temperature, and represents the heat sources. Assuming all quantities , and to be constant and that the heat source is independent on the polar angle and working in cylindrical coordinates , (2.1) becomes
The heat source used in the study () was the electrical power density (W/m3) dissipated into the tissue throughout the temporal interval . Due to the cylindrical geometry of the theoretical model, we employed the same formulation used by Haemmerich et al. [5], and hence we have where is the current density at the conductor surface, is the electrical conductivity, and is the electrode radius.
The initial and boundary conditions in the case of an internally cooled electrode are where is the temperature fixed by the refrigeration of the electrode.
We change to dimensionless variables (where ) which leads us to the problem (where we have defined ) with the following initial and boundary conditions: (Dirichlet’s boundary condition).
3. Preliminary Results about Bessel’s Functions
To solve the initial boundary value problem (2.8), (2.9), (2.10), (2.11), we shall need to use deep well-known properties of Bessel’s and modified Bessel functions of complex argument, which we present now to facilitate reading the paper. First we recall the expression of the modified Bessel functions and of first and second class and order 0 and the expression of the Bessel function of second class and order 0 ( is the Euler-Mascheroni constant) which implies that where and are even holomorphic functions on . Remark that Moreover, we shall need the following relations:
Finally we recall the asymptotic expansions for of and the modified Bessel functions of first class and integer order (see [7], 7.2.3), which will be necessary for delicate computations in the following sections: for every .
4. Resolution of the Initial-Boundary Value Problem
Taking Laplace’s transforms with respect to and using (2.9), (2.10), and (4.2), we obtain The homogeneous equation associated to (4.1) is a modified Bessel equation of order 0 with general solution Since the Wronskian determinant of and is (formula (19) of 3.71 in [7]), the method of variation of parameters gives us
So we take where and are functions independent on to be chosen in such a way that (4.2) and (4.3) hold.
First, we remark the expected but nontrivial fact that In fact, by L’Hôpital’s rule and formula (7) in Section 3.71 of [7], we have Using the asymptotic expansions of , and given in Section 7.23 in [7], we obtain and (4.7) follows easily. Analogously,
As a consequence of (4.7) and (4.10), we need to choose Then, to satisfy (4.3), we obtain obtaining finally
To compute the inverse Laplace transform , we proceed in several steps. All the involved functions have a branch point in , so we will use the Bromwich’s contour of Figure 2 where denotes a fixed positive real number such that all the used functions are holomorphic on the set . We denote for subsequent use .
4.1. Computation of
As in the integral we begin computing, for every ,
As , it follows easily from the asymptotic expansions (3.8) and (3.9) that if and is large enough in order that , and and, moreover, and , we have the estimation because and . A similar result holds if by (3.10). Then we can to apply Bromwich’s formula to find .
The residue in the pole is On the other hand, by (3.1) we obtain Hence, Having in mind (3.7) and (3.3), we obtain because the imaginary part of with and and has limit if approach to the real axis. In definitive, after the change ,
4.2. Computation of
To find we have , and hence, as a consequence of the previous result, we obtain
4.3. Computation of
Now we compute By (3.8), if , and is large enough in order that and , we have (because ) and so Bromwich’s inversion formula cannot be used directly to find . To circumvent this complication we need to proceed in the way finding this last inverse with means of Bromwich’s contour of Figure 2 since if is large enough.
First, we remark that it follows easily from (3.3) that On the other hand, the function has a pole of order 2 in . As the function has continuous partial derivatives of any order in some neighbourhood of and every fixed , we can write
Finally, since is an even function, by (3.3), (3.5), and (3.4), we have
As a consequence, putting , we obtain from (4.26) that
4.4. Computation of
To compute we have and so, in order to find using the convolution theorem and our previous results, we can write
4.5. Complete Solution
Collecting our previous results, by the second translation theorem, Fubini’s theorem, and some natural computations, we obtain finally that can be rewritten in the following way which is more suitable for numerical computations:
It is interesting to check that, as expected, actually we have at every time . In fact, as if , in the computation of , we can use Fubini’s theorem and obtain explicitly the integration with respect to , which is Finally, after the variable change and the application of formula (5) in [7], Section 13.53, we obtain and this gives us .
5. Results and Discussion
Once the solution was achieved, we observed that it was hard to make a direct plot of the temperatures with Mathematica 6.0 software (Wolfram Research Inc., Champaign, IL, USA). The computer took around 24 hours for each plot at a fixed time .
Although we had continuously employed dimensionless variables in the analytical solution of the problem, as here we considered the case of liver RF ablation, we used the following values for the hepatic tissue characteristics: density of 1060 kg/m3, specific heat of 3600 J/kg·K, thermal conductivity of 0.502 W/m·K, and electrical conductivity of 0.33 S/m [8]. Since this was a case of well-perfused tissue, we considered the following blood characteristics: density of 1000 kg/m3, specific heat of 4148 J/kg·K, and a perfusion rate of [9]. Blood temperature, and hence initial tissue temperature , was 37°C. For all the simulations, we used a current density of value 5 mA/mm2.
In order to assess the effect of different coolant temperatures on the temperature profile, we used different boundary temperatures on the electrode surface . Figure 3 shows the temperature profiles at different times for a coolant temperature of 5°C. Likewise Figures 4 and 5 show the same plots for coolant temperatures of 15°C and 25°C, respectively.
The results showed that the maximal temperature in the tissue is reached ≈3 mm from the electrode (see Figures 3–5), which confirms previous experimental findings [10]. We also observed that the temperature was rising until achieving a steady-state at infinite time (thick line in Figures 3–5).
We also observed that the temperature distributions were similar for the three values of coolant temperature (5°C, 15°C, and 25°C). The differences were only significant at temperatures very close to the probe. This finding also agrees with previous experimental results in which little difference was observed in lesion size when coolant temperature was varied [5]. In this respect, if we considered the 50°C line as a thermal lesion mark [11] in Figures 3–5, the lesion diameter would be around 2 cm, which is exactly that observed experimentally in perfused hepatic tissue [11] and slightly smaller than that observed in nonperfused (ex vivo) hepatic tissue [5].
Our results were achieved with a current density of 5 mA/mm2, and it seems obvious that higher current density would cause bigger lesions. However, once the tissue temperature reaches 100°C, the charred tissue around the electrode creates a highly resistive electrical interface, which impedes further power deposition in the tissue [5]. Since our analytical solution is not able to model the nonlinear processes involved in these phenomena, we chose a value of current density in our simulations to keep the tissue temperature always lower than 100°C. This was also used in the previous analytical model by Haemmerich et al. [5]. Temperature dependence of tissue electrical and thermal conductivity was not taken into account. Taking this into account would make the problem nonlinear and most likely impossible to solve analytically. However, in real RF ablation experiments, this effect is present as well as phase transition, charring, and other hard to model effects. Therefore, although comparison of analytical solution to experimental results was favorable, future studies using numerical methods, such as finite element method, should be conducted. However, it is necessary to point out that the solution in (4.37) is exact, while that the solution provided by FEM is an approximation.
Since it is known that the density current pattern around a needle-like electrode is highly heterogeneous [8, 9], the value of 5 mA/mm2 chosen for our simulations cannot be related with the values of current usually employed in RF liver ablation (1-2 A). In spite of this, if we consider a 30 mm long and 0.75 mm radius electrode, the value of 5 mA/mm2 provides a current total of 0.7 A, which is a bit smaller than the values experimentally observed.
6. Conclusion
We have solved analytically the transient time problem of RF ablation with a needle-like internally cooled electrode by using the bioheat equation (i.e., considering the blood perfusion term). The temperature distributions computed from the theoretical model matched the experimental results obtained in previous studies, which suggests the utility of the model and its analytical solution to study the thermal performance of this kind of electrode.
Acknowledgments
This work received financial support from the Spanish “Plan Nacional de I + D + i del Ministerio de Ciencia e Innovación” Grant no. TEC2008-01369/TEC. The translation of this paper was funded by the Universitat Politècnica de València, Spain.