Angiogenesis model with Erlang distributed delays

MATHEMATICAL BIOSCIENCES AND ENGINEERING Volume 14, Number 1, February 2017 doi:10.3934/mbe.2017001 pp. 1–15 ANGIOGENESIS MODEL WITH ERLANG DISTRIBUTED DELAYS Emad Attia Department of Mathematics, Faculty of Science Damietta University New Damietta 34517, Egypt Marek Bodnar and Urszula Foryś Institute of Applied Mathematics and Mechanics Faculty of Mathematics, Informatics and Mechanics University of Warsaw Banacha 2, 02-097 Warsaw, Poland Abstract. We consider the model of angiogenesis process proposed by Bodnar and Foryś (2009) with time delays included into the vessels formation and tumour growth processes. Originally, discrete delays were considered, while in the present paper we focus on distributed delays and discuss specific results for the Erlang distributions. Analytical results concerning stability of positive steady states are illustrated by numerical results in which we also compare these results with those for discrete delays. 1. Introduction. In this paper we focus on the process of angiogenesis, which means formation of new vessels from pre-existing ones. It is a normal and vital process in growth and development of animal organisms, but it is also essential in the transition of avascular forms of solid tumours into metastatic ones. Small tumours (less than 1–2 mm3 of volume) receive nutrients from the outside by simple diffusion. Then necrotic core starts to grow, which is associated with the secretion of angiogenic factors by cancer cells. The supply of nutrients allows cancer to grow further, and moreover vessels form the path for spreading cancer cells in the organism. Mathematical modelling of this process is inextricably linked with the idea of anti-angiogenic treatment first considered by Folkman in 1971 (c.f. [15]). However, as anti-angiogenic agents were not known that time, it took more than 20 years for the specific models of angiogenesis and anti-angiogenic treatment appear (c.f. [17] were one of the best known models of it has been proposed). Various approaches were used to describe the angiogenesis process in mathematical language. The simplest approach is based on ordinary differential equations describing the dynamics of tumour and vasculature, like in the approach of Hahnfeldt et al. [17]. This idea was extended by many authors, including d’Onofrio and Gandolfi [14], Agur et al. [2], Bodnar and Foryś [8], Poleszczuk et al. [21], Piotrowska and Foryś [20, 16]. However, in many cases spatio-temporal dynamics of 2010 Mathematics Subject Classification. Primary: 34K11, 34K13, 34K18, 37N25; Secondary: 92B05. Key words and phrases. Delay differential equations, stability analysis, Hopf bifurcation, angiogenesis, tumour growth. 1 2 EMAD ATTIA, MAREK BODNAR AND URSZULA FORYŚ vessels and tumour structure is important. To reflect such a structure one can use the approach of partial differential equations. First models of that type was based on reaction-diffusion equations with the process of chemotaxis taken into account; cf. [12]. Many papers focusing on chemotaxis modelling was published by M. Chaplain and coauthors (cf. [12], [19] where also haptotaxis process was considered, [3] with the influence of external forces). Yet another approach is based on cellular automata. Models of that type was proposed by Rieger and coauthors (cf. [6, 23]), while Alarcón et al. [1] considered hybrid cellular automaton in the context of multiscale modelling. Very nice review of various types of spatio-temporal models of vasculogenesis and angiogenesis processes could be find in [22], where a list of many other references on that topic is available. In this paper we consider the model of tumour angiogenesis proposed in [7] and studied in [9] in the context of discrete delays. However, discrete delays could be only an approximation of delays present in real life, and therefore, following [10] we decided to incorporate distributed delays and compare the results with those for the discrete case. We mainly focus on the Erlang distributions, however some results for general distributions are also obtained. Partial results for the distributed delay in the vessel formation process was presented in [5]. Here, we incorporate distributed delays both in the vessel formation and tumour growth processes. Thus we consider the following system of the first order differential equations with distributed delays   N (t) Ṅ (t) = αN (t) 1 − , 1 + f1 (h1 (Et )) Ṗ (t) = f2 (E(t))N (t) − δP (t),  Ė(t) = f3 (h2 (Pt ))ds − α 1 − (1.1) N (t) 1 + f1 (h1 (Et )) ! E(t), where N (t), P (t), E(t), α and δ describe the tumour size, the amount of regulating proteins, the effective vessel density, the tumour proliferation rate and the degradation of proteins, respectively. Both parameters could also include the influence of treatment. The functions hi : C((−∞, 0], R), i = 1, 2, are defined as follows Z ∞ hi (φ) = ki (s)φ(−s) ds, 0 and we assume ki (s) : [0, ∞) → [0, ∞) are probability densities with finite expectations, that is Z ∞ Z ∞ ki (s) ds = 1 and 0 < ski (s) ds < ∞, i = 1, 2. 0 0 For any φ ∈ C(R, R), the function φt is defined as φt (s) = φ(t+s) for s ∈ (−∞, 0]. Moreover, we assume that the functions fi are locally Lipschitz and there exist positive constants A2 , A3 , B1 , B3 and m3 such that (A1) f1 is increasing, f1 (0) = 0 and limE→+∞ f1 (E) = B1 > 0, (A2) f2 is decreasing and convex, f2 (0) = A2 > 0 and limE→+∞ f2 (E) = 0, (A3) f3 is increasing, f3 (0) = −A3 < 0, f3 (m3 ) = 0 and limP →+∞ f3 (P ) = B3 . For detailed derivation of the model described by Eqs. (1.1) we refer to [7, 9]. To close the problem we need to define initial data. Let C = C((−∞, 0], R3 ) be the space of continuous functions defined on the interval (−∞, 0] with values in R3 , and as C+ we define the subspace of C that consist of the functions with ANGIOGENESIS MODEL WITH ERLANG DISTRIBUTED DELAYS 3 non-negative values. Because the distributions of delays have infinite supports, we need to control the behaviour of initial functions at infinity (see [18]). To this end, let define a Banach space Φ in the following way ( ) Φ= kφkΦ = φ ∈ C : lim φ(s)es = 0 and s→−∞ sup s∈(−∞,0] sup s∈(−∞,0] |φ(s)es | < ∞ , s |φ(s)e |, and we consider initial functions from the set Φ+ = Φ ∩ C+ . 2. Analysis of the model. In this section we consider basic properties of system (1.1) for general distributions ki (s) and study stability in the case of Erlang distributions. For all results presented here we assume the functions fi , i = 1, 2, 3, fulfil conditions (A1)–(A3). Theorem 2.1. For any arbitrary initial function φ = (φN , φP , φE ) ∈ Φ+ there exists a unique solution of system (1.1) in Φ+ defined on t ∈ [0, +∞). Moreover, for all t ≥ 0 any solution satisfies the following inequalities Nmin ≤ N (t) ≤ Nmax ,   A2 0 ≤ P (t) ≤ max Nmax , φ2 (0) , δ (2.1)
0 ≤ E(t) ≤ φ3 (0) exp (B3 + α(Nmax − 1))t , where Nmin = min {1, φN (0)} , Nmax = max {φN (0), 1 + B1 } . Proof. It is easy to see that the right-hand side of system (1.1) is locally Lipschitz, which yields local existence of the solution of (1.1). Non-negativity follows from the form of this right-hand side. Inequalities (2.1) could be obtained in the same way as in [9]. Then the global existence of the solutions can be proved by the use of Theorem 2.7 from [18, Chapter 2]. Now, we turn to steady states. It is obvious that there are at least two steady states  A  2 A = (0, 0, 0) and B = 1, ,0 , δ compare [9]. Moreover, there can exist positive steady states Di = (N̄i , m3 , Ēi ), with N̄i = 1 + f1 (Ēi ), and Ēi are zeros of the function g(x) = f2 (x)(1 + f1 (x)) − δm3 . (2.2) Stability of the steady states A and B does not depend on time delays, as in the discrete case; c.f. [9]. We recall the results without the proof. Proposition 2.2. The trivial steady state A of system (1.1) exists and is unstable independently of the model parameters. The semi-trivial steady state B of system (1.1) is locally asymptotically stable for A2 < δm3 and unstable for A2 > δm3 . 4 EMAD ATTIA, MAREK BODNAR AND URSZULA FORYŚ 2.1. Stability of positive steady states. In this section, we focus on examining the stability of the positive steady states Di in the case of distributed delays. First, we give some general results regarding stability and later we consider Erlang probability distribution, which is a special type of Gamma distribution with the shape parameter being a natural number. In the general case let us define Ki (λ) = Z ∞ ki (s)e−λs ds. (2.3) 0 Then, the stability matrix of system (1.1) for the steady state Di reads   −α − λ 0 αd1 K1 (λ) , −δ − λ −N̄i d2 M (N̄ , P̄ , Ē) =  f2 (Ēi ) αbĒi d3 Ēi K2 (λ) −αbd1 Ēi K1 (λ) − λ where 1 , d1 = f1′ (Ēi ) > 0, d2 = −f2′ (Ēi ) > 0, 1 + f1 (Ēi ) Hence, the characteristic quasi-polynomial has the form b= d3 = f3′ (m3 ) > 0. W (λ) =λ3 + C1 λ2 + C2 λ + (λ2 + δλ)C3 K1 (λ)+ + (λ + α)C4 K2 (λ) − C3 C5 K1 (λ)K2 (λ), with (2.4) βd2 d3 , C5 = δd3 m3 , β = bĒi . b2 Conditions (A1)–(A3) guarantee positivity of b, β and di , i = 1, 2, 3. Consequently, Ci > 0 for i = 1, . . . , 5. C1 = α + δ, C2 = αδ, C3 = αβd1 , C4 = Theorem 2.3. If g ′ (Ēi ) > 0, where g is given by (2.2), then the positive steady state Di = (N̄i , m3 , Ēi ) is unstable. Proof. We show that the characteristic function (2.4) has at least one positive real root. In the proof of Theorem 3.4 in [9] it is shown that the sign of W (0) = αC4 − C3 C5 is reverse to the sign of g ′ (Ēi ). Therefore, the assumption g ′ (Ēi ) > 0 implies that W (0) < 0. Further, it is easy to see that W (λ) → +∞ as λ → +∞. Then there exists at least one real positive eigenvalue, and this implies that Di is unstable. Now, we focus on the cases when only one of the considered processes is delayed and the other is instantaneous. First, we consider k2 (s) = δD (s), where δD means Dirac-delta distribution, that is only the first delay is present in the system. In the theorem presented below we shall use the following auxiliary polynomials and positive zeros of these polynomials. Let us define w1 (ω) = −ω 3 − C3 ω 2 + (C2 + C4 − δC3 )ω − C3 C5 , 2 and w3 (ω) = −(C1 + C3 )ω − δC3 ω + αC4 − C3 C5 , (2.5) (2.6) w4 (ω) = −(C1 − C3 )ω 2 + δC3 ω + αC4 + C3 C5 . (2.7) To obtain positive zeros of w1 one needs to assume C2 + C4√− δC3 > 0. More−2C + 4C 2 +12(C +C −δC ) 3 2 4 3 3 over, the value at positive extremum, that is at Pmax = , 6 must be positive. Hence, it is necessary that w1 (Pmax ) > 0. In such a case ANGIOGENESIS MODEL WITH ERLANG DISTRIBUTED DELAYS 5 there exist two positive zeros P1 < P2 of w1 . Next, we easily see that if g ′ (Ēi ) < 0, √ that is αC4 − C3 C5 > 0, then w3 has exactly one positive zero P3 = −δC3 + δ 2 C32 −4(C1 +C3 )(αC4 −C3 C5 ) . Similarly, if C1 > C3 , then w4 has exactly one 2(C1 +C3 ) positive zero p δC3 + δ 2 C32 + 4(C1 − C3 )(αC4 + C3 C2 ) P4 = . 2(C1 − C3 ) In the following, we require P1 < P3 and P4 < P2 . Theorem 2.4. Assume that g ′ (Ēi ) < 0 and Di = (N̄i , m3 , Ēi ) is a positive steady state of system (1.1). Assume also that the delay distribution k2 is concentrated at zero, that is k2 (s) = δD (s). If (i) βd1 < 1 and 2α2 b < d3 < min d2 Ēi  


δ b 2 2 , δ 2 − α 2 β 2 d1 2 − 1 , 1 − β d 1 2d2 Ēi 2β 2 d1 2 m3 or (ii) C1 > C3 , w1 (Pmax ) > 0 and w1 (Pi ) > 0 for i = 3, 4, where w1 is defined by (2.5) and P3 , P4 are positive zeros of (2.6) and (2.7), respectively, then Di is locally asymptotically stable for any distribution k1 . Proof. The steady state Di is stable in the case without delay for g ′ (Ēi ) < 0. Hence, we need to show that the characteristic function has no purely imaginary zeros. Therefore, we investigate the behaviour of W (iω). Let K1 (iω) = η1 − iζ1 , η1 = Z∞ 0 k1 (s) cos(ωs) ds, ζ1 = Z∞ k1 (s) sin(ωs) ds. 0 Thus W (iω) = −iω 3 − C1 ω 2 + i(C2 + C4 )ω + ((−ω 2 + iδω)C3 − C3 C5 )(η1 − iζ1 ) + αC4 , and we have Re(W (iω)) = −C1 ω 2 + αC4 − C3 (ω 2 + C5 )η1 + δC3 ωζ1 , (2.8) Im(W (iω)) = −ω 3 + (C2 + C4 )ω + δC3 ωη1 + C3 (ω 2 + C5 )ζ1 . Assume that there exists ω > 0 such that W (iω) = 0, then
2 + −ω 3 + (C2 + C4 )ω =
2
2 −C3 (ω 2 + C5 )η1 + δC3 ωζ1 + δC3 ωη1 + C3 (ω 2 + C5 )ζ1 . −C1 ω 2 + αC4
2 (2.9) For Eq. (2.9) we have

L.H.S = ω 6 + C12 − 2(C2 + C4 ) ω 4 + (C2 + C4 )2 − 2αC1 C4 ω 2 + α2 C42 ,


R.H.S = C32 ω 4 + δ 2 + 2C5 ω 2 + C52 ζ12 + η12 . 6 EMAD ATTIA, MAREK BODNAR AND URSZULA FORYŚ Schwarz inequality yields 2  ∞ ∞ Z Z Z  cos(ωs)k1 (s) ds =  cos(ωs) d 0 0 ≤ = Z∞ 0 Z∞ 2 cos (ωs) d Z 0 s 0 s  2 k1 (u) du   Z∞ Z d k1 (u) du 0 s k1 (u) du 0  cos2 (ωs)k1 (s) ds, 0 2 ∞ Z∞ Z  sin(ωs)k1 (s) ds ≤ sin2 (ωs)k1 (s) ds. 0 0 Consequently, ζ12 + η12 ≤ 1. Then, ω satisfying W (iω) = 0 should fulfil 0 = L.R.S. − R.H.S.    2 6 2 4 ω + C1 − 2C̃2,4 ω + C̃2,4 − 2αC1 C4 ω 2 +

+ α2 C42 − C32 ω 4 + δ 2 + 2C5 ω 2 + C52 ,  where C̃2,4 = C2 + C4 . Let us denote y = ω 2 and define the auxiliary function     2  2 2 2 3 2 2 C̃2,4 − 2αC1 C4 − C3 δ − 2C5 F (y) = y + C1 − 2C̃2,4 − C3 y + y+ + α2 C42 − C32 C52 , Existence of purely imaginary eigenvalue requires F (y) ≤ 0 for some positive y. However, all coefficients of this polynomial are positive under the assumptions, and √ then purely imaginary eigenvalue i y does not exist. Clearly, the free term is positive due to g ′ (Ēi ) < 0. Moreover, inequalities (i) are equivalent to β d 2 d3 βd1 < 1, α2 < and 2b2 ) (r 2β d2 d3 2β 2 d1 2 d3 m3 2 2 2 , α (β d1 − 1) + . δ > max b2 1 − β 2 d1 2 Therefore,
2β d2 d3 C12 − 2(C2 + C4 ) − C32 = δ 2 − α2 β 2 d1 2 − 1 − > 0, b2 (C2 + C4 )2 − 2αC1 C4 − δ 2 C32 − 2C32 C5 = 2  2α (δ + α) β d2 d3 β d 2 d3 − αδ + − δ 2 α2 β 2 d1 2 − 2α2 β 2 d1 2 δ d3 m3 = b2 b2  

β d 2 d3 β d 2 d3 2 1 − β 2 d1 2 α2 δ − 2α2 β 2 d1 2 d3 m3 δ + > 0. − 2α b2 b2 Due to the continuous dependance of eigenvalues on the model parameters this completes the proof of the first part. ANGIOGENESIS MODEL WITH ERLANG DISTRIBUTED DELAYS 7 For the proof of the second part, notice that (2.8) implies Re(W (iω)) ≥ − (C1 + C3 )ω 2 − δC3 ω + αC4 − C3 C5 , Re(W (iω)) ≤ − (C1 − C3 )ω 2 + δC3 ω + αC4 + C3 C5 , Im(W (iω)) ≥ − ω 3 − C3 ω 2 + (C2 + C4 − δC3 )ω − C3 C5 . It is clear that, Re(W (iω)) > Re(W (iω)) < 0 for ω ∈ [0, P3 ), 0 for ω ∈ (P4 , ∞). and Im(W (iω)) > 0 for ω ∈ (P1 , P2 ), Inequalities w1 (Pi ) > 0 for i = 3, 4, imply P1 < P3 and P2 > P4 . Moreover, w3 (0) < w4 (0) and w3′ (ω) < w4′ (ω) for ω > 0. Therefore P3 < P4 and there is no ω ≥ 0 such that W (iω) = 0. This completes the proof of the theorem. Remark 1. If the second condition of assumption (i) of Theorem 2.4 is satisfied, then δ 2 > 3α2 . If this inequality holds, we can choose sufficiently small d1 and later we can take a proper value of d3 such that the condition (i) of Theorem 2.4 holds. In the following, we shall consider Erlang distributed delays. The density of Erlang distribution is given by  m m−1  a (s − σ) e−a(s−σ) , s ≥ σ, (m − 1)! (2.10) k(s) =  0, otherwise, where a, σ > 0, and a is a scaling parameter. To the case σ = 0 we refer as to the non-shifted Erlang distribution, while to the case σ > 0 we refer as to the shifted one. The mean value of the non-shifted Erlang distribution is given by m a and the variance is equal to am2 . The average delay is equal to this mean, obviously, and the deviation measures the concentration of the delay about the average value. For the shifted Erlang distribution the mean value is σ + m a and the variance is the same as for the non-shifted one. On the other hand, taking the limit m → +∞ and with discrete delay τ . By direct keeping m/a = τ constant we R ∞obtain system (1.1) am −λσ . calculation, it is found that 0 k(s)e−λs ds = (a+λ) me Now, let us consider the first distribution to be Erlang, k1 (s) = k(s), while the other is still k2 (s) = δD (s). In this case, if C5 6= a(a − δ), then the characteristic function (2.4) could be replaced by  
WI (λ) = (a + λ)m λ3 + C1 λ2 + C2 + C4 λ + αC4 +   (2.11) + am C3 λ2 + δC3 λ − C3 C5 e−λσ . As the case C5 = a(a − δ) is non-generic, in the following we assume C5 6= a(a − δ), and consider the simplest possible Erlang distribution with m = 1 and σ = 0. Theorem 2.5. If m = 1, σ = 0 and g ′ (Ēi ) < 0, then the positive steady state Di is locally asymptotically stable. Proof. For this case the characteristic function WI simplifies to    
WI (λ) = (a + λ) λ3 + C1 λ2 + C2 + C4 λ + αC4 + a C3 λ2 + δC3 λ − C3 C5 . 8 EMAD ATTIA, MAREK BODNAR AND URSZULA FORYŚ The Routh-Hurwitz Criterion for WI reads q1 q2 q3 > q32 + q12 q4 , (2.12) where q1 = a + C 1 , q2 = a(C1 + C3 ) + η2 , q3 = a(η2 + δC3 ) + αC4 , q4 = aη4 , and η2 = C2 + C4 , η4 = αC4 − C3 C5 , Notice that inequality (2.12) is equivalent to P (a) = a3 a3 + a2 a2 + a1 a + a0 > 0, where a3 = (C1 + C3 )(η2 + δC3 ) − η4 , a2 = αC4 (C1 + C3 ) + C2 (η2 + δC3 ) + C1 (C1 + C3 )(η2 + δC3 ) + C4 (η2 + δC3 ) − (η2 + δC3 )2 − 2C1 η4 , a1 = C1 (η22 + αC3 C4 + δC3 η2 ) − αC4 (η2 + 2δC3 ) + C12 C3 C5 , a0 = αC4 (C1 η2 − αC4 ). Due to the definitions of C1 and ηi we have:    


a0 = αC4 δ + α C2 + C4 − αC4 = αC4 δ + α C2 + δC4 > 0,   a1 = (α + δ) (C2 + C4 )2 + αC3 C4 + δC3 (C2 + C4 ) − αC4 (C2 + C4 )+ − 2αδC3 C4 + C12 C3 C5  
= α η2 C2 + αC3 C4 + δC2 C3 + δη2 η2 + δC3 > 0, a2 = αC4 C3 + C2 (η2 + δC3 ) + (C1 (α + δ) + C1 C3 )(η2 + δC3 ) + C4 (η2 + δC3 )+ − (η2 + δC3 )2 − αC1 C4 + 2C3 C5 , = αC4 C3 + (δC1 + (α + δ)C3 )(η2 + δC3 ) + αC1 (C2 + δC3 )+ − δC2 C3 − δC4 C3 − δ 2 C32 + 2C3 C5 , = αC4 C3 + (δC1 + α)C3 (η2 + δC3 ) + αC1 (C2 + δC3 ) + 2C3 C5 > 0, a3 = (δ + α + C3 )(C2 + C4 + δC3 ) − αC4 + C3 C5 = (δ + C3 )(C2 + C4 + δC3 ) + α(C2 + δC3 ) + C3 C5 > 0. Hence, P (a) > 0 for every positive a. This means that the condition of stability is always satisfied. Remark 2. Although Theorem 2.4 gives condition guaranteeing stability of the positive steady state for any delay distribution, Theorem 2.5 says that the positive steady state, if it is stable for the case without delay, cannot lose stability when the tumour growth process is delayed according to the Erlang distribution. Now, we switch to the case when k1 (s) = δD (s), while the second distribution is Erlang, that is k2 (s) = k(s) described by Eq. (2.10). Note that if αC4 6= C4 a+C3 C5 , then λ is the eigenvalue defined by (2.4) if and only if λ is zero of 

 WII (λ) = (a + λ)m λ3 + C1 + C3 λ2 + C2 + δC3 λ +   (2.13) + am C4 λ + αC4 − C3 C5 e−λσ . ANGIOGENESIS MODEL WITH ERLANG DISTRIBUTED DELAYS 9 Again, because the case αC4 = C4 a + C3 C5 is non-generic, we do not consider it here, and in the following we assume αC4 6= C4 a + C3 C5 . Thus, studying the stability of the positive steady states Di of system (1.1) is equivalent to study zeros of the polynomial WII defined by (2.13). Below, we again study the simplest case of Erlang distribution with m = 1 and σ = 0. Proposition 2.6. If m = 1, σ = 0, g ′ (Ēi ) < 0 and Q1 Q2 Q3 > Q23 + Q21 Q4 , where Q1 = C1 + C3 + a, Q2 = C2 + δC3 + a(C1 + C3 ), Q3 = a(C2 + C4 + δC3 ), Q4 = a(αC4 − C3 C5 ), then the positive steady state Di is locally asymptotically stable. Proof. For m = 1 and σ = 0, the polynomial WII defined by (2.13) reads WII (λ) = λ4 + (C1 + C3 + a)λ3 + (C2 + δC3 + a(C1 + C3 ))λ2 + + a(C2 + C4 + δC3 )λ + a(αC4 − C3 C5 ), and the assertion of the proposition comes directly from the Routh-Hurwitz Criterion. Now, we try to answer the question when the assumptions of Proposition 2.6 are satisfied. To simplify calculations and shorten notation, let us denote η1 = C1 + C 3 , η2 = C2 + δC3 , With this notation we have Q1 = η1 + a, Q2 = η2 + aη1 , η4 = αC4 − C3 C5 . Q3 = a(η2 + C4 ), Q4 = aη4 , Qi > 0 for i = 1, . . . , 4. Now, the Routh-Hurwitz condition is equivalent to    


a2 η1 η2 + C4 − η4 + a η2 + C4 η12 − C4 − 2η1 η4 +  
+ η1 η2 η2 + C4 − η1 η4 > 0. (2.14) Notice that the coefficient of a2 is positive. Clearly, using the definitions of η1 , η4 and C1 we have
η1 η2 +C4 −η4 = η1 η2 +(α+δ+C3 )C4 −αC4 +C3 C5 = η1 η2 +(δ+C3 )C4 +C3 C5 > 0. Now, we have only three possibilities: 1. η2 (η2 + C4 ) < η1 η4 – there exists exactly one critical value of a, below which the steady state is unstable and above which it is stable (average delay is 1/a in this case);

2. η2 + C4 η12 − C4 /2 < η1 η4 < η2 (η2 + C4 ) and the discriminant of the quadratic polynomial is positive – there exist
two critical
values of a; 3. η1 η4 < η2 (η2 + 1) and either η1 η4 < η2 + C4 η12 − C4 /2 or the discriminant of the quadratic polynomial is negative – the steady state is stable for all a. To obtain two changes of stability, we need to have  2   



η2 +C4 η12 −C4 −2η1 η4 > 4η1 η2 η2 +C4 −η1 η4 η1 η2 +C4 −η4 , (2.15) together with (η12 − C4 )(η2 + C4 ) < η1 η4 < η2 (η2 + C4 ). 2 10 EMAD ATTIA, MAREK BODNAR AND URSZULA FORYŚ Inequality (2.15) is equivalent to

2

η2 + C4 η12 − C4 − 4η1 η4 η12 − C4 − 4η12 η2 η2 + C4 + 4η1 η2 η4 + 4η13 η4 > 0, and collecting terms with η12 − C4 we obtain 

2


 η2 +C4 η12 −C4 −4η1 η4 η12 −C4 +4η1 η4 η2 +η12 −η1 η2 η2 +C4 > 0. (2.16) Notice that the free and linear terms of (2.16) are positive under the assumption
η1 η2 η2 + C 4 η4 > and η12 < C4 . η2 + η12 We have η12 − C4 = (α + δ + αδ)2 − βd2 d3 /b2 , and it is negative for sufficiently large d2 d3 or sufficiently small b. Eventually, two stability switches are possible under the assumptions
η 1 η2 η2 + C 4 η2 (η2 + C4 ) < η4 < and η12 < C4 . 2 η2 + η1 η1 Proposition 2.7. If g ′ (Ēi ) 6= 0 and (i) Di is unstable for σ = 0, then it is unstable for all σ > 0; (ii) Di is stable for σ = 0, then there exists σ0 > 0, such that Di is stable for σ < σ0 , and Di is unstable for σ > σ0 . Furthermore, if fi ∈ C2 , i = 1, 2, 3, then Hopf bifurcation occurs at σ0 . Proof. For the characteristic function WII given by (2.13), we define the auxiliary function
m 2
F (y) = y a2 + y y + ((C1 + C3 )2 − 2(C3 δ + C2 ))y + (C3 δ + C2 )2 − − a2m C42 y − a2m (C4 α − C3 C5 )2 , where y = ω 2 . Notice, that (C1 + C3 )2 − 2(C3 δ + C2 ) = α2 (1 + βd1 )2 + δ 2 > 0. Clearly, F has at least one positive zero, since the coefficient of y with the highest power is positive, while the free term is negative. We show that this is the unique simple positive root. Note, that all coefficients of F are positive with the exception of the free term which is negative and the coefficient of y, which sign is not determined. However, in both cases, there exists exactly one change of sign in the coefficients of the polynomial F , and the Descartes’ Rule of Signs implies that F has a unique and simple positive root. This, together with Theorem 1 from [13], completes the proof. Eventually, we consider both distributions to be Erlang with the same parameters. 2 2 Proposition 2.8. If m = 1, σ = 0, g ′ (Ēi ) < 0, q11 q22 q33 > q33 + q11 q44 and where 2 2 2 (q11 q44 − q55 )(q11 q22 q33 − q33 − q11 q44 ) > q55 (q11 q22 − q33 )2 + q11 q55 , q11 = C1 + 2a , q22 = C2 + 2a C1 + a2 + aC3 , q33 = 2a C2 + a2 C1 + a2 C3 + aC3 δ + aC4 , q44 = a2 C2 + a2 C3 δ + a2 C4 + aC4 α, then the positive steady state Di is locally stable. q55 = a2 C4 α − a2 C3 C5 , ANGIOGENESIS MODEL WITH ERLANG DISTRIBUTED DELAYS 11 Proof. For m = 1 and σ = 0, the characteristic function (2.4) reads 
1 W (λ) = λ3 + C1 λ2 + C2 λ (a + λ)2 + a(λ2 + δλ)C3 (a + λ) + (a + λ)2  + a(λ + α)C4 (a + λ) − a2 C3 C5 . Hence, we need to study roots of the polynomial λ5 +(C1 +2a )λ4 +(C2 +2a C1 +a2 +aC3 )λ3 +(2a C2 +a2 C1 +a2 C3 +aC3 δ+aC4 )λ2 + + (a2 C2 + a2 C3 δ + a2 C4 + aC4 α)λ + a2 C4 α − a2 C3 C5 = 0. A direct application of the Routh-Hurwitz Criterion completes the proof. 3. Numerical simulations. For the numerical simulations we choose functions fi and parameters values proposed in [9], that is f1 (E) = b1 E n , c1 + E n f2 (E) = a2 c2 , c2 + E f3 (P ) = b3 (P 2 − m23 ) m23 b3 a3 + P2 , and a2 = 0.4, a3 = 1, b1 = 2.3, b3 = 1, c1 = 1.5, c2 = 1, α = 1, δ = 0.34. (3.1) For these values of parameters there exist three positive steady states: D1 ≈ (1.04, 1.05, 0.17), D2 ≈ (1.37, 1.05, 0.54), D3 ≈ (2.67, 1.05, 1.99). Now, we can influence the model dynamics changing the value of δ. This parameter could reflect an application of some treatment that blocks VEGF activation, which can be modelled by the increase of clearance rate. Then the steady state D1 exists for 0.331 < δ < 0.381, D3 for δ < 0.368 and D2 for 0.331 < δ < 0.368. The steady state D2 is unstable, while D1 and D3 are stable for the case without time delay. In Table 1 we presented critical values of delay at which the steady states D1 Table 1. Critical values of τ at which the positive steady state loses stability. δ discrete m=1 m=2 m=5 steady state D1 0.332 0.346 0.36 0.368 0.378 66.7 33.4 29.3 43.6 182 steady state does 176 54.7 69.1 106.1 460 89.9 29.6 37.4 56.6 234 steady state 0.3 0.332 0.346 4.49 5.89 7.53 not lose stability 5.58 9.36 14.4 4.03 5.97 8.34 D3 0.36 0.368 13.0 94.0 32.2 16.6 284 135 and D3 lose their stability, for some values of parameter δ, in the case of discrete delay, as well as the critical average delay at which the steady states D1 and D3 lose their stability for the case of the Erlang delay distribution, both in the vessel formation process, while the process of tumour growth is instantaneous. The average delay is calculated as m/a. These results are also illustrated in Fig. 1. Numerical simulations suggest that if δ converges to the limits of the interval of the existence of the steady state D1 , (average) critical value of delay tends to +∞. Similar behaviour is observed for critical value of time delay for the steady state D3 , however for this state the interval of existence is bounded only from the right-hand side. It can be also noticed that for m = 1, the positive steady states do not change their 12 EMAD ATTIA, MAREK BODNAR AND URSZULA FORYŚ D3 D1 m=2 m=3 m=5 discrete 40 τcr τcr 100 m=2 m=3 m=5 discrete 20 50 0.33 0.34 0.35 0.36 0.37 0.3 0.32 δ 0.34 0.36 δ Figure 1. Critical average delay, that is m/acr for various values of m in the dependance on δ in the case when only the process of tumour growth is delayed; left – graphs for the steady state D1 , right – graphs for the steady state D3 . effective vessel density stability, and this result is similar to those obtained when the tumour growth is delayed according to the Erlang distribution with m = 1. However, in the latter case the stability does not depend on other model parameters, while in the first case it depends on the model parameters. Numerical simulations suggest that for δ close to 0.346 the critical value of delay for the steady state D1 is the smallest one while it is larger for other δ. The regions below the curves in Fig. 1 are stability regions. For the state D1 stability regions for Erlang distributions is larger than in the discrete case, and decreases with increasing m. On the other hand, for the steady state D3 , the critical value of delay seems to be an increasing function of δ. For some values of δ the region of stability in the case of m = 2 is more than twice as large as in the discrete case. However, numerical simulations suggest, that for δ close to 0.3, the stability region for the Erlang distribution with m = 5 is smaller than for the discrete case, and this is somehow unexpected result. In Fig. 2 we see exemplary solutions of system (1.1) for parameters given by (3.1) and τ = 10. For m = 1 the solution converges to the steady state very fast, while for m = 5 it seems that oscillations are sustained. tumour volume 3 2.5 2 m=1 discrete delay 1.5 m=5 3 2 1 1 0 100 200 time 300 400 0 100 200 time 300 400 Figure 2. Solutions of system (1.1) for parameters given by (3.1) and τ = 10, with time delay present only in the vessel formation term. For comparison, we present solutions of system (1.1) with Erlang distributed delay with parameters m = 1 and σ = 0 or discrete delay present only in the tumour growth process. In this case Theorem 2.5 shows that the steady state is ANGIOGENESIS MODEL WITH ERLANG DISTRIBUTED DELAYS 13 effective vessel density tumour volume stable independently of the model parameters. This is seen in Fig. 3. Eventually, we present the case of both processes delayed with the same Erlang distribution in Fig. 4. 2.7 2.6 m=1 discrete delay 2.5 0 20 40 60 m=5 80 2.1 2 1.9 1.8 100 0 20 40 time 60 80 100 time effective vessel density tumour volume Figure 3. Solutions of system (1.1) for parameters given by (3.1) and τ = 10, with time delay present only in the tumour growth term. Here, for Erlang distribution, the steady state is stable, and solutions for m = 1 and m = 5 are almost identical. 3 2 m=1 discrete delay 1 m=5 15 10 5 0 0 20 40 60 time 80 100 0 20 40 60 80 100 time Figure 4. Solutions of system (1.1) for parameters given by (3.1) and τ = 5, with time delay present in both terms. From our numerical analysis it is clear that the most robust is the model with Erlang distribution with m = 1, while the most sensitive is the model with discrete delays. Moreover, in the case when the delay is present only in the tumour growth process, the influence of the magnitude of m is almost not visible. As we expect, for both delays being present, the behaviour of solutions is much more unpredictable. However, the result for m = 1 seems to be optimistic. Clearly, Erlang distribution with m = 1 is the most reasonable, and therefore this means that in reality oscillatory behaviour should be exception not rule. 4. Conclusions and discussion. In this paper a model of tumour angiogenesis with distributed delays was considered. We proved basic mathematical properties of the model showing that solutions are unique, non-negative and well defined on the whole positive half-line. We formulated conditions on the model parameters that guarantee lack of change of local stability of a steady state for any distribution of delays. On the other hand, we proved condition under which stability change can take place and Hopf bifurcation occurs. We gave more strict conditions in the case when delays are distributed according to Erlang distributions. Our results indicate that the model with distributed delays is more stable than with discrete ones. In 14 EMAD ATTIA, MAREK BODNAR AND URSZULA FORYŚ particular, we observe stabilisation of the solution in a steady state value for some delay distributions while in the same time solutions of the model with discrete delays exhibit oscillations. In the case of Erlang distributions we observe that the behaviour of the solution for large shape parameter is closer to the behaviour of the solution to the model with discrete delays. The model considered in this paper is an extension of the model proposed earlier by Agur et. al. [2]. In this paper Agur et al. tried to simplified more complex computer model of angiogenesis process proposed in [4]. However, this model always exhibits oscillatory dynamics, which is not realistic. On the other hand, according to the data presented in [4], such type of the dynamics should be also present in the model, and therefore they included time delays into their model. Our idea was to combine the properties of the Hahnfeldt et al. model with the properties of the one presented in [2]. Here we decided to use distributed delays instead of discrete ones in order to make the model more realistic comparing to the previous discrete case [7]. Our results show that both type of the model dynamics could be observed for the model with delays distributed according to Erlang distributions, depending on the shape parameter, which is good from the point of view of potential applications. Although we have not validated our model with experimental data, it is done for a small modification of this model and the results should be published shortly; cf. [11]. REFERENCES [1] T. Alarcón, M. R. Owen, H. M. Byrne and P. K. Maini, Multiscale modelling of tumour growth and therapy: The influence of vessel normalisation on chemotherapy, Computational and Mathematical Methods in Medicine, 7 (2006), 85–119. [2] Z. Agur, L. Arakelyan, P. Daugulis and Y. Ginosar, Hopf point analysis for angiogenesis models, Discrete Contin. Dyn. Syst. B , 4 (2004), 29–38. [3] A. R. A. Anderson and M. A. J. Chaplain, Continuous and discrete mathematical models of tumour-induced angiogenesis, Bull. Math. Biol., 60 (1998), 857–899. [4] L. Arakelyan, V. Vainstein and Z. Agur, 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, 5 (2002), 203–214. [5] E. Attia, M. Bodnar and U. Foryś, Angiogenesis model with erlang distributed delay in vessels formation, Proceedings of XIX National Conference on Application of Mathematics in Biology and Medicine (Regietów), ed. Faculty of Mathematics, Informatics and Mechanics, University of Warsaw, September 2015, pp. 9–14. [6] K. Bartha and H. Rieger, Vascular network remodeling via vessel cooption, regression and growth in tumors, Journal of Theoretical Biology, 241 (2006), 903–918. [7] M. Bodnar and U. Foryś, Angiogenesis model with carrying capacity depending on vessel density, J. Biol. Sys., 17 (2009), 1–25. [8] M. Bodnar and U. Foryś, Influence of time delays on the Hahnfeldt et al. angiogenesis model dynamics, Appl. Math. (Warsaw), 36 (2009), 251–262. [9] M. Bodnar, M. J. Piotrowska, U. Foryś and E. Nizińska, Model of tumour angiogenesis – analysis of stability with respect to delays, Math. Biosci. Eng., 10 (2013), 19–35. [10] M. Bodnar and M. J. Piotrowska, Stability analysis of the family of tumour angiogenesis models with distributed time delays, Communications in Nonlinear Science and Numerical Simulation, 31 (2016), 124–142. [11] M. Bodnar, P. Guerrero, R. Perez-Carrasco and M. J. Piotrowska, Deterministic and stochastic study for a microscopic angiogenesis model: Applications to the Lewis Lung carcinoma, PLoS ONE , 11 (2016), e0155553. [12] H. M. Byrne, M. A. J. Chaplain, Growth of non-nerotic tumours in the presence and absence of inhibitors, Math. Biosci., 130 (1995), 151–181. [13] K. L. Cooke and P. van den Driessche, On zeroes of some transcendental equations, Funkcj. Ekvacioj, 29 (1986), 77–90. [14] A. d’Onofrio and A. Gandolfi, Tumour eradication by antiangiogenic therapy: Analysis and extensions of the model by Hahnfeldt et al. (1999), Math. Biosci., 191 (2004), 159–184. ANGIOGENESIS MODEL WITH ERLANG DISTRIBUTED DELAYS 15 [15] J. Folkman, Tumor angiogenesis: therapeutic implications, N. Engl. J. Med., 285, (1971), 1182–1186. [16] U. Foryś and M. J. Piotrowska, Analysis of the Hopf bifurcation for the family of angiogenesis models II: The case of two nonzero unequal delays, Appl. Math. and Comp., 220 (2013), 277–295. [17] P. Hahnfeldt, D. Panigrahy, J. Folkman and L. Hlatky, Tumor development under angiogenic signaling: A dynamical theory of tumor growth, treatment response, and postvascular dormancy, Cancer Res., 59 (1999), 4770–4775. [18] Y. Hino, S. Murakami and T. Naito, Functional Differential Equations with Infinite Delay, Lecture Notes in Mathematics, vol. 1473, Springer-Verlag, New York, 1991. [19] M. E. Orme and M. A. J. Chaplain, Two-dimensional models of tumour angiogenesis and anti-angiogenesis strategies, IMA J. Math. Appl. Med. Biol., 14 (1997), 189–205. [20] M. J. Piotrowska and U. Foryś, Analysis of the Hopf bifurcation for the family of angiogenesis models, J. Math. Anal. Appl., 382 (2011), 180–203. [21] J. Poleszczuk, M. Bodnar and U. Foryś, New approach to modeling of antiangiogenic treatment on the basis of Hahnfeldt et al. model, Math. Biosci. Eng., 8 (2011), 591–603. [22] M. Scianna, C.G. Bell, and L. Preziosi, A review of mathematical models for the formation of vascular networks Journal of Theoretical Biology, 333 (2013), 174–209. [23] M. Welter, K. Bartha and H. Rieger, Vascular remodelling of an arterio-venous blood vessel network during solid tumour growth, Journal of Theoretical Biology, 259 (2009), 405–422. Received November 23, 2015; Accepted April 06, 2016. E-mail address: emadr@du.edu.eg E-mail address: mbodnar@mimuw.edu.pl E-mail address: urszula@mimuw.edu.pl