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

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. )e 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.


Introduction
e 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][2][3][4][5]).
is 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 that N(t) is the size of prey population and P(t) is the size of the predator population at time t, then the Lotka-Volterra model is given by the following equations: with N(0) > 0 and P(0) > 0. Here, β 1 is the per capita maximum filtering rate and c 1 is the death rate of the prey N(t), while the parameter g 1 denotes the strength of intraspecific competition. e predator death rate and predation rate are, respectively, denoted by c 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. e 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. ese responses are used to model the phenomena of predation, which captures the usual properties, for instance, positivity and increase (see also [10][11][12][13]). e 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]).

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) � rN 2 (1 − (N/K)) for a weak Allee effect and (dN/dt) � rN(1 − (N/K))((N/A) − 1) for a strong Allee effect. Figure 1 shows a per capita growth rate (1/N)(dN/dt) of the population with strong and weak Allee effects. e straight line shows the logistic growth, and the red curve displays a weak Allee effect, while the blue curve shows a strong Allee effect. e 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/N)(dN/dt) < 0 (the species will die out) and for N > A, (1/N)(dN/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.
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 form b(N) ≡ (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) e 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 b(0) � 0 and b ′ (N) > 0 if N ∈ [0, ∞), i.e., Allee effect decreases as density increases, and lim N⟶∞ b(N) � 1 means that Allee effect disappears at high densities. erefore, the term b(N) is considered as a weak Allee effect function in a rectangular hyperbola form, known as Michaelis-Menten-like function [21].

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][22][23][24][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. e rest of this paper is organized as follows: Model formulation is presented in Section 2. 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 2 Complexity concludes the study with a summary of the reported findings and future recommendations.

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][27][28][29]). Sen et al. [26] discussed the Allee effect on two-preys' growth function, where the predator is generalized. ey 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. e 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. eir 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 onepredator system by introducing constant periodic releases of predators through periodically spraying a pesticide on the prey. ey were then able to show that the system remains permanent under certain conditions. Herein, we generalize model (2) to a multispecies preypredator system (two-preys one-predator). e model consists of two teams of preys with densities x(t) and y(t), interacting with one team of predator with density z(t). 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. e model takes the general form: with initial conditions Here, ϕ i (θ) (i � 1, 2, 3) are smooth initial functions. e coefficients α and β represent the coefficients of competition of two preys x and y, in the absence of predator. e description of rest of model parameters along with their symbols is presented in Table 1. 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.
It 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. erefore, we incorporated two discrete time delays τ 1 and τ 2 in the reaction response functionals in the predator growth to represent the reaction time. e 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) δy(t)z(t)/(1 + cy(t)), response indicates that it is a hard-to-capture prey compared to the first species (see Figure 2).
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).

Positivity and Boundedness of the Solution.
e 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.

Proof. Model (3) can be represented in a matrix form
where U � (x, y, z) T ∈ R 3 and Let 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 with initial conditions x(0) > 0 and y(0) > 0; we can easily Adding the two equations of (8) yields where κ � min β 1 , β 2 . If we integrate both sides of (9), we get Parameter Description α 1 , α 2 Strength of Allee effect β 1 , β 2 Per capita maximum filtering rate of population g 1 , g 2 Strength of intracompetition c 1 , c 2 Death rate for preys α, β Coefficient of competition e, δ Decrease rate of x(t) and y(t) due to predation by z(t) β 3 Predator death rate c Magnitude of interference between the second type of prey ∈ An equal transformation rate of predator to preys x(t) and y(t) βxy αxy  Complexity Since (x(0) + y(0)) > 0, the solutions are bounded, which clearly shows that lim t⟶∞ sup(x(t) + y(t)) ≤ κ.

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) 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

Existence of Interior Equilibrium Points. System
erefore, the equilibria are the points of intersection of these zero growth isoclines regardless of the parameter values.

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.

Existence of Bistability.
e 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]. e coexistence between two stable attractors can be determined by increasing or decreasing the value of some control parameters. erefore, 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. e 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.

Stability and Bifurcation Analysis of the Interior Equi
librium E * . Now, we study the stability of the interior equilibrium E * ≡ (x * , y * , z * ), at which the Jacobian matrix is Here, e characteristic equation for the interior point E * ≡ (x * , y * , z * ) is then given by Here, such that To gain insight regarding interior equilibrium E * , 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 erefore, the interior equilibrium E * is locally asymptotically stable if us, based on Routh-Hurwitz Criteria, all the roots of (20) have negative real parts.

Moreover, system (3) undergoes a Hopf bifurcation at
Case 3. When τ 2 � 0 and τ 1 > 0, in the same manner of the pervious case, we arrive at the following theorem.

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 E * with parameter values: δ � 1.  Figure 5 shows the numerical simulations of the delayed system (3) around the steady state E * . e interior steady state E * is asymptotically stable when τ 1 < τ * 1 and τ 2 ∈ (0, τ * 2 ). e 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 z(t).
e 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 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 R 3 + converges to one of the interior equilibria.    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.

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. e 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 . e 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.

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.