# The Gestation Delay: A Factor Causing Complex Dynamics in Gause-Type Competition Models

Authors: Zizhen Zhang; Ranjit Kumar Upadhyay; Rashmi Agrawal; Jyotiska Datta
Journal: Complexity
 (2018)

---

## Abstract

In this paper, we consider a Gause-type model system consisting of two prey and one predator. Gestation period is considered as the time delay for the conversion of both the prey and predator. Bobcats and their primary prey rabbits and squirrels, found in North America and southern Canada, are taken as an example of an ecological system. It has been observed that there are stability switches and the system becomes unstable due to the effect of time delay. Positive invariance, boundedness, and local stability analysis are studied for the model system. Conditions under which both delayed and nondelayed model systems remain globally stable are found. Criteria which guarantee the persistence of the delayed model system are derived. Conditions for the existence of Hopf bifurcation at the nonzero equilibrium point of the delayed model system are also obtained. Formulae for the direction, stability, and period of the bifurcating solution are conducted using the normal form theory and center manifold theorem. Numerical simulations have been shown to analyze the effect of each of the parameters considered in the formation of the model system on the dynamic behavior of the system. The findings are interesting from the application point of view.


---

## Body

## 1. Introduction

Persistence, structure, and dynamics of biological species have received great attention from ecologists as well as mathematical biologists. The interspecific interference such as competition and predation dominates the dynamics of the system. Competitive exclusion principle by Gause states that when two species compete for the same resources within the environment, one of them will eventually outcompete and displace the other. The displaced population may become locally extinct by either migration or death, or it may adapt to a sufficiently distinct niche within the environment [1]. Upadhyay and Iyengar [2] studied the extinction and coexistence of two competing prey species of a Gause-type ecological model in which the predator species is influenced by the damage effect caused by crowding from the members of its own population.In ecosystem modeling, the type of predation is very important, and consecutively, the dynamics become richer with more realistic assumptions [3]. The functional response (rate of prey consumption by an average predator) is one of the important features that govern the dynamics of the prey-predator system. Holling type II saturated functional response has been widely used to investigate the behavior of the prey-predator system. In real situation, conversion of prey biomass into predator biomass does not take place instantaneously. There is always a time lag of gestation. To make the system considered in [2] more realistic, we introduced the gestation period as a time delay and studied the effect of delay on the dynamics of the model system. Kar and Batabyal [4] have studied a two-prey one-predator system in the presence of time delay and derived the conditions for persistence and global stability of the model system. Recently, Deka et al. [5] considered a general Gause-type two-prey one-predator model system and studied the effect of predation on two competing prey species. The authors observed Hopf bifurcation by taking the death rate of the predator population as a bifurcation parameter. Jana et al. [6] also studied a two-prey one-predator model system in which the prey species is divided into weaker and stronger classes due to predator’s catching efficiency. The evolutionary effects of selective disturbances on an evolving trait (e.g., body size and maturation age) of the predator individuals in a one-predator two-prey community are investigated by Meng et al. [7]. Using methods of adaptive dynamics, they constructed an invasion fitness function and obtained the conditions for evolutionary branching and stability under selective disturbances in both monomorphic and dimorphic populations. Lu et al. [8] developed a rigorous analysis to investigate how the dual phenomena of prey infection and hunting strategies of the predators influence the disease dynamics in a predator-prey model with infected prey and stage-structured predator. A three-species model system is also studied by Gakkhar and Gupta [9] in which two competing prey species grow logistically but the third species acts as a predator as well as host, and they observed that the possibility of survival of all three species is due to commensalism of a prey and predator population. The dynamic behavior of a two-prey one-predator system with impulsive effects concerning periodic biological and chemical control is investigated by Pei et al. [10], and it was observed that there exist asymptotically stable pest-eradication periodic solutions if the impulsive period is less than some threshold. Wang et al. [11] studied a two-prey one-predator system with Watt-type functional response and impulsive control strategy and established conditions for extinction of prey and predator population. The problem of adaptive control and chaotic behavior of a continuous-time two-prey and one-predator system is studied by El-Gohary and Al-Ruzaiza [12] with nonlinear feedback.Many different food chain models including three species are constructed. Das et al. [13] studied chaotic dynamics of a three-species competition model with noise. Dubey and Upadhyay [14] studied the interaction of two predators competing for a single prey species with ratio-dependent-type functional response and established the conditions for persistence and extinction of the species. Song et al. [15] explored the combined influence of multiple delays with two different scenarios, i.e., the intrinsic growth rate of the resource is less/more than those of the consumer and predator in a three-species model with consumer-eat-resource and predator-eat-consumer. The significance of top predator interference and gestation delay is studied by Jana et al. [16] on a three-species food chain model involving intermediate and top predator populations, and they observed the subcritical Hopf bifurcation phenomena. Wang and his group explored the dynamics of different stochastic epidemic models like SIRS and FIV with different infectious forces under intervention strategies and environmental variability [17–20].Motivated by the above works, we explain the dynamical behavior of the undertaken delayed model system and study the effect of control parameters with time delay on the dynamics of the model system. Since the model represents the very common food chain model of real-world systems, we are taking one of them as an example—a species found most throughout North America and from southern Canada to central Mexico, known as bobcats. Its main prey are rabbits and squirrels, which depend on grass and leaves. Kittens may be taken as prey by adult bobcats when prey levels are low, but this is very rare and does not influence much the dynamics of the population [21]. Some species like golden eagle and coyotes used to predate on bobcats. Diseases, accidents, hunters, and starvation are the leading causes of death of the species. A bobcat is able to survive for long periods without food but will eat heavily when prey is abundant and used to store their hunted prey sometimes to eat later. The gestation period of the species is generally from 60 to 70 days. We study this time lag by introducing the time delay. In our work, interspecific competition coefficient of both prey species along with the intraspecific interference coefficient of the predator is considered. The effect of all of these parameters is studied. We have considered Holling type II functional response as predation of both the prey and studied the effect of predation rates b13,b23 with the conversion rate of prey into predator β31,β32 on the model system. The death rate of predators is also affecting the dynamics of the model system significantly. All these parameters are considered in numerical simulations to show the behavior and changes in the dynamics of the model system.The rest of the organization of this paper is as follows: In Section2, the model system is formulated. Also, positive invariance and boundedness have been shown in this section. Stability analysis for both delayed and nondelayed model systems has been made in Section 3. The conditions for the existence of Hopf bifurcation are also performed here. Next, Section 4 is devoted to establish the formulae to determine the direction, period, and stability of Hopf bifurcation of the bifurcating periodic solution and some numerical simulations are carried out to illustrate the analytical results in Section 5. Finally, we end the paper with a brief conclusion and discussion in Section 6.


## 2. The Mathematical Model

Consider a Gause-type model to draw the scenario of an ecological system where a predator (bobcats) speciesXt predates two competing prey species N1 (squirrels) and N2 (rabbits) (see Figure 1). It is well known that gestation delay, as a time lag, plays a complicated role on the dynamics of the prey-predator system. Thus, the considered system can be governed by the following delay differential equations [2, 5]:
(1)dN1dt=N1r11−N1K1−b13X1+a1N1−α12N2,dN2dt=N2r21−N2K2−b23X1+a2N2−α21N1,dXdt=Xβ31N1t−τ1+a1N1t−τ+β32N2t−τ1+a2N2t−τ−δ−γX,with the initial conditions
(2)N1θ=ϕ1θ≥0,N2θ=ϕ2θ≥0,Xθ=ϕ3θ≥0,θ∈−τ,0,ϕi0>0,i=1,2,3,where ϕ:−τ,0→ℜ3 with the norm
(3)ϕ=sup−τ≤θ≤0ϕ1θ,ϕ2θ,ϕ3θ,such that ϕ=ϕ1,ϕ2,ϕ3.Figure 1
Prey-predator model system.Here, all the parameters are positive.r1,K1,α12, r2,K2,α21 are the intrinsic growth rate, carrying capacity, and interspecific competition coefficient of the prey species N1,N2, respectively, b13,β31, b23,β32 are the feeding rates per predator per unit prey consumed, assimilation, or conversion efficiency of the predator of the prey species N1,N2, respectively, and δ,γ are the death rate and intraspecific interference coefficient of the predator species X in the absence of the prey. The physical interpretation of ai is that ai=1/Ai, where Ai is the half saturation constant of the prey species Ni. The values of all the parameters are positive.


### 2.1. Positive Invariance

It is necessary to insure that the species always survive. For this, we need to show that the system is positively invariant.Theorem 2.1.
All the solutions of (1) are positive with positive initial conditions.Proof.
We haveN10>0,N20>0,X0>0 from the initial conditions and
(4)N1t=N10exp∫0tr11−N1K1−b13X1+a1N1−α12N2ds,N2t=N20exp∫0tr21−N2K2−b23X1+a2N2−α21N1ds,Xt=X0exp∫0tβ31N1t−τ1+a1N1t−τ+β32N2t−τ1+a2N2t−τ−δ−γXds.Hence, the result follows.


### 2.2. Boundedness

Theorem 2.2.
For any positive solution ϕt=N1t,N2t,Xt of system (1), there exists a T>0 such that 0<N1t≤M1,0<N2t≤M2 and 0<Xt≤M3 for t≥T, where
(5)0<M1=K1,0<M2=K2,0<M3=1γβ31K1+β32K2.Proof.
Obviously, from Theorem2.1, N1t,N2t,Xt remains nonnegative. So, we need to show that N1t≤M1,N2≤M2 and Xt≤M3 only.
From the first equation of model system (1), we get
(6)dN1dt≤r1N11−N1K1.Therefore, by using the standard comparison rule, we have(7)limsupt→+∞N1t≤K1.Similarly, from the second equation of model system (1), we get
(8)dN2dt≤r2N21−N2K2.Therefore, by using the standard comparison rule, we have(9)limsupt→+∞N2t≤K2.Now, from the third equation of model system (1), we have
(10)dXdt≤Xβ31K1+β32K2−γX2=β31K1+β32K2X1−XL,L=1γXβ31K1+β32K2.Therefore, by using the standard comparison rule, we have(11)limsupt→+∞Xt≤1γβ31K1+β32K2.


## 2.1. Positive Invariance

It is necessary to insure that the species always survive. For this, we need to show that the system is positively invariant.Theorem 2.1.
All the solutions of (1) are positive with positive initial conditions.Proof.
We haveN10>0,N20>0,X0>0 from the initial conditions and
(4)N1t=N10exp∫0tr11−N1K1−b13X1+a1N1−α12N2ds,N2t=N20exp∫0tr21−N2K2−b23X1+a2N2−α21N1ds,Xt=X0exp∫0tβ31N1t−τ1+a1N1t−τ+β32N2t−τ1+a2N2t−τ−δ−γXds.Hence, the result follows.


## 2.2. Boundedness

Theorem 2.2.
For any positive solution ϕt=N1t,N2t,Xt of system (1), there exists a T>0 such that 0<N1t≤M1,0<N2t≤M2 and 0<Xt≤M3 for t≥T, where
(5)0<M1=K1,0<M2=K2,0<M3=1γβ31K1+β32K2.Proof.
Obviously, from Theorem2.1, N1t,N2t,Xt remains nonnegative. So, we need to show that N1t≤M1,N2≤M2 and Xt≤M3 only.
From the first equation of model system (1), we get
(6)dN1dt≤r1N11−N1K1.Therefore, by using the standard comparison rule, we have(7)limsupt→+∞N1t≤K1.Similarly, from the second equation of model system (1), we get
(8)dN2dt≤r2N21−N2K2.Therefore, by using the standard comparison rule, we have(9)limsupt→+∞N2t≤K2.Now, from the third equation of model system (1), we have
(10)dXdt≤Xβ31K1+β32K2−γX2=β31K1+β32K2X1−XL,L=1γXβ31K1+β32K2.Therefore, by using the standard comparison rule, we have(11)limsupt→+∞Xt≤1γβ31K1+β32K2.


## 3. Stability Analysis and Existence of Hopf Bifurcation

Model system (1) has seven possible nonnegative equilibrium points, namely, E00,0,0,E1K1,0,0,E20,K2,0,E3N~1,N~2,0,E4N^1,0,X^,E50,N⌣2,X⌣,E∗N1∗,N2∗,X∗. Here, we are omitting the existence criteria of all these equilibrium points as it has been shown in Upadhyay and Iyenger [2].


### 3.1. Linear Stability Analysis

In this section, we will discuss the stability criteria of all abovementioned equilibrium points.(i)
The eigenvalues of the variational matrix aroundE0 are r1,r2,−δ. Hence, E0 is a saddle point
(ii)
The eigenvalues of the variational matrix aroundE1 are −r1,r2−α21K1 and K1β31/1+a1K1−δ. If r2<α21K1 and K1β31/1+a1K1<δ, then E1 is asymptotically stable else unstable
(iii)
The eigenvalues of the variational matrix aroundE2 are −r2,r1−α12K2 and K2β32/1+a2K2−δ. If r1<α12K2 and K2β32/1+a2K2<δ, then E2 is asymptotically stable else unstable
(iv)
The predator-free equilibrium point is given asE3N~1,N~2,0 where N~1=r2K1r1−α12K2/r2r1−α21α12K1K2 and N~2=r1K2r2−α21K1/r2r1−α21α12K1K2. Thus, equilibrium point E3 exists if r1−α12K2 and r2−α21K1 are of the same sign that is either
(12)r1>α12K2,r2>α21K1,or
(13)r1<α12K2,r2<α21K1
One eigenvalue of the variational matrix aroundE3 is
(14)λ1=β31N~11+a1N~1+β32N~21+a2N~2−δ,and the other two eigenvalues are roots of the following equation:
(15)λ2+A1λ+A2=0,where
(16)A1=r1N~1K1+r2N~2K2>0,A2=r1r2K1K2−α12α21N~1N~2.We note thatB>0 under condition (12). Thus, the equilibrium point E3N~1,N~2,0 is locally asymptotically stable if λ1<0 and E3 exists under condition (12).Again,E3 is unstable in the positive direction orthogonal to the N1−N2 plane, i.e., the X direction if λ1>0
(v)
The equilibrium point is given asE4N^1,0,X^ where X^=1/γβ31N^1/1+a1N^1−δ and N^1 is a positive solution of the following equation:
(17)b1N^13+b2N^12+b3N^1+b4=0,where
(18)b1=γr1a12,b3=γr11−2K1a1+K1b13β31−δa1,b2=γr1a12−K1a1,b4=−K1γr1+δb13.By Descarte’s rule of sign, (17) has a unique positive solution N1=N^1 if b2>0, i.e., if
(19)K1<2a1.One eigenvalue of the variational matrix aroundE4 is
(20)λ1=r2−b23X^−α21N^1,and the other two eigenvalues are the roots of the following characteristic equation
(21)λ2+B1λ+B2=0,where
(22)B1=γX^+rK1−a1b13X^1+a1N^12N^1,B2=γrK1−a1b13X^1+a1N^12X^N^1+b13β13N^1X^1+a1N^13.Thus, the other two eigenvalues have negative real parts if(23)rK1>a1b13X^1+a1N^12.It follows that the equilibriumE4N^1,0,X^ is locally asymptotically stable if λ1 (given by (20)) is negative and condition (23) holds. Also E4 is unstable in the positive direction orthogonal to the N1−X plane, i.e., in the N2 direction if λ1>0
(vi)
The equilibrium point is given asE50,N⌣2,X⌣ where X⌣=1/γβ32N⌣2/1+a2N⌣2−δ and N⌣2 is a positive solution of the following equation:
(24)l1N⌣23+l2N⌣22+l3N⌣2+l4=0,where
(25)l1=γr2a22,l3=γr21−2K2a2+K2b23β32−δa2,l2=γr2a22−K2a2,l4=−K2γr2+δb23.Similar toE4, it can be checked that the equilibrium E50,N⌣2,X⌣ exists if
(26)K2<2a2,β32N⌣21+a2N⌣2>δ.Again,E5 is locally asymptotically stable if
(27)r2K2>a2b23X⌣1+a2N⌣22,λ∗=r1−b13X⌣−α12N⌣2.Ifλ∗>0, then E5 is unstable in the N1 direction
(vii)
For the equilibrium pointE∗N1∗,N2∗,X∗, we discuss both the cases separatelyCase 1.
Whenτ=0.The eigenvalues of the “J” matrix around E∗ are the roots of the following equation:
(28)λ3+A1λ2+A2λ+A3=0,where
(29)A1=−a11+a22+a33,A2=a22a33+a11a22+a11a33−a23a32−a21a12−a31a13,A3=−a11a22a33−a32a23+a21a13a32−a12a33+a31a12a23−a13a22,and
(30)a11=−r1N1∗K1+a1b13N1∗X∗1+a1N1∗2,a12=−α12N1∗,a13=−b13N1∗1+a1N1∗,a21=−α21N2∗,a22=−r2N2∗K2+a2b23N2∗X∗1+a2N2∗2,a23=−b23N2∗1+a2N2∗,a31=β31X∗1+a1N1∗2,a32=β32X∗1+a2N2∗2,a33=−γX∗.From the Routh-Hurwitz criterion, we can state the following theorem.Theorem 3.1.
The equilibrium point E∗ is locally asymptotically stable if A1,A3>0 and A1A2>A3.Case 2.
Whenτ≠0.We linearize model system (1) at E∗N1∗,N2∗,X∗ and obtain xt=N1t−N1∗,yt=N2t−N2∗,zt=Xt−X∗ as the perturbed variables. Then the linearized system of model (1) is given by
(31)ddtxtytzt=M1xtytzt+M2xt−τyt−τzt−τ,where
(32)M1=−r1N1∗K1+a1b13N1∗X∗1+a1N1∗2−a12N1∗−b13N1∗1+a1N1∗−a21N2∗−r2N2∗K2+a2b23N2∗X∗1+a2N2∗2−b23N2∗1+a2N2∗00−γX∗,and
(33)M2=000000β31X∗1+a1N1∗2β32X∗1+a2N2∗20.


### 3.2. Existence of Hopf Bifurcation

The characteristic equation of the linearized system aroundE∗N1∗,N2∗,X∗ is given by
(34)M1+M2e−λτ−λI=0,i.e., it takes the form given as
(35)λ3+Aλ2+Bλ+C+e−λτDλ+E=0,where
(36)A=−a11+a22+a33,B=a11a22+a11a33+a22a33−a12a21,C=a12a21a33−a11a22a33,D=−a13b31+a23b32,E=a13a22b31+a11a23b32−a12a23b31+a13a21b32.Letλ=iω0,τ=τ0, substituting this into (35), for the sake of simplicity, denoting ω0 and τ0 by ω,τ.
(37)iω3+Aω2−iBω−C−iDωcosωτ−Ecosωτ−Dωsinωτ+iEsinωτ=0.Separating real and imaginary parts, we have(38)ω3−Bω=Dωcosωτ−Esinωτ,Aω2−C=Ecosωτ+Dωsinωτ.The following equation is obtained by some mathematical calculations(39)cosωτ=1E−Dωsinωτ+Aω2−C,(40)sinωτ=DA−Eω3+BE−DCωD2ω2+E2.Assin2ωτ+cos2ωτ=1, from (38), we have
(41)ω6+k1ω4+k2ω2+k3=0,where
(42)k1=A2−2B,k2=B2−2AC−D2,k3=C2−E2.Without loss of generality, assume that it has at least one positive real rootω0, for which from (40), we have
(43)τk=1ω0sin−1DA−Eω3+BE−DCωD2ω2+E2+2kπω0,where k=0,1,2,…, and then ±iω0 are a pair of purely imaginary roots of (35) with τk.Let us define(44)τ0=minτk,i=0,1,2….Also, differentiating (34) with respect to τ, we get
(45)Redλdτ−1τ=τ0=Usinωτ+Vcosωτ−D2ω2E2+D2ω2ω2,where
(46)U=2ADω3+BEω−3Eω3,V=3Dω4+2AEω2−BDω2.Hence, the transversality condition is satisfied ifUsinωτ+Vcosωτ−D2ω2>0. Thus, we can now state the following theorem.Theorem 3.2.
The equilibrium point E∗ of model system (1) is locally asymptotically stable when τ∈0,τ0 if the conditions of Theorem 3.1 hold and are unstable for τ>τ0. Furthermore, the system undergoes Hopf bifurcation at E∗ when the value of τ crosses through τ0 provided that Usinωτ+Vcosωτ−D2ω2>0.


### 3.3. Global Stability Analysis

Theorem 3.3.
The positive equilibrium point E∗ of system (1) is globally asymptotically stable provided that
(i)
ηii>0,i=1,2,3
(ii)
η11η22η33−η12η21η33−η11η23η32>0where
(47)η11=r1K1−b13a1X∗1+a1N1∗1+a1m1,η12=α12,η13=b13X∗1+a1M1,η21=α21,η22=r2K2−b23a2X∗1+a2N2∗1+a2m2,η23=b23X∗1+a2M2,η31=−β311+a1N1∗1+a1m1,η32=−β321+a2N2∗1+a2m2,η33=γ,and m1<N1t<M1, m2<N2t<M2, m3<Xt<M3 for t>0.Proof.
To prove that the equilibrium pointE∗x∗,y∗,z∗ of model system (1) is globally asymptotically stable, the method of Lyapunov function is used. We derive a sufficient condition which guarantees that E∗ is globally asymptotically stable. For mathematical convenience, we make the following transformations of the following variables:
(48)N1t=N1∗ext,N2t=N2∗eyt,Xt=X∗ezt.These coordinate change transforms the positive equilibriumE∗ into the trivial equilibrium xt=yt=zt=0 for all t>0. Due to the above transformations, model system (1) is reduced as follows:
(49)dxdt=−r1N1∗K1+b13a1X∗N1∗1+a1N1∗1+a1N1ext−1−α12N2∗eyt−1−b131+a1N1ezt−1,dydt=−α21N1∗ext−1+−r2N2∗K2+b23a2X∗N2∗1+a2N2∗1+a2N2eyt−1−b231+a2N2ezt−1,dzdt=−γX∗ezt−1+β31N1∗1+a1N1∗1+a1N1t−τext−τ−1+β32N2∗1+a2N2∗1+a2N2t−τeyt−τ−1.LetV1=xt. Computing the upper right derivative of V1t along the solution of (1), we get
(50)D+V1t≤−r1N1∗K1+b13a1X∗N1∗1+a1N1∗1+a1N1ext−1−α12N2∗eyt−1−b131+a1N1ezt−1≤−r1N1∗K1+b13a1X∗N1∗1+a1N1∗1+a1m1ext−1−α12N2∗eyt−1−b131+a1M1ezt−1=−η11N1∗ext−1−η12N2∗eyt−1−η13X∗ezt−1.Again, letV2=yt. Computing the upper right derivative of V2t along the solution of (1), we obtain
(51)D+V2t≤−α21N1∗ext−1+−r2N2∗K2+b23a2X∗N2∗1+a2N2∗1+a2N2eyt−1−b231+a2N2ezt−1≤−α21N1∗ext−1+−r2N2∗K2+b23a2X∗N2∗1+a2N2∗1+a2m2eyt−1−b231+a2M2ezt−1=−η21N1∗ext−1−η22N2∗eyt−1−η23X∗ezt−1.Now, letV3=zt. Computing the upper right derivative of V3t along the solution of (1), we obtain
(52)D+V3t≤−γX∗ezt−1+β31N1∗1+a1N1∗1+a1N1t−τext−τ−1+β32N2∗1+a2N2∗1+a2N2t−τeyt−τ−1.Again, due to the structure of (51), we consider the following functional:
(53)V33t=V3t+β31N1∗1+a1N1∗1+a1N1t−τ∫t−τtexs−1ds+β32N2∗1+a2N2∗1+a2N2t−τ∫t−τteys−1ds,whose upper right derivative along the solution of system (1) is given by
(54)D+V33t≤D+V3t+β31N1∗1+a1N1∗1+a1N1t−τext−1−ext−τ−1+β32N2∗1+a2N2∗1+a2N2t−τeyt−1−eyt−τ−1≤−γX∗ezt−1+β31N1∗1+a1N1∗1+a1N1t−τext−1+β32N2∗1+a2N2∗1+a2N2t−τeyt−1≤−γX∗ezt−1+β31N1∗1+a1N1∗1+a1m1ext−1+β32N2∗1+a2N2∗1+a2m2eyt−1=−η31N1∗ext−1−η32N2∗eyt−1−η33X∗ezt−1.According to the assumptions of the theorem, we know thatη=ηij3×3 is an M matrix; hence, there exist positive constants hii=1,2,3 such that,
(55)η11h1+η12h2+η13h3>0,η21h1+η22h2+η23h3>0,η31h1+η32h2+η33h3>0.Let us consider the Lyapunov functionV defined by
(56)Vt=h1V1t+h2V2t+h3V33t.Then from (50), (51), and (54), we obtain
(57)dVdt≤−η11h1N1∗ext−1−η12h2N2∗eyt−1−η13h3X∗ezt−1−η21h1N1∗ext−1−η22h2N2∗eyt−1−η23h3X∗ezt−1−η31h1N1∗ext−1−η32h2N2∗eyt−1−η33h3X∗ezt−1.It follows from (55) that, for all t≥T+r, we have
(58)dVdt≤−l1N1∗ext−1−l2N2∗eyt−1−l3X∗ezt−1.Since system (1) is positively invariant, one can see that there exist positive constants m1,m2,m3 and T∗>T+r, such that N1∗ext=N1t≥m1, N2∗eyt=N2t≥m2, and X∗ezt=Xt≥m3 for t≥T∗. Also, system (1) is bounded, and hence, there exist positive constants M1,M2,M3, such that N1∗ext=N1t≤M1, N2∗eyt=N2t≤M2, and X∗ezt=Xt≤M3 for t≥T∗.Using the mean-valued theorem, one obtains,(59)N1∗ext−1=N1∗eθ1txt,N2∗eyt−1=N2∗eθ2tyt,X∗ezt−1=X∗eθ3tzt,where N1∗eθ1t lies between N1t and N1∗, N2∗eθ2t lies between N2t and N2∗, and X∗eθ3t lies between Xt and X∗.Lettingα=minm1l1,m2l2,m3l3, then it follows that, for t≥T∗,
(60)D+V≤−αxt+yt+zt.Noting that the Lyapunov function is such that(61)Vt≥minh1,h2,h3xt+yt+zt.Hence, by applying the method of Lyapunov function and (60), we can conclude that the zero solution of system (49) is globally asymptotically stable; therefore, the positive equilibrium of system (1) is globally asymptotically stable. The proof is complete.Note.
The positive equilibrium pointE∗ will also be globally asymptotically stable for nondelayed model system (1). The conditions under which E∗ remains globally asymptotically stable are stated in the following lemma.Lemma 3.1.
If the positive equilibrium point E∗ is locally asymptotically stable, then it is also globally asymptotically stable in the interior of the positive octant under the following conditions:
(62)r1K1>b13a1X∗1+a1N1∗1+a1m1,r2K2>b23a2X∗1+a2N2∗1+a2m2,(63)m122<m11m22,where
(64)m11=r1K1−b13a1X∗1+a1N1∗1+a1m1,m22=c1r2K2−b23a2X∗1+a2N2∗1+a2m2,m12=α12+c1α21,c1=β32c2b231+a2N2∗,c2=b131+a1N1∗β31,and m1<N1t<M1, m2<N2t<M2, and m3<Xt<M3 for t>0.Proof.
The proof of the lemma is referred to in Appendix.


### 3.4. Uniform Persistence of the Model System

Theorem 3.4.
Suppose that system (1) satisfies the following conditions:
(65)Li>0,fori=5,6,7,where
(66)L5=β31N~11+a1N~1+β32N~21+a2N~2−δ,L6=r2−b23X^−α21N^1,L7=r1−b13X⌣−α12N⌣2.
Then system (1) is uniformly persistent.Proof.
We prove this theorem by the method of suitable Lyapunov function [22]. Let the Lyapunov function for system (1) be
(67)σN1t,N2t,Xt=N1p1N2p2Xp3,where p1,p2, and p3 are positive constants. Clearly, σ is a nonnegative C1 function defined in R+3. Then we have
(68)ψN1t,N2t,Xt=σ̇N1t,N2t,XtσN1t,N2t,Xt=p1Ṅ1N1+p2Ṅ2N2+p3ẊX=p1r11−N1K1−b13X1+a1N1−α12N2+p2r21−N2K2−b23X1+a2N2−α21N1+p3β31N1t−τ1+a1N1+β32N2t−τ1+a2N2−δ−γX.By computingψ at the boundary equilibria, we have
(69)ψ0,0,0=p1r1+p2r2−p3δ,ψK1,0,0=p2r2−α21K1+p3β31K11+a1K1−δ=p2L1+p3L2,ψ0,K2,0=p1r1−α12K2+p3β32K21+a2K2−δ=p1L3+p3L4,ψN~1,N~2,0=p3β31N~11+a1N~1+β32N~21+a2N~2−δ=p3L5,ψN^1,0,X^=p2r2−b23X^−α21N^1=p2L6,ψ0,N⌣2,X⌣=p1r1−b13X⌣−α12N⌣2=p1L7.If we choosep1 and p2 large enough so that p1r1+p2r2−p3δ>0, hence, ψ is positive at E0. Also, if conditions (12), (19), and (23) satisfy, then p2L1+p3L2>0 and p1L3+p3L4>0, i.e., ψ is positive at E1 and E2. And, if inequalities in (65) hold, then ψ is positive at all the remaining boundary equilibria. Therefore, system (1) is uniformly persistent, which follows from Theorem 3.12 of [22]. The proof is complete.


## 3.1. Linear Stability Analysis

In this section, we will discuss the stability criteria of all abovementioned equilibrium points.(i)
The eigenvalues of the variational matrix aroundE0 are r1,r2,−δ. Hence, E0 is a saddle point
(ii)
The eigenvalues of the variational matrix aroundE1 are −r1,r2−α21K1 and K1β31/1+a1K1−δ. If r2<α21K1 and K1β31/1+a1K1<δ, then E1 is asymptotically stable else unstable
(iii)
The eigenvalues of the variational matrix aroundE2 are −r2,r1−α12K2 and K2β32/1+a2K2−δ. If r1<α12K2 and K2β32/1+a2K2<δ, then E2 is asymptotically stable else unstable
(iv)
The predator-free equilibrium point is given asE3N~1,N~2,0 where N~1=r2K1r1−α12K2/r2r1−α21α12K1K2 and N~2=r1K2r2−α21K1/r2r1−α21α12K1K2. Thus, equilibrium point E3 exists if r1−α12K2 and r2−α21K1 are of the same sign that is either
(12)r1>α12K2,r2>α21K1,or
(13)r1<α12K2,r2<α21K1
One eigenvalue of the variational matrix aroundE3 is
(14)λ1=β31N~11+a1N~1+β32N~21+a2N~2−δ,and the other two eigenvalues are roots of the following equation:
(15)λ2+A1λ+A2=0,where
(16)A1=r1N~1K1+r2N~2K2>0,A2=r1r2K1K2−α12α21N~1N~2.We note thatB>0 under condition (12). Thus, the equilibrium point E3N~1,N~2,0 is locally asymptotically stable if λ1<0 and E3 exists under condition (12).Again,E3 is unstable in the positive direction orthogonal to the N1−N2 plane, i.e., the X direction if λ1>0
(v)
The equilibrium point is given asE4N^1,0,X^ where X^=1/γβ31N^1/1+a1N^1−δ and N^1 is a positive solution of the following equation:
(17)b1N^13+b2N^12+b3N^1+b4=0,where
(18)b1=γr1a12,b3=γr11−2K1a1+K1b13β31−δa1,b2=γr1a12−K1a1,b4=−K1γr1+δb13.By Descarte’s rule of sign, (17) has a unique positive solution N1=N^1 if b2>0, i.e., if
(19)K1<2a1.One eigenvalue of the variational matrix aroundE4 is
(20)λ1=r2−b23X^−α21N^1,and the other two eigenvalues are the roots of the following characteristic equation
(21)λ2+B1λ+B2=0,where
(22)B1=γX^+rK1−a1b13X^1+a1N^12N^1,B2=γrK1−a1b13X^1+a1N^12X^N^1+b13β13N^1X^1+a1N^13.Thus, the other two eigenvalues have negative real parts if(23)rK1>a1b13X^1+a1N^12.It follows that the equilibriumE4N^1,0,X^ is locally asymptotically stable if λ1 (given by (20)) is negative and condition (23) holds. Also E4 is unstable in the positive direction orthogonal to the N1−X plane, i.e., in the N2 direction if λ1>0
(vi)
The equilibrium point is given asE50,N⌣2,X⌣ where X⌣=1/γβ32N⌣2/1+a2N⌣2−δ and N⌣2 is a positive solution of the following equation:
(24)l1N⌣23+l2N⌣22+l3N⌣2+l4=0,where
(25)l1=γr2a22,l3=γr21−2K2a2+K2b23β32−δa2,l2=γr2a22−K2a2,l4=−K2γr2+δb23.Similar toE4, it can be checked that the equilibrium E50,N⌣2,X⌣ exists if
(26)K2<2a2,β32N⌣21+a2N⌣2>δ.Again,E5 is locally asymptotically stable if
(27)r2K2>a2b23X⌣1+a2N⌣22,λ∗=r1−b13X⌣−α12N⌣2.Ifλ∗>0, then E5 is unstable in the N1 direction
(vii)
For the equilibrium pointE∗N1∗,N2∗,X∗, we discuss both the cases separatelyCase 1.
Whenτ=0.The eigenvalues of the “J” matrix around E∗ are the roots of the following equation:
(28)λ3+A1λ2+A2λ+A3=0,where
(29)A1=−a11+a22+a33,A2=a22a33+a11a22+a11a33−a23a32−a21a12−a31a13,A3=−a11a22a33−a32a23+a21a13a32−a12a33+a31a12a23−a13a22,and
(30)a11=−r1N1∗K1+a1b13N1∗X∗1+a1N1∗2,a12=−α12N1∗,a13=−b13N1∗1+a1N1∗,a21=−α21N2∗,a22=−r2N2∗K2+a2b23N2∗X∗1+a2N2∗2,a23=−b23N2∗1+a2N2∗,a31=β31X∗1+a1N1∗2,a32=β32X∗1+a2N2∗2,a33=−γX∗.From the Routh-Hurwitz criterion, we can state the following theorem.Theorem 3.1.
The equilibrium point E∗ is locally asymptotically stable if A1,A3>0 and A1A2>A3.Case 2.
Whenτ≠0.We linearize model system (1) at E∗N1∗,N2∗,X∗ and obtain xt=N1t−N1∗,yt=N2t−N2∗,zt=Xt−X∗ as the perturbed variables. Then the linearized system of model (1) is given by
(31)ddtxtytzt=M1xtytzt+M2xt−τyt−τzt−τ,where
(32)M1=−r1N1∗K1+a1b13N1∗X∗1+a1N1∗2−a12N1∗−b13N1∗1+a1N1∗−a21N2∗−r2N2∗K2+a2b23N2∗X∗1+a2N2∗2−b23N2∗1+a2N2∗00−γX∗,and
(33)M2=000000β31X∗1+a1N1∗2β32X∗1+a2N2∗20.


## 3.2. Existence of Hopf Bifurcation

The characteristic equation of the linearized system aroundE∗N1∗,N2∗,X∗ is given by
(34)M1+M2e−λτ−λI=0,i.e., it takes the form given as
(35)λ3+Aλ2+Bλ+C+e−λτDλ+E=0,where
(36)A=−a11+a22+a33,B=a11a22+a11a33+a22a33−a12a21,C=a12a21a33−a11a22a33,D=−a13b31+a23b32,E=a13a22b31+a11a23b32−a12a23b31+a13a21b32.Letλ=iω0,τ=τ0, substituting this into (35), for the sake of simplicity, denoting ω0 and τ0 by ω,τ.
(37)iω3+Aω2−iBω−C−iDωcosωτ−Ecosωτ−Dωsinωτ+iEsinωτ=0.Separating real and imaginary parts, we have(38)ω3−Bω=Dωcosωτ−Esinωτ,Aω2−C=Ecosωτ+Dωsinωτ.The following equation is obtained by some mathematical calculations(39)cosωτ=1E−Dωsinωτ+Aω2−C,(40)sinωτ=DA−Eω3+BE−DCωD2ω2+E2.Assin2ωτ+cos2ωτ=1, from (38), we have
(41)ω6+k1ω4+k2ω2+k3=0,where
(42)k1=A2−2B,k2=B2−2AC−D2,k3=C2−E2.Without loss of generality, assume that it has at least one positive real rootω0, for which from (40), we have
(43)τk=1ω0sin−1DA−Eω3+BE−DCωD2ω2+E2+2kπω0,where k=0,1,2,…, and then ±iω0 are a pair of purely imaginary roots of (35) with τk.Let us define(44)τ0=minτk,i=0,1,2….Also, differentiating (34) with respect to τ, we get
(45)Redλdτ−1τ=τ0=Usinωτ+Vcosωτ−D2ω2E2+D2ω2ω2,where
(46)U=2ADω3+BEω−3Eω3,V=3Dω4+2AEω2−BDω2.Hence, the transversality condition is satisfied ifUsinωτ+Vcosωτ−D2ω2>0. Thus, we can now state the following theorem.Theorem 3.2.
The equilibrium point E∗ of model system (1) is locally asymptotically stable when τ∈0,τ0 if the conditions of Theorem 3.1 hold and are unstable for τ>τ0. Furthermore, the system undergoes Hopf bifurcation at E∗ when the value of τ crosses through τ0 provided that Usinωτ+Vcosωτ−D2ω2>0.


## 3.3. Global Stability Analysis

Theorem 3.3.
The positive equilibrium point E∗ of system (1) is globally asymptotically stable provided that
(i)
ηii>0,i=1,2,3
(ii)
η11η22η33−η12η21η33−η11η23η32>0where
(47)η11=r1K1−b13a1X∗1+a1N1∗1+a1m1,η12=α12,η13=b13X∗1+a1M1,η21=α21,η22=r2K2−b23a2X∗1+a2N2∗1+a2m2,η23=b23X∗1+a2M2,η31=−β311+a1N1∗1+a1m1,η32=−β321+a2N2∗1+a2m2,η33=γ,and m1<N1t<M1, m2<N2t<M2, m3<Xt<M3 for t>0.Proof.
To prove that the equilibrium pointE∗x∗,y∗,z∗ of model system (1) is globally asymptotically stable, the method of Lyapunov function is used. We derive a sufficient condition which guarantees that E∗ is globally asymptotically stable. For mathematical convenience, we make the following transformations of the following variables:
(48)N1t=N1∗ext,N2t=N2∗eyt,Xt=X∗ezt.These coordinate change transforms the positive equilibriumE∗ into the trivial equilibrium xt=yt=zt=0 for all t>0. Due to the above transformations, model system (1) is reduced as follows:
(49)dxdt=−r1N1∗K1+b13a1X∗N1∗1+a1N1∗1+a1N1ext−1−α12N2∗eyt−1−b131+a1N1ezt−1,dydt=−α21N1∗ext−1+−r2N2∗K2+b23a2X∗N2∗1+a2N2∗1+a2N2eyt−1−b231+a2N2ezt−1,dzdt=−γX∗ezt−1+β31N1∗1+a1N1∗1+a1N1t−τext−τ−1+β32N2∗1+a2N2∗1+a2N2t−τeyt−τ−1.LetV1=xt. Computing the upper right derivative of V1t along the solution of (1), we get
(50)D+V1t≤−r1N1∗K1+b13a1X∗N1∗1+a1N1∗1+a1N1ext−1−α12N2∗eyt−1−b131+a1N1ezt−1≤−r1N1∗K1+b13a1X∗N1∗1+a1N1∗1+a1m1ext−1−α12N2∗eyt−1−b131+a1M1ezt−1=−η11N1∗ext−1−η12N2∗eyt−1−η13X∗ezt−1.Again, letV2=yt. Computing the upper right derivative of V2t along the solution of (1), we obtain
(51)D+V2t≤−α21N1∗ext−1+−r2N2∗K2+b23a2X∗N2∗1+a2N2∗1+a2N2eyt−1−b231+a2N2ezt−1≤−α21N1∗ext−1+−r2N2∗K2+b23a2X∗N2∗1+a2N2∗1+a2m2eyt−1−b231+a2M2ezt−1=−η21N1∗ext−1−η22N2∗eyt−1−η23X∗ezt−1.Now, letV3=zt. Computing the upper right derivative of V3t along the solution of (1), we obtain
(52)D+V3t≤−γX∗ezt−1+β31N1∗1+a1N1∗1+a1N1t−τext−τ−1+β32N2∗1+a2N2∗1+a2N2t−τeyt−τ−1.Again, due to the structure of (51), we consider the following functional:
(53)V33t=V3t+β31N1∗1+a1N1∗1+a1N1t−τ∫t−τtexs−1ds+β32N2∗1+a2N2∗1+a2N2t−τ∫t−τteys−1ds,whose upper right derivative along the solution of system (1) is given by
(54)D+V33t≤D+V3t+β31N1∗1+a1N1∗1+a1N1t−τext−1−ext−τ−1+β32N2∗1+a2N2∗1+a2N2t−τeyt−1−eyt−τ−1≤−γX∗ezt−1+β31N1∗1+a1N1∗1+a1N1t−τext−1+β32N2∗1+a2N2∗1+a2N2t−τeyt−1≤−γX∗ezt−1+β31N1∗1+a1N1∗1+a1m1ext−1+β32N2∗1+a2N2∗1+a2m2eyt−1=−η31N1∗ext−1−η32N2∗eyt−1−η33X∗ezt−1.According to the assumptions of the theorem, we know thatη=ηij3×3 is an M matrix; hence, there exist positive constants hii=1,2,3 such that,
(55)η11h1+η12h2+η13h3>0,η21h1+η22h2+η23h3>0,η31h1+η32h2+η33h3>0.Let us consider the Lyapunov functionV defined by
(56)Vt=h1V1t+h2V2t+h3V33t.Then from (50), (51), and (54), we obtain
(57)dVdt≤−η11h1N1∗ext−1−η12h2N2∗eyt−1−η13h3X∗ezt−1−η21h1N1∗ext−1−η22h2N2∗eyt−1−η23h3X∗ezt−1−η31h1N1∗ext−1−η32h2N2∗eyt−1−η33h3X∗ezt−1.It follows from (55) that, for all t≥T+r, we have
(58)dVdt≤−l1N1∗ext−1−l2N2∗eyt−1−l3X∗ezt−1.Since system (1) is positively invariant, one can see that there exist positive constants m1,m2,m3 and T∗>T+r, such that N1∗ext=N1t≥m1, N2∗eyt=N2t≥m2, and X∗ezt=Xt≥m3 for t≥T∗. Also, system (1) is bounded, and hence, there exist positive constants M1,M2,M3, such that N1∗ext=N1t≤M1, N2∗eyt=N2t≤M2, and X∗ezt=Xt≤M3 for t≥T∗.Using the mean-valued theorem, one obtains,(59)N1∗ext−1=N1∗eθ1txt,N2∗eyt−1=N2∗eθ2tyt,X∗ezt−1=X∗eθ3tzt,where N1∗eθ1t lies between N1t and N1∗, N2∗eθ2t lies between N2t and N2∗, and X∗eθ3t lies between Xt and X∗.Lettingα=minm1l1,m2l2,m3l3, then it follows that, for t≥T∗,
(60)D+V≤−αxt+yt+zt.Noting that the Lyapunov function is such that(61)Vt≥minh1,h2,h3xt+yt+zt.Hence, by applying the method of Lyapunov function and (60), we can conclude that the zero solution of system (49) is globally asymptotically stable; therefore, the positive equilibrium of system (1) is globally asymptotically stable. The proof is complete.Note.
The positive equilibrium pointE∗ will also be globally asymptotically stable for nondelayed model system (1). The conditions under which E∗ remains globally asymptotically stable are stated in the following lemma.Lemma 3.1.
If the positive equilibrium point E∗ is locally asymptotically stable, then it is also globally asymptotically stable in the interior of the positive octant under the following conditions:
(62)r1K1>b13a1X∗1+a1N1∗1+a1m1,r2K2>b23a2X∗1+a2N2∗1+a2m2,(63)m122<m11m22,where
(64)m11=r1K1−b13a1X∗1+a1N1∗1+a1m1,m22=c1r2K2−b23a2X∗1+a2N2∗1+a2m2,m12=α12+c1α21,c1=β32c2b231+a2N2∗,c2=b131+a1N1∗β31,and m1<N1t<M1, m2<N2t<M2, and m3<Xt<M3 for t>0.Proof.
The proof of the lemma is referred to in Appendix.


## 3.4. Uniform Persistence of the Model System

Theorem 3.4.
Suppose that system (1) satisfies the following conditions:
(65)Li>0,fori=5,6,7,where
(66)L5=β31N~11+a1N~1+β32N~21+a2N~2−δ,L6=r2−b23X^−α21N^1,L7=r1−b13X⌣−α12N⌣2.
Then system (1) is uniformly persistent.Proof.
We prove this theorem by the method of suitable Lyapunov function [22]. Let the Lyapunov function for system (1) be
(67)σN1t,N2t,Xt=N1p1N2p2Xp3,where p1,p2, and p3 are positive constants. Clearly, σ is a nonnegative C1 function defined in R+3. Then we have
(68)ψN1t,N2t,Xt=σ̇N1t,N2t,XtσN1t,N2t,Xt=p1Ṅ1N1+p2Ṅ2N2+p3ẊX=p1r11−N1K1−b13X1+a1N1−α12N2+p2r21−N2K2−b23X1+a2N2−α21N1+p3β31N1t−τ1+a1N1+β32N2t−τ1+a2N2−δ−γX.By computingψ at the boundary equilibria, we have
(69)ψ0,0,0=p1r1+p2r2−p3δ,ψK1,0,0=p2r2−α21K1+p3β31K11+a1K1−δ=p2L1+p3L2,ψ0,K2,0=p1r1−α12K2+p3β32K21+a2K2−δ=p1L3+p3L4,ψN~1,N~2,0=p3β31N~11+a1N~1+β32N~21+a2N~2−δ=p3L5,ψN^1,0,X^=p2r2−b23X^−α21N^1=p2L6,ψ0,N⌣2,X⌣=p1r1−b13X⌣−α12N⌣2=p1L7.If we choosep1 and p2 large enough so that p1r1+p2r2−p3δ>0, hence, ψ is positive at E0. Also, if conditions (12), (19), and (23) satisfy, then p2L1+p3L2>0 and p1L3+p3L4>0, i.e., ψ is positive at E1 and E2. And, if inequalities in (65) hold, then ψ is positive at all the remaining boundary equilibria. Therefore, system (1) is uniformly persistent, which follows from Theorem 3.12 of [22]. The proof is complete.


## 4. Properties of Hopf Bifurcation

In this section, the formulae for determining the direction, stability, and period of the bifurcating periodic solutions of system (1) are presented here. The method that we used is based on the normal form theory and center manifold theorem as described in Hassard et al. [23] to determine the abovementioned properties of Hopf bifurcation. In the previous section, we noted that for the critical value of time delay, τ=τ0, the family of periodic solutions bifurcates from the steady state under certain condition. Without loss of generality, we can say that (35) has a pair of purely imaginary roots ±iω0 at these critical values of τ and the system undergoes Hopf bifurcation. Hence, for any root of (35) of the form λτ=ντ−iωτ,ντk=0 and ωτk=ω0 and ∂ν/∂τ∣τ=τk≠0. Let τ=τ0+μ,μ∈R, so that μ=0 is a Hopf bifurcation value for the system. Define the space of continuous real-valued functions as C=C−1,0,R3. Let x1=N1−N1∗,x2=N2−N2∗,x3=X−X∗, and xit=xiτt for i=1,2,3; delay system (1) then transforms to functional differential equation in C as
(70)ẋ=Lμxt+fμ,xt,where xt=x1t,x2t,x3tT∈R3,xtθ=xt+θ,θ∈−1,0 and Lμ:C→R,f:C×R→R are given by
(71)Lμϕ=τ0+μM1ϕ0+M2ϕ−1,such that
(72)M1=−r1N1K1+α1N1b13X1+α1N12−N1α12−b13N11+α1N1−N2α21−r2N2K2+α2N2b23X1+α2N22−b23N21+α2N200−γX,M2=000000Xβ311+α1N12Xβ321+α2N220,and
(73)fμ,ϕ=τ0+μ−r1K1ϕ120−b13ϕ10ϕ301+α1ϕ10−α12ϕ10ϕ20−r2K2ϕ220−b23ϕ20ϕ301+α2ϕ20−α21ϕ10ϕ20β31ϕ1−1ϕ301+α1ϕ1−1+β32ϕ2−1ϕ301+α2ϕ2−1−γϕ320.By the Riesz representation theorem, there exists a functionηθ,μ whose components are of bounded variation for θ∈−1,0 such that
(74)Lμϕ=∫−10ϕθdηθ,μ,forϕ∈C.In view of (71), we can choose
(75)ηθ,μ=τ0+μM1δθ+M2δθ+1,where δθ is the Dirac delta function.Forϕ∈C1−1,0,R3, define
(76)Aμϕ=dϕθdθ,θ∈−1,0,∫−10ϕsdηs,μ=Lμϕ,θ=0,and
(77)Rμϕ=0,θ∈−1,0,fθ,μ,θ=0.System (70) is then equivalent to
(78)ẋt=Aμxt+Rμxt,where xtθ=xt+θ for θ∈−1,0.Forψ∈C10,1,R3∗, define
(79)A∗ψs=−dψsds,s∈0,1,∫−10ψ−tdηTt,0,s=0,and a bilinear inner product
(80)ψ,ϕ=ψ¯0.ϕ0−∫θ=−10∫ξ=0θψ¯ξ−θdηθϕξdξ,where ηθ = ηθ,0. Then A0 (from here onwards, we shall refer A0 as A) and A∗ are adjoint operators. Since ±iω0τ0 are the eigenvalues of A, they will also be the eigenvalues of A∗. We need to compute eigenvectors of A and A∗ corresponding to +iω0τ0 and −iω0τ0, respectively.Suppose thatqθ=α1,α2,α3Teiω0τ0θ be the eigenvector of A corresponding to eigenvalue iω0τ0, then
(81)Aqθ=iω0τ0qθ,which, for θ=0, by direct computation and choosing α1=1, we get
(82)α2=a21+a23α3iω0−a22,and
(83)α3=iω0−a11iω0−a22−a21a12a23a12+a13iω0−a22,where M1=aiji,j=1,2,3 and M2=biji,j=1,2,3.Similarly, assuming thatq∗s=Dα1∗,α2∗,α3∗eiω0τ0s be the eigenvectors of A∗ corresponding to eigenvalue −iω0τ0, we have
(84)A∗q∗s=−iω0τ0q∗s,where
(85)α1∗=1,α2∗=a13b31e−iω0τ0−iω0+a33iω0+a11a21iω0+a33−a23b31e−iω0τ0,and
(86)α3∗=−a13+α2∗a23iω0+a33.Using the normalization condition, i.e.,q∗s,qθ=1, we get
(87)q∗s,qθ=D¯1,α¯2∗,α¯3∗1,α2,α3T−∫−10∫ξ=0θD¯1,α¯2∗,α¯3∗e−iω0τ0ξ−θdηθ1,α2,α3Teiω0τ0ξdξ=D¯1+α2α¯2∗+α3α¯3∗−∫−101,α¯2∗,α¯3∗θeiω0τ0θdηθ1,α2,α3T=D¯1+α2α¯2∗+α3α¯3∗+τ0α¯3∗b31+b32α2e−iω0τ0,which gives
(88)D¯=1P,where
(89)P=1+α2α¯2∗+α3α¯3∗+τ0α¯3∗b31+b32α2e−iω0τ0.Proceeding in the same manner as Hassard et al. [23], we compute the coordinates to describe the center manifold C0 at μ=0. Let xt be the solution of (78) when μ=0. Define
(90)zt=q∗,xt,Wt,θ=xtθ−2Reztqθ.On the center manifoldC0, we have
(91)Wt,θ=Wz,z¯,θ=W20θz22+W11θzz¯+W02θz¯22+W30θz36+⋯,where z and z¯ are local coordinates for center manifold C0 in the direction of q∗ and q¯∗. We now consider only the real solution xt∈C0 of (78), which gives
(92)ż=iω0τ0z+q¯∗.f0,Wz,z¯,0+2Rezq0=iω0τ0z+gz,z¯,where
(93)gz,z¯=q¯∗0f0z,z¯=g20z22+g11zz¯+g02z¯22+g21z2z¯2+⋯.It follows from (90) and (91) that
(94)xtθ=Wz,z¯,θ+2Rezqθ=W20θz22+W11θzz¯+W02θz¯22+zα1,α2,α3Teiω0τ0θ+z¯α¯1,α¯2,α¯3Te−iω0τ0θ+⋯,so that
(95)x1t0=z+z¯+W2010z22+W1110zz¯+W0210z¯22+⋯,x2t0=α2z+α¯2z¯+W2020z22+W1120zz¯+W0220z¯22+⋯,x3t0=α3z+α¯3z¯+W2030z22+W1130zz¯+W0230z¯22+⋯,x1t−1=ze−iω0τ0+z¯eiω0τ0+W201−1z22+W111−1zz¯+W021−1z¯22+⋯,x2t−1=α2ze−iω0τ0+α¯2z¯eiω0τ0+W202−1z22+W112−1zz¯+W022−1z¯22+⋯,x3t−1=α3ze−iω0τ0+α¯3z¯eiω0τ0+W203−1z22+W113−1zz¯+W023−1z¯22+⋯.Now, from (92), we have
(96)gz,z¯=τ0D¯α¯1∗,α¯2∗,α¯3∗−r1K1ϕ120−b13ϕ10ϕ301+a1ϕ10−α12ϕ10ϕ20−r2K2ϕ220−b23ϕ20ϕ301+a2ϕ20−α21ϕ10ϕ20β31ϕ1−1ϕ301+a1ϕ1−1+β32ϕ2−1ϕ301+a2ϕ2−1−γϕ320=τ0D¯z2−r1K1−α12α2−b13α3+α¯2∗−α21α2−α22r2K2−b23α2α3+α¯3∗−γα32+β31α3e−iω0τ0+β32α2α3e−iω0τ0+zz¯−2r1K1−α12Reα2−b13Reα3+α¯2∗−α21Reα2−2α22r2K2−b23α2α¯3+α¯2α3+α¯3∗−2γα32+β31α¯3e−iω0τ0+α3eiω0τ0+β32α2α¯3e−iω0τ0+α¯2α3eiω0τ0+z¯2−r1K1−α12α¯2−b13α¯3+α¯2∗−α21α¯2−α¯22r2K2−b23α¯2α¯3+α¯3∗−γα¯32+β31α¯3eiω0τ0+β32α¯2α¯3eiω0τ0+z2z¯−r1K12W1110+W2010−α12W1120+α¯2W20102+W20202+α2W1110−b13W1130+α¯3W20102+W20302+α3W1110+α¯2∗−α21W1120+α¯2W20102+W20202+α2W1110−r2K22α2W1120+α¯2W2020−b23α2W1130+α¯3W20202+α¯2W20302+α3W1120+α¯3∗−γ2α3W1130+α¯3W2030+β31e−iω0τ0W1130+α¯3W201−12+eiω0τ0W20302+α3W111−1+β32α2e−iω0τ0W1130+α¯3W202−12+α¯2eiω0τ0W20302+α3W112−1.Comparing the coefficient with (93), we have
(97)g20=2τ0D¯−r1K1−α12α2−b13α3+α¯2∗−α21α2−α22r2K2−b23α2α3+α¯3∗−γα32+β31α3e−iω0τ0+β32α2α3e−iω0τ0,g11=τ0D¯−2r1K1−α12Reα2−b13Reα3+α¯2∗−α21Reα2−2α22r2K2−b23α2α¯3+α¯2α3+α¯3∗−2γα32+β31α¯3e−iω0τ0+α3eiω0τ0+β32α2α¯3e−iω0τ0+α¯2α3eiω0τ0,g02=2τ0D¯−r1K1−α12α¯2−b13α¯3+α¯2∗−α21α¯2−α¯22r2K2−b23α¯2α¯3+α¯3∗−γα¯32+β31α¯3eiω0τ0+β32α¯2α¯3eiω0τ0,g21=2τ0D¯−r1K12W1110+W2010−α12W1120+α¯2W20102+W20202+α2W1110−b13W1130+α¯3W20102+W20302+α3W1110+α¯2∗−α21W1120+α¯2W20102+W20202+α2W1110−r2K22α2W1120+α¯2W2020−b23α2W1130+α¯3W20202+α¯2W20302+α3W1120+α¯3∗−γ2α3W1130+α¯3W2030+β31e−iω0τ0W1130+α¯3W201−12+eiω0τ0W20302+α3W111−1+β32α2e−iω0τ0W1130+α¯3W202−12+α¯2eiω0τ0W20302+α3W112−1.In order to computeg21, we need to compute W20θ and W11θ.From (90) and (92), we have
(98)Ẇ=ẋt−żq−z¯̇q¯=AW−2Req¯∗0f0qθ,θ∈−1,0,AW−2Req¯∗0f0qθ+f0,θ=0,(99)Ẇ=AW+Hz,z¯,θ,with
(100)Hz,z¯,θ=H20θz22+H11θzz¯+H02θz22+⋯.Also, onC0, using the chain rule, we get
(101)Ẇ=Wzż+Wz¯z¯̇.It follows from (92), (99), and (101) that
(102)A−2iω0τ0W20=−H20,(103)AW11=−H11,etc.Now, forθ∈−1,0,
(104)Hz,z¯,θ=−q¯∗0f0qθ−q∗0f¯0q¯θ=−gz,z¯qθ−g¯z,z¯q¯θ=−g20qθ+g¯02q¯θz22−g11qθ+g¯11q¯θzz¯+⋯,which in comparing the coefficients with (100) gives
(105)H20θ=−g20qθ−g¯02q¯θ,(106)H11θ=−g11qθ−g¯11q¯θ.From (102), (105), and the definition of A, we have
(107)Ẇ20θ=2iω0τ0W20θ+g20qθ+g¯02q¯θ.Note thatqθ=q0eiω0τ0θ; hence,
(108)W20θ=ig20ω0τ0qθ+ig¯023ω0τ0q¯θ+E1e2iω0τ0θ.Similarly, from (103), (106), and the definition of A, we have
(109)Ẇ11θ=g11qθ+g¯11q¯θ,(110)W11θ=−ig11ω0τ0qθ+ig¯11ω0τ0q¯θ+E2,where E1=E11,E12,E13 and E2=E21,E22,E23∈R4 are constant vectors, to be determined.It follows from the definition ofA and (102) that
(111)∫−10W20θdηθ=2iω0τkW200−H200,and
(112)∫−10W11θdηθ=−H110,where ηθ=η0,θ.From (98) and (100), we get
(113)H200=−g20q0−g¯02q¯0+2τ0H1,H2,H3T,and
(114)H110=−g11q0−g¯11q¯0+2τ0P1,P2,P3T,where
(115)H1=−r1K1−α12α2−b13α3,H2=−α21α2−α22r2K2−b23α2α3,H3=−γα32+β31α3e−iω0τ0+β32α2α3e−iω0τ0,P1=−2r1K1−α12Reα2−b13Reα3,P2=−α21Reα2−2α22r2K2−b23α2α¯3+α¯2α3,P3=−2γα32+β31α¯3e−iω0τ0+α3eiω0τ0+β32α2α¯3e−iω0τ0+α¯2α3eiω0τ0.Using (108) and (113) in (111) and noticing that
(116)iω0τ0I−∫−10eiω0τ0θdηθq0=0,and
(117)−iω0τ0I−∫−10e−iω0τ0θdηθq¯0=0,we obtain
(118)2iω0τ0I−∫−10e2iω0τ0θdηθE1=2τ0H1,H2,H3T.This leads to(119)2iω0−a11−a12−a13−a212iω0−a22−a23−b31e−2iω0τ0−b32e−2iω0τ02iω0−a33E1=2H1H2H3.Solving this system for K1, we obtain
(120)E11=2A~H1−a12−a13H22iω0−a22−a23H3−b32e−2iω0τ02iω0−a33,E12=2A~2iω0−a11H1−a13−a21H2−a23−b31e−2iω0τ0H32iω0−a33,E13=2A~2iω0−a11−a12H1−a212iω0−a22H2−b31e−2iω0τ0−b32e−2iω0τ0H3,where
(121)A~=2iω0−a11−a12−a13−a212iω0−a22−a23−b31e−2iω0τ0−b32e−2iω0τ02iω0−a33.Similarly, substituting (110) and (114) into (112), we get
(122)−a11−a12−a13−a21−a22−a23−b31−b32−a33E2=2P1P2P3,and hence,
(123)E21=2B~P1−a12−a13P2−a22−a23P3−b32−a33,E22=2B~−a11P1−a13−a21P2−a23−b31P3−a33,E23=2B~−a11−a12P1−a21−a22P2−b31−b32P3,where
(124)B~=−a11−a12−a13−a21−a22−a23−b31−b32−a33.Thus, we can determineW20θ and W11θ from (108) and (110). Furthermore, g21 can be expressed by the parameters and delay term. Thus, we can compute the following values:
(125)c10=i2ω0τ0g20g11−2g112−g0223+g212,μ2=−Rec10Reλ′τ0,β2=2Rec10,T2=−Imc10+μ2Imλ′τ0ω0τ0,which determine the qualities of bifurcating periodic solution in the center manifold at the critical value τ0. Now, we state the results in the following theorem.Theorem 4.
Suppose that μ2 determines the direction of the Hopf bifurcation. If μ2>0, then the Hopf bifurcation is supercritical and the bifurcating periodic solutions exit forτ>τ0. If μ2<0, then the Hopf bifurcation is subcritical and the bifurcating periodic solutions exit forτ<τ0. β2determines the stability of the bifurcating periodic solutions. The bifurcating periodic solutions are stable ifβ2<0and unstable ifβ2>0, andT2determines the period of the bifurcating periodic solutions. The period increases ifT2>0and decreases ifT2<0.


## 5. Numerical Simulation

In order to study the dynamic behavior of model system (1), we consider a set of parameter values as
(126)r1=4.8,r2=4,K1=95,K2=99,b13=0.27,b23=0.35,a1=1.4,a2=0.25,α12=0.02,α21=0.01,β31=1.45,β32=1.25,δ=5.3,γ=0.005.With respect to the considered parameter set, the nonzero equilibrium point is calculated asE∗77.3512,38.3201,50.7297. Also, corresponding with ω0 and the critical value of time delay, τ0 are given as 0.799559 and 0.20757, respectively. In Figure 2, we can see that the nondelay system is stable for this set of parameter values. Also, the system remains stable for τ=0.2<τ0=0.20757 and becomes unstable for τ=0.25>τ0 as τ0 increases. The system bifurcates at the critical value of delay, i.e., at τ0. From Figure 2, we can also see the effect of time delay that destabilizes the system as it becomes larger than the critical value. Now, we are interested to study the effect of the other parameters as well on the dynamics of the model system.Figure 2
(a) Phase plot of nondelayed model system (1). (b) Bifurcation diagram with respect to time delay which shows that the system bifurcates at τ=0.20757=τ0. (c) Time series of system (1) remains stable for τ=0.2<τ0. (d) Time series of system (1) becomes unstable at τ=0.25>τ0. Parameters are taken as given in the text.
(a)
(b)
(c)
(d)


### 5.1. Effect of Carrying Capacities of Prey PopulationsK1,K2

We found that the carrying capacityK1 of the first prey N1 does not affect the dynamics of the system (see Figure 3(a)) whereas the carrying capacity K2 of the second prey N2 changes the nature of the model system significantly. The system shows the limit cycle behavior from the stable state due to increase in K2 (see Figure 3(b)). Now, we have plotted the graph to show the relation between these parameters (K1 and K2) and the critical value of time delay. In both the cases, we observed that the value of τ0 decreases as K1 or K2 increases (see Figures 3(c) and 3(d)). It is interesting to note that the values of τ0 vary in the same ranges of K1 and K2. It is easy to check in Figure 4 that the system changes its dynamics as the value of time delay changes. In Figures 4(a) and 4(c), we can see that for some values of K1 or K2, the system shows stable focus and for some other values, it shows the limit cycle behavior accordingly as the value of τ=0.2 is less or greater than the respective value of τ0. Since τ=0.3 and τ=0.5 are more than the values of all τ0 for the respective K1 and K2 taken in Figures 4(b) and 4(d), the system shows the limit cycle behavior for each cases.Figure 3
Phase plot diagram of nondelayed model system (1) for different (a) K1 and (b) K2. Relation between the critical value of time delay and (c) K1 and (d) K2. Rest of the parameters is taken as given in the text.
(a)
(b)
(c)
(d)Figure 4
Phase plot diagram of system (1) for different (a) K1 at τ0=0.2, (b) K1 at τ0=0.3, (c) K2 at τ0=0.2, and (d) K2 at τ0=0.5. The rest of the parameters is taken as given in the text.
(a)
(b)
(c)
(d)


### 5.2.Effect of Intrinsic Growth Rate of Prey Populations r1,r2

In the same manner, we can see that the growth rate ofN2 only affects the dynamics of the system but not of N1. System dynamics change from the stable state to the limit cycle state as we increase the value of r2 but system remains unchanged when r1 is changed (see Figure 5). Now, we are checking the effect with the time delay and it is easy to see by comparing Figures 5(c) and 5(d) with Figure 5(a). The system remains stable at τ=0.15 while it turns into stable limit cycle when τ=0.3 is taken for all the values of r1 in both cases. While at τ=0.2 and τ=0.5 in Figures 5(e) and 5(f), resp., we observed that the system shows cyclic behavior for some r2 and remains stable for other values of r2 as in Figure 5(b).Figure 5
Phase plot diagram of nondelayed model system (1) for different (a) r1 and (b) r2. Phase plot diagram of system (1) for different r1 at (c) τ=0.15 and (d) τ=0.3. Phase plot diagram of system (1) for different r2 at (e) τ=0.2 and (f) τ=0.5. The rest of the parameters is taken as given in the text.
(a)
(b)
(c)
(d)
(e)
(f)


### 5.3. Effect of Predation Rate ofN1 and Conversion Efficiency of Predators b13,β31

In Figures6(a) and 6(b), we can see that the system changes its behavior from stable focus to limit cycle as we increased the value of β31 from 1.3 to 1.55 whereas b13 does not have effect on the dynamics of the system; it remains at stable focus. Moreover, it is interesting to see the reverse effect of increasing β31 and b13 in a particular range on the critical values of time delay (see Figures 6(c) and 6(d)). The respective critical value increases with an increase of b13 whereas it decreases with an increase of β31. Also, we observed that the decreasing value of β31 maintains the system to the stable state for larger time lag.Figure 6
Phase plot diagram of nondelayed model system (1) for different (a) β31 and (b) b13. Relation between the critical value of time delay and (c) β31 and (d) b13. The rest of the parameters is taken as given in the text.
(a)
(b)
(c)
(d)But, the predation rate ofN2 and conversion efficiency of predator b23,β32 do not follow any particular pattern for the respective critical values in any particular range. So, we study the effect of these parameters on the nondelayed model system only. In Figure 7(a), we can see that the system changes its behavior from stable focus to limit cycle as we are increasing the value of β32 from 1.05 to 1.35. Also, due to the increase in b23 from 0.15 to 0.45, the system loses its stability and shows limit cycle behavior (Figure 7(b)). We observed that the system becomes stable as we increase the death rate of the predator species (Figure 8(a)). As we increase δ from 5.1 to 5.4, the system becomes stable from limit cycle behavior. We also study the effect of intraspecific interference coefficient of the predator on the model system and found that the system becomes stable from limit cycle behavior as we increase the value (see Figure 8(b)). For a particular set of parameter values, the system shows limit cycle behavior for δ<5.3 and becomes stable for δ≥5.3. Similarly, to keep the system stable, we also need γ≥0.005; otherwise, the system shows limit cycle behavior for γ<0.005.Figure 7
Phase plot diagram of nondelayed model system (1) for different (a) β32 and (b) b23. The rest of the parameters is taken as given in the text.
(a)
(b)Figure 8
Phase plot diagram of nondelayed model system (1) for different (a) δ and (b) γ. The rest of the parameters are taken as given in the text.
(a)
(b)


## 5.1. Effect of Carrying Capacities of Prey PopulationsK1,K2

We found that the carrying capacityK1 of the first prey N1 does not affect the dynamics of the system (see Figure 3(a)) whereas the carrying capacity K2 of the second prey N2 changes the nature of the model system significantly. The system shows the limit cycle behavior from the stable state due to increase in K2 (see Figure 3(b)). Now, we have plotted the graph to show the relation between these parameters (K1 and K2) and the critical value of time delay. In both the cases, we observed that the value of τ0 decreases as K1 or K2 increases (see Figures 3(c) and 3(d)). It is interesting to note that the values of τ0 vary in the same ranges of K1 and K2. It is easy to check in Figure 4 that the system changes its dynamics as the value of time delay changes. In Figures 4(a) and 4(c), we can see that for some values of K1 or K2, the system shows stable focus and for some other values, it shows the limit cycle behavior accordingly as the value of τ=0.2 is less or greater than the respective value of τ0. Since τ=0.3 and τ=0.5 are more than the values of all τ0 for the respective K1 and K2 taken in Figures 4(b) and 4(d), the system shows the limit cycle behavior for each cases.Figure 3
Phase plot diagram of nondelayed model system (1) for different (a) K1 and (b) K2. Relation between the critical value of time delay and (c) K1 and (d) K2. Rest of the parameters is taken as given in the text.
(a)
(b)
(c)
(d)Figure 4
Phase plot diagram of system (1) for different (a) K1 at τ0=0.2, (b) K1 at τ0=0.3, (c) K2 at τ0=0.2, and (d) K2 at τ0=0.5. The rest of the parameters is taken as given in the text.
(a)
(b)
(c)
(d)


## 5.2.Effect of Intrinsic Growth Rate of Prey Populations r1,r2

In the same manner, we can see that the growth rate ofN2 only affects the dynamics of the system but not of N1. System dynamics change from the stable state to the limit cycle state as we increase the value of r2 but system remains unchanged when r1 is changed (see Figure 5). Now, we are checking the effect with the time delay and it is easy to see by comparing Figures 5(c) and 5(d) with Figure 5(a). The system remains stable at τ=0.15 while it turns into stable limit cycle when τ=0.3 is taken for all the values of r1 in both cases. While at τ=0.2 and τ=0.5 in Figures 5(e) and 5(f), resp., we observed that the system shows cyclic behavior for some r2 and remains stable for other values of r2 as in Figure 5(b).Figure 5
Phase plot diagram of nondelayed model system (1) for different (a) r1 and (b) r2. Phase plot diagram of system (1) for different r1 at (c) τ=0.15 and (d) τ=0.3. Phase plot diagram of system (1) for different r2 at (e) τ=0.2 and (f) τ=0.5. The rest of the parameters is taken as given in the text.
(a)
(b)
(c)
(d)
(e)
(f)


## 5.3. Effect of Predation Rate ofN1 and Conversion Efficiency of Predators b13,β31

In Figures6(a) and 6(b), we can see that the system changes its behavior from stable focus to limit cycle as we increased the value of β31 from 1.3 to 1.55 whereas b13 does not have effect on the dynamics of the system; it remains at stable focus. Moreover, it is interesting to see the reverse effect of increasing β31 and b13 in a particular range on the critical values of time delay (see Figures 6(c) and 6(d)). The respective critical value increases with an increase of b13 whereas it decreases with an increase of β31. Also, we observed that the decreasing value of β31 maintains the system to the stable state for larger time lag.Figure 6
Phase plot diagram of nondelayed model system (1) for different (a) β31 and (b) b13. Relation between the critical value of time delay and (c) β31 and (d) b13. The rest of the parameters is taken as given in the text.
(a)
(b)
(c)
(d)But, the predation rate ofN2 and conversion efficiency of predator b23,β32 do not follow any particular pattern for the respective critical values in any particular range. So, we study the effect of these parameters on the nondelayed model system only. In Figure 7(a), we can see that the system changes its behavior from stable focus to limit cycle as we are increasing the value of β32 from 1.05 to 1.35. Also, due to the increase in b23 from 0.15 to 0.45, the system loses its stability and shows limit cycle behavior (Figure 7(b)). We observed that the system becomes stable as we increase the death rate of the predator species (Figure 8(a)). As we increase δ from 5.1 to 5.4, the system becomes stable from limit cycle behavior. We also study the effect of intraspecific interference coefficient of the predator on the model system and found that the system becomes stable from limit cycle behavior as we increase the value (see Figure 8(b)). For a particular set of parameter values, the system shows limit cycle behavior for δ<5.3 and becomes stable for δ≥5.3. Similarly, to keep the system stable, we also need γ≥0.005; otherwise, the system shows limit cycle behavior for γ<0.005.Figure 7
Phase plot diagram of nondelayed model system (1) for different (a) β32 and (b) b23. The rest of the parameters is taken as given in the text.
(a)
(b)Figure 8
Phase plot diagram of nondelayed model system (1) for different (a) δ and (b) γ. The rest of the parameters are taken as given in the text.
(a)
(b)


## 6. Conclusions and Discussions

In this work, we have studied a Gause-type two-prey one-predator model system. We have discussed the stability properties of equilibrium points for the model system including positive invariance and boundedness. Nonlinear stability analysis for both delayed and nondelayed model systems has been established. Conditions for the existence of Hopf bifurcation at the nonzero equilibrium point of the delayed model system have been obtained. Criteria for the system to be in uniform persistence have been derived, and it is observed that delay due to gestation does not affect uniform persistence of the system. The formulae for the direction, stability, and period of bifurcating solution have been calculated using the normal form theory and center manifold theorem. Numerical simulations are presented to verify the analytical predictions. It illustrates that the time delay destabilizes the system depending on its strength. We have studied the effect of all the control parameters on the dynamics of the model system.It has been observed in our simulation experiments that some of the parameters change the dynamics of the system either into a stable state or to an unstable state as explained in the following:(i)
The conversion rateβ31 of prey N1 into predator is sufficient to break the stability of the nondelayed model system as it crosses through β31=1.49
(ii)
To maintain the stability of the nondelayed model system, we need to take (for a particular set of parameter) intraspecific interference coefficient,γ≥0.005, death rate of predator species, δ≥5.3, and carrying capacity of the second prey population, K2≤114
(iii)
Time lag increases (up to which the system remains in a stable state) with increase in the predation rate ofN1, 0.01≤b13≤0.5, whereas it decreases with an increase in the conversion rate of prey N1 into predator, 1.2≤β31≤1.48
(iv)
Time lag decreases with an increase in the carrying capacities of both prey populations (K1,K2) in the range of 80–110, but time delay up to which the system remains stable is larger for decreasing K2
(v)
The control parameters (r2,K2,b23) appearing in the growth equation of the second prey population have significant impact on the dynamics of the model system as compared to the control parameters (r1,K1,b13) for the first prey population. Increase in the values of K2 and b13 results in limit cycle behavior of the system, and stability of the system is achieved by increasing the growth rate (r2) of prey N2The rest of the parameters does not change the nature of the dynamical system much. We have also studied the effect of the growth rates of the prey population on the critical value of time delay and observed that the value ofτ0 decreases first up to a certain level for a particular value of r1 and r2, and then it starts to increase due to the increase in the growth rate of the prey population. The population of the bobcats depends primarily on the population of its prey, and it has been found generally stable and healthy [24]. Since the population is leading the wildcat species in the skin trade, the bobcat is legally hunted in small numbers. In some areas, they are quite rare, while in others, they have stable and sometimes dense populations. Hence, some states allow regulated hunting, while in others, they are protected. On a regional level, the bobcat is totally protected in ten USA states, in Canada, hunting and trade is regulated, and in Mexico, hunting is regulated in five states [25]. Using the above results, we can reach to the goal of making a healthy stable ecosystem and plan for regulating the bobcat population in North America and southern Canada.



---

*Source: 1589310-2018-11-07.xml*