This study included 312 HCC patients from Hunan Provincial People’s Hospital/The First Affiliated Hospital of Hunan Normal University (referred to as dataset A) and 86 patients from the Third Affiliated Hospital of Sun Yat-sen University (referred to as dataset B) between February 2018 and October 2023. Institutional review board approval was obtained from each participating center. The inclusion criteria were: (1) a solitary tumor, (2) pathologically confirmed HCC following surgical resection, (3) available information on MVI status and grade, and (4) preoperative Gd-EOB-DTPA MRI conducted within two weeks prior to surgery. The exclusion criteria were: (1) macrovascular invasion, (2) prior HCC treatment before MRI (e.g., radiofrequency ablation, microwave ablation, or transcatheter arterial chemoembolization), (3) tumors larger than 10 cm in maximum diameter (as previous studies (23–25) have shown a greater likelihood of MVI in such cases), (4) inadequate MRI quality, and (5) missing pathological or clinical data. The two datasets were combined to form a total cohort of 398 patients. A randomly selected 20% of this cohort was designated as a fixed test set, and the remaining 80% of the cases were used for 5-fold cross-validation. Clinical data, including variables such as age, gender, etiology, cirrhosis, MVI status, and pathological differentiation, were extracted from electronic medical records and are summarized in Table 1 .

For dataset A, MRI scans were performed using 1.5T or 3.0T MRI machines from GE (Signal Greator, Premier), Philips (Achieva, Ingenia), and Siemens (Magnetom Trio, Magnetom Prisma, Vida). For dataset B, MRI scans were conducted using 1.5T or 3.0T MRI systems from GE (Optima MR360, Signa Excite, Discovery MR750, Signa Architect), Philips (Achieva), Siemens (Magnetom Prisma), and United Imaging (uMR790). All patients underwent fat-saturated T1-weighted pre-contrast scans, followed by scans in the AP, PVP, and HBP. AP images were acquired 20-30 seconds after gadolinium contrast injection, PVP images 60-70 seconds post-injection, and HBP images were taken 20-30 minutes after contrast administration.

Two experienced radiologists, each with over ten years of experience in MRI diagnostics, independently and blindly assessed the radiological features of the tumors. A consensus was reached regarding the following six characteristics: (1) tumor size, (2) non-smooth tumor margin, (3) radiological capsules (26), (4) intratumoral artery (27), (5) arterial peritumoral enhancement (28), and (6) peritumoral hypointensity on HBP (29).

Primary tumors from the HBP images were manually delineated by a senior radiologist, Yuli Zeng, with over 15 years of experience, using the ITK-SNAP 3.4 software platform (www.itksnap.org). The AP and PVP images were then registered to the HBP images, which served as reference images. To correct for low-frequency intensity nonuniformity, N4 bias field correction (30) was applied to all images. All images were resampled to an isotropic voxel size of 1 × 1 × 1 mm³ using B-spline interpolation, while the delineated tumor masks were resampled using nearest neighbor interpolation.

For each sequence (i.e., AP, PVP, and HBP), 105 radiomics features were extracted from the original images using the PyRadiomics package (31). These features included 14 shape features and 91 texture features. Additionally, texture features were extracted using wavelet filters (HHH, HHL, HLH, HLL, LHH, LHL, LLH, LLL) and Laplacian of Gaussian filters with sigmas of 2.0, 3.0, 4.0 and 5.0, which resulted in a comprehensive set of 1197 features for each sequence.

The images of HCC patients typically exhibit AP hyperenhancement followed by washout in the PVP (16). To assess radiomic changes during dynamic contrast enhancement, the features from the AP images were compared to those from the PVP images using the following Equation 1:

Where F D e l t a represents the change in features between AP images and PVP images, which is time-related and predictive to MVI, F A P represents the features extracted from AP images, and F P V P refers to the features generated from the PVP images.

Radiomics offers detailed insights into tumor phenotypes and the tumor microenvironment (32). To capture intra-tumoral heterogeneity, we performed intra-tumoral habitat partitioning in two steps: case-based clustering and subregion feature extraction. Case-based clustering was conducted independently for each tumor using the k-means algorithm with squared Euclidean distances between voxel intensities. The number of clusters was set to five due to the small tumor volumes in this study. The clustering process is performed using the in-house nnFAE software. In the subregion feature extraction step, radiomics features, including histogram, GLCM, GLRLM, NGTDM, GLSZM and GLDM, were extracted from each subregion without applying additional filters.

Subsequently, we applied a Gaussian mixture model to perform unsupervised clustering of radiomic features across all tumor habitats. The optimal number of clusters, representing the diversity of the tumor ecosystem, was determined using the Bayesian Information Criterion (BIC) (33), which generated an intra-tumoral ecological diversity (iTED) feature vector for each patient, which could then be used for further analysis. Each iTED feature reflects the optimal number of clusters corresponding to specific radiomic features. For example, the iTED_entropy feature represents the optimal number of clusters for assessing tumor heterogeneity, using traditional entropy as a metric. While conventional entropy measures the unpredictability or variability of image values, iTED_entropy quantifies the complexity of intra-tumoral heterogeneity by evaluating entropy at the cluster level. This iTED feature vector provides a novel approach to tumor characterization, potentially offering new insights into tumor behavior and structure (34).

The generation of the conventional radiomics features, delta radiomics features and iTED features is shown in Figure 1 . To address potential variability in the radiomics features caused by differences in imaging protocols across the two centers, the ComBat harmonization method (35) was applied.

To preselect features with high stability, we simulated delineation perturbations based on the training cohort. Morphological operations, including dilation and erosion, were applied slice-by-slice using a circular structural element with distances of 1 mm and 2 mm. This process generated four distinct VOIs, labeled D1, D2, E1, and E2. To assess feature stability, we used the inter-class correlation coefficient (ICC) (36), classifying features as having high (ICC ≥ 0.75), moderate (0.75 > ICC ≥ 0.50), or low (ICC < 0.50) stability. Following established guidelines (37), we applied the ICC (2,1) model as defined by Shrout and Fleiss (38) and calculated the ICC using the Pingouin statistical library (https://github.com/raphaelvallat/pingouin).

ICC was calculated for all five ROIs, including the original ROI, and the dilated (1 mm and 2 mm) and eroded (1 mm and 2 mm) ROIs. Only first-order and textural features were evaluated for stability, with features having an ICC greater than 0.75 selected for the next stage of the feature selection pipeline. Shape-related features, however, were directly retained and included in the pipeline without undergoing stability evaluation.

For both conventional radiomic features and delta radiomics features, the selection process began with retaining features that showed significant differences between patients with and without MVI, as determined by the Mann–Whitney U-test. Next, we selected features that achieved an Area Under the Curve (AUC) greater than 0.60 in univariate logistic regression analysis. To further refine the feature set, the minimum redundancy and maximum relevancy (mRMR) method was applied to eliminate redundant and irrelevant features. Finally, the Least Absolute Shrinkage and Selection Operator (LASSO) algorithm was used to reduce the feature set to only the most predictive features. For conventional radiomic features, this process resulted in the selection of 18 features from AP images, 13 from PVP images, and 10 from HBP images. For delta radiomics features, 14 features were retained for model development.

For the iTED features, we first applied z-score normalization to standardize the features and removed those with minimal variance. Next, features with an AUC greater than 0.55 in univariate logistic regression were retained. The LASSO algorithm was then applied, leaving 4 features from AP images, 4 from PVP images, and 9 from HBP images for further analysis.

The details of the selected features and their corresponding ICC values are provided in Supplementary Tables S1 - S7 of the Supplementary Material .

Differences in clinical and radiological characteristics between the training and test cohorts were assessed using the Mann-Whitney U test for continuous variables and the chi-square test for categorical variables. In the training dataset, five-fold cross-validation with stratified sampling was performed to ensure consistent category proportions. A random forest (RF) model was constructed to classify patients with or without MVI, and Bayesian optimization (39) was applied to fine-tune the model’s hyperparameters.

Ultimately, eight model types were constructed using (1) A clinical-radiological model (CR model) using demographic, pathological, and radiological features, (2) conventional radiomic features from AP images (M_CVT-AP), (3) conventional radiomic features from PVP images (M_CVT-PVP), (4) conventional radiomic features from HBP images (M_CVT-HBP), (5) delta radiomics features (M_Delta), (6) iTED features from AP images (M_iTED-AP), (7) iTED features from PVP images (M_iTED-PVP), and (8) iTED features from HBP images (M_iTED-HBP). Then, we developed a fusion_R model to combine the predictions from the above seven radiomics models using a stacking algorithm (40). Furthermore, we constructed a fusion_CR model, which combines the radiomics models with the CR model. The workflow of developing the fusion model by stacking algorithm is shown in Figure 2 .

The performance of the models in predicting MVI was evaluated using the Area Under the Receiver Operating Characteristic Curve (AUC) with 95% confidence intervals, as well as accuracy (ACC), sensitivity, and specificity. Delong’s test was employed to compare the AUCs of different models, with statistical significance set at p < 0.05.

A total of 398 HCC patients were included in the study, with 318 patients in the training dataset (mean age 56 years; 271 males, 47 females) and 80 patients in the testing dataset (mean age 55 years; 67 males, 13 females). As shown in Table 1 , the two sets were well-balanced as there were no statistically significant differences in clinical-radiological characteristics either between the training and testing sets or within each set (p = 0.104-1.000).

Univariate analysis identified one demographic factor (HBV infection), four radiological factors (tumor size, nonsmooth tumor margin, arterial peritumoral enhancement, and peritumoral hypointensity on HBP), and two pathological factors (degree of differentiation and cirrhosis, stage 4 fibrosis) as being associated with MVI in the training set. Multivariable analysis ( Table 2 ) revealed that tumor size (OR = 1.20, 95% CI: 1.06–1.36, p < 0.001), nonsmooth tumor margin (OR = 2.84, 95% CI: 1.63–5.06, p < 0.001), and cirrhosis (stage 4 fibrosis) (OR = 2.02, 95% CI: 1.19–3.50, p = 0.01) were significant predictors of MVI and were incorporated into the clinical-radiological (CR) model. The CR model achieved an AUC of 0.784 (95% CI: 0.766–0.802) in the training dataset, 0.722 (95% CI: 0.661–0.784) in the validation dataset, and 0.677 (95% CI: 0.610–0.744) in the testing dataset.

The performance of models based on conventional radiomics features is presented in Table 3 . The M_CVT-AP model achieved an AUC of 0.723 (95% CI: 0.655–0.789) in the testing cohort, demonstrating superior diagnostic performance compared to the M_CVT-PVP model (AUC = 0.672, 95% CI: 0.611–0.734) and the M_CVT-HBP model (AUC = 0.620, 95% CI: 0.601–0.639). The higher signal contrast within and between tumors on AP images, due to significant enhancement, likely contributed to this improved performance. Conversely, the M_CVT-HBP model exhibited lower sensitivity (0.496), likely due to the minimal signal variation observed in HCC lesions during the HBP phase. Figure 3 illustrates the mean receiver operating characteristic (ROC) curves, the probability distribution of classes, and the confusion matrices for the M_CVT-AP, M_CVT-PVP, and M_CVT-HBP models.

The performance of the M_Delta model is summarized in Table 4 . The M_Delta model achieved an AUC of 0.707 (95% CI: 0.678–0.735), outperforming the M_CVT-PVP model (AUC = 0.672, 95% CI: 0.611–0.734) but falling short of the M_CVT-AP model (AUC = 0.723, 95% CI: 0.655–0.789). Importantly, the M_Delta model exhibited a higher sensitivity (0.672, 95% CI: 0.590–0.753) compared to the M_CVT-AP model (0.635, 95% CI: 0.508–0.763) and the M_CVT-PVP model (0.520, 95% CI: 0.448–0.593), suggesting that delta radiomics features, which capture time-related changes, provide valuable predictive information for MVI beyond what is offered by AP and PVP features alone.

As shown in Table 5 , the M_iTED-HBP model (AUC = 0.727, 95% CI: 0.706–0.749) outperformed the M_iTED-AP model (AUC = 0.613, 95% CI: 0.575–0.651) and the M_iTED-PVP model (AUC = 0.691, 95% CI: 0.676–0.707). Interestingly, the performance of the iTED models in different phases was the opposite of that observed in the conventional radiomics models. However, the sensitivity of the iTED models was relatively low (Sensitivity = 0.460–0.545) in the testing cohort. Representative MRI images are displayed in Figure 4 . MVI-positive HCC cases were found to have a higher proportion of habitat-4 in the tumor center. Figure 5 provides a visual assessment of clustering effectiveness while emphasizing intra-tumor heterogeneity. In the figure, habitat-1 represents regions with high enhancement; habitat-2 corresponds to areas with medium to medium-high enhancement; habitat-3 includes regions with low or no enhancement; habitat-4 highlights areas of cystic degeneration and necrosis; and habitat-5 encompasses regions with medium-low enhancement.

As shown in Table 6 , the fusion_R model demonstrated excellent discriminatory performance, achieving an AUC of 0.823 (95% CI: 0.816–0.831) and an accuracy of 0.775 (95% CI: 0.753–0.796) in the testing cohort. The fusion_R model outperformed the base classifiers (M_CVT-AP, M_CVT-PVP, M_CVT-HBP, M_Delta, M_iTED-AP, M_iTED-PVP, and M_iTED-HBP) across nearly all evaluation metrics in both the validation and testing cohorts. Figures 6a, b shows the ROC and precision-recall (PR) curves for the fusion_R model alongside the best-performing conventional radiomics model (M_CVT-AP), delta radiomics model (M_Delta), and iTED model (M_iTED-HBP).

The performance of the fusion_R model was comparable to that of the fusion_CR model (AUC = 0.823 vs. AUC = 0.830, p = 0.718), suggesting that while clinical-radiological features had predictive value, their contribution to enhancing the radiomics-based prediction was minimal. The ROC and PR curves for the fusion_R and fusion_CR models are displayed in Figures 6c, d . Additionally, we applied sigmoid calibration to the fusion_R and fusion_CR models. However, the calibration resulted in no significant improvement in performance, with the fusion_R model showing a slight change (pre: 0.823 vs. post: 0.825) and the fusion_CR model exhibiting minimal variation (pre: 0.830 vs. post: 0.828). Supplementary Figure S1 in the Supplementary Material presents the model calibration curves both before and after calibration.

The p values of the Delong test between the different models are shown in Figure 7 . The fusion model significantly improved MVI discrimination compared to every other model (p = 0.000–0.050) except the M_CVT-AP model (fusion_R: p = 0.101, fusion_CR: p = 0.054). No significant differences were found in the performance of images from different phases, whether using conventional radiomics features (p = 0.096–0.420) or iTED features (p = 0.106–0.744). Additionally, for images from the same phase, there was no significant difference in performance between conventional radiomics features and iTED features (p = 0.158 for AP images, p = 0.844 for PVP images, and p = 0.157 for HBP images). These findings suggest that, although the predictive power of iTED features is not as strong as conventional radiomics features, iTED features have the potential to replace conventional radiomics features to some extent.