Numerical Simulation of Microwave Ablation in the Human Liver

Abstract

Microwave thermal ablation was developed as an alternative to other forms of thermal ablation procedures. The objective of this study is to numerically model a microwave ablation probe operating at the 2.45 GHz level using the finite element and finite volume methods to provide a comprehensive and repeatable study within a human male approximately 25 to 30 years old. The three-dimensional physical model included a human liver along with the surrounding tissues and bones. Three different input powers (10, 20, and 30 watts) were studied, along with the effect of the probe’s internal coolant flow rate. One of the primary results from the numerical simulations was the extent of affected tissue from the microwave probe. The resulting time and temperature results were used to predict tissue damage using an injury integral method. The numerical approach was validated with available experimental data and was found to be within 6% of the average experimentally measured temperatures.

1. Introduction

One of the most common malignancies globally is liver cancer, which results in more than half a million fatalities every year [1]. Possible treatments for liver cancer are surgical operation, cryosurgery, chemotherapy, radiation therapy, and thermal ablation [2]. Thermal ablation has emerged as a prominent minimally invasive modality for treating tumors in the liver and other major organs. The goal of thermal ablation procedures is to raise the temperature of a tumor and a margin of the surrounding healthy tissue to lethal temperatures, which are above 50 °C [3,4]. Radiofrequency (RF) ablation is the most ubiquitous platform available for clinical use; however, it suffers several limitations, including relatively long treatment times and difficulties in creating uniform zones of tissue heating [5]. More recently, microwave thermal ablation has been developed as an alternative to RF ablation systems.

During Microwave ablation (MWA), a coaxial antenna is used to deliver a high frequency (915 MHz or 2.45 GHz) electromagnetic (EM) field into a tumor [3]. When EM energy is applied to the tissue, the molecules with an intrinsic dipole moment (such as water molecules) are forced to continuously realign with the applied field, thus increasing the kinetic energy and elevating local temperatures [6]. Compared to other ablation modalities, microwaves offer the potential for greater volumetric heating, short treatment times, higher intratumoral temperatures, and simultaneous operation of multiple ablation probes [3]. Treatment times for microwave ablation vary from device to device and on the size or number of tumors being treated; however, the tissues are generally treated within 10 min. Because of the invasive nature of the device, computational modeling of microwave ablations has been used as a precise and repeatable technique that can assist with system design, treatment planning, and procedural analysis.

A literature review found studies comparing different antenna frequencies. Curto et al. compared ablation at 9.15 MHz to 2.45 GHz using the finite element method [3]. Their results suggested that faster heating was observed with a single antenna operating at 2.45 GHz. Sun et al. compared water-cooled 915 MHz and 2.45 GHz slot antennas with experiments in ex vivo and in vivo tissues [7]. Their studies suggested that larger ablation zones are achieved with 915 MHz antennas. Finally, Keangin et al. compared single slot and double slot antennas on microwave power absorbed, specific absorption rate, and temperature distribution [2]. Their results suggested that single slot antennas show higher absorption rates and temperature distributions. While these studies provide some insight into finite element modeling techniques and antenna frequency design, they do not provide enough details to replicate the work/results. This is because certain details, such as material properties, mesh/grid information, governing equations, boundary conditions, or dimensions/sizes of the probe or solution domain were not given.

An experimental comparison of RF and MWA devices was performed in [8] to compare the effectiveness of both types of applicators. This was completed on ex vivo porcine livers using internally cooled devices. In each test, two simultaneous devices were powered at 50, 60, 70, and 80 W for 10 min, and the temperature and ablation volume results were observed. The findings of [8] suggest that MWA devices had larger ablation volumes and faster rises in temperature compared to RF devices. The authors conclude that using two MWA devices working simultaneously may be more advantageous than a multipolar RF device in larger liver tumors.

In the recent work of [9], the authors present numerical results (Finite Element Analysis) for four different novel microwave probes (operating at 2.45 GHz) placed in a human liver model. The simulation was two-dimensional and focused on the relative location of the probe tip to a tumor to determine the best position. The authors of [9] applied input powers ranging from 6 to 50 W on tumors ranging from 20 to 55 mm in diameter for a duration of 10 min. Some details, such as the criterion for tissue ablation, convergence criteria for the simulations, grid details, source of the material properties for the tumor, etc., were not provided in the paper. The main result of the work in [9] is guidance on probe positioning related to a tumor’s location and size.

The work presented in [10] presents a low-cost microwave needle for thermal ablation and addresses design and manufacturing aspects. This includes considerations for probe design, material selection, and the fabrication process. While the paper discusses the use of simulations to investigate design aspects (geometry, materials, coolant, etc.), no details (solution domain, governing equations, material properties, boundary conditions, grid details, convergence criteria, etc.) about how the simulations were preformed was provided, only the simulation results. The microwave needle was supplied with 20 W of input power, and it was applied for 200 s. The authors of [10] compared their simulation results with experiments and found good agreement with the predicted and actual ablation volumes.

Several papers in the published literature, such as [11,12], point out the potential importance of tissue shrinkage or expansion during thermal ablation. Specifically, ref. [11] evaluated three different microwave devices and found the volume change results were device-specific. According to the discussion presented in [12], tissue experiencing microwave ablation can contract by 30%. To adequately incorporate the effects of tissue shrinkage or expansion in a numerical model, it would be necessary to have detailed information on how the tissue’s biological processes (such as blood prefusion) are affected by corresponding volume changes. Tissue is a porous medium, and the pore sizes would change with the overall tissue volume. Unfortunately, such information does not appear to be readily available in the published literature and is a potential limitation on numerical studies.

For a more extensive review of both experimental and numerical aspects related to microwave ablation and other thermal therapies, readers are directed to [12,13].

The objective of this study is to fully characterize microwave ablation at the 2.45 GHz level using the finite element and finite volume methods to provide a comprehensive, repeatable study. The three-dimensional anatomically correct numerical model includes a human liver, along with the surrounding biological media (such as bones, muscles, fat) and the microwave probe used to ablate the tissue. Several design parameters will be investigated, and the resulting temperatures and ablation depths will be presented. While several papers in the published literature provide some insight into the modeling techniques, they lack important modeling details presented in this paper.

Overall Solution Domain Geometry

The microwave power will be transferred into the tissue using an 80 mm long commercially available probe comprised of an inner conductor (copper), a dielectric material (Teflon), a polyimide tube, an outer conductor (copper), a stainless steel tube, a ceramic tip, and a flowing coolant fluid (water). The probe is shown schematically in Figure 1.

The probe is positioned in a numerically generated human body created from geometry files for the human anatomy obtained from [14,15,16]. Using the acquired geometry files as a reference, the present authors created an anatomically-correct numerical geometry corresponding to a human male, approximately 25 to 30 years old. The solution domain was limited to a torso segment to include the tissues surrounding the liver (such as bone, muscle, fat, skin) and microwave probe.

The resulting human model used in the present work is shown in Figure 2. The solution area surrounding the human model represented the air region and was created to give ample distance between the human model and the boundary conditions. Specific details about the boundary conditions for the human model will be discussed in the forthcoming sections.

2. Numerical Methods

The forthcoming discussion is focused on the governing equations, which will be solved for the three-dimensional physical solution domain presented in the previous section. Two separate grids (meshes) will be utilized depending on the set of equations to be solved. The first set of equations to be presented, Equations (1)–(6), will be applied using a finite volume grid (mesh) corresponding to the CFD (Computational Fluid Domain) region. The remaining equations, Equations (7)–(11), are applied to a separate finite element grid. Both grids encompass the same physical geometry, and details about both grids will be given in Section 2.1.

The Reynolds-averaged Navier–Stokes equations (RANS equations), represented in Equation (1), are time-averaged equations of motion for fluid flow. The idea behind these equations was first originated by Osborne Reynolds [17]. The energy equations, Equations (3) and (4), are derived from energy conservation and Fourier’s law [18]. In the case of Equation (4), the generation terms or the source terms are separated into three distinct contributions, namely the metabolic heat generation, blood prefusion, and the energy absorbed by the microwave field. In terms of Maxwell’s equations, Equations (8)–(12), several individuals were instrumental in their formulation. However, Heinrich Hertz is credited for proving the equations experimentally, which eventually led to their widespread use [19].

The fluid flow outside of the human model is expected to be laminar, and the fluid motion is caused by natural convection due to the temperature difference between the human skin and the surrounding air environment. In Einstein notation, the governing equations for unsteady, three-dimensional, constant property laminar flow are

$$\rho \frac{\partial u_j}{\partial t}+\rho \left(u_i\frac{\partial u_j}{\partial x_i}\right)=-\frac{\partial p}{\partial x_j}+\mu \frac{{\partial }^2{u}_j}{\partial x_i^2}+\rho g_j$$

(1)

where $i=1,2,3; j=1,2,3$ and $i\ne j$

$$\frac{\partial \rho }{\partial t}+\frac{\partial u_j}{\partial x_i}=0$$

(2)

where $\rho$ is the fluid density, $u$ is a velocity component, $t$ is time, $x$ is a direction coordinate, $p$ is pressure, $\mu$ is the dynamic viscosity, $g$ is the standard earth gravity, and i and j represent the three coordinate directions.

The energy equation for the fluid domain is

$$\rho c_p\frac{\partial T}{\partial t}+\rho c_p\frac{\partial \left(u_{i}T\right)}{\partial x_i}=k\frac{{\partial }^{2}T}{\partial x_i^2} i=1,2,3;j=1,2,3 \mathrm{and} i\ne j$$

(3)

In this equation, T is the fluid temperature, c_p is the specific heat at constant pressure, $\rho$ is the fluid density, and k is the thermal conductivity.

The energy equation for describing heat transfer inside of the solid domains is

$$\rho c\frac{\partial T}{\partial t}=k\left(\frac{{\partial }^{2}T}{\partial x_i^2}\right)+{\dot{S}}_{\mathrm{met}}+{\dot{S}}_{\mathrm{perfusion}}+{\dot{S}}_{\mathrm{microwave}} i=1,2,3$$

(4)

where ${\dot{S}}_{\mathrm{met}}$ is the metabolic heat generation inside of the human tissue, ${\dot{S}}_{\mathrm{perfusion}}$ represents the contribution of blood perfusion to energy transfer within the tissue, and ${\dot{S}}_{\mathrm{microwave}}$ is the energy absorbed from the microwave probe.

The numerical values of ${\dot{S}}_{\mathrm{met}}$ are prescribed based on location in the human body whereas ${\dot{S}}_{\mathrm{perfusion}}$ comes from the well-known approach by Pennes [20]. The perfusion source term takes the form

$${\left(\rho c_{\mathrm{p}}\right)}_{\mathrm{b}}\psi \left(1-\zeta \right)\left(T_{\mathrm{b}}-T_{\mathrm{t}}\right)$$

(5)

where the subscripts t and b refer to local tissue and blood, respectively, $\psi$ is the tissue-specific perfusion rate, and $\zeta$ is an equilibrium constant between capillary blood and the tissue, which is assumed to be uniform throughout the tissue where $0\le \zeta \le 1$. Based on the work of [21], $\zeta$ was chosen to be equal to 0. The blood properties for Equation (5) were assigned a density of 1204 kg/m³ and a specific heat of 2672 J/kg$·$°C. In the present model, the temperature $T_{\mathrm{b}}$ is the mass-based average temperature of the blood within a specific region of the body, and $T_{\mathrm{t}}$ is the local tissue temperature anywhere within the body.

The source term ${\dot{S}}_{\mathrm{microwave}}$ is computed from the electromagnetic field distribution in the solids and is defined as

$${\dot{S}}_{\mathrm{microwave}}=\frac{\sigma {\left|E\right|}^2}{2}$$

(6)

where $\left|E\right|$ represents the computed electric field strength (V/m) and $\sigma$ is the electric conductivity (S/m) of the material. If $\sigma$ = 0, the material is dielectric (electrically insulating material).

When divided by the tissue density $\rho$, ${\dot{S}}_{\mathrm{microwave}}$ is referred to as the Specific Absorption Rate (SAR), which represents the electromagnetic power (watts) deposited per unit mass (kg) by the microwave energy

$$\mathrm{SAR}=\frac{\sigma {\left|E\right|}^2}{2\rho }$$

(7)

The propagation of the microwaves is governed by Maxwell’s equations expressed as

$$\frac{\partial E_{\mathrm{j}}}{\partial x_{\mathrm{i}}}-\frac{\partial E_{\mathrm{i}}}{\partial x_{\mathrm{j}}}=-\frac{\partial B_{\mathrm{k}}}{\partial t} i=1,2,3; j=1,2,3;k=1,2,3 \mathrm{and} i\ne j\ne k$$

(8)

$$\frac{\partial H_{\mathrm{j}}}{\partial x_{\mathrm{i}}}-\frac{\partial H_{\mathrm{i}}}{\partial x_{\mathrm{j}}}=J_k+\frac{\partial D_{\mathrm{k}}}{\partial t}$$

(9)

$$\frac{\partial D_{\mathrm{i}}}{\partial x_i}={\rho }_{\mathrm{o}}$$

(10)

$$\frac{\partial B_{\mathrm{i}}}{\partial x_{\mathrm{i}}}=0$$

(11)

where $E$ is the electric field strength (V/m), $B$ is the magnetic field (T), $H$ is the magnetic field strength (A/m), $J$ is the electric current density (A/m²), $D$ is the electric flux density (C/m²), and ${\rho }_{\mathrm{o}}$ is the volume charge density (C/m³). The constitutive relationships are

$$\begin{array}{c}J=\sigma E \\ B={\mu }_{\mathrm{o}}{\mu }_{r}H \\ D={\epsilon }_{\mathrm{o}}{\epsilon }_{r}E\end{array}$$

(12)

2.1. Computational Grid and Convergence

The grids were generated using commercially available software ANSYS, and the numerical solutions were generated using solvers from ANSYS.

For the finite volume mesh (CFD), the implicit solver uses a second-order-accurate discretization for the advection terms in the RANS equations, which were solved iteratively using a false-transient approach. A mesh independency study using several different computational meshes was employed by systematically varying the node count (a node is the location where the calculations are performed). For the CFD mesh independence study, it was found that increasing the node number from 4.3 million to 52 million nodes gave rise to temperature variations no greater than 0.5%. A representative image of the mesh used for the CFD simulations in the neighborhood of the microwave probe is shown in Figure 3. The mesh, as seen in Figure 3, is comprised of tetrahedral elements. Each tetrahedral element is vertex-centered and has six integration points. The results reported here correspond to a mesh of 24 million nodes.

The CFD solutions were considered to be converged to sufficient accuracy when the root-mean-square (RMS) residuals for all governing equations were 10⁻⁶ or less. The simulations were performed as a transient using time steps between 0.001 and 0.0001 s. The total computational time varied based on node count and simulation parameters. The average computational time was between 15 to 20 days per simulation.

A separate mesh and convergence criteria for the HFSS (High Frequency Simulation Software) simulations were needed. The HFSS simulations used an adaptative tetrahedral mesh and multiple passes to refine the mesh to reach a target maximum Delta S criterion (convergence criterion for HFSS simulations). Delta S represents the change or difference in the S parameters (reflection and transmission coefficients) as the grid is refined and rerun. The maximum Delta S criteria varied between 0.02 to 0.002 with a maximum of 10 passes to reach convergence. In other words, the simulation would run and refine the grid 10 times (or less) until the change in S parameters from one grid to another was less than 2% or 0.2%. Each adaptive mesh pass added 100,000 to 200,000 elements in the region around the microwave probe, and the resulting mesh was a mixture of first and second order elements. The SAR value, which represents the electromagnetic power deposited per unit mass by the microwave energy, did not change or affect the temperature results when the HFSS mesh had greater than 2 million elements.

2.2. Boundary Conditions

For the air region boundary conditions, all the surfaces far from the human/biological domain were specified as openings (entrainment) at a gauge pressure of 0 Pa and an ambient temperature of 25 °C. These ‘weak’ boundary conditions allowed the air to enter or leave the solution domain at any velocity or direction. The remaining air region boundaries were specified as symmetry boundaries. For the internal human boundaries, the surfaces were specified as a constant temperature wall at 37 °C.

For the microwave probe, a Wave Port boundary condition was applied on the far end of the device (extending outside of the tissue and into the air domain, as seen in Figure 2a, 80 mm from the device’s tip) on the inner and outer conductor. Depending on the case being analyzed, a specific constant wattage was applied as the energy source. On the far-field boundaries, an EM radiation boundary was applied.

The microwave probe has an internal coolant (water) to prevent the temperatures near the device from getting too hot too quickly. Biological tissues that get too hot (temperatures above 100 °C) can lead to tissue charring, which occurs near the near neighborhood of the device [11]. Without active cooling, microwave ablations devices are known to reach high temperatures, up to 120 °C [12]. The charring makes the positioning and/or the removal of the probe challenging and therefore undesirable. The optimal coolant flowrate to maintain reasonable temperatures was not known apriori and therefore was made a design parameter. The mass flow rate of the coolant was parametrically varied, from 0.5 g/s to 2 g/s, and entered the microwave probe through an annulus at 25 °C. After passing through the microwave probe and turning 180 degrees, the coolant exited the probe.

2.3. Material Properties

It is necessary to define both the thermophysical properties (specific heat c_p, density $\rho$, thermal conductivity k, and dynamic viscosity $\mu$) and the electromagnetic properties (relative permittivity ${\epsilon }_{\mathrm{r}}$ and relative permeability ${\mu }_{\mathrm{r}}$) of each of the materials. For the initial set of numerical cases, all of the material properties will be specified as homogeneous and isotropic. For reference, the permittivity of free space ${\epsilon }_{\mathrm{o}}$ $\cong$ 8.8542 × 10⁻¹² F/m and the permeability of free space ${\mu }_{\mathrm{o}}$ $\cong$ 1.2566 × 10⁻⁶ H/m. It was necessary to obtain the material properties from several different references [22,23,24,25,26,27,28] since no single reference contained all of the required properties.

Table 1. Thermophysical properties for the solid materials [2,4,22,23,24,25].

		Density $\mathit{\rho }$ (kg/m3)	Specific Heat ${\mathit{c}}_{\mathbf{p}} (\mathbf{J}/\mathbf{kg}\mathit{·}$°C)	$\mathbf{Thermal} \mathbf{Conductivity} \mathit{k} (\mathbf{W}/\mathbf{m}\mathit{·}$°C)
Tissues	Liver	1100	3350	0.55
	Muscle	1050	3770	0.5
	Fat	850	2500	0.2
	Skin	1000	3770	0.21
	Bone	1300	1590	1.16
Probe	Ceramic	6050	400	2.2
	Copper	8960	386	401
	Stainless steel	8055	500	16.2
	Polyimide	1400	1320	0.471
	Teflon	800	1000	0.25

lists thermophysical properties for all of the solid materials, and

Table 2. Electromagnetic properties for a frequency of 2.45 GHz [26,27,28].

		$\mathbf{Relative} \mathbf{Permittivity} {\mathit{\epsilon }}_{\mathit{r}}$	$\mathbf{Relative} \mathbf{Permeability} {\mathit{\mu }}_{\mathit{r}}$	Electric Conductivity σ (S/m)
Tissues	Liver	44	1	1.79
	Muscle	49.6	1	2.56
	Fat	12	1	0.82
	Skin	44	1	1.85
	Bone	4.8	1	0.21
Probe	Ceramic	29	1	0
	Copper	1	0.999	5.8 × 107
	Stainless steel	1	1	1.1 × 106
	Polyimide	3.5	1	0
	Teflon	1.7	1	0

lists the electromagnetic properties for a frequency of 2.45 GHz used for the present simulations.

The air density obeyed the Ideal Gas Law where the ambient temperature and pressure were 25 °C and 1 atm, respectively. The air had a constant specific heat of 1000 J/kg$·$°C, a thermal conductivity of 0.0261 W/m$·$°C, and a dynamic viscosity of 1.835 × 10⁻⁵ Pa$·$s. The coolant inside of the probe was water and had a density of 998.2 kg/m³, a specific heat of 4182 J/kg$·$°C, a thermal conductivity of 0.6 W/m$·$°C, and a dynamic viscosity of 0.89 mPa$·$s.

2.4. Initial Conditions

Since the simulations were time dependent, it was necessary to specify the starting conditions (i.e., temperatures, velocities, and pressures) for all of the biological tissues, the probe, and the surrounding air environment. A stepwise approach was adopted to obtain more realistic starting conditions. In the first step, all biological tissues were initialized at 37 °C, and the surrounding air was initialized at 25 °C. For this first step, the model was run to simulate 30 min to develop an appropriate temperature gradient in the human tissues and allow natural convection to occur. The results of the first step were then used to initialize the second step, where the probe, at 25 °C, was placed inside the human liver and subsequently activated after 5 s.

3. Results and Discussion

Several different cases were identified to study the effect of the parameters, such as input power and coolant flow rate, on tissue damage depth and time. These investigated cases are listed in

Table 3. Parametric variation of input parameters.

Case #	Input Power ${\mathit{P}}_{\mathbf{input}} \left(\mathbf{W}\right)$	Coolant Flow Rate $\dot{\mathit{m}}$ (g/s)
1	10	1
2	20	1
3	30	0
4	30	0.5
5	30	1
6	30	2

. Each of the listed cases were simulated for 5 min with the microwave probe active. Due to losses in the probe, only a portion of the input power was converted into thermal energy absorbed by the liver. The remaining amounts of the thermal energy not absorbed directly by the liver were absorbed by surrounding tissues or parts of the probe. Case #3 was included, as a reference, to observe the impact of not utilizing a coolant. It is expected that the resulting temperatures without a coolant will exceed the desired temperature range.

Table 4. Distribution of absorbed thermal energy.

Input Power ${\mathit{P}}_{\mathbf{input}} \left(\mathbf{W}\right)$	Total Input Power after Losses ${\mathit{P}}_{\mathbf{input}-\mathbf{losses}} \left(\mathbf{W}\right)$	Power Absorbed by Liver ${\mathit{P}}_{\mathbf{absorbed},\mathbf{liver}} \left(\mathbf{W}\right)$	Power Absorbed by Surrounding Tissues ${\mathit{P}}_{\mathbf{absorbed},\mathbf{surrounding}} \left(\mathbf{W}\right)$
10	7.15	5.96	0.92
20	14.30	11.91	1.83
30	21.44	17.85	2.75

reports the resulting distributions of energy absorbed directly by the liver and the surrounding tissues (muscle and fat). The results show that ~29% of the input energy is lost due to the device, and the remaining ~71% is converted into thermal energy for heating. After losses, ~83% of the energy is absorbed directly into the liver tissue, ~13% entered the tissues surrounding the liver, and the remaining ~4% was absorbed by the probe’s components.

3.1. Qualitative Temperature Results

The human numerical model created for this work was three-dimensional, and it is inherently difficult to present the forthcoming results in three dimensions. Therefore, to simplify the presentation, a single plane of observation that bisects the microwave probe in the axial direction will be used. This plane’s location relative to the human model is illustrated in Figure 4a,b. Figure 4b displays a color contour diagram of temperature on the selected plane. Since the probe is axisymmetric in design, the microwave field and the heating of the liver tissue occur in an axisymmetric shape. An illustration of the axisymmetric shape is illustrated by a typical SAR pattern, shown in Figure 5 for Case 5 (30 W input power and coolant flow rate of 1 g/s). Figure 5 shows two color contour diagrams, (a) and (b), for Case 5 with different legends for clarity. As seen in this figure, the SAR pattern is axisymmetric around the microwave probe. Therefore, the subsequent temperature results (such as in Figure 6) shown on the plane represent the thermal results in any radial cross-section because of the symmetric shape of the microwave field.

Figure 6 displays the characteristic heating pattern of the tissue, as a function of time, due to the microwave probe. The results correspond to an input power of 30 W and a coolant flow rate of 1 g/s. At the start of the procedure, t = 0 s, the probe is at 25 °C, and the surrounding liver tissue is 37 °C. As seen in Figure 6, once the microwave field is activated, the liver tissue quickly increases in temperature, and the heated region of the tissue forms a spherical shape. The highest temperatures occur 15 mm from the probe tip in the axial direction and 2 mm in the radial direction from the outside of the probe. As seen in the color contour diagrams, the coolant fluid circulating through the probe can maintain reasonable temperatures (less than 55 °C) at the interface between the probe and the tissue to prevent undesirable charring of the tissue.

3.2. Quantitative Results

A useful result is quantifying tissue damage caused by heat generated from the microwave EM field. A commonly used metric for cell injury/damage due to heating is the injury integral method [29,30,31]. This method, also known as the Arrhenius damage integral, uses a first-order chemical reaction approach to quantify tissue damage by protein denaturation and has the form

$$\mathsf{\Omega }\left(t\right)=\ln \left(\frac{C\left(0\right)}{C\left(t\right)}\right)=\mathrm{A}\underset{0}{\overset{t}{\int }}{\mathrm{e}}^{\left[-{\mathrm{E}}_{\mathrm{a}}/\mathrm{R}T\left(t\right)\right]}\mathrm{d}t$$

(13)

where $\mathsf{\Omega }\left(t\right)$ is the damage index as a function of time, $C\left(0\right)$ is the number of undamaged cells at time zero, $C\left(t\right)$ is the number of undamaged cells at time $t$, $\mathrm{A}$ is the frequency factor of liver tissue (8 $\times 10$³⁹ s⁻¹ [32]), ${\mathrm{E}}_{\mathrm{a}}$ is the activation energy of the liver tissue (257,700 J/mol [32]), $\mathrm{R}$ is the universal gas constant (8.3143 J/mol$·$K), and $T\left(t\right)$ is the tissue temperature (in Kelvin) at any location and time. For reference, a value of $\mathsf{\Omega }\left(t\right)$ = 0 signifies the tissue is undamaged, $\mathsf{\Omega }\left(t\right)$ = 0.7 indicates 50% of the cells would be still viable, and $\mathsf{\Omega }\left(t\right)$ = 1 corresponds to 37% of the cells being still viable. Within the published literature, different criteria for $\mathsf{\Omega }\left(t\right)$ have been chosen for thermal injury. As discussed in [33], values of $\mathsf{\Omega }\left(t\right)$ = 1 or 0.53 are commonly chosen in the published literature. For instance, ref. [34] choose $\mathsf{\Omega }\left(t\right)$ = 1 because the author believes that is the point when thermal necrosis occurs. However, ref. [33] suggests a threshold value of $\mathsf{\Omega }\left(t\right)$ should be when the first observable change occurs to the tissue and not complete cell death. Therefore, a value of $\mathsf{\Omega }\left(t\right)$ = 0.1 might be a better criterion because it could be the point when tissue damage would likely be detectable [33]. An alternative and potentially advantageous presentation of $\mathsf{\Omega }\left(t\right)$ is in terms of the percent of surviving cells [31], which is defined as

$$\mathrm{Percent} \mathrm{of} \mathrm{surviving} \mathrm{cells}=\frac{1}{{\mathrm{e}}^{\mathsf{\Omega }\left(t\right)}}\times 100$$

(14)

This is a useful presentation of the predicted tissue damage from the injury integral since its values are bound between 0 and 100%, and it is better at revealing what is happening in the tissue compared with values of $\mathsf{\Omega }\left(t\right)$. The alternative to the percent of surviving cells is the damage fraction or necrotic fraction, equal to 100—percent of surviving cells.

To present results using Equation (14) in a meaningful manner, several points relative to the probe’s location in the tissue have been identified in Figure 7. Each of the points shown in Figure 7 have been numbered, and the temperature variation at each point as a function of time was recorded. Figure 8 has been prepared to show the temperature variation as a function of time for several of the points identified in Figure 7 for Case 5 (30 W input power and a coolant mass flow rate of 1 g/s).

As seen in Figure 8, the probe tip (point 1) reaches a maximum temperature of 62 °C over the five-minute interval, and 5 mm from the tip in the axial direction (point 3) does not exceed 50 °C. At the probe surface away from the tip, point 4 shows an initial cooling of the tissue due to the circulating coolant, though it eventually reaches a maximum temperature of 41 °C. The highest temperatures are located in the neighborhood of point 16, which reaches a temperature of 129 °C at five minutes. At point 20, 8 mm from point 16 in the radial direction, the temperature reaches 75 °C. Overall, the figure shows the temperatures within a region located 20 mm radially from the probe surface are sufficiently high to cause tissue damage.

With the time and temperature data, Equation (14) can plot the percentage of undamaged cells at any location. The percentage of undamaged cells corresponding to the locations shown in Figure 8 have been plotted in Figure 9 for Case 5. Not surprisingly, the survival of cells in the regions near points 16 and 5 drops very quickly due to the elevated temperatures at these locations. After 105 s, the liver tissue within 10 mm radially from the probe surface reached a value of $\mathsf{\Omega }\left(t\right)$ = 1.0 (37% of the cells still viable or 63% damaged) except for some locations in direct contact with the probe. At point 4, because of the coolant flow in the probe, 83% of the cells are still surviving ($\mathsf{\Omega }\left(t\right)$ = 0.185 or 17% damaged) after five minutes. Since the microwave field is not very strong beyond the probe tip, point 3 reaches a 30% survival rate ($\mathsf{\Omega }\left(t\right)$ = 1.2 or 70% damaged) at 5 min.

The most influential parameter that was varied was the input power to the microwave probe. In general, it was expected that higher power levels would increase temperatures and reduce the overall time to necrose the liver tissue. Figure 10 has been prepared to show the percentage of undamaged cells for several key locations: at the tip (point 1), at 1 cm deep (point 20), and 2 cm deep (point 25) for the three input powers of 10, 20, and 30 W (Cases 1, 2, and 5). As expected, the 30 W input power necrosed the tissue faster than the 20 and 10 W input power.

Comparing the results for the 20 and 30 W cases at a depth of 1 cm reveal that the time to reach $\mathsf{\Omega }\left(t\right)$ = 1.0 (37% of the cells still viable or 63% damaged) was different by 35 s. However, at 10 W, the time to reach $\mathsf{\Omega }\left(t\right)$ = 1.0 at a depth of 1 cm took an extra 135 s. At the probe tip, the time to reach $\mathsf{\Omega }\left(t\right)$ = 1.0 for the 20 W case required an extra 70 s compared to the 30 W input power, and at 10 W, the tip did not get hot enough in five minutes to reach $\mathsf{\Omega }\left(t\right)$ = 1.0. At a depth of 2 cm, the percentage of surviving cells remained above 60% (40% damaged) over the five-minute time period.

Overall, the time to reach $\mathsf{\Omega }\left(t\right)$ = 1.0 varied based on location and input power. Of the three locations displayed in Figure 9, the probe tip is strongly influenced by the circulating coolant. For the 30 W case, the SAR field at the tip was 200 W/kg compared to the maximum values occurring near point 15, SAR > 100,000 W/kg, as shown in Figure 5. The SAR field created by the probe is relatively strong (~500 W/kg) at a depth of 1 cm in the tissue causing the time to reach $\mathsf{\Omega }\left(t\right)$ = 1.0 for the 20 and 30 W cases to be within 35 s of each other. At a depth of 2 cm from the probe, the SAR field is much weaker (~40 W/kg), and heat conduction through the tissue plays a more important role in the resulting tissue damage.

Based on comparing the three input powers, 10, 20, and 30 W, the 30 W input power is the most effective at damaging the tissue cells in the shortest amount of time. In order to effectively damage the tissue 2 cm radially from the probe surface, the duration of the microwave field at 30 W would need to be increased to 480 s (8 min).

The next set of results to be presented relates to the effect of the coolant flow rate. The purpose of the coolant is to maintain reasonable probe temperatures to prevent the charring of the tissue in contact with the probe. Having too low a flow rate could lead to tissue charring, whereas too high a flow rate could diminish the effectiveness of the intended therapy. Several cases (Cases 3, 4, 5, 6) were examined where the coolant flowrate was parametrically varied from 0 to 2 g/s. The case without a coolant flowrate is a limiting case where the temperatures are expected to be too high. For the non-coolant case, the probe was filled with stagnant air instead of water. The temperatures at three locations along the probe surface, points 1, 4, and 15, are presented in Figure 11 to show the effect of varying the coolant flow rate.

As seen in Figure 11, there is an obvious effect of the internal coolant flowrate. For the case without a flowing coolant (Case 3), the tissue temperatures elevate much faster (points 4 and 15) along the probe surface. However, the temperature results corresponding to coolant flowrates of 1 and 2 g/s are within 1% or less with each other. It is interesting to note that at point 4, there is an initial decrease in temperature (due to the coolant) for the first 6 s. Point 15 is located closer to the region of maximum tissue temperature and does not experience an initial temperature decrease.

3.3. Comparison with Experimental Results

To validate the present modeling approach, a specific simulation was performed to compare with published experimental data presented. Since experimental data for microwave probe tissue ablation performed in vivo within a human liver was not readily available, temperature data reported in [35] was chosen for comparison with the present simulations. The work performed in [35] was done with bovine liver. In general, bovine liver is used as a surrogate for human liver because of the similar physical properties. For example, the permittivity of human liver is ~43 at 2.45 GHz [36] compared with bovine liver at ~44 at 2.45 GHz [35] in the same temperature range.

The experimental work performed in [35] used a microwave probe operating at 2.45 GHz and 30 W, similar to the present study. The specific probe used in [35] was of a different design than the present study, but it utilized an internal cooling system fed by a constant temperature bath at 10 °C. The specific coolant flowrate and physical design of the microwave probe used in [35] were not discussed. The experimental microwave ablation was performed ex vivo in a bovine liver held within a test fixture. The bovine liver started at 15 °C, and the microwave energy was applied for 10 min.

A simulation was performed with the present microwave probe using 30 W of input power and a coolant flow rate of 1 g/s. To match the experiment, the simulation started with the liver tissue at 15 °C, and the microwave energy was applied for 10 min. Since the experimental work of [35] did not use living tissue, the present simulations excluded both the metabolic heat generation and perfusion from the energy equation. The simulation boundary conditions of the tissue sample were assumed to be adiabatic since (a) the tissue sample was large compared to the area affected by the microwave field, and (b) the actual conditions surrounding the sample that was used in [35] were unknown.

The experimental temperature measurements reported in [35] were obtained with a temperature probe located 10 mm axially from the probe tip and 5 mm radially outward from the microwave probe surface. As reported in [35], this specific location corresponded to the maximum power deposited by the microwave antenna. In the present numerical simulation, the maximum power deposited by the microwave antenna occurs approximately 14 to 15 mm from the tip and 2 to 3 mm radially from the microwave probe surface. For comparison purposes, temperatures from the simulation were extracted from the zone of maximum deposited power and plotted in Figure 12. The corresponding experimental temperature results presented in [35], with a 95% confidence interval, are also reproduced in Figure 12.

Examination of Figure 12 shows very good agreement between the simulated reproduction of the work performed in [35] and the temperature results from the experimental measurements. In general, the temperatures predicted from the simulation were generally within 6% of the average experimentally measured temperatures. Slight differences in the temperature results are expected because (a) the probes are of different designs, (b) the tissue used in the experiments will not be as homogenous as in the simulated tissue, (c) it is expected there would be heat loss to the environment during the experiments and this was not modeled in the simulation, and (d) experimentally locating a temperature probe at the point of the maximum microwave field could influence the measured temperature (because of the probe’s thermal heat capacity, electromagnetic material properties, response time, etc.) whereas the temperature probe was not modeled in the simulation.

4. Conclusions

Microwave thermal ablation procedures have become a popular alternative to other existing ablation procedures (such as RF ablation) for dysfunctional tissue treatment. While several papers in the published literature provide some insight into finite element modeling techniques and antenna frequency design, they lack many important modeling details. The objective of this study was to fully characterize microwave ablation at the 2.45 GHz level using the finite element and finite volume methods to provide a comprehensive, repeatable study.

The three-dimensional numerical model included a physically realistic human liver (i.e., realistic dimensions, physical properties, biological processes, etc.) along with the surrounding tissues and bones found in the human body. Three different input powers (10, 20, and 30 watts) were studied, along with the effect of the probe’s internal coolant flowrate. One of the primary results from the numerical simulations was the extent of affected tissue (elevated temperatures) from the microwave probe.

Not unexpectedly, the higher power inputs led to higher tissue temperatures in a shorter duration of time. The resulting time and temperature results were then used to predict tissue damage using an injury integral method discussed in Section 3.2. For instance, with an input power of 10 W, it took more than twice as much time to kill 50% of the tissue cells within a 1 cm radius compared to 30 W (less than 100 s).

The numerical approach was then validated with available experimental data from the literature. Since human-specific experimental data was unavailable, it was necessary to use experimental data generated with ex vivo bovine liver. Bovine liver is often used as a surrogate for human liver because of its similar material properties. However, since the experiments were performed ex vivo, the tissue was not living and did not experience biological processes, such as blood perfusion or metabolic heat generation. The numerical simulation results were compared to the experimentally measured temperatures and were found to be generally within 6%. While slight differences in the temperatures predicted by the numerical simulation and the results reported in the experiment can be expected, the overall agreement is considered satisfactory.