Novel Unsupervised Machine Learning Model Using Multi-Parametric Radiomics for Prognostic Stratification of Bifocal Hepatocellular Carcinoma

Introduction

Hepatocellular carcinoma (HCC) has a 5-year relative survival rate of approximately 18%, making it the third leading cause of cancer-related deaths globally. The grim prognosis is underscored by the nearly identical incidence and mortality rates.¹ Multifocal HCC (mHCC), particularly characterized by its biological heterogeneity, exhibits worse long-term outcomes compared to single-nodular HCC.² Notably, bifocal HCC (bHCC), accounting for more than 80% of mHCC cases, is often treated aggressively with radical resection in many Asian countries.^3–5 Despite curative-intent treatments, the postoperative recurrence rate among bHCC patients remains alarmingly high, which is largely attributed to pronounced tumor heterogeneity.⁶ Given this challenge, effective risk assessment and long-term prognostic stratification are crucial for guiding clinical decision-making and optimizing management of bHCC patients. While previous studies primarily focused on risk stratification for solitary HCC,^7–10 these approaches are inherently limited in bHCC as they fail to account for the comprehensive tumor burden derived from both lesions. Therefore, a comprehensive evaluation should encompass both tumor lesions in bHCC to capture the full extent of tumor characteristics.

Radiomics has demonstrated its capacity to characterize the complex imaging phenotypes of tumors and has shown significant potential in assessing both the tumor heterogeneity and prognosis of HCC.¹¹ Prior studies have predominantly utilized supervised machine learning techniques, including random forests, support vector machines, convolutional neural networks, transfer learning, and other deep learning algorithms.⁷ Recently, the integration of unsupervised machine learning algorithms with radiomics has opened new avenues for uncovering intricate patterns in medical imaging datasets.^12,13 Unlike traditional supervised radiomics analyses, unsupervised machine learning algorithms eliminate the need for manual labeling of clinical outcomes; this approach groups patients based on intrinsic phenotypic similarities rather than predefined outcomes.¹⁴

Similarity network fusion (SNF), an unsupervised machine learning algorithm, constructs networks of patients based on each available data type and then effectively merges these into a single network, representing the full spectrum of underlying data.¹⁵ This approach substantially outperforms single-data-type analyses and established integrative methods in identifying cancer heterogeneity, and it has proven effective in predicting survival outcomes.^16–20 The originality of this study lies in the application of SNF to address the unique challenges of bHCC. Unlike traditional integrative methods that often suffer from noise and loss of information during simple feature concatenation, SNF constructs patient similarity networks for two nodules and MRI multi-sequences, subsequently fusing them into a unified network. This approach is particularly suitable for bHCC as it effectively captures the synergistic radiomics signals from both lesions without discarding the imaging features of either nodule. By computing and fusing these patient similarity networks derived from multi-parametric MRI radiomics separately, SNF leverages the complementary nature of the data to better explore the heterogeneity in bHCC.

Therefore, the purpose of this study was to develop an SNF-based radiomics framework that jointly integrates two-lesion and multi-parametric MRI features to identify Radiomics Image Heterogeneity (RIH) clusters in bHCC. Based on subtype classification, we stratified bHCC patients into different prognostic risk subgroups preoperatively, and developed a clinical-radiological-radiomics model using multivariate Cox regression analysis to predict recurrence-free survival (RFS) and overall survival (OS).

Materials and Methods

The Institutional Review Board of Zhongshan Hospital approved this retrospective study (B2021-682R), which was conducted in accordance with the Declaration of Helsinki. The requirement to obtain written informed consent was waived due to the retrospective design of this study. We stated that patient data was strictly confidential. This study was conducted in accordance with the CLEAR, CLEAR-E3, and METRICS guidelines to ensure transparency, reproducibility, and methodological rigor in radiomics research.^21–23

Study Patients

From January 2015 to January 2019, patients with bHCC reported by MRI were retrospectively and continuously registered by retrieving our institutional picture archive and communication system (PACS). Inclusion criteria were as follows: Patients who received curative hepatic resection as the first-line treatment (R0 resection) and patients without major vascular invasion on preoperative MRI (n = 489). Exclusion criteria were as follows: (a) previous history of treatments before MRI including (chemotherapy, radiotherapy, trans-arterial chemoembolization, radiofrequency ablation, hepatectomy, liver transplantation and systemic therapy) (n = 148); (b) presence of extrahepatic metastasis (n = 51); (c) any lesion pathologically diagnosed as non-HCC after surgery, including combined hepatocellular carcinoma-cholangiocarcinoma (n = 16), intrahepatic cholangiocarcinoma (n = 4); (d) incomplete follow-up data (n = 63); (e) poor image quality due to breathing motion artifacts (n = 8), which assess by radiologist 1 (six years in liver imaging); (f) long time interval (longer than 1 month) between MRI examination and surgical resection (n = 17). Finally, 182 patients were included in our study and were randomly divided into training set and validation set in a ratio of 1:1 (Figure 1).

Figure 1 Patient enrollment flowchart. Abbreviations: HCC, hepatocellular carcinoma; bHCC, bifocal HCC.

Laboratory and Pathological Data

All variables were retrospectively collected through the electronic medical record system in our institute. Clinical information, including sex and age, was obtained during admission through medical history taking. Laboratory data were obtained from blood tests conducted within one week prior to surgery, including status of HBV or HCV infection, preoperative HBV deoxyribonucleic acid (HBV-DNA) load, levels of alpha fetoprotein (AFP), total bilirubin (TBIL), direct bilirubin (DBIL), total protein (TP), albumin (ALB), alanine aminotransferase (ALT), aspartate aminotransaminase (AST), alkaline phosphatase (ALP), g-glutamyl transferase (GGT), total bile acid (TBA), prothrombin time (PT), and platelet count (PLT). Following the surgical procedure, pathological evaluation of the resected specimens was performed by two experienced pathologists within 1–2 weeks. Histopathologic features such as microvascular invasion (MVI), cirrhosis of the noncancerous liver parenchyma, satellite nodule, and Edmondson-Steiner (E-S) grade were evaluated based on these findings.

Clinical Outcome

Follow-up clinical data were obtained by a review of medical records or telephone interviews. RFS was defined as the time from the date of surgery to the date of first recurrence, metastasis, or last follow-up. OS was defined as the time from the date of surgery to the date of death or last follow-up. The data were censored on December 31, 2021.

Power Analysis

Prior to the study, a power analysis was conducted to determine the minimum sample size required to detect significant differences in diagnostic performance, as measured by the area under the receiver operating characteristic (AUC-ROC) curve, with a significance level of 0.05 and a statistical power ranging from 0.80 to 0.95. The results indicated that to detect AUC values of 0.66, 0.70, 0.75, and 0.80, at least 86, 62, 42, and 30 patients, respectively, would be required. Given that 182 patients were included in the study and the observed AUC values ranged from 0.66 to 0.83, the sample size was deemed sufficient to ensure a high level of statistical confidence and validity. To further ensure the robustness of our prognostic modeling, we assessed the adequacy of the training dataset (n = 91) using the Events Per Variable (EPV) criterion. Given the high maturity of our cohort, with a maximum follow-up of 80 months and minimal censoring, the event density (127 for RFS and 52 for OS) provided ample statistical power. In the training cohort, the EPV ratios were maintained at a level that supports stable parameter estimation for both the RFS and OS multivariable Cox models, effectively minimizing the risk of overfitting.

MRI Protocol

All patients were scanned in a 1.5-T MR scanner (MAGNETOM Aera, Siemens Healthcare). Routine liver protocols consisted of transverse T2-weighted imaging (T2WI), T1-weighted imaging (PRE) and diffusion-weighted imaging (DWI). Dynamic imaging was performed with a T1-weighted fat-suppressed sequence. The arterial phase (AP) was acquired when the contrast agent (gadolinium diethylenetriamine pentaacetic acid, Gd-DTPA; Magnevist, Bayer HealthCare) reached the ascending aorta after intravenous administration at a dose of 0.1 mmol/kg at a rate of 2 mL/s. The portal venous phase (VP) and delayed phase sequences (DP) were acquired at 70–90 s and 160–180 s, respectively. Detailed information is available in Table S1. The use of a single scanner was intentional to minimize inter-scanner variability and technical noise, which are known to significantly affect the stability and reproducibility of radiomics features during the model development phase.

MRI Features Evaluation

Three radiologists with 3, 8, and 19 years of experience in abdominal imaging, respectively, blinded to clinical and pathological information independently evaluated the following qualitative imaging features: 1) Larger Diameter > 5 cm; 2) satellite nodule; 3) hemorrhage in mass; 4) fat in mass; 5) arterial rim enhancement; 6) radiological capsule; 7) arterial peritumoral enhancement; 8) mosaic architecture; 9) nodule-in-nodule architecture; 10) unsmooth margin; 11) atypical enhancement pattern (Supplementary material 1). If any HCC lesion showed the qualitative MRI feature described above, the feature was considered present in the patient. The radiological liver cirrhosis was also assessed. After independent interpretation of images, discrepancies between readers were resolved by a consensus discussion. 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 also calculated.

Tumor Segmentation and Radiomics Feature Extraction

A radiologist (X.J., with 3 years of abdominal imaging analysis experience) performed tumor segmentation by ITK-SNAP software, these segmentation results were checked by a senior radiologist (C.Y., with 19 years of abdominal imaging analysis experience). Volumes of interest (VOIs) were manually delineated on 6 sequences of pre-T1WI, AP, PVP, DP, T2WI-FS and DWI with b values of 500 s/mm². In addition, MR images of randomly selected 30 lesions were delineated again after 1 month by X.J. to assess the intra-observer reproducibility, and these 30 MRI images were also delineated by another radiologist (F.W., with 8 years of abdominal imaging analysis experience) independently to evaluate inter-observer reproducibility.

According to the radiomics checklist,^21–23 image preprocessing involved the following steps: 1) Intensity non-uniformity was corrected using N4 bias field correction (the Python-based package “Advanced Normalization Tools”);

(2) Multi-sequence registration was performed using Symmetric Normalization (SyN) in antspyx, with pre-contrast T1WI as the fixed reference;

(3) An adaptive normalizer excluded voxels below the 1% and above the 99% intensity percentiles;

(4) Intensity values were rescaled to a [0, 1] range via min-max normalization (SimpleITK, Python version 3.10).

Radiomic features, including first-order features and textural features using the PyRadiomics package, extracted from multiple imaging sequences: pre-contrast T1WI, AP, PVP, DP, T2WI-FS, and DWI (b = 500 s/mm²). First-order features described the distribution of voxel intensities, and textural features were obtained based on five textural matrices to describe the radiological pattern of the ROI, including the gray level cooccurrence 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). Detailed parameter settings for feature extraction, along with their justifications and supporting references, are provided in Supplementary Material 2 to ensure experimental reproducibility. To account for the collective tumor burden in patients with bHCC, two lesions were treated as a single entity rather than selecting only the largest one. For each patient, radiomics features were extracted independently from both lesions, and a patient-level feature representation was calculated using the arithmetic mean:

where R1 and R2 represent the radiomics features of the first and second lesion, respectively. This averaging-based integration was implemented to capture the comprehensive radiomics phenotype and reflect the total tumor burden of patients with bHCC.

Intra- and inter-class correlation coefficients (ICCs) of radiomics features were calculated for repeatability analysis, with only features having ICC values >0.75 retained for further analysis. A total of 2447 features were, respectively, obtained in each sequence. To further refine the feature set, we excluded features with zero variance across patients and those with a median absolute deviation (MAD) ≤ 1,^20,24 as only features with a MAD > 1 were considered highly variable; the resulting number of highly variable, heterogeneity-related features identified across the six MR sequences is summarized in Table S2.

SNF and Spectral Clustering

Inter-patient similarity matrices were constructed using multi-sequences radiomics features from pre-T1WI, AP, PVP, DP, T2WI-FS and DWI sequences, with Pearson’s correlation employed as the core metric. The similarity matrices, representing patient similarity based on multi-sequence radiomics data, were then fused using the SNF algorithm (R package “SNFtool”¹⁵) separately for the training and testing sets.

Spectral clustering was subsequently performed on the fused similarity network to categorize patients into distinct RIH clusters. Following established methodological frameworks, the optimal number of clusters was determined using the “Spectrum” function within the “R package Spectrum”, which evaluates the algorithm’s internal consistency and the intrinsic structural properties of the patient similarity network. With predefined parameters (t = 10, sigma = 0.4),²⁵ the algorithm identified two clusters as the most robust and stable latent structure for partitioning the cohort, representing the most significant divergence in imaging phenotypes. The association between RIH subtypes and prognosis was assessed using Kaplan–Meier survival curves and Log rank tests.

Statistical Analysis

Statistical analysis was performed using Python (v3.10) with packages “statsmodels” and “scikit-survival”. Continuous variables were presented as mean ± standard deviation or median with interquartile range, depending on statistics distribution determined by Shapiro–Wilk test. Categorical variables were presented as counts with percentages. Statistical significance of intra-group differences was calculated by Chi-squared test, Student’s t test or Mann Whitney U-test. For high-dimensional radiomics data, we integrated the stable and highly variable features (identified via ICC and MAD as described above) into a representative RIH cluster using the SNF algorithm. This dimensionality reduction strategy ensured that the subsequent multivariable analysis remained statistically parsimonious.

Cox proportional hazard regression was applied to assess each variable’s association with OS or RFS. Variables with significant statistical significance in the univariate Cox regression were passed on to multivariate Cox regression, and features shown significance in the multivariate analysis were used to construct prediction model. The selection of variables for the final models was strictly controlled to maintain a balanced EPV ratio, ensuring that the prognostic findings are statistically sound. K-M curve analysis was used to evaluate OS and RFS with Log rank testing for group differences. Concordance index (C-Index), ROC curve analysis, calibration curve were used as evaluation metrics for model performance. 95% confidence interval (CI) was calculated using a bootstrapped method with 1000 random resamples. Models were constructed on the training dataset, then evaluated on the testing set. The p value of 0.05 was considered statistically significant.

Results

Study Population

Figure 1 illustrated the initial enrollment of 489 patients in this retrospective study. Following the application of specific exclusion criteria, 182 patients with bHCC who underwent preoperative multi-parametric MRI were ultimately enrolled. These patients were randomly divided into training and test sets in a 1:1 ratio. Table 1 presented the baseline characteristics of the enrolled patients. For both the training and test cohorts, no significant differences were observed between patients regarding clinical and imaging characteristics.

Table 1 The Baseline Characteristics of Patients with Bifocal HCC

Tables 2 and 3 separately showed the results of Cox regression analysis focusing on clinical and radiologic hallmarks. In the univariable Cox regression analysis of the training cohort, the presence of Unsmooth Margin, pathological satellite nodule, MVI and AFP were independent risk factors for shorter RFS within the clinical and imaging models (all p < 0.05, Table 2). The Cox regression analysis also indicated that larger diameter >5 cm, the presence of Satellite Nodule in imaging features, along with pathological Satellite Nodule, E-S grade, MVI and AFP in clinical features were independent predictors of OS (all p < 0.05, Table 2).

Table 2 Univariate Cox Regression Analyses for RFS and OS in Bifocal HCC Patients

Table 3 Multivariate Cox Regression Analyses for RFS and OS in Bifocal HCC Patients

In the multivariable Cox regression analysis, the predictors of RFS included Radiological Unsmooth Margin (HR, 3.081; 95% CI, 1.386–6.850) in the Imaging model and Pathological Satellite Nodule (HR, 3.310; 95% CI, 1.401–7.817), MVI (HR, 1.646; 95% CI, 0.973–2.787) in the Clinical model (Table 3). For OS, the predictors identified were Radiological Satellite Nodule (HR, 3.556; 95% CI, 1.288–9.821) in Imaging model and MVI (HR, 2.488; 95% CI, 1.081–5.725) in Clinical model (Table 3). Based on these analyses, we constructed the Imaging model and the Clinical model for RFS and OS, respectively. In Clinico-Imaging model, Unsmooth Margin (HR, 2.670; 95% CI, 1.189–6.0), Pathological Satellite Nodule (HR, 3.076; 95% CI, 1.355–6.986) and MVI (HR, 1.713; 95% CI, 1.064–2.758) were independent risk factors for shorter RFS. The Clinico-Imaging model for OS comprised MVI (HR, 1.776; 95% CI, 1.095–2.878) and Radiological satellite nodule (HR, 1.955; 95% CI, 1.070–8.542).

Unsupervised Learning Based on Radiomic Features Among bHCC Patients

In this study, we extracted multi-parametric radiomics from two lesions of bHCC that contained comprehensive radiomics information. The extracted features comprised first-order statistics and textural features derived from the GLCM, GLDM, GLRLM, GLSZM, and NGTDM matrices. The averaged features of the two lesions were used to represent the tumor radiomics characteristics of the patient. Features with zero variance and the MAD less than 1 were filtered to reduce dimensionality and minimize noise in the data (Table S2).

After calculating the similarity between the radiomics of bHCC and constructing the similarity matrix, the patients were classified into different categories using the spectral clustering algorithm (Figure 2). The two diverse subtypes suggested significant differences in the heterogeneity of radiomics features. We found that Radiomics Subtype B is characterized by high radiomics image heterogeneity (high-RIH), whereas Radiomics Subtype A is characterized by low radiomics image heterogeneity (low-RIH) (Figure 3).

Figure 2 Spectral clustering heatmap of similarity network fusion based on radiomics features of bHCC patients (A) in the training cohort and (B) in the test cohort. The color bar on the left of the heatmap identifies the two main subtypes derived from the clustering analysis: Radiomics Subtype A (blue) and Radiomics Subtype B (Orange).

Figure 3 Representative MRI images of high RIH (A and B) and low (C and D) RIH. (A and B) Low-RIH image of a 66-year-old female patient. No recurrence was observed in this patient until the last follow-up. (C and D) High-RIH image of a 65-year-old male patient. Recurrence was detected 12 months after surgery and OS was 32 months.

Characteristics of Radiomics Subtypes

Interestingly, the radiomics-based clustering stratified patients solely based on radiomics features, rather than on any predefined clinical labels. However, the TTD of low-RIH cluster were significantly lower than that of high-RIH (TTD: low-RIH/high-RIH = 5.55 [4.60–7.78] vs. 7.50 [6.50–9.15], p = 0.004) (Table 4). Significant differences were also observed in etiology of liver disease, TP, ALB levels, PLT and the Milan Criteria in training cohort (Hepatitis B virus infection: p = 0.032; TP > 65g/L: p = 0.034; ALB > 35g/L: p = 0.040; PLT > 100 × 10⁹/L: p = 0.019; the Milan Criteria: p < 0.001). Regarding imaging features, the enhancement type also showed a statistically significant difference between the two clusters (p = 0.001).

Table 4 Comparison of Radiomics Subgroups Identified by Unsupervised Analysis

RIH Cluster Reveals Prognostic Value in Predicting Outcomes of bHCC

RIH cluster based solely on robust radiomics features stratified two patient groups with significantly different RFS (log-rank p = 0.0059) and OS (log-rank p = 0.0343) (Figure 4). RIH cluster was identified as an independent predictor of RFS in multivariate Cox regression analyses, after adjusting for the clinical and conventional imaging characteristics (HR = 1.782 [95% CI: 1.189–2.670], p = 0.005; Table 3 and Figure S1). RIH cluster were also significantly associated with OS (HR = 2.008 [95% CI: 1.119–3.605], p = 0.019; Table 3 and Figure S1). The AUC of RIH cluster predicting 1-,3-,5- years OS in training cohort were 0.78, 0.83, 0.81, respectively (Figure S2). The AUC of RIH cluster predicting 1-,3-,5- years RFS in training cohort were 0.75, 0.69, 0.66, respectively.

Figure 4 Kaplan–Meier curves of RFS and OS risk group stratified by RIH clusters. Abbreviations: OS, overall survival; RFS, recurrence-free survival; RIH, radiomics image heterogeneity.

Improved Efficacy of Fusion Model in Predicting Prognosis

The Clinico-Imaging-Radiomics models (RIH-based model) was developed to predict RFS based on Pathological Satellite Nodule (HR = 1.946, 95% CI = 1.094–3.460, p = 0.023), MVI (HR = 1.714, 95% CI = 1.231–2.386, p = 0.001) and RIH Cluster (HR = 1.782, 95% CI = 1.189–2.670, p = 0.005) (Table 3). Additionally, we constructed the RIH-based model to predict OS by incorporating Radiological Satellite Nodule (HR = 1.982, 95% CI = 1.008–3.901, p = 0.047), MVI (HR = 4.350, 95% CI = 2.358–8.028, p < 0.001) and RIH Cluster (HR = 2.008, 95% CI = 1.119–3.605, p = 0.019) (Table 3).

Based on RIH-based model, all patients were stratified into low-risk and high-risk cohorts. Figure S3 illustrated that patients in high-risk group had significantly shorter OS than those in low-risk group (log-rank p < 0.001 in training cohort, long-rank p = 0.001 in test cohort) and RFS (log-rank p = 0.006 in training cohort, long-rank p = 0.007 in test cohort).

The C-index demonstrated that the RIH-based model, with the addition of the RIH cluster, could better predict progression compared to the other three models, particularly the Imaging model (training cohort: form 0.670 [0.557–0.780] to 0.797 [0.711–0.871], p = 0.011; testing cohort: from 0.626 [0.535–0.729] to 0.740 [0.658–0.819], p = 0.030) (Figure 5 and Table S3). Similarly, the RIH-based model enhanced the predictive efficacy for RFS in comparison to others model, particularly the Imaging model (training cohort: form 0.587 [0.527–0.647] to 0.664 [0.589–0.729], p = 0.041; testing cohort: from 0.546 [0.488–0.611] to 0.647 [0.577–0.709], p = 0.034).

Figure 5 Comparison of C-indexes for Clinical, Imaging, Clinco-Imaging and Clinco-Imaging-Rad models in both training and testing sets for (A) OS and (B) RFS.

Discussion

In this study, we innovatively extracted the radiomics features from two lesions in bHCC patients, using the mean radiomics features to represent each patient’s tumor characteristics. We then employed unsupervised machine learning and spectral clustering to establish a non-invasive MRI-based radiomics model (RIH cluster) to assess tumor heterogeneity. The RIH cluster identified two subgroups of bHCC with significantly different RFS and OS. The RIH cluster remained an independent predictor for clinical outcomes even after multivariable Cox regression analysis. Finally, we constructed the RIH-based model to predict outcomes in bHCC patients after resection by integrating RIH clusters with clinical and imaging features.

HCC is characterized by the high tumor heterogeneity, which results the early tumor recurrence and poor clinical outcome. ²⁶ Although radiomics has emerged as an innovative tool for characterizing such heterogeneity ²⁷ and predict outcome ²⁸ in solitary HCC, its application in non-solitary cases has remained constrained. Previous studies have often relied on a single dominant or largest lesion to represent the disease status, potentially overlooking the overall tumor burden and the collective phenotypic state of the patient. To address this limitation, we employed the arithmetic mean to integrate features from both nodules into a unified, patient-level representation. Unlike traditional approaches, this fusion method ensures that the biological contribution of the other lesion is not ignored, thereby providing a more holistic reflection of the total disease burden. While averaging might partially smooth inter-lesion variance, it provides a stable and reproducible estimator of the average aggressiveness across different anatomical locations, reducing the impact of idiosyncratic noise from a single lesion. As an early exploration into fusion-based radiomics for non-solitary HCC, this approach establishes a primary framework for representing multifocal disease; however, we recognize that more sophisticated strategies, such as volume-weighted averaging or multi-instance learning, warrant further investigation in the future. Parallel to this patient-level feature integration, the choice of predictive modeling is crucial for capturing complex biological variations. While previous studies have commonly employed supervised algorithms for prediction, the reliance of supervised machine learning on large training datasets often leaves insufficient data for validation. In contrast, unsupervised machine learning provides an alternative approach that can overcome this limitation.²⁹ SNF, as an unsupervised clustering analysis, comprehensively analyzes data by constructing patient networks from multiple data types and fusing them into a single network that represents the full spectrum of the data.^12,20,30 In the prior study, Su et al²⁰ formulated a noninvasive radiomics methodology using SNF analysis to effectively investigate intratumor heterogeneity in breast cancer. However, to date, there has been a lack of studies focusing on HCC, especially regarding bHCC, which is characterized by a higher degree of tumor heterogeneity.²⁶

Therefore, we utilized the SNF to analyze the tumor heterogeneity of bHCC based on radiomics features of two lesions images. This algorithm capitalized on the complementary nature of multiple MRI sequences, without requiring a priori knowledge, to compute and fuse patient similarity networks derived from each sequence. By utilizing an unsupervised approach like SNF, we reduce the risk of overfitting that often plagues supervised models, as the feature selection and clustering are driven by the intrinsic structure of the data rather than being forced to fit specific clinical labels. Based on the similarity matrix, we constructed a Laplacian matrix of radiomics to identify two distinct subgroups. We observed significant prognostic differences between the subgroups, with low-RIH patients exhibiting better outcomes, and high-RIH patients showing shorter RFS and OS. Unlike biopsy-based methods that are limited by sampling bias and may fail to represent the entire tumor,^31,32 our radiomics framework offers a comprehensive, three-dimensional, and non-invasive evaluation. The RIH clusters were found to be complementary to routine clinical and imaging tests, suggesting that radiomics can capture “invisible” biological variations between multiple lesions.

In the fusion model, Pathological Satellite Nodule and MVI were independent predictors for RFS, while Radiological Satellite Nodule and MVI predicted OS. Our study corroborated the widely recognized role of MVI as a significant risk factor for HCC recurrence following curative resection.^33,34 Satellite Nodule, which represents intrahepatic metastases within the venous drainage area surrounding the tumor, is closely associated with MVI.^35,36 They imply a heightened risk of aggressiveness and recurrence.³⁷ Notably, some predictors that were significant in univariate analysis, such as AFP levels and morphological features like unsmooth margins, lost their statistical significance in the multivariate models. This phenomenon likely stems from the inherent biological and statistical correlations between these variables. For instance, elevated AFP and unsmooth margins are often closely associated with MVI^38,39 and pathological satellite nodules,⁴⁰ which are more direct indicators of tumor invasiveness. When entered into the same multivariate Cox model, the more powerful independent predictors, specifically MVI and the RIH cluster, tend to capture the shared variance, thereby mediating the prognostic value of these overlapping clinical and imaging factors. Our findings further demonstrated that integrating these variables with the RIH cluster significantly enhanced predictive performance. From a clinical perspective, this integrated model may potentially assist clinicians in identifying high-risk patients who might benefit from more aggressive adjuvant therapies or intensified surveillance protocols after surgery. However, its practical implementation in clinical decision-making requires further validation through prospective, multicenter studies to ensure its reliability and generalizability across diverse patient populations.

To ensure the reliability and generalizability of our prognostic models, we implemented several rigorous strategies to address the risk of overfitting. First, we employed a strict dimensionality reduction workflow for radiomics, including repeatability and variability filtering, to distill over a thousand features into a single representative RIH cluster. Second, we adhered to the EPV principle to maintain a balanced ratio between event density and the number of predictors, ensuring that our multivariable Cox models remained statistically parsimonious.

Our study has several limitations. Firstly, our study is an exploratory investigation with a retrospective design, which may introduce potential selection bias. One of the primary limitations is the relatively small sample size from a single center. Although we verified model stability through internal testing and bootstrap resamples to achieve high consistency in performance, the lack of multi-center validation limits the immediate generalizability of our findings to different imaging protocols or diverse patient populations. Specifically, while the use of a single MRI scanner ensured the consistency of radiomics features, the lack of scanner variability in our training set may affect the robustness of the RIH model when applied to images acquired from different manufacturers or hardware configurations. Furthermore, while our inclusion and exclusion criteria were clearly defined to ensure a homogeneous surgical cohort, we excluded patients with macrovascular invasion (eg, portal vein tumor thrombus). Although this was necessary to focus on candidates for radical hepatectomy, it may further limit the generalizability of our findings to the broader HCC population, particularly those in more advanced clinical stages. Secondly, our study focused solely on the relationship between tumor imaging heterogeneity and prognosis. Although previous research has explored the potential links between radiomics and genomic expression levels, the biological significance of imaging heterogeneity remains unclear. Therefore, we plan to conduct integrated studies combining radiomics and genomics to further investigate the potential connections between imaging phenotypes and genomic features, aiming to elucidate the multidimensional information underlying tumor heterogeneity. Thirdly, to address the potential risk of overfitting and the instability of radiomic features, we strictly adhered to radiomics guidelines and checklists to maximize the reliability of radiomic analyses. This included implementing a multi-step feature selection process and maintaining a conservative variable-to-event ratio to ensure the findings were not driven by random noise. Despite these precautions, further validation across different imaging protocols and centers is warranted to ensure the reproducibility of our results.

Conclusion

In this study, we employed unsupervised machine learning to assess radiomics features heterogeneity in patients with bifocal HCC, ultimately constructing a RIH-based model for predicting RFS and OS. While the model demonstrates potential for enhancing prognostic stratification, further external validation in multicenter cohorts is warranted to confirm its clinical applicability as a reliable tool for guiding therapeutic decision-making.