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

Authors: Juan A. López Molina; María J. Rivera; Enrique Berjano
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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∘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.


---

## Body

## 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) 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(TC) which corresponds with the coolant temperature inside the electrode. Accordingly, the theoretical model has one dimension (1D), that is, axis r.
(a)(b)


## 2. Governing Equations

Temperature distribution in the tissue during RF ablation is mathematically obtained by solving the Pennes bioheat equation [4]:(2.1)ηc∂T∂t=∇⋅(k∇T)+S-ηbcbωb(T-Tb),
where η,c, and k are the density, specific heat, and thermal conductivity of the tissue; respectively, ηb,cb, and ωb are the density, specific heat, and perfusion coefficient of the blood, Tb is the blood temperature, and S represents the heat sources. Assuming all quantities η,ηb,c,cb,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(2.2)ηc∂T∂t(r,t)=k(∂2T∂r2(r,t)+1r∂T∂r(r,t))+S(r,t)-ηbcbωb(T(r,t)-Tb).The heat source used in the study (S) was the electrical power density (W/m3) 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(2.3)ηc∂T∂t(r,t)=k(∂2T∂r2(r,t)+1r∂T∂r(r,t))+j02r02σr2-ηbcbωb(T(r,t)-Tb),
where j0 is the current density at the conductor surface, σ is the electrical conductivity, and r0 is the electrode radius.The initial and boundary conditions in the case of an internally cooled electrode are(2.4)T(r,0)=Tb,∀r>r0,(2.5)lim⁡r→∞T(r,t)=Tb,∀t>0,(2.6)T(r0,t)=TC,∀t>0,
where TC is the temperature fixed by the refrigeration of the electrode.We change to dimensionless variables(2.7)ρ=rr0,ξ=αtr02,V(ρ,ξ)=σkj02r02(T(r0ρ,r02ξα)-Tb),
(where α=k/cη) which leads us to the problem(2.8)-(∂2V∂ρ2+1ρ∂V∂ρ)+∂V∂ξ+βV=1ρ2,
(where we have defined β∶=ηbcbωbr02/k) with the following initial and boundary conditions:(2.9)V(ρ,0)=0,∀ρ>1,(2.10)lim⁡ρ→∞V(ρ,ξ)=0,∀ξ>0,ξ>0,(2.11)V(1,ξ)=σk(TC-Tb)j02r02:=(-B),
(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 I0(z) and K0(z) of first and second class and order 0(3.1)I0(z)=∑m=0∞z2m22mm!Γ(m+1),K0(z)=-I0(z)log⁡(z2)+∑m=0∞1+(1/2)+⋯+(1/m)-γ22m(m!)2z2m,
and the expression of the Bessel function of second class and order 0(3.2)Y0(z)=2π(J0(z)log⁡(z2)-∑m=0∞(-1)m1+(1/2)+⋯+(1/m)-γ22m(m!)2z2m),
(γ is the Euler-Mascheroni constant) which implies that(3.3)K0(z)=-I0(z)log⁡(z2)+Z(z),(3.4)Y0(z)=2π(log⁡(z2)J0(z)-R(z)),
where Z(z) and R(z) are even holomorphic functions on ℂ. Remark that(3.5)Z(zi)=R(z).
Moreover, we shall need the following relations:(3.6)K0′(z)=-K1(z),(3.7)∀z∈C,I0(z)=J0(zi).Finally we recall the asymptotic expansions for|z|→∞ of K0(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.8)K0(z)~π2ze-z(1+O(1z)),if|Arg(z)|≤3π2,(3.9)Iν(z)~ez2πz(1+4ν2-12222z+(4ν2-12)(4ν2-32)242!(2z)2+O(1z3))+e-z+(ν+(1/2))πi2πz(1+O(1z)),if-π2<Arg(z)<3π2,(3.10)Iν(z)~ez2πz(1+4ν2-12222z+(4ν2-12)(4ν2-32)242!(2z)2+O(1z3))+e-z-(ν+(1/2))πi2πz(1+O(1z)),if-π2<Arg(z)<3π2,if-3π2<Arg(z)<π2,
for every ν∈ℕ∪{0}.


## 4. Resolution of the Initial-Boundary Value Problem

Taking Laplace’s transformsD(ρ,s,β):=𝔏[V(ρ,ξ)](ρ,s,β) with respect to ξ and using (2.9), (2.10), and (4.2), we obtain(4.1)ρ2d2Ddρ2+ρdDdρ-(s+β)ρ2D=-1s,(4.2)lim⁡ρ→∞D(ρ,s,β)=0,(4.3)D(1,s,β)=-Bs.
The homogeneous equation associated to (4.1) is a modified Bessel equation of order 0 with general solution(4.4)D(ρ,s,β)=C1(s)I0(ρs+β)+C2(s)K0(ρs+β).
Since the Wronskian determinant of I0(z) and K0(z) is W(I0(z)K0(z))=-1/z (formula (19) of  3.71 in [7]), the method of variation of parameters gives us(4.5)C1′(s)=-1sK0(s+βρ)ρ,C2′(s)=1sI0(s+βρ)ρ.So we take(4.6)C1(s)=-1s∫1ρK0(s+βv)vdv+M1(s),C2(s)=1s∫1ρI0(vs+β)vdv+M2(s),
where M1(s) and M2(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(4.7)L1:=lim⁡ρ→∞I0(ρs+β)∫ρ∞K0(vs+β)vdv=0.
In fact, by L’Hôpital’s rule and formula (7) in Section  3.71 of [7], we have(4.8)L1=lim⁡ρ→∞-1ρs+βK0(ρs+β)I02(ρs+β)I1(ρs+β).
Using the asymptotic expansions of I0(x),I1(x), and K0(x) given in Section  7.23 in [7], we obtain(4.9)lim⁡x→∞K0(x)I02(x)I1(x)=lim⁡x→∞2πxe2x(1+O(1/x))(1+O(1/x))22πxexex(1+O(1/x))=0,
and (4.7) follows easily. Analogously,(4.10)L2:=lim⁡ρ→∞K0(ρs+β)∫1ρI0(vs+β)vdv=lim⁡ρ→∞I0(ρs+β)K02(ρs+β)ρs+βK1(ρs+β)=lim⁡x→∞e-x(1+O(1/x))(1+O(1/x))22x2e-x(1+O(1/x))=0.As a consequence of (4.7) and (4.10), we need to choose(4.11)M1(s)=1s∫1∞K0(vs+β)vdv.
Then, to satisfy (4.3), we obtain(4.12)M2(s)=-BsK0(s+β)-I0(s+β)sK0(s+β)∫1∞K0(vs+β)vdv,
obtaining finally(4.13)D(ρ,s,β)=1sI0(ρs+β)∫ρ∞K0(vs+β)vdv+1sK0(ρs+β)∫1ρI0(vs+β)vdv-BK0(ρs+β)sK0(s+β)-I0(s+β)K0(ρs+β)sK0(s+β)∫1∞K0(vs+β)vdv=∫ρ∞L-1[f1](ρ,ξ,v)vdv+∫1ρL-1[f2](ρ,ξ,v)vdv-BL-1[f3](ρ,ξ,v)-∫1∞L-1[f4](ρ,ξ,v)vdv.To compute the inverse Laplace transform𝔏-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 𝒟β:=ℂ∖]-∞,-β].Figure 2
Bromwich’s integration contour.


### 4.1. Computation of𝔏-1[f1](ρ,ξ,v)

Asv≥ρ in the integral(4.14)∫ρ∞L-1[f1](ρ,ξ,v)vdv,
we begin computing, for every v≥ρ,(4.15)F1(ρ,ξ,v):=L-1[f1](ρ,ξ,v):=L-1[I0(ρs+β)K0(vs+β)s](ρ,ξ,v).Asv≥ρ, it follows easily from the asymptotic expansions (3.8) and (3.9) that if s∈𝒟β and |s+β| is large enough in order that |O(1/ρs+β)|≤1/2,|O(1/vs+β)|≤1/2, and |β/s|≤1/2 and, moreover, Re(s)≤γ0 and Im(s)>0, we have the estimation(4.16)|I0(ρs+β)K0(vs+β)s|≤9(|e-(v-ρ)s+β|+|e-(v+ρ)s+β|)8vρ|s||s+β|≤9(|e-(v-ρ)|s+β|cos⁡(Arg(s+β)/2)|+|e-(v+ρ)|s+β|cos⁡(Arg(s+β)/2)|)8vρ|s||s+β|≤92vρ|s|3/2,
because cos(Arg(s+β)/2)≥0 and |1+β/s|≤1-1/2=1/2. A similar result holds if Im(s)<0 by (3.10). Then we can to apply Bromwich’s formula to find 𝔏-1[f1](ρ,ξ,v).The residue in the poles=0 is(4.17)R(s=0)=lim⁡s→0sesξI0(ρs+β)K0(vs+β)s=K0(vβ)I0(ρβ).
On the other hand, by (3.1) we obtain(4.18)lim⁡s→-β,s∈Dβ(s+β)esξI0(ρs+β)K0(vs+β)s=0.
Hence,(4.19)L-1[f1]=K0(vβ)I0(ρβ)-12πi∫-∞-βesξsK0(v-s-βi)I0(ρ-s-βi)ds+12πi∫-∞-βesξsK0(-v-s-βi)I0(-ρ-s-βi)ds.
Having in mind (3.7) and (3.3), we obtain(4.20)L-1[f1]=K0(vβ)I0(ρβ)-12πi∫-∞-βπi2eyξyJ0(ρ-y-β)(J0(v-y-β)+iY0(-v-y-β))dy+12πi∫-∞-βπi2eyξyJ0(ρ-y-β)(J0(v-y-β)+iY0(v-y-β))dy=K0(vβ)I0(ρβ)+12∫-∞-βeyξyJ0(ρ-y-β)J0(v-y-β)dy,
because the imaginary part of log(z) with z∈𝒟β and Re(z)<β and Im(z)>0 has limit πi if z approach to the real axis. In definitive, after the change -y=x,(4.21)F1(ρ,ξ,v)=L-1[f1]=K0(vβ)I0(ρβ)-12∫β∞e-xξxJ0(ρx-β)J0(vx-β)dx.


### 4.2. Computation of𝔏-1[f2](ρ,ξ,v)

To find(4.22)∫1ρL-1[f2](ρ,ξ,v)dv,
we have v≤ρ, and hence, as a consequence of the previous result, we obtain(4.23)F2(ρ,ξ,v):=L-1[f2](ρ,t,v):=L-1[I0(vs+β)K0(ρs+β)s](ρ,ξ,v)=F1(v,ξ,ρ).


### 4.3. Computation of𝔏-1[f3](ρ,ξ,v)

Now we compute(4.24)F3(ρ,ξ):=L-1[f3](ρ,ξ):=L-1[K0(ρs+β)sK0(s+β)](ρ,ξ).
By (3.8), if s∈𝒟β,Re(s)≤γ0, and |s+β| is large enough in order that |O(1/ρs+β)|≤1/2 and |O(1/s+β)|≤1/2, we have(4.25)|K0(ρs+β)sK0(s+β)|≤1ρ1|s||e-(ρ-1)s+β||1+O(1/ρs+β)||1+O(1/s+β)|≤3ρ|s|e-(ρ-1)|s+β|cos⁡(Arg(s+β)/2)≤3ρ|s|,
(because cos(Arg(s+β)/2)≥0) and so Bromwich’s inversion formula cannot be used directly to find F3(ρ,ξ). To circumvent this complication we need to proceed in the way(4.26)F3(ρ,ξ)=L-1[K0(ρs+β)sK0(s+β)]=ddξL-1[K0(ρs+β)s2K0(s+β)],
finding this last inverse with means of Bromwich’s contour of Figure 2 since(4.27)|K0(ρs+β)s2K0(s+β)|≤3ρ|s|2,
if |s+β| is large enough.First, we remark that it follows easily from  (3.3) that(4.28)lim⁡s→-βs+βs2K0(ρs+β)K0(s+β)=0.
On the other hand, the function esξf3(s)/s has a pole of order 2 in s=0.  As the function
(4.29)esξK0(ρs+β)K0(s+β)
has continuous partial derivatives of any order in some neighbourhood of s=0 and every fixed t, we can write(4.30)ddξRes[esξK0(ρs+β)s2K0(s+β)]s=0=ddξ(lim⁡s→0ddss2esξK0(ρs+β)s2K0(s+β))=lim⁡s→0∂∂s∂∂ξesξK0(ρs+β)K0(s+β)=K0(ρβ)K0(β).Finally, sinceR(z) is an even function, by (3.3), (3.5), and (3.4), we have(4.31)lim⁡r→∞,ɛ→012πi(-∫L1esξK0(ρs+β)s2K0(s+β)ds-∫L2esξK0(ρs+β)s2K0(s+β)ds)=-12πi∫-∞-βesξs2-I0(ρ-s-βi)(log⁡(ρ-s-β/2)+(π/2)i)+Z(ρ-s-βi)-I0(-s-βi)(log⁡(-s-β/2)+(π/2)i)+R(-s-βi)ds+12πi∫-∞-βesξs2-I0(ρ-s-βi)(log⁡(ρ-s-β/2)-(π/2)i)+Z(ρ-s-βi)-I0(-s-βi)(log⁡(-s-β/2)-(π/2)i)+R(-s-βi)ds=-12πi∫-∞-βesξs2-J0(ρ-s-β)i-Y0(ρ-s-β)-J0(-s-β)i-Y0(-s-β)ds+12πi∫-∞-βesξs2J0(ρ-s-β)i-Y0(ρ-s-β)J0(-s-β)i-Y0(-s-β)ds=-1π∫-∞-βesξs2J0(ρ-s-β)Y0(-s-β)-Y0(ρ-s-β)J0(-s-β)J02(-s-β)+Y02(-s-β)ds.As a consequence, puttings=-x, we obtain from (4.26) that (4.32)F3(ρ,ξ):=L-1[K0(ρs+β)sK0(s+β)](ρ,ξ)=K0(ρβ)K0(β)+1π∫β∞e-xξxJ0(ρx-β)Y0(x-β)-Y0(ρx-β)J0(x-β)J02(x-β)+Y02(x-β)dx.


### 4.4. Computation of𝔏-1[f4](ρ,ξ,v)

To compute(4.33)∫1∞L-1[f4](ρ,ξ,v)vdv,
we have v≥1 and so, in order to find(4.34)L-1[f4](ρ,ξ,v):=L-1[K0(ρs+β)I0(s+β)K0(vs+β)sK0(s+β)](ρ,ξ,v),
using the convolution theorem and our previous results, we can write(4.35)F4(ρ,ξ,v):=L-1[f4](ρ,ξ,v)=L-1[K0(ρs+β)sK0(s+β)](ρ,ξ,v)*L-1[I0(s+β)K0(vs+β)](ρ,ξ,v)=F3(ρ,ξ)*ddξL-1[I0(s+β)K0(vs+β)s](ρ,ξ,v)=F3(ρ,ξ)*ddξF1(1,ξ,v)=12F3(ρ,ξ)*∫β∞e-xξJ0(x-β)J0(vx-β)dx=12∫0ξF3(ρ,w)(∫β∞e-x(ξ-w)J0(x-β)J0(vx-β)dx)dw.


### 4.5. Complete Solution

Collecting our previous results, by the second translation theorem, Fubini’s theorem, and some natural computations, we obtain finally(4.36)V(ρ,ξ)=∫ρ∞1v(K0(vβ)I0(ρβ)-12∫β∞e-xξxJ0(ρx-β)J0(vx-β)dx)dv+∫1ρ1v(K0(ρβ)I0(vβ)-12∫β∞e-xξxJ0(vx-β)J0(ρx-β)dx)dv-BF3(ρ,ξ)-12∫1∞1v(∫0ξF3(ρ,w)(∫β∞e-x(ξ-w)J0(x-β)J0(vx-β)dx)dw)dv,
that can be rewritten in the following way which is more suitable for numerical computations:(4.37)V(ρ,ξ)=∫1∞1vK0(vβ)I0(ρβ)dv-∫1ρ1vK0(vβ)I0(ρβ)dv+∫1ρ1vK0(ρβ)I0(vβ)dv-BF3(ρ,ξ)-12∫1∞1v(∫β∞e-xξxJ0(ρx-β)J0(vx-β)dx)dv-12∫1∞1v(∫0ξF3(ρ,w)(∫β∞e-x(ξ-w)J0(x-β)J0(vx-β)dx)dw)dv.It is interesting to check that, as expected, actually we haveV(1,ξ)=-B at every time ξ>0. In fact, as F3(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(4.38)∫0ξe-x(ξ-w)dw=1x-e-xξx.
Finally, after the variable change z=x-β and the application of formula (5) in [7], Section  13.53, we obtain(4.39)∫β∞J0(x-β)J0(vx-β)xdx=2∫0∞J0(z)J0(vz)z2+βzdz=2I0(β)K0(vβ),
and this gives us V(1,ξ)=-B.


## 4.1. Computation of𝔏-1[f1](ρ,ξ,v)

Asv≥ρ in the integral(4.14)∫ρ∞L-1[f1](ρ,ξ,v)vdv,
we begin computing, for every v≥ρ,(4.15)F1(ρ,ξ,v):=L-1[f1](ρ,ξ,v):=L-1[I0(ρs+β)K0(vs+β)s](ρ,ξ,v).Asv≥ρ, it follows easily from the asymptotic expansions (3.8) and (3.9) that if s∈𝒟β and |s+β| is large enough in order that |O(1/ρs+β)|≤1/2,|O(1/vs+β)|≤1/2, and |β/s|≤1/2 and, moreover, Re(s)≤γ0 and Im(s)>0, we have the estimation(4.16)|I0(ρs+β)K0(vs+β)s|≤9(|e-(v-ρ)s+β|+|e-(v+ρ)s+β|)8vρ|s||s+β|≤9(|e-(v-ρ)|s+β|cos⁡(Arg(s+β)/2)|+|e-(v+ρ)|s+β|cos⁡(Arg(s+β)/2)|)8vρ|s||s+β|≤92vρ|s|3/2,
because cos(Arg(s+β)/2)≥0 and |1+β/s|≤1-1/2=1/2. A similar result holds if Im(s)<0 by (3.10). Then we can to apply Bromwich’s formula to find 𝔏-1[f1](ρ,ξ,v).The residue in the poles=0 is(4.17)R(s=0)=lim⁡s→0sesξI0(ρs+β)K0(vs+β)s=K0(vβ)I0(ρβ).
On the other hand, by (3.1) we obtain(4.18)lim⁡s→-β,s∈Dβ(s+β)esξI0(ρs+β)K0(vs+β)s=0.
Hence,(4.19)L-1[f1]=K0(vβ)I0(ρβ)-12πi∫-∞-βesξsK0(v-s-βi)I0(ρ-s-βi)ds+12πi∫-∞-βesξsK0(-v-s-βi)I0(-ρ-s-βi)ds.
Having in mind (3.7) and (3.3), we obtain(4.20)L-1[f1]=K0(vβ)I0(ρβ)-12πi∫-∞-βπi2eyξyJ0(ρ-y-β)(J0(v-y-β)+iY0(-v-y-β))dy+12πi∫-∞-βπi2eyξyJ0(ρ-y-β)(J0(v-y-β)+iY0(v-y-β))dy=K0(vβ)I0(ρβ)+12∫-∞-βeyξyJ0(ρ-y-β)J0(v-y-β)dy,
because the imaginary part of log(z) with z∈𝒟β and Re(z)<β and Im(z)>0 has limit πi if z approach to the real axis. In definitive, after the change -y=x,(4.21)F1(ρ,ξ,v)=L-1[f1]=K0(vβ)I0(ρβ)-12∫β∞e-xξxJ0(ρx-β)J0(vx-β)dx.


## 4.2. Computation of𝔏-1[f2](ρ,ξ,v)

To find(4.22)∫1ρL-1[f2](ρ,ξ,v)dv,
we have v≤ρ, and hence, as a consequence of the previous result, we obtain(4.23)F2(ρ,ξ,v):=L-1[f2](ρ,t,v):=L-1[I0(vs+β)K0(ρs+β)s](ρ,ξ,v)=F1(v,ξ,ρ).


## 4.3. Computation of𝔏-1[f3](ρ,ξ,v)

Now we compute(4.24)F3(ρ,ξ):=L-1[f3](ρ,ξ):=L-1[K0(ρs+β)sK0(s+β)](ρ,ξ).
By (3.8), if s∈𝒟β,Re(s)≤γ0, and |s+β| is large enough in order that |O(1/ρs+β)|≤1/2 and |O(1/s+β)|≤1/2, we have(4.25)|K0(ρs+β)sK0(s+β)|≤1ρ1|s||e-(ρ-1)s+β||1+O(1/ρs+β)||1+O(1/s+β)|≤3ρ|s|e-(ρ-1)|s+β|cos⁡(Arg(s+β)/2)≤3ρ|s|,
(because cos(Arg(s+β)/2)≥0) and so Bromwich’s inversion formula cannot be used directly to find F3(ρ,ξ). To circumvent this complication we need to proceed in the way(4.26)F3(ρ,ξ)=L-1[K0(ρs+β)sK0(s+β)]=ddξL-1[K0(ρs+β)s2K0(s+β)],
finding this last inverse with means of Bromwich’s contour of Figure 2 since(4.27)|K0(ρs+β)s2K0(s+β)|≤3ρ|s|2,
if |s+β| is large enough.First, we remark that it follows easily from  (3.3) that(4.28)lim⁡s→-βs+βs2K0(ρs+β)K0(s+β)=0.
On the other hand, the function esξf3(s)/s has a pole of order 2 in s=0.  As the function
(4.29)esξK0(ρs+β)K0(s+β)
has continuous partial derivatives of any order in some neighbourhood of s=0 and every fixed t, we can write(4.30)ddξRes[esξK0(ρs+β)s2K0(s+β)]s=0=ddξ(lim⁡s→0ddss2esξK0(ρs+β)s2K0(s+β))=lim⁡s→0∂∂s∂∂ξesξK0(ρs+β)K0(s+β)=K0(ρβ)K0(β).Finally, sinceR(z) is an even function, by (3.3), (3.5), and (3.4), we have(4.31)lim⁡r→∞,ɛ→012πi(-∫L1esξK0(ρs+β)s2K0(s+β)ds-∫L2esξK0(ρs+β)s2K0(s+β)ds)=-12πi∫-∞-βesξs2-I0(ρ-s-βi)(log⁡(ρ-s-β/2)+(π/2)i)+Z(ρ-s-βi)-I0(-s-βi)(log⁡(-s-β/2)+(π/2)i)+R(-s-βi)ds+12πi∫-∞-βesξs2-I0(ρ-s-βi)(log⁡(ρ-s-β/2)-(π/2)i)+Z(ρ-s-βi)-I0(-s-βi)(log⁡(-s-β/2)-(π/2)i)+R(-s-βi)ds=-12πi∫-∞-βesξs2-J0(ρ-s-β)i-Y0(ρ-s-β)-J0(-s-β)i-Y0(-s-β)ds+12πi∫-∞-βesξs2J0(ρ-s-β)i-Y0(ρ-s-β)J0(-s-β)i-Y0(-s-β)ds=-1π∫-∞-βesξs2J0(ρ-s-β)Y0(-s-β)-Y0(ρ-s-β)J0(-s-β)J02(-s-β)+Y02(-s-β)ds.As a consequence, puttings=-x, we obtain from (4.26) that (4.32)F3(ρ,ξ):=L-1[K0(ρs+β)sK0(s+β)](ρ,ξ)=K0(ρβ)K0(β)+1π∫β∞e-xξxJ0(ρx-β)Y0(x-β)-Y0(ρx-β)J0(x-β)J02(x-β)+Y02(x-β)dx.


## 4.4. Computation of𝔏-1[f4](ρ,ξ,v)

To compute(4.33)∫1∞L-1[f4](ρ,ξ,v)vdv,
we have v≥1 and so, in order to find(4.34)L-1[f4](ρ,ξ,v):=L-1[K0(ρs+β)I0(s+β)K0(vs+β)sK0(s+β)](ρ,ξ,v),
using the convolution theorem and our previous results, we can write(4.35)F4(ρ,ξ,v):=L-1[f4](ρ,ξ,v)=L-1[K0(ρs+β)sK0(s+β)](ρ,ξ,v)*L-1[I0(s+β)K0(vs+β)](ρ,ξ,v)=F3(ρ,ξ)*ddξL-1[I0(s+β)K0(vs+β)s](ρ,ξ,v)=F3(ρ,ξ)*ddξF1(1,ξ,v)=12F3(ρ,ξ)*∫β∞e-xξJ0(x-β)J0(vx-β)dx=12∫0ξF3(ρ,w)(∫β∞e-x(ξ-w)J0(x-β)J0(vx-β)dx)dw.


## 4.5. Complete Solution

Collecting our previous results, by the second translation theorem, Fubini’s theorem, and some natural computations, we obtain finally(4.36)V(ρ,ξ)=∫ρ∞1v(K0(vβ)I0(ρβ)-12∫β∞e-xξxJ0(ρx-β)J0(vx-β)dx)dv+∫1ρ1v(K0(ρβ)I0(vβ)-12∫β∞e-xξxJ0(vx-β)J0(ρx-β)dx)dv-BF3(ρ,ξ)-12∫1∞1v(∫0ξF3(ρ,w)(∫β∞e-x(ξ-w)J0(x-β)J0(vx-β)dx)dw)dv,
that can be rewritten in the following way which is more suitable for numerical computations:(4.37)V(ρ,ξ)=∫1∞1vK0(vβ)I0(ρβ)dv-∫1ρ1vK0(vβ)I0(ρβ)dv+∫1ρ1vK0(ρβ)I0(vβ)dv-BF3(ρ,ξ)-12∫1∞1v(∫β∞e-xξxJ0(ρx-β)J0(vx-β)dx)dv-12∫1∞1v(∫0ξF3(ρ,w)(∫β∞e-x(ξ-w)J0(x-β)J0(vx-β)dx)dw)dv.It is interesting to check that, as expected, actually we haveV(1,ξ)=-B at every time ξ>0. In fact, as F3(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(4.38)∫0ξe-x(ξ-w)dw=1x-e-xξx.
Finally, after the variable change z=x-β and the application of formula (5) in [7], Section  13.53, we obtain(4.39)∫β∞J0(x-β)J0(vx-β)xdx=2∫0∞J0(z)J0(vz)z2+βzdz=2I0(β)K0(vβ),
and this gives us V(1,ξ)=-B.


## 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 timet.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 (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/m3, specific heat (cb) of 4148 J/kg·K, and a perfusion rate (ωb) of 6.410-3s-1 [9]. Blood temperature, and hence initial tissue temperature (Tb), was 37°C. For all the simulations, we used a current density (j0) 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(TC). 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.Figure 3
Temperatures at 60 s, 180 s, and 360 s and limit temperature (thick line) with a temperature of 5°C in the electrode surface and a current intensity 5 mA/mm2.Figure 4
Temperatures at 60 s, 180 s, and 360 s and limit temperature (thick line) with a temperature of 15°C in the electrode surface and a current intensity 5 mA/mm2.Figure 5
Temperatures at 60 s,180 s, and 360 s and limit temperature (thick line) with a temperature of 25°C in the electrode surface and a current intensity 5 mA/mm2.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.



---

*Source: 834807-2012-01-11.xml*