Hopf bifurcation and stability analysis of the Rosenzweig-MacArthur predator-prey model with stage-structure in prey

Abstract: We consider a stage-structure Rosenzweig-MacArthur model describing the predator-prey interaction. Here, the prey population is divided into two sub-populations namely immature prey and mature prey. We assume that predator only consumes immature prey, where the predation follows the Holling type II functional response. We perform dynamical analysis including existence and uniqueness, the positivity and the boundedness of the solutions of the proposed model, as well as the existence and the local stability of equilibrium points. It is shown that the model has three equilibrium points. Our analysis shows that the predator extinction equilibrium exists if the intrinsic growth rate of immature prey is greater than the death rate of mature prey. Furthermore, if the predation rate is larger than the death rate of predator, then the coexistence equilibrium exists. It means that the predation process on the prey determines the growing effects of the predator population. Furthermore, we also show the existence of forward and Hopf bifurcations. The dynamics of our system are confirmed by our numerical simulations.


Introduction
Mathematically, there are some complex phenomena in nature that are appealing to be investigated. One example is the interaction between prey and its predator. The first predator-prey interaction model was developed independently by [1] and [2], and therefore it is known as the Lotka-Volterra predatorprey model. Since then, the model for the interaction between prey and its predator is continuously studied and developed by both mathematicians and biologists [3]. One of important features that determine the characteristics of predator-prey interaction is the functional response which explains how many prey are eaten by predators per unit time. For example, the famous Rosenzweig-MacArthur predator-prey model implements the Holling-Type II functional response p(x) = ax b+x , where a and b are positive parameters related to the predator's attack rate and the handling time, respectively [4,5]. The Rosenzweig-MacArthur model has attracted many scholars to study the effect of various factors on predator-prey interactions, for example the effect of refuge [6,7,8,9], harvesting [10], quiescence [11], disease transmission among predators [12] or among prey [13] and the mutant population [14].
The transition between immature and mature stages from an adaptive dynamics perspective has recently been evaluated in [15]. In the predator-prey interaction, the age or size of the prey certainly determines the level of predation. Even in many cases, predators are only able to attack when the prey is still immature. Indeed, small size or immature prey tend to be more easily captured and eaten by predators. Moreover, there are many factors that significantly influence the dynamics of the predatorprey models with stage-structure. For example, [16] and [17] have studied the harvesting policy in a stage-structure predator-prey model. In [18], a ratio-dependent prey-predator model with stagestructure was investigated, and the uniform persistence and impermanence of the model were carried out through the sufficient conditions. Other features which influence the interaction of predator-prey with stage-structure have also been studied in many references, such as maturation delay [19,20,21], periodic functional responses [22,23], prey refuge [24], combination of prey refuge and additional food for predator [25], anti-predator [26,27,28], as well as diffusive effects [29].
Recently, the authors [30] have studied a Rosenzweig-MacArthur predator-prey model with stagestructure in prey: where x 1 ≡ x 1 (t) and x 2 ≡ x 2 (t) are respectively the densities of immature and mature prey population at time t, and y ≡ y(t) is the density of predator population at time t. The following assumptions are taken in deriving model (1.1) : (H1) The immature prey do not produce offspring and their growth rate depends entirely on the reproduction by the mature prey, where the growth is assumed to be logistically with constant intrinsic rate r > 0 and constant carrying capacity K > 0.
(H2) The immature prey is more vulnerable so predators only consume immature prey. The predation mechanism is assumed to follow the Holling Type II functional response where the maximum predation rate is v 0 > 0 and the environmental protection for immature prey is denoted by n 1 > 0.
(H3) The immature prey grow up and turn into mature prey with a conversion rate α > 0.
(H4) The growth of mature prey population depends only on the conversion of immature prey into mature prey. The mature prey do not have a risk to be attacked by predator. The death rate of the mature prey population is δ 0 > 0.
(H5) The predator consumes only the immature prey and the growth rate of predator population is proportional to the predation rate, where the conversion rate of consumed prey into predator births is c > 0. The death rate of the predator population is δ 1 > 0.
Notice that the growth of immature prey in model (1.1) is inhibited by the intraspecific competition between the mature and immature prey. In nature, the intraspecific competition also occur between mature prey. Such intraspecific competition will be incoporated in this paper and hence we propose the following stage-structure predator-prey model System (1.2) describes a predator-prey interaction where the predator only consumes immature prey. One ecological example of such a predator-prey system is the interaction between Moluccan megapode (Eulipoa wallacei) and rats Rattus sp. Rattus sp. is reported to prey upon Eulipoa eggs and chicks, and it cannot attack adult Eulipoa, see [31]. Another example of predator-prey system which can be described by model (1.2) is Chinese fire-bellied newt, which is unable to feed on the mature Rana chensinensis, can only prey on its immature, see [32]. For the sake of convenience, we describe all variables and parameters in system (1.2) in Table 1. We remark that the intraspecific competition between the mature and immature prey is assumed to have the same strength as that between mature prey. To simplify our analysis, we reduce the number of parameters in system (1.2) by considering the following non-dimensional model As far as we are aware, the dynamics of system (1.3) have not been studied. Hence, a Rosenzweig-MacArthur predator-prey model with stage-structure in prey is introduced in this paper and the dynamical properties of the model are investigated. The organization of this paper is as follows. In Section 2, we analyze the boundedness, existence and uniqueness of the solution of the system (1.3). In Section 3 we determine all possible equilibrium points including the conditions for their existence and stability properties. The existence of Hopf-bifurcation and numerical simulations to illustrate our analytical finding are presented in Section 4 and Section 5, respectively. The conclusion of our study is given in the last section.

Existence and uniqueness of solutions
The existence and uniqueness of solutions of system (1.3) will be investigated in the region [0, ∞) × Ω M where Ω M = (u 1 , u 2 , u 3 ) T ∈ R 3 + : max {|u 1 | , |u 2 | , |u 3 |} ≤ M for sufficiently large M. The existence of M is assured by the boundedness of the solutions, which will be discussed later. Let U = (u 1 , u 2 , u 3 ) T ,Ū = (ū 1 ,ū 2 ,ū 3 ) T and consider a mapping For any U,Ū ∈ Ω M , we can show that Clearly that H(U) satisfies the Lipschitz condition with respect to U, and therefore system (1.3) with any positive initial condition Hence, we have the following theorem related to the existence and uniqueness of the system (1.3).

Boundedness of solutions
It is clear that d dτ (u 1 + u 2 ) ≤ 0, and thus u 1 + u 2 is monotone decreasing. Let's denote Based on Barbalat lemma, we can show that if κ > 1 then This contradiction leads to κ = 1. Since u 1 + u 2 is monotone decreasing function, we have that meaning that there exists τ 0 > 0 such that for τ > τ 0 we get u 1 + u 2 ≤ 1 + or u 1 ≤ 1 + − u 2 for any ≥ 0. Then, from the second equation in system (1.3) we get for any τ ≥ τ 0 that It is obvious that (2.10) By taking = 0, we get Similarly, we can also show that Hence u 1 and u 2 is bounded above. The latter inequality shows that there exists τ 1 such that u 1 ≤ β 1 Using the first and the third equations in system (1.3), we obtain (2.13) By taking = 0 and denoting η = |β 2 − α 1 | β 1 α 1 +β 1 , we get the following inequality from which we have Hence, u 3 is also bounded above because Now we suppose that our assumption that u 1 + u 2 ≥ 1 is violated, i.e. there exists τ 2 > 0 such that we have for the first time u 1 (τ 2 ) + u 2 (τ 2 ) = 1. Then, we easily have Consequently, whenever once a solution with u 1 + u 2 has entered the interval (0, 1) then it remains bounded there for all τ > τ 2 , i.e.
From the second equation in system (1.3) we get Using the previous argument, we can prove that Similarly, we can also show that Again, by defining u = u 1 + u 3 α 2 and using the previous argument, we can easily show that u is bounded above and therefore u 3 is also bounded.

Existence and stability analysis of equilibrium points
It is easy to show that system (1.3) has three non-negative equilibrium points as follows.
From the description above, it can be seen that if the conditions (3.2) are fulfilled then all three equilibrium points exist. The local stability of each equilibrium points of system (1.3) is shown in the following theorem. (ii) The equilibrium point E 1 is locally asymptotically stable if ω > ω.
Proof. The local stability of all equilibrium points can be studied from the linearization of the system (1.3). The Jacobian matrix of the system (1.3) at a point (u 1 , u 2 , u 2 ) is given by By observing the eigenvalues of the Jacobian matrix (3.3) at each equilibrium point, we have the following stability properties.
(iii) The characteristic equation of the Jacobian matrix of the system (1.3) at E * is given by the following cubic equation where , The stability of E * can be determined by the Routh-Hurwitz criterion [33], i.e. E * is locally asymptotically stable if ϕ i > 0, i = 1, 3 and From conditions (3.1) and (3.4), we note that if the intrinsic growth rate of immature prey is less than the death rate of mature prey then the extinction point E 0 will be locally asymptotically stable. Otherwise, the predator-free point E 1 exists whenever the death rate of mature prey is less than the intrinsic growth rate of immature prey. Thus, the existence of predator-free point E 1 and the stability of the extinction point E 0 are dependent on both intrinsic growth rate of immature prey and the death rate of mature prey. Next, condition (3.2) says that the coexistence equilibrium E * exists if the death rate of the predator is less than the predation rate.
Based on the existence and stability results of equilibrium points of the system (1.3), we remark that when ω > ω then the coexistence equilibrium point (E * ) does not exist. In this case, parameter β 1 is a threshold: if it is greater than unity then the asymptotically stable extinction equilibrium point (E 0 ) is the unique equilibrium point; if it is less than unity then E 0 becomes unstable and there appears an asymptotically stable equilibrium point (E 1 ). Hence the system (1.3) undergoes a forward bifurcation at β 1 = 1. In fact a forward bifurcation may also be observed when β 1 < 1, where in this case ω behaves as the threshold parameter. Since β 1 < 1, the extinction point (E 0 ) is unstable. Furthermore, when ω > ω, equilibrium point E 1 is asymptotically stable and E * does not appear. On the contrary, if ω < ω then E 1 becomes unstable and E * exists if α 2 > β 2 . If the Routh-Hurwitz criterion for the characteristics equation (3.6) is fulfilled then E * is asymptotically stable. Consequently, we have the following corollary about the occurrence of bifurcation in the system (1.3). (ii) If β 1 < 1, α 2 > β 2 , ϕ 1 > 0, ϕ 3 > 0 and ϕ 1 ϕ 2 > ϕ 3 then the system (1.3) undergoes a forward bifurcation at ω = ω from E 1 to E * .

Existence of Hopf-bifurcation
In this section, we study the Hopf-bifurcation around the coexistence equilibrium point E * = (u * 1 , u * 2 , u * 3 ) of the system (1.3). We consider ω and α 2 as the bifurcation parameters. ω = n 1 /K and α 2 = cv 0 /r are chosen as the bifurcation parameters because r and K are strongly related to the growth of immature prey, which controls energy input in the predator-prey system. Furthermore a, v 0 and n 1 are important parameters governing the exchange of energy from prey to predator.
Therefore the following equation is obtained .
According to Theorem 4.1, there exists a Hopf bifurcation in the stage-structure predator-prey model (1.3) where the Hopf bifurcation is controlled by ω. In fact, using the same argument as in the proof of Theorem 4.1, we can show that the Hopf bifurcation can also be controlled by parameter α 2 . The possibility of the Hopf bifurcation occurance is stated in the following theorem.

Numerical simulations
For graphical confirmation of the previous analytical results, several numerical simulations have been carried out. The model (1.3) is solved using the fourth-order Runge-Kutta method with some different initial conditions and some hypothetical values of the parameters. We first consider the following parameter values: α 1 = 0.5, α 2 = β 2 = 0.3, ω = 0.03 and β 1 = 1.01. Because β 1 > 1, Theorem 3.1 says that the extinction of population point (E 0 ) is the only equilibrium point which is locally asymptotically stable. This can be understood from the fact that in this case β 1 = δ 0 r > 1, i.e. the death rate of mature prey is larger than the intrinsic growth rate of immature prey. Hence the system will asymptotically be convergent to the extinction equilibrium point. The behavior of this case is depicted in Figure 1.(a). Next, we carry out a simulation using the same previous parameter values, except β 1 = 0.3. In addition to the unstable extinction of population point (E 0 ), system (1.3) also has the predator-free equilibrium E 1 = (0.2625, 0.4375, 0.0). Since α 2 = β 2 , we have ω = 0 and therefore ω > ω. Based on Theorem 3.1, E 1 is asymptotically stable. The case of α 2 = β 2 is equivalent to the case of cv 0 = δ 1 . Hence, from the third equation in the system (1.2), we have that the per capita growth rate of predator is less than the predator death rate. This explains the extinction of predator population. The numerical solutions which describe this situation is plotted in Figure 1.(b). To see a complete view of the dynamics of the system (1.3) with parameter values α 1 = 0.5, α 2 = β 2 = 0.3, and ω = 0.03, in Figure 2 we plot the equilibrium densities of immature and mature prey as function of β 1 . The equilibrium density of predator in this case is zero and thus it is not shown in the picture. Figure 1 shows that if β 1 < 1 then E 0 is unstable while E 1 is stable. On the other hand, when β 1 > 1, E 1 disappears and E 0 becomes stable. Thus, the system (1.3) undergoes a forward bifurcation.
We now choose α 1 = α 2 = 0.5 and β 1 = β 2 = 0.3. Using Theorem 4.1, we can check that the system undergoes a Hopf bifurcation around E * where the bifurcation point is at ω * = 0.03989. To verifiy the  Figure 3. It can be seen in this figure that our numerical solutions are convergent to a periodic solution when ω = 0.03 < ω * and convergent to E * when ω = 0.04 > ω * . This behavior indicates that the coexistence equilibrium point (E * ) changes its stability, i.e. E * is unstable and convergent to a limit cycle when ω < ω * and it becomes asymptotically stable if ω > ω * . In other words, the system (1.3) experiences a Hopf bifurcation driven by parameter ω. Figure 3 gives an indication that the Hopf bifurcation is supercritical because the limit cycle is stable, As stated in Theorem 4.2, the Hopf bifurcation can also be controlled by parameter α 2 . To show this behavior we perform simulation using α 1 = 0.5, β 1 = β 2 = 0.3 and ω = 0.03. The critical α * 2 for the occurance of Hopf bifurcation can be determined using Theorem 4.2, where in this case we get α * 2 = 0.5822. Using Theorem 3.1, it can be verified that the coexistence equilibrium point (E * ) is unstable for the case of α = 0.55, and it is asymptotically stable for α = 0.59. As depicted in Figure 4, such stability properties coincide with our numerical simulations. It is clearly seen in Figure 4.(a) that for α 2 = 0.55 < α * 2 , E * is unstable and the solutions converge to a limit cycle (periodic solution). If we take α 2 = 0.59 > α * 2 , then the solution is convergent to E * , see Figure 4.(b). From Figure 4, we also see the occurance of (supercritical) Hopf bifurcation controlled by α 2 .
To see the detail dynamics of the system (1.3) with parameter α 1 = 0.5, β 1 = 0.3 and β 2 = 0.3, we plot the stability region in (ω, α 2 )−plane, see Figure 5. In this picture we see that there are three different regions: the cyan region represents the stability area of E 1 (E 0 exists but it is unstable, E * does not exist); the green region represents the stability area of E * (E 0 and E 1 exists but they are not stable); and the yellow region corresponds to the stable limit cycle (all the three equilibrium points exists   but they are unstable). We identify the occurrence of a Hopf bifurcation as a change from the green region to the yellow region due to the changes of parameters ω or α 2 . We also notice the occurrence of a forward bifurcation caused by the changes of α 2 , see the changes from cyan region to green region. To provide a clearer illustration, we plot a bifurcation diagram for system (1.3) with parameter α 1 = 0.5, β 1 = β 2 = 0.3 and ω = 0.03, see Figure 6. From this figure we can see the appearance of transcritical bifurcation which is caused by the changing of the value of α 2 . When α 2 < α 1 2 = 0.3343, the equilibrium point E 1 is asymptotically stable. On the contrary, if α 1 2 < α 2 < α 2 2 = 0.3922 then equilibrium point E 1 becomes ustable and equilibrium point E * is asymtotically stable. In Figure 6 we also observe that the system (1.3) undergoes 2-Hopf bifurcation. In this case, the equilibrium point E * is asymptotically stable for α 1 2 < α 2 < α 2 2 = 0.3922 or α 2 > α * 2 = 0.5822, and it is unstable for α 2 2 < α 2 < α * 2 . In the latter case, the system (1.3) is convergent to a periodic solution or limit cycle.

Conclusion
A Rosenzweig-MacArthur predator-prey model with stage-structure in prey has been discussed. It was shown that the proposed model has three equilibrium points, i.e. the extinction of population (E 0 ), the predator-free point (E 1 ), and the coexistence point (E * ). All of these equilibrium points are conditionally asymptotically stable. Our analysis also showed that the proposed model exhibits a Hopfbifurcation which can be driven by (ω) or (α 2 ). Furthermore, the proposed model may also undergo a forward bifurcation for suitable parameter values. The analytical findings have been confirmed by our numerical simulations.