Dynamical Behavior of a Modiﬁed Leslie–Gower One Prey–Two Predators with Competition

: We study the dynamics of a modiﬁed Leslie–Gower one prey–two predators model with competition between predator populations. The model describes complex dynamics in the permanence, global stability and bifurcation. It is shown that there are eight possible equilibrium states. Two equilibrium states, i.e., the extinction of all of the species state and the extinction of both predators state are always unstable, while the other equilibrium states are conditionally locally and globally asymptotically stable. We also analyzed numerically the e ﬀ ect of competition between predators. Our numerical simulations showed that the competition rate of the second-predator may induce the transcritical bifurcation, the saddle-node bifurcation as well as the bi-stability phenomenon. Such numerical results are consistent with the analytical results. and A.; Writing—original draft, D.S.; Writing—review


Introduction
The model of population interaction has attracted great interest of many scholars. One of the important mathematical models is the model of the predator-prey interaction [1]. A well-known model that describes the predator-prey interaction is the Leslie-Gower model [2]. In this model, the predator is assumed to grow logistically, where its carrying capacity is proportional to the density of prey. In the case of prey scarcity, predators can consume other populations but their growth will be bounded by the fact that their main prey is not abundantly available. To account for such a situation, Aziz-Alaoui and Okiye [3] proposed a modified Leslie-Gower model by introducing a constant that measures the environment protection for the predator. Since then, many researchers have studied the modified Leslie-Gower models with various types of functional responses [4,5], harvesting [6][7][8][9], Allee effect [10][11][12][13], etc.
The model of predator-prey interaction is one of the food chain building block models that involves two species. Although predator-prey interaction models have been widely studied or developed, many ecological and modeling problems remain open. One of the developments that has been made includes adding another prey or another predator. For example, some scholars have proposed mathematical models that describe the interaction between two preys and one predator and have studied dynamics such as the stability of equilibrium states to determine conditions of coexistence in ecosystems, species extinction, existence of stable periodic solutions or transitions between prey and predator [14][15][16][17][18]. There are also many mathematical models describing the interaction of one prey with two predators. Sarwandi et al. [19] studied the effects of a prey refuge on the two competing predators. Sayekti et al. [20] considered the effects of prey harvesting on the interaction of one prey and two predators. They found that coexistence can be obtained if the prey harvesting rate is less than the prey intrinsic growth rate. We notice that all above mentioned models are derived based on the Lotka-Volterra predator-prey model. Recently, Alebraheem and Abu-Hasan [21][22][23] introduced mathematical models of one prey-two predators interactions, which are based on the Leslie-Gower model. Hence, they considered that the two predators grow following a logistic model, where the carrying capacity is only dependent on the prey density.
Recently, we [24] revisited the model of Alebraheem and Abu-Hasan [23]. We assumed in [24] that without competition, the interaction of prey with the two predators follows the Leslie-Gower model with a Holling type I functional response. If competition is also considered, then the model is where x(t), y(t) and z(t) represent the population density of the prey, the predator-1 and the predator-2 at time t, respectively and r, r 1 and r 2 are, respectively, the intrinsic growth rate of x, y and z.
The carrying capacity of prey is k. Parameters β 1 and β 2 are the maximum value of the per capita reduction rate of x due to y and z, respectively. Parameters α 1 and α 2 denote the competition rate of predator-1 and predator-2, respectively. All parameters in the model are positive. As opposed to the model in [23], the intrinsic growth rates of both predators in system (1) are constant. In [23], the intrinsic growth rates of both predators are proportional to the density of prey. The dynamical behavior of system (1) was studied by Savitri et al. [24]. They only focused on local stability analysis and show some phase portraits by numerical simulations; transcritical bifurcation and saddle-node bifurcation have not been observed. As explained by Aziz-Alaoui and Okiye [3], the Leslie-Gower term (y/x) in the system (1) describes the reduced predator population due to scarcity of its favorite foods. In this article it is assumed that if there is a severe scarcity of prey, predators can consume other populations, but their growth will be bounded by the fact that their favorite food (x) is not available in large quantities. To take into account such a situation, following Aziz-Alaoui and Okiye [3], we add a constant in the denominator to obtain a modified Leslie-Gower one prey-two predators model with competition between predators.
where k 1 and k 2 are measures of environmental protection for predator-1 and predator-2, respectively.
In this article we investigate the dynamics of system (2), which includes the permanence, the existence of equilibrium states and local and global stability of equilibrium states. We also discuss the bifurcation numerically using MatCont. This article is organized in the following manner. In Section 2, we give some preliminary theorems and show that all solutions of the model are always non-negative, permanent and bounded. We analyze the existence of the equilibrium states, and the local and global stability in Section 3. In Section 4 we show numerical simulations to illustrate the analytical results and explore transcritical bifurcation, saddle-node bifurcation and bistability. Finally, we end with a brief conclusion in Section 5.

Non-Negativity
Notice that model (2) describes the interaction of three populations, and therefore, for biological relevance, we require that the solutions of this model must be non-negative. The non-negativity of solutions are assured by the following theorem: Theorem 1. All solutions of system (2) with initial values x(0) > 0, y(0) > 0 and z(0) > 0 are always non-negative.

Boundedness
In nature, resources for every population are always limited such that the population should also be bounded. The boundedness of population prey and predators as solution of system (2) is demonstrated by utilizing the following comparison lemma: [25].
Proof. From the first equation of system (2), we have Since r > 0 and k > 0, Lemma 1 gives Then, for any sufficiently small ε > 0, there is a t 1 ≥ 0 such that Based on this result and the second equation of system (2) we have From Lemma 1 we obtain lim t→+∞ supy(t) ≤ k + ε + k 1 .

Permanence
In this part, we show the permanence of system (2). System (2) is said to be permanent if there exist positive constants l and L such that Proof. From Theorem 2, we have that for any sufficiently small ε > 0, there are a t 2 , t 3 ≥ 0 such that LetT = max{t 1 , t 2 , t 3 }. Using Theorem 2 and the first equation in system (2), we can show that for any t ≥T, By applying Lemma 1 and taking ε → 0 we obtain Similarly, using the second equation in system (2), we get for any t ≥T, By using Lemma 1 and letting ε → 0 we obtain As in the previous case, the third equation in system (2) gives and therefore (by Lemma 1 and taking ε → 0 ) we get Observe that if (β 1 (k+k 1 )+β 2 (k+k 2 )) r < 1, Hence, by choosing l ≤ min{l 1 , l 2 , l 3 } and L ≥ max{L 1 , L 2 , L 3 }, we can easily see that system (2) is permanent. This completes the proof.

Equilibrium State
By solving the following system of equations, we can show that system (2) has eight feasible positive equilibrium states: All equilibrium states and their existence conditions are shown in Table 1. We remark that x * in the equilibrium state E 8 = (x * , y * , z * ) is the positive solution of the cubic equation where

Equilibrium State Existence Condition
The existence of positive root x * of a cubic equation can be easily derived by Cardano's criteria. The detail of Cardano's criteria can be seen, for example, in [10] and is not discussed here. Furthermore, the values of y * and z * are respectively given by We also require that y * and z * be positive.

Local Stability of Equilibrium States
The local behavior of system (2) is investigated by considering its linear approximation around each equilibrium state. The variational matrix of the linearized system around a point (x, y, z) is given by The stability of an equilibrium state (x,ŷ,ẑ) is determined by the eigenvalues of the variational matrix evaluated at the corresponding equilibrium state, i.e., J(x,ŷ,ẑ).
Proof. See Appendix B.
Proof. See Appendix C.
Proof. See Appendix D.

Remark 1.
(i) We notice that model (1) has only four equilibrium states (see [24]). This is a consequence of the use of the Leslie-Gower model, which does not allow the extinction of prey, and therefore equilibrium states E 1 , E 3 , E 4 and E 7 do not exist in the model (1). In the model (2) we introduced environmental protection for both predators, which are denoted by k 1 and k 2 , respectively. Both constants cause model (2) to have an additional four equilibrium states where each density of prey equilibrium states is zero, i.e., E 1 , E 3 , E 4 and E 7 . (ii) From Theorems 4-8, we can classify the stability of each equilibrium state, as seen in Table 2. The extinction of all of the species state and the extinction of the both predators state are always unstable, while other equilibrium states are conditionally stable. This can be interpreted as meaning that the two predators cannot become extinct at the same time, even though there are no prey in the environment. This phenomenon can be understood from model (2) that predator-1 and predator-2 have environmental protection. We can see in Table 2 that the environmental protection is very important in determining the stability of the equilibrium points.

Global Stability of Equilibrium States
In the following we present the global stability analysis of some equilibrium states. The global stability is performed by choosing an appropriate Lyapunov function for each equilibrium state. Theorem 9. The equilibrium state E 3 = (0, 0, k 2 ) of the system (2) is globally asymptotically stable if α 1 l 3 > r 1 and 1 + k 2 Proof. We assume that α 1 l 3 > r 1 and construct a Lyapunov function The time derivative of this Lyapunov function is It is clear that if 1 + k 2 l 1 +k 2 < l 3 k+k 2 then dV 1 dt ≤ 0. We observe that dV 1 dt = 0 if and only if x = 0, y = 0 and z = k 2 . Hence, according to La Salle's invariant principle, the equilibrium state E 3 = (0, 0, k 2 ) is globally asymptotically stable if α 1 l 3 > r 1 and 1 + k 2 Theorem 10. The equilibrium E 4 = (0, k 1 , 0) of the system (2) is globally asymptotically stable if α 2 l 2 > r 2 and 1 + k 1 Proof. The theorem can be proven using a similar procedure as in the Proof of Theorem 9, i.e., by choosing a Lyapunov function V 2 = y − k 1 − k 1 ln y k 1 + bz. The detail of the Proof can be seen in Appendix E.
Theorem 11. The coexistence equilibrium E 8 = (x * , y * , z * ) of system (2) is globally asymptotically stable if Proof. The theorem can be proven by defining a Lyapunov function Detail of the Proof can be seen in Appendix F.

Numerical Simulations
In this section we solve model (2) using the 4th-order Runge-Kutta method with parameter values given in Table 3. The numerical simulations are performed to study the dynamics of model (2), particularly to investigate the change of dynamical behavior due to the change of the competition rate of predator-2 (α 2 ). To give a view on the dynamical behavior, we consider four different cases (α 2 = 0.2, α 2 = 0.4, α 2 = 0.5, α 2 = 0.6) as the consequence of the previous stability analysis. Such dynamical behavior is depicted in Figure 1, which shows some phase portraits of model (2) for four different values of α 2 . Proof. The theorem can be proven using a similar procedure as in the Proof of Theorem 9, i.e., by Detail of the Proof can be seen in Appendix F. □

Numerical Simulations
In this section we solve model (2) using the 4th-order Runge-Kutta method with parameter values given in Table 3. The numerical simulations are performed to study the dynamics of model (2), particularly to investigate the change of dynamical behavior due to the change of the competition rate of predator-2 ( ). To give a view on the dynamical behavior, we consider four different cases ( = 0.2, = 0.4, = 0.5, = 0.6) as the consequence of the previous stability analysis. Such dynamical behavior is depicted in Figure 1, which shows some phase portraits of model (2) for four different values of .    Case (1). When α 2 = 0.2, it is found that 0.88 = r 1 r 2 > α 1 α 2 k 1 k 2 = 0.252, and therefore according to Theorem 7, E 7 = (0.0, 0.31529, 1.31974) is asymptotically stable. In this case, no other equilibrium state is stable (see Figure 1a).
Case (2). If we increase the competition rate of predator-2 such that α 2 = 0.4, then we have a unique coexistence equilibrium state, namely E 8 = (0.03305, 0.47462, 1.18577). E 8 is asymptotically stable since the Routh-Hurwitz criterion in Theorem 8 is satisfied. Our numerical simulation shown in Figure 1b confirms this stability property.
Case (4). Finally, we take α 2 = 0.6 and show the result of our numerical simulation in Figure 1d. From Figure 1d, it is observed that only E 6 = (0.7875, 2.5875, 0.0) is stable. This can be understood from the fact that only the stability condition in Theorem 7 is satisfied (3.52 = r 2 (kb 1 + r) < rα 2 (k + k 2 ) = 4.96), while the stability conditions for other equilibrium states are not fulfilled.
From previous numerical simulations, we observe that the dynamical property of model (2) changes as we vary the value of α 2 . To understand the dynamics of model (2) with parameters values in Table 3, we check the existence of equilibrium states and their local stability properties using the previous theoretical results. From the existence conditions, it is found that our model has the following equilibrium states: E 1 = (0, 0, 0), E 2 = (2.8, 0, 0), E 3 = (0, 0, 1.4), E 4 = (0, 1.8, 0) and E 6 = (0.7875, 2.5875, 0). The extinction of the predator-1 equilibrium state (E 5 ) does not exist, while the prey extinction equilibrium state (E 7 ) exists if α 2 < 0.6111. Based on Theorems 4-7, we can check that E 1 , E 2 , E 3 and E 4 are unstable. E 7 is stable if α 2 < α * 1 2 = 0.330499, and the stability of E 6 is achieved whenever α 2 > α * 2 2 . The density of the prey population in the coexistence equilibrium state is determined by the cubic Equation (7). We applied the well-known Cardano formula in MAPLE Software to verify all possible real roots of Equation (7) as well as to check the existence of the coexistence equilibrium state. For biological reasons, the coexistence equilibrium state exists if the densities of prey (x * ), predator-1 (y * ) and predator-2 (z * ) are positive. Based on the cubic Equation (7), there are three possible roots. We denote those possible roots as x * 8 , x * 8a and x * 8b and the respective coexistence equilibrium states as E 8 = x * 8 , y * 8 , z * 8 , E 8a = x * 8a , y * 8a , z * 8a and E 8b = x * 8b , y * 8b , z * 8b . We found that one of roots of Equation (7), say x * 8 , is a real number if α 2 < α * 3 2 = 0.560599. Under this condition, y * 8 and z * 8 are positive, but x * 8 is positive only if α 2 > α * 1 2 = 0.330499. Hence, E 8 exists if α * 1 2 < α 2 < α * 3 2 . The second root x * 8a is a positive real number if α 2 < α * 3 2 . However, y * 8a > 0 and z * 8a > 0 can be achieved only if α 2 > α * 2 2 = 0.425121, so that E 8a exists if α * 1 2 < α 2 < α * 3 2 . Finally, it is found that the third root x * 8b is also a real number whenever α * 3 2 < α 2 < 13.7915378. However, x * 8b < 0 and therefore E 8b is not biologically feasible. The stability of the coexistence equilibrium point can be checked by applying the Routh-Hurwitz criteria (see Theorem 8). Due to the explicit forms of the coexistence equilibrium points (x * 8 and x * 8a ) being very complicated, we investigated the stability of coexistence equilibrium points numerically using a numerical continuation package MatCont. The numerical continuation was performed not only to study the stability of the coexistence equilibrium points, but also to examine all possible equilibrium states and their stability properties. In Figure 2 we plot the equilibrium state density of predator-1 (y) as the value of α 2 is changed. It is noticed that the equilibrium state density of prey (x) and that of predator-2 (z) are similar to this graph. It appears that our theoretical results are confirmed by the results of numerical continuation. The detail of the dynamics shown in Figure 2 could be explained as follows.

Transcritical Bifurcation
From Figure 1a,b we observe the exchange of stability between the prey extinction equilibrium state ( ) and the coexistence equilibrium state ( ) , indicating there occurs a transcritical bifurcation caused by . The transcritical bifurcation is more clearly seen in Figure 2. For a relatively small , only the equilibrium state is stable. As increases, there appears a branch point (BP1) at * . This BP1 indicates the appearance of transcritical bifurcation, i.e., there is a change of stability between the extinction of the prey equilibrium state ( ) and the coexistence equilibrium state ( ). When < * , is asymptotically stable while is unstable. On the other hand, if * < < * , then is unstable while becomes asymptotically stable. This shows that when the competition rate of predator-2 is relatively low, the predator-2 can grow well, while if the competition rate of predator-2 is large enough, then the growth of predators is reduced and therefore predation is decreased. As a result, prey can survive, which is shown by the stability of the coexistence equilibrium state ( ).
Another transcritical bifurcation may also occur between and , where the bifurcation point is at BP2 (or = * ). We remark that those phenomena of transcritical bifurcations are only possible mathematically. Indeed, does not biologically exist for < * , since in this case the value of * is negative. Similarly, does not exist for biological reasons for < * . Thus, the transcritical bifurcations in our case are not biologically feasible.

Saddle-Node Bifurcation
As seen in Figure 2, if * < < * , then there exists a unique coexistence equilibrium state ( ), which is asymptotically stable (see also Figure 1b). A further observation shows that for * < < * there are two coexistence equilibrium states. One of them is asymptotically stable while the other one is unstable. The stability of one of these coexistence equilibrium states can also be seen in Figure 1c. Then, these two coexistence equilibrium states collide at the limit point (LP), i.e., at = * . If > * , then the two coexistence equilibrium states disappear and only the extinction of predator-2 equilibrium state is stable. Hence there occurs a saddle-node bifurcation driven by .

Bi-Stability
Another interesting dynamics that can be observed in Figure 2 is the appearance of two different stable equilibrium states. This interesting occurrence is known as the bi-stability phenomenon. The bi-stability phenomenon occurs when * < < * . In this case, both extinction of the predator-2

Transcritical Bifurcation
From Figure 1a,b we observe the exchange of stability between the prey extinction equilibrium state (E 7 ) and the coexistence equilibrium state (E 8 ), indicating there occurs a transcritical bifurcation caused by α 2 . The transcritical bifurcation is more clearly seen in Figure 2. For a relatively small α 2 , only the equilibrium state E 7 is stable. As α 2 increases, there appears a branch point (BP 1 ) at α * 1 2 . This BP 1 indicates the appearance of transcritical bifurcation, i.e., there is a change of stability between the extinction of the prey equilibrium state (E 7 ) and the coexistence equilibrium state (E 8 ). When α 2 < α * 1 2 , E 7 is asymptotically stable while E 8 is unstable. On the other hand, if α * 1 2 < α 2 < α * 2 2 , then E 7 is unstable while E 8 becomes asymptotically stable. This shows that when the competition rate of predator-2 is relatively low, the predator-2 can grow well, while if the competition rate of predator-2 is large enough, then the growth of predators is reduced and therefore predation is decreased. As a result, prey can survive, which is shown by the stability of the coexistence equilibrium state (E 8 ).
Another transcritical bifurcation may also occur between E 6 and E 8a , where the bifurcation point is at BP 2 (or α 2 = α * 2 2 ). We remark that those phenomena of transcritical bifurcations are only possible mathematically. Indeed, E 8 does not biologically exist for α 2 < α * 1 2 , since in this case the value of x * 8 is negative. Similarly, E 8a does not exist for biological reasons for α 2 < α * 2 2 . Thus, the transcritical bifurcations in our case are not biologically feasible.

Saddle-Node Bifurcation
As seen in Figure 2, if α * 1 2 < α 2 < α * 2 2 , then there exists a unique coexistence equilibrium state (E 8 ), which is asymptotically stable (see also Figure 1b). A further observation shows that for α * 2 2 < α 2 < α * 3 2 there are two coexistence equilibrium states. One of them is asymptotically stable while the other one is unstable. The stability of one of these coexistence equilibrium states can also be seen in Figure 1c. Then, these two coexistence equilibrium states collide at the limit point (LP), i.e., at α 2 = α * 3 2 . If α 2 > α * 3 2 , then the two coexistence equilibrium states disappear and only the extinction of predator-2 equilibrium state is stable. Hence there occurs a saddle-node bifurcation driven by α 2 .

Bi-Stability
Another interesting dynamics that can be observed in Figure 2 is the appearance of two different stable equilibrium states. This interesting occurrence is known as the bi-stability phenomenon.
The bi-stability phenomenon occurs when α * 2 2 < α 2 < α * 3 2 . In this case, both extinction of the predator-2 equilibrium state (E 6 ) and the coexistence equilibrium (E 8 ) are asymptotically stable. Here the dynamics of model (2) are very sensitive to the initial values. For some relatively large initial values of predator-2, we have a situation where the solutions converge to the coexistence equilibrium state (E 8 ), as shown by the black-trajectories in Figure 1c. Hence, all three species (prey, predator-1 and predator-2) can coexist. On the other hand, whenever the initial predator-2 is relatively small, then the solutions converge to the extinction of the predator-2 equilibrium state (E 6 ) (see blue-trajectories in Figure 1c).

Conclusions
We have discussed a model that describes the interaction of one prey-two predators with competition between the two predators. The model is derived from the modified Leslie-Gower predator-prey model. The model is proven to be permanent and all solutions of the model are always non-negative and bounded. We found that the proposed model has eight equilibrium states. There are two equilibrium states that are always unstable, namely the extinction of all of the species and the extinction of the two predators. Other equilibrium states are conditionally stable. We have also shown the global stability of some equilibrium states by defining suitable Lyapunov functions. The results of dynamical analysis were confirmed by our numerical simulations. Moreover, our numerical simulations have shown that the proposed model exhibits rich dynamics as we observed the occurrence of transcritical bifurcation, saddle-node bifurcation and bi-stability phenomenon, which are driven by the competition rate of predator-2 species.

Conflicts of Interest:
The authors declare that there are no conflicts of interest regarding the publication of this paper.