This retrospective study was approved by the institutional review boards of the two hospitals (2023(145)), and the need for informed consent was waived. A total of 499 patients with hepatocellular carcinoma (HCC) who underwent treatment with DEB-TACE at the Second Affiliated Hospital of Chongqing Medical University, and 235 patients from the First Affiliated Hospital of Army Medical University between January 2019 and February 2023. The inclusion criteria were as follows: 1) DEB-TACE as first-line treatment except the cases that previous treatment was 1 month ago and the target lesions treated were different from the current DEB-TAC; 2) Enhanced MRI within 4 weeks before DEB-TACE; 3) Enhanced MRI or enhanced CT within 3 months after DEB-TACE; 4) Nodular or Massive HCC. First-line treatment was defined as the initial treatment administered to a patient who had not received any prior therapy at the time of their HCC diagnosis. The exclusion criteria were as follows: 1) preoperative use of other examination methods or lack of preoperative imaging; 2) lack of postoperative imaging; 3) Diffuse HCC or lesions with a diameter <10 mm; and 4) poor image quality or lack of clinical data. Ultimately, 114 HCC patients from center 1 served as the training cohort, and 43 patients from center 2 constituted an independent external validation cohort and were included in the study (Figure 1).

Routine preoperative clinical characteristics and descriptions of the DEB-TACE procedure are provided in the Supplementary Material. All patients were evaluated using multiphase enhanced MRI or CT within 1–3 months after DEB-TACE, with the follow-up endpoint of early peritumoral recurrence, defined as the internal or marginal portion of the lesion that was enhanced in the arterial phase and faded in the venous and delayed phases.

All MRI examinations included arterial phase (AP), portal venous phase (PVP), and delayed phase (DP) images, which were obtained at 15–25 seconds, 50–60 seconds, and 150–180 seconds after contrast injection, respectively. Details regarding the MR acquisition parameters of the two centers are presented in (Supplementary Table S1).

The imaging data were collected using the picture archiving and communication system (PACS) at the Second Affiliated Hospital of Chongqing Medical University and the First Affiliated Hospital of the Army Military Medical University, and patients’ preoperative multiphase enhanced MRI were exported in DICOM format. The MRI underwent resampling through linear interpolation to achieve a voxel size of 1×1×1 mm3, thereby standardizing the voxel spacing. The region of interest (ROI) was delineated by two radiologists using the 3D-Slicer software (version 4.10.2, https://download.slicer.org), which provides a powerful function for semi-automatic segmentation. Contrast-enhanced magnetic resonance imaging (AP, PVP, and DP) was performed to segment the tumor and avoid the surrounding tumor vessels. To capture features from the 3 mm peritumoral area (ROI-external) in the PVP, which has a higher potential for microvascular invasion, a dilation algorithm was applied to obtain the 3-mm wide area. The combined intratumoral and peritumoral areas (ROI-plus) were generated simultaneously. Importantly, non-hepatic regions within the ROI were subtracted either semi-automatically or manually on a slice-by-slice basis, as appropriate. Ultimately, five ROIs were identified from these three phases after the segmentation process for each patient. (Figure 2a).

All handcrafted features were extracted utilizing the PyRadiomics package. The voxel intensity values were discretized using a fixed bin width of 5. A total of 960 quantitative features were calculated from each ROI in accordance with the guidelines set by the Image Biomarker Standards initiative (IBSI) (20), including first-order statistical features, shape features, textural features, and transformation features. Statistics-based textural features can reflect the homogeneity of the images and the arrangement of properties that change slowly or periodically on the body surface, including gray-level co-occurrence matrix (GLCM), gray-level run length matrix (GLRLM), gray-level size zone matrix (GLSZM), neighboring gray tone difference matrix (NGTDM), and gray-level dependence matrix (GLDM) features. There were advanced filters applied using the Laplacian of Gaussian (LoG, sigma 1.0 mm) and wavelet decompositions with all possible combinations of high (H) or low (L) pass filters in each of the three dimensions (HHH, HHL, HLH, LHH, LLL, LLH, LHL, and HLL). Detailed information on these features is available in PyRadiomics (http://PyRadiomics.readthedocs.io/en/latest/).

Each slice of the tumor was bound by a cubic bounding box during data preprocessing to ensure that the entire tumor was contained within the bounding box. Then, an area of 224 × 224 pixels containing the tumor was cropped as the final image input for the DL models. ResNet-34 was used to build a signature on MRI (21). Due to the millions of learnable parameters in DL models, training it is computationally expensive and requires a large number of images. With transfer-learning technology, a DL model can be trained on less data. A DL model with ResNet34 architecture was trained using the ImageNet dataset with PyTorch (version 1.4.1; PyTorch. org). In the ResNet34 model, the fully connected and softmax layers were removed, and the output values of the nodes in the last layer were used as DL features. The bounding boxes of the images on three adjacent MRI slices were combined into a three-channel image as DL model input. To achieve robust prediction, all three-channel images of each tumor were fed into the DL model (Figure 2b). Based on the DL model, a total of 512 DL features were extracted and selected (Figure 2). The average value of the prediction probability of multiple tumor images was calculated as the DL signature.

For the training cohort, a four-step procedure was used for feature reduction. First, interobserver correlation coefficients (ICCs) were calculated to explore the stability and reproducibility of features, and only the features with both inter- and intra-ICCs > 0.80 were considered to have agreeable reproducibility and were chosen for further analysis. The abnormal values were replaced by the median, all features were standardized, and z-score normalization of MRI signal intensities was performed to eliminate the variance of features before selection. Second, we removed redundant and less-relevant features using minimum redundancy and maximum relevance (mRMR). Then, the optimized feature subsets were selected using the least absolute shrinkage and selection operator (LASSO) method with 10-fold cross-validation (Supplementary Figure S1). Finally, a multivariate logistic regression analysis was performed to build the signature. The radiomic signature was generated using a linear combination of selected features weighted by their respective regression coefficients. The cut-off value was then identified using Youden’s index to divide the patients into non-EPR and EPR subgroups.

There are five radiomics models: the arterial phase model (AP), portal venous phase model (PVP), delay phase model (DP), peritumoral 3 mm portal venous phase model (PVP_Pri3mm), tumor plus peritumoral portal venous model (PVP_Plus3mm), and one deep learning radiomics model based on the portal venous phase (PVP_DLR). Moreover, a clinical model and related combined models, such as the deep learning radiomics and radiomics model (PVP_DLRR), deep learning radiomics nomogram (PVP_DLRN), and deep learning radiomics and radiomics nomogram (PVP_DLRRN) were established. The performance of all established models for HCC recurrence was measured using receiver operating characteristic (ROC) analysis, and the area under the ROC curve (AUC) was calculated and compared among cohorts using the DeLong test. In addition, sensitivity and specificity were measured. The net reclassification index (NRI) was calculated to compare the performance of the models.

Statistical analyses were performed using SPSS (version 26.0, https://www.ibm.com/spss) and R (version 4.2.1, available at http://www.R-project.org). The chi-square test or Fisher’s exact test was used for nominal variables. A logistic regression analysis was performed using the “glmnet” package. The diagnostic performances of the models were compared using ROC analysis, and the differences in the AUCs between these models were compared using the Delong test. Receiver operating characteristic (ROC) curves were plotted using the “pROC” package. All statistical tests were two-sided, and statistical significance was set at P <0.05.

A flowchart of the study is shown in (Figure 2). The baseline clinical characteristics and demographics of the training and validation groups are summarized in (Table 1). The enrolled patients were allocated to a training set (n=114) or an external validation set (n=43). The efficacy of DEB-TACE was balanced for patients in the two cohorts, with early peritumoral recurrence rates of 40.3% (n=46) and 46.5% (n=20) for the training and independent external validation cohorts, respectively. Images of the two patients are shown in (Figure 3) and (Figure 4). No significant differences were detected in sex, age, ALT, ChildPugh, HBsAg, cirrhosis, portal hypertension, tumor number, tumor margin, rim enhancement, or peritumoral enhancement between the recurrence and non-recurrence groups (P>0.05). Moreover, the AFP levels (P = 0.031) were significantly different between the two groups in the training cohort. AST level (P = 0.015), tumor size (P = 0.002), and BCLC stage (P = 0.016) also showed statistically significant differences in the external validation cohort. AFP in the training cohort was constructed for the clinical model (Table 1).

A total of 960 radiomics features were extracted from each phase of the MRI, and ICC was used to select 922 (96% remaining) features with high robustness (ICC>0.8). Next, mRMR and Lasso were applied to further select features, and 1, 2, 8, 5, and 1 features with rich information remained in the AP, DP, PVP, PVP_Plus3mm, and PVP_Pri3mm, respectively. Multivariate logistic regression analysis was performed using weighted summation to obtain the radiomic signature. The selected features and their relative coefficients are presented in (Supplementary Table S2). The distribution of radiomic signatures has shown good separability in early peritumoral recurrence. Moreover, the features of the last fully connected layer of ResNet34 were weighted to obtain a deep-learning signature (PVP_DLR). Fifteen features were selected to construct the PVP PVP-DLR model (Supplementary Table S2).

Five radiomics models and one deep learning radiomics model were established (Figure 5), and the efficacy of each model in the training and validation sets is tabulated in (Table 2). Among the radiomics models, the PVP model had a higher efficacy, with an AUC of 0.751 (95% CI, 0.659–0.843) in the training set and 0.691 (95%CI, 0.530–0.853) in the external validation cohort compared to the AP model, with an AUC of 0.724 (95%CI, 0.625–0.823), and the DP model with an AUC of 0.749 (95%CI, 0.659–0.839). Therefore, PVP was selected to explore intratumoral and peritumoral information, and the PVP_Pri3mm and PVP_Plus3mm models were established. The PVP_Plus3mm with an AUC of 0.754 (95%CI, 0.672–0.855) and PVP_Pri3mm model with an AUC of 0.727 (95%CI, 0.630-0.824) had comparable efficacy to the PVP model, but no significant difference (P = 0.916, P = 0.325) remained. In addition, the deep learning radiomics model based on the venous phase (PVP_DLR) had the best and most stable efficacy in the training set, with an AUC of 0.802 (0.717-0.887), and a higher efficacy in the validation set, with an AUC of 0.774 (95% CI: 0.700-0.783). The distribution of the prediction results for each model is shown in (Supplementary Figure S2).

In the training cohort, hand-craft based signature, DL-based signature, and AFP level were independent factors for EPR prediction using backward stepwise multivariable analysis. However, only PVP_DLR was significant, so we combined the signatures to build the PVP_DLRRN model (Table 3, Figure 6a). As shown in (Table 2), the performance of the combined models PVP_DLRR (AUC, 0.804), PVP_DLRN (AUC, 0.797), and PVP_DLRRN (AUC, 0.802) was not higher than that of PVP_DLR (AUC, 0.802), which was further confirmed in the external validation cohorts. There were no significant differences between the combined models (P>0.05). Furthermore, all the combined models showed significantly higher AUCs than the clinical model in the two cohorts, which also outperformed the handcrafted and DL signatures (P<0.05). NRI and IDI analyses revealed that the integration of image signatures into the DLR performed satisfactorily in the two cohorts, indicating an improved classification accuracy for the EPR prediction clinical model. The calibration curves of PVP_DLR, PVP_DLRR, PVP_DLN, and PVP_DLRRN demonstrated that the model-predicted EPR was well-calibrated with the actual observations in the cohorts(P>0.05) (Figure 6b). Additionally, DCA graphically indicated that the DLRRN provided a net benefit over other models over the relevant threshold range in the entire cohort (Figure 6c). (Figures 5, 6) shows correctly classified examples from the EPR and non-EPR, respectively.