Our study received approval from the Ethics Committee of Ruijin Hospital, Shanghai Jiao Tong University School of Medicine. We retrospectively reviewed all HCC patients who received HAIC from January 2021 to December 2024 at Ruijin Hospital and Ruijin Hospital North.

Patients were included if (1) HCC diagnosed according to the American Association for the Study of Liver Diseases (AASLD) guidelines; (2) age ≥18 years; (3) at least one measurable unresectable lesion; (4) ECOG performance status of 0–1; (5) Child–Pugh class A or B; (6) available contrast-enhanced MRI within 1 month before HAIC and within 1 month after the final treatment cycle. Exclusion criteria included (1) receipt of additional treatments during follow-up; (2) severe cardiac, pulmonary, or renal insufficiency; (3) history of other malignancies; (4) poor-quality imaging or incomplete clinical data.

We included 105 patients who met the eligibility criteria. A flowchart of our patient selection process is provided in Figure 1A.

All patients underwent HAIC using a standardized FOLFOX regimen, repeated every 3 weeks. Under local anesthesia, femoral arterial access was obtained, and a conventional catheter with a microcatheter was selectively advanced into the tumor-feeding hepatic artery. The FOLFOX regimen consisted of oxaliplatin (85 mg/m²) administered via arterial infusion on day 1, followed by leucovorin (400 mg/m²) and fluorouracil (5-FU, 400 mg/m²) delivered intra-arterially, and continuous arterial infusion of 5-FU (2400 mg/m²) over 46 hours (days 2–3). No embolic agents were used.

Follow-up was conducted until the last recorded date (July 18, 2025). Treatment response was evaluated according to the modified Response Evaluation Criteria in Solid Tumors (mRECIST) criteria (12) based on MRI obtained within 1 month after completion of HAIC. Responses were categorized as complete response (CR), partial response (PR), stable disease (SD), or progressive disease (PD). Patients were subsequently grouped as having objective response (OR = CR + PR) or no objective response (nOR = SD + PD). Two radiologists specializing in abdominal interventional radiology (8 and 3 years of experience) independently evaluated all cases while blinded to clinical data; disagreements were resolved through consensus with a highly experienced radiologist (25 years of experience).

The workflow of our study is showed in Figure 1B. Four T1-weighted imaging (T1WI) sequences were included: pre-contrast (PRE), arterial late phase (LAP), portal venous phase (PVP), and delayed phase (DP). MRI scanners consisted of the uMR 790 3.0T, Philips Ingenia 3.0T, and GE SIGNA Architect AIR 3.0T. Detailed imaging parameters are provided in the Supplementary Material S1. Raw DICOM files were obtained from the Picture Archiving and Communication System (PACS) system for analysis.

PVP images obtained prior to HAIC were imported into a software named 3D Slicer (version 5.2.2; www.slicer.org), and tumors were manually segmented by two radiologists (3 and 7 years of experience). Regions of Interest (ROIs) were delineated slice-by-slice.

To minimize variability across scanners and sequences, all images underwent N4 bias-field correction and were resampled to 1×1×1 mm. Registration was then conducted using ANTs (https://antsx.github.io/ANTs/), followed by min–max intensity normalization. All image processing was performed using Python (version 3.8.16; https://www.python.org). Following inter-sequence registration, all MRI sequences were spatially aligned; therefore, the tumor ROIs delineated on the portal venous phase were automatically mapped and correspondingly synchronized to the same anatomical locations on the remaining sequences.

Habitat subregion partitioning was performed using the open-source Python-based HABIT toolkit (https://github.com/lichao312214129/HABIT.git). For each tumor, unsupervised K-means clustering was applied to the voxel-wise intensity features within the ROI. The optimal number of clusters was then determined using the Elbow method by computing a value called the inertia that represents the sum of squared distances from each sample to its nearest cluster center. The inertia was calculated for k-means with k ranging from 2 to 10; we identified k as the point where the reduction of inertia became steady in the inertia–k curve. In our analyses, we selected k = 3. (Figure 2F). The intratumoral heterogeneity of HCC, as characterized by distinct habitat subregions, is shown in Figures 2A–E.

Three types of features were considered: clinical features, radiomic features, and habitat features. Clinical features included demographic information, laboratory test data, and imaging features of tumors and were screened by univariate logistic regression analyses followed by multivariable logistic regression analyses. Radiomic features were extracted from whole-tumor ROIs whereas habitat features were extracted from tumor subregions, using PyRadiomics (https://pypi.org/project/pyradiomics/).

The stability of the extracted features was tested in 30 randomly selected patients by performing three repetitive extractions to calculate its coefficient of variation (CV), which was less than 15%. Features showing CV < 15% were considered stable. To further assess the inter-observer reproducibility of feature extraction, intraclass correlation coefficient (ICC) was used, and radiomic features with ICC < 0.8 were excluded. After normalization process was completed, highly correlated features (|r| > 0.9) were considered redundant and removed. Feature selection was done using recursive feature elimination (RFE) process followed by LASSO regression analysis (13) with 10-fold cross-validation technique after which feature selection was done based on iterative model building approach.

The three types of features were incorporated sequentially to build three models: a clinical model, a clinical-radiomics (CR) model, and a clinical-radiomics-habitat (CRH) model with corresponding risk scores. We internally validated the models using 1,000 times bootstrap resampling. Given the small sample size in our study, bootstrap resampling was implemented to make sure that the evaluation of the models was robust (14, 15).

To assure feature independence and complementary information, we performed Spearman correlation and variance inflation factor (VIF) analyses with the combined features of the CR and CRH models. Severe redundancy indicated by |r|>0.8, and significant multicollinearity was indicated by VIF >10. Here, |r| is the absolute value of the Spearman correlation coefficient measuring the linear association between variables and VIF reflects how much a parameter estimate is increased due to multicollinearity.

Python (version 3.10.19; https://www.python.org) was used for all statistical analyses, with key functions from the scikit-learn, lifelines, matplotlib, and scipy libraries. For continuous variables, we reported those with a normal distribution as mean ± standard deviation (SD), while non-normally distributed ones were summarized as median (interquartile range). Categorical variables were presented as counts and percentages. For inter-group comparisons, independent t-test or Mann–Whitney U test according to distributions was used for continuous variables, and chi-squared test was used for categorical variables.

We assessed model discrimination with receiver operating characteristic (ROC) curves, area under the curve (AUC), and confusion matrices. Calibration and reliability were checked using calibration plots, and decision curve analysis was used to explore clinical utility. When comparing AUCs across models, we applied DeLong’s test. A two-sided p-value < 0.05 was considered statistically significant.

Of 218 HCC patients who received HAIC, 47 were excluded due to concomitant malignancies, 21 for poor image quality or missing clinical data, and 45 for incomplete follow-up (Figure 1). Ultimately, 105 patients were included (mean age 57.3 ± 11.7 years; 97 male [92.4%], 8 female [7.6%]). According to treatment response, 9 achieved CR, 48 PR, 34 SD, and 14 PD. Overall, 57 patients (54.3%) were classified as OR, and 48 (45.7%) as nOR. Baseline characteristics are summarized in Table 1. Compared with the nOR group, the OR group showed notable differences in vascular invasion and tumor diameter; other baseline variables were similar (Table 2). Here, vascular invasion was defined as macrovascular invasion, specifically referring to the involvement of the portal vein and hepatic vein.

A total of 428 radiomics features and 1,140 habitat features were initially extracted per patient. For each sequence, 93 features were extracted from both the whole-tumor ROI and each habitat subregion, including 18 first-order features, 24 GLCM features, 14 GLDM features, 16 GLRLM features, 16 GLSZM features, and 5 NGTDM features. Additionally, 14 shape features were calculated for each ROI, and habitat volume and proportion were recorded. Details of the feature extraction process are provided in Supplementary Material S2. After feature selection, six radiomics features (DP_Original_firstorder_10Percentile, DP_Original_glrlm_ShortRunEmphasis, DP_Original_firstorder_90Percentile, DP_Original_firstorder_Maximum, LAP_Original_firstorder_10Percentile, PRE_Original_firstorder_Maximum) and six habitat features (DP_Habitat3_firstorder_Kurtosis, PVP_Habitat3_firstorder_RootMeanSquared, PVP_Habitat1_firstorder_Minimum, PRE_Habitat1_firstorder_Range, PRE_Habitat2_firstorder_Skewness, DP_Habitat2_firstorder_Minimum) were retained as the most predictive (Figure 3).

Univariate and multivariate logistic regression identified Vascular Invasion (p=0.015) and d_tumor (p=0.037) as independent clinical predictors (Table 3). These two factors were used to construct the clinical model, while the CR model was built by combining clinical features with the selected radiomics features. The CRH model was further developed by incorporating habitat analysis features. Risk scores were calculated using binary logistic regression as follows:

No significant redundancy or multicollinearity was observed (CR: 8 features, mean VIF = 2.445, max=4.930; CRH: 14 features, mean VIF = 2.802, max=5.990), supporting feature independence and model stability (Figure 4).

ROC analysis showed that the clinical, CR, and CRH models achieved AUCs of 0.633 (95% CI: 0.525–0.742), 0.747 (95% CI: 0.653–0.841), and 0.771 (95% CI: 0.682–0.860), respectively. DeLong tests confirmed that the CR model showed significantly better performance than the clinical model (Z=–1.972, p=0.05), and the CRH model performed even better (Z=–2.386, p=0.02). Although the difference between the CR and CRH models did not reach statistical significance (Z = –1.071, p = 0.28), bootstrap validation indicated stable model performance (SD of AUC: CR = 0.048, CRH = 0.046). By incorporating six additional habitat features, the CRH model enhanced the representation of tumor microenvironment heterogeneity, increasing feature dimensionality and providing more biologically meaningful prognostic information (Figure 5).

We validated the performance of the CRH model in classifying patients. Figure 6A depicts a confusion matrix with satisfactory overall accuracy. The difference in predicted risk scores between OR and nOR groups is statistically significant (p<0.001), with a median risk score of 0.626 in the OR group compared with that of 0.486 in the nOR group (Figure 6B). This indicated that our model can distinguish between patients with different clinical outcomes. Figure 6C shows the Sankey diagram analysis, which revealed good agreement between our model’s prediction of early response and patients’ actual outcomes.