A Race between Tumor Immunoescape and Genome Maintenance Selects for Optimum Levels of (epi)genetic Instability

The human immune system functions to provide continuous body-wide surveillance to detect and eliminate foreign agents such as bacteria and viruses as well as the body's own cells that undergo malignant transformation. To counteract this surveillance, tumor cells evolve mechanisms to evade elimination by the immune system; this tumor immunoescape leads to continuous tumor expansion, albeit potentially with a different composition of the tumor cell population (“immunoediting”). Tumor immunoescape and immunoediting are products of an evolutionary process and are hence driven by mutation and selection. Higher mutation rates allow cells to more rapidly acquire new phenotypes that help evade the immune system, but also harbor the risk of an inability to maintain essential genome structure and functions, thereby leading to an error catastrophe. In this paper, we designed a novel mathematical framework, based upon the quasispecies model, to study the effects of tumor immunoediting and the evolution of (epi)genetic instability on the abundance of tumor and immune system cells. We found that there exists an optimum number of tumor variants and an optimum magnitude of mutation rates that maximize tumor progression despite an active immune response. Our findings provide insights into the dynamics of tumorigenesis during immune system attacks and help guide the choice of treatment strategies that best inhibit diverse tumor cell populations.


Introduction
In 1909, Paul Ehrlich was the first to propose the idea that the immune system scans for and eliminates nascent transformed cells in the human body [1]. This hypothesis received much interest from both immunologists and cancer researchers and led to experiments with tumors transplanted into mice; these studies suggested the existence of tumor-associated antigens and formed the basis of the idea of immune surveillance [2]. Since these landmark studies in the 1950s, the model of cancer immune surveillance has gained widespread acceptance, and the central role of immune effector cells, such as B, T, and natural killer (NK) cells, have been elucidated [3,4,5,6,7]. NK cells and CD8 z cd T cells were found to recognize and kill tumor cells through the interaction of specific cell surface receptors with tumor cell ligands [3,8,9,10,11,12]. Similarly, CD4 z and CD8 z ab T cells recognize MHC class II and class I molecules on tumor cells, respectively, and B cells produce antibodies against tumor antigens [3,6,13]. When the immune system fails to eliminate all tumor cells, then the malignant cell population continues to grow -a phenomenon termed ''tumor immunoescape''. The interaction with the immune system, however, may significantly decimate the abundance of tumor cells and select for those phenotypes with relative resistance against immune system attacks. The ''cancer immunoediting'' hypothesis then predicts that, while one outcome is complete eradication of a developing tumor, another is the generation of a sculpted tumor cell population that either displays reduced immunogenicity [4] or an increased ability to inhibit antitumor immune responses [6,14,15,16]. The latter capacity may be imparted via diverse mechanisms [17,18]: (i) tumor cells can lose their MHC class I molecules, enabling them to evade CTL attacks [19]; (ii) while the immunodominant epitope becomes the main target of immune responses, cells with other phenotypes may continue proliferating in the ''shadow'' of the dominant clone [20]; (iii) furthermore, tumor cell secretion of immunosuppressive cytokines such as TGF-b and IL-10 can reduce the efficiency of the immune response [21], and (iv) a modification of death signaling may prevent cells from undergoing apoptosis [22].
Tumor immunoescape is driven by the generation of tumor cell variants [17,23]. Frequent genetic and epigenetic alterations enable tumor cells to lose MHC class I molecules, produce immunosuppressive cytokines, and generate other phenotypes that are selected to escape immunosurveillance. Although cells with normal rates of accumulating such alterations may also manage to evade the immune response, this process is accelerated by the evolution of genomic instabilities [24,25]. Genomic instabilities are common in most cancer types [26], and two main categories have been identified: in the majority of tumors, chromosomal instability (CIN) leads to an increased rate of losing or gaining whole chromosomes or parts of chromosomes during cell division [24]; in a smaller fraction of cancers, a mismatch-repair deficiency leads to microsatellite instability (MIN) at the nucleotide level [27]. Similar to genomic instabilities, epigenetic instabilities were also recently found to contribute to tumorigenesis by modulating the production of oncogenic proteins [28].
An increased chance of accumulating (epi)genetic alterations during cell divisions enhances the rate of generating tumor cell variants that may evade the immune response; however, high rates of alterations may in turn lead to an error catastrophe in that a functioning genome cannot be sustained when error-prone replication produces excess damage [29]. The concept of an error catastrophe was first introduced to describe the behavior of RNA viruses [30], and numerous observations about the extinction of such viruses due to excess error have been reported [31,32,33]. These findings imply that a mutator phenotype does not serve as an unequivocal benefit for tumor cells, but also harbors a risk of extinction if the extent of variability in the population crosses a threshold. A delicate balance between the cost of a potential error catastrophe and the benefits of outracing the immune response enables tumor cells to survive and expand despite immune system attacks.
Several mathematical models have been designed to provide insights into the dynamics of tumorigenesis under immunosurveillance or an error catastrophe of tumor cells [34,35,36,37,38,39,40,41]. Most studies of the effects of immunosurveillance on tumor evolution considered a homogeneous population of tumor cells and concentrated on phenomena such as tumor dormancy and immunoescape. Studies of the tumor error catastrophe, in contrast, investigated quasispecies models in simplified situations without an immune response. The dynamics of tumor immunoescape and error thresholds, however, result from an interaction between both components; such studies are still lacking from the literature. In this paper, we investigate an integrated model of both concepts during tumor progression -the effects of tumor immunosurveillance and the consequences of a mutator phenotype of tumor cells. We introduce specific immune responses to a formulation of the quasispecies model and study the balance between evasion of immunosurveillance and prevention of an error catastrophe. This study reveals the effects of various tumor antigens on specific immune responses from the viewpoint of evolutionary dynamics, and provides new perspectives on optimum treatment strategies of tumors subjected to immunosurveillance and -editing.

Results
During the early phases of tumorigenesis, immune system cells such as NK and CD8+ T cells attack tumor cells and may succeed in suppressing their expansion; this outcome is referred to as ''tumor dormancy''. However, if the immune system cannot successfully eradicate a tumor, then eventually a subset of tumor cells will acquire the phenotypes necessary for immunoescape. Depending on the magnitude of the mutation rate of these cells, the tumor cell population may then be at risk of going extinct due to the generation of excess damage -the event of an error catastrophe. To investigate the dynamics, conditions, and likelihood of these events, we designed a mathematical model of tumor and immune system cells.
In the context of our mathematical model, initially there is only a single type of tumor cells -those cells that originally founded the tumor. Denote the abundance of these original tumor cells by x 0 . They divide at rate f 0 and die at rate d. During each division of such a cell, a new variant tumor cell is produced with probability u. Each tumor variant may have evolved a phenotype which allows it to evade eradication by immune system cells. The different tumor variants are enumerated as cell types i, and the abundance of each type is given by x i . Let us first assume that all tumor variants divide at rate a and die at rate d. These assumptions ensure that intra-variant competition is stronger than inter-variant competition, such that each tumor cell only competes with cells of its own type for oxygen, nutrients, and space; these assumptions will be relaxed later on. The total number of tumor variants is denoted by N. Since we consider a maximum number of N tumor variants, the growth rate of the original tumor cells reduces to f 0 (1{Nu).
In addition to tumor cells, we also consider immune system cells that launch a specific immune response against each particular tumor variant. Denote the abundance of immune system cells specific to the original tumor clone by y 0 , and the abundance of those specific to tumor variant i by y i , for i~1,:::,N. These immune system cells inhibit tumor variants at rate p and are generated by interactions with the tumor cells at rate c. Figure 1 displays a schematic representation of this framework. Specific immune responses such as CTLs recognize their target tumor cells through random interactions and identification of antigens presented on the cell surface [42,43]. In the context of the mathematical model, we assume that immune system cells encounter tumor cells at a rate proportional to the latter cells' frequency, x i =(szX ); the parameter s represents the coefficient of interactions between immune and tumor cells, and X~P i x i represents the total number of tumor cells. Expansion and differentiation of specific immune responses (e.g., precursor CTL proliferation and their differentiation into effector CTLs) are also regulated through interactions with tumor cells [44]. Therefore,

Author Summary
Immunologic surveillance is a function of the immune system which serves to constantly monitor the body for microorganisms, foreign tissue, and cancer cells. To evade this surveillance and subsequent elimination, cancer cells evolve strategies to prevent being recognized and killed by immune system cells; one mechanism is to increase the rate at which genetic and/or epigenetic variability is generated. The benefits of an increased variability of cancer cells to counteract immune surveillance, however, stands in contrast to the costs associated with such heightened mutation rates: the risk of an inability to maintain essential genome structure and functions. To study such situations arising in tumorigenesis, we designed a novel mathematical framework of tumor immunosurveillance and the evolution of mutation rates. We then utilized this framework to study how increased mutation rates and immunologic surveillance affect the abundance of tumor and immune system cells. We found that there exists an optimum number of tumor variants and an optimum magnitude of mutation rates that maximize tumor progression despite the presence of actively proliferating and functioning immune system cells. Our study contributes to an understanding of cancer development during immune system attacks and also suggests treatment strategies for heterogeneous tumor cell populations.
immune system elimination of tumor cells and the proliferation of CTLs are described in a frequency-dependent way, with maximum elimination and proliferation rates of p and c, respectively. The lifespan of CTLs is considered to be exponentially distributed with mean 1=d days.

The basic mathematical model
With these considerations, we define the basic mathematical model including tumor variants and their specific immune responses by ð1Þ Baseline values of model parameters and their respective ranges used for simulations are presented in Table 1. A subset of these parameter values were estimated in [39,45,46]. Since the original tumor cell clone has been suggested to proliferate faster than variant tumor cells [37,47], the division rate of variant tumor cells is a~ef 0 where 0vev1; we assume a default value of e~0:05 but also perform sensitivity analyses (see a later section). We found that, although some of these values are rough estimates and might deviate when measured by other groups or in other systems, our main results are qualitatively preserved within broad ranges around our baseline values.
Let us now discuss the possible outcomes of interactions between the immune system and the tumor cell population: there may be tumor dormancy, partial immunoescape, complete immunoescape, and the event of an error catastrophe. In the dormancy state, immunosurveillance serves to effectively suppress the tumor cell population. In the partial immunoescape state, some tumor variants (but not all) achieve immunoescape while in the complete immunoescape state, the immune response is completely unsuccessful. Finally, in the error catastrophe state, the original tumor clone, which has the highest division rate, goes extinct due to the accumulation of excess alterations. We now outline how the original tumor clone, the tumor cell variants, and the specific immune system cells behave during the accumulation of alterations and the evolution of higher mutation rates.

Tumor immunoescape and error catastrophe
The four qualitative outcomes of the interaction between tumor cells and the immune system -dormancy, partial and complete immunoescape as well as error catastrophe -are most significantly influenced by two systems parameters: the mutation rate generating tumor variants (u) and the maximum number of tumor variant types that can emerge (N). We therefore investigated the dynamics of tumor evolution in dependence of these parameters, and identified three analytical thresholds (N Ã 1 , N Ã 2 , and N Ã 3 ) separating the potential outcomes ( Figure 2). The formulas and detailed mathematical analyses of these thresholds are provided in the Methods section.
As long as the number of tumor variants is less than the first threshold, N Ã 1 , immune responses suppress all tumor variants (tumor dormancy). When the number of variants exceeds this threshold, however, then some tumor cells escape from the specific immune response (partial immunoescape). Once the number of variants passes the second threshold, N Ã 2 , all tumor cells escape from immune responses (complete immunoescape). This finding implies that tumor cells can evade immune surveillance by accumulating a sufficiently large extent of intratumor heterogeneity. However, if the number of variants exceeds the third threshold, N Ã 3 , then an error catastrophe of tumor cells occurs, in which the original tumor clone can no longer maintain an expanding population and the original tumor cells therefore go extinct. We also found that, as the mutation increases, the threshold N Ã 1 increases while N Ã 2 and N Ã 3 decrease. In all scenarios, however, tumor eradication is unlikely -although the tumor cell burden may shrink by a large amount -when the growth rate of the original tumor clone is negligibly small as compared to their death rate by apoptosis and/or interactions with the immune system.
We have thus established that although high rates of accumulating alterations allow tumor cells to reach a state of complete immunoescape, those cells with an excessively high mutation rate suffer an error catastrophe as the number of tumor variant types increases. These systems dynamics suggest that there is an optimum amount of instability that optimizes tumor evolution (i.e. maximizes the number of tumor cells) while maintaining a functioning genome. given by x i . Tumor variants divide at rate a and die at rate d. The total number of tumor variant types is denoted by N. Since we consider N tumor variant types, the growth rate of the original tumor cells reduces to f 0 (1{Nu). In addition to tumor cells, we also consider immune system cells that launch a specific immune response against each particular tumor variant. Denote the abundance of immune system cells specific to the original tumor clone by y 0 , and those specific to tumor variant i by y i for i~1,:::,N. These immune system cells inhibit tumor variants at rate p and are generated by interactions with the tumor cells at rate c. Immune system cells encounter tumor cells at a rate proportional to the latter cells' frequency, x i =(szX ); the parameter s represents the coefficient of interactions between immune and tumor cells, and X~P i x i represents the total number of tumor cells. The lifespan of immune system cells is exponentially distributed with mean 1=d days. doi:10.1371/journal.pcbi.1002370.g001 The optimal rate of tumor evolution Let us now investigate the system dynamics for varying mutation rates and identify those regimes in which the total tumor cell number is maximized. Every time a new tumor variant arises, the dynamics of tumor evolution rapidly converges to its steady state; we therefore analyze the dynamics in steady state. The total number of tumor cells depends on the number of variants as well as the mutation rate, and an optimum combination of these parameter values exists that maximizes the total tumor cell number. In Figure 3, we demonstrate how the total number of tumor cells is affected by the number of tumor variant types for three different cases in which the mutation rate is u~0:004, u~0:01, and u~0:025, respectively. Detailed mathematical analyses of those equilibria are provided in the Methods section.   In situations in which all tumor cells are effectively suppressed by the immune response (tumor dormancy), the total number of tumor cells increases with the number of variant types. In situations in which some tumor cell types manage to escape from immune surveillance, the total number of tumor cells increases as both the number of variant types and the mutation rate increase (Figure 3). However, in situations in which all tumor cell types completely escape from their specific immune responses, there exist two thresholds regarding the total number of tumor cells: N c and u c (see Figure 2). In this scenario, the total number of tumor cells increases until the number of variant types and the mutation rate, respectively, exceed the values of N c and u c ; once crossing these thresholds, the tumor cell number decreases as N and u further increase. Therefore, tumor cells with an excessively high mutation rate cannot continue to become more abundant as the number of variant types increases (Figure 3), but there is an optimum, nontrivial parameter regime that maximizes the number of tumor cells.
Our results demonstrate that there are two strategies to maximize the rate of tumor evolution so that the total tumor cell mass is maximally large: one is to maintain a low mutation rate, since then the tumor cell population can increase the number of variant types along the threshold u c (see Figure 2); another is to keep the number of variant types relatively small, since then the tumor cell population can increase the mutation rate along the threshold N c (see Figure 2). When both the mutation rate and the number of variant types are large, then the tumor cell population cannot maintain its maximum number without decreasing one of the two parameters.

Fitness of variant tumor cells
Let us now investigate how the division rate of variant tumor cells affects the evolution of tumor cells during their interaction with immune system cells. Recall that in the basic model, the division rate is a~ef 0 , and that the threshold for an error catastrophe to occur (N Ã 3 ) is independent of the division rate. To investigate the dependence of the system behavior on this division rate, we chose four different e and performed a sensitivity analysis for the thresholds N Ã 1 and N Ã 2 . Figure 4 displays how the division rate of variant tumor cells influences the outcome of tumor immunoescape in the plane of mutation rates (u) and the number of tumor variants (N). The four panels of the figure represent cases with different values of e. The higher the fitness of variant tumor cells becomes, the more easily they escape from immunosurveillance. However, the qualitative profiles of the system dynamics are preserved; that is, tumor cells with high mutation rates tend to reach a complete immunoescape while tumor cells with low mutation rates effectively produce a diverse population and thus increase in number.

Extensions of the mathematical model
Let us next consider additional effects arising during tumor progression such as competition among tumor cells of different variant types, the presence of an innate immune response such as NK cells, which non-specifically target all tumor variants, and differential growth rates among tumor cell variants. In order to investigate the conditions for outcomes such as tumor immunoescape and error catastrophe in these more complex scenarios, we established an extended model, given by where i~1,:::,N. The parameter a i represents the division rate of tumor variant i. We now assume that each tumor cell competes with all other tumor cells so that the death term becomes dx i X . Furthermore, the variable z describes innate immune responses, for instance by NK cells which attack tumor cells without antigen specificity. The parameters g, r, k, and s represent, respectively, the maximum proliferation rate of NK cells, the maximum elimination rate of tumor cells by NK cells, the decay rate of NK cells, and a constant source of NK cells.
The dynamics of tumor progression considering these situations are shown in Figure 5. We investigated how intervariant tumor cell competition ( Figure 5A), incorporation of an innate immune response ( Figure 5B), growth rates which differ between individual tumor variants ( Figure 5C), and all three effects simultaneously modulate the thresholds between outcomes as well as the optimum parameter regimes for maximizing tumor cell numbers. Competition among tumor cells of the same variant type renders it difficult for the tumor cell population to completely escape from immune surveillance and to increase the total cell number beyond a relatively small value, irrespective of the mutation rate ( Figure 5A). However, when only an innate immune response is present without intervariant competition, then there is a larger parameter regime in which complete immunoescape is possible. Furthermore, the total number of tumor cells is larger in this situation as compared to the above case ( Figure 5B). Similar to this scenario, the presence of different growth rates for individual variant clones allows for the existence of a large number of tumor cells as well as a large regime in which complete immune escape can be achieved ( Figure 5C). Finally, when all three aspects are combined in the mathematical model, then the region of complete immune escape becomes very small; this effect is mainly driven by the incorporation of interal competition. The total number of tumor cells also remains below a rather small threshold for this case ( Figure 5D).

Optimum treatment strategies for diverse tumor cell populations
Finally, let us discuss the effects of different treatment modalities on the rates of cancer progression and the chance of immunoescape. Since the behavior of tumor cells and thus patient outcomes are to a considerable extent driven by the interactions between tumor and immune system cells, we considered both traditional chemotherapy and treatment options that stimulate the immune system to launch or sustain an attack against the tumor cell population. In general, immune therapies have not been proven to be very effective against many tumor types; one of the few exceptions is represented by adoptive cell therapy, which is used in the treatment of metastatic melanoma and causes regressions in about 50% of patients [48]. Recently, however, synergistic effects of immunotherapy in combination with chemotherapy have been reported in both human and animal trials [49,50,51], and several mechanisms were identified that may explain these synergistic effects [52].
To study the effects of chemotherapy, immune therapy, and combination therapy on the dynamics of tumor evolution, we introduced a series of different treatment types into the mathematical framework and identified optimal treatment strategies for diverse tumor cell populations ( Figure 6). These different treatment modalities were tested in situations in which tumor cells had previously achieved complete immunoescape and consisted of a large number of tumor cells. The number of tumor variants and the mutation rate were considered to be N~50 and u~0:01 at the time of treatment initiation. Chemotherapy then reduces the number of tumor variants and kills tumor cells proportional to the tumor cell number present. We also considered the case in which chemotherapy reduces the growth rate of tumor cells. Then the model after treatment initiation is given by Here chemotherapy reduces the number of tumor variants to N' and either kills the tumor cells at rate b or reduces the growth rates to f ' 0 and a'. Immunotherapy increases the number of specific immune cells (y i for i~0,:::,N) during treatment. We then utilized this system to investigate optimum treatment strategies. First, let us consider the effects of chemotherapeutic agents which reduce the number of tumor cells by inducing cell deaths at a rate proportional to the cell number present within the tumor. Administration of such treatments decreases the total cell number, but may not be capable of leading to complete eradication of all tumor cells ( Figure 6A) unless its effects are sufficiently (and maybe unrealistically) strong ( Figure 6B). Second, consider chemotherapeutic drugs which reduce the number of tumor variant types as well as the growth rates of tumor cells. Again, administration of such treatments decreases the total cell number but is incapable of achieving complete eradication of tumor cells ( Figure 6C). Third, consider the administration of immunotherapy which increases the population of tumor-specific immune system cells. Such therapy alone is not able to decrease the abundance of tumor cells by a large extent ( Figure 6D). However, when combining chemotherapy and immunotherapy, an effective decrease of the tumor cell population can be achieved, which may ultimately lead to tumor eradication and a cure ( Figure 6E). Notably, in situations in which the mutation rate is small, the administration of combination therapy is more successful in eradicating all tumor cells as compared to situations in which the mutation rate is high ( Figure 6F).
In conclusion, our mathematical model predicts successful outcomes of combination therapy when (i) chemotherapy is administered which induces tumor cell death at a significantly large rate, or (ii) combination therapy is administered which reduces the number of tumor variants, induces tumor cell death, and replenishes immune cell populations. When the mutation rate of tumor cells is small, combination therapy is more effective than when variations arise at a large rate. An explanation of these findings can be found in Figure 2 -activation of the immune response alone does not change the state of the tumor cell population once it has reached complete immunoescape; in that case, the number of tumor cells does not decrease ( Figure 6D). A reduction of the number of tumor variant types and tumor cells by administering chemotherapy alone allows for partial immunoescape or dormancy states, but there is an insufficient abundance of immune system cells to effectively control the tumor cell population ( Figure 6A). However, combination therapy which enables immune cells to be activated in the states of partial immunoescape or dormancy is capable of eradicating the tumor ( Figure 6E). Thus, our mathematical framework is capable of identifying those treatment modalities that have the potential to lead to a cure of the tumor.

Discussion
In this paper, we have investigated the dynamics of tumor progression under immune system surveillance while considering the effects of increasing rates at which (epi)genetic alterations are generated. We defined specific situations that can arise due to the interactions of immune system cells and tumor cells. When the tumor cell population is able to persist under immunosurveillance without leading to tumor growth, then a state of tumor dormancy ensues. Should the immune system not be capable of efficiently suppressing the tumor cell population, then partial or complete immunoescape is possible, depending on whether some or all tumor clones evade immune system inhibition. Finally, an error catastrophe occurs when the tumor cells evolve mutation rates that are incompatible with the maintenance of a functioning genome due to excess error.
The dynamics of the system and likelihood of these different states depend on the rate at which variability emerges in the population (denoted by the mutation rate u per cell division) as well as the number of distinct tumor clones (given by N) that are distinguished by their capabilities of generating a specific immune response (see Figure 2). If both quantities are excessively large, then an error catastrophe occurs and the original tumor cell population goes extinct. In intermediate regimes, states of dormancy and partial or complete immunoescape are possible. We also investigated the extent to which the total number of tumor cells depends on these parameters and identified regimes in which the maximum number of tumor cells is attained. Moreover, we relaxed the model assumptions to consider more complex scenarios such as growth competition among tumor variants, innate immune responses that non-specifically recognize and kill tumor cells, and different growth rates of tumor variants. These studies revealed that the patterns of states do not vary significantly as the assumptions of competition, growth, and innate immune responses are altered; however, internal competition among tumor variants renders it difficult for tumor cells to achieve complete immune escape.
Finally, we investigated the effects of different treatment modalities on the rates of tumor progression and found that administration of both chemotherapy and immunotherapy leads to optimum response rates, thereby confirming recent experimental findings [49,50,51]. These investigations have direct implications for the clinical management of cancers since they incorporate both mutator phenotype and the interactions between tumor cells and the immune system. A consideration of these factors is essential for an understanding of the dynamics of tumor cell populations evolving during immune system attacks. Our results thus suggest that combination therapy incorporating both chemotherapeutic and immunostimulatory agents would lead to best results in the clinic.
Our mathematical framework represents only one possibility of modeling the system of tumor and immune system cells. This modeling choice was made for reasons of mathematical simplicity as well as availability of parameter estimates; however, several model extensions are conceivable. For instance, the spatial components of the system could be incorporated such that the spatio-temporal aspects can be considered. Also, we have neglected stochasticity in our formulation of the mathematical model since both population sizes and mutation rates are large, and therefore deterministic dynamics dominate. However, for more detailed investigations -such as a determination of the probability that a certain phenotype arises -the stochasticity of the system cannot be neglected. Such studies will be the topic of future contributions. Furthermore, interactions with the microenvironment such as with stromal cells and other factors could be considered. The incorporation of these extensions are complicated by the fact that few quantitative estimates are available. The determination of system parameters necessary for including such factors into a mathematical framework is an important goal of the field.

Mathematical analysis
Let us first consider the basic model, equation (1), in detail. The basic model can be considered qualitatively as a 4-dimensional ODE system (although equation (1) is a 2(Nz1)-dimensional ODE system) when analyzing the equilibria, since all parameters in the equations describing x i and y i are the same, and therefore x i and y i have the same properties at the equilibria. Hence we can consider that X~x 0 zNx i with respect to the equilibria, where N is a model parameter.

Existence conditions of equilibria
We investigated the existence conditions of the equilibria of model (1). The model has seven possible equilibria: E 0~( 0,0,0,0); and x x i is the positive root of the following equation: The existence condition of E p is ff 0 (1{Nu)(c{d){ddsg=ddNw x x i . In addition, we calculate the existence conditions of E z as follows.
Note that x z 0~x z i w0 and y z i w0 are equivalent to Nv(c{d)=d and Nv(c{d)=d{ds=(azf 0 u), respectively. If Nv(c{d)=d, then we have the following relations: because G(1=u)v0, G(N)~0 always has real roots N c and N N c (N c v N N c ). Here we assume that f 0 (c{d){ddsw0: otherwise E z never exists and E c becomes stable. Therefore, y z 0 w0 is equivalent to NvN c . Note that the roots of G(N)j d~0~0 are (c{d)=d and 1=u, which implies that N c v(c{d)=d, 1=u. Furthermore, when u~0, G(N)w0 is equivalent to Nv(c{d)=d{ds=f 0 . Since (c{d)=d{ds=av(c{d)=d{ds=f 0 in the context of f 0 wa (which is a suitable assumption), we can roughly estimate that (c{d)=d{ds=(azf 0 u)vN c if uvv1. Consequently, if Nv(c{d)=d{ds=(azf 0 U), f 0 (c{d){ddsw0, f 0 wa and uvv1, then E z exists.

Derivation of thresholds
Consider the situation in which E z exists, all tumor cells are suppressed by their specific immune responses, and the number of tumor variants is 1 (i.e., N~1). Let us investigate the transversal eigenvalue in the y i {direction at E p . Here we define Since sz x x 0 zN x x i~c x x 0 =d, the transversal eigenvalue is evaluated as follows: This result implies that, if N Ã 1 vN, then E z disappears from R 4 z and y i approaches 0 near E p . Hence only the original tumor cells are suppressed by their specific immune responses while tumor variants escape from immune surveillance once the number of tumor variant exceeds N Ã 1 . Next, we investigate the transversal eigenvalue in the y 0 {direction at E c . Here we define N Ã 2 to satisfy The transversal eigenvalue is evaluated as follows: This result implies that N Ã 2 vN is equivalent to Ly' 0 =Ly 0 v0 near E c ; that is, if N Ã 2 vN, then E p disappears from R 4 z and y 0 converges to 0 near E c (E c is stable). Thus, the immune response against the original tumor clone also becomes inactivated if the number of tumor variants exceeds N Ã 2 . Note that we consider a restricted region of u and N, where 2(1{Nu)f 0 (c{d) 2 {d(c{d)(aNz2ds){f 0 ud 2 N 2 w0. Consequently, all tumor cells escape from their immune responses once the number of tumor variants exceeds N Ã 2 . Furthermore, the original tumor cell clone is no longer sustainable (i.e., an error catastrophe occurs) as soon as the number of the variants exceeds N Ã 3~1 =u.

The total number of tumor cells
Let us now calculate the total number of tumor cells at equilibrium. Note that the dynamics of the basic model, equation (1), might not converge to an equilibrium, but may oscillate if NvN Ã 2 . When the number of tumor variants is 0vNvN Ã 1 , then the total number of tumor cells at E z is X z~( Nz1)ds=fc{(Nz1)dg. Since we evaluate LX z =LNw0 and LX z =Lu~0, the total number of tumor cells increases as the number of variants grows. When the number of variants is N Ã 1 vNvN Ã 2 , then the total number of tumor cells at E p is X X~(cN x x i zds)=(c{d). Again, as we evaluate L X X =LNw0 and L X X =Luw0, the total number increases as the number of variants and the mutation rate increase. When the number of variants is N Ã 2 vNvN Ã 3 , then the total number of tumor cells at E c isX X~f2f 0 (1{Nu)zaNzN ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a 2 z4f 2 0 (1{Nu)u q g=2d.
Here we find critical thresholds N c and u c such as LX X =LNj N~Nc~0 and LX X =Luj u~uc~0 . Therefore, the total number of tumor cells increases as long as the number of variants and the mutation rate are N Ã 3 vNvN c and uvu c , respectively. However, once the number of variants and the mutation rate exceed N c and u c , respectively, the total number of tumor cells decreases. Eventually, when the number of variants is N Ã 3 vN, the total number of tumor cells at E e is X e~a N=d.