Skip to main content

Mathematical modeling of solid cancer growth with angiogenesis

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 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 α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, k1(E), as k2 depends on C, k2(C), where k1 and k2 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 k3 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 (k4). 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 γ, 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 α3 and ε, which depend respectively on angiogenesis and cancer cells, and carrying capacities k3 and k4, 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 d t C= α 1 C [ 1 - C k 1 ( E ) ] - μ 1 C. 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

d d t C = α 1 C 1 - C k 1 ( E ) - β 1 T k 1 ( E ) - μ 1 C - C m δ D ( t ) d d t E = α 2 E 1 - E k 2 ( C ) - γ E T - μ 2 E d d t T = C m δ D ( t ) + α 3 A T 1 - T k 3 + α 3 T 1 - T k 1 ( E ) - β 2 C k 1 ( E ) - μ 3 T d d t P = γ E T - δ P - μ 4 P d d t A = δ P + ε T A 1 - A k 4 - μ 5 A .

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 k1 and k2 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 α 3 AT ( 1 - T / k 3 ) α 3 T ( 1 - T / k 1 ) . The above impulsive system taking into account these simplifications can be written as, for t > 0,

d d t C = α 1 C 1 - C k 1 - β 1 C T - μ 1 C d d t E = α 2 E 1 - E k 2 - γ E T - μ 2 E d d t T = α 3 A T 1 - T k 3 - β 2 C T - μ 3 T d d t P = γ E T - δ P - μ 4 P d d t A = δ P + ε T A 1 - A k 4 - μ 5 A ,
(1)

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:

( C ( 0 ) = C 0 C m , E ( 0 ) = E 0 , T ( 0 ) = C m , P ( 0 ) = 0 , A ( 0 ) = 0 ) ,
(2)

where C0 and E0 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.

Table 1 Symbols and definitions

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 Q ̄ = ( 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 d t C=0. The first is the equilibrium where there is not any cells, that is,

Q ̄ a b s = ( 0 , 0 , 0 , 0 , 0 ) .

For obvious reason, this point is not considered here.

The trivial equilibrium point, given by

Q ̄ 0 = ( C 0 , E 0 , 0 , 0 , 0 ) ,

where C0 and E0 are

C 0 = k 1 α 1 ( α 1 - μ 1 ) E 0 = k 2 α 2 ( α 2 - μ 2 ) ,
(3)

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 k1(E0) and k2(C0).

The coordinates of the non-trivial equilibrium point Q ̄ * are

C ̄ = k 1 α 1 [ ( α 1 - μ 1 ) - β 1 T ̄ ] Ē = k 2 α 2 [ ( α 2 - μ 2 ) - γ T ̄ ] P ̄ = γ k 2 ( μ 4 + δ ) α 2 [ ( α 2 - μ 2 ) - γ T ̄ ] T ̄ Ā = k 3 k 3 - T ̄ μ 3 α 3 + β 2 k 1 α 3 α 1 ( α 1 - μ 1 ) - β 1 β 2 k 1 α 3 α 1 T ̄ ,
(4)

where T ̄ is the positive solution of the equation

f ( T ) = g ( T ) g 1 ( T ) + g 2 ( T ) ,
(5)

with the fourth degree polynomial f(T), and the third degree polynomials g1(T) and g2(T) being given by

f ( T ) = γ δ k 2 ( μ 4 + δ ) μ 5 α 2 k 3 [ ( α 2 - μ 2 ) - γ T ] × ( k 3 - T ) 2 T g 1 ( T ) = 1 - ε μ 5 T ( k 3 - T ) × μ 3 α 3 + β 2 k 1 α 3 α 1 ( α 1 - μ 1 ) - β 1 β 2 k 1 α 3 α 1 T g 2 ( T ) = ε k 3 μ 5 k 4 T × μ 3 α 3 + β 2 k 1 α 3 α 1 ( α 1 - μ 1 ) - β 1 β 2 k 1 α 3 α 1 T 2 .
(6)

By inspecting C ̄ , Ē and Ā, the positive solution T ̄ , in order to Q ̄ * be biologically feasible, must satisfy, respectively, the constraints T ̄ < T C , T ̄ < T E and T ̄ < T A , where

T C = α 1 - μ 1 β 1 T E = α 2 - μ 2 γ T A = k 3 .
(7)

However, the constraint for the equilibrium Ā is T ̄ < T g , where

T g = α 1 α 3 β 1 β 2 k 1 μ 3 α 3 + β 2 k 1 α 3 α 1 ( α 1 - μ 1 ) = α 1 μ 3 β 1 β 2 k 1 + α 1 - μ 1 β 1 ,
(8)

but it is always satisfied when T ̄ < T C , since T g > T C .

With respect to the constraints, we must have T ̄ < 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 β1 and γ, and directly with the carrying capacity k3. Hence, for higher values of removing parameters β1 and γ, and lower values of carrying capacity k3, the constraints are decreased, that is, the steady state cancer cells T ̄ must assume lower values. Let us suppose that k3 is high, which implies that cancer cells 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, T ̄ must assume lower values, because it must be lower than both T E and k3. Therefore, cancer cells will achieve higher values for high carrying capacity k3, and lower sprouting rate γ. Further angiogenesis cells must be originated by the clonal division from a small amount of pre-angiogenesis cells. In appendix A we show the analysis of equation (5). Summarizing, the number of solutions T ̄ depends on the values assigned to the model's parameters, which can be 0, 2 or 4. Notwithstanding, the solution 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 T ̄ as small (T < , and the corresponding equilibrium point Q ̄ < * ) and big (T>, and the corresponding equilibrium point 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 non-linear 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 T ̄ of equation (5) situates below T C . However, as β1 increases, T ̄ can surpass T C , after intercepting T C at β 1 = β 1 c . As T ̄ increases with β1, C ̄ , given by the first equation of (4), decreases and reaches zero at β 1 = β 1 c , and C ̄ <0 sinceafter. At β 1 = β 1 c we have T ̄ = T C c = ( α 1 - μ 1 ) / β 1 c . Let us suppose that T C c < min { T E , T A } , the minimum between T E and T A . When C ̄ =0, there arises another feasible equilibrium point Q ̄ c , with the coordinates

C ̄ = 0 Ē = k 2 α 2 ( α 2 - μ 2 ) - γ T ̄ P ̄ = γ k 2 ( μ 4 + δ ) α 2 ( α 2 - μ 2 ) - γ T ̄ T ̄ Ā = k 3 μ 3 α 3 ( k 3 - T ̄ ) ,
(9)

where T ̄ = T C c , and the equilibrium T ̄ does not depend anymore on the interacting parameter β1. Hence for β 1 β 1 c we have a constant value T ̄ = ( α 1 - μ 1 ) / β 1 c , which implies in constant values for Ē, P ̄ and Ā.

See more discussion below, in the stability analysis.

The equilibrium point Q ̄ c corresponds to the case where the cancer cells displaced normal cells. By inspecting the equation for 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 Q ̄ c can be avoided as the equilibrium point when k3 assumes small values, that is, if k3 < min{T E , T C }, minimum between T E and T C . In this situation, the equilibrium value T ̄ is always smaller than k3, 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

T ̄ = k 3 Ā α 3 Ā α 3 - μ 3 - β 2 C ̄ .

When β2 increases, T ̄ decreases, and at a certain value of β2, say β 2 t h , T ̄ assumes zero value. When T ̄ =0, the unique feasible equilibrium point is the trivial Q ̄ 0 . When the influence of the normal cells and the environment is higher ( β 2 β 2 t h ), 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 C ̄ =0, the unique feasible equilibrium point is the trivial Q ̄ 0 . Hence, we expect that Ē > 0 for γ ≥ 0, because of T ̄ < T E = ( α 2 - μ 2 ) /γ 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 T ̄ is solution of the polynomial of degree two

0 = 1 - ε μ 5 T ( k 3 - T ) + ε k 3 μ 5 k 4 T × μ 3 α 3 + β 2 k 1 α 3 α 1 ( α 1 - μ 1 ) - β 1 β 2 k 1 α 3 α 1 T ,
(10)

with another being given by T ̄ = T g . Moreover, from equation (4), we have P ̄ =0 (due to γ = 0, the role of the parameter δ does not matter). In this case, there arises another feasible equilibrium point Q ̄ p , with the coordinates

C ̄ = k 1 α 1 [ ( α 1 - μ 1 ) - β 1 T ̄ ] Ē = k 2 α 2 ( α 2 - μ 2 ) P ̄ = 0 Ā = k 3 k 3 - T ̄ μ 3 α 3 + β 2 k 1 α 3 α 1 ( α 1 - μ 1 ) - β 1 β 2 k 1 α 3 α 1 T ̄ . .
(11)

Equation (10) has two positive solutions (not shown) when T g > k3 and ε > μ5/k3, where T g is given by equation (8), while a unique positive solution occurs when T g = k3. Let us introduce a threshold of ε, named εth. We remark that this threshold value appears only for γ = 0. Hence, when ε > εth, with εth> μ5/k3, we have two positive solutions, named T < ε (and the corresponding equilibrium point Q ̄ < p ) and T > ε (and the corresponding equilibrium point Q ̄ > p ), which collapse to one at ε = εth, and for ε < εththere is not positive solution. The equilibrium point 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 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 ( T ( 0 ) > T < ε ) .

Now, for lower values of γ (γ 0), as a consequence of the appearing of εthfor γ = 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 γ < γththe unique solution is T ̄ = 0 (this fact is shown numerically). Note that as ε increases, γthdecreases. 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 T ̄ , which must be substituted in the variables given in equation (4), must be solution of f(T) = g(T), where

f ( T ) = γ δ k 2 ( μ 4 + δ ) μ 5 α 2 k 3 [ ( α 2 - μ 2 ) - γ T ] × ( k 3 - T ) T g ( T ) = μ 3 α 3 + β 2 k 1 α 3 α 1 ( α 1 - μ 1 ) - β 1 β 2 k 1 α 3 α 1 T .
(12)

This equation clearly shows that a positive solution T ̄ is feasible for sufficiently higher values of γ, that is, γ> γ 0 t h , where γ 0 t h corresponds to ε = 0. Due to the fact that ε increases as γthdecreases, γ 0 t h is the upper bound of the threshold of γ. Even the cancer cells do not promote the proliferation of the sprouted endothelial (pre-angiogenesis) 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 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 γ) 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 pre-angiogenesis parameter δ, when δ , Q ̄ appears, similar to the equilibrium point Q ̄ p , because lim δ P ̄ =0. However, from the third equation of (4), we have lim δ δ P ̄ =γĒ T ̄ . Hence, the coordinates of Q ̄ are those given in equation (11), but T ̄ is solution of the equation f(T) = g1(T) + g2(T), where

f ( T ) = γ k 2 μ 5 α 2 k 3 [ ( α 2 - μ 2 ) - γ T ] ( k 3 - T ) 2 T ,

and g1(T) + g2(T) is given be equation (6). Notice that f(T) differs from that in equation (6) by the factor δ/(μ4 + δ). The equilibrium point Q ̄ represents the absence of the intermediate pre-angiogenesis, that is, angiogenesis occurs directly and instantaneously from epithelial cells.

Summarizing, the non-trivial equilibrium points Q ̄ * ( Q ̄ < * and Q ̄ > * , when they exist), Q ̄ c (for β 1 β 1 c ), and Q ̄ p (when γ = 0) are biologically feasible, but Q ̄ c is not a real cancer equilibrium point, since C ̄ =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 β 2 t h and γth(and, eventually, εth, for γ = 0), and a critical value β 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 k1, k2, k3 and k4 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 ( J ̄ - λ I ) =0, where J ̄ is the Jacobian J evaluated at the equilibrium point under analysis,

J ̄ = - α 1 k 1 C ̄ 0 - β 1 C ̄ 0 0 0 - α 2 k 2 Ē - γ Ē 0 0 - β 2 T ̄ 0 - α 3 k 3 Ā T ̄ 0 j 1 0 γ T ̄ γ Ē - j 2 0 0 0 j 3 δ - j 4

where j1, j2, j3 and j4 are

j 1 = α 3 T ̄ 1 - T ̄ k 3 j 2 = μ 4 + δ j 3 = ε Ā 1 - Ā k 4 j 4 = δ P ̄ Ā + ε k 4 Ā T ̄ .
(13)

The Jacobian J ̄ was obtained using equilibrium equations given in (4) assuming that C ̄ 0. When C ̄ =0, the element in the first line and first column is a 11 = ( α 1 - μ 1 ) - β 1 T ̄ , and in the third column, a13 = 0.

We present analytical results corresponding to the equilibrium points Q ̄ 0 and Q ̄ c .

The eigenvalues (from Jacobian J ̄ 0 ) corresponding to the trivial equilibrium point Q ̄ 0 are: λ1 = -μ5, λ2 = - (μ4 + δ), λ3 = - (μ3 + β2C0), λ4 = - (α1 - μ1), and λ5 = - (α2 - μ2). In the last two eigenvalues we used the coordinates given by equation (3). Hence, the trivial equilibrium Q ̄ 0 is locally asymptotically stable whenever this point exists (α1 > μ1 and α2 > μ2), because all eigenvalues become negative.

The non-trivial equilibrium point Q ̄ 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 Jacobian J ̄ c . Once C ̄ =0, one of the eigenvalues is λ 1 = - β 1 ( T ̄ - T C ) , arising one of the conditions to the not real cancer equilibrium Q ̄ 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-matrix J ̄ 1 c , which comes from J ̄ c excluding first row and first column. Let us define matrix A as A = - J ̄ 1 c , or,

A = α 2 k 2 Ē γ Ē 0 0 0 α 3 k 3 Ā T ̄ 0 - j 1 - γ T ̄ - γ Ē j 2 0 0 - j 3 - δ j 4 ,
(14)

where j1, j2, j3 and j4 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 to J ̄ 1 c have negative real part. Hence, Q ̄ c , which is biologically feasible but not real cancer, is locally asymptotically stable only if T ̄ > T C . Therefore, Q ̄ c appears when one of the constraints to Q ̄ * be biologically feasible is violated: in the equation (4), instead of C ̄ <0, we have C ̄ =0 for β 1 β 1 c .

The local stability of the non-trivial equilibrium Q ̄ * is performed numerically. However, we stress the fact that the equilibrium Q ̄ c is an extension of Q ̄ * for β 1 β 1 c . Hence, the big solution T> of equation (5) for β 1 β 1 c must be stable. For this reason we showed the local stability of 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 4thorder 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 point Q ̄ 0 has the cancer free concentrations: C0 = 9 cells/uv and E0 = 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 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 Q ̄ 0 , that is, at time t = 0 the variables assume C(0) = C0 - C m , E(0) = E0, 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) = C0. The initial amount of mutations C m can either grow up to cancer, or fades out.

Table 2 Values of the parameters

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 non-trivial equilibrium points, corresponding to small (T < ) and big (T > ) equilibrium values of cancer cells obtained from equation (5), are, respectively, Q ̄ < * = ( 8 . 92 , 9 . 85 , 0 . 08 , 0 . 07 , 0 . 71 ) and Q ̄ > * = ( 5 . 32 , 2 . 64 , 3 . 68 , 0 . 88 , 1 . 96 ) , 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 Q ̄ 0 is always stable.

Table 3 Eigenvalues

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 Q ̄ 0 ; while for T(0) = 1.024, the dynamical trajectories attain the non-trivial equilibrium 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 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 Q ̄ 0 when T(0) = (1 - ϵ) × 0.08, and to the non-trivial Q ̄ > * 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.

Figure 1
figure 1

Dynamical trajectories using values in 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: trivial Q ̄ 0 for T(0) = 1.023 (a), or non-trivial Q ̄ > * (b). The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.

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 Q ̄ 0 , and not those of the unstable Q ̄ < * . Hence, T(0) is not comparable with T < , one of the coordinates of the break-point Q ̄ < * . The cancer cells proliferate above the subclinical threshold of 103 cells and reach 109 cells which is the X-ray 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 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 b1), 0.1 (b2), 1.0 (b3) and 10.0 (b4). In Figure 2(b) a zoom near lower values of β1 is shown to enhance the small solution T < .

Figure 2
figure 2

Two Positive solutions for four values of β 1 . Two positive roots for the equation f(T) = g(T), for four values of β1. In (a) we show both roots, while in (b), a zoom near lower values of β1. Legends of curves stand for: b1 = 0.01 (Table 2), b2 = 0.1, b3 = 1.0 and b4 = 10.0.

Let us vary β1 and compute the corresponding equilibrium value 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 non-trivial equilibrium point Q ̄ * . In Figure 3 we show the coordinates of the equilibrium points Q ̄ < * (a) and 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 β 1 = β 1 c =2.1273×1 0 - 2 , we have C ̄ =0, at which the big non-trivial Q ̄ > * disappears and arises an another equilibrium Q ̄ > c , with C ̄ =0 and other coordinates given by equation (9), which has fixed value T > = T C c = ( α 1 - μ 1 ) / β 1 c =4.2306 for β 1 β 1 c . Coordinates of Q ̄ > c are same for all β 1 β 1 c . Mathematically, there is another value of β1, which does not change the existing equilibrium point Q ̄ c , at which we have T> = T g , that is, β 1 = β 1 A = 2 . 8001 × 1 0 - 2 . At this value we have Ā = 0, with T > = T g A = α 1 μ 3 β 1 A β 2 k 1 + α 1 - μ 1 β 1 A = 4 . 9996 .

Figure 3
figure 3

Positive equilibrium values varying β 1 . The coordinates of the positive equilibrium points varying β1. In (a) we show the coordinates of the small equilibrium point Q ̄ < * , and in (b), of the big equilibrium point Q ̄ > * . The coordinates of Q ̄ < * is quite insensitive with variation in β1, and in (b) we also show the curves of T C ( C ̄ =0 at T ̄ = T C ) and T g (Ā = 0 at T ̄ = T g ). The small root T<decreases very slowly, while the big one T> increases up to T> = T g , and then, decreases. The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.

Figure 3 shows that trivial Q ̄ 0 and small non-trivial Q ̄ < * exist for all β1, but big non-trivial is Q ̄ > * ( β 1 < β 1 c ) or Q ̄ > c ( β 1 β 1 c ) . In Figure 4, we show the bifurcation diagram, considering T ̄ as a function of β1 (curve T in (a) and (b) of Figure 3). For all values of β1 the small equilibrium point Q ̄ < * is the break-point, which divides two regions where trivial Q ̄ 0 or big non-trivial Q ̄ > * (or Q ̄ > c ) is attracting point. Initial conditions set in a small region marked with I and Ia are attracted to the trivial equilibrium point Q ̄ 0 . However, initial conditions set in regions marked with II and III are attracted to big non-trivial equilibrium point Q ̄ > * for β 1 < β 1 c , while for β 1 β 1 c (regions marked with IIa and IIIa), to 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.

Figure 4
figure 4

Bifurcation diagram varying β 1 . The bifurcation diagram of T ̄ varying β1, showing the attracting regions. With respect to the coordinates of small equilibrium point Q ̄ < * for all β1, the trivial Q ̄ 0 is attractor for initial conditions in a very small region I (and Ia). Initial conditions in regions II and III are attracted to Q ̄ > * for β 1 < β 1 c , and initial conditions in regions IIa and IIIa are attracted to Q ̄ > c for β 1 β 1 c . The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.

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 Q ̄ c , which is biologically feasible, but does not describe evolution of cancer. Dynamical trajectories are attracted to Q ̄ 0 (a) when T(0) = 0.151, while for T(0) = 0.152, to 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 C ̄ (a), or remains in the horizontal axis (b). Increasing β1 in 100-fold, T(0) decreased 7-fold. For another β1 = 100.0, when T(0) = 0.142, trajectories go to Q ̄ 0 ; while for T(0) = 0.143, to Q ̄ > c , which coordinates have exactly the same values found in the previous case (β1 = 10.0). Both cases (figures not shown) correspond to β 1 > β 1 c , hence Q ̄ > c are the same, while Q ̄ < c decrease little bit.

Figure 5
figure 5

Dynamical trajectories using β 1 = 10.0. Dynamical trajectories of the system (1) considering the values of parameters given in Table 2, except β1 = 10.0. The initial conditions determine the region of attraction: trivial Q ̄ 0 for T(0) = 0.151 (a), or non-trivial Q ̄ > * for T(0) = 0.152 (b). The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.

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 Q ̄ < * (a) and Q ̄ > * (b) by varying β2. Positive solutions for equation (5) disappeared for higher β2.

Figure 6
figure 6

Positive equilibrium varying β 2 . The coordinates of the positive equilibrium points varying β2. In (a) we show the coordinates of the small equilibrium point Q ̄ < * , and in (b), of the big equilibrium point Q ̄ > * . The small root T<increases, while the big one T> decreases up to β 2 c , at which T<= T<, and then f(T) = g(T) does not have positive solution. The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.

In Figure 7, we show the bifurcation diagram, considering T ̄ as a function of β2 (curve T in (a) and (b) of Figure 6). For β 2 < β 2 t h , initial conditions set in region marked with I are attracted to the trivial equilibrium point Q ̄ 0 , while for those set in regions II and III are attracted to big non-trivial equilibrium point Q ̄ > * . However, for β 2 > β 2 t h all initial conditions are attracted to 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 β 2 t h , above which all trajectories go to trivial equilibrium. At the threshold value β 2 t h =4.8376×1 0 - 2 both roots assume same value, that is, T < ( β 2 t h ) = T > ( β 2 t h ) = 0 . 7053 .

Figure 7
figure 7

Bifurcation diagram varying β 2 . The bifurcation diagram of T ̄ varying β2, showing the attracting regions. With respect to the coordinates of small equilibrium point Q ̄ < * for β 2 < β 2 t h , the trivial Q ̄ 0 is attractor for initial conditions in region I, and for regions II and III, trajectories are attracted to Q ̄ > * . For β 2 > β 2 t h , all initial conditions (region Ia) are attracted to Q ̄ 0 . The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.

For β 2 < β 2 t h 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 Q ̄ 0 ; while for T(0) = 10.61, to Q ̄ > * . Notice that increasing β2 in 4-fold, T(0) increased 10-fold, showing that negative influence on cancer cells by normal cells affects strongly (very sensitive) in the cancer dynamics. For β 2 > β 2 t h , 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 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) = C0 = 9.0, because the amount C m in the initial condition for T is higher than C0.

Figure 8
figure 8

Dynamical trajectories using β 2 = 0.04. Dynamical trajectories of the system (1) considering the values of parameters given in Table 2, except β2 = 0.04 (near β 2 t h ). The initial conditions determine the region of attraction: trivial Q ̄ 0 for T(0) = 10.60 (a), or non-trivial Q ̄ > * for T(0) = 10.61 (b). The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.

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, β 1 > β 1 c . With respect to β2, we observed a threshold for β 2 , β 2 t h , 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 Q ̄ < * (a) and 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 pre-existing 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 non-trivial equilibrium point Q ̄ * for sufficiently higher values of γ.

Figure 9
figure 9

Positive equilibrium values varying γ. The coordinates of the positive equilibrium points varying γ. In (a) we show the coordinates of the small equilibrium point Q ̄ < * , and in (b), of the big equilibrium point Q ̄ > * . 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.

In Figure 10, we show the bifurcation diagram, considering 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 γ > γthinitial conditions set in regions marked with I, Ia and Ib are attracted to the trivial equilibrium point Q ̄ 0 . However, initial conditions set in regions marked with II and III are attracted to Q ̄ > * ( γ t h < γ < γ 1 c ) ; to a limit cycle circulating Q ̄ > * in regions IIa and IIIa ( γ 1 c < γ < γ 2 c ) ; and to Q ̄ 0 in regions IV and V ( γ > γ 2 c ) , where E decreases and, then, increases to equilibrium value E0. The Hopf bifurcation occurs at γ = γ 1 c (supercritical) and γ= γ 2 c (subcritical) [17]. The special values are: γth= 5.450 × 10-4 and T > ( γ t h ) =1.74, γ 1 c =2.3196×1 0 - 2 and T > ( γ 1 c ) = 2 . 1029 , and γ 2 c =2.5161×1 0 - 2 and T > ( γ 2 c ) = 1 . 9443 . Figure 9(b) showed that at γ γ 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 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.

Figure 10
figure 10

Bifurcation diagram varying γ. The bifurcation diagram of T ̄ varying γ, showing the attracting regions. In (a) we show wide range of variation of γ, and in (b), a zoom near origin. For γ <γth(b) the trivial Q ̄ 0 is attractor for all initial conditions (region Ic). For γ > γth(a), with respect to the coordinates of small equilibrium point Q ̄ < * , the trivial Q ̄ 0 is attractor for regions I, Ia and Ib; and we have three possibilities: (1) for γ t h <γ< γ 1 c , the non-trivial Q ̄ > * is attractor for initial conditions in regions II and III; (2) for γ 1 c <γ< γ 2 c , stable limit cycle circulating unstable Q ̄ > * , in regions IIa and IIIa; and (3) for γ> γ 2 c , trivial Q ̄ 0 is attractor for regions IV and V. In the latter case, the way to reaching the trivial equilibrium is different for initial conditions in Ib and IV or V. The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.

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 time-delay 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 Q ̄ < * (a) and 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.

Figure 11
figure 11

Positive equilibrium values varying ε. The coordinates of the positive equilibrium points varying ε. In (a) we show the coordinates of the small equilibrium point Q ̄ < * , and in (b), of the big equilibrium point Q ̄ > * . The coordinates of Q ̄ < * and Q ̄ > * decrease. The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.

In Figure 12, we show the bifurcation diagram, considering 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 Q ̄ 0 , and for initial conditions set in regions II and III, trajectories are attracted to 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).

Figure 12
figure 12

Bifurcation diagram varying ε. The bifurcation diagram of T ̄ varying ε, showing the attracting regions. With respect to the coordinates of small equilibrium point Q ̄ < * , the trivial Q ̄ 0 is attractor for initial conditions in a very small region I, and for initial conditions in regions II and III Q ̄ > * is the attractor (a). There is not a critical value due to influx γP in equation for A. However, for γ = 0, we have bifurcation diagram similar to the Figure 10.b: for ε <εth, trivial Q ̄ 0 is attractor in I (and for all initial conditions in Ia), and Q ̄ > * 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.

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 Q ̄ 0 (a) when T(0) = 0.416, while for T(0) = 0.417, to Q ̄ > * (b). The dynamical trajectories in (a) are similar to Figure 1(a). Increasing ε in 100-fold, T(0) decreased 2.5-fold.

Figure 13
figure 13

Dynamical trajectories using ε = 1.0. Dynamical trajectories of the system (1) considering the values of parameters given in Table 2, except ε = 1.0. The initial conditions determine the region of attraction: trivial Q ̄ 0 for T(0) = 0.416 (a), or non-trivial Q ̄ > * for T(0) = 0.417 (b). The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.

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 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, 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 k3 and k4 values comparable to k1 and k2 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 k3 and k4 must be lower than k1 and k2, 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 Q ̄ 0 (a) when T(0) = 9.51 × 10-3, while for T(0) = 9.52 × 10-3, to Q ̄ > * (b). The cancer is triggered at around 400 days. Decreasing k3 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 k3 = 5.0, while in Figure 14, near k3 = 0.1. Then all previous simulations can be translated to real world by appropriate scaling factors.

Figure 14
figure 14

Dynamical trajectories using alternative values. Dynamical trajectories of the system (1) considering the fixed values of parameters given in Table 2, except those in the column marked with **. The initial conditions determine the region of attraction: trivial Q ̄ 0 for T(0) = 9.51 × 10-3 (a), or non-trivial Q ̄ > * for T(0) = 9.52 × 10-3 (b). The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.

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, β 2 t h 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 k3 is very low in comparison with constraint T E ). With respect to γ, there are two critical values, named γ 1 c and γ 2 c . When γ increases, the epithelial cells E decrease, and damped oscillations appear. When γ approaches to γ 1 c , the oscillations are less damped, and when surpasses γ 1 c , regular oscillations occur. However, the amplitude of regular oscillations increases as γ approaches to γ 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 γ 2 c . Again both γ 1 c and γ 2 c 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 γ = γ 1 c (the supercritical Hopf bifurcation) and cease at γ= γ 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 Q ̄ * for γ = 0.75; (2) regular oscillations around Q ̄ * in the interval [0.76,0.80]; and (3) to the trivial equilibrium 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 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 β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 Q ̄ 0 (not shown) when T(0) = 2.492 × 10-1, while for T(0) = 2.493 × 10-1, to Q ̄ > * (a); and (2) for δ = 10.0, to Q ̄ 0 (not shown) when T(0) = 6.608 × 10-3, while for T(0) = 6.609 × 10-3, to 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.

Figure 15
figure 15

Dynamical trajectories using δ = 0.001 and δ = 10.0. Dynamical trajectories of the system (1) considering the fixed values of parameters given in Table 2, except those in the column marked with **. The initial conditions determine the region of attraction. For δ = 0.001 non-trivial Q ̄ * is attractor when T(0) = 2.493 × 10-1 (a), while for δ = 10.0, when T(0) = 6.609 × 10-3 (b). The scales of vertical and horizontal axes must be multiplied by the factors shown in the legends to obtain the actual values.

Conclusions

Many models [5] [19] have already been proposed to describe cancer growth, and some of those models have explicitly considered the spatial dimension [2022], 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 = t0 = 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 (γ > γ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 k1 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 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 β1 and γ must be restricted. When β 1 > β 1 c we have an equilibrium point Q ̄ c with C ̄ =0, while for γ 1 c <γ< γ 2 c , we have limit cycle (sustained oscillations) with E ~ 0, but for γ> γ 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 time-delays 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 γ, 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

d d t E = α 2 E 1 - E k 2 - Φ ( t ) = μ 2 E d d t P = Φ ( t ) - δ P - μ 4 P ,

where Φ(t) is the total intermittent sproutings rate. One form of Φ(t) is

Φ ( t ) = i = 0 n γ i θ ( t - τ i ) θ ( τ i + 1 - t ) E T ,

where γ 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 γ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

Φ ( t ) = i = 0 n s i δ ( t - τ i ) E ,

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,

d d t T = α 3 A T 1 - T k 3 - β 2 C T - μ 3 T - μ Q T Q k q + Q ,

and adding an equation for the drug administration as

d d t Q = q ( t ) - l μ Q T Q k q + Q - μ q Q ,

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

q ( t ) = i = 0 n u i δ ( t - τ i ) ,

where u i is, for i = 0, ···, n, the i-th concentration of drug administered at time τ i . Another is

q ( t ) = i = 0 n q i θ ( t - τ i ) θ ( τ i + 1 - t ) ,

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 q2i> 0 and q2i+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 γthand β 2 t h , 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 γthto surpass γ, and to decrease β 2 t h 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 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 f(T) is such that

f ( - ) = - f ( 0 ) = 0 f ( + ) = - ,

and has three non-negative roots: 0, (α2 - μ2) = T E , and k3 = 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.

Figure 16
figure 16

Qualitative behavior of f ( T ). Qualitative behavior of f(T), with k3 <T E (a) and k3 > T E (b).

The function g1(T) is such that

g 1 ( - ) = + g 1 ( 0 ) = g 0 = k 3 μ 3 α 3 + β 2 k 1 α 3 α 1 ( α 1 - μ 1 ) > 0 g 1 ( + ) = - ,
(A.1)

and has three positive roots: ε/μ5, k3 = T A , and α 1 μ 3 β 1 β 2 k 1 + α 1 - μ 1 β 1 = T g (see equations (7) and (8)). The function g2 (T) is such that

g 2 ( - ) = - g 2 ( 0 ) = 0 g 2 ( + ) = + ,

and has two non-negative roots: 0, and α 1 μ 3 β 1 β 2 k 1 + α 1 - μ 1 β 1 = T g (double root). The dominant terms of the third degree polynomials g1(T) and g2(T), when T± ∞, are

g 1 ( T ) = v 1 T 3 ; w i t h v 1 = - ε β 1 β 2 k 1 μ 5 α 3 α 1 g 2 ( T ) = v 2 T 3 ; w i t h v 2 ε ( β 1 β 2 k 1 ) 2 k 3 μ 5 ( α 3 α 1 ) 2 k 4 .

The third degree polynomial g(T) is sum of g1(T) and g2(T). Notice that, for T ≥ 0, we have g2(T) ≥ 0, while g1(T) changes signal. Hence, the effect of g2(T) in the sum (g(T)) is the increasing in the values of g1(T), except at T = T g , at which we have g1(T g ) = g2(T g ). The asymptotic behavior of the polynomial g(T) depends on the velocities of increasing v2 and v2 (or coefficients) of T3. When β1 β2 < α3 α1 k4/(k1 k3) (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

g ( - ) = + g ( 0 ) = g 0 g ( + ) = - .

It can be shown that there is one solution or three positive solutions. When β1 β2 > α3 α1 k4/(k1 k3) (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

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 g0 is given by equation (A.1). The existence of other positive solutions depends on the relative position of the roots of g1(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 g1(T), with three positive solutions (a), and one solution (b); when T g is between the roots ε/μ5 and k3, with three positive solutions (c); and when T g is the smallest root of g1(T), with three positive solutions (d). T1 and T2, the roots of g1(T), stand for ε/μ5 and k3, depending on the relative positions between them. Figure 17(a) and 17(b) shows clearly the effect of g2(T) increasing g1(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 g2(T) is sufficiently higher.

Figure 17
figure 17

Qualitative behavior of g ( T ). Qualitative behavior of g(T), for small β1 β2 (weak interaction): T g is the greatest root of g1(T), with three positive solutions (a), and one solution (b); T g is between the roots ε/μ5 and k3, with three positive solutions (c); and T g is the smallest root of g1(T), with three positive solutions (d). T1 and T2, the roots of g1(T), stand for ε/μ5 and k3.

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. 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) = -g0 < 0.

  2. 2.

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

  3. 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. 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 {(α2 - μ2)/γ, k3}, the minimum between (α2 - μ2)/γ and k3. 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 T ̄ < T A , T ̄ < T E and T ̄ < T g . Additionally, the constraint T ̄ < T C are satisfied, once 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 = k3.

Figure 18
figure 18

Qualitative behavior of f ( T ) - g ( T ). Qualitative behavior of f(T) - g(T): none (a), two (b) and four (c) solutions. In (d), we show two positive solutions not biologically feasible, because they are greater than the constraint T A = k3.

When there is not any positive solution for f(T) - g(T) = 0, the unique equilibrium point is the trivial Q ̄ 0 . In the case of two positive solutions, the small root (T < ) forms the unstable equilibrium Q ̄ < * , and the big one (T>) forms a possibly stable (let us for simplicity say that it is stable) equilibrium Q ̄ > * . Between two stable equilibria Q ̄ > * and Q ̄ 0 (this always exists and is stable), we have an unstable point Q ̄ < * . We call as the 'break-point' [28, 29] the unstable equilibrium point Q ̄ < * , which separates two attracting regions containing one of the equilibrium points Q ̄ > * and Q ̄ 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 point Q ̄ < * such that one of the equilibrium points Q ̄ > * and Q ̄ 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].

Appendix B: The local stability of the non-trivial equilibrium Q ̄ c

Let us show that the matrix A, given by equation (14), is an M-matrix [31] for big non-trivial equilibrium Q ̄ > c , where T> is one of the coordinates.

Definition. We say that the n × n matrix A = [a ij ] is a non-singular M-matrix if a ij ≤ 0, ij, and there exists a matrix B ≥ 0 and a real number s > 0 such that

A = s I B a n d s > ( ρ ( B ) ,

where I is the identity matrix and ρ is the spectral radius [32].

Or, equivalently:

Proposition 1. A is a non-singular M-matrix if and only if the real part of its eigenvalues is greater than zero.

Proposition 2. A (elements a ij ) is a non-singular M-matrix 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,

a i i d i > j i a i j d j ,

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

α 2 k 2 Ē d 1 > γ Ē d 2 α 3 k 3 Ā T ̄ d 2 > α 3 k 3 T ̄ ( k 3 - T ̄ ) d 1 ( μ 4 + δ ) d 3 > γ T d 1 + γ Ē d 2 δ P ̄ Ā + ε k 4 Ā T ̄ d 4 > ε k 4 Ā k 4 - Ā d 2 + δ d 3 ,
(B.1)

because Ā can be greater than k4. The equilibrium values correspond to the point Q ̄ c , given by equation (9).

Let us define

d 1 = γ k 2 + φ α 2 d 2 = 1 d 3 = γ 2 k 2 T ̄ + γ α 2 Ē + ( γ T ̄ α 2 ) φ α 2 ( μ 4 + δ ) d 4 = Ā - φ k 3 - T ̄ ,

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

0 < φ < φ n φ d ,
(B.2)

where the numerator φ n and denominator φ d > 0 are

φ n = δ P ̄ Ā + ε k 4 Ā T ̄ Ā k 3 - T ̄ - ε k 4 Ā k 4 - Ā - δ γ k 2 ( α 2 - μ 2 ) α 2 ( μ 4 + δ ) φ d = δ P ̄ Ā + ε k 4 Ā T ̄ 1 k 3 - T ̄ + δ ( γ T ̄ + α 2 ) α 2 ( μ 4 + δ ) ,

with the last term of the numerator being obtained using the relation

γ 2 k 2 T ̄ + γ α 2 Ē = γ k 2 ( α 2 - μ 2 ) ,

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 Q ̄ > c , with coordinate T> as solution of (5), is stable. Additionally, we assume that Ā > k4, which seems reasonable since the equilibrium Q ̄ c is stable for T ̄ > ( α 1 - μ 1 ) / β 1 = T C . Substituting the coordinates of the equilibrium Q ̄ c , the numerator of φ can be written as

φ n = Ā ( k 3 - T ̄ ) 2 Φ ( T ̄ ) ,

with

Φ ( T ̄ ) = ε k 4 α 3 Φ 1 ( T ̄ ) - γ δ k 2 α 3 α 2 ( μ 4 + δ ) μ 3 k 3 × k 3 - T ̄ 2 Φ 2 ( T ̄ ) ,
(B.3)

where

Φ 1 ( T ̄ ) = k 4 α 3 T ̄ 2 - 2 k 3 ( k 4 α 3 - μ 3 ) T ̄ + k 3 2 ( k 4 α 3 - μ 3 ) Φ 2 ( T ̄ ) = γ T ̄ 2 - 2 ( α 2 - μ 2 ) T ̄ + ( α 2 - μ 2 ) ,

and we must show that Φ ( T ̄ ) >0, for higher value of T ̄ , that is, T ̄ = T > . Let us assume that k4α3 > μ3. Then, we have Φ ( T ̄ ) >0, for T ̄ 0, because the discriminant of Φ 1 ( T ̄ ) is

Δ 1 = - 4 k 3 2 μ 3 ( k 4 α 3 - μ 3 ) < 0 .

Hence, Φ 2 ( T ̄ ) determines the existence of positive φ n . The discriminant of Φ 2 ( T ̄ ) is

Δ 2 = 4 γ ( α 2 - μ 2 ) ( T E - T A ) ,

where T E and T A are given in equation (7). We have two possibilities. First, when T E < T A , we have Φ 2 ( T ̄ ) >0, for all T ̄ 0, and Φ ( T ̄ ) >0 for

ε > ε min ,

where

ε min = γ δ k 2 α 3 k 4 α 3 α 2 ( μ 4 + δ ) μ 3 k 3 ( k 3 - T > ) 2 Φ 2 ( T > ) Φ 1 ( T > ) .
(B.4)

Hence, when ε > εmin, φ n > 0 and we have a positive number φ obeying equation (B.2). Second, when T E > T A , Φ 2 ( T ̄ ) has two positive solutions Φ2<and Φ2>given by

Φ 2 < = α 2 - μ 2 γ - α 2 - μ 2 γ α 2 - μ 2 γ - k 3 = T E - T E ( T E - T A ) Φ 2 > = α 2 - μ 2 γ + α 2 - μ 2 γ α 2 - μ 2 γ - k 3 = T E + T E ( T E - T A ) ,

and we have Φ 2 ( T ̄ ) 0 , for Φ 2 < T ̄ Φ 2 > ; otherwise, Φ 2 ( T ̄ ) >0. Notice that Φ2> is out of the range of feasibility, once Φ2> >T E . We know that ( α 1 - μ 1 ) / β 1 < T ̄ < k 3 , but the bigger solution under consideration is T > = ( α 1 - μ 1 ) / β 1 c = T E c , for β 1 β 1 c . It is easy to show that Φ2(k3) < 0; but, Φ 2 ( T E c ) 0, if T E c Φ 2 < , and Φ 2 ( T E c ) >0, if T E c < Φ 2 < . Hence, Φ 2 ( T ̄ ) 0 if T E c < Φ 2 < . Hence, Φ ( T ̄ ) >0, for T ̄ 0, when T E c Φ 2 < , or when T E c Φ 2 < and T ̄ Φ 2 < ; and Φ ( T ̄ ) >0, when T E c < Φ 2 < and T E c < T ̄ < Φ 2 < for ε > εmin.

Summarizing, φ n > 0 occurs, in order to have positive number φ:

1. T A < T E or γ < (α2 - μ2)/k3 - weak capacity of recruitment of the normal epithelial cells by cancer cells. We have:1.a If T E c Φ 2 < : φ n > 0 without restriction about T>.

1.b If T E c < T ̄ < Φ 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.

  1. 2.

    T A > T E or γ > (α2 - μ2)/k3 - strong capacity of coopting normal epithelial cells by cancer cells: φ n > 0, if ε > εmin, where T> satisfies εmin.

When γ is small, the big equilibrium 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 Q ̄ > c can be unstable (case 2), in which case Hopf bifurcation can occur (see Appendix C).

In the case of the small equilibrium Q ̄ < c , with one coordinate T<, it is unstable. Assuming that Ā < k4, we show here that A corresponding to small T< is not an M-matrix (we are not proving that Q ̄ < c is unstable). In this case, equation (B.3) becomes

Φ ( T ̄ ) = - ε k 4 α 3 Φ 1 ( T ̄ ) - γ δ k 2 α 3 α 2 ( μ 4 + δ ) μ 3 k 3 × ( k 3 - T ̄ ) 2 Φ 2 ( T ̄ ) ,

where

Φ 1 ( T ̄ ) = k 4 α 3 ( T ̄ - k 3 ) 2 - μ 3 k 3 2 Φ 2 ( T ̄ ) = γ T ̄ 2 - 2 ( α 2 - μ 2 ) T ̄ + ( α 2 - μ 2 ) .

The bigger roots of Φ 1 ( T ̄ ) and Φ 2 ( T ̄ ) are not biologically feasible, because they are, respectively, higher than k3 and (α2 - μ2)/γ. 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 k4 ≥ 1, and T < 2 = T E 1 - 1 - T A T E , assuming that T E T A . Hence, if T ̄ < 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 satisfy 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.

Figure 19
figure 19

Dynamical trajectories using γ near γ 1 c Dynamical trajectories of the system (1) considering the values of parameters given in Table 2, except γ near γ 1 c . When T(0) = 0.45 (for T(0) = 0.44 both cases go to Q ̄ 0 ), Q ̄ > * is the attracting point for γ = 2.31 × 10-2 (a), and limit cycle with small amplitude circulating unstable Q ̄ > * appears for γ = 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.

Dynamical trajectories are shown in Figure 19 for γ near γ 1 c : γ 1 = 2 . 31 × 1 0 - 2 γ c 1 and γ 2 = 2 . 32 × 1 0 - 2 1 0 - 2 γ 1 c . When T(0) = 0.44, trajectories of both cases go to Q ̄ 0 (trajectories similar to Figure 1(a), not shown). When T(0) = 0.45, dynamical trajectories for γ1 go to stable Q ̄ > * (a), while for γ2, they oscillate around unstable Q ̄ > * (b), in which case limit cycle arises.

Now, we show the dynamical trajectories for γ near γ 2 c . Again, for γ 3 = 2 . 519 × 1 0 - 2 γ 2 c and γ 4 = 2 . 520 × 1 0 - 2 γ 2 c , when T(0) = 0.41, trajectories of both cases go to 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 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 γ 2 c . In both cases we have γ> γ 2 c , and Figure 21 shows trajectories going to 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 γ> γ 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.

Figure 20
figure 20

Dynamical trajectories using γ near but lower than γ 2 c Dynamical trajectories of the system (1) considering the values of parameters given in Table 2, except γ near but lower than γ 2 c . When T(0) = 0.42 ( Q ̄ 0 is attracting for T(0) = 0.41), limit cycle with large amplitude circulating unstable Q ̄ > * occurs for γ = 2.519 × 10-2. Regular oscillations are observed in all variables: C and E (a), T (b), P (c) 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.

Figure 21
figure 21

Dynamical trajectories using γ near but greater than γ 2 c Dynamical trajectories of the system (1) considering the values of parameters given in Table 2, except γ near but greater than γ 2 c . When T(0) = 0.42 ( Q ̄ 0 is attracting for T(0) = 0.41), limit cycle disappears for γ = 2.520 × 10-2. Being Q ̄ > * unstable, the dynamical trajectories go to trivial Q ̄ 0 after one oscillation (a). However, the number of oscillations increases if γ is very close to γ 2 c . For γ = 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.

Comparing Figures 19, 20 and 21, we observe that, as γ increases, the real part of the complex eigenvalues decreases (see Table 3), and damped oscillations persist for longer times. At γ = γ 1 c , real part is zero. For instance, a pair of complex number has real part -2.3 × 10-7 at γ = 2.3196 × 10-2, and +9.0 × 10-7 at γ = 2.3197 × 10-2. In both cases, two of them are complex number with negative real part and one negative number. As γ increases sinceafter γ 1 c , amplitude of the limit cycle increases, and at γ= γ 2 c 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 γ showed that there is a threshold for γ, plus two critical values. For small values, γ < γth, cancer cells can not induce new blood vessels, and cancer fades out. As γ 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 k1 and k2.

References

  1. Jemal A, Siegel R, XJ Q, Ward E: Cancer Statistics, 2010. CA Cancer J Clin. 2010, 60: 277-300. 10.3322/caac.20073.

    Article  PubMed  Google Scholar 

  2. Parkin DM, Pisani PJF: Global Cancer Statistics. CA Cancer J Clin. 1999, 49: 33-64. 10.3322/canjclin.49.1.33.

    Article  CAS  PubMed  Google Scholar 

  3. Sompayrac L: How Cancer Works. 2004, Boston: Jones and Bartlett Publishers

    Google Scholar 

  4. Byrne HM: Using mathematics to study solid tumour growth. Proceedings of the 9th General Meetings of European Women in Mathematics. 1999, 81-107.

    Google Scholar 

  5. Araujo RP, McElwain DLS: A History of the Study of Solid Tumour Growth: The Contribution of Mathematical Modelling. Bull Math Biol. 2004, 66: 1039-1091. 10.1016/j.bulm.2003.11.002.

    Article  CAS  PubMed  Google Scholar 

  6. Michelson S, Leith JT: Positive Feedback and Angiogenesis in Tumor Growth Control. Bull Math Biol. 1997, 59: 233-254. 10.1007/BF02462002.

    Article  CAS  PubMed  Google Scholar 

  7. d'Onofrio A, Gandolfi A: Tumour Eradication by Antiangiogenic Therapy: Analysis and Extensions of the Model by Hahnfeldt et al. (1999). Math Biosc. 2004, 191: 159-184. 10.1016/j.mbs.2004.06.003.

    Article  Google Scholar 

  8. Risau W: Mechanisms of Angiogenesis. Nature. 1997, 386: 671-674.

    Article  CAS  PubMed  Google Scholar 

  9. Chaplain MAJ, McDougall SR, Anderson ARA: Mathematical Modeling of Tumor-Induced Angiogenesis. Annu Rev Biomed Eng. 2006, 8: 233-257. 10.1146/annurev.bioeng.8.061505.095807.

    Article  CAS  PubMed  Google Scholar 

  10. Edelstein-Keshet L: Mathematical Models in Biology. 1988, New York: McGraw Hill, Inc

    Google Scholar 

  11. Nani F, Freedman HI: A mathematical Model of Cancer Treatment by Immunotherapy. Math Biosc. 2000, 163: 159-199. 10.1016/S0025-5564(99)00058-9.

    Article  CAS  Google Scholar 

  12. Ruggiero RA, Bustoabad OD: The Biological Sense of Cancer: A Hypothesis. Theoret Biol Med Modelling. 2006, 3: 43:1-14.

    Article  Google Scholar 

  13. Kitagawa M, Utsuyama M, Kurata M, Yamamoto K, Yuasa Y, Ishikawa Y, Arai T, Hirokawa K: Cancer and Aging: Symposium of the 27th Annual Meeting of the Japanese Society for Biomedical Gerontology. Tokyo Cancer Immunol Immunother. 2005, 54: 623-634. 10.1007/s00262-004-0622-9.

    Article  PubMed  Google Scholar 

  14. Komarova NL, Mironov V: On The Role of Endothelial Progenitor Cells in Tumor Neovascularization. J Theoret Biol. 2005, 235: 338-349. 10.1016/j.jtbi.2005.01.014.

    Article  CAS  Google Scholar 

  15. Murray JD: Mathematical Biology. 1989, New York: Springer-Verlag

    Book  Google Scholar 

  16. Press WH, Flannery BP, Teukolsky SA, Vetterling WT: Numerical Recipes The Arts of Scientifc Computing (FORTRAN Version). 1989, Cambridge: Canbridge University Press

    Google Scholar 

  17. Kuznetsov YA: Elements of Applied Bifurcation Theory. 1995, New York: Springer-Verlag

    Book  Google Scholar 

  18. Agur Z, Larakelyan L, Daugulis P, Ginosar Y: Hopf Point Analysis for Angiogenesis Model. Discr Contin Dynam Syst. 2004, 4 (1): 29-38.

    Article  Google Scholar 

  19. Peirce SM: Computational and Mathematical Modeling of Angiogenesis. Microcirculation. 2008, 15 (8): 739-751. 10.1080/10739680802220331.

    Article  PubMed Central  PubMed  Google Scholar 

  20. Chaplain MAJ: Mathematical Modelling of Angiogenesis. J Neuro-Oncology. 2000, 50: 37-51. 10.1023/A:1006446020377.

    Article  CAS  Google Scholar 

  21. Ribba B, Colin T, Schnell S: A Multiscale Mathematical Model of Cancer, And Its Use in Analyzing Irradiation Therapies. Theoret Biol Med Modelling. 2006, 3: 7:1-19.

    Article  Google Scholar 

  22. Budu-Grajdeanu P, Schugart RC, Friedman A, Valentine C, Agarwal AK, Rovin BH: A Mathematical Model of Venous Neointimal Hyperplasia Formation. Theoret Biol Med Modelling. 2008, 5: 2:1-9.

    Article  Google Scholar 

  23. Arakelyan L, Vainstein V, Agur Z: A Computer Algorithm Describing The Process of Vessel Formation and MAturation, and Its Use for Predicting The Effects of Anti-angiogenic and Anti-maturation Therapy on Vascular Tumor Growth. Angiogenesis. 2002, 5: 203-214. 10.1023/A:1023841921971.

    Article  CAS  PubMed  Google Scholar 

  24. Stamatakos GS, Kolokotroni EA, Dionysiou DD, Geordiadi EC, Desmedt C: An Advanced Discrete State-Discrete Event Multiscale Simulation Model of The Response of A Solid Tumor to Chemotherapy: Mimicking a Clinical Study. J Theoret Biol. 2010, 266: 124-139. 10.1016/j.jtbi.2010.05.019.

    Article  CAS  Google Scholar 

  25. Yang HM: Modeling Directly Transmitted Infections in a Routinely Vaccinated Population - The Force of Infection Described by Volterra Integral Equation. Appl Math Comput. 2001, 122 (1): 27-58. 10.1016/S0096-3003(00)00011-4.

    Article  Google Scholar 

  26. Yang HM: Modelling Vaccination Strategy Against Directly Transmitted Diseases Using a Series of Pulses. J Biol Syst. 1998, 6 (2): 187-212. 10.1142/S0218339098000145.

    Article  Google Scholar 

  27. Kenner J, Sneyd J: Mathematical Physiology. 1998, New York: Springer

    Google Scholar 

  28. Bradley DJ, May RM: Consequences of Helminth Aggregation for the Dynamics of Schistosomiasis. Trans R Soc Trop Med Hyg. 1978, 73: 262-273.

    Article  Google Scholar 

  29. May RM: Togetherness Amongst Schistosome: Its Effects on the Dynamics of the Infection. Math Biosc. 1977, 35: 301-343. 10.1016/0025-5564(77)90030-X.

    Article  Google Scholar 

  30. Esteva L, Yang HM: Mathematical Model to Assess the Control of Aedes aegypti Mosquitoes by The Sterile Insect Technique. Math Biosc. 2005, 198: 132-147. 10.1016/j.mbs.2005.06.004.

    Article  Google Scholar 

  31. Raimundo SM, Massad E, Yang HM: Modelling Congenital Transmission of Chagas' Disease. BioSystems. 2010, 99: 215-222. 10.1016/j.biosystems.2009.11.005.

    Article  PubMed  Google Scholar 

  32. Berman A, Plemmons RJ: Nonnegtive Matrices in the Mathematical Sciences. 1979, New York: Academic Press

    Google Scholar 

Download references

Acknowledgements

This work was supported by a grant from FAPESP (Projeto Temático). HMY gratefully acknowledges a Fellowship awarded by CNPq. HMY thanks PFA Mancera for careful reading of the manuscript and providing corrections.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Hyun M Yang.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Authors’ original file for figure 3

Authors’ original file for figure 4

Authors’ original file for figure 5

Authors’ original file for figure 6

Authors’ original file for figure 7

Authors’ original file for figure 8

Authors’ original file for figure 9

Authors’ original file for figure 10

Authors’ original file for figure 11

Authors’ original file for figure 12

Authors’ original file for figure 13

Authors’ original file for figure 14

Authors’ original file for figure 15

Authors’ original file for figure 16

Authors’ original file for figure 17

Authors’ original file for figure 18

Authors’ original file for figure 19

Authors’ original file for figure 20

Authors’ original file for figure 21

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Yang, H.M. Mathematical modeling of solid cancer growth with angiogenesis. Theor Biol Med Model 9, 2 (2012). https://doi.org/10.1186/1742-4682-9-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1742-4682-9-2

Keywords