All mouse studies were approved by the Pfizer Institutional Animal Care and Use Committee according to established guidelines.

To obtain intra-tumoral T cell profiling data, four immunocompetent Balb/c mice were inoculated with CT26 tumor cells. The tumors were harvested 7 days after inoculation. Frequencies of CD8⁺ inactive T cells [CD8⁺CD4-TCF7+ (Kurtulus et al., 2019)], CD8⁺ active T cells [CD8⁺CD4-PD1+CD44⁺ (Baaten et al., 2010)], CD8⁺ proliferating T cells [CD8⁺CD4-Ki67+ (Scholzen and Gerdes, 2000)], CD8⁺ cytotoxic T-lymphocytes [CTLs, CD8⁺CD4-Granzyme B+ (Pardo et al., 2009)] and CD8⁺ exhausted CTLs [CD8⁺CD4-PD1+Tim3+Lag3+ (Maruhashi et al., 2020; Yang et al., 2020)] out of CD45⁺ T cells were quantified by cell surface markers using flow cytometry.

Tumor growth inhibition data are from five studies of immunocompetent Balb/c mice inoculated with CT26 tumor cells. The mice were treated with phosphate buffer saline (PBS, control), 0.5% methylcellulose (control), water (control), anti-mouse-PD-1 antibody (CD279, RMP1-14, BioXCell®), or anti-mouse PD-L1 antibody (B7-H1, 10F.9G2, BioXCell®). Tumor sizes were measured using calipers. Tumor volumes were calculated as (width x width x length)/2. Every study had a control group and 1-2 anti-PD-(L)1 treatment groups. Details of the five studies are shown in Table 1.

A QSP model was originally designed to investigate sources of variability in response to anti-CTLA4 (Qiao, 2022) and the model was updated to incorporate the mechanism of action of anti-PD-(L)1. The QSP model consists of three compartments, including the plasma and peripheral compartments to capture the pharmacokinetics of anti-PD-(L)1 antibodies, and the tumor compartment to describe the pharmacodynamics within the tumor microenvironment.

In the tumor compartment the concentrations of tumor cells, PD-1+CD8⁺ T cells, and other PD-1+ T cells were modeled. Tumor cells are damaged through interactions with CTLs or through anti-PD-L1 ADCC. Damaged tumor cells undergo apoptosis at a constant rate. Other PD-1+ T cells serve as an alternate binding site for PD-L1 on tumor cells. PD-1+CD8⁺ T cells include CD8⁺ T cells in the inactive, active, proliferating, CTL, and exhausted CTL stages depicted in Figure 1A. PD-1+CD8⁺ T cells are synthesized as inactive cells ( E I ). The cells go through eight proliferation stages ( E A , and E P i for i = 1 . . 7 ) before becoming CTLs ( E D i for i = 1 . . . 10 ) (Qiao, 2022). A total of eight proliferation stages was chosen to match the magnitude of T cell was chosen to match the magnitude of T cell expansion observed in (Yoon et al., 2010). Each CTL can damage ten tumor cells before becoming exhausted ( E D 0 ) (Halle et al., 2016). All CTLs are assumed to be reactive towards the tumor. All PD-1+CD8⁺ T cells express the same number of PD-1 receptors per cell. Inactive, active, and proliferating CD8⁺ T cells degrade at a rate of k d , C D 8 while CTLs degrade at a rate of k d , C T L .

The binding of PD-L1 expressed on tumor cells to PD-1 expressed on inactive PD-1+CD8⁺ T cells ( P I : L V ) controls T cell behavior. This binding inhibits the influx ( f i n f l u x of inactive PD-1+CD8⁺ T cells into the tumor microenvironment, and the activation ( f a c t i v a t e ) of PD-1+CD8⁺ T cells in the tumor compartment (Eqs 1–3). In the absence of antibody treatments, PD-1:PD-L1 per inactive CD8⁺ T cell is at its steady state concentration (approximately 6.6 PD-1:PD-L1 per inactive CD8⁺ T cell) causing inactive PD-1+CD8⁺ T cells to enter the tumor compartment at the rate of k s , I V t (Eq. 2). The influx rate ensures that the CD8⁺ T cells reach a steady state concentration during control simulations. The maximum influx and activation rates are denoted E max i n f l u x and E max a c t i v a t e , respectively. Both functions reach half maximal values at E C 50 . d E I d t = f i n f l u x P I : L V E I − k d , C D 8 E I − f a c t i v a t e P I : L V E I E I + k A 2 I b a s a l E A (1) f i n f l u x P I : L V E I = max 0 , E max i n f l u x 1 + E C 50 6.6 − E max i n f l u x 1 + E C 50 P I : L V E I + k s , I V t (2) f a c t i v a t e P I : L V E I = E max a c t i v a t e − E max a c t i v a t e 1 + E C 50 P I : L V E I (3)

Equations describing the subsequent T cell stages follow the set up in (Qiao, 2022) and can be found with the full model equations in Supplementary Material.

Parameters describing receptor expression, antibody half-lives, binding constants, and plasma concentrations were sourced from relevant literature (Supplementary Table S1). The unknown model parameters were estimated in three stages. First, pharmacokinetic parameters of RMP1-14 (anti-PD-1) and 10F.9G2 (anti-PD-L1) were fit from the published single dose data (Zalba et al., 2020) using two-compartment models with linear elimination (Supplementary Figure S1).

Second, T cell kinetic parameters including maximum CD8⁺ T cell activation rate, T cell proliferation half-life, CD8⁺ T cell half-life, CTL half-life, CTL-mediated tumor cell killing rate, and baseline influx rate of inactive CD8⁺ T cells were constrained using average concentrations of CD8⁺ T cells from untreated CT26 tumors in Balb/c mice that were derived from flow cytometry data. Full details can be found in Supplementary Material.

Third, the full QSP model was fit to tumor volume profiles from five in vivo mouse studies. Non-linear mixed effects (NLME) model fitting and local sensitivity analysis were performed iteratively to identify the key parameters to capture inter-individual variability. Local sensitivity analysis was centered on the population parameter estimates. In each iteration, sensitive parameters were included as random effects in the NLME model. Based on the results of fitting the inter-individual variability, study site ID of the five in vivo mouse studies was identified as a significant covariate. An exponential covariate model was included to capture study site level effects. Goodness of fit was evaluated using Akaike information criteria, Bayesian information criteria and visual inspection of individual fits. The control and treatment groups of all the studies were fit simultaneously.

Model fitting was performed in Monolix 2020R1 (Lixoft SAS, a Simulations Plus company). Local sensitivity analysis was performed in MATLAB version R2020b (MathWorks software). Plots were generated in R version 4.0.3.

Local sensitivity analysis was performed using either the population fit or a chosen fitted animal as the starting parameter set. Parameters were perturbed ± 30% one at a time. The percent change in the area under the tumor volume curve (AUC) was calculated for control, anti-PD-1, and anti-PD-L1 treatment simulations. Three doses of 10 mg/kg given every 3 days (Q3D) were simulated with the first dose administered on day 7. The AUC was chosen to summarize the changes to tumor volume over the entire 30 days.

Local sensitivity analysis was performed in MATLAB version R2020b (MathWorks software). Plots were generated in R version 4.0.3.

The QSP model describes the pharmacokinetics of anti-PD-(L)1, interactions between the antibodies and their molecular targets, tumor cell dynamics, intra-tumoral T cell dynamics, and CTL-mediated tumor cell killing (Figure 1A). Five stages of T cell development were tracked in the tumor compartment, described as inactive, active, proliferating, cytotoxic, and exhausted T cells. The concentration of each of these T cell populations depends on the T cell kinetic parameters including baseline CD8⁺ T cell influx rate, CD8⁺ T cell half-life, CTL half-life, maximum CD8⁺ T cell activation rate, CD8⁺ T cell proliferation half-life, and CTL-mediated tumor cell killing rate. These parameters were constrained using the average concentrations of CD8⁺ T cells from untreated CT26 tumors in Balb/c mice (Figure 1B, Supplementary Table S1) that were derived from flow cytometry data where the five CD8⁺ T cell populations were defined by cell surface markers (Supplementary Material).

The maximum CD8⁺ T cell activation rate limits the maximum amount of change in the CTL concentration post antibody treatments. Specifically, simulations indicate that anti-PD-1 treatments administered Q3Dx3 induced a maximum increase of 2.8% in CTL concentration and 0.017% reduction in tumor volume 10 days after treatment initiation compared to control simulations, whereas anti-PD-L1 treatment administered Q3Dx3 induced a maximum increase of 56% in CTL concentration and 52% reduction in tumor volume (Figures 1C,D). Differences in response to anti-PD-1 and anti-PD-L1 treatments is due to differences in receptor occupancy and the ADCC effect of anti-PD-L1. These simulations indicated that activation and expansion of resident CD8⁺ T cells in the tumor microenvironment alone are insufficient for modeling the anti-tumor responses of anti-PD-(L)1 treatments.

When a PD-1:PD-L1 interaction dependent CD8⁺ T cell influx was incorporated, the increase in CTL concentration due to anti-PD-(L)1 treatments significantly increased (from 2.8% to 65000% for anti-PD-1, and from 56% to 208000% for anti-PD-L1), as did the reduction in tumor volume (from 0.017% to 35% for anti-PD-1, and from 52% to 75% for anti-PD-L1) (Figures 1C,D).

In summary, intratumor T cell profiling data were used to characterize kinetic parameters of CD8⁺ T cells in CT26 + tumor microenvironment. Model simulations indicate that dependence of the influx of CD8⁺ T cells to the tumor microenvironment on the amount of PD-1:PD-L1 interactions between T cells and tumor cells is necessary for anti-PD-(L)1-mediated anti-tumor responses.

Significant variability in tumor volume profiles and responses to anti-PD-(L)1 treatments was observed in the in vivo studies with the CT26 syngeneic mouse tumor model. The dataset consisted of 136 mice with CT26 tumors treated with PBS, 0.5% MC or water (control group), anti-PD-1 antibody (RMP1-14) or anti-PD-L1 antibody (10F.9G2) (treatment details in Table 1). Representative data are shown in Figure 2A and the full dataset is shown in Supplementary Figure S4.

To quantify the differences in individual tumor volume profiles, the observed rates of tumor volume change between day 7 and day 13 after treatment start were calculated for each mouse. The observed rates of tumor volume change varied significantly between studies, as identified through a Kruskal-Wallis rank sum test (p = 0.0054).

The average observed rates of tumor volume change between the treatment group and the control group within each study were tested for significance using a Wilcoxon rank sum test. As shown in Figures 2B,C, anti-PD-1 treatment caused a statistically significant decrease in the observed rates of tumor volume change only at the highest dose (study 3, 10 mg/kg Q3Dx3). Anti-PD-L1 treatment caused statistically significant decreases in the observed rates of tumor volume change in studies 3 and 4, but not in study 5 (p-values are unadjusted, significance threshold is 0.01/6 = 0.0016 with Bonferroni correction) despite the same dosing regimen being used in the studies and similar tumor sizes at treatment start in studies 3 and 5 (Supplementary Figure S2, Supplementary Table S2).

Tumors with observed rates of change below 25 mm³ per day were classified as having a slow rate of change, while the others were classified as having a fast rate of change. The threshold of 25 mm³ per day was chosen as it identifies all mice that experience tumor regression after treatment, or a period of regression followed by progression (Supplementary Figure S3). While mice with fast tumor growth are present in all the dose groups, the number of tumors with a slow observed rate of tumor volume change increased as the anti-PD-1 dose increased (Figure 2D), indicating dose response.

In summary, substantial inter-study and inter-individual variability is present in the tumor volume profiles. Anti-PD-1 led to dose response in tumor growth inhibition, whereas dose ranging data were not available for anti-PD-L1 (Figure 2E).

Next, we parameterized the QSP model to understand the variabilities observed in the tumor volume profiles post anti-PD-(L)1 treatments. PK parameters for RMP1-14 and 10F.9G2 were estimated based on literature data (Zalba et al., 2020) (Supplementary Figure S1). The tumor volume-time profiles were used to fit the CD8⁺ T cell proliferation half-life, maximum rate of ADCC, tumor cell apoptosis rate, maximum rate of CD8⁺ T cell influx, CTL-mediated tumor cell killing rate, tumor cell proliferation rate, and tumor carrying capacity (Supplementary Table S1). Model fitting was approached iteratively, alternating between local sensitivity analysis and NLME modeling to identify a minimum number of parameters to describe the inter-study and inter-individual variabilities in the tumor volume profiles. The final model included random effects for CD8⁺ T cell proliferation half-life, tumor cell proliferation rate and tumor carrying capacity, and included study site ID as a covariate for population tumor cell proliferation rate and tumor carrying capacity (Table 2). The model was sufficient to replicate the variation in the tumor volume profiles (Figure 3A, Supplementary Figure S4, and Table 2). Removing the random effects for any of the three parameters diminished the goodness of fit (results not shown).

The individual predictions recapitulate the observations, except for very small tumor sizes (Figure 3B). The distribution of CD8⁺ T cell proliferation half-life is clustered around the average value, indicating difficulty estimating the distribution given the available data (Figure 3C). The distributions of individual tumor cell proliferation rates and tumor carrying capacities show a bimodal distribution (Figures 3D,E). The tumors in studies 1-3 tend to have lower tumor cell proliferation rates and higher carrying capacities, while the tumors in studies 4 and 5 have higher tumor cell proliferation rates and lower carrying capacities (Figure 3F).

In summary, the model replicates the variability in CT26 tumor volume-time profiles by incorporating inter-individual variations in three parameters, namely, CD8⁺ T cell proliferation half-life, tumor cell proliferation rate, and tumor carrying capacity.

We utilized the individual parameter estimates to identify the sources of variability in anti-tumor responses. Anti-tumor responses to anti-PD-(L)1 treatments were still categorized based on the observed rates of tumor volume change between day 7 and day 13 after treatment initiation. Mice with an observed rate of tumor volume change below 25 mm³/day were considered responders, while all others were considered non-responders.

Estimated values for individual parameters including CD8⁺ T cell proliferation half-life, tumor cell proliferation rate, and tumor carrying capacity were compared between responders and non-responders (Figure 4A). The CD8⁺ T cell proliferation half-lives and tumor carrying capacities do not significantly differ between responders and non-responders across studies (p-values are unadjusted, significance threshold is 0.05/3 = 0.016 with Bonferroni correction), potentially due to the high shrinkage values (Table 2). The average tumor cell proliferation rates significantly differ between responders and non-responders in study 5, and the tumor cell proliferation rate correlates with the observed rate of tumor volume change between day 7 and day 13 after treatment start for most studies and treatments (Supplementary Figure S5). The study-specific findings remain exploratory due to small sample sizes.

Furthermore, we sought to identify additional parameters that affect treatment response through local sensitivity analysis centered on parameter sets fitted to two mice with different responses to anti-PD-L1 treatment. The two mice were selected from the same study to have similar tumor cell proliferation rates and tumor carrying capacities and to respond to anti-PD-L1 10 mg/kg Q3Dx3 treatment differently. In simulations, the responder shows complete tumor regression in response to both 10 mg/kg Q3Dx3 anti-PD-1 and 10 mg/kg Q3Dx3 anti-PD-L1, while the non-responder shows continued tumor growth when treated with either anti-PD-1 or anti-PD-L1 (data not shown). For the local sensitivity analysis, each parameter was perturbed by 30% in both directions from the nominal parameter value one at a time. The new parameter set was used to simulate tumor volume-time profiles under the same conditions. The simulated AUCs were compared to the simulated AUCs with the nominal parameter values (Figure 4B).

In the control simulations, increasing the PD-1:PD-L1 dissociation constant (KD) or the fraction of tumor volume not occupied by cells (tumor void fraction) caused the tumor to completely regress. Decreasing the tumor cell proliferation rate, PD-L1 per tumor cell, or PD-1 per CD8⁺ T cell also caused the tumor to regress. The results suggest these intrinsic tumor properties determine the growth of the tumor in the absence of treatments (Figure 4B1).

The intrinsic tumor cell proliferation rate is also an important determinant of treatment response. A 30% decrease in the tumor cell proliferation rate caused the non-responder to experience complete tumor regression in both anti-PD-1 and anti-PD-L1 treatment simulations. Similarly, a 30% increase in the tumor cell proliferation rate caused the responder to experience a large increase in tumor volume (Figures 4B2,B3). This is consistent with the difference in average tumor cell proliferation rate between responders and non-responders identified in Figure 4A.

Increasing PD-L1 expression has a minimal impact on treatment response due to high receptor engagement during either treatment. At the nominal PD-L1 expression level, during anti-PD-1 simulations, receptor engagement can drop from approximately 85% to as low as 53% between doses. Consequently, increasing PD-L1 expression induced a 0.9% and 6.1% increase in tumor volume AUC for the non-responder and responder, respectively (Figure 4B2). In anti-PD-L1 simulations, increasing PD-L1 had minimal impact on AUC (Figure 4B3) due to PD-L1 receptor occupancy being maintained at >97% during treatment. For both treatments, 10 mg/kg Q3Dx3 is near the top of the dose response curve (Figures 1C,D), and a small increase in receptor expression is insufficient for anti-PD-(L)1 treatments to lose efficacy.

Properties of the immune system such as CD8 half-life, CD8 proliferation half-life, max CD8⁺ T cell influx rate, and the CTL-mediated tumor cell killing rate influenced response to anti-PD-(L)1 treatments (Figures 4B2,B3) but had no impact on tumor growth during control simulations (Figure 4B1). Decreasing the CD8⁺ T cell proliferation half-life causes the non-responder AUC to decrease 89% and 79% during anti-PD-1 and anti-PD-L1 treatments respectively while the opposite change causes the responder to experience a 24% and 405% increase in AUC, respectively. Similarly, in the non-responder, increasing the CD8⁺ T cell half-life causes a 54% and 75% decrease in AUC during anti-PD-1 and anti-PD-L1 treatments respectively while for the responder, decreasing the CD8⁺ T cell half-life causes a 26% and 460% increase in AUC. Smaller percent changes are seen in the responses to perturbations in the max CD8⁺ T cell influx rate and the CTL-mediated tumor cell killing rate.

In summary, while random effects on tumor carrying capacity, tumor cell proliferation rate, and CD8⁺ T cell proliferation half-life explain the variability in the profiles of CT26 tumors post anti-PD-(L)1 treatments, baseline tumor growth and treatment response are sensitive to both intrinsic tumor properties and immune system properties. Variability beyond the three parameters fit to individuals in this analysis should be considered to improve the prediction intervals in future analyses.

PD-1 and PD-L1 expression levels are potential sources of variability in clinical populations. In analysis of clinical trial results, patients are often stratified by PD-L1 expression (Doroshow et al., 2021). The preclinical model used here assumes all mice had the same PD-(L)1 expression levels. However, the preclinical QSP model can be used to examine how the modeled tumor volume for each mouse changes in response to different PD-1 and PD-L1 expression levels. For each fitted animal parameter set, the PD-1 or PD-L1 expression level was varied over a range of values and the model was used to predict tumor growth under the control, anti-PD-1, or anti-PD-L1 treatment. Figure 5 shows the percent changes in tumor volume from treatment start to 12 days after treatment start. Mice are ordered along the x-axis by estimated tumor cell proliferation rates. Day 12 was chosen as it is late enough to see differences between responders and non-responders, yet early enough to show differences in the strength of response. Experimentally quantified expression levels are marked by black arrows. Results stratified by study can be found in Supplementary Figures S6, S7.

The model simulations indicate that low PD-(L)1 expression was associated with spontaneous tumor regression in mice with low estimated tumor cell proliferation rates (Figures 5A,B). Tumor volume increases as PD-(L)1 expression levels or estimated tumor cell proliferation rates increase. Anti-PD-(L)1 treatments suppress tumor growth for higher expression levels compared to the control group. Mice with high estimated tumor cell proliferation rates show treatment resistance at all PD-(L)1 expression levels. Mice with a tumor cell proliferation rate below 0.33 per day show an inverted bell shape response curve as receptor expression increases, shown in the transition from light blue to dark blue and returning to light blue as receptor expression increases in each column. This indicates that each tumor cell proliferation rate has a different range of receptor expression levels where the treatment has the highest effect suggesting that dosing levels could be adjusted to each individual for the maximum effect.

In summary, PD-1 and PD-L1 expression levels determine treatment response, but not independent of the tumor cell proliferation rate. Mice with higher receptor expression levels or faster growing tumors may benefit from higher dose levels to saturate the receptors.