Radiomics Model Based on Gadoxetic Acid Disodium-Enhanced MR Imaging to Predict Hepatocellular Carcinoma Recurrence After Curative Ablation

Introduction

As the fourth most common cause of cancer-related death worldwide,¹ hepatocellular carcinoma (HCC) is a heavy disease burden, especially in Eastern Asia and sub-Saharan Africa where 85% of HCC cases are estimated to occur.² The five-year survival rate of patients with HCC is approximately 18%,³ but significant regional differences in the prognosis of HCC are observed due to differences in the availability of surveillance programs and effective treatments.^4,5 In China, although the period between 2000 and 2014 saw a decline in the incidences and mortality rates of HCC, there were still 364,800 new HCC cases and 318,800 liver cancer deaths in 2014.⁶

Various treatment modalities have been developed for HCC based on the extent of disease and the severity of underlying liver disease. For patients with localized HCC that are ineligible for resection or transplantation without macrovascular invasion, ablation is an important therapy option. Ablation includes microwave ablation (MWA) and radiofrequency ablation (RFA), which have been used to treat HCC for years in China. These techniques have some theoretical advantages.^7,8 Results from randomized controlled trials and meta-analyses suggest that the results of MWA are comparable with those of RFA.^9,10

However, recurrence after curative ablation remains a concern and influences subsequent therapy.^11,12 Given that an early diagnosis and effective treatment are associated with survival of patients with HCC,^13,14 a tool to predict the prognosis of patients with HCC who underwent curative ablation is needed. Although various models have been developed to predict prognosis of patients with HCC,^15–17 few of these models were specifically designed for patients with HCC who have undergone curative ablation, and most available models are generally of a low quality.¹⁸ A possible reason for this is that most of the prediction models are based on clinical and initial pathological predictors, such as alpha-fetoprotein (AFP), tumor size, vascular invasion, tumor number, tumor differentiation, and neutrophil–lymphocyte ratio.¹⁸

To further improve the performance of prediction models for HCC prognosis, novel predictors should be explored, such as gene signatures¹⁹ and non-coding RNA signatures.²⁰ Radiomics is a promising predictor for developing HCC prognosis prediction models.^21,22 One of the advantages of radiomics is that it can extract quantitative radiology information and assess prognosis in a non-invasive manner. However, radiomics-based prognosis prediction models are limited for patients with HCC after MWA and RFA. Currently, available radiomics-based models are of poor quality; thus, they have not yet reached clinical utility.²³ Therefore, the present study aimed to develop a magnetic resonance imaging (MRI)-based radiomics model to predict recurrence of HCC after curative ablation.

Materials and Methods

Study Design and Study Population

The study adopted a retrospective cohort study design. A flowchart of the study design is presented in Figure 1. A total of 132 patients with HCC who were admitted to Sun Yat-sen Cancer Center between June 2012 and April 2018 and underwent curative ablation after being evaluated by a multidisciplinary team were enrolled. Curative ablation (either RFA or MWA) was conducted by experienced surgeons or interventional doctors guided by ultrasound and computed tomography. Detailed procedures are described below.

Figure 1 Flowchart of study design. Three-dimensional (3D) segmentation was performed on 3D serial raw magnetic resonance imaging scans within 1 month before microwave ablation. Radiomic features were extracted from regions of interest, which were outlined around the lesion. A random survival forest algorithm was employed for radiomic feature selection and model construction. The predictive performance of constructed models was evaluated by receiver operating characteristic (ROC) curves, time-dependent ROC curves, decision curves, and Kaplan–Meier survival curves. Abbreviations: GLCM, gray-level co-occurrence matrix; GLRLM, gray-level run-length matrix; GLSZM, gray-level size zone matrix; NGTDM, neighboring gray tone difference matrix; GLDM, gray-level dependence matrix; LoG, logarithm.

Detailed inclusion criteria are as follows: 1) clinically diagnosed HCC following non-invasive criteria defined by the American Association for the Study of Liver Disease based on specific imaging features;²⁴ 2) the longest diameter of the largest lesion >10 mm; 3) no capsular invasion or invasion of adjacent organs; 4) no extrahepatic metastasis at initial diagnosis; 5) no tumors in the hepatic veins or main portal vein; 6) data from contrast-enhanced MRI (contrast agent, gadoxetic acid disodium) within 1 month before ablation. Patients who met the following exclusion criteria were excluded: 1) follow-up time <2 years after ablation; 2) radiotherapy, chemotherapy, and interventional therapy before ablation; 3) coexistence of other malignant tumors; 4) severe images artifacts.

Patients were randomly divided into a training cohort (n = 92) and a validation cohort (n = 40) at a ratio of 7:3. Data from the training cohort were used to construct the prediction models, and data from the validation cohort were used for internal validation.

This study was approved by the institutional review board of Sun Yat-sen University Cancer Center, and the requirement for informed consent was waived.

Ablation Procedure

Patients remained in a supine position on the computed tomography table, and 5–15 mL of 1% lidocaine was used to induce local infiltration. An antenna was percutaneously inserted into the tumor at the end of expiration under computed tomography guidance. Axial computed tomography was performed to ensure that the antenna passed through the marked center of the tumor’s maximal circle in the axial plane after puncture, and the antenna tip reached the superior-most axial computed tomography slice that contained the tumor. When necessary, the antenna was adjusted slightly. After puncture, 1–2 mg/kg of intravenous propofol (Diprivan, Zeneca, Macclesfield, UK) and 50–100 μg fentanyl (Fentaini, Renfu, Yichang, China) was used to induce general anesthesia.

For MWA, a power of 45 W for 5–10 minutes per MWA application was used, and the time of MWA was based on lesion necrosis observed by computed tomography during the procedure (emission facility: frequency, 2450 MHz; output power, 0–120 W; precision of temperature control, ±0.1°C; discharge waveform, continuous wave; produced by Qi Ya Medical Treatment Facility Limited Company, Nanjing, China). Overlapping ablation was performed for tumors with a relatively irregular shape or after failure to completely cover the tumor by increasing the ablation time.

For RFA, an alternating current generator (RF3000, Boston Scientific, Boston, MA, USA) and an electrode needle with an insulated 15-G outer needle, which housed 10 solid retractable curved electrodes with a diameter of 3.5 cm when expanded like an umbrella, were used. A second application of ablation was performed when a marked increase in impedance was not observed. Multiple overlapping ablation rounds were performed when the greatest tumor dimension was >3.0 cm.

Study Outcome

Recurrence-free survival (RFS) was the study outcome. To determine the occurrence of RFS, all patients were followed up for at least 2 years from the date of curative ablation with regular examination of serum AFP and contrast-enhanced computed tomography or MRI every three months, until the first occurrence of RFS. Patients who were lost to follow-up or who died due to other reasons were censored.

Clinical Characteristics and Selection of Predictors

We retrospectively collected the baseline clinical characteristics of the study population, including age, sex, hepatitis B virus surface antigen (HBsAg) positivity, albumin (ALB), alkaline phosphatase, alanine aminotransferase, aspartate aminotransferase, direct bilirubin, gamma-glutamyl transferase (GGT), indirect bilirubin, total bilirubin, AFP, tumor size, and number of tumors.

To select factors associated with RFS for model construction, univariable Cox regression was employed to assess the association of each clinical factor with RFS. Factors with a p-value of <0.1 were further included in the multivariable Cox regression analysis, and a backward stepwise selection strategy was used based on the Akaike information criterion.

MRI Examination and Extraction of Radiomic Features

All patients underwent MRI within 1 month before curative ablation. Detailed MRI parameters, including acquisition sequence and type of scanner, are presented in Supplementary Table 1. In brief, all transverse images were obtained with a field of view of 44 × 33 cm or 40 × 30 cm. Gadoxetic acid disodium (Primovist, Bayer Schering Pharma, Berlin, Germany, 0.1 mL/kg) and 0.9% saline chaser (20 mL) were administered as a bolus at a flow rate of 1 mL/s. After injection of contrast agent, dynamic enhanced images in arterial, portal venous, and transitional phases were obtained, while additional hepatobiliary phase images were obtained with a 20-minute delay. T2-weighted fat suppression images were obtained using a navigator-triggered technique. The imaging appearance of the lesions is listed in Supplementary Table 2.

To extract tumoral radiomic features, a radiologist (L.Z.) with 10 years of experience in diagnostic liver imaging manually segmented the lesion volume of interest (VOI) in a three-dimensional manner using ITK-SNAP (version 3.8.0, http://www.itksnap.org/pmwiki/pmwiki.php) in T2 fat suppression, T1 fat suppression, hepatic arterial phase, portal venous phase, and hepatocellular phase sequences (20-minute delay after administration of contrast agent). In cases of multiple lesions, the largest lesion was chosen as the representative lesion to segment.

To extract peritumoral radiomics, the peritumoral VOI was also segmented using a Python morphology dilation algorithm, which automatically dilated the boundaries of each lesion by 5 mm. A peripheral ring was also generated with automated dilatation and shrinkage of the tumoral boundaries by 5 mm on each side, namely outside and inside the boundary, resulting in a ring with a thickness of 10 mm (5 + 5 mm). The portion beyond the liver parenchyma was removed manually if the VOI was beyond the liver parenchyma after expansion.

Before radiomic feature extraction, to reduce potential image differences caused by different MRI systems, all images were resampled to a voxel size of 1 × 1 × 1 mm³. We used the open-source package Pyradiomics (version 2.12; https://pyradiomics.readthedocs.io/en/2.1.2/) to extract radiomic features, including tumoral radiomic (T), peritumoral radiomic (PT 5 mm), and peritumoral radiomic (PT 5 + 5 mm) features.

For each radiomic feature, the following five sequences were analyzed, respectively: T2 fat suppression, T1 fat suppression, hepatic arterial phase, portal venous phase, and hepatocellular phase (20-minute delay after administration of contrast agent). Within each sequence, 1316 radiomic features were extracted, including 14 shape features, 18 first-order intensity statistics features, 75 texture features, 465 logarithmic features, and 744 wavelet features. In total, 19,740 radiomic features were extracted for each patient. Extraction was conducted by one single reader (reader 1). However, to evaluate the stability of the extracted radiomic features, another reader (reader 2) randomly selected 30 patients and repeated the extraction. The intraclass correlation coefficient (ICC) was calculated to test the stability of the extracted radiomic features.

All extracted radiomic features were standardized with Z-scores before further analysis, and only radiomic features with excellent stability (ICC > 0.75) were included.

Model Construction and Evaluation

We built the following five models first: clinical model, radiomics model-T, radiomics model-PT (5 mm), radiomics model-PT (5 + 5 mm), and radiomics model-T+PT. A combined model incorporating clinical characteristics and radiomic features was then developed.

Before model construction, we performed hyperparameter optimization, including the number of trees, the number of random variables for splitting, the minimum number of events in terminal nodes, and the maximum depth to which a tree should be grown, to increase model generalizability. Hyperparameters were tuned using five-fold cross-validation based on the training cohort. Hyperparameter optimization was performed using a Bayesian model-based optimization algorithm search through a predefined hyperparameter space.

We used a random survival forest with the log-rank splitting method for model learning based on the above-mentioned selected features and hyperparameters. RSF ranks candidate features based on the variable importance (VIMP), which is calculated by comparing out-of-bag prediction performance for the permuted feature to the original feature. A forest of 1000 trees was grown using log-rank splitting, and VIMP was recorded for each feature. The analysis was repeated 100 times independently, and VIMP was averaged over the runs. We then aggregated the top 20 engineered features from the RSF algorithm.

The performance of the models built based on the training cohort was further examined in the validation cohort using the concordance index (C-index). The average performance in the validation cohort was presented as the final model performance, which was generated by bootstrap resampling (n = 1000).

Statistical Analyses

Baseline characteristics between the training cohort and the validation cohort were compared using the t-test, the Chi-squared test, or the Mann–Whitney U-test, as appropriate. Receiver operating characteristic (ROC) curves and time-dependent ROC curve analyses were further performed to evaluate the discriminative abilities of the models. A decision curves analysis (DCA) was conducted to assess the utility of models for decision-making. The median risk value derived from the combined model was used as a cutoff for risk categorization, and the distribution of RFS was compared between groups by Kaplan–Meier survival curve analysis.

Statistical analyses were performed using R software (v. 3.6.3; http://www.Rproject.org). The following packages were used in the analyses: package “survival” for Cox regression; package “RandomForestSRC” for radiomic feature selection; packages “mlr,” “mlrMBO,” “randomForestSRC,” and “boot” for hyperparameter optimization, model building, and validation; package “timeROC” for time-dependent ROC curve analysis; package “dca” for DCA; and package “survminer” for Kaplan–Meier survival curve analysis.

Results

Baseline Characteristics of the Study Population

The study population had a median age of 54 years (46–62 years), and 86.4% of patients (114/132) were male. A total of 92.4% of patients (122/132) were HBsAg positive. The median tumor size was 20.00 mm (15.75–24.25 mm), and 19.7% of patients (26/132) had more than two tumor sites. One hundred patients (75.8%) underwent RFA, while 32 patients (24.2%) underwent MWA. As presented in Table 1, there were no significant differences observed in the baseline characteristics between the two cohorts.

Table 1 Baseline Characteristics of the Study Population

Follow-Up and Survival of the Study Population

The median follow-up times for the training and validation cohorts were 818 days (733–903 days) and 1029 days (794–1264 days), respectively. In the training cohort, 60.9% of patients (56/92) had confirmed tumor recurrence during follow-up at time points ranging from 33–1911 days. Median RFS was 592 days, and average RFS was 878 days. Overall 1-, 2-, and 3-year cumulative RFS rates were 69.5%, 42.5%, and 28.0%, respectively. Survival of the validation cohort was similar to the training cohort (Table 1). Median RFS was 654 days, and average RFS was 1152 days. Overall 1-, 2-, and 3-year cumulative RFS rates were 75.0%, 46.6%, and 46.6%, respectively. Survival of the validation cohort was similar to the training cohort (Table 1).

Clinical Characteristics Associated with RFS

Results of univariable and multivariable Cox regression analyses are presented in Supplementary Table 3. The following clinical characteristics were significantly associated with RFS: ALB (hazard ratio [HR] 0.66, 95% confidence interval [CI] 0.39–1.12, compared with ALB ≤40 g/L), GGT (HR 1.65, 95% CI 1.02–2.66, compared with GGT ≤60 U/L), tumor size (HR 3.25, 95% CI 1.79–5.90, compared with tumor size ≤30 mm), and AFP (HR 1.64, 95% CI 1.04–2.60, compared with AFP ≤25 ng/mL).

Radiomic Features Associated with RFS

The top 20 best ranking radiomic features were identified from 19,740 radiomic features using a random survival forest algorithm. Among them were 8 tumoral radiomic features, 6 peritumoral radiomic features (5 mm), and 6 peritumoral radiomic features (5 + 5 mm). Detailed information about radiomic features is presented in Supplementary Table 4.

Performance of Different Models

Based on the radiomic features selected from specific segmentation methods and combined with clinical characteristics, eleven models were developed, including radiomics model-T, radiomics model-PT (5 mm), radiomics model-PT (5 + 5 mm), radiomics model-T+PT-T1, radiomics model-T+PT-T2, radiomics model-T+PT-A, radiomics model-T+PT-V, radiomics model-T+PT-HBP, radiomics model-T+PT, the clinical model, and the combined model. The C-index for each model in the training and validation cohorts is presented in Table 2. Based on the C-index in the validation cohort, radiomics model-T+PT had better predictive performance (C-index 0.698, 95% CI 0.640–0.755) compared with the clinical model (C-index 0.614, 95% CI 0.499–0.695), while the combined model had the best predictive performance (C-index 0.706, 95% CI 0.638–0.763). Similar results were also observed in the ROC analysis and the time-dependent ROC analysis (Figures 2 and 3). In the validation cohort, the AUC of the clinical model was 0.672, with a sensitivity of 0.5 and a specificity of 0.767, the AUC of the radiomics signature was 0.817, with a sensitivity of 0.9 and a specificity of 0.7, the AUC of the combined model was 0.833, with a sensitivity of 0.9 and a specificity of 0.7. The DCA analysis (Figure 4) suggested that the combined model had a better net benefit compared with the other two models. When stratified based on the median of the combined model, distinct RFS distributions were observed with the Kaplan–Meier analysis (Figure 5), which indicated that the low-risk group had a higher RFS probability compared with the high-risk group (Log rank test: training cohort, p < 0.0001; validation cohort, p = 0.007). In the training cohort, the median time to recurrence was 322 days (IQR=260 days) of the high-risk group and 1857 days (IQR=1574 days) of the low-risk group. In the validation cohort, the median time to recurrence was 404 days (IQR=420 days) of the high-risk group and more than 2037 days of the low-risk group.

Table 2 Concordance Index of Different Models

Figure 2 Receiver operating characteristic curves of the models for predicting recurrence-free survival. (A) Training cohort; (B) validation cohort.

Figure 3 Time-dependent receiver operating characteristic curves of the models for predicting recurrence-free survival. (A) Training cohort; (B) validation cohort.

Figure 4 Decision curves of the models for predicting recurrence-free survival. (A) Training cohort; (B) validation cohort.

Figure 5 Kaplan–Meier survival analysis of recurrence-free survival according to risk strata defined by the combined model. (A) Training cohort; (B) validation cohort.

Subgroup Analysis

The study population was further categorized by AFP concentration (>25 ng/mL or ≤25 ng/mL), tumor size (>30 mm or ≤30 mm), GGT concentration (>60 U/L or ≤60 U/L), and ALB concentration (>40 g/L or ≤40 g/L). As presented in Supplementary Figures 1–4, within each subgroup, the median derived from the combined model showed consistent performance for risk stratification (p<0.001 for all).

Discussion

Given the lack of a practical prognosis predictive model specifically for HCC patients who undergo curative ablation, our study extracted radiomic features from MRI before curative ablation in patients with HCC and developed MRI-based radiomics models using a random survival forest algorithm to predict recurrence after curative ablation. Several models were built and compared, and the combined model, which incorporated 20 radiomic features (including tumoral, peritumoral [5 mm], and peritumoral [5 + 5 mm] radiomic features) and 4 clinical characteristics, had the best predictive performance. Based on this combined model, patients with HCC who underwent curative ablation were categorized as high-risk and low-risk for recurrence and showed distinct RFS distributions after curative ablation. This approach helped to identify patients with HCC who underwent curative ablation but were at a high risk of recurrence. Thus, this approach may have the potential to improve prognosis by more frequent monitoring²⁵ or by combined use with other therapies,^26,27 although this is yet to be proven.

Prognosis prediction models for HCC were traditionally developed using clinical characteristics and pathological characteristics at a macro-level, such as binary variables indicating microscopic vascular invasion.²⁸ However, considering that HCC has heterogeneous genotypes and phenotypes,^29,30 it is difficult to develop an ideal predictive model that could be applied to all HCC patients. Incorporating predictors strongly associated with events is one of the important methods to improve the performance of a predictive model. Compared with models developed using novel predictors, such as genes, RNA, and serum biomarkers,^19,20,31 our study used predictors extracted from MRI enhanced with gadoxetic acid disodium. Tumor segmentation was manually performed by an experienced radiologist. This increased the feasibility of the model since MRI data were more accessible compared with the above-mentioned novel predictors in patients with HCC. To some extent, radiomic features extracted from medical imaging can be seen as tumor micro-level features, which may help to distinguish the heterogeneous features of HCC. In addition, compared with pathological information obtained via biopsies, radiomic features could provide information on the whole tumor in a non-invasive way. Numerous studies support the association of radiomic features with tumor features and HCC prognosis,^32,33 making radiomic features potential predictors of prognosis. As expected, we noted that the model based on radiomic features only (radiomics model-T+PT) showed a nearly fair predictive performance (C-index 0.698, 95% CI 0.640–0.755), and the performance of the clinical model improved after incorporating the identified radiomic features (C-index increased from 0.614 [95% CI 0.499–0.695] to 0.706 [95% CI 0.638–0.763]). However, it is worth noting that models based on a single radiomic feature (tumoral or peritumoral) show poor performance, suggesting that a combination of tumoral and peritumoral radiomic features is necessary. As far as we know, there are no MRI-based radiomics models that have been specifically developed to predict recurrence after MWA for patients with HCC; thus, we were unable to compare the results of our study with those of other studies. Yuan et al³⁴ developed a contrast-enhanced computed tomography-based radiomics model to predict early recurrence after curative ablation for patients with HCC, and similar results were observed. Specifically, incorporating radiomic predictors increased the predictive performance of the clinical predictor-based model. Compared with the results of our study, the combined model in Yuan’s study had better performance (C-index  0.755 [95% CI 0.651–0.860] vs 0.706 [95% CI 0.638–0.763], respectively); however, given the different study population, a direct comparison may not be suitable. Moreover, our research was based on MR images rather than computed tomography images, which are recommended in clinical practice.

Our study had some strengths that should be highlighted. First, given that reproductivity is a major concern in the field of radiomics,²³ we only selected radiomic features with excellent stability (ICC > 0.75) for model development. Second, MR images used in the study were enhanced with gadoxetic acid disodium, which allow to identify lesion hypointensity that is one of the key steps during hepatocarcinogenesis,³⁵ and therefore may help radiologists to segment lesions more precisely by better delineating lesion boundaries. Third, we extracted different radiomic features, including tumoral, peritumoral (5 mm), and peritumoral (5 + 5 mm), and in total, 19,740 radiomics features were identified for each patient. This might provide more information and therefore be more representative of the heterogeneous features of HCC. Although this was a retrospective single-center study without external cohort validation, four different MR scanners were used in the study, which indicate that our radiomics model may be widely performed in multiple centers. Finally, a random forest algorithm was employed for radiomic feature selection and model construction, which enabled a high degree of predictive accuracy.³⁶

There were several limitations in our study that should be highlighted. First, the study used a retrospective study design, which is prone to information bias. Second, images used in this study were obtained using several MRI systems, which made it difficult to standardize the acquisition parameters. This may have resulted in bias to the post-imaging workflow. In addition, the model development was based on limited patients from a single center, and only internal validation was performed. External validity should be considered in further studies to examine the generalizability of our model.

Conclusion

The radiomics model, which combined radiomic features extracted from gadoxetic acid disodium-enhanced MR images with clinical characteristics, could be useful to predict HCC recurrence after curative ablation.