Critical point analysis for three-variable cancer angiogenesis

We perform critical-point analysis for three-variable systems that represent essential processes of the growth of the angiogenic tumor, namely, tumor growth, vascularization, and generation of angiogenic factor (protein) as a function of effective vessel density. Two models that describe tumor growth depending on vascular mass and regulation of new vessel formation through a key angiogenic factor are explored. The first model is formulated in terms of ODEs, while the second assumes delays in this regulation, thus leading to a system of DDEs. In both models, the only nontrivial critical point is always unstable, while one of the trivial critical points is always stable. The models predict unlimited growth, if the initial condition is close enough to the nontrivial critical point, and this growth may be characterized by oscillations in tumor and vascular mass. We suggest that angiogenesis per se does not suffice for explaining the observed stabilization of vascular tumor size.

1. Introduction.The crucial role of blood supply in the tumor development is presently well recognized (see, e.g., [11,16,27]).Most solid tumors start their growth in avascular form, where the cells obtain oxygen, nutrients, and growth factors by diffusion from the outer environment.Such a form has a size limitation, shown both experimentally (e.g., [29,34]) and theoretically (e.g., [2,1,19,12,35]).To continue its growth, the tumor needs to initiate the process of angiogenesisformation of new blood vessels inside the tissue.This newly created vascular network enables tumor expansion and growth by carrying all the needed nutrients into the tumor tissue.Angiogenesis normally occurs in healthy tissues, for example, in wound healing.It consists of several stages, from endothelial cells proliferation and migration to establishment of the full vascular network of mature vessels covered by the smooth muscle cells.This process is regulated by specific factors that are secreted by the cells in the tissue.The regulation of angiogenesis is complex and includes several feedback processes balancing each other at different stages.During tumor development, angiogenesis is initiated by factors secreted by tumor cells, because of nutrient deprivation [15,30].Ineffective angiogenesis can lead to insufficient blood supply and thus limit tumor growth and invasion.In recent decades, different models of tumor angiogenesis and tumor growth dependence on angiogenesis has been proposed.Those models concentrate on different key aspects of the angiogenic process.Some models focus on the formation of new vasculature in the presence of stimulating factors (e.g., [13]).Others define the whole process of tumor growth as depending on the vascular mass (e.g., [9,31,24]).More detailed models can describe some specific aspects of the interplay between tumor dependence on the vasculature and the regulation of angiogenesis by normal and cancerous cells.
In [3], a family of models of tumor angiogenesis was presented.These models are new in assuming that the effective vessel density (rather than vessel mass per se, and vessel maturation) may be significant in modulating vascular tumor dynamics.The work was aimed at explaining the experimentally observed oscillations in tumor growth ( [4,22]) and relating these oscillations to the intrinsic properties of the angiogenic process.It was shown that in all the proposed models, Hopf points exist whenever a delay is assumed in the regulation of tumor proliferation and vessels formation.The existence of Hopf points is associated with periodic behavior in the system-oscillations of tumor size and vessel volume.Time delays and their role in tumor growth and oscillations were also considered in connection with other processes, for example, proliferation and apoptosis of cells (cf.[7,8,10,18,20]).
In the present work, we further investigate the general behavior of two models considered in [3] with and without delay.These models assume that tumor dynamics are controlled by angiogenesis and, in particular, by the effective vessel density.This assumption will be examined by analyzing the models and comparing them to real-life dynamics.To this end, we examine stability and asymptotic behavior of solutions.This integrates with the previous work to yield effective characterization of those models.In the future, we will extend this investigation to more complicated models.
2. Model description.We present short description of the models under discussion.For a more detailed description see [3].We consider a system including three time-dependent variables: tumor size, N (t); amount of regulating protein, P (t); and vessel volume V (t).
The process is described by a system of ODEs (with or without delays), and the interdependence between variables is described using "sigmoid-shaped" functions, assumed to be smooth.Such functions are used to represent, in a generic way, the response of biological systems to different stimuli.Their form reflects experimental observations that the system reactions change only within a certain range of stimuli values.Inside this range, the dependence is monotonic, but possibly with a nonconstant derivative.We do not propose any specific formulae for those functions; rather, we describe their crucial properties.We also use the variable for effective vessel density, which is calculated from tumor size and vessel volume: E = V /N (see [3,5]).
The equations for the first model are In (1), tumor growth rate is described by f 1 , which depends on vessel density, E. It is an increasing sigmoid function with f 1 (0) = −a 1 < 0 (when the vessel density is insufficient, tumor cells die) and lim x→∞ f 1 (x) = b 1 > 0 (maximal possible growth rate is limited by maximal proliferation rate of the cells).The rate of protein secretion by tumor cells in (1), which depends on vessel density is described by f 2 .This function satisfies f 2 (0) = a 2 > 0, lim x→∞ f 2 (x) = 0, and decreases to zero convexly (the secretion is the response to nutrient deprivation that results from reduced vessel density).The growth rate of vasculature regulated by the protein amount is described by f 3 .This function has properties similar to f 1 : f 3 (0) = −a 3 < 0 (so vessels regress if not stimulated) and lim x→∞ f 3 (x) = b 3 > 0, and the function is increasing.
Another form of the same system is obtained by a substitution V → E: The second model uses similar equations but assumes the existence of time delays in the growth rates' dependence on E and P : 3. Three-variable model without delay.In this section, we study system (2).We focus on the analysis of (2), instead of (1) in original variables, for two main reasons.The first one is of biological nature.We focus on the importance of effective vessel density in the process of angiogenesis.Therefore, the properties of (2) are of our main interest.On the other hand, the system after substitution has better mathematical structure-it is naturally defined for every nonnegative values of variables N , P , E, while (1) is not defined for N = 0 because of the definition of E. Hence, we study the behavior of solutions to (2) and then interpret it in the original variables N , P , V .
Remark 1. Assume that f i (x), i = 1, 2, 3 are of class C 1 for nonnegative x.Then for every nonnegative initial values N 0 , P 0 , E 0 ≥ 0 and every t ≥ 0, there exists a unique nonnegative solution to (2).
Proof.Local existence and uniqueness of solutions to (2) is obvious.We show that such a solution is nonnegative.For the first and the third equations of (2), nonnegativity is an easy consequence of the form of these equations.Hence, for the second equation, we have Ṗ ≥ −δP, which implies P (t) ≥ P 0 e −δt ≥ 0. Global existence of solutions is also an easy consequence of the form of the right-hand side.Specifically, it is linearly bounded, which guarantees existence of solutions for every t ≥ 0.

Critical points.
In this section, we study the existence and local stability of the critical points in system (2).
Obviously, there exists at least one critical point (0, 0, 0).Let ( N , P , Ē) be a nontrivial critical point.In this case, f 1 ( Ē) = 0.The definition of f 1 implies that there is only one point Ē > 0. From the third equation, we obtain either Ē = 0 or f 3 ( P ) = f 1 ( Ē) = 0. Hence, f 3 ( P ) = 0.Because of the definition of f 3 , the point P > 0 is also unique.The second equation of system (2) implies N = δ P f 2 ( Ē) .The analysis above shows that for every parameter value there exist two critical points-the trivial one and the nontrivial Looking for other (semi-trivial) critical points, we obtain the following result.If N = 0, then P = 0 and the third equation implies f This means that such a critical point exists only when f 1 (0) < −a 3 (for f 1 (0) = −a 3 , the trivial and semi-trivial points coincide).Hence, for a 1 > a 3 , we have the third critical point, (0, 0, Ẽ) with f 1 ( Ẽ) = −a 3 .In this case Ẽ < Ē.
For a 1 ≤ a 3 , we have only two critical points (0, 0, 0), and ( N , P , Ē).Now we study the stability of the critical points.
Lemma 3.1.The nontrivial critical point ( N , P , Ē) is always unstable.Stability of the trivial critical point depends on the existence of a semi-trivial point, i.e., the magnitude of a 1 .If a 1 < a 3 , then the trivial point is stable; if a 1 > a 3 , then it is unstable and the semi-trivial point that appears is stable.
Proof.The Jacobi matrix for (2) has the following form: Hence, we have the following:   , and therefore, the trivial critical point is a saddle for a 1 > a 3 , or it is a stable node for a 1 < a 3 .There is a saddle-node bifurcation when (0, 0, Ẽ) appears.
is a stable node when it exists.
In this case, the characteristic polynomial is equal to W (λ) = −P (λ), where It is obvious that there exists at least one positive characteristic value λ 0 > 0. Hence, the nontrivial critical point is unstable independent of the parameters.Now, we study the behavior of solutions near ( N , P , Ē) in more details.Inequality P (λ) > 0 for λ > 0 guarantees the uniqueness of the positive zero of the polynomial P (λ).Hence, if other zeros of P (λ) are real, then they are negative.We show that if they are complex, then their real part is negative.Assume λ 1 = α + iβ and λ 2 = λ1 are complex conjugate characteristic values.Then Finally, the nontrivial critical point is either a saddle or unstable spiral.The direction of instability is described by the characteristic vector for λ 0 ; in the plane perpendicular to this characteristic vector, we have a stable focus.This characteristic vector is calculated to be equal to 3.2.The model in original variables.Coming back to the original variables N , P , V and the system described by ( 1) with E = V N , we see that (1) is not properly defined for N = 0. On the other hand, the whole right-hand side can be defined also for N = 0, taking Hence, we reformulate the model in the following way: where It is easy to see that both functions g i are continuous for nonnegative N , V. Therefore, the right-hand side of ( 4) is properly defined and continuous for every nonnegative value of N , P , V .Moreover, the right-hand side of ( 4) is locally a Lipschitz function.It is obvious for N = 0 because in this case the right-hand side is of class C 1 .Hence, it is enough to check this property at the neighborhood of points (0, V ), for arbitrary V ≥ 0. It is also enough to check this property for g i (N, V ), i = 1, 2, separately.
We have This implies that there exists a unique solution to (4) for every nonnegative initial condition.This also implies that (0, 0, 0) is the critical point for our system.Moreover, this is the only trivial critical point.In fact, ).If N = 0, then g 2 (N, V ) = 0, and this yields P = 0. From the third equation of (4) we obtain V = 0.
For N = 0, we have V = ĒN, and therefore, g 1 (N, ĒN ) = δP and f 3 (P ) ĒN = 0. Hence, f 3 (P ) = 0, and this implies P = P , N = N = δ P f 2 ( Ē) , and V = V = Ē N .We interpret the trivial critical points of (2) in the original variables.Both of them form the same trivial solution (0, 0, 0) to (4).The difference lies in the way the vessel volume depends on the tumor size when tumor size tends to 0. The trivial point (0, 0, 0) means that the ratio V N → 0 as V, N → 0, while the semitrivial one implies V N → Ẽ as V, N → 0. Hence, vessel volume depends either more than linearly (V ≈ N 1+γ , γ > 0, for the trivial case) or linearly (V ≈ ĒN , for the nontrivial one) on the tumor size.
We can conclude the analysis above with the following statement: Corollary 3.2.The trivial solution (0, 0, 0) to ( 4) is stable independent of the parameters.The nontrivial solution ( N , P , V ) is unstable independent of the parameters.
For the trivial solution (0, 0, 0), we obtain the following estimation.Assuming P → 0 and E → 0 from the first and third equations of (1), we get for every positive N , V .For P and E near 0, the following inequalities hold: of the negativity of V .Integrating the inequalities above with respect to t and changing variables in every integral, we obtain for every > 0 and some constants c 1 , c 2 > 0. This shows that E → 0 with the rate V 1− a 1 a 3 in the case a 3 > a 1 .

4.
Examples of the asymptotic behavior of (2).In this section, we present some examples of the behavior of solutions to (2) for some particular functions f i , i = 1, 2, 3, and discuss the dependence on these particular forms.For simulations, we use the Michaelis-Menten-type functions as well as the exponential type proposed in [3].In both cases, we obtained similar results.In Figures 1-5 the graphs correspond to the functions  ).The nontrivial one is (2, 1, 1).The initial-condition coordinates in this graph are (2, 1, 0.9), (1, 1, 1), (1, 1, 0.1), (1, 0.1, 0.1) and (1, 0.1, 1).All these points belong to the basin of attraction of the semi-trivial critical point.The asymptotic behavior of the system observed in simulations depends on the initial conditions and model parameters in the following way: if the initial condition belongs to the basin of attraction of the trivial or semi-trivial (if it exists) critical point, then the solutionobviously tends to this point (see Figures 1 and 2).Otherwise, N (t) → +∞.In such a case, we observe three types of the behavior: • N, P → ∞, E stabilizes at some level (Fig. 3), • N, E → ∞, P stabilizes at some level (Fig. 4), • all variables N, P, E → ∞ (Fig. 5).
The question arises about the connection between the asymptotic behavior and the coefficients of the model in general and the type of functions f i in the particular case used in our simulations.
If not, then there exists t such that E( t) = Ȇ and Ė( t) ≥ 0. On the other hand, the model implies Ė( t) < b 3 − f 1 ( Ȇ) Ȇ = 0, which contradicts the definition of t.Now, we show that if )) E 0 < 0, the function E(t) decreases at the beginning.But the same inequality is true for every point t such that E(t) > Ȇ.Hence, E(t) decreases until the inequality E(t) > Ȇ holds.Therefore, either E(t) > Ȇ for every t > 0 and then E(t) is decreasing, the inequality E(t) < E 0 holds, and E(t) has a limit, obviously; or there exists t such that E( t) = Ȇ, which means the previous estimation holds.
Proof.If P → ∞, then f 3 (P ) → b 3 .Hence, for every small > 0 there exists t such that for t > t we have Let Ȇ be such that f 1 ( Ȇ ) = b 3 − .Hence, if E < Ȇ , then E is increasing and for sufficiently large t we have E > Ȇ − .On the other hand, if E > Ȇ, then E is decreasing, so for sufficiently large t we have E < Ȇ + .Therefore, E → Ȇ.  3. The times for the corresponding graphs are t = 5 for the upper graph and t = 50 for the lower ones.As noted previously, we see three different orbits in the upper graph.In the lower-left graph, the second and third orbits join, while in the lower-right graph, the first and second orbits do not join.However, for all orbits, we observe an increase in all variables.
Remark 3. If E → Ȇ, then N, P → +∞ and P N → c, where c is some constant.
Proof.For E → Ȇ, the asymptotic behavior of N is described by Finally, we see that the behavior described in Lemma 4.1 and Remarks 2 and 3 does not depend on the particular choice of the functions f i .In two other cases, this does not hold true.
Hence, to obtain If the inverse inequality holds (i.e., 2b 1 < b 3 ,) then formula (5) implies P → 0, contrary to the assumption.Stabilization of P at some level P is possible when P is such that , and to obtain stabilization, one needs b 3 > b = 2b 1 .Another example.In the case analyzed above, if so this expression is a double exponential.Hence, stabilization is possible for every b 3 > b > b 1 .This shows that the results depend on the particular choice of f 2 .
Summarizing, in the case of f 2 (E) = a 2 E+1 , we observe the following three type of asymptotic behavior of large tumors:

5.
The model with time delays.In this section we focus on the model with delay (3).To solve equations (3), we define nonnegative continuous initial function (N (s), P (s), E(s)) on the interval [−τ, 0] where τ = max{τ i , i = 1, 2}, (see, e.g., [25]).The basic properties (such as existence, uniqueness, and nonnegativity) of (3) for such an initial function are the same as for (2) with nonnegative initial data.We are mainly interested in stability properties and their dependence on the delays.
Lemma 5.1.Stability and instability of the trivial and nontrivial critical points do not depend on the magnitude of both delays.Stability of the semi-trivial point depends on the magnitude of τ 1 .
Proof.The form of the right-hand side of (3) implies that stability properties of the trivial solution (0, 0, 0) do not change-expanding f 1 (E(t − τ 1 )) and f 3 (P (t − τ 2 )) at E = 0 and P = 0 (see e.g., [25]), we obtain the same free terms as for the case without delay; namely, f 1 (0) = −a 1 and f 3 (0) = −a 3 .These free terms form the linear part of the right-hand side for the trivial solution.Hence, the trivial solution is stable or unstable (depending on the parameters) independent of the delay.
To calculate the characteristic equation (see, e.g., [25]), we use three matrices: which is connected with the terms without delay; which is connected with the terms with τ 1 ; and which is connected with the term with τ 2 .Using M , M 1 , M 2 , we form the following characteristic equation: where I is the identity matrix.
For the semi-trivial solution, we obtain that the matrix M +M 1 e −λτ1 +M 2 e −λτ2 − λ I is a triangular one.Hence, its determinant is equal to the product of the terms on the main diagonal; that is, −(a 3 + λ)(δ + λ)(f 1 ( Ẽ) Ẽe −λτ 1 + λ) = 0.This yields that the semi-trivial solution is stable for f 1 ( Ẽ) Ẽτ 1 < π 2 and unstable for greater τ 1 .Stability does not depend on the magnitude of τ 2 .For a detailed derivation, see [17].
The nontrivial solution ( N , P , Ē) is unstable independent of the magnitude of both delays-the characteristic equation always has a real positive zero.This equation reads as Looking for the real solutions to the equation above, we see that the right-hand side is positive a decreasing function of λ ∈ R + , while the left-hand side is equal to 0 for λ = 0 and tends to +∞ as λ → +∞.This implies the existence of a real positive solution.
It should be noted that if the critical point is unstable in the case without delay, then it is also unstable for every discrete non-zero delays (it can be shown using, for example, Mikhailov criterion; for details, see [17,23,28]).
In [3], it is shown that the Hopf bifurcation in the nontrivial critical point appears in the model with nonzero delays (3), while for the model without delay (2), it is not possible.6. Discussion.We analyzed the stability properties of the simplest models of vascular tumor dynamics introduced in [3]-the difference between them being the assumption of delay in the changes in cells growth rate.We have shown the existence and uniqueness of a positive solution for any positive initial conditions in both models.In both models, there exist two or three critical points, depending on the parameters.There is always the trivial point (0,0,0) and the nontrivial point; the semi-trivial point can appear, depending on parameter values.The nontrivial critical point is always unstable.For the model without delay (2), if the semi-trivial critical point exists then, it is always stable; otherwise, the trivial critical point is stable.
We have also shown that the existence of the semi-trivial critical point and its bifurcation from the trivial point are the result of substitution E = V N .System (2) can be reformulated in the original variables (4), having only one trivial critical point.This point is always stable.The trivial and semi-trivial critical points for (2) both correspond to this trivial critical point for (4)-the difference being in the rate of convergence to zero for V and N .
For the system without delay (2), the nontrivial point is either a saddle or unstable spiral, so the solution starting around it will grow to infinity, possibly in the asymptotic direction; see section 3.
For the system with delay (3), the trivial and the nontrivial critical points have the same stability properties as in the model without delay.The semi-trivial point (if it exists) will be stable for τ 1 small enough.In this system, the Hopf bifurcation is possible for the nontrivial critical point.However, this point will always be unstable.Typically, the stable Hopf bifurcation appears when the periodic orbit bifurcates from the stable critical point.This is not the case, but stability of the Hopf bifurcation is not excluded.We suspect that the solution around the nontrivial critical point will grow to infinity, possibly oscillating in the plane orthogonal to the direction of growth.
Both systems predict stability of the trivial or the semi-trivial critical point.This means that starting close enough to this point, the solution will converge to it.For solutions starting close enough to the nontrivial critical point, the models suggest unlimited growth (for both cell number and vessel volume nearly proportionally), possibly characterized by oscillations in the delayed model.This corresponds to the behavior of experimentally observed growing vascular tumors that show oscillations both in tumor and vascular mass during their growth [22].
The characteristics of the critical points in the discussed models rule out the possibility of tumor-load stabilization, which is sometimes observed experimentally [14,32,33].Admittedly, cancer growth law is still a controversial phenomenon [6,21,26].Nevertheless, it can be suggested by this work that angiogenesis per se is not sufficient for describing asymptotic tumor growth.Some additional processes (such as vessel maturation and regression) should be introduced to improve this description.Another limitation of these models is that they do not describe the initiation of tumor growth and the onset on angiogenesis, since for small tumors with small vessel volume, the solution will tend to zero.This results from our model assumptions that neglect both the ability of small tumors to survive and grow in the absence of blood vessels; and the migration of endothelial cells, which is important at the onset of the angiogenic process.Such a neglection is justified, since the models aim to describe the angiogenic phase, in which both tumor size and vascular mass are relatively large.
We believe that a more comprehensive description can be obtained by including models that represent additional regulation pathways for angiogenesis.At present, we are extending our model to include more biological processes involved in angiogenesis as well as tumor growth at the early stages.

3 .
These are the simplest functions of the Michaelis-Menten type that fulfill the assumptions.For all figures, we have δ = a 2 = 1 and different values of d i and a i , i = 1, 3.

Figure 2 .
Figure 2.For this graph, we used d 1 = 1, a 1 = 0.5, d 3 = 2, and a 3 = 1.The semi-trivial critical point does not exist, because a 3 > a 1 .The nontrivial one has the same coordinates as in Fig. 1.The initial-condition coordinates are also the same.All these points belong to the basin of attraction of the trivial critical point.

Figure 3 .
Figure 3.The parameters of the model are the same as in Figure 1.We have b 1 > b 3 .The initial-condition coordinates are (20, 20, 11), (10, 5, 30), (0.1, 10, 20).The upper-left graph corresponds to t = 5-we see three different orbits.The upper-right graph corresponds to t = 15-now the orbits join and appears as one orbit.The lower-left graph shows the stabilization of E, and the lower-right graph describes the linear asymptotic dependence between N and P .

Figure 4 .
Figure 4.The parameters of the model are d 1 = 1, a 1 = 0.5, d 3 = 4, and a 3 = 1.We have b 3 > 2b 1 .The initial-condition coordinate are the same as in Figure 3.The times for the corresponding graphs are also the same.As noted previously, we see three different orbits in the upper-left graph, while in the upper-right graph the orbits join and appeares as one orbit.The lower-left graph shows the stabilization of P and the lower-right graph describes the linear asymptotic dependence between E and N .

Figure 5 .
Figure 5.The parameters of the model are d 1 = 1, a 1 = 0.5, d 3 = 1.6, and a 3 = 1.We have b 1 < b 3 < 2b 1 .The initialcondition coordinates are the same as in Figure3.The times for the corresponding graphs are t = 5 for the upper graph and t = 50 for the lower ones.As noted previously, we see three different orbits in the upper graph.In the lower-left graph, the second and third orbits join, while in the lower-right graph, the first and second orbits do not join.However, for all orbits, we observe an increase in all variables.

Remark 4 .
If b 3 > b 1 and E, P → ∞, then E → ∞ exponentially with the exponent equal to b 3 − b 1 .Proof.It is obvious because of the asymptotic growth of E described by Ė ≈ (b 3 − b 1 )E.Consider the function f 2 then we expect the unlimited growth of N , P and the stabilization of E. • If b 1 < b 3 < 2b 1 , then we expect the unlimited growth of every variable N , P , E. • If b 3 > 2b 1 , then we expect the unlimited growth of N , E and the stabilization of P .