# Stability Regions and Bifurcation Analysis of a Delayed Predator-Prey Model Caused from Gestation Period

Authors: Nikenasih Binatari; Fajar Adi-Kusumo; Lina Aryati
Journal: International Journal of Differential Equations
 (2022)

---

## Abstract

We consider a saturated predator-prey system with delayed time. The system is motivated by natural disease management using the interactions between the original and the treated species populations, such asAedes aegypti and Wolbachia mosquitoes, fertile and infertile pests as a pesticide’s effect, uninfected and infected cancer cells by an oncolytic virus, and so forth. The delayed time shows the gestation effect of the treated populations where the impact on the stability of the unique positive equilibrium point of the system will be studied. We obtain the exact formula of the equilibrium point where it is asymptotically stable for the nondelay case. The stability region of the nonzero solution is given in parameter space following the Pontryagin criteria. Furthermore, some conditions, such that for delay case this solution is conditionally stable, are also provided in this study.


---

## Body

## 1. Introduction

The need for food continues to increase in line with the increasing population and the higher level of domestic consumption. However, agricultural and plantation production yields can decrease or fail if the presence of pests around the land is not controlled. One technique to reduce the proportion of pests is through the application of Integrated Pest Management (IPM). Some techniques applied to suppress pest attacks are good agricultural practices, biological control, breeding and growth of resistant varieties, and the use of chemical insecticides [1].The most common pest control technique used by farmers is the use of chemical insecticides. It is due to the fact that the pesticides can be applied easily and quickly and can reach a large area. However, the use of pesticides can also reduce the quality of crop yields. In particular, in excessive use, the pesticides can harm farmers and the environment. On the other hand, biological control is the safest and the most effective technique to control the pest. This control method utilizes the predatory behavior as the natural pest control. However, this method is time-consuming and inefficient. Based on these facts, the analysis related to the dynamics of the prey population by considering its predator and the presence of pesticides needs to be studied.The interaction of predator and prey populations can be presented in a mathematical model, which was introduced by Lotka in a simple manner, in which the growth of the prey-predator population is assumed to be influenced only by the birth and the interaction of both populations [2]. Recently, one of the developments in modeling has been to involve the delayed time, such as on the predator density [3, 4], on both prey and predator density [5, 6], on the predator death rate [7], and on the harvesting time [8], where the existence of the delayed time could change the stability behavior of the system.The impact of the delayed time in pest controlling in farming has been studied by Abraha et al. [9], where the authors considered the delays in the level of awareness as a saturated term. In [10], the authors adopt the predator prey system to study the interaction of two populations of one species, which is original and treated populations with saturation term for nondelay case. Meanwhile, the gestation period also affects the growth rate of the original population. Mathematically, the gestation period is represented by the delayed time of the growth rate; see [11–13]. Therefore, in this paper, we include the gestation period of original population by adding a delayed time parameter to a saturated model in [10]. The objective of this research is to study existence and stability of the equilibrium point of the delayed time system. The study is important to understand the effect of the gestation period to have successful treatment.The stability of an equilibrium point in the system of ordinary or delay differential equations is usually determined by calculating the roots of the characteristic equation of the linearized system near an equilibrium point. An equilibrium point is called asymptotically stable if and only if the real parts of the roots of characteristic equation are negative [14]. The difference between the characteristic equation of the ordinary differential equations (ODE) and delay differential equation (DDE) systems is as follows. The characteristic equation of an ODE system is a polynomial, which usually has finite number of characteristic roots. Meanwhile, in DDE system, the characteristic equation is a quasi-polynomial [15] that has infinitely many characteristic roots [16, 17].Due to the complexity of the equation, the stability behavior of the delayed time system in some references was observed without computing the roots. The fundamental criterion for calculating the roots of quasi-polynomial equations was based on the work of Pontryagin [18] that gave necessary and sufficient conditions for a quasi-polynomial which has a principle term. Some other methods have been proposed to study the stability of exponential polynomials, such as Lyapunov’s direct method, Yesipovich-Svirskii criterion, τ-decomposition method, Chebotarev criterion, D-subdivision method, Nyquist criterion, Bode criterion, Nichols criterion, and Michailov criterion. In this paper, we consider a quasi-polynomial equation which has a principle term, so that the stability regions of the equation will be derived by applying Pontryagin Theorem.The stability regions of some cases of the second-order delay differential equations that use Pontryagin Theorem have been studied in [19–23]. In this paper, we obtain a new type of quasi-polynomial characteristic equation, that is, λ2−T¯λ+D¯−E¯e−λτ=0, where the methods to find the roots of the quasi-polynomial equation will be studied. It is important to determine the stability regions in the parameter space and the bifurcation analysis in terms of the delay time parameter.The remainder of this paper is organized as follows. In Section2, we start with model formulations and then continue with the main results in Section 3. In Subsection 3.1, we show the existence conditions of the positive equilibrium point. By applying the Pontryagin Theorem, we provide the sufficient and necessary condition in terms of the parameters and the delay parameter such that the positive equilibrium point is asymptotically stable; we show this in Section 3.2. In the next sections, we identify the bifurcation point and its conditions, and then finally some simulations are given to illustrate the results.


## 2. Problem Formulation

We consider a two-dimensional delay differential equation as follows:(1)dxdt=r1x−r1xx+yK−bxyx+y+a,dydt=r2y−r2yxτ+yK+bxyx+y+a−βy.System (1) is based on the predator prey model with logistic growth. The system is motivated by a treatment model of populations with one species which was studied in [10]. The model in [10] showed the interaction between the original population x and the treated population y for the delay parameter τ=0. The densities of the populations were assumed to follow the logistic growth with the intrinsic growth rates r1 and r2, respectively, and the carrying capacity K. The other assumption is that the interaction between both populations is able to reduce the original population that turns into the treated population. The existence of the nontrivial positive equilibrium point due to the parameter space and the stability criteria were studied in [10], but the explicit expression of the equilibrium point was not presented.The generalized model of System (1) in the sense of parameter space has been presented in the work of Adi-Kusumo et al. [24]. The authors in [24] introduced two new parameters that show the endurance of the populations and efficacy of the treatment was introduced. Moreover, although the expression of a nontrivial positive equilibrium point is provided, the stability behavior of the equilibrium point has not yet been studied.The existence of the nontrivial positive equilibrium points of the nondelay case of System (1) depends on the parameter space. In [10], there is a region in the parameter space where the nontrivial positive equilibrium point is unique. The appearance of this equilibrium point is important in the sense of applications. The stability of the equilibrium point shows the situation where the original population cannot be removed completely by the treatment.In this paper, we add the delay timeτ to the system in [10], where xτ=xt−τ; see System (1). The delayed term can be interpreted as the gestation period of the original population that affects the interaction with the treated populations. Our work focuses on understanding the effect of the delay term on the behavior of the nontrivial positive equilibrium point of System (1). Moreover, we will obtain the stability region and some conditions for the appearance of Hopf bifurcations.


## 3. Stability Analysis

### 3.1. The Equilibrium Points

Letℛ=ℝ+∪0 be a nonnegative real number set. The domain of System (1) is denoted by(2)Ω=x,y∈ℛ2|0<x+y≤K,where all the parameters r1,r2,K,b,a, and β are nonnegative. The equilibrium solutions of System (1) must satisfy(3)r1x−r1xx+yK−bxyx+y+a=0,(4)r2y−r2yxτ+yK+bxyx+y+a−βy=0.In the absence of delay, System (1) has four possible types of equilibrium points, that is,E00,0 as a trivial equilibrium, E1K,0 and E20,y¯ with y¯=β−r2/r2 as semitrivial equilibria, and E3x∗,y∗ as the nontrivial equilibrium point. The equilibrium points E0,E1, and E2 indicate the extinction of one or both populations.There is more than one positive nontrivial equilibrium point due to the parameter space. However, in this paper, we restrict our discussion to the parameter space that the system has a unique positive equilibrium point. The sufficient conditions of the uniqueness were given in [10]. Here, we provide the expression of the unique positive equilibrium point in Theorem 1.Suppose thaty=fx is the nonzero solution of (4), M=−βK/2r2+K−a/2, N=K2β2−2K2βr2+K2r22−2Kaβr2+2Kar22+a2r22, m1=−r1Kβ/r2+bK, m0=KMr1+Kar1−M2r1−Mar1, V=4r22m0m1−4Km12br2/4r22, and W=4r22m02−m12N/4r22. In Theorem 1, we show the feasibility and the uniqueness conditions of the nontrivial positive equilibrium point x∗,y∗.Theorem 1.
If the conditions that(5)r1β−r2b>0,r2−β>0,β<bKK+a,−r11+f′x2x+2fx+a−K−Kbf′x<0,
hold, then there is a unique positive equilibriumE3x∗,y∗ of System (1), where(6)x∗=−V+V2−4m12W2m12,y∗=M−x∗+12r2N+4Kbr2x∗.Proof.
Since we are looking for nonzero solution of (3) and (4), we can simplify these equations into(7)r1−r1x+yK−byx+y+a=0,(8)r2−r2x+yK+bxx+y+a−β=0.
We solve (8) by acting it as a quadratic equation of y. The formula of positive y∗ is given in (6). Next, substituting this solution to (7) yields(9)m1x∗+m0=12r2m1N+4Kbr2x∗.Here we obtain nonlinear form of equation. By squaring both sides to neglect the square root, we obtain a quadratic equation. Since we only consider a unique positive equilibrium, the positive roots are given in (6).
Talking about the stability behavior of the equilibrium point, the following theorem shows the type of stability for positive equilibrium point if it exists.Theorem 2.
For System (1), if the delay does not appear in the system, then we have that equilibriumE3 is locally asymptotically stable if it exists.Proof.
The characteristic equation of the linearized system atE3x∗,y∗ is(10)λ2+r1x∗+r2y∗Kλ+bx∗y∗r1−r2x+y+aK+abx+y+a3=0.
By the conditions in Theorem1, we can derive that(11)bβ>K+aK>1.
This implies that(12)r1−r2>0.
Therefore, both coefficients of the characteristic equation are positive. In other words, the roots are two negative real numbers or a pair of complex number with a negative real part. Here, we can conclude that equilibriumE3 is locally asymptotically stable if it exists.


### 3.2. The Stability Region

By the fact that the nontrivial positive equilibrium point (see (6)) is locally asymptotically stable for the absence of delay case, the solutions near the equilibrium points will converge to the equilibrium. However, the existence of a delay in System (1) could bring different behavior. In this section, we will identify the delay effect on the stability behavior of the solution near the nontrivial positive equilibrium point.If the property holds for all values of the delays, it is named delay-independent. Meanwhile, if the stability is preserved for some values of delays but become unstable for other values, it is named delay-dependent.Let us assume that the conditions in (5) hold such that the positive equilibrium point exists. Since we only consider the solution around the positive equilibrium point, we linearize System (1) nearE3x∗,y∗. Using transformations Z1=x−x∗ and Z2=y−y∗, the linearized system can be written as(13)Z˙1Z˙2=A∗B∗C∗D∗Z1Z2+00E∗0Z1τZ2τ.We have(14)A∗=−r1x∗K+bx∗y∗x∗+y∗+a2,B∗=−r1x∗K−bx∗x∗+y∗+a+bx∗y∗x∗+y∗+a2,C∗=by∗x∗+y∗+a−bx∗y∗x∗+y∗+a2,D∗=−r2y∗K−bx∗y∗x∗+y∗+a2,E∗=−r2y∗K,Z1τ=xτ−x∗,Z2τ=yτ−y∗.Suppose that the solution is a nontrivial exponential solutionZ1Z2T=eλtu1u2T, where u1 and u2 are constants. To obtain a nontrivial solution, the value of λ should satisfy(15)detλI−A∗B∗C∗D∗−00E∗0e−λτ=0.I is the identity matrix. According to the definition defined in [25], (15) is later called the characteristic equation associated with (13). The characteristic equation of a delayed linear system is in transcendent form and for System (1) the characteristic equation is given by(16)λ2−T¯λ+D¯−E¯e−λτ=0,where(17)D¯=A∗D∗−B∗C∗,T¯=A∗+D∗,danE¯=E∗B∗.Moreover, since we deal with positive equilibrium point and all parameters are nonnegative, we obtain thatT<0 and E>0.As applied for ODEs, the zero solution of (16) is asymptotically stable if and only if it has no zeros in the right half plane λ|Reλ≥0 [26]. Beforehand, we study the stability sets in its parameter spaces. To investigate the stability region of this transcendent equation, we consider the Pontryagin Theorem. Next, we are going to provide the sufficient and necessary condition such that the zero solution of System (13) is asymptotically stable.If a polynomialhz,λ does not possess a principle term, then it definitely has many number of zeros with positive real part. In other words, the zero solution of this polynomial is not asymptotically stable. However, if it has a principal term, then the Pontryagin Theorem could be applied.Theorem 3 (the Pontryagin Theorem [18]).
LetHz=hz,ez, where hz,w is a polynomial with a principal term. Function Hiy is now separated into real and imaginary parts; that is, we set Hiy=Fy+iGy. If all the zeros of function Hz lie in the open left half plane, then the zeros of functions Fy and Gy are real and interlacing and(18)Δy=G′yFy−GyF′y>0.It is for all realy. Moreover, in order that all the zeros of function Hz lie in the open left half plane, it is sufficient that one of the following conditions is satisfied:(1)
All of the zeros of functionsFy and Gy are real and interlacing, and inequality (18) is satisfied for at least one value of r.(2)
All the zeros of functionFy are real and for each of these zeros y=r the condition in (18) is satisfied; that is, F′rGr<0.(3)
All the zeros of functionGy are real and for each of these zeros inequality (18) is satisfied; that is, G′rFr>0.Multiplying both sides of characteristic (16) by eλτ and letting z=λτ will transform that equation into(19)z2ez−Tzez+Dez−E=0,where T=τT¯,D=τ2D¯, and E=τ2E¯. We obtain that the left side of this equation has a principal term, z2ez, so the Pontryagin Theorem could be applied.Let us define(20)Hz=z2ez−Tzez+Dez−E.Here we can restate(21)Hiy=Fy+iGy,where(22)Fy=−y2cosy+Tysiny+Dcosy−E,(23)Gy=−y2siny−Tycosy+Dsiny.In the following theorem, we show that, forT≤D≤E, the zero solution of (16) is not asymptotically stable. Therefore, for the next, we restrict our discussion for D<T or D>E.Theorem 4.
GivenT<0 and E>0, if T≤D≤E, then the zero solution of (16) is not asymptotically stable.Proof.
It is clear thaty=0 is a root for function Gy. Here we obtain that(24)Δ0=G′0F0=−T+DD−E≤0.
Thus, by Theorem3, we obtain that the zero solution of (16) is not asymptotically stable for T≤D≤E.
It is obvious that all roots of functionF and/or G are real. In Theorem 5, one shows sufficient and necessary conditions such that all the roots of a function are real numbers. To use the theorem, let us consider a polynomial in the form(25)fz,u,v=∑m=0p∑n=0qbmnzmϕnu,v,where ϕmnu,v is a homogeneous polynomial of degree n in u and v. Denote brszrϕrsu,v as the principal term of fz,u,v; then(26)fz,u,v=zrϕ∗su,v+∑m=0r−1∑n=0sbmnzmϕmnu,v,where ϕ∗su,v=∑n=0sbrnϕrnu,v. Also we let(27)Φ∗sz=ϕ∗scosz,sinz.Theorem 5.
Letfz,u,v be a polynomial with principal term zrϕrsu,v. If ε is such that Φ∗sε+iy≠0 for all real y, then, in the strip −2πk+ε≤x≤2πk+ε, z=x+iy, function Fz=fz,cosz,sinz will have for all sufficiently large values of k exactly 4sk+r zeros. Thus, in order for function Fz0 to have only real roots, it is necessary and sufficient that, in the interval −2πk+ε≤x≤2πk+ε, it has exactly 4sk+r real roots for all sufficiently large k.Letgy,u,v=−y2v−Tyu+Dv; then Gy=gy,cosy,siny. The principal term of this function is −y2v, so we obtain that r and s in Theorem 5 are 1 and 2, respectively. Moreover,(28)Φ∗1ε+iy=−sinε+iy,=−sinεcoshy−icosεsinhy.Takingε=π/2 yields(29)Φ∗1π2+iy=−coshy<0.Hence, Theorem5 implies that the zeros of Gy are all real if and only if Gy has exactly 4k+2 roots in the interval −2πk+π/2≤y≤2πk+π/2.Theorem 6.
GivenT<0 and E>0, for D>0, the zero solution of System (13) is to be asymptotically stable if and only if(1)
D>E,(2)
Fr∗>0.Herer∗ is the root of function G and r∗∈π,2π.Proof.
⇒
As mentioned before thaty=0 is a root of function Gy, the zero solution of this function is asymptotically stable; therefore we obtain that Δ0>0. Since D>0, this implies D−T>0. By the condition that Δ0>0, we have D>E. Next, we will show that the second condition holds. Suppose that rj is a nonzero solution of Gy=0. From equation (22) we can write(30)D=rj2sinrj+Trjcosrjsinrj.
By substituting (30) into the first derivative function G′y, we can rewrite this function as(31)G′rj=−2rjsinrj−Tcosrj+Trjsinrj.
The value ofG′rj can be seen in Figure 1.
Since the minimum value ofTy/siny is greater than the maximum value of Tcosy in interval π,2π, we can conclude that(32)Tysiny>Tcosy,forT<0.
Moreover, the Pontryagin Theorem yields(33)FrjG′rj>0,for allj∈ℤ.
Forrj=r∗∈π,2π, from (32), we obtain that G′r∗>0. From (33) and G′r∗>0, we can conclude that(34)Fr∗>0S.
⇐
To show that the zero solution is asymptotically stable, we are going to prove that the Pontryagin conditions hold by showing that all roots ofGy are real numbers and Δy0>0 for all y0 real roots of function Gy. It is clear that Gy is an odd function; meanwhile, G′y and Fy are even functions. So, we restrict our discussion to the nonnegative domain and use the property of odd and even functions to extend the domain.(1)
IfD>E, then all roots of function Gy are real roots.Fory=jπ, we get Gjπ=−1j+1Tjπ. This implies that Gjπ<0 for all j positive odd numbers and Gjπ>0 for all j positive even numbers. According to Intermediate Value Theorem, there exists at least one real root. Moreover, there is an odd number of root. For y=y0 nonzero root of Gy, it satisfies(35)wy0=ξy0,wherewy=D/y−y and ξt=Tcoty. For D>E>0, we obtain that w′y<0 and, for T<0, we obtain that ξ′y>0. Because ξy is a periodic function with period π and both functions are monotone, we can say that there is only one real root in each interval jπ,j+1π. This implies in interval −2kπ,2kπ that there exist 4k+1 roots for all k∈ℤ, including y=0. Taking 0<ε<π/2 and for k sufficiently large, function Gy has one root and no roots in intervals 2k,2k+ε and −2kπ,−2kπ+ε, respectively. Therefore, the number of roots in −2kπ+ε,2kπ+ε is 4k+2. In other words, function Gy has all real roots.(2)
IfFr∗>0, then Δrj>0.Fory=0, it is clear that Δ0>0. Suppose that y=rj is the root of function Gy in the interval jπ,j+1π for all j positive integers. By substituting (30) into function Fy, we can rewrite this function as(36)Frj=Trj1sinrj−E.
Now, let us consider that forj=2k where k∈ℤ+∪0, it is obvious that Fr2k<0. Since G′r2k<0, we can conclude that it is satisfying the Pontryagin Condition:(37)G′r2kFr2k>0.
Next, forj=2k+1 where k∈ℤ+∪0, using the similar argument as given in (32) yields G′r2k+1>0. Based on the monotonicity of r2k+1 converging to 2k+1π, limy⟶2k+1π+Fy=∞, and we can conclude that Fr2k+1<F2k+3 for all k∈ℤ∪0. The illustration of Frj can be seen in Figure 2. Since Fr∗>0, we obtain that F2k+1>0. According to this, the Pontryagin Condition also holds for all k∈ℤ∪0.(38)G′r2k+1Fr2k+1>0.
We have shown that the Pontryagin Condition is satisfied for nonnegative domain. Next, using the properties thatFy and G′y are both even functions, we obtain that(39)F−rjG′−rj=FrjG′rj>0.
The proof is complete.
Next, whether the stability behavior of the equilibrium is changed by the existence of the delay or not will be discussed.Figure 1
The graph ofG′y and its value at the roots of Gy, rj in positive domain.Figure 2
The graph ofFy and its value at the roots of Gy, rj in positive domain.


## 3.1. The Equilibrium Points

Letℛ=ℝ+∪0 be a nonnegative real number set. The domain of System (1) is denoted by(2)Ω=x,y∈ℛ2|0<x+y≤K,where all the parameters r1,r2,K,b,a, and β are nonnegative. The equilibrium solutions of System (1) must satisfy(3)r1x−r1xx+yK−bxyx+y+a=0,(4)r2y−r2yxτ+yK+bxyx+y+a−βy=0.In the absence of delay, System (1) has four possible types of equilibrium points, that is,E00,0 as a trivial equilibrium, E1K,0 and E20,y¯ with y¯=β−r2/r2 as semitrivial equilibria, and E3x∗,y∗ as the nontrivial equilibrium point. The equilibrium points E0,E1, and E2 indicate the extinction of one or both populations.There is more than one positive nontrivial equilibrium point due to the parameter space. However, in this paper, we restrict our discussion to the parameter space that the system has a unique positive equilibrium point. The sufficient conditions of the uniqueness were given in [10]. Here, we provide the expression of the unique positive equilibrium point in Theorem 1.Suppose thaty=fx is the nonzero solution of (4), M=−βK/2r2+K−a/2, N=K2β2−2K2βr2+K2r22−2Kaβr2+2Kar22+a2r22, m1=−r1Kβ/r2+bK, m0=KMr1+Kar1−M2r1−Mar1, V=4r22m0m1−4Km12br2/4r22, and W=4r22m02−m12N/4r22. In Theorem 1, we show the feasibility and the uniqueness conditions of the nontrivial positive equilibrium point x∗,y∗.Theorem 1.
If the conditions that(5)r1β−r2b>0,r2−β>0,β<bKK+a,−r11+f′x2x+2fx+a−K−Kbf′x<0,
hold, then there is a unique positive equilibriumE3x∗,y∗ of System (1), where(6)x∗=−V+V2−4m12W2m12,y∗=M−x∗+12r2N+4Kbr2x∗.Proof.
Since we are looking for nonzero solution of (3) and (4), we can simplify these equations into(7)r1−r1x+yK−byx+y+a=0,(8)r2−r2x+yK+bxx+y+a−β=0.
We solve (8) by acting it as a quadratic equation of y. The formula of positive y∗ is given in (6). Next, substituting this solution to (7) yields(9)m1x∗+m0=12r2m1N+4Kbr2x∗.Here we obtain nonlinear form of equation. By squaring both sides to neglect the square root, we obtain a quadratic equation. Since we only consider a unique positive equilibrium, the positive roots are given in (6).
Talking about the stability behavior of the equilibrium point, the following theorem shows the type of stability for positive equilibrium point if it exists.Theorem 2.
For System (1), if the delay does not appear in the system, then we have that equilibriumE3 is locally asymptotically stable if it exists.Proof.
The characteristic equation of the linearized system atE3x∗,y∗ is(10)λ2+r1x∗+r2y∗Kλ+bx∗y∗r1−r2x+y+aK+abx+y+a3=0.
By the conditions in Theorem1, we can derive that(11)bβ>K+aK>1.
This implies that(12)r1−r2>0.
Therefore, both coefficients of the characteristic equation are positive. In other words, the roots are two negative real numbers or a pair of complex number with a negative real part. Here, we can conclude that equilibriumE3 is locally asymptotically stable if it exists.


## 3.2. The Stability Region

By the fact that the nontrivial positive equilibrium point (see (6)) is locally asymptotically stable for the absence of delay case, the solutions near the equilibrium points will converge to the equilibrium. However, the existence of a delay in System (1) could bring different behavior. In this section, we will identify the delay effect on the stability behavior of the solution near the nontrivial positive equilibrium point.If the property holds for all values of the delays, it is named delay-independent. Meanwhile, if the stability is preserved for some values of delays but become unstable for other values, it is named delay-dependent.Let us assume that the conditions in (5) hold such that the positive equilibrium point exists. Since we only consider the solution around the positive equilibrium point, we linearize System (1) nearE3x∗,y∗. Using transformations Z1=x−x∗ and Z2=y−y∗, the linearized system can be written as(13)Z˙1Z˙2=A∗B∗C∗D∗Z1Z2+00E∗0Z1τZ2τ.We have(14)A∗=−r1x∗K+bx∗y∗x∗+y∗+a2,B∗=−r1x∗K−bx∗x∗+y∗+a+bx∗y∗x∗+y∗+a2,C∗=by∗x∗+y∗+a−bx∗y∗x∗+y∗+a2,D∗=−r2y∗K−bx∗y∗x∗+y∗+a2,E∗=−r2y∗K,Z1τ=xτ−x∗,Z2τ=yτ−y∗.Suppose that the solution is a nontrivial exponential solutionZ1Z2T=eλtu1u2T, where u1 and u2 are constants. To obtain a nontrivial solution, the value of λ should satisfy(15)detλI−A∗B∗C∗D∗−00E∗0e−λτ=0.I is the identity matrix. According to the definition defined in [25], (15) is later called the characteristic equation associated with (13). The characteristic equation of a delayed linear system is in transcendent form and for System (1) the characteristic equation is given by(16)λ2−T¯λ+D¯−E¯e−λτ=0,where(17)D¯=A∗D∗−B∗C∗,T¯=A∗+D∗,danE¯=E∗B∗.Moreover, since we deal with positive equilibrium point and all parameters are nonnegative, we obtain thatT<0 and E>0.As applied for ODEs, the zero solution of (16) is asymptotically stable if and only if it has no zeros in the right half plane λ|Reλ≥0 [26]. Beforehand, we study the stability sets in its parameter spaces. To investigate the stability region of this transcendent equation, we consider the Pontryagin Theorem. Next, we are going to provide the sufficient and necessary condition such that the zero solution of System (13) is asymptotically stable.If a polynomialhz,λ does not possess a principle term, then it definitely has many number of zeros with positive real part. In other words, the zero solution of this polynomial is not asymptotically stable. However, if it has a principal term, then the Pontryagin Theorem could be applied.Theorem 3 (the Pontryagin Theorem [18]).
LetHz=hz,ez, where hz,w is a polynomial with a principal term. Function Hiy is now separated into real and imaginary parts; that is, we set Hiy=Fy+iGy. If all the zeros of function Hz lie in the open left half plane, then the zeros of functions Fy and Gy are real and interlacing and(18)Δy=G′yFy−GyF′y>0.It is for all realy. Moreover, in order that all the zeros of function Hz lie in the open left half plane, it is sufficient that one of the following conditions is satisfied:(1)
All of the zeros of functionsFy and Gy are real and interlacing, and inequality (18) is satisfied for at least one value of r.(2)
All the zeros of functionFy are real and for each of these zeros y=r the condition in (18) is satisfied; that is, F′rGr<0.(3)
All the zeros of functionGy are real and for each of these zeros inequality (18) is satisfied; that is, G′rFr>0.Multiplying both sides of characteristic (16) by eλτ and letting z=λτ will transform that equation into(19)z2ez−Tzez+Dez−E=0,where T=τT¯,D=τ2D¯, and E=τ2E¯. We obtain that the left side of this equation has a principal term, z2ez, so the Pontryagin Theorem could be applied.Let us define(20)Hz=z2ez−Tzez+Dez−E.Here we can restate(21)Hiy=Fy+iGy,where(22)Fy=−y2cosy+Tysiny+Dcosy−E,(23)Gy=−y2siny−Tycosy+Dsiny.In the following theorem, we show that, forT≤D≤E, the zero solution of (16) is not asymptotically stable. Therefore, for the next, we restrict our discussion for D<T or D>E.Theorem 4.
GivenT<0 and E>0, if T≤D≤E, then the zero solution of (16) is not asymptotically stable.Proof.
It is clear thaty=0 is a root for function Gy. Here we obtain that(24)Δ0=G′0F0=−T+DD−E≤0.
Thus, by Theorem3, we obtain that the zero solution of (16) is not asymptotically stable for T≤D≤E.
It is obvious that all roots of functionF and/or G are real. In Theorem 5, one shows sufficient and necessary conditions such that all the roots of a function are real numbers. To use the theorem, let us consider a polynomial in the form(25)fz,u,v=∑m=0p∑n=0qbmnzmϕnu,v,where ϕmnu,v is a homogeneous polynomial of degree n in u and v. Denote brszrϕrsu,v as the principal term of fz,u,v; then(26)fz,u,v=zrϕ∗su,v+∑m=0r−1∑n=0sbmnzmϕmnu,v,where ϕ∗su,v=∑n=0sbrnϕrnu,v. Also we let(27)Φ∗sz=ϕ∗scosz,sinz.Theorem 5.
Letfz,u,v be a polynomial with principal term zrϕrsu,v. If ε is such that Φ∗sε+iy≠0 for all real y, then, in the strip −2πk+ε≤x≤2πk+ε, z=x+iy, function Fz=fz,cosz,sinz will have for all sufficiently large values of k exactly 4sk+r zeros. Thus, in order for function Fz0 to have only real roots, it is necessary and sufficient that, in the interval −2πk+ε≤x≤2πk+ε, it has exactly 4sk+r real roots for all sufficiently large k.Letgy,u,v=−y2v−Tyu+Dv; then Gy=gy,cosy,siny. The principal term of this function is −y2v, so we obtain that r and s in Theorem 5 are 1 and 2, respectively. Moreover,(28)Φ∗1ε+iy=−sinε+iy,=−sinεcoshy−icosεsinhy.Takingε=π/2 yields(29)Φ∗1π2+iy=−coshy<0.Hence, Theorem5 implies that the zeros of Gy are all real if and only if Gy has exactly 4k+2 roots in the interval −2πk+π/2≤y≤2πk+π/2.Theorem 6.
GivenT<0 and E>0, for D>0, the zero solution of System (13) is to be asymptotically stable if and only if(1)
D>E,(2)
Fr∗>0.Herer∗ is the root of function G and r∗∈π,2π.Proof.
⇒
As mentioned before thaty=0 is a root of function Gy, the zero solution of this function is asymptotically stable; therefore we obtain that Δ0>0. Since D>0, this implies D−T>0. By the condition that Δ0>0, we have D>E. Next, we will show that the second condition holds. Suppose that rj is a nonzero solution of Gy=0. From equation (22) we can write(30)D=rj2sinrj+Trjcosrjsinrj.
By substituting (30) into the first derivative function G′y, we can rewrite this function as(31)G′rj=−2rjsinrj−Tcosrj+Trjsinrj.
The value ofG′rj can be seen in Figure 1.
Since the minimum value ofTy/siny is greater than the maximum value of Tcosy in interval π,2π, we can conclude that(32)Tysiny>Tcosy,forT<0.
Moreover, the Pontryagin Theorem yields(33)FrjG′rj>0,for allj∈ℤ.
Forrj=r∗∈π,2π, from (32), we obtain that G′r∗>0. From (33) and G′r∗>0, we can conclude that(34)Fr∗>0S.
⇐
To show that the zero solution is asymptotically stable, we are going to prove that the Pontryagin conditions hold by showing that all roots ofGy are real numbers and Δy0>0 for all y0 real roots of function Gy. It is clear that Gy is an odd function; meanwhile, G′y and Fy are even functions. So, we restrict our discussion to the nonnegative domain and use the property of odd and even functions to extend the domain.(1)
IfD>E, then all roots of function Gy are real roots.Fory=jπ, we get Gjπ=−1j+1Tjπ. This implies that Gjπ<0 for all j positive odd numbers and Gjπ>0 for all j positive even numbers. According to Intermediate Value Theorem, there exists at least one real root. Moreover, there is an odd number of root. For y=y0 nonzero root of Gy, it satisfies(35)wy0=ξy0,wherewy=D/y−y and ξt=Tcoty. For D>E>0, we obtain that w′y<0 and, for T<0, we obtain that ξ′y>0. Because ξy is a periodic function with period π and both functions are monotone, we can say that there is only one real root in each interval jπ,j+1π. This implies in interval −2kπ,2kπ that there exist 4k+1 roots for all k∈ℤ, including y=0. Taking 0<ε<π/2 and for k sufficiently large, function Gy has one root and no roots in intervals 2k,2k+ε and −2kπ,−2kπ+ε, respectively. Therefore, the number of roots in −2kπ+ε,2kπ+ε is 4k+2. In other words, function Gy has all real roots.(2)
IfFr∗>0, then Δrj>0.Fory=0, it is clear that Δ0>0. Suppose that y=rj is the root of function Gy in the interval jπ,j+1π for all j positive integers. By substituting (30) into function Fy, we can rewrite this function as(36)Frj=Trj1sinrj−E.
Now, let us consider that forj=2k where k∈ℤ+∪0, it is obvious that Fr2k<0. Since G′r2k<0, we can conclude that it is satisfying the Pontryagin Condition:(37)G′r2kFr2k>0.
Next, forj=2k+1 where k∈ℤ+∪0, using the similar argument as given in (32) yields G′r2k+1>0. Based on the monotonicity of r2k+1 converging to 2k+1π, limy⟶2k+1π+Fy=∞, and we can conclude that Fr2k+1<F2k+3 for all k∈ℤ∪0. The illustration of Frj can be seen in Figure 2. Since Fr∗>0, we obtain that F2k+1>0. According to this, the Pontryagin Condition also holds for all k∈ℤ∪0.(38)G′r2k+1Fr2k+1>0.
We have shown that the Pontryagin Condition is satisfied for nonnegative domain. Next, using the properties thatFy and G′y are both even functions, we obtain that(39)F−rjG′−rj=FrjG′rj>0.
The proof is complete.
Next, whether the stability behavior of the equilibrium is changed by the existence of the delay or not will be discussed.Figure 1
The graph ofG′y and its value at the roots of Gy, rj in positive domain.Figure 2
The graph ofFy and its value at the roots of Gy, rj in positive domain.


## 4. Bifurcation Analysis for Delay Parameter

A bifurcation analysis for the expansion of the model with parameterized logistic growth has been carried out in [24]. At the beginning, we have discussed that the positive equilibrium point is asymptotically stable for nondelay case. Due to the continuity of the delay, the following theorem shows that there exists τ^>0 such that interior equilibrium point will remain asymptotically stable for every τ∈[0,τ^). However, one should check whether the stability behavior is preserved or changed by the existence of the delay. Therefore, we will also investigate whether there is τ0 such that E3 changes stability to τ>τ0. If there is, then this equilibrium point is called conditionally stable. In this chapter, we will look for when the Hopf bifurcation occurs with respect to the delay parameter and then analyze the characteristic if there is a small change at the bifurcation point.Theorem 7.
For delay case, there existsτ^ such that the E3 equilibrium point of System (1) is asymptotically stable for every τ∈[0,τ^).Proof.
By the continuous dependence of roots of (16) and Theorem 2, there exists τ^ such that Reλτ<0 for τ∈[0,τ^).
It is well known from the stability theory of linear time-invariant time-delayed systems (LTI-TDS) that the pure imaginary characteristic roots are the only possible transition points from stable to unstable behavior, and vice versa [27]. First, we consider looking for the existence of these bifurcation points. The bifurcation point is reached when v=0 and the first partial derivative with respect to w of the characteristic function is not equal to zero. Suppose that λ=iww>0 is the root of (16).(40)−e−iwτE¯+D¯−iT¯w−w2=0.
Separating the real and imaginary parts yields(41)−w2+D¯=E¯coswτ,(42)T¯w=E¯sinwτ.
Squaring the two sides and adding them together will give a quartic polynomial given in(43)−w2+D¯2+T¯2w2=E¯2.
Expanding this formula yields(44)w4+T¯2−2D¯w2+D¯2−E¯2=0.
By considering (43) as a quadratic equation for a new variable w2, the solution to this equation can be expressed in terms of(45)w1,22=−12T¯2+D¯±12T¯4−4D¯T¯2+4E¯2.
In the following lemma, we present sufficient and necessary conditions on which there is only one positive root forw2. Consequently the characteristic in (16) has exactly two imaginary roots, ±iw0.Lemma 1.
Suppose thatD¯2−E¯2<0, and T¯4−4DT¯2+4E¯2>0. The characteristic in (16) has only two pure imaginary roots, ±iw0, where(46)w0=−12T¯2+D¯+12T¯4−4D¯T¯2+4E¯21/2,for(47)τn=1w0cos−11E¯D¯−w02+2nπw0,n∈ℕ.Proof.
Use the properties of real numbers. Forp>0, s2+p>s≥±s. According to these properties, by taking s=T¯2−2D¯ and p=4E¯2−D¯2, we can derive that(48)w12=−12T¯2+D¯+12T¯4−4D¯T¯2+4E¯2>0,(49)w22=−12T¯2+D¯−12T¯4−4D¯T¯2+4E¯2<0.
Suppose thatw0 is the characteristic root of (16). Since we only consider for w0 real number, the only positive real root is given in (47), w0=w12. Moreover, by substituting w=w0 into (40) one obtains that the value of τ is given by (1). This completes the proof.
From these results, we can conclude that when the delay time reachesτ=τn, we get a pure imaginary characteristic root. Consequently, the solution of the system is a periodic solution. Subsequently, the smallest bifurcation value is achieved for τ=τ0 where(50)τ0=1w0cos−11E¯D¯−w02.
Sinceτ0 is the smallest bifurcation point, according to Theorem 7, we can conclude that the zero solution of the characteristic equation is asymptotically stable for τ∈[0,τ0). In the following theorem, we will show that there is a change in the stability properties before and after the bifurcation point by referring to the first derivative of the real part of the characteristic root at the equilibrium point. Suppose that the eigenvalues are expressed as a function of delay, λτ=vτ+iwτ.Theorem 8.
Letλτ=vτ+iwτ be the root of (16) satisfying vτ0=0,wτ0=w0 defined by (45); then the following transversality condition holds:(51)dvλdτ|v=0,w=w0,τ=τ0>0.Proof.
Differentiating the characteristic in (16) with respect to τ and using the implicit function theorem yield(52)2λdλdτ−T¯dλdτ−E¯e−λτ−τdλdτ−λ=0,ordλdτ=−E¯λe−λτ2λ−T¯+E¯τe−λτ.
Next, it will be shown that the derivative at the bifurcation value is nonzero. According to the characteristic equation, we use(53)e−λτ=λ2−T¯λ+D¯E¯.
We obtain that the first derivative of the real part of eigen value with respect toτ at bifurcation point τ0,w0 is(54)dvdτ|0,ω0,τ0=w02T¯2+2w04−2w02D¯Z0,where Z0=τw02+τD¯−T¯2+2w0−τT¯w02. From (47), since w02=w12, one obtains that w02T¯2+2w04−2w02D¯=w02T¯2+2w02+2E¯2=w02T¯4−4DT¯2+4E¯2>0. Therefore, both the numerator and the denominator of the first derivative are positive. In other words,(55)dvλdτ|v=0,w=w0,τ=τ0>0.
Therefore, the transversality condition holds.
Another possibility for the roots of (42) is that both roots are positive. The sufficient and necessary condition on which there are two positive roots for w2 is given in Lemma 2.Lemma 2.
Suppose thatD¯2−E¯2>0, T¯2−2D¯<0, and T¯4−4D¯T¯2+4E¯2>0. The characteristic in (16) has two pure positive imaginary roots, iw0±, where(56)w0±=−12T¯2+D¯±12T¯4−4D¯T¯2+4E¯21/2.This is for(57)τn=1w0±cos−11E¯D¯−w0±2+2nπw0±,n∈ℕ.Proof.
The third condition implied that the characteristic roots are real number. Moreover, note that(58)w0++w0−=−T¯2+2D¯>0,w0+×w0−=D¯2−E¯2>0.
The multiplication of two real numbers is positive; therefore both numbers are positive or negative. Since the summation is also positive, we can conclude that both numbers are positive.
In the following theorem, we will show that if the conditions are satisfied, then there is a change in stability for the given delay value.Theorem 9.
Letλτ=vτ+iwτ be the root of (16) satisfying vτ0=0,wτ0± defined in (55); the following transversality condition holds:(59)dvλdτv=0,w=w0+,τ=τ0+>0,dvλdτ|v=0,w=w0−,τ=τ0−<0.Proof.
Similar to Theorem8, we obtain that(60)dvdτ|v=0,w=w0−,τ=τ0−=w0+T¯2+2w0+2−2D¯Z0+=w0+2T¯4−4D¯T¯2+4E¯2Z0+>0,dvdτ|v=0,w=w0−,τ=τ0−=w0−T¯2+2w0−2−2D¯Z0−=−w0−2T¯4−4D¯T¯2+4E¯2Z0−<0.Z0±=τw0±2+τD¯−T¯2+2w0±−τT¯w0±2.
From the result in Theorems8 and 9, we can conclude that the positive equilibrium of System (1) is conditionally stable if it is satisfying sufficient conditions given in Lemma 1 or Lemma 2. Otherwise, it is absolutely stable equilibrium point. Next, the results of the numerical analysis show the effect of the presence of delay on the system solution for the initial value around the positive equilibrium point.


## 5. Numerical Simulation

In this section, we would like to verify the results from the previous sections by showing some numerical simulations. Here, we consider the parameter values given in [10].(61)r1=40,K=100,r2=2,a=0.05,b=0.02,β=0.003.The equilibrium points for this set of parameters are(62)E00,0,E1100,0,E20,99.85,E310.5295,89.4258.Here we consider the positive equilibrium point E3 which later in the figures will be depicted as star∗. Using these parameter values, we obtain that(63)D¯=7.604830378,T¯=−6.000317366,E¯=7.533271847.Since this satisfies the first condition of Theorem6, we can see that the unstable condition holds whenever the second condition is not satisfied. However, the coefficient of function Gy is also dependent on the delay and the formula of the first positive root of Gy has not been provided. Therefore, we analyze the stability properties of this model using the bifurcation analysis of the delay term given in Lemmas 1 and 2.Here, we takex0=10 and y0=88 as the constant initial population of prey and predator which is close to the positive equilibrium point. For τ=0, as given in the phase portrait in Figure 3(a) and the solution graph for the predator in Figure 3(b), we can see that the solution converges to E3. Let us perform comparison to the models with delay, for τ=2 and τ=40.Figure 3
Nondelay case,τ=0. (a) The phase portrait; (b) the number of predators converging to the equilibrium point.
(a)(b)Forτ=2, we have D¯2−E¯2=1.08326036, T¯2−2D¯=20.79414773, and T¯4−4D¯T¯2+3E¯=428.063539. Since the conditions in both Lemmas 1 and 2 are not satisfied, we can conclude that this equilibrium point is absolutely stable. However, the populations reach oscillation near the equilibrium point and as the time goes they converge to their equilibrium value as seen in Figure 4. Though different property is shown for τ=40, it also converges to equilibrium point. This behavior can be seen in Figure 5.Figure 4
Delay case,τ=2. (a) The phase portrait; (b) the number of predators with short delay.
(a)(b)Figure 5
Delay case,τ=40. (a) The phase portrait; (b) the number of predators with long delay.
(a)(b)


## 6. Concluding Remarks

In this paper, we have provided the formula of the positive equilibrium point and the condition of its existence. For nondelay case, the stability of this point is asymptotically stable and, due to the continuity property, this behavior remains stable for some interval delay time. However, the existence of delay could change the stability behavior of the equilibrium point. As can be seen in this paper that we have not provided the stability conditions forD<0, according to Theorem 4 one of the conditions for this case is that D should be less than T. Nevertheless, the analytical explanation according to Pontryagin Theorem for this case has not been discussed before.



---

*Source: 3711158-2022-04-22.xml*