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

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., [15]). 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 that is the size of prey population and is the size of the predator population at time t, then the Lotka–Volterra model is given by the following equations:with and . Here, is the per capita maximum filtering rate and is the death rate of the prey , while the parameter 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 [1013]).

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: for a weak Allee effect and for a strong Allee effect.

Figure 1 shows a per capita growth rate 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 , (the species will die out) and for , , 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 in the prey growth function of model (1), which is considered as the probability of finding a mate (see [19]), so that

The parameter is the strength of Allee effect, , where R is the average area that can be searched by an individual [20]. We may notice that and if , i.e., Allee effect decreases as density increases, and means that Allee effect disappears at high densities. Therefore, the term 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 [2125]).

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 and 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 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 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., [2629]). 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 and , interacting with one team of predator with density . 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:with initial conditions

Here, () 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.

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. Therefore, we incorporated two discrete time delays and 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) , 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).

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 remains positive for all , which satisfieswhere and .

Proof. Model (3) can be represented in a matrix formwhere andLet , since the right-hand side of system (3) is locally Lipschitz on , such that , where , , and . According to [30], the solutions of (6) with initial conditions (4) exist uniquely on the interval , where ; therefore, all solutions exist on the first quadrant of the plane.
To prove the boundedness of solutions for system (3), let us first consider the case when the predator is absent, so thatwith initial conditions and ; we can easily show that for and , such that and for and , where . Adding the two equations of (8) yieldswhere . If we integrate both sides of (9), we getSince , the solutions are bounded, which clearly shows that .
To extend the analysis to (3), let us consider , . Also assume that and choose . By considering the derivative of , for for some fixed positive time T, we haveSince x and y are nonnegative and bounded by κ,Due to the positivity of z and since the parametric condition exists for ρ, the differential inequality is bounded above, such that , i.e., there exists N where for all , which implies the boundedness of z, such that

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 and 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 , , and , in . Therefore, the equilibria are the points of intersection of these zero growth isoclines regardless of the parameter values.

An interior equilibrium point exists with , and such that and , where is the root(s) the following equation:

The coefficients , , are defined by

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 , which leads to . 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 (), 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.

3.2. Stability and Bifurcation Analysis of the Interior Equilibrium

Now, we study the stability of the interior equilibrium , at which the Jacobian matrix is

Here,

The characteristic equation for the interior point is then given by

Here,such that

To gain insight regarding interior equilibrium , we discuss the stability of interior steady states and Hopf bifurcation conditions of the threshold parameters and by considering the following different cases.

Case 1. When , equation (17) becomesTherefore, the interior equilibrium is locally asymptotically stable if(i) holdsThus, based on Routh–Hurwitz Criteria, all the roots of (20) have negative real parts.

Case 2. For , then equation (17) becomesWe assume for some values of (), there exists a real number ω such that is a root of (21); then, we getSquaring and adding both of the equations yieldwhereBy Descartes’ rule of signs, equation (22) has at least one positive root if(i) holdsEliminating from (22) yieldswhere .
By differentiating (21) with respect to such that and , the transversality condition can be obtained in this form:Here,Then, a Hopf bifurcation occurs for if ; i.e., . We arrive at the following theorem.

Theorem 1. Let and hold, where ; then, there exists such that remains stable for and unstable for , where defined by (25). Moreover, system (3) undergoes a Hopf bifurcation at when .

Case 3. When , in the same manner of the pervious case, we arrive at the following theorem.

Theorem 2. For system (3), with , there exists a positive number , such that the equilibrium point is locally asymptotically stable for and unstable for , where . Furthermore, Hopf bifurcation occurs at :where .

Case 4. When and , we assume that as a variable parameter and as fixed on its stable interval. Let as a root of (17); separating real and imaginary parts impliesThus, eliminating the trigonometric functions ( and ) from (29) and (30) yieldswhereEquation (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 if(i) therefore, we haveHere,To study the Hopf bifurcation analysis, we fix in its stable interval and differentiate equations (29) and (30) with respect to . Then, substitute and , we havewhereFrom (35), we getAs , then Hopf bifurcation occurs for .
Therefore, we arrive at the following theorem.

Theorem 3. If exists, such that and hold, with , then there exists a positive threshold parameter such that the interior equilibrium is locally asymptotically stable for and unstable , where as in (38). Additionally, system (3) undergoes Hopf bifurcation at when .

Remark 2. Similarly, for , there exists a threshold parameter such that the interior equilibrium is locally asymptotically stable for and unstable . Also, Hopf bifurcation occurs for system (3) as , where is given bywithThe 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 .

Theorem 4. If and , then system (3) is globally asymptotically stable at the interior equilibrium point .

Proof. We suggest the Lyapunov function at of the form:where are nonnegative constants. Take derivative with respect to time t, yielding toThus, based on the assumptions: , , and , we can getHence, 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 and (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 are denoted by

Hence, sensitivity functions due to small perturbations in Allee parameter are given by a system of DDEs:

To estimate the sensitivity functions , , and , 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 satisfy the system of DDEs:

We then solve (45) along with (3) to evaluate , , and (see Figure 4).

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:

Figure 5 shows the numerical simulations of the delayed system (3) around the steady state . The interior steady state is asymptotically stable when and . The model undergoes a Hopf bifurcation when and . Figure 6 displays the Hopf bifurcation diagrams of and which are obtained numerically by maximum and minimum amplitudes of . The left banner displays the threshold parameter with , while right banner shows that the threshold parameter with .

Figure 5 displays a bistability of two interior equilibrium points, for the DDEs model (3), when parameter α varies from to . If the interior equilibria exist, any trajectory starting from the interior of converges to one of the interior equilibria.

Figure 7 shows the sensitivity of the dynamics of system (3) due to small changes in the severity of Allee effects and . The left banners show the numerical simulations with different values of () and fixed value of , while right banners display the simulations with different values of () and fixed value of . The phase portrait gets stretched over time as is reduced, while low values of 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: , , and to evaluate the sensitivity of the state variables due to a small perturbation in and . The oscillation behaviour indicates that the species population is very sensitive to small changes in the parameter. It is clear that and 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.

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

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by PhD project, United Arab Emirates University (UAE).