Radiofrequency ablation for liver tumors abutting complex blood vessel structures: treatment protocol optimization using response surface method and computer modeling

Abstract

Objective

To achieve a result of a large tumor ablation volume with minimal thermal damage to the surrounding blood vessels by designing a few clinically-adjustable operating parameters in radiofrequency ablation (RFA) for liver tumors abutting complex vascular structures.

Methods

Response surface method (RSM) was employed to correlate the ablated tumor volume ($R_a$) and thermal damage to blood vessels ($D_t$) based on RFA operating parameters: ablation time, electrode position, and insertion angle. A coupled electric-thermal-fluid RFA computer model was created as the testbed for RSM to simulate RFA process. Then, an optimal RFA protocol for the two conflicting goals, namely (1) large tumor ablation and (2) small thermal damage to the surrounding blood vessels, has been achieved under a specific ablation environment.

Results

Linear regression analysis confirmed that the RFA protocol significantly affected $R_a$ and $D_t$ (the adjusted coefficient of determination $R_{\mathit{adj}}^2$ = 93.61% and 95.03%, respectively). For a proposed liver tumor scenario (liver tumor with a dimension of 4$\times$3$\times$2.9 cm³ abutting a complex vascular structure), an optimized RFA protocol was found based on the regression results in RSM. Compared with a reference RFA protocol, in which the electrode was centered in the tumor with a 12-min ablation time, the optimized RFA protocol has increased $R_a$ from 98.1% to 99.6% and decreased $D_t$ from 4.1% to 0.4%, achieving nearly the complete ablation of proposed liver tumor and ignorable thermal damages to vessels.

Conclusion

This work showed that it is possible to design a few clinically-adjustable operating parameters of RFA for achieving a large tumor ablation volume while minimizing thermal damage to the surrounding blood vessels.

Keywords:

• Radiofrequency ablation
• liver tumor
• blood vessels
• heat sink
• treatment protocol optimization

1. Introduction

Radiofrequency ablation (RFA) is one of the energy-based minimally invasive tumor ablation modalities [1,2]. RFA has been widely used in the treatment of various tumors with favorable clinical outcomes [3–8]. However, when treatment is applied to tumors abutting to large blood vessels, the heat generated by the electrode is dissipated by those large vessels (‘heat-sink’ effect), leading to incomplete tumor ablation [9–11]. The incomplete ablation usually leads to the treatment failure and local recurrence [12,13]. Currently, there are no effective, safe, and easy ways available clinically to eliminate or minimize the ‘heat-sink’ effect of large vessels when RFA is used. It is still a challenge for RFA to achieve the complete ablation of liver tumors abutting large vessels and minimal thermal damages to surrounding vessels. Some methods attempted clinically include a temporary interruption of the blood flow (i.e., Intermittent Pringle maneuver [14] or Balloon occlusion [15]) and the development of new ablation electrodes [16–19]. However, due to cooling protection of blood fluid, the method of vessel occlusion could usually cause thermal damages to vessels, leading to severe post-ablation complications, such as hemorrhage. Mathematical modeling and simulation have been applied widely in both the design and development of novel clinical protocols and the optimization and improvement of existing clinical protocols [20,21]. A comprehensive and validated mathematical model can provide reliable information not only on the treatment performance, but the spatiotemporal temperature and electric field distributions in the treatment area, which are difficult to attain in experiments [19,22]. Shao et al. [19,23] used mathematical models to investigate the thermal ablation of para-macrovascular tumors and the influence of the ‘heat-sink’ effect on the ablation zone. They designed a novel ablation electrode and a specific ablation protocol to reduce the ‘heat-sink’ effect of large blood vessels. They concluded that the distance between the RFA electrode and the blood vessel plays a major role in the incomplete ablation during the heating process and a new bipolar electrode matrix system can achieve a larger tissue ablation zone while arriving at a lower peak temperature outcome. Huang et al. [24] established a finite difference method to calculate the ablation zone formation via RFA and investigated the influence of the location, diameter, and orientation of the blood vessel to the electrode on the ablation zone formation. They found that a vessel that is orthogonal to the RF electrode has less cooling impact in the ablation zone than a parallel vessel. González-Suárez et al. [25] established a computer model to simulate in vivo RFA experiments in which the liver surface was heated by an internally cooled electrode placed as close as possible to a vessel. The computational findings verified that the heat-sink effect could protect the portal vein wall.

It is worth mentioning that the computer models currently used in literature only considered a simplified vascular system within homogenous liver tissues [19,21,24,25]. However, a growing body of evidence suggests that a liver tumor with different thermal and electrical properties should be included in the computer model. Furthermore, given that the liver is very vascular, the model should include a liver tumor abutting to a major vascular structure. Moreover, the size of the tumor and the distance between the tumor and the large blood vessels are also required to be considered in the model as they can affect the complete ablation of tumors and the thermal damage to the vessels. A model incorporating all these elements would be more representative of the treatment planning of RFA.

This paper aims to combine an engineering optimization method, response surface method (RSM) with computer modeling to enhance RFA treatments for a liver tumor scenario mentioned-above. Further, RFA treatment protocols with the goal of completely ablating the liver tumor while minimizing the thermal damage to the vessels were determined by optimizing only operating parameters that can be adjustable by radiologists or clinician during treatment in clinic. These parameters included the ablation time, the electrode position, and the angle of RF electrode. By considering the real liver anatomy, a 3 D RFA liver tumor model was built and placed in proximity to a complex vascular system including the hepatic vein, portal vein, and hepatic artery. In the such tumor position, the heat-sink effect of vessels has a significant effect on treatment outcome. The heat dissipation was modeled by using the Navier-Stokes equation combined with the Pennes bioheat transfer equation. The proposed model was then verified by comparing the tissue impedance change, ablation zone, and tissue temperature between the current model and ex vivo experiments in literature.

2. Materials and methods

2.1. Geometrical modeling

As shown in Figure 1(a), a 3 D liver tumor with a complex vascular system was created in a cubic liver tissue with a dimension of 120$\times$120$\times$130 mm³. The dimension of the liver tissue was determined by a convergence test with a criterion of less 0.1% change in the ablated volume between two contiguously increasing dimensions of liver tissue domain. An ellipsoidal liver tumor with dimensions of 40, 30, and 29 mm along the X-, Y- and Z-axis, respectively was generated within the liver tissue. Three tubular structures (i.e., hepatic vein, portal vein, and hepatic artery) in an anatomical structure were also modeled around the liver tumor. A commercial cool-tip RFA electrode (Cool-tip^TM RF Ablation System E Series, Medtronic, Minneapolis, MN, USA) with a diameter of 15 mm and the active tip of 30 mm was inserted into the liver tumor. The inside of the electrode was cooled with chilled saline at the temperature of 10 $°C.$ The relative distribution of the elements in the 3 D model is shown in Figure 1(b).

Figure 1. The geometric model of RFA for a liver tumor abutting a vascular structure: (a) a 3 D view and (b) X-Z cross-sectional view (with the variables ($x,z,\mathrm{and}\mathrm{ }\theta$) of electrode placement; dimension in mm and out of scale).

The diameters of the main hepatic vein and branch veins were taken as 9 and 7.2 mm, respectively, consistent with measurements from human liver anatomy [26] (as illustrated in Table 1). The velocity of blood flow in the hepatic vein was set as 39.5 cm/s. Similarly, the diameter of the portal vein and hepatic artery were 10 and 4 mm, respectively, and the flow velocities of the portal vein and hepatic artery were 19.8 and 44.3 cm/s, respectively. The thickness of all vessel walls was assumed set at 0.5 mm [27–29]. The direction of blood flow in the three vascular structures are shown in the Figure 1(a). Note that the flow directions of the portal vein and the hepatic artery are opposite.

Table 1. The velocity and diameter of three vessels in the hepatic vascular system [27].

Vessels	Velocity (cm/s)	Diameter (mm)
Portal vein	19.8 ± 5.1	8–12
Hepatic vein	39.5 ± 12.2	5–9
Hepatic artery	44.3 ± 15.8	2–6

2.2. Coupled electric-thermal-fluid modeling of RFA

A commercial finite element software (COMSOL Multiphysics 5.5, COMSOL inc., Burlington, MA, USA) was used to simulate RFA procedure of the coupled multi-physics. During ablation process, RF power was applied to the electrode, and the Laplace equation was used to describe the electric field generated in the tissue of interest: (1) $$\begin{array}{c}\begin{array}{c}\nabla ·\sigma \nabla V=0\end{array}\end{array}$$(1) where $\sigma$ is the temperature-dependent tissue electrical conductivity and $V$ is the applied voltage. Electric field intensity ($E$) and current density ($J$) can be calculated by (2) $$\begin{array}{c}\begin{array}{c}E=-\nabla V\end{array}\end{array}$$(2) (3) $$\begin{array}{c}\begin{array}{c}J=E·\sigma \end{array}\end{array}$$(3)

The local power density that causes the tissue heating is the product of the current density and the electric field strength as follows (4) $$\begin{array}{c}\begin{array}{c}{Q}_{hs}=J\cdot E={\sigma \left|\nabla V\right|}^2\end{array}\end{array}$$(4)

The Pennes bioheat transfer Equation [29] was used to govern the heat transfer in the tissues, as shown in Equation (5). Meanwhile, for the blood vessel, the fluid dynamic equation and the heat transfer equation were combined to simulate the ‘heat-sink’ effect of vessels. Given that, the Pennes equation was modified by adding a convection term to the right side of the equation, as given by Equation (6). (5) $$\begin{array}{c}\begin{array}{c}\rho c_p\frac{\partial T}{\partial t}=\nabla ·k\nabla T+{\rho }_b{c}_{bp}{\omega }_b\left(T_b-T\right)+Q_{hs}+Q_m\end{array}\end{array}$$(5) (6) $$\begin{array}{c}\begin{array}{c}\rho c_p\frac{\partial T}{\partial t}=\nabla ·k\nabla T+Q_{hs}\end{array}-{\rho }_b{c}_b{\mathit{u}}_b·\nabla T\end{array}$$(6) where $\rho$ is the density of the tissue, $c_p$ is the specific heat, $T$ is the temperature, $k$ is the thermal conductivity. ${\rho }_b$ is the blood density, $c_{bp}$ is the specific heat of blood, ${\omega }_b$ is the blood perfusion rate, $T_b$ is the temperature of blood entering the tissue, $Q_m$ is the heat generated in the process of tissue metabolism, which is usually negligible due to its small magnitude compared with other items, and $Q_{hs}$ is an external heat power generated by RF energy and it can be solved in Equation (4). In the convection term, ${\mathit{u}}_b$ is the blood flow velocity field in the blood vessel, which is derived by the incompressible Navier–Stokes equation [30], which consists of the momentum conservation equation (i.e., Equation. (7)) and mass conservation equation (i.e., Equation. (8)(8) $$\begin{array}{c}\begin{array}{c}\nabla ·{\mathit{u}}_b=0\end{array}\end{array}$$(8) ) (7) $$\begin{array}{c}\begin{array}{c}{\rho }_b\frac{\partial {\mathit{u}}_b}{\partial t}+{\rho }_b{\mathit{u}}_b·\nabla {\mathit{u}}_b=-\nabla P+\mu {\nabla }^2{\mathit{u}}_b+F\end{array}\end{array}$$(7) (8) $$\begin{array}{c}\begin{array}{c}\nabla ·{\mathit{u}}_b=0\end{array}\end{array}$$(8) where $P$ is the pressure, $\mu$ is the viscosity of the blood 0.0021(kg/m·s), and $\mathit{F}$ is the volume force, which was not considered owing to its relatively small magnitude in this model.

In this study, the water vaporization of liver tissues was also considered as temperature-dependent [30] (9) $$\begin{array}{c}\begin{array}{c}\rho c\frac{\partial T}{\partial t}=\frac{\partial T}{\partial t}\{\begin{array}{c}{\rho }_l{c}_l 0\le T\le 99°C\\ H_{fg}C 99<T\le 100°C\\ {\rho }_g{c}_g T>100°C\end{array}\end{array}\end{array}$$(9) where ${\rho }_i$ and $c_i$ are the density and specific heat of the tissue that change due to the phase change, that is, at a temperature below 100 °C ($i$=$l$ for liquid tissue phase) and a temperature above 100 °C ($i$=$g$ for gaseous tissue phase), $H_{fg}$ is the latent heat, that is, the product of the latent heat of vaporization of water and the water density at 100 °C (2.162 × 109 J/m³), and $C$ is the water ratio inside the liver tissue (68%). The thermal conductivity and electrical conductivity of liver tissues and liver tumor were also considered as temperature dependent in the study [30,31] (10) $$\begin{array}{c}\begin{array}{c}k\left(T\right)=\{\begin{array}{c}{k}_0+\Delta k\left(T-T_{\mathit{ref}}\right) T\le 100°C\\ k_0+\Delta k\left(100°\mathrm{C}-T_{\mathit{ref}}\right) T\le 100°C\end{array}\end{array}\end{array}$$(10) (11) $$\begin{array}{c}\begin{array}{c}\sigma \left(T\right)=\{\begin{array}{c}\sigma \left(37\right)e^{0.015\left(T-37\right)} 0\le T<99°C\\ \sigma \left(37\right)2.5345 99°C\le T\le 100°C\\ \sigma \left(37\right)(2.5345-0.50183(T-100) 100°C<T<105°C\\ \sigma \left(37\right)2.5345\times {10}^{-2} T>105°C\end{array}\end{array}\end{array}$$(11) where $k_0$ and $\sigma \left(37\right)$ are the reference thermal conductivity and electrical conductivity, which are measured at a reference temperature of 37 °C, respectively. The blood perfusion of liver tissue was given by Equation (12)(12) $$\begin{array}{c}\begin{array}{c}{\omega }_b=\{\begin{array}{c}{\mathrm{\omega }}_{b0} T<60°C\\ 0 T\ge 60°C\end{array}\end{array}\end{array}$$(12) (12) $$\begin{array}{c}\begin{array}{c}{\omega }_b=\{\begin{array}{c}{\mathrm{\omega }}_{b0} T<60°C\\ 0 T\ge 60°C\end{array}\end{array}\end{array}$$(12)

To evaluate the ablation zone, the Arrhenius model was used to determine the damaged tissue, which was given by Equation (13)(13) $$\begin{array}{c}\begin{array}{c}\Omega =A{\int }_0^T{e}^{\frac{-\Delta E}{RT\left(\tau \right)}}d\tau \end{array}\end{array}$$(13) [19] (13) $$\begin{array}{c}\begin{array}{c}\Omega =A{\int }_0^T{e}^{\frac{-\Delta E}{RT\left(\tau \right)}}d\tau \end{array}\end{array}$$(13) where $\Omega$ is the degree of tissue death, $A$ is the frequency factor, $\Delta E$ is the activation energy for the irreversible damage reaction, $R$ is the universal gas constant (8.314 J mol⁻¹ K⁻¹), and $T$ is the tissue temperature inside the computational domain at a specified time. For the liver tissue and liver tumor, $A$= 7.390 × 10³⁹s⁻¹ and 3.247 × 10⁴³s⁻¹, respectively, and $\Delta E$= 2.577 × 10⁵ J mol⁻¹and 2.814 × 10⁵ J mol⁻¹, respectively [32]. Regarding to tissue damage, $\Omega \left(t\right)=1$ corresponds to a probability of 63% cell death and $\Omega \left(t\right)=4.6$ corresponds to a probability of 99% cell death. In the study, we defined the $\Omega \left(t\right)=4.6$ as the threshold of tissue ablation for a higher probability of tumor tissue death and $\Omega \left(t\right)=1$ as the threshold of thermal damage to blood vessels for less thermal damage [24].

Except for the electrical conductivity, thermal conductivity, and blood perfusion rate, other parameters of the tissues pertaining to the model were assumed to be constants, summarized in Table 2.

Table 2. Characteristics of the materials used in the theoretical model [32–34].

Modeling element	$\rho$(kg m^-3)	$C$(J kg^-1 K^-1)	$k_0$(W m^-1 K^-1)	$\delta (37)$(S m^-1)	${\omega }_{b0}$(kg m^-3 s^-1)
Liver tissue	1080^a	3455^a	0.515^c	0.203^c	17.03
	370^b	2156^b
Liver tumor	1045^a	3760^a	0.6^c	0.5^c	2.09
	370^b	2156^b
RF electrode	6450	840	18	1.0 × 10^8
Insulated shaft	70	1045	0.026	1.0 × 10^-5
Blood vessel	1000	4180	0.49	0,667

^aat temperature below 100 $°C;$ ^bat temperature above 100 $°C;$ ^cassessed at 37 $°C.$

2.3. Mesh and boundary conditions

A mesh convergence test was performed in the study with the criterion of the difference in the temperature of two points between two contiguous meshes smaller than 0.1%. The two points were chosen at 5 mm away perpendicular to the center of the active tip and 5 mm away along the active tip when electrode is placed at the central position (i.e., $x$ = 65 mm, $z$ = 58.5, $\theta$ = 0$°$). Eventually, the free tetrahedral elements in a total number of 214091 were used in the model.

The initial temperature of the finite element model of RFA was assumed to be uniform and the same as the internal temperature of the human body (14) $$\begin{array}{c}\begin{array}{c}T\left(t=0\right)=37°C\end{array}\end{array}$$(14)

Adopting the power delivery method from Ben-David et al. [35], a time-dependent power was applied on the active tip of the electrode in the study (15) $$\begin{array}{c}P\left(t\right)=P_{\mathit{app}}\left(t\right)\end{array}$$(15) where $P_{\mathit{app}}(t)$ is the applied power and has an initial value of 1.5 W. The transient variation in $P_{\mathit{app}}(t)$ is governed by an impedance-controlled algorithm.

The electrical boundary between the external healthy liver tissue and the electrode cavity was electrically insulated (16) $$\begin{array}{c}\begin{array}{c}n·J=0\end{array}\end{array}$$(16)

All other boundaries were assumed to be of continuity (17) $$\begin{array}{c}\begin{array}{c}n·\left(J_i-J_j\right)=0\end{array}\end{array}$$(17) The thermal boundary condition across the active part of the electrode was (18) $$\begin{array}{c}\begin{array}{c}-n·\left(-k\nabla T\right)=h_c\left(T_c-T\right)\end{array}\end{array}$$(18) where $h_c$ is the convection heat transfer coefficient and $T$ is the temperature of the internally chilled saline. Equation (18)(18) $$\begin{array}{c}\begin{array}{c}-n·\left(-k\nabla T\right)=h_c\left(T_c-T\right)\end{array}\end{array}$$(18) describes the convective heat transfer between the internally chilled saline and the surrounding tissue. In the present study, the flow rate of the internally chilled saline was chosen to be 45 mL min⁻¹, which results in a value of $h_c$ of 4416.3 W m⁻²K⁻¹ [34]. The temperature of the internally chilled saline was chosen to be 10 °C. The selected flow rate and temperature combination was sufficient to maintain the temperature along the surface of the active tip around 10 − 15 $°C$ during RFA. All the other thermal boundary conditions were considered as continuity (19) $$\begin{array}{c}\begin{array}{c}n·\left(k_i\nabla T_i-k_j\nabla T_j\right)=0\end{array}\end{array}$$(19)

Then, all outer surfaces of the liver tissue (i.e., the external boundaries) were set as a Dirichlet boundary condition ($T=37°C$) and a ground terminal (V = 0). Meanwhile, the inlet and outlet boundaries of blood flow were also defined as blood average velocities in the hepatic vascular system, as shown in Table 1, and zero pressure boundary ($P_0=0$ Pa), respectively.

2.4. Impedance-controlled ablation protocol

The RFA protocol adopted in this study was the optimized impedance-controlled ablation protocol proposed by Ben-David et al. [35], which was introduced as a way of avoiding high impedance of the tissues surrounding the electrode as a result of complete water vaporization of the tissue surrounding the RF electrode (i.e., ‘roll-off’). To maximize the efficiency of energy deposition, the generator algorithm first endeavored to achieve a pre-set power based upon the active tip length over 30 s, and then pulsed the power in response to changes in tissue impedance, as shown in Figure 2. While the impedance was stable, power was increased in an increment of 5 W every 30 s. When the impedance increased over a pre-set trigger point (i.e., 120 $\mathrm{\Omega }$), the input power was halted for at least 15 s before it was resumed. When impedance rose so rapidly that a power could not be maintained for at least 10 s, the input power was then reduced by 5 W when energy deposition resumed.

Figure 2. The demonstration of impedance changes along with the energy supply under the impedance-controlled ablation protocol.

2.5. Thermal damage to vessels

The closest edge of thermal damage (i.e., $\Omega \left(t\right)=1$) to inner vessel walls was always in the cross-section ($\mathcal{R}$) perpendicular to the blood vessel and passing through the electrode. The thermal damage of the vessel wall could be evaluated by $D_t,$ which was the ratio of damaged wall area ($S_t$) and the intact wall area ($S$) in the$\mathcal{ }\mathcal{R},$ as shown in Figure 3. (20) $$\begin{array}{c}\begin{array}{c}{D}_t=\frac{S_{\mathrm{t}}}{S}\times 100\%\end{array}\end{array}$$(20)

Figure 3. The quantification method of vessel thermal damage.

2.6. Treatment protocol optimization

In this study, RSM was used to do RFA protocol optimization, which was divided into three main steps: (1) to determine optimizing variables, (2) to design experiments, and (3) to build a statistical model for finding optimal results. As mentioned in Introduction, the purpose of this study is to achieve a larger tumor ablation volume for complete ablation of liver tumors and minimal thermal damages to vessels simultaneously by optimizing a few clinically adjustable RFA operating parameters. Therefore, two targets of treatment optimization could be defined: (1) the ablated ratio ($R_a$), which is the ratio of the ablated tumor volume ($V_{\mathit{tumor}}$) to the whole tumor volume ($V_{\mathit{all}}$); (2) the thermal damage to vessels$,\mathrm{ }{D}_t$ in Equation (20)(20) $$\begin{array}{c}\begin{array}{c}{D}_t=\frac{S_{\mathrm{t}}}{S}\times 100\%\end{array}\end{array}$$(20) . (21) $$\begin{array}{c}\begin{array}{c}{R}_a=\frac{V_{\mathit{tumor}}}{V_{\mathit{all}}}\times 100\%\end{array}\end{array}$$(21)

As shown in Figure 1(b), four operating parameters taken as variables were chosen to be optimized, including the ablation time ($t$), insertion position of electrode along X-axis ($x$), Z-axis ($z$), and the angle ($\theta$) [23–25]. The adjustable range of the ablation time (i.e., 660–780 s) was determined by considering a 12-min RFA protocol for a live tumor of 3 cm diameter in clinic [35]. Therefore, the insertion of RF electrode was determined by the three parameters (i.e., $x,$ $z,$ and $\theta ,$ as shown in Figure 1(b)). Typically, the electrode was centered in the tumor for a best match between the predicated ablation zone and tumor, i.e., $x$ = 65 mm, $z$ = 58.5 mm, and $\theta$ = 0° in the study [36]. Therefore, we chose this central position of electrode and the clinical ablation time (i.e., 12 min) as a reference protocol. Based on this central position, the range of the $x,$ $z$ and $\theta$ for optimization were adopted as 60 − 70 mm, 53.5 − 63.5 mm, and −30° − 30°, respectively. In general, this multi-objective treatment optimization problem can be defined as (22) $$\begin{array}{c}\{\begin{array}{c}\begin{array}{l}V\mathit{ariable} t,x,z,\theta \\ E\mathit{valuate} R_a=f\left(t,x,z,\theta \right) \mathit{and} D_t=f\left(t,x,z,\theta \right) \\ M\mathit{aximize} R_a \\ M\mathit{inimize} D_t \\ S\mathit{ubject} to: \\ 660 s\le t\le 780 s\\ 60 mm\le x\le 70 mm\\ 53.5 mm\le z\le 63.5 mm\\ -30°\le \theta \le 30°\end{array}\end{array}\end{array}$$(22)

2.7. Statistical analysis

A commercial statistical software of Minitab (Version 20, Minitab Inc., PA, USA) was used in the study to perform the statistical modeling of RSM, of which Box-behnken method in RSM was chosen to design the experiment matrix for simulations. Based on the simulation results, a statistical predictive model was built by using linear regression. The adjusted coefficient of determination ($R_{\mathit{adj}}^2$) was used to evaluate the accuracy of model, and each variable in the model were considered as statistical significance at p < 0.05.

3. Results

3.1. Validation of the RFA model

The validation of the RFA model was performed by comparing the temperature results of the simulation and an ex vivo experiment in literature [27]. The geometry and initial conditions of RFA model was adjusted to be consistent with the conditions in the ex vivo experiment. The coolant volume flow rate inside the electrode was set at 25 ml/min and the blood perfusion was ignored in the $ex \mathit{vivo}$ liver. A good agreement in the temperature at the end of RFA treatment (i.e., 720 s) between the simulation and experiment results was recognized at 5, 10, 15, and 20 mm from the center of RF electrode, as shown in Figure 4 and Table 3. Meanwhile, similar ablation diameters also can be found with a minor difference of 2.83% (4.47 mm in the simulation vs. 4.6 ± 0.2 mm in the experiment). A similar acceptable accuracy of the proposed model was also validated by comparing the temperatures at three locations (10, 15, and 30 mm from the center of RF electrode) between the computer model and an$\mathit{\text{ ex }}\mathit{vivo}$ experiment study under a constant power delivery method (i.e., 50 V) for 200 s [37], as illustrated in Table 3. All these comparisons support the accuracy of the proposed RFA model as a reasonable testbed for treatment optimization.

Figure 4. (a) the temperature change at 5, 10, 15, and 20 mm from the center of RF electrode, and (b) the temperature distribution at the ending of RFA treatment.

Table 3. Comparison of the temperature between the proposed RFA computer model and the ex vivo experiments.

Location (mm)	Temperature (°C)		Temperature (°C)
	Simulation	Experiment [27]	Simulation	Experiment* [37]
5	98.9	100.4 ± 1.4	–	–
10	78.3	79.5 ± 4.7	53.5	52 ± 4
15	59.6	57.3 ± 4.2	37.0	33 ± 2
20	45.5	46.2 ± 4.4	30.1	27 ± 1

*An approximate data derived from Figure 4 in [37].

3.2. RSM results

Table 4 demonstrates the designed simulation matrix including 25 cases corresponding with the responses (i.e., $R_a$ and $D_t$) for four variables (i.e., $t,$ $x,$ $z,$ and $\theta$), and extra data on the simulation results were provided in Table A1 of Appendix. Meanwhile, based on the simulated results, Table 5 shows the linear regression results using a full quadratic polynomial model presented in Equation (23). (23) $$Y={\beta }_0+{\beta }_1{t}_n+{\beta }_2{x}_n+{\beta }_3{z}_n+{\beta }_4{\theta }_n+{\beta }_{12}{t}_n{x}_n+{\beta }_{13}{t}_n{z}_n+{\beta }_{14}{t}_n{\theta }_n+{\beta }_{23}{x}_n{z}_n+{\beta }_{24}{x}_n{\theta }_n+{\beta }_{34}{z}_n{\theta }_n+{\beta }_{11}{t_n}^2+{\beta }_{22}{x_n}^2+{\beta }_{33}{z_n}^2+{\beta }_{44}{{\theta }_n}^2$$(23) where $Y$ represents the responses (i.e., $R_a$ or $D_t$), $t_n,$ $x_n,$ $z_n,$ and ${\theta }_n$ are rescaled variables of $t,$ $x,$ $z,$ and $\theta$ using min-max normalization method. Generally, two regression results showed a significant influence of four variables on both two responses (both p < 0.001). The influence of linear, square and two-way interaction were considered for each variable in the regression model. The ablation time ($t$) did not present a significant effect on both $R_a$ and $D_t$ in the scope of 660–780 s. Otherwise, other variables has showed the significance on the response of $R_a$ or $D_t.$ Especially, the significances of $z$ and $\theta$ were presented in the square part of the regression models for $R_a$ (p = 0.003 and < 0.001, respectively) and $D_t$ (p = 0.009 and < 0.001 respectively). Meanwhile, only $z$ and $\theta$ has a significant effect in the way of two-way interaction (p < 0.001 for $R_a,$ and p = 0.043 for $D_t$). In the linear part, the significant effects of $x,$ $z,$ and $\theta$ are only presented in the response of $D_t$ with p < 0.001, 0.039, and < 0.001, respectively.

Table 4. Results of ablation and thermal damage for each designed simulation.

No.	Variables				Responses
	t (s)	x (mm)	z (mm)	$\theta$ (°)	$R_a$ (%)	$D_t$ (%)
1	720	60	53.5	0	86.0%	5.2%
2	660	65	53.5	0	88.6%	6.1%
3	720	65	53.5	30	90.3%	5.5%
4	720	65	53.5	−30	39.5%	27.3%
5	780	65	53.5	0	91.6%	8.4%
6	720	70	53.5	0	84.8%	15.0%
7	660	60	58.5	0	97.5%	0.0%
8	720	60	58.5	30	70.8%	0.0%
9	720	60	58.5	−30	70.4%	11.4%
10	780	60	58.5	0	99.3%	0.3%
11	660	65	58.5	30	65.4%	4.2%
12	660	65	58.5	−30	64.2%	21.5%
13	720	65	58.5	0	98.1%	4.1%
14	780	65	58.5	30	67.8%	4.2%
15	780	65	58.5	−30	66.2%	22.2%
16	660	70	58.5	0	96.7%	10.8%
17	720	70	58.5	30	60.7%	14.2%
18	720	70	58.5	−30	60.4%	30.6%
19	780	70	58.5	0	98.3%	11.6%
20	720	60	63.5	0	84.2%	0.0%
21	660	65	63.5	0	89.6%	7.1%
22	720	65	63.5	−30	42.1%	6.5%
23	720	65	63.5	30	90.3%	19.9%
24	780	65	63.5	0	92.5%	7.1%
25	720	70	63.5	0	83.5%	11.9%

Table 5. The value of coefficient in RSM regression model.

	$R_a$	p Value	$D_t$	p Value
${\beta }_0$	0.342	<0.001	0.203	0.003
${\beta }_1$	0.030	0.385	0.002	0.530
${\beta }_2$	0.093	0.143	0.086	<0.001
${\beta }_3$	0.777	0.928	−0.166	0.039
${\beta }_4$	1.660	0.695	−0.516	<0.001
${\beta }_{12}$	−0.002	0.982	0.004	0.910
${\beta }_{13}$	−0.001	0.991	−0.023	0.543
${\beta }_{14}$	0.004	0.964	−0.007	0.854
${\beta }_{23}$	0.005	0.956	0.020	0.594
${\beta }_{24}$	−0.001	0.991	−0.050	0.203
${\beta }_{34}$	−0.990	<0.001	0.084	0.043
${\beta }_{11}$	0.008	0.918	0.018	0.587
${\beta }_{22}$	−0.134	0.105	0.055	0.111
${\beta }_{33}$	−0.282	0.003	0.100	0.009
${\beta }_{44}$	−1.156	<0.001	0.338	<0.001
$R_{\mathit{adj}}^2$	93.61%		95.03%

3.3. Multi-objective optimization

As shown in Figure 5, multi-objective optimization was performed (i.e., maximal $R_a$ and minimal $D_t$) and the change tendency of $R_a$ and $D_t$ along with the variables (i.e., $t,$ $x,$ $z,$ and $\theta$) were found. A monotonical increase relationship can be observed in the two pairs of $t-R_a$ and $x-D_t,$ and however, a non-monotonic change was found in other pairs. For achieving a maximal $R_a$ and minimal $D_t,$ an optimal set for all variables was predicted and set as $t$ = 780 s, $x$ = 62 mm, $z$ = 57.8 mm, and $\theta$ = 3.6°. The predictive optimal results of $R_a$ and $D_t$ are 99.0% with a 95% Confidence Interval (CI) of 93.4 − 100% and 0.1% with a 95% CI of 0 − 3.4%, respectively. Correspondingly, the computer model results showed in agreement with the predictive results on $R_a$ and $D_t$ using the optimal set, which are 99.6% and 0.4%, respectively, as shown in Figure 6(a–c). The results of the reference protocol are illustrated in Figure 6(d–f) with $R_a$ = 98.1% and $D_t$ = 4.1%. The demonstrations of the computer model results in other views were also provided in Figures A1 and A2 of Appendix.

Figure 5. The multi-objective optimization results of RSM (the red solid line and blue dashed line represent the optimal solutions and predictive results of the regression models, respectively).

Figure 6. The RFA treatment results in X-Z cross-sectional view: (a–c) ablation, thermal damage, and temperature distribution, respectively for the optimized protocol; (d–f) ablation, thermal damage, and temperature distribution, respectively for the reference protocol. (The black line is the contour line of tumor).

4. Discussion

Using a coupled electric-thermal-fluid RFA model, we demonstrated the feasibility of optimizing a few clinically adjustable RFA operating parameters for a liver tumor in an anatomically sensitive location (tumor abutting a complex vascular structure). As shown in Figure 5, although a similar change tendency of $R_a$ along with four variables could be observed for $D_t,$ it is difficult to find a consistent variable set for both $R_a$ and $D_t$ simultaneously. A possible explanation is the nonuniform distribution of blood vessels in liver compared to the tumor zone, which caused the following two issues for RFA treatment. First, the thermal distribution during RFA treatment was distorted by the heat-sink effect of blood vessels, as shown the temperature distribution in Figure 6. Thus, it is difficult RFA treatment to generate a uniform ablation zone that matches the tumor shape (i.e., ellipsoid in this study). Second, avoiding thermal damages to the blood vessels is another constraint. Thus, it is difficult for radiologists to find an appropriate RFA treatment protocol, which could ablate the tumor completely and avoid thermal damages to blood vessels for the tumor scenario presented in this study. Instead, what is needed is an optimized protocol plan before RFA treatment. Such a personalized plan is necessary to achieve a more efficient and safe treatment rather than using a constant protocol (such as, the reference protocol).

We believe that the calculation of an optimized protocol could be pursued for each patient, especially for those tumors abutting to major vascular structure. In this study, maximizing the completeness of tumor ablation ($R_a$) while minimizing thermal damage to vessels ($D_t$) have been simultaneously achieved by optimizing four clinically adjustable RFA operating parameters, including the ablation time (i.e., $t$) and placement of RF electrode (i.e., $x,$ $z,$ and $\theta$) using an engineering optimization method (i.e., RSM). Compared with the reference protocol, the optimized protocol achieved a better simulation result, of which $R_a$ increased from 98.1% to 99.6% and $D_t$ decreased from 4.1% to 0.4%. As shown in Figure A3 of Appendix, the optimized protocol successfully avoided the thermal damage penetrating the hepatic vein wall while increasing the ablated tumor volume (See Figure A1 of Appendix). Thus, it is possible for such a liver tumor abutting complex blood vessel structure to achieve a complete tumor ablation without thermal damages to blood vessels by the proposed RFA protocol optimization method, i.e., RSM. The benefit to the patient is a minimally invasive procedure that yields a lower risk of recurrence, while minimizing complications.

In addition, a regression model has been established to predict the relationship between variables and their responses. The optimized protocol agreed nicely with the predicted simulation results. In this study, three conclusions can be made from the RSM regression model (Table 5 and Figure 5). They are that 1) ablation time ($t$) in such range is not a significant variable affecting $R_a$ and $D_t.$ Two possible reasons are that the thermal equilibrium was reached at the end of RFA treatment [32], and the 2-min range of $t$ chosen in the study is not long enough to introduce a significant influence on ablation results; 2) between the electrode and blood vessel, two distances ($x$ and $z$) were considered of which $x$ represents the distance along the electrode direction and $z$ represents the distance perpendicular to electrode. A more significant effect for the treatment outcome could be observed in $z$ rather than $x.$ In other words, RFA treatment is more sensitive to the blood vessel on the side of the electrode compared with it at the bottom of electrode; and 3) both $R_a$ and $D_t$ are very sensitive to insertion angle ($\theta$) by which it dominates the direction of ablation zone. It is worth mentioning that the insertion of electrode along with the long axis of tumor is not always the best choice, especially for tumors surrounding the blood vessels, like the one in this study.

However, there are still two limitations with this study needed to be overcome in the future. First, the vascular structure in the study was built based on data from the literature. An image-based 3 D model should provide a more accurate and patient-specific vascular structure. Second, only simulation was performed to validate the efficiency of engineering optimization on enhancing RFA treatment outcome. Future studies of in vivo experiments are warranted for the application of engineering optimization in clinical RFA treatment.

5. Conclusion

The combination of computer modeling and engineering optimization technique can determine a few operating parameters in the process of the ablation of liver tumors abutting complex blood vessel structures. A large tumor ablation with minimal thermal damages to surrounding blood vessels can be achieved simultaneously by using the optimized RFA treatment protocol. The proposed method and finding can contribute to the application of RFA in the treatment of liver tumors in critical anatomical locations.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This article was supported by the National Natural Science Foundation of China through a funding grant (Grant No. 81801795) to the corresponding author (BZ), RUH Foundation Research Grant [Grant No. R210211] to WZ and the Program of Production, Research and Development of Minhang District of Shanghai through a funding grant [Grant No. 2019MHC066] to the corresponding author (ZQ).