# Stability and Hopf Bifurcation of Three-Species Prey-Predator System with Time Delays and Allee Effect

Authors: F. A. Rihan; H. J. Alsakaji; C. Rajivganthi
Journal: Complexity
 (2020)

---

## Abstract

Allee effect is one of the important factors in ecology, and taking it into account can cause significant impacts in the system dynamics. In this paper, we study the dynamics of a two-prey one-predator system, where the growth of both prey populations is subject to Allee effects, and there is a direct competition between the two-prey species having a common predator. Two discrete time delaysτ1 and τ2 are incorporated into the model to represent reaction time of predators. Sufficient conditions for local stability of positive interior equilibrium and existence of Hopf bifurcations in terms of threshold parameters τ1∗ and τ2∗ are obtained. A Lyapunov functional is deducted to investigate the global stability of positive interior equilibrium. Sensitivity analysis to evaluate the uncertainty of the state variables to small changes in the Allee parameters is also investigated. Presence of Allee effect and time delays in the model increases the complexity of the model and enriches the dynamics of the system. Some numerical simulations are provided to illustrate the effectiveness of the theoretical results. The model is highly sensitive to small changes in Allee parameters at the early stages and with low population densities, and this sensitivity decreases with time.


---

## Body

## 1. Introduction

The dynamical relationship between prey and their predators has long been and will continue to be one of the dominant themes in ecology due to its universal existence and importance (see, e.g., [1–5]). This relationship/interaction between two or more species has been essential in theoretical ecology since the famous Lotka–Volterra equations [6, 7], which are a system of first order, nonlinear differential equations that describe the dynamics and interactions between two or more species of biological systems. Of course, the qualitative properties of a prey-predator system such as stability of the steady states, bifurcation analysis, and oscillation of the solutions usually depend on the system parameters (see [8]).Suppose thatNt is the size of prey population and Pt is the size of the predator population at time t, then the Lotka–Volterra model is given by the following equations:(1)dNtdt=Ntβ1−γ1−g1Nt−eNtPt,dPtdt=Pt−γ+eNt,with N0>0 and P0>0. Here, β1 is the per capita maximum filtering rate and γ1 is the death rate of the prey Nt, while the parameter g1 denotes the strength of intraspecific competition. The predator death rate and predation rate are, respectively, denoted by γ and e. In the above model, it is assumed that prey population is subjected to logistic growth rate and the exponential type functional response.We should also mention here that one important component of prey-predator relationships is the functional response of predators to their prey(s)’ densities. The response of predators to different prey densities depends on the feeding behavior of individual predators. In [9], Holling discussed three different types of functional responses: Holling type I (linear), type II, and type III, etc. These responses are used to model the phenomena of predation, which captures the usual properties, for instance, positivity and increase (see also [10–13]).The authors believe that Allee effect and time delays greatly increase the likelihood of local and global extinction and can produce a rich variety of dynamic effects. It is a natural question that how the introduction of Allee effect in the prey growth function changes the system dynamics of the prey-predator system. However, before we introduce the final model, we give brief preliminaries about Allee effects and time delays in the prey-predator model (see [14, 15]).


### 1.1. Allee Effect

Allee effect was firstly reported by the American ecologist Allee [16], when he asked “what minimal numbers are necessary if a species is to maintain itself in nature?” Allee, in [16], shows that the growth rate is not always positive for small densities, and it may not be decreasing as in the logistic model either. In general, Allee effect mechanisms arise from cooperation or facilitation among individuals in the species [17]. A population is said to have an Allee effect if the growth rate per capita is initially an increasing function for the low density. It can be classified into two types: strong and weak. A strong Allee effect takes place when the population density is less than the specified threshold population considered, resulting in the species dying out. However, if the population density is greater than the threshold, the growth rate will remain positive [18], while a weak Allee effect means that the per capita growth rate cannot go below zero and remains positive.Now, we show how an Allee effect can be modelled, and how the per capita growth rate is affected with a weak Allee effect or a strong Allee effect throughout the simple examples:dN/dt=rN21−N/K for a weak Allee effect and dN/dt=rN1−N/KN/A−1 for a strong Allee effect.Figure1 shows a per capita growth rate 1/NdN/dt of the population with strong and weak Allee effects. The straight line shows the logistic growth, and the red curve displays a weak Allee effect, while the blue curve shows a strong Allee effect. The negative density dependence at low population sizes is described as a strong Allee effect, where there exists a threshold population level A, such that for N<A, 1/NdN/dt<0 (the species will die out) and for N>A, 1/NdN/dt>0, the growth will remain positive [18]. However, when the growth rate remains positive at low population densities, it is considered as a weak Allee effect.Figure 1
(a) The per-capita growth rate1/NdN/dt vs. population Nt with logistic (black dashes), strong (blue curve), and weak (red curve) Allee effects. (b) The population growth rate dN/dt vs. population Nt. For the strong Allee effect, the y-intercept of the per capita growth rate is less than zero at zero density, while in weak Allee effect, the y-intercept cannot go below zero.
(a)
(b)For multispecies models, there are flexible ways to formulate the Allee effects. For example, due to difficulties in finding mates when prey population density becomes low, Allee effect takes place in prey species. Herein, we propose and incorporate an additive Allee effect of the formbN≡N/α1+N in the prey growth function of model (1), which is considered as the probability of finding a mate (see [19]), so that(2)dNtdt=Ntβ1Ntα1+Nt−γ1−gNt−eNtPt,dPtdt=Pt−γ+eNt.The parameterα1 is the strength of Allee effect, α1=1/R, where R is the average area that can be searched by an individual [20]. We may notice that b0=0 and b′N>0 if N∈0,∞, i.e., Allee effect decreases as density increases, and limN⟶∞bN=1 means that Allee effect disappears at high densities. Therefore, the term bN is considered as a weak Allee effect function in a rectangular hyperbola form, known as Michaelis–Menten-like function [21].


### 1.2. Time Delays

Usually, the individuals of the prey and predator species usually pass through various life stages during their entire life span and the involved morphology differs from one stage to another. Construction of delay differential equation models is a well-known modelling strategy to take care of the stage-specific activities which are responsible for significant change in the dynamics of interacting populations. In various existing literature studies, the biological processes like incubation, gestation, maturation, and reaction time, are taken care of by introducing relevant time-delay parameters to the models for prey predator and other types of interacting populations. Incorporating time lags (or time delays) in biological models makes the systems much more realistic, as it can destabilize the equilibrium points and give rise to a stable limit cycle, causing oscillations to grow and enriching the dynamics of the model. Time delays have been considered and extremely studied by many authors in prey-predator models and biological systems (see [21–25]).Motivation to what we mentioned above, it is interesting and important to study the impact of time delays and Allee effect on the dynamics of three-species prey-predator models. In this paper, we extend the work in [26] and study the dynamics of a two-prey one-predator system, where the growth of both prey populations is subject to Allee effects, and there exists a direct competition between the two-prey species having a common predator. Two discrete time delays τ1 and τ2 are incorporated into the predator growth equation to represent the reaction time with each prey. Sensitivity analysis to evaluate the uncertainty of the state variables to small changes in the Allee parameters is also considered.The rest of this paper is organized as follows: Model formulation is presented in Section2. Local stability and bifurcation analysis of the steady states are discussed in Section 3, and global asymptotic stability of interior steady state is discussed in Section 4. We also utilize sensitivity functions to evaluate the sensitivity (uncertainty) of the state variables (preys and predator populations) to small changes in the severity Allee parameters through Section 5. Some numerical simulations are presented in Section 6 to show the effectiveness of the theoretical results. Finally, Section 7 concludes the study with a summary of the reported findings and future recommendations.


## 1.1. Allee Effect

Allee effect was firstly reported by the American ecologist Allee [16], when he asked “what minimal numbers are necessary if a species is to maintain itself in nature?” Allee, in [16], shows that the growth rate is not always positive for small densities, and it may not be decreasing as in the logistic model either. In general, Allee effect mechanisms arise from cooperation or facilitation among individuals in the species [17]. A population is said to have an Allee effect if the growth rate per capita is initially an increasing function for the low density. It can be classified into two types: strong and weak. A strong Allee effect takes place when the population density is less than the specified threshold population considered, resulting in the species dying out. However, if the population density is greater than the threshold, the growth rate will remain positive [18], while a weak Allee effect means that the per capita growth rate cannot go below zero and remains positive.Now, we show how an Allee effect can be modelled, and how the per capita growth rate is affected with a weak Allee effect or a strong Allee effect throughout the simple examples:dN/dt=rN21−N/K for a weak Allee effect and dN/dt=rN1−N/KN/A−1 for a strong Allee effect.Figure1 shows a per capita growth rate 1/NdN/dt of the population with strong and weak Allee effects. The straight line shows the logistic growth, and the red curve displays a weak Allee effect, while the blue curve shows a strong Allee effect. The negative density dependence at low population sizes is described as a strong Allee effect, where there exists a threshold population level A, such that for N<A, 1/NdN/dt<0 (the species will die out) and for N>A, 1/NdN/dt>0, the growth will remain positive [18]. However, when the growth rate remains positive at low population densities, it is considered as a weak Allee effect.Figure 1
(a) The per-capita growth rate1/NdN/dt vs. population Nt with logistic (black dashes), strong (blue curve), and weak (red curve) Allee effects. (b) The population growth rate dN/dt vs. population Nt. For the strong Allee effect, the y-intercept of the per capita growth rate is less than zero at zero density, while in weak Allee effect, the y-intercept cannot go below zero.
(a)
(b)For multispecies models, there are flexible ways to formulate the Allee effects. For example, due to difficulties in finding mates when prey population density becomes low, Allee effect takes place in prey species. Herein, we propose and incorporate an additive Allee effect of the formbN≡N/α1+N in the prey growth function of model (1), which is considered as the probability of finding a mate (see [19]), so that(2)dNtdt=Ntβ1Ntα1+Nt−γ1−gNt−eNtPt,dPtdt=Pt−γ+eNt.The parameterα1 is the strength of Allee effect, α1=1/R, where R is the average area that can be searched by an individual [20]. We may notice that b0=0 and b′N>0 if N∈0,∞, i.e., Allee effect decreases as density increases, and limN⟶∞bN=1 means that Allee effect disappears at high densities. Therefore, the term bN is considered as a weak Allee effect function in a rectangular hyperbola form, known as Michaelis–Menten-like function [21].


## 1.2. Time Delays

Usually, the individuals of the prey and predator species usually pass through various life stages during their entire life span and the involved morphology differs from one stage to another. Construction of delay differential equation models is a well-known modelling strategy to take care of the stage-specific activities which are responsible for significant change in the dynamics of interacting populations. In various existing literature studies, the biological processes like incubation, gestation, maturation, and reaction time, are taken care of by introducing relevant time-delay parameters to the models for prey predator and other types of interacting populations. Incorporating time lags (or time delays) in biological models makes the systems much more realistic, as it can destabilize the equilibrium points and give rise to a stable limit cycle, causing oscillations to grow and enriching the dynamics of the model. Time delays have been considered and extremely studied by many authors in prey-predator models and biological systems (see [21–25]).Motivation to what we mentioned above, it is interesting and important to study the impact of time delays and Allee effect on the dynamics of three-species prey-predator models. In this paper, we extend the work in [26] and study the dynamics of a two-prey one-predator system, where the growth of both prey populations is subject to Allee effects, and there exists a direct competition between the two-prey species having a common predator. Two discrete time delays τ1 and τ2 are incorporated into the predator growth equation to represent the reaction time with each prey. Sensitivity analysis to evaluate the uncertainty of the state variables to small changes in the Allee parameters is also considered.The rest of this paper is organized as follows: Model formulation is presented in Section2. Local stability and bifurcation analysis of the steady states are discussed in Section 3, and global asymptotic stability of interior steady state is discussed in Section 4. We also utilize sensitivity functions to evaluate the sensitivity (uncertainty) of the state variables (preys and predator populations) to small changes in the severity Allee parameters through Section 5. Some numerical simulations are presented in Section 6 to show the effectiveness of the theoretical results. Finally, Section 7 concludes the study with a summary of the reported findings and future recommendations.


## 2. Delayed Model with Allee Effect for the Two-Prey One-Predator System

Many studies have been done on multispecies prey-predator systems, including local and global bifurcations and different types of chaos (see, e.g., [26–29]). Sen et al. [26] discussed the Allee effect on two-preys’ growth function, where the predator is generalized. They explained how the Allee effect can suppress the chaotic dynamics and the route to chaos in prey growth by comparing it with a model without the Allee effect. In [27], the authors studied dynamics of three species (two preys and one predator) delayed prey-predator model with cooperation among the preys against predation. The growth rate for preys is thought to be logistic. Delays are taken just in the growth components for each of the species. Takeuchi and Adachi [28] considered two preys with logistic growth rates and an exponential functional response, where the predator survives on two-prey populations. Their results showed that the apparent chaotic behavior is a result of the periodic solution when one of the two-prey has greater competitive strength than the other. Song and Li [29] explored the dynamic behaviors of a Holling II two-prey one-predator system by introducing constant periodic releases of predators through periodically spraying a pesticide on the prey. They were then able to show that the system remains permanent under certain conditions.Herein, we generalize model (2) to a multispecies prey-predator system (two-preys one-predator). The model consists of two teams of preys with densities xt and yt, interacting with one team of predator with density zt. We also incorporate Allee effects in the growth functions of the two-prey populations, and there exists a direct competition between the two-prey species having a common predator.The model takes the general form:(3)dxtdt=xtβ1xtα1+xt−γ1−g1xt−αxtyt−extzt,dytdt=ytβ2ytα2+yt−γ2−g2yt−βxtyt−δytzt1+cyt,dztdt=−β3zt+εext−τ1zt−τ1+εδyt−τ2zt−τ21+cyt−τ2,with initial conditions(4)xθ=ϕ1θ>0,yθ=ϕ2θ>0,zθ=ϕ3θ>0,θ∈−τ,0,τ=maxτ1,τ2.Here,ϕiθ (i=1,2,3) are smooth initial functions. The coefficients α and β represent the coefficients of competition of two preys x and y, in the absence of predator. The description of rest of model parameters along with their symbols is presented in Table 1.Table 1
One biological meaning for the parameters of model (3).
Parameter
Description
α
1, α2
Strength of Allee effect
β
1, β2
Per capita maximum filtering rate of population
g
1, g2
Strength of intracompetition
γ
1, γ2
Death rate for preys
α, β
Coefficient of competition
e, δ
Decrease rate ofxt and yt due to predation by zt
β
3
Predator death rate
c
Magnitude of interference between the second type of prey
∈
An equal transformation rate of predator to preysxt and ytIt is reasonable to assume that the death (predation) of preys is instantaneous when attacked by their predator but their contribution to the growth of predator population must be delayed by some time delay. Therefore, we incorporated two discrete time delaysτ1 and τ2 in the reaction response functionals in the predator growth to represent the reaction time. The interaction between first species of prey and predator is assumed to be governed by Holling type I. While the interaction between the second species of prey and predator is assumed to be governed by Holling type II (cyrtoid functional) δytzt/1+cyt, response indicates that it is a hard-to-capture prey compared to the first species (see Figure 2).Figure 2
Mathematic scheme of the three-species predator prey system (3).To investigate the role of time delay and Allee effect on the dynamics of the system, we first discuss the boundedness and positivity of the solutions of system (3) with the given positive initial conditions (4).


### 2.1. Positivity and Boundedness of the Solution

The positivity of the solutions indicates the existence of the population, while the boundedness explains the natural control of growth due to the restriction of resources. We arrive at the following lemma.Lemma 1.
Every solution of system (3) corresponding to initial conditions (4) defined on 0,∞ remains positive for all t≥0, which satisfies(5)limt⟶∞supxt+yt≤κ,limt⟶∞supzt≤N,where κ=minβ1,β2 and N>0.Proof.
Model (3) can be represented in a matrix form(6)U˙t=FU,where U=x,y,zT∈ℝ3 and(7)FU=F1UF2UF3U=xβ1xα1+x−γ1−g1x−αxy−exzyβ2yα2+y−γ2−g2y−βxy−δyz1+cy−β3z+εext−τ1zt−τ1+εδyt−τ2zt−τ21+cyt−τ2.
Letℝ+3=0,∞3, since the right-hand side of system (3) is locally Lipschitz on C:ℝ+3+1⟶ℝ3, such that FiUuit=0,U∈ℝ+3≥0, where u1=x, u2=y, and u3=z. According to [30], the solutions of (6) with initial conditions (4) exist uniquely on the interval 0,ξ, where 0<ξ≤∞; therefore, all solutions exist on the first quadrant of the xyz plane.
To prove the boundedness of solutions for system (3), let us first consider the case when the predator is absent, so that(8)dxdt=xβ1xα1+x−γ1−g1x−αxy≡G1x,y,dydt=yβ2yα2+y−γ2−g2y−βxy≡G2x,y,with initial conditions x0>0 and y0>0; we can easily show that G1x,y≥0 for y=0 and x<β1−γ1/g1, such that β1>γ1 and G2x,y≥0 for x=0 and y<β2−γ2/g2, where β2>γ2. Adding the two equations of (8) yields(9)ddtx+y=xβ1xα1+x−γ1−g1x+yβ2yα2+y−γ2−g2y−xyα+β≤xβ1−γ1−g1x+yβ2−γ2−g2y≤β1x+β2y≤κx+y,where κ=minβ1,β2. If we integrate both sides of (9), we get(10)xt+yt≤x0+y0e−κt.
Sincex0+y0>0, the solutions are bounded, which clearly shows that limt⟶∞supxt+yt≤κ.
To extend the analysis to (3), let us consider 0<ϕ1θ+ϕ2θ+ϕ30<M, θ∈−τ,0. Also assume that Wt=εxt−τ1+εδyt−τ2+z and choose 0<ρ<β3. By considering the derivative of W, for t>T+τ for some fixed positive time T, we have(11)dWdt+ρW≤εxt−τ1β1+ρ−xt−τ1+εδyt−τ2β2+ρ−yt−τ2+ρ−β3z.
Sincex and y are nonnegative and bounded by κ,(12)dWdt+ρW≤ε+εδκ+ρ−β3z≤M.
Due to the positivity ofz and since the parametric condition exists for ρ, the differential inequality is bounded above, such that dW/dt≤M−ρW, i.e., there exists N where 0<Wt<N for all t>T, which implies the boundedness of z, such that limt⟶∞supzt≤N.


## 2.1. Positivity and Boundedness of the Solution

The positivity of the solutions indicates the existence of the population, while the boundedness explains the natural control of growth due to the restriction of resources. We arrive at the following lemma.Lemma 1.
Every solution of system (3) corresponding to initial conditions (4) defined on 0,∞ remains positive for all t≥0, which satisfies(5)limt⟶∞supxt+yt≤κ,limt⟶∞supzt≤N,where κ=minβ1,β2 and N>0.Proof.
Model (3) can be represented in a matrix form(6)U˙t=FU,where U=x,y,zT∈ℝ3 and(7)FU=F1UF2UF3U=xβ1xα1+x−γ1−g1x−αxy−exzyβ2yα2+y−γ2−g2y−βxy−δyz1+cy−β3z+εext−τ1zt−τ1+εδyt−τ2zt−τ21+cyt−τ2.
Letℝ+3=0,∞3, since the right-hand side of system (3) is locally Lipschitz on C:ℝ+3+1⟶ℝ3, such that FiUuit=0,U∈ℝ+3≥0, where u1=x, u2=y, and u3=z. According to [30], the solutions of (6) with initial conditions (4) exist uniquely on the interval 0,ξ, where 0<ξ≤∞; therefore, all solutions exist on the first quadrant of the xyz plane.
To prove the boundedness of solutions for system (3), let us first consider the case when the predator is absent, so that(8)dxdt=xβ1xα1+x−γ1−g1x−αxy≡G1x,y,dydt=yβ2yα2+y−γ2−g2y−βxy≡G2x,y,with initial conditions x0>0 and y0>0; we can easily show that G1x,y≥0 for y=0 and x<β1−γ1/g1, such that β1>γ1 and G2x,y≥0 for x=0 and y<β2−γ2/g2, where β2>γ2. Adding the two equations of (8) yields(9)ddtx+y=xβ1xα1+x−γ1−g1x+yβ2yα2+y−γ2−g2y−xyα+β≤xβ1−γ1−g1x+yβ2−γ2−g2y≤β1x+β2y≤κx+y,where κ=minβ1,β2. If we integrate both sides of (9), we get(10)xt+yt≤x0+y0e−κt.
Sincex0+y0>0, the solutions are bounded, which clearly shows that limt⟶∞supxt+yt≤κ.
To extend the analysis to (3), let us consider 0<ϕ1θ+ϕ2θ+ϕ30<M, θ∈−τ,0. Also assume that Wt=εxt−τ1+εδyt−τ2+z and choose 0<ρ<β3. By considering the derivative of W, for t>T+τ for some fixed positive time T, we have(11)dWdt+ρW≤εxt−τ1β1+ρ−xt−τ1+εδyt−τ2β2+ρ−yt−τ2+ρ−β3z.
Sincex and y are nonnegative and bounded by κ,(12)dWdt+ρW≤ε+εδκ+ρ−β3z≤M.
Due to the positivity ofz and since the parametric condition exists for ρ, the differential inequality is bounded above, such that dW/dt≤M−ρW, i.e., there exists N where 0<Wt<N for all t>T, which implies the boundedness of z, such that limt⟶∞supzt≤N.


## 3. Local Stability and Hopf Bifurcation

In this section, we investigate the qualitative behaviour of system (3) by studying the local stability of positive equilibrium points and Hopf bifurcation analysis, which provides a deeper insight into the model to address the behavioral change of solutions as a response to changes in a particular parameter. Since time lags τ1 and τ2 have a significant impact in the complexity and dynamics of the model, we consider them as bifurcation parameters.


### 3.1. Existence of Interior Equilibrium Points

System (3) has some boundary and interior equilibrium points. However, we only focus on the dynamic analysis of interior equilibrium points. In order to obtain the attainable equilibrium points for system (3), the zero growth isoclines of the system are given by xβ1x/α1+x−γ1−g1x−αxy−exz=0, yβ2y/α2+y−γ2−g2y−βxy−δyz/1+cy=0, and −β3z+εexz+εδyz/1+cy=0, in ℝ+3=x,y,z∈ℝ3:x,y,z≥0. Therefore, the equilibria are the points of intersection of these zero growth isoclines regardless of the parameter values.An interior equilibrium pointℰ∗≡x∗,y∗,z∗ exists with β1x∗/α1+x∗−γ1−g1x∗−αy∗−ez∗=0, β2y∗/α2+y∗−γ2−g2y∗−βx∗−δz∗/1+cy∗=0, and −β3+εex∗+εδy∗/1+cy∗=0 such that x∗=1/εeβ3−εδy∗/1+cy∗>0 and z∗=1/eβ1β31+cy∗−εδy∗/1+cy∗εe+β3+g1β3−εδy∗/1+cy∗−γ1−αy∗>0, where y∗ is the root(s) the following equation:(13)Gy=σ1y4+σ2y3+σ3y2+σ4y+σ5=0.The coefficientsσi, i=1,…,5, are defined by(14)σ1=α2c2,σ2=c2α2+α2cg2+ββ2cεe−βδe−δα+cγ2,σ3=β1δ2ε−cδβ1β2εe+β3+α2βδεc−ββ2α2c2−βδεεe−δ2εg1+cβ2+c2β2+cg1δβ3+cδγ1+α2−δα,σ4=ββ2−α2βδεεe+δβ1β2+cδβ1β3α2+β1α2δ2εe+β3−g1α2δ2ε+cα2δβ3g1+cα2δγ1+cα2γ2−β2−β2c+α2g2−δβ2g1−δγ1+γ2,σ5=δβ1β3α2εe+β3−α2β2βεe−α2δβ3g1−α2δγ1+α2γ2.The nature of the roots for (13) is determined by noting the sign of its discriminant [31]. Therefore, a sufficient condition that guarantees that (13) has at least one positive root is σ5<0, which leads to δβ1β3α2/εe+β3+α2γ2<α2β2β/εe+α2δβ3g1+α2δγ1. Thus, system (3) can have at most four interior equilibria in the presence of the Allee effect. However, in the absence of Allee effect, we arrive at the following remark.Remark 1.
In the absence of the Allee effect (α1=α2=0), the interior equilibria for system (3) are reduced to at most three interior equilibria. Consequently, Allee effect can generate or eradicate interior equilibria. It may stabilize or destabilize the system.


#### 3.1.1. Existence of Bistability

The phenomenon of bistability has been recognized experimentally in some biological situations but much more commonly in theoretical models, such as the dynamics of animal populations [32]. The coexistence between two stable attractors can be determined by increasing or decreasing the value of some control parameters. Therefore, the system pursues one branch of equilibrium points when increasing a control parameter until a threshold limit point is reached at which the system jumps to another branch of stable equilibrium points. Bistability occurs when the system can converge to two different equilibrium points, depending on the variation of the initial conditions in the same parametric region. Or the system is able to evolve into either one of two equilibrium points by increasing or decreasing the level of one of the system’s parameters.The underlying model (3) displays bistability of two interior equilibria, which is based on the variation of the coefficient of competition α (see Figure 3). Both equilibria are locally asymptotically stable.Figure 3
Bistability of two interior equilibria for the delayed system (3), with α=0.9 and α=0.5. Both equilibria are locally asymptotically stable; other parameter values are given in (45).


### 3.2. Stability and Bifurcation Analysis of the Interior Equilibriumℰ∗

Now, we study the stability of the interior equilibriumℰ∗≡x∗,y∗,z∗, at which the Jacobian matrix is(15)J=F1F2F3F4F5F6I1e−λτ1I2e−λτ2F7+I3e−λτ1+I4e−λτ2.Here,(16)F1=β1x∗α1+x∗1+α1α1+x∗−2g1x∗−γ1−αy∗−ez∗<0,F2=−αx∗,F3=−ex∗,F4=−βy∗,F5=β2y∗α2+y∗1+α2α2+y∗−2g2y∗−γ2−βx∗−δz∗1+cy∗2<0,F6=−δy∗1+cy∗,F7=−β3,I1=εez∗,I2=εδz∗1+cy∗2,I3=εex∗,I4=εδy∗1+cy∗.The characteristic equation for the interior pointℰ∗≡x∗,y∗,z∗ is then given by(17)Aλ+Bλe−λτ1+Cλe−λτ2=0.Here,(18)Aλ=λ3+R1λ2+R2λ+R3,Bλ=N1λ2+N2λ+N3,Cλ=M1λ2+M2λ+M3,such that(19)R1=−F1−F5−F7,R2=F1F5+F1F7+F5F7−F2F4,R3=F2F4F7−F1F5F7,N1=−I3,N2=F1+F5I3−F3I1,N3=F2F4I3+F3F5I1−F2F6I1−F1F5I3,M1=−I4,M2=F1+F5I4−F6I2,M3=F2F4I4+F1F6I2−F3F4I2−F1F5I4.To gain insight regarding interior equilibriumℰ∗, we discuss the stability of interior steady states and Hopf bifurcation conditions of the threshold parameters τ1 and τ2 by considering the following different cases.Case 1.
Whenτ1=τ2=0, equation (17) becomes(20)λ3+R1+N1+M1λ2+R2+N2+M2λ+R3+N3+M3=0.
Therefore, the interior equilibriumℰ∗ is locally asymptotically stable if(i)
H1R1+N1+M1>0,R3+N3+M3>0andR1+N1+M1R2+N2+M2>R3+N3+M3 holds
Thus, based on Routh–Hurwitz Criteria, all the roots of (20) have negative real parts.Case 2.
Forτ1=0,τ2>0, then equation (17) becomes(21)λ3+R1+N1λ2+R2+N2λ+R3+N3+M1λ2+M2λ+M3e−λτ2=0.
We assume for some values of (τ2>0), there exists a real number ω such that λ=iω is a root of (21); then, we get(22)−R1+M1ω2+R3+N3=M1ω2−M3cosωτ2−M2ωsinωτ2,−ω3+R2+N2ω=M3−M1ω2sinωτ2−M2ωcosωτ2.
Squaring and adding both of the equations yield(23)ω6+a1ω4+a2ω2+a3=0,where(24)a1=R1+M12−2R2+N2−M12,a2=R2+N22−2R1+M1R3+N3+2M1M3−M22,a3=R3+N32−M32.
By Descartes’ rule of signs, equation (22) has at least one positive root ω1 if(i)
H2R12+2R1M1>2R2+N2andR3+N32<M32 holds
Eliminatingsinω1τ2 from (22) yields(25)τ2,j=1w1arccosR3+N3−R1+N1w12M1w12−M3+M2R2+N2w12−M2w14M3−M1w122−M2w12+2jπw1,where j=0,1,2,….
By differentiating (21) with respect to τ2 such that ω=ω1 and τ2=τ2,j, the transversality condition can be obtained in this form:(26)ℛedλdτ2−1=A1A4−A2A3A2A4.
Here,(27)A1=R2+N2−3ω12R2+N2ω12−ω14+2R1+N1ω1R3+N3ω1−R1+N1ω13,A2=ω14−R2+N2ω122+R3+N3ω1−R1+N1ω132,A3=M22ω12+2M1ω13−M3ω1M1ω1,A4=M2ω122+M3ω1−M1ω132.
Then, a Hopf bifurcation occurs forτ2 if ℛedλ/dτ2−1>0; i.e., A1A4>A2A3. We arrive at the following theorem.Theorem 1.
LetH1 and H2 hold, where τ1=0; then, there exists τ2>0 such that ℰ∗ remains stable for τ2<τ2′ and unstable for τ2>τ2′, where τ2′=minτ2,j defined by (25). Moreover, system (3) undergoes a Hopf bifurcation at ℰ∗ when τ2=τ2′.Case 3.
Whenτ2=0andτ1>0, in the same manner of the pervious case, we arrive at the following theorem.Theorem 2.
For system (3), with τ2=0, there exists a positive number τ1, such that the equilibrium point ℰ∗ is locally asymptotically stable for τ1<τ1′ and unstable for τ1>τ1′, where τ1′=minτ1,j. Furthermore, Hopf bifurcation occurs at τ1=τ1′:(28)τ1,j=1w0arccosR3+M3−R1+M1w22N1w22−N3−N2R2+M2w22+N2w24N1w22−N32+N2w22+2jπw2,where j=0,1,2,….Case 4.
Whenτ1>0 and τ2>0, we assume that τ1 as a variable parameter and τ2 as fixed on its stable interval. Let λ=iw as a root of (17); separating real and imaginary parts implies(29)−w3+R2w+M1w2−M3sinwτ2+M2wcoswτ2=N3−N1w2sinwτ1−N2wcoswτ1,(30)−R1w2+R3+M3−M1w2coswτ2+M2wsinwτ2=N1w2−N3coswτ1−N2wsinwτ1.
Thus, eliminating the trigonometric functions (sinωτ1 and cosωτ1) from (29) and (30) yields(31)w6+ξ4w5+ξ3w4+ξ2w3+ξ1w2+ξ0=0,where(32)ξ4=−2M1sinwτ2,ξ3=R1+M12−2R2−N12−2M2coswτ2,ξ2=2M1R2+M3sinwτ2−2M3R1coswτ2,ξ1=−2M3R2sinwτ2,ξ0=R32+M32−N32+2M3R3+R1M1coswτ2.
Equation (31) is a preternatural equation in a complicated form; it is quite difficult to predict the nature of its roots. Thus, by applying Descartes rule of signs, we can say that (31) has at least one positive root ω0 if(i)
H3ξ4>0sinceM1<0andξ0<0; therefore, we have(33)τ1,j=1w0arccosAD+CBA2+C2+2jπw0,j=0,1,2,….
Here,(34)A=N1w02−N3,B=−w03+R2w0+M3−M1w02sinw0τ2+cosw0τ2,C=N2w0,D=−R1w02+R3+M1w02−M3cosw0τ2+M2w0sinw0τ2.
To study the Hopf bifurcation analysis, we fixτ2 in its stable interval and differentiate equations (29) and (30) with respect to τ1. Then, substitute τ1=τ1,0 and w=w0, we have(35)Q2dℜλdτ1τ1=τ1,0+Q1dwdτ1τ1=τ1,0=Q3,−Q1dℜλdτ1τ1=τ1,0+Q2dwdτ1τ1=τ1,0=Q4where(36)Q1=−3w02+R2+2N1w0−N2w0τ1,0sinw0τ1,0+N2+N1τ1w02−N3τ1,0cosw0τ1,0+2w0M1−M2τ2w0sinw0τ2+M1τ2w02−M3τ2+M2cosw0τ2,Q2=−2R1w0+N1w02τ1,0−N3τ1,0+N2sinw0τ1,0+N2w0τ1−2N1w0cosw0τ1,0+M2+M1w02τ2−M3τ2sinw0τ2+M2w0τ2−2M1w0cosw0τ2,Q3=N2w02sinw0τ1,0+N3w0−N1w03cosw0τ1,0,Q4=N2w02cosw0τ1,0+N1w03−N3w0sinw0τ1,0.
From (35), we get(37)dℜλdτ1τ1=τ1,0=Q2Q3−Q1Q4Q22+Q12.
AsQ2Q3>Q1Q4, then Hopf bifurcation occurs for τ1=τ1,0.
Therefore, we arrive at the following theorem.Theorem 3.
Ifℰ∗ exists, such that H1 and H3 hold, with τ2∈0,τ2′, then there exists a positive threshold parameter τ1∗ such that the interior equilibrium ℰ∗ is locally asymptotically stable for τ1<τ1∗ and unstable τ1>τ1∗, where τ1∗=minτ1,j as in (38). Additionally, system (3) undergoes Hopf bifurcation at ℰ∗ when τ1=τ1∗.Remark 2.
Similarly, forτ1∈0,τ1′, there exists a threshold parameter τ2∗ such that the interior equilibrium ℰ∗ is locally asymptotically stable for τ2<τ2∗ and unstable τ2>τ2∗. Also, Hopf bifurcation occurs for system (3) as τ2=τ2∗, where τ2∗=minτ2,j is given by(38)τ2,j=1w3arccosA1D1+C1B1A12+C12+2jπw3,j=0,1,2,…,with(39)A1=M1w32−M3,B1=w33−R2w3+N3−N1w32sinw3τ1−N2w3cosw3τ1,C1=M2w3,D1=−R1w32+R3+cosw3τ1+N2w3sinw3τ1.
The proofs are obtained in the same manner of the above analysis.


## 3.1. Existence of Interior Equilibrium Points

System (3) has some boundary and interior equilibrium points. However, we only focus on the dynamic analysis of interior equilibrium points. In order to obtain the attainable equilibrium points for system (3), the zero growth isoclines of the system are given by xβ1x/α1+x−γ1−g1x−αxy−exz=0, yβ2y/α2+y−γ2−g2y−βxy−δyz/1+cy=0, and −β3z+εexz+εδyz/1+cy=0, in ℝ+3=x,y,z∈ℝ3:x,y,z≥0. Therefore, the equilibria are the points of intersection of these zero growth isoclines regardless of the parameter values.An interior equilibrium pointℰ∗≡x∗,y∗,z∗ exists with β1x∗/α1+x∗−γ1−g1x∗−αy∗−ez∗=0, β2y∗/α2+y∗−γ2−g2y∗−βx∗−δz∗/1+cy∗=0, and −β3+εex∗+εδy∗/1+cy∗=0 such that x∗=1/εeβ3−εδy∗/1+cy∗>0 and z∗=1/eβ1β31+cy∗−εδy∗/1+cy∗εe+β3+g1β3−εδy∗/1+cy∗−γ1−αy∗>0, where y∗ is the root(s) the following equation:(13)Gy=σ1y4+σ2y3+σ3y2+σ4y+σ5=0.The coefficientsσi, i=1,…,5, are defined by(14)σ1=α2c2,σ2=c2α2+α2cg2+ββ2cεe−βδe−δα+cγ2,σ3=β1δ2ε−cδβ1β2εe+β3+α2βδεc−ββ2α2c2−βδεεe−δ2εg1+cβ2+c2β2+cg1δβ3+cδγ1+α2−δα,σ4=ββ2−α2βδεεe+δβ1β2+cδβ1β3α2+β1α2δ2εe+β3−g1α2δ2ε+cα2δβ3g1+cα2δγ1+cα2γ2−β2−β2c+α2g2−δβ2g1−δγ1+γ2,σ5=δβ1β3α2εe+β3−α2β2βεe−α2δβ3g1−α2δγ1+α2γ2.The nature of the roots for (13) is determined by noting the sign of its discriminant [31]. Therefore, a sufficient condition that guarantees that (13) has at least one positive root is σ5<0, which leads to δβ1β3α2/εe+β3+α2γ2<α2β2β/εe+α2δβ3g1+α2δγ1. Thus, system (3) can have at most four interior equilibria in the presence of the Allee effect. However, in the absence of Allee effect, we arrive at the following remark.Remark 1.
In the absence of the Allee effect (α1=α2=0), the interior equilibria for system (3) are reduced to at most three interior equilibria. Consequently, Allee effect can generate or eradicate interior equilibria. It may stabilize or destabilize the system.


### 3.1.1. Existence of Bistability

The phenomenon of bistability has been recognized experimentally in some biological situations but much more commonly in theoretical models, such as the dynamics of animal populations [32]. The coexistence between two stable attractors can be determined by increasing or decreasing the value of some control parameters. Therefore, the system pursues one branch of equilibrium points when increasing a control parameter until a threshold limit point is reached at which the system jumps to another branch of stable equilibrium points. Bistability occurs when the system can converge to two different equilibrium points, depending on the variation of the initial conditions in the same parametric region. Or the system is able to evolve into either one of two equilibrium points by increasing or decreasing the level of one of the system’s parameters.The underlying model (3) displays bistability of two interior equilibria, which is based on the variation of the coefficient of competition α (see Figure 3). Both equilibria are locally asymptotically stable.Figure 3
Bistability of two interior equilibria for the delayed system (3), with α=0.9 and α=0.5. Both equilibria are locally asymptotically stable; other parameter values are given in (45).


## 3.1.1. Existence of Bistability

The phenomenon of bistability has been recognized experimentally in some biological situations but much more commonly in theoretical models, such as the dynamics of animal populations [32]. The coexistence between two stable attractors can be determined by increasing or decreasing the value of some control parameters. Therefore, the system pursues one branch of equilibrium points when increasing a control parameter until a threshold limit point is reached at which the system jumps to another branch of stable equilibrium points. Bistability occurs when the system can converge to two different equilibrium points, depending on the variation of the initial conditions in the same parametric region. Or the system is able to evolve into either one of two equilibrium points by increasing or decreasing the level of one of the system’s parameters.The underlying model (3) displays bistability of two interior equilibria, which is based on the variation of the coefficient of competition α (see Figure 3). Both equilibria are locally asymptotically stable.Figure 3
Bistability of two interior equilibria for the delayed system (3), with α=0.9 and α=0.5. Both equilibria are locally asymptotically stable; other parameter values are given in (45).


## 3.2. Stability and Bifurcation Analysis of the Interior Equilibriumℰ∗

Now, we study the stability of the interior equilibriumℰ∗≡x∗,y∗,z∗, at which the Jacobian matrix is(15)J=F1F2F3F4F5F6I1e−λτ1I2e−λτ2F7+I3e−λτ1+I4e−λτ2.Here,(16)F1=β1x∗α1+x∗1+α1α1+x∗−2g1x∗−γ1−αy∗−ez∗<0,F2=−αx∗,F3=−ex∗,F4=−βy∗,F5=β2y∗α2+y∗1+α2α2+y∗−2g2y∗−γ2−βx∗−δz∗1+cy∗2<0,F6=−δy∗1+cy∗,F7=−β3,I1=εez∗,I2=εδz∗1+cy∗2,I3=εex∗,I4=εδy∗1+cy∗.The characteristic equation for the interior pointℰ∗≡x∗,y∗,z∗ is then given by(17)Aλ+Bλe−λτ1+Cλe−λτ2=0.Here,(18)Aλ=λ3+R1λ2+R2λ+R3,Bλ=N1λ2+N2λ+N3,Cλ=M1λ2+M2λ+M3,such that(19)R1=−F1−F5−F7,R2=F1F5+F1F7+F5F7−F2F4,R3=F2F4F7−F1F5F7,N1=−I3,N2=F1+F5I3−F3I1,N3=F2F4I3+F3F5I1−F2F6I1−F1F5I3,M1=−I4,M2=F1+F5I4−F6I2,M3=F2F4I4+F1F6I2−F3F4I2−F1F5I4.To gain insight regarding interior equilibriumℰ∗, we discuss the stability of interior steady states and Hopf bifurcation conditions of the threshold parameters τ1 and τ2 by considering the following different cases.Case 1.
Whenτ1=τ2=0, equation (17) becomes(20)λ3+R1+N1+M1λ2+R2+N2+M2λ+R3+N3+M3=0.
Therefore, the interior equilibriumℰ∗ is locally asymptotically stable if(i)
H1R1+N1+M1>0,R3+N3+M3>0andR1+N1+M1R2+N2+M2>R3+N3+M3 holds
Thus, based on Routh–Hurwitz Criteria, all the roots of (20) have negative real parts.Case 2.
Forτ1=0,τ2>0, then equation (17) becomes(21)λ3+R1+N1λ2+R2+N2λ+R3+N3+M1λ2+M2λ+M3e−λτ2=0.
We assume for some values of (τ2>0), there exists a real number ω such that λ=iω is a root of (21); then, we get(22)−R1+M1ω2+R3+N3=M1ω2−M3cosωτ2−M2ωsinωτ2,−ω3+R2+N2ω=M3−M1ω2sinωτ2−M2ωcosωτ2.
Squaring and adding both of the equations yield(23)ω6+a1ω4+a2ω2+a3=0,where(24)a1=R1+M12−2R2+N2−M12,a2=R2+N22−2R1+M1R3+N3+2M1M3−M22,a3=R3+N32−M32.
By Descartes’ rule of signs, equation (22) has at least one positive root ω1 if(i)
H2R12+2R1M1>2R2+N2andR3+N32<M32 holds
Eliminatingsinω1τ2 from (22) yields(25)τ2,j=1w1arccosR3+N3−R1+N1w12M1w12−M3+M2R2+N2w12−M2w14M3−M1w122−M2w12+2jπw1,where j=0,1,2,….
By differentiating (21) with respect to τ2 such that ω=ω1 and τ2=τ2,j, the transversality condition can be obtained in this form:(26)ℛedλdτ2−1=A1A4−A2A3A2A4.
Here,(27)A1=R2+N2−3ω12R2+N2ω12−ω14+2R1+N1ω1R3+N3ω1−R1+N1ω13,A2=ω14−R2+N2ω122+R3+N3ω1−R1+N1ω132,A3=M22ω12+2M1ω13−M3ω1M1ω1,A4=M2ω122+M3ω1−M1ω132.
Then, a Hopf bifurcation occurs forτ2 if ℛedλ/dτ2−1>0; i.e., A1A4>A2A3. We arrive at the following theorem.Theorem 1.
LetH1 and H2 hold, where τ1=0; then, there exists τ2>0 such that ℰ∗ remains stable for τ2<τ2′ and unstable for τ2>τ2′, where τ2′=minτ2,j defined by (25). Moreover, system (3) undergoes a Hopf bifurcation at ℰ∗ when τ2=τ2′.Case 3.
Whenτ2=0andτ1>0, in the same manner of the pervious case, we arrive at the following theorem.Theorem 2.
For system (3), with τ2=0, there exists a positive number τ1, such that the equilibrium point ℰ∗ is locally asymptotically stable for τ1<τ1′ and unstable for τ1>τ1′, where τ1′=minτ1,j. Furthermore, Hopf bifurcation occurs at τ1=τ1′:(28)τ1,j=1w0arccosR3+M3−R1+M1w22N1w22−N3−N2R2+M2w22+N2w24N1w22−N32+N2w22+2jπw2,where j=0,1,2,….Case 4.
Whenτ1>0 and τ2>0, we assume that τ1 as a variable parameter and τ2 as fixed on its stable interval. Let λ=iw as a root of (17); separating real and imaginary parts implies(29)−w3+R2w+M1w2−M3sinwτ2+M2wcoswτ2=N3−N1w2sinwτ1−N2wcoswτ1,(30)−R1w2+R3+M3−M1w2coswτ2+M2wsinwτ2=N1w2−N3coswτ1−N2wsinwτ1.
Thus, eliminating the trigonometric functions (sinωτ1 and cosωτ1) from (29) and (30) yields(31)w6+ξ4w5+ξ3w4+ξ2w3+ξ1w2+ξ0=0,where(32)ξ4=−2M1sinwτ2,ξ3=R1+M12−2R2−N12−2M2coswτ2,ξ2=2M1R2+M3sinwτ2−2M3R1coswτ2,ξ1=−2M3R2sinwτ2,ξ0=R32+M32−N32+2M3R3+R1M1coswτ2.
Equation (31) is a preternatural equation in a complicated form; it is quite difficult to predict the nature of its roots. Thus, by applying Descartes rule of signs, we can say that (31) has at least one positive root ω0 if(i)
H3ξ4>0sinceM1<0andξ0<0; therefore, we have(33)τ1,j=1w0arccosAD+CBA2+C2+2jπw0,j=0,1,2,….
Here,(34)A=N1w02−N3,B=−w03+R2w0+M3−M1w02sinw0τ2+cosw0τ2,C=N2w0,D=−R1w02+R3+M1w02−M3cosw0τ2+M2w0sinw0τ2.
To study the Hopf bifurcation analysis, we fixτ2 in its stable interval and differentiate equations (29) and (30) with respect to τ1. Then, substitute τ1=τ1,0 and w=w0, we have(35)Q2dℜλdτ1τ1=τ1,0+Q1dwdτ1τ1=τ1,0=Q3,−Q1dℜλdτ1τ1=τ1,0+Q2dwdτ1τ1=τ1,0=Q4where(36)Q1=−3w02+R2+2N1w0−N2w0τ1,0sinw0τ1,0+N2+N1τ1w02−N3τ1,0cosw0τ1,0+2w0M1−M2τ2w0sinw0τ2+M1τ2w02−M3τ2+M2cosw0τ2,Q2=−2R1w0+N1w02τ1,0−N3τ1,0+N2sinw0τ1,0+N2w0τ1−2N1w0cosw0τ1,0+M2+M1w02τ2−M3τ2sinw0τ2+M2w0τ2−2M1w0cosw0τ2,Q3=N2w02sinw0τ1,0+N3w0−N1w03cosw0τ1,0,Q4=N2w02cosw0τ1,0+N1w03−N3w0sinw0τ1,0.
From (35), we get(37)dℜλdτ1τ1=τ1,0=Q2Q3−Q1Q4Q22+Q12.
AsQ2Q3>Q1Q4, then Hopf bifurcation occurs for τ1=τ1,0.
Therefore, we arrive at the following theorem.Theorem 3.
Ifℰ∗ exists, such that H1 and H3 hold, with τ2∈0,τ2′, then there exists a positive threshold parameter τ1∗ such that the interior equilibrium ℰ∗ is locally asymptotically stable for τ1<τ1∗ and unstable τ1>τ1∗, where τ1∗=minτ1,j as in (38). Additionally, system (3) undergoes Hopf bifurcation at ℰ∗ when τ1=τ1∗.Remark 2.
Similarly, forτ1∈0,τ1′, there exists a threshold parameter τ2∗ such that the interior equilibrium ℰ∗ is locally asymptotically stable for τ2<τ2∗ and unstable τ2>τ2∗. Also, Hopf bifurcation occurs for system (3) as τ2=τ2∗, where τ2∗=minτ2,j is given by(38)τ2,j=1w3arccosA1D1+C1B1A12+C12+2jπw3,j=0,1,2,…,with(39)A1=M1w32−M3,B1=w33−R2w3+N3−N1w32sinw3τ1−N2w3cosw3τ1,C1=M2w3,D1=−R1w32+R3+cosw3τ1+N2w3sinw3τ1.
The proofs are obtained in the same manner of the above analysis.


## 4. Global Stability of Interior Steady Stateℰ∗

In this section, we study the global stability of system (3) around interior steady state ℰ∗≡x∗,y∗,z∗.Theorem 4.
Ifβ1α1<g1α1+x∗α1+x and β2α21+cy∗1+cy+δcz∗α2+y∗α2+y<g2α2+y∗α2+y1+cy∗1+cy, then system (3) is globally asymptotically stable at the interior equilibrium point ℰ∗.Proof.
We suggest the Lyapunov function atℰ∗≡x∗,y∗,z∗ of the form:(40)Vt=ϱ1xt−x∗−x∗lnxtx∗+ϱ2yt−y∗−y∗lnyty∗+ϱ3zt−z∗−z∗lnztz∗,where ϱ1,ϱ2,andϱ3 are nonnegative constants. Take derivative V with respect to time t, yielding to(41)V˙t=ϱ1x−x∗xx˙t+ϱ2y−y∗yy˙t+ϱ3z−z∗zz˙t=ϱ1x−x∗β1xα1+x−γ1−g1x−αy−ez+ϱ2y−y∗β2yα2+y−γ2−g2y−βx−δyzy1+cy+ϱ3z−z∗−β3+εext−τ1zt−τ1z+εδyt−τ2zt−τ2z1+cyt−τ2≤ϱ1x−x∗β1xα1+x−β1x∗α1+x∗−g1x−x∗−αy−y∗−ez−z∗+ϱ2y−y∗β2yα2+y−β2y∗α2+y∗−g2y−y∗−βx−x∗+δy∗z∗y∗1+cy∗−δyzy1+cy+ϱ3z−z∗εext−τ1zt−τ1z+εδyt−τ2zt−τ2z1+cyt−τ2−εex∗+εδy∗1+cy∗≤−ϱ1g1x−x∗2−ϱ2g2y−y∗2−ϱ1α+ϱ2βx−x∗y−y∗+εeϱ3−eϱ1x−x∗z−z∗+ϱ1x−x∗β1xα1+x−β1x∗α1+x∗+ϱ2y−y∗β2yα2+y−β2y∗α2+y∗+ϱ2y−y∗δy∗z∗y∗1+cy∗−δyzy1+cy+ϱ3z−z∗ϵδy1+cy−ϵδy∗1+cy∗≤−ϱ1g1x−x∗2−ϱ2g2y−y∗2−ϱ1α+ϱ2βx−x∗y−y∗+εeϱ3−eϱ1x−x∗z−z∗+β1ϱ1x−x∗2α1α1+x∗α1+x+β2ϱ2y−y∗2α2α2+y∗α2+y+δϱ2y−y∗−z−z∗1+cy+cz∗y−y∗1+cy∗1+cy+εδϱ3y−y∗z−z∗11+cy−cy1+cy∗1+cy.
Thus, based on the assumptions:β1α1<g1α1+x∗α1+x, β2α21+cy∗1+cy+δcz∗α2+y∗α2+y<g2α2+y∗α2+y1+cy∗1+cy, and εϱ3<maxϱ1,ϱ2, we can get(42)V˙t≤ϱ1α1β1α1+x∗α1+x−ϱ1g1x−x∗2+εδϱ3−δϱ21+cyy−y∗z−z∗+δϱ2cz∗1+cy∗1+cy+ϱ2α2β2α2+y∗α2+y−ϱ2g2y−y∗2+εeϱ3−ϱ1ex−x∗z−z∗−ϱ1α+ϱ2βx−x∗y−y∗−εδϱ3cy1+cy∗1+cyz−z∗y−y∗≤0.
Hence, the proof is complete.


## 5. Sensitivity Analysis to Severity of Allee Effect

Here, we study the sensitivity of model solution of (3), with respect to the parameters α1 and α2 (strength Allee effect). It is quite common for a model to exhibit high sensitivity to small variations in some parameters, while showing robustness to variation in other parameters. There are different ways to find the sensitivity functions of DDEs [33]. Nevertheless, we utilize the so-called “direct approach” to find sensitivity functions of model (3). The sensitivity functions with respect to Allee parameters αii=1,2 are denoted by(43)Sxαit:=∂∂αixt,Syαit:=∂∂αiyt,Szαit:=∂∂αizt.Hence, sensitivity functions due to small perturbations in Allee parameterα1 are given by a system of DDEs:(44)Sxα1′t=Sxα1tβ1xtα1+xt−γ1−2g1xt−αyt−ezt−αSyα1txt−eSzα1txt+β1xtα1Sxα1t−xtα1+xt2,Syα1′t=Syα1tβ2ytα2+yt−γ2−2g2yt−βxt+yt−βSxα1t+α2β2Syα1tα2+yt2−δSyα1tzt1+cyt2+Szα1tyt1+cyt,Szα1′t=−β3Szα1t+εeSxα1t−τ1zt−τ1+Szα1t−τ1xt−τ1+εδSyα1t−τ2zt−τ21+cyt−τ22+Szα1t−τ2yt−τ21+cyt−τ2.To estimate the sensitivity functionsSxα1t, Syα1t, and Szα1t, we have to solve the system of sensitivity equations (44) together with the original system (3).Similarly, the sensitivity functions due to small changes in Allee coefficientα2 satisfy the system of DDEs:(45)Sxα2′t=Sxα2tβ1xtα1+xt−γ1−2g1xt−αyt−ezt+β1xtα1Sxα2tα1+xt2−αSyα2txt−eSzα2txt,Syα2′t=Syα2tβ2ytα2+yt−γ2−2g2yt−βxt−βSxα2tyt+β2ytα2Syα2t−ytα2+yt2−δSyα2tzt1+cyt2+Szα2tyt1+cyt,Szα2′t=−β3Szα2t+εeSxα2t−τ1zt−τ1+Szα2t−τ1xt−τ1+εδSyα2t−τ2zt−τ21+cyt−τ22+Szα2t−τ2yt−τ21+cyt−τ2.We then solve (45) along with (3) to evaluate Sxα2t, Syα2t, and Szα3t (see Figure 4).Figure 4
Sensitivity functions of model solution of system (3) with respect to Allee parameters α1 and α2. The top banners show the sensitivity functions for xt, yt and zt with respect to small changes in Allee parameter α1. However, the bottom banners display the sensitivity with respect to α2. They show that the model is very sensitive to the small perturbations of Allee parameters in early time intervals and the sensitivity decreases by time. The two parameters α1 and α2 are significant in the model and cause high impact in early stages of interactions.


## 6. Numerical Simulations

Some numerical simulations of system (3) are carried out, using Matlab package DDE23, to confirm our theoretical results. We first investigate the behavior of the model around ℰ∗ with parameter values:(46)α=0.9,α1=0.001,α2=0.001,β=1.35,γ2=1,γ1=1,β1=2,β2=2,β3=1,ε=0.5,e=5,δ=1.Figure5 shows the numerical simulations of the delayed system (3) around the steady state ℰ∗. The interior steady state ℰ∗ is asymptotically stable when τ1<τ1∗ and τ2∈0,τ2∗. The model undergoes a Hopf bifurcation when τ1=τ1∗=4.34 and τ2<τ2∗=5.33. Figure 6 displays the Hopf bifurcation diagrams of τ1 and τ2 which are obtained numerically by maximum and minimum amplitudes of zt. The left banner displays the threshold parameter τ1∗=4.34 with τ2<τ2∗, while right banner shows that the threshold parameter τ2∗=5.54 with τ1<τ1∗.Figure 5
Numerical simulations of system (3) around the steady state ℰ∗. Top banners show that ℰ∗ is asymptotically stable when τ1=5.54<τ1∗ and τ2∈0,τ2∗. Bottom banners display a Hopf bifurcation when τ1=τ1∗=4.34 and τ2<τ2∗=5.34; the other parameter values are given in (45).
(a)
(b)
(c)
(d)Figure 6
The Hopf bifurcation diagrams ofτ1 and τ2 which are obtained numerically by maximum and minimum amplitudes of zt. (a) The threshold parameter τ1∗=4.34 with τ2<τ2∗. (b) The threshold parameter τ2∗=5.54 with τ1<τ1∗.
(a)
(b)Figure5 displays a bistability of two interior equilibrium points, for the DDEs model (3), when parameter α varies from α=0.5 to α=0.9. If the interior equilibria exist, any trajectory starting from the interior of ℛ+3 converges to one of the interior equilibria.Figure7 shows the sensitivity of the dynamics of system (3) due to small changes in the severity of Allee effects α1 and α2. The left banners show the numerical simulations with different values of α1 (0.001≤α1≤0.02) and fixed value of α2=0.001, while right banners display the simulations with different values of α2 (0.01≤α2≤0.02) and fixed value of α1=0.01. The phase portrait gets stretched over time as α1 is reduced, while low values of α2 increases the oscillations over time. The presence of Allee effect in the model enriches the dynamics of the system, while Figure 4 exhibits the absolute values of sensitivity functions: ∂xt/∂α1,2, ∂yt/∂α1,2, and ∂zt/∂α1,2 to evaluate the sensitivity of the state variables due to a small perturbation in α1 and α2. The oscillation behaviour indicates that the species population is very sensitive to small changes in the parameter. It is clear that α1 and α2 are important in the model and have a significant impact in the dynamics, specially in the early stages of time. However, the sensitivity to these parameters decreases with time.Figure 7
The sensitivity of the dynamics of system (3) due to small changes in the severity of Allee effect α1 and α2. The left banners show the numerical simulations with different values of α1(0.001≤α1≤0.02) and fixed value of α2=0.001. The right banners display the simulations with different values of α2(0.01≤α2≤0.02) and fixed value of α1=0.01. The phase portrait gets stretched over time as α1 reduced, while low values of α2 increase the oscillations over time. The presence of Allee effect in the model enriches the dynamics of the system.


## 7. Conclusion

In this paper, we established the two-prey one-predator mode with time delays and a weak Allee effect in the preys’ growth functions, where there is a direct competition between prey populations. Although the model is simple, the system exhibits rich dynamic behaviour such as bistability of equilibria, Hopf bifurcation, and period doubling chaos. Nonnegativity and boundedness of the solutions have been investigated. Some new sufficient conditions for local and global asymptotic stability of interior steady states have been deduced. In addition, Hopf bifurcation with respect to time delay threshold parametersτ1∗ and τ2∗ has been studied. The model undergoes a Hopf bifurcation when time delays pass through its critical values. We also investigated the sensitivity of model solutions to small perturbations in the severity Allee effects α1 and α2. The obtained results confirm that Allee effect has a significant impact in the dynamics in the early stages of interaction. It has been seen by the numerical simulations that time delay and Allee effects play an important role in the dynamics of prey-predator systems. Introduction of time delay and Allee effects, in the model, improves the stability results, enriches the dynamics of the system, keeps the population densities in balance, and makes the model closer to reality.As part of future work, more sophisticated models with harvesting terms, control variables, and effect of environmental noise can be studied. Fractional derivatives, instead of integer-order derivatives in the same or similar model, to consider the long-term memory effect, with a saturating functional response, will also be our future goal.



---

*Source: 7306412-2020-01-31.xml*