Temperature simulation of microwave ablation based on improved specific absorption rate method compared to phantom measurements

Abstract

Purpose: To propose an improved specific absorption rate (SAR) computation method by incorporating a thermal conductivity item, develop mathematical simulation models based on the SAR, and enhance the accuracy of temperature simulation for 2450-MHz microwave fields. Materials and methods: Experimental data of 2450-MHz microwave antenna were obtained with a phantom and then new SAR equations were calculated according to Pennes bio-heat transfer equation. Next, the SAR equations were regarded as heating sources and SAR-derived temperature changes were simulated numerically; Finally, temperature measurement results were compared with simulation data to validate the accuracy of the new SAR method. Results: The simulated and measured temperature changes were generally in good agreement for the phantom. According to comparison results, the maximum error was less than 6 °C, the average error was less than 2 °C, and the standard deviation was less than 1 °C. Conclusions: This new SAR-derived simulation method can significantly simplify the simulation process and improve the prediction accuracy of microwave ablation temperature field.

Keywords:

• Microwave ablation
• specific absorptionrate
• temperature distribution simulation

1. Introduction

Microwave ablation (MWA) has emerged as a well-targeted strategy for tumor treatment and has made marked advancements in recent years. Its principle is to heat tumor tissues by percutaneously inserted microwave (MW) antenna and induce cell protein coagulation or even denaturation, thereby killing tumor cells and achieving treatment purposes [1]. Compared with traditional therapeutic methods, MWA has fewer complications and minimally invasive advantages [2–4]. Although MWA has capacious application prospects, the temperature field is still difficult to predict and monitor precisely. Therefore, how to control ablation regions has become a long recognized difficulty and continued to challenge scientific researchers and clinicians [5,6].

To this end, numerous researches have paid considerable attention to temperature field distributions and ablation zone sizes of MW antenna. A wide variety of ablation simulation methods based on the bio-heat transfer equation have been proposed [7,8]. Practical coagulative outcomes of MWA are expected to be obtained with physical experiments and computer-simulation methods. In this way, patient-specific therapy of tumors may be achieved by further referring to the patient's medical images [9]. The most critical parameter affecting temperature field distribution is the specific absorption rate (SAR). Due to constant power during MW treatment, the SAR distributions could be obtained by experiment methods and applied to different tissues by modifying thermal parameter values in simulation process. However, the simulation process of SAR tends to ignore the thermal conductivity item, and the time selected for computing temperature-rising rate is not very reasonable. While SAR modeling is maturing rapidly, the accuracy of temperature predictions is lagging behind [10]. Thus, the process of temperature field distribution simulation is still to be improved [11,12].

In this study, we focused on the precise computation of SAR and the simulation of temperature field based on the computed SAR. SAR computation and temperature modeling were separated to reduce the computational complexity. We firstly performed tissue phantom experiments to obtain actual temperature data, and then employed an improved method to derive SAR equations of MW antenna. Next, finite element analysis method was used to simulate MWA temperature field distribution, where the SAR was regarded as the source item of heating. Finally, simulation results were compared with measurement values to determine the feasibility of the SAR method. In brief, this research could be used in a situation in which the heating power was constant. This SAR-derived method could provide tissue-specific simulation, simplified modeling, ablative margin and far field temperature prediction, and feedback strategy application.

2. Materials and methods

2.1 Experimental setup

The experimental system consisted of a water-cooling MWA device (KY-2000; Kangyou Microwave Energy Sources Institute, Nanjing, China), a data acquisition device (34970 A; Agilent Technologies Inc., Santa Clara, CA, USA), some copper-constantan thermocouple thermometers, and a MW phantom (Figure 1).

Figure 1. Experimental setup of MWA.

The MWA device had a MW frequency of 2450 MHz and an output power of 5–120 W. It was equipped with coaxial cables of low energy consumption and was attached to hard MW antenna by these cables. The MWA device was also equipped with a peristaltic pump for water cooling. The MW antenna had a diameter of 1.9 mm and a length of 180 mm. The antenna employed aperture emission, and the aperture was 1.0 mm in width and 11 mm from the tip. The copper-constantan thermocouple thermometers had a diameter of 0.003 inch and could measure the temperature changes in substantial real-time. The MW phantom was prepared according to Jiang et al. [13] and was composed of 3.5 weight percent sodium carboxymethyl cellulose, 27.15 weight percent polyvinyl chloride, 0.35 weight percent sodium chloride, and 69 weight percent distilled water.

2.2 Acquisition of actual temperature data

MW antenna 3 D temperature field was axisymmetric around the antenna in a homogeneous medium, but was asymmetric between the aperture’s front side (z > 0) and rear side (z < 0) due to water-cooling effects (Figure 2). The ordinate axis z was oriented to coincide with the antenna direction; the abscissa axis R was perpendicular to the antenna direction; and the front end of the antenna aperture was set as the origin (Figure 2). Consequently, 3 D spatial modeling of the MW antenna’ SAR distribution could be actually transformed into 2 D planar modeling of SAR distribution along R > 0 direction in the RZ plane. Based on MW antenna characteristics and phantom properties, 22 thermocouples were arranged in the RZ plane to characterize the temperature field distribution. Table 1 showed coordinates of the thermocouples.

Figure 2. Spatial distribution view of 22 thermocouples.

Table 1. Coordinates of 22 thermocouples.

Rear side (r,z) (cm)		Front side (r,z) (cm)
No.	Coordinate	No.	Coordinate
1	(1.0, −4.0)	13	(1.0, 0.5)
2	(1.0, −3.5)	14	(1.0, 1.5)
3	(1.0, −3.0)	15	(1.0, 2.5)
4	(1.0, −2.0)	16	(1.0, 3.0)
5	(1.0, −1.0)	17	(1.5, 0.5)
6	(1.5, −3.0)	18	(1.5, 1.5)
7	(1.5, −2.0)	19	(1.5, 2.5)
8	(1.5, −1.0)	20	(2.0, 0.5)
9	(2.0, −2.0)	21	(2.0, 1.0)
10	(2.0, −1.5)	22	(2.0, 1.5)
11	(2.0, −1.0)
12	(2.5, −1.0)

Temperature measurements were made in the MW phantom. The MW antenna was horizontally passed through antenna hole of antenna fixing plate, and the thermocouples were vertically passed through corresponding thermometer holes of thermometer fixing plate (Figure 2). The MW antenna and thermocouples were disposed at 90 degree in the phantom, and the insertion depth of the antenna and thermocouples depended on the fixing plates. The power and heating duration of the MWA device were set at 40-60 W (40 W, 45 W, 50 W, 55 W, 60 W) and 600 s respectively. The sampling interval of the data acquisition device was 2 s. During MW heating, the data acquisition device recorded the temperature data of each thermocouple in real time and plotted corresponding temperature curves. After the heating process, the collected temperature data were stored in a text file for subsequent processing.

2.3 Computation of SAR distribution

SAR represents MW energy acting on one point per unit time, which can lead to temperature rising and thermal conduction of the point. A new computing method of SAR was proposed in this study. Thermal diffusion and conduction in biological tissue could be described by the Pennes bio-heat equation [14]: (1) where ρ is the density (kg/m³); c is the specific heat (J/kg·°C); T is the temperature as a function of time (°C); k is the thermal conductivity (W/m·°C); W_b is the blood perfusion rate (kg/m³·s); C_b is the specific heat of blood (J/kg·°C); T_b is the temperature of blood (°C); q_r is the microwave power density, i.e., SAR (W/kg); q_m is the metabolism heat (W/kg); ▽² is the Laplacian operator and represents the divergence of temperature gradient. Because there is no metabolism or blood perfusion in the phantom, q_m and W_b in the Pennes bio-heat equation is ignored for the phantom heating. Thus, the left-hand side of the equation comprises the temperature-rising item and the thermal conductivity item, and the right-hand side of the equation is the SAR. So the Pennes bio-heat equation can be simplified as follows: (2)

The actual temperature variation and thermal conductivity value could be obtained from the experiments and then the SAR distribution could be derived based on Equation (2)(2) . SAR distribution varied with location, and originated from some characteristic points’ temperature data, so not all the temperature changes could be analyzed [15]. To obtain SAR distribution, some characteristic points were selected firstly, and then their temperature curves were recorded during MW ablation. In the process of SAR calculation, the sum of temperature variation rate and thermal conductivity value was obtained, and was then used to simulate SAR distribution with data fitting method [16,17]. Studies have shown that SAR distribution was associated with properties of the phantom and antenna, yet did not have a linear relationship with application power [18], so SAR at each power must be determined separately.

Calculation of temperature-rising item. The temperature-rising item represented the energy for raising the temperature, i.e. the product of the temperature-rising rate (slope of temperature curve) dT/dt, phantom density ρ, and specific heat c. Because all the temperature curves exhibited the characteristics from a rapid ramp-up to a substantial steady state, the initial 100 s temperature data (the first fifty measured values) exhibited better linear relationship and were therefore chosen as valid data. This procedure was repeated twenty times to give average values for reducing the error. Then average temperature variation curves were plotted and fitted by a linear trendline. According to the linear trendline, the slope of the curve was directly extracted. It was then multiplied by density and specific heat to get the temperature-rising item at corresponding measuring points.

Calculation of thermal conductivity item. The thermal conductivity item represented the energy transferred to other ablation zones from one point, i.e. the product of thermal conductivity k and ▽²T. Consequently, second order partial derivatives of the temperature T in the R direction and z direction ( and ) needed to be calculated respectively. Because the temperature at some point was known, its second order partial derivative could be calculated by the forward-difference method. The first and second order partial derivatives of the temperature in the R direction were computed respectively by (3) and (4) where h was 0.5 cm, depending on the collected data. The first and second order partial derivatives of the temperature in the z direction were calculated in the same way, as expressed in Equation (5)(5) : (5)

Finally, the second order partial derivatives of some points in the R direction and Z direction were derived using Equations (4)(4) and (5)(5) , and were then multiplied by thermal conductivity k to obtain the thermal conductivity item k▽²T at corresponding measuring points.

2.4 Finite element modeling

The COMSOL Multiphysics analysis software (COMSOL Inc., Palo Alto, CA, USA) was used to carry out pre-processing, simulation computation, and post-processing of finite element models for validating the SAR distribution. Based on the antenna structure and the simulation computation simplification principle, the 2 D planar model of the MW ablation phantom was established to represent the actual 3 D axisymmetric model, where the antenna central axis was the symmetry axis and the center of the aperture front end was set to be the origin. Further, the antenna length was 180 mm; the antenna diameter was 1.9 mm; and the distance from the aperture front end to the tip of the antenna was 11 mm. The tip had a length of 3 mm and was composed of inner core, aperture and outer shell. The phantom was disposed to have a length of 200 mm and a width of 100 mm, and the positions of thermocouples were marked in the phantom.

The modeling procedure employed free meshing and implemented mesh refinement in the range from the antenna tip to the aperture, to ensure sufficiently fine meshes and accurate thermal simulation [19]. The modeling procedure employed bio-heat transfer model to simulate the thermal ablation. The source item was SAR distribution function. The ablation time was 600 s. The initial temperature of the phantom and the antenna was set to be 25 °C. The phantom boundary was set to be insulation boundary. The thermal injury was computed by Arrhenius equation [20].

3. Results

3.1 SAR distribution

The sum of the temperature-rising item and the thermal conductivity item led to SAR. Depending on antenna design features and phantom properties, SAR distributions were fitted by an exponential function in the r direction and were fitted by a cubic polynomial in the z direction [21]. Further, the SAR distributions in the z direction were fitted separately on the front side (z > 0) without water-cooling and the rear side (z < 0) with water-cooling. The experimental results showed that this regressed form was optimal in fitting the SAR distributions and the correlation coefficients were greater than 0.96. The fitting equation was represented by (6)

Finally, the SAR distributions were fitted with different weights by the coefficient of variation method. The SAR distributions of 2450 MHz MW antenna at 40 W, 45 W, 50 W, 55 W, 60 W power groups were shown in Table 2.

Table 2. The SAR functions of 2450 MHz microwave antenna.

Microwave power	SAR distribution functions (W/kg)		Correlation coefficient (R-square)
	r direction	z direction
40 W	f(r) = 4419.63e^−3.509r	f(z)Front = 0.376z^3–1.433z^2–0.68z + 4.816	0.979
	f(r) = 4419.63e^−3.509r	f(z)Rear =0.241z^3–0.103z^2–0.197z + 4.816	0.966
45 W	f(r) = 3796.26e^−3.070r	f(z)Front = 0.307z^3–1.166z^2–0.566z + 3.92	0.981
	f(r) = 3796.26e^−3.070r	f(z)Rear = 0.191z^3–0.975z^2–0.106z + 3.92	0.975
50 W	f(r) = 3114.95e^−2.984r	f(z)Front = 0.205z^3–0.709z^2–0.773z + 3.181	0.989
	f(r) = 3114.95e^−2.984r	f(z)Rear = 0.146z^3–0.794z^2+0.0716z + 3.181	0.971
55 W	f(r) = 4266.36e^−3.408r	f(z)Front = 0.364z^3–1.309z^2–1.017z + 5.044	0.981
	f(r) = 4266.36e^−3.408r	f(z)Rear = 0.237z^3–1.263z^2+0.0356z + 5.044	0.985
60 W	f(r) = 3776.64e^−3.081r	f(z)Front = 0.277z^3–1.005z^2–0.748z + 3.842	0.978
	f(r) = 3776.64e^−3.081r	f(z)Rear = 0.175z^3–0.941z^2+0.0351z + 3.842	0.972

3.2 Simulation results of temperature field

The computed SAR source loads were applied to the finite element model to obtain corresponding simulation results. Take the SAR distribution of the 60 W power group for example (Figure 3), the ablation zone was similar to an oval shape. The black line was 54 °C isotherm (i.e. ablation zone boundary). The highest temperature was 123.6 °C and the lowest temperature was 25 °C. In order to highlight the ablation area, the area with the temperature more than 54 °C was considered as the thermal coagulation zone and the boundary was displayed in a black line (54 °C). The temperature field isotherms took oval like shapes with the antenna as an axis of symmetry, and the temperature was gradually lowered in the outward direction. In the central region around the aperture, the isotherms were much denser; the temperature was much higher; and the temperature gradient was larger.

Figure 3. Temperature distributions for 60 W at 600 s.

Figure 4 shows the measurement values and simulation values of three measuring points labeled 13, 17, 20 respectively for 40 W, 45 W, 50 W, 55 W, 60 W power groups at 600 s. As can be seen from the plots in Figure 4, the more consistent temperature variations can be obtained for the measured temperatures and simulation values, especially in far field areas. The early temperature at point (1.0,0.5) quickly rose due to the action of MW energy, and then gradually became stable. The temperature-rising curves at points (1.5,0.5) and (2.0,0.5) were relatively consistent, because these points were further from the antenna aperture so that their temperature variations could be mainly attributed to thermal conductivity. To indicate the errors between simulation data (T_i^s) and measurement data (T_i^m), maximum error (α), average error (β) and standard deviation (δ) were given by Equation (7)(7) , (8)(8) and (9)(9) , respectively. (7) (8) (9) where i represented sequence number of sampling points; n denoted total number of sampling points.

Figure 4. The measurement values and simulation values of labeled points 13,17,20 respectively for 40 W, 45 W, 50 W, 55 W, 60 W power group.

Table 3 gave the errors of three measuring points labeled 13, 17, 20 for 40 W, 45 W, 50 W, 55 W, 60 W power groups. As listed in Table 3, the maximum error was generally less than 6 °C, the average error was generally less than 2 °C, and the standard deviation was less than 1 °C. The simulation errors could meet the essential clinical requirements.

Table 3. The error between simulation data and measurement data respectively for 40 W, 45 W, 50 W, 55 W, 60 W power.

	Power (W)
Point (r,z) (cm)	40 W	45 W	50 W	55 E	60 W
(1.0,0.5)
α (°C)	7.34	5.89	5.39	5.04	3.97
β (°C)	2.54	2.20	1.89	2.22	1.73
γ (°C)	2.36	1.78	1.69	1.21	1.33
(1.5,0.5)
α (°C)	1.19	1.26	1.29	2.63	2.31
β (°C)	0.56	0.74	0.86	0.96	1.43
γ (°C)	0.25	0.36	0.57	0.57	0.32
(2.0,0.5)
α (°C)	0.61	0.85	1.71	1.20	1.03
β (°C)	0.37	0.52	0.54	0.59	0.46
γ (°C)	0.19	0.27	0.35	0.29	0.28

3.3 Simulation results of ablation zone sizes

The ablated lesions could be divided into three zones: carbonization zone, coagulation zone, and inflammatory reaction zone [22]. In this research, the ablation coagulation zone’s effective boundary temperature was set as 54 °C. The MW ablation zones were characterized by a long diameter, a short diameter and a prepositive distance, as shown in Figure 5. The long diameter represented the maximum length of ablation zones along the z axis. The short diameter represented the maximum length of ablation zones along the r axis. The prepositive distance represented the length from the antenna tip to the front side of the ablation zone boundary. Table 4 shows the simulation and measurement values of ablation zones for the 60 W power group at 600 s. It could be seen that these characteristic lengths had excellent consistency.

Figure 5. The characteristic dimensions of the ablation zone in RZ plane.

Table 4. The measured and simulated values of ablation zones at 60 W.

60 W & 600 s	Long diameter (cm)	Short diameter (cm)	Prepositive distance (cm)
Reference value	5.5 ± 0.3	3.4 ± 0.2	0.3 ± 0.05
Simulation value	5.191	3.362	0.341

4. Discussion

It should be noted that SAR was the integral effects of the electrical parameters, but these parameters were generally not available especially for tumor tissues and thus took fixed values or empirical data in the simulation process. This will greatly affect the simulation precisions. Furthermore, directly measuring these electrical parameters is impracticable in the MWA procedure. Therefore, the precise modeling of SAR is key to the temperature simulation and prediction. During the simulation process, experimentally measured SAR can be used as the heat source to the Pennes bio-heat equation, rather than employing complex, computationally intensive electromagnetic solvers. Because the experimentally measured SAR is much closer to the actual characteristics of the tissue, the simulation results have higher accuracy. The simulation speed is also greatly improved.

During the computation of SAR, the thermal conductivity item was taken into consideration and was approximately computed by FDM. Next, finite element models were built to simulate the temperature field distribution for various power groups. Finally, simulation results were verified based on the phantom experiments. Previous studies [21] proposed the SAR equation fitting method, but didn’t consider thermal conduction effects. Therefore, the novelty of this research is to minimize the errors due to thermal conduction in their experimental measurement of SAR.

The SAR functions were experimentally obtained by the thermocouple measurements. The thermocouples consisted of metallic materials, which might lead to RF interference and produce some amount of temperature errors. However, the errors were also related to the size of the thermocouples. The study tried to reduce the influence of the thermocouples by the following ways. Firstly, The thermocouples we used were very thin and the temperature resolution was about 0.1 °C, which greatly reduced the interference between the electromagnetic field and the thermocouples; Secondly, in the solving process, we mainly employed the first 100 seconds of heating time to calculate SAR, so thermocouples’ heat transfer effects were limited; Thirdly, inserting directions between the thermocouples and the MW antenna were perpendicular, in order to minimize the contact area between the thermocouple and the phantom. As such, the simulation accuracy could be guaranteed.

As compared with traditional simulation methods, the SAR method will enable tissue-specific simulation, simplified modeling of treatment outcome and improved temperature prediction accuracy. The static SAR over the course of the ablation was obtained and was suitable for the liver ablation cases because of the constant powers. Due to the changes in tissue electrical properties, SAR would likely change with temperature during an ablation. So during the ablation process, a thermocouple can be inserted into living tissues to provide feedback for SAR changes. So the SAR method can be applied for tissue-specific simulation by incorporating real-time measurement with a thermocouple. Figure 6 depicts the comparison plots of measurement values, simulation values by SAR method, and simulation results by traditional electromagnetic coupling method for 60 W power group. As can be seen in this Figure, better agreement is achieved between the measured temperatures and simulation values by SAR method. The advantages of the SAR method are also apparent from the error comparisons between Tables 3 and 5. By using 60 W power for example, the errors between simulation data by SAR method and measurement data are significantly lower.

Figure 6. The comparison plots of measurement values, simulation values by SAR method, and simulation results by traditional method for 60 W power group.

Table 5. The error between traditional simulation data and measurement data for 60 W power.

	Point (r,z) (cm)
Error (°C)	(1.0,0.5)	(1.5,0.5)	(2.0,0.5)
α	5.89	2.35	1.56
β	2.51	1.09	1.01
γ	2.52	0.50	0.51

More importantly, the temperatures of far field points are of particular interest in clinical hyperthermia applications, because clinicians need to know the ablative margins. According to the research, the SAR method can better predict far field temperatures and the average errors for far field points are less than or near 1 °C. Therefore, this method has a relatively higher value of practical application.

Furthermore, the exact derivation of SAR function has especially scientific significance. According to Equation (2)(2) , the 54 °C isotherms over time can be easily obtained without having to carry out finite element analysis and computation as long as the SAR functions have sufficient accuracy. This can greatly reduce temperature prediction time of MWA and thus improve treatment effects of MWA.

In spite of the above advances, there is continuing desire to enhance the simulation accuracy. Blood vessels are intensive in human tissues and thermal-physical properties vary with temperature, which can affect the simulation results. The size of the ablation lesion is inversely related to the amount of blood flow and then the models with no blood flow overestimate the size of the ablation [23]. Therefore, experiments considering blood flow should be performed in the future, and temperature-dependent biological thermal parameters should be obtained. Since SAR seems to dominate in temperature simulation, further research is needed by actual measurement feedback to improve SAR accuracy.

5. Conclusions

This paper has proposed a temperature simulation method based on an improved SAR numerical model. The SAR distribution functions for different powers were obtained by incorporating the thermal conductivity item. Results show that the SAR-derived simulation method can greatly improve the simulation precision and speed. These researches can provide some useful reference to the temperature predictions in the clinic MWA procedures.

Disclosure statement

The authors declared they have no competing interests. The authors alone are responsible for the content and writing of the paper.

Additional information

Funding

This work was supported by Beijing Natural Science Foundation [Grant No. 7174279]; National Natural Science Foundation of China [Grant No. 71661167001]; Beijing Education Committee Foundation Program for the Top Young Innovative Talents [Grant No. CIT&TCD201404053]; and National Key Technology R&D Program [Grant No. 2015BAI02B03].