 Research
 Open access
 Published:
Mathematical modeling of solid cancer growth with angiogenesis
Theoretical Biology and Medical Modelling volume 9, Article number: 2 (2012)
Abstract
Background
Cancer arises when within a single cell multiple malfunctions of control systems occur, which are, broadly, the system that promote cell growth and the system that protect against erratic growth. Additional systems within the cell must be corrupted so that a cancer cell, to form a mass of any real size, produces substances that promote the growth of new blood vessels. Multiple mutations are required before a normal cell can become a cancer cell by corruption of multiple growthpromoting systems.
Methods
We develop a simple mathematical model to describe the solid cancer growth dynamics inducing angiogenesis in the absence of cancer controlling mechanisms.
Results
The initial conditions supplied to the dynamical system consist of a perturbation in form of pulse: The origin of cancer cells from normal cells of an organ of human body. Thresholds of interacting parameters were obtained from the steady states analysis. The existence of two equilibrium points determine the strong dependency of dynamical trajectories on the initial conditions. The thresholds can be used to control cancer.
Conclusions
Cancer can be settled in an organ if the following combination matches: better fitness of cancer cells, decrease in the efficiency of the repairing systems, increase in the capacity of sprouting from existing vascularization, and higher capacity of mounting up new vascularization. However, we show that cancer is rarely induced in organs (or tissues) displaying an efficient (numerically and functionally) reparative or regenerative mechanism.
Background
A total of 562, 875 cancer deaths were recorded in the United States in 2007, and it is estimated that approximately 570, 000 died from cancer in 2010. The overall estimate of approximately 1.53 million new cases does not include carcinoma in situ of any site except urinary bladder, nor does it include basal cell and squamous cell cancers of the skin. Greater than 2 million unreported cases of basal cell and squamous cell skin cancer, approximately 54, 010 cases of breast carcinoma in situ, and 46, 770 cases of melanoma in situ are expected to be newly diagnosed in 2010 [1]. In the world, the estimate for 1990 suggested a total of 8.1 million new cases, divided almost exactly between developed and developing countries, and 5.2 million cancer deaths, about 55% of which occurred in developing countries [2].
Cells of some organs, as the heart, stop proliferation after reaching their size, but others, like skin cells and cells that line body cavities, must proliferate almost continuously to replenish cells that are lost. Cancer arises when within a single cell multiple malfunctions of control systems occur, which are, broadly, the system that promote cell growth and the system that protect against erratic growth (for instance, tumor suppressor gene p53, which name comes from protein of molecular weight 53,000, and RB, from retinoblastoma). Additional systems within the cell must be corrupted so that a cancer cell, to form a mass of any real size, produces substances (such as VEGF  Vascular Endothelial Growth Factor) that promote the growth of new blood vessels. Cellular growthcontrol systems can be corrupted either when cellular genes are mutated or when proteins produced during a viral infection interfere with control system function. Multiple mutations are required before a normal cell can become a cancer cell by corruption of multiple (about five) growthpromoting systems [3].
With respect to the control systems that protect against cancer, they are classified in two general types: systems that prevent mutations, and systems that deal with mutations once they occurred. For instance, there are mechanisms in cells that can convert toxic byproducts of cellular metabolism or carcinogenic substances from the environment into harmless chemicals, preventing DNA damage. In addition, cells have a number of different systems that can detect damaged DNA and fix it, avoiding mutations. However, if the DNA of a cell has been damaged that repair would be impossible, there are systems which trigger the cell to die. Hence, even though there are trillions of cells, only one in three humans will get cancer during his lifetime, and cancer is mainly a disease of the elderly [3].
In order to develop effective treatments, it is important to identify the mechanisms to controlling cancer growth, how they interact, and how they can most easily be manipulated to eradicate (or manage) the disease. Aiming to gain such insight, it is usually necessary to perform large numbers of timeconsuming and intricate experiments  but not always. Through the development and solution of mathematical models that describe different aspects of solid tumor growth, applied mathematics has the potential to prevent excessive experimentation and also to provide biologists with complementary and valuable insight into the mechanisms that may control the development of solid tumors [4]. In [5] it is provided a concise history of the study of tumor growth, discussing some of the most influential mathematical models and their relationship to experimental studies, and illustrating how the field of cancer research has evolved due to these interactions between theoretical and experimental approaches.
In this paper we develop a mathematical model to understand the solid cancer growth inducing angiogenesis. We determine the steady states of the model and the stability of the equilibrium points are performed. We simulate the model to assess the appearance of cancer. Treatments aiming the control of cancer are not dealt here (see [6] and [7] for tumor control).
Methods
Our objective is the development of a simple mathematical model to describe the cancer growth dynamics inducing angiogenesis. The model does not include important physiological processes like the diffusion of oxygen into a solid where it is consumed by metabolic processes, the outward diffusion of lactic acid from a solid which produces it by metabolic processes and the diffusion of oxygen away from a blood vessel into a region with an oxygen debt. First, we describe an organ of human body free of cancer and, then, cancer cells appear in this organ.
The normal cells of an organ (the concentration at time t is designed as C) are feed by an existing vasculature formed by endothelial cells (the concentration at time t is designed as E). Adult endothelial cells and cells of an organ of human body are normally quiescent apart from certain developmental processes (e.g., embryogenesis), and proliferation of these cells aims to replenish losses (for instance, cells dying due to aging and wound). The normal cells and surrounding vascular networks (epithelial cells) are governed by the logistic growth, with intrinsic growth rates α_{1} and α_{2}, respectively, and they are under mortality rates μ_{1} and μ_{2}, respectively. Normal cells produce substances to promote the growth of blood vessels in a suitable network to attend their needs. The extension of vascularization (blood vessels network) depends on the size and function of the organ, and conversely, the blood vessels network determines the activity of the organ. In other words, the size of an organ could be determined by the surrounding network of blood vessels, that is, k_{1}(E), as k_{2} depends on C, k_{2}(C), where k_{1} and k_{2} are the carrying capacities of the organ (normal cells) and the extension of blood vessels network (epithelial cells), respectively. In a cancer free state, normal cells (C) and blood vessels (E) are in a steady state, where the amount of normal cells depends on the feeding network, as well as the amount of epithelial cells depends on the size and function of organs that must feed.
Let us suppose that a series of accumulating mutations in normal cells (one single cell may be sufficient) corrupted systems that promote cell growth and that protect against erratic growth. Also systems that control the production of substances which induce new blood vessels are corrupted in the cancer cells. These sequential events initiate a rapid growth of cancer due to angiogenesis. Angiogenesis is the process by which new blood vessels develop from an existing vasculature, through endothelial cell sprouting, proliferation, and fusion [8]. The angiogenesis links the relatively harmless avascular and the potentially fatal vascular growth phases of the tumor [9].
Let us consider what occurs in cancer growth: (1) there is an initial amount of cancer cells (the model does not deal with the origin of cancer cells) appearing in a completely health organ; (2) these cells induce the formation of prevascular cells from existing vascular network; (3) after a period of time, the formation of new blood vessels is initiated; and (4) both normal and cancer cells compete for nutrients and space (proteins, oxygen, etc.).
For cancer cells (the concentration at time t is designed as T), the growth, which is likely to be limited by energetic constraints, is limited by the carrying capacity k_{3} and the surrounding new vessels originated by angiogenesis (the concentration at time t is designed as A). The new vessels, represented by angiogenesis cells A, have growth depending on tumor cells T and are limited by available space (k_{4}). The growth of blood vessels (E and A) are constrained by physical restrictions, even they are never subject to either oxygen or nutrient deprivation. We assume the existence of an early stage of angiogenesis, called preangiogenesis (the concentration of cells at time t is designed as P). This stage corresponds to the release by tumoral cells of VEGF, which starts to diffuse into the surrounding tissue and approach the endothelial cells of nearby blood vessels. Endothelial cells subsequently respond to the VEGF concentration gradient by forming sprouts. From these sprouts, there occur migration and proliferation toward the tumor, and these new vessels are called as angiogenesis (endothelial cells A).
The tumoral cells induce new vascular network from the existing one to feed cancer cells. The preangiogenesis cells are formed by the mass action law [10] at a constant rate γ, the rate of sprouts formation, from which new vascularizations occur resulting in angiogenesis cells after an average period of time δ^{1}, where δ is transfer rate from preangiogenesis to angiogenesis cells. The cancer, preangiogenesis and angiogenesis cells are under the mortality rates μ_{3}, μ_{4} and μ_{5}, respectively. The cancer and angiogenesis cells are governed by the logistic growth: intrinsic growth α_{3} and ε, which depend respectively on angiogenesis and cancer cells, and carrying capacities k_{3} and k_{4}, respectively. Cancer cells can also grow logistically using existing vascularization at rate {\alpha}_{3}^{\prime}. The rates at which normal and cancer cells compete themselves for resources and space are {\beta}_{1}^{\prime} and {\beta}_{2}^{\prime}. Both parameters also take into account effects of the environment.
Preangiogenesis cells deserve some words. New blood vessels to nourish cancer cells must sprout from preexisting blood vessels. However, there is an elapse of time from the appearance of sprouts until the complete formation of new blood vessels that effectively feed cancer cells. By doing this, we take into account a timedelay between the formation of new network of blood vessels (angiogenesis cells A) and tumor growth. Notice that the probability of sproutings becoming angiogenesis cells is δ/(δ + μ_{4}). If we let δ → ∞, this phase is negligible (probability one).
Let us consider the normal cells C taking into account the above descriptions. In the absence of cancer cells T, normal cells obey \frac{d}{dt}C={\alpha}_{1}C\left[1\frac{C}{{k}_{1}\left(E\right)}\right]{\mu}_{1}C. The cancer cells influence negatively normal cells, and this term is described by \frac{{\alpha}_{1}C{{\beta}^{\prime}}_{1}T}{{k}_{1}\left(E\right)} with negative signal. Due to mutation, some normal cells become cancer cells, which can occur continuously. However, we will assume that the mutations occur instantaneously at a given time (cancer cells grow exponentially in the beginning, hence further mutations, being in small amount, can be disregarded), when an amount of normal cells are transferred to cancer cells compartment. For other types of cells, the dynamical equations can be obtained similarly.
Based on the above definitions of variables and model's parameters, the dynamics of the cancer growth is described by the following system of equations
In this model, an amount C_{ m }of normal cells are mutated to cancer cells at time t = 0, which is described by the Dirac delta function δ_{ D }(t): it assumes ∞ at t = 0, and 0, otherwise.
Let us simplify the model assuming that k_{1} and k_{2} are constant, depending only on the status corresponding to cancer free. Under this simplification, we can define the interaction parameters as {\beta}_{1}={\beta}_{1}^{\prime}{\alpha}_{1}/{k}_{1} and {\beta}_{2}={\beta}_{2}^{\prime}{\alpha}_{3}^{\prime}/{k}_{1}. Another simplification assumes that the rate at which cancer cells increase due to new vascularization is much higher than due to existing one, or we let {\alpha}_{3}AT\left(1T/{k}_{3}\right)\gg {\alpha}_{3}^{\prime}T\left(1T/{k}_{1}\right). The above impulsive system taking into account these simplifications can be written as, for t > 0,
in order to describe a perturbation introduced at t = 0 in a cancer free status. Hence, the initial conditions supplied to the system of equations become:
where C_{0} and E_{0} are the cancer free steady state values (see below), and C_{ m }is the instantaneous mutation of normal cells to initiate cancer growth. The vital dynamics of normal cells is similar than that presented in [11]. In Table 1 we present the summary of variables, which depend on time, for instance C(t), and parameters, which do not depend on time, of the model.
We next obtain the steady states of the system of equations (1), and the corresponding characteristic equations in order to establish the stability of the steady states.
Steady states
We determine the steady states of the system (1), denoted by \stackrel{\u0304}{Q}=\left(\stackrel{\u0304}{C},\u0112,\stackrel{\u0304}{T},\stackrel{\u0304}{P},\u0100\right). The coordinates of the equilibrium point are obtained by letting the time derivatives in equation (1) equal to zero, for instance, \frac{d}{dt}C=0. The first is the equilibrium where there is not any cells, that is,
For obvious reason, this point is not considered here.
The trivial equilibrium point, given by
where C_{0} and E_{0} are
corresponds to the absence of cancer. The trivial equilibrium point exists for α_{1} > μ_{1} and α_{2} > μ_{2}. In this model the normal and epithelial cells are not dependent on each other. One possibility of dependence can be done by the carrying capacities k_{1}(E_{0}) and k_{2}(C_{0}).
The coordinates of the nontrivial equilibrium point {\stackrel{\u0304}{Q}}^{*} are
where \stackrel{\u0304}{T} is the positive solution of the equation
with the fourth degree polynomial f(T), and the third degree polynomials g_{1}(T) and g_{2}(T) being given by
By inspecting \stackrel{\u0304}{C}, Ē and Ā, the positive solution \stackrel{\u0304}{T}, in order to {\stackrel{\u0304}{Q}}^{*} be biologically feasible, must satisfy, respectively, the constraints \stackrel{\u0304}{T}<{T}_{C},\stackrel{\u0304}{T}<{T}_{E} and \stackrel{\u0304}{T}<{T}_{A}, where
However, the constraint for the equilibrium Ā is \stackrel{\u0304}{T}<{T}_{g}, where
but it is always satisfied when \stackrel{\u0304}{T}<{T}_{C}, since T_{ g }> T_{ C }.
With respect to the constraints, we must have \stackrel{\u0304}{T}<\text{min}\left\{{T}_{C},{T}_{E},{T}_{A}\right\}, where min{T_{ C }, T_{ E }, T_{ A }} is the minimum among T_{ C }, T_{ E }and T_{ A }. The constraints depend inversely with removing parameters β_{1} and γ, and directly with the carrying capacity k_{3}. Hence, for higher values of removing parameters β_{1} and γ, and lower values of carrying capacity k_{3}, the constraints are decreased, that is, the steady state cancer cells \stackrel{\u0304}{T} must assume lower values. Let us suppose that k_{3} is high, which implies that cancer cells \stackrel{\u0304}{T} can achieve higher values. If the sproutings rate γ from existing blood vessels is high, in principle it seems to be beneficial to cancer cells. However, higher values of γ result in lower values for the constraint T_{ E }, and, consequently, \stackrel{\u0304}{T} must assume lower values, because it must be lower than both T_{ E }and k_{3}. Therefore, cancer cells will achieve higher values for high carrying capacity k_{3}, and lower sprouting rate γ. Further angiogenesis cells must be originated by the clonal division from a small amount of preangiogenesis cells. In appendix A we show the analysis of equation (5). Summarizing, the number of solutions \stackrel{\u0304}{T} depends on the values assigned to the model's parameters, which can be 0, 2 or 4. Notwithstanding, the solution \stackrel{\u0304}{T} must obey the constraints given in equation (7), which reduces the number of solutions up to 2. In the case of 2 solutions, we denote the solutions \stackrel{\u0304}{T} as small (T_{ < }, and the corresponding equilibrium point {\stackrel{\u0304}{Q}}_{<}^{*}) and big (T_{>}, and the corresponding equilibrium point {\stackrel{\u0304}{Q}}_{>}^{*}). Hence, the constraints play an important role in the number of equilibrium points. Notice that the interacting parameter ε does not appear in the constraints. Let us analyze possible scenarios resulting by varying the nonlinear term parameters β_{1}, β_{2}, γ and ε. First, let us vary β_{1}. Notice that T_{ C }depends on β_{1}, but T_{ E }and T_{ A }do not. For small β_{1}, the solution \stackrel{\u0304}{T} of equation (5) situates below T_{ C }. However, as β_{1} increases, \stackrel{\u0304}{T} can surpass T_{ C }, after intercepting T_{ C }at {\beta}_{1}={\beta}_{1}^{c}. As \stackrel{\u0304}{T} increases with β_{1}, \stackrel{\u0304}{C}, given by the first equation of (4), decreases and reaches zero at {\beta}_{1}={\beta}_{1}^{c}, and \stackrel{\u0304}{C}<0 sinceafter. At {\beta}_{1}={\beta}_{1}^{c} we have \stackrel{\u0304}{T}={T}_{C}^{c}=\left({\alpha}_{1}{\mu}_{1}\right)/{\beta}_{1}^{c}. Let us suppose that {T}_{C}^{c}<\text{min}\left\{{T}_{E},{T}_{A}\right\}, the minimum between T_{ E }and T_{ A }. When \stackrel{\u0304}{C}=0, there arises another feasible equilibrium point {\stackrel{\u0304}{Q}}^{c}, with the coordinates
where \stackrel{\u0304}{T}={T}_{C}^{c}, and the equilibrium \stackrel{\u0304}{T} does not depend anymore on the interacting parameter β_{1}. Hence for {\beta}_{1}\ge {\beta}_{1}^{c} we have a constant value \stackrel{\u0304}{T}=\left({\alpha}_{1}{\mu}_{1}\right)/{\beta}_{1}^{c}, which implies in constant values for Ē, \stackrel{\u0304}{P} and Ā.
See more discussion below, in the stability analysis.
The equilibrium point {\stackrel{\u0304}{Q}}^{c} corresponds to the case where the cancer cells displaced normal cells. By inspecting the equation for \stackrel{\u0304}{C}, the first in equation (4), the normal cells drop out to zero for a sufficiently higher values of β_{1}. In another words, when cancer cells have very higher fitness than normal cells and overcome the competition for resources, they can lead to the exclusion of the normal cells (of course, the cancer diseased person dies before reaching this equilibrium). Notice that {\stackrel{\u0304}{Q}}^{c} can be avoided as the equilibrium point when k_{3} assumes small values, that is, if k_{3} < min{T_{ E }, T_{ C }}, minimum between T_{ E }and T_{ C }. In this situation, the equilibrium value \stackrel{\u0304}{T} is always smaller than k_{3}, and all coordinates of the equilibrium point are positive (Appendix A).
Second, let us now vary β_{2}. The third equation of (1) in the equilibrium can be written as
When β_{2} increases, \stackrel{\u0304}{T} decreases, and at a certain value of β_{2}, say {\beta}_{2}^{th},\stackrel{\u0304}{T} assumes zero value. When \stackrel{\u0304}{T}=0, the unique feasible equilibrium point is the trivial {\stackrel{\u0304}{Q}}_{0}. When the influence of the normal cells and the environment is higher ({\beta}_{2}\ge {\beta}_{2}^{th}), then cancer can not be maintained.
Cancer onset depends on the interplay among the initial condition T(0) = C_{ m }and interacting parameters β_{1} and β_{2}. Once cancer cells are originated by mutation from several sources, these cells interact with healthy cells and the surrounding environment. While the parameter β_{1} is practically insensitive in the behavior of the dynamics, however, there is a threshold in the parameter β_{2}, above which cancer can not be established. This result shows that cancer is rarely induced in organs (or tissues) displaying an efficient (numerically and functionally) reparative or regenerative mechanism [12], which is the reason for the incidence of cancer rising exponentially with age [13].
Finally, let us analyze the effects of varying the sprouting rate γ and angiogenesis increasing rate ε.
When γ increases, Ē decreases according to the second equation of (4). Let us assume that Ē reaches zero at a certain value named γ^{c}, and Ē < 0 sinceafter. When Ē = 0, differently from \stackrel{\u0304}{C}=0, the unique feasible equilibrium point is the trivial {\stackrel{\u0304}{Q}}_{0}. Hence, we expect that Ē > 0 for γ ≥ 0, because of \stackrel{\u0304}{T}<{T}_{E}=\left({\alpha}_{2}{\mu}_{2}\right)/\gamma is always obeyed, and there is not a critical value γ^{c}. When we decrease γ, we observe that equation (5) may not have biologically feasible solutions.
Initially, let us analyze γ = 0. In this case, we have f(T) = 0, and the equilibrium value \stackrel{\u0304}{T} is solution of the polynomial of degree two
with another being given by \stackrel{\u0304}{T}={T}_{g}. Moreover, from equation (4), we have \stackrel{\u0304}{P}=0 (due to γ = 0, the role of the parameter δ does not matter). In this case, there arises another feasible equilibrium point {\stackrel{\u0304}{Q}}^{p}, with the coordinates
Equation (10) has two positive solutions (not shown) when T_{ g }> k_{3} and ε > μ_{5}/k_{3}, where T_{ g }is given by equation (8), while a unique positive solution occurs when T_{ g }= k_{3}. Let us introduce a threshold of ε, named ε^{th}. We remark that this threshold value appears only for γ = 0. Hence, when ε > ε^{th}, with ε^{th}> μ_{5}/k_{3}, we have two positive solutions, named {T}_{<}^{\epsilon} (and the corresponding equilibrium point {\stackrel{\u0304}{Q}}_{<}^{p}) and {T}_{>}^{\epsilon} (and the corresponding equilibrium point {\stackrel{\u0304}{Q}}_{>}^{p}), which collapse to one at ε = ε^{th}, and for ε < ε^{th}there is not positive solution. The equilibrium point {\stackrel{\u0304}{Q}}^{p} represents the origin of angiogenesis from other sources, not from existing vascular system feeding an organ of the human body. Two facts must occur to appearing of {\stackrel{\u0304}{Q}}^{p}: sufficiently higher capacity of cancer cells to originate new vascular system from external sources (ε > ε^{th}), and higher initial amount, at time 0, of vascular cells \left(T\left(0\right)>{T}_{<}^{\epsilon}\right).
Now, for lower values of γ (γ ≳ 0), as a consequence of the appearing of ε^{th}for γ = 0, equation (5) does not have positive solution and there is not biologically feasible equilibrium point. Let us introduce a threshold of γ, named γ^{th}. For γ < γ^{th}the unique solution is \stackrel{\u0304}{T}=0 (this fact is shown numerically). Note that as ε increases, γ^{th}decreases. Hence, higher the capacity of building up new vascularization ε, less the amount of sprouts originated from existing blood vessels needed to angiogenesis. In the special case ε = 0 (cancer cells do not promote growth in the new vascularization), the equilibrium \stackrel{\u0304}{T}, which must be substituted in the variables given in equation (4), must be solution of f(T) = g(T), where
This equation clearly shows that a positive solution \stackrel{\u0304}{T} is feasible for sufficiently higher values of γ, that is, \gamma >{\gamma}_{0}^{th}, where {\gamma}_{0}^{th} corresponds to ε = 0. Due to the fact that ε increases as γ^{th}decreases, {\gamma}_{0}^{th} is the upper bound of the threshold of γ. Even the cancer cells do not promote the proliferation of the sprouted endothelial (preangiogenesis) cells (ε = 0), if the capacity of endothelial cell sprouting is higher (γ > γ^{th}), then cancer cells can be maintained.
There are two mechanisms by which the tumor's vasculature is built: (1) as a tumor grows, it excretes tumor angiogenesis factor, which help activate endothelial cells of nearby blood vessels, initiating angiogenesis; (2) bone marrow derived endothelial progenitor cells are mobilized into the blood, by means of longrange signaling, and they are recruited from blood by the tumor by shortrange signaling, resulting in vasculogenesis (the vascular endothelium of the nearby vessels become activated and allows the endothelial progenitor cells to extravase and start a cycle of differentiation/division) [14]. In both cases, the newly stimulated endothelial cells (described by parameter γ) and the recruited endothelial progenitor cells (this phenomenon is not considered in the model; however, the case γ = 0 can in some extent be understood as vasculogenesis) enter a stage of clonal expansion and continue to form blood vessels, which is described by the parameter ε.
With respect to the linear term of the preangiogenesis parameter δ, when \delta \to \infty ,{\stackrel{\u0304}{Q}}^{\infty} appears, similar to the equilibrium point {\stackrel{\u0304}{Q}}^{p}, because \underset{\delta \to \infty}{\text{lim}}\stackrel{\u0304}{P}=0. However, from the third equation of (4), we have \underset{\delta \to \infty}{\text{lim}}\delta \stackrel{\u0304}{P}=\gamma \u0112\stackrel{\u0304}{T}. Hence, the coordinates of {\stackrel{\u0304}{Q}}^{\infty} are those given in equation (11), but \stackrel{\u0304}{T} is solution of the equation f(T) = g_{1}(T) + g_{2}(T), where
and g_{1}(T) + g_{2}(T) is given be equation (6). Notice that f(T) differs from that in equation (6) by the factor δ/(μ_{4} + δ). The equilibrium point {\stackrel{\u0304}{Q}}^{\infty} represents the absence of the intermediate preangiogenesis, that is, angiogenesis occurs directly and instantaneously from epithelial cells.
Summarizing, the nontrivial equilibrium points {\stackrel{\u0304}{Q}}^{*} ({\stackrel{\u0304}{Q}}_{<}^{*} and {\stackrel{\u0304}{Q}}_{>}^{*}, when they exist), {\stackrel{\u0304}{Q}}^{c} (for {\beta}_{1}\ge {\beta}_{1}^{c}), and {\stackrel{\u0304}{Q}}^{p} (when γ = 0) are biologically feasible, but {\stackrel{\u0304}{Q}}^{c} is not a real cancer equilibrium point, since \stackrel{\u0304}{C}=0. Especially, {\stackrel{\u0304}{Q}}^{p} is a cancer equilibrium point only if vasculogenesis can act alone to generate new blood vessels. The analysis of the model showed the existence of thresholds {\beta}_{2}^{th} and γ^{th}(and, eventually, ε^{th}, for γ = 0), and a critical value {\beta}_{1}^{c}.
We analyzed the steady states of the model described by the system of ordinary differential equations (1) assuming that all parameters of the linear terms (μ_{1}, μ_{2}, μ_{3}, μ_{4}, μ_{5} and δ) plus the intrinsic growth rates α_{1}, α_{2} and α_{3} and carrying capacities k_{1}, k_{2}, k_{3} and k_{4} are fixed values. But, the parameters regarded to the interaction between different cells (β_{1}, β_{2}, γ and ε) are generally unknown, hence they were allowed to vary to analyze broad range of variation in these parameters.
Local stability analysis
The local stability of the equilibrium points [15] is determined by the eigenvalues of the characteristic equation det \left(\stackrel{\u0304}{J}\lambda I\right)=0, where \stackrel{\u0304}{J} is the Jacobian J evaluated at the equilibrium point under analysis,
where j_{1}, j_{2}, j_{3} and j_{4} are
The Jacobian \stackrel{\u0304}{J} was obtained using equilibrium equations given in (4) assuming that \stackrel{\u0304}{C}\ne 0. When \stackrel{\u0304}{C}=0, the element in the first line and first column is {a}_{11}=\left({\alpha}_{1}{\mu}_{1}\right){\beta}_{1}\stackrel{\u0304}{T}, and in the third column, a_{13} = 0.
We present analytical results corresponding to the equilibrium points {\stackrel{\u0304}{Q}}^{0} and {\stackrel{\u0304}{Q}}^{c}.
The eigenvalues (from Jacobian {\stackrel{\u0304}{J}}^{0}) corresponding to the trivial equilibrium point {\stackrel{\u0304}{Q}}^{0} are: λ_{1} = μ_{5}, λ_{2} =  (μ_{4} + δ), λ_{3} =  (μ_{3} + β_{2}C_{0}), λ_{4} =  (α_{1}  μ_{1}), and λ_{5} =  (α_{2}  μ_{2}). In the last two eigenvalues we used the coordinates given by equation (3). Hence, the trivial equilibrium {\stackrel{\u0304}{Q}}^{0} is locally asymptotically stable whenever this point exists (α_{1} > μ_{1} and α_{2} > μ_{2}), because all eigenvalues become negative.
The nontrivial equilibrium point {\stackrel{\u0304}{Q}}^{c}, with coordinates given in equation (9), has one of them zero (\stackrel{\u0304}{C}=0). The eigenvalues corresponding to the trivial equilibrium point {\stackrel{\u0304}{Q}}^{c} are obtained from the Jacobian {\stackrel{\u0304}{J}}^{c}. Once \stackrel{\u0304}{C}=0, one of the eigenvalues is {\lambda}_{1}={\beta}_{1}\left(\stackrel{\u0304}{T}{T}_{C}\right), arising one of the conditions to the not real cancer equilibrium {\stackrel{\u0304}{Q}}^{c} being stable, that is, \stackrel{\u0304}{T}>{T}_{C}, where T_{ C }is given by equation (7). The remaining four eigenvalues are obtained from the submatrix {\stackrel{\u0304}{J}}_{1}^{c}, which comes from {\stackrel{\u0304}{J}}^{c} excluding first row and first column. Let us define matrix A as A={\stackrel{\u0304}{J}}_{1}^{c}, or,
where j_{1}, j_{2}, j_{3} and j_{4} were defined in equation (13). In appendix B we show that, for big solution T_{>}, A is an Mmatrix, and, hence, all eigenvalues associated to {\stackrel{\u0304}{J}}_{1}^{c} have negative real part. Hence, {\stackrel{\u0304}{Q}}^{c}, which is biologically feasible but not real cancer, is locally asymptotically stable only if \stackrel{\u0304}{T}>{T}_{C}. Therefore, {\stackrel{\u0304}{Q}}^{c} appears when one of the constraints to {\stackrel{\u0304}{Q}}^{*} be biologically feasible is violated: in the equation (4), instead of \stackrel{\u0304}{C}<0, we have \stackrel{\u0304}{C}=0 for {\beta}_{1}\ge {\beta}_{1}^{c}.
The local stability of the nontrivial equilibrium {\stackrel{\u0304}{Q}}^{*} is performed numerically. However, we stress the fact that the equilibrium {\stackrel{\u0304}{Q}}^{c} is an extension of {\stackrel{\u0304}{Q}}^{*} for {\beta}_{1}\ge {\beta}_{1}^{c}. Hence, the big solution T_{>} of equation (5) for {\beta}_{1}\ge {\beta}_{1}^{c} must be stable. For this reason we showed the local stability of {\stackrel{\u0304}{Q}}^{c}, even though this is not a real cancer equilibrium.
Results and Discussion
In this section we present numerical simulations of the model, which are discussed. The numerical methods used are bisection (to find zeros of polynomials) and 4^{th}order RungeKutta (to solve system of ordinary differential equations) [16].
Numerical analysis of the model
Numerical simulations are performed taking into account the values and units of the model's parameters given in Table 2. These values have the purpose of illustrating the outcomes of the model. Cancer free equilibrium point {\stackrel{\u0304}{Q}}^{0} has the cancer free concentrations: C_{0} = 9 cells/uv and E_{0} = 10 cells/uv, where uv stands for an arbitrarily unit of volume. Hereafter we will omit units of all variables and parameters. In this section we deal with the steady state of the system (1), determining the equilibrium points and bifurcation diagrams. We also study the dynamical trajectories of the system (1) assuming that mutations occurred and cancer cells arise suddenly. For this reason the initial conditions are given by equation (2). The reason behind this is the fact that the trivial equilibrium point {\stackrel{\u0304}{Q}}^{0} always exists and is stable. Hence, the initial conditions supplied to the system of equations correspond to an appearance of C_{ m }number of cancer cells in a cancer free situation {\stackrel{\u0304}{Q}}^{0}, that is, at time t = 0 the variables assume C(0) = C_{0}  C_{ m }, E(0) = E_{0}, T(0) = C_{ m }, P(0) = 0 and A(0) = 0. The initial amount of cancer cells is originated by mutation from several sources. In some special situations we will use C(0) = C_{0}. The initial amount of mutations C_{ m }can either grow up to cancer, or fades out.
Based on the values given in Table 2, the constraints, using equations (7) and (8), are T_{ C }= 9.0 (at β_{1} = 0.01), T_{ E }= 5.0, T_{ A }= 5.0 and T_{ g }= 14.0. The coordinates of the nontrivial equilibrium points, corresponding to small (T_{ < }) and big (T_{ > }) equilibrium values of cancer cells obtained from equation (5), are, respectively, {\stackrel{\u0304}{Q}}_{<}^{*}=\left(8.92,9.85,0.08,0.07,0.71\right) and {\stackrel{\u0304}{Q}}_{>}^{*}=\left(5.32,2.64,3.68,0.88,1.96\right), from equation (4). The corresponding five eigenvalues (λ_{ i }, i = 1, ..., 5) are given in Table 3, showing that the small is unstable (λ_{1} > 0), while the big is stable (all λ_{ i }'s are negative). The big equilibrium T_{ > }= 3.68 satisfies the constraints given in equation (7). The trivial equilibrium {\stackrel{\u0304}{Q}}^{0} is always stable.
In Figure 1 we illustrate the dynamical trajectories of system (1) taking into account the values of parameters given in Table 2. The dynamical trajectories depend on the initial conditions. In Figure 1(a), due to T(0) = 1.023, the dynamical trajectories are attracted to the trivial equilibrium {\stackrel{\u0304}{Q}}^{0}; while for T(0) = 1.024, the dynamical trajectories attain the nontrivial equilibrium {\stackrel{\u0304}{Q}}_{>}^{*} (Figure 1(b), which shows a sudden relapse to cancer onset at around 400 days). The exact critical amount of cancer cells T* that divides two regions of attraction situates in the open interval (1.023,1.024). Notice that T* is much higher than T_{ < }= 0.08, which comes out due to other initial conditions. However, if the initial conditions are the coordinates of {\stackrel{\u0304}{Q}}_{<}^{*}, except T(0), that is, C(0) = 8.92, E(0) = 9.85, P(0) = 0.07 and A(0) = 0.71, then, for an arbitrary value ϵ > 0, the dynamical trajectories go to trivial {\stackrel{\u0304}{Q}}^{0} when T(0) = (1  ϵ) × 0.08, and to the nontrivial {\stackrel{\u0304}{Q}}_{>}^{*} when T(0) = (1 + ϵ) × 0.08 (for instance, ϵ = 0.001, figures not shown). In this case, we have T* = T_{ < }. Hence, {\stackrel{\u0304}{Q}}_{<}^{*} is the breakpoint, and its coordinates generate a hypersurface that divides two attracting regions (see Appendix A). Notice that T(0) = 1.023 which plays the role of separating two attracting regions is 12fold higher than T_{<} = 0.08.
The initial amount of normal cells that suffer mutation C_{ m }is crucial to trigger cancer disease. Figure 1 shows that if this amount is below a critical value, that is, T(0) = C_{ m }< T*, then repairing systems act efficiently and cancer does not settle in an organ. However, if the mutated cells surpass the critical value, the repairing systems do not avoid the onset of cancer.
Next, numerical simulations are performed in order to assess the effects of varying parameters β_{1}, β_{2}, γ and ε. In all dynamical trajectories, remember that unstable solution T_{ < }is very small, but T(0) is higher due to the initial conditions supplied to system (1) correspond to the coordinates of the trivial equilibrium {\stackrel{\u0304}{Q}}^{0}, and not those of the unstable {\stackrel{\u0304}{Q}}_{<}^{*}. Hence, T(0) is not comparable with T_{ < }, one of the coordinates of the breakpoint {\stackrel{\u0304}{Q}}_{<}^{*}. The cancer cells proliferate above the subclinical threshold of 10^{3} cells and reach 10^{9} cells which is the Xray detectable threshold [11].
Interaction between normal and cancer cells  β_{1} and β_{2}
Direct competition between normal and cancer cells for resources and space occur in order to grow. But there are many factors that affects both populations, like indirect interaction as excretion released by cells, changes in the environment and control mechanisms. The parameters β_{1} and β_{2} take into account these factors also.
Let us assess the cancer cells affecting negatively in the normal cells, by varying β_{1}.
In Figure 2 we show \stackrel{\u0304}{T}, the solutions of f(T)  g(T) = 0. Figure 2(a), which corresponds to the figure shown in Appendix A with T_{ E }= T_{ A }= 5.0, shows the existence of two positive solutions: small (denoted by T_{ < }) and big (denoted by T_{>}) for four values of β_{1}: 0.01 (labelled b_{1}), 0.1 (b_{2}), 1.0 (b_{3}) and 10.0 (b_{4}). In Figure 2(b) a zoom near lower values of β_{1} is shown to enhance the small solution T_{ < }.
Let us vary β_{1} and compute the corresponding equilibrium value \stackrel{\u0304}{T} using equation (5). As shown in Figure 2, we have two positive solutions: small, T_{ < }, and big, T_{>}. Substituting this solution into equation (4), we obtain the coordinates of the nontrivial equilibrium point {\stackrel{\u0304}{Q}}^{*}. In Figure 3 we show the coordinates of the equilibrium points {\stackrel{\u0304}{Q}}_{<}^{*} (a) and {\stackrel{\u0304}{Q}}_{>}^{*} (b). In (b) we also show the curves of T_{ C }and T_{ g }, which intercept the curve of big solution T_{>}. When T_{>} = T_{ C }, which occurs at {\beta}_{1}={\beta}_{1}^{c}=2.1273\times 1{0}^{2}, we have \stackrel{\u0304}{C}=0, at which the big nontrivial {\stackrel{\u0304}{Q}}_{>}^{*} disappears and arises an another equilibrium {\stackrel{\u0304}{Q}}_{>}^{c}, with \stackrel{\u0304}{C}=0 and other coordinates given by equation (9), which has fixed value {T}_{>}={T}_{C}^{c}=\left({\alpha}_{1}{\mu}_{1}\right)/{\beta}_{1}^{c}=4.2306 for {\beta}_{1}\ge {\beta}_{1}^{c}. Coordinates of {\stackrel{\u0304}{Q}}_{>}^{c} are same for all {\beta}_{1}\ge {\beta}_{1}^{c}. Mathematically, there is another value of β_{1}, which does not change the existing equilibrium point {\stackrel{\u0304}{Q}}^{c}, at which we have T_{>} = T_{ g }, that is, {\beta}_{1}={\beta}_{1}^{A}=2.8001\times 1{0}^{2}. At this value we have Ā = 0, with {T}_{>}={T}_{g}^{A}=\frac{{\alpha}_{1}{\mu}_{3}}{{\beta}_{1}^{A}{\beta}_{2}{k}_{1}}+\frac{{\alpha}_{1}{\mu}_{1}}{{\beta}_{1}^{A}}=4.9996.
Figure 3 shows that trivial {\stackrel{\u0304}{Q}}^{0} and small nontrivial {\stackrel{\u0304}{Q}}_{<}^{*} exist for all β_{1}, but big nontrivial is {\stackrel{\u0304}{Q}}_{>}^{*}\left({\beta}_{1}<{\beta}_{1}^{c}\right) or {\stackrel{\u0304}{Q}}_{>}^{c}\left({\beta}_{1}\ge {\beta}_{1}^{c}\right). In Figure 4, we show the bifurcation diagram, considering \stackrel{\u0304}{T} as a function of β_{1} (curve T in (a) and (b) of Figure 3). For all values of β_{1} the small equilibrium point {\stackrel{\u0304}{Q}}_{<}^{*} is the breakpoint, which divides two regions where trivial {\stackrel{\u0304}{Q}}^{0} or big nontrivial {\stackrel{\u0304}{Q}}_{>}^{*} (or {\stackrel{\u0304}{Q}}_{>}^{c}) is attracting point. Initial conditions set in a small region marked with I and Ia are attracted to the trivial equilibrium point {\stackrel{\u0304}{Q}}^{0}. However, initial conditions set in regions marked with II and III are attracted to big nontrivial equilibrium point {\stackrel{\u0304}{Q}}_{>}^{*} for {\beta}_{1}<{\beta}_{1}^{c}, while for {\beta}_{1}\ge {\beta}_{1}^{c} (regions marked with IIa and IIIa), to {\stackrel{\u0304}{Q}}_{>}^{c}. Notice that T_{<} always decreases very slowly with β_{1} (see also Figure 3(a)), showing that as β_{1} increases, less amount of initial cancer cells is needed to trigger a cancer. But this variation is quite insensitive.
In Figure 5 we show the dynamical trajectories, changing only β_{1} in 1000 times the values of parameters given in Table 2, or β_{1} = 10.0, being attracted to the special equilibrium point {\stackrel{\u0304}{Q}}^{c}, which is biologically feasible, but does not describe evolution of cancer. Dynamical trajectories are attracted to {\stackrel{\u0304}{Q}}^{0} (a) when T(0) = 0.151, while for T(0) = 0.152, to {\stackrel{\u0304}{Q}}_{>}^{c} (b). The dynamical trajectories in Figure 5 are similar to Figure 1, except C, which decreases abruptly to zero and situates practically in the vertical axis, and after 650 days increases quickly to \stackrel{\u0304}{C} (a), or remains in the horizontal axis (b). Increasing β_{1} in 100fold, T(0) decreased 7fold. For another β_{1} = 100.0, when T(0) = 0.142, trajectories go to {\stackrel{\u0304}{Q}}^{0}; while for T(0) = 0.143, to {\stackrel{\u0304}{Q}}_{>}^{c}, which coordinates have exactly the same values found in the previous case (β_{1} = 10.0). Both cases (figures not shown) correspond to {\beta}_{1}>{\beta}_{1}^{c}, hence {\stackrel{\u0304}{Q}}_{>}^{c} are the same, while {\stackrel{\u0304}{Q}}_{<}^{c} decrease little bit.
Let us assess the normal cells affecting negatively in the cancer cells, by varying β_{2}.
In Figure 6 we show the coordinates of the equilibrium points {\stackrel{\u0304}{Q}}_{<}^{*} (a) and {\stackrel{\u0304}{Q}}_{>}^{*} (b) by varying β_{2}. Positive solutions for equation (5) disappeared for higher β_{2}.
In Figure 7, we show the bifurcation diagram, considering \stackrel{\u0304}{T} as a function of β_{2} (curve T in (a) and (b) of Figure 6). For {\beta}_{2}<{\beta}_{2}^{th}, initial conditions set in region marked with I are attracted to the trivial equilibrium point {\stackrel{\u0304}{Q}}^{0}, while for those set in regions II and III are attracted to big nontrivial equilibrium point {\stackrel{\u0304}{Q}}_{>}^{*}. However, for {\beta}_{2}>{\beta}_{2}^{th} all initial conditions are attracted to {\stackrel{\u0304}{Q}}^{0}, which is the unique equilibrium, because there is not any positive solution for equation (5). Hence, there is a threshold of the parameter β_{2}, denoted by {\beta}_{2}^{th}, above which all trajectories go to trivial equilibrium. At the threshold value {\beta}_{2}^{th}=4.8376\times 1{0}^{2} both roots assume same value, that is, {T}_{<}\left({\beta}_{2}^{th}\right)={T}_{>}\left({\beta}_{2}^{th}\right)=0.7053.
For {\beta}_{2}<{\beta}_{2}^{th} the dynamical trajectories are similar to that shown in Figure 1. For instance, when β_{2} = 0.04, for T(0) = 10.60 trajectories are attracted to {\stackrel{\u0304}{Q}}^{0}; while for T(0) = 10.61, to {\stackrel{\u0304}{Q}}_{>}^{*}. Notice that increasing β_{2} in 4fold, T(0) increased 10fold, showing that negative influence on cancer cells by normal cells affects strongly (very sensitive) in the cancer dynamics. For {\beta}_{2}>{\beta}_{2}^{th}, all initial conditions are attracted to trivial disregarding initial conditions. Figure 8 shows an extreme example, changing only β_{2} in 100 times the values of parameters given in Table 2, or β_{2} = 1.0. In this simulation we consider a very high (unrealistic) initial conditions T(0) = 100.0, even though {\stackrel{\u0304}{Q}}^{0} is the attractor: we observe that T decreases quickly (practically in the vertical axis), and A decreases slowly. We stress the fact that we used as the initial condition for C, the equilibrium value of normal cells, that is, C(0) = C_{0} = 9.0, because the amount C_{ m }in the initial condition for T is higher than C_{0}.
As β_{1} increases, the initial number of cancer cells (for instance, originated by mutation) needed to trigger a cancer T(0) = C_{ m }is decreased, but smoothly. Also, the normal cells can be displaced by cancer cells for values of β_{1} higher than its critical value, that is, {\beta}_{1}>{\beta}_{1}^{c}. With respect to β_{2}, we observed a threshold for {\beta}_{2},{\beta}_{2}^{th}, above which cancer can not be settled. Also, as β_{2} increases, the initial number of cancer cells needed to trigger a cancer C_{ m }increases quickly, avoiding the process of cancer disease. Therefore, cancer can be settled in an organ if the following combination matches: better fitness of cancer cells (β_{1} increases), and decrease in the efficiency of the repairing systems (β_{2} decreases).
Recruitment of existing epithelial cells  γ
New network of blood vessels is created by cancer cells to provide nutrients and oxygen to support their growth. This new network depends on the capacity of originating sprouts from the existing blood vessel network, and is described by the parameter γ, which is varied.
In Figure 9 we show the coordinates of the equilibrium points {\stackrel{\u0304}{Q}}_{<}^{*} (a) and {\stackrel{\u0304}{Q}}_{>}^{*} (b) by varying γ. The small T_{<} decreases, while the big T_{>} decreases after an initial increase. Notice that T_{>} decreases due to the fact that preexisting blood vessels E decreases quickly. The new blood vessels A also decrease. The maximum of the variables occurs at around γ = 6.57 × 10^{3} (for C is minimum). In Figure 9(b) we show the curve T_{ E }, which depends on γ and situates always above the equilibrium value T_{>}. Hence, we have nontrivial equilibrium point {\stackrel{\u0304}{Q}}^{*} for sufficiently higher values of γ.
In Figure 10, we show the bifurcation diagram, considering \stackrel{\u0304}{T} as a function of γ (curve T in (a) and (b) of Figure 9, which appear at γ = γ^{th}). In (b), a zoom near zero is shown. In Figures 10(a) and 10(b), for γ > γ^{th}initial conditions set in regions marked with I, Ia and Ib are attracted to the trivial equilibrium point {\stackrel{\u0304}{Q}}^{0}. However, initial conditions set in regions marked with II and III are attracted to {\stackrel{\u0304}{Q}}_{>}^{*}\left({\gamma}^{th}<\gamma <{\gamma}_{1}^{c}\right); to a limit cycle circulating {\stackrel{\u0304}{Q}}_{>}^{*} in regions IIa and IIIa \left({\gamma}_{1}^{c}<\gamma <{\gamma}_{2}^{c}\right); and to {\stackrel{\u0304}{Q}}^{0} in regions IV and V \left(\gamma >{\gamma}_{2}^{c}\right), where E decreases and, then, increases to equilibrium value E_{0}. The Hopf bifurcation occurs at \gamma ={\gamma}_{1}^{c} (supercritical) and \gamma ={\gamma}_{2}^{c} (subcritical) [17]. The special values are: γ^{th}= 5.450 × 10^{4} and {T}_{>}\left({\gamma}^{th}\right)=1.74,{\gamma}_{1}^{c}=2.3196\times 1{0}^{2} and {T}_{>}\left({\gamma}_{1}^{c}\right)=2.1029, and {\gamma}_{2}^{c}=2.5161\times 1{0}^{2} and {T}_{>}\left({\gamma}_{2}^{c}\right)=1.9443. Figure 9(b) showed that at \gamma \simeq {\gamma}_{1}^{c}, E is very low. For γ < γ^{th}(b), all initial conditions set in region marked with Ic are attracted to the trivial equilibrium point {\stackrel{\u0304}{Q}}^{0}, which is the unique equilibrium, because there is not any positive solution for equation (5). Hence, there is a threshold of the parameter γ, denoted γ^{th}, below which all trajectories go to trivial equilibrium.
In cancer disease, it is expected that a small number of epithelial cells must be recruited in order to build up new vascularization from the sproutings. For small values of γ but higher than the threshold (γ > γ^{th}), the dynamical trajectories follow those shown in Figure 1. In Appendix C we illustrate the Hopf bifurcation, which behavior is not compatible with real cancer. Considering number of tumor cells, concentration of growth factor and volume of blood vessels feeding the tumor, Agur et al. [18] showed that Hopf bifurcation can not occur if ordinary differential equations are used. But, Hopf bifurcation can occur if timedelay is encompassed. Our model presented Hopf bifurcation, however only in a range of values of parameter γ which is not compatible with biological findings.
Solid tumors need extra source of resources to attend the quick growth of cancer cells. Hence, cancer can be settled in an organ if the capacity of sprouting from existing vascularization is sufficiently higher (γ > γ^{th}). However, it must not be so higher in order to avoid the death of the cancer diseased person due to normal cells being displaced by cancer cells quickly.
Capacity of building up new vascularization  ε
The appearance of shunts from existing blood vessels to initiate new vascularizations was described by the parameter γ. New network of blood vessels is created by cancer cells to provide nutrients and oxygen to support their growth. After a period of time δ^{1}, new vessels are built up from the shunts. The capacity of mounting up new vessels by cancer cells is analyzed by varying ε.
In Figure 11 we show the coordinates of the equilibrium points {\stackrel{\u0304}{Q}}_{<}^{*} (a) and {\stackrel{\u0304}{Q}}_{>}^{*} (b) by varying ε. Both small T_{ < }and big T_{ > }decrease. In (b), due to the solution T_{ > }(ε = 0) = 5.0, we used equation (12) for f(T) and g(T) to obtain Ā(0) = 4.457.
In Figure 12, we show the bifurcation diagram, considering \stackrel{\u0304}{T} as a function of ε (curve T in (a) and (b) of Figure 11). When γ > 0, there are not neither special nor threshold values: when initial conditions are set in a small region I, trajectories are attracted to the trivial equilibrium point {\stackrel{\u0304}{Q}}^{0}, and for initial conditions set in regions II and III, trajectories are attracted to {\stackrel{\u0304}{Q}}_{>}^{*} (a). This behavior results by the existence of influx in equation for A, given by the term γP. Hence, if we let γ = 0, a different bifurcation arises (b).
There is a threshold of ε, denoted ε^{th}, below which all trajectories go to trivial equilibrium. For ε > ε^{th}, where ε^{th}= 5.250 × 10^{2}, we have similar behavior than that observed in (a), but T_{>} is increasing. In Figure 13 we show the dynamical trajectories depending on the initial conditions for γ = 0, changing only ε in 100 times the values of parameters given in Table 2, or ε = 1.0. Dynamical trajectories are attracted to {\stackrel{\u0304}{Q}}^{0} (a) when T(0) = 0.416, while for T(0) = 0.417, to {\stackrel{\u0304}{Q}}_{>}^{*} (b). The dynamical trajectories in (a) are similar to Figure 1(a). Increasing ε in 100fold, T(0) decreased 2.5fold.
Due to the influx γP, which is a linear term, in the equation for A, the parameter ε that describes the mounting up of new vascularization does not present any special behavior, except by the dependency with the initial conditions. However, when γ = 0, there arises a threshold of ε, called ε^{th}, above which new vascularizations promoted by cancer cells can occur. The case γ = 0 means that cancer cells are supported exclusively by the new network of blood vessels originated from surrounding tissues for instance, and the preexisting one maintains its function of nourishing exclusively the normal cells.
Angiogenesis is the process by which new blood vessels develop from an existing vasculature, through endothelial cell sprouting, proliferation, and fusion. Hence, angiogenesis create new vascularization from sprouting originated in existing vasculature, which was called preangiogenesis. Due to the endothelial cell sprouting promoted by the vascular endothelial growth factor, cancer cells can grow even in the absence of new vascularization (ε = 0), being the size of cancer cells big with higher capacity of mounting up new vascularization (increasing ε).
Discussion
In foregoing section we have used for k_{3} and k_{4} values comparable to k_{1} and k_{2} in order to enhance the results. In other words, cancer related cells are allowed to grow comparable to the size of normal cells, which is not true. In real world, cancer related cells are found in much smaller size, hence k_{3} and k_{4} must be lower than k_{1} and k_{2}, being the relative size depending on the organ of the body.
In Figure 14 we show dynamical trajectories using the initial conditions given in equation (2). The values of the parameters are given in the column marked with ** of Table 2, and other parameters are those given in fixed values. Dynamical trajectories are attracted to {\stackrel{\u0304}{Q}}^{0} (a) when T(0) = 9.51 × 10^{3}, while for T(0) = 9.52 × 10^{3}, to {\stackrel{\u0304}{Q}}_{>}^{*} (b). The cancer is triggered at around 400 days. Decreasing k_{3} in 50fold (and some values of other parameters were also changed), the initial value that divides two attracting regions T(0) is around 100fold lower than that shown in Figure 1. In Figure 1, T reaches the asymptotic value near k_{3} = 5.0, while in Figure 14, near k_{3} = 0.1. Then all previous simulations can be translated to real world by appropriate scaling factors.
The initial amount of cancer cells originated from normal cells by mutations plays an important role in the dynamics of cancer growth. When the interacting parameters β_{1}, γ and ε increase, the small solution T_{<} decreases. This decrease in the initial cancer cells necessary to trigger the cancer is mediated by cancer cells that: (1) inhibit and occupy normal cells habitat (β_{1}); (2) produce higher amount of substances which cause new blood vessels to grow (γ); and (3) construct new blood networks to nourish themselves (ε). In opposite way, when the interacting parameter β_{2} increases, the small solution T_{<} also increases. Hence, when normal cells inhibit the growth of cancer cells by many factors as better fitness, increasing repairing action, and induction of apoptosis, then the onset of cancer is avoided.
Another important aspect in the cancer growth is the appearance of thresholds. The interacting parameters β_{2} and γ present threshold values, respectively, {\beta}_{2}^{th} and γ^{th}. Usually the systems that control the production of substances that induce the formation of new blood vessels to grow operate normally, which have as consequence that cancer cells are unable to recruit the blood to supply their need to continue to proliferate, and they fade out at this early stage. However, cancer cells may begin to produce substances which cause new blood vessels to grow. This phenomenon is characterized by the threshold of γ. The threshold of β_{2} can be understood as the well functioning repairing mechanisms and the low fitness of cancer cells in comparison with normal cells. The threshold of ε arises only for γ = 0, which mimics cancer cells being nourished only by the new network of blood vessels originated from surrounding tissues by vasculogenesis [14].
There are also critical values for parameters β_{1} and γ. Critical value for β_{1} is a mathematical artefact, because it is meaningless biologically (in general k_{3} is very low in comparison with constraint T_{ E }). With respect to γ, there are two critical values, named {\gamma}_{1}^{c} and {\gamma}_{2}^{c}. When γ increases, the epithelial cells E decrease, and damped oscillations appear. When γ approaches to {\gamma}_{1}^{c}, the oscillations are less damped, and when surpasses {\gamma}_{1}^{c}, regular oscillations occur. However, the amplitude of regular oscillations increases as γ approaches to {\gamma}_{2}^{c}, resulting for lowest values of E, T and A reaching zero values (see figures in Appendix C). When the lowest values are incapable to trigger new burst of cancer cells, the oscillations cease and the trivial is the attractor. This occurs when γ surpasses {\gamma}_{2}^{c}. Again both {\gamma}_{1}^{c} and {\gamma}_{2}^{c} do not bear any biological meanings, because the cure of cancer is due to the elimination of preexisting network of blood vessels. The biological meaningless sustained oscillations begin at \gamma ={\gamma}_{1}^{c} (the supercritical Hopf bifurcation) and cease at \gamma ={\gamma}_{2}^{c} (the subcritical Hopf bifurcation).
For instance, considering values of Table 2, the column with **, we observe the same behavior than that found in figures shown in Appendix C as γ increases: (1) damped oscillations around {\stackrel{\u0304}{Q}}^{*} for γ = 0.75; (2) regular oscillations around {\stackrel{\u0304}{Q}}^{*} in the interval [0.76,0.80]; and (3) to the trivial equilibrium {\stackrel{\u0304}{Q}}^{0} for γ = 0.81. In the region of limit cycle, we observed that: (1) C oscillates between (8.9, 9.0) for all γ; E, P and A oscillate between (0, M), where M increases as γ increases; and (3) T oscillates between (m, 0.1), where m decreases as γ increases. When m becomes very small, there is not burst of cancer cells, limit cycle is destroyed, and cancer fades out for higher values of γ. Notice that the amplitude of oscillations of normal cells C is not affected, in opposite way of epithelial cells E which drops to near 0.
Agur et al. [18] showed that ordinary differential equations admit Hopf bifurcation if and only if at least one timedelay is introduced in the tumor growth modeling. They then concluded that an appropriate candidate for describing the cancer growth is the alternative that includes timedelays in the tumor proliferation or angiogenesis process. They also concluded that further mathematical research is warranted for exploring timedelays in the biologically realistic domains in the parameter spaces. Our model, however, showed Hopf bifurcation without timedelay (in fact, there is an elapse of time between preangiogenesis and angiogenesis cells), and sustained oscillations occur only in biologically not realistic domains in the parameter space.
When the initial conditions, especially T(0) = C_{ m }, are such that the nontrivial equilibrium is attractor, the cancer cells reach the level T_{>}. C_{ m }increases with increasing β_{1}, decreases with β_{2} and ε. Cancer cells can grow and reach higher levels when they affect negatively normal cells (β_{1}), but reach lower levels when normal cells acts as a barrier against them (β_{2}). When γ varies, T_{>} increases in the initial phase (E decreases), and then decreases (E is practically zero). The parameter γ plays an equivalent role of β_{1}, but, restricted only to E, which decreases it dramatically. For lower γ, there is sufficient number of E to increase A, but epithelial cells are exhausted as γ increases, and A decreases (see equation (4)). Finally, T_{>} decreases with ε, which comes out due to relative higher value of γ. The behavior of ε is strongly dependent on γ due to the influx γP in equation for A: for lower γ (also γ = 0), T_{>} increases with ε (see Figure 12(b)).
We introduced in the model an intermediate phase between epithelial cells and angiogenesis cells. The purpose was to consider a delay in new blood vessel formation (angiogenesis) by the period of time δ^{1}. This can be suppressed by letting δ → ∞. In Figure 15 we show dynamical trajectories using the initial conditions given in equation (2), and values of the parameters given in the column marked with ** of Table 2, and other parameters are those given in fixed values, except δ. Dynamical trajectories are attracted: (1) for δ = 0.001, to {\stackrel{\u0304}{Q}}^{0} (not shown) when T(0) = 2.492 × 10^{1}, while for T(0) = 2.493 × 10^{1}, to {\stackrel{\u0304}{Q}}_{>}^{*} (a); and (2) for δ = 10.0, to {\stackrel{\u0304}{Q}}^{0} (not shown) when T(0) = 6.608 × 10^{3}, while for T(0) = 6.609 × 10^{3}, to {\stackrel{\u0304}{Q}}_{>}^{*} (b). The cancer is triggered at around 900 and 360 days, respectively for δ = 0.001 and 10.0. Including Figure 14(b), the cancer trigger is delayed and initial cancer formation due to mutation must be increased as δ decreases.
Conclusions
Many models [5] [19] have already been proposed to describe cancer growth, and some of those models have explicitly considered the spatial dimension [20–22], which has been shown to play a key role in the understanding of various tumor growth processes. Other models considered computational approach [23, 24]. Spatial and computational modelings were analyzed numerically. However, we developed a zerodimensional model for the initial stages of tumor angiogenetic growth in order to obtain analytical results.
In this model, various cell species are supposed to satisfy LotkaVolterra growth laws. Neither metastasis nor geometry of the solid cancer were taken into account. The initial conditions (2), supplied to the dynamical systems (1), describe an impulsive system: in a steady state, a perturbation is introduced at a time t = t_{0} = 0 as a form of pulse (Dirac delta function). This pulse mimics normal cells mutating to cancer cells. Cancer cells then promote the mounting of new network of blood vessels to nourish them after an elapse of time δ^{1}. We introduced the preangiogenesis cells P to include the time delay in order to completing functional angiogenesis.
From the model, we conclude that the dynamical trajectories depend on the initial conditions supplied to the system, and also on interacting parameters. The cumulative effects of mutation is essential to originate a cancer cell. This effect is captured by the initial amount of cancer cells originating from normal cells, denoted by C_{ m }. A sufficient number of cells must suffer mutation in order to a concentration of C_{ m }cells really bear all necessary mutations to become effectively cancer cells. Our model is spatially homogeneous, hence the initial number of cancer cells is C_{ m }× V, where V is the volume of an organ of human body. In extremely favorable environmental and individual conditions, this initial number can be one.
Initially, cancer cells always can grow. But they fade out if they are unable to build up new blood vessels in order to supply their needs. The capacity of inducting new vascularization from existing blood vessel network must be efficient (γ > γ^{th}), and better fitness (increasing capacity of proliferation and capturing nutrients, decreasing mortality, etc.) of cancer cells (β_{1} increases) in comparison with normal cells. Another aspect of cancer growth is corruption of the repairing systems, and in some extend we can think of that normal cells influencing negatively cancer cells play the role of repairing (fixing mutated DNA and inducting apoptosis). The parameter β_{2} measures the efficient action of repairing system, and the effect of decreasing this value result in a higher level of corruption in the repairing system. Since cancer cells can recruit epithelial cells to form new blood vessels, the capacity of proliferation of new vessels (ε) does not present threshold value.
The parameter γ can be thought of deviation of nutrients and oxygen from normal cells to feed cancer cells. This deviation does not affect normal cells because the carrying capacity of normal cells k_{1} does not depend on the size of network of blood vessels. As γ increases, the deviated network plus the new one mounted by the action of angiogenesis effectors nourish the cancer cells. Hence, γ = 0 means that preexisting network of blood vessels feeds normal cells, while the new network nourishes cancer cells. In this case, the capacity of cancer cells in promoting new vascularization is essential, which must surpass the threshold value ε^{th}. We obtained biologically feasible nontrivial equilibrium points, that is, the coordinates of the equilibrium points are positively defined. However, due to the simplifications assumed by model, the range of variations of parameters like β_{1} and γ must be restricted. When {\beta}_{1}>{\beta}_{1}^{c} we have an equilibrium point {\stackrel{\u0304}{Q}}^{c} with \stackrel{\u0304}{C}=0, while for {\gamma}_{1}^{c}<\gamma <{\gamma}_{2}^{c}, we have limit cycle (sustained oscillations) with E ~ 0, but for \gamma >{\gamma}_{2}^{c}, we have abrupt increase of E, and the trivial equilibrium is attained. Both results can not be acceptable for cancer growth description.
Some results obtained here can be understood as vasculogenesis (when γ = 0) [14]. The dynamics of normal and cancer cells are similar than that presented in the model proposed by Nani and Freedman [11]. However, we did not take into account the action of immune system, while they did not take into account the angiogenesis. Agur et al. [18] proposed to examine the occurrence of Hopf bifurcation in the clinical context, that is, to check whether or not one can contain tumor growth by imposing timedelays in the processes of neovascularization. Our results showed that Hopf bifurcation occurs in biologically not realistic domains in the parameter space.
In a future work we will analyze a model in which the sizes of the normal and cancer cells are allowed to depend on the overall network of blood vessels: normal and cancer cells compete for nutrients provided by the preexisting blood vessels, while cancer cells have additional source originated from angiogenesis. If we take these effects into account in the model, maybe Hopf bifurcation can be avoided. There are several ways to improving model given in equation (1). One is the dependency of normal cells with the size, which decreases with increasing γ, of the existing vasculature.
For instance, we can deal with intermittent process instead of a continuous process of sprouting from existing blood vessels. We can change the second and fourth equations by
where Φ(t) is the total intermittent sproutings rate. One form of Φ(t) is
where γ_{ i }and τ_{ i }(τ_{0} = 0) are, for i = 0, ···, n, respectively, the ith sprouting rate and the time interval during which sproutings occur. The Heaviside function θ(x) is such that θ(x) = 1 if x > 0, and θ(x) = 0, otherwise. To be intermittent, we must have γ_{2i}> 0 and γ_{2i+1}= 0, that is, sproutings occur in the time interval [0, τ_{1}] with rate γ_{0}, do not occur in the interval [τ_{1}, τ_{2}], and so on. Another possibility is
where s_{ i }and τ_{ i }are, for i = 0, ···, n, respectively, the ith proportion of epithelial cells generating sproutings and the time at which sproutings occur. The proportion s_{ i }can depend on T. Vaccination campaigns against viral infections were analyzed considering age interval vaccination (Heaviside function) [25] or a series of pulses (Dirac function) [26].
Another improvement of the model (1) is the introduction of the immune response. This can be done by introducing lymphocyte cells action as Nani and Freedman dealt with [11].
In the model (1), chemotherapy that acts specifically against tumor cells can be introduced easily. An intermittent chemotherapy can be introduced in the model by adding one term in the first equation, that is,
and adding an equation for the drug administration as
where Q is the concentration of drug at time t, and μ_{ Q }and μ_{ q }are the rates of, respectively, intake of drug by cancer cells and elimination of drug by body. The parameter l is the amount of drugs intake by one cancer cells, and q(t) is the drug administration rate. The kinetics of drug intake follows MichaelisMenten [27]. If q(t) = q, a fixed value, we have a continuous regimen of administration. Intermittent drug administration can be considered. First, we can define
where u_{ i }is, for i = 0, ···, n, the ith concentration of drug administered at time τ_{ i }. Another is
where q_{ i }and τ_{ i }are, for i = 0, ···, n, the ith drug administration rate and the time interval during which drug is administrated. To be intermittent, we must have q_{2i}> 0 and q_{2i+1}= 0.
In this paper we obtained threshold values, but we did not deal with the effects of controlling mechanisms, which are left to a further work. Briefly, for instance, two parameters can be used in order to control cancer growth. The thresholds of the parameters γ and β_{2}, named γ^{th}and {\beta}_{2}^{th}, can be managed as follows.
Let us suppose that γ and β_{2} have values above and below their respective thresholds. To control cancer, parameters of the model (including dynamics of controls) must be managed to increase γ^{th}to surpass γ, and to decrease {\beta}_{2}^{th} to fall below β_{2}. Another way to control cancer is acting on the probability of sproutings becoming angiogenesis cells, in order to decrease δ/(δ + μ_{4}), which can be done by decreasing the rate of transformation from preangiogenesis to angiogenesis cells δ and/or increasing the mortality rate μ_{4} of preangiogenesis cells.
Appendix A: Nontrivial equilibrium point
The nontrivial equilibrium value of cancer cells corresponding to model (1), \stackrel{\u0304}{T}, is the positive solution of the equation (5), that is, f(T) = g(T). The fourth degree polynomial f(T) is such that
and has three nonnegative roots: 0, (α_{2}  μ_{2})/γ = T_{ E }, and k_{3} = T_{ A }is a double root (see equation (7)). Figure 16 shows the qualitative behavior of f(T): T_{ A }< T_{ E }(a), and T_{ A }> T_{ E }(b). Notice that Figure 16(a) corresponds to higher carrying capacity for cancer cells.
The function g_{1}(T) is such that
and has three positive roots: ε/μ_{5}, k_{3} = T_{ A }, and \frac{{\alpha}_{1}{\mu}_{3}}{{\beta}_{1}{\beta}_{2}{k}_{1}}+\frac{{\alpha}_{1}{\mu}_{1}}{{\beta}_{1}}={T}_{g} (see equations (7) and (8)). The function g_{2} (T) is such that
and has two nonnegative roots: 0, and \frac{{\alpha}_{1}{\mu}_{3}}{{\beta}_{1}{\beta}_{2}{k}_{1}}+\frac{{\alpha}_{1}{\mu}_{1}}{{\beta}_{1}}={T}_{g} (double root). The dominant terms of the third degree polynomials g_{1}(T) and g_{2}(T), when T → ± ∞, are
The third degree polynomial g(T) is sum of g_{1}(T) and g_{2}(T). Notice that, for T ≥ 0, we have g_{2}(T) ≥ 0, while g_{1}(T) changes signal. Hence, the effect of g_{2}(T) in the sum (g(T)) is the increasing in the values of g_{1}(T), except at T = T_{ g }, at which we have g_{1}(T_{ g }) = g_{2}(T_{ g }). The asymptotic behavior of the polynomial g(T) depends on the velocities of increasing v_{2} and v_{2} (or coefficients) of T^{3}. When β_{1} β_{2} < α_{3} α_{1} k_{4}/(k_{1} k_{3}) (generically referred to as weak interaction between normal and cancer cells when values of β_{1} and β_{2} are proportional), in which case T_{ g }is higher, we have
It can be shown that there is one solution or three positive solutions. When β_{1} β_{2} > α_{3} α_{1} k_{4}/(k_{1} k_{3}) (generically referred to as strong interaction between normal and cancer cells when β_{1} and β_{2} are proportional), in which case T_{ g }is lower, we have
and there is one negative solution and two positive solutions. In both cases T_{ g }is always a positive root of g(T), and g_{0} is given by equation (A.1). The existence of other positive solutions depends on the relative position of the roots of g_{1}(T) and the coefficients of the polynomials.
Figure 17 shows the qualitative behavior of g(T) for small β_{1} β_{2} (weak interaction): when T_{ g }is the greatest root of g_{1}(T), with three positive solutions (a), and one solution (b); when T_{ g }is between the roots ε/μ_{5} and k_{3}, with three positive solutions (c); and when T_{ g }is the smallest root of g_{1}(T), with three positive solutions (d). T_{1} and T_{2}, the roots of g_{1}(T), stand for ε/μ_{5} and k_{3}, depending on the relative positions between them. Figure 17(a) and 17(b) shows clearly the effect of g_{2}(T) increasing g_{1}(T) and changing the roots of g(T), named {T}_{1}^{g} and {T}_{2}^{g}:{T}_{1}^{g}>{T}_{1} and {T}_{2}^{g}<{T}_{2}, and both {T}_{1}^{g} and {T}_{2}^{g} disappear when g_{2}(T) is sufficiently higher.
The positive solution of the equation (5) is the intersection between the curves f(T) and g(T), or roots of f(T)g(T). Notice that:

1.
We do not have negative solutions, because f(T) < 0 and g(T) > 0, for T < 0, which implies that f(T)  g(T) < 0. At T = 0 we have f(0)  g(0) = g_{0} < 0.

2.
We have f(∞)  g(∞) → ∞, even when g(∞) → +∞, because f(T) is fourth degree polynomial, and g(T), third degree.

3.
f(T)  g(T) is a continuous function, a fourth degree polynomial, where there is not solution for T < 0, f(0)  g(0) < 0 and f(∞)  g(∞) < 0. Hence, by the intermediate value theorem, we have 0, 2 or 4 positive solutions in the interval (0, ∞).

4.
We are searching biologically feasible equilibrium, hence, according to equation (7), \stackrel{\u0304}{T} must be lower than the lowest root T_{ m }of f(T), that is, \stackrel{\u0304}{T}\in \left(0,{T}_{m}\right), where T_{ m }= min {(α_{2}  μ_{2})/γ, k_{3}}, the minimum between (α_{2}  μ_{2})/γ and k_{3}. Therefore, in the interval (0, T_{ m }) we have at most 2 positive solutions.
Figure 18 shows qualitative behavior of f(T)  g(T). Figures 18(b) and 18(c) (first two lower roots) show two positive solutions satisfying \stackrel{\u0304}{T}<{T}_{A},\stackrel{\u0304}{T}<{T}_{E} and \stackrel{\u0304}{T}<{T}_{g}. Additionally, the constraint \stackrel{\u0304}{T}<{T}_{C} are satisfied, once \stackrel{\u0304}{T}<{T}_{g}. Hence, biologically feasible solutions are at most 2. However, Figures 18(c) (last two higher roots) and 18(d) show two positive solutions not biologically feasible, because they are greater than the constraint T_{ A }= k_{3}.
When there is not any positive solution for f(T)  g(T) = 0, the unique equilibrium point is the trivial {\stackrel{\u0304}{Q}}^{0}. In the case of two positive solutions, the small root (T_{ < }) forms the unstable equilibrium {\stackrel{\u0304}{Q}}_{<}^{*}, and the big one (T_{>}) forms a possibly stable (let us for simplicity say that it is stable) equilibrium {\stackrel{\u0304}{Q}}_{>}^{*}. Between two stable equilibria {\stackrel{\u0304}{Q}}_{>}^{*} and {\stackrel{\u0304}{Q}}^{0} (this always exists and is stable), we have an unstable point {\stackrel{\u0304}{Q}}_{<}^{*}. We call as the 'breakpoint' [28, 29] the unstable equilibrium point {\stackrel{\u0304}{Q}}_{<}^{*}, which separates two attracting regions containing one of the equilibrium points {\stackrel{\u0304}{Q}}_{>}^{*} and {\stackrel{\u0304}{Q}}^{0}. In another words, we have a hypersurface obeying, e.g., h\left(\stackrel{\u0304}{C}\left({T}_{<}\right),\u0112\left({T}_{<}\right),{T}_{<},\stackrel{\u0304}{P}\left({T}_{<}\right),\u0100\left({T}_{<}\right)\right)=0, generated by the coordinates of the equilibrium point {\stackrel{\u0304}{Q}}_{<}^{*} such that one of the equilibrium points {\stackrel{\u0304}{Q}}_{>}^{*} and {\stackrel{\u0304}{Q}}^{0} is an attractor depending on the relative position of the initial conditions supplied to the dynamical system (1) with respect to the hypersurface h [30].
Appendix B: The local stability of the nontrivial equilibrium {\stackrel{\u0304}{Q}}^{c}
Let us show that the matrix A, given by equation (14), is an Mmatrix [31] for big nontrivial equilibrium {\stackrel{\u0304}{Q}}_{>}^{c}, where T_{>} is one of the coordinates.
Definition. We say that the n × n matrix A = [a_{ ij }] is a nonsingular Mmatrix if a_{ ij }≤ 0, i ≠ j, and there exists a matrix B ≥ 0 and a real number s > 0 such that
where I is the identity matrix and ρ is the spectral radius [32].
Or, equivalently:
Proposition 1. A is a nonsingular Mmatrix if and only if the real part of its eigenvalues is greater than zero.
Proposition 2. A (elements a_{ ij }) is a nonsingular Mmatrix if and only if the diagonal entries are positive, and there exists a positive diagonal matrix D (diagonal elements d_{ i }> 0), such that AD is strictly diagonal dominant, that is,
for i = 1, 2, ..., n.
According to the first part of Proposition 2, the matrix A has positive diagonal elements, see equation (14).
The second part of Proposition 2 is written as
because Ā can be greater than k_{4}. The equilibrium values correspond to the point {\stackrel{\u0304}{Q}}^{c}, given by equation (9).
Let us define
where φ > 0. By these definitions, the first three inequalities of equation (B.1) hold. To prove the last inequality, we substitute above definitions, and we obtain
where the numerator φ_{ n }and denominator φ_{ d }> 0 are
with the last term of the numerator being obtained using the relation
which is valid in the equilibrium. Since φ_{ d }> 0, if we show that φ_{ n }> 0, then there exists a positive number φ.
Let us consider the case where there are two positive solutions for the equation f(T) = g(T). In this case, we show that the bigger equilibrium {\stackrel{\u0304}{Q}}_{>}^{c}, with coordinate T_{>} as solution of (5), is stable. Additionally, we assume that Ā > k_{4}, which seems reasonable since the equilibrium {\stackrel{\u0304}{Q}}^{c} is stable for \stackrel{\u0304}{T}>\left({\alpha}_{1}{\mu}_{1}\right)/{\beta}_{1}={T}_{C}. Substituting the coordinates of the equilibrium {\stackrel{\u0304}{Q}}^{c}, the numerator of φ can be written as
with
where
and we must show that \Phi \left(\stackrel{\u0304}{T}\right)>0, for higher value of \stackrel{\u0304}{T}, that is, \stackrel{\u0304}{T}={T}_{>}. Let us assume that k_{4}α_{3} > μ_{3}. Then, we have \Phi \left(\stackrel{\u0304}{T}\right)>0, for \stackrel{\u0304}{T}\ge 0, because the discriminant of {\Phi}_{1}\left(\stackrel{\u0304}{T}\right) is
Hence, {\Phi}_{2}\left(\stackrel{\u0304}{T}\right) determines the existence of positive φ_{ n }. The discriminant of {\Phi}_{2}\left(\stackrel{\u0304}{T}\right) is
where T_{ E }and T_{ A }are given in equation (7). We have two possibilities. First, when T_{ E }< T_{ A }, we have {\Phi}_{2}\left(\stackrel{\u0304}{T}\right)>0, for all \stackrel{\u0304}{T}\ge 0, and \Phi \left(\stackrel{\u0304}{T}\right)>0 for
where
Hence, when ε > ε_{min}, φ_{ n }> 0 and we have a positive number φ obeying equation (B.2). Second, when {T}_{E}>{T}_{A},{\Phi}_{2}\left(\stackrel{\u0304}{T}\right) has two positive solutions Φ_{2<}and Φ_{2>}given by
and we have {\Phi}_{2}\left(\stackrel{\u0304}{T}\right)\le 0, for {\Phi}_{2<}\le \stackrel{\u0304}{T}\ge {\Phi}_{2>}; otherwise, {\Phi}_{2}\left(\stackrel{\u0304}{T}\right)>0. Notice that Φ_{2>} is out of the range of feasibility, once Φ_{2>} >T_{ E }. We know that \left({\alpha}_{1}{\mu}_{1}\right)/{\beta}_{1}<\stackrel{\u0304}{T}<{k}_{3}, but the bigger solution under consideration is {T}_{>}=\left({\alpha}_{1}{\mu}_{1}\right)/{\beta}_{1}^{c}={T}_{E}^{c}, for {\beta}_{1}\ge {\beta}_{1}^{c}. It is easy to show that Φ_{2}(k_{3}) < 0; but, {\Phi}_{2}\left({T}_{E}^{c}\right)\le 0, if {T}_{E}^{c}\ge {\Phi}_{2<}, and {\Phi}_{2}\left({T}_{E}^{c}\right)>0, if {T}_{E}^{c}<{\Phi}_{2<}. Hence, {\Phi}_{2}\left(\stackrel{\u0304}{T}\right)\le 0 if {T}_{E}^{c}<{\Phi}_{2<}. Hence, \Phi \left(\stackrel{\u0304}{T}\right)>0, for \stackrel{\u0304}{T}\ge 0, when {T}_{E}^{c}\ge {\Phi}_{2<}, or when {T}_{E}^{c}\ge {\Phi}_{2<} and \stackrel{\u0304}{T}\ge {\Phi}_{2<}; and \Phi \left(\stackrel{\u0304}{T}\right)>0, when {T}_{E}^{c}<{\Phi}_{2<} and {T}_{E}^{c}<\stackrel{\u0304}{T}<{\Phi}_{2<} for ε > ε_{min}.
Summarizing, φ_{ n }> 0 occurs, in order to have positive number φ:
1. T_{ A }< T_{ E }or γ < (α_{2}  μ_{2})/k_{3}  weak capacity of recruitment of the normal epithelial cells by cancer cells. We have:1.a If {T}_{E}^{c}\ge {\Phi}_{2<}:{\phi}_{n}>0 without restriction about T_{>}.
1.b If {T}_{E}^{c}<\stackrel{\u0304}{T}<{\Phi}_{2<}: we have two possibilities
1.b.1 If T_{>} > Φ_{2<}: φ_{ n }> 0 without restriction about T_{>}.
1.b.2 If T_{>} < Φ_{2<}: φ_{ n }> 0 if ε > ε_{min}, where T_{>} satisfies ε_{min}, equation (B.4). Higher proliferation of angiogenesis cells must occur.

2.
T_{ A }> T_{ E }or γ > (α_{2}  μ_{2})/k_{3}  strong capacity of coopting normal epithelial cells by cancer cells: φ_{ n }> 0, if ε > ε_{min}, where T_{>} satisfies ε_{min}.
When γ is small, the big equilibrium {\stackrel{\u0304}{Q}}_{>}^{c}, with one coordinate T_{>}, is stable without conditions (cases 1.a and 1.b.1). However, for sufficiently higher values of γ, the big equilibrium point {\stackrel{\u0304}{Q}}_{>}^{c} can be unstable (case 2), in which case Hopf bifurcation can occur (see Appendix C).
In the case of the small equilibrium {\stackrel{\u0304}{Q}}_{<}^{c}, with one coordinate T_{<}, it is unstable. Assuming that Ā < k_{4}, we show here that A corresponding to small T_{<} is not an Mmatrix (we are not proving that {\stackrel{\u0304}{Q}}_{<}^{c} is unstable). In this case, equation (B.3) becomes
where
The bigger roots of {\Phi}_{1}\left(\stackrel{\u0304}{T}\right) and {\Phi}_{2}\left(\stackrel{\u0304}{T}\right) are not biologically feasible, because they are, respectively, higher than k_{3} and (α_{2}  μ_{2})/γ. Let us define {T}_{m}=\text{min}\left\{\underset{<}{\overset{1}{T}},{T}_{<}^{2}\right\}, where \text{min}\left\{\underset{<}{\overset{1}{T}},{T}_{<}^{2}\right\} is the minimum between the small roots of, respectively, {\Phi}_{1}\left(\stackrel{\u0304}{T}\right) and {\Phi}_{2}\left(\stackrel{\u0304}{T}\right) given by {T}_{<}^{1}={k}_{3}\left(1\sqrt{1\frac{{\mu}_{3}}{{k}_{4}{\alpha}_{3}}}\right), assuming that k_{4} ≥ 1, and {T}_{<}^{2}={T}_{E}\left(1\sqrt{1\frac{{T}_{A}}{{T}_{E}}}\right), assuming that T_{ E }≥ T_{ A }. Hence, if \stackrel{\u0304}{T}<{T}_{m}, we have \Phi \left(\stackrel{\u0304}{T}\right)>0 and there is not a positive number φ. Notice that, when T_{ E }< T_{ A }, we have {\Phi}_{2}\left(\stackrel{\u0304}{T}\right)>0, and it is enough to satisfy \stackrel{\u0304}{T}<{T}_{<}^{1}.
Appendix C: Hopf bifurcation
In Figures 19, 20 and 21 we illustrate the Hopf bifurcation (see Figure 10 in the main text), using values of parameters given in Table 2, except γ assuming higher values. The following figures are mathematical results, not cancer in an organ.
Dynamical trajectories are shown in Figure 19 for γ near {\gamma}_{1}^{c}:{\gamma}_{1}=2.31\times 1{0}^{2}\lesssim {\gamma}_{c}^{1} and {\gamma}_{2}=2.32\times 1{0}^{2}\gtrsim 1{0}^{2}\gtrsim {\gamma}_{1}^{c}. When T(0) = 0.44, trajectories of both cases go to {\stackrel{\u0304}{Q}}^{0} (trajectories similar to Figure 1(a), not shown). When T(0) = 0.45, dynamical trajectories for γ_{1} go to stable {\stackrel{\u0304}{Q}}_{>}^{*} (a), while for γ_{2}, they oscillate around unstable {\stackrel{\u0304}{Q}}_{>}^{*} (b), in which case limit cycle arises.
Now, we show the dynamical trajectories for γ near {\gamma}_{2}^{c}. Again, for {\gamma}_{3}=2.519\times 1{0}^{2}\lesssim {\gamma}_{2}^{c} and {\gamma}_{4}=2.520\times 1{0}^{2}\gtrsim {\gamma}_{2}^{c}, when T(0) = 0.41, trajectories of both cases go to {\stackrel{\u0304}{Q}}^{0} (trajectories similar to Figure 1(a), not shown). When T(0) = 0.42, Figure 20 shows dynamical trajectories for γ_{3} oscillating around unstable {\stackrel{\u0304}{Q}}_{>}^{*}. Notice that the amplitude of regular oscillations of the variables is very high: C and E (a), T (b), P (c) and A (d). For the same T(0) = 0.42, we show the dynamical trajectories for γ_{4} and γ_{5} = 2.519675 × 10^{2}, which is slightly higher than {\gamma}_{2}^{c}. In both cases we have \gamma >{\gamma}_{2}^{c}, and Figure 21 shows trajectories going to {\stackrel{\u0304}{Q}}^{0} with different number of oscillations: for γ_{4} (a), we have one oscillation, while for γ_{5}, three oscillations, showing C and E (b), T (c), P and A (d). When \gamma >{\gamma}_{2}^{c}, we observe that E goes initially to zero, and, after 400 days, increases abruptly, according to Figure 21(a). In Figure 5(a) we showed that C has similar behavior.