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.
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
Figure 1 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 and outer radius . In some studied cases, here we also consider a two-compartment model, with cancerous tissue located in the zone , while the healthy tissue is in the region , represented by Figure 1(b). Inside this annular section of tissue is an active RF electrode of radius and the dispersive electrode is considered to be at the periphery of the tissue.
(a)
(b)
Consider then the case of a fragment of tissue in the shape of a concentric annulus with inner radius and outer radius . The tissue is originally at the basal temperature , where is the metabolic heat generation of the tissue [W ], C is the blood arterial temperature, and and are the blood density [kg ] and specific heat [J ], while is the blood perfusion rate []. The inner radius of this section is an RF active electrode held at a low temperature while the outer radius is at the unaltered temperature . The inner radius also has an applied voltage , 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 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 ] is considered for the range of temperatures C. The domain of the solution is . Table 1 presents the set of properties and geometrical parameters used throughout, unless specified otherwise when required. For most cases was taken as C, unless specified otherwise.
The governing equation for this problem is Pennes’ Bioheat equation [12]:where , , and are the tissue thermal conductivity [W ], density [kg ], and specific heat [J ], is the metabolic heat generation per unit volume within the tissue, and is the blood perfusion heat transport term. is the volumetric heat generation due to circulation of electrical currents in the RF range [W ].
The calculation of requires solution of the electrical problem. This is stated in terms of Laplace’s equation:where is the electrical conductivity of tissue.
The heat generation per unit volume from RF currents circulation, , that appears in (1) is determined fromThe electrical potential [W ] and electrical current density [A ] can be determined fromso that
Considering that in this problem the voltage is only a function of the radius and is constant, (2) can be stated assubject to boundary conditions
With these equation and set of boundary conditions, the solution can be stated as
The electrical potential and electrical current density can be determined by introducing (8) into (4):
And the heat generation per unit volume is determined fromwhere .
2.2. Analytical Solution
Considering a temperature of tissue , (1) can be written assubject to
Using a change of variable , (11) can be written assubject to
This problem is nonhomogeneous in the equation itself, because of term and also because of the nonhomogeneous boundary condition at . The problem is solved here using Green’s functions [13]. The solution in terms of Green’s functions can be stated as where is Green’s function. In this case , , and . is the tissue thermal diffusivity [ ].
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 equationsubject towhere 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 variablethe equation can be written aswhere is a constant that later will be related to the eigenvalues of the problem.
By defining and dividing every term by two separated equations are produced:subject towhere has been introduced. The solution to this equation is where the eigenvalues are the positive roots of equation
Equation (22) is called Bessel’s differential equation of order zero, whose solutions and are called Bessel’s functions of order zero of the first and second kind, respectively.
The equation for isThe general solution to this equation is
Writing again and considering summation over all possible eigenvalues , we get
Applying the initial condition expressed in (19), can be written as
The solution to the problem in terms of Green’s function can be stated as
A comparison of terms, after replacing by , provide Green’s function for this problem
Once Green’s function for the problem is obtained, the solution to (13) and (14) can be stated as whereThe first term of (32) accounts for the source term , while the last term accounts for the nonhomogeneous boundary condition at . This last term is evaluated aswhere
Since, in this case , the integral to account for the nonhomogeneous boundary condition at can be evaluated easily, then
The final form of the solution is
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 in the neighborhood where the termswitches from positive to negative or vice versa. The search of regions where roots are located consisted of calculating , starting with and with increments , to find zones of two consecutive values where switches sign. Once these regions are identified, Newton’s method was run at each one of those regions, trying to find the values of that produce a value of (38) closer to zero. The iterative procedure uses the following equation: The iterations proceed as long as , where is the convergence criteria. When the convergence criteria is met, is taken as . In this investigation a value was used. Table 2 presents a list of the first 25 eigenvalues obtained by this procedure for and . 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).
The integrals that appear in the solution are evaluated using Simpson’s rule, using a fine mesh of 397 divisions of the region . 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 () 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 and do not change with or .
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 as 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 for the case where and . A very good agreement between these two solutions is obtained.
(a)
(b)
(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] ( , J , W , S , , J , , and C). In [9] the domain is m, while in this work it is , 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
Figure 3 shows temperature profiles at different times for three different blood perfusion rates. Figure 3(a) is for W , (no metabolic heat generation, no blood perfusion), and V. The voltage was adjusted so that the highest temperature as is 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 is arguably incorrect for the case , while for shorter times it seems to be justified by the nature of the temperature profiles. Figure 3(b) is for W m, , and V. The voltage here was also adjusted so that the highest temperature as is 100C. 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 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 in (37). Steady state is reached when the previous term approaches unity, and since 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 ), 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 and V. The voltage here was also adjusted so that the highest temperature as is C. Figure 3(c) shows the accentuation of the effects described before as blood perfusion is further increased.
Figure 4 presents temperature profiles at different times for W, and different blood perfusion rates, V and C. The voltage was not adjusted so that the highest temperature of 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 is justified, because all temperature profiles approach asymptotically as for .
(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 , and made equal to zero.
The previous equation defines the maximum temperature within the tissue. It can be solved to determine the critical radius, , 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 is determined from the transcendental equation, it can be introduced into (37) evaluated at to obtain the maximum temperature, ,
Considering the case C and also recalling , the value of that produces C is determined from
Figure 5(a) shows the effect of perfusion on voltage that is needed to produce roll-off (C) and on volume of tissue damage (defined as the volume of tissue above 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.
(a)
(b)
Figure 5(b) shows evolution of volume of tissue damage (defined as the volume of tissue above a threshold temperature of 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.
Figure 6(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 values used to generate the graph, which was already found by Haemmerich et al. [8].
(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 term that can be considered. could be any function of 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:where 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 of (3) can be obtained:where
Expression (44) can be introduced into (37) and analytical solution is obtained for that case. A new parameter , where S , was introduced ( is a value usually reported for healthy hepatic tissue).
Figure 7 shows steady state temperature distributions for different values of and . 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 V, which is the voltage that produces roll-off at in the case and . The effect of different electrical conductivity in a first layer of tissue is to sensibly modify the temperature distribution within the tissue. It is seen that for large values of a second hump is forming around , modifying significantly the volume of tissue whose temperature is above C. In Figures 7(a) and 7(b) it is seen that the case has a higher maximum temperature than the case , but the temperature far from the electrode is higher for . The explanation we found for this behavior is that as is increased, less heat is produced in the zone , but more heat is generated in the region 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 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 the current density has decreased as to not produce a large second hump.
(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 C 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 Section 3.3).
Results from Figure 3 cannot be compared to experimental results since they were computed by adjusting applied voltage to avoid reaching 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 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 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, Figure 5(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, Figure 6(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].
Figure 7 suggests that when the tissue layer closest to the electrode surface is more electrically conductive (, i.e., ), 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, , 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.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work was supported by the Universidad Autónoma de San Luis Potosí (México), which granted Ricardo Romero-Méndez a sabbatical leave to do research in the field of biomedical engineering, and the Government of Spain through the “Plan Estatal de Investigación, Desarrollo e Innovación Orientada a los Retos de la Sociedad” (Grants TEC 2014-52383-C3-R and TEC 2014-52383-C3-1-R).