Conformal coverage of liver tumors by the thermal coagulation zone in 2450-MHz microwave ablation

Abstract

Purpose: To optimize treatment schemes using 2450-MHz microwave ablation (MWA), a novel conformal coverage method based on bipolar-angle mapping is proposed that determines whether a liver tumor is completely encompassed by thermal coagulation zones.

Materials and methods: Firstly, three-dimensional (3-D) triangular mesh data of liver tumors were reconstructed from clinical computed tomography (CT) slices using the Marching cubes (MC) algorithm. Secondly, characterization models of thermal coagulation zones were established based on finite element simulation results of 40, 45, 50, 55, and 60 W ablations. Finally, coagulation zone models and tumor surface data were mapped and fused on a two-dimensional (2-D) plane to achieve conformal coverage of liver tumors by comparing the corresponding polar radii.

Results: Optimal parameters for ablation treatment of liver tumors were efficiently obtained with the proposed conformal coverage method. Fifteen liver tumors were obtained with maximal diameters of 12.329–78.612 mm (mean ± standard deviation, 39.094 ± 19.447 mm). The insertion positions and orientations of the MWA antenna were determined based on 3-D reconstruction results of these tumors. The ablation patterns and durations of tumors were planned according to the minimum mean standard deviations between the ablative margin and tumor surface.

Conclusion: The proposed method can be applied to computer-assisted MWA treatment planning of liver tumors, and is expected to guide clinical procedures in future.

Keywords:

• Conformal coverage
• three-dimensional reconstruction
• microwave ablation
• coagulation zone
• treatment planning
• liver tumor

1. Introduction

Primary liver tumor is the 6th most common malignancy worldwide with approximately 745,000 deaths annually [1]. Microwave ablation (MWA) has emerged as an alternative treatment for liver tumors and is increasingly utilized in the inactivation of tumor cells [2–6]. The effects of MWA treatment depend on therapeutic temperature. Generally, coagulation necrosis will occur instantly at temperatures of 55 °C or higher [7,8]. For liver tumors, a temperature of 54 °C is often regarded as the thermal coagulation threshold [9,10]. In tumor ablation procedures, it is necessary that thermal coagulation zones are sufficiently wide to cover all the tumor tissue as well as a 5–10 mm safety margin [11]. A safety margin over 5 mm could reduce local tumor progression (LTP) significantly. Therefore, in this work, we analyze the conformal treatment of liver tumors by an isothermal surface of 54 °C (IS-54).

Image-guided MWA technology is a precise method of selectively and locally inducing coagulation necrosis of tumors inside the body [12,13]. Preoperative MWA planning based on two-dimensional (2-D) medical images (computed tomography (CT) scans, magnetic resonance imaging (MRI) scans, ultrasonic images) of patients, and temperature distribution simulation models is critical to the optimal success of ablative therapies [14]. The core mission of tumor MWA treatment planning is to ensure the conformal ablation of liver tumors and the minimum damage to healthy tissue. Preoperative planning typically involves the proper insertion positions of the antenna into the targeted tumor, as well as the requisite ablation power and ablation duration (typically 5–20 min) to achieve a thermal coagulation zone. However, MWA clinical applications remain hindered by a lack of reliable, automatic, and personalized treatment planning [15].

Three-dimensional (3-D) surface results of liver tumors are advantageous for displaying spatial relationships. In this study, 3-D surface results were quantitatively analyzed to obtain the coordinates of the gravity center (GC), the length of the maximum diameter (MD), and the orientation of the maximum diameter through the gravity center. Using the evolving properties of MWA coagulation zones over time, IS-54 surfaces for 40, 45, 50, 55, and 60 W were characterized by characteristic length growth models and shape variation factors. The 3-D surface results of liver tumors and IS-54 were then projected and mapped onto a 2-D plane. The purpose of this study was to simply and efficiently achieve conformal coverage of liver tumors by IS-54 and output the optimal MWA treatment parameters for clinicians.

2. Materials and methods

2.1. Surface reconstruction of liver tumors

Surface images of liver tumors were produced to obtain the accurate dimensional measurements required for treatment planning. Among various medical imaging modalities, CT is often used to acquire patient-specific anatomical data of liver tumors due to its higher signal-to-noise ratio and better spatial resolution [16]. The experimental data of this study include two datasets. The first dataset consists of 3 CT scans from the General Hospital of PLA (Beijing, China), where each CT series contains 512 × 512 × 60 voxels and the size of each voxel is 0.703 mm × 0.703 mm × 2.5 mm. The second dataset consists of 12 CT scans from the public dataset 3Dircadb1 (http://www.ircad.fr/research/3dircadb), where each CT series contains 512 × 512 × (74–260) voxels and the size of each voxel is (0.57–0.87) mm×(0.57–0.87) mm×(1–4) mm.

Liver tumors in these datasets were manually segmented by clinical experts. 3-D surface results of liver tumors were reconstructed by the MC algorithm based on these segmented data and then output as STL files. These STL files were read and analyzed using a VC++ 6.0 platform to derive the position and size information of liver tumors.

2.2. Finite element method (FEM) simulation of MWA and experimental validation

Ex vivo porcine livers were employed in this study because they exhibited similar properties to human livers. 2-D planar axis-symmetry FEM models for 40, 45, 50, 55, and 60 W ablations were built. The FEM models were designed to have a length of 200 mm and a width of 100 mm and were meshed by the free meshing method. The temperature distributions within ex vivo porcine livers were obtained by solving the Helmholtz harmonic wave equation (Equation (1)(1) $$\nabla \times {\mu }_r^{-1}(\nabla \times E)-k_0^2({\epsilon }_r-\frac{j\sigma }{\omega {\epsilon }_0})E=0;$$(1) ) [17] and the Pennes bio-heat transfer equation (Equation (2)(2) $$\rho \mathit{SAR}=\frac{1}{2}\sigma |E{|}^2;$$(2) ) [18]. Because there was no metabolism or blood perfusion in the ex vivo porcine liver, the Pennes bio-heat transfer equation was simplified to Equation (4)(4) $$\rho \cdot \mathit{Cp}\cdot \frac{\mathit{dT}}{\mathit{dt}}=\nabla \cdot k\nabla T+\rho \mathit{SAR},$$(4) . (1) $$\nabla \times {\mu }_r^{-1}(\nabla \times E)-k_0^2({\epsilon }_r-\frac{j\sigma }{\omega {\epsilon }_0})E=0;$$(1) (2) $$\rho \mathit{SAR}=\frac{1}{2}\sigma |E{|}^2;$$(2) (3) $$\rho \cdot \mathit{Cp}\cdot \frac{\mathit{dT}}{\mathit{dt}}=\nabla \cdot K\nabla T+\rho Q+\rho \mathit{SAR}+{\rho }_b{c}_{b}w(T-T_b);$$(3) (4) $$\rho \cdot \mathit{Cp}\cdot \frac{\mathit{dT}}{\mathit{dt}}=\nabla \cdot k\nabla T+\rho \mathit{SAR},$$(4) where ${\epsilon }_r$ is the relative dielectric constant of biological tissue; ${\epsilon }_0$ is the relative dielectric constant of vacuum; $\sigma$ is the tissue conductivity (S/m); $E$ is the electric field intensity (V/m); ${\mu }_r$ is the relative permeability (${\mu }_r$=1); $k_0$ is the free space wave number; $\rho$ is the tissue density (kg/m³); $\mathit{SAR}$ is the specific absorption rate (W/kg); T is temperature (°C); $t$ is time (s); $\mathit{Cp}$ is the specific heat capacity (J/kg/°C); $K$ is the thermal conductivity (W/m/°C); $w$ is the blood perfusion rate (1/s); $Q$ is the metabolic heat (W/kg); and the subscript $b$ represents the blood property.

The initial temperature and the boundary conditions of ex vivo porcine livers in these models were set to 21 °C. The power-feeding point was modeled using a power port-boundary condition and the water-cooled effects of MW antenna were characterized by the convective heat transfer coefficient (set to 2400 W/m²/°C). The temperature-dependent forms of the characteristic parameters [19] were employed in the FEM models, as represented by Equations (5)–(9): (5) $$\mathit{Cp}\left(T\right)=\{\begin{array}{cc}3628+9.812\cdot (T-25)\hfill & T\le 70°\mathrm{C}\hfill \\ 4069.540-2.152\cdot \frac{\partial W\left(T\right)}{\partial T}\hfill & T\ge 70°\mathrm{C}\hfill \end{array};$$(5) (6) $$W\left(T\right)=\{\begin{array}{cc}0.778\cdot (1-\hspace{0.17em}\text{exp}\hspace{0.17em}(\frac{T-106}{3.420}\left)\right)\hfill & 70°\mathrm{C}\le T<100°\mathrm{C}\hfill \\ 7.053-0.0640\cdot T\hfill & 100°\mathrm{C}\le T<104°\mathrm{C}\hfill \\ 0.778\cdot \hspace{0.17em}\text{exp}\hspace{0.17em}(-\frac{T-80}{34.320})\hfill & 104°\mathrm{C}\le T\hfill \end{array};$$(6) (7) $$K\left(T\right)=0.565+0.00111\cdot (T-25);$$(7) (8) $$\epsilon \left(T\right)=44.3\cdot \left[1-\frac{1}{1+\hspace{0.17em}\text{exp}\hspace{0.17em}(5.2-0.0519\cdot T)}\right];$$(8) and (9) $$\sigma \left(T\right)=1.8\cdot \left[1-\frac{1}{1+\hspace{0.17em}\text{exp}\hspace{0.17em}(5.324-0.0598\cdot T)}\right],$$(9) where $T$ is temperature (°C); $\mathit{Cp}(T),$ $W(T),$ $K(T),$ $\epsilon (T),$ and $\sigma (T)$ represent the temperature-dependent functions of specific heat capacity (J/kg/°C), tissue water content (kg/m³), thermal conductivity (W/m/°C), relative dielectric constant, and electrical conductivity (S/m), respectively.

To verify the accuracy of FEM models, 20 sets of ex vivo pig liver experiments were conducted under the same conditions. The experimental setup included a cooled-shaft microwave unit of 2450-MHz (KY-2000; Kangyou Microwave Energy Sources Institute, Nanjing, China), a peristaltic pump (BTDI-100; Baoding Longer Precision Pump, Inc, Hebei, China), thermometers (YWY-2; Kangyou Microwave Energy Sources Institute, Nanjing, China), a steel ruler, and ex vivo pig livers (Figure 1). This study used MWA power settings of 40, 45, 50, 55, and 60 W. The positions of the thermometers were disposed at V₁ (10 mm, −5 mm), V₂ (10 mm, −15 mm), V₃ (10 mm, 0 mm), V₄ (20 mm, 5 mm), V₅ (25 mm, 5 mm), and V₆ (10 mm, 10 mm). These thermometers were used to measure the temperature changes to validate the temperature prediction accuracy of the FEM models. In addition, a steel ruler was used to measure the size of thermal coagulation zones after MWA.

Figure 1. MWA experimental setup for validation.

2.3. Characterization of IS-54 evolution

In previous studies, IS-54 on the r-z plane was often expressed by characteristic lengths, whose growing functions with time take logarithmic forms [20]. These characteristic lengths generally include the long-axis and short-axis diameters of thermal coagulation zones. The ablated zone was sectioned along the antenna shaft, and the diameter of the zone along the antenna insertion axis was measured as the long-axis diameter. The short-axis diameter was orientated perpendicular to the long-axis diameter. According to the water-cooled properties of this MW antenna, the forward long-axis length $L_f,$ backward long-axis length $L_b,$ and short-axis radius $L_s$ are taken as the characteristic lengths $L$ (Figure 2) to more accurately characterize the change of coagulation zone size. Previous experiments have demonstrated that $L$ can be represented by power functions (similar to logarithmic forms), as illustrated in Equation (10)(10) $$L=p_1+p_2\cdot t^{p_3},$$(10) : (10) $$L=p_1+p_2\cdot t^{p_3},$$(10) where $p_1-p_3$ are the coefficients of the characteristic length growth models.

Figure 2. The schematic view of characteristic lengths. $L_f$ is the forward long-axis length, $L_b$ is the backward long-axis length, $L_s$ is the short-axis radius, ${\phi }_R$ is the rotation polar-angle, ${\phi }_O$ is the offset polar-angle, and $\lambda$ is the polar radius. The origin is disposed at the MW emission point. The blue line indicates the thermal coagulation zone.

The changing size of the thermal coagulation zones over time was obtained from the FEM models and $p_1-p_3$ was then fitted by 1stOpt (7 D-Soft High Technology, Inc., Beijing, China), which is a mathematical optimization analysis toolbox. The Levenberg-Marquardt (LM) optimization algorithm was used to fit functions. Because the growth patterns of thermal coagulation zones in the initial MWA stage were irregular, this study predominantly analyzed the growth model after 100 s of ablation. Based on the properties of the MW antenna, the forward IS-54 and the backward IS-54 were separately modeled using the selected characteristic lengths [21]: (11) $$z_f(t,r)=(L_f-(L_f/L_s{}^2)\cdot r^2)/(1+\alpha \cdot (k_1\cdot r^2+k_2\cdot r^4\left)\right);$$(11) (12) $$z_b(t,r)=(L_b-(L_b/L_s{}^2)\cdot r^2)/(1+\beta \cdot (q_1\cdot r^2+q_2\cdot r^4)+q_3\cdot r^6+q_4\cdot r^8+q_5\cdot r^{10}),$$(12) where $k_1-k_2$ and $q_1-q_5$ are the coefficients of the forward IS-54 and the backward IS-54, respectively; $\alpha$ and $\beta$ are the forward shape variation factor and the backward shape variation factor, respectively. These shape variation factors are used to deal with the minor variability of IS-54 shapes over time and also take the form of power functions: (13) $$\alpha =m_1+m_2\cdot t^{m_3};$$(13) (14) $$\beta =n_1+n_2\cdot t^{n_3},$$(14) where $m_1-m_3$ and $n_1-n_3$ are the coefficients of the shape variation factors.

Based on the FEM data at some discrete time (100 s, 150 s, 200 s, etc), the IS-54 models were precisely obtained from Equations (10)–(14). Then, to validate the effectiveness and applicability of this modeling method, the derived IS-54 models were compared with the FEM data at some other time (206 s, 302 s, 408 s, 466 s, 523 s, and 552 s).

2.4. Conformal coverage algorithm for treating physical liver tumors

Any point, M, on the 3-D thermal coagulation zone models or physical liver tumors can be characterized by a polar presentation: polar radius $\lambda ,$ rotation polar-angle ${\phi }_R$[0, 2π], and offset polar-angle ${\phi }_O$[-π/2, π/2], as shown in Figure 2. Conformal coverage of liver tumors by IS-54 can be achieved by comparing the corresponding polar radii in the case of ${\phi }_R$ alignment and ${\phi }_O$ alignment. As such, a novel conformal coverage algorithm based on bipolar-angle mapping is proposed in this study (Figure 3).

Figure 3. Work flowchart of the conformal coverage algorithm.

2.4.1. 2-D Projection and mapping of physical liver tumors

Preoperative CT slices of liver tumors after manual segmentation were loaded to form 3-D volume data. Then, 3-D surface reconstructions of ROI (region of interests) were conducted using the MC algorithm. The surface rendering results of liver tumors were quantified to derive the length of the tumor maximum diameter ($l_{\mathit{MD}}$), the coordinate of the tumor gravity center ($P_0(x_0, y_0 , z_0)$), and the orientation vector of the maximum diameter through the gravity center ($\overrightarrow{v}$), where $l_{\mathit{MD}}$ indicates the tumor’s size, $P_0(x_0, y_0 , z_0)$ is the insertion position of the MW antenna [22], and $\overrightarrow{v}$ represents the insertion direction of the MW antenna. If the 3-D surface result of the tumor was composed of $n$ points ($\{P_i(x_i, y_i , z_i)|1\le i\le n\}$), then $P_0(x_0, y_0 , z_0),$ $l_{\mathit{MD}},$ and $\overrightarrow{v}$ were computed as follows: (15) $$P_0(x_0,y_0,z_0)=\frac{1}{n}\sum _{i=1}^n{P}_i=\frac{1}{n}\sum _{i=1}^n(x_i,y_i,z_i);$$(15) (16) $$l_{\mathit{MD}}=\max \left(\sqrt{{(x_i-x_j)}^2+{(y_i-y_j)}^2+{(z_i-z_j)}^2}\right)\hspace{1em}(i\ne j);$$(16) (17) $$\overrightarrow{v}=(x_m-x_0,y_m-y_0,y_m-y_0),$$(17) where $i$ represents the index of the tumor points, and point $P_m(x_m, y_m , z_m)$ is the endpoint of the maximum diameter through the gravity center.

To ensure that the gravity center of the tumor was located at the origin (i.e. the microwave emission point), coordinate translation was performed as follows: (18) $$P_i=P_i-P_0.$$(18)

To reduce the computational loads, 3-D coordinates $P_i(x_i, y_i , z_i)$ of the tumor were projected into 2-D coordinates $P_i(r_i, z_i)$ on the r-z plane, as illustrated in Equation (19)(19) $$\begin{array}{c}{r}_i=\sqrt{x_i^2+y_i^2}\hfill \\ z_i=z_i\hfill \end{array}.$$(19) : (19) $$\begin{array}{c}{r}_i=\sqrt{x_i^2+y_i^2}\hfill \\ z_i=z_i\hfill \end{array}.$$(19)

Because the maximum diameter was required to align with the insertion direction of the MW antenna, coordinate rotation of the tumor was performed according to Equation (20)(20) $$\begin{array}{c}{r}_i\prime =r_i\cdot \hspace{0.17em}\text{cos}\hspace{0.17em}(\theta )+z_i\cdot \hspace{0.17em}\text{sin}\hspace{0.17em}(\theta )\hfill \\ z_i\prime =z_i\cdot \hspace{0.17em}\text{cos}\hspace{0.17em}(\theta )-r_i\cdot \hspace{0.17em}\text{sin}\hspace{0.17em}(\theta )\hfill \end{array},$$(20) : (20) $$\begin{array}{c}{r}_i\prime =r_i\cdot \hspace{0.17em}\text{cos}\hspace{0.17em}(\theta )+z_i\cdot \hspace{0.17em}\text{sin}\hspace{0.17em}(\theta )\hfill \\ z_i\prime =z_i\cdot \hspace{0.17em}\text{cos}\hspace{0.17em}(\theta )-r_i\cdot \hspace{0.17em}\text{sin}\hspace{0.17em}(\theta )\hfill \end{array},$$(20) where $\theta$ is the intersection angle between the maximum diameter and the z axis.

After coordinate transformation of the tumor surface data, the polar radius ${\lambda }_{\mathit{tumor}}$ was mapped based on bipolar-angles ${\phi }_R$ and ${\phi }_O$ (Equation (21)(21) $$\begin{array}{c}{\lambda }_{\mathit{tumor}}=\sqrt{r_i\prime \cdot r_i\prime +z_i\prime \cdot z_i\prime }\hfill \\ {\phi }_O=\arcsin (z_i\prime /{\lambda }_{\mathit{tumor}})\hfill \\ {\phi }_R=\arcsin (y_i/\mathit{abs}(r_i\left)\right)\hfill \end{array}.$$(21) ), where ${\phi }_R$[0, 2π] was designated as the horizontal axis and ${\phi }_O$[-π/2, π/2] was designated as the vertical axis. (21) $$\begin{array}{c}{\lambda }_{\mathit{tumor}}=\sqrt{r_i\prime \cdot r_i\prime +z_i\prime \cdot z_i\prime }\hfill \\ {\phi }_O=\arcsin (z_i\prime /{\lambda }_{\mathit{tumor}})\hfill \\ {\phi }_R=\arcsin (y_i/\mathit{abs}(r_i\left)\right)\hfill \end{array}.$$(21) 2-D projection and mapping of physical liver tumors was then achieved through the above process.

2.4.2. Bipolar-angle mapping of IS-54

IS-54 models for 40, 45, 50, 55, and 60 W power levels were built according to Equations (10)–(14). The offset angle ${\phi }_O$[-π/2, π/2] of IS-54 was discretized into 1800 equal parts in order to align with discrete values of 3-D tumor surface results. However, discretization of ${\phi }_R$ was not required due to the rotational symmetry of IS-54. 2-D bipolar-angle mapping of IS-54 models was achieved based on ${\phi }_R$ and ${\phi }_O,$ where the corresponding polar radius ${\lambda }_{\mathit{IS}-54}$ was represented by color brightness. The color brightness was enhanced with increasing ${\lambda }_{\mathit{IS}-54}.$

2.4.3. Superposition and fusion of liver tumors and IS-54

2-D projection views of the liver tumors and IS-54 on the r-z plane were superimposed and their bipolar-angle maps were fused together in order to intuitively learn the conformal coverage conditions of the liver tumors.

2.4.4. Conformal coverage of liver tumors by IS-54

IS-54 gradually grows with MWA duration. As shown in Figure 4, when the minimum difference of corresponding polar radii was equal to only 5 mm (this value could be changed based on the specific safety margin requirement, which was set to 5 mm in this study), conformal coverage was accomplished; i.e., the safety margin was guaranteed and the optimal ablation time was obtained. Therefore, the objective function was defined as follows: (22) $$\min ({\lambda }_{\mathit{IS}-54}-{\lambda }_{\mathit{tumor}})=5.$$(22)

Figure 4. Conformal coverage illustration of liver tumor by IS-54. The blue dots indicate the mesh nodes of the tumor surface. $L_f$ is the forward long-axis length, $L_b$ is the backward long-axis length, $L_s$ is the short-axis radius, ${\phi }_O$ is the offset polar-angle, ${\lambda }_{\mathit{IS}-54}$ represents the polar radius of IS-54, and ${\lambda }_{\mathit{tumor}}$ represents the polar radius of the tumor surface model. The origin is set at the MW emission point.

When the objective function was satisfied, the corresponding lost function was computed. The lost function was defined as the standard deviation (SD) between the ablative margin (AM) and the tumor surface: (23) $$\delta =\sqrt{\frac{1}{n}\sum _{i=1}^n{\left({\lambda }_{\mathit{IS}-54}\right(i)-{\lambda }_{\mathit{tumor}}(i\left)\right)}^2}\hspace{1em}i=1,2,\dots ,n.$$(23)

By changing power levels, the corresponding SD was computed. The suitable MWA parameters (ablation power, ablation duration) were output based on the minimum SD. In addition, the gravity center of the tumor $P_0(x_0, y_0 , z_0)$ was designated as the insertion position of the MW antenna and the orientation vector of the maximum diameter through the gravity center ($\overrightarrow{v}$) determined the insertion orientation of the MW antenna. According to this algorithm, optimal MWA treatment schemes could be determined.

3. Results

3.1. 3-D Reconstruction results of liver tumors

In order to test the applicability of the conformal coverage algorithm, 15 liver tumors with different shapes and different sizes were obtained and reconstructed, as shown in Figure 5. The 3-D surface data as STL files were quantitatively analyzed using the VC++ 6.0 platform to output the gravity center coordinates, the maximum diameter lengths, and the maximum diameter vectors through the gravity centers of these tumors, where the upper-left corner of the first CT slice was designated as the origin.

Figure 5. Surface reconstruction results of liver tumors.

The maximal diameter of these tumors ranged from 12.329 mm to 78.612 mm (mean ± SD, 39.094 ± 19.447 mm). The detailed data of these surface results are listed in Table 1. Based on this detailed information, the physical tumors were precisely characterized.

Table 1. Detailed data of tumor reconstruction results.

Case No.	Coordinate of GC (mm)	Length of MD (mm)	Vector of MD through GC
1	(78.963, 178.249, 94.768)	23.288	(2.599, 1.047, 94.768)
2	(126.804, 150.569, 83.236)	78.612	(−6.526, 23.102, −31.986)
3	(57.083, 147.263, 47.089)	34.177	(−0.271, −1.576, −17.089)
4	(101.743, 191.876, 125.884)	33.914	(0.287, −8.336, −13.878)
5	(129.562, 163.574, 109.898)	37.447	(2.596, −6.392, −18.698)
6	(113.836, 109.864, 172.459)	39.087	(1.789, 4.511, −13.694)
7	(101.536, 138.872, 159.527)	70.266	(−9.405, −7.362, −25.108)
8	(99.055, 146.384, 91.701)	30.309	(−4.055, −9.821, −9.201)
9	(162.682, 74.405, 92.662)	37.537	(2.252, −5.964, −22.249)
10	(86.428, 194.483, 68.258)	73.064	(−11.346, −10.271, −32.235)
11	(113.189, 137.676, 158.726)	27.855	(19.799, −10.109, −39.435)
12	(94.541,168.543, 96.074)	28.631	(−5.196, 2.467, −12.861)
13	(141.204, 89.013, 186.047)	22.721	(10.416, −1.233, −3.641)
14	(60.457, 98.598, 176.138)	12.329	(−1.747, −3.408, −4.932)
15	(120.144, 126.804, 185.016)	37.171	(−1.014, −11.664, −15.409)
Average	N/A	39.094	N/A

N/A: not applicable.

3.2. FEM simulation results and experimental verification

Twenty sets of experiments were conducted on ex vivo porcine livers under 40, 45, 50, 55, and 60 W power levels, respectively. The FEM simulation results of MWA temperature distributions and coagulation zones were then validated using these experimental data.

Figure 6 shows the temperature comparison plots of the measurements and the simulations at several verification points. The average of the maximum errors, the average errors, and the average standard deviations between simulation temperatures and measurements at V₁, V₂, V₃, V₄, V₅, and V₆ were 2.879 °C, 1.126 °C, and 0.836 °C, respectively.

Figure 6. The temperature comparison plots of the measurements and the simulations at (a) the point V₁, (b) the point V₂, (c) the point V₃, (d) the point V₄, (e) the point V₅, and (f) the point V₆.

According to gross pathologic examination, the ablation zones included a central “white coagulation zone” surrounded by a “red hemorrhagic zone” [23]. The white coagulation zone was generally accepted to represent coagulated tissue and was consistent with irreversible cell death [24]. The characteristic lengths were measured to validate the simulation size of the thermal coagulation zones. The errors between the measurements and the simulation data were within ±5%. The comparison of characteristic lengths for the coagulation zones demonstrates that the FEM models can accurately simulate the actual ablation effects of ex vivo porcine livers and effectively indicate the changes of thermal coagulation zones. Figure 7 shows some sectional images of porcine livers after 600 s of MWA at 45, 50, and 60 W power levels, where the blue line represents IS-54 obtained from FEM models. It is clear that FEM data can effectively reflect the changes of thermal coagulation zones. Therefore, this study successfully derived the characterization functions of IS-54 over time based on FEM models.

Figure 7. The sectional images of porcine livers after 600 s at 45, 50, and 60 W power levels. IS-54 indicated by the blue line was overlain onto the pathology images of the ablated livers.

3.3. Characterization of IS-54 evolution

According to FEM model data and the function forms of IS-54 (Equations (10)–(14)), the corresponding coefficients of IS-54 models were fitted by 1stOpt, as summarized in Tables 2 and 3. It was demonstrated in this study that the shape variation characteristics of IS-54 at 40, 45, 50, 55, and 60 W were extremely similar. Therefore, $\alpha$ and $\beta$ for different power levels employed the same power function forms. The expressions of $\alpha$ and $\beta$ were fitted by 1stOpt, as shown by Equations (24)(24) $$\alpha =0.1049+583.8803\cdot t^{-1.0124}\hspace{1em}({\mathrm{R}}^2 = 0.9993);$$(24) and (25)(25) $$\beta =1.023-6072.0134\cdot t^{-1.9404}\hspace{1em}({\mathrm{R}}^2 = 0.9885).$$(25) : (24) $$\alpha =0.1049+583.8803\cdot t^{-1.0124}\hspace{1em}({\mathrm{R}}^2 = 0.9993);$$(24) (25) $$\beta =1.023-6072.0134\cdot t^{-1.9404}\hspace{1em}({\mathrm{R}}^2 = 0.9885).$$(25)

Table 2. Characteristic length growth models at various ablation powers.

Power levels	Characteristic length (mm)	p1	p2	p3	Correlation coefficient (R^2)
40 W	Lf	−0.0132	3.586	0.257	0.9992
	Lb	−28.860	75.070	−0.393	0.9989
	Ls	−8.223	5.743	0.215	0.9999
45 W	Lf	−0.0333	3.938	0.247	0.9992
	Lb	−33.745	59.735	−0.282	0.9995
	Ls	−8.959	6.407	0.207	0.9999
50 W	Lf	0.0139	3.916	0.252	0.9986
	Lb	−31.800	65.540	−0.349	0.9987
	Ls	−9.322	6.797	0.203	0.9999
55 W	Lf	−0.00219	4.233	0.242	0.9989
	Lb	−33.438	63.865	−0.330	0.9980
	Ls	−9.763	7.240	0.199	0.9999
60 W	Lf	−3.134	6.698	0.195	0.9986
	Lb	−34.978	59.289	−0.290	0.9989
	Ls	−10.502	7.934	0.190	0.9998

Table 3. Corresponding coefficients of IS-54.

Powers
Coefficients	40 W	45 W	50 W	55 W	60 W
k1	−0.00169	−0.00159	−0.00150	−0.00145	−0.00141
k2	−1.104 × 10^−5	−9.152 × 10^−6	−8.180 × 10^−6	−6.981 × 10^−6	−6.497 × 10^−6
q1	−0.0149	−0.0131	−0.0118	−0.0108	−0.0114
q2	0.000439	0.000359	0.000303	0.000265	0.000266
q3	−5.503 × 10^−6	−4.186 × 10^−6	−3.320 × 10^−6	−2.747 × 10^−6	−2.634 × 10^−6
q4	2.904 × 10^−8	2.050 × 10^−8	1.523 × 10^−8	1.192 × 10^−8	1.097 × 10^−8
q5	−5.560 × 10^−11	−3.643 × 10^−11	−2.537 × 10^−11	−1.876 × 10^−11	−1.660 × 10^−11

Therefore, the precise mathematical functions of IS-54 curves could be obtained, where thermal coagulation zone sizes were determined by the characteristic length growth models and the shape changes of these zones were dependent on the shape variation factors.

The shapes and sizes of IS-54 simulation isotherms and the prediction curves were generally in good agreement. At 206 s, 302 s, 408 s, 466 s, 523 s, and 552 s, the standard deviations between the simulations and the predictions were 0.777 mm, 0.772 mm, 0.799 mm, 0.817 mm, 0.818 mm, and 0.881 mm, respectively. Furthermore, the paired student t-test in SPSS (SPSS Inc., Chicago, USA) was used to identify differences between these results, and the p values were .0835, .670, .635, .894, .110, and .725, respectively. All p values were greater than .05, so there was no statistically significant difference. In summary, these IS-54 characterization functions were able to perfectly predict changes in the thermal coagulation zones over time. 3 D thermal coagulation zone models were then obtained based on the rotation operation of IS-54.

3.4. Determining liver tumor treatment schemes using the conformal coverage algorithm

Selection of the treatment scheme in tumor MWA depends on correct localization of the MW antenna and suitable energy applications by KY-2000. The insertion position and direction of the MW antenna is determined by the gravity center of the liver tumor and the orientation vector of the maximum diameter through the gravity center, respectively. The energy applications rely on ablation power and ablation duration.

3.4.1. R-z projection results of physical liver tumors

According to Equations (18)–(20), the surface reconstruction results of liver tumors were processed to achieve the r-z projection, coordinate translation, and coordinate rotation operations. Figure 8 illustrates r-z projection views for cases 1–4, in which the gravity centers of the tumors have been translated to the origin and the maximum diameters through the gravity centers of the liver tumors are aligned with the insertion direction of the MW antenna. The projection results of physical liver tumors were then overlain on r-z plan views of IS-54, as shown in Figure 4. Through the combined 2-D image, we qualitatively observed whether IS-54 completely covered the original tumor for each tumor node in all directions.

Figure 8. R-z projection views of liver tumors. The blue dots indicate the mesh nodes of the tumor surfaces. The gravity centers of the tumors have been placed at the origin.

3.4.2. Mapping results of liver tumors and IS-54

In order to accurately quantify the conformal coverage conditions of liver tumors, 2-D mapping of tumor surface results and IS-54 were performed based on bipolar-angles ${\phi }_O$ and ${\phi }_R,$ as shown in Figure 9, where the red regions represent the backward color maps of IS-54, the green regions represent the forward color maps of IS-54, and the color brightness reflects the lengths of polar radii of IS-54. The color brightness of these maps increased over time due to the increasing polar radii of IS-54. The maximum brightness in forward and backward IS-54 were close to ${\phi }_O$ of π/2 and -π/2, respectively. On the condition of bipolar-angle alignment of IS-54 and liver tumors, the tumor nodes protruding from the IS-54 safety margin are displayed as black dots on the IS-54 maps. As such, tumor maps and IS-54 maps (using case 4 under 60 W as an example) were fused together. IS-54 gradually evolved and increased with time, so the black dots on the fusion maps gradually decreased in number.

Figure 9. Fusion maps of IS-54 and liver tumors. The red regions represent the backward color maps of IS-54, and the green regions represent the forward color maps of IS-54. The black dots indicate the tumor nodes protruding from IS-54 safety margin.

3.4.3. Treatment schemes based on the conformal coverage algorithm

As shown in Figure 10, when the black dots on the fusion maps completely disappeared, conformal coverage was accomplished. At different heating powers, reasonable ablation parameters of liver tumors were derived based on Equation (22)(22) $$\min ({\lambda }_{\mathit{IS}-54}-{\lambda }_{\mathit{tumor}})=5.$$(22) and the SD of errors were calculated based on Equation (23)(23) $$\delta =\sqrt{\frac{1}{n}\sum _{i=1}^n{\left({\lambda }_{\mathit{IS}-54}\right(i)-{\lambda }_{\mathit{tumor}}(i\left)\right)}^2}\hspace{1em}i=1,2,\dots ,n.$$(23) . Finally, the optimal treatment schemes of liver tumors were output according to the minimum SD.

Figure 10. MWA planning interface based on conformal coverage algorithm. The upper-left part of the figure shows the plan projection view of the tumor model and IS-54, the upper-right part shows the fusion map of the tumor model and IS-54, the lower-left part is the selection module of thermal ablation modes and the input module of liver tumor models, and the lower-right part is the output module of optimal treatment parameters.

Because KY-2000 is suitable for the treatment of single tumors within 50 mm, some tumors in this study could not be treated by the ablator and hence the corresponding ablation parameters were not calculated. The specific ablation parameters of cases 1–15 are shown in Table 4.

Table 4. Ablation parameters of liver tumors.

Case No.	Ablation power (W)	Ablation duration (s)	Minimum SD (mm)
1	60	176	5.059
2	N/A	N/A	N/A
3	60	567	5.035
4	60	783	7.933
5	50	1089	6.851
6	55	1155	8.355
7	N/A	N/A	N/A
8	55	573	7.377
9	45	819	8.468
10	N/A	N/A	N/A
11	40	367	5.508
12	50	596	6.444
13	45	326	5.091
14	40	123	5.493
15	60	1012	7.151

N/A: not applicable.

In addition, details of the insertion position and direction of the MW antenna are listed in Table 1. The cross-sectional views of cases 1–5 are illustrated in Figure 11, which indicate that the MWA antenna was inserted to the gravity center positions and was aligned with the maximum diameters through the gravity centers.

Figure 11. The cross-sectional views of liver tumors and MW antenna. The blue lines indicate the insertion positions and orientations of the MW antenna.

4. Discussion

For thermal ablation treatments, patients with a safety margin of greater than 5 mm have a lower probability of developing LTP [25]. For a highly precise estimate, specific IS-54 models of ex vivo porcine livers were established in this study and CT images of patients were used to derive quantitative data of tumors. 2-D projection images and fusion maps of physical tumors and IS-54 were then obtained. Projection images enabled a qualitative comparison to assist in determining the approximate ablation duration and reducing the computational loads. Fusion maps provided a quantitative comparison of physical tumors and IS-54. The major contribution of this study is the new conformal coverage method based on these fusion maps that can determine the treatment parameters in real time and enhance the effectiveness of clinical MWA. Moreover, point-to-point rigid registration of tumor data and IS-54 is not necessary.

An advantage of this bipolar-angle mapping method is that it can achieve conformal ablations of all tumors. In previous studies using combined 3-D images of physical tumors and ablation zones [25,26], not all cases could achieve a 5-mm safety margin. Some conformal coverage methods have been proposed based on comparisons of specific lengths [27,28]; thus, the precision requires improvement. The largest 3-D sphere fitting thermal coagulation zone has often been used for planning purposes [29]. Moreover, some tumor maps have been derived to achieve parameterization of the tumor surface, including spherical map [30] and cylindrical map [22]. To reduce the high computational load, a tumor pseudo-cylindrical map has been proposed for reliable therapy assessment and as an interactive navigation tool; however, this map is not suitable for concave tumors and deskewing is compulsory [22].

Based on the proposed conformal coverage algorithm, it is easy to obtain the optimal MWA treatment parameters when the maximum diameters of the liver tumors are much smaller than the treatment limit of KY-2000 (<50 mm); for example, the 1st, 13th and 14th liver tumors shown in Table 1. However, if the size of the liver tumor is close to 50 mm, then conformal treatment is difficult to achieve. This is due to the considerations of the 5-mm safety margin and the specific morphology of the liver tumors (e.g. the 5th, 6th, 9th and 15th cases). When the size of the liver tumors is larger than 50 mm, the effective treatment parameters cannot be obtained (e.g. the 2nd, 7th and 10th cases). In addition, if the shapes of the liver tumors are similar to that of IS-54 (a pear-like shape), there are smaller ablation errors; i.e. less damage to normal tissue (e.g. the 1st, 3rd, and 13th cases). For some irregular and larger tumors, MWA conformal treatments can be achieved by applying higher powers [31], infusing saline into the tissue surrounding the MW antenna [32], or using multiple applicators [33] to increase the ablation zone sizes.

In this study, the coagulation zone sizes were measured solely on the 54 °C isotherm, because the isotherm could be simply and quickly derived. In order to verify the validity of the temperature threshold method, IS-54 was compared with existing coagulation zone models based on thermal equivalent doses, including the Arrhenius model and CEM₄₃ model [34–36].

The Arrhenius model is the most widely used method for evaluating thermal coagulation zones, which is given as follows: (26) $$\Omega (\tau )={\int }_0^{\tau }A\cdot \hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{\Delta E}{R\cdot T\left(t\right)}\right)dt;$$(26) (27) $${\theta }_d=1-\hspace{0.17em}\text{exp}\hspace{0.17em}(-\Omega (\tau \left)\right),$$(27) where ${\theta }_d$ is the tissue death fraction (%), $\Omega (\tau )$ is the tissue death index, $\tau$ is time (s), $A$ is the frequency factor (s⁻¹), $\Delta E$ is the activation energy for the irreversible coagulation reaction (J/mol), $R$ is the universal gas constant (8.314 J/mol/K), and $T(t)$ is the absolute temperature (K). For healthy liver tissue,$A=7.390\times {10}^{39}$ s⁻¹, $\Delta E=2.577\times {10}^5$ J/mol/K [37]. When $\Omega (\tau )=1$ (i.e. ${\theta }_d=63.212\%$), the biological tissue is considered to have undergone complete thermal coagulation.

Figure 12 shows a comparison between the IS-54 model and the Arrhenius thermal damage model for different MW durations. The red region represents the thermal coagulation zone, the black line represents the IS-54 model and the white line represents the Arrhenius thermal damage model. The Arrhenius thermal damage model has a non-smooth outline and is slightly smaller than IS-54. However, with time, the Arrhenius thermal damage model becomes approximately the same as IS-54; thus, the shape and size changes of the Arrhenius thermal damage model are similar to those of IS-54 for different MWA durations. Therefore, the modeling method of IS-54 is also suitable for the Arrhenius thermal damage model. This further verifies the effectiveness of the proposed method.

Figure 12. A comparison between the IS-54 and the Arrhenius thermal damage model at (a) 206 s, (b) 302 s, (c) 408 s, (d) 466 s, (e) 523 s, and (f) 552 s. The red region represents the thermal coagulation zone, the black line represents the IS-54 model, and the white line represents the Arrhenius model.

It is hypothesized that the CEM₄₃ model is effective for the in vivo MWA condition. To this end, IS-54 was compared with the CEM₄₃ model, which is defined as: (28) $$\mathit{CE}{M}_{43}={\int }_0^{\tau }{\left[R_{\mathit{CEM}}\right]}^{43-T}dt,$$(28) where CEM₄₃ is the cumulative equivalent minutes at 43 °C (min); $R_{\mathit{CEM}}$ is a dimensionless factor, and is set to 0.25 and 0.5 for temperatures below and above 43 °C, respectively; $\tau$ is time (s); and $T$ is temperature (°C).

The lethal thermal dose threshold for livers is 240 CEM₄₃; i.e., CEM₄₃=240 min [38]. Figure 13 compares the IS-54 model and the 240 CEM₄₃ thermal damage model under different MWA durations. The aurantius region represents the thermal coagulation zone, the black line represents the IS-54 model, and the white line represents the 240 CEM₄₃ thermal damage model. IS-54 is significantly larger than the 240 CEM₄₃ thermal damage model, which is attributable to the fact that actual blood perfusion of liver is not considered in the FEM simulation models. Therefore, in future work, the IS-54 model based on FEM simulation data of in vivo livers should be further derived for an effective comparison with the 240 CEM₄₃ model. Despite this, the 240 CEM₄₃ thermal damage model also exhibits time-evolution characteristics. Therefore, in the evaluation of the 240 CEM₄₃ thermal damage model, the characteristic length growth model and the shape variation factor can also be used to characterize the shape and size changes of the thermal coagulation zones, although the modeling forms may be different from IS-54. In addition, compared with the 240 CEM₄₃ thermal damage model, the IS-54 model is relatively simple and requires a significantly shorter calculation time. Hence, the IS-54 model was employed in the proposed conformal coverage algorithm.

Figure 13. A comparison between the IS-54 and the 240 CEM₄₃ thermal damage model at (a) 206 s, (b) 302 s, (c) 408 s, (d) 466 s, (e) 523 s, and (f) 552 s. The aurantius region represents the thermal coagulation zone, the black line represents the IS-54 model, and the white line represents the 240 CEM₄₃ model.

This study only considered the homogeneous tissue model instead of healthy and tumorous tissue models, which enables a direct comparison with existing experimental studies. If the tumor is located very near to the liver boundary or the blood vessel, a geometrically-accurate surface model of the liver or blood vessel should be reconstructed and truncated around the tumor to reduce the high computational demand [39]. The geometrically-truncated model can then be mapped on the bipolar-angle map of the tumor and IS-54. In order to ensure the liver boundary or the blood vessel is not ablated, the polar radius of IS-54 must be less than the corresponding polar radius of the liver boundary or blood vessel. As such, the liver capsule and the blood vessel will not be destroyed according to the proposed conformal coverage algorithm.

Nevertheless, this study has several limitations. Firstly, post-interventional CT images are not acquired and, therefore, a comparison to post-interventional data cannot be performed. Secondly, the physical constraints of antenna insertion are not considered [40–42]. Thirdly, the characterization models of IS-54 were derived based on ex vivo porcine livers and the proposed correlation might change if there is a large blood vessel in close proximity to the target. Heat sinks caused by large blood vessels during treatment would lead to incomplete thermal coagulation [43,44]. In future work, an FEM simulation model including large vessels can be established to derive the growth models of specific characteristic lengths (e.g. the distance from the MW antenna to large vessels) and the functional expressions of thermal coagulation zones. In addition, tumor deformation during microwave ablation was not considered, which is another limitation of this study. The IS-54 model would, therefore, overestimate the actual treated tissue and the safety margin would be larger than the fixed 5 mm due to tissue contraction [45,46]. However, the tissue contraction during MWA could be evaluated using CT images [47] or measured using small specimens of ex vivo liver [48]. The precise evaluation results were then employed to amend the characteristic length growth models and shape variation factors of IS-54. In addition, by changing the objective function (Equation (22)(22) $$\min ({\lambda }_{\mathit{IS}-54}-{\lambda }_{\mathit{tumor}})=5.$$(22) ), the proposed conformal coverage method can be extended to cover different conditions, not a fixed 5 mm pretreatment margin.

5. Conclusion

This study presents a conformal coverage method for liver tumors during MWA. Using 2-D bipolar-angle mapping of IS-54 models and tumor surface data, MWA treatment schemes can be quickly and precisely formulated. Therefore, the proposed method can be applied to computer-assisted MWA treatment planning of liver tumors. In future study, in vivo properties of tumor tissues will be considered to obtain more accurate MWA treatment parameters. Furthermore, the effectiveness of this conformal coverage method should be validated through a comparison with clinical MWA results of liver tumors.

Acknowledgements

The authors would like to thank the anonymous reviewers for their constructive comments and suggestions.

Disclosure statement

The authors declare that they have no conflict of interest.

Additional information

Funding

This work was funded by Beijing Natural Science Foundation [Grant No. 7174279] and National Natural Science Foundation of China [Grant No. 71661167001, 61801312 and 61871005].