Automatic control of finite element models for temperature-controlled radiofrequency ablation

Background The finite element method (FEM) has been used to simulate cardiac and hepatic radiofrequency (RF) ablation. The FEM allows modeling of complex geometries that cannot be solved by analytical methods or finite difference models. In both hepatic and cardiac RF ablation a common control mode is temperature-controlled mode. Commercial FEM packages don't support automating temperature control. Most researchers manually control the applied power by trial and error to keep the tip temperature of the electrodes constant. Methods We implemented a PI controller in a control program written in C++. The program checks the tip temperature after each step and controls the applied voltage to keep temperature constant. We created a closed loop system consisting of a FEM model and the software controlling the applied voltage. The control parameters for the controller were optimized using a closed loop system simulation. Results We present results of a temperature controlled 3-D FEM model of a RITA model 30 electrode. The control software effectively controlled applied voltage in the FEM model to obtain, and keep electrodes at target temperature of 100°C. The closed loop system simulation output closely correlated with the FEM model, and allowed us to optimize control parameters. Discussion The closed loop control of the FEM model allowed us to implement temperature controlled RF ablation with minimal user input.


Background
Radio-frequency (RF) ablation has become of considerable interest as a minimally invasive treatment for primary and metastatic liver tumors. Hepatocellular carcinoma is one of the most common malignancies, worldwide with an estimated annual number of 500,000 deaths [1]. Surgi-cal resection offers the best chance of long-term survival, but is rarely possible. In many patients with cirrhosis or with multiple tumors, hepatic reserve is inadequate to tolerate resection and alternative means of treatment are necessary [2]. In RF ablation, RF current of 450 to 500 kHz is delivered to the tissue via electrodes inserted percutaneously or during surgery. Different modes of controlling the electromagnetic power delivered to tissue can be utilized. Power-controlled mode (P = constant), temperature-controlled mode (T = constant) and impedancecontrolled mode (Z < constant) are commonly used. The electromagnetic energy is converted to heat by resistive heating. Tissue damage can occur at temperatures above 43°C with heating times of several hours; at 50°C cell necrosis occurs after ~3 min [12]. A commonly used mode is temperature-controlled ablation, where the tip temperature of the electrodes is kept at a predetermined value, usually around 100°C. The hepatic electrodes (Fig. 1) used for temperature-controlled ablation have temperature sensing elements (thermistors or thermocouples) embedded in the prong tips. The sensors report temperature back to the generator, which then applies an appropriate amount of power to the electrode to keep the temperature constant.
Researchers have been using the finite element method (FEM) to simulate both cardiac and hepatic RF ablation [3][4][5][6][7][8][9][10][11]. When using the FEM to model temperature-controlled ablation, the applied voltage has to be adjusted to keep the tip temperature constant. Previously, most researchers used manual adjustment of applied voltage and trial-and-error methods to perform temperature-controlled ablation [3,[5][6][7][8]. Implementation of temperaturecontrolled feedback in models is of importance to obtain results comparable to clinical devices that use this type of control. We modeled a clinically used hepatic ablation electrode (15 gauge, RITA medical systems model 30) as described previously [6]. We implemented a control algorithm for a PI controller -a commonly used controller type -in a C++ program to change the applied voltage between the time steps. We used a closed loop system simulation to optimize control parameters for the PI controller.
Geometry of fully deployed Rita model 30 umbrella probe used in FEM model

Finite element method
RF ablation destroys tissue by thermal energy, which is converted from electric energy. The current flows from the conductive electrode through the tissue to a surface dispersive electrode. Tissue in close vicinity of the electrode tip is heated by resistive heating.
The heating of tissue during RF ablation is governed by the bioheat equation: where ρ is the density (kg/m 3 ), c is the specific heat (J/ (kg·K)), and k is the thermal conductivity (W/(m·K)). J is the current density (A/m 2 ) and E is the electric field intensity (V/m). T bl is the temperature of blood, ρ bl is the blood density (kg/m 3 ), c bl is the specific heat of the blood (J/(kg·K)), and w bl is the blood perfusion (1/s). h bl is the convective heat transfer coefficient accounting for the blood perfusion. Q m (W/m 3 ) is the energy generated by metabolic processes and was neglected since it is small compared to the other terms [13]. Equation (1) defines the solution in the spatial domain encompassing electrodes and tissue.
Tissue properties were assumed temperature independent and are described in more detail in [6]. We used the commercial software ABAQUS/Standard 6.3 (Hibbitt, Karlsson & Sorensen, Inc., Pawtucket, RI) for solving the coupled thermo-electrical analysis. All analysis was performed on a SUN Blade-1000 workstation equipped with 2.5 GB of RAM and 80 GB of hard disk space. Fig. 1 shows the geometry of the RITA model 30, 4-prong electrode we used in our models. The electrode was placed within a cylinder (80 mm diameter, 50 mm length) of liver tissue. The outer surfaces of the cylinder were set to 37°C (thermal boundary condition), and 0 V (electrical boundary condition). We performed quasi-static analysis. Due to the symmetry of the arrangement, we could reduce computing time by only modeling a quarter of the cylinder. We used the same model as in a previous study [3]; the model consisted of ~35,000 tetrahedral elements and 7,000 nodes. The node spacing was small next to the electrode (0.2 mm) and larger at the model boundary (2 Step response of the original dynamic system (FEM model) and of the approximation Figure 2 Step response of the original dynamic system (FEM model) and of the approximation.
. Perfusion was included in the model according to the Pennes model [14]. The blood perfusion wbl used in this model was 6.4·10 -3 1/s [15].
An input file was submitted to ABAQUS, in which geometry, material properties, step time and boundary conditions were specified. The applied voltage was one of the boundary conditions and remained constant during each step. ABAQUS created a results file in which temperatures of all nodes of the FEM model at the end of the step were written. We created a C++ program that reads the temperature at the electrode tip from the results file and creates a new input file where the applied voltage is set according to a control algorithm; the C++ program then calls ABAQUS and waits until the next step is completed. We implemented a PI controller in our control algorithm.

Closed loop system simulation
Since the FEM model takes several hours to complete, using the FEM model to optimize control parameters was not feasible. Therefore we used a closed loop simulation to adjust parameters of the controller, and then used these parameters in the FEM model. Before implementing the controller, we analyzed the dynamic system (i.e. the FEM model) to be controlled. This system consisted of the ablation electrode, tissue, and dispersive electrode. The input variable of the system was the voltage applied to the electrode. The output variable was the temperature measured at the tip of the electrode. Initially, we determined the step response by applying a constant voltage of 25 V for 180 s. We then approximated the transfer function of this system by a time-discrete transfer function of the following form: We used the control system simulation software ANA 2.52 (Freeware, Dept. of Control Engineering, Tech. Univ. Vienna/Austria) to analyze the control system. This software allowed us to approximate the system from its step response by a recursive least square algorithm and gave the parameters a 0 , a 1 , a 2 , b 1 and b 2 of (2). Fig. 2 shows the step responses of the original dynamic system (i.e. the FEM model) and of the approximation according to (2). The parameters used for the approximation in (2) were: a 0 = 1.0, a 1 = -0.959, a 2 = 0.127, b 1 = 1.734, b 2 = -1.150. The sampling time of the approximated system was 10 s, since that was also the sampling time used later in the digital PI controller. This ensured a more accurate simulation model of the closed-loop control system.
Once we had identified an approximation of the dynamic system (FEM model), we designed a feedback control system. There are different ways of controlling the tip temperature such as PID control, adaptive control, neural network prediction control and fuzzy logic control. We chose the relatively simple PI controller for our control system. Fig. 3 shows the complete closed-loop control system. T s is the desired set tip temperature and T t is the current tip temperature. The input of the PI controller e = T s -T t . The output of the PI controller u (corresponds to the applied voltage) is fed into the dynamic system.
The transfer function of the PI controller is: Closed loop system incorporating a PI controller and the dynamic system (FEM model) Figure 3 Closed loop system incorporating a PI controller and the dynamic system (FEM model). To implement this controller in software, we have to use the discrete time-domain version of (4), The second term of (5) represents the approximation of the integral term in (4) by trapezoidal numerical integration. The behavior of the PI controller is determined by the two parameters K p and K i .
We simulated the behavior of the closed-loop system with the software ANA. From in-vivo experiments in pigs with the RITA 500 generator and the Model 30 electrode we determined that it takes between 1 and 2 min for the tip temperature to reach a target temperature of 100°C, which is the temperature used in clinical cases at our insti-tution. The RITA 500 generator does not allow recording of the tip temperature signal, so we don't have exact data of the tip temperature over time. We empirically chose parameters for the PI controller to minimize overshoot and obtain similar temporal behavior as in the in-vivo experiments.
The parameters we used for the PI controller were K p = 0.02, K i = 0.0064. Fig. 4 shows the results of the closed-loop control system simulation. The target tip temperature of 100°C was reached 100 s after start of the ablation. The maximum overshoot was 11%, which was reached after 148 s. The maximum voltage of 24.5 V was applied after 100 s.

Controlled FEM model
We implemented the PI controller with the parameters resulting from the control system simulation in a control program written in C++. The software control program determined the applied voltage and the step time. An input file for ABAQUS was created and ABAQUS solved Temporal behavior of tip temperature and applied voltage of the closed loop system The tip temperature was read from the result file by the control program. The program then determined the applied voltage and step size for the next step, and ABAQUS was restarted with the modified input file. With decreasing step size (i.e. time that ABAQUS simulated the model with constant voltage) total computation time gets longer. As a compromise, we chose 10 s as the initial step size. Note that the FEM solver divided each step into smaller increments, starting at 0.05s. Also, the FEM solver performed convergence tests to ensure that the increment size was sufficiently small.
This scheme was repeated until the ablation had been simulated for the desired time. Fig. 6 is a flow chart of the algorithm implemented in the control program. The initial step size was 10 s. Subsequently, the step size increased once the tip temperature change between the steps decreased below a certain value.
We simulated ablation for 12 min using the FEM model and the control program. The simulation took 250 min to run and was divided into 35 steps. The individual steps took between 6 and 11 min to run in ABAQUS. Without using the control program, manual interaction would be necessary after each step to manually change the input file for ABAQUS and apply a new input voltage and set the step size. In this case the operator would have to check the results file after each step (i.e. in 6 to 11 min intervals), modify the input file and restart ABAQUS.  Fig. 4 also shows the results of the control system simulation. There is good correlation between the control simulation and the closed loop system consisting of the FEM model and the control software. We explain deviations between the two by inaccuracies of the approximation (see (2)) and differences in Temporal behavior of tip temperature and applied voltage of the closed loop system step sizes. In the control system simulation a constant sampling time (i.e. step time) of 10 s was used. However, the algorithm implemented in the control program did not use constant step times.

Results and Discussion
Step time was increased if the change in tip temperature between the steps was below a certain value. Note that the parameters of the controller must be modified when boundary conditions (e.g. perfusion), electrode geometry etc. are changed. Otherwise there will be changes in the dynamic behavior of the closed-loop system. Fig. 5 shows the temperature and Flowchart of the control software implementing a PI control algorithm Figure 6 Flowchart of the control software implementing a PI control algorithm.

Initial conditions:
Step  Fig. 4). Both the overshoot and the settling time of tip temperature for the case with perfusion were higher. Also, the voltage required to keep the tip temperature at the set temperature value was higher because the perfusion carried heat away. To obtain the same performance during the initial period in both cases, different controllers had to be used, e.g. the parameters of the PI controller had to be modified. Our model only included a quarter of the actual electrode. In models where non-uniform heating due to vessels occurs, the electrodes will obtain different temperatures [6]. In this case, the temperature of the hottest electrode should be used for control to avoid tissue overheating.
In hepatic RF ablation, ablation times clinically used go up to 35 min. Since the heat-up period is comparably much smaller (1 to 2 min.), it is not of essential importance that the temporal behavior of the control algorithm reproduces the control algorithm used in clinical devices during the heat-up period. From our experience, the temperature distribution in the FEM model reaches close to steady state at the end of the simulation due to the long simulation times. As long as the tip temperature is kept within a small range around the target temperature after the initial heat-up period, the model results (i.e. final temperature distribution) should not differ significantly. However, with knowledge of the actual control parameters and algorithms of commercial devices (e.g. obtained from measurements), accurate simulations of these devices is possible using our methods.

Conclusion
We implemented a closed loop control system to a FEM model, to automate simulation of temperature-controlled RF ablation. We further used a closed loop control system simulation to optimize control parameters. Previously, researchers often applied constant power, or used timeconsuming trial-and-error methods to determine required voltage. Furthermore, if control parameters and algorithms of commercial devices are known or can be measured, an accurate simulation of commercial devices is possible.

Authors' contributions
DH carried out computer simulations. JG participated in design of the study. All authors read and approved the final manuscript.