Characterizing immune profiles in hepatocellular carcinoma patients benefiting from pembrolizumab and lenvatinib using machine learning.

Introduction
Hepatocellular carcinoma (HCC) represents the predominant form of primary liver cancer, ranking as the sixth most common cancer and the third leading cause of cancer-related mortality worldwide [1]. More than half of HCC cases are diagnosed at an advanced stage, necessitating systemic therapy. Combination immunotherapy utilizing immune checkpoint inhibitors (ICIs) has emerged as a promising therapeutic strategy for unresectable HCC (uHCC) [2–4]; however, the objective response rate to ICI-based immunotherapy remains modest, ranging from 20 to 30%. Although the phase III LEAP-002 study did not meet expectations [5], the combination of pembrolizumab and lenvatinib (PL) demonstrated comparable tumor responses and prolonged survival benefits in responders, particularly in Asian real-world data [6, 7]. Nevertheless, the identification of an effective biomarker for predicting tumor response to combination immunotherapy remains an unmet need.
Given that changes in the composition of circulating immune cell profiles (ICPs) may be associated with clinical responses to ICIs in patients with advanced cancer, these changes may be valuable for predicting the patients’ clinical responses to combination immunotherapies [8, 9]. Our previous studies have demonstrated that circulating ICPs can be used to predict the efficacy of nivolumab monotherapy in patients with advanced HCC [10, 11]. In this study, we investigated the differences in circulating ICPs between patients with uHCC who achieved objective responses to PL and those who did not. Additionally, we adopted a machine learning model that determines the contribution of each immune cell subset to ICPs. Considering that the LEAP-002 trial did not achieve its original design endpoints, the generalizability of its findings to uHCC patients remains limited. Moreover, real-world data have shown that patients with stable disease (SD) following PL therapy derived only modest clinical benefit compared to objective responders, with significantly shorter progression-free survival (9.0 vs. 15.3 months) and overall survival (12.0 months vs. not reached) [6]. To better elucidate the immune checkpoint profiles (ICPs) associated with effective tumor response, we focused our analysis on objective responders and non-responders, excluding patients with SD. This approach aimed to identify immune signatures with clearer clinical relevance.

Materials and methods

Patients, treatment, and outcome measurement
From July 2019 to July 2023, 51 patients who received front-line combination therapy of pembrolizumab plus lenvatinib for uHCC in Taipei Veterans General hospital were prospectively enrolled in a biomarker study prior to the initiation of ICI treatment [12]. During this period, the ICIs approved by Taiwan Food and Drug Administration for HCC included nivolumab, pembrolizumab and atezolizumab. Unlike sorafenib or lenvatinib, the combination of atezolizumab and bevacizumab was not reimbursed by the national health insurance until August 2023, leading to substantial out-of-pocket costs for patients [13]. Consequently, based on the meaningful tumor response and survival benefit demonstrated by pembrolizumab plus lenvatinib in real-world data [6, 7], this combination was selected through shared decision-making during this period in Taiwan [14]. Peripheral blood mononuclear cells (PBMCs) were collected from patients prior to the initiation of immunotherapy, provided they were free of active infection or had not taken antibiotics, antifungal drugs, steroids, non-steroidal anti-inflammatory drugs, or immunomodulatory agents within four weeks before starting ICI treatment. Immunophenotyping of PBMCs and its association with tumor response were analyzed. Additionally, patients with a radiology-confirmed objective response (responders, R) and those with primary progressive disease (non-responders, NR) were identified based on RECIST version 1.1 criteria. The study algorithm is presented in Supplementary Fig. 1. The study was approved by the Institutional Review Board of Taipei Veterans General Hospital (approval numbers: 2017-09-007CC, 2019-07-007AC, and 2019-08-006B). All blood collection were performed after participants signing the informed consent form.
Lenvatinib was administered according to the standard body weight (BW)-based recommendations: 12 mg/day for patients with a BW ≥ 60 kg or 8 mg/day for BW < 60 kg. Pembrolizumab was administered intravenously at a dose of 200 mg or 2–3 mg/kg every three weeks [5, 6]. Safety assessments and grading were performed using the National Cancer Institute Common Terminology Criteria for Adverse Events (NCI CTCAE, version 5.0). Tumor imaging was conducted by computed tomography or magnetic resonance imaging at screening and every 9 weeks during treatment [5, 6]. Tumor response was assessed using RECIST version 1.1 criteria based on contrast-enhanced radiologic images [15].

Reagents and antibodies
The reagents and antibodies used in this study are listed in Supplementary Tables 1 and 2. All materials were aliquoted and stored under conditions recommended by the manufacturer.

Preparation of PBMCs
PBMCs were prepared following the manufacturer’s protocol. Briefly, blood was prepared with PBS in a 1:1 dilution and separated by Ficoll–Paque (1.077 g/mL; Cytiva, Marlborough, MA, USA) density gradient centrifugation. PBMCs were collected from the buffy coat, Suspended in 10% DMSO in PBS, and cryopreserved at − 80 °C until analysis.

Immunostaining and flow cytometry
Immunostaining of PBMCs was conducted, as previously described [16]. Briefly, frozen PBMCs were thawed, and debris and dead cells were removed. Cells were sequentially stained with antibodies against surface and intracellular markers. Fluorescence was analyzed using a Navios flow cytometer (Beckman Coulter, Brea, CA, USA), and the fluorescence data were acquired and evaluated with Kaluza analysis software (version 1.3; Beckman Coulter).

Data acquisition, pairwise comparison, and machine learning
Immune cell Subsets are defined in Supplementary Table 3, for which the filtration pedigrees are depicted in Supplementary Fig. 2. Following data acquisition, myeloid cells (CD45+) were initially isolated. These myeloid cells were then divided into two populations: CD14/19+ cells (monocytes and B cells) and CD14/19− cells. Subsequently, the CD14/19− cells were gated using CD3 and CD56 markers to identify T cells (CD3+CD56−), NKT cells (CD3+CD56+), and NK cells (CD3−CD56+). Within these T, NK, and NKT cell populations, CD4+ and CD8+ subpopulations were further delineated using CD4 and CD8 markers, respectively. Finally, each identified population was analyzed for the expression of CD25, PD-1, and PD-L1 to identify regulatory cells (CD25+), PD-1+ cells, and PD-L1+ cells. Because pembrolizumab is an anti-PD-1 antibody that inhibits the interaction of PD-1 with PD-L1 [17], clinical responses to PL can be determined based on the proportion of PD-1-related immune cells. Hence, we constructed a series of pedigrees for examining the proportion of T-cell and NK-cell lineages, including CD4+, CD8+, FoxP3+CD25+ (regulatory cells), and their PD-1+/PD-L1+ subpopulations [18]. Notably, a study reported that patients with biliary cancer who did not demonstrate benefits from anti-PD-1 therapy had a higher proportion of circulating monocytes than those who demonstrated benefits from the therapy [19]. Consequently, in this study, we included monocytes in our pedigree. We conducted pairwise comparisons of the proportions of immune cell subsets, expressed as the proportions of each subset and the total CD45+ cells, between the R and NR groups.
This study adopted the machine learning (ML) model from a previous study [16], with the detailed parameter settings were provided in Supplementary Table 4. Briefly, the data of patients were randomly divided into training and validation datasets (3:1 ratio). We excluded immune cell subsets with highly correlated features (|Pearson’s coefficient| ≥ 0.8), and features were selected through recursive feature elimination. This model operates using a random forest classifier with repeated stratified k-fold cross-validation, and Optuna [20] is used for hyperparameter optimization. Receiver operating characteristic (ROC) analysis and a confusion matrix were used to evaluate the ML model. In addition, the contribution of each immune cell subset to ICPs was analyzed using SHapley Additive exPlanations (SHAP) method [21]. The ML model was modified and implemented in Python.

Statistical analysis
Continuous variables were expressed as the median (interquartile ranges [IQR]), while categorical variables were analyzed as frequencies and percentages. Pearson chi-square analyses or Fisher’s exact tests were used to compare categorical variables, and the Mann-Whitney U tests or Kruskal-Wallis tests with Dunn’s multiple comparisons were applied for continuous variables. The optimal cutoff values of bacterial abundance to predict tumor response were assessed using the area under receiver operating characteristic (ROC) curves (AUC). For all analyses, p < 0.05 was considered statistically significant. All statistical analyses were performed using the Statistical Package for Social Sciences (SPSS 26.0 for Windows, SPSS Inc, Chicago, IL, USA) and Prism 9.5.1 (GraphPad Software, Huston, CA, USA). ML model was developed in a Python environment using the scikit-learn [22], Optuna, and SHAP packages.

Results

Demographics of participants
Participant characteristics are detailed in Table 1. The participants were predominantly male (76.5%), with a median age of 62.6 years, and the majority had chronic hepatitis B virus infection (60.8%). At baseline, most participants had Child-Pugh class A liver function (70.6%), although 60.8% had albumin–bilirubin grade greater than 1. In addition, 45 patients (88.2%) exhibited Barcelona Clinic Liver Cancer stage C cancer. Following PL, 16 participants exhibited objective responses, with a complete response in one subject and partial responses in 15 participants. In addition, 22 participants exhibited stable disease, and 11 participants exhibited progressive disease. We stratified 16 participants with objective responses and 11 with progressive diseases into the R and NR groups, respectively. The two groups exhibited comparable liver function, incidence of portal vein invasion, extrahepatic metastasis, tumor characteristics, and serologic profiles.

Peripheral ICP differs between R and NR groups
Compared with the NR group, the R group exhibited significantly higher proportions of peripheral total T lymphocytes (38.94% vs. 28.95%, p = 0.026) and CD8 T lymphocytes (10.78% vs. 8.74%, p = 0.014). By contrast, the proportions of monocytes were marginally lower in R than NR (3.72% vs. 8.91%, p = 0.397) (Fig. 1).
Regarding immune tolerance, the two response groups exhibited similar percentages of regulatory T cells, although a significant decrease in CD4 NKTreg cells was noted in the NR groups (Fig. 2A). Additionally, the R group exhibited higher proportions of PD-1+ CD4 T cells (0.57% vs. 0.41%), CD8 NK cells (0.05% vs. 0%), and PD-L1+ CD8 T cells (3.34% vs. 2.15%) compared with the NR group (Fig. 2B). By contrast, the proportions of PD-L1+ monocytes were higher in the NR group than in the R group (1.35% vs. 8.25%, Fig. 2C).

The trained ML model accurately discriminates between R and NR groups in terms of their ICPs
To determine the contribution of immune cell subsets with distinct ICPs to tumor response, ICP data from the R and NR groups were used to train an ML model to differentiate R from NR patients. Before model training, the ICP data were divided into training and validation datasets in a 3:1 ratio. These datasets were similar in all demographic characteristics except for gender (Supplementary Table 5). After training, the contribution of each immune cell subset in the discrimination was evaluated by SHAP method. As shown in Figs. 3A & B, CD8+ T cells was the highest contributed feature, following by PD-1+ CD8 NK, PD-L1+ monocyte, PD-L1+ CD4 T, PD-L1+ CD8 T, and total T cells. While applying validating dataset in testing the performance of the trained ML model, the trained ML model could correctly determine whether a patient was in the R (precision = 1.00) or NR (precision = 0.80) group on the basis of their ICP (Figs. 3C & D). The ROC analysis revealed an area under the curve (AUC) of 0.9167 (Fig. 3E). This AUC value indicated that the model could effectively differentiate between the ICPs of the R versus NR groups. In the SHAP analysis, six previously identified features were highly contributed to discrimination between ICPs from the R and NR groups as well. CD8 T cells exerted the most influence on the discrimination ability of the model (Fig. 3F). In summary, the ML model successfully discriminated between ICPs in the R versus NR groups; information on CD8 T cells was most instrumental to this discrimination.

Discussion
The present study aimed to characterize the patterns of ICPs for T- and NK-based cell lineages in uHCC patients who achieved objective tumor response following PL therapy. In pairwise comparisons, responders exhibited higher proportions of T (particularly CD8 T cells), CD4 NKTreg cells, and PD-1+ subpopulations of CD4 T and CD8 T cells than non-responders. The ML model discriminated between ICPs between responders and non-responders to PL therapy. The results indicated that alterations in CD8 T cells, T cells, PD-1 subpopulations of CD8 NK and CD4 T cells, and PD-L1 subpopulations of monocytes and CD8 T cells were instrumental to this discrimination, with alterations in CD8 T being the most instrumental.
The efficacy of combination immunotherapy for tumor control in uHCC patients varies from 20.1 to 34%, which is due to the heterogeneity of the tumor microenvironment [4, 5, 23, 24]. Various biomarkers, including circulating miRNA, tumor DNA, and serum α-fetoprotein, have been proposed for predicting clinical responses to targeted ICIs [25]. Reliable biomarkers such as intratumoral PD-L1 expression, dMMR, and MSI-H have been applied in advanced solid tumors, but these biomarkers, namely dMMR and MSI-H, have been detected in < 2% of uHCC patients [25–27]. Furthermore, a study reported that only 26% of uHCC patients with high PD-L1 expression exhibited clinical responses to ICIs [28]. Thus, an effective approach should be developed to identify patients likely to benefit from these therapies.
This study is the first to characterize the patterns of ICPs in uHCC patients who exhibited objective responses to PL. Previous research on other cancers has revealed similar trends. Rubio et al. noted that patients with metastatic solid tumors who exhibited persistent survival benefits from ICIs had higher proportions of activated CD8 T cells and dendritic cells and lower proportions of myeloid-derived suppressor cells and monocytes than those without survival benefits [29, 30]. Our previous studies on uHCC patients treated with nivolumab suggested that compared with nonresponders, responders exhibited higher proportions of PD-1+ or PD-L1+ CD8 T cells and PD-L1+ monocytes and lower proportions of PD-1+ CD8 NKT cells and PD-1+ B cells compared to NRs [10, 11]. Using these findings, an effective approach can be developed for identifying uHCC patients likely to benefit from ICI therapy.
Our study results suggest that uHCC patients who responded to PL exhibited higher proportions of CD8 T cells and their PD-1+ subpopulations. This finding, along with existing literature, underscores that the correlation between peripheral CD8 T cell proportions and clinical response to immunotherapies is contingent upon the specific type of ICI. For instance, a study on anti-PD-L1 therapy, such as atezolizumab plus bevacizumab (A + B) therapy, reported that responders often present with a lower baseline proportion of peripheral CD8 T cells compared to non-responders [31]. Given that the effect of A + B therapy on CD8 T cells primarily involves promoting their migration to the tumor site [32], a higher baseline proportion of peripheral CD8 T cells in non-responders may indicate an inherent defect in the tumor infiltration capacity of these cells [33]. In contrast, anti-PD-1 therapy directly enhances CD8 T cell proliferation and activation [34, 35]. Therefore, a higher baseline proportion of peripheral CD8 T cells suggests a larger pool of readily responsive, antigen-specific T cells capable of benefiting from anti-PD-1 blockade. However, our previous work reported a lower, albeit not statistically significant, baseline proportion of peripheral CD8 T cells in nivolumab responders among HCC patients [10]. This apparent discrepancy may be attributed to the concurrent use of lenvatinib, which is known to enhance T-cell infiltration [36]. In conclusion, the proportion of peripheral CD8 T cells serves as a significant correlate of clinical response to ICIs; however, the directional trend of this correlation between responders and non-responders is fundamentally influenced by the specific type of ICI administered.
The negative correlation between the proportion of CD8 T cells and the clinical response to PL suggests that CD8 T-cell dysfunction may hinder clinical responses. ICIs combined with HCC-specific cytotoxic T cells may provide enhanced treatment efficacy, as demonstrated in rodent models [37]. Smith et al. reported a metastatic nasopharyngeal carcinoma patient who received autologous Epstein–Barr virus-specific cytotoxic T cells followed by nivolumab; the patient exhibited a complete response and disease-free progression until treatment completion [38]. This result supports that combining ICIs with adoptive T-cell therapy enhances tumor control. Currently, the clinical evaluation of three cell products—two chimeric antigen receptor T-cell products targeting CD133 and glypican-3 and one TCR-modified T-cell product targeting HBV-specific antigens—is ongoing [39–41]. As an alternative strategy, a novel ex vivo CD8 T-cell expansion technique targeting unknown tumor neoantigens can achieve tumor control in HCC patients without CD133 or glypican-3 overexpression (approximately 25%) [42].
The present study found that uHCC patients who did not demonstrate benefits from PL exhibited significantly lower proportions of CD4 NKTreg cells than those who demonstrated benefits from PL. NKT cells, which are activated by CD1d-restricted antigens, secrete pro- or anti-inflammatory cytokines, such as IFN-γ and TGF-β, influencing the reactions of T-helper cells and the production of PD-1+ CD8 T cells in the liver [43–45]. NKT cells play roles in the pathogenesis of non-alcoholic steatohepatitis (NASH) and NASH-related HCC [46]. Anti-PD-1 therapy can enhance the activity of NKT cells, which increases the proportion of PD-1+ CD8 T cells, thus increasing the incidence of NASH-related HCC [46, 47]. TGF-β can promote the transformation of NKT cells into NKTreg cells, which do not secrete IFN-γ, thus preventing CD8 T-cell proliferation and reducing HCC progression related to anti-PD-1 therapy [48]. These findings suggest that higher proportions of NKTreg cells in patients demonstrating benefits from PL may reflect improved immune modulation, and that anti-inflammatory therapies can enhance uHCC patients’ clinical responses to immunotherapy.
The present study has three key limitations: its small sample size, its sole focus on PL, and the Limited applicability of its findings to non-viral-related HCC. Although we identified key factors predicting clinical responses in uHCC patients before PL initiation, the small sample size potentiates the risk of overfitting and reduces the robustness of the statistical power. Additionally, our study only examined immune checkpoint proteins in the context of PL; thus, the results may not be applicable to other combination therapies. Further studies involving patients receiving different immunotherapy regimens are needed to assess the broader applicability of immune cell profiles. Furthermore, our findings are primarily based on patients with virus-related HCC, which accounted for 77.8% of analyzed cases. Thus, their generalizability to non-viral HCC is limited. Although a prior study suggested that ICIs may be less effective in non-viral HCC [46], including MASLD-related cases, subsequent data have been inconclusive, and this remains an area requiring further investigation [49]. Due to the small number of non-viral HCC patients in our analyzed cohort (n = 6), subgroup analysis was not feasible. Future studies with larger non-viral cohorts are warranted to validate the applicability of our findings.
This study presents novel patterns of ICPs in uHCC patients who demonstrated benefits from PL. This study identified three key features for predicting clinical responses to other combination immunotherapies, and these critical features can be applied for developing new approaches to enhance the efficacy of tumor control in uHCC patients.

Supplementary Information