A Cellular Automaton Model for Tumor Dormancy: Emergence of a Proliferative Switch

Malignant cancers that lead to fatal outcomes for patients may remain dormant for very long periods of time. Although individual mechanisms such as cellular dormancy, angiogenic dormancy and immunosurveillance have been proposed, a comprehensive understanding of cancer dormancy and the “switch” from a dormant to a proliferative state still needs to be strengthened from both a basic and clinical point of view. Computational modeling enables one to explore a variety of scenarios for possible but realistic microscopic dormancy mechanisms and their predicted outcomes. The aim of this paper is to devise such a predictive computational model of dormancy with an emergent “switch” behavior. Specifically, we generalize a previous cellular automaton (CA) model for proliferative growth of solid tumor that now incorporates a variety of cell-level tumor-host interactions and different mechanisms for tumor dormancy, for example the effects of the immune system. Our new CA rules induce a natural “competition” between the tumor and tumor suppression factors in the microenvironment. This competition either results in a “stalemate” for a period of time in which the tumor either eventually wins (spontaneously emerges) or is eradicated; or it leads to a situation in which the tumor is eradicated before such a “stalemate” could ever develop. We also predict that if the number of actively dividing cells within the proliferative rim of the tumor reaches a critical, yet low level, the dormant tumor has a high probability to resume rapid growth. Our findings may shed light on the fundamental understanding of cancer dormancy.

as cellular dormancy, angiogenic dormancy and immunosurveillance have been proposed, a comprehensive understanding of cancer dormancy and the "switch" from a dormant to a proliferative state still needs to be strengthened from both a basic and clinical point of view.Computational modeling enables one to explore a variety of scenarios for possible but realistic microscopic dormancy mechanisms and their predicted outcomes.The aim of this paper is to devise such a predictive computational model of dormancy with an emergent "switch" behavior.Specifically, we generalize a previous cellular automaton (CA) model for proliferative growth of solid tumor that now incorporates a variety of cell-level tumor-host interactions and different mechanisms for tumor dormancy, for example the effects of the immune system.Our new CA rules induce a natural "competition" between the tumor and tumor suppression factors in the microenvironment.This competition either results in a "stalemate" for a period of time in which the tumor either eventually wins (spontaneously emerges) or is eradicated; or it leads to a situation in which the tumor is eradicated before such a "stalemate" could ever develop.We also predict that if the num-

Introduction
Cancer dormancy, the phenomena that the tumor's volume or the number of tumor cells stays at a very low level for a certain period of time before the tumor begins to grow rapidly, has been an outstanding issue in cancer research for many years [1,2].Currently, the mechanisms responsible for the "switch" from a dormant state to a rapid growth state for different tumors are not well understood, although it is well known that such a "switch" in secondary metastatic tumors can be triggered by the removal of the primary tumor.This could eventually lead to failure of tumor treatment and fatal outcomes for the patient.Therefore, a comprehensive understanding of the "switch" from a dormant to a proliferative state is crucial to our fundamental understanding of cancer progression and recurrence and might lead to the development of novel treatments for cancer.Dormancy has been observed in many types of cancer.This includes tumor dormancy before any metastases take place and the latency of cancer recurrence after therapy.In some cases of pancreatic cancer, the tumor can remain in a benign dormant state for about 20 years [3].During this time, it is undetectable by conventional clinical methods, and it is only afterwards that the tumor becomes highly malignant and grows aggressively with highly fatal outcomes after about a year.In the cases of breast and prostate cancer, it is reported that 20%-45% of patients will relapse years or decades later after the resection of the primary tumor [4][5][6].In addition, recurrence has been observed in brain tumors, which indicates the existence of a large number of micrometastases that are dormant in the presence of the primary tumor [7,8].
Extensive studies over years have revealed three major cancer dormancy mechanisms: cellular dormancy, angiogenic dormancy and immunosurveillance [1,2].On the cellular level, a tumor cell could be arrested at a certain stage of the cell cycle and unable to complete the cell division process successfully, resulting in a dormant solitary cell [9][10][11].
On the cell population level, when the population does not gain enough ability to recruit blood vessels and promote neovascularization, the tumor cannot obtain sufficient nutrients necessary for its proliferation and as a result, angiogenic dormancy occurs [12,13].On other hand, immunosurveillance operates when the immune system suppresses the proliferation of tumor cell population and leads to the dormancy of the tumor [14][15][16][17][18]. Figure 1 A comprehensive understanding of cancer dormancy and the "switch" from a dormant to a proliferative state still needs to be strengthened.This is mainly due to the fact that efficient and accurate experimental or clinical approaches to track the states of individual cells in a dormant tumor in vivo throughout the entire dormancy period are still under development [19][20][21].
Given the current need for further understanding of dormancy, computational modeling provides a powerful means to probe various scenarios for the underlying mechanisms.Specifically, modeling enables one to probe a variety of different dormancy scenarios by examining different combinations of mechanisms in order to see which ones provide possible explanations for experimental and clinical observations.Over the past few decades, computational modeling has played an important role in the study of the progression of solid tumors [22]; a variety of models based on different mathematical schemes have been developed, including continuum models [23][24][25][26], discrete cell models [27,28] and hybrid models [29,30].Various models have been used to investigate cancer dormancy caused by cancer-immune interactions and other mechanisms, including ordinary differential equation-based models [31,32], stochastic differen-tial equation-based models [33], models based on kinetic theory for active particles [34][35][36], and cellular automaton models [37].However, the aforementioned studies neither explicitly demonstrated how the dynamic process of active proliferation after a certain period of dormancy emerges from various microscopic mechanisms nor showed the associated growth dynamics of the "switch" phenomenon.
Therefore, predictive computational models that incorporate cellular-level microscopic mechanisms are needed to address these important issues.

Materials and Methods
We divide the two-dimensional square simulation box into different polygonal units (i.e., automaton cells).Our model is coarse-grained, allowing us to grow the tumor from a very small size with a cross section of roughly 1000 real cells through to a fully developed tumor with a cross section consisting of 2.0 × 10 6 cells.Specifically, the innermost automaton cells represent roughly 100 real tumor cells or or a region of host microenvironment of similar size, while the outermost automaton cells represent roughly 10 4 real tumor cells or a region of host microenvironment of similar size.To generate the automaton cells in the simulation box, we first fill the simulation box with non-overlapping circular disks (or spheres in 3D) using random-sequentialaddition packing method [43] until there is no void space left for additional circular disks (or spheres in 3D).Periodic boundary conditions are used for generating the packing.Then we divide the simulation box into polygons (or polyhedra in 3D), each polygon (or polyhedron in 3D) associated with a particle center, such that any point within a polygon (or polyhedron in 3D) [i.e., a Voronoi polygon (or polyhedron in 3D)] is closer to its associated particle center than to any other particle centers.The resulting Voronoi polygons (or polyhedra in 3D) are referred to as automaton cells.In this paper, we focus on the two-dimensional case, but our model should be readily generalized to three dimensions.In this section, we will introduce our CA dormancy model, which modifies our previous basic CA models of tumor growth [41,42,44]  • Quiescent cells more than a certain distance δ n from the tumor's edge are turned necrotic.
The tumor's edge, which is assumed to be the source of nutrients, consists of all ECMassociated cells that border the tumor.The critical distance δ n for quiescent cells to turn necrotic is computed as follows: where a is the necrotic thickness controlled by nutritional needs, d is the Euclidean spatial dimension and L t is the distance between the geometric centroid x c of the tumor (i.e., , where N is the total number of cells in the tumor) and the tumor edge cell that is closest to the quiescent cell under consideration.
• Proliferative cells more than a certain distance δ p from the tumor's edge are turned quiescent.
The critical distance δ p is given by where b is the proliferative thickness controlled by nutritional needs, d is the spatial dimension and L t is the distance between the geometric tumor centroid x c and the tumor edge cell that is closest to the proliferative cell under consideration.
• The probability of division for a proliferative cell used in our model is where p 0 = 0.192 is the base probability of division linked to cell-doubling time, ρ ECM is the local ECM density, ω * = 2ρ 0 ECM is a parameter taking into account the effect of pressure, ξ = ρ ECM /ρ 0 ECM is the ratio of current average ECM density over the initial density, and ℓ and w are, respectively, the length and width of local protrusion tips.

Interactions between the tumor and the microenvironmental suppression factors
Here, we specify the additional interaction rules between the tumor and the microenvironmental suppression factors beyond the aforementioned ones for noninvasive proliferative growth, which were not included in our previous CA models.We assume that there are two possible states of proliferative cells, dormant or actively dividing, depending on their interactions with the microenvironmental suppression factors.
• Initially, we assume that all proliferative cells are kept in dormant states by the microenvironmental suppression factors, which means that they are not able to divide.For example, mutated leukaemic cells in acute myeloid leukaemia acquire resistance to cytotoxic T lymphocytes-mediated cell lysis, whose degree is related to the level of the cell's expression of B7-H1 or B7.1 [45].For simplicity, we divide the phenotypic changes into two different types: weak changes and strong changes with respect to their resistance to the suppression factors in the microenvironment (i.e.their ability to actively proliferate).Henceforth, we will refer to these phenotypic changes as "transformations" and the cells that undergo these changes as "transformed" cells for simplicity.Strong-type "transformed" cells gain a larger competition advantage and thus have a greater ability to divide actively.The quantities x W and x S are the fractions of weak-type "transformations" and strong-type "transformations".Henceforth, we set x W 0.99 and x S 0.01.
• At each subsequent day, the microenvironmental suppression factors will counteract the weak-type "transformed" and strong-type "transformed" cells with probabilities α W and α S .The suppression factors in the microenvironment will either kill the "transformed" cells or turn them back into dormant cells [14,15].[46,47].Therefore, when the microenvironmen-tal factors counteract the "transformed" cells, the fraction of the cells that are killed can be coupled with the growth rate of the tumor by where k 0 is a constant characterizing the viding tumor cells and k is set to be 0 [46,47].
• Due to the "competition" between the tumor and the suppression factors in the microenvironment, the ratio of the number of actively dividing proliferative cells over the total number of proliferative cells n pro acti /n pro changes with time.The larger is this ratio n pro acti /n pro , the larger is the amount of nutrients the tumor tissue consumes.As a result, the nutrient concentration around the tumor depends on the ratio n pro acti /n pro .Therefore, we make the necrotic thickness a and proliferative thickness b functions of n pro acti /n pro : where a 0 = 0.58 mm 1/2 and b 0 = 0.30 mm • At each day, each dormant proliferate cell is checked to see if it enters the active state according to the interaction rules.Once it begins to actively divide, it proliferates according to the CA rules for proliferative tumor growth.
• At each day, each active proliferative cell is checked to see if it is killed or turned back into dormant cell according to the interaction rules.
• Quiescent cells and necrotic cells act accord-ing to CA rules for proliferative tumor growth.
However, the values of parameters a and b determining the transitions from necrotic cells to quiescent cells and from proliferative cells to quiescent cells respectively are influenced by interaction rules, as mentioned above.
Note that our CA dormancy model should be readily generalized to angiogenic dormancy by explicitly considering the angiogenic process and vascular tumor growth.This is beyond the scope of this work and will be addressed in future work.

Results
In this section, we apply our CA model and show that it produces a dormancy period of the tumor that can lead to a subsequent emergent "switch" behavior to a proliferative state.A homogeneous distribution of ECM density is used for simplicity [41].A circular growth permitting region containing

Statistics of tumor growth
Here we consider the growth of a proliferative tumor

Suppression rate vs transformation rate
Here we investigate growth dynamics of the tumor under different suppression rates α and phenotypic transformation rates γ.The suppression rate α is defined as the following weighted average: proliferative and dormant regions.By "proliferative", we mean that the tumor resumes rapid proliferation after a period of dormancy and the length of the dormancy period is less than a virtual patient's life; by "dormant", we mean that the tumor remains in a dormant state during the whole life of a person and undetected by conventional clinical methods (usually clinicians call such tumors "benign" [48][49][50]).The solid line separates the two regions, and crossing this boundary line is associated with a "phase transition".

Rigidity of the microenvironment
Various mechanical cues in the microenvironment could influence the growth dynamics of the tumor from those reported here [51].

Strength of the suppression factors
Here we investigate how tumor growth dynamics changes with the strength of the microenvironmental suppression factors.As shown in Figure 6, increasing the fraction of actively dividing tumor cells that are killed [i.e., increasing k 0 in the equation ( 4)] when the microenvironmental suppression factors (which we recall could either kill the "trans-  Specifically, our work shows that even if the immune system maintains its strength throughout the tumor growth process, there is still a high possibility that the immune system could eventually fail, which is to be contrasted with the simple explanations that it becomes weaker as the tumor develops [15][16][17][18]. (a)    shows an image of tumor tissue surrounded by im-mune cells.Figure1(b) compares the morphology and vascular structure of dormant and fast-growing tumors.
in order to investigate tumor dormancy.Our goal is to formulate a dynamical model in which the "switch" to a proliferative state spontaneously emerges by incorporating additional interactions between the tumor and the microenvironment, for example the effects of immune system, which were not included in our previous CA model.The new rules of our CA model induce a "competition" between the tumor's propensity to proliferate and the microenvironmental factors that suppress its growth.In our model, a fraction of the dormant cells undergo phenotypic transformations triggered by intracellular factors or external stimulus and acquire the ability to actively proliferate.Subsequently, those microenvironmental factors act to suppress the growth of these transformed tumor cells either by killing some of these cells or turning these actively dividing proliferative cells back into dormant cells.The "competition" between the tumor and the microenvironmental suppression factors either results in a "stalemate" for a period of time in which the tumor either eventually wins (spontaneously emerges) or is eradicated; or it leads to a situation in which the tumor is eradicated before such a "stalemate" could ever develop.Since we are mainly interested in the situations in which tumor growth involves a period of dormancy, we will henceforth focus on those situations in which a "stalemate" between the tumor and the microenvironmental suppression factors develops.Our model demonstrates that a variety of parameters characterizing the tumor-host interactions may greatly alter the growth dynamics of the tumor.These parameters include the rate of phenotypic transformation, by which the tumor cells gain the ability to proliferate against those suppression factors, the suppression rate imposed by suppression factors on individual tumor cells, and the mechanical rigidity of the microenvironment.The growth dynamics influenced by these parameters include the existence of a dormant period in tumor's growth, the length of the dormant period (if there exists one) and the existence of a sudden "switch" to a highly proliferative state.We also demonstrate that if the number of actively dividing cells within the proliferative rim reaches a critical, yet low level, the tumor has a high probability to begin rapid proliferation.While we study a 2D CA model for simplicity in this paper, our model can be easily generalized to three dimensions (3D).
The microenvironment surrounding a tumor is mainly composed of stroma cells and extracellular matrix (ECM).In the current model, we explicitly take into account the effects of the ECM macromolecule density, ECM degradation by the proliferative cells, and the pressure built up due to the ECM deformation by tumor growth.The effects of the stroma cells are not explicitly considered in our current model.Henceforth, we will refer to the re-gions of microenvironment as ECM-associated cells for simplicity.In addition, we consider the interactions between the tumor and the various suppression factors in the microenvironment, for example the immune system.Since we consider development of primary tumor or local recurrences of micrometastases under microenvironmental suppression, invasive tumor growth is not a mechanism relevant for our purposes and hence is not included in our dormancy model.In our model of noninvasive proliferative tumor growth, tumor cells can be in one of the three possible states: proliferative, quiescent or necrotic, depending on their nutrient supply.Proliferative cells are tumor cells that have enough nutrients and possess the ability to divide.Quiescent (or arrested) cells are tumor cells that are alive, but do not have enough nutrient supply to support cell division.Quiescent cells can eventually become inert, necrotic (dead) cells due to an insufficient nutrient supply.In this paper, we focus on avascular tumor growth and assume that there is no explicit angiogenesis during the growth process (although this assumption can be relaxed).The nutrients available to tumor cells are those that diffuse into the tumor region through tumor edge.As the tumor grows, the amount of nutrient supply, which is proportional to the perimeter of the tumor interface (or surface area of the tumor interface in 3D), cannot meet the needs of all of the tumor cells.As a result, quiescent and necrotic regions emerge near the center of the tumor.The state of a tumor cell is determined by its distance to the tumor edge (i.e., the source of nutrients).We assume that proliferative cells more than δ p away from the tumor edge become quiescent and quiescent cells more than δ n away from the tumor edge become necrotic (see details below).
by introducing several additional parameters to incorporate the interactions between tumor cells and the microenvironmental suppression factors.This dynamical model is capable of producing situations in which a "switch" from a dormant state to a proliferative state spontaneously emerges.Noninvasive proliferative tumor growth We now specify the cellular automaton rules used in our model for noninvasive proliferative tumor growth.Each ECM-associated automaton cell is assigned a specific density ρ ECM , representing the density of the ECM molecules within the automaton cell.If a proliferative cell divides, its daughter cell occupies a nearby ECM-associated cell.The daughter cell pushes away or degrades the ECM within the ECM-associated cell it occupies.Initially, a tumor is introduced by designating several automaton cells at the center of the growth permitting region as proliferative tumor cells.Then time is discretized into units, with each time step representing one day.At each time step, the tumor grows according to the following cellular automaton rules.

•
At each day, beyond the aforementioned CA rules for proliferative noninvasive growth, each dormant proliferative cell has a certain probability γ to change in their phenotypes due to intracellular factors or external stimulus.The cell with phenotype change gains different degrees of resistance to the suppression factors in the microenvironment, depending on the specific phenotype change the cell undergoes.
When the number of tumor cells reaches a certain threshold N T , strong reactions of the microenvironmental factors are triggered and those factors start to kill the "transformed" cells.The parameter N T is introduced to ensure that the tumor is not completely removed by the microenvironmental suppression factors.Note that the particular choice of N T barely has any effect on the simulation results within a relatively wide range of N T values.In this work N T is set to be 50, a sufficiently small value that leads to biophysically realistic outcomes.As the tumor grows, the microenvironmental factors are weakened by the tumor, resulting in weaker suppression of the tumor cells of the suppression factors in the microenvironment, dA/dt is the daily area change of the tumor (i.e. the growth rate of the tumor), and △r C is the critical value of the tumor's growth rate.In this work, △r C is chosen as half of the tumor's maximum growth rate under suppression, but our numerical tests have revealed that the simulation results are insensitive to the choice of △r C as long as △r C is smaller than the tumor's maximum growth rate.When the growth rate of the tumor reaches this critical value, the suppression factors become too weak to kill any actively di-

∼ 2 ×
10 4 automaton cells is employed.Simulating the growth of a 2D tumor from several cells (representing roughly 1000 real cells) to a macroscopicsize tumor (a cross section of 5 cm 2 consisting of ∼ 2 × 10 6 real cells) with a period of dormancy up to a person's life (∼ 80 years) generally takes no more than a few minutes on a standard Dell Workstation (Precision T3400).Initially, a few automaton cells at the center of the growth-permitting region are designated as pro-liferative cells.Then the initial tumor is allowed to grow according to our CA model incorporating the additional interaction rules between the tumor and the suppression factors in the microenvironment.Certain geometrical characteristics of the tumor (e.g., tumor area, areas of different tumor cell populations) and its morphology (e.g., the geometrical positions of the tumor cells) are collected every T c days.We set γ = 0.005, α W = 0.75, α S = 0.15, k 0 = 0.8 and use these parameter values throughout this paper, except where otherwise stated.
in a confined space with ρECM = 0.30.As shown in Figure2(a), with the interactions between the tumor and the microenvironmental suppression factors incorporated, there exists a period of dormancy in the tumor's growth.Specifically, for the initial approximate 900 days, the tumor stays in a dormant state.Suddenly at approximately day 900, the tumor switches its behavior and begins rapid proliferation.The virtual patient would die 100 days after this critical point in time.Figure 2(b) shows the areas A of different populations normalized by the area of the growth-permitting area A 0 .For purposes of comparison, Figure 2(c) and Figure 2(d) show the statistics of the tumor growth without the suppression of microenvironmental factors.Moreover, by comparing Figure 2(a) and Fig-ure 2(c), it is seen that the interactions between the tumor and the microenvironmental suppression factors lead to the existence of a dormancy period and a subsequent emergent "switch" behavior of the tumor from a dormant state to a proliferative state.Also, from the comparison of Figure 2(b) and Figure 2(d), one can see that the additional interaction rules alter the fractions of necrotic cell population and proliferative cell population within the tumor.When suppression of the tumor growth is present, the necrotic region decreases and the proliferative region increases relatively; the area of the quiescent region remains almost unchanged.

Figure 3
Figure 3 shows snapshots of the simulated 2D tumor.It can be clearly seen that the tumor develops a highly aspherical morphology due to the interactions between the tumor and the microenvironmental factors.Figure 3 also demonstrates that the tumor hardly grows during the period of dormancy, but once the "switch" occurs, the tumor expands very rapidly.Henceforth, we will use the "CA dormancy model" to investigate the effects of the various parameters characterizing the tumor-host interactions on the growth dynamics of the tumor.These parameters include the rate of phenotypic transformation, by which the tumor cells gain the ability to proliferate against those suppression factors, the suppression rate imposed by suppression factors on individual tumor cells, and the mechanical rigidity of the microenvironment.

Figure 2 .
Figure 2. (Color online) Upper panel: statistics of a simulated noninvasive tumor growing in the ECM with ρECM = 0.3 and microenvironmental suppression factors, as predicted by the "CA dormancy model".(a) Tumor area A T normalized by the area A 0 of the growth permitting region.(b) Areas of different cell populations normalized by the area A 0 of the growth permitting region.Lower panel: statistics of a simulated noninvasive tumor growing in the ECM with ρECM = 0.3 without suppression.(c) Tumor area A T normalized by the area A 0 of the growth permitting region.(d) Areas of different cell populations normalized by the area A 0 of the growth permitting region.

Figure 3 .
Figure 3. (Color online) Snapshots of a simulated noninvasive tumor growing in the ECM with ρECM = 0.3 on different days given by the CA dormancy model.Upper panel: Dormancy period.Lower panel: Regrowth period.

Figure 4 .
Figure 4. (Color online) (a) The "critical" point at which the noninvasive tumor growing in the ECM with ρECM = 0.3 switches from a dormant state to a proliferative state as functions of α and γ.(b) A schematic phase diagram that characterizes the growth dynamics of a noninvasive tumor growing in the ECM with ρECM = 0.3 under different α and γ.

[ 51 ]
. Here we only consider the effects of the ECM macromolecule density, ECM degradation by the proliferative cells, and the pressure built up due to the ECM deformation by tumor growth.As shown in Figure 5, when ECM rigidity increases, the time at which the switch occurs gets delayed significantly and the final size of the tumor when it plateaus appreciably decreases.For example, with all the other parameters fixed, when the tumor grows in a soft ECM with ρECM = 0.15, the "switch" point occurs approximately on day 250.This is to be contrasted with growth in a rigid ECM with ρECM = 0.45 where the dormancy period could last for 2,000 days before a "switch" to a proliferative state occurs.Also, the plateau size of the tumor growing in a soft ECM with ρECM = 0.15 is five times as large as that of one growing in a rigid ECM with ρECM = 0.45.In other words, when the tumor grows in a harsher microenvironment, it's harder for the tumor to break out of a dormant state and potential proliferative growth is largely suppressed.Note that the rigidity of the microenvironment could also affect tumor growth via various intracellular signaling processes [51].Those mechanotransduction effects will be incorporated into our CA dormancy model in future work, which could result in different scenarios

Figure 5 .
Figure 5. (Color online) Simulated tumor area A T normalized by the area A 0 of the growth permitting region of a noninvasive tumor growing in the ECM with different ρECM .

Figure 6 .
Figure 6.(Color online) Tumor area A T normalized by the area A 0 of the growth permitting region of a simulated noninvasive tumor growing in the ECM under different killing rates by microenvironmental suppression factors.The parameter k 0 is the fraction that the suppression factors from the microenvironment kill the actively dividing proliferative cells when the suppression factors counteract these cells.
Also, for a tumor of a specific type, we can extract the parameter values in our model by fitting our simulation results to the statistics of a real in-vitro or in-vivo tumor of this type.Then we can utilize our model to explore optimal treatment strategies for the tumors of this specific type.In addition, once we determine the effects of a specific microenvironmental factor (e.g., specific integrins[59]) on the parameter values in our model, we could then study the effects of this microenvironmental factor on the tumor growth dynamics.Our current CA dormancy model is still preliminary, and to achieve our ultimate goal of understanding cancer dormancy and progression, we need to develop robust models that incorporate appropriate cell-level tumor-host interactions that are informed by experiments.For example, by explicitly considering angiogenesis and using more realistic distribution of "transformed" tumor cells' resistance to microenvironmental suppression factors (currently we just divide the "transformed" cells into two types with respect to their resistance to microenvironmental suppression factors: weak and strong), our model might be able to yield more realistic results and improve our understanding of cancer dormancy and progression.Also, the effects of tumor cell competition, cooperation and the microenvironmental changes caused by tumor cell activities could be incorporated to further strengthen our CA dormancy model[60].It is noteworthy that although we employ interaction rules based on a discrete cell model to describe "competition" between the tumor and the microenvironmental suppression factors, alternatives such as evolutionary game theory implemented by partial-differential equations are also available to address the interplay between the tumor and the microenvironment[61].