Cancer-induced immunosuppression can enable effectiveness of immunotherapy through bistability generation: A mathematical and computational examination.

Cancer immunotherapies rely on how interactions between cancer and immune system cells are constituted. The more essential to the emergence of the dynamical behavior of cancer growth these interactions are, the more effectively they may be used as mechanisms for interventions. Mathematical modeling can help unearth such connections, and help explain how they shape the dynamics of cancer growth. Here, we explored whether there exist simple, consistent properties of cancer-immune system interaction (CISI) models that might be harnessed to devise effective immunotherapy approaches. We did this for a family of three related models of increasing complexity. To this end, we developed a base model of CISI, which captures some essential features of the more complex models built on it. We find that the base model and its derivates can plausibly reproduce biological behavior that is consistent with the notion of an immunological barrier. This behavior is also in accord with situations in which the suppressive effects exerted by cancer cells on immune cells dominate their proliferative effects. Under these circumstances, the model family may display a pattern of bistability, where two distinct, stable states (a cancer-free, and a full-grown cancer state) are possible. Increasing the effectiveness of immune-caused cancer cell killing may remove the basis for bistability, and abruptly tip the dynamics of the system into a cancer-free state. Additionally, in combination with the administration of immune effector cells, modifications in cancer cell killing may be harnessed for immunotherapy without the need for resolving the bistability. We use these ideas to test immunotherapeutic interventions in silico in a stochastic version of the base model. This bistability-reliant approach to cancer interventions might offer advantages over those that comprise gradual declines in cancer cell numbers.

Introduction

Mathematical modeling of cancer-immune system interactions (CISI) can reveal the fundamental mechanisms that govern the dynamics of tumor growth (Altrock, Liu, Michor, 2015, Eftimie, Bramson, Earn, 2010), and represent and important tool to devise and test new forms of immunotherapy in silico (Talkington et al., 2018). The modeling relies on the appropriate integration of how cancer and immune cells affect one another (De Boer, Hogeweg, Dullens, De Weger, Den Otter, 1985, de Pillis, Radunskaya, Wiseman, 2005, Goldstein, Faeder, Hlavacek, 2004, Kronik, Kogan, Vainstein, Agur, 2008, Kuznetsov, Makalkin, Taylor, Perelson, 1994). Recent studies have uncovered a plethora of interactions by which cancer cells affect immune cells, and vice versa (Mellman, Coukos, Dranoff, 2011, Eftimie, Bramson, Earn, 2010). For instance, cancer cells elicit immune responses by a variety of effector cells (Parish, 2003, Smyth, Godfrey, Trapani, 2001, Mellman, Coukos, Dranoff, 2011). These effector cells, in particular white blood cells, natural killer cells (NKs) and cytotoxic T lymphocytes (CTLs) can lyse cancer cells (Quesnel, 2008), inhibiting tumor growth or even eliminating microscopic tumors altogether — a process termed immunosurveillance (Burnet, 1957, Burnet, 1967). However, cancers have also been shown to be able to suppress the proliferation of effector cells, which typically target cancer cells with specific biochemical signatures (Kooi, Zhang, Patenia, Edwards, Platsoucas, Freedman, 1996, Hamanishi, Mandai, Iwasaki, Okazaki, Tanaka, Yamaguchi, Higuchi, Yagi, Takakura, Minato, Honjo, Fujii, 2007). Cancer cells accrue mutations that, by changing these signatures, enable them to partially evade immune recognition (Altrock, Liu, Michor, 2015, Parsa, Waldron, Panner, Crane, Parney, Barry, Cachola, Murray, Tihan, Jensen, Mischel, Stokoe, Pieper, 2007, Hanahan, Weinberg, 2011). Furthermore, cancers may actively downregulate immune responses elicited against them (Keir, Butte, Freeman, Sharpe, 2008, Mellor, Munn, 2004, Aggarwal, Pittenger, 2005, Munn, Mellor, 2004, Marigo, Dolcetti, Serafini, Zanovello, Bronte, 2008), for example by recruiting the action of T regulatory cells (Mellman, Coukos, Dranoff, 2011, Ohta, Gorelik, Prasad, Ronchese, Lukashev, Wong, Huang, Caldwell, Liu, Smith, Chen, Jackson, Apasov, Abrams, Sitkovsky, 2006, Facciabene, Peng, Hagemann, Balint, Barchetti, Wang, Gimotty, Gilks, Lal, Zhang, Coukos, 2011), leading to immunosupression. A summary of these interactions shows that all combinations of stimulation and suppression on growth between cancer and immune cells may act simultaneously (see Fig. 1A). These interactions direct the interplay between the cancer and the immune system. Thus, their integration into mathematical models can reveal how immunotherapeutic approaches may be employed with maximum efficiency. The main immunotherapy approaches today work by impairing mechanisms that allow cancers to suppress immune action or by the administration of effector cells to the host (Dougan, Dranoff, 2009, Mellman, Coukos, Dranoff, 2011) (see Fig. 1B).

Fig. 1

Cancer-immune system interactions and effects of immunotherapy. A) Interactions governing the dynamics between cancer cells (T) and immune cells targeting the cancer cells —termed effector cells— (E). A complex web of interactions has been identified (Mellman, Coukos, Dranoff, 2011, Eftimie, Bramson, Earn, 2010, Altrock, Liu, Michor, 2015), with both cell types capable to both stimulate and suppress each others’ growth. B) Immunotherapy acts by either increasing the killing rate of effector cells, for example by administrating new effector cells into the host (adoptive T cell transfer), or by impairing the escape mechanisms cancer cells adopt to avoid being cleared, for example by monoclonal antibody therapy.

Due to their importance for immunotherapy, mathematical models of cancer-immune interactions have been the focus of intense theoretical efforts over the last decade (Altrock, Liu, Michor, 2015, Eftimie, Bramson, Earn, 2010). At the heart of this effort lies the task of identifying what features of cancer-immune dynamics are most effectively used to achieve immunotherapeutic success, i.e. cancer eradication. Following this tradition, we focus on deterministic, population based non-spatial models for analysis. We do this for two reasons. First, their mathematical simplicity is better suited to unambiguously reveal those model properties that are crucial for the modification of a cancer’s state, since the number of factors influencing the dynamics is manageable and understandable (Eftimie et al., 2010). Second, we can draw from a broad body of theory about immune system pressures developed within cancer research (Kuznetsov, Makalkin, Taylor, Perelson, 1994, Kuznetsov, Knott, 2001, d’Onofrio, 2005, d’Onofrio, 2008, Eftimie, Dushoff, Bridle, Bramson, Earn, 2011, Eftimie, Bramson, Earn, 2010, Lopez, Seoane, Sanjuan, 2014) as well as in virus dynamics, and in particular, human immunodeficiency virus (HIV) dynamics research (Nowak, May, 2000, Wodarz, 2007, De Boer, 2007, Althaus, De Boer, 2008, Conway, Perelson, 2014, Garcia, Regoes, 2015, Graw, Regoes, 2009, Graw, Richter, Oxenius, Regoes, 2011, Yates, Graw, Barber, Ahmed, Regoes, Antia, 2007, Althaus, De Boer, 2008, Elemans, Florins, Willems, Asquith, 2014, Garcia, Richter, Graw, Oxenius, Regoes, 2015, Regoes, Barber, Ahmed, Antia, 2007). To explore immunotherapeutic interventions in silico, we use a stochastic version of one of the models analyzed.

In this study, we have explored whether there exist consistent properties of cancer-interaction models that facilitate effective immunotherapy approaches. We investigated this question in a family of three related models of increasing complexity. To this end, we first developed a base model of cancer immune system interaction that captures some essential features of more complex models. The purpose of the model is to gain qualitative insights and then to serve as a guide for treatment therapies based on immune action. We analyzed under what parameter regimes the model produces biologically plausible behavior, and investigated how steady states are affected by changes in these parameters. We then successively extended the model to first include more complex features of cancer-immune system interactions (Conway and Perelson, 2014), such as saturating proliferative stimulation and exhaustion and in a second step, to include the action of NKs and CTLs. The common properties identified in these models were then used to study how combinations of immunotherapeutic treatments may work together to achieve eradication. To this end, we implemented stochastic simulations of the base model and analyzed how the dynamics are affected by adoptive T cell transfer (Rosenberg, 1991, Rosenberg, Yang, Restifo, 2004, Dudley, Wunderlich, Robbins, Yang, Hwu, Schwartzentruber, Topalian, Sherry, Restifo, Hubicki, Robinson, Raffeld, Duray, Seipp, Rogers-Freezer, Morton, Mavroukakis, White, Rosenberg, 2002, Rosenberg, Restifo, Yang, Morgan, Dudley, 2008), as well as by the disruption of immune evasion mechanisms of the cancer through for example monoclonal antibody therapy (Mellman et al., 2011; Brahmer, Drake, Wollner, Powderly, Picus, Sharfman, Stankevich, Pons, Salay, McMiller, Gilson, Wang, Selby, Taube, Anders, Chen, Korman, Pardoll, Lowy, Topalian, 2010).

The R code used to produce the figures of this manuscript, as well as to run the computations and stochastic simulations, is publicly available under https://doi.org/10.6084/m9.figshare.11536824.v1.

Materials and methods

To analyze the algebraic properties of a system of equations involving cancer-immune interactions, we used the program Mathematica (Wolfram Research, 2011). To find equilibrium points in situations where this was not algebraically possible, we used the rootSolve package in R (Soetaert, Herman, 2008, Soetaert, Soetaert, Petzoldt, Setzer, 2010).

Since all ordinary differential equations (ODEs) here described are deterministic, the time course of the decline of cancer cell numbers will always follow the same continuous trajectory given identical initial conditions. However, when small cancer cell numbers are reached, the temporal order at which the discrete events occur that underpin the dynamics will become important. Such events include the replenishment of immune cells and cancer cell deaths. Thus, at small cell numbers, accounting for the stochasticity of these events will add realism to the simulation, and help decide when eradication has effectively been achieved. To this end, we employed the Gillespie algorithm, where the interactions between cell types are explicitly simulated. Stochastic simulations of all ODEs were run in the R language for statistical computing (Team, 2012) by using the Gillespie algorithm (Gillespie, 1977) with tau leaping in the adaptivetau package (Johnson, 2011). If not stated otherwise, simulations were run with the set of parameter values given in Table 1. For alternative strategies to account for the stochasticity of CISI at the temporal mesoscale see (d’Onofrio, 2010).

Table 1

Standard parameter values for the base model.

Parameter	Description	Value
a	maximum replication rate of cancer cells (de Pillis et al., 2005)	0.514 day^−1
b	inverse carrying capacity of tumor (de Pillis et al., 2005)	1.02·10^−9cells^−1
k	killing efficacy of immune cells (Ganusov, Goonetilleke, Liu, Ferrari, Shaw, McMichael, Borrow, Korber, Perelson, 2011, Wick, Yang, Corey, Self, 2005)	10^−4−10cells^−1day^−1
σ	replenishment rate of immune cells (Althaus and De Boer, 2008)	10 cells·day^−1
d	immune cell death rate (Ogg, Kostense, Klein, Jurriaans, Hamann, McMichael, Miedema, 1999, Casazza, Betts, Picker, Koup, 2001)	1·10^−2day^−1
m	maximum immune cell proliferation rate	−10^−6day^−1

To model treatment, two possible procedures were considered. First, an increase in the ability of immune cells to detect and eliminate cancer cells — their killing rate or efficacy. Second, adoptive immune cell transfer, corresponding to the injection of immune cells into the system (Rosenberg, 1991, Rosenberg, Yang, Restifo, 2004). Both of these mechanisms enhance the suppression of cancer growth by the immune system and can be applied in concert as combination immunotherapy. The time at which the killing rate is first enhanced, the treatment time τ, can differ from the time at which cancer-specific immune cells are first injected into the system τ. Also, the time period during which each of these treatment approaches are administered can vary, with Δτ the treatment period for killing efficacy enhancement, and Δτ the period for immune cell transfer.

We assume that treatment always consists of the administration of either immunoactivating compounds or immune cells into the host system, and we denote the amount of compound delivered as the administered dose. In the increased killing efficacy approach, we assume that the alteration induced by the administration of the compound is permanent, which is reflected in a change of the killing efficacy parameter of NK or CTLs, Δc, or Δk, respectively. The change occurs gradually over the time course of the treatment. For example, an initial CTL killing efficacy k before treatment initiation will by increased by Δk/Δτ every day, leading to a final efficacy of .

In the immune cell transfer approach the change in immune cell numbers is not permanent. Prior to the transfer, immune cells are assumed to be rendered ineffective at a predefined rate (for example by the shedding of NKG2D ligands such as MIC-A, MIC-B (Mellman et al., 2011) or alternatively, by the upregulation the ligands PD-L1 or PD-L2 (Kooi, Zhang, Patenia, Edwards, Platsoucas, Freedman, 1996, Hamanishi, Mandai, Iwasaki, Okazaki, Tanaka, Yamaguchi, Higuchi, Yagi, Takakura, Minato, Honjo, Fujii, 2007)). Once transferred into the system, immune cell numbers will be affected by already present cancer cells. Thus, immune cell numbers will change depending on the state of the cancer, because the cancer exerts a suppressive effect on immune cell proliferation. Conversely, cancer cell numbers will vary due to immune cell killing. Analogously to the dosage of killing efficacy increasing compounds, effector cells are administered at daily doses of ΔE/Δτ cells, until the full dose of ΔE has been dispensed.

Results

Base-model of cancer-immune system interactions

We develop a base model of cancer-immune cell interaction. The model follows a large body of theory that uses two-equation deterministic ODEs to describe the interaction between cancer tumor cells and immune system cells (Kuznetsov, Makalkin, Taylor, Perelson, 1994, Kuznetsov, Knott, 2001, d’Onofrio, 2005, d’Onofrio, 2008, Eftimie, Dushoff, Bridle, Bramson, Earn, 2011, Eftimie, Bramson, Earn, 2010, Lopez, Seoane, Sanjuan, 2014). This model aims to replicate some basic features of cancer dynamics with a minimum of added complexity. With such an approach, qualitative insights about the behavior of the system can be obtained by relatively simple mathematical analysis. To this end, we make four fundamental assumptions. First, we assume the existence of immune cells, which are able to detect and kill tumor cells (Parish, 2003, Smyth, Godfrey, Trapani, 2001). These cells may eliminate microscopic tumors before they grow to endanger the organism; a process termed immunosurveillance (Burnet, 1957, Burnet, 1967). Second, these immune cells comprise the action of all cells that control tumor growth by antigen recognition and subsequent elimination (Eftimie et al., 2010), including natural killer (NK) cells (Khar, 1997) and CTLs (Boon and van der Bruggen, 1996) and are termed effector cells. The process by which the cancer cells are neutralized is called lysis (Quesnel, 2008). A background level of effector cells is present at all times (de Pillis and Radunskaya, 2003). Third, we assume that tumor growth is well described by a logistic growth in the absence of immune cells (Eftimie et al., 2010). Fourth, the interactions between tumor and cancer cells are governed by mass-action kinetics (Kuznetsov et al., 1994).

The third assumption of logistic growth warrants special discussion. The dynamics of tumor growth remain a debated issue in the literature. Benzekry et al. compared different theoretical growth dynamics in lung and breast tumor data of mice (Benzekry et al., 2014). They concluded that Gompertzian growth typically best predicts data. Gompertzian growth dynamics are motivated by the observation that tumor growth prior to detection appears to be faster than after detection (Eftimie et al., 2010). This suggests that the initial unbounded growth may be limited by the exhaustion of growth resources or cancer cell’s mutual growth impairment. Gompertzian growth dynamics, developed on the basis of birth and death processes, account for this behavior, and yield a sigmoidal type of growth curve for cancer cell numbers. However, the Gompertz model suffers from serious drawbacks, while other models did almost equally as well as the Gompertzian in predicting data (Benzekry et al., 2014). Because an upper bound of the proliferation rate is imposed by cell division time, the Gompertz model cannot adequately describe the dynamics of very small tumors (d’Onofrio, 2008, Eftimie, Bramson, Earn, 2010). Furthermore, theoretical analysis reveals that Gompertzian growth is at odds with the immuno-surveillance hypothesis, because the immune response is unable to eradicate cancers that grow in a Gompertzian fashion (Eftimie, Bramson, Earn, 2010, d’Onofrio, 2005). Thus, given these assessments, we chose a growth model that retains a sigmoidal cancer growth curve shape, namely logistic growth. With this, we preserve the notion of an exhaustion of growth resources. We note however, that alternative growth models may also successfully capture tumor growth patterns.

Our base model differs from most models in the literature in that it combines proliferative and suppressive effects of cancer cells on immune cells in one single term that describes its net effect. In this way, we can analyze how the systems behaves depending on the net effect of these opposing forces.

The equations for the base model are:

Tumor or cancer cells T, grow at a maximal rate a in a logistic fashion. The population density of the cancer cells is regulated by the coefficient b, which acts as an inverse carrying capacity. The cancer cells are detected and killed by effector cells E at a net rate k (Kuznetsov et al., 1994). We assume a constant supply of effector cells at a rate σ (Althaus, De Boer, 2008, De Boer, 2007), and a death rate d per effector cell (Wodarz, 2007). Effector cells can either be stimulated to proliferate or be impaired in their growth at a rate m, the net growth increase or decrease due to the presence of cancer cells. In other words, m can attain positive as well as negative values.

We proceed to analyze the possible equilibria of this system. We start by observing that is always a fixed point of the ODEs above. Two further solutions for T* can be represented by the quadratic formula (see Appendix A):where

A closer inspection of the properties of the dynamics of (1), (2) reveals that all biologically relevant cases, namely those in which are consistent with m < 0 (see Appendix A). This corresponds to a net immune cell proliferation suppression by cancer, which can arise by various mechanisms (Conway, Perelson, 2014, Keir, Butte, Freeman, Sharpe, 2008, Mellor, Munn, 2004, Aggarwal, Pittenger, 2005, Munn, Mellor, 2004, Marigo, Dolcetti, Serafini, Zanovello, Bronte, 2008, Kooi, Zhang, Patenia, Edwards, Platsoucas, Freedman, 1996, Hamanishi, Mandai, Iwasaki, Okazaki, Tanaka, Yamaguchi, Higuchi, Yagi, Takakura, Minato, Honjo, Fujii, 2007). m < 0 is also where a bistability pattern in the steady states of (1), (2) emerges. The alternative m > 0, leads to scenarios that are at odds with well established concepts of cancer modeling, and produce incomplete dynamics (see Appendix A). In particular, they conflict with the well-established notion of an immunological barrier (Kuznetsov et al., 1994); the idea that tumors have to grow above a critical threshold to reach a large size close to carrying capacity (Eftimie et al., 2010). Temporary changes in the activity of the immune system can lead to fluctuations in tumor size which place its size above the barrier, which then gives rise to cancer.

There are three possible cases of sign arrangements of the roots of (3) under m < 0: i) both negative, ii) one positive and one negative, iii) both positive. Scenario i) has as only biologically plausible solution. Only the cancer-free state exists. Case ii) signifies a single attractive equilibrium at a non-zero tumor size. The emergence of a cancerous cell suffices to ignite a replicative process that induces the establishment of a tumor close to carrying capacity 1/b. Thus iii), where is the only case which admits stable equilibria compatible with an existing immunological barrier.

For to be satisfied, and while assuming that a > 0 and b > 0 for biological reasons, we obtain the following conditions (see Appendix A):

These results reveal a bistability pattern that is mediated by the killing efficacy k. Fig. 2 shows how the increase in the parameter k leads to a bifurcation in the stable states of T and to bistability for k. At k < k the system will reside in the aforementioned case ii). By increasing k above k but below k, the system will enter case iii), and move to i) as k > k. In line with our expectations, values of k below the threshold value k represent a similar situation as would be expected in the absence of immune cells: unchecked cancer growth. If k is gradually increased above k, the tumor cells would have to begin replicating at increasingly large initial sizes in order to avoid being absorbed by the attractor at that is, to be suppressed by the immune system. This is where the bifurcation appears, and now three equilibria, of which two are stable, dominate the dynamics. The parameter range spanned between k and k, is highly dependent on the ratio of b and m, as well as σ. Lastly, large values of k above k entail that even few effector cells are able to clear the tumor, and even large tumors are eliminated with certainty.

Fig. 2

Bistability of immune control of tumor growth. Increasing values of k lead to the emergence of a bifurcation in equilibrium tumor cell numbers. The emergence of the bifurcations is indicated by black, vertical dashed lines. The upper blue line shows the largest stable steady states and the lower dashed blue line shows the smallest of the steady states, while red lines show unstable steady states for a given k value. If bT is one, the cancer has reached carrying capacity. At no cancer cells exist. Thus, above a threshold value k, the cancer is cleared. The parameter values are (Althaus, De Boer, 2008, De Boer, 2007), (Ogg, Kostense, Klein, Jurriaans, Hamann, McMichael, Miedema, 1999, Casazza, Betts, Picker, Koup, 2001), (de Pillis et al., 2005), (de Pillis et al., 2005). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

At intermediate k (case iii) the base level of effector cells crucially affects the dynamics. In the absence of cancer, the equilibrium level of effector cells is at σ/d. The appearance of cancer cells will induce a killing process mediated by k. High enough effector cell levels will suppress tumor growth. In reality, the tumor might fluctuate above some threshold within which effector cells control it, and the suppression of further T cell proliferation (m < 0) by the cancer will dominate over the killing. Thus, the conditions (5), (6), (7), (8) mark the region of killing efficacy values that constitute an immunological barrier to the cancer. Increasing k values imply that this barrier is heightened: existing effector cells improve their capacity to completely eliminate the cancer.

This analysis may serve as a model to think about how to efficiently combine immunotherapy approaches. These results suggest that one mechanism to generate biologically plausible bistability is consistent with situations in which the cancer’s immunosuppressive effects outweigh its immunoproliferative effects (m < 0). The existence of a bifurcation implies that increasing killing efficacies will not simply gradually diminish cancer cell numbers. Instead, an increase above a critical value k can terminally clear the cancer. Equally, if k is located in the range k < k < k, the adoptive transfer of effector cells may work by perturbing cancer cell numbers below the unstable equilibrium point after which it will be absorbed into and be cleared.

Model variations

The here presented base-model, despite being strongly simplified, can help us intuitively understand more complex models of cancer-immune system interactions. In the following, we demonstrate that the basic bistability phenomenon can be replicated in two related, but more complex variations of the base-model. When possible, we also analyze whether multistability can arise, since it is relevant to cancer dormancy, and also because it might mitigate the effect of treatment.

Incorporating saturation effect of tumor cells on immune response

A natural way to extend the base-model is to include more biologically plausible assumptions about the behavior of effector cells. Here, we retain the basic assumptions going into the base model, but refine the way that effector cells are reacting to the presence of tumor cells. In particular, we follow the approach by Conway et al., which allows for a very general set of behaviors of effector cells to arise and also provides estimates for the parameters (Conway and Perelson, 2014). This study is set in the context of HIV, where a great body of theory has been devoted to the particulars of CTL behavior (Wodarz, 2007). The equations are:

Here, the interactions governing the rate of change in tumor cells have remained intact. Effector cell growth can now be stimulated at a maximum rate b (Davenport, Ribeiro, Chao, Perelson, 2004, Conway, Perelson, 2014). The proliferation rate saturates with the number of tumor cells T, and is half-maximal at the Michaelis-Menten constant κ (Conway and Perelson, 2014). In contrast, effector cells can be exhausted by contact with tumor cells, and die from its consequences at a maximal rate d (Johnson, Kochin, McAfee, Stromnes, Regoes, Ahmed, Blattman, Antia, 2011, Conway, Perelson, 2014). As with proliferation, the effects of exhaustion – mediated by κ – saturate. For typical values of these parameters see Table S1 in Appendix B. Models with a saturation term in have been analyzed before (Kirschner and Panetta, 1998), and have been thoroughly discussed in for example (Talkington et al., 2018).

As with the base model, is always a fixed point. The rest of the fixed points are determined by the roots of a cubic equation (see Appendix B). For the system of equations to generate bistability, that is, to give rise to exactly two positive equilibria in T, the following set of conditions needs to be satisfied (see Appendix B): where

Analogously to the base model, a closer inspection of this result reveals that biologically plausible dynamics are consistent with (see Appendix B.1). Interestingly, this expression is independent of κ or κ. Since a, b > 0, this implies that which is analogous to the situation where m < 0 in the base model. The treatment rationale identified in the base model may therefore also be applicable for the extended base model with saturation (10), (11).

In Conway et al’s work (Conway and Perelson, 2014), this condition is satisfied by b < d, which is also functionally equivalent to m < 0 in our base-model: the effector cell population decreases due to cancer-mediated exhaustion. Importantly, the parameter choice in Conway et al. is also consistent with the emergence of bistable and multistable equilibira.

Besides bistable patterns, the system (10), (11) can also generate a pattern of multistability. The conditions to obtain four equilibria in T for (10), (11) reads:where A, B, C and D are as in (15). Again, a biologically plausible arrangement of equilibrium points is in agreement with A < 0. Figure A1 shows that bistability becomes common for d > 1 and with the elimination of cancer ensuing after a critical kthreshold is surpassed.

Fig. A1

Number of possible non-negative steady states ofTof the model(10), (11). Shown is the region of parameter space spanned by killing efficacy k and effector cell exhaustion d divided into regions of different non-negative fixed points of the model (10), (11). The number of roots were identified numerically by means of the rootSolve package in R (Soetaert, Herman, 2008, Soetaert, Soetaert, Petzoldt, Setzer, 2010). Numbers indicate the number of non-negative fixed points in each parameter region. The region that supports multistability is colored in violet. The region that admits bistability (3 fixed points) is adjacent to the multistability region. The other used paramaters are as in Table 2.

The existence of multiple stable equilibria may be interpreted as cancer dormancy (see Discussion). Cancer dormancy is the phenomenon of a period of non-growth of tumors. Often, this occurs in small, nearly undetectable tumors residing within body tissues (Wilkie, 2013). Fig. A1 shows that multistability is possible in (10), (11). The existence of multistability in two-equation models has been predicted and shown in other work (d’Onofrio, 2008, de Vladar, González, 2004, Kuznetsov, Makalkin, Taylor, Perelson, 1994). Here, we give precise analytical conditions for its emergence under (10), (11). When in a multistable regime, an increase in killing efficacy k might not directly lead to cancer eradication if treatment is started when the cancer is near carrying capacity 1/b. Instead, a new microscopic steady state (MISS, (d’Onofrio, 2008)) might be attained before a further increase in k leads to cancer clearance.

Incorporating natural killer (NK) cells and tumor-specific CTL response

In a next step, we incorporated a further level of complexity by distinguishing between two types of effector cells: natural killer or NK cells, and cytotoxic T lymphocytes or CTLs. Models that account of the different roles between NK and CTL can be highly complex (de Pillis et al., 2005). To better understand where possible bistabilities originate from, we restrict ourselves to an extension of the base model with saturation, following an approach inspired by de Pillis et al. (2005) and Conway and Perelson (2014). Effector cells are now split into NK cells N and CTLs E:

Here, cancer cells T are killed at rates c by NK cells, and at rates k by CTLs. The dynamics of the NK is now analogous to the dynamics of immune cells in (10), (11). We assume a constant supply of NK cells σ stemming from the host’s hematopoesis. NK cells die naturally at a rate μ. The maximum NK proliferation rate induced by the presence of cancer cells is b, and the saturation coefficient is κ. Again, exhaustion occurs at a maximum rate d and the saturation in T is half-maximal at κ. For CTLs, following (de Pillis et al., 2005), we assume that there is no constant supply of cells. Instead, CTL growth is stimulated by the NK-cancer cell interactions at a rate ω. This modeling approach ensures that CTLs are activated only after the emergence of the NK immune response (de Pillis et al., 2005). Several biological mechanisms appear to exist by which NKs can stimulate CTL growth (Pallmer and Oxenius, 2016). The work of Fan et al. suggests that already activated NK cells can facilitate the priming of CTLs by means of IFN-γ (Fan et al., 2006). The proliferation and exhaustion terms are as in the extended base model.

The model (21), (22), (23) can generate substantially more complex behavior than the two previously analyzed. As in the previous models (1), (2) and (10), (11), is always a fixed point. However, the NK-CTL model may allow for up to five other fixed points (see Appendix C). This is because finding the steady states of (21), (22), (23) can be reduced to finding the roots of the rate of change of T, which is a polynomial of sixth order. The five fixed points besides are the roots of a fifth order polynomial, for which no general solutions exist. Thus, we cannot draw similar conclusions for the existence of real fixed points as in the previous, relying on analogous conditions on discriminants Δ > 0 or Δ > 0. However, similarly to the base model (1), (2), a similar condition to m < 0 is compatible with biologically plausible arrangements of fixed point’s stabilities. The expression analogous to m < 0 is (see Appendix C):

The analogy to m < 0 arises from k, σ, ω > 0, valid bounds in most biological contexts. If tumor cells are able to exhaust NKs and CTLs, that is, if both and then the system can display biologically reasonable behavior. Note that unlike the previous models, this condition can be satisfied by other means as well, such as increasing k. Again, the condition is independent of the four saturation coefficients κ.

Due to the analytical unfeasability of the model (21), (22), (23), we resorted to numerical methods to prove the existence of basic bistability patterns (see Appendix C) (Soetaert, Herman, 2008, Soetaert). We found that the system is able to display bistability similar to that found in the extended base model with saturation (see Fig. A2).

Fig. A2

Bistability of cancer cell fixed points induced by immune control of tumor growth in model(21), (22), (23). Panels A) to C) show the function dT/dt(T) (blue lines) for different parameter values of c. c is varied from (A) to (B) 100 (C). The blue points indicate the numerically found roots of dT/dt(T). Bistability emerges at . The parameter values were as specified in Table 3.

Bistability-based Strategies of Cancer Immunotherapy

The existence of bistability patterns in simple non-spatial cancer models as well as its variations, can be informative to the assessment of immunotherapeutic options and of their efficacy. Taking the base model (1), (2) as a foundation, three intervention approaches seem apparent. First, the elimination of the exhaustive effects of cancer on the immune cells (m < 0 → m > 0). Second, the increase of the killing efficacy of effector cells above some threshold (k < k → k > k). Third, the administration of effector cells (). In terms of the dynamics, this represents pushing the state of the system into the attraction basin of . Combinations of these therapy approaches have previously been explored in simulations (Kirschner and Panetta, 1998).

In current immunotherapy, the two main available tools for cancer cell reduction correspond to the second (antibody therapy) and third (adoptive T cell transfer) options (Dougan, Dranoff, 2009, Mellman, Coukos, Dranoff, 2011). In antibody therapy, an increase in killing efficacy is attained by disrupting cancer cells' mechanisms for impairing T cell action. This impairment occurs by the acquisition of mutations in cancer cells that, for example, lead to the expression of the PD-L1 and PD-L2 ligands (Kooi, Zhang, Patenia, Edwards, Platsoucas, Freedman, 1996, Hamanishi, Mandai, Iwasaki, Okazaki, Tanaka, Yamaguchi, Higuchi, Yagi, Takakura, Minato, Honjo, Fujii, 2007). These ligands are known to bind to the PD-1 receptors on T cell surfaces, thereby downregulating the activation of the T cells. Thus, a direct effect of PD-1 immune checkpoint blockade is to halt cancer-induced T cell anergy and exhaustion (Mellman and Steinman, 2001). In this work we assume that ultimately, monoclonal antibodies binding to the ligands effectively increase k by interrupting this cancer escape mechanism. Possible increases in T cell recruiting and proliferation are neglected. In adoptive T cell transfer (Rosenberg, 1991, Rosenberg, Yang, Restifo, 2004), T cells are pre-programmed to kill host cells that carry particular biochemical signatures, for example certain peptides on their surface. The signatures are chosen such that they match characteristic features of cancer cells. Subsequently, these specific T cells are grown and injected into the blood stream of the patient (Mellman et al., 2011).

These two novel treatment methods can be combined to take advantage of the bistability phenomenon in cancer. We used the base model (1), (2) to investigate how a combination of both approaches could be used to clear the tumor, while accounting for the stochastic effects arising from singular cell-to-cell interactions. Starting from an already established tumor, increasing the efficacy of killing of effector cells will not by itself necessarily lead to the elimination of the tumor, unless very high levels of killing efficacy can be attained. Instead, increasing the killing efficacy by two orders of magnitudes will lead the system to equilibrate at tumor cell numbers lower than the carrying capacity (see Fig. 2). If the treatment with monoclonal antibodies has been sufficiently effective, it will have shifted the system into a regime with two stable equilibria, out of which one is the cancer-free state. If now the system is perturbed further with adoptive T cell transfer into the attraction basin of the cancer-free equilibrium, the waning of the effects of the first treatment will not lead to the reemergence of the cancer. Thus, the generation of a temporary bistability in the system can be exploited to perturb it into a cancer-free state.

To model the combined effects of killing efficacy increases by PD-1 specific monoclonal antibody and adoptive T cell transfer treatments, we used stochastic version of the base cancer immune interaction model (1), (2) (see Materials and Methods). Fig. 3 shows the time courses of the dynamics with and without treatment. Without treatment, a cancer that has surpassed the immunological barrier will grow unrestricted up to levels very close to carrying capacity (Fig. 3A). The combined administration of effector cells and killing efficacy-increasing compounds at first only gradually reduces cancer numbers (Fig. 3B). Daily administered effector cells ΔE/Δτ (Materials and Methods) can only temporarily affect the dynamics before they are rapidly suppressed and exhausted by the cancer (m < 0). When the state of the system is pushed into the attraction basin of the cancer-free state, cancer numbers rapidly go to zero. Further injections of effector cells become unnecessary, and effectors build up.

Fig. 3

Simulation runs of the base model without A) and with B) immunotherapeutic treatment. A) Cancer replication begins at and the natural equilibrium of the effectors is at . The cancer grows to carrying capacity in a time frame of around 45 days. B) Combined immunotherapeutic treatment is initiated at 50 days after the cancer has begun to grow. It lasts for 21 days in antibody therapy, and 28 days in adoptive cell transfer. Killing efficacies are increased to while a total of cells are injected in a gradual fashion. Immune cells are eliminated rapidly by the powerful immune exhaustion effects exerted by the cancer cells. Parameter values are as in Table 1. In particular, .

We then investigated whether an increase in k and the administration of effector cells E work together in a synergistic or antagonistic fashion to remove the tumor. Fig. 4 shows the outcome of combination immunotherapy, initiated simultaneously for antibody and T cell injections. Each immunotherapeutic approach may clear the cancer on its own. A marked frontier between cancer presence and clearance emerges. A lowering of killing efficacies along this frontier will lead to insufficient pressure to clear the cancer, but can be compensated by an increase in adoptive transfer doses. The linear shape of the frontier indicates that the two approaches do not mutually impair each others’ function.

Fig. 4

Cancer cell numbers after combination immunotherapy after 120 days. Immunotherapy begins at day 50, after the cancer has grown to full size, and lasts for a single day. Combination immunotherapy is implemented by increasing the killing efficacy of the effector cells and adoptive immune transfer of effector cells. Darker colors indicate high cancer cell numbers, implying that the cancer persists, while lighter colors indicate low cancer cell numbers. Parameter values are as in Table 1.

Discussion

In this study, we have shown that the base model can only reproduce biologically plausible behavior if the suppressive effects exerted by cancer cells on immune cells dominate their proliferative effects. Under these circumstances, the base model displays a conspicuous pattern of bistability: The cancer-immune interaction dynamics gives rise to two distinct, stable states (a cancer-free, and a full-grown tumor state). Under bistability, the modification of the killing efficacy can lead to a bifurcation in cancer cell numbers, where the system may abruptly be tipped into a new, cancer-free state. Furthermore, in situations where exhaustion prevails over proliferation in immune cells, all analyzed models can produce bistability patterns that are biologically plausible. If this condition is not satisfied, the base model cannot produce biologically plausible behavior across a wide range of k values. We also formulated more complex extensions of the base model, which can generate multistability, a dynamic behavior that can be interpreted as stable microscopic cancers (d’Onofrio, 2005), or cancer dormancy (Kuznetsov, Makalkin, Taylor, Perelson, 1994, Wilkie, 2013). We gave the exact conditions under which multistability might arise in one model extension inspired by Conway and Perelson (2014). We also examined how bistability may be used for effective combination immunotherapy. We tested a combination of two different immunotherapeutic approaches in stochastic simulations of the base model. We found that the combination of treatment interventions is able to clear the cancer, and that the different treatment approaches do not impair one another.

How both treatment approaches investigated here would work in isolation has been studied in other work (d’Onofrio, 2008, Kuznetsov, Knott, 2001, Eftimie, Bramson, Earn, 2010). However, it is less apparent that they may be employed to work in concert without mutual impairment. This is important when considering side effects of these therapies. For example, the administration of PD-1 antibodies in mice have resulted in lung inflammation and cardiomyopathy (Nishimura, Okazaki, Tanaka, Nakatani, Hara, Matsumori, Sasayama, Mizoguchi, Hiai, Minato, Honjo, 2001, Nishimura, Nose, Hiai, Minato, Honjo, 1999, Mellman, Coukos, Dranoff, 2011). Thus, combination immunotherapy may help minimize the risks associated with standalone approaches.

Another advantage of this study lies in that it specifically shows how multistability may arise from standard assumptions about cancer-immune system interactions (saturation terms), and in that it gives precise conditions for its emergence. Although the first extended model with saturation (10), (11) is not intended to admit multiple steady states, these arise as a consequence of the basic assumptions about interactions.

Intermediate-sized cancers in multi-stable regimes may be interpreted as cancer dormancy (d’Onofrio, 2008, Kuznetsov, Makalkin, Taylor, Perelson, 1994, de Vladar, González, 2004), but the model (10), (11) does not explicitly explain how they may escape immune control. For immune escape to arise, additional processes must be assumed, such as stochastic perturbations or immunoediting (Wilkie and Hahnfeldt, 2013a). A discussion of how tumors might rise to large numbers in models very similar to the base model with saturation has been given in Kuznetsov et al. (1994) and by Wilkie and Hahnfeldt (2013b). In the Kuznetsov et al. model, a separatrix between the two main attraction basins passes close by the trivial, cancer-free steady state. When a cancer arises and starts to replicate at low numbers, it should follow a trajectory into a dormant, stable steady state. However, stochastic fluctuations can push the system’s state into the other attraction basin. In the Wilkie and Hahnfeldt model, the saddle point is the dormant state itself. Another of the main mechanisms hypothesized to drive the transition from dormancy to large tumors is immunoediting (Wilkie, 2013): The prolonged growth suppression of the tumor by the immune response leads to the selection of cancer mutations that escape immune pressure, effectively reducing the immune killing efficacy k Wilkie and Hahnfeldt (2013a). The model (10), (11) can also offer an intuitive explanation for this process, whereby a smaller, undetectable and stable equilibrium of cancer cells is maintained by a relatively weak immune response. The further decrease of the immune response efficacy by means of immune escape processes leads to the establishment of a full grown tumor. This is exemplified by the fact that decreasing k in Fig. A1 pushes the system into a region of parameter space where there exists a stable steady state for the tumor at carrying capacity, that is the full grown cancer state.

Most other models of cancer-immune interaction so far have attributed the phenomenon of dormant states to the existence of an additional compartment: quiescent cancer cells (Page, Uhr, 2005, Wilkie, 2013). These cells are assumed to replicate at a slower rate than normal cancer cells, and can revert back to a fast growing state by means of phenotypic switching or by acquiring further mutations (Wilkie, 2013). In line with other work, the model (10), (11) explains the existence of dormancy by a specific balance of cancer cell growth and killing attained in cancer-immune interactions, without relying on any additional compartments or assumptions. A notable example of how dormancy can emerge from cancer-immune interactions alone is given in Kuznetsov et al. (1994). A third mechanism for the emergence of dormancy has been given by Wilkie and Hahnfeldt (2013b). In this mechanism, dormant states are represented by saddle nodes traversed by a separatrix demarcating the adjacent attractor regions of either growth progression or tumor clearance.

This point emphasizes a last advantage of non-spatial ODE models: Understanding cancer growth requires an appropriate description of cancer-immune system interactions at multiple scales (Altrock et al., 2015). These scales range from cancer microenvironments to large numbers of already systemic cancers. ODE models offer useful tools to combine the behavior of both into a single framework (Eftimie et al., 2010), accounting for the frequency-dependent growth at early stages as well as the dominant immunosuppressive effects achieved by cancer when approaching carrying capacity levels.

A major caveat of the base model is that it cannot elicit an immune response. The elicitation of an immune response by a cancer, with a subsequent rise in effector cells, is an important aspect of cancer-immune system interactions. To achieve this, the base model would have to include a proliferation term for effectors that takes on different values than the suppression in the T, E-plane. The base model is thus more useful to study the aspects of how immunotherapy can be deployed to return the CISI system below an immunological barrier by external perturbation.

While simple modeling frameworks offer greater possibilities for an in-depth understanding, this approach also has its limitations. For instance, we have largely neglected stochastic attributes of cancer-immune system interactions in our mathematical analysis. These may mostly arise from the discreteness of cell-to-cell interactions, and are well captured by simulating the models by a Gillespie algorithm (Gillespie, 1977). We have addressed this shortcoming by adopting a stochastic simulation framework to implement immunotherapy. Other forms of stochasticity —for example the random accrual of malignant mutations in cancer cells— are also not explictly modeled. Instead, they are assumed to be captured by model parameter values. Environmentally based fluctuations (Eftimie, Bramson, Earn, 2010, Bose, Trimper, 2009), or changes in the exerted immune pressure due to, for example, disease (Mina et al., 2015), are also neglected.

The models here analyzed do also not account for spatial structure (discussed in more detail in Araujo, McElwain, 2006, Roose, Chapman, Maini, 2007, Chaplain, 2008). Spatial structure may change the way that effector cell killing affects cancer growth, as well as how the presence of cancer cells may mediate immune cell proliferation. In particular, we did not explore the fractional cell kill laws as introduced by de Pillis et al. (2005), hypothesized to account for some of the geometrical features of tumors (Lopez et al., 2014). In this approach, the total killing exerted by effectors K(E, T) is governed by the de Pillis-Radunskaya-Wiseman (PRW) law Lopez et al. (2014), where and . Structurally, how (de Pillis et al., 2005) implement the recruitment of NK cells differs only slightly from our implementation in (21). However, how CTLs are recruited differs in structure from the simpler terms analyzed here in our models. With λ < 1 (obtained from model fits to mouse data (Lopez et al., 2014)), the behavior of the recruiting is qualitatively similar to that studied here: the recruiting of CTLs would then saturate with increasing cancer cell numbers T, but continue growing with increasing effector cell numbers E.

We thus assume that while the PRW law introduces an advantageous new concept in the modeling in tumor-immune system interactions, our deviating from it will not yield marked qualitative differences. Instead of attempting to capture tumor geometry behavior, the models analyzed here are rooted in the tradition of virus dynamics —especially HIV— which assumes well-mixed cell types (Nowak, May, 2000, Conway, Perelson, 2014). Thus, the PRW law seems to mainly address problems arising from tumor geometry, while this study focuses mainly on systemic cancer types —cancer types that do not manifest in single tumors only (Eftimie et al., 2010).

We have also not included the action of cytokines in our analysis, which are typically accounted for with a separate, additional equation (Arciero, Jackson, Kirschner, 2004, Kirschner, Panetta, 1998, Eftimie, Bramson, Earn, 2010). We have therefore not been able to assess the effectiveness of cytokine-based immunotherapy approaches in combination with the ones studied here. Models with cytokines display features like the persistence of large tumors, tumor dormancy, and tumor clearance upon immunotherapeutic treatment, as well as oscillations between these states (Eftimie, Bramson, Earn, 2010, Kirschner, Panetta, 1998). Including cytokines into a more comprehensive modeling framework would be an interesting topic for future work.

Our study has to be interpreted in the context of other models of cancer-immune system interactions. The most comprehensive mathematical analysis of two-equation models has been put forth by (d’Onofrio, 2005, d’Onofrio, 2008). d’Onofrio analyzed a generalized mathematical model in two variables —x denotes cancer and y effector cell densities—, deriving some general results on the existence of steady states and cancer eradication given some broad mathematical conditions on the interactions between these two cell types. Solutions were provided for the generalized model, but except for the rate of adoptive transfers θ(t) (where t is time), no specific dependence was given on how steady states change with parameter value modifications. Our base model corresponds to a special case of his general model (1), (2) in d’Onofrio (2008)), with and . As in the models investigated in this study, d’Onofrio has observed that his generalized model admits a cancer-free state, and also predicted that it may attain multiple stable equilibria, which he interpreted as microscopic steady states (MISS) and which we interpret in the context of cancer dormancy. Our own results thus confirm some of d’Onofrios, but go further to explore how specific cancer-immune system interaction models are concretely affected by changing dynamical properties that may be tailored for immunotherapy. In particular, we wanted to explore some of the properties of CISI models that underpin the mechanisms that may give rise to bistability patterns (the aforementioned dominance of immunosuppressive effects). In our models, we were therefore more interested in explicit analytical results, which would allow us to study how bistability patterns depend on effector cell killing efficacy. We also extended this approach to include how adoptive transfer might function under conditions with stochasticity.

A similarly comprehensive analysis of how CISI models may give rise to successful adoptive immunotherapy treatments has recently been put forth by Talkington et al. (2018). Similarly to our own conclusions, Talkington et al. identify bistability as a major prerequisite for successful adoptive immunotherapy. Their approach is also to review a series of models of increasing complexity, whereby complexity is understood to represent incorporations of additional aspects of the immune system, such as helper cells, interleukin and naïve T cells into a base model. The base model is Kuznetsov et al.’s early model from 1994 (Kuznetsov et al., 1994), which, akin to our base model, assumes only two compartments: tumor and effector cells. The Kuznetsov model allows for bistability, with a stable state of the tumor close to cancer eradication. For all other models, Talkington et al. show that when they can give rise to bistability, adoptive immunotherapy leads to successful outcomes in simulations.

Future work could address whether the mechanism for bistability emergence identified in this study, the dominance of immunosuppression by cancer over immune cell proliferation, is also the one that gives rise to bistability in the models examined by Talkington et al. To this end, it will be useful to embark on a more comprehensive analysis of how mathematical models as the ones put forth here, bring about some behavior of interest, such as bistability. One approach that could be taken in this direction follows the axiomatic modeling framework pioneered by Komarova et al. (2003) and (d’Onofrio, 2008, d’Onofrio, 2005).

In a departure from the work of de Pillis et al. (2005), we have concentrated on model features typically used in HIV modeling, borrowing in particular from the study of (Conway and Perelson, 2014). The reason for this choice is that ODE-based modeling has a long tradition in HIV and virus dynamics modeling (Wodarz, 2007, Nowak, May, 2000, Perelson, Ribeiro, 2013). A great wealth of data have helped to validate interaction terms of different cell types in mixtures, and particularly, how CTLs kill. In our view, these advantages can be fruitfully employed in cancer modeling as well. A more inter-disciplinary integration, in particular with respect to CTL behavior, will benefit both fields, and allow for the analysis of structural similarities between models that might be harnessed for immunotherapy design.

Our models show that biologically plausible cancer-immune system interactions may be utilized to induce cancer-free states. Increases in the killing efficacy of immune effector cells can destroy the bistability pattern inherent in those models, abruptly removing the basis for cancer growth.

Funding

This work was supported by the European Research Council Advanced Grant [grant number PBDR 268540]; the Swiss National Science Foundation [grant number P2EZP3_162257]; and SystemsX the Swiss Initiative for Systems Biology [grant number 51FSP0_163566]. F.F would like to acknowledge the NIH COBRE Program(grant no. 1P20GM130454) and the Bill & Melinda Gates Foundation (award no. OPP1217336) for generous financial support.

CRediT authorship contribution statement

Victor Garcia: Conceptualization, Investigation, Methodology, Formal analysis, Writing - original draft, Writing - review & editing. Sebastian Bonhoeffer: Conceptualization, Funding acquisition, Project administration, Supervision, Resources, Writing - original draft, Writing - review & editing. Feng Fu: Conceptualization, Investigation, Methodology, Formal analysis, Funding acquisition, Project administration, Supervision, Resources, Writing - original draft, Writing - review & editing.

Table 2

Parameter values for the extended base model with saturation (10), (11).

Parameter	Description	Value
a	maximum replication rate of unmutated tumor cells (de Pillis et al., 2005)	0.514 day^−1
b	inverse carrying capacity of tumor (de Pillis et al., 2005)	1.02·10^−9cells^−1
k	killing efficiency of effector cells or CTLs (Ganusov, Goonetilleke, Liu, Ferrari, Shaw, McMichael, Borrow, Korber, Perelson, 2011, Wick, Yang, Corey, Self, 2005, Elemans, Florins, Willems, Asquith, 2014, de Pillis, Radunskaya, Wiseman, 2005)	10^−5−10^−3cells^−1day^−1
σ	replenishment rate of effector cells or CTLs E (Althaus, De Boer, 2008, De Boer, 2007)	5 · 103cells·day^−1
d	effector or CTL loss rate (Conway and Perelson, 2014)	2 day^−1
de	effector CTL immune impairment rate	10^−3−10^3day^−1
be	maximum CTL proliferation rate (Davenport, Ribeiro, Chao, Perelson, 2004, Conway, Perelson, 2014)	1 day^−1
κe	saturation coefficient of effector immune growth	5 · 102 cells
κd	saturation coefficient of effector immune exhaustion	2.5 · 104 cells

Table 3

Parameter values for the extended model with saturation in NK and CTL (21), (22), (23).

Parameter	Description	Value
a	maximum replication rate of unmutated tumor cells (de Pillis et al., 2005)	0.514 day^−1
b	inverse carrying capacity of tumor (de Pillis et al., 2005)	1.02·10^−9cells^−1
c	killing efficacy natural killer cells	10^−7−10^0cells^−1day^−1
k	killing efficiency of effector cells or CTLs (Ganusov, Goonetilleke, Liu, Ferrari, Shaw, McMichael, Borrow, Korber, Perelson, 2011, Wick, Yang, Corey, Self, 2005, Elemans, Florins, Willems, Asquith, 2014)	10^−5−10cells^−1day^−1
σ	replenishment rate of effector cells or CTLs E (Althaus, De Boer, 2008, De Boer, 2007)	10−10^3cells·day^−1
bn	proliferation rate of natural killer cells N	1.5 cells^−1day^−1
dn	exhaustion rate of natural killer cells N	1 cells^−1day^−1
μ	death rate of natural killer cells N	1 day^−1
d	effector loss rate (Conway and Perelson, 2014)	2 day^−1
de	CTL immune impairment rate (Johnson, Kochin, McAfee, Stromnes, Regoes, Ahmed, Blattman, Antia, 2011, Conway, Perelson, 2014)	2 day^−1
be	maximum CTL proliferation rate (Davenport, Ribeiro, Chao, Perelson, 2004, Conway, Perelson, 2014)	1 day^−1
ω	growth stimulation rate of cytotoxic T lymphocytes E by natural killer cells (de Pillis et al., 2005)	1.1·10^−7cells^−1day^−1
κbe	saturation coefficient of CTL immune growth	5 · 102 cells
κde	saturation coefficient of CTL immune exhaustion	2.5 · 104 cells
κbn	saturation coefficient of NK immune growth	104 cells
κdn	saturation coefficient of NK immune exhaustion	103 cells