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 ( ) and thermal damage to blood vessels ( ) 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 and (the adjusted coefficient of determination = 93.61% and 95.03%, respectively). For a proposed liver tumor scenario (liver tumor with a dimension of 4 3 2.9 cm3 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 from 98.1% to 99.6% and decreased 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.


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][4][5][6][7][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][10][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 'heatsink' 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][17][18][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 'heatsink' 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 alez-Su arez 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.

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Â120Â130 mm 3 . 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 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][28][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.

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: where r 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 The local power density that causes the tissue heating is the product of the current density and the electric field strength as follows 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).
where q is the density of the tissue, c p is the specific heat, T is the temperature, k is the thermal conductivity. q b is the blood density, c bp is the specific heat of blood, x 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, 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)) where P is the pressure, l is the viscosity of the blood 0.0021(kg/mÁs), and 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] qc where q 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 3 ), 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] where k 0 and r 37 ð Þ 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) To evaluate the ablation zone, the Arrhenius model was used to determine the damaged tissue, which was given by Equation (13)  where X is the degree of tissue death, A is the frequency factor, DE is the activation energy for the irreversible damage reaction, R is the universal gas constant (8.314 J mol À1 K À1 ), 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 39 s À1 and 3.247 Â 10 43 s À1 , respectively, and DE¼ 2.577 Â 10 5 J mol À1 and 2.814 Â 10 5 J mol À1 , respectively [32].
Regarding to tissue damage, X t ð Þ ¼ 1 corresponds to a probability of 63% cell death and X t ð Þ ¼ 4:6 corresponds to a probability of 99% cell death. In the study, we defined the X t ð Þ ¼ 4:6 as the threshold of tissue ablation for a higher probability of tumor tissue death and X t ð Þ ¼ 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.

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, h ¼ 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 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 where P app ðtÞ is the applied power and has an initial value of 1.5 W. The transient variation in P 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 All other boundaries were assumed to be of continuity The thermal boundary condition across the active part of the electrode was where h c is the convection heat transfer coefficient and T is the temperature of the internally chilled saline. Equation (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 À1 , which results in a value of h c of 4416.3 W m À2 K À1 [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 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.

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 X), 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.

Thermal damage to vessels
The closest edge of thermal damage (i.e., X t ð Þ ¼ 1) to inner vessel walls was always in the cross-section (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 R, as shown in Figure 3.

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 tumor ) to the whole tumor volume (V all ); (2) the thermal damage to vessels, D t in Equation (20).
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 (h) [23][24][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 h, 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 h ¼ 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 h 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

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 2 adj ) was used to evaluate the accuracy of model, and each variable in the model were considered as statistical significance at p < 0.05.

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 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 ex 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. 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 h), 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).

RSM results
where Y represents the responses (i.e., R a or D t ), t n , x n , z n , and h n are rescaled variables of t, x, z, and h 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   Figure 4 in [37]. Table 4. Results of ablation and thermal damage for each designed simulation. No.

Variables
Responses significances of z and h 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 h 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 h are only presented in the response of D t with p < 0.001, 0.039, and < 0.001, respectively.

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 h) 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 h ¼ 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.

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 h) 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 (h) 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.

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.