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 mm$^3$ 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 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; c.f.[12]. Many papers focusing on chemotaxis modelling was published by M. Chaplain and coauthors (c.f.[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 (c.f.[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

$\begin{split} \dot{{N}}(t) &= \alpha N(t) \biggl( 1-\frac{N(t)}{1+f_1(h_1(E_t)) } \biggr) ,\\ \dot{{P}}(t) &= f_2(E(t)) N(t) - \delta P(t),\\ \dot{{E}}(t) &= \Biggl(f_3(h_2(P_t)) ds - \alpha \biggl( 1-\frac{N(t)}{1+f_1(h_1(E_t))} \biggr) \Biggr) \ E(t), \end{split}$

(1.1)

where $N(t)$, $P(t)$, $E(t)$, $\alpha$ and $\delta$ 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 $h_i:\mathbf{C}((-\infty,0] , \mathbb{R})$, $i=1$, $2$, are defined as follows

$h_i(\phi) = \int_0^{\infty} k_i(s)\phi(-s){\rm{d}} s,$

and we assume $k_i(s): [0,\infty) \rightarrow [0,\infty)$ are probability densities with finite expectations, that is

$\int_0^\infty k_i(s){\rm{d}} s=1 \;\; \text{and}\;\; 0<\int_0^\infty s k_i(s) {\rm{d}} s<\infty, i=1, \, 2.$

For any $\phi\in \mathbf{C}(\mathbb{R} , \mathbb{R})$, the function $\phi_t$ is defined as $\phi_t(s) = \phi(t+s)$ for $s\in (-\infty,0]$.

Moreover, we assume that the functions $f_i$ are locally Lipschitz and there exist positive constants $A_2,\ A_3,\ B_1,\ B_3$ and $m_3$ such that

(A1) $f_1$ is increasing, $f_1(0)=0$ and $\lim_{E\rightarrow+\infty}f_1(E) =B_1>0$,

(A2) $f_2$ is decreasing and convex, $f_2(0)=A_2>0$ and $\lim_{E\rightarrow+\infty}f_2(E) =0$,

(A3)$f_3$ is increasing, $f_3(0)=-A_3<0$, $f_3(m_3)=0$ and $\lim_{P\rightarrow+\infty}f_3(P)=B_3$.

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 $\mathcal{C}=\mathbf{C}((-\infty,0] , \mathbb{R}^3)$ be the space of continuous functions defined on the interval $(-\infty,0]$ with values in $\mathbb{R}^3$, and as $\mathcal{C}_+$ we define the subspace of $\mathcal{C}$ that consist of the functions with 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 $\Phi$ in the following way

$\begin{split} \Phi&=\left \{ \phi \in \mathcal{C}: \lim\limits_{s \rightarrow -\infty} \phi(s){\rm{e}}^s=0 \mbox{and} \sup\limits_{s \in (-\infty,0]} |\phi(s) {\rm{e}}^s|<\infty \right \}, \\ \|\phi\|_\Phi &=\sup\limits_{s \in (-\infty,0]} |\phi(s) {\rm{e}}^s|, \end{split}$

and we consider initial functions from the set $\Phi_+=\Phi\cap\mathcal{C}_+$.

---

2. Analysis of the model

In this section we consider basic properties of system (1.1) for general distributions $k_i(s)$ and study stability in the case of Erlang distributions. For all results presented here we assume the functions $f_i$, $i=1,2,3$, fulfil conditions (A1)-(A3).

Theorem 2.1.  For any arbitrary initial function $\phi=(\phi_N,\phi_P,\phi_E) \in \Phi_+$ there exists a unique solution of system (1.1) in $\Phi_+$ defined on $t \in [0,+\infty)$. Moreover, for all $t \geq 0$ any solution satisfies the following inequalities

$\begin{split} N_{\min{}} \leq N(t) &\leq N_{\max{}}, \\ 0 \leq P(t) &\leq \max \left\{\frac{A_2}{\delta}N_{\max{}},\phi_2(0) \right\},\\ 0 \leq E(t) &\leq \phi_3(0) \exp\big((B_3+\alpha(N_{\max{}}-1))t\big), \end{split}$

(2.1)

where

$N_{\min{}} = \min \left\{1,\phi_N(0)\right\}, N_{\max{}}=\max \left\{\phi_N(0),1+B_1\right\}.$

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~= (0, 0, 0) \text{and} B = \Bigl(1, \frac{A_2}{\delta}, 0\Bigr) ,$

compare [9]. Moreover, there can exist positive steady states $D_i = (\bar{N_i}, m_3, \bar{E_i})$, with $\bar{N_i}=1+f_1(\bar{E_i})$, and $\bar{E_i}$ are zeros of the function

$g(x)=f_2(x)(1+f_1(x))-\delta m_3 .$

(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 $A_2 < \delta m_3$ and unstable for $A_2 > \delta m_3$.

---

2.1. Stability of positive steady states

In this section, we focus on examining the stability of the positive steady states $D_i$ 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

$K_i(\lambda) = \int_0^\infty k_i(s) {\rm{e}}^{-\lambda s} {\rm{d}} s .$

(2.3)

Then, the stability matrix of system (1.1) for the steady state $D_i$ reads

$M(\bar{N},\bar{P},\bar{E})= \left(\begin{matrix} -\alpha - \lambda&0 &\alpha d_1 K_1(\lambda)\\ f_2(\bar{E_i})&-\delta-\lambda&- \bar{N_i}\ d_2\\ \alpha b \bar{E_i}&d_3 \bar{E_i} K_2(\lambda) & -\alpha b d_1 \bar{E_i} K_1(\lambda) -\lambda \end{matrix}\right),$

where

$b=\frac{1}{1+f_1(\bar E_i)}, d_1 = f_1'(\bar E_i)>0, d_2 = -f_2'(\bar E_i)>0, d_3 = f_3'(m_3)>0 .$

Hence, the characteristic quasi-polynomial has the form

$\begin{split} W(\lambda)=&\lambda^3+C_1\lambda^2+C_2\lambda+(\lambda^2+\delta \lambda)C_3 K_1(\lambda)+\\ &+(\lambda+\alpha)C_4K_2(\lambda) -C_3C_5K_1(\lambda)K_2(\lambda), \end{split}$

(2.4)

with

$C_1=\alpha+\delta, C_2=\alpha \delta, C_3=\alpha \beta d_1, C_4=\frac{\beta d_2 d_3}{b^2}, C_5=\delta d_3 m_3, \beta= b\bar{E_i}.$

Conditions (A1)-(A3) guarantee positivity of $b$, $\beta$ and $d_i$, $i=1$, $2$, $3$. Consequently, $C_i>0$ for $i=1,\ldots,5$.

Theorem 2.3.  If $g'(\bar{E_i}) > 0$, where $g$ is given by (2.2), then the positive steady state $D_i=(\bar{N_i}, m_3, \bar{E_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)=\alpha C_4-C_3 C_5$ is reverse to the sign of $g'(\bar{E_i})$. Therefore, the assumption $g'(\bar{E_i}) > 0$ implies that $W(0)<0$. Further, it is easy to see that $W(\lambda)\rightarrow +\infty$ as $\lambda \rightarrow +\infty$. Then there exists at least one real positive eigenvalue, and this implies that $D_i$ 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 $k_2(s)=\delta_D(s)$, where $\delta_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

$w_1(\omega) = -\omega^3-C_3 \omega^2+(C_2+C_4-\delta C_3 )\omega- C_3C_5,$

(2.5)

$w_3(\omega) = -(C_1+C_3)\omega^2-\delta C_3 \omega+\alpha C_4-C_3C_5 ,$

(2.6)

and

$w_4(\omega) = -(C_1-C_3)\omega^2+\delta C_3 \omega+\alpha C_4+C_3C_5.$

(2.7)

To obtain positive zeros of $w_1$ one needs to assume $C_2+C_4-\delta C_3 >0$. Moreover, the value at positive extremum, that is at $P_{\max}=\frac{-2C_3 + \sqrt{4C_3^2 +12 (C_2+C_4-\delta C_3)}}{6}$, must be positive. Hence, it is necessary that $w_1(P_{\max})>0$. In such a case there exist two positive zeros $P_1<P_2$ of $w_1$. Next, we easily see that if $g'(\bar E_i)<0$, that is $\alpha C_4-C_3C_5>0$, then $w_3$ has exactly one positive zero $P_3 = \frac{-\delta C_3+\sqrt{\delta^2C_3^2-4(C_1+C_3)(\alpha C_4-C_3C_5)}}{2(C_1+C_3)}$. Similarly, if $C_1>C_3$, then $w_4$ has exactly one positive zero

$P_4 = \frac{\delta C_3+ \sqrt{\delta^2C_3^2+4(C_1-C_3)(\alpha C_4+C_3C_2)}}{2(C_1-C_3)}.$

In the following, we require $P_1<P_3$ and $P_4<P_2$.

Theorem 2.4.   Assume that $g'(\bar{E_i}) < 0$ and $D_i = (\bar{N_i}, m_3, \bar{E_i})$ is a positive steady state of system (1.1). Assume also that the delay distribution $k_2$ is concentrated at zero, that is $k_2(s)=\delta_D(s)$. If

(ⅰ)${\beta}{d_{{1}}}<1$ and

$\frac{2 \alpha^2 b}{d_2 \bar{E_i} }<d_3< \min\left \{\frac{b}{2d_2 \bar{E_i}}\left(\delta^2-{\alpha}^{2}\left({\beta}^{2}{d_{{1}}}^{2}-1\right)\right),\ \frac{\delta}{2{\beta}^{2}{d_{{1}}}^{2} m_{{3}}} \left(1-{\beta}^{2}{d_{{1}}}^{2}\right)\right\},$

or

(ⅱ) $C_1>C_3$, $w_1(P_{\max})>0$ and $w_1(P_i)>0$ for $i=3$, $4$, where $w_1$ is defined by (2.5) and $P_3$, $P_4$ are positive zeros of (2.6) and (2.7), respectively,

then $D_i$ is locally asymptotically stable for any distribution $k_1$.

Proof. The steady state $D_i$ is stable in the case without delay for $g'(\bar E_i)<0$. Hence, we need to show that the characteristic function has no purely imaginary zeros. Therefore, we investigate the behaviour of $W(i\omega)$.

Let

$K_1(i \omega)=\eta_1-i\zeta_1, \eta_1=\int\limits_0^\infty k_1(s) \cos(\omega s){\rm{d}} s, \ \ \zeta_1= \int\limits_0^\infty k_1(s) \sin(\omega s){\rm{d}} s.$

Thus

$W(i\omega)=-i\omega^3-C_1\omega^2+i(C_2+C_4)\omega+((-\omega^2+i\delta \omega)C_3 -C_3C_5)(\eta_1-i\zeta_1)+\alpha C_4 ,$

and we have

$\begin{array}{r@{}l} \mathfrak{Re}(W(i\omega))&{}=-C_1\omega^2+\alpha C_4 - C_3 ( \omega^2+ C_5 ) \eta_1+\delta C_3 \omega \zeta_1 ,\\\\ \mathfrak{Im}(W(i\omega))&{}=-\omega^3+(C_2+C_4)\omega +\delta C_3 \omega \eta_1 + C_3 ( \omega^2+ C_5)\zeta_1. \end{array}$

(2.8)

Assume that there exists $\omega>0$ such that $W(i\omega)=0$, then

$\begin{split} & \left(-C_1\omega^2+\alpha C_4\right)^2+\left(-\omega^3+(C_2+C_4)\omega\right)^2 = \\ & \left( -C_3 ( \omega^2+ C_5 ) \eta_1+\delta C_3 \omega \zeta_1 \right)^2 +\left( \delta C_3 \omega \eta_1+ C_3 ( \omega^2+ C_5)\zeta_1 \right)^2. \end{split}$

(2.9)

For Eq. (2.9) we have

$\begin{split} L.H.S& = \omega^6+ \left( C_1^2 -2(C_2+C_4) \right) \omega^4 +\left( (C_2+C_4)^2 -2 \alpha C_1 C_4 \right) \omega^2 +\alpha^2 C_4^2 , \\ R.H.S&= C_3^2 \left( \omega^4+\left( \delta^2 +2 C_5 \right) \omega^2 +C_5^2 \right) \left(\zeta_1^2+\eta_1^2\right). \end{split}$

Schwarz inequality yields

$\begin{split} \left(\int\limits_0^{\infty} \cos(\omega s) k_1(s) {\rm{d}} s\right)^2 & = \left(\int\limits_0^{\infty} \cos(\omega s) {\rm{d}} \left(\int_0^s k_1(u) {\rm{d}} u\right) \right)^2 \\ &\leq \int\limits_0^{\infty} \cos^2(\omega s) {\rm{d}} \left(\int_0^s k_1(u) {\rm{d}} u\right) \int\limits_0^{\infty} {\rm{d}} \left(\int_0^s k_1(u) {\rm{d}} u\right)\\ &= \int\limits_0^{\infty} \cos^2(\omega s) k_1(s) {\rm{d}} s,\\ \left(\int\limits_0^{\infty} \sin(\omega s) k_1(s) {\rm{d}} s\right)^2 & \leq \int\limits_0^{\infty} \sin^2(\omega s) k_1(s) {\rm{d}} s. \end{split}$

Consequently, $\zeta_1^2+\eta_1^2\leq 1$. Then, $\omega$ satisfying $W(i\omega)=0$ should fulfil

$0=L.R.S.-R.H.S.$

$\begin{split} \omega^6+ &\left( C_1^2 -2\tilde C_{2,4} \right) \omega^4 +\left( \left(\tilde C_{2,4} \right)^2 -2 \alpha C_1 C_4 \right) \omega^2 + \\ & +\alpha^2 C_4^2 - C_3^2 \left( \omega^4+\left( \delta^2 +2 C_5 \right) \omega^2 +C_5^2 \right) , \end{split}$

where $\tilde C_{2,4} = C_2+C_4$. Let us denote $y=\omega^2$ and define the auxiliary function

$\begin{split} F(y) &= y^3+\left(C_1^2-2\tilde C_{2,4}-C_3^2\right) y^2+ \left(\left(\tilde C_{2,4}\right)^2-2 \alpha C_1 C_4 - C_3^2\Bigl(\delta^2 -2 C_5\Bigr) \right)y +\\ &+\alpha^2 C_4^2-C_3^2C_5^2 , \end{split}$

Existence of purely imaginary eigenvalue requires $F(y) \leq 0$ for some positive $y$. However, all coefficients of this polynomial are positive under the assumptions, and then purely imaginary eigenvalue $i\sqrt{y}$ does not exist.

Clearly, the free term is positive due to $g'(\bar E_i)<0$. Moreover, inequalities (i) are equivalent to

${\beta}{d_{{1}}}<1, \alpha^2 < {\frac {\beta\,d_{{2}}d_{{3}}}{2{b}^{2}}} \;\;\; \mbox{and }$

${\delta}>\max \left\{\sqrt{{\alpha}^{2}({\beta}^{2}{d_{{1}}}^{2}-1)+{\frac {2 \beta\,d_{{2}}d_{{3}}}{{b}^{2}}}},\frac{2{\beta}^{2}{d_{{1}}}^{2}d_{{3}}m_{{3}}}{1-{\beta}^{2}{d_{{1}}}^{2} }\right\}.$

Therefore,

$C_1^2-2(C_2+C_4)-C_3^2={\delta}^{2}- {\alpha}^{2} \left( {\beta}^{2}{d_{{1}}}^{2}-1 \right) -{\frac {2 \beta\,d_{{2}}d_{{3}}}{{b}^{2}}} >0 ,$

$\;\;\;\;\;(C_2+C_4)^2-2 \alpha C_1 C_4 -\delta^2 C_3^2 -2 C_3^2 C_5= \\ \left( \alpha\,\delta+{\frac {\beta\,d_{{2}}d_{{3}}}{{b}^{2}}} \right) ^{2} -{\frac {2 \alpha\, \left( \delta+\alpha \right) \beta\, d_{{2}}d_{{3}}}{{b}^{2}}} - {\delta}^{2}{\alpha}^{2}{\beta}^{2}{d_{{1}}}^{2}-2{\alpha}^{2}{\beta}^{2}{d_{{1}}}^{2}\delta\,d_{{3}}m_{{3}} = \\\left( \left( 1-{\beta}^{2}{d_{{1}}}^{2} \right) {\alpha}^{2} \delta-2 {\alpha}^{2}{\beta}^{2}{d_{{1}}}^{2}d_{{3}}m_{{3}}\right) \delta+\frac{\beta\,d_{{2}}d_{{3}}}{b^2} \left( {\frac {\beta\,d_{{2}}d_{{3}}}{{b}^{2}}}- 2{\alpha}^{2} \right) >0 .$

Due to the continuous dependance of eigenvalues on the model parameters this completes the proof of the first part.

For the proof of the second part, notice that (2.8) implies

$\begin{split} \mathfrak{Re}(W(i\omega)) \geq&-(C_1+ C_3) \omega^2-\delta C_3 \omega+\alpha C_4-C_3C_5,\\ \mathfrak{Re}(W(i\omega)) \leq&-(C_1- C_3) \omega^2+\delta C_3 \omega+\alpha C_4+C_3C_5,\\ \mathfrak{Im}(W(i\omega)) \geq&-\omega^3-C_3 \omega^2+(C_2+C_4-\delta C_3 )\omega- C_3C_5. \end{split}$

It is clear that,

$\begin{array}{l} \mathfrak{Re}(W(i\omega))&>&0 \ \mbox{for } \omega \in [0,P_3), \mbox{and } \mathfrak{Im}(W(i\omega))>0 \ \mbox{for } \omega \in (P_1,P_2),\\ \mathfrak{Re}(W(i\omega))&<&0 \ \mbox{for } \omega \in (P_4,\infty). \end{array}$

Inequalities $w_1(P_i)>0$ for $i=3$, $4$, imply $P_1<P_3$ and $P_2>P_4$. Moreover, $w_3(0)<w_4(0)$ and $w_3'(\omega)<w_4'(\omega)$ for $\omega>0$. Therefore $P_3<P_4$ and there is no $\omega\geq 0$ such that $W(i\omega)=0$. This completes the proof of the theorem.

Remark 1.  If the second condition of assumption (ⅰ) of Theorem 2.4 is satisfied, then

$\delta^2>3\alpha^2 .$

If this inequality holds, we can choose sufficiently small $d_1$ and later we can take a proper value of $d_3$ 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

$k(s)=\begin{cases} \displaystyle\frac{{a}^{m}(s-\sigma)^{m-1}}{(m-1)!} {\rm{e}}^{-{a}(s-\sigma)},& s\geq \sigma,\\ 0,&\textrm{otherwise}, \end{cases}$

(2.10)

where $a$, $\sigma>0$, and $a$ is a scaling parameter. To the case $\sigma=0$ we refer as to the non-shifted Erlang distribution, while to the case $\sigma>0$ we refer as to the shifted one. The mean value of the non-shifted Erlang distribution is given by $\frac{m}{a}$ and the variance is equal to $\frac{m}{a^2}$. 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 $\sigma+\frac{m}{a}$ and the variance is the same as for the non-shifted one. On the other hand, taking the limit $m\to+\infty$ and keeping $m/a=\tau$ constant we obtain system (1.1) with discrete delay $\tau$. By direct calculation, it is found that $\int_0^\infty k(s) {\rm{e}}^{-\lambda s} {\rm{d}} s=\frac{a^{m}}{(a+\lambda)^{m}} {\rm{e}}^{-\lambda \sigma}$.

Now, let us consider the first distribution to be Erlang, $k_1(s)=k(s)$, while the other is still $k_2(s)=\delta_D(s)$. In this case, if $C_5 \ne a(a-\delta)$, then the characteristic function (2.4) could be replaced by

$\begin{split} W_I(\lambda) &= (a+\lambda)^m\Bigl(\lambda^3+C_1\lambda^2 +\bigl(C_2+C_4\bigr)\lambda +\alpha C_4\Bigr)+ \\& +a^m\Bigr(C_3\lambda^2+\delta C_3\lambda -C_3C_5\Bigl){\rm{e}}^{-\lambda\sigma}. \end{split}$

(2.11)

As the case $C_5 = a(a-\delta)$ is non-generic, in the following we assume $C_5 \ne a(a-\delta)$, and consider the simplest possible Erlang distribution with $m=1$ and $\sigma=0$.

Theorem 2.5.  If $m=1$, $\sigma=0$ and $g'(\bar{E_i}) < 0$, then the positive steady state $D_i$ is locally asymptotically stable.

Proof. For this case the characteristic function $W_I$ simplifies to

$W_I(\lambda) = (a+\lambda)\Bigl(\lambda^3+C_1\lambda^2 +\bigl(C_2+C_4\bigr)\lambda +\alpha C_4\Bigr) +a\Bigr(C_3\lambda^2+\delta C_3\lambda -C_3C_5\Bigl).$

The Routh-Hurwitz Criterion for $W_I$ reads

$q_1 q_2 q_3>q_3^2+q_1^2 q_4 ,$

(2.12)

where

$q_1 = a+C_1, q_2 = a(C_1+C_3)+\eta_2, q_3 = a(\eta_2+\delta C_3) + \alpha C_4, q_4 = a\eta_4, \ $

and

$\eta_2 = C_2+C_4, \eta_4 = \alpha C_4-C_3C_5 ,$

Notice that inequality (2.12) is equivalent to

$P(a) = a_3 a^3 + a_2 a^2 + a_1 a~+ a_0 >0,$

where

$\begin{split} a_3 &= (C_1 + C_3) (\eta_2 + \delta C_3) - \eta_4 ,\\ a_2 &= \alpha C_4 (C_1 + C_3) + C_2 (\eta_2 + \delta C_3) + C_1 (C_1 + C_3) (\eta_2 + \delta C_3)\\ & \ \ \ \ + C_4 (\eta_2 + \delta C_3 ) - (\eta_2 +\delta C_3 )^2 - 2 C_1 \eta_4, \\ a_1 &= C_1 (\eta_2^2 + \alpha C_3 C_4 + \delta C_3\eta_2 ) -\alpha C_4 (\eta_2 + 2\delta C_3 ) + C_1^2 C_3C_5, \\ a_0 &= \alpha C_4 (C_1\eta_2 - \alpha C_4). \end{split}$

Due to the definitions of $C_1$ and $\eta_i$ we have:

$\begin{split} a_0&= \alpha C_4 \Bigl(\bigl(\delta+\alpha\bigr)\bigl(C_2+C_4\bigr) - \alpha C_4\Bigr) = \alpha C_4 \Bigl( \bigl(\delta+\alpha\bigr) C_2 + \delta C_4 \Bigr)>0, \\ a_1&= (\alpha+\delta) \Bigl((C_2+C_4)^2 + \alpha C_3 C_4 + \delta C_3(C_2+C_4) \Bigr) -\alpha C_4 (C_2+C_4) +\\ & -2\alpha \delta C_3 C_4 + C_1^2 C_3C_5 \\ &= \alpha\Bigl(\eta_2C_2+\alpha C_3C_4+\delta C_2C_3\Bigr) + \delta\eta_2\bigl(\eta_2+\delta C_3\bigr)>0 , \\ a_2 &= \alpha C_4 C_3 + C_2 (\eta_2 + \delta C_3) + (C_1(\alpha+\delta) + C_1C_3) (\eta_2 + \delta C_3) + C_4 (\eta_2 + \delta C_3 )+\\ & - (\eta_2 +\delta C_3 )^2 - \alpha C_1 C_4+2C_3C_5, \\ &=\alpha C_4 C_3 + (\delta C_1 + (\alpha+\delta)C_3) (\eta_2 + \delta C_3) +\alpha C_1 (C_2+ \delta C_3) +\\ & - \delta C_2 C_3- \delta C_4 C_3 -\delta^2 C_3^2 +2C_3C_5, \\ &=\alpha C_4 C_3 + (\delta C_1 + \alpha )C_3 (\eta_2 + \delta C_3) +\alpha C_1 (C_2+ \delta C_3) +2C_3C_5>0,\\ a_3&= (\delta+\alpha + C_3) (C_2+C_4 + \delta C_3) - \alpha C_4+C_3C_5 \\ &= (\delta+ C_3) (C_2+C_4 + \delta C_3) +\alpha (C_2+\delta C_3) +C_3C_5 > 0. \\ \end{split}$

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 $k_1(s)=\delta_D(s)$, while the second distribution is Erlang, that is $k_2(s)=k(s)$ described by Eq. (2.10). Note that if $\alpha C_4\neq C_4a+C_3C_5$, then $\lambda$ is the eigenvalue defined by (2.4) if and only if $\lambda$ is zero of

$\begin{split} W_{II}(\lambda) &= (a+\lambda)^{m}\Bigl(\lambda^3+\bigl(C_1+C_3\bigr)\lambda^2 +\bigl(C_2+\delta C_3\bigr)\lambda \Bigr) + \\ & +a^m\Bigr(C_4\lambda+\alpha C_4 -C_3C_5\Bigl){\rm{e}}^{-\lambda \sigma}. \end{split}$

(2.13)

Again, because the case $\alpha C_4= C_4a+C_3C_5$ is non-generic, we do not consider it here, and in the following we assume $\alpha C_4\neq C_4a+C_3C_5$. Thus, studying the stability of the positive steady states $D_i$ of system (1.1) is equivalent to study zeros of the polynomial $W_{II}$ defined by (2.13). Below, we again study the simplest case of Erlang distribution with $m=1$ and $\sigma=0$.

Proposition 2.6.  If $m=1$, $\sigma=0$, $g'(\bar{E_i}) < 0$ and $Q_1 Q_2 Q_3>Q_3^2+Q_1^2 Q_4$, where

$\begin{aligned} Q_1&=C_1+C_3+a,& Q_2&=C_2+\delta C_3+a (C_1+C_3), \\ Q_3&=a(C_2+C_4+\delta C_3), & Q_4&=a (\alpha C_4 -C_3C_5), \end{aligned}$

then the positive steady state $D_i$ is locally asymptotically stable.

Proof. For $m=1$ and $\sigma=0$, the polynomial $W_{II}$ defined by (2.13) reads

$\begin{split} W_{II}(\lambda)&= \lambda^4+(C_1+C_3+a)\lambda^3+(C_2+\delta C_3+a (C_1+C_3))\lambda^2 +\\& +a(C_2+C_4+\delta C_3)\lambda+a (\alpha C_4 -C_3C_5) , \end{split}$

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

$\eta_1 = C_1+C_3, \eta_2 = C_2+\delta C_3, \eta_4 = \alpha C_4-C_3C_5.$

With this notation we have

$Q_1=\eta_1+a, Q_2=\eta_2+a \eta_1, Q_3=a(\eta_2+C_4), Q_4 =a \eta_4,$

$Q_i>0 \ \ \text{for} \ \ i=1,\ldots ,4.$

Now, the Routh-Hurwitz condition is equivalent to

$\begin{split} &a^2\Bigl(\eta_1 \bigl(\eta_2+C_4\bigr)-\eta_4\Bigr) + a \Bigl(\bigl(\eta_2+C_4\bigr)\bigl(\eta_1^2-C_4\bigr)-2\eta_1\eta_4\Bigr)+\\ & +\eta_1\Bigl(\eta_2\bigl(\eta_2+C_4\bigr)-\eta_1\eta_4\Bigr) >0 . \end{split}$

(2.14)

Notice that the coefficient of $a^2$ is positive. Clearly, using the definitions of $\eta_1$, $\eta_4$ and $C_1$ we have

$\eta_1 \bigl(\eta_2+C_4\bigr)-\eta_4 = \eta_1\eta_2 + (\alpha+\delta+C_3)C_4 - \alpha C_4+C_3C_5 =\eta_1\eta_2+(\delta+C_3)C_4+C_3C_5>0.$

Now, we have only three possibilities:

1. $\eta_2(\eta_2+C_4)< \eta_1\eta_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.$\bigl(\eta_2+C_4\bigr)\bigl(\eta_1^2-C_4\bigr)/2 <\eta_1\eta_4 < \eta_2(\eta_2+C_4)$ and the discriminant of the quadratic polynomial is positive -there exist two critical values of $a$;

3.$\eta_1\eta_4 < \eta_2(\eta_2+1)$ and either $\eta_1\eta_4< \bigl(\eta_2+C_4\bigr)\bigl(\eta_1^2-C_4\bigr)/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

$\Bigl(\bigl(\eta_2+C_4\bigr)\bigl(\eta_1^2-C_4\bigr)-2\eta_1\eta_4\Bigr)^2>4 \eta_1\Bigl(\eta_2\bigl(\eta_2+C_4\bigr)-\eta_1\eta_4\Bigr)\Bigl(\eta_1 \bigl(\eta_2+C_4\bigr)-\eta_4\Bigr) ,$

(2.15)

together with

$\frac{(\eta_1^2-C_4)(\eta_2+C_4)}{2}<\eta_1\eta_4< {\eta_2(\eta_2+C_4)}.$

Inequality (2.15) is equivalent to

$\bigl(\eta_2+C_4\bigr)\bigl(\eta_1^2-C_4\bigr)^2 -4\eta_1\eta_4\bigl(\eta_1^2-C_4\bigr) -4 \eta_1^2\eta_2\bigl(\eta_2+C_4\bigr)+4\eta_1\eta_2\eta_4 +4\eta_1^3\eta_4 >0 ,$

and collecting terms with $\eta_1^2-C_4$ we obtain

$\bigl(\eta_2+C_4\bigr)\bigl(\eta_1^2-C_4\bigr)^2 -4\eta_1\eta_4\bigl(\eta_1^2-C_4\bigr) +4\eta_1 \Big( \eta_4 \bigl(\eta_2 +\eta_1^2\bigr) - \eta_1\eta_2\bigl(\eta_2+C_4\bigr)\Big) >0 .$

(2.16)

Notice that the free and linear terms of (2.16) are positive under the assumption

$\eta_4> \frac{ \eta_1\eta_2\bigl(\eta_2+C_4\bigr)}{\eta_2+\eta_1^2} \text{and} \eta_1^2<C_4 .$

We have $\eta_1^2-C_4=(\alpha+\delta+\alpha\delta)^2-\beta d_2d_3/b^2$, and it is negative for sufficiently large $d_2d_3$ or sufficiently small $b$.

Eventually, two stability switches are possible under the assumptions

$\frac{ \eta_1\eta_2\bigl(\eta_2+C_4\bigr)}{\eta_2+\eta_1^2} <\eta_4<\frac{\eta_2(\eta_2+C_4)}{\eta_1} \text{and} \eta_1^2<C_4 .$

Proposition 2.7.  If $g'(\bar E_i) \ne 0$ and

(ⅰ)$D_i$ is unstable for $\sigma=0$, then it is unstable for all $\sigma>0$;

(ⅱ)$D_i$ is stable for $\sigma=0$, then there exists $\sigma_{0}>0$, such that $D_i$ is stable for $\sigma<\sigma_{0}$, and $D_i$ is unstable for $\sigma> \sigma_{0}$. Furthermore, if $f_i \in \mathbf{C}^2, i~= 1, 2, 3$, then Hopf bifurcation occurs at $\sigma_{0}$.

Proof.  For the characteristic function $W_{II}$ given by (2.13), we define the auxiliary function

$\begin{split} F(y)&=y\bigl(a^2+y\bigr)^{m} \left(y^2+((C_1+C_3)^2-2(C_3\delta+C_2))y+(C_3\delta+C_2)^2\right)-\\ & -a^{2m} C_4^2y- a^{2m} (C_4\alpha -C_3C_5)^2 , \end{split}$

where $y=\omega^2$. Notice, that

$(C_1+C_3)^2-2(C_3\delta+C_2)=\alpha^2(1+\beta d_1)^2+\delta^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.

Proposition 2.8.   If $m=1$, $\sigma=0$, $g'(\bar{E_i}) < 0$, $q_{11} q_{22} q_{33}>q_{33}^2+q_{11}^2q_{44}$ and

$(q_{11}q_{44}-q_{55})(q_{11}q_{22}q_{33}-q_{33}^2-q_{11}^2q_{44})>q_{55}(q_{11}q_{22}-q_{33})^2+q_{11}q_{55}^2,$

where

$\begin{split} & q_{11}= C_1+2 a~, q_{22}=C_2+2 a~C_1+a^2+a C_3,\\ & q_{33}=2 a~C_2+a^2 C_1+a^2 C_3+a C_3 \delta+a C_4,\\ & q_{44}=a^2 C_2+a^2 C_3 \delta +a^2 C_4+a C_4 \alpha, q_{55}=a^2 C_4 \alpha-a^2C_3C_5 , \end{split}$

then the positive steady state $D_i$ is locally stable.

Proof. For $m=1$ and $\sigma=0$, the characteristic function (2.4) reads

$\begin{array}{l} W(\lambda)&=&\frac{1}{(a+\lambda)^2}\Big(\left(\lambda^3+C_1\lambda^2+C_2\lambda\right)(a+\lambda)^{2}+a(\lambda^2+\delta \lambda)C_3 (a+\lambda)+\\ &+&a(\lambda+\alpha)C_4(a+\lambda) -a^2C_3C_5\Big). \end{array}$

Hence, we need to study roots of the polynomial

$\lambda^5+(C_1+2 a~)\lambda^4+(C_2+2 a~C_1+a^2+a C_3)\lambda^3+(2 a~C_2+a^2 C_1+a^2 C_3+a C_3 \delta+a C_4)\lambda^2+\\ +(a^2 C_2+a^2 C_3 \delta +a^2 C_4+a C_4 \alpha) \lambda+a^2 C_4 \alpha-a^2C_3C_5=0.$

A direct application of the Routh-Hurwitz Criterion completes the proof.

---

3. Numerical simulations

For the numerical simulations we choose functions $f_i$ and parameters values proposed in[9], that is

$\begin{align*} f_1(E) &= \frac{b_1 E^n}{c_1+E^n}\,, & f_2(E) &= \frac{a_2c_2}{c_2+E}\,, & f_3(P) &= \frac{b_3(P^2-m_3^2)}{\frac{m_3^2b_3}{a_3}+P^2}\,, \end{align*}$

and

$a_2 = 0.4, \ \ a_3=1, \ \ b_1 = 2.3, \ \ b_3=1, \ \ c_1 = 1.5, \ \ c_2 = 1, \ \ \alpha = 1, \ \ \delta = 0.34.$

(3.1)

For these values of parameters there exist three positive steady states:

$D_1\approx(1.04, 1.05, 0.17), D_2\approx (1.37, 1.05, 0.54), D_3\approx (2.67, 1.05, 1.99) .$

Now, we can influence the model dynamics changing the value of $\delta$. 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 $D_1$ exists for $0.331<\delta<0.381$, $D_3$ for $\delta<0.368$ and $D_2$ for $0.331<\delta<0.368$. The steady state $D_2$ is unstable, while $D_1$ and $D_3$ are stable for the case without time delay.

In Table 1 we presented critical values of delay at which the steady states $D_1$ and $D_3$ lose their stability, for some values of parameter $\delta$, in the case of discrete delay, as well as the critical average delay at which the steady states $D_1$ and $D_3$ 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 $\delta$ converges to the limits of the interval of the existence of the steady state $D_1$, (average) critical value of delay tends to $+\infty$. Similar behaviour is observed for critical value of time delay for the steady state $D_3$, 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 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 $\delta$ close to $0.346$ the critical value of delay for the steady state $D_1$ is the smallest one while it is larger for other $\delta$. The regions below the curves in Fig. 1 are stability regions. For the state $D_1$ 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 $D_3$, the critical value of delay seems to be an increasing function of $\delta$. For some values of $\delta$ 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 $\delta$ 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 $\tau=10$. For $m=1$ the solution converges to the steady state very fast, while for $m=5$ it seems that oscillations are sustained.

For comparison, we present solutions of system (1.1) with Erlang distributed delay with parameters $m=1$ and $\sigma=0$ or discrete delay present only in the tumour growth process. In this case Theorem 2.5 shows that the steady state is 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.

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 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; c.f.[11].

---