The combination of camrelizumab and apatinib represents a promising treatment strategy for patients with advanced hepatocellular carcinoma (aHCC). However, the specific patient populations that may benefit from this combination therapy, as well as the changes in efficacy after adjusting the medication regimen to avoid serious adverse reactions, remain uncertain.
Methods
We employ a quantitative systems pharmacology (QSP) approach to address these significant clinical issues. A QSP model is established by integrating pharmacokinetic data of camrelizumab and apatinib, generating a virtual patient cohort for rapid and reliable virtual clinical studies.
Results
Ultimately, our model identifies the pre-treatment CD8+/Treg ratio, CD4+/Treg ratio, and the density of myeloid-derived suppressor cells (MDSCs) as key predictive biomarkers. Furthermore, through computer-simulated clinical trials, we find that reducing the dose of apatinib in combination therapy to 125 mg can still achieve therapeutic effects comparable to the original dose.
Discussion
These findings provide valuable insights for future drug development and clinical trial design.
advanced hepatocellular carcinomaapatinibbiomarkerscamrelizumabdose optimizationquantitative systems pharmacologyThe author(s) declared that financial support was received for this work and/or its publication. This work was supported by the Joint Funds for the Innovation of Science and Technology, Fujian Province (Grant No. 2024Y9595), the Research Fund for Clinical Studies of Innovative Drugs after Market Launch (Grant No. WKZX2024CX501208), and Fujian Provincial Natural Science Foundation of China (Grant Nos. 2025J08268 and 2024J011068).section-at-acceptanceSystems ImmunologyIntroduction
Hepatocellular carcinoma (HCC) is the most prevalent form of liver cancer, with both incidence and mortality rates steadily rising worldwide (1). Chronic infections with hepatitis B virus (HBV) and hepatitis C virus (HCV) are the primary risk factors contributing to the development of HCC (2). Despite significant advancements in disease screening technologies, the majority of HCC patients are diagnosed at an advanced stage, referred to as advanced hepatocellular carcinoma (aHCC) (3). At this advanced stage, traditional therapeutic approaches—including liver transplantation, surgical resection, and conventional chemotherapy—exhibit limited efficacy, thereby rendering systemic therapy the primary treatment option for patients with aHCC (4).
Systemic treatment for aHCC has long posed a significant therapeutic challenge. Prior to 2017, treatment predominantly relied on anti-angiogenic tyrosine kinase inhibitors (TKIs) such as sorafenib (5). However, the approval of novel first- and second-line therapies has introduced immune checkpoint inhibitor (ICI)-based combination strategies, which have reinvigorated treatment possibilities for aHCC, thereby expanding the range of available therapeutic options (6, 7). Clinical trials, including RESCUE and CARES-310, have demonstrated that the combination of camrelizumab (an ICI) and apatinib (a TKI) significantly prolongs overall survival (OS) and progression-free survival (PFS) in patients with aHCC (8, 9). In light of these promising findings, the Chinese Society of Clinical Oncology (CSCO) Hepatocellular Carcinoma Expert Committee has recommended the camrelizumab-apatinib combination as one of the first-line treatment options for aHCC in its clinical guidelines (10).
Despite the significant efficacy of this combination therapy, the objective response rates (ORR) in the first-line treatment arms of the RESCUE and CARES-310 trials were 34.3% and 25%, respectively. This indicates that a substantial proportion of patients with aHCC do not respond optimally to the camrelizumab-apatinib combination. Therefore, identifying predictive biomarkers to better select patients likely to benefit from this combination therapy has become crucial. The current study investigates a limited range of biomarkers with a relatively small sample size, which affects the reliability of the selection results (11, 12). Additionally, a high proportion of patients required dose adjustments or changes in treatment frequency due to adverse reactions. In the RESCUE and CARES-310 trials, 60.6% (20 patients) and 72.1% (137 patients) of patients experienced grade 3 or higher adverse drug reactions, leading to treatment interruptions or dose modifications. However, the impact of dose reductions on treatment efficacy remains inadequately explored. Given these limitations, exploring reliable biomarkers and optimizing treatment strategies in larger cohorts has become urgent.
Virtual clinical trials, which enable rapid, cost-effective, and efficient drug evaluation in large populations, appear to be a promising solution (13). Quantitative systems pharmacology (QSP), leveraging demographic data, disease mechanisms, and drug properties, has emerged as a preferred approach for conducting virtual clinical trials (14). According to the International Society of Pharmacology, QSP is predominantly applied in the development of cancer immunotherapies to address critical issues such as dose and schedule optimization, combination therapy evaluation, biomarker discovery, and clinical trial design (15). In 2020, Sové et al. developed a QSP platform for immuno-oncology (QSP-IO) that incorporates detailed mechanisms of key immune interactions (16). This platform facilitates the creation of QSP models of varying complexity for IO research. Utilizing this platform, researchers have developed QSP models that account for anti-cancer drug characteristics, tumor immune microenvironment dynamics, and spatial multi-omics data, thereby identifying efficacy-related biomarkers and optimizing therapeutic strategies for diseases such as breast cancer, lung cancer, and liver cancer, which contribute to the standardization of clinical drug use (17–19).
This study addresses two major challenges in the combination therapy of camrelizumab and apatinib for patients with aHCC: the relatively low ORR and the high incidence and severity of adverse reactions. To overcome these therapeutic limitations, we enhanced an existing QSP model of aHCC based on the QSP-IO platform developed by Sové et al., enabling simulation of the therapeutic effects of camrelizumab and apatinib under various dosing regimens. Our primary objective is to identify predictive biomarkers that can guide patient selection and maximize treatment benefit. Additionally, we aim to evaluate the impact of dose reduction on therapeutic efficacy and to optimize dosing strategies for improved clinical outcomes.
Materials and methods
In this study, we further refine a QSP model for aHCC based on the QSP-IO platform developed by Sové et al. to simulate the therapeutic effects of camrelizumab and apatinib. Our objective is to identify predictive biomarkers that can guide precise patient selection and ultimately enhance treatment success rates. Additionally, we aim to investigate the impact of dose reduction on treatment efficacy and optimize the dosing regimen for improved clinical outcomes.
Model development
The QSP model was developed using QSP-IO, a modular platform based on MATLAB (MathWorks, version R2023b, Natick, MA, USA) SimBiology, specifically designed for creating QSP models in IO research (16). The model comprises four compartments that represent different biological environments: the tumor, tumor-draining lymph nodes (LNs), central (blood), and peripheral (all other organs) compartments. This structure aims to comprehensively capture the complete virtual patient. As a modular framework, the model incorporates ten distinct modules that describe various biological components and processes. These include cancer cells, effector T cells (Teff), regulatory T cells (Treg), antigen-presenting cells (APCs), tumor-specific neoantigens, tumor-associated self-antigens, immune checkpoints, bone marrow-derived suppressor cells (MDSCs), macrophages, and therapeutic drug pharmacokinetics.
The model workflow consists of five key steps: (1) Initializing Parameters: The load_parameters function loads the initial set of model parameters; (2) Model Initialization: The simbio_init function creates a simulated biological object, integrating modules that define species, events, and rules; (3) Baseline Simulation: Initial conditions are simulated using baseline parameters and the expected initial tumor diameter; (4) Defining the Dosing Regimen: The regimen is specified, including drug dosage, dosing interval, and treatment duration; (5) Conducting Virtual Clinical Trials: Multiple parameter sets are generated to represent virtual patients, enabling simulations and analyses of clinical outcomes (20). Global uncertainty and sensitivity analyses were conducted using Latin Hypercube Sampling (LHS) and the Partial Rank Correlation Coefficient (PRCC) methods to evaluate the impact of various parameters on the observed outcomes of the model (21). In the PRCC analysis, outcome indicators such as tumor volume, Treg cell density, and macrophage density are particularly noteworthy. Parameters that exhibited significant correlations with the model outputs (PRCC values greater than 0.05 or less than -0.05, with statistical significance) were identified and further analyzed (22).
Model evaluation
The efficacy of the virtual patient model was assessed based on the Response Evaluation Criteria in Solid Tumors (RECIST) v.1.1 (23). RECIST v.1.1 is a classification system utilized to evaluate the response of solid tumors to treatment, primarily based on the percentage change in tumor diameter relative to the baseline or minimum diameter. Tumor diameter is defined as the sum of the longest diameters of all measurable lesions. In this study, tumors were assumed to be spherical, and tumor diameter was calculated from tumor volume.
The ORR and median duration of response (DOR) predicted by the model were calculated, with the median and 95% confidence intervals (CIs) reported. These results were subsequently compared to those from the CARES-310 clinical trial, facilitating the refinement and adjustment of the model’s parameter distributions (9). To align with the clinical setup of the CARES-310 study, the treatment duration was established at 700 days, with tumor volume assessed every 8 weeks. In the CARES-310 trial, the combination therapy of camrelizumab and apatinib yielded an ORR of 25.4% (95% CI: 20.3-31.0%) and a median DOR of 14.8 months (95% CI: 8.4-NR). Considering the significant differences in tumor progression and body tolerance between patients receiving first-line therapy and those receiving second-line therapy or higher, only the ORR from the clinical study involving non-first-line treatments with camrelizumab and apatinib was utilized to preliminarily assess the predictive performance of the model (24, 25).
In the safety analysis, we observed that apatinib monotherapy (750 mg) in aHCC patients was associated with hypertension incidences of 48% (all grades) and 28% (grade ≥3) (25). Notably, when combined with camrelizumab (250 mg), the incidence increased to 70% and 38%, respectively, suggesting a synergistic effect (9). To capture this key safety signal, we added a blood pressure dynamics module based on drug exposure levels of apatinib and camrelizumab, enabling quantitative prediction of hypertension risk. The model parameters were optimized using hypertension incidence data from apatinib monotherapy, and data from the combination with camrelizumab were used to validate the model’s predictive capability for hypertension occurrence.
Model applicationExploration of biomarkers
The developed QSP model was employed to conduct a virtual clinical trial involving 500 patients. The treatment regimen comprised camrelizumab (200 mg, administered biweekly) combined with apatinib (250 mg, taken daily). According to the RECIST v.1.1 criteria, patients were classified as responders or non-responders. Informed by biomarker findings from other QSP models and subgroup analyses from clinical trials of the camrelizumab-apatinib combination, we identified seven key biomarkers relevant to the virtual patient cohort: initial tumor diameter, PD-L1 expression, CD8+ T cell density, Treg cell density, Macrophage density, CD4+/Treg ratio, and CD8+/Treg ratio (18, 26, 27). These biomarkers were investigated to evaluate their potential as indicators of treatment efficacy. The ORR and its 95% confidence intervals for various subgroups were plotted using R (version 4.2.3) and RStudio software (Posit, version 1.4.1717, Boston, MA, USA).
Subgroup analysis has inherent limitations in capturing complex interactions and nonlinear effects among variables. To address these limitations, we integrated machine learning techniques for covariate selection. The process involved several key steps: First, we applied least absolute shrinkage and selection operator (LASSO) regression for initial screening to identify potential biomarkers (28). Subsequently, we utilized these biomarkers to construct a random forest model, further quantifying their importance and predictive performance regarding treatment outcomes. We specifically focused on biomarkers identified in both the subgroup analysis and machine learning selection, considering them as key covariates. Receiver operating characteristic (ROC) curve analysis was performed on the selected key biomarkers to determine optimal cut-off values for predicting treatment response. This approach offers a more comprehensive understanding of the relationship between biomarkers and treatment outcomes, thereby providing stronger support for clinical decision-making.
LASSO regression incorporates an L1 regularization term in the loss function, which compels many regression coefficients to shrink towards zero. This feature selection method reduces model complexity and enhances predictive accuracy, particularly when addressing datasets with numerous irrelevant features. Furthermore, LASSO aids in alleviating issues of multicollinearity, which is prevalent in high-dimensional datasets. Furthermore, we employed Random Forest, a supervised machine learning technique, to evaluate the significance of the selected biomarkers in predicting patient responses to treatment (29). The model was implemented using Python (Python Software Foundation, version 2024.20.0, San Francisco, CA, USA). A parameter grid was established to optimize critical model parameters, including the number of decision trees and the maximum depth of the trees, thereby defining the hyperparameter search space.
Optimization of dosing regimen
In real-world clinical studies on advanced hepatocellular carcinoma (HCC), camrelizumab and apatinib are frequently employed, either as monotherapy or in combination, at fixed doses. The standard dosing regimens are as follows: (1) Camrelizumab (200 mg every 2 weeks) plus apatinib (250 mg once daily); (2) Apatinib monotherapy: 750 mg once daily; (3) Camrelizumab monotherapy: 3 mg/kg every 3 weeks. During treatment, the dosing frequency of apatinib may be adjusted in response to adverse effects. Possible adjustments include: (4) Apatinib taken every other day; (5) Apatinib taken daily for 5 consecutive days, followed by a 2-day break; (6) Apatinib taken daily for 2 consecutive days, followed by a 2-day break. These modifications aim to mitigate adverse effects while preserving therapeutic efficacy. For other solid tumors, the combination regimen is as follows: (7) Camrelizumab (200 mg every 3 weeks) plus apatinib (250 mg once daily). To explore the feasibility of dose reduction and sequential therapy, we also investigated the following regimens: (8) No treatment; (9) Camrelizumab (200 mg every 2 weeks) plus apatinib (125 mg once daily); (10) Apatinib (250 mg daily for 7 consecutive days), followed by camrelizumab (200 mg every 2 weeks).
In addition to the virtual patients utilized in the CARES-310 clinical study, five additional cohorts, each comprising 500 patients, were generated. These cohorts were employed to simulate the aforementioned dosing regimen in order to evaluate its efficacy compared to the standard treatment plan. Efficacy was assessed every 8 to 12 weeks using RECIST v.1.1 criteria, with the aim of identifying the dosing regimen that offered optimal efficacy and minimal adverse effects (hypertension).
ResultsModel development
Figure 1 illustrates the molecular and cellular interactions that govern the dynamics of various species within the model. The modular structure of the model in this study builds upon the QSP model previously established by Wang et al (30). In addition to incorporating the pharmacokinetic parameters of camrelizumab into the original model, we integrated the pharmacokinetics of apatinib, its tumor-killing effect, and its inhibition of tumor angiogenesis to develop a new apatinib module (31–33).
The dynamics of the major species in the quantitative systems pharmacology model of patients with advanced hepatocellular carcinoma. C, central compartment; P, peripheral compartment; LN, lymph node compartment; T, tumor compartment; kcl, elimination rate constant from central compartment; q_LN, transfer rate constant between central and LN compartments; q_T, between central and tumor compartments; q_P, between central and peripheral compartments; q_LD, between tumor and LN compartments; APC, antigen-presenting cell; mAPC, mature antigen presenting cell; nT, naïve T cell; aT, activated T cell; Teff, effector T cell; Texh, exhausted T cell; Th, T helper cell; Treg, regulatory T cell; MDSC, myeloid-derived suppressor cell; NO, nitric oxide; Arg-I, arginase I; camr, camrelizumab; apat, apatinib.
Complex flowchart depicting immune response pathways in tumor, central, peripheral, and lymphatic system regions, illustrating interactions among cancer cells, immune cells, cytokines, and treatments apatinib and camrelizumab, with symbol legend and pharmacokinetic compartments below.
In the cancer compartment, it depicts processes such as cancer cell proliferation, the release of self- and cancer-associated antigens from dying cancer cells, macrophage recruitment and polarization, and the elimination of cancer cells by Teff. The activation of naïve CD8+ and CD4+ T cells by mature APCs (mAPCs) is shown in the lymph node compartment, followed by their differentiation into Teff and Treg, after which they are transported to both the central and peripheral compartments. Following intravenous administration of camrelizumab, the drug is transported to the cancer compartment, where it enhances the cytotoxic activity of Teff cells against cancer cells by binding to the PD-1 receptor. Concurrently, after oral administration, apatinib is absorbed into the central compartment, from where it is transported to the cancer compartment, where it exerts direct cytotoxic effects on cancer cells. Additionally, apatinib inhibits vascular endothelial growth factor receptor (VEGFR) activity, thereby blocking tumor vasculature growth. This QSP model comprehensively accounts for the impact of these drugs on the tumor cell population, with the corresponding computational equations provided in Equation 1:
Here, C represents the number of cancer cells, and kgrowth is the growth rate constant for cancer cells. Ctotal denotes the total number of cancer cells in the tumor compartment, while Cmax is the maximum carrying capacity for cancer cells in the tumor. kinnate is the death rate constant for cancer cells induced by the innate immune system, and kapat represents the rate constant for the inhibition of cancer cell proliferation by apatinib. The concentration of apatinib in the tumor compartment is denoted as [APAT], and IC50apat refers to the concentration of apatinib required to inhibit 50% of cancer cell proliferation. fapat is the uptake rate of apatinib by cancer cells, and kTcell represents the rate constant for cancer cell death induced by Teff. kC is the maximum killing rate of cancer cells by Teff cells. kTreg is the rate constant for the inhibition of cancer cell proliferation by Treg, and CTreg denotes the number of Treg in the tumor compartment. T represents the number of effector T cells. The model also incorporates the inhibitory effects of PD-1 immune checkpoint (HPD1), TGF-β (HTGF), and MDSC(HMDSC) on T cell activity. kM1 represents the rate constant for the phagocytosis of cancer cells by M1 macrophages, and M1 is the number of M1 macrophages. kMac quantifies the dependence of the phagocytosis rate on the ratio of M1 macrophages to cancer cells. Finally, the inhibitory effects of signal regulatory protein alpha (SIRPA, HMac) and IL-10 (HIL10) on macrophage phagocytosis are included in the model. Cmax is calculated using Equation 2.
kK,g represents the baseline rate of tumor angiogenesis driven by angiogenic factors. Cvas denotes the effective concentration of these factors during the tumor vascularization process. EC50vas is the half-maximal concentration required to elicit a vascular response mediated by angiogenic factors. kK,d characterizes the endogenous inhibition rate of vascular growth, which is regulated by mechanisms such as endothelial cell apoptosis. Vcell corresponds to the volume of an individual aHCC cell. kK,apat quantifies the inhibitory effect of apatinib on tumor angiogenesis.
In the blood pressure dynamics module for apatinib and camrelizumab exposure, the relationship between drug concentration and blood pressure changes is described by Equation 3:
Here, EBPbase represents the baseline blood pressure of the patient before treatment, BP_Emax is the maximum blood pressure increase induced by apatinib, BP_EC50 is The half maximal effective concentration (EC50) of apatinib-induced hypertension, C_Capat is the central compartment concentration of apatinib, BP_Hill is the Hill coefficient for apatinib-induced hypertension, C_Ccamr is the central compartment concentration of camrelizumab, and f_BPcamr is the coefficient describing the synergistic effect of camrelizumab on apatinib-induced blood pressure elevation. Blood pressure of 120–159 mmHg was classified as grade 1/2 hypertension, while ≥160 mmHg was classified as grade 3/4 hypertension. To avoid the occurrence of outliers, the maximum predicted blood pressure value was set to 220 mmHg.
To accurately reflect the distribution of real clinical patient characteristics, the model parameter subsets of the virtual patient queue are varied (refer to Supplementary Table S1). The baseline values and ranges for the population characteristic parameters are sourced from both experimental and clinical data, as outlined in Supplementary Table S2. Additional parameters, species, and the rules governing the model are detailed in the supplementary materials of Wang et al.’s article (34). For certain parameters without reported experimental values (Supplementary Table S2), fitting was performed based on clinical trial results in aHCC patients treated in CARES-310 or with apatinib/camrelizumab monotherapy. In this study, sensitivity analysis was conducted only for variables with assigned variability.
In this study, we developed a virtual patient model comprising 500 patients to simulate tumor treatment responses. During initialization, 10.4% of patients failed to reach the expected baseline tumor size, and 6.4% showed parameter non-convergence. During model execution, we noted that the parameter values for a subset of virtual patients (10 out of 416) exceeded the pre-established reasonable range (refer to Supplementary Table S3). The characteristics of these non-target virtual patients (94/500) are presented in Supplementary Table S4, with t-tests conducted to compare them with the complete dataset.
To further investigate the correlation between modeling parameters and key outcome indicators, such as tumor volume, we conducted a sensitivity analysis using the PRCC (Figure 2). The results indicated a significant positive correlation between tumor growth rate, neoantigen-MHC binding affinity, and tumor angiogenesis rate with tumor volume. Additionally, the tumor’s initial diameter, the recruitment rate of MDSCs, and the level of C-C motif chemokine ligand 2 (CCL2) demonstrated a significant positive correlation with MDSC cell density in the tumor. Notably, a strong negative correlation was identified between tumor growth rate and the cell density of MDSCs and macrophages.
Global uncertainty and sensitivity analysis were conducted using Latin Hypercube Sampling (LHS) to assign thirty-one parameters based on our estimated distribution. The Partial Rank Correlation Coefficient (PRCC) between selected post-treatment observations and input parameters is presented in the form of a heatmap. ***, extremely significant; **, highly significant, *, Significance.
Heatmap illustrating correlation values between various biological parameters (rows) and eight tumor microenvironment features (columns), with color indicating strength and direction of correlation; significance is denoted by asterisks for p-values as defined in the accompanying legend.Model evaluation
In Table 1, we present a comparison of the simulation results derived from our model for 406 reasonable virtual patients against actual data obtained from the CARES-310 clinical trial. The ORR following treatment with standard doses of camrelizumab (200 mg, intravenously administered biweekly) in conjunction with apatinib (250 mg, orally administered once daily) was found to be 25%, while the median DOR was recorded at 14.9 months. The difference between these predictive values and the results of the CARES-310 clinical trial is within 5%, and it closely aligns with the ORR observed in clinical studies of non-first-line treatments using camrelizumab or apatinib monotherapy. These findings validate the accuracy and reliability of our model’s predictions. Additionally, Figure 3 illustrates the time-varying response and the waterfall plot of tumor dynamics, which align with the findings observed in the CARES-310 study, further corroborating the reliability of the QSP model we have established. As shown in Supplementary Table S3, the deviation between our model predictions and the observed hypertension incidence in the CARES-310 clinical trial was 6.14%, demonstrating the model’s robust predictive performance.
Efficacy prediction for the virtual patient cohort generated based on calibrated parameter distribution.
Efficacy measures
Camrelizumab Plus Apatinib
Camrelizumab monotherapy (3 mg/kg q2w)
Apatinib monotherapy (750 mg qd)
Observation (95% CI)
Prediction (95% CI)
Observation (95% CI)*
Prediction (95% CI)
Observation (95% CI)*
Prediction (95% CI)
n
272
406
109
436
261
418
ORR(%)
25.4(20.3 - 31.0)
25.4(21.1 - 29.6)
11.9(6.5 - 19.5)
13.5(10.3 - 16.8)
10.7(7.0 - 15.0)
10.8(7.8 - 13.8)
Median DOR(month)
14.8(8.4 - NR)
14.9(11 - NR)
NE
NE
NE
NE
ORR, objective response rate; DOR, duration of response; NR, not reached; NE, not evaluated; *, second-line or later therapy; q2w, biweekly; qd, Once a day.
Rate of response (A) and the best overall response (B) in model-predicted tumor diameter of 406 randomly virtual patients. Response is assessed by RECIST V.1.1. Median (thick lines) and individual (thin line) rate of response are shown in PD (red), SD (yellow), and PR/CR (blue) subgroups. Dashed lines correspond to +20% and −30%. CR, complete response; PD, progressive disease; PR, partial response; SD, stable disease.
Panel A is a line chart showing individual percent tumor size change over time, color-coded by response category: red for progressive disease (PD), yellow for stable disease (SD), and blue for partial or complete response (PR/CR); dotted horizontal lines separate the response groups. Panel B is a waterfall bar chart with blue bars showing percent tumor size change for individual cases and gray lines indicating standard deviations, with the same response group divisions as Panel A.Model applicationExploration of biomarkers
This study further investigates seven measurable biomarkers through subgroup analysis to identify potential predictors of patient response to treatment. Figure 4 presents the ORR and the corresponding confidence intervals for patient groups stratified by the median values of each biomarker. This analysis identifies several key factors influencing the ORR in combination therapy, including the Teff/Treg ratio, MDSC density, macrophage density, interleukin-10 (IL-10) concentration, interferon-gamma (IFN-γ) concentration, CCL2 concentration, CD8+/Treg ratio, and CD4+/Treg ratio.
Subgroup analysis of the combination therapy. Virtual patients (N = 406) are divided into 14 subgroups based on the pretreatment values of selected biomarkers. Objective response rates (ORR) for each group are given along with the 95% CI estimated. The black solid line indicated the ORR for the total population. Groups were divided by median values; Interaction p value indicates whether two or more variables show significant interaction (p<0.05); initial_tumour_diameter, baseline tumor diameter; PD1 expression, average PD1 per T cell; PDL1 expression, baseline PDL1 per tumor cell; TGFb, transforming growth factor-β; MDSC, myeloid-derived suppressor cell; M, macrophage; IL10, interleukin 10; IFNg, interferon-γ; CCL2, C-C motif chemokine ligand 2.
Forest plot comparing overall response rates (ORR) with 95 percent confidence intervals across multiple immune-related and tumor characteristics subgroups in a cohort of 203 patients, showing interaction P values, with significant values marked by asterisks.
The results of the LASSO regression analysis are depicted in Figure 5A, which illustrates the relationship between the binomial deviance curve and the logarithm of the tuning parameter λ (log(λ)). The optimal λ value was identified through the analysis of the deviance curve, yielding an optimal value of 0.012, which corresponds to log(λ) = -4.44. Utilizing this optimal λ, we subsequently generated the coefficient profile, as presented in Figure 5B and Table 2. LASSO regression analysis identified the CD8+/Treg ratio, CD4+/Treg ratio, MDSC density, CD8+ density, CCL2 concentration, IL-10 concentration, Treg density, and initial tumor diameter as variables of interest. The first six parameters were consistent with the results of the subgroup analysis and are therefore considered biomarkers closely associated with therapeutic efficacy. ROC analysis of these parameters, shown in Figure 6, indicated that the CD8+/Treg ratio had the highest AUC, with a cut-off value of 0.86.
Selection of significant parameters in the clinicopathologic from variables in the training set. (A) Ten-fold cross-validation for tuning parameter selection in the LASSO model. The solid vertical lines represented the binomial deviance ± standard error (SE). The dotted vertical lines were drawn at the optimal values using the minimum and 1-SE criteria. (B) LASSO coefficient profiles. The LASSO shrank some regression coefficients to exactly zero. The binomial deviance curve was plotted versus log (λ), where λ was the tuning parameter. LASSO, Least Absolute Shrinkage and Selection Operator; coef_min, λ giving minimum cross-validation error; coef_1se, largest λ within one standard error of the minimum.
Panel A shows a line chart with red dots representing mean squared error across log-transformed lambda values during model regularization, with error bars and two vertical dashed lines indicating optimal lambda points. Panel B displays a coefficient path plot with multiple colored lines illustrating coefficient estimates as a function of log-transformed lambda, and two dashed vertical lines labeled with log coefficients marking tuning parameter selections.
ROC curve chart comparing six biomarker metrics for prediction performance, with individual colored curves for CD8+/Treg ratio, CD4+/Treg ratio, MDSC density, CD8+ density, CCL2 concentration, and IL-10 concentration; associated AUC and cutoffs for each are provided in the legend. Sensitivity is plotted against one minus specificity.
Using simulation results from 406 virtual patients in the combination therapy group, we developed a random forest model comprising 5000 trees, utilizing 9 biomarkers (derived from LASSO screening results) as predictor variables. The model exhibited robust predictive performance, achieving accuracies of 0.94 on the training set and 0.80 on the test set. Figure 7A presents the Receiver Operating Characteristic (ROC) curve of the random forest classifier, further validating the model’s strong classification capability. Additionally, Figure 7B illustrates the importance of each biomarker within the model. The results indicate that, prior to treatment, the CD4+/Treg ratio and initial tumor diameter are crucial biomarkers closely associated with patient response to treatment.
Exploring biomarkers using a random forest classifier. (A) Receiver operation characteristic curve (ROC) of random forest classifier. The dashed line represents the performance of the random classifier, which is the ideal classifier without classification ability. (B) The importance of features in random forest classifier. AUC, areas under curve; TGFb, transforming growth factor-β; MDSC, myeloid-derived suppressor cell; M, macrophage; IL10, interleukin 10; IFNg, interferon-γ; CCL2, C-C motif chemokine ligand 2; initial_tumour_diameter, baseline tumor diameter.
Panel A shows a receiver operating characteristic (ROC) curve with an area under the curve (AUC) of zero point eight five, indicating good model performance. Panel B is a horizontal bar chart ranking eight features by their importance in the model, with CD4 plus Treg ratio and CD8 plus Treg ratio as the highest, followed by initial tumor diameter, MDSC density, CCL2, CD8 plus density, Treg density, and IL10.Optimization of dosing regimen
To evaluate the potential of the current QSP platform for dose optimization, we designed and conducted ten distinct virtual clinical trials utilizing existing virtual patient datasets. The detailed results of the efficacy evaluations conducted at 8-week and 12-week intervals for each clinical trial group are presented in Table 3. Our analysis considered both ORR and median DOR, revealing that synchronous combination therapy demonstrated a slight advantage over sequential therapy. To optimize the dosing frequency of apatinib, we compared various dosing regimens. Research has demonstrated that, in addition to daily dosing regimens, a continuous daily dosing regimen lasting five days followed by a two-day cessation also exhibits comparable ORR. Additionally, we observed comparable ORR at doses of 250 mg and 125 mg of apatinib, while the dose reduction from 250 mg to 125 mg led to a 25% (0%–56%) decrease in all-grade hypertension and a 27% (0%–61%) decrease in grade ≥3 hypertension (Supplementary Table S6), highlighting its potential to mitigate apatinib-associated toxicity.
Prediction of therapeutic efficacy under different administration strategies.
Medication strategy
Efficacy prediction
Camrelizumab
Apatinib
ORR (%)
Median DOR (month)
*200 mg q2w
250 mg qd
26.0 (23.8-28.1)
14.0 (11.2-15.9)
NA
750 mg qd
11.4 (10.3-12.5)
13.5 (12.1-14.9)
3 mg/kg q3w
NA
13.0 (12.4-13.5)
18.4 (16.8-18.7)
200 mg q2w
250 mg qod
22.2 (19.6-25.2)
13.8 (12.1-14.9)
200 mg q2w
250 mg qd × 5, then 2-day hiatus
24.4 (22.4-28.0)
12.8 (11.2-13.1)
200 mg q2w
250 mg qd × 2, then 2-day hiatus
21.9 (19.2-24.6)
14.8 (14.0-14.9)
200 mg q3w
250 mg qd
25.2 (23.1-27.6)
13.7 (11.2-14.9)
200 mg q2w
125 mg qd
24.4 (21.5-26.8)
13.1 (11.2-14.9)
200 mg q2w began on the seventh day
250 mg qd
25.6 (23.8-27.6)
13.1 (11.2-14.9)
NA
NA
NA
NA
ORR, objective response rate; DOR, duration of response; NA, Not Applicable; *, standard treatment plan; q2w, biweekly; q3w, once every three weeks; qd, once a day; qod, once every other day.
Discussion
In the past decade, although systemic treatment for aHCC has progressed, new therapeutic strategies continue to encounter challenges such as low response rates and drug resistance (35). Utilizing the QSP-IO platform developed by Wang et al. and previously reported pharmacodynamic and pharmacokinetic parameters, we established a QSP model for aHCC patients undergoing the camrelizumab plus apatinib regimen. By generating virtual patients to simulate the CARES-310 trial, the model’s predicted ORR and median DOR were highly consistent with clinical observations. This consistency demonstrates the potential of virtual clinical trials in identifying key biomarkers, optimizing treatment plans, and providing robust support for the screening of advantageous populations and the customization of individualized dosing regimens for camrelizumab combined with apatinib treatment. Compared with a previous QSP model of aHCC patients, our model incorporated more parameters (279 vs. 140), which partly explains the higher rate of non-convergence (6% vs. 3%) (19). Each virtual patient was simulated starting from a single tumor cell, and the initial parameter set required to achieve the preset initial tumor volume was recorded as that of a standard virtual patient. We found no strong correlation between tumor volume and initial tumor diameter, contrary to the findings of Wang et al.’s study (36). Building upon the original QSP model, we incorporated an apatinib module to characterize the pharmacokinetic process of apatinib in the body and its effects on tumor cells. As demonstrated in the sensitivity analysis, the inclusion of the apatinib module did not result in an underestimation of the effect of immunotherapy. This is attributed to the fact that the efficacy of the combination therapy was significantly related not only to the parameters associated with apatinib but also to other immune parameters. Another interesting observation was the negative correlation between tumor growth rate and MDSC/macrophage density in the sensitivity analysis (Figure 2). Mechanistically, although a higher tumor growth rate enhances CCL2 secretion and promotes immune cell recruitment, the rapid expansion of tumor volume outpaces this recruitment, ultimately leading to a decrease in relative immune cell density. This finding underscores the importance of distinguishing between absolute numbers and relative density when analyzing the tumor microenvironment. Constructing virtual cohorts that closely match target populations is essential before conducting virtual clinical trials. However, developing reliable oncology virtual patients remains challenging due to limited clinical data in public databases and the heterogeneity of trial-recruited cohorts (37). To address this, researchers have designed screening methods based on probability density matching or parameter space mapping algorithms to align with clinical data distributions (38, 39). However, the complexity of QSP models and strong parameter correlations can introduce biases or significantly increase computational costs. To balance model applicability and computational efficiency, this study employed LHS for parameter generation and filtered virtual patients (VPs) using predefined feasibility thresholds.
In the biomarker analysis, we conducted a subgroup analysis to investigate significant variables associated with treatment response. For the combination therapy of camrelizumab and apatinib, our model confirmed that the CD8+/Treg ratio, CD4+/Treg ratio, and Teff/Treg ratio were the most predictive biomarkers, as they exhibited a higher ORR in the high-expression subgroup. Studies have demonstrated that the CD8+/Treg ratio in tumors of elderly mice is lower than that in their younger counterparts (40). This observation may elucidate the findings of the CARES-310 clinical study, which indicates that patients under the age of 65 derive greater benefits from combination therapy compared to those aged 65 and older (9). Furthermore, the model identified that the high CD8+ cell subgroup could achieve superior efficacy, aligning with the biomarkers reported for this combination therapy regimen in other solid tumors (41). However, discrepancies were noted between the model-predicted biomarkers and clinical observations. As depicted in Figure 4, the pretreatment PD-L1 expression level was not as closely correlated with efficacy as indicated in other studies, potentially due to Xu et al.’s choice of progression-free survival (PFS) as the efficacy endpoint (26). A similar situation was observed with APC cell concentration, where the disparity may be attributed to differences in immune biology between patients with resectable and unresectable hepatocellular carcinoma (42). In recent years, machine learning algorithms have gained significant traction in clinical decision-making (43). Unlike traditional methods, machine learning does not require the assumption of key drivers to infer or identify unknown associations. However, the ‘black box’ nature of these algorithms often results in outcomes that lack mechanistic explanations (44). The QSP model, grounded in complex mechanisms, addresses this limitation by providing a more comprehensive understanding. Furthermore, the availability of clinical and preclinical data remains limited. The QSP model can generate extensive datasets that are challenging to obtain using conventional technologies. Currently, an increasing number of studies are exploring the integration of QSP models with machine learning techniques (45). For instance, Coletti et al. developed a QSP model for prostate cancer immunotherapy and constructed a decision tree to analyze the results, thereby predicting the potential causality of various treatment regimens (46). In the realm of liver cancer, Lasso regression has emerged as a commonly employed algorithm, particularly for predicting patients’ recurrence risks and identifying biomarkers related to drug treatment efficacy (47).
To further identify individuals who may benefit from tailored treatment plans and to develop personalized diagnostic and therapeutic strategies for advanced liver cancer patients, we conducted LASSO regression analysis on 17 parameters of virtual patients prior to treatment administration. This analysis revealed 8 biomarkers closely associated with treatment response status, including the CD8+/Treg ratio and CCL2 concentration in tumors (see Figure 5, Table 2). Research indicates that CCL2 is overexpressed in human liver cancer and is closely linked to patient prognosis (48). Previous work by Xia et al. demonstrated that combining PD-L1 blockade with FGFR4 or MAPK inhibition significantly reduces CCL2 levels and suppresses metastatic progression (49). Informed by these findings and our QSP modeling results, we hypothesize that incorporating CCL2/CCR2 blockade into camrelizumab plus apatinib therapy may alleviate myeloid-driven immunosuppression and enhance antitumor immune activity, thereby providing a mechanistically grounded strategy for personalized combination treatment in aHCC. Subsequently, the variables selected through LASSO regression were utilized to construct a random forest model. Typically, an area under the curve (AUC) of ≥ 0.7 signifies sufficient predictive capability (50). Our findings suggest that the random forest model established with the 8 biomarkers demonstrates strong predictive performance for the test set (AUC = 0.85). Notably, our research identified that the CD8+/Treg ratio, CD4+/Treg ratio, and initial tumor diameter significantly contribute to the random forest model. The six biomarkers identified through the combination of subgroup analysis and machine learning in this study show promising clinical applicability. In 2025, the European Society for Medical Oncology (ESMO) guidelines recommended histopathological confirmation of aHCC prior to initiating systemic therapy (51). This recommendation provides an excellent opportunity for biomarker testing, as residual biopsy samples can be quantitatively analyzed using immunohistochemistry, flow cytometry, and other techniques (51–53). Integrating the quantitative cut-off values identified in this study into existing clinical workflows may offer a direct and feasible reference for developing individualized treatment strategies, thereby facilitating the clinical translation of these biomarkers. Biomarker exploration using machine learning on QSP model outputs in aHCC has been reported previously, but unlike that study, our random forest model did not include PD-L1 expression (19). Our model, based on eight biomarkers, achieved an AUC of 0.85, indicating strong predictive performance and potential clinical utility. When PD-L1 expression was incorporated into the model, the AUC decreased to 0.79. This difference reflects variations in drug mechanisms and calibration trials. Notably, in the CARES-310 trial, PD-L1 was also not significantly associated with combination therapy efficacy, supporting our findings. Rather than acting as a determinant of response, PD-L1 appears to reflect adaptive immune evasion; its expression is both temporally dynamic and spatially heterogeneous (54, 55). In hepatocellular carcinoma, multiple concurrent immunosuppressive circuits—including MDSCs, tumor-associated macrophages, chemokine signaling, and angiogenic pathways—can attenuate the effects of PD-1/PD-L1 blockade, thereby limiting the predictive utility of PD-L1 as a standalone biomarker (56). Consistently, QSP-IO simulations indicate that PD-L1 levels do not capture overall immune competence, whereas mechanistic indices integrating effector and suppressive components (such as CD8+/Treg or CD4+/Treg ratios) better account for inter-patient heterogeneity in response to combination therapy. Our model achieved 80% accuracy in the test set, slightly lower than the 91.3% reported previously, likely due to differences in sample sizes (training set: n=325 vs. n=1365; test set: n=81 vs. n=50) (19). Nevertheless, the clinical utility of these biomarkers requires confirmation in large-scale real-world cohorts. Future studies will focus on collecting patient data from our center to prospectively or retrospectively evaluate their predictive performance and robustness, thereby providing a foundation for their integration into clinical decision-making.
During the dose simulation process, we explored not only the existing clinical trial regimen but also the feasibility of reducing the dose or frequency of administration. The simulated 750 mg apatinib monotherapy yielded an ORR consistent with that observed in the AHELP clinical trial (25). Furthermore, the ORR of camrelizumab monotherapy in the second-line treatment clinical trial (11.9% vs. 13.5%) closely aligned with our simulation results (24). Reducing apatinib to 125 mg achieved efficacy comparable to the standard dose (ORR 25.9% vs. 24.4%; median DOR 14.0 vs. 13.0 months) while lowering the incidence of hypertension by more than 25%. Prior evidence suggests that anti-angiogenic treatment at an appropriate dose range is sufficient to achieve vascular normalization and facilitate immune infiltration (57). Increasing the dose beyond this range may fail to enhance immunologic synergy and instead intensify hypoxia and toxicity. In contrast, lower dosing has been proposed to extend the duration of vascular normalization, thereby augmenting the therapeutic benefit of immunotherapy (58). Consistently, QSP simulations suggest that the biologically effective dose of apatinib required for immune modulation in combination regimens is below currently used levels, supporting dose de-escalation as a strategy to maintain efficacy while improving tolerability.
The QSP model developed in this study demonstrated promising results in identifying biomarkers and optimizing doses for the camrelizumab and apatinib combination regimen. However, it has certain limitations. Firstly, this study used tumor size dynamics in aHCC patients as the core efficacy output, without further evaluation of progression-free or overall survival. Moreover, HBV/HCV status and Child-Pugh classification were not incorporated into the model framework, limiting subgroup analyses. Secondly, we simplified apatinib’s anti-angiogenic effect, which reduces the maximum tumor carrying capacity, into a single parameter (kK_apat). While this provides a clearer representation of its inhibitory effect on tumor cell proliferation, it reduces mechanistic interpretability. Future studies will incorporate detailed signaling pathways to better evaluate efficacy, adverse effects, and potential resistance mechanisms (59). Thirdly, due to limited computational resources and the lack of systematically integrated spatial multi-omics data for calibration and validation, intratumoral spatial heterogeneity was not incorporated. Future improvements could leverage single-cell sequencing and spatial transcriptomics to better capture the complexity of the tumor microenvironment. Emerging evidence indicates that, compared with the tumor core, macrophage recruitment and M2 polarization are more pronounced at the tumor margin, leading to localized immunosuppression and accelerated progression of hepatocellular carcinoma (60). These spatially distinct immune features suggest that future QSP models for advanced HCC may benefit from explicitly separating the tumor into core and margin compartments, thereby better capturing the invasive behavior and therapeutic resistance observed in clinical practice. Finally, during the VP generation process, we can utilize the machine learning surrogate modeling approach proposed by Myers et al. to pre-screen parameter combinations while maintaining full model validation, significantly improving the efficiency of VP generation (61).
Conclusions
This study expands the functionality of the QSP-IO platform to make it applicable to patients with aHCC receiving combination therapy with camrelizumab and apatinib. The model was calibrated and validated using real-world clinical trial data. A virtual clinical trial was conducted based on the established QSP model to simulate the therapeutic outcomes of camrelizumab and apatinib administration. For biomarker identification, we combined subgroup analysis with machine learning approaches, and identified the CD8+/Treg ratio, CD4+/Treg ratio, and the density of MDSCs as potential predictive biomarkers for treatment efficacy. Additionally, the feasibility of reducing the apatinib dose from 250 mg to 125 mg was evaluated, providing a reference for future dose optimization in clinical practice. Overall, this study highlights the significant potential of QSP in biomarker discovery and optimization of therapeutic strategies.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.
Author contributions
WH: Writing – review & editing, Data curation, Writing – original draft, Methodology, Visualization, Validation. GT: Writing – review & editing, Data curation, Writing – original draft, Supervision, Resources. DL: Visualization, Writing – original draft, Formal Analysis, Writing – review & editing. CW: Writing – original draft, Methodology, Writing – review & editing, Supervision. JZ: Writing – review & editing, Writing – original draft. ZJ: Methodology, Software, Writing – review & editing, Writing – original draft, Data curation, Validation. LY: Funding acquisition, Supervision, Writing – review & editing, Writing – original draft, Investigation, Project administration.
Acknowledgments
We thank Figdraw (www.figdraw.com) for generation of figures in this article.
Conflict of interest
The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2026.1617227/full#supplementary-material
ReferencesBrayFLaversanneMSungHFerlayJSiegelRLSoerjomataramI.
Global cancer statistics 2022: globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. (2024) 74:229–63. doi: 10.3322/caac.21834, PMID: 38572751LlovetJMKelleyRKVillanuevaASingalAGPikarskyERoayaieS.
Hepatocellular carcinoma. Nat Rev Dis Primers. (2021) 7:6. doi: 10.1038/s41572-020-00240-3, PMID: 33479224CelsaCGiuffridaPStornelloCGrovaMSpatolaFRizzoGEM.
Systemic therapies for hepatocellular carcinoma: the present and the future. Recenti Prog Med. (2021) 112:110–6. doi: 10.1701/3559.35371, PMID: 33624623El DikaIMakkiIAbou-AlfaGK.
Hepatocellular carcinoma, novel therapies on the horizon. Chin Clin Oncol. (2021) 10:12. doi: 10.21037/cco-20-113, PMID: 32527116SangroBSarobePHervas-StubbsSMeleroI.
Advances in immunotherapy for hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol. (2021) 18:525–43. doi: 10.1038/s41575-021-00438-0, PMID: 33850328ChenYHuHYuanXFanXZhangC.
Advances in immune checkpoint inhibitors for advanced hepatocellular carcinoma. Front Immunol. (2022) 13:896752. doi: 10.3389/fimmu.2022.896752, PMID: 35757756VogelAMeyerTSapisochinGSalemRSaborowskiA.
Hepatocellular carcinoma. Lancet. (2022) 400:1345–62. doi: 10.1016/S0140-6736(22)01200-4, PMID: 36084663XuJShenJGuSZhangYWuLWuJ.
Camrelizumab in combination with apatinib in patients with advanced hepatocellular carcinoma (Rescue): A nonrandomized, open-label, phase ii trial. Clin Cancer Res. (2021) 27:1003–11. doi: 10.1158/1078-0432.CCR-20-2571, PMID: 33087333QinSChanSLGuSBaiYRenZLinX.
Camrelizumab plus rivoceranib versus sorafenib as first-line therapy for unresectable hepatocellular carcinoma (Cares-310): A randomised, open-label, international phase 3 study. Lancet. (2023) 402:1133–46. doi: 10.1016/S0140-6736(23)00961-3, PMID: 37499670LiuXLuYZhouWPengTZhouJBiH.
Chinese multidisciplinary expert consensus on immune checkpoint inhibitor-based combination therapy for hepatocellular carcinoma (2023 edition). Liver Cancer. (2024) 13:355–75. doi: 10.1159/000535496, PMID: 39114757ZengZJiangYLiuCZhuGMaFYangL.
Efficacy and biomarker exploration of camrelizumab combined with apatinib in the treatment of advanced primary liver cancer: A retrospective study. Anti-cancer Drugs. (2021) 32:1093–8. doi: 10.1097/CAD.0000000000001127, PMID: 34232941LiangLWenYHuRWangLXiaYHuC.
Safety and efficacy of pd-1 blockade-activated multiple antigen-specific cellular therapy alone or in combination with apatinib in patients with advanced solid tumors: A pooled analysis of two prospective trials. Cancer Immunol Immunother. (2019) 68:1467–77. doi: 10.1007/s00262-019-02375-z, PMID: 31451841AzerKKaddiCDBarrettJSBaiJPFMcQuadeSTMerrillNJ.
History and future perspectives on the discipline of quantitative systems pharmacology modeling and its applications. Front Physiol. (2021) 12:637999. doi: 10.3389/fphys.2021.637999, PMID: 33841175WangHArulrajTIppolitoAPopelAS.
From virtual patients to digital twins in immuno-oncology: lessons learned from mechanistic quantitative systems pharmacology modeling. NPJ Digit Med. (2024) 7:189. doi: 10.1038/s41746-024-01188-4, PMID: 39014005LemaireVBassenDReedMSongRKhaliliSLienYTK.
From cold to hot: changing perceptions and future opportunities for quantitative systems pharmacology modeling in cancer immunotherapy. Clin Pharmacol Ther. (2023) 113:963–72. doi: 10.1002/cpt.2770, PMID: 36282521SoveRJJafarnejadMZhaoCWangHMaHPopelAS.
Qsp-io: A quantitative systems pharmacology toolbox for mechanistic multiscale modeling for immuno-oncology applications. CPT Pharmacometr Syst Pharmacol. (2020) 9:484–97. doi: 10.1002/psp4.12546, PMID: 32618119ArulrajTWangHEmensLASanta-MariaCAPopelAS.
A transcriptome-informed qsp model of metastatic triple-negative breast cancer identifies predictive biomarkers for pd-1 inhibition. Sci Adv. (2023) 9:eadg0289. doi: 10.1126/sciadv.adg0289, PMID: 37390206WangCYDaiHRTanYPYangDHNiuXMHanL.
Development and evaluation of a quantitative systems pharmacology model for mechanism interpretation and efficacy prediction of atezolizumab in combination with carboplatin and nab-paclitaxel in patients with non-small-cell lung cancer. Pharm (Basel). (2024) 17:238. doi: 10.3390/ph17020238, PMID: 38399453SoveRJVermaBKWangHHoWJYarchoanMPopelAS.
Virtual clinical trials of anti-pd-1 and anti-ctla-4 immunotherapy in advanced hepatocellular carcinoma using a quantitative systems pharmacology model. J Immunother Cancer. (2022) 10:e005414. doi: 10.1136/jitc-2022-005414, PMID: 36323435WangHMaHSoveRJEmensLAPopelAS.
Quantitative systems pharmacology model predictions for efficacy of atezolizumab and nab-paclitaxel in triple-negative breast cancer. J Immunother Cancer. (2021) 9:e002100. doi: 10.1136/jitc-2020-002100, PMID: 33579739MarinoSHogueIBRayCJKirschnerDE.
A methodology for performing global uncertainty and sensitivity analysis in systems biology. J Theor Biol. (2008) 254:178–96. doi: 10.1016/j.jtbi.2008.04.011, PMID: 18572196ZhouYTChuJHZhaoSHLiGLFuZYZhangSJ.
Quantitative systems pharmacology modeling of her2-positive metastatic breast cancer for translational efficacy evaluation and combination assessment across therapeutic modalities. Acta Pharmacol Sin. (2024) 45:1287–304. doi: 10.1038/s41401-024-01232-9, PMID: 38360930EisenhauerEATherassePBogaertsJSchwartzLHSargentDFordR.
New response evaluation criteria in solid tumours: revised recist guideline (Version 1.1). Eur J Cancer (Oxford England: 1990). (2009) 45:228–47. doi: 10.1016/j.ejca.2008.10.026, PMID: 19097774QinSRenZMengZChenZChaiXXiongJ.
Camrelizumab in patients with previously treated advanced hepatocellular carcinoma: A multicentre, open-label, parallel-group, randomised, phase 2 trial. Lancet Oncol. (2020) 23:571–80. doi: 10.1016/S1470-2045(20)30011-5, PMID: 32112738QinSLiQGuSChenXLinLWangZ.
Apatinib as second-line or later therapy in patients with advanced hepatocellular carcinoma (Ahelp): A multicentre, double-blind, randomised, placebo-controlled, phase 3 trial. Lancet Gastroenterol Hepatol. (2021) 6:559–68. doi: 10.1016/S2468-1253(21)00109-6, PMID: 33971141XuJZhangYJiaRYueCChangLLiuR.
Anti-pd-1 antibody shr-1210 combined with apatinib for advanced hepatocellular carcinoma, gastric, or esophagogastric junction cancer: an open-label, dose escalation and expansion study. Clin Cancer Res. (2019) 25:515–23. doi: 10.1158/1078-0432.CCR-18-2484, PMID: 30348638JuWTXiaRHZhuDWDouSJZhuGPDongMJ.
A pilot study of neoadjuvant combination of anti-pd-1 camrelizumab and vegfr2 inhibitor apatinib for locally advanced resectable oral squamous cell carcinoma. Nat Commun. (2022) 13:5378. doi: 10.1038/s41467-022-33080-8, PMID: 36104359TibshiraniR.
Regression Shrinkage and Selection Via the Lasso. J R Stat Soc Ser B: Stat Methodol. (1996) 58:267–88. doi: 10.1111/j.2517-6161.1996.tb02080.xIshwaranHKogalurUBBlackstoneEHLauerMS.
Random survival forests. Ann Appl Stat. (2008) 2:841–60. doi: 10.1214/08-AOAS169WangHArulrajTKimkoHPopelAS.
Generating immunogenomic data-guided virtual patients using a qsp model to predict response of advanced nsclc to pd-L1 inhibition. NPJ Precis Oncol. (2023) 7:55. doi: 10.1038/s41698-023-00405-9, PMID: 37291190YangCQinS.
Apatinib targets both tumor and endothelial cells in hepatocellular carcinoma. Cancer Med. (2018) 7:4570–83. doi: 10.1002/cam4.1664, PMID: 30109780YangXMYePPLiuXLGuoZXKanMLiQ.
Population pharmacokinetics and pharmacogenetics of apatinib in adult cancer patients. Br J Clin Pharmacol. (2023) 89:1862–72. doi: 10.1111/bcp.15665, PMID: 36662574WangCYShengCCMaGLXuDLiuXQWangYY.
Population pharmacokinetics of the anti-pd-1 antibody camrelizumab in patients with multiple tumor types and model-informed dosing strategy. Acta Pharmacol Sin. (2021) 42:1368–75. doi: 10.1038/s41401-020-00550-y, PMID: 33154554WangHZhaoCSanta-MariaCAEmensLAPopelAS.
Dynamics of tumor-associated macrophages in a quantitative systems pharmacology model of immunotherapy in triple-negative breast cancer. iScience. (2022) 25:104702. doi: 10.1016/j.isci.2022.104702, PMID: 35856032YangXYangCZhangSGengHZhuAXBernardsR.
Precision treatment in advanced hepatocellular carcinoma. Cancer Cell. (2024) 42:180–97. doi: 10.1016/j.ccell.2024.01.007, PMID: 38350421WangHSoveRJJafarnejadMRahmehSJaffeeEMStearnsV.
Conducting a virtual clinical trial in her2-negative breast cancer using a quantitative systems pharmacology model with an epigenetic modulator and immune checkpoint inhibitors. Front Bioeng Biotechnol. (2020) 8:141. doi: 10.3389/fbioe.2020.00141, PMID: 32158754WangHArulrajTIppolitoAPopelAS.
Quantitative systems pharmacology modeling in immuno-oncology: hypothesis testing, dose optimization, and efficacy prediction. Handb Exp Pharmacol. (2024) 289:261–84. doi: 10.1007/164_2024_735, PMID: 39707022AllenRJRiegerTRMusanteCJ.
Efficient generation and selection of virtual populations in quantitative systems pharmacology models. CPT Pharmacometr Syst Pharmacol. (2016) 5:140–6. doi: 10.1002/psp4.12063, PMID: 27069777DerippeTFouliardSDeclevesXMagerDE.
Accelerating robust plausible virtual patient cohort generation by substituting ode simulations with parameter space mapping. J Pharmacokinet Pharmacodyn. (2022) 49:625–44. doi: 10.1007/s10928-022-09826-8, PMID: 36316606SitnikovaSIWalkerJAPrickettLBMorrowMValge-ArcherVERobinsonMJ.
Age-induced changes in anti-tumor immunity alter the tumor immune infiltrate and impact response to immuno-oncology treatments. Front Immunol. (2023) 14:1258291. doi: 10.3389/fimmu.2023.1258291, PMID: 37920465LiuJLiYLiQLiangDWangQLiuQ.
Biomarkers of response to camrelizumab combined with apatinib: an analysis from a phase ii trial in advanced triple-negative breast cancer patients. Breast Cancer Res Treat. (2021) 186:687–97. doi: 10.1007/s10549-021-06128-4, PMID: 33634417XiaYTangWQianXLiXChengFWangK.
Efficacy and safety of camrelizumab plus apatinib during the perioperative period in resectable hepatocellular carcinoma: A single-arm, open label, phase ii clinical trial. J Immunother Cancer. (2022) 10:e004656. doi: 10.1136/jitc-2022-004656, PMID: 35379737WangQQiaoWZhangHLiuBLiJZangC.
Nomogram established on account of lasso-cox regression for predicting recurrence in patients with early-stage hepatocellular carcinoma. Front Immunol. (2022) 13:1019638. doi: 10.3389/fimmu.2022.1019638, PMID: 36505501AghamiriSSAminRHelikarT.
Recent applications of quantitative systems pharmacology and machine learning models across diseases. J Pharmacokinet Pharmacodyn. (2022) 49:19–37. doi: 10.1007/s10928-021-09790-9, PMID: 34671863ZhangTAndroulakisIPBonatePChengLHelikarTParikhJ.
Two heads are better than one: current landscape of integrating qsp and machine learning: an isop qsp sig white paper by the working group on the integration of quantitative systems pharmacology and machine learning. J Pharmacokinet Pharmacodyn. (2022) 49:5–18. doi: 10.1007/s10928-022-09805-z, PMID: 35103884ColettiRLeonardelliLParoloSMarchettiL.
A qsp model of prostate cancer immunotherapy to identify effective combination therapies. Sci Rep. (2020) 10:9063. doi: 10.1038/s41598-020-65590-0, PMID: 32493951ChenDLiuJZangLXiaoTZhangXLiZ.
Integrated machine learning and bioinformatic analyses constructed a novel stemness-related classifier to predict prognosis and immunotherapy responses for hepatocellular carcinoma patients. Int J Biol Sci. (2022) 18:360–73. doi: 10.7150/ijbs.66913, PMID: 34975338LiDZhangTGuoYBiCLiuMWangG.
Biological impact and therapeutic implication of tumor-associated macrophages in hepatocellular carcinoma. Cell Death Dis. (2024) 15:498. doi: 10.1038/s41419-024-06888-z, PMID: 38997297XieMLinZJiXLuoXZhangZSunM.
Fgf19/fgfr4-mediated elevation of etv4 facilitates hepatocellular carcinoma metastasis by upregulating pd-L1 and ccl2. J Hepatol. (2023) 79:109–25. doi: 10.1016/j.jhep.2023.02.036, PMID: 36907560DongBZhangHDuanYYaoSChenYZhangC.
Development of a machine learning-based model to predict prognosis of alpha-fetoprotein-positive hepatocellular carcinoma. J Trans Med. (2024) 22:455. doi: 10.1186/s12967-024-05203-w, PMID: 38741163LiMWangLCongLWongCCZhangXChenH.
Spatial proteomics of immune microenvironment in nonalcoholic steatohepatitis-associated hepatocellular carcinoma. Hepatology. (2024) 79:560–74. doi: 10.1097/HEP.0000000000000591, PMID: 37733002TanJFanWLiuTZhuBLiuYWangS.
Trem2(+) macrophages suppress cd8(+) T-cell infiltration after transarterial chemoembolisation in hepatocellular carcinoma. J Hepatol. (2023) 79:126–40. doi: 10.1016/j.jhep.2023.02.032, PMID: 36889359MagenAHamonPFiaschiNSoongBYParkMDMattiuzR.
Intratumoral dendritic cell-cd4(+) T helper cell niches enable cd8(+) T cell differentiation following pd-1 blockade in hepatocellular carcinoma. Nat Med. (2023) 29:1389–99. doi: 10.1038/s41591-023-02345-0, PMID: 37322116WangYZhouYYangLLeiLHeBCaoJ.
Challenges coexist with opportunities: spatial heterogeneity expression of pd-L1 in cancer therapy. Adv Sci (Weinh). (2024) 11:e2303175. doi: 10.1002/advs.202303175, PMID: 37934012TopalianSLTaubeJMAndersRAPardollDM.
Mechanism-driven biomarkers to guide immune checkpoint blockade in cancer therapy. Nat Rev Cancer. (2016) 16:275–87. doi: 10.1038/nrc.2016.36, PMID: 27079802FuYGuoXSunLCuiTWuCWangJ.
Exploring the role of the immune microenvironment in hepatocellular carcinoma: implications for immunotherapy and drug resistance. Elife. (2024) 13:e95009. doi: 10.7554/eLife.95009, PMID: 39146202Dianat-MoghadamHNedaeiniaRKeshavarzMAziziMKazemiMSalehiR.
Immunotherapies targeting tumor vasculature: challenges and opportunities. Front Immunol. (2023) 14:1226360. doi: 10.3389/fimmu.2023.1226360, PMID: 37727791FukumuraDKloepperJAmoozgarZDudaDGJainRK.
Enhancing cancer immunotherapy using antiangiogenics: opportunities and challenges. Nat Rev Clin Oncol. (2018) 15:325–40. doi: 10.1038/nrclinonc.2018.29, PMID: 29508855ZhangYWangHOliveiraRHMZhaoCPopelAS.
Systems biology of angiogenesis signaling: computational models and omics. WIREs Mech Dis. (2022) 14:e1550. doi: 10.1002/wsbm.1550, PMID: 34970866WuLYanJBaiYChenFZouXXuJ.
An invasive zone in human liver cancer identified by stereo-seq promotes hepatocyte-tumor cell crosstalk, local immunosuppression and tumor progression. Cell Res. (2023) 33:585–603. doi: 10.1038/s41422-023-00831-1, PMID: 37337030MyersRCAugustinFHuardJFriedrichCM.
Using machine learning surrogate modeling for faster qsp vp cohort generation. CPT Pharmacometr Syst Pharmacol. (2023) 12:1047–59. doi: 10.1002/psp4.12999, PMID: 37328956
Edited by: M. Cristina Vega, Spanish National Research Council (CSIC), Spain
Reviewed by: Xiangliang Liu, The First Hospital of Jilin University, China