Mathematical modeling of solid cancer growth with angiogenesis

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 growth-promoting 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 growth-control 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) growth-promoting 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 by-products 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 time-consuming 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 a 1 and a 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 pre-vascular 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 pre-angiogenesis (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 pre-angiogenesis cells are formed by the mass action law [10] at a constant rate g, 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 pre-angiogenesis to angiogenesis cells. The cancer, pre-angiogenesis 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 a 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 α 3 . The rates at which normal and cancer cells compete themselves for resources and space are β 1 and β 2 . Both parameters also take into account effects of the environment.
Pre-angiogenesis cells deserve some words. New blood vessels to nourish cancer cells must sprout from pre-existing 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 time-delay 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 d The cancer cells influence negatively normal cells, and this term is described by α 1 Cβ 1 T k 1 (E) 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 β 1 = β 1 α 1 /k 1 and β 2 = β 2 α 3 /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 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 byQ = (C,Ē,T,P,Ā) . The coordinates of the equilibrium point are obtained by letting the time derivatives in equation (1) equal to zero, for instance, d dt C = 0 . The first is the equilibrium where there is not any cells, that is, Q abs = (0, 0, 0, 0, 0).
For obvious reason, this point is not considered here. The trivial equilibrium point, given bȳ corresponds to the absence of cancer. The trivial equilibrium point exists for a 1 >μ 1 and a 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 non-trivial equilibrium pointQ * are whereT is the positive solution of the equation

Rate of inhibition of cancer cells by normal cells
The summary of the variables and parameters of the model.
with the fourth degree polynomial f(T), and the third degree polynomials g 1 (T) and g 2 (T) being given by By inspectingC , Ē and Ā, the positive solutionT , in order toQ * be biologically feasible, must satisfy, respectively, the constraintsT < T C , However, the constraint for the equilibrium Ā isT < T g , where but it is always satisfied whenT < T C , since T g >T C .
With respect to the constraints, we must haveT < min{T C , T E , T A } , 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 b 1 and g, and directly with the carrying capacity k 3 . Hence, for higher values of removing parameters b 1 and g, and lower values of carrying capacity k 3 , the constraints are decreased, that is, the steady state cancer cellsT must assume lower values. Let us suppose that k 3 is high, which implies that cancer cellsT can achieve higher values. If the sproutings rate g from existing blood vessels is high, in principle it seems to be beneficial to cancer cells. However, higher values of g result in lower values for the constraint T E , and, consequently,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 g. 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 solutionsT depends on the values assigned to the model's parameters, which can be 0, 2 or 4. Notwithstanding, the solutionT 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 solutionsT as small (T < , and the corresponding equilibrium pointQ * < ) and big (T > , and the corresponding equilibrium pointQ * > ). 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 non-linear term parameters b 1 , b 2 , g and ε. First, let us vary b 1 . Notice that T C depends on b 1 , but T E and T A do not. For small b 1 , the solu-tionT of equation (5) situates below T C . However, as b 1 increases,T can surpass T C , after intercepting T C at β 1 = β c 1 . AsT increases with b 1 ,C , given by the first equation of (4), decreases and reaches zero at Let us suppose that T c C < min{T E , T A } , the minimum between T E and T A . WhenC = 0 , there arises another feasible equilibrium pointQ c , with the coordinates whereT = T c C , and the equilibriumT does not depend anymore on the interacting parameter b 1 . Hence for β 1 ≥ β c 1 we have a constant valueT = (α 1 − μ 1 )/β c 1 , which implies in constant values for Ē,P and Ā.
See more discussion below, in the stability analysis.
The equilibrium pointQ c corresponds to the case where the cancer cells displaced normal cells. By inspecting the equation forC , the first in equation (4), the normal cells drop out to zero for a sufficiently higher values of b 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 thatQ 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 valueT is always smaller than k 3 , and all coordinates of the equilibrium point are positive (Appendix A).
Second, let us now vary b 2 . The third equation of (1) in the equilibrium can be written asT When b 2 increases,T decreases, and at a certain value of b 2 , say β th 2 ,T assumes zero value. WhenT = 0 , the unique feasible equilibrium point is the trivialQ 0 . When the influence of the normal cells and the environment is higher ( β 2 ≥ β th 2 ), then cancer can not be maintained.
Cancer onset depends on the interplay among the initial condition T(0) = C m and interacting parameters b 1 and b 2 . Once cancer cells are originated by mutation from several sources, these cells interact with healthy cells and the surrounding environment. While the parameter b 1 is practically insensitive in the behavior of the dynamics, however, there is a threshold in the parameter b 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 g and angiogenesis increasing rate ε.
When g increases, Ē decreases according to the second equation of (4). Let us assume that Ē reaches zero at a certain value named g c , and Ē < 0 sinceafter. When Ē = 0, differently fromC = 0 , the unique feasible equilibrium point is the trivialQ 0 .
Hence, we expect that Ē > 0 for g ≥ 0, because ofT < T E = (α 2 − μ 2 )/γ is always obeyed, and there is not a critical value g c . When we decrease g, we observe that equation (5) may not have biologically feasible solutions.
Initially, let us analyze g = 0. In this case, we have f(T) = 0, and the equilibrium valuē T is solution of the polynomial of degree two with another being given byT = T g . Moreover, from equation (4), we haveP = 0 (due to g = 0, the role of the parameter δ does not matter). In this case, there arises another feasible equilibrium pointQ 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 g = 0. Hence, when ε >ε th , with ε th >μ 5 /k 3 , we have two positive solutions, named T ε < (and the corresponding equilibrium pointQ p < ) and T ε > (and the corresponding equilibrium pointQ p > ), which collapse to one at ε = ε th , and for ε <ε th there is not positive solution. The equilibrium pointQ 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 ofQ 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 (T(0) > T ε < ). Now, for lower values of g (g ≳ 0), as a consequence of the appearing of ε th for g = 0, equation (5) does not have positive solution and there is not biologically feasible equilibrium point. Let us introduce a threshold of g, named g th . For g <g th the unique solution isT = 0 (this fact is shown numerically). Note that as ε increases, g 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 equili-briumT , which must be substituted in the variables given in equation (4), must be This equation clearly shows that a positive solutionT is feasible for sufficiently higher values of g, that is, γ > γ th 0 , where γ th 0 corresponds to ε = 0. Due to the fact that ε increases as g th decreases, γ th 0 is the upper bound of the threshold of g. 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 (g >g 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 long-range signaling, and they are recruited from blood by the tumor by short-range 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 g) and the recruited endothelial progenitor cells (this phenomenon is not considered in the model; however, the case g = 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 pre-angiogenesis parameter δ, when δ → ∞,Q ∞ appears, similar to the equilibrium pointQ p , because lim δ→∞P = 0 .
However, from the third equation of (4), we have lim δ→∞ δP = γĒT . Hence, the coordinates ofQ ∞ are those given in equation (11), butT is solution of the equation f(T) 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 pointQ ∞ represents the absence of the intermediate pre-angiogenesis, that is, angiogenesis occurs directly and instantaneously from epithelial cells.
Summarizing, the non-trivial equilibrium pointsQ * (Q * < andQ * > , when they exist), Q c (for β 1 ≥ β c 1 ), andQ p (when g = 0) are biologically feasible, butQ c is not a real cancer equilibrium point, sinceC = 0 . Especially,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 β th 2 and g th (and, eventually, ε th , for g = 0), and a critical value β c 1 . 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 a 1 , a 2 and a 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 (b 1 , b 2 , g 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 (J − λI) = 0 , whereJ is the Jacobian J evaluated at the equilibrium point under analysis, The JacobianJ was obtained using equilibrium equations given in (4) assuming that C = 0 . WhenC = 0 , the element in the first line and first column is , and in the third column, a 13 = 0.
We present analytical results corresponding to the equilibrium pointsQ 0 andQ c .
In the last two eigenvalues we used the coordinates given by equation (3). Hence, the trivial equilibriumQ 0 is locally asymptotically stable whenever this point exists (a 1 >μ 1 and a 2 >μ 2 ), because all eigenvalues become negative.
The non-trivial equilibrium pointQ c , with coordinates given in equation (9), has one of them zero (C = 0 ). The eigenvalues corresponding to the trivial equilibrium point Q c are obtained from the JacobianJ c . OnceC = 0 , one of the eigenvalues is , arising one of the conditions to the not real cancer equilibriumQ c being stable, that is,T > T C , where T C is given by equation (7). The remaining four eigenvalues are obtained from the sub-matrixJ c 1 , which comes fromJ c excluding first row and first column. Let us define matrix A as A = −J c 1 , 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 M-matrix, and, hence, all eigenvalues associated toJ c 1 have negative real part. Hence,Q c , which is biologically feasible but not real cancer, is locally asymptotically stable only ifT > T C . Therefore,Q c appears when one of the constraints toQ * be biologically feasible is violated: in the equation (4), instead of C < 0 , we haveC = 0 for β 1 ≥ β c 1 . The local stability of the non-trivial equilibriumQ * is performed numerically. However, we stress the fact that the equilibriumQ c is an extension ofQ * for β 1 ≥ β c 1 . Hence, the big solution T > of equation (5) for β 1 ≥ β c 1 must be stable. For this reason we showed the local stability ofQ 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 Runge-Kutta (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 pointQ 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 pointQ 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 situationQ 0 , that is, at time t = 0 the variables assume C 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 b 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) (4). The corresponding five eigenvalues (l i , i = 1, ..., 5) are given in Table 3, showing that the small is unstable (l 1 > 0), while the big is stable (all l i 's are negative). The big equilibrium T > = 3.68 satisfies the constraints given in equation (7). The trivial equilibriumQ 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 The values for model's parameters. The unity of [•] is number of cells of type • per unit of volume, •/uv. Parameters regarded to the non-linear terms of the dynamical system are allowed to vary (indicated by the symbol*). The values of the parameters that differ from the fixed values are in the column indicated with **. Eigenvalues (l i , i = 1,..., 5) corresponding to the non-trivial equilibrium pointsQ * < (unstable) andQ * > (stable) using values given in Table 2.
initial conditions. In Figure 1(a), due to T(0) = 1.023, the dynamical trajectories are attracted to the trivial equilibriumQ 0 ; while for T(0) = 1.024, the dynamical trajectories attain the non-trivial equilibriumQ * > (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 ofQ * < , 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 trivialQ 0 when T(0) = (1 -) × 0.08, and to the non-trivialQ * > when T(0) = (1 + ) × 0.08 (for instance, = 0.001, figures not shown). In this case, we have T* = T < . Hence,Q * < is the break-point, and its coordinates generate a hyper-surface that divides two attracting regions (see Appendix A). Notice that T(0) = 1.023 which plays the role of separating two attracting regions is 12-fold 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 b 1 , b 2 , g 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 equilibriumQ 0 , and not those of the unstableQ * < . Hence, T(0) is not comparable with T < , one of the coordinates of the break-pointQ * < . The cancer cells proliferate above the subclinical threshold of 10 3 cells and reach 10 9 cells which is the X-ray detectable threshold [11].
Interaction between normal and cancer cells -b 1 and b 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  Table 2. Dynamical trajectories of the system (1) considering the values of parameters given in Table 2. The initial conditions determine the region of attraction: trivialQ 0 for T(0) = 1.023 (a), or non-trivialQ * > (b). The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values. interaction as excretion released by cells, changes in the environment and control mechanisms. The parameters b 1 and b 2 take into account these factors also.
Let us assess the cancer cells affecting negatively in the normal cells, by varying b 1 .
In Figure 2 we showT , the solutions of f(T)g(T) = 0. Figure  Let us vary b 1 and compute the corresponding equilibrium valueT 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 non-trivial equilibrium pointQ * . In Figure 3 we show the coordinates of the equilibrium points Q * < (a) andQ * > (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 β 1 = β c 1 = 2.1273 × 10 −2 , we haveC = 0 , at which the big non-trivialQ * > disappears and arises an another equili-briumQ c > , withC = 0 and other coordinates given by equation (9), which has fixed value T > = T c Mathematically, there is another value of b 1 , which does not change the existing equilibrium pointQ c , at which we have T > = T g , that is, Figure 3 shows that trivialQ 0 and small non-trivialQ * < exist for all b 1 , but big nontrivial isQ * > (β 1 < β c 1 ) orQ c > (β 1 ≥ β c 1 ) . In Figure 4, we show the bifurcation diagram, consideringT as a function of b 1 (curve T in (a) and (b) of Figure 3). For all values of b 1 the small equilibrium pointQ * < is the break-point, which divides two regions where trivialQ 0 or big non-trivialQ * > (orQ c > ) is attracting point. Initial conditions set in a small region marked with I and Ia are attracted to the trivial equilibrium pointQ 0 .  However, initial conditions set in regions marked with II and III are attracted to big non-trivial equilibrium pointQ * > for β 1 < β c 1 , while for β 1 ≥ β c 1 (regions marked with IIa and IIIa), toQ c > . Notice that T < always decreases very slowly with b 1 (see also Figure 3(a)), showing that as b 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 b 1 in 1000 times the values of parameters given in Table 2 In Figure 6 we show the coordinates of the equilibrium pointsQ * < (a) andQ * > (b) by varying b 2 . Positive solutions for equation (5) disappeared for higher b 2 .
In Figure 7, we show the bifurcation diagram, consideringT as a function of b 2 (curve T in (a) and (b) of Figure 6). For β 2 < β th 2 , initial conditions set in region marked with I are attracted to the trivial equilibrium pointQ 0 , while for those set in regions II and III are attracted to big non-trivial equilibrium pointQ * > . However, for β 2 > β th 2 all initial conditions are attracted toQ 0 , which is the unique equilibrium, because there is not any positive solution for equation (5). Hence, there is a threshold of the parameter b 2 , denoted by β th 2 , above which all trajectories go to trivial equilibrium. At the threshold value β th 2 = 4.8376 × 10 −2 both roots assume same value, that is, T < (β th 2 ) = T > (β th 2 ) = 0.7053 . For β 2 < β th 2 the dynamical trajectories are similar to that shown in Figure 1.

that negative influence on cancer cells by normal cells affects strongly (very sensitive)
in the cancer dynamics. For β 2 > β th 2 , all initial conditions are attracted to trivial disregarding initial conditions. Figure 8 shows an extreme example, changing only b 2 in 100 times the values of parameters given in Table 2  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 b 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 b 1 higher than its critical value, that is, β 1 > β c 1 . With respect to b 2 , we observed a threshold for β 2 , β th 2 , above which cancer can not be settled. Also, as b 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 (b 1 increases), and decrease in the efficiency of the repairing systems (b 2 decreases).

Recruitment of existing epithelial cells -g
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 g, which is varied.
In Figure 9 we show the coordinates of the equilibrium pointsQ * < (a) andQ * > (b) by varying g. The small T < decreases, while the big T > decreases after an initial increase. Notice that T > decreases due to the fact that pre-existing blood vessels E decreases quickly. The new blood vessels A also decrease. The maximum of the variables occurs at around g = 6.57 × 10 -3 (for C is minimum). In Figure 9(b) we show the curve T E , which depends on g and situates always above the equilibrium value T > . Hence, we have non-trivial equilibrium pointQ * for sufficiently higher values of g.
In Figure 10, we show the bifurcation diagram, consideringT as a function of g (curve T in (a) and (b) of Figure 9, which appear at g = g th ). In (b), a zoom near zero is shown. In Figures 10(a) and 10(b), for g >g th initial conditions set in regions marked with I, Ia and Ib are attracted to the trivial equilibrium pointQ 0 . However, initial conditions set in regions marked with II and III are attracted toQ * > (γ th < γ < γ c 1 ) ; to a limit cycle circulatingQ * > in regions IIa and IIIa (γ c 1 < γ < γ c 2 ) ; and toQ 0 in regions IV and V (γ > γ c 2 ), where E decreases and, then, increases to equilibrium value E 0 .  Figure 9 Positive equilibrium values varying g. The coordinates of the positive equilibrium points varying g. In (a) we show the coordinates of the small equilibrium pointQ * < , and in (b), of the big equilibrium pointQ * > . In (b) we also show the curve of T E , which is always greater than T > . The small root T < decreases, while the big one T > decreases after an initial increase. The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.  The Hopf bifurcation occurs at γ = γ c 1 (supercritical) and γ = γ c 2 (subcritical) [17]. The special values are: g th = 5.450 × 10 -4 and T > (γ th ) = 1.74, γ c 1 = 2.3196 × 10 −2 and T > (γ c 1 ) = 2.1029 , and γ c 2 = 2.5161 × 10 −2 and T > (γ c 2 ) = 1.9443. Figure 9(b) showed that at γ γ c 1 , E is very low. For g < g th (b), all initial conditions set in region marked with Ic are attracted to the trivial equilibrium pointQ 0 , which is the unique equilibrium, because there is not any positive solution for equation (5). Hence, there is a threshold of the parameter g, denoted g 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 g but higher than the threshold (g > g 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 time-delay is encompassed. Our model presented Hopf bifurcation, however only in a range of values of parameter g 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 (g >g 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 g. 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 12, we show the bifurcation diagram, consideringT as a function of ε (curve T in (a) and (b) of Figure 11). When g >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 pointQ 0 , and for initial conditions set in regions II and III, trajectories are attracted toQ * > (a). This behavior results by the existence of influx in equation for A, given by the term gP. Hence, if we let g = 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 g = 0, changing only ε in 100 times the values of parameters given in Table 2 Due to the influx gP, 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 g = 0, there arises a threshold of ε, called ε th , above which new vascularizations promoted by cancer cells can occur. The case g = 0 means that cancer cells are supported exclusively by the new network of blood vessels originated from surrounding tissues for instance, and the pre-existing 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,   Figure 12 Bifurcation diagram varying ε. The bifurcation diagram ofT varying ε, showing the attracting regions. With respect to the coordinates of small equilibrium pointQ * < , the trivialQ 0 is attractor for initial conditions in a very small region I, and for initial conditions in regions II and IIIQ * > is the attractor (a). There is not a critical value due to influx gP in equation for A. However, for g = 0, we have bifurcation diagram similar to the Figure 10.b: for ε <ε th , trivialQ 0 is attractor in I (and for all initial conditions in Ia), andQ * > is attractor in II and III (b) for ε > ε th . The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values. angiogenesis create new vascularization from sprouting originated in existing vasculature, which was called pre-angiogenesis. 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 toQ 0 (a) when T(0) = 9.51 × 10 -3 , while for T(0) = 9.52 × 10 -3 , toQ * > (b). The cancer is triggered at around 400 days. Decreasing k 3 in 50-fold (and some values of other parameters were also changed), the initial value that divides two attracting regions T(0) is around 100-fold 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 b 1 , g 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 (b 1 ); (2) produce higher amount of substances which cause new blood vessels to grow (g); and (3) construct new blood networks to nourish themselves (ε). In opposite way, when the interacting parameter b 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 b 2 and g present threshold values, respectively, β th 2 and g 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 g. The threshold of b 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 g = 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 b 1 and g. Critical value for b 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 g, there are two critical values, named γ c 1 and γ c 2 . When g increases, the epithelial cells E decrease, and damped oscillations appear. When g approaches to γ c 1 , the oscillations are less damped, and when surpasses γ c 1 , regular oscillations occur. However, the amplitude of regular oscillations increases as g approaches to γ c 2 , 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 g surpasses γ c 2 . Again both γ c 1 and γ c 2 do not bear any biological meanings, because the cure of cancer is due to the elimination of pre-existing network of blood vessels. The biological meaningless sustained oscillations begin at γ = γ c 1 (the supercritical Hopf bifurcation) and cease at γ = γ c 2 (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 g increases: (1) damped oscillations aroundQ * for g = 0.75; (2) regular oscillations aroundQ * in the interval [0.76,0.80]; and (3) to the trivial equilibriumQ 0 for g = 0.81. In the region of limit cycle, we observed that: (1) C oscillates between (8.9, 9.0) for all g; E, P and A oscillate between (0, M), where M increases as g increases; and (3) T oscillates between (m,0.1), where m decreases as g 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 g. 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 time-delay is introduced in the tumor growth modeling. They then concluded that an appropriate candidate for describing the cancer growth is the alternative that includes time-delays in the tumor proliferation or angiogenesis process. They also concluded that further mathematical research is warranted for exploring time-delays in the biologically realistic domains in the parameter spaces. Our model, however, showed Hopf bifurcation without time-delay (in fact, there is an elapse of time between pre-angiogenesis 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 non-trivial equilibrium is attractor, the cancer cells reach the level T > . C m increases with increasing b 1 , decreases with b 2 and ε. Cancer cells can grow and reach higher levels when they affect negatively normal cells (b 1 ), but reach lower levels when normal cells acts as a barrier against them (b 2 ). When g varies, T > increases in the initial phase (E decreases), and then decreases (E is practically zero). The parameter g plays an equivalent role of b 1 , but, restricted only to E, which decreases it dramatically. For lower g, there is sufficient number of E to increase A, but epithelial cells are exhausted as g increases, and A decreases (see equation (4)). Finally, T > decreases with ε, which comes out due to relative higher value of g. The behavior of ε is strongly dependent on g due to the influx gP in equation for A: for lower g (also g = 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, toQ 0 (not shown) when T(0) = 2.492 × 10 -1 , while for T (0) = 2.493 × 10 -1 , toQ * > (a); and (2) for δ = 10.0, toQ 0 (not shown) when T(0) = 6.608 × 10 -3 , while for T(0) = 6.609 × 10 -3 , toQ * > (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][21][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 zero-dimensional 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 Lotka-Volterra 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 pre-angiogenesis 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 (g > g th ), and better fitness (increasing capacity of proliferation and capturing nutrients, decreasing mortality, etc.) of cancer cells (b 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 b 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 g 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 g increases, the deviated network plus the new one mounted by the action of angiogenesis effectors nourish the cancer cells. Hence, g = 0 means that pre-existing 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 non-trivial 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 b 1 and g must be restricted. When β 1 > β c 1 we have an equilibrium pointQ c withC = 0 , while for γ c 1 < γ < γ c 2 , we have limit cycle (sustained oscillations) with E~0, but for γ > γ c 2 , 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 g = 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 neo-vascularization. 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 pre-existing 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 g, 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 g i and τ i (τ 0 = 0) are, for i = 0, ···, n, respectively, the i-th 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 g 2i > 0 and g 2i+1 = 0, that is, sproutings occur in the time interval [0, τ 1 ] with rate g 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 i-th 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 Michaelis-Menten [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 i-th concentration of drug administered at time τ i . Another is where q i and τ i are, for i = 0, ···, n, the i-th 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 g and b 2 , named g th and β th 2 , can be managed as follows.
Let us suppose that g and b 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 g th to surpass g, and to decrease β th 2 to fall below b 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 pre-angiogenesis to angiogenesis cells δ and/or increasing the mortality rate μ 4 of pre-angiogenesis cells.

Appendix A: Non-trivial equilibrium point
The non-trivial equilibrium value of cancer cells corresponding to model (1),T , is the positive solution of the equation (5), that is, f(T) = g(T). The fourth degree polynomial and has three non-negative roots: 0, (a 2 -μ 2 )/g = 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 α 1 μ 3 β 1 β 2 k 1 + α 1 − μ 1 β 1 = T g (see equations (7) and (8)). The function g 2 (T) is such that ⎧ ⎨ ⎩ g 2 (−∞) = −∞ g 2 (0) = 0 g 2 (+∞) = +∞, Yang Theoretical Biology and Medical Modelling 2012, 9:2 http://www.tbiomed.com/content/9/1/2 and has two non-negative roots: 0, and 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 b 1 b 2 < a 3 a 1 k 4 /(k 1 k 3 ) (generically referred to as weak interaction between normal and cancer cells when values of b 1 and b 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 b 1 b 2 > a 3 a 1 k 4 /(k 1 k 3 ) (generically referred to as strong interaction between normal and cancer cells when b 1 and b 2 are proportional), in which case T g is lower, we have ⎧ ⎨ ⎩ g(−∞) = −∞ g(0) = g 0 g(+∞) = +∞, 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 b 1 b 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. 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: 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),T must be lower than the lowest root T m of f(T), that is,T ∈ (0, T m ), where T m = min {(a 2 -μ 2 )/g, k 3 }, the minimum between (a 2 -μ 2 )/g and k 3 . Therefore, in the interval (0, T m ) we have at most 2 positive solutions. simplicity say that it is stable) equilibriumQ * > . Between two stable equilibriaQ * > and Q 0 (this always exists and is stable), we have an unstable pointQ * < . We call as the 'break-point' [28,29] the unstable equilibrium pointQ * < , which separates two attracting regions containing one of the equilibrium pointsQ * > andQ 0 . In another words, we have a hyper-surface obeying, e.g., h(C(T < ),Ē(T < ), T < ,P(T < ),Ā(T < )) = 0, generated by the coordinates of the equilibrium pointQ * < such that one of the equilibrium points Q * > andQ 0 is an attractor depending on the relative position of the initial conditions supplied to the dynamical system (1) with respect to the hyper-surface h [30].
The bigger roots of 1 (T) and 2 (T) are not biologically feasible, because they are, respectively, higher than k 3 and (a 2 -μ 2 )/g. Let us define T m = min{T 1 < , T 2 < }, where min{T 1 < , T 2 < } is the minimum between the small roots of, respectively, 1 (T) and 2 (T) given by T 1 < = k 3 1 − 1 − μ 3 k 4 α 3 , assuming that k 4 ≥ 1, and , assuming that T E ≥ T A . Hence, ifT < T m , we have (T) > 0 and there is not a positive number . Notice that, when T E <T A , we have 2 (T) > 0, and it is enough to satisfyT < 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 g assuming higher values.
The following figures are mathematical results, not cancer in an organ.
Dynamical trajectories are shown in Figure 19 for g near γ c 1 : γ 1 = 2.31 × 10 −2 γ 1 c and γ 2 = 2.32 × 10 −2 10 −2 γ c 1 . When T(0) = 0.44, trajectories of both cases go tō Q 0 (trajectories similar to Figure 1(a), not shown). When T(0) = 0.45, dynamical trajectories for g 1 go to stableQ * > (a), while for g 2 , they oscillate around unstableQ * > (b), in which case limit cycle arises. Now, we show the dynamical trajectories for g near γ c 2 . Again, for γ 3 = 2.519 × 10 −2 γ c 2 and γ 4 = 2.520 × 10 −2 γ c 2 , when T(0) = 0.41, trajectories of both cases go toQ 0 (trajectories similar to Figure 1(a), not shown). When T(0) = 0.42, Figure 20 shows dynamical trajectories for g 3 oscillating around unstableQ * > . 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 g 4 and g 5 = 2.519675 × 10 -2 , which is slightly higher than γ c 2 . In both cases we have  Table 2, except g near γ c 1 . When T(0) = 0.45 (for T(0) = 0.44 both cases go toQ 0 ),Q * > is the attracting point for g = 2.31 × 10 -2 (a), and limit cycle with small amplitude circulating unstableQ * > appears for g = 2.32 × 10 -2 (b). The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values. Q * > , and Figure 21 shows trajectories going toQ 0 with different number of oscillations: for g 4 (a), we have one oscillation, while for g 5 , three oscillations, showing C and E (b), T (c), P and A (d). When γ > γ c 2 , 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.
Comparing Figures 19, 20 and 21, we observe that, as g increases, the real part of the complex eigenvalues decreases (see Table 3), and damped oscillations persist for longer times. At γ = γ c 1 , real part is zero. For instance, a pair of complex number has real part -2.3 × 10 -7 at g = 2.3196 × 10 -2 , and +9.0 × 10 -7 at g = 2.3197 × 10 -2 . In both cases, two of them are complex number with negative real part and one negative number. As g increases sinceafter γ c 1 , amplitude of the limit cycle increases, and at γ = γ c 2 disappears. Observe that a finite number of oscillations occurs before reaching the trivial equilibrium. The increasing in the amplitude of regular oscillations resulted in an unsustainable value of T > , and, hence, this value (due to being lower than a critical value) can not trigger new burst of cancer cells. Numerical results regarding to g showed that there is a threshold for g, plus two critical values. For small values, g <g th , cancer cells can not induce new blood vessels, and cancer fades out. As g increases, amplitude of damped oscillations increases and appears stable limit cycle. The limit cycle separates cancer state (Q * ) to cure (Q 0 ). However, the cure occurs at the expense of death of cancer cells due to elimination of the pre-existing network of blood vessels. This phenomenon occurs due to the absence of dependency between normal cells (C) and pre-existing epithelial cells (E) in the carrying capacities k 1 and k 2 .  Table 2, except g near but greater than γ c 2 . When T(0) = 0.42 (Q 0 is attracting for T(0) = 0.41), limit cycle disappears for g = 2.520 × 10 -2 . BeingQ * > unstable, the dynamical trajectories go to trivialQ 0 after one oscillation (a). However, the number of oscillations increases if g is very close to γ c 2 . For g = 2.519675 × 10 -2 and T(0) = 0.42,Q 0 is attained after three oscillations: C and E (b), T (c), P and A (d). The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.