Computational simulation of the predicted dosimetric impact of adjuvant yttrium-90 PET/CT-guided percutaneous ablation following radioembolization

Background 90Y PET/CT post-radioembolization imaging has demonstrated that the distribution of 90Y in a tumor can be non-uniform. Using computational modeling, we predicted the dosimetric impact of post-treatment 90Y PET/CT-guided percutaneous ablation of the portions of a tumor receiving the lowest absorbed dose. A cohort of fourteen patients with non-resectable liver cancer previously treated using 90Y radioembolization were included in this retrospective study. Each patient exhibited potentially under-treated areas of tumor following treatment based on quantitative 90Y PET/CT. 90Y PET/CT was used to guide electrode placement for simulated adjuvant radiofrequency ablation in areas of tumor receiving the lowest dose. The finite element method was used to solve Penne’s bioheat transport equation, coupled with the Arrhenius thermal cell-death model to determine 3D thermal ablation zones. Tumor and unablated tumor absorbed-dose metrics (average dose, D50, D70, D90, V100) following ablation were compared, where D70 is the minimum dose to 70% of tumor and V100 is the fractional tumor volume receiving more than 100 Gy. Results Compared to radioembolization alone, 90Y radioembolization with adjuvant ablation was associated with predicted increases in all tumor dose metrics evaluated. The mean average absorbed dose increased by 11.2 ± 6.9 Gy. Increases in D50, D70, and D90 were 11.0 ± 6.9 Gy, 13.3 ± 10.9 Gy, and 11.8 ± 10.8 Gy, respectively. The mean increase in V100 was 7.2 ± 4.2%. All changes were statistically significant (P < 0.01). A negative correlation between pre-ablation tumor volume and D50, average dose, and V100 was identified (ρ < − 0.5, P < 0.05) suggesting that adjuvant radiofrequency ablation may be less beneficial to patients with large tumor burdens. Conclusions This study has demonstrated that adjuvant 90Y PET/CT-guided radiofrequency ablation may improve tumor absorbed-dose metrics. These data may justify a prospective clinical trial to further evaluate this hybrid approach.


Background
Tumor targeting in yttrium-90 ( 90 Y) radioembolization differs from other radionuclide therapies that are infused systemically and find their targets through high affinity to cellular receptors. Instead, the distribution of 90 Y within a tumor depends strongly on the catheter position during infusion, downstream fluid dynamics, and arterial perfusion in both tumor and uninvolved liver [1,2]. Due to the mechanical nature of 90 Y microsphere trapping, dose nonuniformities within tumor have been demonstrated through histological analyses [3][4][5][6] and 90 Y PET/CT imaging [7][8][9][10]. While several metrics have been used to predict response following radioembolization [11], average tumor absorbed dose (D avg ) is the most common metric. However, non-uniform deposition of 90 Y may result in a sub-therapeutic dose to portions of the tumor, which has been shown to correlate with poor response even if D avg is favorable [9].
Like many liver-directed therapies, 90 Y radioembolization is commonly classified as a palliative treatment due to the relatively poor prognosis of patients with liver cancer. However, this does not suggest that there is no utility in treatment optimization. Response following locoregional hepatic therapy has been shown to correlate with improved patient survival, prompting the use of multimodality therapies [12,13] to improve tumor response. One such example of multimodality therapy is combined trans-arterial chemoembolization (TACE) and radiofrequency ablation (RFA) for the treatment of hepatocellular carcinoma (HCC) [12]. Since the widespread use of post-radioembolization 90 Y PET/CT, there has been an interest in utilizing this imaging data to improve patient therapy. While multimodality treatment utilizing 90 Y radioembolization has not been widely studied, several authors have attempted the use of 90 Y PET/CT to provide multistage patient-specific 90 Y treatments with positive outcomes [14,15].
Image-guided percutaneous ablation is the standard minimally invasive treatment for eligible patients with small (<3 cm) liver tumors [16]. However, the efficacy of ablation decreases with increasing tumor size. For patients with HCC, RFA generally results in poor outcomes for tumors greater than 5 cm in size, where alternative treatments such as 90 Y radioembolization are commonly employed [17]. While combined RFA and TACE has demonstrated improved response in treating larger tumors [12], the utility of combining percutaneous ablation with 90 Y radioembolization has not yet been evaluated. 3D 90 Y PET/CT-based dosimetry and dose-response thresholds [18,19] may be used in principle as a guide for percutaneous ablation of under-treated regions of tumor, potentially improving therapy.
In this work, 90 Y PET/CT has been used as a guide for simulated modeling of adjuvant RFA. Improvement in tumor absorbed-dose metrics in the unablated tumor has been calculated to determine if this hybrid technique warrants further evaluation in a prospective clinical trial.

Methods
The study was carried out in compliance with the Declaration of Helsinki and has been approved by the institutional review board at each participating site. Written consent was obtained from patients at site A (University of Tennessee IRB #3502), while a waiver of informed consent was obtained for data collected at site B (Wright State University IRB #SC6291). Patients treated with 90 Y radioembolization using resin [site A, SIR-Spheres®, SIR-Tex Medical Ltd, North Sydney, Australia] or glass [site B, Therasphere® BTG, London, UK] microspheres for primary or secondary liver cancer and who received post-treatment 90 Y PET/CT imaging were reviewed. Fourteen patients (age 49-80 years) met the inclusion criteria of tumor D70 < 100 Gy [9] for resin microspheres (D70 < 150 Gy for glass) with visualized areas of decreased microsphere uptake >2 cm in size identified on 90 Y PET/CT, where D70 is the minimum absorbed dose (Gy) to 70% of the tumor volume. These patients received radioembolization for treatment of hepatocellular carcinoma (HCC, n = 10), intrahepatic cholangiocarcinoma (n = 2), or liver dominant metastatic disease (n = 2). The median model for end-stage liver disease score was 8, with 10 Child-Pugh class A patients and 4 class B patients. Additional data are available in Table 1. , and no filter. 90 Y imaging performance of these PET/CT scanners at activity concentrations common in radioembolization with resin microspheres has been evaluated previously [20]. All patients received post-radioembolization 90 Y PET/CT imaging on the day of treatment.
Under the supervision of a dual-board-certified nuclear medicine radiologist, tumors were contoured in three dimensions on 90 Y PET/CT referencing appropriate pretreatment hepatic protocol CT or MRI in patients with HCC or 18 FDG PET/CT images in patients with metastatic disease. Contours were drawn using Osirix 7.5 [Pixmeo, Bernex, CH], exported to Matlab 2014b [Mathworks, Natick, MA], and converted to a 3D mask with the same pixel size as 90 Y PET/CT data. For each tumor, 90 Y PET/CT-based 3D dosimetry was performed using the local deposition method [21] and dose-volume histograms were generated. The D50, D70, D90, maximum and average tumor dose (D avg ), and V100 were calculated for each tumor since these metrics were previously validated as prognostic indices for radioembolization [9], where D50 and D90 are the minimum doses to 50 and 90% of tumor volume and V100 is the percentage of tumor volume receiving more than 100 Gy.

Radiofrequency ablation simulation
Because RFA is the most common treatment in the USA for small tumor ablation in situ [16], it was selected as the simulated ablation modality in this study. While in vivo and ex vivo RF temperature profiles and ablation zones have been published previously, we elected to use biophysical modeling which allowed for flexible calculation of ablation zones for varied ablation times, electrode angles, positions, and tumor tissue characteristics (HCC or metastatic disease) which could be defined in a 3D voxel space matching that of 90 Y PET/CT. This model was compared with data previously reported from similar models [22] and in vivo measurements [23] to establish its accuracy under a fixed set of conditions.
The model consisted of a water-cooled 17-gauge straight RF electrode with a 30-mm active tip. The electrode model was placed inside a 20 × 20 × 20 cm block of liver tissue, with electrothermal properties discussed later in this section. A 3D tetrahedral mesh was applied to the model for finite element method (FEM) analysis of the heat transport problem. The maximum mesh element size was restricted to 8 mm in tissue, 1 mm in the insulating portion of the electrode, and 0.25 mm in the active electrode tip. The electrode model and tetrahedral mesh are shown in Fig. 1. The FEM was used to solve Penne's bioheat transport equation [24], which consists of two partial differential equations: where ∇ is the Laplace operator, σ is the electrical conductivity (S/m), V is the electric potential (volts), ρ is tissue material density (Kg/m 3 ), c is tissue specific heat (J/Kg/°C), k is the electrical conductivity of tissue (W/m/°C), w b is the perfusion rate of blood per unit volume of tissue (m 3 /m 3 /s), and c b is the specific heat of blood (J/Kg/°C). T a is the temperature of arterial blood, held constant at 37°C, and T is the temperature as a function of time, t. The initial conditions for all geometries at all nodes (tissue and electrode) were set to 37°C, with liver boundary conditions maintained at 37°C throughout the simulation. The potential (V) was modulated using a simulated feedback circuit to maintain the temperature at the hottest tissue node to <105°C. The perfusion rate of blood per unit volume in tissue (w b ) was varied during simulation to account for vascular coagulation [25]. The first-order kinetic Arrhenius model [25] was employed to estimate fractional vascular coagulation  (Table 2) with RF ablation cell killing boundary (Ω = 6.9) axial and lateral range indicated (right) as a function of time and temperature during ablation, as described in Eq. 3.
C f is the fraction of capillary flow that has been occluded, and R is the ideal gas constant (8.314 J/mol/°K). A is the Arrhenius pre-exponential factor, a measure of molecular collision frequency (1.98e106 1/s), and E a is the activation energy, the amount of energy required to transform molecules from their original state to a damaged state (6.67e5 J/mol). T(t) is the temperature as a function of time. Following a previous investigation [25], W b was linearly decreased with increasing C f (Eq. 4) to simulate decreasing blood perfusion as RFA creates localized coagulation.
In Eq. 4, w b,nc is the perfusion rate of blood in the absence of any heat-induced vascular coagulation. The Penne's heat transport equation [24] and similar approximations have been previously used for modeling of radiofrequency ablation in tissue, with several examples [22,26] reviewing additional relevant details.
The RFA electrode model, tetrahedral mesh, and FEM analysis of the Penne's bioheat equation were simulated as described using COMSOL (FEMLAB) 5.0 [COMSOL, Stockholm, Sweden]. Electrode angle, depth, ablation time, and tissue parameters were predefined through a Matlab interface with COMSOL for each location simulated. After each simulation, continuous temperature data as a function of time was re-binned in Matlab into a 3D spatial map with voxel sizes matching radioembolization data from quantitative 90 Y PET/CT. Calculations were performed on a 6th-generation Intel Core-i7 system with 32 Gb of memory.

Tissue parameters
Tissue parameters (σ, ρ, c, k, c b , w b,nc ) were obtained from the literature for normal liver, cirrhotic liver, and HCC (Table 2). Tissue parameters for normal liver were used for the four patients treated for metastatic disease or cholangiocarcinoma (Table 1) since no data specific to these tumor types were available. Data for HCC were used for simulation in the remaining patients. For reference, data for cirrhotic liver, with characteristically decreased w b,nc , is also shown in Table 2.

Y PET/CT-Guided RFA
The tumor destruction boundary corresponding to each ablation location was calculated using the Arrhenius model for cell killing from FEM calculated timedependent temperature (Eq. 5).
When applied to the prediction of thermal cell death, the Arrhenius pre-exponential factor (A) is 2.984e80 1/s, the activation energy (E a ) is 5.06e5 J/mol, and SF is the surviving fraction of cells. An Ω value of 6.9 was selected, corresponding to a SF of 0.001 or 99.9% cell killing within the ablation zone.
Under the supervision of a board-certified interventional radiologist, up to four electrode locations were selected in each patient, with ablation time up to 600 s per location. Electrode locations were selected using realistic percutaneous access paths, as would be performed using conventional clinical techniques [16] under CT or CT fluoroscopy guidance, with the acquired 90 Y PET/CT as a guide. Time and position were varied to limit the ablation zone to areas of tumor receiving <100 Gy when possible and, in cases where the ablation zone was close to the liver capsule or gallbladder, to spare these tissues. The ablation zone boundary was maintained at least 1 cm from nearby gastrointestinal tissue.

Dosimetry with and without RFA
Following calculation of 3D ablation zones for each location, the Ω = 6.9 threshold was used to create a 3D mask of the absolute ablation boundary. This mask was subtracted from the pre-defined tumor mask, leaving unablated tumor, treated only with 90 Y radioembolization. Dose-volume histograms (DVH), D50, D70, D90, D avg , and V100 were recomputed using the same 3D dosimetric dataset (less ablated tissue) to compare 90 Y radioembolization with adjuvant PET/CT-guided ablation to radioembolization alone in this patient cohort.

Statistical analysis
The Kolmogorov-Smirnov test was used to assess the assumption of normality for all datasets evaluated. Differences in tumor dose metrics between 90 Y radioembolization and 90 Y radioembolization with PET/CT-guided ablation were evaluated using a paired-sample T test. The potential correlation between both tumor size and number of ablation sites on post 90 Y PET/CT-guided ablation dose metrics was assessed using the Spearman's correlation statistic. Linear regression was used to explore this relationship if Spearman's statistic suggested a correlation (|ρ| > 0.5). A P value <0.05 was considered statistically significant through all analysis.

Results
Thirty-three RFA simulations were performed for the adjuvant simulated therapy of 15 tumors. The mean computation time per simulation was 7.8 min. Simulated 90 Y PET/CT-guided RFA resulted in a calculated decrease in active tumor volume when cell killing inside the Ω = 6.9 boundary was assumed. The mean pre-and post-ablation tumor volume was 197.2 ± 102.3 cm 3 and 179.3 ± 98.7 cm 3 , respectively. The average percent decrease in tumor volume was 12.1 ± 7.9%. Spearman's statistic showed a moderate positive correlation between the number of ablation sites and the absolute change in tumor volume (ρ = 0.59, P = 0.02). Ablation and volumetric data for each tumor is shown in Table 3. Lack of consistency in the volumetric effect of ablation is secondary to varying overlap of ablation zones or ablation zones not fully contained within tumor boundaries.  Table 1) are shown in Figs. 2 and 3. Adjuvant simulated ablation had the effect of modifying the shape of the DVH, reducing the percentage of tumor volume receiving lower absorbed doses in all cases. However, the effect was more pronounced for smaller tumors than larger tumors, as illustrated in Fig. 4a, b for patients 1 and 5, respectively.
In addition to affecting DVH shape, simulated 90 Y PET/ CT-guided RFA resulted in statistically significant increases in the D avg , D50, D70, D90, and V100 compared  to 90 Y radioembolization alone. The mean absolute increase in D avg was 11.2 ± 6.9 Gy (P < 0.001). Mean increases in D50, D70, and D90 were 11.0 ± 6.9 Gy (P < 0.001), 13.3 ± 10.9 Gy (P < 0.001), and 11.8 ± 10.8 Gy (P < 0.01), respectively. The mean increase in V100 was 7.2 ± 4.2% (P < 0.001). Additional calculated dosimetric details comparing radioembolization alone and radioembolization with simulated RFA are given in Table 4. The limited volume of tumor that can be affected by RFA is reflected in the difference in the percent change in D50, D70, and D90. For example, D50 was increased by an average of 11.2%, D70 by 18.1%, and D90 by 43.8%. Since D90 is the minimum dose to 90% of the tumor volume, it shows the most drastic difference since adjuvant RFA targets a relatively small fraction of tumor volume receiving the lowest dose.

Discussion
The distribution of 90 Y radioembolization within the tumor can only be partially controlled by the treating physician, i.e., by careful microcatheter placement and through modification of downstream fluid dynamics. It is logical, therefore, to combine 90 Y radioembolization with an adjuvant treatment modality that can be more precisely targeted, such as percutaneous ablation. Compared to RFA alone, there is substantial evidence that multimodality treatments such as combined TACE-RFA can result in improved survival in patients with intermediate size HCC lesions [12,27]. While there are obvious similarities between TACE-RFA and combined radioembolization-RFA, it is the differences that warrant discussion within the context of the findings in this manuscript. The primary utility of 90 Y PET/CT is that it allows for proactive planning of alternative or adjuvant therapies immediately following radioembolization [14]. While not fully established, published absorbed-dose thresholds [18,19] and quantitative 90 Y PET/CT can be used to plan secondary therapies, potentially including the percutaneous ablation of areas of tumor receiving low absorbed dose. By the same token, areas of tumor receiving high absorbed doses can be spared superfluous secondary intervention. Unlike combined TACE-RFA, strict quantification of the change in absorbed dose, a physical quantity related to tumor response, can be predicted prior to adjuvant ablation following radioembolization-as has been done in this work. To this end, this study has demonstrated statistically significant improvements in both unablated tumor absorbed-dose metrics and DVH shape after simulated ablation following 90 Y treatment in a 14-patient cohort.
The primary endpoint of this effort was successful in elucidating that improvement in unablated tumor dose metrics following adjuvant ablation appear possible; however, the clinical significance of this is unknown. In addition, neither the technical feasibility nor the complication rate associated with these two therapies in tandem has been evaluated. Because percutaneous ablation is highly targeted and spares normal liver tissue, however, it is likely that complication rates will be low. Radioembolization-RFA can be contrasted in theory to patients who received both external beam radiation therapy and 90 Y radioembolization-a tandem treatment with a notable increase in hepatic toxicity [28].
The relationship between improved tumor absorbeddose metrics and improved response has been demonstrated in 90 Y radioembolization [9,18]. However, there is sufficient variability in the dose-response thresholds reported in the literature that a prediction of improved treatment efficacy in the patient cohort analyzed in this manuscript cannot be made. The effect of combined radioembolization and percutaneous ablation on tumor  (Table 1), pre-treatment tumor volume = 124.1 cm 3 . b Patient 5 (Table 1), pre-treatment tumor volume = 374.6 cm 3 response must, therefore, be evaluated as part of a future clinical trial. Such a trial focusing on outcomes would also allow exploration into situations where the tumor as a whole was grossly under-treated, but without substantial inhomogeneity. While adjuvant ablation in the case of homogeneous tumor dose would not result in a substantial change in dose metrics, outcomes may still be affected due to the concomitant biological effects of both treatment modalities.
One limitation of this study is that the only percutaneous ablation modality simulated was RFA. While RFA was selected due to its wide use in hepatic tumor ablation [16], radioembolization combined with other ablation modalities may be associated with different results. For example, the mean difference in pre-and post-RFA tumor volume among all patients was 18.9 cm 3 (Table 3). This relatively small volume is not only due to lower ablation times used to target only under-treated tumor in some patients but also due to the relatively small RFA ablation zone size. In addition, the correlation identified between tumor volume and absolute change in several tumor absorbed-dose metrics (ΔD50, ΔD70, ΔD avg , ΔV100) can be at least partially attributed to RFA ablation zone size. Alternatives such as microwave ablation are associated with larger ablation volumes may result in further improvement in unablated tumor absorbed-dose metrics if used following radioembolization.
A second limitation is that electrothermal tissue characteristics and in particular, w b,nc , can vary substantially and have a marked impact on calculated and actual RF ablation zone size [22]. To this end, it is likely that the model employed in this work, utilizing w b,nc for normal liver and HCC (Table 2), resulted in an underestimation of the ablation zone size. One clinically well-known phenomenon in the thermal ablation of small HCC tumors in patients with cirrhosis is the "oven effect" [29] in which the decreased blood perfusion rate (w b,nc ) in cirrhotic liver (Table 2) results in an increased ablation zone size. Since blood perfusion carries heat away from the electrode, it ultimately limits the ablation zone size regardless of the modality of thermal ablation used. However, in 90 Y PET/CT-guided RFA, patients receive radioembolization with glass or resin microspheres prior to adjuvant ablation. The embolic nature of the microspheres, resin greater than glass, may decrease tumor blood perfusion and, therefore, could increase ablation zone size beyond that simulated in this work. The oven effect may be contributory to increased ablation zone size in RFA-TACE [12]; however, the effects following 90 Y radioembolization cannot be determined at this preclinical stage.
Finally, this study has not considered the effects of respiratory motion on the 90 Y PET/CT acquisition, tumor boundary delineation, calculated dose metrics, or the potential effect of respiratory motion on clinical electrode placement. Several studies [30,31] have described changes in tumor volume and SUV in 18 FDG PET/CT when respiratory gating was used. However, investigation into the feasibility of respiratory gating in 90 Y PET/CT has been limited to a few examples [32,33]. The usefulness of 90 Y PET respiratory gating may be limited by the small branching ratio of positron emission in 90 Y [34] and the clinical difficulty of increasing an already long acquisition time. Respiratory motion undoubtedly had an effect on the tumor dose metrics reported for each patient in this manuscript; however, the extent of which could not be quantified. As this effort proceeds into clinical investigation, more work into the potential use of amplitude based gating in 90 Y PET [32] may be warranted and, at the very least, physicians must be aware of the potential effects of motion on 90 Y PET data when planning adjuvant therapy. Although this work precedes a clinical trial, the timing of radioembolization and subsequent ablation will be of concern and warrants brief discussion. 90 Y radioembolization is a permanent implant, delivering 97.5% of its absorbed dose to tissue in the first 2 weeks after infusion. However, because there is no biological removal of 90 Y, the absorbed dose is committed immediately after the microspheres are infused. In theory, the timing, therefore, between 90 Y treatment and ablation is not critical so long as ablation is performed before structural tumor changes occur creating deviation with the ablation guidance model, i.e., the post-treatment 90 Y PET/CT. One additional concern is whether heating of microspheres in situ could result in the release of free 90 Y-a complication that must be investigated with in vitro experiments and patient bioassay before and during a future clinical trial. However, even if systemic release is found to be a potential issue, waiting two or more weeks before adjuvant ablation will allow nearly all the infused radioactivity to decay, minimizing safety concerns.

Conclusions
This study has demonstrated that adjuvant 90 Y PET/ CT-guided ablation may improve tumor absorbed-dose metrics. These data may justify a prospective clinical trial to further evaluate this hybrid approach.
Authors' contributions ASP, LKF, and YCB were involved in the study design, implementation, analysis, and manuscript preparation. AL designed and prepared the software for the 90 Y/RFA treatment-planning simulator and assisted in the manuscript preparation. CG assisted in the simulated RFA probe placement for ablation and the manuscript preparation. SK assisted in the manuscript preparation and patient identification. All authors read and approved the final manuscript.  Linear regression for dose metrics vs tumor volume. a ΔD50 vs tumor volume (Slope fit = −0.046 Gy/cm 3 , R 2 = 0.505, P = 0.003). b ΔD avg vs tumor volume (Slope fit = −0.04 Gy/cm 3 , R 2 = 0.354, P = 0.02). c ΔV100 vs tumor volume (Slope fit = −0.034%/cm 3 , R 2 = 0.615, P = 5.4e−4)