Mathematical modeling identifies LAG3 and HAVCR2 as biomarkers of T cell exhaustion in melanoma

Summary Cytotoxic T lymphocytes (CTLs) control tumors via lysis of antigen-presenting targets or through secretion of cytokines such as interferon-γ (IFNG), which inhibit tumor cell proliferation. Improved understanding of CTL interactions within solid tumors will aid the development of immunotherapeutic strategies against cancer. In this study, we take a systems biology approach to compare the importance of cytolytic versus IFNG-mediated cytostatic effects in a murine melanoma model (B16F10) and to dissect the contribution of immune checkpoints HAVCR2, LAG3, and PDCD1/CD274 to CTL exhaustion. We integrated multimodal data to inform an ordinary differential equation (ODE) model of CTL activities inside the tumor. Our model predicted that CTL cytotoxicity played only a minor role in tumor control relative to the cytostatic effects of IFNG. Furthermore, our analysis revealed that within B16F10 melanomas HAVCR2 and LAG3 better characterize the development of a dysfunctional CTL phenotype than does the PDCD1/CD274 axis.


INTRODUCTION
Immunotherapy is an emerging strategy for treatment of cancer, with an ever growing number of immunotherapies having reached clinical trials or been approved already. 1 Blood cancers were among the first to be successfully treated with immunotherapy 2 ; to date solid tumors have proved to be more challenging. Despite this, several treatments are already available for solid tumors and many more are under trial. 3 Despite some success with immunotherapy so far, there remains a pressing need for greater mechanistic understanding of immune cell interactions within solid tumors. Such understanding may help expand the scope of immunotherapies to different cancers, identify biomarkers to predict patients benefiting from immunotherapy, 4 optimize dosing schedules for immunotherapies, 5,6 or identify combination therapeutic strategies. 7 Computational models are a useful tool to develop such understanding since they can link data from different sources and make quantitative predictions for what we should expect under different conditions. CD8 + cytotoxic T lymphocytes (CTLs) are key players in the anti-cancer immune response, and many immunotherapy strategies are focused on them. Prominent examples are blockade of inhibitory receptors such as programmed cell death protein 1 PDCD1 (common alias: PD-1, encoded by the Pdcd1 gene) expressed on CTLs to ''remove the brakes'' on the immune response 8,9 or adoptive transfer of engineered T cells. 10,11 Understanding the functioning of CTLs inside tumors is of foundational importance for rational immunotherapy design. Secretion of the cytokine interferon-g (IFNG) is a hallmark of activated CTLs, yet due to its pleiotropic effects the exact effects of this cytokine in solid tumors remain poorly understood. Some have even noted the ''paradoxical'' role of IFNG in tumor progression, 12 paradoxical in the sense that IFNG can have both pro-tumor and anti-tumor effects. Among the pro-tumor effects, IFNG can lead to recruitment of suppressive cells like regulatory T cells or myeloid-derived suppressor cells (MDSCs) or induce expression of immune checkpoint ligands like programmed death-ligand 1CD274 (common alias: PD-L1, encoded by the Cd274 gene) on tumor cells. 13,14 Among the anti-tumor effects, IFNG can aid in the recruitment of innate immune effectors, kill tumor cells, or exert antiproliferative effects on tumor cells. [15][16][17] Despite the recognition that IFNG has both pro-tumor and anti-tumor effects on cancer cells, these effects have not been quantified in great detail. With respect to anti-tumoral effects, we have previously used computational models to demonstrate how an antiproliferative effect mediated by cytokines could potently limit tumor progression since through cytokine signaling CTLs can control many tumor cells, stalling tumor growth and buying time for killing of tumor cells by CTLs. 18,19 However, in our previous modeling work, no direct data linking the proliferation of tumor cells to cytokine levels inside the tumors were available. Indeed, antiproliferative effects of IFNG are mediated by inhibitors of cyclin-dependent kinases, which result in arrest of tumor cells at the G 1 phase of the cell cycle, as shown in a variety of cell lines. [15][16][17] With respect to quantification of pro-tumoral effects of IFNG, it is currently unclear what is the relative importance of various T cell exhaustion markers that limit the immune response against cancer. Although the clinical success of immune checkpoint blocking antibodies demonstrates the significance of T cell exhaustion in tumor progression, there remains an urgent need for mechanistic and quantitative insight into exhaustion pathways in order to identify effective treatment combinations and biomarkers predicting response.
In the current study we sought to gain a quantitative understanding of the role of both the dynamics of IFNG-mediated antiproliferative effects and the exhaustion marker dynamics. Therefore, we utilize data from Matshushita and coworkers where they explicitly explored the antiproliferative effects of IFNG following in vivo T cell adoptive transfer by means of a cell cycle sensor in B16F10 melanoma. 15 Moreover, we exploited dynamic mRNA expression data from that study, allowing us to longitudinally quantify the expression of various exhaustion markers throughout adoptive transfer treatments. Although mathematical models of T cell exhaustion have been developed, 20,21 these have so far focused only on PDCD1-or cytotoxic T-lymphocyte-associated protein 4 (CTLA-4)-mediated inhibition and not on other checkpoints available in our data, thus necessitating development of a new mathematical model here. As a first step to quantify the importance of IFNG in tumor control relative to the canonical killing functions of CTLs, we developed an ODE model, which integrated data from Matshushita et al. 15 to arrive at a coherent, quantitative description of the intratumoral activities of CTLs and their interactions with the tumor following adoptive transfer. Consistent with our previous work on this topic, 19 our model predicted that the cytotoxic effects of CTLs were small and that cytostatic effects of IFNG were responsible for almost all of the observed difference in tumor growth between CTL-treated versus untreated tumors. As a second step, we extended our analysis and modeling to integrate the dynamics of various exhaustion markers as potential explanation for the short window of IFNG production, with CTLs losing the ability to produce IFNG within days of CTLs infiltrating the tumor. Markers of CTL exhaustion such as HAVCR2 (encoded by the Havcr2 gene) and LAG3 (encoded by the Lag3 gene) increased over this period, suggesting that CTLs had become exhausted inside the tumor. In contrast to HAVCR2 and LAG3, the dynamics of PDCD1 and CD274 did not coincide with the dynamics of CTL exhaustion, suggesting a relatively minor role for PD-1/PD-L1 as determinants of CTL exhaustion in the B16F10 melanoma model, at least at late stages of anti-tumour immune responses.

Presence of CTLs correlates with cell-cycle arrest in tumor cells
A previous study employed a Fluorescence Ubiquitin Cell Cycle Indicator (FUCCI) cell-cycle sensor to show that adoptive transfer of CTLs induced G1-phase cell-cycle arrest of B16F10 tumor cells in an IFNG-dependent manner; 15 however, the temporal evolution of this arrested state and correlation with the number of tumor-infiltrating CTLs was not explicitly quantified. Therefore, we exploited previously unquantified images from the same study, taken at multiple time points after CTL transfer, to estimate the number of tumor-infiltrating CTLs and B16F10 tumor nuclei. To quantify the number of B16F10 nuclei in either the G 1 phase or in the S-G 2 -M phases at different time points after CTL transfer, we developed automated pipelines using the ilastik 22 cell density estimation tool (see STAR Methods). Comparison of the ilastik predictions for small subregions of sample images ( Figure 1A) to manual counts made for the same images demonstrated that our pipeline was reliable ( Figure 1B). Moreover, our estimated densities of G 1 phase ( Figure 1C) or S-G 2 -M phase ( Figure 1D) on day 3 were comparable to those in the study of Matsushita et al., 15 as were the ratios of cells in G 1 :S-G 2 -M phases ( Figure 1E). Due to their irregular morphology, detection of CTLs was difficult to automate, so we instead performed a manual count of the number of CTLs across all images ( Figure 1F). We found a strong negative correlation between the number of CTLs and the G 1 :S-G 2 -M ratio in the sample images ( Figure 1G), with a Pearson's correlation coefficient of À0.60 (95% confidence interval between À0.78 and À0.33) allowing us to reject the null hypothesis of no correlation (p = 0.00015). In summary, the G 1 cell-cycle arrest following CTL transfer lasted for up to 5 days, and its temporal dynamics were closely linked to the presence of CTLs inside the tumor.

Tumor cell-cycle arrest correlates with tumor growth reduction
To check if the temporary G 1 cell-cycle arrest was consistent with tumor volume progression, we also incorporated tumor volume measurements into our analysis. iScience Article sufficient to describe tumor progression over the studied interval (Figure 2A), i.e., within the observed range of tumor sizes, there was not yet any indication for a potential carrying capacity limiting tumor growth. Note that although there was considerable variability between tumor growth rates of individual mice, the population estimate (0.4 day À1 ) was similar to the mean of the growth rates estimated from  iScience Article individual mice ( Figure S1). Volume estimates were available from three separate experiments with CTL treatment ( Figure 2B). We noted some minor yet apparent systematic differences between experiments. For instance, almost all volumes recorded on day À1 were larger in one of the biological replicates (compare red and green points in Figures 2A and 2B, day À1). Despite these minor discrepancies, the broad pattern of tumor progression was similar across replicates, with substantially arrested growth between days 3-7. Nevertheless, such systematic differences between experiments could potentially distort our results, for example, because the switching from 12 mice to 2 mice between measurements going from days 7-10 ( Figure 2B) would artificially introduce a period of tumor growth above even the untreated growth rate into our data. To avoid this issue, we converted the data into estimates of the tumor growth rate between measurement intervals for both the data without CTL transfer ( Figure 2C) and those with CTL transfer ( Figure 2D). For the experiments where CTLs were transferred, this resulted in consistent values between experiments and allowed us to safely incorporate the additional measurements from the 2 mice that were recorded up until day 14. From this analysis, reduced tumor growth was apparent between days 3-7 ( Figure 2D; points centered on days 4 and 6) but growth recovery in the measurement interval between days 7-10 ( Figure 2D; point centered on day 8.5). Therefore the period of tumor growth reduction was coincident with the period of G 1 -phase tumor cell-cycle arrest ( Figure 1E) and by extension also coincident with the presence of CTLs within the tumor ( Figure 1F) Loss of IFNG production precedes loss of CTLs from tumors IFNG secreted by CTLs was the putative agent which led to cell-cycle arrest and the transient reduction of tumor progression in our studied data. 15 As a proxy for IFNG levels inside the tumor, we used mRNA expression data recorded within the same experiments as the previously analyzed image ( Figure 1) and volume progression data ( Figure 2). We found that Cd8a transcription dynamics ( Figure 3A, row 1) matched the CTL dynamics measured in the images ( Figure 1F), indicating agreement between the transcriptomics data and the imaging data with respect to CTL abundance. However, the dynamics of Ifng transcription appeared much different from those of the CTLs ( Figure 3A, row 2). Ifng transcription peaked sharply on iScience Article day 3 after CTL transfer but had dropped sharply by day 5 and returned to basal levels on day 7, when CTLs still remained inside the tumor. To verify the dynamics of Ifng, we also checked Stat1 and Socs1 ( Figure 3A, rows 3-4), which are downstream of the IFNG receptor 23 in the IFNG signaling pathway. These followed very similar dynamics to Ifng mRNA, lending support to the idea that the Ifng mRNA expression data were a suitable proxy for IFNG signaling dynamics inside the tumor.
We hypothesized that the difference in dynamics between CTLs and Ifng transcription was due to a gradual CTL exhaustion inside the tumor, leading to a loss of their effector functions. Exhausted T cells display hierarchical loss of effector functions including proliferative ability, capacity to kill target cells, and secretion of cytokines such as IFNG. 24,25 Several genes are associated with the exhausted T cell state, 26,27 and as T cells become progressively more exhausted, they express a greater diversity of inhibitory receptors. 25 Indeed, we could identify transcripts for a number of well-described immune checkpoint molecules in the mRNA dataset, including PDCD1, its ligand CD274, LAG3, and HAVCR2 (Figure 3B). Overall, our analysis suggests that the pulse of IFNG production remains brief despite CTLs still being present within the tumor and might be due to development of an exhausted phenotype among the transferred CTLs.
IFNG transcription dynamics are compatible with G 1 -phase tumor cell-cycle arrest Due to the early reduction in IFNG signaling, it is unclear whether IFNG can be entirely responsible for the G 1 -phase tumor cell-cycle arrest, which followed highly similar dynamics to the CTLs. To test the compatibility of the Ifng transcription data with the dynamics of the CTLs and the tumor cell-cycle dynamics, we developed an ODE model. Our ODE model describing the interactions between CTLs and the tumor (Figure 4A) features an explicit description of the cell cycle of tumor cells, in which they cycle from G 1 phase into S-G 2 -M phases at rate k gs and then back into G 1 phase at rate k sg . The model also features CTLs, which kill tumor cells at rate k e and produce IFNG, which precludes tumor cells from transferring from G1 phase to S phase. The sensitivity of tumor cell-cycle arrest to IFNG is determined by the parameter k i . To test the iScience Article contribution of the two CTL effector functions to tumor control (i.e., killing and antiproliferative effect), we linearly interpolated between the experimental data for the number of CTLs ( Figure 4B) and for Ifng expression ( Figure 4C) and used these interpolations directly as inputs to our model. Subsequently, we tested different combinations of the parameters k e and k i ( Figure 4D) to find the best fit to the tumor growth rate ( Figure 4E, red line) and the S-G 2 -M:G 1 ratios ( Figure 4F, red line) determined from the experimental data. Our best-fitting parameter set ( Figure 4D; marked with black circle) had a value of k e = 0.9 (CTL À1 day À1 ) although other values for k e in the range 0-3 (CTL À1 day À1 ) led to relatively low errors, consistent with killing rates we have previously estimated for CTLs against B16F10 melanoma tumors. 19 The bestfitting value for the antiproliferative effect (k i = 8.1 IFN À1 mm 3 ) led to sharp reductions in the transition rate of tumor cells out of the G 1 phase for the Ifng expression levels found in our data ( Figure 4G). At the peak of Ifng expression on day 3, the transition rate from G 1 to S-G 2 -M phases (k gs ) was reduced to 7% of its original value, and even at the lower Ifng expression levels measured on other days, k gs was  When we took the best-fitting parameters and disabled killing by setting k e = 0 ( Figure 4D; marked with black square), most of the tumor growth reduction was preserved ( Figure 4E, blue dashed line). In contrast, taking our best-fitting parameters and disabling the antiproliferative effect of IFNG ( Figure 4D; marked with triangle) resulted in only a very small reduction in the net growth rate of the tumors ( Figure 4E, green dashed line). These results show that for the best-fitting parameter combination, the antiproliferative effect rather than CTL cytotoxicity played a dominant role in tumor control. We also tested whether only cytotoxicity, or only an antiproliferative effect, was compatible with the data by fitting a reduced model with k e or k i set equal to zero ( Figure S2). The Akaike information criterion (AIC) of the full model was 8.0; for the model without killing, the AIC was 6.7, and for the model without antiproliferative effect, the AIC was 131.3. Overall, these results support our previous analysis showing that an antiproliferative effect of IFNG is more important than CTL cytotoxicity to control B16F10 tumors. 19 Moreover, these results show that the dynamics of IFNG are compatible with the dynamics of the tumor cell-cycle arrest, despite the apparently short duration of IFNG production.

CTL exhaustion quantitatively explains Ifng transcription dynamics
In order to explain the dynamics of CTLs and IFNG and to quantify the importance of different immune checkpoints in these dynamics, we extended our ODE model ( Figure 5A). In this extended model, CTLs infiltrate the tumor at a basal rate s 0 , expand within the tumor at rate s e , and die with rate d e . In addition to their killing of tumor cells, CTLs inside the tumor produce IFNG. Finally our model includes the immune checkpoints LAG3, HAVCR2, PDCD1, and its ligand CD274, which decrease the activity of CTLs (see STAR Methods). We fit this ODE model simultaneously to all the experimental data discussed in Figures 1, 2, and 3 (see STAR Methods; best-fitting parameter sets provided in Table 1), resulting in a model, which nicely explained the observed dynamics ( Figures 5B and 5C). Among the best-fitting parameter sets, we observed a tendency for either HAVCR2 or LAG3 to make the most significant contribution toward CTL exhaustion ( Figure 5D). The dominant immune checkpoint tended to decay more slowly in the model fit ( Figure 5C), but regardless of which inhibitor was dominant, CTL dynamics and tumor growth ( Figure 5B) as well as total CTL activity ( Figure 5E) were similar. Importantly, without using any checkpoints, we could not obtain a good fit to any of the experimental measurements ( Figure S3 green lines), demonstrating that a regulatory mechanism like T cell exhaustion is required to explain the T cell anti-tumour response.
Since none of our best-fitting parameter sets predicted a significant role for PDCD1-CD274 inhibition in the in vivo experimental setting with B16F10 tumors, we asked how important the small early contribution of PDCD1-CD274 to CTL exhaustion was for the model dynamics. To address this, we disabled each inhibitor in turn. In the LAG3-dominant model ( Figure S4, red lines), only disabling LAG3 had a significant effect on model dynamics, but disabling PDCD1/CD274 or HAVCR2 had negligible effect. Similar results were obtained for the HAVCR2-dominant model ( Figure S4, blue lines) in which only disabling HAVCR2 disrupted model dynamics, but disabling PDCD1/CD274 or LAG3 had no effect. These results suggest that inhibition by either LAG3 or HAVCR2 is required to explain the experimental data, but PDCD1/CD274 inhibition is dispensable. However, it could be that our stochastic optimization procedure simply did not find any parameter sets where PDCD1/CD274 inhibition was important, so we also tried fitting our model with each immune checkpoint separately. Including only PDCD1-CD274 as an inhibitor (purple lines; Figure S3) resulted in a poor fit to the number of CTLs counted inside the tumor, as well as the dynamics of PDCD1 and CD274 themselves. In contrast, with LAG3 ( Figure S3, blue dashed lines) or HAVCR2 ( Figure S3, red dotted lines) as sole inhibitors, we achieved an equally good fit as with the full model including all checkpoints together ( Figure S3, black lines). Application of the AIC to our fits confirmed that having PDCD1/CD274 as the only inhibitors led to significantly worse performance than the other tested models ( Figure 5F), while the models with only HAVCR2 or only LAG3 had a slightly improved score than the full model because of the need to estimate fewer parameters. Note that these model variants gave similar predictions to the simplified model, where CTLs and IFNG were used as inputs to the model, in terms of the relative importance of cytotoxic effects versus antiproliferative effects ( Figure S5).
To further verify the relative importance of each checkpoint, we utilized a variant of the Metropolis Hastings algorithm to assess the practical identifiability of our estimated parameters. Four of our estimated parameters (d e , d i , k i , s e ) were fully identifiable ( Figure S6).  Figure S7A). Even extremely high or low values of d l , d p , or d pl led to relatively moderate decreases in log likelihood, whereas the log likelihood decreased more strongly for large values of d t . This occurs because, in this setting, only HAVCR2 contributes significantly to exhaustion, so modifying d l , d p , or d pl has hardly any effect on model dynamics except for those of the inhibitor in question. Thus, the uncertainty in the currently used immune checkpoint data is too great to accurately estimate the values of the immune checkpoint disappearance rates. Moreover, parameters controlling the contribution of each immune checkpoint toward CTL exhaustion (k p , k t , k l ) were not fully identifiable iScience Article because the minimum credible value for all of these parameters was unbounded on a logarithmic scale. Importantly, however, k p (contribution of PDCD1/CD274 toward CTL exhaustion) was predicted to be close to zero ( Figure S6), while either k t or k l (contributions of respectively HAVCR2 and LAG3 toward CTL exhaustion) had to have a high value ( Figure 5G), which generally co-occurred with a low value for the disappearance rate for the respective checkpoint ( Figures S7B and S7C). Simulations with parameter sets drawn from the estimated posterior distributions described the data quite well but revealed 2-3 main trajectories suggesting the possibility of multiple local optima (Figures S7D and S7E). Coloring iScience Article by parameter values revealed that two trajectories were enriched for either low or high values of the parameter d p ( Figure S8A). The other two trajectories that stood out were enriched for either low or high values of d l ( Figure S8B) and d t ( Figure S8C) reflecting the previously identified HAVCR2dominant and LAG-dominant fits. Otherwise, there was little correlation between fitted parameters ( Figure S9).
Overall, these results show that the brief window of IFNG production is quantitatively consistent with development of an exhausted state among CTLs and that expression of the immune checkpoint molecules LAG3 and HAVCR2 correlates best with the development of this exhausted state. Our analysis suggests that, of the three immune checkpoints considered, PDCD1/CD274 is the least important determinant of the exhausted CTL state; however, our model remains compatible with PDCD1/CD274 playing a role in CTL exhaustion at early time points after CTL infiltration of the tumor.

DISCUSSION
In a previous study on which our work was built, Matsushita et al. found that tumor control of B16F10 melanoma by CTLs was mediated by a combination of cytotoxic and cytostatic effects, with cytostatic effects being due to IFNG-mediated cell-cycle arrest of the melanoma cells. 15 However, the progression of cytostatic and antiproliferative effects over time was not explicitly explored. Here, we analyzed image data, tumor volume measurements, and transcriptomics data from the study by Matsushita et al., 15 using data acquired at multiple time points after CTL transfer. We found that the presence of CTLs inside the tumor strongly correlated with tumor cell-cycle arrest, as well as with the inhibition of volumetric tumor growth. However, IFNG signaling within tumors followed early dynamics, with CTLs primarily producing IFNG early after arrival in the tumors. Since IFNG loss preceded the recovery of tumor cell proliferation, it was unclear whether IFNG signaling could completely account for the observed tumor cell-cycle arrest and what role T cell exhaustion had in these processes. Therefore, we developed an ODE model to describe tumor growth, CTL infiltration, CTL production of IFNG and subsequent interference with cell-cycle progression, and also tumor cell killing by CTLs. Using this model we could describe all the experimental data, which led us to conclude that IFNG-mediated tumor cell-cycle arrest, together with killing of tumor cells by CTLs, was sufficient mechanism to account for the experimental data. We also used our models to compare the contribution of CTL-mediated cytotoxic or cytostatic effects toward tumor control. Our models predicted only a minor contribution of CTL killing toward tumor control compared to the IFNG-mediated cell-cycle arrest, consistent with our prior findings. 19 As part of our study, we developed a model describing the dynamics and effector functions of tumor-infiltrating CTLs. Based on mRNA expression data, Ifng transcription peaked on day 3, had fallen sharply by day 5, and was virtually zero on day 7. This was in contrast to the number of CTLs, which remained present in similar numbers on days 3 and 5 and were still observable in reasonable numbers at late time points. Therefore, our model included the development of CTL exhaustion in order to account for this loss in ability to produce IFNG. Another alternative, which might explain reduced CTL function, is a reduction in stimulation via the T cell receptor, which could also account for the decreasing PDCD1 expression on the CTLs since PDCD1 is normally expressed during CTL activation. 28 However, we consider this unlikely since the tumors remained large during the period of the experiments, meaning that CTLs should remain in contact with tumor cells and thus remain stimulated. CTL exhaustion is identified by a progressive increase in the number and diversity of inhibitory receptors expressed by CTLs. 25,29,30 We identified four well-known inhibitory molecules among the available transcriptomics data: LAG3, HAVCR2, PDCD1, and the PDCD1 ligand CD274. With our model we were able to obtain good fits if the exhausted state was correlated with HAVCR2 or LAG3 but not with PDCD1/CD274, which was due to the early peak of Pdcd1 and Cd274 transcription that was already well in decline on day 5 while CTL numbers in the tumor remained high. This early peak was not compatible with the idea that CTLs were becoming gradually more exhausted over time. On the other hand, Lag3 and Havcr2 increased relative to the CTLs over time and therefore correlated most with the loss of Ifng transcription. Consistent with our model prediction, LAG3 and HAVCR2 have been previously shown to have high correlation with a dysfunctional ''exhausted'' phenotype in CD8 + CTLs in melanoma. 26 Our model was not compatible with the dynamics of PDCD1/CD274 as sole correlates of the exhausted state, which appears at first sight to contradict reports indicating that PDCD1/CD274 signaling is relevant for immunosuppression in melanoma, 31 iScience Article points, it may be that CD274/PDCD1 plays only an initial role in immune suppression and that this role is taken over later by other checkpoints. This is consistent with findings that blockade of LAG3 as well as PDCD1 receptors is required to prevent relapse in melanoma 31 and with recent clinical results demonstrating that co-administration of PDCD1 and LAG3 inhibitors improves survival compared to administration of PDCD1 inhibitors alone. 35 One caveat for the data employed to fit our model is that only one probe was available per checkpoint. Therefore future experiments should confirm the dynamics of the expression of these immune checkpoint molecules and further investigate their contribution to T cell exhaustion. Another potential criticism of our approach is that it may be unfair to consider LAG3 and HAVCR2 separately from PDCD1/CD274-mediated immunosuppression since PDCD1 upregulation usually precedes LAG3 and HAVCR2 upregulation following T cell activation. 25 Indeed, our model is phenomenological with respect to checkpoint dynamics, so any causal relationships, which may exist between the expression of different checkpoints, are not taken into account because additional data would be required to describe such causal interactions. Therefore, the immune checkpoint component of our model should be viewed as a means to explore hypothetical scenarios where only single checkpoints have a role in T cell exhaustion.
Our model implies that the reduced activity of CTLs and in particular the apparent reduction in IFNG, which preceded the disappearance of CTLs in the tumor by several days, could be explained by the development of an exhausted phenotype in the tumor-infiltrating CTLs. Moreover, in our model IFNG played an important role in driving this exhausted phenotype. For exhaustion related to the CD274/PDCD1 axis, this is clearly justified because IFNG can induce upregulation of CD274 on tumor cells. 13 Moreover, IFNG induces increased antigen presentation on tumor cells, 12 which should lead to increased stimulation of CTLs via their T cell receptors. This could explain the contribution of IFNG toward upregulation of the other immune checkpoints included in our model, which are more commonly associated with excessive and prolonged exposure to antigen. 30 In order to further study the dynamics of the CTL population in the tumor, it would be useful to perform a second transfer of CTLs, which may help elucidate the extent to which the mechanisms of decline in CTL function are due to transferred CTLs becoming exhausted (and therefore a second transfer of ''fresh'' CTLs should result in similar anti-tumour effects) or are due to resistive mechanisms deployed by the tumor (in which case a second transfer of CTLs would be expected to provide only limited benefit).
In our analysis, we estimated several parameters, and these can be compared to values in the literature. A wide range of estimates (0.02-1 day À1 ) for the death rate of CTLs have been reported, 20,36-38 which are compatible with our own estimate of 0.8 days À1 . We also estimated that the infiltration rate of CTLs into the tumor was very low, with our estimate hitting the lower parameter bound of s 0 = 0.001 CTLs mm À3 day À1 . The lower bound of s o was set on the basis that tumors were around 100 mm 3 while CTLs were infiltrating, and one could interpret the value s 0 = 0.001 CTLs mm À3 day À1 as a 10% chance of a single CTL infiltrating a tumor of size 100 mm 3 per day. This value seems implausibly low, thus our modeling of CTL infiltration into the tumor as a constant rate process may be an oversimplification. Since s 0 was low, the large numbers of CTLs in the tumor were instead explained by a significant contribution from CTL expansion inside the tumors (estimated as the high value of s e = 8-9 days À1 ). CTLs have been reported to divide up to five times per day in vivo, 39 which is in a similar range albeit somewhat lower than the 8-9 divisions per day we estimated. While the variability in tumor growth rates estimated for individual untreated mice ( Figure S1) does increase uncertainty regarding estimates of other parameters, the early volume growth rates in the presence of CTLs matched the population growth rates in the absence of CTLs quite well, so we expect the impact of this variability to be low. In conclusion, there remains some uncertainty surrounding our parameter estimates for CTL infiltration and intratumoral expansion, underscoring the need for future experimental studies to measure the rates of both processes following adoptive CTL transfer. Importantly, we do not expect these uncertainties to affect our main conclusions regarding CTL killing and antiproliferative effects or the importance of various checkpoints. This is because those conclusions require the model to accurately describe the number of CTLs inside the tumors over time, which was the case.
In our analysis, we used mRNA expression as a substitute for protein expression. Previous studies report that mRNA levels are substantially predictive of protein expression levels. 40,41 Moreover, although some delay should be expected between mRNA expression and protein expression, this delay has been estimated to last for only a few hours 41  iScience Article of measurements made across several days. For the in vivo setting we studied, the rapid decline in the S-G 2 -M:G 1 ratio after transfer of CTLs indeed suggests that protein expression rapidly follows mRNA expression. Conversely, the S-G 2 -M:G 1 ratio does not appear to recover immediately upon downregulation of Ifng mRNA. One explanation could be that the effect of IFNG lasts longer than the protein due to downstream signaling. Another is that tumor cells are very sensitive to low levels of IFNG, therefore the effect could persist even after IFNG synthesis has substantially declined. The latter explanation seems to be in line with a study, which found that bystander sensing of IFNG could occur at distances of over 40 cell lengths, 13 implying high sensitivity of tumor cells to this cytokine.
In summary, we have presented a mathematical model that can successfully predict inhibition of tumor growth following adoptive T cell transfer. We used this model to quantify the contribution of IFNG and cytotoxicity to the anti-tumor activity of CTLs, which led to the conclusion that IFNG contributes most to tumor growth blockade by CTLs. Our model also includes anti-tumorigenic (antiproliferative, enhancing recruitment of CTLs) and pro-tumorigenic (driver of CTL exhaustion) effects of IFNG. The presence of opposing effects of IFNG has led to descriptions of an ''IFNG paradox''. 12 Our model, by including these different effects associated with IFNG, can serve as a quantitative baseline to be augmented in future and may help guide further experimental work.

Limitations of the study
Our study was limited by a lack of direct data concerning several important aspects of the CTL dynamics within the tumor. First, we had no direct data on the killing rate (k e ) of the CTLs inside the tumor, so this parameter was allowed to vary freely during the fitting process. Depending on which immune checkpoint molecules were included in the fitting process, we recovered a range of different values for k e , although these parameters were all plausible and comparable with other values for the killing rate of tumor cells by CTLs in vivo reported elsewhere, 42 including that for attack of B16F10 cells. 19 Interestingly, the higher killing rate estimated using the model with immune checkpoints implies that the low predicted impact of killing on the tumors in these experiments might have been due to checkpoint-mediated suppression of killing. This shows that it is difficult to obtain proper ''intrinsic'' killing rate estimates without direct measurement. Importantly, however, the choice of immune checkpoint molecules and the resultant values of k e did not impact our conclusion that IFNG-mediated cell-cycle arrest was the main determinant of tumor control. A second limitation surrounds our model of CTL exhaustion inside the tumors. Two specific questions we could not address due to the whole-tumour microarray data we used were: 1) which cells were expressing inhibitory molecules and 2) whether our results would have been different had we included ''missing'' relevant molecules from the transcriptomics data, e.g., the immune checkpoint CTLA-4. Unbiased gene expression data generated at the single-cell level, using single-cell RNA sequencing techniques, would therefore be interesting to incorporate into similar modeling strategies in future. A third limitation is the possible presence of other tumor-infiltrating cells, such as MDSCs which are recruited to B16F10 tumors after adoptive transfer of CTLs and exert suppressive effects on the tumor-infiltrating CTLs. 43 The frequency of various immune cell types can be inferred from either single cell or bulk transcriptomic data using computational methods, 44 and it would be interesting to extend our modeling approach to include other relevant immune cells using such methods in future. Finally, our estimates of the numbers of tumor cells and CTLs inside the tumor are based on thin 2D tumor sections. Since CTLs are smaller than tumor cells, they are more likely to remain undetected within the thin tissue slices, possibly resulting in an underestimate of the CTL to tumor cell ratio. It should be noted though that with image-based approaches typically more T cells are detected compared to approaches involving cell isolation procedures. 45 Nevertheless, using images based on 3D image stacks to determine the densities of cells inside tumors would likely result in more accurate estimates in future.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

Analysis of gene expression data
Microarray data were downloaded from the Gene Expression Omnibus (GEO) database (series GSE57304; samples GSM1379331-GSM1379344). These data correspond to the same set of experiments as the image and tumor volume progression data we have used, and the methodology for acquisition of these data has been described previously. 15 Briefly, tumor tissues from mice were harvested on different days (1, 3, 5, 7) after CTL transfer, or on the same days in the untreated (without CTLs) condition. Each sample contained 500 ng of pooled RNA from 3 to 4 different tumors, and microarray analysis was performed with 45,018 probes to quantify expression levels of the targeted genes. We performed similar data processing steps to the original publication: probes were discarded when their gIsWellAboveBG flag was zero at all samples, and we normalised different samples at the 75th percentile.

Basic ODE model
We The resulting tumors grow exponentially when the ratio of cells in G and S states is at its steady state value. When CTLs (E) are introduced into the tumors, our basal model of tumor growth (Equations 1 and 2) is modified to include two possible effects CTLs can have on the tumor. The first of these effects is direct killing of tumor cells, which occurs with a constant rate k e (per CTL). As we have done previously, 18,19 we take the total killing activity of CTLs to be directly proportional to the number of CTLs inside the tumor, such that the total killing activity of the CTLs is given by ak e E. Here, a is a scalar used to modify the effector functions of CTLs if their activity is reduced due to being exhausted (see . Note that in general the total amount of CTL killing is expected to saturate with CTL and/or target cell density, and can be described with a double saturation (DS) function. However, in a situation where the target cell density greatly exceeds the CTL density, as is the case in our modeled scenarios, this DS function reduces to the here employed (simpler) killing term. [18][19][20] We consider that killing is directed equally toward cells in G 1 or S-G 2 -M phases, so that the fraction of tumor cells in either state (i.e. G ðS + GÞ À 1 or S ðS + GÞ À 1 , respectively) determines the fraction of the total killing activity that each subset of tumor cells receives.
The second effect that CTLs can have on the tumor is an antiproliferative effect, mediated by IFNG, which results in an arrest of the cell cycle in the G 1 phase. To include this effect in our model we reduce the transition of cells out of the G 1 phase by scaling with the term ð1 + k i I=V Þ À 1 . Here, the variable I represents the total quantity of IFNG inside the tumor and V is the variable representing tumor volume. Thus, the term I/V represents the concentration of IFNG inside the tumors, and k i determines the concentration dependence of the IFNG dependent reduction in the rate at which tumor cells can leave the G 1 phase. The equations to describe the evolution of the number of tumor cells in G 1 or S-G 2 -M phases become: dS dt = k gs G ð1 + k i I=V Þ À 1 À k sg S À ak e E S ðS + GÞ À 1 ; iScience Article dG dt = À k gs G ð1 + k i I=V Þ À 1 + 2k sg $ S À ak e E G ðS + GÞ À 1 :

(Equation 4)
To test the compatibility of Ifng transcription dynamics with G 1 tumor cell-cycle arrest Equations 3 and 4 were used directly. I and E were estimated by linearly interpolating between the mean of the experimental data at each available time point, and these linear interpolations were used as inputs to the model.

ODE model with CTL dynamics
To describe the dynamics of the CTL population and their production of IFNG, we extended our basic ODE model with further equations. We consider that after transfer, CTLs would begin to arrive in any given region of the tumor at a constant rate s 0 . We take a constant rate of CTL arrival per unit volume of tumor, hence the rate at which CTLs can find the tumor scales with tumor volume. Moreover, we consider that CTLs expand within the tumor at rate s e , and die at a constant rate, d E . CTL expansion inside the tumor is also reduced according to the level of CTL exhaustion (a). Thus, CTL dynamics is described by the following equation: We consider CTLs to be the major source of IFNG inside the tumors, therefore IFNG production is proportional to the number of CTLs, but is reduced according to their level of exhaustion (a), and IFNG disappears from the system with a rate d i : Finally, we include a mechanism whereby CTLs become exhausted inside the tumor. T cell exhaustion is characterised by a loss of effector functions along with a progressive increase in the amount and diversity of inhibitory receptors expressed by T cells. 25,29,30 We used the well described PDCD1, CD274, LAG3 & HAVCR2 inhibitory molecules as indicators of exhausted T cells, 26 which in our model appear with variable names P, P L , L, and H (respectively): ) We tested several model variants, one with no immune checkpoints, three in which we consider one checkpoint at a time, and one considering all checkpoints simultaneously. All inhibitory molecules follow similar dynamics, increasing in proportion to the number of CTLs inside the tumor and disappearing from the system with different rate constants d p ,d pl ,d l , and d t (respectively). Production is increased proportional to the term ð1 + k A I =V Þ, which allows for a contribution of IFNG to the exhausted state. Note that the IFNG induction of exhaustion may be direct or indirect, e.g. by increasing antigenicity of tumor cells and thereby increasing stimulation of T cells via the T cell receptor. In our model, we consider PDCD1, LAG3, and HAVCR2 as expressed on the membrane of CTLs, so the ratio of each of these checkpoint molecules to the number of CTLs determines the overall level of exhaustion of the CTLs in our model. CD274 is modeled differently, being a ligand for the PDCD1 receptor, and the concentration of CD274 in our model is multiplied together with the membrane density of PDCD1 expressed on CTLs to determine the contribution from PDCD1-CD274 signaling (see Equation 11 in ref. 20 ). To describe the joint effect of these inhibitory molecules on the CTLs, we consider a weighted sum R: iScience Article toward the level of exhaustion of CTLs. In absence of detailed information about the impact of exhaustion level on CTL functions (killing, IFNG production, expansion), we take all these functions to be equally reduced with the level of CTL exhaustion: Here, k ex is the level of exhaustion at which all effector functions are half of their maximum value. Thus, Equation 12 scales the exhaustion level to a scaled term a which can range from 0 to 1 and is applied to the relevant rate constants in Equations 3, 4, 5, and 6.

Parameter estimation
Parameter estimation for the two parameters in the basal tumor growth model was performed separately from the other parameters. We consider that the density of the tumor cells remained constant over time, so that volumetric tumor growth rate could be taken as a proxy for the expansion rate of the tumor cell population. We obtained an estimate for the untreated tumor growth rate (g) from fitting an exponential model of tumor growth to the volumetric growth data for the untreated tumors. Moreover, estimates of the ratio of tumor cells in the S-G 2 -M: G 1 phase gave a second measurement allowing the two-parameter (k gs , k sg ) basic model of tumor growth to be completely defined, considering that the ratio of S-G 2 -M: G 1 phase tumor cells has reached a steady state (which is reasonable since we deal with data two weeks after tumor inoculation). Then, using the equation for exponential growth: dðS + GÞ dt = g ðS + GÞ;

(Equation 13)
where g is the tumor growth rate, one can substitute the left hand side of Equation 13 with k sg S, noting that k sg S is the total rate of new tumor cell production obtained from summing Equations 1 and 2. Following the substitution, an expression for g can be found in terms of S and G: The resulting Equations 17 and 18 can be combined to remove the common terms (ge gt ). This leads to the following: where we have omitted the subscript with the understanding that Equation 19 is valid only when the tumor is growing exponentially with the ratio S/G at a steady state value. Therefore, by substituting the expressions for dG dt and dS dt given in Equations 1 and 2, the Equations 14 and 19 can be rearranged to express the k sg and k gs parameters for the basal tumor growth model as a function of tumor growth rate (g) and the ratio of S-G 2 -M: G 1 phase tumor cells at steady state: k sg = g ð1 + S = GÞðS=GÞ À 1 ; (  46 This revealed that the parameters k A and k ex were not identifiable, so we set both values equal to 1. The remaining model parameters which relate to the dynamics of the CTLs and their effects on the tumors were obtained together, by fitting to all available data over time simultaneously (i.e. the CTL counts; the S/G ratios; the volumetric growth data; the Ifng expression data; the immune checkpoint expression data). To determine the optimal parameters we maximised the log likelihood: Where y m ðt j Þ is the model estimate for observable m output at time point t j , y m;j and s 2 m;j are respectively the mean and standard deviation of the corresponding experimental measurements. For the immune checkpoints there was only one sample per time point corresponding to the CTL treated condition, precluding calculation of the standard deviation. Therefore we estimated the standard deviation of the checkpoints based on the standard deviation of a list of housekeeping genes 47 (most stable mouse transcripts). A linear regression of the mean and standard deviation expression of the housekeeping genes in our data gave the relation: s = 0:054 + 0:29y, with adjusted R-squared 0.986. This relation was used to estimate the standard deviation for each individual immune checkpoint. Growth rate measurements between days 10-14 had an outsized effect on our optimization due to their very low standard deviation originating from only two mice in a single experiment, thus were excluded from fitting. For the basic model with I and E used as inputs, there were only two parameters (k i and k e ) to be estimated, therefore we tested all combinations of the parameters in the range 0-20 at intervals of 0.1 and selected the combination with the highest log likelihood. For the model including CTL dynamics we used a multi-start strategy with 20,000 initial parameter guesses selected using latin hypercube sampling, then optimised using the L-BFGS-B algorithm. We set relatively permissive bounds on the parameters as follows: all protein degradation parameters (d i , d t , d l , d p , d pl ) between 0.01 and 100 days À1 , corresponding to half-lives between 10 min and 70 days. Bounds of 0.01-100 days À1 were used for the killing rate of CTLs, and bounds of 0.01-10 days À1 were used for CTL death rate. All other parameters were bounded in the range 10 À3 to 10 3 . Best-fitting parameter sets are provided in Table 1 in main text.
To approximate the posterior distribution of the fitted parameters, we used the Adaptive Metropolis Parallel Hierarchical Sampling algorithm. 48 We ran 20 auxiliary chains alongside the main chain, using the best 21 parameter sets derived from the multi-start optimisation procedure for initialisation. Chains were run for 100,000 steps and the final 50,000 samples of the main chain were used for sampling.