A Validated Mathematical Model of Cell-Mediated Immune Response to Tumor Growth

Therapeutic vaccines are being developed and tested as a promising new approach to treatment for cancer patients. There are still many unanswered questions about how the immune system interacts with a growing tumor, and which components of the immune system play significant roles in responding to immunotherapy. Mathematical models provide an analytic framework in which to address such questions, and these models can be used both descriptively and predictively. It is important to develop models of tumor growth that include a representation of an immune response. The ultimate goal is to create models that can reflect a system's response to emerging biological therapies, such as vaccine therapy. Mathematical modeling of tumor growth and treatment has been approached by a number of researchers using a variety of models over the past decades (for overviews, see for example refs. 1–5).

Our new mathematical model of tumor-immune interactions sheds light on the differing roles of the natural killer (NK) and CD8⁺ T cells in suppressing various tumor cell lines in mice and humans. NK cells are large granular lymphocytes that do not express markers of either T- or B-cell lineage. A constituent of innate immunity, they recognize and destroy tumor cells, among others, independent of prior exposure. NK cells are thought to play a key role in preventing the development of clinical cancer by killing abnormal cells before they multiply and grow.

T cells, which carry the CD3⁺ marker, are morphologically small lymphocytes in the peripheral blood. They develop in the thymus and orchestrate the immune system response to infected or malignant cells. CD3⁺CD8⁺ T cells (also called CD8⁺ T cells) are a critical subpopulation of T-lymphocytes which can be cytotoxic to tumor cells provided previous sensitization has occurred.

Our model is based on and validated by recent experimental studies by Diefenbach et al. (6), in which mouse tumor cell lines are modified to express higher levels of immune stimulating NKG2D ligands. We further validate our model using additional human data provided by Dudley et al. (7), in which subjects with metastatic melanoma are treated with highly selected tumor-reactive T cells. Both the mouse and the human studies provide vital experimental information about tumor growth rates and effector to target lysis rates. The model is used to explore the dynamics of tumor rejection, the specific role of the NK and CD8⁺ T cells, and the development of protective immunity to subsequent tumor rechallenge.

The mathematical structure of the model is built on earlier modeling work from refs. 8 and 9 in which tumor growth, an immune response, and chemotherapy treatment are represented by a system of differential equations. Our current model also extends other lower-dimensional models, such as that described in ref. 10, in which different cell populations are represented as interacting species. Other mathematical models that include an immune interaction with a tumor are described in refs. 3, 11–24.

The dynamics of our mathematical model, as well as certain parameter values, are borrowed from assertions and data provided in the following two publications. A brief summary of each publication follows.

Mouse data (6). In this study, the authors show that ectopic expression of certain murine NKG2D ligands in several tumor cell lines resulted in rejection of the tumors by syngeneic mice, and that this rejection was mediated by both NK cells and CD8⁺ T cells. A retrovirus expression system was used to ectopically express high levels of these ligands in thymoma, T-cell lymphoma, and melanoma. Groups of mice were inoculated s.c. with syngeneic tumor cell transductants. The experimental results in this work showed that if expression of these ligands was low, tumor establishment, growth, and successful evasion of immune responses in mice took place. However, sufficiently high levels of ligand expression created a significant barrier to tumor establishment in the mice. The authors stated that there is hope that these treatments with ligand-transduced tumor cells might eventually lead to the ability to effectively eliminate preestablished tumors, as well as naturally arising tumors. The main conclusions of the experiments in ref. 6 can be summarized as follows: (a) NKG2D ligand expression (at sufficiently high levels) causes activation of CD8⁺ T cells and NK cells. (b) Ligand-transduced tumor cells can stimulate protective immunity to tumor rechallenge with ligand-negative tumor cells. (c) Tumor cells engineered to up-regulate ligand expression are much less likely to establish tumors than wild-type low-ligand tumor cells. (d) Typical levels of NKG2D ligands naturally found on most tumor cell lines are low and thus insufficient to trigger host response. The hope is that tumor immunity can be boosted by engineering cells with higher ligand levels.

Human data (7). Patients with metastatic melanoma were treated with the adoptive transfer of highly selected tumor-reactive T cells directed against overexpressed self-derived differentiation antigens. In previous trials, increasing the number of circulating CD8⁺ CTL precursor cells showed no correlation with tumor regression. However, lympho-depleting chemotherapy has been shown to have a positive effect on the efficacy of T-cell transfer therapy, perhaps because lympho-depletion disrupts the regulating mechanisms of the immune system. The patients in this study received immuno-depleting chemotherapy for 7 days before being treated with tumor-infiltrating lymphocytes (TIL) expanded in vitro. Six of the thirteen patients in the study showed objective clinical responses, with significant shrinkage of metastatic deposits. Two of these patients were selected as exhibiting particularly strong immune responses. After validating the functional forms of our mathematical model against the mouse data in ref. 6, we tested the generality of our model by comparing mathematical predictions to the human data in this published study. We used data collected from the two strongly responding patients to validate our model because careful measurements of lymphocytosis were made in both cases (see Fig. 2). Data on the immune responses of less strongly responding patients were not available. An intriguing question to ask in future work when such data do become available is whether the dynamics predicted in this model can account for the nonresponding or low-responding cases.

The specific biological assumptions we took into account when developing our model equations are based on both accepted knowledge of immune system function and conclusions stated in refs. 6 and 7. The assumptions include:

• (1) The tumor cells grow logistically in the absence of an immune response. This assumption is based on previous work, and is also suggested by graphs in ref. 6 of tumor growth in mice without an immune response. Other mathematical forms, such as Gompertzian growth and power-law growth (see ref. 25), have been considered in other contexts, and can also be made to fit these data. However, we chose here the logistic form because it was the simplest form for which predictions also agreed with the experimental data.

• (2) Both NK cells and CD8⁺ T cells can kill tumor cells.

• (3) Tumor cells have the potential to engender cytocidal activity in previously naive and noncytotoxic cells. Both NK cells and CD8⁺ T cells respond to tumor cells by increasing metabolic activity and releasing various lymphokines, including IFN-γ. Note that the level of effector cell “effectiveness” depends on both the number of cells present as well as the cytotoxicity of individual cells. In the model, we do not separate the measures of high effectiveness per cell from an increase in cell population, but measure the combined overall increase in effectiveness in response to tumor.

• (4) As part of innate immunity, NK cells are always present and active in the system, even in the absence of tumor cells.

• (5) As part of specific immunity, tumor-specific CD8⁺ T cells are recruited once tumor cells are present.

• (6) Each NK cell and CD8⁺ T cell will eventually become inactivated after some number of encounters with tumor cells.

In the equations, we denote the three cell populations by

• T(t), tumor cell population at time t;

• N(t), total level of NK cell effectiveness at time t;

• and L(t), total level of tumor-specific CD8⁺ T-cell effectiveness at time t.

Using the list of assumptions from above, we describe the system as three coupled differential equations, where each equation gives the rate of change of the particular cell population in terms of growth and death, cell-cell kill, cell recruitment, and cell inactivation. In particular,

The mathematical forms of the growth and death terms for tumor and immune cell populations will reflect assumptions (1), (4), and (5). Assumption (2) is reflected in the cell-cell kill term; assumption (3) gives rise to the effector cell recruitment terms; and assumption (6) is incorporated through the effector inactivation terms.

Immune recruitment terms are generally assumed to be of a Michaelis-Menten form (see, e.g., ref. 10 in which Michaelis-Menten dynamics are derived for immune cell recruitment by cancer cells). These dynamics are commonly used in mathematical tumor models that include an immune component because they allow for a saturation effect (see, e.g., ref. 19). In the case of the CD8⁺ T cells, in addition to being recruited by interactions with T-cell processed tumor cells through a Michaelis-Menten dynamic, additional CD8⁺ T cells are stimulated by the interaction of NK cells with tumor cells. This NK stimulation is represented by the rNT term in Eq. (C). The term rNT, representing interactions between NK cells and tumor cells, is the vehicle through which we model the fact that the specific immune response of the CD8⁺ T cells is activated only after the activation of the earlier response of innate immunity.

To determine the fractional cell kill dynamics, data from chromium release assays published in refs. 6 and 7 were used. Chromium release assays determine the ability of CD8⁺ T cells to lyse target cells expressing specific ligands. The assays in both refs. 6 and 7 were standard 4-hour ⁵¹Cr release assays. Standard techniques exist for collecting, storing, and coculturing patients' immune cells with tumor cells, a procedure which can be implemented before the onset of treatment, or anytime thereafter. The lytic activity of these cells can then be analyzed with the assay (see, e.g., approaches referenced in ref. 7). The fractional cell kill term for the NK cells was assumed to be proportional to the size of the NK cell population. This assumption was consistent with all the data we examined. However, the same assumption was not consistent with the data for tumor-specific CD8⁺ T cells.

We therefore introduced a new functional form for the (CD8⁺ T)-tumor kill term, represented by D in Eq. (A), given explicitly in Eq. (D) below. This term we say is of a “rational” form. Because this term is an innovation in this model, its derivation is fully discussed in “Fractional cell kill,” and results using this form are presented in “New functional forms.”

Substituting specific mathematical forms for each of the growth, death, recruitment, and inactivation terms yields the following system of equations:

where

Table 1 provides a detailed listing of the parameters in this model, along with their units, descriptions, numerical values for the simulations, and reference sources from which these values were taken. Detailed development of all terms, except for the new fractional cell kill term D discussed below, can be found in ref. 8.

A crucial component of the model equations is the set of terms describing the interaction or “competition” between the tumor cells and either the NK cells or the CD8⁺ T cells. It is common to assume that that the action of the effector cells is to reduce the tumor cell population by a term proportional to both the effector and tumor populations (see, for example, refs. 1, 26, 20, 16, 10, 8). Mathematically, this is expressed using the simplest product form for the competition between effector cells and tumor cells. In our case, the effect of the NK cells on the tumor cell population would be expressed with the term −cNT, whereas the form of the CD8⁺ T-cell effect would be −dLT, where c and d are proportionality parameters determined through experiment. As discussed below, we found that although this simple linear product term was sufficient to reproduce in simulations the experimental interaction of NK cells with tumor cells, it did not accurately reflect the action of the CD8⁺ T cells.

We did data fitting experiments using a generalization of the product term: a power term that allows for exponential growth. Using cell lysis data from ref. 6, we employed an iterative process to find the parameters c and v in cN^vT, and d and λ in dL^λT. The procedure involved minimizing the distance between the data points and the predicted percent lysis curves generated by the model over a range of c, v, d, and λ values. For each [c, v] and [d, λ] pair, a prediction was made by solving a system of differential equations up to time T_final = 4 hours, with initial values from the effector/target ratio data in ref. 6. The optimization was done using an iterative procedure to determine the least-squares fit to the data, again through repeated solutions of the differential equations. We found that the best-fit exponent for the NK kill term was v ≈ 1. Because a good mathematical model will be one in which the desired behaviors of the system are captured using the simplest mathematics possible, we chose to keep the product form −cNT to describe the effect of the NK cells on tumor cells. In fact, the optimal value of c, determined using our algorithm, reproduced the lysis rate data extremely well (see Fig. 1, top).

However, when fitting for variables d and λ for the CD8⁺ T-cell kill term, we found that the power form produced growth curves for T that were not particularly good fits to the data provided in ref. 6. Instead, we found that we could produce curves that better fit the data by allowing this term to have the rational form given in Eq. (D), for which we also had to determine parameter s. In Eq. (D), the exponent λ represents how the lysis rate depends on the effector/target ratio, the parameter s affects the steepness of the curve, and parameter d gives the maximum lysis rate. Further comment on the fits to the data with this rational law dynamic is given in Results. We note that the additional parameter in Eq. (D) gives three degrees of freedom, so that a better fit to the data should be expected using the rational form. However, because the observations in ref. 6 give five data points for each cell type considered, the closeness of fit to the data supports the idea that the form of this term is correct. In particular, both in vitro and in vivo experiments indicate that percent lysis seems to be a function of the ratio of CD8⁺ T cells to tumor cells, explaining the appearance of (L/T). Furthermore, the data indicate that the percent of cells lysed never exceeds a maximum, a saturation effect that is reflected by the rational form given in Eq. (D).

This saturation effect highlights the fact that the NK cells and CD8⁺ T cells are interacting with tumor cells in a qualitatively different way because there is no saturation level for the NK cell competition term. It may be that the NK cell-kill rate could achieve saturation as well in theory, but in practice this does not occur. On the other hand, it may be that the antigen-specific T cells follow this curve to saturation because they are targeting a specific tumor type, and are therefore more effective in terms of cell-cell interactions.

From ref. 6 we were able to get data on the growth curves in the absence of an immune response, which allowed us to estimate parameters a and b. These model parameters were estimated from the data in ref. 6 by minimizing a least-squares distance using optimization software built into MATLAB6. Data measuring the percent of IFN-γ producing immune cells as a function of ligand expression allowed us roughly to estimate immune recruitment rates stimulated both by ligand-transduced and control-transduced tumor cells. Other parameters, such as the background source rate for NK cells (σ) and death rates for immune cells (f and m), were taken from the literature (e.g., refs. 10 and 27). Although some of these variables are rough estimates, and may deviate from other specific data, the model as a whole qualitatively describes the observed data both in the mouse and in the human experiments.

The model can be used to simulate the effect of enhancing ligand expression on tumor cells by allowing the relevant parameters to depend on the tumor cell type. The relevant parameters in this model are c and d, the effectiveness of the immune cells, along with g and j, the recruitment variables.

Figure 1 (top) plots the effector/target lysis data from ref. 6 for NK cells, along with our simulated model curves. The ligand-transduced tumor cells are lysed at a higher rate by NK cells than those that are control transduced. The two values of NK-lysis parameter c estimated from the two sets of data accurately reproduce the effects of this ligand transduction.

In Fig. 1 (middle), effector/target lysis data and simulations for the CD8⁺ T cells are presented. For our experiments, four CD8⁺ T-cell lysis parameters were determined through fitting to the four ligand transduction data sets of ref. 6, and these are all able to capture the different experimental outcomes. For brevity, only the two cases representing priming and rechallenge with control-transduced cells and priming and rechallenge with ligand-transduced cells are presented in Fig. 1. Figure 1 (middle) plots the experimental data against the mathematical model prediction using the best-fit parameter values for both the power form and the new rational form of the competition term. Note that in Fig. 1 (middle, left) in which we compare fits to data for non-ligand-transduced cells, although the difference between the fit achieved by the traditional power kill law and that achieved by the new rational kill law is not clearly visible, the numerical difference in the error term is present. This can be seen in Fig. 1 (bottom, left). Here, we plot the numerical errors between the predictions and the data, allowing a comparison between the goodness-of-fit of the power form and the goodness-of-fit of the rational form of the competition term. In the right panel of the center row, the superiority of the fit achieved by the rational kill law over the power kill law is visible and striking. Similarly, the numerical error bars of the bottom panel, right side, reflect the much smaller error achieved by the rational kill law. It seems that it is critical to employ the rational law to fit ligand-transduced cell data, whereas the use of either the rational or the power law for non-ligand-transduced data will give us an acceptable fit. This may indicate that the more effective the immune cells are at lysing their target cells, the more they follow a rational law dynamic.

New functional forms. In contrast to NK cells, CD8⁺ T cells have to be primed to be activated. Therefore, for the CD8⁺ T cell to tumor cell competition term, D, there were four different data sets to be fitted. These came from first priming with either control-transduced cells or ligand-transduced cells, then challenging with either control-transduced or ligand-transduced cells. Experimental data indicated that the least effective combination was both to prime and to challenge with control-transduced cells, in which case there was at most 10% lysis. On the other hand, the most effective combination was to prime and then challenge with ligand-transduced cells, resulting in a maximum of nearly 70% lysis. In the model, these cases are distinguished by the estimates of the parameters d and λ, which changed by a factor of up to 5. Thus, the estimated values of these parameters as components of the rational cell kill term D may be highly significant factors distinguishing between tumors that are weakly and strongly antigenic.

The form for the rational competition term D is mainly phenomenological in the sense that it models observable outcomes, not direct underlying mechanisms. It is not immediately clear what the individual components of D biologically represent. The use of phenomenological dynamics in modeling biological processes is quite common and can serve to provide predictive capabilities in the model. Such descriptive (as opposed to explanatory) dynamics are frequently used as a foundation on which to build models of tumor development. For example, see the comparison of several phenomenological tumor growth models presented in ref. 28 (p. 239). Perhaps future investigations may elucidate the underlying mechanisms that give rise to the form of D.

Validation of model with human data. To validate the fundamental model dynamics with respect to the new rational form of the tumor-specific cell lysis term, we did another comparison of power-law versus rational law predictions, this time using human (CD8⁺ T)-tumor lysis data. Figure 2 shows the results of this comparison. The top graph shows the power-law predictions plotted against (CD8⁺ T)-tumor lysis data for two separate patients. It is clear that the power-law prediction does not fit the data particularly well. On the other hand, the bottom graph shows the prediction using our newly introduced rational law. In this case, the model can predict cell lysis quite accurately, even when applied to this human data set.

For this particular set of data, effector cells are fairly efficient at lysing tumor cells, with a maximum lysis rate of around 60%. Note that, as with the ligand-transduced mouse data, the difference between the power law and rational law models is quite pronounced, once again indicating that the rational law model is particularly well suited to simulating cases in which effector cell lysis rates are relatively strong.

It is necessary in each case to find the parameters which will describe the particular type of tumor-immune interaction under study. The two data sets pictured here underline a feature inherent in the modeling process: there is a wide variety of cell behavior between any two different patients. Care must therefore be taken in making sweeping statements about specific responses to treatments, and any quantitative information must be interpreted as one possibility, and not as a firm predictor in any given case. However, a large set of simulations, along with some analysis of the sensitivity of the model to parameter fluctuations, can certainly provide a general picture of possible behaviors under certain conditions. Further comparisons may lead to new insights into the nature of the differences between different tumor types, as well as different immunotherapeutic protocols.

Sensitivity analysis using human variables. To discover which components of the model contribute most significantly to determining final tumor size, we did a sensitivity analysis. Model sensitivity was assessed by measuring the effect of small variable changes on the final volume of the tumor as represented by a simulation of system evolution over 25 days. Because ultimately we are interested in predicting a patient's response to immunotherapy treatment, we used human data for the sensitivity study. In particular, the variable set from patient 9, available in ref. 7, and for whom lysis data are plotted with squares in Fig. 2, was used as the base point. Each parameter was perturbed from its estimated value by 1%, and the corresponding percent change in final tumor volume was calculated.

The results of this parameter sensitivity analysis for the mathematical model are shown in Fig. 3. The system in this case is found to be most sensitive to the exponent in the CD8⁺ T lysis term, λ, as well as to the tumor growth variable a. This suggests that, in addition to the aggressiveness of the tumor, as represented by growth variable a, even very small changes in the cytolytic effectiveness of tumor-specific T cells, as represented by shifts in the value of λ, can affect clinical outcome. This would indicate that any treatment which might enhance this effectiveness should aggressively be pursued. By contrast, the size of the tumor after 25 days is not very sensitive to the NK cell competition variable c. According to this model then, the cytolytic activity of the NK cell population alone is not a determining factor in the eventual size of the tumor, and should be considered in conjunction with CD8 cell activity.

Model simulations with mouse variables: vaccination and thresholds for immune efficacy. The simulations show what this model would predict under three different experimental scenarios similar to those reported in ref. 6. These simulations explain some of the reported experimental observations (see ref. 6, Figs. 2 and 3, p. 167-8). Ligand-transduced cells stimulate the immune response sufficiently to control tumor growth (Fig. 4, top right), whereas control-transduced tumor cells escape immune defenses (Fig. 4, top left). In Fig. 4 (top left), the immune system is rechallenged at day 10 after priming with control-transduced cells, and the tumor escapes surveillance. In Fig. 4 (top right), the immune system is again rechallenged at day 10 with control-transduced cells, but the primary challenge was with ligand-transduced cells. This simulation shows that the tumor is controlled, indicating the development of immunity. Changing ligand levels on the cells requires changes to the model parameters d, λ, and s (all the parameters involved in the rational T-cell kill term D), as well as c (strength of NK cell kill), g (NK cell recruitment rate), and j (CD8⁺ T-cell recruitment rate). Numerical values for these parameters with varying ligand levels are provided in Table 1.

Simulations generated by a validated mathematical model can be used to detect thresholds for immune efficacy. In Fig. 4 (bottom), we reproduce with a computational solution of our mathematical model the qualitative results of three sets of experiments that were presented in Fig. 2 of ref. 6 (p. 167). For the experiments in ref. 6, groups of mice were challenged with either 10⁴, 10⁵, or 10⁶ ligand-transduced tumor cells, then tumor establishment was tracked. For our in silico simulations, we also challenge the mathematical system with these three levels of tumor cells. Figure 4 (bottom left) shows simulated tumor cell growth over time in response to these three initial levels of tumor burden in the absence of CD8⁺ T-cell activity, reflecting the experiments in which the mice were depleted of CD8⁺ T cells. This simulation represents a system lacking a strong antigen-specific immune response. The system can control a small tumor, but tumor challenges of 10⁵ cells or more escape immune system control.

Figure 4 (bottom center) shows simulated tumor growth outcomes for the same three experiments done in the absence of NK cells, reflecting the experiments with mice depleted of NK cells. The system is now able to control initial tumor burdens of up to 10⁵ cells, but a larger challenge of 10⁶ cells escapes immunosurveillance.

Figure 4 (bottom right) shows simulated results with both NK and CD8⁺ T cells active, reflecting the experiments on mice with intact immune systems. With both the NK cells and the CD8⁺ T cells working together, initial tumor burdens of up to 10⁶ cells are controlled.

This model incorporates tumor-immune interaction terms of a form that is qualitatively different from those commonly used. This form provides a good fit with experimental data resulting from priming and rechallenge with different combinations of tumor cell types. The mathematical model suggests the value in continuing to research the mechanisms by which NK cells and CD8⁺ T cells induce tumor cell lysis. The functional descriptions presented in this article suggest that further laboratory tests may aid in determining why the two types of immune cells give rise to such different cell interaction dynamics. We hypothesize that the more effective the immune cell kill is, the more closely it follows a rational law dynamic.

The model currently includes no self-regulatory terms in the equations or down-regulation of an activated immune response. We have yet to resolve this issue in our model, but it is currently receiving active attention. Although the current model does not address all processes of the immune system, it does, despite its simplicity, fit the empirical data.

These experimental and simulated results, observed together with the cell lysis data presented and the parameter sensitivity analysis, highlight the importance of CD8⁺ T cell activation on final outcome. Model results seem to indicate that to promote tumor regression, it may be necessary (although perhaps not sufficient) to focus on increasing CD8⁺ T cell activity. In fact, we propose that there may be a direct positive correlation between the patient-specific efficacy of the CD8⁺ T cell response, as measured by cytotoxicity assays, and the likelihood of a patient favorably responding to immunotherapy treatments.

Here we list all of the parameters used in the model, their meaning and the estimated values. The first set of data in Table 1 reflects the experiments run to simulate the mouse experiments from ref. 6. The second set of data in Table 2 applies to the human data from ref. 7.

Grant support: L.G. de Pillis is supported in part by the W.M. Keck Foundation.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.