A New Conformal Penetrating Heating Strategy for Atherosclerotic Plaque

(1) Background: A combination of radiofrequency (RF) volumetric heating and convection cooling has been proposed to realize plaque ablation while protecting the endothelial layer. However, the depth of the plaque and the thickness of the endothelial layer vary in different atherosclerotic lesions. Current techniques cannot be used to achieve penetrating heating for atherosclerosis with two targets (the specified protection depth and the ablation depth). (2) Methods: A tissue-mimicking phantom heating experiment simulating atherosclerotic plaque ablation was conducted to investigate the effects of the control parameters, the target temperature (Ttarget), the cooling water temperature (Tf), and the cooling water velocity (Vf). To further quantitatively analyze and evaluate the ablation depth and the protection depth of the control parameters, a three-dimensional model was established. In addition, a conformal penetrating heating strategy was proposed based on the numerical results. (3) Results: It was found that Ttarget and Tf were factors that regulated the ablation results, and the temperatures of the plaques varied linearly with Ttarget or Tf. The simulation results showed that the ablation depth increased with the Ttarget while the protection depth decreased correspondently. This relationship reversed with the Tf. When the two parameters Ttarget and Tf were controlled together, the ablation depth was 0.47 mm–1.43 mm and the protection depth was 0 mm–0.26 mm within 2 min of heating. (4) Conclusions: With the proposed control algorithm, the requirements of both the ablation depth and the endothelium protection depth can be met for most plaques through the simultaneous control of Ttarget and Tf.


Introduction
Restenosis remains a challenge in treating atherosclerosis with percutaneous transluminal angioplasty (PTA) [1][2][3][4]. The main reason for restenosis is the endothelial injuries caused by the mechanical force of balloon or stent expansion and the alteration of hemodynamics due to stent struts as well as the consequent proliferation of smooth muscle cells (SMCs) [1,5]. Thermal ablation can safely and efficiently remove unwanted tissue [6,7] and has also been used in the treatment of atherosclerosis [8][9][10]. Studies have demonstrated that radiofrequency (RF) balloon angioplasty can effectively ablate atherosclerotic plaque, open occluded vessels, and fuse dissected vascular media [11]. However, restenosis persists since current thermal techniques inevitably overheat the thin endothelium [10,12].
New thermal balloon angioplasty has been proposed to realize the selective heating of atherosclerotic plaque [13]. A radiofrequency microelectrode array is embedded on the surface of a balloon for penetrating volume heating of the plaque tissue, while the circulating cooling agent inside the balloon provides surface convection cooling to protect the endothelial layer [14]. The feasibility of the selective heating strategy was confirmed, but how to ablate plaques of different thicknesses remains unknown. To achieve a successful treatment for plaques of different thicknesses, it is important to precisely control the two modalities at the same time. When the RF volumetric heating power is excessive, overheating occurs and damages the endothelium and the surrounding nerves, which causes restenosis and other secondary injuries [2,15], while overprotection of the endothelium and surrounding tissues results in insufficient plaque ablation.
Successful heating strategies to achieve the desired ablation depth while protecting the thin endothelial layer are quite different from the conformal radiofrequency ablation of tumors. Tumor ablation needs to precisely control the far end of the ablation and sufficiently cover the distal margin of the tumor (several centimeters away from the probe) while sparing the surrounding tissues or important organs. This was achieved through the optimization of the heating power or electrode configuration design [16,17]. For the ablation of atherosclerotic plaques, two parameters, the penetrating ablation depth and the protection depth, need to be satisfied at the same time.
In addition, the thickness of the endothelial layer and the depths of the atherosclerotic plaques vary greatly in different atherosclerotic lesions due to different levels of stenosis and media hyperplasia [18][19][20]. As shown in Figure 1, a typical plaque includes a necrotic core, proliferative SMCs and an overlying intima [3,21]. The necrotic core and the proliferative SMCs are the unwanted tissue and need to be ablated, while the intima, containing a thin endothelial layer, needs to be protected. The intima-media thickness (IMT), which describes the plaque thickness, is approximately 0.55-0.95 mm in peripheral arteries [22] and 0.81-1.01 mm in the femoral artery [18], and the intimal thickness ranges from 0.05 mm to 0.12 mm [23]. The required ablation depth (A) and the protection depth (P), as defined in Figure 1, are different in every lesion. A conformal ablation strategy for atherosclerotic plaques with various thicknesses that protects the endothelium, which also has different dimensions, is necessary. The target temperature for proliferated SMCs is over 50 • C [24,25], while the temperature in the whole layer of the intima should be under 50 • C, which means that different protection depth requirements and different thermal gradients are needed in the intima. endothelial layer [14]. The feasibility of the selective heating strategy was confirmed, but how to ablate plaques of different thicknesses remains unknown. To achieve a successful treatment for plaques of different thicknesses, it is important to precisely control the two modalities at the same time. When the RF volumetric heating power is excessive, overheating occurs and damages the endothelium and the surrounding nerves, which causes restenosis and other secondary injuries [2,15], while overprotection of the endothelium and surrounding tissues results in insufficient plaque ablation. Successful heating strategies to achieve the desired ablation depth while protecting the thin endothelial layer are quite different from the conformal radiofrequency ablation of tumors. Tumor ablation needs to precisely control the far end of the ablation and sufficiently cover the distal margin of the tumor (several centimeters away from the probe) while sparing the surrounding tissues or important organs. This was achieved through the optimization of the heating power or electrode configuration design [16,17]. For the ablation of atherosclerotic plaques, two parameters, the penetrating ablation depth and the protection depth, need to be satisfied at the same time.
In addition, the thickness of the endothelial layer and the depths of the atherosclerotic plaques vary greatly in different atherosclerotic lesions due to different levels of stenosis and media hyperplasia [18][19][20]. As shown in Figure 1, a typical plaque includes a necrotic core, proliferative SMCs and an overlying intima [3,21]. The necrotic core and the proliferative SMCs are the unwanted tissue and need to be ablated, while the intima, containing a thin endothelial layer, needs to be protected. The intima-media thickness (IMT), which describes the plaque thickness, is approximately 0.55-0.95 mm in peripheral arteries [22] and 0.81-1.01 mm in the femoral artery [18], and the intimal thickness ranges from 0.05 mm to 0.12 mm [23]. The required ablation depth (A) and the protection depth (P), as defined in Figure 1, are different in every lesion. A conformal ablation strategy for atherosclerotic plaques with various thicknesses that protects the endothelium, which also has different dimensions, is necessary. The target temperature for proliferated SMCs is over 50 °C [24,25], while the temperature in the whole layer of the intima should be under 50 °C, which means that different protection depth requirements and different thermal gradients are needed in the intima. Therefore, in this study, the influence of the control parameters on heat generation and heat transfer in the newly proposed radiofrequency balloon angioplasty was investigated both experimentally and numerically. Based on the relationship between the ablation results and the control parameters, a control algorithm was proposed to determine the feasibility of bidirectional control for realizing the conformal heating of plaques. Therefore, in this study, the influence of the control parameters on heat generation and heat transfer in the newly proposed radiofrequency balloon angioplasty was investigated both experimentally and numerically. Based on the relationship between the ablation results and the control parameters, a control algorithm was proposed to determine the feasibility of bidirectional control for realizing the conformal heating of plaques.

Experimental Setup
A tissue-mimicking phantom, which consisted of 5.5% gelatin, 0.46% NaCl powder, and deionized water [26], was used to mimic arterial tissue. As shown in Figure 2, a polyimide catheter (inner diameter: 2.5 mm, outer diameter: 2.94 mm) simulating an angioplasty balloon was inserted into the tissue-mimicking phantom, which was placed in a 37 • C thermostat water bath (DK-8A, Shanghai Jinghong Experimental Equipment Co., Ltd., Shanghai, China). To achieve the bidirectional depth control of the penetrating thermal treatment, the bipolar model was chosen. Therefore, each time, only one pair of electrodes ( Figure 3a) was used to deliver the 460 kHz radiofrequency energy to the phantom, though there were 3 × 4 electrodes (width and spacing: 1 mm, length: 4 mm) attached to the outer surface of the catheter. The water inside the low-temperature thermostat bath was pumped into the catheter through a soft silicone tube by an adjustable-speed pump (DC40D-2480A, Zhongke Electromechanical Co., Ltd., Shenzhen, China). The host computer controlled the homemade radiofrequency source [27], which was calibrated with an oscilloscope (Agilent InfiniiVision, MSO-X 2022A, Keysight Technologies, Inc., Santa Rosa, CA, USA), and the adjustable-speed pump via a data acquisition device (DAQ, USB3102A, ART Technology, Beijing, China).

Experimental Setup
A tissue-mimicking phantom, which consisted of 5.5% gelatin, 0.46% NaCl powder, and deionized water [26], was used to mimic arterial tissue. As shown in Figure 2, a polyimide catheter (inner diameter: 2.5 mm, outer diameter: 2.94 mm) simulating an angioplasty balloon was inserted into the tissue-mimicking phantom, which was placed in a 37 °C thermostat water bath (DK-8A, Shanghai Jinghong Experimental Equipment Co., Ltd, Shanghai, China). To achieve the bidirectional depth control of the penetrating thermal treatment, the bipolar model was chosen. Therefore, each time, only one pair of electrodes ( Figure 3a) was used to deliver the 460 kHz radiofrequency energy to the phantom, though there were 3×4 electrodes (width and spacing: 1 mm, length: 4 mm) attached to the outer surface of the catheter. The water inside the low-temperature thermostat bath was pumped into the catheter through a soft silicone tube by an adjustable-speed pump (DC40D-2480A, Zhongke Electromechanical Co., Ltd, Shenzhen, China). The host computer controlled the homemade radiofrequency source [27], which was calibrated with an oscilloscope (Agilent InfiniiVision, MSO-X 2022A, Keysight Technologies, Inc., Santa Rosa, CA, USA), and the adjustable-speed pump via a data acquisition device (DAQ, USB3102A, ART Technology, Beijing, China).

Figure 2.
Illustration of the tissue-mimicking phantom heating experiment. The radiofrequency catheter system consisted of an RF source and a microelectrode array, The cooling unit (an adjustable pump and a low-temperature thermostat bath), the control unit (DAQ and host computer), and the tissue-mimicking phantom were immersed in a 37 °C thermostat bath. A temperature measurement device connected to the host computer was used to measure, display, and record the temperature.
Three calibrated temperature sensors (TC1, TC2, and TC3) (K-Type thermocouple, 44AWG, Sandvik, Sweden) were inserted into the tissue-mimicking phantom parallel to the catheter to monitor the phantom temperature in the depth direction (see Figure 3a). They were positioned vertically above the middle of the two working electrodes using the holder with holes. The phantom with thermocouples was photographed from the top and side views before and after 2 minutes of heating. The accurate positions of the Figure 2. Illustration of the tissue-mimicking phantom heating experiment. The radiofrequency catheter system consisted of an RF source and a microelectrode array, The cooling unit (an adjustable pump and a low-temperature thermostat bath), the control unit (DAQ and host computer), and the tissue-mimicking phantom were immersed in a 37 • C thermostat bath. A temperature measurement device connected to the host computer was used to measure, display, and record the temperature.
Three calibrated temperature sensors (TC1, TC2, and TC3) (K-Type thermocouple, 44AWG, Sandvik, Sweden) were inserted into the tissue-mimicking phantom parallel to the catheter to monitor the phantom temperature in the depth direction (see Figure 3a). They were positioned vertically above the middle of the two working electrodes using the holder with holes. The phantom with thermocouples was photographed from the top and side views before and after 2 min of heating. The accurate positions of the thermocouples were determined using the image processing software Image J. The distances between TC1, TC2, and TC3 and the outer surface of the catheter were 0.85 ± 0.09 mm, 2.46 ± 0.17 mm, and 4.08 ± 0.13 mm. Therefore, TC1, TC2, and TC3 could represent the locations of plaque, adventitia, and surrounding tissue, respectively [28]. Another thermocouple, TC0 (K-Type thermocouple, 44AWG, Sandvik, Sweden), was located on the inner surface of the tissue phantom and represented the inner surface of the blood vessel. TC0 was positioned in the middle of the two working electrodes to measure the temperature at the control point T target (Figure 3b), and TC0 was isolated with polyimide film (thickness: 0.05 mm) to avoid electromagnetic interference. Similarly, to reduce electromagnetic interference and minimize interference in temperature measurement, TC1, TC2, and TC3 were insulated with thin polyimide tubes (inner diameter: 0.2 mm, thickness: 0.02 mm), and their hot junction tips were covered with thermally conductive silicone grease (STARS-922, Balance Stars, China), which has good electrical insulation and thermal conductivity (>0.671 W/(m·K)). mm, 2.46 ± 0.17 mm, and 4.08 ± 0.13 mm. Therefore, TC1, TC2, and TC3 could represent the locations of plaque, adventitia, and surrounding tissue, respectively [28]. Another thermocouple, TC0 (K-Type thermocouple, 44AWG, Sandvik, Sweden), was located on the inner surface of the tissue phantom and represented the inner surface of the blood vessel. TC0 was positioned in the middle of the two working electrodes to measure the temperature at the control point Ttarget (Figure 3b), and TC0 was isolated with polyimide film (thickness: 0.05 mm) to avoid electromagnetic interference. Similarly, to reduce electromagnetic interference and minimize interference in temperature measurement, TC1, TC2, and TC3 were insulated with thin polyimide tubes (inner diameter: 0.2 mm, thickness: 0.02 mm), and their hot junction tips were covered with thermally conductive silicone grease (STARS-922, Balance Stars, China), which has good electrical insulation and thermal conductivity (>0.671 W/(m•K)).

Heating Conditions
To protect the endothelial layer from overheating, the temperature-controlled mode is chosen, and the adjustable parameters included the target temperature (Ttarget), the cooling water temperature (Tf), and the cooling water velocity (Vf). Specifically, Ttarget is a parameter that adjusts the RF heating power, while Tf and Vf are parameters that adjust the convective cooling power. To determine the key parameters for ablation control, the influence of these parameters was investigated. Three groups of experiments were designed to study the influence of the parameters on the temperature profile inside the phantom:  [29] (Appendix A, Equations (A1)-(A5)), which was also used in the numerical simulation.
Considering the short duration of several minutes in the real procedure of conventional balloon angioplasty [25,30], the whole experiment included two stages: a 45-second expansion of the balloon using cooling water and 2 minutes of convection cooling combined with RF heating. Temperatures (T1, T2, and T3) were measured every 500 ms using

Heating Conditions
To protect the endothelial layer from overheating, the temperature-controlled mode is chosen, and the adjustable parameters included the target temperature (T target ), the cooling water temperature (T f ), and the cooling water velocity (V f ). Specifically, T target is a parameter that adjusts the RF heating power, while T f and V f are parameters that adjust the convective cooling power. To determine the key parameters for ablation control, the influence of these parameters was investigated. Three groups of experiments were designed to study the influence of the parameters on the temperature profile inside the phantom: Group one: T target was set at 35 • C, 38 • C, 40 • C, 43 • C, 45 • C, and 48 • C. T f was 15 • C, and V f was 2.85 m/s.
Group two: T f was adjusted from 10 • C to 30 • C in 5 • C intervals with T target set at 38 • C and V f set at 2.85 m/s.
Group three: V f ranged from 1.51 m/s to 4.57 m/s in steps of 0.76 m/s, while T target was 38 • C and T f was 20 • C. The corresponding volume flow (Q) ranged from 7.4 ml/s to 22.43 ml/s, and the convective heat transfer coefficient ranged from 7221 W/(m 2 ·K) to 21,745 W/(m 2 ·K), with an average interval of 3631 W/(m 2 ·K). The convective heat transfer coefficient (h) was calculated according to the Gnielinski formula [29] (Appendix A, Equations (A1)-(A5)), which was also used in the numerical simulation.
Considering the short duration of several minutes in the real procedure of conventional balloon angioplasty [25,30], the whole experiment included two stages: a 45-s expansion of the balloon using cooling water and 2 min of convection cooling combined with RF heating. Temperatures (T1, T2, and T3) were measured every 500 ms using thermocouples (TC1, TC2, and TC3) with Keysight DAQ970A. Temperature T0 served as the control point. Each experiment was repeated three times.

Model Geometry
To further analyze the relationship between the ablation dimensions and the treatment parameters, an applicable three-dimensional model was established to obtain the temperature and SAR (specific absorption rate) distribution. Different from the radiofrequency catheter ablation in cardiac arrhythmias, such as in [31], the bioheat trans fer model coupled with an RF electric field and cooling water circulating inside a catheter was studied. At the same time, the temperature dependences of the tissue's electrical conductivity and thermal conductivity were considered in our model, which was different from previous studies [32,33].
The three-dimensional (3D) geometry was established in COMSOL Multiphysics (COMSOL, Inc., Burlington, MA, USA). A 50 × 50 × 24 mm 3 block and a hollow cylinder with an inner diameter of 2.5 mm and an outer diameter of 2.94 mm in the middle were used to simulate the tissue and the catheter ( Figure 4a). As the fibrotic plaque was mainly composed of proliferating SMCs, which have similar properties to the arterial wall, the fibrotic plaque and the arterial wall were regarded as homogenous [34]. The paired copper electrodes (thickness: 25 µm, width: 1 mm, length: 4 mm) with a spacing of 1 mm were distributed circumferentially on the outer layer of the polyimide catheter. To ensure the precision of the calculation, a finer mesh (minimum size: 0.075 mm, maximum size: 1 mm) was used in the fan-shaped zone close to the working electrodes and the tissue (Figure 4b).
In total, there were 380,825 tetrahedral mesh elements, 50,862 triangular mesh elements, and 1650 edge mesh elements in our model. More detailed mesh information is provided in Appendix A (Table A1 and Figure A4).

Model Geometry
To further analyze the relationship between the ablation dimensions and the treatment parameters, an applicable three-dimensional model was established to obtain the temperature and SAR (specific absorption rate) distribution. Different from the radiofrequency catheter ablation in cardiac arrhythmias, such as in [31], the bioheat transfer model coupled with an RF electric field and cooling water circulating inside a catheter was studied. At the same time, the temperature dependences of the tissue's electrical conductivity and thermal conductivity were considered in our model, which was different from previous studies [32,33].
The three-dimensional (3D) geometry was established in COMSOL Multiphysics (COMSOL, Inc., Burlington, MA, USA). A 50 × 50 × 24 mm 3 block and a hollow cylinder with an inner diameter of 2.5 mm and an outer diameter of 2.94 mm in the middle were used to simulate the tissue and the catheter (Figure 4a). As the fibrotic plaque was mainly composed of proliferating SMCs, which have similar properties to the arterial wall, the fibrotic plaque and the arterial wall were regarded as homogenous [34]. The paired copper electrodes (thickness: 25 μm, width: 1 mm, length: 4 mm) with a spacing of 1 mm were distributed circumferentially on the outer layer of the polyimide catheter. To ensure the precision of the calculation, a finer mesh (minimum size: 0.075 mm, maximum size: 1 mm) was used in the fan-shaped zone close to the working electrodes and the tissue ( Figure  4b). In total, there were 380825 tetrahedral mesh elements, 50862 triangular mesh elements, and 1650 edge mesh elements in our model. More detailed mesh information is provided in Appendix A (Table A1 and Figure A4).

Mathematical Modeling
For low frequencies, as used in RFA (460 kHz), Maxwell's equation can be simplified as the quasi-static approach, Equation (1): where is the electrical conductivity (S/m) and V is the electrical potential (V) in the tissue. For the two working electrodes, one was applied with voltage (V=V(k)), and the other one was set to 0. Electrical insulation conditions were applied to the surface in contact

Mathematical Modeling
For low frequencies, as used in RFA (460 kHz), Maxwell's equation can be simplified as the quasi-static approach, Equation (1): where σ is the electrical conductivity (S/m) and V is the electrical potential (V) in the tissue. For the two working electrodes, one was applied with voltage (V = V(k)), and the other one was set to 0. Electrical insulation conditions were applied to the surface in contact with the catheter and the outer surface of the phantom, as it was exposed to pure water with low electrical conductivity (2 × 10 −4 S/m) [35].
To quickly reach T target , avoid overshooting, and achieve a small steady-state error, the maximum voltage (30.9 (V)) was used when the temperature difference (e(k)) between the target temperature and the measured temperature was larger than 4 • C. When the difference is less than 4 • C, the voltage V(k) was calculated with the PID algorithm, which ranged from 13.5 (V) to 30.9 (V). The lower and upper limits of the voltage were manually adjusted. The PID-controlled voltage (V(k)) was calculated using the following formulas, Equations (2)- (4): where T target and T 0 (k) are the target temperature and the measured temperature in each sample, respectively. u(k) (V) is the control function in each sample. The proportional gain (K p ), integral gain (K i ), and derivative gain (K d ) were equal to 1.30 (V/K), 0.19 (V/K), and 0 (V/K), respectively. K p , K i , K d and the range of u(k) in Equation (4) were set to be the same as in the experiment, which was manually adjusted to ensure reaching the target temperature quickly (<10 s) and to keep steady with a control precision of ±0.5 • C. In the experiment (Section 2.1), K d was set to 0.56 (V/K) to avoid overshooting, given the response times of the thermocouples. The volumetric heat generation rate (Q ec ) (W/m 3 ) was calculated using Equation (5): where E is the root-mean-square value of the electric field intensity (V/m) and σ is the electrical conductivity (S/m). The heat transfer governing equation is given as follows, Equation (6): where ρ, c, and k are the density (kg/m 3 ), the specific heat capacity (J/(kg·K)), and the thermal conductivity (W/(m·K)), respectively. The electrical and thermal properties of the phantom, the polyimide catheter, and the copper electrodes are listed in Table 1.
The convective thermal boundary was applied to the inner surface of the catheter to simulate the flow of cooling water in the catheter, Equation (7): where h (W/(m 2 ·K)) is the convective heat transfer coefficient and T w1 and T f are the temperature of the inner wall of the balloon catheter and the cooling water temperature.
To simulate the in vivo thermal treatment, the initial temperature was set to 37 • C and the outer surface of the phantom was set as thermal insulation.

Optimal Heating Strategy
To further determine the heating strategy for the conformal treatment of the fibrotic plaques, an optimization algorithm was proposed based on the quantitative relationships between the ablation dimensions and different combinations of the parameters, which were obtained from the numerical model. Figure 5a illustrates the temperature at three monitoring points with the target temperature (T target ), the cooling water temperature (T f ), and the flow velocity (V f ) set to be 38 • C, 20 • C, and 2.85 m/s, respectively. It can be seen that the temperature (T0) at the control point reached the setpoint (38 • C) in less than 10 s. As seen in Figure 5a, T1, which represented a point in the plaque, had the highest temperature, while T0, which was located on the inner surface of the phantom, was the lowest. The temperatures of T2 (46.11 ± 0.01 • C) and T3 (40.89 ± 0.02 • C) were much lower than that in the plaque (T1: 62.12 ± 0.05 • C), over 16-22 • C. As temperatures T1-T3 reached equilibrium in less than 45 s, the final equilibrium temperatures under different conditions were used and are shown in Figure 5b-d.

Numerical Modeling
The numerical model was validated using the tissue-mimicking phantom heating experiments. The comparison of the numerical calculation results with the experimental results is presented in Figures A1-A3 (Ttarget = 40 °C, Tf = 15 °C, and Vf = 2 m/s). A one-way sensitivity analysis was performed to verify the robustness of the model [39]. It was found that when the electrical conductivity increased by 20% from 0.30 S/m to 0.36 S/m, the maximum value of temperature T1 (the point located in the plaque) changed from 59.92 °C to 59.93 °C, while when the thermal conductivity increased by 20% from 0.51 W/(m•K) to 0.63 W/(m•K), the maximum T1 changed from 58.96 °C to 61.13 °C, which was less than 3%. The convergence of the model was verified, as shown in Figure A5.
The relationships between the lesion dimensions and the two main control parameters were investigated. Moreover, to evaluate the evolution of absorbed RF energy during the two-minute heating process, the SAR (specific absorption rate) distribution was also calculated.  Figure 5b illustrates the steady-state temperatures at the measured points (T1-T3) with different target temperatures (T target ). It was found that, for all three points, the final temperature increased linearly with T target (R 2 > 99.5%). The steady-state temperatures at these points were found to decrease linearly with the cooling water temperature (T f ) ( Figure 5c) and increase slightly with the flow rate (V f ) (Figure 5d).
The temperature of TC1, which represented the plaque area, changed with T target , T f , and V f . In particular, the range of the T1 temperature could be adjusted by more than 25 • C by adjusting T target or T f (Figure 5b,c), while it could be adjusted by less than 5 • C by adjusting V f (Figure 5d). T target and T f were the two major factors affecting the steady temperature in the phantom (T1-T3, especially T1).
However, only the temperatures at certain points were obtained through the experimental measurement. To determine the temperature distribution, especially to define the ablation range, numerical results are obtained.

Numerical Modeling
The numerical model was validated using the tissue-mimicking phantom heating experiments. The comparison of the numerical calculation results with the experimental results is presented in Figures A1-A3 (T target = 40 • C, T f = 15 • C, and V f = 2 m/s). A one-way sensitivity analysis was performed to verify the robustness of the model [39]. It was found that when the electrical conductivity increased by 20% from 0.30 S/m to 0.36 S/m, the maximum value of temperature T1 (the point located in the plaque) changed from 59.92 • C to 59.93 • C, while when the thermal conductivity increased by 20% from 0.51 W/(m·K) to 0.63 W/(m·K), the maximum T1 changed from 58.96 • C to 61.13 • C, which was less than 3%. The convergence of the model was verified, as shown in Figure A5.
The relationships between the lesion dimensions and the two main control parameters were investigated. Moreover, to evaluate the evolution of absorbed RF energy during the twominute heating process, the SAR (specific absorption rate) distribution was also calculated.

Numerical Results
A 50 • C isothermal curve was used to determine the ablation region, as 50 • C has been found to be the lethal temperature for arterial-related cells [24]. The ablation depth and protection depth were defined as the maximum and minimum distances of a 50 • C isothermal away from the outer surface of the catheter.
(a) The evolution of temperature/SAR distribution The evolution of the ablation zone during a two-minute heating process is shown in Figure 6a-d, with T target = 35 • C, T f = 18.5 • C, and V f = 2 m/s (h = 9303 W/(m 2 *K)). The ablation zone (marked with a green line) appeared approximately 0.25 mm away from the inner surface of the phantom after 12 seconds of heating, and it was initiated from the two corners between the paired electrodes with the shapes of raindrops (Figure 6a). Then, the two raindrops fused quickly and extended further (Figure 6b,c). The ablation depth increased with the heating time, while the protection depth decreased with time. After 30 s, the 50 • C isothermal regions, including the ablation depth and the protection depth, reached a steady state (Figure 6d).
To analyze the deposition of RF energy in tissue during the heating process, the SAR distribution is shown in Figure 6e-h. The distribution of SAR was found to mainly concentrated in the region close to the working electrodes, and it hardly changed with time though the electrical conductivity changes with the temperature. A comparison of Figure 6a-d with Figure 6e-f shows the influence of the surface convection effect. Volumetric heating (SAR) was established quickly in the tissue. It was the propagation of surface cooling that changed the temperature profile in the tissue.
(b) The effects of parameters on the dimensions of the lesions The target temperature (T target ) and the cooling agent temperature (T f ) were the two critical parameters regulating the temperature profile based on the previous experiment ( Figure 5). The local temperature in the tissue increased with T target and decreased with the flow temperature (T f ). To determine how the ablation zone varied with these two control parameters, different values of T target and T f were simulated, and the results are shown in Figure 7. With an increase in T target , the ablation zone expanded from the intermediate region, while with an increase in T f the ablation zone shrank. The two parameters had contrary effects on the ablation depth (A) and the protection depth (P). With a single control parameter, the ablation depth (A) and protection depth (P) changed synchronously, and it was impossible to realize the conformal control target for the diseased blood vessel. The two control parameters may need to be used simultaneously. distribution is shown in Figure 6e-h. The distribution of SAR was found to mainly concentrated in the region close to the working electrodes, and it hardly changed with time though the electrical conductivity changes with the temperature. A comparison of Figure  6a-d with Figure 6e-f shows the influence of the surface convection effect. Volumetric heating (SAR) was established quickly in the tissue. It was the propagation of surface cooling that changed the temperature profile in the tissue.

(b) The effects of parameters on the dimensions of the lesions
The target temperature (Ttarget) and the cooling agent temperature (Tf) were the two critical parameters regulating the temperature profile based on the previous experiment ( Figure 5). The local temperature in the tissue increased with Ttarget and decreased with the flow temperature (Tf). To determine how the ablation zone varied with these two control parameters, different values of Ttarget and Tf were simulated, and the results are shown in Figure 7. With an increase in Ttarget, the ablation zone expanded from the intermediate region, while with an increase in Tf the ablation zone shrank. The two parameters had contrary effects on the ablation depth (A) and the protection depth (P). With a single control parameter, the ablation depth (A) and protection depth (P) changed synchronously, and it was impossible to realize the conformal control target for the diseased blood vessel. The two control parameters may need to be used simultaneously.  Therefore, to further quantitatively establish the relationship between the ablation zone geometry and the two parameters, a parameter sweep analysis was performed. The target temperature (Ttarget) ranged from 29 °C to 43 °C, and the coolant temperature circulating inside the RF balloon catheter (Tf) ranged from 0.5 °C to 36.5 °C at an interval of 1 °C. Vf was kept at 2 m/s in all conditions. The ranges of the parameters were chosen to protect the endothelium from excessive heating.
Among the sweeping cases, if the temperature of the endothelial layer was below 50 °C and there were regions above 50 °C, it was defined as a feasible case, and the combined condition was labeled in green, as shown in Figure 8. The cases with excessive heating, where temperatures close to the endothelial layer were all above 50 °C, were labeled in red. The cases with no temperature higher than 50 °C were labeled in blue. The green diagram then clearly labeled all possible combined conditions that could be used to realize the control of both the protection depth and the ablation depth. It could also be seen that the green diagonal area with feasible cases, which accounted for 16.2% (90/555) of the total Therefore, to further quantitatively establish the relationship between the ablation zone geometry and the two parameters, a parameter sweep analysis was performed. The target temperature (T target ) ranged from 29 • C to 43 • C, and the coolant temperature circulating inside the RF balloon catheter (T f ) ranged from 0.5 • C to 36.5 • C at an interval of 1 • C. V f was kept at 2 m/s in all conditions. The ranges of the parameters were chosen to protect the endothelium from excessive heating.
Among the sweeping cases, if the temperature of the endothelial layer was below 50 • C and there were regions above 50 • C, it was defined as a feasible case, and the combined condition was labeled in green, as shown in Figure 8. The cases with excessive heating, where temperatures close to the endothelial layer were all above 50 • C, were labeled in red. The cases with no temperature higher than 50 • C were labeled in blue. The green diagram then clearly labeled all possible combined conditions that could be used to realize the control of both the protection depth and the ablation depth. It could also be seen that the green diagonal area with feasible cases, which accounted for 16.2% (90/555) of the total cases, was relatively narrow. zone geometry and the two parameters, a parameter sweep analysis was performed. The target temperature (Ttarget) ranged from 29 °C to 43 °C, and the coolant temperature circulating inside the RF balloon catheter (Tf) ranged from 0.5 °C to 36.5 °C at an interval of 1 °C. Vf was kept at 2 m/s in all conditions. The ranges of the parameters were chosen to protect the endothelium from excessive heating.
Among the sweeping cases, if the temperature of the endothelial layer was below 50 °C and there were regions above 50 °C, it was defined as a feasible case, and the combined condition was labeled in green, as shown in Figure 8. The cases with excessive heating, where temperatures close to the endothelial layer were all above 50 °C, were labeled in red. The cases with no temperature higher than 50 °C were labeled in blue. The green diagram then clearly labeled all possible combined conditions that could be used to realize the control of both the protection depth and the ablation depth. It could also be seen that the green diagonal area with feasible cases, which accounted for 16.2% (90/555) of the total cases, was relatively narrow.  To determine the control ranges of these two parameters, the relationships between the two parameters and the ablation depth (A) as well as the protection depth (P) in all feasible cases were analyzed, as shown in Figure 9. It can be seen that the ablation depth (A) ranged from 0.47 mm to 1.43 mm, while the protection depth (P) ranged from 0 mm to 0.26 mm, which covered the range for most arterial plaques [40,41], such as in peripheral [22], brachial, and radial arteries [23]. (Ttarget) ranged from 29 °C to 43 °C, and the cooling water temperature (Tf) ranged from 0.5 °C to 36.5 °C.
To determine the control ranges of these two parameters, the relationships between the two parameters and the ablation depth (A) as well as the protection depth (P) in all feasible cases were analyzed, as shown in Figure 9. It can be seen that the ablation depth (A) ranged from 0.47 mm to 1.43 mm, while the protection depth (P) ranged from 0 mm to 0.26 mm, which covered the range for most arterial plaques [40,41], such as in peripheral [22], brachial, and radial arteries [23].

Algorithm of Optimal Heating Strategy
With the obtained quantitative relationships between the ablation depth (A), the protection depth (P), and the two parameters (Ttarget and Tf), an optimal heating strategy was then established to investigate the possibility of the bidirectional control of ablation. Two implicit discrete functions describing the relationship, A = M (Tf, Ttarget) and P = N (Tf, Ttarget), were obtained from the above parameter sweep analysis. For the excessive ablation cases and insufficient ablation cases, the values of the protection depth or ablation depth were set to zero in M (Tf, Ttarget) and N (Tf, Ttarget), respectively. Then, the control problem was

Algorithm of Optimal Heating Strategy
With the obtained quantitative relationships between the ablation depth (A), the protection depth (P), and the two parameters (T target and T f ), an optimal heating strategy was then established to investigate the possibility of the bidirectional control of ablation. Two implicit discrete functions describing the relationship, A = M (T f , T target ) and P = N (T f , T target ), were obtained from the above parameter sweep analysis. For the excessive ablation cases and insufficient ablation cases, the values of the protection depth or ablation depth were set to zero in M (T f , T target ) and N (T f , T target ), respectively. Then, the control problem was simplified to an optimization problem with boundaries. The bilinear interpolation and iteration from these two functions were proposed based on the linear relationship between ablation depth/protection depth and T target /T f in the region. The targeted ablation depth (A t ) and the protection thickness (P t ) were used as input variables. The two existing discrete functions (M and N) were used to calculate the objective function { A t − M + (A t /P t ) P t − N }. Before starting iterative interpolation, the initial values T f,0 and T target,0 were obtained by calculating the function min{ A t − M + (A t /P t ) P t − N }. Then, A 0 and P 0 were obtained from the functions M (T f,0 , T target,0 ) and N (T f,0 , T target,0 ). A 0 and P 0 were compared with the targeted values of A t and P t . If the relative errors were less than 5% for both the protection depth and the ablation depth, the iteration continued. Two terminal conditions were used to stop the iteration: (1) if the objective function { A t − A k + (A t /P t ) P t − P k } was less than 0.01 mm and (2) if the iteration number (k) reached 10. During the iteration, T target,k and (or) T f,k were adjusted according to the distances between the current values (the ablation depth (A k ) and the protection depth (P k )) and the targeted values (A t and P t ) to obtain new T target, k+1 and T f, k+1 values. M and N were linearly interpolated bidirectionally from M (T target,k, T f,k ) and N (T target,k, T f,k ) to calculate the new corresponding ablation depth (A k+1 ) and protection depth (P k+1 ). Then, k, A k , and P k , as well as M and N, were updated, and the above procedures were repeated. A detailed algorithm flowchart is shown in Figure 10. well as M and N, were updated, and the above procedures were repeated. A detailed algorithm flowchart is shown in Figure 10. Two cases with different treatment requirements were used to illustrate the above algorithm. For case 1, the objective ablation depth and protection depth were set to be 0.78 mm and 0.06 mm (Figure 11a). According to the algorithm, the control parameters were then determined to be Ttarget = 41.9 °C and Tf = 33.6 °C. With these two parameters, the ablation geometry was calculated using the model described above. The ablation depth was found to be 0.78 mm, and the protection was found to be 0.06 mm. As the intimal thickness increased (1.02 mm in ablation depth and 0.1 mm in protection depth), as shown in case 2 (Figure 11b), Ttarget = 29.9 °C and Tf = 8.6 °C were obtained with the algorithm, and there were errors of 0.01 mm in both ablation depth and protection depth. Two cases with different treatment requirements were used to illustrate the above algorithm. For case 1, the objective ablation depth and protection depth were set to be 0.78 mm and 0.06 mm (Figure 11a). According to the algorithm, the control parameters were then determined to be T target = 41.9 • C and T f = 33.6 • C. With these two parameters, the ablation geometry was calculated using the model described above. The ablation depth was found to be 0.78 mm, and the protection was found to be 0.06 mm. As the intimal thickness increased (1.02 mm in ablation depth and 0.1 mm in protection depth), as shown in case 2 (Figure 11b), T target = 29.9 • C and T f = 8.6 • C were obtained with the algorithm, and there were errors of 0.01 mm in both ablation depth and protection depth.
Two cases with different treatment requirements were used to illustrate the above algorithm. For case 1, the objective ablation depth and protection depth were set to be 0.78 mm and 0.06 mm (Figure 11a). According to the algorithm, the control parameters were then determined to be Ttarget = 41.9 °C and Tf = 33.6 °C. With these two parameters, the ablation geometry was calculated using the model described above. The ablation depth was found to be 0.78 mm, and the protection was found to be 0.06 mm. As the intimal thickness increased (1.02 mm in ablation depth and 0.1 mm in protection depth), as shown in case 2 (Figure 11b), Ttarget = 29.9 °C and Tf = 8.6 °C were obtained with the algorithm, and there were errors of 0.01 mm in both ablation depth and protection depth.

Discussion
Based on the experimental study, it was found that the temperatures at the three monitoring points (T1, T2, and T3) inside the tissue, representing the locations of plaque, adventitia, and surrounding tissue [28], changed linearly with the three control parameters. Moreover, the two parameters Ttarget and Tf were found to be the main factors that adjusted the ablation results, including the ablation depth and protection depth.

Discussion
Based on the experimental study, it was found that the temperatures at the three monitoring points (T1, T2, and T3) inside the tissue, representing the locations of plaque, adventitia, and surrounding tissue [28], changed linearly with the three control parameters. Moreover, the two parameters T target and T f were found to be the main factors that adjusted the ablation results, including the ablation depth and protection depth. However, the RF volumetric heating power and the convection cooling power were simultaneously affected by T target and T f . The ablation depth and protection depth changed synchronously with a single control parameter. Only through the simultaneous control of the two parameters could bidirectional control of confined penetrating heating be realized.
Current RF electrodes and the cooling agent setup were found to be able to cover an ablation depth ranging from 0.47 mm to 1.43 mm, with a protection depth ranging from 0 mm to 0.26 mm, which is sufficient for most arterial plaques [40,41]. Further explorations into increasing this range would help broaden its application in clinical settings.
The proposed algorithm combining bilinear interpolation and iteration was used to obtain the optimal combination of T target and T f for the targeted ablation depth and protection depth, which is essentially an inverse problem of heat transfer [42,43]. In our study, this was simplified as a multivariable optimization problem, and the results of the numerical calculations were used to simplify the process of determining physical relations. The algorithm illustrated the feasibility of using the two proposed parameters for the bidirectional control of the arterial plaque treatment. Though there were errors of 0.01 mm in control for both the ablation depth and the protection depth, as long as the resulting protection region is larger than the objective, with a sufficient ablation depth, this should be clinically acceptable.

Limitation
There are certain limitations to this study. Though a one-way sensitivity analysis verified the robustness of the current model, a model that takes into account the differences in the electrical and thermal properties of heterogeneous plaques would give more precise predictions, especially for advanced atherosclerotic plaques, which include lipids and calcified regions [44,45]. The feasibility of the bidirectional control strategy using two main control parameters, T target and T f , was proven with numerical calculations. Further precision imaging of the ablation, together with careful experimental design, is needed.

Conclusions
The temperature of the control point on the inner surface of the endothelial layer (T target ) and the temperature of the cooling flow (T f ) were found to be two effective factors for adjusting the ablation results in the newly proposed radiofrequency balloon angioplasty.
The two parameters have contrary effects on the control of the ablation depth and protection depth when used alone. An ablation range covering most clinical plaque treatment requirements can be attained through the simultaneous use of the two parameters. The feasible combinations of the conditions were found. A bidirectional penetrating heating strategy can be realized with the proposed bilinear interpolating and iterating algorithm. Future studies to widen the control range while considering the complexity of real plaque properties would help with the technique's further application.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The datasets created and/or analyzed during the current investigation are available upon reasonable request from the corresponding author. All figures in this paper are original.

Acknowledgments:
The authors would like to thank to their colleagues in the Bioheat and Mass Transfer Laboratory for their valuable discussions and advice in improving the writing.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: RF radiofrequency PID proportional-integral-derivative T target the target temperature T f the temperature of cooling water  V f  the velocity of cooling water  PTA  percutaneous transluminal angioplasty  SMCs  smooth muscle cells  IMT  intima-media thickness  DAQ  data acquisition card  TC0, TC1, TC2, and TC3 thermocouples on points TC0, TC1, TC2, and TC3  T0, T1, T2, and T3  temperatures measured by thermocouples TC0, TC1, TC2, and TC3  SAR specific absorption rate A ablation depth, unit: millimeter P protection depth, unit: millimeter A k calculated ablation depth, unit: millimeter P k calculated protection depth, unit: millimeter where Nu f is the Nusselt number of cooling water, λ is the thermal conductivity of cooling water, and D is the catheter diameter. Re is the Reynolds number of cooling water; f is the Darcy resistance coefficient; Pr f is the Prandtl number; and ρ, u, µ, v, and a are the density, velocity, coefficient of viscosity, kinematic viscosity, and thermal diffusion coefficient of the cooling water in the catheter.      ± 0.11V) and 70.5 Ω (615.9 ± 3.7 Ω vs. 545.4 ± 0.4 Ω) for the applied voltage and measured impedance. These differences may have been due to the following reasons: 1) There was a difference in the applied voltage between the simulation and the experiment. The voltage adjusted by the PID controller had a linear relationship with e(k) in the simulation, while the RF voltage output by the RF generator may have had a slightly nonlinear relationship with e(k) in the experiment. 2) There was a difference in the position of the control point between the simulation and experiment. Due to measurement precision, the position of the temperature control point (TC0) in the simulation and the experiment may have had some differences. 3) There was a difference in electrical and thermal conductivity between the experiment and simulation. In the model, the properties of the artery wall were used, and there were some differences between the artery and the tissue-mimicking phantom used in the experiment, though the electrical conductivity of the tissue-mimicking phantom was adjusted to match the electrical conductivity of the arterial wall. At the same time, the impedance interface was neglected in our model but was present in the experiment and brought a higher impedance than that in the simulation.
Considering the real output of the RF generator, positioning precision, the difference in the electrical conductivity, and the similar trend of the matching curve between the numerical and experimental results, these deviations are acceptable.  The deviations between the numerical results and the experimental results were 4.26 • C (59.95 ± 0.09 • C vs. 55.69 ± 0.14 • C), 3.99 • C (46.68 ± 0.12 • C vs. 42.69 ± 0.07 • C), and 1.61 • C (40.98 ± 0.18 • C vs. 39.37 ± 0.08 • C) for T1, T2, and T3, and 0.98 V (24.24 ± 0.23 V vs. 25.22 ± 0.11 V) and 70.5 Ω (615.9 ± 3.7 Ω vs. 545.4 ± 0.4 Ω) for the applied voltage and measured impedance. These differences may have been due to the following reasons: (1) There was a difference in the applied voltage between the simulation and the experiment. The voltage adjusted by the PID controller had a linear relationship with e(k) in the simulation, while the RF voltage output by the RF generator may have had a slightly nonlinear relationship with e(k) in the experiment. (2) There was a difference in the position of the control point between the simulation and experiment. Due to measurement precision, the position of the temperature control point (TC0) in the simulation and the experiment may have had some differences. (3) There was a difference in electrical and thermal conductivity between the experiment and simulation. In the model, the properties of the artery wall were used, and there were some differences between the artery and the tissue-mimicking phantom used in the experiment, though the electrical conductivity of the tissue-mimicking phantom was adjusted to match the electrical conductivity of the arterial wall. At the same time, the impedance interface was neglected in our model but was present in the experiment and brought a higher impedance than that in the simulation.
Considering the real output of the RF generator, positioning precision, the difference in the electrical conductivity, and the similar trend of the matching curve between the numerical and experimental results, these deviations are acceptable.     T1_finer mesh T2 T2_finer mesh T3 T3_finer mesh Figure A5. Heating curves at different points (T0, T1, T2, and T3) to compare the convergence of the models with different mesh densities. The mesh elements increased from 380,825 to 1,381,825.