Bifurcation Behavior Analysis in a Predator-Prey Model

A predator-prey model is studied mathematically and numerically. The aim is to explore how some key factors influence dynamic evolutionary mechanism of steady conversion and bifurcation behavior in predator-prey model. The theoretical works have been pursuing the investigation of the existence and stability of the equilibria, as well as the occurrence of bifurcation behaviors (transcritical bifurcation, saddle-node bifurcation, and Hopf bifurcation), which can deduce a standard parameter controlled relationship and in turn provide a theoretical basis for the numerical simulation. Numerical analysis ensures reliability of the theoretical results and illustrates that three stable equilibria will arise simultaneously in the model. It testifies the existence of Bogdanov-Takens bifurcation, too. It should also be stressed that the dynamic evolutionary mechanism of steady conversion and bifurcation behavior mainly depend on a specific key parameter. In a word, all these results are expected to be of use in the study of the dynamic complexity of ecosystems.


Introduction
The dynamical behaviors between different populations and their complex properties have been given close attention by biologists and ecologists.Since the pioneering work of Lotka and Volterra, the research interest in predator-prey dynamics has achieved constant attention.It is well known that these models can directly reflect changes in the size of populations.Considerable improvements are that the relevant theories become more and more complete in this category in recent years [1][2][3][4][5].
The predator-prey models have extensive applicability in the field of biological problems.The biologist can use them to study the relationship between species in different domains [6][7][8][9][10].Yang and Zhao [11] have established a fish-algae consumption model to explore how to apply the complex dynamics between fish-algae populations to expound the mechanism of algae blooms; these results will be helpful in controlling algae bloom.González-Olivares and Rojas-Palma [12] have established a Gause type predator-prey model with Allee effects and considered three standard functional responses, respectively.They found that different types of functional responses will lead to a model's dynamic behavior change.Their results perfectly explore that the expression of one interaction term has a significant effect on the stability and persistence of the population model, which is meaningful in establishing biological mathematical model.Dhar et al. [13] have proposed a mathematical model to study how instantaneous nutrient recycling affects the dynamic characteristic of aquatic ecosystem.The nutrient supply rate was found to affect the local stability of equilibrium; this work was significant for models involving flowing waters.Luo [14] has considered a mathematical model to study how the periodic environment influences the internal operating characteristics of the aquatic ecosystem.He pointed out that the temperature has a certain time periodicity, the fluctuation of temperature can lead to changes in intrinsic carrying capacity, and the growth rate of prey populations can also be influenced by time periodicity; these results in [14] are more accordant with the actual situation.
It is a common phenomenon in nature that one predator lives on multiple prey species.To address this issue, some researchers have begun to consider alternative prey to depict the dynamic predator mechanism [15][16][17][18].It is easy to find out that alternative prey can variously affect the dynamic capture feature.On one hand, alternative prey can increase the predation quantity for the focal prey, because more prey biomass may result in higher predation rates for both prey items.On the other hand, the alternative prey population can also lower predation on the focal prey because of predator preference for the alternative prey resources [19].At the same time, many studies show that the alternative population can intensively influence the dynamic behavior of the aquatic ecosystem [20,21].Based on this mechanism, Kar and Chattopadhyay [22] have developed a two-species predator-prey model which includes the effect of alternative prey.The model can be depicted as where  and  represent prey and predator population densities or biomass at time , respectively.Here,  1 stands for the intrinsic growth rate of the prey without any environment limitations,  1 is the environmental carrying capacity of the prey in the absence of the predator,  1 /( 1 +) is the Holling type-II functional response [23], which is used to depict the average feeding rate of the predator when the predator spends time seeking prey, where  1 is the half saturation constant for the Holling type-II,  1 is the grazing rate of the predator population,  1 and  1 are the conversional rate and mortality rate of predator, respectively, and  1 (1 − / 1 ) indicates the portion of biomass of predator increments from the alternative prey, where  1 represents the growth rate of the predator on account of the alternative prey.From the formula, we can see that when the quantity of focal prey  approaches the environmental carrying capacity  1 , the amount of alternative prey consumed by the predator will tend to be zero [24].The concept of Allee effect was firstly derived from the research of Allee and Bowen [25].Since then, the Allee effect received the attention of many researchers [26][27][28][29].Allee effects can be roughly classified into two types: strong and weak [30].For these two forms, there is a critical value that is referred to as the Allee threshold, respectively.The fist form means that if the population size is below the threshold, the species will become extinct.When the growth rate gradually decreases but remains positive with a low population size, the Allee effect is described as weak.Aulisa and Jang [31] have established a continuous-time predator-prey model to study influences in dynamical behaviors when the prey population possesses Allee effect.They pointed out that both species will become extinct if the prey population size falls below a certain threshold.Pan et al. [32] have considered a reaction-diffusion phytoplankton-zooplankton model with double Allee effects on prey population.They pointed out that the Allee effects can make the dynamical behaviors of a system increasingly complex.
The strong Allee effect can be depicted by the following form: where  1 is Allee effect threshold.If population density or size is below the threshold, this population is doomed to extinction.Here 0 <  1 <  1 , the other parameters' significance is the same as model (1).Now we will establish a predator-prey model with strong Allee effects: The parameters are all greater than zero and have the same significance as above.For simplicity, we write the above model in dimensionless form as follows: we take the scaling  =  1 ,  =  1 / 1 , and  = / 1 ; then, model (3) can be simplified as where  =  1  1 / 1 ,  =  1 / 1 ,  =  1 / 1 ,  =  1 / 1 ,  =  1 / 1 , and apparently  > 1.
In the succedent subsections of this paper, we introduce the problem of determining the number of equilibria.Stability analysis of equilibrium point is also presented.Then, we provide a demonstration of several types of bifurcation and give a set of parameter values to prove the existence of Bogdanov-Takens (BT) bifurcation.In Section 3, numerical simulations are given to illustrate the theoretical analysis results, followed by conclusions in Section 4.
The Jacobian matrix of model (4a) and (4b) at the equilibrium point   is Theorem 1.
Proof.According to (6), Imitating the proof of Theorem 2, we can know then, det ( ( 4 , 4 ) ) = −( 4 )( 4 ) > 0, which means two eigenvalues have the same sign.Since Solving () = 0, we can get that there exists a set  ⊂ (1/, 1), such that () < 0 for all  ∈ , where  = ((1 +  −  + √1 −  +  +  Considering the existence and stability conditions integrated, we know that there will appear three stable equilibria synchronously if the scope of those parameters satisfies the above conditions.Based on the above demonstration, some cases where the stability of equilibria may occur are cited in Table 1, where the existence of equilibria is classified according to the magnitude of  and .

Local Bifurcation.
From Table 1 we know that the variation of parameter value will lead to the number change of interior equilibria.We take  as variable parameter and find that changing the value of  will vary equilibrium's number; when the value of  increases across the threshold  TC1 = /( + 1), the interior equilibrium  4 bifurcates from  1 , and when the value of  is across  TC2 = /( + 1) + (1 − 1/), another interior equilibrium  3 can bifurcate from  2 .Then, there exist two transcritical bifurcations, which are denoted by TC1 and TC2, respectively.Theorem 4. (1) Model (4a) and (4b) undergoes transcritical bifurcation at  1 (1, 0) when the value of parameter  equals the transcritical bifurcation threshold  1 = /( + 1).
Proof.It can be easily seen that  1 (1, 0) coincides with  4 ( 4 ,  4 ) when  = /( + 1).According to the theorems in [33,34], it can be found that one interior equilibrium point branches off from  1 when  passes the threshold  TC1 = /( + 1) and also conforms with the transversality condition for transcritical bifurcation.The same notation we followed in this paper is mentioned in [33].
Proof.Because the interior equilibrium point  3 is always a saddle, then Hopf bifurcation can only take place at  4 .Parameter  can drive equilibrium  4 into an unstable state when  >   , so  =   is the critical value where the stability of  4 changes.Next we will prove the necessary condition for Hopf bifurcation to occur.To testify the transversality condition of the Hopf bifurcation of the model's solution, we take  = () + (), where  is an eigenvalue of the Jacobian matrix  (( Using ( 16), we know Det ( ( 4 , 4 ,  ) ) > 0 and Tr ( ( 4 , 4 ,  ) )/ ̸ = 0.Then, the transversality condition of a Hopf bifurcation is satisfied [35].

Numerical Results
In order to verify the correctness and feasibility of the theoretical results, a series of numerical simulations will be depicted in detail.Many phase diagrams are given to display the dynamics properties of model ( 4a) and (4b), which are based on pplane8 routines [38] in Matlab 7.1.In fact, according to [39], we can know pplane8 is a powerful tool for studying planar autonomous systems of differential equations, which can rapidly and accurately draw trajectories of each phase plane, count each critical point, and correctly characterize each equilibrium point of the studied systems.It is easy to find from Table 1 that model (4a) and (4b) only has an interior equilibrium point  3 if  ≥ ; then, the premise of parametric ranges in Figure 1 is  > .It should be pointed out from Figure 1(a) that two vertical lines passing through points (/( + 1), 0) and (/( + 1) + (1 − /), 0) are the transcritical bifurcation curves, which have been named TC1 and TC2, respectively.Furthermore, if  < /(+1) is feasible, that is to say, model (4a) and (4b) does not have any interior equilibrium point when the value of  is positioned within region I of Figure 1(a), three boundary equilibria  0 (0, 0),  1 (1, 0), and  2 (2/3, 0) exist. 0 (0, 0) is asymptotically stable, while  1 (1, 0) and  2 (2/3, 0) are unstable, the results of which have been shown in Figure 2(a) ( = 0.131).We know those three boundary equilibria  0 (0, 0),  1 (1, 0), and  2 (2/3, 0) always exist no matter what the value of  is.Thus, in the following cases we do not introduce the existence of boundary equilibria separately.However, if the value of  is gradually increasing across the transcritical bifurcation TC1 and finally enters into region II, model (4a) and (4b) has an interior equilibrium point  3 .It is worth emphasizing that the critical value of transcritical bifurcation is  TC1 = /( + 1) = 0.13167 and the interior equilibrium point  3 (0.97064, 0.039763) is a saddle point, boundary equilibria  0 and  1 are stable, and  2 is an unstable node point, which has been shown in Figure 2(b) ( = 0.132).As the value of  gradually increases to  TC2 = /( + 1) + (1 − 1/) = 0.132083, which is a critical value of another transcritical bifurcation and finally enters into region III, model (4a) and (4b) has two interior equilibria  3 (0.91328, 0.093458) and  4 (0.71172, 0.052831), which are unstable;  0 and  1 are stable, and  2 is a saddle point (see Figure 2(c),  = 0.1325).
In order to clearly explore the steady characteristic of the interior equilibrium point  4 , Figure 1(c) will be given, which is the partially enlarged view of Figure 1(b).It should be stressed that the interior equilibrium  4 will change the stable state if the value of  increased beyond the line L1 and finally enters into region VI, which suggests that a Hopf bifurcation can lead to the appearance of a limit cycle in the vicinity of the interior equilibrium  4 (0.83816, 0.11816) if  = 0.1357878256 59 (see Figure 2(f)).From Figure 2(f), we can observe that this limit cycle around  4 is stable as it attracts two neighboring trajectories: the trajectory (bottle green curve) lying inside the limit cycle and the trajectory (blue curve) lying outside; these two trajectories move ectad and entad, respectively, and converge on the limit cycle.But at this time the interior equilibrium point  3 (0.85397, 0.11709) and boundary equilibrium point  2 (2/3, 0) are unstable, and  0 (0, 0) and  1 (1, 0) are stable.When the value of  enters into domain VI,  4 (0.86298, 0.042889) becomes stable,  3 (0.92702, 0.032122) remains unstable, and  0 and  1 are all stable in this domain.Thus, it is interesting to know that model (4a) and (4b) will show three stationary phenomenon (see Figure 2(g),  = 0.13579).Due to the parameter values, the nature of  3 and  4 only can be seen clearly on the enlarged view (see Figure 2(g)).Therefore, we take another group of values to exhibit the three stable states in overall view in Figure 3(b).If the value of  gradually increases and reaches  SN =  +  + 2 √  +  = 0.1357900212, it is easy to find that  3 (0.84605, 0.11789) and  4 (0.84605, 0.11789) coincide at line L2 and the coincident point is called saddlenode point, which has features of both saddle and node points (see Figure 2(h)).However, it will disappear if the value of  is increased higher than  SN , which signifies that model (4a) and (4b) will undergo a saddle-node bifurcation if the value of  increases across the threshold value  SN = ++2 √ +.In a word, it is worthy of our summary that model (4a) and (4b) can show a complex dynamic evolutionary process of steady conversion and bifurcation behavior with increase of key parameter .
In addition, model (4a) and (4b) has three boundary equilibria  0 (0, 0),  1 (1, 0), and  2 (2/3, 0) and an interior equilibrium point  3 , which is unstable if  >  or  = .However, it is worth stressing that the boundary equilibrium point  0 (0, 0) is a saddle point if  >  and is an unstable high-order singularity if  =  (see Figures 4(a Based on the above analysis, the key parameter  can impose influence on dynamic evolutionary mechanism of steady conversion and bifurcation behavior and lead to model (4a) and (4b) having three stationary phenomena and multiple bifurcation behaviors, which can in turn prove that the theoretical results are correct and the complex dynamics of model (4a) and (4b) mainly depend on some key parameters.Moreover, these results show that the method of using mathematical model to study the ecological problems is feasible.

Conclusions
On the basis of the theories and methods of ecology, a predator-prey model is studied numerically and analytically in this paper.The aim is to probe how some key factors influence dynamic evolutionary mechanism of steady conversion and bifurcation behavior in predator-prey model.The theoretical works have been promoting the investigation of the existence and stability of the equilibria, as well as some conditions of some bifurcations behaviors, such as transcritical bifurcation, saddle-node bifurcation, and Hopf bifurcation, which can deduce a standard parameter controlled relationship and in turn provide a theoretical basis for the numerical simulation.Numerical analysis indicates that the dynamic evolutionary mechanism of steady conversion and bifurcation behavior mainly depend on a specific key parameter .Within this framework, the direct and indirect effects caused by the specific key parameter  are investigated by means of bifurcation analysis and phase diagram.It is obvious to find that the existence and stability of interior equilibria  3 and  4 mainly depend on a key parameter .These results suggest that the key parameter  plays an important role in the prey-predator model.In addition, when the value of the key parameter  is larger than some critical value, the model can possess multiple bifurcation behaviors, such as transcritical bifurcation, saddle-node bifurcation, and Hopf bifurcation.Thus, it is worthwhile to remark that the key parameter  has a profound effect on the population bifurcation dynamical behaviors.In a word, some key parameters can alter population dynamics and features in prey-predator model.In addition, it is our hope that all these results can be applied in the study of the dynamic complexity of ecosystems.

Figure 1 :
Figure 1: Bifurcation diagram of model (4a) and (4b) in  −  plane constructed for  = 1.5,  = 2.0,  = 0.1,  = 0.395 (a), and  = 0.405 (b, c).Red perpendicular dashed lines represent the critical value for different bifurcation occurrence.Three horizontal lines stand for three boundary equilibria:  0 (golden),  1 (gray), and  2 (black).Two interior equilibria  3 and  4 are presented by blue and green curves, respectively.The solid lines indicate equilibrium in stable state and dotted line indicates unstable state, where, according to the stability theorem, these results are calculated and drawn based on Maple 14 platform.

Figure 2 :Figure 3 :Figure 4 :
Figure2: Phase portraits for different concomitant case of the model's equilibria.The horizontal axis is prey population  and the vertical axis is predator population .The green curves are stable or unstable orbits; the red point is the equilibrium point.We take  = 1.5,  = 2.0,  = 0.1,  = 0.395 in (a-c), and  = 0.405 in (d-h).