Multiparametric MRI-Based Machine Learning Radiomics Prognostic Models for Multifocal Hepatocellular Carcinoma Beyond Milan Criteria: A Retrospective Study.

Introduction
Hepatocellular carcinoma (HCC) is the fifth most common cancer and third leading cause of cancer-related deaths,1 with a 5-year survival rate of 21%.2 Multifocal HCC (MHCC), present in 41–75% of patients at diagnosis,3 is biologically and anatomically heterogeneous4 and generally has a poorer prognosis than solitary HCC.5 According to the Barcelona Clinic Liver Cancer staging system, liver transplantation (LT), resection, or ablation is recommended for MHCC within the Milan criteria (2–3 nodules, ≤ 3 cm), while transarterial chemoembolization (TACE) is advised for cases beyond these limits (2–3 nodules, > 3 cm or ≥ 4 nodules).6 Recently, there has been growing interest in expanding resection indications for selected patients with MHCC beyond Milan criteria.7–9
While LT is associated with superior long-term survival, limited organ availability and strict criteria restrict its use. Consequently, surgical resection remains the preferred curative option for resectable HCC.10 Prognostic factors such as microvascular invasion (MVI), satellite lesions, and poor differentiation are associated with recurrence and survival,9,11 but are typically assessed postoperatively. This highlights the need for accurate preoperative models to guide surgical decision-making in patients with MHCC beyond the Milan criteria.
Multiparametric magnetic resonance imaging (mpMRI) is widely used for HCC diagnosis and staging. Several studies have linked imaging characteristics, including tumor margin, peritumoral hyperenhancement, peritumoral hepatobiliary phase hypointensity, and diffusion restriction to histological aggressiveness.12–16 Nonetheless, the subjectivity and interobserver variability of these features limit their reliability.17
Radiomics, a quantitative imaging technique, enables the extraction of high-dimensional features that reflect tumor heterogeneity and has shown promise in predicting clinical outcomes.18 Despite its potential, radiomics has rarely been applied to mpMRI for MHCC beyond the Milan criteria. This study addresses this gap by: (1) identifying radiomics-based subtypes of MHCC beyond the Milan criteria using unsupervised clustering, and (2) developing and validating mpMRI-based radiomics models to predict postoperative outcomes. The potential improvement in predictive performance by integrating radiomics features with clinical variables was also evaluated.

Materials and Methods
This retrospective, single-center study was approved by the institutional review board of the Zhongshan hospital (B2021-682R). The requirement to obtain informed consent was waived because of the retrospective nature of the study. All methods were implemented in accordance with the approved regulations and the Declaration of Helsinki. Patient confidentiality was strictly protected in compliance with applicable laws and institutional policies. All personal identifiers were removed or anonymized before analysis, and the data were stored on secure, access-restricted servers.

Patient Cohorts
Patients suspected of MHCC based on MRI reports and who underwent curative hepatectomy between January 2015 and January 2019 were retrospectively registered (Figure 1). A lesion with a satellite nodule (a smaller nodule within 2 cm of the main tumor) was considered one HCC lesion. Inclusion criteria were: (a) pathologically confirmed HCC; (b) available clinicopathological data; and (c) no vascular invasion or biliary tract thrombosis on preoperative MRI. Exclusion criteria were as follows: (a) history of antitumor treatment; (b) presence of other malignant tumors; (c) missing or poor-quality MRI data; (d) MRI-surgery interval >1 month; and (e) patients within the Milan criteria.

MRI Study Protocols
MR images were acquired using a 1.5 T MR scanner (MAGNETOM Aera, Siemens Healthcare) equipped with a phased-array body coil. The protocol included T2-weighted, diffusion-weighted (b = 0 and 500 sec/mm2), and T1-weighted dynamic contrast-enhanced MRI (DCE-MRI) sequences (Supplementary Table E1). DCE-MRI was performed after administering gadolinium diethylenetriamine pentaacetic acid (Gd-DTPA; Magnevist, Bayer HealthCare) at 0.1 mmol/kg of body weight and 2.0 mL/sec flow rate. Arterial (20–30 s), portal venous (70–90 s), and delayed (160–180 s) phase images were obtained after the contrast injection.

Clinical Characteristics and Radiologic Assessment
Patients information from electronic health records included: (a) demographic and clinical characteristics, namely, age, sex, hepatitis B virus (HBV) infection, HBV deoxyribonucleic acid (HBV-DNA) load, levels of alpha fetoprotein (AFP), total bilirubin (TBIL), direct bilirubin (DBIL), total protein (TP), glutamic-pyruvic transaminase (ALT), glutamic oxaloacetic transaminase (AST), alkaline phosphatase (ALP), g-glutamyl transferase (GGT), total bile acid (TBA), and platelet (PLT); (b) histopathologic features, namely, ki67, MVI, Glypican-3 (GPC3), and Edmondson-Steiner (ES) grade; (c) follow-up information, namely, overall survival (OS) and recurrence-free survival (RFS). OS was defined as the time from surgery to death or last follow-up. RFS was the time from surgery to the first recurrence, metastasis, or last follow-up. Data were censored on December 31, 2021, which was the predefined cutoff at the study design stage to ensure consistency and integrity of survival analyses.
Three readers (F. W., X. D. Z. and C. Y., with 9, 16 and 20 years of experience in abdominal imaging, respectively) were blinded and independently reviewed all images to evaluate radiological features. Discrepancies were resolved through a majority vote. The evaluated radiological features included tumor number, tumor diameter, liver cirrhosis, satellite nodule, hemorrhage in mass, fat in mass, arterial rim enhancement, radiological capsule, arterial peritumoral enhancement, mosaic architecture, and nonsmooth tumor margin. If any HCC lesion displayed these MRI features, it was considered present in the patient. The radiological total tumor diameter (TTD), and ratio of the largest to the smallest tumor diameter (RLSD = diameter of the largest tumor/diameter of the smallest tumor) were calculated. The tumor burden score (TBS) was defined as TBS2 = (largest tumor diameter)2 + (tumor number)2.

Image Preprocessing and Radiomics Feature Extraction
Image preprocessing included N4 bias field correction, voxel spacing resampling, and intensity normalization using SimpleITK (Python version 3.10).19 Before segmentation, arterial phase (AP), portal venous phase (VP), and delayed phase (DP) images were aligned with pre-contrast images through non-linear registration using Advanced Normalization Tools (ANTs, version 2.3.5). For MHCC patients, only the largest tumor was contoured. Tumors were manually segmented slice-by-slice using 3D Slicer software (5.2.2) on T2-weighted images, apparent diffusion coefficient (ADC) map (generated from the combination of b = 0 and 500 sec/mm2), and pre-contrast DCE-MRI images by two readers. A 5-mm width band was generated around the tumor boundaries, excluding areas beyond the liver parenchyma. The volumes of interest (VOIs) included tumor, peritumor and tumor-peritumor regions.
Radiomics features of the three VOIs were extracted using Python (version 3.10) with the “Pyradiomics” package from the original, Laplacian, wavelet, local binary pattern, and gradient-transformed mpMRI. The features included shape, histogram-based and textural features defined by gray-level co-occurrence matrix (GLCM), gray-level dependence matrix (GLDM), gray-level run length matrix (GLRLM), gray-level size-zone matrix (GLSZM) and neighboring gray tone difference matrix (NGTDM). Dynamic changes in radiomics features between DCE-MRI phases were quantified by calculating the rate of change in features between sequential phases (Supplementary S1). Features were normalized using Z-scores, and those with zero variance or a mean absolute deviation (MAD) less than 1 were excluded from further analysis.

Unsupervised Clustering
To ensure feature-level repeatability, intra- and inter-class correlation coefficients (ICCs) were calculated by extracting radiomics features twice from the same images under consistent preprocessing conditions. Features with ICC values > 0.75 were retained for further analysis. Spectral clustering analysis was performed on mpMRI-based radiomics features from three VOIs. Patient-by-patient similarity matrices were constructed using Pearson’s correlation of multi-view radiomics features. Similarity network fusion was independently performed for both training and validation sets. In the validation set, cluster labels were assigned using the nearest centroid method based on the training results. The unsupervised spectral clustering algorithm (K=1, t=10)20 was applied as a purely data-driven approach to identify imaging-based subtypes without priori clinical or pathological assumptions. Kaplan-Meier (KM) survival curves and Log rank tests assessed associations between subtypes and outcomes.
The heterogeneity index (HI) combined with Z-score-normalized features was positively and negatively related to tumor heterogeneity (Supplementary S2 and S3).

Machine Learning Based Radiomics Risk Score
Initial radiomics features were selected using univariate Cox proportional hazard regression analysis. The Machine Learning (ML) based radiomics risk score (RRS) for OS and RFS was derived using a supervised extreme gradient boosting (XGBoost)-LASSO Cox proportional hazard regression model in the training cohort.21 This model integrated XGBoost with L1-penalized Cox regression, optimizing a partial log-likelihood loss function across the radiomics feature space. Hyperparameters were evaluated using stratified 5-fold cross-validation on the internal validation set, selecting the penalizing factor that maximized the averaged partial log-likelihood. Patients were stratified into high- and low-risk groups based on the median RRS. KM analysis and Log rank tests were used to compare OS and RFS between the risk-defined groups in both the training and validation sets. The model implementation utilized the scikit-learn Python library with the XGBoost-LASSO framework.

Development of Prognostic Models
Three prognostic model types were developed: (a) PreOp (preoperative clinico-radiological features), (b) radiomics (including unsupervised subtypes, HI, and RRS), and (c) combined. Uni- and multivariable Cox proportional hazard regression analyzed clinico-radiological and radiomics characteristics. Variables association with OS or RFS (P < 0.1) in univariable analysis were included in the multivariable models. Prediction models were constructed using scikit-learn and statsmodels in Python. The potential of radiomics predictors to improve prognostication beyond conventional risk factors in MHCC beyond Milan criteria was evaluated.

Statistics Analysis
Statistics analysis was performed with medcalc (version 20.0.4) and Python (version 3.10) using the “snfpy”, “statsmodels”, “sklearn” and “scikit-survival” packages. Continuous variables were reported as mean ± SD or median (IQRs), based on the Shapiro–Wilk test distribution. Categorical variables were presented as numbers and percentages. Intra-group differences were assessed using the Chi-squared test, Student’s T test, or Mann–Whitney U-test. Uni- and multivariable Cox proportional hazard regression identified risk factors, reporting hazard ratios (HR), and 95% confidence intervals (CI). KM survival curve analysis was used to evaluate OS and RFS with log-rank validation for group differences. Model performance was evaluated using the concordance index (C-Index) and time-dependent ROC analysis, while calibration curves assessed alignment between nomogram predictions and observed outcomes. To statistically compare the prognostic performance, C-index values of the PreOp model and radiomics model were compared with those of the combined model using the Wilcoxon signed-rank test, separately in both the training and validation sets. A two-sided P value < 0.05 was considered statistically significant.

Results

Patient Demographic Characteristics
A total of 592 patients who underwent hepatectomy for MHCC (without vascular invasion and biliary tract thrombosis on preoperative MRI) were consecutively enrolled. Based on the inclusion and exclusion criteria, 436 patients were excluded: 154 for receiving preoperative anticancer treatment, 37 for other malignancies, 33 for missing or poor-quality MRI, 15 for undergoing MRI more than one month before surgery, and 197 for meeting the Milan criteria (Figure 1). Consequently, 156 patients were included for further analysis and randomly allocated to the training (n = 78) and validation (n = 78) groups using a random number generator in Python (version 3.10). There were no statistically significant differences in the baseline characteristics between these sets (Supplementary Table E2).

Identification of Unsupervised MHCC Subtypes
In the training set, spectral clustering analysis based on the Pearson correlation coefficients of each subject’s similarity network fusion (Figure 2a) identified two distinct clusters (Figure 2b, subtypes A and B). Patient characteristics of these subtypes are shown in Table 1. Subtype B was significantly associated with a higher tumor number (P = 0.002), the presence of GPC3 (P = 0.028), and higher Ki67 (P = 0.003). No significant differences were observed in other demographic, clinico-radiological, or pathological characteristics (P = 0.101 to P = 0.929). In the validation set, notable differences included DBIL > 6.8 mmol/L (P = 0.037), AST > 40 U/L (P = 0.030), arterial rim enhancement (P = 0.050) and the presence of MVI (P = 0.015).

Identification of Heterogeneity Index
The HI in MHCC was evaluated using heterogeneity-related first-order and textual radiomic features based on arithmetic formula and standardized feature definitions (Supplementary S1). Considering the three VOIs and mpMRI, 21615 features were initially analyzed. Subtype B demonstrated significantly higher HI than subtype A in both the training (P < 0.001) and validation (P = 0.048) sets, indicating greater imaging heterogeneity (Figure 2c).

Prognostic Value of Unsupervised MHCC Subtypes
Unsupervised subtypes showed significant prognostic disparities in the training set, with subtype B associated with poorer OS (HR, 2.496 [95% CI: 1.129, 5.518]; P = 0.024; Figure 3a) and RFS (HR, 2.345 [95% CI: 1.341, 4.102]; P = 0.003; Figure 3b) compared with subtype A. This was validated in the validation set, in which subtype B also showed worse OS (HR, 2.420 [95% CI: 1.114, 5.258]; P = 0.026; Figure 3c) and RFS (HR, 1.744 [95% CI: 1.023, 2.973]; P = 0.041; Figure 3d). These findings underscore the clinically relevant prognostic information encoded within the unsupervised subtypes.

Radiomics Risk Score Predicts Prognosis
Based on radiomics-defined prognostic subtypes, the clinical utility of a ML-based RRS was examined for outcome prediction in MHCC beyond the Milan criteria. From mpMRI, 21615 features were extracted from the three VOIs. LASSO-based feature selection method identified 17 features for OS (Supplementary Table E3) and 13 features for RFS (Supplementary Table E4). Patients were stratified into low- and high-risk groups based on the mean value of RRS. Distinct OS (HR, 4.290 [95% CI: 1.927. 9.552]; P < 0.001; Figure 4a) and RFS (HR, 2.992 [95% CI: 1.709, 5.238]; P < 0.001; Figure 4c) were observed within the training sets. Similarly, high-risk patients in the validation set demonstrated poorer OS (HR, 5.273 [95% CI: 2.438, 11.403]; P < 0.001; Figure 4b) and RFS (HR, 3.972 [95% CI: 2.198. 7,175]; P < 0.001; Figure 4d) using the same stratification. Based on Kaplan–Meier survival analysis, the actual 1-, 3-, and 5-year overall survival (OS) rates were 82.4%, 58.3%, and 42.1% in the training set, and 85.7%, 62.1%, and 45.0% in the validation set, respectively.

Variables Associated with Prognosis
For MHCC patients beyond Milan criteria, univariable analysis indicated that OS was significantly associated with the presence of satellite nodule, fat in mass, unsupervised subtype, and RRS. Tumor number, tumor diameter > 5 cm, TTD, RLSD, TBS, presence of unsmooth tumor margin, unsupervised subtype, HI and RRS were significant predictors of worse RFS (Table 2). Multivariable Cox regression analysis revealed that OS could be independently predicted by the unsupervised subtype (HR, 1.697 [95% CI:1.019, 2.825]; P = 0.042) and RRS (HR, 44.278 [95% CI: 3.094, 633.667]; P = 0.005), along with the presence of satellite nodule (HR, 1.909 [95% CI: 1.047, 3.480]; P = 0.035) (Table 3). For RFS, the RRS (HR, 1.248 [95% CI: 1.152–1.351]; P < 0.001) and unsupervised subtype (HR, 1.493 [95% CI: 1.045–2.132]; P = 0.028) remained independent prognostic factors after adjusting for confounders (Table 4). These results indicate that the unsupervised radiomics-based subtype, though derived in a data-driven exploratory step, contributes meaningfully to prognostication and was retained in the final combined models.

Combined Model Improve Performance
Integration of the proposed radiomics predictors enhanced the PreOp model for OS prediction, as shown by an increased C-index (training set: from 0.616 [95% CI: 0.515, 0.731] to 0.712 [95% CI: 0.643, 0.788]; validation set: from 0.522 [95% CI: 0.423, 0.630] to 0.710 [95% CI: 0.634, 0.768]; Figure 5a). Similarly, RFS prediction models also demonstrated an improved C-index with the combined model compared to the PreOp model in Figure 5b (training set: from 0.653 [95% CI: 0.572, 0.717] to 0.735 [95% CI: 0.687, 0.771]; validation set: from 0.574 [95% CI: 0.481, 0.643] to 0.698 [95% CI: 0.629, 0.758]), and the radiomics model alone (training set: 0.714 [95% CI: 0.650, 0.773]; validation set: 0.635 [95% CI: 0.558, 0.700]), confirming the added prognostic performance of radiomics.

In the validation set, the combined model predicted 1-, 3-, and 5-year OS with corresponding area under the receiver operating characteristic curves (AUCs) of 0.70 (95% CI: 0.52, 0.88), 0.75 (95% CI: 0.66, 0.83), and 0.75 (95% CI: 0.66,0.83), respectively (Figure 6a). It also accurately predicted RFS at 1 year (AUC, 0.80 [95% CI: 0.70, 0.89]), 3 years (AUC, 0.71 [95% CI: 0.58, 0.83]), and 5 years (AUC, 0.75 [95% CI: 0.62, 0.87]) (Figure 6a) in the validation set. The combined model variables are illustrated in nomogram plots for OS and RFS, providing individualized probability estimates and highlighting the significance of each characteristic (Figure 6b and c).

Discussion
The variability in prognosis of intermediate-stage HCC, particularly those beyond the Milan criteria, highlights the need for accurate survival risk stratification.6,11,22 The implementation of survival risk stratification in HCC can be enhanced through a radiomics approach using preoperative prognostic models.23,24 However, no studies have identified personalized prognostic radiomics biomarkers for predicting outcomes in MHCC patients beyond Milan criteria after hepatectomy. Our study found a statistically significant association between unsupervised subtypes and prognosis in these patients. Two quantitative combined models were developed to predict OS and RFS, demonstrating enhanced performance when incorporating radiomics and preoperative risk factors, as shown by improved C-index values. These findings demonstrate the potential of mpMRI-based radiomics for noninvasive outcome prediction in MHCC post-hepatectomy.
Tumor burden is the primary predictor for stratifying HCC patients beyond Milan criteria.12,15 Tumor size and number remain key indicators of aggressiveness, highlighting the need for more precise markers. Radiomics provides detailed tumor lesion information, reflecting intratumoral heterogeneity and aiding precision management.16,25 Although the clustering was performed in an unsupervised, data-driven approach, the unsupervised subtypes showed distinct clinico-pathological characteristics, including differences in MVI status, Ki67, and GPC3, indicating possible biological underpinnings. Notably, these subtypes were not only statistically associated with survival outcomes but also emerged as independent prognostic factors in multivariable analysis, and were integrated into the final predictive models. Previous studies have shown that MRI-derived radiomics can predict the response to lenvatinib monotherapy in unresectable HCC.26 Radiomic features provide a radiological representation of tumor phenotypes. Integrating radiomics with clinical and pathological data could improve treatment response and prognosis prediction in precision oncology.
Clinical applicability was demonstrated by accurately predicting OS and RFS with a new XGBoost-lasso model. Integrating radiomics predictors (including the unsupervised subtype, HI, and RRS) into preoperative combined models significantly improved the predictive power, as measured by the C-index, compared to the conventional preoperative or radiomics-only models. Adding radiomics factors to the models showed their crucial role in enhancing prognostic performance. This finding suggests that incorporating mpMRI-derived radiomics improves risk assessment and patient stratification in MHCC patients beyond Milan criteria.
Radiomics provides detailed information on whole-tumor heterogeneity. Spatial tissue heterogeneity, captured by texture features, identifies image patterns based on voxel intensities and spatial arrangement. In this study, most RRS features, such as “GLDM”, “GLCM” and “GLSZM”, were texture features or wavelet transformations features related to both OS and RFS. Our findings align with previous studies showing that texture features outperformed first-order features for prognosis,27,28 supporting the inclusion of texture features is vital for enhancing prognostic performance.
Sequential radiomics captures relative net change in image features, including arterial phase hyperenhancement and washout, for HCC diagnosis. In this work, all RRS features of OS were derived from sequential radiomics, indicating that changes in DCE-MRI may have prognostic value. Similarly, 10 of 13 RFS features were sequential. Additionally, the VP_from_AP feature set, reflecting changes between the AP and VP, correlates with the “washout” pattern, which is a key criteria in HCC guidelines.29,30 These sequential features provide a dynamic view of tumor characteristics, enhancing the predictive power of radiomics models.
Numerous studies have emphasized the importance of analyzing the peritumoral region, as the tumor microenvironment might contain valuable prognostic information. Peritumoral features have been reported to predict MVI status, prognosis, and treatment response in HCC.23,31,32 In this study, nearly half of the RRS features were derived from peritumor or peritumor-tumor areas. Our combined model, incorporating tumor, peritumor, and tumor-peritumoral features, proved effective for HCC management. These findings confirm that peritumoral features add prognostic value in HCC.23,33 Integrating multiple VOI data emphasizes the role of intratumoral heterogeneity and surrounding tissue in patient stratification. Thus, combing multidimensional and multiregional radiomics features enhances prognostic prediction.
Several preoperative radiologic features, such as the presence of satellite nodules, fat in mass, unsmooth margin, and TTD, were identified as important prognostic indicators. These findings are consistent with previous studies.13,16,34–36 Although preoperative AFP has been widely evaluated,37,38 our findings and others39–41 suggest it is not robust, and thus it was excluded from our final model. Notably, our current work builds on our earlier study,42 which focused on predicting MVI with conventional imaging. Here, advanced radiomics and machine learning extend prognostic evaluation to overall and recurrence-free survival in a more advanced HCC.
The combined models offer a noninvasive tool to aid preoperative risk stratification in MHCC patients, potentially supporting surgical decision-making and personalized postoperative surveillance. To facilitate clinical adoption, future work should focus on prospective validation and development of user-friendly decision support tools that seamlessly incorporate radiomics findings into routine practice.
This study has several limitations. First, it is retrospective single-center study with a limited sample size and no external validation, which may limit generalizability. Future multicenter prospective studies are needed to confirm model robustness across institutions and imaging protocols. Second, MRI scans were acquired with a single 1.5-T scanner, which may restrict applicability to other imaging systems and field strengths. Third, the study cohort exhibited a gender imbalance, with significantly more men patients, reflecting the higher incidence of HCC in East Asia but potentially introducing selection bias. Fourth, manual tumor segmentation is time-consuming and may hinder clinical scalability. Although ICC filtering was applied to ensure test–retest stability of radiomics features, variability across different scanners, protocols, and institutions remains a known limitation and requires further multicenter validation. Future advances in automated segmentation are expected to substantially reduce the workload. Finally, the clinical implementation of the model faces several challenges, including standardization of imaging protocols, reproducibility of radiomics features across centers, and integration into clinical decision-making pathways.
In conclusion, this study demonstrates radiomics as a noninvasive method for stratifying MHCC patients beyond Milan criteria after hepatectomy by decoding clinical and prognostic differences. We present two radiomics-based models that enhance conventional preoperative models, offering a reliable tool for prognostication. It may provide new insights into clinical decision-making.