Analytical Model Based on a Cylindrical Geometry to Study RF Ablation with Needle-Like Internally Cooled Electrode

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 ≈3mm from the electrode, which confirms previous experimental findings. We also observed that the temperature distributions were similar for three coolant temperature values 5◦C, 15◦C, and 25◦C . The differences were only notable in temperature very close to the probe. Finally, considering the 50◦C 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.


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  Figure 1: a Cluster of three internally cooled needle-like electrodes used to ablate biological tissues by means of RF currents. b Analytical model representing a simplified scenario of an ideal conductor with infinite length totally immersed in homogeneous tissue. The internal cooling was modeled by means of a Dirichlet's thermal boundary condition with a constant temperature T C which corresponds with the coolant temperature inside the electrode. Accordingly, the theoretical model has one dimension 1D , that is, axis r.
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 wellperfused 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.

Governing Equations
Temperature distribution in the tissue during RF ablation is mathematically obtained by solving the Pennes bioheat equation 4 : Mathematical Problems in Engineering   3 where η, c, and k are the density, specific heat, and thermal conductivity of the tissue; respectively, η b , c b , and ω b are the density, specific heat, and perfusion coefficient of the blood, T b is the blood temperature, and S represents the heat sources. Assuming all quantities η, η b , c, c b , k, and ω b to be constant and that the heat source is independent on the polar angle θ and working in cylindrical coordinates r, θ , 2.1 becomes The heat source used in the study S was the electrical power density W/m 3 dissipated into the tissue throughout the temporal interval 0, ∞ . 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 j 0 is the current density at the conductor surface, σ is the electrical conductivity, and r 0 is the electrode radius. The initial and boundary conditions in the case of an internally cooled electrode are where T C is the temperature fixed by the refrigeration of the electrode. We change to dimensionless variables where α k/cη which leads us to the problem Mathematical Problems in Engineering where we have defined β : η b c b ω b r 2 0 /k with the following initial and boundary conditions: Dirichlet's boundary condition .

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 I 0 z and K 0 z of first and second class and order 0 Finally we recall the asymptotic expansions for |z| → ∞ of K 0 z and the modified Bessel functions of first class and integer order I ν z see 7 , 7.2.3 , which will be necessary for delicate computations in the following sections: 3.9 3.10 for every ν ∈ N ∪ {0}.

Resolution of the Initial-Boundary Value Problem
Taking Laplace's transforms D ρ, s, β : L V ρ, ξ ρ, s, β 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 I 0 z and K 0 z is W I 0 z K 0 z −1/z formula 19 of 3.71 in 7 , the method of variation of parameters gives us Mathematical Problems in Engineering So we take where M 1 s and M 2 s 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 I 0 x , I 1 x , and K 0 x given in Section 7.23 in 7 , we obtain and 4.7 follows easily. Analogously,

4.10
As a consequence of 4.7 and 4.10 , we need to choose Mathematical Problems in Engineering 7 Then, to satisfy 4.3 , we obtain

4.13
To compute the inverse Laplace transform L −1 D ρ, s , we proceed in several steps. All the involved functions have a branch point in s −β, so we will use the Bromwich's contour of Figure 2 where γ 0 denotes a fixed positive real number such that all the used functions are holomorphic on the set Re s ≥ γ 0 . We denote for subsequent use D β : C\ − ∞, −β .

4.19
Having in mind 3.7 and 3.3 , we obtain

4.20
because the imaginary part of log z with z ∈ D β and Re z < β and Im z > 0 has limit πi if z approach to the real axis. In definitive, after the change −y x,

Computation of
we have v ≤ ρ, and hence, as a consequence of the previous result, we obtain

4.25
because cos Arg s β /2 ≥ 0 and so Bromwich's inversion formula cannot be used directly to find F 3 ρ, ξ . 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

4.30
Finally, since R z is an even function, by 3.3 , 3.5 , and 3.4 , we have

12 Mathematical Problems in Engineering
As a consequence, putting s −x, we obtain from 4.26 that 4.32

Computation of
we have v ≥ 1 and so, in order to find using the convolution theorem and our previous results, we can write

4.35
Mathematical Problems in Engineering 13

Complete Solution
Collecting our previous results, by the second translation theorem, Fubini's theorem, and some natural computations, we obtain finally

4.36
that can be rewritten in the following way which is more suitable for numerical computations:

4.37
It is interesting to check that, as expected, actually we have V 1, ξ −B at every time ξ > 0. In fact, as F 3 1, ξ 1 if ξ > 0, in the computation of V 1, ξ , we can use Fubini's theorem and obtain explicitly the integration with respect to w, which is Finally, after the variable change z x − β and the application of formula 5 in 7 , Section 13.53, we obtain and this gives us V 1, ξ −B.

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 t.
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/m 3 , specific heat c of 3600 J/kg·K, thermal conductivity k 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 η b of 1000 kg/m 3 , specific heat c b of 4148 J/kg·K, and a perfusion rate ω b of 6.410 −3 s −1 9 . Blood temperature, and hence initial tissue temperature T b , was 37 • C. For all the simulations, we used a current density j 0 of value 5 mA/mm 2 .
In order to assess the effect of different coolant temperatures on the temperature profile, we used different boundary temperatures on the electrode surface T C . Figure 3 shows the temperature profiles at different times for a coolant temperature of 5 • C. Likewise The results showed that the maximal temperature in the tissue is reached ≈3 mm from the electrode see  , 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  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/mm 2 , 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/mm 2 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/mm 2 provides a current total of 0.7 A, which is a bit smaller than the values experimentally observed.

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.