# An Analytical Solution for Radiofrequency Ablation with a Cooled Cylindrical Electrode

Authors: Ricardo Romero-Méndez; Enrique Berjano
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

We present an analytical solution to the electrothermal mathematical model of radiofrequency ablation of biological tissue using a cooled cylindrical electrode. The solution presented here makes use of the method of separation of variables to solve the problem. Green’s functions are used for the handling of nonhomogeneous terms, such as effect of electrical currents circulation and the nonhomogeneous boundary condition due to cooling at the electrode surface. The transcendental equation for determination of eigenvalues of this problem is solved using Newton’s method, and the integrals that appear in the solution of the problem are obtained by Simpson’s rule. The solution obtained here has the possibility of handling different functional dependencies of the source term and nonhomogeneous boundary condition. The solution provides a tool to understand the physics of the problem, as it shows how the solution depends on different parameters, to provide mathematical tools for the design of surgical procedures and to validate other modeling techniques, such as the numerical methods that are frequently used to solve the problem.


---

## Body

## 1. Introduction

Radiofrequency (RF) ablation is a surgical technique considered as marginally invasive, used, for instance, to destroy malignant tissue in bland organs, such as liver, kidney, and lungs. The objective of this treatment is to heat the tumor to a temperature above a threshold where the tissue is damaged. Initially, this procedure was limited to small volumes of tissue damage. Rossi et al. [1] and Goldberg [2] reported that the diameter of tissue damage attained by a single electrode was of only 1.5 cm, reason for which RF ablation was only considered viable for treatment of small tumors. A problem of RF ablation is tissue charring around the electrode, which produces a phenomenon called roll-off, in which the tissue surrounding the electrode increases its electrical resistance and with it the circuit impedance, with which current circulation stops and no further damage to tissue is produced [3]. Several alternatives have been used to avoid early appearance of roll-off and to increase the volume of tissue damaged: internally cooled electrodes, expandable electrodes, hybrid electrodes, and so forth [4]. The first alternative is to use internal cooling of the active electrode. Internally cooled electrode (known as cool-tip electrode) is one of the most frequently used alternatives for hepatic tissue damage [5, 6]; it consists of a needle shaped electrode with an active tip, inside which there is a circulation of chilled water, which cools down the tissue that surrounds the electrode and avoids carbonization of the tissue [7].Much of the studies that have been performed in relation to RF current ablation have relied on experimental techniques and numerical methods. Very few analytical solutions have been reported in relation to this phenomenon. Haemmerich et al. [8] presented analytical results of the modeling of RF ablation for steady state conditions, where blood perfusion and metabolic heat generation were neglected. López Molina et al. [9] presented a solution to the transient problem for an infinite cylindrical geometry that was solved by Laplace transform. Here we present an alternative solution to the time dependent RF ablation heating process from a cooled cylindrical RF electrode in a finite domain using Green’s function. This solution allows including any kind of source term (time dependent and with any spatial dependence) and boundary conditions that may be time dependent. The techniques presented here provide a closed form solution of this phenomenon, with which it is possible to understand better the influence of the different factors influencing tissue heating, and it is useful also as a tool for the design of surgical protocols and for validating numerical solutions obtained for the purpose of modeling RF ablation.


## 2. Methods

### 2.1. Problem Description and Mathematical Model

Figure1 shows how the geometry of the problem was defined. RF ablation is performed by placing the active electrode in the middle of the tumor. The problem is illustrated by Figure 1(a), where the tissue region is modeled as a cylindrical volume, with the active electrode placed inside the tissue and the dispersive electrode being placed on the external surface. The tissue consists of a region of cancerous tissue (tumor) of spherical geometry, surrounded by healthy tissue, represented by the cylindrical geometry. The active electrode is illustrated in Figure 1(a) as the part of the RF applicator placed inside the sphere that represents the cancerous tissue. The axial mid-plane of the conductive part of the active electrode, represented by the circle placed at the mid-height of the domain, is modeled as a symmetry plane, where the electric and thermal problems behave as axisymmetric one dimensional problems, depending only on the radial coordinate. For simplicity, the region to be analyzed is this plane: the annulus of inner radius ri and outer radius ro. In some studied cases, here we also consider a two-compartment model, with cancerous tissue located in the zone ri<r<rδ, while the healthy tissue is in the region rδ<r<ro, represented by Figure 1(b). Inside this annular section of tissue is an active RF electrode of radius ri and the dispersive electrode is considered to be at the periphery of the tissue.Figure 1
(a) Physical situation modeled: an RF applicator is inserted in the tissue and the active electrode is placed in the center of a spherical tumor. (b) Analytical model has one dimension (r) and a domain from ri (active electrode radius) to ro (where the dispersive electrode is assumed to be placed). The active electrode is in the middle of a tumor of radius rδ, and consequently the healthy tissue occupies a region between rδ and ro.
(a)
(b)Consider then the case of a fragment of tissue in the shape of a concentric annulus with inner radiusri and outer radius ro. The tissue is originally at the basal temperature Tbasal=Tb+q˙m/(ρbωcb), where q˙m is the metabolic heat generation of the tissue [W m-3], Tb=37∘C is the blood arterial temperature, and ρb and cb are the blood density [kg m-3] and specific heat [J kg-1K-1], while ω is the blood perfusion rate [s-1]. The inner radius of this section is an RF active electrode held at a low temperature Tc while the outer radius is at the unaltered temperature Tbasal. The inner radius also has an applied voltage Vi, while the outer radius is considered held at a zero voltage (i.e., mimicking the dispersive electrode). Because of voltage potential, electrical current circulates through the tissue, and the temperature gradually increases within the tissue. This phenomenon is relevant because exposing the tissue to high temperature produces tissue damage. The procedure will be considered done when the tissue reaches a temperature of 100∘C at some point, because at that temperature tissue charring around the electrode occurs, which produces a phenomenon called roll-off, in which the tissue surrounding the electrode increases suddenly its electrical resistivity and with it the circuit impedance, with which current circulation stops and no further damage to tissue is produced. A constant electrical conductivity σ [S m-1] is considered for the range of temperatures Tc≤T≤100∘C. The domain of the solution is 0.75mm<r<10cm. Table 1 presents the set of properties and geometrical parameters used throughout, unless specified otherwise when required. For most cases Tc was taken as 5∘C, unless specified otherwise.Table 1
Properties of liver tumor tissue used in this study [10, 11].
Property
Nomenclature
Magnitude
Units
Tissue density
ρ
1000
kgm-3
Tissue thermal conductivity
k
0.512
Wm-1 K-1
Tissue specific heat
c
4200
Jkg-1 K-1
Blood density
ρ
b
1000
Kgm-3
Blood perfusion
ω
0.0015
s
-
1
Blood specific heat
c
b
3600
Jkg-1 K-1
Tissue electrical conductivity
σ
0.188
Sm-1
Metabolic heat generation
q
˙
m
700
Wm-3The governing equation for this problem is Pennes’ Bioheat equation [12]:(1)ρc∂T∂t=∇·k∇T-ρbωcbT-Tb+q˙m+q˙RF,where k, ρ, and c are the tissue thermal conductivity [W m-1K-1], density [kg m-3], and specific heat [J kg-1K-1], q˙m is the metabolic heat generation per unit volume within the tissue, and ρbωcb(T-Tb) is the blood perfusion heat transport term. q˙RF is the volumetric heat generation due to circulation of electrical currents in the RF range [W m-3].The calculation ofq˙RF requires solution of the electrical problem. This is stated in terms of Laplace’s equation:(2)∇·σ∇V=0,where σ is the electrical conductivity of tissue.The heat generation per unit volume from RF currents circulation,q˙RF, that appears in (1) is determined from(3)q˙RF=E·J.The electrical potential E [W m-1] and electrical current density J [A m-2] can be determined from(4)E=-∇V,J=σEso that(5)q˙RF=σ∇V2.Considering that in this problem the voltage is only a function of the radius andσ is constant, (2) can be stated as(6)ddrrdVdr=0subject to boundary conditions(7)Vri=Vi,Vro=0.With these equation and set of boundary conditions, the solution can be stated as(8)Vr=Viln⁡ro/riln⁡r-ln⁡ro.The electrical potentialE and electrical current density J can be determined by introducing (8) into (4):(9)E=-1rViln⁡ro/rie^r,J=-σrViln⁡ro/rie^r.And the heat generation per unit volume is determined from(10)q˙RF=E·J=σVi2r2ln⁡ro/ri2=κr2,where κ=σVi2/[ln⁡(ro/ri)]2.


### 2.2. Analytical Solution

Considering a temperature of tissueT=T(r,t), (1) can be written as(11)1r∂∂rkr∂T∂r-ρbωcbT-Tb+q˙m+q˙RF=ρc∂T∂t,subject to(12)Tri,t=Tc,Tro,t=Tbasal,Tr,0=Tbasal.Using a change of variableθ=T-Tbasal, (11) can be written as(13)1r∂∂rkr∂θ∂r-ρbωcbθ+q˙RF=ρc∂θ∂t,subject to(14)θri,t=θc,θro,t=0,θr,0=0.This problem is nonhomogeneous in the equation itself, because of termq˙RF and also because of the nonhomogeneous boundary condition at r=ri. The problem is solved here using Green’s functions [13]. The solution in terms of Green’s functions can be stated as (15)θr,t=∫ϑGr,t∣r∗,ττ=0Frdv∗+αk∫τ=0tdτ∫ϑGr,t∣r∗,τgr∗,τdv∗-α∫τ=0tr∗∂G∂r∗r=rifτdτ,where G(r,t∣r∗,τ) is Green’s function. In this case F(r)=0,  g(r∗,τ)=q˙RF=κ/r∗2, and  f(τ)=θc.  α=k/(ρc) is the tissue thermal diffusivity [m2 s-1].To complete the solution, suitable Green’s function must be obtained. The solution of an associated homogeneous equation and boundary conditions provides a way to identify Green’s function for this problem. We take the auxilliary equation(16)1r∂∂rkr∂ψ∂r-ρbωcbψ=ρc∂ψ∂t,subject to(17)ψri,t=0,(18)ψro=0,(19)ψr,0=Fir,where Fi(r) is a generic initial condition function.The method of separation of variables can be used to solve (16). The method is based on the concept of expressing a multivariable dependent function as the product of several functions that are each dependent on one of the different independent variables. In this particular case, by using the new variable(20)ψr,t=RrΓtthe equation can be written as(21)Γd2Rdlr2+ΓrdRdr-ρbωcbkΓR=RαdΓdt=-λ2,where λ is a constant that later will be related to the eigenvalues of the problem.By definingm2=ρbωcb/k and dividing every term by ΓR two separated equations are produced:(22)d2Rdr2+1rdRdr+β2R=0subject to(23)Rri=0,Rro=0,where β2=λ2-m2 has been introduced. The solution to this equation is(24)Rr=J0βnrY0βnro-J0βnroY0βnr, where the eigenvalues βn are the positive roots of equation(25)J0βnriY0βnro-J0βnroY0βnri=0.Equation (22) is called Bessel’s differential equation of order zero, whose solutions J0(βnr) and Y0(βnr) are called Bessel’s functions of order zero of the first and second kind, respectively.The equation forΓ is(26)dΓdt+αλ2=0.The general solution to this equation is(27)Γt=Cnexp⁡-αλ2t.Writing againψ=R(r)Γ(t) and considering summation over all possible eigenvalues βn, we get(28)ψr,t=∑j=1∞Cnexp⁡-αβn2+m2tJ0βnrY0βnro-J0βnroY0βnr.Applying the initial condition expressed in (19), ψ(r,t) can be written as (29)ψi=π22∑n=1∞βn2J02βnriJ02βnri-J02βnro∫rirorJ0βnrY0βnro-J0βnroY0βnrFir·exp⁡-αβn2+m2tJ0βnrY0βnro-J0βnroY0βnrdr.The solution to the problem in terms of Green’s function can be stated as(30)ψr,t=∫riror∗Gr,t∣r∗,ττ=0Fir∗dr∗.A comparison of terms, after replacingt by t-τ, provide Green’s function for this problem (31)Gr,t∣r∗,τ=π22∑n=1∞βn2J02βnriJ02βnri-J02βnroJ0βnr∗Y0βnro-J0βnroY0βnr∗·exp⁡-αβn2+m2t-τJ0βnrY0βnro-J0βnroY0βnr.Once Green’s function for the problem is obtained, the solution to (13) and (14) can be stated as (32)θr,t=π22k∑n=1∞1-e-αβn2+m2tβn2+m2βn2J02βnriJ02βnri-J02βnroRoβnr∫riror∗Roβnr∗q˙RFr∗dr∗-α∫τ=0tr∗∂G∂r∗r=riθcdτ,where(33)Roβnr=J0βnrY0βnro-J0βnroY0βnr.The first term of (32) accounts for the source term q˙RF, while the last term accounts for the nonhomogeneous boundary condition at r=ri. This last term is evaluated as(34)αr∗∂G∂r∗r=ri=απ22∑n=1∞e-αβn2+m2t-τβn2J02βnririJ02βnri-J02βnroRoβnrRo′βnri,where(35)Ro′βnri=βn-J1βnriY0βnro+J0βnroY1βnri.Since, in this casef(t)=θc=constant, the integral to account for the nonhomogeneous boundary condition at r=ri can be evaluated easily, then (36)-α∫τ=0tr∗∂G∂r∗riθcdτ=∑n=1∞θcπ2ri2βn2+m2βn2J02βnriJ02βnri-J02βnroe-αβn2+m2t-1·RoβnrRo′βnri.The final form of the solution is(37)Tr,t=π22k∑n=1∞1-e-αβn2+m2tβn2+m2βn2J02βnriJ02βnri-J02βnroRoβnr∫riror∗Roβnr∗q˙RFr∗dr∗+∑n=1∞θcπ2ri2βn2+m2βn2J02βnriJ02βnri-J02βnroe-αβn2+m2t-1RoβnrRo′βnri+Tb+q˙mρbωcb.Equation (25) for determination of eigenvalues of the solution is a transcendental equation. In order to determine the roots a numerical procedure was used. The procedure consisted of identifying values of βn in the neighborhood where the term(38)Reiβn=J0βnriY0βnro-J0βnroY0βnriswitches from positive to negative or vice versa. The search of regions where roots are located consisted of calculating Rei(β), starting with β=0 and with increments Δβ=0.1, to find zones of two consecutive β values where Rei switches sign. Once these regions are identified, Newton’s method was run at each one of those regions, trying to find the values of βn that produce a value of (38) closer to zero. The iterative procedure uses the following equation:(39)βni+1=βni-Rei′βniReiβni. The iterations proceed as long as Rei(βni+1)>ϵ, where ϵ is the convergence criteria. When the convergence criteria is met, βni+1 is taken as βn. In this investigation a value ϵ=10-14 was used. Table 2 presents a list of the first 25 eigenvalues obtained by this procedure for ri=0.75mm and ro=10cm. It is important to see that eigenvalues depend only on geometric conditions of the problem. 625 eigenvalues were used in obtaining numerical values of (37).Table 2
First 20 eigenvalues of transcendental equation (25) for this example.
Eigenvalue number
β
n
1
27.7412491590629
2
59.7215157304870
3
91.6538083942138
4
123.5335630443243
5
155.3754428957278
6
187.1895966699011
7
218.9827041355604
8
250.7593126035128
9
282.5226419505083
10
314.2750529359010
11
346.0183280786956
12
377.7538469185017
13
409.4826995869618
14
441.2057629309977
15
472.9237530699340
16
504.6372626347768
17
536.3467877607103
18
568.0527480482156
19
599.7555015887882
20
631.4553564558279The integrals that appear in the solution are evaluated using Simpson’s rule, using a fine mesh of 397 divisions of the regionri<r<ro. The integral obtained with 397 divisions was compared to a calculation of the integral with 794 divisions. The comparison of the integral of (37) for the leading eigenvalue (beta1) gave a difference of 0.0058% between the cases of calculation of the integral with 397 divisions and 794 divisions. It is important to mention that the integrals that appear in (37) are evaluated once for every eigenvalue βn and do not change with r or t.


## 2.1. Problem Description and Mathematical Model

Figure1 shows how the geometry of the problem was defined. RF ablation is performed by placing the active electrode in the middle of the tumor. The problem is illustrated by Figure 1(a), where the tissue region is modeled as a cylindrical volume, with the active electrode placed inside the tissue and the dispersive electrode being placed on the external surface. The tissue consists of a region of cancerous tissue (tumor) of spherical geometry, surrounded by healthy tissue, represented by the cylindrical geometry. The active electrode is illustrated in Figure 1(a) as the part of the RF applicator placed inside the sphere that represents the cancerous tissue. The axial mid-plane of the conductive part of the active electrode, represented by the circle placed at the mid-height of the domain, is modeled as a symmetry plane, where the electric and thermal problems behave as axisymmetric one dimensional problems, depending only on the radial coordinate. For simplicity, the region to be analyzed is this plane: the annulus of inner radius ri and outer radius ro. In some studied cases, here we also consider a two-compartment model, with cancerous tissue located in the zone ri<r<rδ, while the healthy tissue is in the region rδ<r<ro, represented by Figure 1(b). Inside this annular section of tissue is an active RF electrode of radius ri and the dispersive electrode is considered to be at the periphery of the tissue.Figure 1
(a) Physical situation modeled: an RF applicator is inserted in the tissue and the active electrode is placed in the center of a spherical tumor. (b) Analytical model has one dimension (r) and a domain from ri (active electrode radius) to ro (where the dispersive electrode is assumed to be placed). The active electrode is in the middle of a tumor of radius rδ, and consequently the healthy tissue occupies a region between rδ and ro.
(a)
(b)Consider then the case of a fragment of tissue in the shape of a concentric annulus with inner radiusri and outer radius ro. The tissue is originally at the basal temperature Tbasal=Tb+q˙m/(ρbωcb), where q˙m is the metabolic heat generation of the tissue [W m-3], Tb=37∘C is the blood arterial temperature, and ρb and cb are the blood density [kg m-3] and specific heat [J kg-1K-1], while ω is the blood perfusion rate [s-1]. The inner radius of this section is an RF active electrode held at a low temperature Tc while the outer radius is at the unaltered temperature Tbasal. The inner radius also has an applied voltage Vi, while the outer radius is considered held at a zero voltage (i.e., mimicking the dispersive electrode). Because of voltage potential, electrical current circulates through the tissue, and the temperature gradually increases within the tissue. This phenomenon is relevant because exposing the tissue to high temperature produces tissue damage. The procedure will be considered done when the tissue reaches a temperature of 100∘C at some point, because at that temperature tissue charring around the electrode occurs, which produces a phenomenon called roll-off, in which the tissue surrounding the electrode increases suddenly its electrical resistivity and with it the circuit impedance, with which current circulation stops and no further damage to tissue is produced. A constant electrical conductivity σ [S m-1] is considered for the range of temperatures Tc≤T≤100∘C. The domain of the solution is 0.75mm<r<10cm. Table 1 presents the set of properties and geometrical parameters used throughout, unless specified otherwise when required. For most cases Tc was taken as 5∘C, unless specified otherwise.Table 1
Properties of liver tumor tissue used in this study [10, 11].
Property
Nomenclature
Magnitude
Units
Tissue density
ρ
1000
kgm-3
Tissue thermal conductivity
k
0.512
Wm-1 K-1
Tissue specific heat
c
4200
Jkg-1 K-1
Blood density
ρ
b
1000
Kgm-3
Blood perfusion
ω
0.0015
s
-
1
Blood specific heat
c
b
3600
Jkg-1 K-1
Tissue electrical conductivity
σ
0.188
Sm-1
Metabolic heat generation
q
˙
m
700
Wm-3The governing equation for this problem is Pennes’ Bioheat equation [12]:(1)ρc∂T∂t=∇·k∇T-ρbωcbT-Tb+q˙m+q˙RF,where k, ρ, and c are the tissue thermal conductivity [W m-1K-1], density [kg m-3], and specific heat [J kg-1K-1], q˙m is the metabolic heat generation per unit volume within the tissue, and ρbωcb(T-Tb) is the blood perfusion heat transport term. q˙RF is the volumetric heat generation due to circulation of electrical currents in the RF range [W m-3].The calculation ofq˙RF requires solution of the electrical problem. This is stated in terms of Laplace’s equation:(2)∇·σ∇V=0,where σ is the electrical conductivity of tissue.The heat generation per unit volume from RF currents circulation,q˙RF, that appears in (1) is determined from(3)q˙RF=E·J.The electrical potential E [W m-1] and electrical current density J [A m-2] can be determined from(4)E=-∇V,J=σEso that(5)q˙RF=σ∇V2.Considering that in this problem the voltage is only a function of the radius andσ is constant, (2) can be stated as(6)ddrrdVdr=0subject to boundary conditions(7)Vri=Vi,Vro=0.With these equation and set of boundary conditions, the solution can be stated as(8)Vr=Viln⁡ro/riln⁡r-ln⁡ro.The electrical potentialE and electrical current density J can be determined by introducing (8) into (4):(9)E=-1rViln⁡ro/rie^r,J=-σrViln⁡ro/rie^r.And the heat generation per unit volume is determined from(10)q˙RF=E·J=σVi2r2ln⁡ro/ri2=κr2,where κ=σVi2/[ln⁡(ro/ri)]2.


## 2.2. Analytical Solution

Considering a temperature of tissueT=T(r,t), (1) can be written as(11)1r∂∂rkr∂T∂r-ρbωcbT-Tb+q˙m+q˙RF=ρc∂T∂t,subject to(12)Tri,t=Tc,Tro,t=Tbasal,Tr,0=Tbasal.Using a change of variableθ=T-Tbasal, (11) can be written as(13)1r∂∂rkr∂θ∂r-ρbωcbθ+q˙RF=ρc∂θ∂t,subject to(14)θri,t=θc,θro,t=0,θr,0=0.This problem is nonhomogeneous in the equation itself, because of termq˙RF and also because of the nonhomogeneous boundary condition at r=ri. The problem is solved here using Green’s functions [13]. The solution in terms of Green’s functions can be stated as (15)θr,t=∫ϑGr,t∣r∗,ττ=0Frdv∗+αk∫τ=0tdτ∫ϑGr,t∣r∗,τgr∗,τdv∗-α∫τ=0tr∗∂G∂r∗r=rifτdτ,where G(r,t∣r∗,τ) is Green’s function. In this case F(r)=0,  g(r∗,τ)=q˙RF=κ/r∗2, and  f(τ)=θc.  α=k/(ρc) is the tissue thermal diffusivity [m2 s-1].To complete the solution, suitable Green’s function must be obtained. The solution of an associated homogeneous equation and boundary conditions provides a way to identify Green’s function for this problem. We take the auxilliary equation(16)1r∂∂rkr∂ψ∂r-ρbωcbψ=ρc∂ψ∂t,subject to(17)ψri,t=0,(18)ψro=0,(19)ψr,0=Fir,where Fi(r) is a generic initial condition function.The method of separation of variables can be used to solve (16). The method is based on the concept of expressing a multivariable dependent function as the product of several functions that are each dependent on one of the different independent variables. In this particular case, by using the new variable(20)ψr,t=RrΓtthe equation can be written as(21)Γd2Rdlr2+ΓrdRdr-ρbωcbkΓR=RαdΓdt=-λ2,where λ is a constant that later will be related to the eigenvalues of the problem.By definingm2=ρbωcb/k and dividing every term by ΓR two separated equations are produced:(22)d2Rdr2+1rdRdr+β2R=0subject to(23)Rri=0,Rro=0,where β2=λ2-m2 has been introduced. The solution to this equation is(24)Rr=J0βnrY0βnro-J0βnroY0βnr, where the eigenvalues βn are the positive roots of equation(25)J0βnriY0βnro-J0βnroY0βnri=0.Equation (22) is called Bessel’s differential equation of order zero, whose solutions J0(βnr) and Y0(βnr) are called Bessel’s functions of order zero of the first and second kind, respectively.The equation forΓ is(26)dΓdt+αλ2=0.The general solution to this equation is(27)Γt=Cnexp⁡-αλ2t.Writing againψ=R(r)Γ(t) and considering summation over all possible eigenvalues βn, we get(28)ψr,t=∑j=1∞Cnexp⁡-αβn2+m2tJ0βnrY0βnro-J0βnroY0βnr.Applying the initial condition expressed in (19), ψ(r,t) can be written as (29)ψi=π22∑n=1∞βn2J02βnriJ02βnri-J02βnro∫rirorJ0βnrY0βnro-J0βnroY0βnrFir·exp⁡-αβn2+m2tJ0βnrY0βnro-J0βnroY0βnrdr.The solution to the problem in terms of Green’s function can be stated as(30)ψr,t=∫riror∗Gr,t∣r∗,ττ=0Fir∗dr∗.A comparison of terms, after replacingt by t-τ, provide Green’s function for this problem (31)Gr,t∣r∗,τ=π22∑n=1∞βn2J02βnriJ02βnri-J02βnroJ0βnr∗Y0βnro-J0βnroY0βnr∗·exp⁡-αβn2+m2t-τJ0βnrY0βnro-J0βnroY0βnr.Once Green’s function for the problem is obtained, the solution to (13) and (14) can be stated as (32)θr,t=π22k∑n=1∞1-e-αβn2+m2tβn2+m2βn2J02βnriJ02βnri-J02βnroRoβnr∫riror∗Roβnr∗q˙RFr∗dr∗-α∫τ=0tr∗∂G∂r∗r=riθcdτ,where(33)Roβnr=J0βnrY0βnro-J0βnroY0βnr.The first term of (32) accounts for the source term q˙RF, while the last term accounts for the nonhomogeneous boundary condition at r=ri. This last term is evaluated as(34)αr∗∂G∂r∗r=ri=απ22∑n=1∞e-αβn2+m2t-τβn2J02βnririJ02βnri-J02βnroRoβnrRo′βnri,where(35)Ro′βnri=βn-J1βnriY0βnro+J0βnroY1βnri.Since, in this casef(t)=θc=constant, the integral to account for the nonhomogeneous boundary condition at r=ri can be evaluated easily, then (36)-α∫τ=0tr∗∂G∂r∗riθcdτ=∑n=1∞θcπ2ri2βn2+m2βn2J02βnriJ02βnri-J02βnroe-αβn2+m2t-1·RoβnrRo′βnri.The final form of the solution is(37)Tr,t=π22k∑n=1∞1-e-αβn2+m2tβn2+m2βn2J02βnriJ02βnri-J02βnroRoβnr∫riror∗Roβnr∗q˙RFr∗dr∗+∑n=1∞θcπ2ri2βn2+m2βn2J02βnriJ02βnri-J02βnroe-αβn2+m2t-1RoβnrRo′βnri+Tb+q˙mρbωcb.Equation (25) for determination of eigenvalues of the solution is a transcendental equation. In order to determine the roots a numerical procedure was used. The procedure consisted of identifying values of βn in the neighborhood where the term(38)Reiβn=J0βnriY0βnro-J0βnroY0βnriswitches from positive to negative or vice versa. The search of regions where roots are located consisted of calculating Rei(β), starting with β=0 and with increments Δβ=0.1, to find zones of two consecutive β values where Rei switches sign. Once these regions are identified, Newton’s method was run at each one of those regions, trying to find the values of βn that produce a value of (38) closer to zero. The iterative procedure uses the following equation:(39)βni+1=βni-Rei′βniReiβni. The iterations proceed as long as Rei(βni+1)>ϵ, where ϵ is the convergence criteria. When the convergence criteria is met, βni+1 is taken as βn. In this investigation a value ϵ=10-14 was used. Table 2 presents a list of the first 25 eigenvalues obtained by this procedure for ri=0.75mm and ro=10cm. It is important to see that eigenvalues depend only on geometric conditions of the problem. 625 eigenvalues were used in obtaining numerical values of (37).Table 2
First 20 eigenvalues of transcendental equation (25) for this example.
Eigenvalue number
β
n
1
27.7412491590629
2
59.7215157304870
3
91.6538083942138
4
123.5335630443243
5
155.3754428957278
6
187.1895966699011
7
218.9827041355604
8
250.7593126035128
9
282.5226419505083
10
314.2750529359010
11
346.0183280786956
12
377.7538469185017
13
409.4826995869618
14
441.2057629309977
15
472.9237530699340
16
504.6372626347768
17
536.3467877607103
18
568.0527480482156
19
599.7555015887882
20
631.4553564558279The integrals that appear in the solution are evaluated using Simpson’s rule, using a fine mesh of 397 divisions of the regionri<r<ro. The integral obtained with 397 divisions was compared to a calculation of the integral with 794 divisions. The comparison of the integral of (37) for the leading eigenvalue (beta1) gave a difference of 0.0058% between the cases of calculation of the integral with 397 divisions and 794 divisions. It is important to mention that the integrals that appear in (37) are evaluated once for every eigenvalue βn and do not change with r or t.


## 3. Results and Discussion

### 3.1. Comparison with Previous Studies

In order to validate the analytical solution presented in the previous section, a comparison of that solution ast→∞ was made with a steady state solution. A possible validation of the solution obtained here is by comparing it to the analytical solution presented by Haemmerich et al. [8]. Solution presented in [8] is the steady state condition of RF heating of tissue without blood perfusion and metabolic heat generation. Figure 2(a) presents a comparison of steady state solution of Haemmerich et al. [8] versus the solution obtained from (37) as t→∞ for the case where q˙m=0 and ω=0. A very good agreement between these two solutions is obtained.Figure 2
Comparison of solution obtained from (37) (a) as t→∞ for κ=13.25, q˙m=0, and ω=0versus solution obtained from model proposed by Haemmerich et al. [8] and (b) versus a numerical solution by the finite volume method for the conditions of Figure 3 in López Molina et al. [9].
(a)
(b)Figure 3
Temperature distribution of tissue subject to different conditions of blood perfusion during RF ablation for (a)κ=13.35, q˙m=0, and ω=0; (b) κ=25.5, q˙m=700 W m-1, and ω=0.0005 s-1; (c) κ=34.55, q˙m=700 W m-1, and ω=0.0015 s-1. In each case values of κ are adjusted so that the maximum temperature as t→∞ does not surpass 100∘C.
(a)
(b)
(c)The solution was also validated by comparing it to the solution obtained by López Molina et al. [9]. Figure 2(b) presents temperature distributions obtained in this solution for the conditions of Figure 3 in López Molina et al. [9] (ρ=1060kg m-3, C=3600 J kg-1K-1, k=0.502 W m-1K-1, σ=0.33 S m-1, ρb=1000kgm-3, cb=4148 J kg-1K-1, ω=6.41×10-3 s-1, and Tc=5∘C). In [9] the domain is 0.75mm<r<∞ m, while in this work it is 0.75mm<r<10cm, but the voltage was adapted so that a similar current density was used in both cases. The very good agreement between the solution of [9] and the solution obtained here is seen. A numerical code based on the control volume method [14], with an implicit method and 200 uniform finite volumes, was also written; results of a simulation using this code for the case of conditions of Figure 3 of López Molina et al. [9] are shown in Figure 2(b). We found a very good agreement between the analytical solution obtained here and this approximate numerical solution.


### 3.2. Influence of Blood Perfusion

Figure3 shows temperature profiles at different times for three different blood perfusion rates. Figure 3(a) is for q˙m=0 W m-3, ω=0s-1 (no metabolic heat generation, no blood perfusion), and Vi=41.23 V. The voltage was adjusted so that the highest temperature as t→∞ is 100∘C. In this case it is clear that the effect of not having blood perfusion is to produce a large zone of high temperature in the steady state case, but it takes a very long time to reach this steady state; after 1 hour the temperature distribution is still very far from the steady state temperature distribution. It is important to notice also that the application of boundary condition at ro=10cm is arguably incorrect for the case t=∞, while for shorter times it seems to be justified by the nature of the temperature profiles. Figure 3(b) is for q˙m=700 W m-3, ω=0.0005 s-1, and Vi=56.98 V. The voltage here was also adjusted so that the highest temperature as t→∞ is 100∘C. From the figure it is seen that blood perfusion has the effect of narrowing the region where the tissue is at a steady state temperature above a 50∘C threshold. It is also seen that the effect of blood perfusion is to accelerate reaching a steady state, which can also be seen from factor (1-e-α(β2+m2)t) in (37). Steady state is reached when the previous term approaches unity, and since m2 grows with blood perfusion, less time is required to achieve steady state when blood perfusion is high. In the absence of blood perfusion, all heat is transported by diffusion, but, in this case, diffusion is a very inefficient way of transporting heat, because the thermal diffusivity is very low (α=1.219×10-6m2s-1), so it takes a long time for heat generated by electrical currents to redistribute in the tissue. When blood perfusion appears, blood plays the role of removing heat from hotspots in an efficient manner, accelerating in this way the achievement of steady state. Figure 3(c) is for ω=0.0015s-1 and Vi=66.32 V. The voltage here was also adjusted so that the highest temperature as t→∞ is 100∘C. Figure 3(c) shows the accentuation of the effects described before as blood perfusion is further increased.Figure4 presents temperature profiles at different times for q˙m=700Wm-3, and different blood perfusion rates, Vi=75 V and Tc=5∘C. The voltage was not adjusted so that the highest temperature of 100∘C is reached at a finite time, which is indicated in the figures for each case. A good selection of applied voltage is important, a practical time must be chosen, as the surgery times should be reasonable, usually no more than 12 minutes. By comparing the three different blood perfusion cases, we observe that blood perfusion delays the appearance of roll-off, but the temperature profile that is reached at roll-off is very similar in all the three cases shown. It is also observed that for all times, application of boundary condition at r=ro=10cm is justified, because all temperature profiles approach asymptotically Tbasal as r→ro for ro=10cm.Figure 4
Temperature distribution of tissue subject to different conditions of blood perfusion during RF ablation forVi=75 V. (a) q˙m=0 and ω=0; (b) q˙m=700 W m-1 and ω=0.0005 s-1; (c) q˙m=700 W m-1 and ω=0.0015 s-1. Simulation times are chosen so that the maximum temperature does not surpass 100∘C.
(a)
(b)
(c)


### 3.3. Use of Solution as Design Tool

The solution provided by (37) can be manipulated to use it for calibrating the parameters of a RF ablation procedure. For instance, a procedure to determine the voltage that is needed to produce roll-off at a particular time is described here.The first derivative of (37) is obtained, evaluated at a chosen time t=ts, and made equal to zero. (40)0=π22k∑n=1∞1-e-αβn2+m2tsβn2+m2βn2J02βnriJ02βnri-J02βnro∂Roβnr∂r∫riror∗Roβnr∗q˙RFr∗dr∗+∑n=1∞θcπ2ri2βn2+m2βn2J02βnriJ02βnri-J02βnroe-αβn2+m2ts-1Ro′βnrRo′βnri.The previous equation defines the maximum temperature within the tissue. It can be solved to determine the critical radius,rcr, where the maximum temperature occurs. Equation (40) is a transcendental equation, which can also be solved by Newton’s method, described in the previous section. When rcr is determined from the transcendental equation, it can be introduced into (37) evaluated at t=ts to obtain the maximum temperature, Tmax,(41)Tmax=π22k∑n=1∞1-e-αβn2+m2tsβn2+m2βn2J02βnriJ02βnri-J02βnroRoβnrcr∫riror∗Roβnr∗q˙RFr∗dr∗+∑n=1∞θcπ2ri2βn2+m2βn2J02βnriJ02βnri-J02βnroe-αβn2+m2ts-1RoβnrcrRo′βnri.Considering the caseTmax=100∘C and also recalling q˙RF=κ/r2, the value of κ that produces Tmax=100∘C is determined from(42)κs=100-∑n=1∞θcπ2ri/2βn2+m2βn2J02βnri/J02βnri-J02βnroe-αβn2+m2ts-1RoβnrcrRo′βnriπ2/2k∑n=1∞1-e-αβn2+m2ts/βn2+m2βn2J02βnri/J02βnri-J02βnroRoβnrcr∫riror∗Roβnr∗1/r∗2dr∗.Figure5(a) shows the effect of perfusion on voltage that is needed to produce roll-off (Tmax=100∘C) and on volume of tissue damage (defined as the volume of tissue above 50∘C) at different times. This is relevant because by adjusting the voltage it is possible to modify the time to roll-off and with it it is possible to modify the volume of tissue damage. In that figure it is seen that higher voltages are required to produce roll-off at a given time if blood perfusion is increased. Blood perfusion removes heat generated in tissue in an efficient manner, so a higher power is required to raise the tissue temperature to roll-off at a given time.Figure 5
Voltage needed to produce roll-off and volume of tissue damage at that time for different blood perfusion rates.
(a)
(b)Figure5(b) shows evolution of volume of tissue damage (defined as the volume of tissue above a threshold temperature of 50∘C) at roll-off time, which occurs at the voltages of Figure 5(a). For small roll-off times, where the volume of tissue damage is small and the tissue damaged is located where the higher current densities occur (very close to the electrode), blood perfusion is not relevant, because heat generated in that zone outperforms blood perfusion potential to remove heat. Blood perfusion becomes relevant for larger roll-off times, because as we move to regions away from the electrode, not much heat is generated there and blood perfusion is capable of significantly lowering tissue temperature, so reducing the volume of tissue damage when blood perfusion is increased.Figure6(a) shows how the voltage required to reach roll-off at a given time is affected by electrode cooling temperature. In this case, a higher voltage is required to reach roll-off at a given time, if cooling temperature is decreased. Cooling the electrode helps to eliminate heat generated, so a higher power is required to raise the tissue temperature to roll-off at a given time if the electrode temperature is lowered. Figure 6(b) presents a comparison of volume of tissue damage at roll-off times. As expected, a lower electrode cooling temperature produces a larger volume of tissue damage, but the difference is not very significative, at least in the range of Tc values used to generate the graph, which was already found by Haemmerich et al. [8].Figure 6
Voltage needed to produce roll-off and volume of tissue damage at that time for different electrode temperatures.
(a)
(b)


### 3.4. Extension to Multilayer Electrical Conductivity

One advantage of the method of solution presented here is the generality of the method in terms of the type of source termq˙RF that can be considered. q˙RF could be any function of r and could also be a function of time. A case of a spatial variation of electrical conductivity of tissue is presented here, which would correspond with a model based on two compartments of different tissues, employed for instance by Zorbas and Samaras [10] to study RF ablation in a tumor surrounded by healthy tissue.Consider the case where the tissue electrical conductivity is not uniform but is layered in two regions:(43)σ=σ1forri≤r≤rδ,σ=σ2forrδ≤r≤ro,where rδ is the external radius of the first layer of tissue, which has an electrical conductivity that is different to the healthy tissue electrical conductivity. In that case, from (2) and the appropriate boundary conditions, expressions for q˙RF of (3) can be obtained:(44)q˙RF=C2σ1r2forri≤r≤rδ,q˙RF=C2σ2r2forrδ≤r≤ro,where (45)C=σ1σ2Viσ2ln⁡ri+σ2-σ1ln⁡rδ-σ1ln⁡ro.Expression (44) can be introduced into (37) and analytical solution is obtained for that case. A new parameter γ=σ1/σ2, where σ2=0.188 S m-1, was introduced (σ2 is a value usually reported for healthy hepatic tissue).Figure7 shows steady state temperature distributions for different values of γ and rδ. Values of γ larger than 1 would model RF ablation of a tumor, since its electrical conductivity is higher than that of healthy tissue [15, 16]. In contrast, values of γ smaller than 1 would model an RF ablation where the target tissue has an electrical conductivity smaller than that of remote tissue, as observed in bone tumor ablation [17]. In all cases shown in Figure 7  Vi=56.87 V, which is the voltage that produces roll-off at t→∞ in the case γ=3 and rδ=3cm. The effect of different electrical conductivity in a first layer of tissue ri≤r≤rδ is to sensibly modify the temperature distribution within the tissue. It is seen that for large values of γ a second hump is forming around r=rδ, modifying significantly the volume of tissue whose temperature is above 50∘C. In Figures 7(a) and 7(b) it is seen that the case γ=2 has a higher maximum temperature than the case γ=3, but the temperature far from the electrode is higher for γ=3. The explanation we found for this behavior is that as γ is increased, less heat is produced in the zone r<rδ, but more heat is generated in the region r>rδ because the overall impedance of the tissue decreases, and the current density is larger. That behavior does not present in the case of Figure 7(c), where the case γ=3 has a higher overall temperature, compared to cases with smaller γ; the second hump is also less notorious in Figure 7(c), because in that case at r>rδ the current density has decreased as to not produce a large second hump.Figure 7
Temperature distribution of tissue with two different layers subject to RF ablation. (a)rδ=1cm, (b) rδ=2cm, and (c) rδ=3cm. For all cases Vi=56.87 V. rδ is the external radius of the first layer of tissue and γ is the ratio of electrical conductivities of first to second layers of tissue, where σ2=0.188 S m-1.
(a)
(b)
(c)In a similar way, the solution expressed by (37) can be used for determination of the temperature distribution in tissue for different types of source terms and dependencies of source term with time and space.


### 3.5. Discussion

The present model is limited to constant properties of tissue and for that reason does not consider temperature dependent properties and does not work for the region where the tissue vaporizes. Our previous experience in computer modeling in relation to the temperature dependence of thermal conductivity is that its effect is really negligible before tissue vaporization [18]. The electrical conductivity increases its value 2%/°C in the range 37∘C ≤T≤100∘C, producing a kind of positive-feedback which basically speeds up the heating.The 1-D model used here does not model the particularities of the electrode tip and connection with the plastic catheter; tissue at those zones absorbs higher RF power (edge effect) and hence is desiccated during the first seconds. In contrast, in the 1-D model, the RF power is equally distributed along all the tissue next to the electrode surface. For this reason the thermal dynamics of a 1-D model is different to that of 2-D/3-D models. The utility of the 1-D model proposed here is to have a fast and easy-to-use tool to study the effect of the design parameters involved in RF ablation with internally cooled needle-like electrodes (see Section3.3).Results from Figure3 cannot be compared to experimental results since they were computed by adjusting applied voltage to avoid reaching 100∘C, and this is not the procedure employed by the RF generators. Likewise, although results from Figure 4 correspond with conditions similar to those employed by RF generators (since they were computed by fixing the applied voltage of 75 V), the simulations were stopped once 100∘C was reached (roll-off). RF generators really employ a mode called “constant impedance” which consists of applying RF pulses once roll-off has occurred. Thermal lesion continues growing after that instant [19]. For these reasons, the diameters of the thermal lesion (assessed from plots of Figures 3 and 4 by the 50∘C isotherm) cannot be directly compared to the experimental results.However, computer results from the 1-D model are quantitatively in agreement with experimental findings. For instance, Figure5(a) indicates that regardless of the perfusion rate, when a high voltage of 95 V is used, roll-off is reached after the first minute (around 60 s). This coincides exactly with what is observed in the clinical practice, where it is usual to work at maximum power, which implies an applied voltage of 100 V. When lower voltages are used, the roll-off importantly delays, as shown in Figure 5(a) and experimentally observed in [20].The model obtained here is 1-D, while most of the experimental and numerical results available are for 2-D/3-D cases. The plane at the mid-height of the active part of the electrode represents the only part of the electrode that behaves quasi one dimensionally, because of symmetry. It has been observed [3] that the last part of the tissue to reach desiccation is precisely the plane at the mid-height of the electrode, so it determines time to roll-off. The numerical results presented by [21] show that for an applied voltage of 80 V, when no blood perfusion is considered, roll-off occurs at 130 seconds. This is qualitatively similar to the results observed in Figure 5(a) of our paper for the case of no blood perfusion and an applied voltage of 80 V, where time to roll-off is 150 seconds.Similarly, Figure6(b) suggest that the effect of the electrode temperature is of little importance in terms of thermal lesion volume, which has also been observed in experiments [8].Figure7 suggests that when the tissue layer closest to the electrode surface is more electrically conductive (σ1>σ2, i.e., γ>1), the temperature profiles for a fixed applied voltage are higher, which means that larger thermal lesions could be created. In a real scenario, tumor tissue is more conductive than the healthy tissue [16]. Overall, all these facts suggest that, due to the higher σ of the tumor, larger lesions can be created.


## 3.1. Comparison with Previous Studies

In order to validate the analytical solution presented in the previous section, a comparison of that solution ast→∞ was made with a steady state solution. A possible validation of the solution obtained here is by comparing it to the analytical solution presented by Haemmerich et al. [8]. Solution presented in [8] is the steady state condition of RF heating of tissue without blood perfusion and metabolic heat generation. Figure 2(a) presents a comparison of steady state solution of Haemmerich et al. [8] versus the solution obtained from (37) as t→∞ for the case where q˙m=0 and ω=0. A very good agreement between these two solutions is obtained.Figure 2
Comparison of solution obtained from (37) (a) as t→∞ for κ=13.25, q˙m=0, and ω=0versus solution obtained from model proposed by Haemmerich et al. [8] and (b) versus a numerical solution by the finite volume method for the conditions of Figure 3 in López Molina et al. [9].
(a)
(b)Figure 3
Temperature distribution of tissue subject to different conditions of blood perfusion during RF ablation for (a)κ=13.35, q˙m=0, and ω=0; (b) κ=25.5, q˙m=700 W m-1, and ω=0.0005 s-1; (c) κ=34.55, q˙m=700 W m-1, and ω=0.0015 s-1. In each case values of κ are adjusted so that the maximum temperature as t→∞ does not surpass 100∘C.
(a)
(b)
(c)The solution was also validated by comparing it to the solution obtained by López Molina et al. [9]. Figure 2(b) presents temperature distributions obtained in this solution for the conditions of Figure 3 in López Molina et al. [9] (ρ=1060kg m-3, C=3600 J kg-1K-1, k=0.502 W m-1K-1, σ=0.33 S m-1, ρb=1000kgm-3, cb=4148 J kg-1K-1, ω=6.41×10-3 s-1, and Tc=5∘C). In [9] the domain is 0.75mm<r<∞ m, while in this work it is 0.75mm<r<10cm, but the voltage was adapted so that a similar current density was used in both cases. The very good agreement between the solution of [9] and the solution obtained here is seen. A numerical code based on the control volume method [14], with an implicit method and 200 uniform finite volumes, was also written; results of a simulation using this code for the case of conditions of Figure 3 of López Molina et al. [9] are shown in Figure 2(b). We found a very good agreement between the analytical solution obtained here and this approximate numerical solution.


## 3.2. Influence of Blood Perfusion

Figure3 shows temperature profiles at different times for three different blood perfusion rates. Figure 3(a) is for q˙m=0 W m-3, ω=0s-1 (no metabolic heat generation, no blood perfusion), and Vi=41.23 V. The voltage was adjusted so that the highest temperature as t→∞ is 100∘C. In this case it is clear that the effect of not having blood perfusion is to produce a large zone of high temperature in the steady state case, but it takes a very long time to reach this steady state; after 1 hour the temperature distribution is still very far from the steady state temperature distribution. It is important to notice also that the application of boundary condition at ro=10cm is arguably incorrect for the case t=∞, while for shorter times it seems to be justified by the nature of the temperature profiles. Figure 3(b) is for q˙m=700 W m-3, ω=0.0005 s-1, and Vi=56.98 V. The voltage here was also adjusted so that the highest temperature as t→∞ is 100∘C. From the figure it is seen that blood perfusion has the effect of narrowing the region where the tissue is at a steady state temperature above a 50∘C threshold. It is also seen that the effect of blood perfusion is to accelerate reaching a steady state, which can also be seen from factor (1-e-α(β2+m2)t) in (37). Steady state is reached when the previous term approaches unity, and since m2 grows with blood perfusion, less time is required to achieve steady state when blood perfusion is high. In the absence of blood perfusion, all heat is transported by diffusion, but, in this case, diffusion is a very inefficient way of transporting heat, because the thermal diffusivity is very low (α=1.219×10-6m2s-1), so it takes a long time for heat generated by electrical currents to redistribute in the tissue. When blood perfusion appears, blood plays the role of removing heat from hotspots in an efficient manner, accelerating in this way the achievement of steady state. Figure 3(c) is for ω=0.0015s-1 and Vi=66.32 V. The voltage here was also adjusted so that the highest temperature as t→∞ is 100∘C. Figure 3(c) shows the accentuation of the effects described before as blood perfusion is further increased.Figure4 presents temperature profiles at different times for q˙m=700Wm-3, and different blood perfusion rates, Vi=75 V and Tc=5∘C. The voltage was not adjusted so that the highest temperature of 100∘C is reached at a finite time, which is indicated in the figures for each case. A good selection of applied voltage is important, a practical time must be chosen, as the surgery times should be reasonable, usually no more than 12 minutes. By comparing the three different blood perfusion cases, we observe that blood perfusion delays the appearance of roll-off, but the temperature profile that is reached at roll-off is very similar in all the three cases shown. It is also observed that for all times, application of boundary condition at r=ro=10cm is justified, because all temperature profiles approach asymptotically Tbasal as r→ro for ro=10cm.Figure 4
Temperature distribution of tissue subject to different conditions of blood perfusion during RF ablation forVi=75 V. (a) q˙m=0 and ω=0; (b) q˙m=700 W m-1 and ω=0.0005 s-1; (c) q˙m=700 W m-1 and ω=0.0015 s-1. Simulation times are chosen so that the maximum temperature does not surpass 100∘C.
(a)
(b)
(c)


## 3.3. Use of Solution as Design Tool

The solution provided by (37) can be manipulated to use it for calibrating the parameters of a RF ablation procedure. For instance, a procedure to determine the voltage that is needed to produce roll-off at a particular time is described here.The first derivative of (37) is obtained, evaluated at a chosen time t=ts, and made equal to zero. (40)0=π22k∑n=1∞1-e-αβn2+m2tsβn2+m2βn2J02βnriJ02βnri-J02βnro∂Roβnr∂r∫riror∗Roβnr∗q˙RFr∗dr∗+∑n=1∞θcπ2ri2βn2+m2βn2J02βnriJ02βnri-J02βnroe-αβn2+m2ts-1Ro′βnrRo′βnri.The previous equation defines the maximum temperature within the tissue. It can be solved to determine the critical radius,rcr, where the maximum temperature occurs. Equation (40) is a transcendental equation, which can also be solved by Newton’s method, described in the previous section. When rcr is determined from the transcendental equation, it can be introduced into (37) evaluated at t=ts to obtain the maximum temperature, Tmax,(41)Tmax=π22k∑n=1∞1-e-αβn2+m2tsβn2+m2βn2J02βnriJ02βnri-J02βnroRoβnrcr∫riror∗Roβnr∗q˙RFr∗dr∗+∑n=1∞θcπ2ri2βn2+m2βn2J02βnriJ02βnri-J02βnroe-αβn2+m2ts-1RoβnrcrRo′βnri.Considering the caseTmax=100∘C and also recalling q˙RF=κ/r2, the value of κ that produces Tmax=100∘C is determined from(42)κs=100-∑n=1∞θcπ2ri/2βn2+m2βn2J02βnri/J02βnri-J02βnroe-αβn2+m2ts-1RoβnrcrRo′βnriπ2/2k∑n=1∞1-e-αβn2+m2ts/βn2+m2βn2J02βnri/J02βnri-J02βnroRoβnrcr∫riror∗Roβnr∗1/r∗2dr∗.Figure5(a) shows the effect of perfusion on voltage that is needed to produce roll-off (Tmax=100∘C) and on volume of tissue damage (defined as the volume of tissue above 50∘C) at different times. This is relevant because by adjusting the voltage it is possible to modify the time to roll-off and with it it is possible to modify the volume of tissue damage. In that figure it is seen that higher voltages are required to produce roll-off at a given time if blood perfusion is increased. Blood perfusion removes heat generated in tissue in an efficient manner, so a higher power is required to raise the tissue temperature to roll-off at a given time.Figure 5
Voltage needed to produce roll-off and volume of tissue damage at that time for different blood perfusion rates.
(a)
(b)Figure5(b) shows evolution of volume of tissue damage (defined as the volume of tissue above a threshold temperature of 50∘C) at roll-off time, which occurs at the voltages of Figure 5(a). For small roll-off times, where the volume of tissue damage is small and the tissue damaged is located where the higher current densities occur (very close to the electrode), blood perfusion is not relevant, because heat generated in that zone outperforms blood perfusion potential to remove heat. Blood perfusion becomes relevant for larger roll-off times, because as we move to regions away from the electrode, not much heat is generated there and blood perfusion is capable of significantly lowering tissue temperature, so reducing the volume of tissue damage when blood perfusion is increased.Figure6(a) shows how the voltage required to reach roll-off at a given time is affected by electrode cooling temperature. In this case, a higher voltage is required to reach roll-off at a given time, if cooling temperature is decreased. Cooling the electrode helps to eliminate heat generated, so a higher power is required to raise the tissue temperature to roll-off at a given time if the electrode temperature is lowered. Figure 6(b) presents a comparison of volume of tissue damage at roll-off times. As expected, a lower electrode cooling temperature produces a larger volume of tissue damage, but the difference is not very significative, at least in the range of Tc values used to generate the graph, which was already found by Haemmerich et al. [8].Figure 6
Voltage needed to produce roll-off and volume of tissue damage at that time for different electrode temperatures.
(a)
(b)


## 3.4. Extension to Multilayer Electrical Conductivity

One advantage of the method of solution presented here is the generality of the method in terms of the type of source termq˙RF that can be considered. q˙RF could be any function of r and could also be a function of time. A case of a spatial variation of electrical conductivity of tissue is presented here, which would correspond with a model based on two compartments of different tissues, employed for instance by Zorbas and Samaras [10] to study RF ablation in a tumor surrounded by healthy tissue.Consider the case where the tissue electrical conductivity is not uniform but is layered in two regions:(43)σ=σ1forri≤r≤rδ,σ=σ2forrδ≤r≤ro,where rδ is the external radius of the first layer of tissue, which has an electrical conductivity that is different to the healthy tissue electrical conductivity. In that case, from (2) and the appropriate boundary conditions, expressions for q˙RF of (3) can be obtained:(44)q˙RF=C2σ1r2forri≤r≤rδ,q˙RF=C2σ2r2forrδ≤r≤ro,where (45)C=σ1σ2Viσ2ln⁡ri+σ2-σ1ln⁡rδ-σ1ln⁡ro.Expression (44) can be introduced into (37) and analytical solution is obtained for that case. A new parameter γ=σ1/σ2, where σ2=0.188 S m-1, was introduced (σ2 is a value usually reported for healthy hepatic tissue).Figure7 shows steady state temperature distributions for different values of γ and rδ. Values of γ larger than 1 would model RF ablation of a tumor, since its electrical conductivity is higher than that of healthy tissue [15, 16]. In contrast, values of γ smaller than 1 would model an RF ablation where the target tissue has an electrical conductivity smaller than that of remote tissue, as observed in bone tumor ablation [17]. In all cases shown in Figure 7  Vi=56.87 V, which is the voltage that produces roll-off at t→∞ in the case γ=3 and rδ=3cm. The effect of different electrical conductivity in a first layer of tissue ri≤r≤rδ is to sensibly modify the temperature distribution within the tissue. It is seen that for large values of γ a second hump is forming around r=rδ, modifying significantly the volume of tissue whose temperature is above 50∘C. In Figures 7(a) and 7(b) it is seen that the case γ=2 has a higher maximum temperature than the case γ=3, but the temperature far from the electrode is higher for γ=3. The explanation we found for this behavior is that as γ is increased, less heat is produced in the zone r<rδ, but more heat is generated in the region r>rδ because the overall impedance of the tissue decreases, and the current density is larger. That behavior does not present in the case of Figure 7(c), where the case γ=3 has a higher overall temperature, compared to cases with smaller γ; the second hump is also less notorious in Figure 7(c), because in that case at r>rδ the current density has decreased as to not produce a large second hump.Figure 7
Temperature distribution of tissue with two different layers subject to RF ablation. (a)rδ=1cm, (b) rδ=2cm, and (c) rδ=3cm. For all cases Vi=56.87 V. rδ is the external radius of the first layer of tissue and γ is the ratio of electrical conductivities of first to second layers of tissue, where σ2=0.188 S m-1.
(a)
(b)
(c)In a similar way, the solution expressed by (37) can be used for determination of the temperature distribution in tissue for different types of source terms and dependencies of source term with time and space.


## 3.5. Discussion

The present model is limited to constant properties of tissue and for that reason does not consider temperature dependent properties and does not work for the region where the tissue vaporizes. Our previous experience in computer modeling in relation to the temperature dependence of thermal conductivity is that its effect is really negligible before tissue vaporization [18]. The electrical conductivity increases its value 2%/°C in the range 37∘C ≤T≤100∘C, producing a kind of positive-feedback which basically speeds up the heating.The 1-D model used here does not model the particularities of the electrode tip and connection with the plastic catheter; tissue at those zones absorbs higher RF power (edge effect) and hence is desiccated during the first seconds. In contrast, in the 1-D model, the RF power is equally distributed along all the tissue next to the electrode surface. For this reason the thermal dynamics of a 1-D model is different to that of 2-D/3-D models. The utility of the 1-D model proposed here is to have a fast and easy-to-use tool to study the effect of the design parameters involved in RF ablation with internally cooled needle-like electrodes (see Section3.3).Results from Figure3 cannot be compared to experimental results since they were computed by adjusting applied voltage to avoid reaching 100∘C, and this is not the procedure employed by the RF generators. Likewise, although results from Figure 4 correspond with conditions similar to those employed by RF generators (since they were computed by fixing the applied voltage of 75 V), the simulations were stopped once 100∘C was reached (roll-off). RF generators really employ a mode called “constant impedance” which consists of applying RF pulses once roll-off has occurred. Thermal lesion continues growing after that instant [19]. For these reasons, the diameters of the thermal lesion (assessed from plots of Figures 3 and 4 by the 50∘C isotherm) cannot be directly compared to the experimental results.However, computer results from the 1-D model are quantitatively in agreement with experimental findings. For instance, Figure5(a) indicates that regardless of the perfusion rate, when a high voltage of 95 V is used, roll-off is reached after the first minute (around 60 s). This coincides exactly with what is observed in the clinical practice, where it is usual to work at maximum power, which implies an applied voltage of 100 V. When lower voltages are used, the roll-off importantly delays, as shown in Figure 5(a) and experimentally observed in [20].The model obtained here is 1-D, while most of the experimental and numerical results available are for 2-D/3-D cases. The plane at the mid-height of the active part of the electrode represents the only part of the electrode that behaves quasi one dimensionally, because of symmetry. It has been observed [3] that the last part of the tissue to reach desiccation is precisely the plane at the mid-height of the electrode, so it determines time to roll-off. The numerical results presented by [21] show that for an applied voltage of 80 V, when no blood perfusion is considered, roll-off occurs at 130 seconds. This is qualitatively similar to the results observed in Figure 5(a) of our paper for the case of no blood perfusion and an applied voltage of 80 V, where time to roll-off is 150 seconds.Similarly, Figure6(b) suggest that the effect of the electrode temperature is of little importance in terms of thermal lesion volume, which has also been observed in experiments [8].Figure7 suggests that when the tissue layer closest to the electrode surface is more electrically conductive (σ1>σ2, i.e., γ>1), the temperature profiles for a fixed applied voltage are higher, which means that larger thermal lesions could be created. In a real scenario, tumor tissue is more conductive than the healthy tissue [16]. Overall, all these facts suggest that, due to the higher σ of the tumor, larger lesions can be created.


## 4. Conclusions

This study presents an analytical solution for the RF ablation of tissue from a cooled cylindrical electrode. The solution matches well with other available analytical solutions and with approximate numerical solutions to the problem. Unlike other available time dependent analytical solutions, this is for a finite spatial domain and is more general, because of the possibility of handling different functional dependencies of the source term and nonhomogeneous boundary condition. Application of boundary condition at the far field,r=ro=10cm, is justified for all cases where perfusion is large and for short times in cases where perfusion is zero. The solution provides a tool for better understanding RF ablation, for validating codes used for numerical modeling of the problem, and can be used as a tool for determination of the set of input parameters that would produce some expected outcome from the RF ablation procedure.



---

*Source: 9021616-2017-03-06.xml*