Habitat imaging on contrast-enhanced magnetic resonance imaging predicts early response to transarterial chemoembolization in hepatocellular carcinoma

Abstract

BACKGROUND

Transcatheter arterial chemoembolization (TACE) is a core treatment for hepatocellular carcinoma (HCC), but significant variability in therapeutic responses exists due to intratumoral heterogeneity (ITH). Conventional predictive methods ignore subregional tumor heterogeneity, limiting their accuracy and generalizability. Habitat imaging enables non-invasive quantification of ITH by characterizing tumor subregional diversity, yet few studies have integrated this technique with clinical and radiomic features to enhance TACE response prediction.

AIM

To develop and validate a combined model integrating clinical, radiomics, and ITH features for predicting early response to TACE in HCC patients.

METHODS

A total of 223 HCC patients were divided into training (n = 107), internal validation (n = 46), and external validation (n = 70) cohorts. Preoperative magnetic resonance imaging images were processed via simple linear iterative clustering superpixel segmentation and Gaussian mixture model clustering (Bayesian information criterion-optimized) to quantify tumor ecological diversity. The ITH and radiomic features were extracted and selected by recursive feature elimination with cross-validation. A combined model was constructed and evaluated using receiver operating characteristic curves, calibration curves, and decision curve analysis. We employed SHapley Additive exPlanations analysis to quantify the contributions of each feature to the combined model.

RESULTS

The combined model achieved the highest area under the curve (AUC) of 0.97 (95% confidence interval: 0.93-0.99) in the training set, outperforming clinical (AUC = 0.86) and radiomics (AUC = 0.71) models. It maintained robust performance in internal (AUC = 0.91) and external (AUC = 0.93) validation cohorts, with good calibration and superior clinical net benefit in decision curve analysis.

CONCLUSION

Integrating ITH features significantly enhances the prediction of TACE response in HCC. The combined model exhibits excellent discriminative power and generalizability, fulfilling the objective of providing a reliable tool for personalized treatment decision-making.

Key Words: Hepatocellular carcinoma; Transcatheter arterial chemoembolization; Intratumoral heterogeneity; Habitat imaging; Magnetic resonance imaging

Core Tip: This study builds a combined model integrating magnetic resonance imaging-based intratumoral heterogeneity habitat imaging with clinical/radiomics features to predict early transcatheter arterial chemoembolization response in hepatocellular carcinoma. Multi-center validation (153 cases vs 70 cases) confirms it outperforms single models, supporting personalized transcatheter arterial chemoembolization decision-making.

INTRODUCTION

Hepatocellular carcinoma (HCC) is one of the common pathological subtypes of primary liver cancer. It ranks 6^th among the most prevalent cancers worldwide and 3^rd among the causes of cancer-related deaths[1]. Due to the insidious onset of HCC, approximately 80.0% of patients are diagnosed at an advanced stage[2], which significantly limits the indications for radical resection[3]. According to the Barcelona Clinic Liver Cancer (BCLC) staging system, transcatheter arterial chemoembolization (TACE) is a key treatment modality for patients with advanced-stage HCC[4,5]. Additionally, for early-stage HCC patients who do not benefit from recommended treatments, the BCLC staging system recommends TACE as an alternative treatment option. Although “TACE combined with targeted therapy and immunotherapy” has become the standard treatment paradigm in the current clinical context, “TACE-only” remains the standard treatment option for specific populations: For patients with small tumor burden, no vascular invasion, and good liver function reserve, “TACE-only” is still the first-line treatment. For some patients with limited economic conditions or comorbidities, “TACE-only” is their only feasible treatment choice.

However, due to the high heterogeneity of HCC, there is significant variability in the early therapeutic responses to TACE among different patients in clinical practice[6]. Intratumoral heterogeneity (ITH) refers to the heterogeneity in genetics, phenotypes, and biological behaviors within a tumor. It manifests as significant differences in cell composition, vascular distribution, metabolic activity, and invasive potential among different subregions or “habitats” within the tumor tissue. This heterogeneity leads to differential responses to treatment[7] and influences clinical treatment decisions. Nevertheless, the current clinical assessment of ITH in HCC and the prediction of therapeutic responses to TACE still face significant challenges[8].

Previous literature reports indicate that common methods for predicting TACE efficacy include routine radiological features of tumors, serological markers, clinical staging, and traditional radiomics models based on computed tomography (CT) or magnetic resonance imaging (MRI)[9,10]. However, the assessment of routine radiological features has certain subjectivity; the predictive value of serological markers and clinical staging is limited; and traditional radiomics models based on CT or MRI extract features from the entire tumor while ignoring the heterogeneity of subregions within the tumor[11]. These factors collectively restrict their clinical application.

Dynamic contrast-enhanced MRI is a sensitive method for assessing the efficacy of TACE treatment in HCC. Compared to CT, dynamic contrast-enhanced MRI provides higher soft-tissue resolution, allowing for clearer visualization of the tumor’s internal structure and regional perfusion differences, which is crucial for quantifying ITH. Habitat imaging technology is a radiomics approach based on the biological microenvironment. Building on traditional radiomics, it divides the tumor into multiple subregions, explores the potential biological information of the tumor, and captures the tumor’s heterogeneous information as comprehensively as possible[12,13]. It clusters tumor cell populations with similar characteristics using quantitative imaging biomarkers, thereby visualizing and quantifying the ITH. Currently, habitat analysis technology has been widely used in research on efficacy and prognosis prediction models for diseases such as glioblastoma and breast cancer[14,15]. Additionally, there have been literature reports on applying this technology to predict HCC prognosis-related indicators (e.g., microvascular invasion), achieving promising results[16].

To our knowledge, no previous study has compared clinical and routine radiological feature models, traditional radiomics models, and ITH-based habitat analysis models for HCC, nor combined these models to develop a new integrated model for predicting the early efficacy of TACE in HCC patients. Therefore, the purpose of this study is to construct a predictive model based on ITH habitat analysis using regional radiomics analysis of preoperative MRI images of HCC patients undergoing TACE, and to examine whether the new model, proposed by combining this method with clinical and routine radiological features as well as traditional radiomics, can effectively predict the early efficacy of TACE in HCC patients.

MATERIALS AND METHODS

Study design and patient criteria

The study was designed in accordance with the Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis guidelines. It strictly adheres to the guidelines’ standardized requirements for data sources, model construction, and validation procedures, and is classified as a type III study based on a clinical retrospective cohort. The dataset used in this study includes data from patients diagnosed with HCC and treated with TACE, collected from two tertiary medical centers between January 2018 and June 2025. A total of 153 patients were enrolled from institution A, and they were randomly divided into an internal training set and an internal validation set at a ratio of 7:3. Institution B served as the external validation set, with 70 patients enrolled.

Inclusion and exclusion criteria

Inclusion criteria: (1) Patients diagnosed with HCC either by histopathological confirmation or in accordance with the latest clinical practice guidelines for the diagnosis and management of HCC issued by the American Association for the Study of Liver Diseases; (2) Patients undergoing TACE for the first time, with no prior receipt of other treatments for HCC; (3) Patients who underwent contrast-enhanced MRI of the upper abdomen within one month before TACE, and had follow-up contrast-enhanced MRI or contrast-enhanced CT images within 4-6 weeks after TACE; and (4) Patients with complete clinical data.

Exclusion criteria: (1) Patients who received other treatments (e.g., surgical resection, radiofrequency ablation, chemotherapy) prior to TACE; (2) Missing or poor-quality MRI images; (3) Missing follow-up data after TACE; and (4) Patients with concomitant life-threatening malignant tumors. The process of patient enrollment and exclusion is shown in Figure 1.

Figure 1 Enrollment and exclusion flowchart of patients with hepatocellular carcinoma. The left side shows the screening process for institution A, where 1631 initial patients with hepatocellular carcinoma who underwent transcatheter arterial chemoembolization were screened. A total of 153 patients were finally enrolled and divided into an internal training set (n = 107) and an internal validation set (n = 46). The right side shows the screening process for institution B, where 743 initial patients with hepatocellular carcinoma who underwent transcatheter arterial chemoembolization were screened. A total of 70 patients were finally enrolled, serving as the external validation set. HCC: Hepatocellular carcinoma; AASLD: American Association for the Study of Liver Diseases; TACE: Transcatheter arterial chemoembolization; MRI: Magnetic resonance imaging.

Patient characteristics

Additionally, a data table was designed in this study to collect routine MRI radiological features and clinical information that may predict the therapeutic response to TACE, specifically including: (1) Patient demographic characteristics and clinical information, such as gender, age, history of hepatitis B, Child-Pugh classification, Eastern Cooperative Oncology Group (ECOG) performance status, and BCLC staging system; (2) Laboratory indicators, namely alkaline phosphatase (ALP), alpha-fetoprotein, alanine aminotransferase, aspartate aminotransferase, total bilirubin, and albumin; (3) Routine MRI radiological features, like the number of tumors, maximum tumor diameter, venous tumor thrombus invasion, and radiological evidence of cirrhosis; and (4) Patient follow-up status and the assessment results of therapeutic response.

This study retrospectively reviewed data from 1631 and 743 patients (from two institutions, respectively) who were diagnosed with HCC by pathology or imaging and had undergone TACE. After exclusion based on predefined criteria, 153 and 70 patients were finally enrolled in the study. Among these, 153 patients from institution A were divided into an internal training set (n = 107) and an internal validation set (n = 46), while 70 patients from institution B served as the external validation set.

The baseline characteristics of the patients are presented in Table 1. The demographic and clinical characteristics were well-balanced across all cohorts, with no statistically significant differences (P > 0.05). Specifically, the proportion of male patients exceeded 75.0% in each cohort, the proportion of cirrhosis ranged from 60.0% to 72.0%, and over 80.0% of patients had a history of hepatitis B. Additionally, Child-Pugh classification was predominantly class A and B in all cohorts. These characteristics were consistent with the common clinical phenotypes of HCC, further confirming baseline consistency across the cohorts.

Table 1 Baseline characteristics of all cohorts, n (%).

Features	Training cohort (n = 107)	Validation cohort (n = 46)	External validation cohort (n = 70)	P value
Age (years), mean ± SD	63.24 ± 11.72	61.80 ± 13.07	65.74 ± 9.59	0.160
Sex				0.791
Male	87 (81.3)	37 (80.4)	54 (77.1)
Female	20 (18.7)	9 (19.6)	16 (22.9)
Size > 5 cm				0.427
Yes	44 (41.1)	17 (37.0)	22 (31.4)
No	63 (58.9)	29 (63.0)	48 (68.6)
Tumor number				0.168
Single	39 (36.4)	17 (37.0)	35 (50.0)
Multiple	68 (63.6)	29 (63.0)	35 (50.0)
Vascular tumor thrombus				0.154
Yes	35 (32.7)	11 (23.9)	14 (20.0)
No	72 (67.3)	35 (76.1)	56 (80.0)
Cirrhosis				0.187
Yes	77 (72.0)	28 (60.9)	42 (60.0)
No	30 (28.0)	18 (39.1)	28 (40.0)
ALP > 125 U/L				0.785
Yes	43 (40.2)	19 (41.3)	25 (35.7)
No	64 (59.8)	27 (58.7)	45 (64.3)
AFP > 400 ng/mL				0.609
Yes	24 (22.4)	11 (23.9)	12 (17.1)
No	83 (77.6)	35 (76.1)	58 (82.9)
AST > 40 U/L				0.178
Yes	63 (58.9)	30 (65.2)	34 (48.6)
No	44 (41.1)	16 (34.8)	36 (51.4)
ALT > 50 U/L				0.836
Yes	31 (29.0)	14 (30.4)	18 (25.7)
No	76 (71.0)	32 (69.6)	52 (74.3)
TBIL < 19 μmol/L				0.562
Yes	56 (52.3)	24 (52.2)	42 (60.0)
No	51 (47.7)	22 (47.8)	28 (40.0)
ALB > 40 g/L				0.087
Yes	21 (19.6)	13 (28.3)	24 (34.3)
No	86 (80.4)	33 (71.7)	46 (65.7)
History of hepatitis B				0.432
Yes	89 (83.2)	41 (89.1)	56 (80.0)
No	18 (16.8)	5 (10.9)	14 (20.0)
ECOG				0.673
0	88 (82.2)	41 (89.1)	60 (85.7)
1	14 (13.1)	4 (8.7)	9 (12.9)
2	5 (4.7)	1 (2.2)	1 (1.4)
BCLC				0.369
0	4 (3.7)	1 (2.2)	6 (8.6)
A	37 (34.6)	18 (39.1)	29 (41.4)
B	24 (22.4)	14 (30.4)	17 (24.3)
C	40 (37.4)	11 (23.9)	17 (24.3)
D	2 (1.9)	2 (4.3)	1 (1.4)
Child-Pugh				0.102
A	61 (57.0)	29 (63.0)	53 (75.7)
B	44 (41.1)	15 (32.6)	16 (22.9)
C	2 (1.9)	2 (4.3)	1 (1.4)

ALP: Alkaline phosphatase; AFP: Alpha-fetoprotein; AST: Aspartate aminotransferase; ALT: Alanine aminotransferase; TBIL: Total bilirubin; ALB: Albumin; ECOG: Eastern Cooperative Oncology Group; BCLC: Barcelona Clinic Liver Cancer.

MRI protocol

All MRI examinations were performed using a 3.0 T magnetic resonance scanner (Signa, Architect, GE Healthcare, IL, United States) with a 16-channel abdominal coil and additional respiratory gating for the abdomen. All patients were placed in the supine position and fasted (no food or water intake) for 4-6 hours before the examination. The scanning sequences were as follows: (1) Axial spoiled gradient-echo sequences for liver acquisition with volume acceleration (LAVA) based on the FLEX fat-water separation technique, including in-phase T1-weighted imaging and out-of-phase T1-weighted imaging; (2) Axial fat-saturated T2-weighted imaging using a fast spin-echo sequence; (3) Axial diffusion-weighted imaging-echo planar imaging with b-values of 50 seconds/mm² and 1000 seconds/mm²; and (4) Axial dynamic contrast-enhanced scans for LAVA: Gadodiamide (GE Healthcare, IL, United States) was administered as a bolus via an antecubital vein using a high-pressure injector at a rate of 3 mL/second with a dose of 0.1 mmol/kg. Scans were performed at 30 seconds, 50 seconds, and 90 seconds after contrast agent injection, respectively, to obtain images of the arterial phase, portal venous phase, and venous phase. All radiomic features in this study were exclusively extracted from arterial-phase images of the dynamic LAVA sequence. MRI acquisition parameters are detailed in Supplementary Table 1.

TACE treatment protocol

All enrolled patients successfully underwent TACE treatment, including conventional TACE (cTACE) and drug-eluting bead TACE (DEB-TACE). Physicians selected either the cTACE or DEB-TACE protocol based on the patient’s liver function, tumor size and number, whether the tumor was diffuse or focal, and the tumor’s location in the liver. The TACE protocol recommended by physicians required patient consent. All treatment procedures were performed in accordance with current clinical practice guidelines and carried out by interventional radiologists with more than 10 years of clinical experience to ensure standardization of the treatment process; the criteria for successful embolization were defined as no contrast medium staining in the tumor area observed on post-embolization angiography.

For cTACE, a mixture of lipiodol (Hengrui Medicine, Jiangsu, China) and 50 mg of epirubicin, gelatin sponge particles (Ailikang, Hangzhou, China), and blank microspheres (Mairuitong, Beijing, China) were used as embolic agents. Under the monitoring of digital subtraction angiography (Philips Allura FD20, Philips Healthcare, Netherlands; or Siemens Artis Q Ceiling, Siemens Healthineers AG, Germany), 2.4Fr microcatheters (Boston Scientific, MA, United States) or 2.6Fr microcatheters (Hengrui Medicine, Jiangsu, China) were used to inject the embolic agents into the hepatic arterial system. For DEB-TACE, microspheres (Callispheres; Hengrui Medicine, Jiangsu, China) with a diameter of 100-300 μm or 300-500 μm were used as carriers (selected based on tumor size) and loaded with 50 mg of epirubicin; if necessary, 100-300 μm gelatin sponge particles and blank microspheres were additionally used for embolization. Meanwhile, during both cTACE and DEB-TACE procedures, 150 mg of oxaliplatin dissolved in glucose solution was alternately injected with the above-mentioned embolic materials. The treatment process was similar for both protocols, with the complete occlusion of the tumor-feeding artery as the treatment endpoint. All patients were hospitalized for postoperative supportive care after TACE and received routine management, including fluid replacement, antiemetic administration, pain management, and monitoring of liver function changes.

Clinical efficacy assessment and image quality evaluation

In this study, patients who underwent TACE underwent follow-up contrast-enhanced CT or MRI within 4-6 weeks after treatment to assess early therapeutic efficacy. The evaluation time point (4-6 weeks) is based on the recommended standard in the Guidelines for Transarterial Chemoembolization in Hepatocellular Carcinoma issued by the Chinese College of Interventionalists[17]. Formulated on the basis of global evidence-based data and clinical consensus, this guideline explicitly identifies 4-6 weeks as the optimal time point for early efficacy evaluation after TACE. Efficacy of TACE for HCC was evaluated using the modified Response Evaluation Criteria in Solid Tumors issued by American Association for the Study of Liver Diseases, by comparing preoperative contrast-enhanced MRI images (obtained within 1 month before TACE) with postoperative follow-up contrast-enhanced MRI or CT images (obtained within 4-6 weeks after TACE). The specific criteria were as follows: (1) Complete response: Disappearance of arterial-phase contrast enhancement in target lesions; (2) Partial response (PR): ≥ 30% reduction in the sum of diameters of target lesions (assessed by arterial-phase contrast enhancement); (3) Stable disease: Reduction not meeting PR criteria or increase not meeting progressive disease (PD) criteria; and (4) PD: ≥ 20% increase in the sum of diameters of target lesions (assessed by arterial-phase contrast enhancement) or the appearance of new lesions. Based on the evaluation results, patients were divided into the TACE response group (including complete response and PR) and TACE non-response group (including stable disease and PD).

The therapeutic response in this study was assessed by two radiologists with more than 5 years of image interpretation experience who were blinded to patients’ clinical information; they evaluated preoperative and postoperative follow-up images, and also conducted image quality evaluation on preoperative contrast-enhanced MRI images. A 5-point Likert scale was used to score the image series in the following aspects: Overall image quality (scored from 1 = poor quality with extensive image noise to 5 = excellent quality with no perceptible noise), lesion depiction (scored from 1 = no visual depiction to 5 = perfect contour depiction), image sharpness (scored from 1 = significant blurriness to 5 = no obvious blurriness), and image noise (scored from 1 = extensive noise to 5 = no noise). An image score of ≥ 3 was considered to meet diagnostic requirements; only images that met diagnostic quality standards were used for subsequent analysis. In case of inconsistent assessments, the two physicians discussed the discrepancies until a consensus was reached.

Image preprocessing and tumor segmentation

As illustrated in Figure 2, the methodological workflow is divided into two parts, encompassing preprocessing and feature acquisition as well as model construction and application. Arterial-phase images from contrast-enhanced MRI scans were exported in Digital Imaging and Communications in Medicine format from the hospital’s Picture Archiving and Communication System. Image preprocessing was performed using the Python SimpleITK toolkit, including grayscale normalization [normalizing image grayscale values to the range (0, 255)] and isotropic resampling (unifying spatial resolution to 1 mm × 1 mm × 1 mm). These steps aimed to reduce the impact of differences in scanning equipment and parameters, ensuring consistency and reliability for subsequent analyses.

Figure 2 Schematic diagram of the methodological workflow. The workflow is divided into two parts: The first part is preprocessing and feature acquisition. First, through region of interest segmentation, masks are obtained from the original images of hepatocellular carcinoma patients; then, via habitat clustering, subregion clustering is performed using simple linear iterative clustering superpixels to obtain habitats; subsequently, feature extraction and feature selection are conducted in sequence. The second part is model construction and application: A combined model is built based on the intratumoral heterogeneity score, radiomics score, and clinical data; model evaluation is completed using the receiver operating characteristic curve, calibration curve, SHapley Additive exPlanation analysis, and decision curve; finally, the model is applied in clinical practice to distinguish between responders (treatment responders) and non-responders (non-treatment responders). ROI: Region of interest; HCC: Hepatocellular carcinoma; GMM: Gaussian mixture model; BIC: Bayesian information criterion; ROC: Receiver operating characteristic; SHAP: SHapley Additive exPlanation.

Preprocessed arterial-phase MRI images of lesions were imported into the open-source software 3DSlicer (version 5.7.0; https://www.slicer.org) for lesion segmentation. Two radiologists, with reference to the lesion’s T2-weighted images and post-contrast venous phase and delayed phase images, manually outlined the region of interest (ROI) slice by slice along the tumor boundary on axial arterial-phase enhanced images of the lesion. The outlined ROI was required to include intratumoral necrosis, hemorrhage, and the peritumoral capsule, but exclude peritumoral satellite nodules. Subsequently, 3DSlicer automatically fused the slice-by-slice outlined ROIs to form a three-dimensional volume of interest. When multiple HCCs were present in a patient’s images, the largest HCC was selected for outlining.

To ensure the stability and reproducibility of lesion segmentation results, the intraclass correlation coefficient (ICC) was used to evaluate segmentation consistency: 30 cases were randomly selected, and two radiologists independently performed ROI outlining to assess inter-rater consistency; meanwhile, rater 1 re-outlined the same 30 cases after 4 weeks to assess intra-rater consistency. The results showed that the ICC for inter-rater segmentation consistency was 0.92 [95% confidence interval (CI): 0.85-0.96], and the ICC for intra-rater segmentation consistency was 0.94 (95%CI: 0.88-0.97). According to ICC evaluation criteria, both values exceeded 0.85, indicating excellent consistency in ROI outlining in this study, with reliable segmentation results that meet the requirements for subsequent feature extraction and model construction.

Tumor subregion segmentation and feature extraction

To further characterize the spatial heterogeneity of intratumoral structure and signal, this study employed the simple linear iterative clustering algorithm to divide the entire tumor into multiple supervoxels, thereby generating objective intratumoral subregions. Subsequently, to capture potential radiological pattern differences among subregions, the Gaussian mixture model (GMM) was used to cluster these radiomic features extracted from subregions. The optimal number of clusters (ranging from 1 to 5) was determined via the Bayesian information criterion (BIC), which quantified the “ecological diversity” within each tumor, i.e., the number of subregions and the complexity of signal distribution. This segmentation and clustering process constituted a modeling framework for tumor ecological heterogeneity, providing a computational basis for the subsequent construction of the ITH index.

For feature extraction, the research team used PyRadiomics (v3.1) to extract a total of 1132 radiomic features from the entire tumor and each of its subregions, including histogram statistical features, morphological features, and texture features (such as gray-level co-occurrence matrix, gray-level run-length matrix, gray-level size zone matrix, gray-level dependence matrix, and neighborhood gray-tone difference matrix). These features strictly adhered to the definitions of the image biomarker standardization initiative. Through the above process, each patient’s tumor yielded a high-dimensional “ecological diversity feature vector”, which reflects the number of tumor subregions, intensity distribution, and structural complexity. This vector was ultimately converted into a continuous ITH index via the model. The index can be understood as a quantitative indicator of “tumor ecosystem complexity” based on MRI images, aiming to serve as a non-invasive surrogate for tumor biological heterogeneity.

Feature selection and model development

In the stage of feature selection and model development, this study designed a multi-level dimensionality reduction process to ensure the stability, independence, and predictive relevance of features in the training dataset. First, low-information features with variance below 75% were excluded, followed by the elimination of multicollinearity (r > 0.75) using the Pearson correlation coefficient. After preliminary screening at the feature level, t-tests or Mann-Whitney U tests were further conducted to compare the feature distribution between patients in the TACE response group and non-response group, retaining only variables with statistically significant differences (P < 0.05).

At the clinical variable level, univariate logistic regression was used to screen clinicopathological variables significantly associated with TACE response (e.g., ECOG performance status, ALP, BCLC staging). These variables were then incorporated into stepwise backward multivariate regression based on the Akaike information criterion to obtain the optimal variable combination. Combining the above features, seven model systems were constructed: Standalone clinical model, ITH model, traditional radiomics model, two-model combinations (ITH-radiomics, clinical + ITH, clinical + radiomics), and a three-model integrated combined model. The clinical model was implemented using multivariate logistic regression, the ITH and radiomics models using random forest based on scikit-learn, and the final combined model using multivariate logistic regression. The probability values output by the models were defined as the ITH index and radiomics score, respectively, with the optimal Youden’s index used as the threshold for classification.

Model performance was evaluated using the internal and external validation sets, including the area under the curve (AUC) with 95%CIs, Hosmer-Lemeshow calibration curves, and decision curve analysis (DCA). Multiple metrics such as sensitivity, specificity, accuracy, precision, and F1-score were also calculated. Results showed that the combined model achieved an AUC of 0.966 in the training set, with AUCs of 0.907 and 0.933 in the internal and external validation sets, respectively, significantly outperforming single-modal models. To enhance the interpretability of the model and verify the relative importance of each feature, we conducted SHapley Additive exPlanations (SHAP) analysis on the final combined model to further validate the incremental value of ITH-related features. The relevant results are presented in Supplementary Figure 1.

Statistical analysis

Continuous variables and categorical variables were presented as median (interquartile range) and n (%), respectively. Fisher’s exact test, Pearson’s χ² test, or Kruskal-Wallis test were used to compare the clinical characteristics of HCC patients after TACE among the training set, internal validation set, and external validation set. Univariate and multivariate logistic regression analyses were conducted to evaluate the association between different patient features (e.g., baseline clinical indicators, radiomic features) and early postoperative efficacy.

The AUC with 95%CI was used to assess the performance of the clinical model, traditional radiomics model, ITH model, and combined model in predicting the early efficacy of TACE for HCC. The Hosmer-Lemeshow goodness-of-fit test was applied to calibrate the agreement between the model-predicted probabilities and actual efficacy. DCA was performed to evaluate the clinical utility of the models. Sensitivity, specificity, accuracy, precision, and F1-score were calculated to quantify the clinical application value of the models. During model construction, variables with a P value < 0.05 in univariate analysis were included in the multivariate logistic regression analysis. All statistical analyses were conducted using Python (3.10.13), and two-tailed P < 0.05 was considered to indicate statistically significant differences.

RESULTS

Tumor microhabitat imaging and feature extraction

For tumor microhabitat imaging, this study used the simple linear iterative clustering algorithm to perform superpixel segmentation on the entire tumor, generating objective intratumoral subregions. Subsequently, the GMM was applied to cluster these subregions, and the BIC was used to determine the optimal number of clusters (ranging from 1 to 5). This successfully quantified the “ecological diversity” of each patient’s tumor (i.e., the number of subregions and the complexity of signal distribution), laying the foundation for subsequent ITH feature extraction and ITH index construction.

Supplementary Table 2 presents the feature selection process and corresponding results for each data source. For feature extraction, PyRadiomics (v3.1) was used, adhering to the definitions of the image biomarker standardization initiative. A total of 1132 radiomic features were extracted from the entire tumor and each subregion, including histogram statistical features, morphological features, and texture features (such as gray-level co-occurrence matrix, gray-level run-length matrix, gray-level size zone matrix, and gray-level dependence matrix). The final combined feature list is detailed in Supplementary Table 3. After screening, the following results were obtained. (1) ITH features: A total of 10 features, including morphological features (e.g., ITH_original_shape_Flatness), multi-scale texture features (e.g., ITH_log-sigma-2-0-mm-3D_glcm_Imc1, ITH_log-sigma-3-0-mm-3D_glrlm_LongRunEmphasis), and wavelet transform features (e.g., ITH_wavelet-LLH_glszm_GrayLevelNonUniformity, ITH_wavelet-LHL_firstorder_Median); (2) Traditional radiomics features: A total of 10 features, including morphological features (original_shape_Elongation, original_shape_Flatness), histogram features (original_firstorder_10Percentile, original_firstorder_90Percentile, original_firstorder_Entropy), and texture features (original_glcm_ClusterShade, original_glcm_DifferenceEntropy); and (3) Combined feature set: Clinical indicators (tumor size, ALP, number of tumors), all 10 ITH features, all 10 radiomics features were integrated to form the final feature matrix for model construction. This matrix covers tumor morphology, signal distribution, structural complexity, and clinical baseline information, enabling comprehensive characterization of tumor biological heterogeneity.

Logistic regression analysis and modeling of factors associated with treatment response

Table 2 presents the performance of various clinical factors in univariate and multivariate logistic regression. Results of univariate analysis showed that tumor size, ALP level, Child-Pugh classification, ECOG performance status, number of tumors, BCLC staging, and tumor thrombus were all significantly associated with the early efficacy of TACE in patients (all P < 0.05). Further multivariate logistic regression analysis indicated that tumor size [odds ratio (OR) = 0.37, 95%CI: 0.18-0.76, P = 0.0066], ALP level (OR = 0.28, 95%CI: 0.13-0.60, P = 0.0009), and number of tumors (OR = 0.39, 95%CI: 0.18-0.82, P = 0.0136) were independent influencing factors for the early efficacy of TACE. Other clinical factors, including tumor thrombus, were not included in the multivariate model due to P values > 0.05.

Table 2 Univariable and multivariable logistic regression analysis of clinical features.

Features	Univariate		Multivariate
	OR (95%CI)	P value	OR (95%CI)1	P value
Sex	0.886 (0.390-2.010)	0.771
Age	1.002 (0.976-1.029)	0.879
Tumor size	0.177 (0.088-0.355)	< 0.001	0.370 (0.180-0.760)	0.0066
Tumor number (multiple)	0.411 (0.212-0.805)	0.009	0.390 (0.180-0.820)	0.0136
Tumor thrombus	0.503 (0.266-0.971)	0.044
Liver cirrhosis	0.686 (0.341-1.379)	0.290
ALP (> 125 U/L)	0.120 (0.058-0.251)	< 0.001	0.280 (0.130-0.600)	< 0.001
AFP (> 400 ng/mL)	0.779 (0.365-1.659)	0.517
Hepatitis B history	0.642 (0.255-1.618)	0.347
ECOG PS	0.352 (0.162-0.741)	0.007
BCLC stage	0.690 (0.488-0.975)	0.035
AST (> 40 U/L)	0.692 (0.358-1.340)	0.275
ALT (> 50 U/L)	1.095 (0.542-2.211)	0.801
Total bilirubin (< 19 μmol/L)	1.382 (0.728-2.623)	0.323
Albumin (> 40 g/L)	0.844 (0.393-1.813)	0.663
Child-Pugh class	0.311 (0.151-0.667)	0.002

¹Empty cells indicate the variable was not included in the multivariate regression model or had no statistical significance in multivariate analysis.

OR: Odds ratio; CI: Confidence interval; ALP: Alkaline phosphatase; AFP: Alpha-fetoprotein; ECOG PS: Eastern Cooperative Oncology Group performance status; BCLC: Barcelona Clinic Liver Cancer; AST: Aspartate aminotransferase; ALT: Alanine aminotransferase.

Performance evaluation and comparison of predictive models

Results of the performance evaluation and comparison of predictive models showed that the predictive efficacy of each model in Table 3 varied across the training set, internal validation set, and external validation set. In the training set, the combined model achieved an AUC of 0.97 (95%CI: 0.93-0.99), with an accuracy of 0.90, sensitivity of 0.88, and specificity of 0.92, significantly outperforming the clinical model (AUC = 0.86, 95%CI: 0.78-0.93), radiomics model (AUC = 0.71, 95%CI: 0.60-0.82), and other pairwise combined models.

Table 3 Prediction performance of different models.

Models	Accuracy	AUC	95%CI	Sensitivity	Specificity	Precision1	F12	Cohort
Clinical model	0.79	0.86	0.78-0.93	0.77	0.82	0.84	0.80	Train
Radiomics model	0.69	0.71	0.60-0.82	0.73	0.63	0.71	0.72	Train
ITH model	0.85	0.89	0.81-0.96	0.79	0.92	0.93	0.85	Train
ITH + clinical model	0.90	0.95	0.90-0.99	0.90	0.90	0.92	0.91	Train
ITH + radiomics model	0.87	0.93	0.87-0.98	0.85	0.90	0.91	0.88	Train
Radiomics + clinical model	0.85	0.90	0.82-0.96	0.92	0.76	0.83	0.87	Train
Combined model	0.90	0.97	0.93-0.99	0.88	0.92	0.93	0.90	Train
Clinical model	0.80	0.78	0.54-0.98	0.75	0.89	0.90	0.82	Validation
Radiomics model	0.67	0.66	0.39-0.89	0.58	0.78	0.78	0.67	Validation
ITH model	0.87	0.90	0.74-1.00	0.75	1.00	1.00	0.86	Validation
ITH + clinical model	0.91	0.88	0.64-1.00	1.00	0.78	0.86	0.92	Validation
ITH + radiomics model	0.87	0.91	0.75-1.00	0.83	0.89	0.91	0.87	Validation
Radiomics + clinical model	0.76	0.68	0.42-0.92	0.75	0.78	0.82	0.78	Validation
Combined model	0.91	0.91	0.71-1.00	0.92	0.89	0.92	0.92	Validation
Clinical model	0.41	0.54	0.33-0.75	0.37	0.75	0.92	0.53	External validation
Radiomics model	0.87	0.78	0.57-0.96	0.90	0.63	0.95	0.93	External validation
ITH model	0.71	0.87	0.77-0.96	0.68	1.00	1.00	0.81	External validation
ITH + clinical model	0.71	0.89	0.78-0.97	0.68	1.00	1.00	0.81	External validation
ITH + radiomics model	0.86	0.94	0.88-0.99	0.84	1.00	1.00	0.91	External validation
Radiomics + clinical model	0.89	0.78	0.56-0.98	0.92	0.63	0.95	0.93	External validation
Combined model	0.83	0.93	0.87-0.99	0.81	1.00	1.00	0.89	External validation

¹Positive predictive value (proportion of true positives among predicted positives).

²Harmonic mean of precision and sensitivity, balancing both metrics.

ITH: Intratumoral heterogeneity; AUC: Area under the curve; CI: Confidence interval.

In the internal and external validation sets, the combined model still maintained good generalization ability: The AUC was 0.91 (95%CI: 0.71-1.00) in the internal validation set and 0.93 (95%CI: 0.87-0.99) in the external validation set, showing a significant advantage in overall performance. The varying performances of different models on the internal training set, internal validation set, and external validation set are visually presented in Figure 3. Further DCA (Figure 4) revealed that within a wide range of threshold probability intervals in the internal test set (Figure 4A), internal validation set (Figure 4B), and external validation set (Figure 4C), the net benefit of the combined model was significantly higher than that of single models (clinical model, radiomics model, ITH model) and other pairwise combined models, indicating good guiding value for clinically screening patients who may benefit from TACE.

Figure 3 Performance and comparison of the seven models in the internal training set, internal validation set, and external validation set. This figure presents the receiver operating characteristic curves of the seven models across three cohorts. A: Internal training set; B: Internal validation set; C: External validation set. The results show that the combined model, integrating clinical variables, traditional radiomics features, and the intratumoral heterogeneity index, exhibits the best performance. Its area under the curve values reach 0.97, 0.91, and 0.93 in the internal training set, internal validation set, and external validation set, respectively, which are the highest among all evaluated models. AUC: Area under the curve; CI: Confidence interval; ITH: Intratumoral heterogeneity.

Figure 4 Decision curve analysis of the seven models in the internal training set, internal validation set, and external validation set. This figure presents the decision curve analysis results of the seven models across three cohorts. A: Internal training set; B: Internal validation set; C: External validation set. As shown in the figure, the combined model, integrating clinical variables, traditional radiomics features, and the intratumoral heterogeneity index, achieves significantly higher net benefits than other single models or pairwise combined models within a wide range of threshold probability intervals in all three panels. This indicates that the combined model has good guiding value for clinical decision-making. ITH: Intratumoral heterogeneity.

Calibration curve analysis (Supplementary Figure 2) further verified the predictive consistency of the combined model: In the training set and internal validation set, the calibration curve of the combined model showed a high degree of fit with the ideal diagonal line, suggesting that the predicted probability of non-response to TACE had a small deviation from the actual outcome. In the external validation set, the curve still maintained a good fit, indicating that the consistency between the predicted values and actual efficacy of the model was high across cohorts from different institutions, further supporting its credibility for clinical application.

DeLong test results in Supplementary Table 4 confirmed that the combined model achieved significantly higher AUC values than the single clinical model, single radiomics model, and radiomics + clinical combined model across the training, internal validation, and external validation sets (all P < 0.05). This statistical evidence further validates the superior discriminative power of the integrated model that incorporates ITH features for predicting TACE response in HCC patients.

To further clarify the relative contribution of each feature to the predictive performance of the combined model, we conducted SHAP analysis on the final combined model (Supplementary Figure 1). Results from the SHAP summary plot showed that ITH-related features ranked among the core predictors, directly confirming the incremental value of ITH features in the model and explaining the core reason why the combined model outperforms single-modality models from the perspective of feature contribution.

DISCUSSION

TACE is a key clinical treatment for patients with HCC. However, HCC exhibits high heterogeneity, leading to significant variations in clinical responses to TACE among different patients and thus affecting clinical treatment decisions[18]. Preoperative prediction of TACE response in HCC patients is therefore crucial[19]. This study developed and validated an ITH model based on preoperative MRI images of HCC patients undergoing TACE; this habitat model demonstrated strong performance in predicting early post-TACE treatment response. When integrating HCC’s ITH features with clinical, conventional radiological, and traditional radiomics features, the combined model achieved the highest predictive performance. This finding highlights that ITH features play a non-negligible role in HCC treatment response and can effectively enhance the predictive accuracy of the combined model, thereby providing important support for clinical efficacy prediction and the development of personalized treatment strategies. SHAP analysis further validates this conclusion: The results show that ITH features rank high in their contribution to prediction outcomes, clarifying that ITH features are key incremental factors for boosting the model’s predictive performance.

Extensive research has been conducted on models predicting post-TACE efficacy in HCC. Early studies mostly relied on clinical models, such as scoring systems constructed using indicators like tumor size, number, and alpha-fetoprotein levels[9,20], but their predictive efficacy was limited. In recent years, traditional radiomics and deep learning methods have been widely used to build predictive models. For example, Kong et al[21] developed and validated a nomogram based on preoperative MRI radiomics, which effectively predicts tumor response to TACE in patients with intermediate-advanced HCC. Combining a radiomics score (Rad-score) with clinical indicators, the model achieved AUCs of 0.861 and 0.884 in the training and validation sets, respectively. Kuang et al[22] retrospectively collected MRI images and clinical data from 153 patients, constructing an MRI-based radiomics nomogram to predict short-term response to TACE in HCC patients with tumors smaller than 5 cm; the AUCs in the training and validation sets were 0.83 and 0.81, respectively. In terms of deep learning, Dong et al[23] applied six machine learning algorithms to build CT image-based models predicting early treatment response to initial cTACE for HCC, and found that the random forest model exhibited the best performance, with an AUC of 0.802. A literature report[24] indicated that models constructed using non-radiomics approaches had a median AUC of 0.73 (95%CI: 0.61-0.77), while radiomics models had a median AUC of 0.79 (95%CI: 0.75-0.82); the median AUCs of models built with different deep learning algorithms mostly clustered around 0.80.

Despite progress in previous studies, both traditional radiomics and deep learning models treat tumors as homogeneous entities when extracting features[25,26], ignoring ITH and limiting their application in HCC with diverse biological behaviors. Prior research has shown that tumor heterogeneity exerts a significant impact on treatment response[27,28], yet few studies have integrated ITH with traditional radiomics and clinical variables to construct models for predicting TACE response. This study is among the first to address this research gap: By segmenting tumor subregions and extracting features to characterize tumor heterogeneity, the combined model built with clinical factors and traditional radiomics demonstrated impressive predictive efficacy. The clinical-traditional radiomics model constructed in this study achieved AUCs of 0.896, 0.676, and 0.782 in the training, internal validation, and external validation sets, respectively, consistent with the efficacy of similar models reported in previous literature. After incorporating ITH into model construction, the clinical-traditional radiomics-ITH combined model showed further improved AUCs of 0.966, 0.907, and 0.933 in the three sets, outperforming most non-radiomics, radiomics, and deep learning models reported in previous studies[24]. Additionally, in DCA curves, the combined model exhibited higher net benefits than other models. This indicates that incorporating ITH as a predictive factor into model construction is of great value; the resulting combined model has clinical utility and provides new insights and important references for clinical treatment decisions. Meanwhile, the model’s strong performance in the external validation set (AUC = 0.933) confirms its good generalization ability. Overall, methodologically, this study combines radiomics with ecosystem modeling concepts, achieving MRI-based heterogeneity characterization through multi-level quantification of intratumoral subregions and statistical learning. Its innovation lies in the automatic definition of tumor “ecological complexity” via GMM clustering and BIC optimization, as well as the validation of predictive generalization in external validation and multi-modal data. This heterogeneity quantification model established from the perspective of “ecological diversity” not only enriches the radiomics framework for TACE efficacy prediction but also provides a scalable technical pathway for future tumor habitat analysis based on multi-modal imaging.

The clinical-traditional radiomics-ITH combined model constructed in this study offers important reference value for preoperative precise decision-making in HCC patients undergoing TACE. First of all, the multi-center study design with independent external validation from institution B stands as a major strength of our research. This design effectively addresses the prevalent single-center bias in radiomics studies, as the external cohort from a distinct clinical institution consistently validated the model’s robust generalization ability, achieving an impressive AUC of 0.93 in the external validation set. Such cross-institutional reproducibility not only enhances the credibility of our findings but also lays a solid foundation for its clinical translation across different healthcare settings. Currently, traditional methods for preoperative prediction of TACE efficacy mainly rely on imaging findings, laboratory test indicators, and patients’ systemic status[29]. However, due to significant heterogeneity in patients’ TACE responses, the predictive ability of these traditional indicators is limited[30,31]. Some patients meet TACE indications but experience early progression postoperatively due to high tumor heterogeneity and aggressive biological behavior; they not only fail to benefit but also bear treatment-related risks such as liver function damage and post-embolization syndrome[32], while causing inefficient consumption of medical resources. By integrating ITH, a key biological feature, the model enables preoperative stratification of patient benefits: For patients predicted to have high benefits, TACE can be clearly prioritized as first-line treatment, avoiding delays in tumor control due to empirical judgment biases; for patients predicted to have low benefits, TACE intensity can be reduced, with priority given to strengthening systemic therapy (e.g., replacing targeted drugs, increasing immunotherapy dosage) or adjusting local treatment methods (e.g., switching to hepatic arterial infusion chemotherapy, radiofrequency ablation)[33]. Furthermore, this model serves as a reference framework for combined therapeutic efficacy to distinguish the contributions of TACE and systemic therapy in the combined effect.

However, there are also some challenges in the clinical translation of this model. Firstly, differences in MRI with different field strengths and scanning parameters may affect feature consistency; to address this, we have provided a standardized scanning parameter template for this study (Supplementary Table 1), and expanded studies involving multiple field strengths and devices will be conducted in the future. Secondly, manual layer-by-layer ROI segmentation is time-consuming, and we will develop an automatic segmentation model subsequently. Finally, MRI examination is more expensive than CT, which may increase the burden on some patients.

This study has several limitations that require objective consideration. First, we acknowledge that ICC-based robustness filtering for radiomic and ITH features was not performed prior to variance and correlation filtering. Although the intra- and inter-rater consistency of ROI segmentation was excellent, ensuring the stability of extracted features against minor segmentation variations via ICC filtering is a standard practice in modern radiomics research. In future studies, we will supplement this step to further enhance the robustness and clinical translatability of the predictive model. Second, it adopts a retrospective cohort design, which inevitably introduces selection bias. Although external validation was performed to improve model reliability, prospective cohort studies are still needed to further validate model efficacy. Third, manual delineation of three-dimensional volumes of interest was used to extract global tumor features; while this maximizes the accuracy and completeness of ITH features, the process is time-consuming and labor-intensive, and susceptible to the operator’s subjective experience. Future research should explore deep learning-based automatic tumor segmentation algorithms to improve segmentation efficiency while reducing subjective errors, laying a technical foundation for clinical translation of the model. Finally, the study only quantified tumor ITH via habitat analysis and did not establish a deep link between features and biological essence, a common challenge in current radiomics research. On one hand, there is a lack of histopathological validation, with no analysis of correlations between imaging ITH features and pathological heterogeneity indicators using postoperative multi-site biopsy specimens. On the other hand, molecular-level data (e.g., circulating tumor DNA mutations, gene expression profiles) were not included, making it impossible to confirm whether the heterogeneity reflected by imaging features corresponds to specific molecular subtypes or treatment resistance mechanisms. A multi-dimensional validation system is needed in future studies to address this limitation.

CONCLUSION

In conclusion, the combined predictive model developed in this study, integrating clinical, conventional radiological, traditional radiomics, and ITH features, shows significant value in predicting early efficacy of initial TACE treatment for HCC patients. It is expected to provide quantitative references for preoperative clinical decisions regarding TACE in HCC patients, assist clinicians in developing personalized treatment plans for different patients, and thereby promote the implementation of precise HCC diagnosis and treatment models in clinical practice.