Quantitative systems pharmacology modeling provides insight into inter-mouse variability of Anti-CTLA4 response.

Clinical responses of immuno-oncology therapies are highly variable among patients. Similar response variability has been observed in syngeneic mouse models. Understanding of the variability in the mouse models may shed light on patient variability. Using a murine anti-CTLA4 antibody as a case study, we developed a quantitative systems pharmacology model to capture the molecular interactions of the antibody and relevant cellular interactions that lead to tumor cell killing. Nonlinear mixed effect modeling was incorporated to capture the inter-animal variability of tumor growth profiles in response to anti-CTLA4 treatment. The results suggested that intratumoral CD8+ T cell kinetics and tumor proliferation rate were the main drivers of the variability. In addition, simulations indicated that nonresponsive mice to anti-CTLA4 treatment could be converted to responders by increasing the number of intratumoral CD8+ T cells. The model provides a mechanistic starting point for translation of CTLA4 inhibitors from syngeneic mice to the clinic.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

There are several quantitative system pharmacology (QSP) models to understand the response of checkpoint inhibitors in humans. Comparable models and analyses for preclinical animal models are lacking.

WHAT QUESTION DID THIS STUDY ADDRESS?

This study integrated QSP and nonlinear mixed effect modeling to understand the variability in the tumor growth response upon treatment of anti‐CTLA4 antibody in CT26 tumor bearing syngeneic mice.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

This study seeks to develop a quantitative model to capture vast variability in the in vivo data of syngeneic mouse studies. The result pinpoints the importance of cellular composition and kinetics in deriving variability of tumor growth response upon anti‐CTLA4 antibody treatment.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

Establishment of a mouse QSP model for immuno‐oncology can improve preclinical study design. In addition, understanding the driver for response variability will benefit preclinical‐to‐clinical translational research.

INTRODUCTION

Cancer treatment has been revolutionized with the advent of immuno‐oncology (IO) therapies, such as ipilimumab, an anti‐CTLA4 antibody. This has led to a rapid increase in the number of IO targets and therapies entering clinical trials. , Despite the promise of IO therapy, the pioneer checkpoint inhibitor drugs only resulted in an overall response rate of ~15% in the United States in 2018. More systematic approaches are required to enable more patients to benefit and fewer clinical failures, including understanding of mechanisms driving variability to IO drug response.

IO drugs require to be studied in immune‐competent mice bearing syngeneic tumors so that an intact immune system is retained and the full interplay among the immune, stromal, and tumor cell systems can be studied. Specifically, upon treatment with an IO agent, a series of immune responses are initiated. , Because no two animals have the same immunogenic state, the immune responses can lead to significant variability upon interaction with tumor cells. In addition, animals’ immunogenicity states are influenced by the environment that ostensibly leads to larger inter‐experiment variation in responses. , The amount of variation in animals’ responses diminishes the utility of typical exposure–response modeling for clinical translation. Furthermore, many IO drug targets have low homology in mice and humans, and therefore mouse surrogate antibodies are required to conduct preclinical studies. Translational strategies based on the mechanism of action (MoA) of IO therapeutic agents are desired for translational research.

Herein, we present an integrated quantitative systems pharmacology (QSP) and nonlinear mixed effect (NLME) model to describe individual tumor growth profiles in CT26 syngeneic mice following anti‐CTLA4 antibody treatment. The model was used to (1) investigate the plausible biological processes driving observed response variability following administration of an anti‐CTLA4 antibody to CT26 syngeneic mice; (2) explore possible explanations for the observed differences between responders and nonresponders; and (3) generate mechanistic hypotheses for increasing the response rate to anti‐CTLA4 treatment.

Anti‐CTLA4 antibodies modulate the immune system through multiple avenues. CTLA4 is predominantly expressed on the surface of regulatory T cells (Tregs), activated CD8+ T cells, and a small percentage of monocytes. Binding of anti‐CTLA4 antibodies to Tregs abrogates the immune suppressive function of Tregs, and depletes Tregs via antibody‐dependent cellular phagocytosis. , , On the other hand, binding of anti‐CTLA4 antibodies to activated CD8+ T cells promotes CD28 co‐stimulatory signaling and subsequent cell activation, proliferation, and differentiation into cytotoxic T lymphocytes (CTLs), which act as serial killers for tumor cells. The QSP model consists of antibody pharmacokinetics (PK) and interactions with key immune cell populations, as well as interactions between the immune cell populations and tumor cells. The QSP model served as the foundation for an NLME model to probe the mechanistic rationale driving inter‐animal variability in tumor growth response upon anti‐CTLA4 antibody treatment in syngeneic mice. Importantly, the model provides a mechanistic starting point for future study of the translation of CTLA4 inhibitors from syngeneic mouse models to the clinic.

MATERIALS AND METHODS

Mouse datasets

All mouse studies were approved by the Pfizer Institutional Animal Care and Use Committee according to established guidelines. For the model training dataset, the CT26 cell line was inoculated into immunocompetent Balb/c mice. Tumor sizes were measured using calipers. Tumor volume was calculated as (width*width*length)/2. Treatments were initiated on day 0 when the tumor volume reached ~116 mm3. Mice were treated with phosphate–buffered saline (PBS, control) or anti‐mouse CTLA4 antibody (BioXell clone 9H10) given at 0.625, 1.25, 2.5, 5, or 10 mg/kg intravenously every 3 days for a total of three doses (Q3dx3) (n = 10/group). Tumor sizes were measured every 2–4 days. Animals were euthanized when tumor volumes reached 2000 mm3. Data from an independent experiment using the same syngeneic tumor model and antibody were used for model validation. In this second experiment, the average tumor volume of each group at day 0 before start of treatment was ~92 mm3. Mice were treated with PBS or 9H10 antibody given at 10 mg/kg intravenously Q3dx3 (n = 10/group).

Anti‐CTLA4 antibody PK in the mouse was from an in‐house study (data not shown). The PK data was well‐described by a two‐compartment model. The same PK parameters were used for all the mice because no individual animal PK data was available.

Figure 1a illustrates the key species and reactions in the QSP model. Mass action reactions were used in the model. At the molecular level, both membrane CTLA4 (mCTLA4) and soluble CTLA4 (sCTLA4) are present in the central, peripheral, and tumor compartments. The sCTLA4 can distribute between the central and the other two compartments. It was assumed that the steady‐state concentrations of sCTLA4 were the same in all the compartments. Distribution, elimination, and binding of anti‐CTLA4 antibody to sCTLA4 and mCTLA4 were described in all the compartments. The endogenous ligand B7 is included in the tumor compartment only because the tumor was assumed to be the site of action, and its interactions with mCTLA4 and sCTLA4 were included here. It was also assumed that the antibody and sCTLA4 complex could distribute between the central and the other two compartments, and the complexes were eliminated from the central compartment at the same rate as the antibody, whereas the complexes of antibody and mCTLA4 were eliminated and cleared at the same rate as mCTLA4. Ordinary differential equations (ODEs) that describe the PK model and molecular interactions are shown in Appendix S1. Receptor occupancy (RO) of antibody binding to CTLA4+ Tregs, CTLA4+ CD8+ T cells and other CTLA4+ cells in the tumor compartment were calculated and linked to tumor cell killing via depletion of Tregs and activation of CTLA4+ CD8+ T cells.

FIGURE 1

Structure of anti‐CTLA4 quantitative system pharmacology (QSP) model developed to capture the response variability in CT26 tumor volume profiles. (a) Schematics of the QSP model. ActCD8, activated CD8+ T cell; InactCD8, inactive CD8+ T cell; ProlifCD8, proliferating CD8+ T cell; CTL, cytotoxic T lymphocyte; TC, tumor cell; TCd, damaged tumor cell; Treg, regulatory T cell. (b) Individual tumor volume profiles in CT26 syngeneic mice in response to anti‐CTLA4 antibody 9H10 treatment given introvenously Q3dx3. (c) Number of mice exhibiting progressive disease (PD), partial response (PR), and complete response (CR) per dose group

At the cellular level, inactive CD8 (E ) is produced with a zero‐order production rate and is converted to activated CD8 (E ) as a reversible first order reaction. Initial concentration of E was derived from literature, and it was assumed to be the same for all the mice during model fitting. The inhibitory effect of Tregs on CD8+ T cell activation is represented by increasing the conversion rate of activated to inactivated CD8+ T cells. Binding of the antibody to CTLA4 on inactivated CD8+ T cells increases the activation rate of CD8+ T cells, and binding of the antibody to CTLA4 on Tregs increases the death rate of Tregs. Activated CD8+ T cells go through seven divisions (E , where i = 1…7) prior to differentiation into CTLs. Tumor cells undergo endogenous proliferation and second order death. CTLs can convert tumor cells into damaged tumor cells (TC ), which are removed with first order elimination. It was assumed that each CTL can kill up to 10 tumor cells before its exhaustion. There are 10 state variables representing CTLs at different exhaustion stages (C , where i = 0…10). Parameter annotations and ODEs describing the cellular kinetics are listed in Appendix S2 and S3, respectively. Parameter values are listed in Table 1.

TABLE 1

Anti‐CTLA4 model parameters

Description	Parameter	Value	Unit	Reference, Note
Target information
Steady state CTLA4 concentration in the plasma	C CTLA4,plasma	9.84E−05	nM	Fitted to PK data given the initial value of 4.23E−05 nM based on bottom‐up estimate
Steady‐state CTLA4 concentration in the peripheral compartment	C CTLA4,peripheral	9.84E−05	nM	Assumed to be the same as in the plasma
Number of CTLA4 receptors per CD8+ T cell	N CTLA4,CD8	1000	#/cell	6
Number of CTLA4 receptors per other CTLA4+ immune cell	N CTLA4,other	1000	#/cell	6
Number of CTLA4 receptors per Treg	N CTLA4,Treg	10,000	#/cell	6
Half‐life of CTLA4	THCTLA4	0.0833	Day	36
Soluble CTLA4 concentration in the plasma (central) compartment	C sCTLA4,plasma	0.4348	nM	Literature value of 10 ng/ml 20 , 22 , 23
Soluble CTLA4 concentration in the peripheral compartment	C sCTLA4,peripheral	0.4348	nM	Assumed to be the same as in the plasma
Soluble CTLA4 concentration in the tumor compartment	C sCTLA4,tumor	0.4348	nM	Assumed to be the same as in the plasma
Distribution half‐life of soluble CTLA4 from the central compartment to the peripheral compartment	TH12sCTLA4	0.0208	Day	Typical half‐life of soluble receptors
Partition coefficient of soluble CTLA4 from the central compartment to the peripheral compartment	P12sCTLA4	1	Unitless	Assumed
Distribution half‐life of soluble CTLA4 from the central compartment to the tumor compartment	TH13sCTLA4	0.0208	Day	Assumed
Partition coefficient of soluble CTLA4 from the central compartment to the tumor compartment	P13sCTLA4	1	Unitless	Assumed
Half‐life of soluble CTLA4	THsCTLA4	0.0417	Day	37
Average number of B7 receptors per tumor cell	N B7,TC	100	#/cell	38, 39
Half‐life of B7	THB7	0.0833	Day	Typical membrane protein half‐life
Association rate between B7 and CTLA4	konB7:CLTA4	86.4	nM−1 day−1	Typical binding rate for protein–protein interaction 40
Dissociation rate between B7 and CTLA4	KDB7:CLTA4	310	nM	Literature value 41 , 42
Cellular information
Concentration of CTLA4+ CD8 cells in the tumor	C CD8	2.4E−06	nM	Fitted to PK given the initial value of 2.4E−07 nM estimated from literature 4 , 43
Concentration of other CTLA4+ cells in the tumor	C other	1.83E−07	nM	Estimated from literature data 4 , 43
Concentration of CTLA4+ Tregs in the tumor	C Treg	1.32E−07	nM	Estimated from literature data 4 , 43
Concentration of tumor cells in the tumor	C TC	8.30E−04	nM	Estimated from literature data  19
Half‐life of inactive CD8+ T cells	THCD8	365	Day	Assume CD8 influx is very small (i.e., anti‐CTLA4 does not recruit new CD8 to the tumor)
Half‐life of Other cells	THother	3.5	day	Assumed
Basal activation of CD8+ T cells	k0 I,A	0	day	Assumed
Antibody information
Association rate between antibody and CTLA4	konAb:CTLA4	86.4	nM−1 day−1	Typical binding rate for protein–protein interaction 40
Dissociation rate between antibody and CTLA4	KDAb:CTLA4	6.7	nM	Estimated EC50 of 9H10 antibody from literature 44
Half‐life of antibody‐CTLA4 complexes	THAb:CTLA4	0.0833	Day	Assumed to be the same as CTLA4
Killing rate of Treg driven by antibody‐CTLA4 per Treg cell	kkillTreg	3.03E−4	Day−1	Estimated from initial fitting; insensitive parameter
PK
Volume of the central compartment	V1	0.001	L	Fitted from PK data
Volume of the peripheral compartment	V2	0.001	L	Fitted from PK data
Half‐life of antibody	THAb	3.2549	Day	Fitted from PK data
Distribution half‐life of antibody from the central compartment to the peripheral compartment	TH12Ab	0.0433	Day	Fitted from PK data
Partition coefficient of antibody from the central compartment to the peripheral compartment	P12Ab	0.3383	Unitless	Fitted from PK data
Distribution half‐life of antibody from the central compartment to the tumor compartment	TH13Ab	0.2897	Day	45
Partition coefficient of antibody from the central compartment to the peripheral compartment	P13Ab	1	Unitless	45

Abbreviations: EC50, half‐maximal effective concentration; PK, pharmacokinetics; Treg, regulatory T cell.

It was assumed that the concentrations of Tregs and CD8+ T cells in the tumor microenvironment are constant in the absence of antibody treatment as estimated from literature , so the number of T cells increases as tumor volume increases. The numbers of tumor cells before antibody treatment were calculated from the measured tumor volumes of individual mice. In the absence of treatment, tumor cells follow a logistic growth curve so that simulated tumor volumes reach a plateau at the prespecified carrying capacity. The changes in the numbers of immune cells and tumor cells together are reflected in the form of tumor volume. Further, immune cells produce sCTLA4 to a detectable circulating level in humans. , , , , The same circulating level of sCTLA4 was assumed for the mouse model.

QSP model development and local sensitivity analysis

The anti‐CTLA4 model was implemented using KroneckerBio version 0.5.2.3. Initial parameter estimation and simulations were performed using MATLAB R2018b (The MathWorks).

For model simulations, the model was run without treatment until all the states reached steady‐state (except tumor cells), and then the time was reset to zero and different dose levels of antibody were administered. Nominal parameter values of each animal were obtained by a combination of manual tuning and fitting the model to individual tumor volume profiles using the maximum likelihood method (fmincon) in MATLAB. Local sensitivity analysis was performed in the presence of a treatment of 10 mg/kg antibody given intravenously Q3dx3. Each parameter value was either increased or decreased two‐fold at a time while keeping all the other parameters at the nominal values. Tumor volumes at day 80 were used as the model output. Sensitivity of parameter i (S ) was quantified by the average percentage differences comparing to the model output of the nominal parameter values (Equation 1).

where , and represent the simulated tumor volumes at day 80 using nominal, doubling and halving the nominal parameter i, respectively.

NLME model development and validation

The model was transferred from KroneckerBio to Monolix for NLME modeling. Model fitting was performed in Monolix version 2019R1 (Lixoft, a Simulation Plus company). Model simulation was performed using mlxR 4.1.0 (Lixoft, a Simulation Plus company) in R 3.6.1.

To train the NLME model, nine parameters without literature values were estimated to capture the tumor volume profiles of the control and treatment groups simultaneously. The parameters were the half‐life of Tregs (THTreg), half‐life of CTLs (THCTL), half‐life of CD8+ T cell proliferation (THproliferation,CD8), proliferation rate of tumor cells (k proliferation,TC), secondary tumor cell death rate (k death,TC), first order elimination rate of damaged tumor cells (k apoptosis,TCd), killing rate of tumor cells driven by CTLs (k kill,CD8), activation rate of CD8+ T cells driven by antibody (k ) and deactivation rate of CD8+ T cells driven by Tregs (k Treg,CD8). Interindividual variability (IIV) was implemented for five of the nine parameters, specifically THproliferation,CD8, THCTL, k proliferation,TC, k death,TC, and k apoptosis,TCd that were selected based on the results of the local sensitivity analysis. Tumor volumes may not be measurable in mouse studies because mice need to be euthanized when tumor volumes reach 2000 mm or tumor sizes are too small to be measured by a caliper. Datapoints of unmeasurable tumors were censored during model fitting according to Monolix instruction. Goodness of model fitting was evaluated by diagnostic plots.

To validate the NLME model, k proliferation,TC was fitted to the tumor volume profiles of the control group from the model validation dataset while fixing the other parameters as fitted from the model training dataset. Median tumor volume profiles were then simulated using the updated k proliferation,TC.

RESULTS

A QSP model captured diverse tumor volume response to an anti‐CTLA4 antibody

A QSP model (Figure 1a) was developed to integrate the molecular and cellular MoA of a murine anti‐CTLA4 antibody, 9H10, in driving tumor regression in mice. In the model training dataset, the tumor volume profiles of the control group showed much less variability in contrast to the treatment groups, ranging in some cases from no response to complete regression for a given dose (Figure 1b). Individual tumor response was manually classified into three categories, including progressive disease (PD) where tumor volumes increased post‐treatment start, partial response (PR) where tumor volumes regressed initially post‐treatment start but then rebound, and complete response (CR) defined as at least 60% reduction in tumor volumes from treatment start to the time that a mouse was euthanized. Although tumor volume dynamics showed significant variability, the number of mice with a CR increased as the dose level increased (Figure 1c). Interestingly, the rates of tumor progression and regression were insensitive to the dose level (Figure S1).

The workflow for model parameterization for the training dataset is illustrated in Figure S2. The individual animal tumor volume profiles were independently fitted in MATLAB first, using the KronekerBio QSP model. Next, a local sensitivity analysis was performed in the presence of 10 mg/kg antibody treatment to identify sensitive parameters regulating tumor volume response. Finally, NLME modeling was performed in Monolix by applying IIV to the sensitive parameters identified from the local sensitivity analysis. Representative results from one mouse in each response category are shown in Figure S3. The local sensitivity analysis indicated that varying PK parameters of the antibody, including elimination rate (kelAb), intercompartment clearance (Q Ab) and intercompartment partition coefficient (P Ab) had minor impact on tumor volume response especially for the mice who did not respond well to treatment. Further, the model was used to predict antibody‐mediated RO on CD8+ T cells and Tregs at 10 mg/kg dose. Receptor occupancy was projected to be high (Figure S4), suggesting that drug exposure at the site of action and target binding may not limit the pharmacological response. Instead, local sensitivity analysis indicated that variation of CD8+ T cell and tumor cell‐related parameters, including THproliferation,CD8, k proliferation,TC, k death,TC, k apoptosis,TCd, and k kill,CD8 could potentially distinguish the different tumor volume profiles in response to treatment (Figure S3). Therefore, IIV of these five parameters were incorporated in an NLME model to capture inter‐animal variability. A series of diagnostic plots for the NLME model are shown, including prediction interval of tumor volume versus time profiles (Figure 2a), the observed versus predicted tumor volumes (Figures 2b and S5), and the distribution of random effect (Figure 2c). The model parameters are shown in Table 1. The fitted population and individual parameters are shown in Table 2 and Appendix S4, respectively.

FIGURE 2

An NLME model captures inter‐animal variability in tumor volume profiles. (a) Model fit of tumor volume profiles of 9H10 treated CT26 tumors. Black dots represent observed data. Red lines represent predicted median tumor volume profiles. Blue shaded areas represent predicted 90% confidence intervals. (b) Goodness‐of‐fit plot showing observed versus model predicted individual data points. (c) Distribution of the standardized random effects of model estimates. NLME, nonlinear mixed effect; PBS, phosphate–buffered saline

TABLE 2

Fitted parameter values for the anti‐CTLA4 NLME QSP model

Description	Parameter	Unit	Value	RSE (%)	Shrinkage
Fixed effect
Half‐life of CD8+ T cell proliferation	THproliferation,CD8	day	0.351	6.10	–
Tumor cell proliferation rate	k proliferation,TC	1/day	0.255	2.39	–
Secondary tumor cell death rate in a logistic growth model	k death,TC	1/nmol/day	1.95e+5	4.40	–
First order elimination rate of damaged tumor cells	k apoptosis,TCd	1/day	0.013	8.85	–
Killing rate of tumor cells driven by CTLs	k kill,CD8	1/nM/day	4.66e+3	11.8	–
Half‐life of Tregs	THTreg	day	7.37	–	–
Deactivation rate of CD8+ T cells driven by Tregs	k Treg,CD8	1/nM/day	3.46e+6	–	–
Half‐life of CTLs	THCTL	day	7.30	–	–
Activation rate of CD8+ T cells driven by antibody‐CTLA4 per CD8+ T cell	k I,A	1/day	8.64e−4	–	–
Random effect
Standard deviation of the random effect on THproliferation,CD8	ω_THproliferation,CD8	day	0.390	10.7	8.65%
Standard deviation of the random effect on k proliferation,TC	ω_kproliferation,TC	1/day	0.173	9.63	−3.98%
Standard deviation of the random effect on k death,TC	ω_kdeath,TC	1/nmol/day	0.219	19.7	12.6%
Standard deviation of the random effect on k apoptosis,TCd	ω_kapoptosis,TCd	1/day	0.520	11.6	−26.2%
Standard deviation of the random effect on k kill,CD8	ω_kkill,CD8	1/nM/day	0.688	11	0.473%
Error model
Proportional error	b	–	0.194	2.72	–

Abbreviations: CTL, cytotoxic T lymphocyte; NLME, nonlinear mixed effect; QSP, quantitative system pharmacology; RSE, relative standard error; Treg, regulatory T cell.

In summary, the CT26 syngeneic mouse data illustrated significant inter‐animal variability in the tumor volume profiles post treatment with 9H10. This variability was not explained by antibody exposure or CTLA4 RO. Instead, allowing five parameters, describing tumor cell and CD8+ T cell kinetics, to vary among mice captured inter‐animal variability in the tumor volume profiles using an NLME model.

Validation of the integrated QSP and NLME model

After the population parameters and IIVs were determined from the model training dataset, an independent dataset was used to validate the model. As shown in Figure 3, the fitted parameter values from the training dataset overestimated the tumor volume profiles of the validation dataset for both the control and the treatment groups. One possibility is that tumors in the validation dataset grew slower than tumors in the training dataset. We re‐estimated the population parameter k proliferation,TC based on the control group from the validation dataset while fixing the other population and IIV parameters as fitted from the training dataset. By updating k proliferation,TC from 0.255 to 0.194/day, model simulations captured the observations of the validation dataset within 95% confidence interval.

FIGURE 3

Validation of the NLME model using dataset from an independent experiment with CT26 syngeneic tumors. Validation dataset consists of a PBS (control) group (left) and a 9H10 treatment group at 10 mg/kg (right). Blue: Simulations with the model fitted from the training dataset. Red: Simulations with the updated tumor proliferation rate fitted from the validation dataset while keeping the other parameter values the same as the model fitted from the training dataset. CI, confidence interval; NLME, nonlinear mixed effect; PBS, phosphate–buffered saline

In summary, the intrinsic tumor growth rate of syngeneic tumor models can vary from experiment to experiment. The target‐related molecular and cellular kinetic model is predictive when the intrinsic tumor growth rate is adjusted, indicating validity of the fitted anti‐CTLA4 model. Further, these results suggested that variability of tumor growth and tumor killing parameters need not be linked to explain inter‐animal or inter‐experiment variability of anti‐CTLA4 response.

No single kinetic parameter can predict variability in anti‐CTLA4 treatment responses

Upon validation of the anti‐CTLA4 model, we tested to see if any of the parameters that were allowed to vary among individuals (THproliferation,CD8, k proliferation,TC, k death,TC, k apoptosis,TCd, and k kill,CD8) correlated with the type of tumor volume response. It was found that while a higher rate of tumor cell proliferation (k proliferation,TC) can lead to a worse response, and separate nonresponders (PD) and responders (CR and PR), no individual parameter can significantly separate the CR and PR groups (Figure 4a). Furthermore, no combination of two (Figure 4b) or more parameters (Figure S6) could distinguish the treatment responses. These results suggested that variability in the intrinsic properties of both immune cells and tumor cells contribute to the treatment response.

FIGURE 4

Single or combination of two parameters cannot distinguish tumors with progressive disease (PD), partial response (PR), and complete response (CR). (a) Box plots summarizing single parameters by response category. *Denotes level of significance; *p value < 0.05; **p value < 0.005; ***p value < 0.0001 from Wilcoxon test. (b) Paired plots for combination of any two parameters and distribution histogram for any single parameter. Colors indicate response groups

Simulations suggest the importance of CD8+ T cell activation in CT26 tumor growth response to anti‐CTLA4 treatment

Multiple parameters could potentially contribute to treatment response among animals. Because anti‐CTLA4 is known to enhance both T cell activation rate and Treg death, we investigated which mechanism was more important to tumor growth response. The local sensitivity analysis (Figure S3) suggested that tumor growth response is more sensitive to anti‐CTLA4 treatment‐mediated CD8+ T cell activation rate (k ) than Treg depletion rate (k kill,Treg). In addition, we tested the impact of baseline levels of CD8+ T cells and Tregs on tumor growth response. There are 12 nonresponsive mice in the training dataset. The model assumed that the initial ratio of CD8+ T cells to Tregs in tumors is ~18:1. , Using the parameter combinations of these 12 mice, simulations predicted no tumor regression regardless of the dose level (0.01–30 mg/kg; Figure S7A), suggesting that these tumors are intrinsically nonresponsive to anti‐CTLA4 treatment. However, increasing the ratio of CD8+ T cells to Tregs by doubling the number of CD8+ T cells, predicted dose dependent tumor regression for most of the mice (Figure S7B), suggesting that therapies boosting CD8+ T cell numbers can potentially convert a tumor from nonresponsive to responsive to 9H10 treatment, with the extent and dose at which tumor regression occurs depending on the combination of parameters. In contrast, increasing the ratio of CD8+ T cells to Tregs by halving the number of Tregs did not result in a similar change in response phenotype (Figure S7C). The simulations indicated that enhancement of CD8+ T cell activation is more important than enhancement of Treg death in determining tumor growth response to anti‐CTLA4 treatment. Results of a representative mouse are shown in Figure 5.

FIGURE 5

Baseline level of CD8+ T cells in the tumor microenvironment can alter treatment response. (a) Simulated dose response of a representative non‐responding mouse from 10 mg/kg treatment group of the training dataset. (b) Simulated dose response using the same parameter values as (a) except doubling the CD8+ T cell concentration at the time of treatment start

DISCUSSION

Although immunotherapy is a promising treatment for many forms of cancers, patients’ responses to immunotherapy are highly variable. Variability in individual animal responses to immunotherapy treatments was also observed in preclinical studies in syngeneic mouse tumor models. There are a few QSP models available describing checkpoint inhibitor treatments in patients. , , , However, comparable models and analyses for preclinical animal models are lacking. Understanding drivers for the variable responses in preclinical animal models will benefit preclinical‐to‐clinical translational research. Herein, we developed a QSP model to describe the tumor growth response to a murine anti‐CTLA4 antibody, 9H10, in mice bearing CT26 syngeneic tumors. The QSP model was then integrated with NLME modeling to describe the variability in individual tumor volume profiles.

During the development of the QSP model, literature was consulted for the MoA of an anti‐CTLA4 antibody, the biology of its molecular and cellular targets, and the interplay among the targeting cell populations in the mouse system. Given the goal of this work was to use the model to explain the response variabilities in the experimental data which were limited to tumor volume profiles, the principle of parsimony was applied by constraining the site of action to the tumor microenvironment, where the interactions between anti‐CTLA4 antibody, CD8+ T cells, Tregs, and tumor cells were the focus. Five parameters were chosen to explain individual mouse data, specifically half‐life of CD8+ T cell proliferation upon activation, tumor cell proliferation rate, tumor cell death rate, elimination rate of damaged tumor cells, and CTL‐mediated killing rate of tumor cells. These parameters determine the numbers of CTLs and tumor cells over time. The association of anti‐CTLA4 treatment response with CD8+ T cell numbers and tumor cell proliferation rate is supported by both preclinical and clinical data. , , , , Although most of the fitted parameters cannot currently be validated due to lack of experimental data, the population estimate of CD8+ T cell proliferation half‐life (8 h) is close to that reported by Yoon et al.

In general, QSP models have more parameters than PK/pharmacodynamic models to capture relevant biological mechanistic details, and as such issues of parameter identifiability are frequently experienced when available data are limited. , , There is a similar issue with the presented model. For each tumor volume profile, more than one combination of the model parameters can likely fit the observed data. The best method to resolve this issue while capturing individual data is an area of active research. One potential computational method is Metropolis‐Hastings sampling that can identify the probability distribution of parameters. Alternatively, an experimental solution is to supplement the available dataset by collecting time course data of inactive, proliferating and cytotoxic CD8+ T cells, Tregs, and tumor cells in tumor samples. While the relevance of collecting those data is obvious, it is challenging to design an optimal study that provides sufficient dynamic information of all cell populations of interest in practice. Having a QSP model available for early exploratory simulations prior to an experimental study can help select optimal dosing regimen and timepoints for data collection, and therefore maximize the chance of obtaining an information‐rich dataset.

Although caution needs to be taken when drawing conclusions from the predictions on cellular kinetics in the presented model, the model was able to predict tumor volume profiles upon 9H10 treatment and to understand the impact of parameter combinations. Comparison of our model training and validation datasets indicated a difference in the intrinsic tumor growth rates between experiments. This can be caused by factors such as cell line passage number, housing environment, human handling, or animal age. In addition, tumor cell growth rate was identified as a key parameter that determines responsiveness to 9H10 treatment. This is not surprising because tumors grow exponentially, and small variations in the tumor growth rate can lead to significant differences in tumor sizes over time post‐treatment. The impact of tumor growth rate on treatment response is not unique to immune‐competent mice. Although outside the scope of the present study, more data from independent experiments can be used to estimate the inter‐experiment variability.

Our results suggested that rather than variability in drug exposure, variabilities in cell composition in the tumor microenvironment at the time of treatment start, and cellular kinetics are the main driver for inter‐animal variability in response to anti‐CTLA4 treatment. Activation of CD8+ T cells may have a pivotal role in determining the responsiveness of a tumor to anti‐CTLA4 treatment. As such, there may be a potential therapeutic advantage of combining an anti‐CTLA4 treatment with therapeutics that expand CD8+ T cells to treat nonresponsive and partially responsive tumors. These results are supported by the observed variability in intratumor cell compositions found in screening studies of syngeneic tumor models. ,

In summary, the presented anti‐CTLA4 QSP mouse model accounts for known pharmacological and biological mechanisms relevant to anti‐CTLA4 treatment and can capture individual animal responses in tumor volume profiles post treatment by adding IIV to parameters related to key tumor and immune cell kinetics. It provides biologically plausible explanations for the observed differences between responders and nonresponders to anti‐CTLA4 treatment. Although the model prediction on intratumoral cellular kinetics remains to be validated, the study demonstrates the power of QSP modeling in answering complex mechanistic questions and generating testable hypotheses. As future studies investigate tumor cell kinetic responses to immunotherapy, this model can help to guide study designs and identify key gaps in our understanding by comparing model predictions to emerging data. The platform nature of the model makes it amenable to reuse and repurposing to support diverse decisions from early drug discovery to clinical studies. More importantly, the model could be translated to human by incorporating predicted human PK, systems parameters such as tumor growth rate, and knowledge of tumors and immune cell profiles. The human model could be used for virtual patient simulations to predict optimal dosing regimens and likelihood of response given a patient’s individual characteristics. Early simulations based on virtual patient modeling can influence clinical data collections and underscore the need for quantification of tumor sizes and immune profiles in individual patient's blood and tumor samples, in alignment with precision medicine strategies.

CONFLICT OF INTEREST

W.Q., J.N., A.H., and A.B. are (or were at the time this work was conducted) employees of Pfizer Inc. L.L., C.Y., F.H., A.M., and L.G. are (or were at the time this work was conducted) employees of Applied BioMath, Inc. and conducted part of this work under contract with Pfizer Inc. L.L. is currently an employee of Biogen Inc. C.Y. is currently an employee of Takeda Inc. A.H. is currently an employee of Regeneron Pharmaceuticals Inc. L.G. is currently an employee of The Leukemia & Lymphoma Society. A.B. is currently an employee of Applied BioMath Inc.

AUTHOR CONTRIBUTIONS

W.Q., F.H., and A.B. wrote the manuscript. L.G. and A.B. designed the research. W.Q., L.L., C.Y., F.H., J.N., A.M., and A.H. performed the research. W.Q., L.L., C.Y., and F.H. analyzed the data.

Figure S1

Click here for additional data file.

Figure S2

Click here for additional data file.

Figure S3

Click here for additional data file.

Figure S4

Click here for additional data file.

Figure S5

Click here for additional data file.

Figure S6

Click here for additional data file.

Figure S7

Click here for additional data file.

Appendix S1

Click here for additional data file.

Appendix S2

Click here for additional data file.

Appendix S3

Click here for additional data file.

Appendix S4

Click here for additional data file.