Noninvasive imaging analysis of vascular phenotypes improves prognostic stratification in primary liver cancer: a multi-cohort study.

Introduction
Primary liver cancers (PLC), principally hepatocellular carcinoma (HCC) and intrahepatic cholangiocarcinoma (ICC), pose a major global health burden due to rising incidence and marked heterogeneity in clinical outcomes1. Precise preoperative subtyping is imperative, given stark differential responses to curative-intent resection, loco-regional therapies, and systemic agents between HCC and ICC2. However, the diagnostic accuracy of the Liver Imaging Reporting and Data System (LI-RADS) remains suboptimal (60-75%), frequently failing to resolve tumors with overlapping imaging phenotypes, particularly combined hepatocellular-cholangiocarcinoma (CHC) and poorly differentiated subtypes3,4. This diagnostic ambiguity critically impedes recurrence risk stratification: HCC recurrence is driven by early intrahepatic dissemination via tumor-associated hyperpermeable sinusoids, whereas ICC recurrence manifests through late-stage lymphatic metastasis facilitated by fibrotic vascular stroma4,5. Although machine learning (ML) models exist for PLC diagnosis or prognosis, current approaches treat these tasks in isolation, neglecting shared pathobiological drivers.
Multi-task learning (MTL) offers a paradigm to overcome this limitation. By concurrently interrelated clinical objectives within a unified architecture, MTL exploits inherent task correlations to enhance data efficiency, mitigate over-fitting, and ultimately improved model generalizability and predictive performance6. Its efficacy hinges on identifying common biological determinants across tasks.
The tumor vascular microenvironment (TVME) represents one such pivotal orchestrator. HCC demonstrates disorganized, hyperpermeable sinusoidal vasculature driven by VEGF-mediated angiogenesis, which fosters hypoxia-induced microvascular invasion (MVI)7. Conversely, ICC exhibits dense fibrotic capillary networks regulated by FGFR/PDGF signaling, promoting stromal barrier formation that enables chemotherapy resistance and lymphatic dissemination8. Critically, the absence of noninvasive tools to decode TVME heterogeneity and its molecular regulators necessitates reliance on invasive biopsy or surgical resection for accurate risk assessment. While conventional radiomics has shown promise in prognostic modeling, its application to deconvolve the specific contributions of the TVME remains challenging. This is primarily because conventional approaches typically rely on agnostic, manually engineered features that describe general tumor texture but are anatomically non-specific. These features fail to distinguish the signal of the vascular network from other tissue components (e.g., cellular regions, necrosis), thereby limiting their capacity to specifically characterize the spatial topology and functional dynamics inherent to tumor vasculature9. Consequently, there is a pressing need for biology-aware feature extraction strategies that can directly target and quantify vascular properties. Recent innovations in Hessian matrix-based Frangi vessel filtering enable quantification of biologically relevant vascular signatures directly from routine CT scans, including vessel density, caliber distribution, tortuosity, and spatial anisotropy10–12. These quantitative vascular (QV) features bridge radiological phenotypes with angiogenic and stromal remodeling biology, enhancing both biological interpretability and potential for generalizability across tasks.
We therefore propose a Multi-Task Vessel computational Network (MTV-Net) integrating PLC reclassification with recurrence prediction using QV features. Our approach comprises three core components: (1) Vessel Phenotyping: Leveraging Hessian-Frangi filtering to quantify spatiotemporal TVME characteristics (vessel richness, enhancement kinetics) from contrast-enhanced CT, thereby capturing subtype-specific vascular signatures; (2) MTL Architecture: Employing established RMTL pipeline with shared feature embedding, facilitating joint optimization of clinically interdependent tasks to enhance robustness and biological fidelity; (3) Radiogenomics Validation: Correlating QV features and bulk RNA sequencing data, using pathways enrichment analyses to identify transcriptomic programs underpinning TVME imaging phenotypes.

Results

The characteristics of enrolled cohorts
The overall workflow of this study is illustrated in Fig. 1, and implemented on six clinical cohorts (n = 1020) across HCC, ICC, and CHC. Demographic and clinic-pathologic characteristics of patients included in diagnostic and prognostic analyses are summarized in Table 1, presented separately for each study cohort. Patient characteristics in the radiogenomics cohort, which was obtained from TCIA database, are detailed in Table S1.

Derivation and validation of the two vessel-imaging biomarkers
In the training cohort, the MTV-Net was optimized through 1500 iterative epochs and 10-fold cross-validation, with a regularization parameter (λ1) set at 0.001 (Fig. S1). The feature coefficients contributing to TAVSPHE and TAVSRE are detailed in Table S2. TAVSPHE demonstrated notable accuracy for differentiating HCC from ICC, yielding AUCs of 0.949, 0.879, and 0.821 in three independent validation cohorts, respectively (Fig. 2A). Notably, a significantly higher TAVSPHE was observed in HCC than that in ICC within each cohort (Fig. 2B). Using Youden’s index to determine optimal cut-off values, patients were stratified into “HCC-like” (scores above the threshold) and “ICC-like” (scores below the threshold). The confusion matrix illustrated strong agreement between model predictions and pathological diagnoses (Fig. 2C), with an overall classification accuracy of 82.0%, sensitivity of 82.1%, and specificity of 82.1%. Class-specific performance metrics are summarized in Table S3.
To assess the stability of TAVSPHE, three independent training-validation datasets combinations were randomly generated from HCC/ICC samples at NFHSMU center, maintaining an 8:2 training-to-validation ratio. Random seeds were predefined to ensure reproducibility of the sampling process. Retrained models were subsequently evaluated against the original model, and no significant differences in AUC values were observed across validation datasets (DeLong’s test, all P > 0.05; Fig. S2), indicating the robust generalizability of TAVSPHE.
TAVSRE demonstrated promising performance in predicting recurrence, with AUC values of 0.667, 0.644, and 0.710 in three validation cohorts (Fig. 2D). Significantly higher TAVSRE were observed in patients with recurrence compared to those without (all P < 0.05; Fig. 2E). Using the ‘survminer’ R package, cohort-specific cut-off values were determined to stratify patients into “High-risk” and “Low-risk” groups. High-risk patients exhibited a significantly elevated recurrence rate (all P < 0.001; Fig. 2F), underscoring the clinical utility of TAVSRE in predicting recurrence across PLC subtypes.

The two vessel-imaging biomarkers accurately predict prognosis
Patients with PLC were stratified into distinct subgroups based on the two vessel-imaging biomarkers. In three independent validation cohorts, both TAVSPHE and TAVSRE were significantly associated with DFS (all P < 0.05), with similar associations observed for OS (Fig. 3A–C). Adjusting for clinic-pathologic risk factors, multivariable Cox regression analyses confirmed that TAVSPHE and TAVSRE remained independent prognostic factors for DFS and OS in each cohort (all P < 0.05; Fig. S3).
The combined use of TAVSPHE and TAVSRE generated four imaging subtypes with distinct prognosis. Kaplan-Meier curves demonstrated that the integrated biomarker was an independent predictor of DFS (P < 0.0001) and OS (P < 0.0001) across all validation cohorts (Fig. 4A, B). In the Meta-cohort comprising 527 patients from three validation datasets, imaging subtype 1 exhibited the highest 3-year DFS (61.8%) and OS (83.1%) rates, whereas subtype 4 had the poorest outcomes (13.9% and 9.1%, respectively; all P < 0.0001; Figure. S4). Conditional stratification was performed through stepwise application of the two imaging biomarkers (TAVSPHE → TAVSRE), while integrating with clinical staging systems (BCLC for “HCC-like” and TNM for “ICC-like”; Fig. 4C). The combined index (TAVSRE/Clinical Stage) effectively stratify patients into four prognostic groups, with significant survival differences within each category (BCLC0-A High vs. Low, P < 0.001, HR = 3.19, 95% CI: 2.18–4.67; BCLCB High vs. Low, P = 0.012, HR = 2.17, 95% CI: 1.19–3.98; TNMI-II High vs. Low, P = 0.0003, HR = 2.50, 95% CI: 1.51–4.15; TNMIII High vs. Low, P = 0.016, HR = 2.21, 95% CI: 1.16–4.21; Fig. 4D and Fig. S5A). Compared to TAVSRE or Clinical Stage alone, the integrated approach improved DFS prediction accuracy, as evidenced by higher C-index values (“HCC-like”: 0.693 vs. 0.649 vs. 0.596; “ICC-like”: 0.651 vs. 0.616 vs. 0.615), AUCs for predicting 2-, 3-, 5-year DFS (“HCC-like”: 0.710-0.721; “ICC-like”: 0.773-0.817; Fig. S5B), highlighting the complementary value of vessel-imaging biomarkers and conventional staging systems.
Notably, some discrepancies in prognostic stratification were observed between imaging subtype 2 and 3, underscoring the need for enhanced risk stratification strategies. To address this, nomograms integrating the two imaging biomarkers and clinic-pathologic variables were developed using Cox proportional hazards regression for predicting DFS and OS in the Meta cohort (Fig. S6). Encouragingly, the integrated model effectively stratified patients into new “High-risk” (44.4%) and “Low-risk” (55.6%) groups, with 3-year DFS rates of 14.3% and 56.2%, and 5-year DFS rates of 12.5% and 46.0%, respectively (all P < 0.0001; Fig. 4E, F). Compared to TAVSRE or clinic model alone (Fig. 4G–J), the integrated approach exhibited superior predictive performance, as indicated by higher AUCs for 2-, 3-, 5-year DFS (0.765–0.815) and OS (0.805–0.821; Fig. S7). The C-index for the integrated model was significantly higher than individual model (DFS: 0.731 vs. 0.634–0.648; OS: 0.789 vs. 0.661–0.674). Calibration curves showed favorable agreement between predicted and observed survival probabilities (Fig. S8A, B). Quantitatively, the integrated model achieved a NRI of 0.189 (95% CI: 0.094–0.285, P < 0.001) and IDI of 0.090 (95% CI: 0.065–0.114, P < 0.001) compared to the clinic model (Fig. S8C), with similar improvements observed against the TAVSRE (NRI: 0.107, 95% CI: 0.045–0.169, P < 0.001; IDI: 0.067, 95% CI: 0.047–0.088, P < 0.001).

The integrated model improves risk stratification in CHC
CHC presents with histological features of both hepatocellular and biliary origin, posing diagnostic and therapeutic challenges. Given that the TAVSPHE and TAVSRE are fundamentally built on vascular patterns in HCC and ICC, CHC serves as an ideal independent validation cohort for assessing the generalizability of these biomarkers. Using MTV-Net, patients initially diagnosed with CHC were reclassified into “HCC-like” or “ICC-like” categories, enabling subsequent risk stratification (Fig. 5A). Representative CHC cases, including their contrast-enhanced CT images, tumor segmentation mask, predicted phenotype, and clinic follow-up data, are illustrated in Fig. 5B.
Historically, ICC exhibits greater biological aggressiveness and higher postoperative recurrence rates compared to HCC. Consistent with this, CHC patients reclassified as “ICC-like” by TAVSPHE demonstrated a recurrence rate of 69.7%, comparable to that of “pure” (pathologically confirmed) ICC (Fig. 5C). These patients also experienced significantly shorter DFS and OS compared to those reclassified as “HCC-like” (Fig. 5D). The distribution of TAVSRE differed significantly between TAVSPHE-predicted subtypes, highlighting the interdependence of these biomarkers (Fig. 5E). TAVSRE exhibited strong predictive performance for postoperative recurrence in CHC, with AUC values ranging from 0.722 to 0.786 for 2-, 3-, 5-year DFS. Notably, recurrence-positive patients had significantly higher TAVSRE than recurrence-free counterparts (Fig. S9). Patients stratified as “High-risk” by TAVSRE consistently exhibited elevated recurrence rate and inferior survival outcomes (Fig. 5F and Fig. S9), reinforcing the association between higher TAVSRE and poor prognosis in PLC. Multivariable Cox regression analysis confirmed that both TAVSPHE and TAVSRE were independent prognostic factors for DFS and OS in CHC (Fig. S10). Similar to HCC and ICC cohorts, the combined use of TAVSPHE and TAVSRE reclassified CHC patients into four imaging subtypes with distinct survival trajectories (Fig. 5G).
To optimize prognostic prediction, nomograms integrating the two imaging biomarkers and clinic-pathologic variables were developed for CHC (Fig. S11). The integrated model achieved high predictive accuracy, with AUC values of 0.856-0.939 for DFS and 0.854-0.949 for OS at 2-, 3-, 5-year time points (Fig. S12). Compared to the vessel model (TAVSRE/TAVSPHE), individual biomarkers (TAVSRE or TAVSPHE) and the clinic model, the integrated approach significantly improved survival prediction, as evidenced by higher C-index values (DFS: 0.785 vs. 0.648-0.710; OS: 0.823 vs. 0.658-0.746). For predicting postoperative recurrence, the integrated model outperformed single biomarker and clinic model, with AUCs of 0.839 vs. 0.683–0.736, respectively (DeLong’s test, all P < 0.05; Fig. 5H, I), underscoring the incremental value of combined vessel-imaging and clinical data in CHC management.

Model interpretability and biologic basis
We conducted a comprehensive assessment of the relative contribution of QV features to “ICC-like” and “High-risk” predictions using SHAP values in the Meta cohort. The top-10 influential variables, as determined by cumulative SHAP values, revealed distinct temporal patterns: seven of these features, primarily derived from the portal-venous phase, were associated with tumor phenotyping, while six arterial-phase features were pivotal for recurrence risk quantification (Fig. 6A, B). Subgroup analyses further demonstrated significant divergence between TAVSPHE subtypes: portal-venous phase features were strongly correlated with recurrence in “HCC-like” cases, whereas arterial-phase features dominated risk assessment in “ICC-like” cases (Fig. 6C, D). These findings suggest that TAVSPHE and TAVSRE effectively capture subtype-specific vascular heterogeneity, aligning with established radiological principles linking multiphasic CT enhancement patterns to tumor phenotypes and prognoses.
Radiogenomics analysis in the TCIA cohort revealed pathway-level associations underlying the prognostic significance of TAVSPHE and TAVSRE. GSEA demonstrated significant up-regulation of tumor aggressiveness-related pathways, including focal adhesion and ECM receptor interaction signaling, in both “ICC-like” and “High-risk” subgroups (Fig. 6E, F). These findings mechanistically explain the adverse prognosis associated with these phenotypes and validate the shared biological basis of the MTL framework. WGCNA identified a key gene module consisting of 683 genes (Table S4), with enrichment for angiogenesis-related pathways. Notably, this module was strongly associated with the formation of imaging subtype 4 (ICC/High), further supporting the role of aberrant angiogenesis in driving poor clinical outcomes (Fig. 6G, H).

Discussion
This study demonstrates that quantifying TVME heterogeneity via multiphasic CT imaging could refine PLC subtyping and serve as a robust prognostic biomarker across PLC subtypes. Leveraging interpretable QV features that systematically encode radiologists’ expertise in tumor vascular spatial distribution, temporal dynamics, and architectural complexity13, we developed two biologically grounded imaging biomarkers: TAVSPHE and TAVSRE. These biomarkers directly correspond to distinct vascular phenotypes reflecting tumor angiogenesis and microenvironmental remodeling14. Notably, the demonstrated cross-subtype generalizability highlights the clinical relevance of this approach, particularly as a practical advancement for risk stratification in resource-limited settings.
Translating artificial intelligence into clinical practice necessitates both predictive accuracy and robust interpretability15. Although deep learning (DL) has demonstrated exceptional performance in specific oncological tasks, such as tumor segmentation and classification, its inherent opacity as a “black-box” model hinders clinical adoption, especially in critical decision-making scenarios16. Conversely, our QV feature-based MTL framework directly addresses this challenge by establishing explicit links between imaging biomarkers and underlying biological mechanisms. This interpretability is fundamental, as it allows clinicians to understand the vascular features driving each prediction, fostering trust and paving the way for clinical integration. The primary contribution of this work is not merely a predictive model, but the establishment of a novel computational framework that simultaneously addresses two distinct but clinically paramount tasks—histological subtyping and prognosis stratification from routine CT scans.
Our framework leverages the inherent advantages of contrast-enhanced CT for volumetric vascular analysis. Moving beyond the two-dimensional, qualitative assessment of Doppler ultrasound, CT enables the quantification of three-dimensional vascular architecture and perfusion dynamics17. The MTV-Net translates these advantages into a non-invasive scoring system, providing a reliable alternative to invasive biopsies for tumor classification and risk assessment. This is particularly significant given the inherent limitations of biopsy, including sampling error and inability to capture the spatial heterogeneity of the entire tumor mass18. By providing a holistic, volumetric characterization of the tumor vasculature, our method mitigates the risk of sampling bias and offers a more comprehensive view of the tumor’s vascular phenotype, which is a recognized key determinant of aggressiveness.
The clinical management of PLC has traditionally been guided by the cell-of-origin classification, treating HCC and ICC as fundamentally distinct19. Mounting evidence now indicates that such rigid categorizations inadequately capture the biological plasticity and heterogeneity within and across these subtypes5,20. Our MTV-Net framework introduces a paradigm shift towards phenotype-driven reclassification. By learning subtype-specific vascular signatures (via TAVSPHE) and shared prognostic drivers of recurrence (via TAVSRE), our model move beyond histological labels to group patients based on their in vivo vascular phenotype. The resulting four-group classification system, integrating both subtyping and risk, demonstrates superior prognostic accuracy across all PLC subtypes. This suggests that vascular phenotypes, as captured non-invasively by our model, may reflect critical underlying biology more accurately than rigid histological categories, offering a more nuanced basis for prognostic stratification and potentially informing future therapeutic strategies.
A pivotal finding of our study, enabled by radiogenomic analysis, is the elucidation of the shared and distinct biological pathways that underpin our TVME-informed biomarkers. The enrichment of the ECM remodeling signaling pathway emerged as a unifying biological theme for both TAVSPHE and TAVSRE. This finding underscores that the ECM is far more than a passive structural scaffold21; it is a dynamically active component of the TVME that critically regulates both vascular function and tumor cell behavior22,23. We hypothesize that aberrant ECM remodeling creates a pathological mechanical and biochemical microenvironment. This dysregulated milieu alters integrin-mediated mechanical cues and growth factor signaling networks, which in turn govern two key, interrelated cancer hallmarks. First, it disrupts the normal processes of vascular morphogenesis, leading to the formation of the abnormal, tortuous vascular architectures that are captured by the phenotype-similarity biomarker TAVSPHE. Second, it activates pro-invasive signaling cascades in tumor cells, fostering the migratory and metastatic potential that is reflected in the recurrence-risk biomarker TAVSRE. Thus, our two imaging biomarkers, though designed for different clinical tasks, can be viewed as non-invasive readouts of the same overarching pathophysiological process of stromal dysregulation.
Delving deeper into the high-risk imaging phenotypes (“ICC-like” and TAVSRE-High), we found their radiological characteristics to be underpinned by a state of hyperactivation in angiogenesis-related pathways. To move beyond a simple correlation, we propose a mechanistic model centered on the concept of ineffective angiogenesis24. While pro-angiogenic signals are highly active, the resulting neovasculature is structurally and functionally incompetent—characterized by immaturity, high permeability, chaotic architecture, and poor perfusion. This model compellingly explains the well-established clinical paradox in oncology where high angiogenic signaling is frequently associated with poor patient prognosis. The Frangi-filter based vascular features extracted by our MTV-Net framework are exquisitely sensitive to vessel morphology. We posit that these features are, in fact, capturing the architectural chaos and dysfunction inherent to ineffective angiogenesis. The “ICC-like” phenotype, with its presumably more disorganized and less functional vasculature, mirrors the known aggressive biology of ICC, which is often associated with a dense, desmoplastic stroma that hampers proper vascularization and promotes treatment resistance25,26.
Several limitations inherent to this retrospective study design warrant acknowledgment. Although multi-center data minimized selection bias, prospective validation in larger, independent cohorts remains essential to confirm clinical utility. Furthermore, while TAVSPHE and TAVSRE demonstrate preliminary yet mechanistically plausible associations with ECM remodeling and angiogenesis, robust functional validation is imperative. This could be achieved through experimental approaches like endothelial lineage tracing or computational hemodynamic modeling. Critical next steps include prospective clinical trials evaluating the predictive and prognostic value of these biomarkers specifically in patients receiving targeted therapy (e.g., VEGF/FGFR inhibitors). Additionally, integrating MTV-Net outputs with dynamic contrast-enhanced MRI and liquid biopsy data holds significant promise for developing comprehensive, multi-parametric diagnostic and monitoring strategies. Despite the promising results of this study, scanner heterogeneity remains a challenge for clinical generalizability. To mitigate this technical variability, advanced harmonization techniques could be implemented in future work. For instance, statistical methods such as ComBat harmonization could be employed to explicitly model and adjust for scanner-specific effects in the feature space, potentially further improving model stability. More sophisticated DL approaches, such as image-to-image translation networks (e.g., CycleGAN), offer the potential to harmonize raw image data from different sources into a standardized style prior to analysis, thereby explicitly decoupling scanner-related variance from biologically relevant signal.
Rather than viewing DL and QV feature-based approaches as mutually exclusive entities, we advocate for a synergistic integration paradigm. DL excels at uncovering complex, high-order imaging features, such as texture patterns that imperceptible to human observers. Conversely, QV features offers mechanistic interpretability by linking imaging phenotypes to established biological knowledge; for instance, alterations in vessel caliber can reflect the response to anti-angiogenic therapies. The proposed integration strategy encompasses two key directions. First, QV features can be incorporated into DL architectures as regulatory elements, such as through the implementation of vascular-attention mechanism. These mechanisms enable the model to prioritize regions of interest based on known vascular pathophysiology, thereby enhancing feature extraction specificity. Second, QV features serve as validation framework for DL-derived biomarkers. By cross-referencing DL-identified imaging signatures against established pathophysiological pathways related to tumor angiogenesis and microenvironmental remodeling, this approach facilitates the transition from data-driven correlative models to biologically interpretable causative systems.
In summary, MTV-Net represents the first-of-its-kind TVME-driven imaging framework capable of simultaneously enabling dual cancer subtyping and prognostic prediction in PLC. By quantitatively translating intratumoral vascular heterogeneity into clinically actionable biomarkers, this approach significantly enhances preoperative risk stratification, offering particular utility for histologically ambiguous cases like CHC. Ultimately, MTV-Net establishes a foundation for personalized therapeutic strategies targeting TVME vulnerabilities.

Methods

Study population
The five in-house cohorts comprised 978 patients underwent hepatectomy at three independent medical centers from June 2013 to July 2022, which were approved by the institutional review board (IRB) of each participating center (NFEC-2022-119; 2023-021; 2025KS-YX-18-01). All procedures performed in this study were accordance with the Declaration of Helsinki and its amendments. Informed consent was waived for the utilization of anonymous data in this retrospective analysis. Overall cohorts were divided into discovery and validation datasets, as depicted in Fig. S13. Specifically, the training cohort (273 HCC, 112 ICC), internal validation cohort (68 HCC, 28 ICC) and expandability validation cohort (66 CHC) were retrieved from the Nanfang Hospital of Southern Medical University (NFHSMU, Guangzhou, China), the external validation-1 cohort (133 HCC, 56 ICC) was recruited from the Second Affiliated Hospital of University of South China (SAHUSC, Hengyang, China), and the external validation-2 cohort (167 HCC, 75 ICC) was recruited from the First Affiliated Hospital of University of South China (FAHUSC, Hengyang, China). The radiogenomics cohort with paired CT scans and Bulk-seq data (42 HCC) is publicly accessible through The Cancer Imaging Archive (TCIA) database (https://cancerimagingarchive.net). Follow-up data of these five in-house cohorts were collected from the initiation of surgery until 31 October 2023. Preoperative contrast-enhanced CT scans as well as complete clinic-pathological information were accessible for all enrolled patients.
The inclusion criteria were as follows: (1) being aged 18 years or older, (2) with a confirmed pathological diagnosis of liver tumor either through surgery or biopsy, (3) underwent a liver CT scan at the initial diagnosis, and (4) available clinical and follow-up data can be retrieved. Exclusion criteria included: (1) CT scans lacking complete images of the entire liver necessary for the reconstruction of three-dimensional lesion; (2) CT scans exhibiting suboptimal image quality that may compromise accurate analysis, including (i) the presence of artifacts, such as motion artifacts resulting from patient movement during scanning and metal artifacts arising from implants within the patient’s body and (ii) inadequate contrast defined by insufficient differentiation in gray levels between various tissues; (3) cases where liver tumors underwent treatments prior to CT scan; and (4) follow-up of less than three months.
This diagnostic/prognostic study followed the Standards for Reporting of Diagnostic Accuracy (STARD) checklist27, and the Transparent Reporting of a multivariable prediction model for Individual Prognosis Or Diagnosis (TRIPOD) checklist28.

CT data acquisition and annotation
Contrast-enhanced CT examinations followed a standardized abdominal imaging protocol to ensure consistency across diverse scanners, with arterial and portal-venous phase acquisitions timed at 30 s and 60 s post-contrast injection, respectively. All images, initially reconstructed at 1.5-5 mm slice thicknesses, underwent isotropic resampling (1 × 1 × 1 mm3 voxels) and intensity normalization to a [-75, 175 HU] soft-tissue window for enhanced vascular-tumor contrast.
Tumor segmentation was executed through a hierarchical workflow as previously described29. In brief, a primary 3D U-Net trained on nnU-Net framework30, performed automated whole-liver segmentation on arterial-phase scans, followed by a secondary 3D U-Net that precisely delineated the region of interest (ROI) of index tumor within the liver mask. To mitigate respiratory motion artifacts, non-rigid registration via ‘Elastix’ toolbox (https://github.com/SuperElastix/ElastixModelZoo) aligned arterial-phase tumor regions to portal-venous phase scans, constrained to the liver volume for spatial coherence. This end-to-end protocol ensured robust, reproducible tumor localization and multiphase vessel features extraction across heterogeneous imaging protocols.

Quantitative vascular features extraction
Initially, the segmented tumor ROI was preprocessed via a multi-scale Hessian-Frangi filter31, enhancing intratumoral vascular structures, as illustrated in Fig. S14. The filter selectively amplifies tubular structures (e.g., blood vessel) by computing second-order intensity derivatives (Hessian matrix eigenvalues), enabling robust detection of vascular networks across 2D axial slices10. From the filtered images, 13 first-order histogram features were then extracted using the ‘PREDICT’ open-source package (https://github.com/Svdvoort/PREDICTFastr) in Python (Version 3.7.1). These features, calculated from intensity distributions partitioned into 50-bin histograms, included: (1) Minimum; (2) Maximum; (3) Range; (4) Interquartile range; (5) Standard deviation; (6) Skewness; (7) Kurtosis; (8) Peak value; (9) Peak position; (10) Energy; (11) Entropy; (12) Mean; and (13) Median. To capture spatial heterogeneity, features were computed across three concentric tumor subregions: (1) Full tumor ROI: Entire index tumor (largest lesion in multifocal cases), (2) Inner region: Central zone (radius ≤ 5 mm from tumor centroid), and (3) Outer region: Peripheral annular zone (radius > 5 mm). This process was independently applied to arterial and portal-venous phase CT images, yielding 39 features per phase (13 features × three regions). The final QV features set comprised 78 parameters per patient, collectively encoding intratumoral vasculature richness, caliber distribution, and spatiotemporal dynamics.

Multi-task learning-based MTV-Net development
MTV-Net was implemented using ‘RMTL’ package in R environment32, integrating two clinically aligned tasks: subtype identification and recurrence prediction. Prior to analysis, all QV features underwent Z-score standardization to ensure scale invariance. RMTL’s automated pipeline was configured to combine five classification and five regression algorithms, unified via cross-task regularization. The model aimed to optimizes the following objective function:In Eq. 1, is task-specific loss (logistic loss for binary classification), is cross-task regularization ( regularization was selected in this study), enforcing shared sparsity patterns across tasks to identify QV features jointly predictive of cancer subtype and survival status, and is Frobenius norm penalty for generalization. and are sets of predictor matrices and corresponding responses of all t tasks, respectively, while each task i contains subjects and p predictors. is coefficient matrix, where is the th column of refers to the coefficient vector of task . Knowledge transfer among tasks is achieved via that jointly modulates MTL model according to the specific prior structure. The and are positive regularization parameters. aims to control the effect of cross-task regularization and could be estimated by cross-validation procedure, whereas is set by investigators in advance ( was set as zero in this study). The model was trained over 1500 iterations with 10-fold cross-validation on the training cohort. The predicted probabilities for cancer subtyping (Tumor-Associated Vessel Score for Phenotyping PLC, referred to as TAVSPHE hereafter) and probabilities for cancer recurrence (Tumor-Associated Vessel Score for Recurrence Estimation, referred to as TAVSRE hereafter) were calculated. Final MTV-Net coefficients were fixed and applied to internal and external validation cohorts to compute TAVSPHE and TAVSRE without retraining.

Evaluation of the model accuracy for PLC subtyping
The optimal cut-off value of TAVSPHE was determined using Youden’s index across validation cohorts, a method that maximizes the sum of sensitivity and specificity for binary classification. Diagnostic performance of TAVSPHE in differentiating HCC from ICC was evaluated using receiver operating characteristic (ROC) curve analysis. Evaluation metrics included the area under the ROC curve (AUC), diagnostic accuracy, sensitivity, specificity, and F1-score. 95% confidence intervals (CI) for all performance metrics were computed via bootstrapping to account for sampling variability. Confusion matrix was constructed to elucidate divergences between TAVSPHE predictions (‘HCC-like’ vs. ‘ICC-like’) and pathological diagnoses.

Evaluation of the model accuracy for prognosis prediction
The optimal cut-off value of TAVSRE was determined across validation cohorts using the surv_cutpoint function in the ‘survminer’ R package, which maximizes the separation of disease-free survival (DFS) distributions. This approach involved calculating cohort-specific thresholds by correlating TAVSRE with DFS outcomes, stratifying patients into high- and low-risk subgroups to mitigate inter-cohort computational batch effects. The same threshold was subsequently applied for overall survival (OS) analysis to ensure consistency. DFS was defined as the interval from surgery to the first documented tumor relapse or death, while OS was defined as the time from surgery to all-cause mortality. Prognostic performance of TAVSRE was evaluated using discrimination and calibration metrics for DFS and OS endpoints. Calibration curves were constructed to assess agreement between model-predicted survival probabilities and observed event outcomes, with bootstrapping employed to validate consistency across resampled datasets. Additionally, associations between TAVSPHE-predicted phenotypes (‘HCC-like’ vs. ‘ICC-like’) and clinical outcomes were evaluated using proportional hazards models, adjusting for TAVSRE and established clinicopathological covariates.

Evaluation of the integrated model accuracy for risk stratification
The TAVSPHE and TAVSRE were integrated into a “hybrid” imaging biomarker framework, generating four combinatorial imaging subtypes by cross-classifying TAVSPHE-predicted subtypes (HCC/ICC) and TAVSRE risk subgroups (High/Low): 1-HCC/Low, 2-HCC/High, 3-ICC/Low, and 4-ICC/High. Survival differences across these subtypes were evaluated using Kaplan-Meier estimates and log-rank tests. To simulate clinical decision-making workflows, conditional stratification was implemented by aligning scoring systems with TAVSPHE-predicted subtypes: the BCLC stage for “HCC-like” cases and the TNM stage for “ICC-like” cases. Two experienced radiologists independently assigned baseline staging, with discrepancies resolved through consensus review by a third senior radiologist. Clinicopathological variables, including age, gender, tumor size, tumor number, differentiation grade, and lymphovascular invasion, were reviewed and collected. Integrated models combining the two vessel-imaging biomarkers and clinical predictors were developed using Cox proportional hazards regression for survival outcomes and logistic regression for binary endpoints. Model performance was evaluated via prediction error curves, with calibration assessed using Hosmer-Lemeshow tests. Relative improvements in predictive accuracy were evaluated through decision curve analysis (DCA) and quantified using net reclassification improvement (NRI) and integrated discrimination improvement (IDI), accounting for reclassification of risk categories.

Model interpretability and transcriptome analysis
The Shapley Additive exPlanations (SHAP) framework was applied to quantify the relative importance of QV features in subtype identification and recurrence prediction, using an XGBoost model implemented via the ‘shapviz’ R package33. This approach decomposes model predictions into feature-level contributions, enabling identification of critical vascular metrics. Gene Set Enrichment Analysis (GSEA) was performed to identify pathways differentially activated between TAVSPHE-predicted subtypes and TAVSRE risk subgroups, with significance defined by a false discovery rate (FDR) less than 0.0534. Weighted gene co-expression network analysis (WGCNA) was employed to construct co-expressed gene modules correlated with TAVSPHE and TAVSRE phenotypes35. Following this, the module most strongly associated with the poorest-prognosis imaging subtype was subjected to functional annotation using ClueGo, a Cytoscape software (version 3.7.2) plug-in36. This tool aggregated non-redundant Gene Ontology (GO) terms into functionally grouped networks, visualizing enriched biological process and molecular pathways underlying high-risk imaging phenotypes.

Statistical analysis
Continuous variables are presented as mean value along with their corresponding standard deviation (SD) or median with interquartile (IQR). Categorical variables are reported as counts (%). For comparisons of subgroups, statistical significance was estimated by Mann-Whitney U test (also known as Wilcoxon rank-sum test) or two-sided Fisher exact test, as appropriate. Survival curves for subgroups in each cohort were generated using Kaplan-Meier method, and log-rank (Mantel-Cox) test was employed to ascertain statistical significance, along with corresponding hazard ratio (HR). The presentation of results from multivariate Cox regression analysis was facilitated via ‘forestplot’ R package. Visualization of Calibri curves was achieved through the implementation of ‘rms’ R package. NRI and IDI were computed using ‘PredictABEL’ R package. The ‘pROC’ and ‘timeROC’ R packages were utilized for plotting ROC curves. DeLong’s test was applied for comparison of AUCs. Statistical analyses were carried out using SPSS (version 22.0, Chicago, IL, USA) and R (version 4.5.0). P value less than 0.05 was considered as statistically significant in two-sided analyses.

Supplementary information