Prediction of prognosis of immune checkpoint inhibitors combined with anti-angiogenic agents for unresectable hepatocellular carcinoma by machine learning-based radiomics

Abstract

Objectives

This study aims to develop and validate a novel radiomics model utilizing magnetic resonance imaging (MRI) to predict progression-free survival (PFS) in patients with unresectable hepatocellular carcinoma (uHCC) who are receiving a combination of immune checkpoint inhibitors (ICIs) and antiangiogenic agents. This is an area that has not been previously explored using MRI-based radiomics.

Methods

111 patients with uHCC were enrolled in this study. After performing univariate cox regression and the least absolute shrinkage and selection operator (LASSO) algorithms to extract radiological features, the Rad-score was calculated through a Cox proportional hazards regression model and a random survival forest (RSF) model. The optimal calculation method was selected by comparing the Harrell’s concordance index (C-index) values. The Rad-score was then combined with independent clinical risk factors to create a nomogram. C-index, time-dependent receiver operating characteristics (ROC) curves, calibration curves, and decision curve analysis were employed to assess the forecast ability of the risk models.

Results

The combined nomogram incorporated independent clinical factors and Rad-score calculated by RSF demonstrated better prognosis prediction for PFS, with C-index of 0.846, 0.845, separately in the training and the validation cohorts. This indicates that our model performs well and has the potential to enable more precise patient stratification and personalized treatment strategies. Based on the risk level, the participants were classified into two distinct groups: the high-risk signature (HRS) group and the low-risk signature (LRS) group, with a significant difference between the groups (P < 0.01).

Conclusion

The effective clinical-radiomics nomogram based on MRI imaging is a promising tool in predicting the prognosis in uHCC patients receiving ICIs combined with anti-angiogenic agents, potentially leading to more effective clinical outcomes.

Introduction

Hepatocellular carcinoma (HCC) was the sixth most commonly observed cancer and the third most common cause of cancer-related deaths worldwide in 2020 [1]. Numerous patients suffer from late-stage disease as a result of delayed identification. The median overall survival is typically ranges from only 6 to 8 months [2].

Recently, systematic treatments for patients with advanced HCC include tyrosine kinase inhibitors (TKIs) and immune checkpoint inhibitors (ICIs). TKIs including sorafenib, regorafenib, lenvatinib, bosutinib, and ramucirumab have been used as first-line and second-line treatments for HCC patients [3,4,5,6,7,8]. Additionally, ICIs, which encompass the inhibitors of programmed death 1 (PD-1) and programmed death ligand 1 (PD-L1), have exhibited remarkable therapeutic advantages for patients with advanced HCC [9]. PD-1 monoclonal antibodies, such as nivolumab, pembrolizumab, camrelizumab, and sintilimab, have received approval for their use in advanced HCC following successful clinical trials [10,11,12,13]. Studies have demonstrated this combination therapy has a synergistic effect to improve overall survival (OS) and PFS [14,15,16]. However, some patients may not experience benefits from the combination of TKIs and ICI therapy, which is also associated with significant adverse events, including gastrointestinal hemorrhage, pneumonia, sepsis, and cardiac arrest. This highlights the urgent need for novel biomarkers to inform clinical decision-making [15, 17].

Radiomics is an emerging field of image analysis that transforms images into quantifiable data, which can be analyzed to construct models that support decision-making [18]. This approach can serve as a guideline for the diagnosis, staging, treatment planning, prediction of treatment response, and determination of patient prognosis in HCC [19, 20]. Dong et al. conducted a study where they used a machine learning (ML) algorithm on noninvasive CT imaging to determine the effectiveness of TKI plus anti-PD-1 antibodies treatment in advanced HCC patients [21]. Using radiomics based on MRI is a valuable method for forecasting the reaction to transarterial chemoembolization (TACE) in HCC [22, 23]. It has also been shown to be effective in assessing the effect of combined targeted therapy [24]. PFS has historically been the most commonly used surrogate endpoint in late-stage clinical trials, and early assessment of PFS is a robust surrogate endpoint for OS in immunotherapy trials for HCC [25, 26]. However, there is currently no published research focusing on the use of MRI-based radiomics to prognosticating progression-free survival (PFS) among patients diagnosed with uHCC.

As an emerging technology, artificial intelligence (AI) radiomics has profoundly changed the field of imaging analysis. It has important application in diagnosis, individualized treatment and prognosis of HCC. ML is a classical technique of AI used to generate predictive models that are widely used in HCC research [27]. In this study, the objective is to develop an MRI-based radiomics model using machine learning ML. We also hypothesize that a comprehensive model, which combines radiomic features with clinical features, will enhance the accuracy of predicting PFS in patients with uHCC.

Methods and materials

Patient selection

From September 2019 to August 2022, 111 patients diagnosed with uHCC according to the national comprehensive cancer network guidelines [28] were enrolled in the study. These patients were treated with ICIs and anti-angiogenic agents at the First Affiliated Hospital of Wenzhou Medical University. The study was conducted in accordance with the Declaration of Helsinki and approved by the Ethics Committee of the First Affiliated Hospital of Wenzhou Medical University. Given that this investigation was retrospective, and involved regular examinations utilizing MRI and the analysis of clinical data, obtaining individual informed consent was not considered mandatory. The eligibility criteria for inclusion in this study were as follows: (1) individuals who had been diagnosed with HCC either clinically or pathologically; (2) patients who received complete treatment with PD-1 inhibitors such as toripalimab (Coherus Biosciences, China), sintilimab (Innovent Biologics, China), or camrelizumab (Hengrui Medicine, China), along with anti-angiogenic agents therapy using sorafenib (Bayer Pharma AG, Germany), lenvatinib (Eisai Co., Ltd., Japan), or regorafenib (Bayer Pharma AG, Germany); (3) patients who underwent enhanced MRI before starting the treatment to obtain imaging data; (4) patients who had comprehensive clinical and treatment information accessible. The exclusion criteria for this study are as follows: (1) patients who underwent alternative anticancer therapies, like interventional therapy, were not included; (2) patients diagnosed with other forms of cancer, like cholangiocarcinoma, were not included; (3) patients with liver function classified as Child-Pugh class C were also not included. (4) patients with MRI data of low quality that was challenging to visualize or quantify were not included. The patients were allocated randomly into two cohorts (at a ratio of 2:1): training cohort (n = 74) and validation cohort (n = 37). The flowchart of patient enrollment was illustrated in Fig. 1.

Fig. 1

Flowchart of patient enrollment. HCC, hepatocellular carcinoma; PD-1 Programmed death 1. *For training and validation of the model, use about 2/3 to 4/5 of the sample for training, and the rest for validation [49]. The 2:1 ratio is a commonly used approach in machine learning and radiomics studies to balance the need for a comprehensive training dataset with the requirement for a substantial validation dataset

Treatment

The treatment protocol was determined by an experienced multidisciplinary team. The therapeutic regiments included the administration of ICIs and TKIs. Toripalimab, at a dose of 3 mg/kg body weight or a fixed amount of 240 mg, was intravenously given every two weeks. Camrelizumab and Sintilimab, at a fixed dose of 200 mg, were intravenously administered every three weeks. Sorafenib, given orally, was provided at a fixed dosage of 400–800 mg per day. For patients with a body weight less than 60 kg, Lenvatinib was orally administered at a fixed dose of 8 mg per day, while patients with a body weight of 60 kg or more received a dosage of 12 mg per day. Regorafenib, taken orally, was administered at a fixed dose of 160 mg per day, every 3 weeks within a 4-week cycle.

Follow-up and clinical endpoint

The subjects underwent regular monthly follow-up appointments throughout the initial year via the outpatient service. Subsequently, they were reassessed every 3–6 months. The endpoint of the study was to determine the PFS of patients. PFS refers to the duration between the introduction of ICIs and TKIs and the occurrence of disease progression or death. Disease progression is determined by RECIST guideline (version 1.1) [29]. The deadline for observation was established until March 30, 2023. The mean follow-up time for both the training and validation cohorts was 6 months.

MRI data acquisition

The images, compliant with DICOM standards, were obtained from the Picture Archiving and Communication System (PACS) records. Acquiring the images involved the utilization of an 8-channel phased-array body coil-equipped Siemens 3.0 T Prisma MR scanner from Siemens Medical Solutions in Erlangen, Germany. The imaging protocol included T1-weighted (T1w), T2-weighted (T2w), diffusion weighted imaging (DWI), and multiphase dynamic contrast-enhanced (DCE) imaging. For the enhanced sequences, the following parameters were used: repetition time (TR) = 2 or 3 ms, echo time (TE) = 1 ms, field of view (FOV) = 420 × 420 mm, matrix (Mat) = 320 × 208, flip angle (FA) = 11°, and slice thickness = 3.0 mm. The MRIs of all patients were conducted using 3.0-T machines within a period of 3 weeks prior to treatment. For the purpose of this study, only the enhanced sequences from the baseline MRI data were selected. Gadodiamide Injection (Omniscan^®,0.2 mL/kg body weight) was intravenously administered with a flow rate of 1 mL/s. MR scans were performed in the arterial phase, portal venous phase, and delayed phase at 20–30 s, 60–70 s, and 180 s after injection, respectively.

Tumor segmentation

The feature data in this study were obtained from specific regions of interest (ROIs). Each patient underwent multiparameter MRIs using a 3.0-T MR scanner. The images were then imported into a 3D Slicer software (version 4.11; https://www.slicer.org/) for manual segmentation. Tumor segmentation was conducted by a radiologic oncologist and a proficient radiologist (5 years of radiological experience) on each transverse slice of 3 enhancement phases. The two reviewers were blinded to clinical data and then independently assessed the images. In cases where patients had multiple liver lesions, the largest lesion was assessed to extract radiomics features. Additionally, a senior radiologist (20 years of experience in abdominal tumor radiology) validated and evaluated each segmentation. Any discrepancies or issues that arose during the segmentation process were resolved by seeking the opinion of another senior radiologist.

Feature extraction

This study used the PyRadiomics platform (version 3.0.1, https://pyradiomics.readthedocs.io/en/latest/index.html) to employ various feature algorithms to quantify information on first order features, shape-based features, texture features (Fig. 2). MR images were resampled into isotropic resolution of 2.0 × 2.0 × 2.0 mm³ to eliminate the inconsistency of spatial resolution. This study utilized this platform to extract 971 features from each phase of the medical images. To ensure consistency and comparability, all radiomics features were normalized using z-scores.

Fig. 2

Radiomics workflow in this study. a Image segmentation was performed in arterial phase, portal venous phase, and delayed phase of Contrast-enhanced MR images. b Features extracted from tumors were classified into shape, first-order statistics, texture, and the features were extracted by least absolute shrinkage and selection operator (LASSO) regression and random survival forest. c clinical nomogram and combined nomogram were established. d Two nomograms were compared by employing time-dependent receiver operating characteristics (ROC) curves, calibration curves, and decision curve analysis. Shape features capture the geometric properties of the tumor, which can be associated with tumor growth and invasiveness. First-order statistics describe basic image characteristics like mean intensity and variance. Texture features reflect the spatial distribution of pixel intensities, indicating tumor heterogeneity and complexity

Development and validation of PFS-Predictive models

All radiomics features were analyzed using a univariable Cox proportional hazard regression approach. Further investigation only considered variables that demonstrated a significant level (P < 0.05) in the univariate analysis. The initial feature set consisted of 2736 features extracted from the MRI images. While this high-dimensional data provides a comprehensive representation of the imaging information, it also increases the risk of overfitting. Therefore, the optimal radiomics features were selected using the LASSO algorithms. The Rad-score for three phase of enhanced MRI was calculated using a linear combination of selected features, which are weighted by their respective LASSO coefficients. Based on Cox regression, a prediction model including three phases of rad-score was established.

Random survival forest (RSF) is an ensemble tree method used to analyze right-censored data, which is a comprehensive method of random forest and survival analysis [30]. RSF can rank the importance of variables, and is known for its high predictive ability and interpretability of the relationship between variables. The RSF method was also used to incorporate the selected features to obtain the final Rad-score, and use Rad-score to build a prediction model.

Subsequently, the prediction ability of the two models was compared by C-index. The more predictive Rad-score enters later studies. Using the Rad-score obtained through the computational model, the participants were categorized into two groups based on their risk level: high-risk signature (HRS) group and low-risk signature (LRS) group. The threshold for categorization was determined by taking the median value. To analyze the relationship between the Rad-score and PFS, Kaplan-Meier survival analysis was conducted in both the training cohort and the validation cohort. In addition to the radiomics model, clinical characteristics were also considered in predicting PFS. Univariate and multivariate logistic regression analyses were conducted to select independent risk factors. Independent clinical factors were utilized to build a Clinical nomogram, while both the Rad-score and clinical risk factors were employed to create a Combined nomogram. These predictive models were then employed to forecast PFS in two cohorts.

Statistical analysis

This study utilized the t test or Mann‒Whitney U test for continuous data analysis for between-group comparisons. The chi-square test or Fisher exact test were used to analyze categorical data. The Kaplan–Meier method was employed for survival curve analysis, which was later compared using the Log rank test. Additionally, Harrell’s concordance index (C-index), time-dependent receiver operating characteristics (ROC) curves were employed to assess the forecast ability of the risk models. Moreover, calibration curves, decision curve analysis, subgroup analysis and N-fold cross-validation were used to measure the agreement between the predicted probabilities generated by our model and the actual observed outcomes and to evaluate the net benefit of using our predictive model across a range of probability thresholds. R software (version 4.0.3, http://www.rproject.org/) was utilized for conducting statistical analyses. The packages employed in the analyses include ‘survival’, ‘rms’, ‘ggplot2’, ‘RMS’, ‘survminer’, ‘randomForestSRC’, ‘dplyr’, ‘timeROC’, and ‘interp’. The code for the modeling described in this paper is available on GitHub at https://github.com/NI-design56/effective-funicular.

Results

Baseline characteristics

In this study, 111 patients who were treated with ICIs and anti-angiogenic agents were included. The mean age was 56.1 years, with 94.6% being male. Tumor size was significant, with 66.7% having tumors larger than 5 cm, and 49.5% having multiple tumors. Portal vein invasion was present in 47.7% of patients. The detailed flowchart of the study was shown in Fig. 1. Subsequently, 74 patients were designated for the training cohort, while 37 patients were allocated to the validation cohort. No significant variations in the baseline characteristics of the patients were observed between the two cohorts (Table 1).

Table 1 Patient baseline characteristics in training and validation cohort

Feature extraction

Each phase consisted of 912 features, with a total of 2736 features initially selected. Subsequently, univariable Cox proportional hazard regression analysis and the LASSO algorithm were applied in each phase, resulting in the retention of 32 features (12, 10, and 10 in the arterial, portal venous, and delayed phase images).

Rad-score calculation

Rad-score calculated by lasso-cox regression

The computational formulae for determining the Rad-score for each patient across the three distinct imaging phases is detailed in the supplementary materials provided with this study. The Rad-score of three phases was used to predict PFS by cox regression. The C-index of the training cohort was 0.760 (95% confidence interval (CI) 0.723–0.789), and the C-index of the Validation cohort was 0.748 (95%CI 0.692–0.763).

Rad-score calculated by RSF

The Grid Search (GS) method was used to determine the combination of mtry and nodesize values when the out-of-bag (OOB) error rate was lowest. As shown in Fig. 3a, under the parameter combination of mtry = 20 and nodesize = 9, the OOB error rate of the RSF model was the lowest (29.2%). Through parameter debugging, the error rate of the model tends to be stable when ntree was 1000 (Fig. 3b). Using the parameters obtained by the above method, the RSF model was constructed, and the importance of the independent variables was ranked using the VIMP method, as shown in the Table 2. Ultimately, 32 features were identified that were considered relatively valuable in all three phases. Considering the potential interrelationships of texture features, the established RSF model algorithm was used to predict PFS by calculating the Rad-score. The C-index of the training cohort was 0.837, and the C-index of the validation cohort was 0.830. Since the C-index of the RSF model was slightly higher than that of the Lasso-Cox regression model, the RSF model was finally used in this study to calculate the Rad-score.

Fig. 3

a Tuning parameter of the RSF model. The OOB error rate is a measure of model accuracy, with lower rates indicating better predictive performance. The black sign in the figure is the parameter combination of mtry = 20 and nodesize = 9 and the corresponding OOB error rate is 29.2%. b Curve of the OOB error rate for RSF model

Table 2 Imaging features selected for the construction of combined model

Prediction models Establishing

In the training cohort, a total of 11 clinical variables were found to be predictive by univariate analysis (Table 3). Among these variables, 5 clinical factors consisting of BMI, tumor size, ALB, AFP, and ascites were independently selected as predictors for recurrence through the application of multivariate analysis. The clinical nomogram was constructed using these variables in the training cohort through Cox regression (Fig. 4a). However, when the combined nomogram was established, it was observed that the value of ALB and AFP for predicting survival was relatively insignificant. Therefore, the combined nomogram was modified to include BMI, tumor size, albumin, and Rad-score (Fig. 4b).

Table 3 The results of univariate and multivariate Cox proportional hazard analysis in this study

Fig. 4

Established clinical nomogram for combination regimen of unresectable hepatocellular carcinoma. (a) The Clinic nomogram for progression-free survival (PFS) including selected clinical risk factors of Table 3. For instance, consider a patient with a BMI of 27 (25 points), high albumin levels (23 points), large tumor size (25 points), presence of ascites (10 points), and an AFP score of 6 (6 points), totaling 89 points. Referring to the nomogram, this patient would fall into the high-risk stratum, with a 3-month PFS rate of approximately 0.5, indicating a moderate to high risk of disease progression. (b) The Combined nomogram for RFS including clinical risk factors and rad-score

This study further assessed the predictive effectiveness of the aforementioned two models by comparing the C-index value and ROC curve. Results indicated that the composite clinical-and-radiomics models, which combined Rad-score with clinical risk factors, exhibited superior ROC curves compared to models using only clinical variables (Fig. 5). In terms of the C-index, the combined nomogram (0.846, 95% CI 0.804–0.879) demonstrated better prognosis prediction for PFS compared to the clinical nomogram (0.752, 95%CI 0.692–0.789). A parallel performance was observed in the validation cohort (C-index, 0.709, 95%CI 0.586–0.772 vs. 0.845, 95%CI0.767-0.893) (Table 4).

Fig. 5

a The time-dependent ROC curves of clinical nomogram in the training cohort. b The time-dependent ROC curves of clinical nomogram in validation cohort. c The time-dependent ROC curves of clinical and radiomics nomogram in the training cohort. d The time-dependent ROC curves of clinical and radiomics nomogram in the validation cohort. AUC Area under the curve

Table 4 The performances of the models in the training and validation cohort

Verification of the radiomics nomogram

This study plotted ROC curves for both the clinical nomograms and the combined nomogram separately in the two cohorts. ROC curves evaluate the model’s ability to distinguish patients with different outcomes. The higher the AUC, the stronger the discriminative power. In our study, the AUC of the combined model (0.966 for 9-month PFS in the training cohort) was higher than that of the clinical nomogram (0.850), indicating better risk stratification. A higher AUC helps identify patients at high risk of disease progression, thereby enabling more targeted treatment plans. The calibration curves (Fig. 6) showed a strong agreement between the predicted and observed PFS in both cohorts. Moreover, Calibration curves assess the consistency between predicted and observed outcomes. Well-calibrated predictions ensure that the model’s outputs can be confidently used for treatment planning and resource allocation. Our combined model showed good calibration in both cohorts (Fig. 7).

Fig. 6

Calibration curves of clinical nomogram and clinicalandradiomics nomogram in the training (a, c) and validation cohort (b, d)

Fig. 7

Decision curve analysis for clinical and clinicalandradiomics nomogram models in the training (a, b, c, d) and validation cohorts (e, f, g, h)

Furthermore, to validate the robustness of our model, we conducted subgroup analyses based on body mass index (BMI). Patients were categorized into high BMI (BMI ≥ 25 kg/m²) and low BMI (BMI < 25 kg/m²) subgroups. The combined nomogram maintained high predictive accuracy within both BMI subgroups, with C-index values of 0.809 (95% CI 0.756–0.848) in the high BMI subgroup and 0.874 (95% CI 0.804–0.931) in the low BMI subgroup, confirming the model’s generalizability across different BMI categories. This subgroup validation using BMI subgroups further reinforces the clinical utility of our model in predicting PFS in patients with uHCC.

Additionally, we used cross-validation to verify the survival prediction model. The cross-validation results showed that the combined nomogram had a significantly lower prediction error compared to a reference model, with a prediction error of 0.071 for the combined nomogram versus 0.129 for the reference model (Fig. 8). This indicates that the combined nomogram not only performs well in terms of discrimination and calibration but also has a robust predictive performance, further validating its reliability and utility in clinical settings.

Fig. 8

Combined model validation: prediction error curves

Prognostic validation of rad‑score

The Rad-score has been validated as a standalone prognostic determinant for individuals diagnosed with uHCC. The classification of Rad-score was determined by establishing a threshold using the median value. This approach effectively segregated patients into either HRS or LRS groups. Survival analysis was conducted in this study using the Kaplan-Meier method, with the main focus being on monitoring disease progression as the primary endpoint (Fig. 9).

Fig. 9

Kaplan–Meier diagrams of HCC patients stratified according to rad-score a PFS survival curve of patients in the training cohort. b PFS survival curve of patients in the validation cohort

Discussion

HCC is one of the most common malignant tumors worldwide, with a particularly poor prognosis for patients. The prognosis of HCC is characterized by significant heterogeneity, and the disease is challenging to treat. Therefore, there is an urgent need to characterize the disease status of individual patients to enable patient stratification and personalized treatment, ultimately maximizing prognosis. In this study, we employed machine learning-based radiomics to predict the prognosis of patients receiving immune checkpoint inhibitors in combination with anti-angiogenic drugs for uHCC. This approach could assist in identifying high-risk patients who would benefit from appropriate treatment.

At present, there are multiple treatment alternatives accessible for advanced and uHCC, involving local-regional therapy, immunotherapy, and molecular targeted medications [31]. The emergence of ICIs has significantly enhanced the prognosis of different solid tumors, including HCC. Moreover, it has displayed synergistic impacts when merged with anti-angiogenic agents [32]. A study was carried out by Shigeta et al. to examine the effects of a combination therapy consisting of anti-PD-1/VEGFR-2 on survival rates in HCC using orthotopic transplantation or induced mouse models [33]. By concurrently targeting VEGFR-2 and PD-1, this dual blockade can promote the normalization of blood vessels, alter the immune microenvironment, and augment the anti-tumor immune response in HCC. Since selecting the most suitable treatment approach for advanced HCC patients from the available options can be challenging, it is crucial to identify potential HCC patients who will respond positively to ICI and anti-angiogenic agents. Optimal treatment strategies and the attainment of therapeutic success heavily rely on the accurate identification of respective conditions. Regrettably, there is a dearth of studies that have prioritized the implementation of ML radiomics techniques in prognosticating the responsiveness of HCC to combination therapy consisting of ICIs and anti-angiogenic agents. In this investigation, researchers adeptly devised radiomics model utilizing ML to forecast the reaction to ICIs in concurrence with anti-angiogenic agents.

Radiomics analysis involves the use of high-throughput feature extraction algorithms to quantitatively measure macroscopic disease features within and between tumors [18, 34], extending beyond the field of HCC to explore its clinical value in different types of tumors. ML is a well-established AI technique employed to develop predictive models, and it is extensively utilized in hepatology research [35]. Previous studies have demonstrated that AI radiomics holds significant promise in forecasting the prognosis of patients with HCC. Akai et al. included 127 hepatocellular carcinoma (HCC) patients who underwent initial hepatectomy and developed a radiomic model utilizing the Random Survival Forest (RSF) algorithm to effectively distinguish between high-risk and low-risk groups for survival risk, demonstrating high consistency [36]. Wang et al. included 201 HCC patients who underwent curative hepatectomy in a retrospective, multicenter study. The 30 radiomics features most associated with survival were selected from preoperative MRI, and combined with preoperative AFP and AST as independent clinical risk factors, the model was constructed by RSF. The model demonstrated a good ability to predict 5-year survival, with a mean AUC of 0.9804 and 0.7578 in the training and validation cohorts, respectively [37]. Bo Z et al. introduced a valuable machine learning radiomics model that showed excellent performance in predicting the response to lenvatinib monotherapy in uHCC, achieving the highest AUC of 0.97 [38].

In the current investigation, researchers constructed a combined nomogram that integrates the Rad-score and three clinical characteristics (body mass index, tumor dimension, and albumin levels) to anticipate the efficacy of ICIs coupled with anti-angiogenic agents in managing uHCC. In this study, the combined model achieved a C-index of 0.846, which is significantly higher than predictive models that rely solely on clinical factors [39]. The models formulated in this study can provide valuable guidance in the judicious administration of ICIs and anti-angiogenic agents in real-world medical practice. The models formulated in this study can provide valuable guidance in the judicious administration of ICIs and anti-angiogenic agents in real-world medical practice. The combined nomogram can identify patients at higher risk of disease progression, allowing for the customization of treatment strategies. For instance, a patient classified as high risk might be considered for more aggressive or alternative therapies, potentially leading to improved patient outcomes. Furthermore, the combined nomogram provides a quantitative risk assessment, which can be used to inform patients about their prognosis and potential treatment outcomes. This information is crucial for shared decision-making, allowing patients to understand their disease trajectory and participate in choosing the most suitable treatment path.

In some radiomics studies, a large number of radiological features have been employed without appropriate feature selection methods, leading to overfitting and consequently poor performance on new, unseen data [40]. This study adopted a rigorous feature selection process using univariate Cox regression and the LASSO algorithm. This approach helps to reduce overfitting by selecting only the most predictive features. Additionally, this study utilized RSF to construct our model, which is less prone to overfitting compared to other machine learning algorithms [41]. Moreover, some radiomics studies have focused solely on imaging features without integrating them with clinical data, which can limit the model’s predictive power as clinical factors often provide important additional information [42]. This study integrates radiological features with clinical variables such as BMI, tumor size, and albumin levels. This integrated approach provides a more comprehensive view of the patient’s condition and enhances the model’s predictive capability.

The use of MRI in this study is a significant differentiator from other imaging modalities such as CT. MRI provides superior soft tissue contrast compared to CT, which is crucial for imaging the liver and detecting small tumors or subtle changes in tumor morphology. This enhanced contrast allows for more accurate delineation of tumor boundaries and the identification of intra-tumoral heterogeneity, which is particularly important in the management of HCC. Furthermore, MRI offers multi-parametric imaging capabilities, including T1-weighted, T2-weighted, diffusion-weighted imaging (DWI), and dynamic contrast-enhanced (DCE) sequences. Each sequence provides unique information about the tumor’s characteristics, such as cellularity, vascularity, and diffusion properties. This multi-parametric approach provides a more comprehensive view of the tumor’s biology and behavior, enhancing the predictive power of our model [43,44,45]. Variations in MRI imaging protocols and scanner settings can significantly impact the reproducibility of radiomic features. To mitigate this issue, the study used a standardized imaging protocol across all patients in our study. This included the use of a consistent MRI scanner model and settings, ensuring that the images were acquired under the same conditions.

Nomograms are commonly used predictive models in evaluating tumor prognosis. Since a majority of hepatocellular carcinoma (HCC) instances are linked to chronic hepatitis and cirrhosis of the liver, the prognostic outcome for liver cancer relies not solely on the tumor, but also on various factors encompassing the patient’s overall nutritional condition, the localized inflammatory microenvironment, and the functioning of the liver reserve. This study found that BMI, tumor size, and albumin levels are associated with the prognosis of uHCC patients receiving combination therapy. BMI serves as a straightforward and intuitive indicator to measure individual obesity levels. Previous research findings support the observation that patients with a higher BMI exhibit an extended period of PFS subsequent to treatment for HCC [46]. This could be attributed to the fact that advanced HCC is a highly consumptive disease, and patients with higher BMI may have better nutritional, physical, and mental statuses, along with greater energy reserves to combat tumor consumption. In addition, tumor size and albumin are important variables in predicting the prognosis of HCC. Many staging systems use a tumor size of 5 cm as a cutoff point for treatment and prognosis of HCC [47, 48]. However, this study did not yield any statistically significant findings concerning ALBI grading. This outcome might be attributed to the inclusion of patients with adequate liver function and the lack of substantial variations in bilirubin levels among them. Furthermore, AFP was omitted from the final integrated model, likely because the HCC patients encompassed in our study were predominantly in the intermediate and advanced stages, with the majority exhibiting high AFP values. When other clinical and radiological features were incorporated into the model, the hazard ratio (HR) of AFP failed to reach statistical significance, and its exclusion did not notably impact the overall performance of the nomogram.

This study has several limitations that warrant acknowledgment. Firstly, the retrospective nature of the study design introduces potential selection biases in the collection and retrieval of MRI and clinical data. This inherent limitation highlights the necessity for additional prospective studies to validate the findings of this investigation. Secondly, there is variability in the use of ICIs and TKIs among patients. However, no studies have yet conducted comparative analyses to establish differences in effectiveness between various PD-1 inhibitors. Thirdly, the short-term follow-up, averaging six months, may restrict our ability to obtain long-term outcomes and fully assess the durability of PFS predictions. To address this limitation, we performed internal validation using N-fold cross-validation to evaluate the model’s performance on a subset of the data, thereby providing a more conservative estimate of its predictive power. Lastly, the study was conducted at a single institution, which may limit the generalizability of the findings to patients in different regions and institutions. In light of these limitations, we plan to conduct subgroup analyses by drug class in future studies to assess the impact of different ICIs and TKIs on treatment outcomes. We also intend to extend the follow-up period in our ongoing research to obtain long-term results, which will enable us to more accurately assess the durability of PFS predictions and validate the model’s performance over a longer time horizon. Furthermore, we plan to explore the use of automatic segmentation techniques, such as deep learning-based methods, to improve the accuracy and repeatability of tumor segmentation. Given that the combination of TKIs and ICIs is a relatively new approach, large-scale, multicenter subgroup analysis studies should be prioritized for future research.

Conclusion

This study successfully developed an MRI-based radiomics model that accurately predicts progression-free survival in patients with unresectable hepatocellular carcinoma. The combined nomogram, which integrates clinical factors and radiomic features, demonstrated superior predictive accuracy compared to clinical nomogram alone. The model’s effectiveness across BMI subgroups highlights its broad applicability. Future research will focus on external validation and the potential inclusion of additional biomarkers to further enhance model performance.