Abstract

Objectives:

To develop and validate radiomics and deep learning models based on the hepatobiliary phase (HBP) of gadoxetic acid-enhanced MRI (EOB-MRI) for the noninvasive prediction of treatment response and prognosis following transarterial chemoembolization (TACE) in hepatocellular carcinoma (HCC).

Materials and methods:

From April 2018 to September 2024, 160 patients with Barcelona Clinic Liver Cancer (BCLC) stage A or B HCC (>3 cm) were retrospectively enrolled and randomly divided into training (n = 112) and test (n = 48) sets. An independent cohort of 38 HCC patients was used for external validation. Twenty-six radiomic features were extracted using LASSO to construct a machine learning model using eXtreme gradient boosting (XGBoost). A deep convolutional neural network (DCNN) based on the ResNet50 architecture was used to develop a deep learning model. A clinical model was built via multivariate logistic regression. Model performance was evaluated by the area under the curve (AUC), calibration curves, and decision curve analysis (DCA). Kaplan–Meier analysis of combined deep learning and radiomics(DLR) scores was used to estimate overall survival in the follow-up cohort (n = 117).

Results:

BCLC stage (P = 0.035), tumor number (P = 0.015), and tumor size (P = 0.013) were independent clinical predictors. In the training set, the AUCs (95% CI) for the clinical, radiomics (XGBoost), and deep learning (DCNN) models were 0.79, 0.84, and 0.96, respectively. In the test set, the AUCs were 0.70, 0.80, and 0.92, respectively. In the external validation set, the AUCs were 0.77, 0.80, and 0.86, respectively. The DCNN model showed superior calibration and the highest net clinical benefit in DCA. In addition, multivariable Cox regression revealed that DLR model output was an independent risk factor for the overall survival (hazard ratio: 15.9, 95% CI: 4.49-56.33; p < 0.001).

Conclusion:

HBP-based AI models effectively predicted TACE response and prognosis in HCC patients, with the DCNN model showing the best performance. The integrated DLR model demonstrated high predictive accuracy and may serve as a reliable tool for individualized treatment planning in precision oncology.

Introduction

Primary liver cancer poses a significant global health burden and is the sixth most commonly diagnosed cancer and the third leading cause of cancer-related death (1). The 5-year survival rate ranges from 5% to 30% (2). Approximately 80% of primary liver cancers are hepatocellular carcinoma (HCC) (3). Surgical resection and liver transplantation are the first-line recommended treatments for curing HCC, with 5-year survival rates above 60%. However, HCC often develops without symptoms and progresses rapidly, and more than half of patients are diagnosed at intermediate or advanced stages, making them ineligible for curative treatments (4). According to the Barcelona Clinic Liver Cancer (BCLC) guidelines, transarterial chemoembolization (TACE) is recommended as the first-line treatment for unresectable HCC, either alone or in combination with other therapies such as ablation, targeted therapy or immunotherapy. TACE achieves an objective response rate of approximately 52.5% and can effectively control tumor progression and prolong survival (5, 6). Additionally, it is commonly used as a bridging or downstaging strategy to reduce the tumor burden and help patients become eligible for surgical resection or liver transplantation. Nevertheless, the treatment response to TACE remains highly variable due to tumor heterogeneity, and nearly half of patients fail to achieve objective remission (5). Early identification of nonresponders is essential, as timely initiation of targeted therapy or immunotherapy has been shown to improve clinical outcomes (7–10).

Current predictive models integrate clinical variables (BCLC stage and subclassification, tumor burden, elevated serum alpha-fetoprotein (AFP), etc.) and morphological features (central location, bilobar distribution, segment I/IV involvement) (11–13). With advances in artificial intelligence (AI), radiomics and deep learning have enabled the extraction and analysis of high-dimensional data from medical images, improving disease diagnosis, prediction of treatment response, and assessment of prognosis. Radiomics involves the quantification of large sets of imaging features, whereas deep learning employs neural networks to capture complex spatial patterns within image data. Recent studies have reported pooled AUCs of 0.82–0.89 for machine learning models in predicting TACE response, highlighting their promising clinical potential (14–25). However, most previous studies have focused on extracellular MR contrast agents, with only a limited number of studies exploring gadoxetic acid (Gd-EOB-DTPA)-enhanced magnetic resonance imaging (EOB-MRI) (11, 26).

The hepatobiliary phase (HBP) of EOB-MRI provides functional insights into tumor pathology and the microenvironment, both of which are closely associated with treatment response. Previous research has demonstrated its superiority in predicting tumor differentiation, microvascular invasion (MVI) and invasiveness, all of which are critical for predicting treatment outcome and prognosis (27–30). Despite its clinical value, the application of AI to EOB-MRI remains limited, with studies focused primarily on small HCCs and conducted in small cohorts (26). The potential of EOB-MRI for predicting TACE response in larger tumors or those at BCLC B stage has yet to be thoroughly explored.

The aim of this study is to develop and validate radiomics and deep learning models based on HBP sequences from EOB-MRI for the noninvasive prediction of TACE response in patients with HCC. In addition, the prognostic performance of the deep learning and radiomics (DLR) models were further evaluated. This approach may facilitate the identification of patients most likely to benefit from TACE and support personalized treatment planning.

Materials and methodsPatients and MRI protocol

This study was approved by the hospital’s Ethics Committee (Approval No. [2022]788). The data were obtained from a prospective project, and this work represents a retrospective analysis of those data. Therefore, the requirement for additional informed consent was waived. A total of 346 consecutive patients who were diagnosed with BCLC stage A or B HCC; confirmed by pathology, clinical-radiologic criteria, or digital subtraction angiography; and treated with TACE between April 2018 and September 2024 were initially enrolled. The inclusion criteria were as follows: (1) confirmed diagnosis of HCC; (2) underwent EOB-MRI within 2 weeks prior to TACE; (3) complete clinical data, including demographic information, hepatitis status, AFP levels, and liver function tests; and (4) follow-up of at least 3 months after TACE. The exclusion criteria were as follows: (1) underwent other treatments, such as radiofrequency ablation, surgical resection, radiotherapy, or systemic chemotherapy; (2) poor image quality; (3) the presence of other malignancies; (4) a diagnosis of advanced HCC (BCLC stage C); and (5) a tumor diameter less than 3 cm. A total of 160 HCC patients were ultimately included. The patient recruitment process is illustrated in Figure 1. External validation was performed on 38 HCC patients treated with TACE in another branch of the institution (The second affiliated hospital of Chongqing Medical University, Yuzhong branch). The inclusion and exclusion criteria for these patients were consistent with those of the previous cohort.

Flowchart of patient enrolment and study design. HCC, hepatocellular carcinoma; TACE, transcatheter arterial chemoembolization; EOB-MRI, gadoxetate disodium-enhanced magnetic resonance imaging; DSA, digital subtraction angiography; RFA, radiofrequency ablation; BCLC, Barcelona Clinic Liver Cancer; OR, objective response; nOR, non-objective response.

MRI examinations were performed on a 3.0T (Magnetom Prisma, Siemens Healthineers, Erlangen, Germany) or 1.5T (Signa HDxt, Milwaukee, WI, USA) scanner with an 18-channel body phased array coil. All patients underwent upper abdominal EOB-MRI. Gd-EOB-DTPA (Primovist, Bayer Schering Pharma, Berlin, Germany) was administered intravenously at a dose of 0.025 mmol/kg at a rate of 1.0 mL/s, followed by a 20 mL saline flush at the same rate. The detailed scanning parameters are described in Supplementary Table 1.

TACE treatment and reference standards for TACE response

All enrolled patients received TACE treatment, including conventional TACE (cTACE), drug-eluting bead TACE (DEB-TACE) and combined TACE. The choice between cTACE and DEB-TACE was determined on the basis of liver function and tumor characteristics. With the exception of embolic agents, the basic treatment processes of cTACE and DEB-TACE are similar. cTACE uses Lipiodol (Yantai, China), gelatin sponge particles, and polyvinyl alcohol as embolic agents. For DEB-TACE, CalliSpheres® microspheres (Merit Medical, USA) with a diameter of 100–300 μm or 300–500 μm were used as carriers and were loaded with 60–80 mg of epirubicin or pirarubicin. All procedures were performed by interventional radiologists with over 10 years of clinical experience who adhered to current practice guidelines. Following TACE, all patients were admitted for postoperative supportive care and routine management.

For patients with multiple lesions, we selected the largest lesion for both radiomic feature extraction and mRECIST evaluation. Tumor remission was assessed on the basis of the modified Response Evaluation Criteria in Solid Tumors (mRECIST) 3–4 months after TACE. The therapeutic response to TACE was classified into four grades: (1) complete response (CR), defined as complete disappearance of any intratumoral arterial enhancement in the targeted lesion; (2) partial response (PR), defined as a decrease of at least 30% in the sum of the diameters of viable target lesions (enhancement in the arterial phase); (3) progressive disease (PD), defined as an increase of at least 20% in the sum of the diameters of viable (enhancing) target lesions; and (4) stable disease (SD), defined as any case that did not qualify for PR or PD. According to mRECIST, patients with CR or PR were categorized as having an objective response (OR), whereas those with SD or PD were classified as having a nonobjective response (nOR) (Figure 2). Tumor response was independently assessed by two certified abdominal radiologists (J.Q. Ma and Y.Y. Liu) with 5 and 10 years of expertise in abdominal imaging on the basis of follow-up enhanced MRI or CT. Among the 160 enrolled patients, 112 were assigned to the OR group, and 48 were assigned to the nOR group. The patients were randomly divided into a training set (112 patients; OR = 78, nOR = 34) and a test set (48 patients; OR = 34, nOR = 14). In the external validation cohort, 32 patients were classified as OR and 6 as nOR.

Representative examples of post-TACE response in HCC patients. (a) A 57-year-old male patient presented with a 43 mm mass in the right hepatic lobe, showing slightly hyperintense on T2WI and DWI, marked arterial enhancement, and hypointensity on the HBP. After TACE, the lesion significantly reduced in size with no residual enhancement, corresponding to OR according to mRECIST criteria. (b) A 51-year-old male patient showed a 38 mm lobulated mass with irregular margins in the right hepatic lobe, hyperintense on T2WI and DWI, exhibiting marked arterial enhancement and hypointensity on HBP. After TACE, the lesion exhibited partial necrosis, but the reduction in tumor size was less than 30%, consistent with SD. TACE, transcatheter arterial chemoembolization; HCC, hepatocellular carcinoma; T2WI, T2-weighted imaging; DWI, diffusion-weighted imaging; AP, arterial phase; PVP, portal venous phase; HBP, hepatobiliary phase; CR, complete response; SD, stable disease; mRECIST, modified Response Evaluation Criteria in Solid Tumors.

Post-TACE follow-up was performed every 2–3 months with CECT/MRI, AFP, and liver function tests. OS (Overall survival) was defined as the time from the first TACE to death from study-related causes (as of September 1, 2024). Those alive or lost to follow-up by that date were censored.

Construction of the radiomics model

The model construction workflow is illustrated in Figure 3. HBP images from EOB-MRI were exported from the Picture Archiving and Communication System in Digital Imaging and Communications in Medicine format. The volume of interest of the tumor areas was manually delineated independently using ITK-SNAP version 3.8.0 (www.itk-snap.org) by two radiologists with 5 and 10 years of expertise in abdominal imaging. The intraclass correlation coefficient (ICC) was used to assess interobserver agreement, with an ICC value greater than 0.75 considered indicative of good reproducibility.

Model development and validation workflow. DCNN, deep convolutional neural network; ICC, intraclass correlation coefficient; HBP, hepatobiliary phase.

To enhance model stability and generalizability, all images were standardized and normalized to minimize interindividual variations in imaging parameters using min-max scaling. A harmonization step was used to correct for batch effect between scanners, ensuring consistency across different scanner models. High-throughput radiomic features were extracted from the preprocessed hepatobiliary phase images using the PyRadiomics module (version 3.0.1) in Python. A total of 1050 radiomics features were extracted for further analysis, including first-order statistical features, shape-based features, the gray level co-occurrence matrix (GLCM), the gray-level run-length matrix (GLRLM), the gray level zone matrix (GLSZM), the neighboring gray tone difference matrix (NGTDM), the gray level dependence matrix (GLDM) and wavelet features. These features were classified into three groups: geometric, intensity, and texture features. Texture features were extracted using four methods: GLCM, GLRLM, GLSZM, and NGTDM.

A total of 112 and 48 target lesions were included in the training cohort and test cohort, respectively. Feature selection was performed using the Mann–Whitney U test to identify significantly different radiomic features between response groups. Least absolute shrinkage and selection operator (LASSO) regression, combined with 10-fold cross-validation, was then used to select the features most relevant to the treatment response (Supplementary Figure 1). A final set of 26 radiomic features was retained. These features were linearly combined using their corresponding LASSO coefficients to calculate a radiomics score (Radscore) for each patient, which was implemented using the Python scikit-learn library (version 1.0.2). In addition, Shapley additive explanation (SHAP) value analysis was utilized to interpret the model and identify key features (Supplementary Figure 1).

Construction of the deep learning model

To improve the stability of the model, all the input images had their pixel intensity ranges normalized to [0, 1]. Data augmentation techniques, including random rotation, horizontal and vertical flipping, and contrast adjustment, were applied to increase image diversity and reduce the risk of overfitting.

A deep convolutional neural network (DCNN) based on the ResNet50 architecture was used to construct a classification model for the HBP of EOB-MRI. ResNet50 employs a residual learning framework, which facilitates deeper network training and improved feature extraction. The fully connected layer was stabilized using batch normalization and regularized with a dropout layer (dropout rate = 0.5) to prevent overfitting. A Softmax activation function was used in the output layer to generate a probability distribution for binary classification.

Cross-entropy loss was adopted as the loss function to quantify the discrepancy between the predicted and actual class labels. Model parameters were optimized using the Adam optimizer combined with a cosine annealing learning rate decay schedule to promote efficient convergence and stable training. The dataset was randomly split into training (n = 112) and test (n = 48) sets at a 7:3 ratio. During training, both loss and accuracy were monitored in real time, and the model’s generalization performance was continuously evaluated using the test set. The model that performed best on the test set was retained.

Nomogram construction

Clinical data preprocessing involved a stratified approach to handle missing values on the basis of variable type and missingness rate. Variables or samples with a missing data rate of 30% or higher were excluded to ensure data quality. For the remaining data, the distribution of continuous variables was assessed using the Shapiro–Wilk test. Normally distributed continuous variables were imputed using the mean, whereas nonnormally distributed variables were imputed using the median. Ordinal categorical variables were imputed with the median, and nominal categorical variables were imputed with the mode. Sensitivity analyses, including complete-case analysis and multiple imputation comparisons, were conducted to assess the robustness of the imputation.

After imputation, categorical variables were subsequently encoded using one-hot encoding. Continuous variables were standardized with Z score normalization if normally distributed or normalized with min–max scaling otherwise to reduce the influence of outliers. To address class imbalance, the synthetic minority oversampling technique (SMOTE) was applied to oversample the minority class and improve model generalizability.

Following data preprocessing, statistical analyses were performed. Clinical variables were categorized on the basis of standard reference ranges and initially analyzed using univariate analysis. Variables with statistical significance were then included, along with the Radscore, in a multivariate logistic regression model. A nomogram was constructed on the basis of the significant predictors identified through multivariate analysis. The final model was selected using a backward stepwise selection method, with the Akaike information criterion (AIC) used as the stopping rule. Stepwise logistic regression was used to calculate odds ratios (ORs) and 95% confidence intervals (CIs).

Evaluation of different models

The predictive performance of different models was evaluated in terms of classification accuracy, sensitivity, specificity, and the area under the receiver operating characteristic curve (AUC) in both the training and test cohorts. Calibration curves were plotted to compare the predicted probabilities against the observed outcomes. Decision curve analysis (DCA) was performed to estimate the net clinical benefit over a range of threshold probabilities and evaluate the clinical utility of each model.

Statistical analyses

All data processing and statistical analyses were conducted using SPSS (version 25.0) and Python (version 3.7.1). For continuous variables, the statistical tests included Student’s t test, Welch’s t test, and the Mann–Whitney U test, depending on the variance assumption. For categorical variables, the chi-square test or Fisher’s exact test was used. Statistical significance was defined as a two-sided p value < 0.05. The deep learning model was constructed using the PyTorch platform (version 1.13.1), whereas the radiomics model was built using an eXtreme gradient boosting (XGBoost) classifier. Tenfold cross-validation and ROC curve analysis were applied to validate the predictive performance of both models. Unless otherwise specified, all hyperparameters and configurations were set to the default values provided by the respective software platforms. In addition, a combined DLR model was constructed using independent factors selected from deep learning and radiomics features. Kaplan–Meier curves were used to assess the association between DLR scores and OS in the follow-up cohort. The Youden index determined the optimal threshold, dividing patients into high- and low-risk groups. Potential prognostic variables for OS were identified through univariate and multivariate Cox regression analyses, with statistical significance set at p < 0.05.

ResultsBaseline characteristics of the patients

Table 1 shows the results of univariate analysis comparing baseline demographic, clinical, and MRI characteristics between the training cohort (n = 112) and test cohort (n = 48). No statistically significant differences were observed between the two groups.

VariablesTraining set (n=112)Validation set (n=48)p valueAge (years)58.2 [10.2]58.6 [11.8]1.000Sex0.818 Male97 (86.6)41 (85.4) Female15 (13.4)7 (14.6)Hypertension0.653 No55 (49.1)39 (81.2) Yes57 (50.9)9 (18.8)Diabetes0.411 No93 (83.0)43 (89.6) Yes19 (17.0)5 (10.4)Etiology HBV105 (93.8)46 (95.8) HCV3 (2.7)1 (2.1) HBV and HCV1 (0.9) Others3 (2.7)1 (2.1)Child–Pugh classification0.918 A83 (74.1)37 (77.1) B/C29 (25.9)11 (22.9)BCLC stage0.890 A60 (53.6)27 (56.2) B52 (46.4)21 (43.8)Treatment modality0.716 c-TACE48 (42.9)21 (43.8) DEB-TACE20 (17.8)8 (16.6) Combined TACE44 (39.3)19 (39.6)AFP (ng/ml)0.470 ≤40087 (77.7)34 (70.8) >40025 (22.3)14 (29.2)PIVAK-II (mAU/ml)0.944 ≤4049 (43.8)20 (41.7) >4063 (56.2)28 (58.3)CEA (ng/ml)1.000 ≤4.5109 (97.3)47 (97.9) >4.53 (2.7)1 (2.1)ALT (U/L)0.381 ≤5077 (68.8)37 (77.1) >5035 (31.2)11 (22.9)AST (U/L)0.469 ≤4060 (53.6)22 (45.8) >4052 (46.4)26 (54.2)GGT (U/L)0.730 ≤6054 (48.2)21 (43.8) >6058 (51.8)27 (56.2)Albumin (g/L)1.000 >4023 (47.9)52 (46.4) ≤4025 (52.1)60 (53.6)Total bilirubin (µmol/L)0.735 ≤17.194 (83.9)42 (87.5) >17.118 (16.1)6 (12.5)Prothrombin time (s)0.852 ≤1332 (66.7)78 (69.6) >1316 (33.3)34 (30.4)INR0.556 ≤1.3100 (89.3)45 (93.8) >1.312 (10.7)3 (6.2)Radiological cirrhosis0.186 Present34 (30.4)9 (18.8) Absent78 (69.6)39 (81.2)Tumor number0.825 <390 (80.4)40 (83.3) ≥322 (19.6)8 (16.7)Tumor size (mm)51.8 [25.6]56.0 [27.2]0.369Tumor margin0.147 Smooth margin60 (53.6)19 (39.6) Nonsmooth margin52 (46.4)29 (60.4)Position0.204 Center31 (27.7)7 (14.6) Peripheral71 (63.4)36 (75.0)Shape0.228 Round62 (55.4)24 (50.0) Irregular34 (30.4)19 (39.6)Peritumoral hyperenhancement0.162 Yes94 (83.9)35 (72.9) No18 (16.1)13 (27.1)HBP intensity0.734 Iso- to hyper7(6.3)1(2.1) hypo89 (79.5)40 (83.4) heterogenous16 (14.3)7 (14.6)LR-M0.823 No101 (90.2)42 (87.5) Yes11 (9.8)6 (12.5)

Clinical characteristics of the training and validation sets.

Data are presented as medians [interquartile ranges] or numbers (percentages). HBV, hepatitis B virus; HCV, hepatitis C virus; BCLC, Barcelona Clinic Liver Cancer; c-TACE, conventional transcatheter arterial chemoembolization; DEB-TACE, drug-eluting beads–transcatheter arterial chemoembolization; AFP, alpha-fetoprotein; PIVKA-II, protein induced by vitamin K absence or antagonist-II; CEA, carcinoembryonic antigen; ALT, alanine amino-transferase; AST, aspartate amino-transferase; GGT, γ-glutamyl transpeptidase; HBP, hepatobiliary phase; LR-M, LI-RADS M (definite or probable malignancy, not specific for hepatocellular carcinoma).

Bold values indicate variable headings and subgroup labels in the table.

Table 2 summarizes the univariate and multivariate analyses comparing baseline characteristics between the OR and nOR groups. In the univariate analysis, significant differences were found in BCLC stage (P = 0.023), γ-glutamyl transpeptidase (GGT, P = 0.038), tumor number (P = 0.001), tumor size (P = 0.001), and tumor location (P = 0.012). Multivariate analysis identified BCLC stage (P = 0.035), tumor number (P = 0.015), and tumor size (P = 0.013) as independent predictors of treatment response.

VariablesUnivariate analysisMultivariate analysisOR (n=112)nOR (n=48)p valueOR (95% CI)p value