Investigation of radiofrequency ablation process in liver tissue by finite element modeling and experiment

Summary. Background. The character of ablation processes with high-frequency electrical current is similar in most biological tissues; however, quantitative characteristics are very different. Consequently, mathematical models of the process have a lot of specific aspects. In this study, we developed mathematical model of radiofrequency ablation in liver tissues with experimental validation of model in ex vivo porcine liver. Methods. The finite element nonlinear computational model for the simulation of the radiofrequency ablation processes and taking into account coupled electrical and thermal phenomena has been developed. The radiofrequency electric current processes are dominated by the active electric conductivity. The heat generation in biological tissues is determined by the electric current density. Simultaneously, the conductivity of the tissue is nonlinearly dependent upon the temperature of the tissue. The model has been implemented in COMSOL Multiphysics computational environment. Tests on physical characteristics of the thermal effect in ex vivo liver tissue have been performed and results compared. Results. Two oval-shaped zones of total and relative tissue destruction were highlighted. The principal distribution of the thermal effect is congruous with the theoretical model; however, the discrepancy of temperatures in experimental and theoretical models increases distally from active perfusion electrode. Conclusions. Distribution of the thermal effect is congruous in the theoretical and experimental model; however, discrepancies of temperatures imply certain inadequacies of the mathematical models. Differences of computed and actual temperatures should be regarded predicting tissue ablation in clinical setting.


Background
Several kinds of tissue ablation by means of highfrequency electrical currents are employed in biomedical practice. The majority of patients with primary or metastatic malignancies confined to the liver are not candidates for resection because of tumor size, location, multifocal character, or inadequate functional hepatic reserve. Local application of heat is tumoricidal (1). Radiofrequency ablation (RFA) is a new local thermal ablative technique for the treatment of unresectable hepatic tumors including hepatocellular carcinoma (2).
Radiofrequency ablation takes place at frequencies of 460-550 kHz. Electrical voltage is being applied to the ablation zone by means of an active probe (needle). The flow of electrical current transfers the energy to the tissue producing temperature reaching 100°C during the interval of 10-15 min (3,4). It is known that at frequency of 500 kHz, the electrical conductivity of liver tissue is approximately 0.148 1/Wm (5) and is essentially defined by the active electrical conductance. The basic physical processes taking place during RFA are electrical current flowing in the tissue and causing volumetric heat generation, as well as heat exchange in thermally conductive tissues. The electrical conductivity is dependent on the temperature of the tissue, thus fully coupled electro-thermal problem is formulated (6). Real RFA processes are even more complex because of the flow of sodium chloride solution injected through the tip of the probe. The filtration of the solution through the tissue and the blood flow in neighboring vessels create additional thermal losses. The analysis can be reasonably simplified by assuming the quasi-static model of electrical conductance. This implies the consideration of alternating voltage and electrical current in terms of their effective values obtained by solving the equations of direct current in conductive media. The lesion of the tissue as a consequence of heating can be evaluated by using the Arrhenius formula (7).
An alternative method of ablation by using highfrequency electric currents is microwave hyperthermia, which takes place at frequencies from 27 MHz to 2.45 GHz. The corresponding values of electrical conductivity at such frequencies are from 0.382 1/Wm to 1.687 1/Wm. Though based on the same physical principle, the two processes are very different from the mathematical point of view. During microwave hyperthermia process, the capacity-dominated electrical conductance is prevailing, and electrical process is governed by Maxwell's equations. In this work, we considered the RFA processes only.
The character of ablation processes is similar in most biological tissues and structures. However, quantitative characteristics are very different when applying the technique in heart (8), liver (6,7), and other biological structures (9). Different applications imply different biomedical accents, different physical constants are necessary, and biological reaction of different tissues to the thermal field may also appear to be different. Consequently, mathematical models of the process have many specific aspects in each application.
The idea of this study is closest to recent publications (6,7) considering the distribution of multiphysical fields in the vicinity of the tip of the ablation electrode. State-of-the-art finite element models generated in COMSOL Multiphysics computational environment are used for theoretical investigations. In this study, the mathematical model of RFA ablation in liver tissues was followed by the experimental validation of model in ex vivo porcine liver.

Mathematical model of RF ablation
Governing equations. The scheme of the investigated domain is presented in Fig. 1A.
Electric current in tissues is described by using the quasi-static model of electrical conductance, the equation of which for axially symmetric formulation reads as [1] Where j -electric potential, [V]; s -conductivity of the tissue, ; r, z -radial and axial coordinates, [m].
Multiplier term 2pr appears in the equation due to axial symmetry of the domain and implies that integrals are calculated in the cylindrical system of coordinates. Therefore, we cancel out neither multiplier 2p from equation [1] nor multiplier term 2pr from subsequent relations, which describe the boundary conditions. The expression 2prs could be rather interpreted as equivalent electrical conductivity coefficient used in the axially symmetric formulation if integrals are calculated over flat triangle elements.
Cauchy (Neuman) boundary condition [2] is used at the insulated boundary of the probe, at the free surface of the tissue and the line of axial symmetry inside of the tissue. Dirichlet boundary condition [3] applies to the surface at the electrically active part of the probe [j E given] and to the external boundary of the investigated domain situated "far enough" from the probe [j E = 0], Fig. 1B. Thermal field in tissues is described by the heat conductance equation as [4] Where T is temperature, Cauchy boundary condition defines the thermal exchange with the environment of the investigated domain: Where a -film coefficient, [W/(Km 2 )]; T¥ -the temperature of the environment, [K]; q -heat flux density, [W/m 2 ]. Equation [5] includes two physically different boundary conditions. Term 2prq indicates a prescribed thermal flux density across the boundary, and term 2pra (T -T¥) describes thermal flux density across the boundary surface of the domain due to heat exchange across the interface between two environments. Therefore, q and a cannot be nonzero at the Medicina (Kaunas) 2007; 43 (4) same part of the surface.
The Dirichlet boundary condition [6] defines the known temperature of the tissue at the distant boundary, where T p = 310 K = 37°C (Fig. 1B).
The coupling of the electric and thermal fields takes place due to the dependency of the electrical conductivity in equation [1] on the temperature of the tissue and due to the volumetric heat source, the magnitude of which is dependent on the electric current strength.
Values of electrical conductivity of biological tissues at different temperatures are not comprehensively investigated at the time being. In practical calculations, the electrical conductivity of the tissue is usually  (4) identified with the electrical conductivity s (T, N) of the sodium chloride solution of normality N and temperature T (°C), which has been expressed by A. Strogyn (10) as: [7] Where is the electrical conductivity of the normality N solution at temperature of 25°C.
The expression of the volumetric heat source reads as , [8] Where s | j| is the power generated by the flow of electric current through the tissues, r b c b w(T-T p )heat power losses due to blood flow perfusion, c b , r b -mass thermal capacity and mass density of the blood, w -perfusion coefficient, Q m -heat power generated due to metabolic processes.
In our analysis of the process in liver tissues, the following physical constants have been used (6) ( Table 1). The damage of the tissue due to heating is obtained by Arrhenius formula (7): , [9] Where c(t) -concentration of live cells; R -universal gas constant; A -"frequency" coefficient, (s -1 ); DE -energy of initiation of irreversible ablation reaction, (J/mol).

Table 1. Physical constants used in the analysis of radiofrequency ablation processes in liver tissues
Temperature of surrounding tissues Parameters A and DE have been characterized for liver tissues by Jacques et al. (11). Damage integral value W=1 corresponds to 63% probability of cell death at a specific point, and value W=4.6 corresponds to 99% probability of cell death at a point in the model. The significance of W=1 has been reported as the point at which tissue coagulation first occurs (12).
General form of equations in COMSOL Multiphysics computational environment. The partial differential equations [1], [4] with boundary conditions [2], [3] and [5], [6] and relation [10] are presented in the "general" form as [10] with boundary conditions , The model of the RFA process presented in COMSOL Multiphysics by symbolic representations [10,11] fully describes the electrical-thermal behavior of the tissue in time by taking into account the nonlinear of the two fields. In this study, our theoretical and experimental investigations were restricted to the analysis of steady distributions of thermal and electrical fields; therefore, the transient terms have been omitted by setting to zero the coefficient matrix d a and using solvers for the steady problem.
Experimentally steady values of temperatures are obtained after the electric current flows through tissues for a very long time. The minimum duration of the time interval is determined experimentally by fixing the measured values at the time moment, after which they do not change in time any more.

Experimental model of RF ablation in ex vivo tissue
Six cadaveric porcine livers (average weight 2300 g) were utilized for the experiment. RFA procedure was performed at the room temperature. The RF delivery system was an ELEKTROTOM HiTT 106 generator system (Berchtold, Tuttlingen, Germany) in mono-polar mode, used in conjunction with a single needle 16-gauge cooled-tip ablation electrode that is used in clinical practice for liver tumor ablation. The probe is 12.0 cm in length with 2 cm uninsulated tip. The distal 2.0 cm of this probe is an electrically conductive metal (i.e. stainless steel), and the proximal 4.0 cm of the probe is covered with an electrically insulating material to prevent cauterization along the shaft. A grounding pad was applied to the cadaveric liver.
The RF generator supplies voltage at a frequency of 375 kHz. The RF generator display indicates power, tissue impedance, and procedure time. Ablation treatments were sequentially performed in five predetermined locations in each porcine liver. The probe was flushed with normal saline, inserted into the liver under Medicina (Kaunas) 2007; 43 (4) direct vision to a depth of approximately 4-5 cm. After insertion of the electrode, RFA was started. Saline at room temperature was infused at a constant speed of 105 mL/h to maintain low needle-tip temperature. After maintaining the baseline power output at 20 W for 1 min, the output was increased to 50 W, and RFA was continued for 10 min.
Both the temperature and tissue impedance were monitored during ablation. Temperature data were acquired by type K thermocouple and TC-08 thermocouple data logger connected to PC supplied with PicoLog data acquisition software (Pico Technology Limited, Cambridgeshire, UK). The active tips of thermocouples were placed 6 and 11 mm apart from the exposed tip of RFA needle. Temperature data were acquired repeatedly every 1 min.
All the experiments were repeated at least five times with identical results obtained.

Numerical results Investigation of the steady electro-thermal RF ablation process
Here we investigate a balanced heat exchange process, which could be obtained if the electric current flows through tissues for a very long time. Experimentally such values are obtained after about 20 min of RFA, when the measured temperature values do not change in time any more. The values of electric potential, electric field strength, temperatures, and other quantities describing the process are obtained at all nodes of the finite element mesh presented in Fig. 2. The mesh of varying density has been selected in order to achieve a compromise between accuracy of computation and reasonable dimensionality of the model. Dense mesh zone has been generated in the vicinity of the tip of the probe, where steep gradients of field variables have been observed. The solution of 19 700 nodes and 39 400 degree of freedom model is obtained in several minutes on 2 MHz Pentium IV computer. The convergence of the solution has been checked by analyzing several characteristic situations in the mesh presented in Fig. 2 and in the mesh of twice smaller element linear dimension.
Along with contour plots, the graph relationships are presented along four selected lines (Fig. 3): linv -line along electrically active part of the probe; linh1, linh2, linh3 -horizontal radial segments of the axisymmetric domain directed outwards from the surface of the probe.
The electrical conductivity of liver tissue is assumed equal to the electrical conductivity of normality N»0.0105 sodium chloride solution and for different temperature values is obtained from formula [7]. For   Fig. 2. Finite element mesh containing 19 700 nodes, 38 640  example, at T=37 o C electrical conductivity value is s =0.134 1/Wm. The temperature field obtained by taking into account the blood perfusion is presented in Fig. 4A. The results are very close to those presented by Mirza et al. (4). As a characteristic feature of temperature distribution, greater values at the ends of electrically active parts of the electrode may be mentioned. Experimentally such effects may be observed only in vivo, when active blood flow is present.
The results closer to ex vivo experiments are obtained at w = 0, i.e., without blood perfusion term taken into account in [8], Fig. 4B. Theoretically, tissue temperature values are even greater than 100 o C. Practically the zero-perfusion situation is not realistic because of the sodium chloride solution injected through the tip of electrically active part of the probe; therefore, the temperatures remain below 100 o C (Fig.  5). However, the character of temperature distribution in ex vivo experiments is closer to simulation results with strongly diminished influence of perfusion. As characteristic features, the oval contours of temperatures and equal values of temperature along the electrode edge can be mentioned (Fig. 6).
The plots of the electric field strength values along the edge of the electrode are presented in Fig. 7. Satisfactory correlation of the results with similar calculations presented by Mirza et al. (4) is observed. It may be reasonably expected that our results are more accurate, as theoretically the values at sharp edges of the electrode profile should be infinite. In calculations, infinity is not obtained because of finite dimensions of elements. In reality, the values are always finite because of geometrical tolerances of the shape of the electrode.

Results of RFA experiments
After the RFA procedure, two oval-shaped zones of tissue destruction were developed: the inner zone of complete lesion and the external zone of relative lesion. The longitudinal dimension of complete lesion zone was 24±2 mm, and transverse dimension was 12±1 mm. The longitudinal dimension of relative lesion zone was 32±2 mm and transverse dimension -22±1 mm (Fig. 8). The zones were assessed according to the temperature distribution and macroscopic changes in the ablated tissue. Temperature gradually increased for 5 minutes and stabilized thereafter in complete lesion zone, whereas temperature was gradually increasing for 10 minutes in relative lesion zone (Fig. 9,  Fig. 10).

Discussion
The clinical and experimental studies of RFA are the matter of special interest in numerous research centers, as the clinical relevance of the method regarding tumor recurrence is still the matter of debate. Research is directed towards optimization of RFA procedure, applying various mathematical models, estimating factors influencing distribution of temperature in different tissues and environments.
Our experimental data of RFA in ex vivo porcine liver are similar to the published series. Minor variations may be dependent on different RF generators and different types of ablation needles used ( Table 2).
The results reported in the literature exhibit very good correspondence of simulation against the experiment performed on liver-like gels (6). In this study, we developed the experimental validation of mathematical model in ex vivo porcine liver, which resulted in more realistic comprehension of sources of some inadequacies still prevailing in the state-of-the-art models when compared with the behavior of real tissues. However, it should be noted that spread of tissue lesion distally up to 10 mm is not explicable by mathematical model of processes taking places in liver tissue during RFA, developed by our group. Similar observations are published (13); however, this phenomenon is still not explained.
We compared computational heat distribution and occurrence of actual zone of tissue necrosis in ex vivo liver tissue as well as alternation of temperature profile during thermal effect in various ranges from active probe. Obtained absolute and relative zones of tissue destruction, assessed according to the temperature distribution and macroscopic changes in the ablated tissue, are in tolerable correlation with theoretically obtained temperature distribution characteristics (Fig. 11). Existence of such zones is reported and characterized in the literature and is of primary clinical importance (14). Our data showed that the differences of numerical results and experimental measurements were in the range of 5%. Comprehension and definiteness of these zones is necessary for proper planning of clinical treatment. Modeling of the processes in 3D space enables to establish the most appropriate points for insertion of the probe ensuring that all volume of the tumor is inside of the absolute destruction zone.
Certain discrepancies between computational model and experimental data were observed. Increasing distance from active electrode resulted in incrementing difference between the computed and actual temperatures (Fig. 12). The latter observation was neither mentioned nor discussed in previous reports. Assumption may be drawn that experiment in ex vivo tissue at room temperature necessitated higher heat   losses due to relatively cooler liver tissue, encircling the zone of thermal destruction. This might have directly influenced increasing difference of temperatures and production of relative tissue destruction area, were computed temperatures should reach 60-65°C. At temperatures of greater than 60°C, there is denaturation of intracellular protein and lipid bilayers, which lead to irreversible cell death (14). At lower temperatures, enzyme inactivation and mitochondrial injury within the cell occur; however, no physical tissue damage is present (15). Inadequacy between the computational model and experiment in cadaveric liver tissue indicate the necessity to get a more precise insight into the influence of the environment. The coupled electric current and heat exchange model used in this work for representing the RFA processes provides a good comprehension of coupled electrical and thermal phenomena. However, it describes the process as taking place in homogenous environment, simulating characteristics of liver tissue by means of specific constants. Structure of real liver, however, is much more complex and non-homogenous, mainly due to gross architecture of triple vascular network and biliary system. This feature is of utmost importance even in cadaveric liver tissue when cooled-tip ablation electrode with permanent saline infusion is applied. It might be speculated that partial drainage of infused solution into vascular struc- Temperatures, no perfusion tures should distort expected uniform distribution of heat within the liver tissue. Obviously, it is not the case in the tumor tissue, where no major vessels are encountered. However, this feature should be seriously regarded when ablation takes place at the boundary of tumor and surrounding liver tissue, as relative destruction zones are anticipated in this location and may be responsible for local tumor recurrence.
Additional efforts should be made to increase the level of adequacy of the model to reality when mathematically describing the thermal exchange mechanisms. Saline at room temperature was infused at a constant speed of 105 mL/h to maintain low probe-tip temperature during the experiments. This is essential to prevent charring of proximate tissues and to maintain optimal propagation of electric current. However, the infused solution actively participates in heat distribution within the tissue and it was not fully regarded in our computational model. Similar inadequacies have been noted by others.
In the present model, the process of heat off-take is described by using certain equivalent values of thermal conductivity of sodium chloride solution, which simultaneously are applied to quantify the thermal conductivity of the tissue. Moreover, current models em-ploy the generalized perfusion coefficients. Further directions for the development of the computational model should include the modeling of the heat advection by means of the filtrating fluid injected through the channel of the needle electrode.

Conclusions
The principal distribution of the thermal effect, as well as zones of total and relative destruction, as delineated by tissue temperature distribution and actual macroscopic changes, is congruous in the theoretical and the experimental models.
Experimental study revealed some characteristic inadequacies observed by comparing the theoretical against the experimental results -discrepancies of computed and actual temperatures, increasing more distally from active ablation electrode, resulting in generation of relative tissue destruction zone. These differences should be seriously regarded predicting tissue ablation in clinical setting.

Acknowledgement
The Lithuanian Science and Studies Foundation sponsored the study.