# Rich Dynamics of a Predator-Prey System with Different Kinds of Functional Responses

Authors: Kankan Sarkar; Subhas Khajanchi; Prakash Chandra Mali; Juan J. Nieto
Journal: Complexity
 (2020)

---

## Abstract

In this study, we investigate a mathematical model that describes the interactive dynamics of a predator-prey system with different kinds of response function. The positivity, boundedness, and uniform persistence of the system are established. We investigate the biologically feasible singular points and their stability analysis. We perform a comparative study by considering different kinds of functional responses, which suggest that the dynamical behavior of the system remains unaltered, but the position of the bifurcation points altered. Our model system undergoes Hopf bifurcation with respect to the growth rate of the prey population, which indicates that a periodic solution occurs around a fixed point. Also, we observed that our predator-prey system experiences transcritical bifurcation for the prey population growth rate. By using normal form theory and center manifold theorem, we investigate the direction and stability of Hopf bifurcation. The biological implications of the analytical and numerical findings are also discussed in this study.


---

## Body

## 1. Introduction

The interaction between prey species and predator species constitutes biological processes of major ecological importance that are ubiquitous in nature. The dynamical interaction between prey-predator models with different kinds of response function is a dominant theme in applied mathematics and theoretical ecology [1]. Mathematical modeling for prey-predator interplays has drawn attention to the mathematicians and scientists since the pioneering worked by Lotka and Volterra in the year 1920s, and there has been a large investigation for their rich dynamics [2–6]. When we go through the interaction between prey-predator system through mathematical modeling, initially it appears to be very straightforward, but at the end, it becomes very challenging tasks from the view point of mathematicians. The investigation and validation for the proper ecological behavior of the proposed model under consideration becomes most challenging and crucial. Since the complicated and rich dynamics for the interactive populations are usual throughout the world, many ecologists have studied the processes that influence the kinetics of predator-prey system and are interested to know which model would be best to represent the interaction among species.The dynamical behavior for the prey-predator system takes part in a significant role in theoretical ecology. There are many ecologists feeling that the unique nonnegative steady state for a predator-prey model is globally asymptotically stable if it becomes stable asymptotically although it is not always true. It can be shown that, under suitable condition, a unique nonnegative locally asymptotically stable steady state has at least one limit cycle around the steady state. Therefore, many mathematicians attempt to utilize some well-studied methods to obtain situations of global stability of the steady state for the prey-predator model [5–10].Generally, the prey-predator system obeys two basic principles: (i) prey species has natural birth rate and both the species has natural death rate; (ii) predator species can only grow by consuming the prey species [11]. In order to explain the predation phenomenon, Holling [12] has proposed different types of functional responses such as type-I, type-II, and type-III. The response functions are the rate of prey consumption by the predator per unit time. The functional responses proposed by Holling [12] are the functions of prey density only. The response functions are independent of predator density. One of the widely studied functional responses is Holling type-II response function, which is characterized by decelerating intake rate. The rate of predation rises as prey density increases, but after a certain stage, the rate of predation remains constant although prey density increases. In the nature, predator population death rate is likely to be maximum when the prey population is minimum due to the lack of food. Later, it is recognized that the response function is also dependent on predator density. In the year 1975, Beddington and DeAngelis [13, 14] introduce a response function, which is alike to Holling type-II response function but holds an additional term, which represents the mutual interaction between the predators. Hence, Beddington–DeAngelis-type response function is a function of both prey and predator density; that is, the rate of predation is function of both prey and predator density. However, depending upon the structural complexity of the prey habitat, the response function for the prey-predator interaction may be different [15].Recently, a series of mathematical models has been developed to better understand the dynamics of the predator-prey system with different kinds of functional responses [6, 16–24]. Dalziel et al. [16] studied a modified Holling type-II predator-prey model, based on the premise that the search rate of predators is dependent on the prey density, rather than constant. They performed a complete analysis of the global behavior of their model and the dichotomy properties. Wang et al. [17] investigated an improved predator-prey model that incorporates seasonality and a predator maturation delay, leading to a system of periodic differential equations with a time delay. They also defined the basic reproduction number R0 and show that it is a threshold parameter determining whether the predators can coexist with the prey species. Kohnke et al. [18] studied a predator-prey model and developed a functional response using timescale separation. In their study, they showed that the resulting functional response contains a single parameter that controls whether the group defense functional response is saturating or dome-shaped. Their model exhibits bifurcation analysis. Antwi-Fordjour et al. [19] investigated the interplay of the mutual interference exponent, and prey refuge, on the behavior of a predator-prey model with a generalized Holling-type functional response considering in particular the “nonsmooth” case. The authors studied the dynamical properties of their model and derived conditions for the occurrence of saddle-node, transcritical, and Hopf bifurcations. Also, they derived the sufficient condition for finite time extinction of the prey population.When we look through literature review for the published research papers on the dynamics of the prey-predator model, we notice that there are multiple numbers of good research works on the dynamics of prey and predators with different kinds of functional responses. Some researchers investigated and elevated some open quires of various kinds of response function. Seo and Kot [2] have compared two prey-predator models with Holling type-I response function. In [25], the author investigated a predator-prey model with different kinds of functional responses by incorporating the role of constant prey refuge. Aiello and Freedman [26] studied a single species model with the effect of discrete time lag, in which the species was divided into mature and immature. Their model demonstrates nonnegative equilibria as the global attractor. Kon et al. [27] established the criteria for the permanence of stage structure prey-predator interactions. These are some existing articles which focused on Beddington–DeAngelis-type functional responses [28–31]. Beddington–DeAngelis-type response function confesses plentiful and biologically relevant prey-predator dynamics [29], which attract researchers to further study the predator-prey model including Beddington–DeAngelis response function. Khajanchi [28] studied a predator-prey model with age-structured incorporating Beddington–DeAngelis-type response function, in which the author performed global stability analysis. Xie et al. [20] studied a predator-prey system by using a system of fractional order differential equations with Holling type-III response function and incorporate the effect of discontinuous harvesting. Predator-prey system has been investigated by introducing discrete time delay with disease in the predator by Huang et al. [21]. Also, nonautonomous Leslie–Gower predator-prey system with nonlinear prey harvesting was studied by Alam et al. [22]. Sarkar and Khajanchi [23] studied a predator-prey system by incorporating the effect of fear with Holling type-II response function, where the authors preformed the bifurcation analysis of the system to better understand the dynamics of the model. Seasonally, varying prey-predator system with the Allee effect model was investigated by Rebelo and Soresina [24], and they studied the coexistence of all the species. Time delay plays an important role in stability switching and may change drastically the qualitative behavior of solutions for the model system under consideration [32–35].The main aim of this manuscript is to study a very simple two-tier continuous predator-prey model with different kinds of functional responses. Although the model is very simple, it gives rich and biologically meaningful dynamics. A comparative study of the stability property was done by considering Holling type-I, Holling type-II, and Beddington–DeAngelis-type response functions. Our predator-prey system experiences Hopf bifurcation with respect to the prey birth rater0, which gives more complicated dynamics of the system. The predator-prey system also experiences transcritical bifurcation with respect to the prey birth rate r0. The biological implications of the theoretical and numerical results are also investigated.This article is organized in the following way: we present a predator-prey model with different kinds of functional responses in the next section. In Section3, we discuss the qualitative behavior for the model, and positivity for the model, boundedness for the solutions, and the position of all possible feasible equilibria and the existence criteria of the equilibria are presented. The criteria for stability of the predator-prey system around those equilibria are studied. In the same section, we performed a comparative study by considering different kinds of functional responses. We performed global stability analysis by constructing suitable Lyapunov function. By using normal form theory and center manifold theorem, we investigate direction and stability of Hopf bifurcation in Section 4. In the next section, we give some numerical simulations to support our analytical findings. Section 6 presents the discussion about key findings for our problem.


## 2. The Model

The interaction between prey and predator species in simplest form is represented by the following system of coupled differential equations [2, 36]:(1)dudt=ufu−gu,vv,dvdt=v−δ2+θgu,v,where vt and ut represent the predator and prey populations, respectively, at any time t and fu stands the birth rate of prey species in exclusion of any predator population as well as death rate and the intraspecific competition between the prey species. Here, δ2 denotes the natural death rate for predator species and θ stands the conversion coefficient of prey’s biomass to predator’s biomass, and gu,v indicates the prey-dependent response function.Our proposed model consists of two population, namely, prey population and predator population. Concentration of prey and predator species at timet can be represented by ut and vt, respectively. Both the prey and predator populations are continuous in time t. The predator-prey interaction with response function gu,v can be expressed by the system of coupled differential equations:(2)dudt=r0u−δ1u−γu2−gu,vv,dvdt=v−δ2+θgu,v.With initial valuesu0=u0≥0 and v0=v0≥0, all the model parameters r0, δ1, γ, δ2, and θ are positive constants and can be interpreted as follows: r0 represents the birth rate of prey species in absence of predators, δ1 represents the death rate for prey population, γ denotes the intraspecific competition of prey species, δ2 represents the death rate for predator population, and θ represents the conversion coefficient of prey’s biomass to predator’s biomass. The dimension of prey and predator biomass is Mass. Also, the dimensions of r0, δ1, γ, and δ2 are Time−1, Time−1, Mass−1⋅Time−1, and Time−1, respectively, and the parameter θ is dimensionless. In exclusion of predator species, the prey population obeys the logistic growth with an intrinsic growth rate r0−δ1>0.


### 2.1. Functional Response

In ecology, the functional responses are the intake rate of predator population as a function of prey population. Holling [12] proposed three types of response function, namely, Holling types-I, II, and III. In our model formulation, we consider Holling Type-II response function in the following form:(3)gu,v=αu1+βu,where α represents the attack rate and β is the handling time to capture the prey population; both α and β are positive constants. The dimensions of α and β are Mass−1Time−1 and Mass−1, respectively. The predation rate increases as prey density increases, but after a certain stage, rate of predation becomes constant although prey density increases. The functional response depends on prey density but not on predator density. By using the functional form of gu,v given by (3), system (2) becomes(4)dudt=r0u−δ1u−γu2−αuv1+βu≡F1u,v,dvdt=αθuv1+βu−δ2v≡F2u,v.


## 2.1. Functional Response

In ecology, the functional responses are the intake rate of predator population as a function of prey population. Holling [12] proposed three types of response function, namely, Holling types-I, II, and III. In our model formulation, we consider Holling Type-II response function in the following form:(3)gu,v=αu1+βu,where α represents the attack rate and β is the handling time to capture the prey population; both α and β are positive constants. The dimensions of α and β are Mass−1Time−1 and Mass−1, respectively. The predation rate increases as prey density increases, but after a certain stage, rate of predation becomes constant although prey density increases. The functional response depends on prey density but not on predator density. By using the functional form of gu,v given by (3), system (2) becomes(4)dudt=r0u−δ1u−γu2−αuv1+βu≡F1u,v,dvdt=αθuv1+βu−δ2v≡F2u,v.


## 3. Qualitative Behavior of the Model

### 3.1. Positivity, Boundedness, and Uniform Persistence

Right side of the system (4) is a continuous function for dependent variable t. Hence, after integration on both sides of each equation of system (4), we get(5)ut=u0exp∫0tr0−δ1−γus−αvs1+βusds>0,∀t≥0,vt=v0exp∫0tαθus1+βus−δ2ds>0,∀t≥0.From the above expressions ofut and vt, it is clear that both ut and vt remain nonnegative for all finite time, that is, t≥0. Hence, ℜ2+=ut,vt:ut≥0,vt≥0 is positively invariant of system (4). In order to prove the boundedness of the prey-predator model (4), we consider the function(6)ϕt=ut+1θvt.Differentiation with respect tot gives(7)dϕtdt=dutdt+1θdvtdt=r0−δ1ut−γu2−δ2θvt.Now, for eachϕt>0, we get(8)ϕtdt+ρϕt≤r0u−γu2,where ρ=minδ1,δ2/θ. The maximum value of r0u−γu2 is r02/4γ. Therefore,(9)ϕtdt+ρϕt≤r024γ=msay.Due to the theory of differential inequality, we have(10)0≤ϕt≤mρ1−e−ρt+ϕ0e−ρtfor t⟶∞, and we have 0≤ϕt≤m/ρ. Therefore, ϕt is bounded. Thus, all the solution of (4) are confined in(11)Ω=u,v:0≤ut+1θvt≤mρ+ɛ,∀ɛ≥0.This implies that the solution of system (4) is bounded. Thus, we have the following theorem.Theorem 1.
For any initial valueu0,v0∈ℜ2+, all the solutions for (4) are nonnegative and bounded.From the ecological perspective, the positivity of the prey-predator system (4) interprets that prey and the predators biologically well-behaved. The boundedness for the solutions of the prey-predator model (4) implies that neither the prey population nor the predator population will grow unboundedly or abruptly for a large period.Definition 1.
System (4) is said to be uniformly persistent if there exist two constants κ and τ, where 0≤κ≤τ such that maxlimsupt⟶∞ut,limsupt⟶∞vt≤τ and minliminft⟶∞ut,liminft⟶∞vt≥κ, where κ and τ are independent of the initial conditions.Theorem 2.
System (4) is uniformly persistent if r0>δ1+αv¯, where v¯ is the upper bound of vt.Proof.
From the proof of Theorem1, we have limsupt⟶∞ut≤m/ρ and limsupt⟶∞vt≤mθ/ρ. We choose τ=maxm/ρ,mθ/ρ. Hence, there exists τ>0 such that maxlimsupt⟶∞ut,limsupt⟶∞vt≤τ.
Now, from the first equation of system (4), we have(12)dudt=r0u−δ1u−γu2−αuv1+βu≥r0u−δ1u−γu2−αuv¯=r0−δ1−αv¯u−γu2,where v¯ is the upper limit of v in Ω.
Thus, we can conclude thatliminft⟶∞ut≥r0−δ1−αv¯/γ with r0>δ1+αv¯. Again, from the second equation of (4), we have(13)dvdt=αθuv1+βu−δ2v≥−δ2v.
Hence, we haveliminft⟶∞vt≥0. Thus, if we choose κ=min0,r0−δ1−αv¯/γ=0 with r0>δ1+αv¯, then minliminft⟶∞ut,liminft⟶∞vt≥κ. Therefore, system (4) will be uniformly persistent if r0>δ1+αv¯.


### 3.2. Equilibria and Their Existence

The singular points for the model (4) are nonnegative solutions for F1u,v=0 and F2u,v=0 in ℜ+2; that is, the equilibrium points of system (4) are the positive point of intersection for zero growth of prey and predator isoclines. The prey isoclines are given by u=0, and v=r0−δ1−γu1+βu/α with predator isoclines is stated as v=0 and the vertical straight line u=δ2/αθ−βδ2. System (4) possesses three biologically meaningful singular points, as follows:(i)
Trivial steady stateE00,0.(ii)
Boundary singular pointE1u¯,0, where u¯=r0−δ1/γ, which is feasible only when r0>δ1; that is, growth of prey population is higher than the death rate for prey population.(iii)
The positive interior equilibrium pointE∗u∗,v∗, where u∗=δ2/αθ−βδ2 and v∗=r0−δ1−γu∗1/α+β/αu∗.The shape and position of the curvev=r0−δ1−γu1+βu/α depends upon multiple parameters. For different values of β, the zero growth prey isoclines are convex curve joining the boundary equilibrium E1u¯,0 and 0,r0−δ1/α (see Figure 1(a)). As β increases, the vertical predator isocline shifted towards increasing u. For different values of α, zero growth prey isocline becomes a convex curve, which meets at boundary equilibrium E1u¯,0 (see Figure 1(b)). As α increases, the vertical predator isocline shifted towards decreasing u and the prey isocline shifted towards decreasing v. Also, for different values of r0, the zero growth prey isoclines are convex curve, but the zero growth predator isoclines do not change position. As r0 increases, the zero growth prey isocline shifted towards v increasing. In all these situations, the interior equilibrium point E∗u∗,v∗ exists only if αθ−βδ2>0 (see Figure 1(c)).Figure 1
Isoclines of the predator-prey system (4) for three different cases with the following set of parameter values: r0=0.04,δ1=0.0125,γ=0.01,α=0.8,β=0.5,θ=0.4, and δ2=0.05. (a) For different β; (b) for different α; (c) for different values of r0, as shown in the inset.
(a)(b)(c)


### 3.3. Local Stability Analysis

In this section, we analyzed the linear stability for the predator-prey system (4). Using linearization techniques, the local behavior of the nonlinear predator-prey system can be investigated. Generally, the predator-prey system is linearized around every biologically feasible singular point, and the model is perturbed by a small quantity and observed whether the model returns to that state of singular point converges to another state of singular point. Analysis of local stability leads to investigate qualitative dynamics for the model system under consideration. For local stability analysis of the considered model (4) around every singular point, first calculate the Jacobian matrix obtained from each of the equilibrium point. The variational matrix of the model (4) at u,v is(14)J=FuFvGuGv,where(15)Fu=r0−δ1−2γu−αv1+βu2,Fv=−αu1+βu,Gu=αθv1+βu2,Gv=αθu1+βu−δ2.Theorem 3.
Ifr0<δ1, then the trivial singular point E00,0 of (4) is locally asymptotically stable and otherwise unstable.Proof.
The eigenvalues of the Jacobian matrix (14) of the system (4) at the trivial equilibrium point E00,0 are r0−δ1 and −δ2. Here, δ2>0; therefore, the trivial equilibrium point E00,0 of the model (4) will be locally asymptotically stable if r0<δ1; that is, the growth rate of prey species is less than the natural death rate of prey species. It can be noted that if the trivial equilibrium point E00,0 is asymptotically stable, then the boundary singular point E1u¯,0 does not exist. For r0=δ1, the trivial equilibrium point E00,0 has eigenvalues 0 and −δ2. So E00,0 is a nonhyperbolic equilibrium point [6].Theorem 4.
Ifδ1<r0<δ1+γδ2/αθ−βδ2, then the boundary equilibrium point E1u¯,0 of the model (4) is locally asymptotically stable and otherwise unstable.Proof.
The eigenvalues of the Jacobian matrix (14) of the system (4) at the equilibrium point E1u¯,0 are λ1=−r0−δ1 and λ2=αθr0−δ1/γ+βr0−δ1−δ2. The boundary equilibrium point E1u¯,0 will be locally asymptotically stable if λ1<0 and λ2<0, and this condition implies δ1<r0<δ1+γδ2/αθ−βδ2. Biologically, r0>δ1 implies that the rate of growth of prey population is higher than the natural rate of death of the prey. Also, r0<δ1+γδ2/αθ−βδ2 implies that the prey growth rate cannot be unbounded. It is to be noted that if the boundary equilibrium state E1u¯,0 is locally asymptotically stable, then the trivial equilibrium state E00,0 is unstable.Theorem 5.
Ifr0−δ1−2γu∗−αv∗/1+βu∗2<0, then the positive interior singular point E∗u∗,v∗ of (4) will be asymptotically stable and otherwise unstable.Proof.
The variational matrix around the coexisting singular pointE∗u∗,v∗ for the system (4) is given by(16)J∗=r0−δ1−2γu∗−αv∗1+βu∗2−αu∗1+βu∗αθv∗1+βu∗20.
The characteristic equation of the above matrixJ∗ can be expressed as(17)λ2−r0−δ1−2γu∗−αv∗1+βu∗2λ+αδ2v∗1+βu∗2=0,orλ2−K11λ+K22=0,where K11=r0−δ1−2γu∗−αv∗/1+βu∗2 and K22=αδ2v∗/1+βu∗2.
Due to well-known Routh-Hurwitz criterion, the characteristic (17) will have negative real roots if K11<0 and K22>0. Hence, r0−δ1−2γu∗−αv∗/1+βu∗2<0 and αδ2v∗/1+βu∗2>0⟹v∗>0 since α,δ2,1+βu∗2>0. It is obvious that v∗>0 because E∗u∗,v∗ will be positive interior singular point of (4). Therefore, the positive interior steady state E∗u∗,v∗ of (4) will be asymptotically stable if r0−δ1−2γu∗−αv∗/1+βu∗2<0 and otherwise unstable.
Figure2 shows the stability region of different singular points of (4) in r0−β parameter-space. The growth rate of prey species r0 varies from 0 to 0.1, and the half saturation constant β varies from 0 to 1; the other parameters are fixed as δ1=0.0125, γ=0.01, α=0.8, θ=0.4, and δ2=0.05. For a very small value of r0, that is, when r0<δ1, the trivial equilibrium E0 is locally asymptotically stable. The stability region of E0 is shown by blue-shaded region in the r0−β plane. In this blue-shaded region, the boundary equilibrium E1 and the interior equilibrium E∗ are unstable. The boundary equilibrium E1 is locally asymptotically stable in the green-shaded region. But the trivial equilibrium E0 and the interior equilibrium E∗ are unstable in this green-shaded region. The red-shaded region of Figure 2 shows the locally asymptotically stable behavior of the interior equilibrium E∗; that is, both the prey and predator populations can exist together in the red-shaded region. For higher value of r0 (such as r0>0.06) and for higher value of β (such as β>0.5), the coexisting equilibrium E∗ loses stability. The trivial equilibrium E0 and the boundary equilibrium E1 are unstable in the red-shaded region. All the equilibria lose stability in the white-shaded region.Figure 2
The stability region of the trivial singular stateE0, boundary singular state E1, and the coexisting singular state E∗ of (4), in the r0−β parameter-space where r0∈0,0.1 and β∈0,1.0. Other parameter values are δ1=0.0125, γ=0.01, α=0.8, θ=0.4, and δ2=0.05. In the figure, (i) within the blue colored region, the equilibrium point E0 is locally asymptotically stable, but boundary singular point E1 and the interior equilibrium E∗ are unstable; (ii) within the green colored region, the equilibrium point E1 is locally asymptotically stable, but trivial equilibria E0 and interior equilibrium E∗ are unstable; and (iii) within the red colored region, the interior equilibrium point E∗ is locally asymptotically stable, but the trivial equilibrium E0 and boundary equilibrium E1 are unstable.


### 3.4. Comparative Study with Different Kinds of Functional Responses

Mainly, three types of response functions (Holling types-I, II, and III) and Beddington–DeAngelis are widely investigated by the researchers in the predator-prey system [29, 37, 38]. Here, we also investigate the role of different kinds of functional responses on the dynamics for the predator-prey system (2). To exemplify this matter, we performed a relative study for (2) with reference to three kinds of response functions, namely, Beddington–DeAngelis, Holling type-I, and Holling type-II (already used in the model (4)). Holling-type response functions are the function of prey population only; that is, the response functions are independent of predator population. But, in the nature, response functions are also the function of predator population. Widely studied Beddington–DeAngelis [13, 14] response function considered the density of predators also. In our study, we have also studied stability of the system by considering Beddington–DeAngelis functional responses and two types of response function. The Beddington–DeAngelis response function is of the following form:(18)gu,v=αu1+βu+ξv,where the parameters α,β, and ξ are nonnegative and have their usual biological meanings. The Beddington–DeAngelis response function is a similar kind of Holling type-II response function but contains an additional term ξv at the denominator. The term ξv indicates the mutual interaction between the predator population. Thus, when the prey density ut is constant and predator density vt is minimum, then the value of the response function gu,v is maximum; when predator density vt is maximum, then the value of the response function gu,v is minimum. This is more realistic because when the number of predator species are less than their prey, consumption due to predators per unit time will be more, but when predator density is more, then the prey consumption due to predators per unit time will be less due to mutual interaction of predator population. But, when ξ=0, then the Beddington–DeAngelis functional response is the same as Holling type-II functional response. Again, for β=0andξ=0, then Beddington–DeAngelis response function is the same as Holling type-I response function. Hence, we may conclude that Beddington–DeAngelis response function is a more generalized form of Holling type-II and type-I functional responses. By considering Beddington–DeAngelis response function, the system (2) becomes(19)dudt=r0u−δ1u−γu2−αuv1+βu+ξv,dvdt=αθuv1+βu+ξv−δ2v.The predator-prey system (19) with Beddington–DeAngelis functional response has three ecologically relevant singular points, namely, trivial singular point E00,0, boundary singular point E1u^,0, and an interior steady state E∗u^∗,v^∗. Here, u^=r0−δ1/γ, and u^∗ is the positive root for the following second-degree equation:(20)θγξu^∗2−r0−δ1θξ−αθ−βδ2u^∗−δ2=0,v^∗=1ξδ2αθ−βδ2u^∗−δ2.We have made a relative investigation of the predator-prey model (2), by considering three different functional responses, which are type-I (model (19) with β=0,ξ=0), Holling type-II (already used in model (4) with ξ=0), and Beddington–DeAngelis. We are not showing detailed analysis of the model (2) for the different functional responses; rather, we have presented the stability scenarios for these three types of response function in Table 1. The parameters r0, θ, and γ are varied to observed the stability behavior of the system. The detailed numerical study and graphical illustrations have been shown in numerical Section 5.Table 1
Stability conditions for different equilibrium points of system (19).
EquilibriaHolling type-I (β=0,ξ=0)Holling type-II (ξ=0)Beddington–DeAngelisE00,0r0<δ1r0<δ1r0<δ1E1u^,0r0>δ1 and αθr0−δ1<δ2r0>δ1 and αθγr0−δ1/γ+βr0−δ12<δ2r0>δ1 and αθγr0−δ1/γ+βr0−δ12<δ2E∗u^∗,v^∗M1+M4<0 and M1M4−M2M3>0 with β=0andξ=0M1+M4<0 and M1M4−M2M3>0 with ξ=0M1+M4<0 and M1M4−M2M3>0M1,M2,M3, and M4 are the elements of the variational matrix around interior singular point E∗u^∗,v^∗, where M1=−γu^∗−αβu^∗v^∗/1+βu^∗+ξv^∗2, M2=−αu^∗1−ξv^∗/1+βu^∗+ξv^∗2, M3=αθv^∗1−βu^∗/1+βu^∗+ξv^∗2, and M4=αθu^∗1−ξv^∗/1+βu^∗+ξv^∗2−δ2.


### 3.5. Global Stability Analysis aroundE∗u∗,v∗

In this subsection, we investigate the sufficient condition for global asymptotic stability of the positive coexisting singular pointE∗u∗,v∗ by establishing appropriate Lyapunov function.Theorem 6.
The positive interior steady stateE∗u∗,v∗ of the system (4) is globally asymptotically stable if γ/αβ>v∗.Proof.
The interior steady stateE∗u∗,v∗ of the prey-predator model (4) fulfills the following conditions:(21)r0−δ1=γu∗+αv∗1+βu∗,δ2=αθu∗1+βu∗.
With the help of these above two equations, the model (4) leads to(22)dudt=u−γu−u∗+αv∗1+βu∗−αv1+βu,dvdt=vαθu1+βu−αθu∗1+βu∗.
We define the Lyapunov functionalVu,v:ℜ2+⟶ℜ such that(23)Vu,v=V1u,v+V2u,v,where V1u,v=u−u∗−u∗lnu/u∗ and V2u,v=Kv−v∗−v∗lnv/v∗. Here, K is a positive constant and is defined as K=1/θ1+βu∗. This particular type of Lyapunov function to study the stability of the interior equilibrium has been used widely [29, 30, 39, 40]. The Lyapunov function Vu,v is continuous on Intℜ2+ and vanishes at E∗u∗,v∗. Also, Vu,v is positive definite in Intℜ2+. Furthermore, ∂V1/∂u>0 when u>u∗, ∂V1/∂u<0 when u<u∗, ∂V2/∂u>0 when v>v∗, and ∂V2/∂v<0 when v<v∗. Hence, E∗u∗,v∗ is a global minimum of Vu,v. The time derivative of V1u,v and V2u,v along the solution for (4) is(24)dV1dt=u−u∗ududt=u−u∗−γu−u∗+αv∗1+βu∗−αv1+βu,dV2dt=Kv−v∗vdvdt=Kαθv−v∗u1+βu−u∗1+βu∗.
Now, differentiating (23) with respect to time t and using the values of dV1/dt and dV2/dt, we obtain(25)dVdt=dV1dt+dV2dt=−γu−u∗2−αu−u∗v−v∗+βu∗v−uv∗1+βu1+βu∗+Kαθu−u∗v−v∗1+βu1+βu∗.
Using the relationu∗v−uv∗=u∗v−v∗−v∗u−u∗, we obtain(26)dVdt=−γu−u∗2+αβv∗u−u∗21+βu1+βu∗−α1+βu∗−Kθ1+βu1+βu∗u−u∗v−v∗.
Putting the value ofK, we obtain(27)dVdt=−γu−u∗2+αβv∗u−u∗21+βu1+βu∗≤−γ−αβv∗u−u∗2.
Hence,dV/dt<0 along the trajectories in ℜ2+ excluding at E∗u∗,v∗ if γ−αβv∗>0, that is, if γ/αβ>v∗. Thus, if γ/αβ>v∗, then the coexisting singular point E∗u∗,v∗ is globally asymptotically stable.


### 3.6. Hopf Bifurcation Analysis

In our study, the growth rater0 of prey population is one of the most important parameters. In this section, we investigate how predator-prey model behavior alters with respect to r0 by using the analysis of Hopf bifurcation. We start with a theorem to study the analysis of Hopf bifurcation.Theorem 7.
When the prey birth rater0 crosses a threshold value r0=r0∗=δ1+2γu∗1+1/βu∗−1/βu∗δ1+γu∗, the predator-prey model (4) experiences Hopf bifurcation around the interior singular point E∗. Necessary and sufficient conditions of Hopf bifurcation to occur are that there exists r0=r0∗ such that the transversality condition holds:(28)ddr0Reλr0r0=r0∗=βu∗1+βu∗≠0.Proof.
A necessary condition for the change in dynamics of the coexisting singular pointE∗u∗,v∗ of the prey-predator model (4) is that the characteristic polynomial (17) should have purely complex eigenvalues. The characteristic roots of the equation (17) will have purely complex roots if and only if K11=0 and K22>0. For K11=0 gives, r0=r0∗=δ1+2γu∗1+1/βu∗−1/βu∗δ1+γu∗, where u∗=δ2/αθ−βδ2. Due to Theorem 5, we can conclude that the positive coexisting singular point E∗u∗,v∗ of the predator-prey model (4) will be stable asymptotically if r0−δ1−2γu∗−αv∗/1+βu∗2<0 and the model (4) experiences Hopf bifurcation at r0=r0∗.
Let the roots of the characteristic (17) are of the form λ1,2=ϕr0±iψr0. Then, λ1+λ2=2ϕr0=K11; that is, Reλr0=ϕr0=1/2K11. To prove the transversality condition, we need to verify the following condition:(29)ddr0Reλr0r0=r0∗=βu∗1+βu∗≠0,ifβ≠0.
Hence, the transversality condition holds, and the system (4) experiences Hopf bifurcation at r0=r0∗.


## 3.1. Positivity, Boundedness, and Uniform Persistence

Right side of the system (4) is a continuous function for dependent variable t. Hence, after integration on both sides of each equation of system (4), we get(5)ut=u0exp∫0tr0−δ1−γus−αvs1+βusds>0,∀t≥0,vt=v0exp∫0tαθus1+βus−δ2ds>0,∀t≥0.From the above expressions ofut and vt, it is clear that both ut and vt remain nonnegative for all finite time, that is, t≥0. Hence, ℜ2+=ut,vt:ut≥0,vt≥0 is positively invariant of system (4). In order to prove the boundedness of the prey-predator model (4), we consider the function(6)ϕt=ut+1θvt.Differentiation with respect tot gives(7)dϕtdt=dutdt+1θdvtdt=r0−δ1ut−γu2−δ2θvt.Now, for eachϕt>0, we get(8)ϕtdt+ρϕt≤r0u−γu2,where ρ=minδ1,δ2/θ. The maximum value of r0u−γu2 is r02/4γ. Therefore,(9)ϕtdt+ρϕt≤r024γ=msay.Due to the theory of differential inequality, we have(10)0≤ϕt≤mρ1−e−ρt+ϕ0e−ρtfor t⟶∞, and we have 0≤ϕt≤m/ρ. Therefore, ϕt is bounded. Thus, all the solution of (4) are confined in(11)Ω=u,v:0≤ut+1θvt≤mρ+ɛ,∀ɛ≥0.This implies that the solution of system (4) is bounded. Thus, we have the following theorem.Theorem 1.
For any initial valueu0,v0∈ℜ2+, all the solutions for (4) are nonnegative and bounded.From the ecological perspective, the positivity of the prey-predator system (4) interprets that prey and the predators biologically well-behaved. The boundedness for the solutions of the prey-predator model (4) implies that neither the prey population nor the predator population will grow unboundedly or abruptly for a large period.Definition 1.
System (4) is said to be uniformly persistent if there exist two constants κ and τ, where 0≤κ≤τ such that maxlimsupt⟶∞ut,limsupt⟶∞vt≤τ and minliminft⟶∞ut,liminft⟶∞vt≥κ, where κ and τ are independent of the initial conditions.Theorem 2.
System (4) is uniformly persistent if r0>δ1+αv¯, where v¯ is the upper bound of vt.Proof.
From the proof of Theorem1, we have limsupt⟶∞ut≤m/ρ and limsupt⟶∞vt≤mθ/ρ. We choose τ=maxm/ρ,mθ/ρ. Hence, there exists τ>0 such that maxlimsupt⟶∞ut,limsupt⟶∞vt≤τ.
Now, from the first equation of system (4), we have(12)dudt=r0u−δ1u−γu2−αuv1+βu≥r0u−δ1u−γu2−αuv¯=r0−δ1−αv¯u−γu2,where v¯ is the upper limit of v in Ω.
Thus, we can conclude thatliminft⟶∞ut≥r0−δ1−αv¯/γ with r0>δ1+αv¯. Again, from the second equation of (4), we have(13)dvdt=αθuv1+βu−δ2v≥−δ2v.
Hence, we haveliminft⟶∞vt≥0. Thus, if we choose κ=min0,r0−δ1−αv¯/γ=0 with r0>δ1+αv¯, then minliminft⟶∞ut,liminft⟶∞vt≥κ. Therefore, system (4) will be uniformly persistent if r0>δ1+αv¯.


## 3.2. Equilibria and Their Existence

The singular points for the model (4) are nonnegative solutions for F1u,v=0 and F2u,v=0 in ℜ+2; that is, the equilibrium points of system (4) are the positive point of intersection for zero growth of prey and predator isoclines. The prey isoclines are given by u=0, and v=r0−δ1−γu1+βu/α with predator isoclines is stated as v=0 and the vertical straight line u=δ2/αθ−βδ2. System (4) possesses three biologically meaningful singular points, as follows:(i)
Trivial steady stateE00,0.(ii)
Boundary singular pointE1u¯,0, where u¯=r0−δ1/γ, which is feasible only when r0>δ1; that is, growth of prey population is higher than the death rate for prey population.(iii)
The positive interior equilibrium pointE∗u∗,v∗, where u∗=δ2/αθ−βδ2 and v∗=r0−δ1−γu∗1/α+β/αu∗.The shape and position of the curvev=r0−δ1−γu1+βu/α depends upon multiple parameters. For different values of β, the zero growth prey isoclines are convex curve joining the boundary equilibrium E1u¯,0 and 0,r0−δ1/α (see Figure 1(a)). As β increases, the vertical predator isocline shifted towards increasing u. For different values of α, zero growth prey isocline becomes a convex curve, which meets at boundary equilibrium E1u¯,0 (see Figure 1(b)). As α increases, the vertical predator isocline shifted towards decreasing u and the prey isocline shifted towards decreasing v. Also, for different values of r0, the zero growth prey isoclines are convex curve, but the zero growth predator isoclines do not change position. As r0 increases, the zero growth prey isocline shifted towards v increasing. In all these situations, the interior equilibrium point E∗u∗,v∗ exists only if αθ−βδ2>0 (see Figure 1(c)).Figure 1
Isoclines of the predator-prey system (4) for three different cases with the following set of parameter values: r0=0.04,δ1=0.0125,γ=0.01,α=0.8,β=0.5,θ=0.4, and δ2=0.05. (a) For different β; (b) for different α; (c) for different values of r0, as shown in the inset.
(a)(b)(c)


## 3.3. Local Stability Analysis

In this section, we analyzed the linear stability for the predator-prey system (4). Using linearization techniques, the local behavior of the nonlinear predator-prey system can be investigated. Generally, the predator-prey system is linearized around every biologically feasible singular point, and the model is perturbed by a small quantity and observed whether the model returns to that state of singular point converges to another state of singular point. Analysis of local stability leads to investigate qualitative dynamics for the model system under consideration. For local stability analysis of the considered model (4) around every singular point, first calculate the Jacobian matrix obtained from each of the equilibrium point. The variational matrix of the model (4) at u,v is(14)J=FuFvGuGv,where(15)Fu=r0−δ1−2γu−αv1+βu2,Fv=−αu1+βu,Gu=αθv1+βu2,Gv=αθu1+βu−δ2.Theorem 3.
Ifr0<δ1, then the trivial singular point E00,0 of (4) is locally asymptotically stable and otherwise unstable.Proof.
The eigenvalues of the Jacobian matrix (14) of the system (4) at the trivial equilibrium point E00,0 are r0−δ1 and −δ2. Here, δ2>0; therefore, the trivial equilibrium point E00,0 of the model (4) will be locally asymptotically stable if r0<δ1; that is, the growth rate of prey species is less than the natural death rate of prey species. It can be noted that if the trivial equilibrium point E00,0 is asymptotically stable, then the boundary singular point E1u¯,0 does not exist. For r0=δ1, the trivial equilibrium point E00,0 has eigenvalues 0 and −δ2. So E00,0 is a nonhyperbolic equilibrium point [6].Theorem 4.
Ifδ1<r0<δ1+γδ2/αθ−βδ2, then the boundary equilibrium point E1u¯,0 of the model (4) is locally asymptotically stable and otherwise unstable.Proof.
The eigenvalues of the Jacobian matrix (14) of the system (4) at the equilibrium point E1u¯,0 are λ1=−r0−δ1 and λ2=αθr0−δ1/γ+βr0−δ1−δ2. The boundary equilibrium point E1u¯,0 will be locally asymptotically stable if λ1<0 and λ2<0, and this condition implies δ1<r0<δ1+γδ2/αθ−βδ2. Biologically, r0>δ1 implies that the rate of growth of prey population is higher than the natural rate of death of the prey. Also, r0<δ1+γδ2/αθ−βδ2 implies that the prey growth rate cannot be unbounded. It is to be noted that if the boundary equilibrium state E1u¯,0 is locally asymptotically stable, then the trivial equilibrium state E00,0 is unstable.Theorem 5.
Ifr0−δ1−2γu∗−αv∗/1+βu∗2<0, then the positive interior singular point E∗u∗,v∗ of (4) will be asymptotically stable and otherwise unstable.Proof.
The variational matrix around the coexisting singular pointE∗u∗,v∗ for the system (4) is given by(16)J∗=r0−δ1−2γu∗−αv∗1+βu∗2−αu∗1+βu∗αθv∗1+βu∗20.
The characteristic equation of the above matrixJ∗ can be expressed as(17)λ2−r0−δ1−2γu∗−αv∗1+βu∗2λ+αδ2v∗1+βu∗2=0,orλ2−K11λ+K22=0,where K11=r0−δ1−2γu∗−αv∗/1+βu∗2 and K22=αδ2v∗/1+βu∗2.
Due to well-known Routh-Hurwitz criterion, the characteristic (17) will have negative real roots if K11<0 and K22>0. Hence, r0−δ1−2γu∗−αv∗/1+βu∗2<0 and αδ2v∗/1+βu∗2>0⟹v∗>0 since α,δ2,1+βu∗2>0. It is obvious that v∗>0 because E∗u∗,v∗ will be positive interior singular point of (4). Therefore, the positive interior steady state E∗u∗,v∗ of (4) will be asymptotically stable if r0−δ1−2γu∗−αv∗/1+βu∗2<0 and otherwise unstable.
Figure2 shows the stability region of different singular points of (4) in r0−β parameter-space. The growth rate of prey species r0 varies from 0 to 0.1, and the half saturation constant β varies from 0 to 1; the other parameters are fixed as δ1=0.0125, γ=0.01, α=0.8, θ=0.4, and δ2=0.05. For a very small value of r0, that is, when r0<δ1, the trivial equilibrium E0 is locally asymptotically stable. The stability region of E0 is shown by blue-shaded region in the r0−β plane. In this blue-shaded region, the boundary equilibrium E1 and the interior equilibrium E∗ are unstable. The boundary equilibrium E1 is locally asymptotically stable in the green-shaded region. But the trivial equilibrium E0 and the interior equilibrium E∗ are unstable in this green-shaded region. The red-shaded region of Figure 2 shows the locally asymptotically stable behavior of the interior equilibrium E∗; that is, both the prey and predator populations can exist together in the red-shaded region. For higher value of r0 (such as r0>0.06) and for higher value of β (such as β>0.5), the coexisting equilibrium E∗ loses stability. The trivial equilibrium E0 and the boundary equilibrium E1 are unstable in the red-shaded region. All the equilibria lose stability in the white-shaded region.Figure 2
The stability region of the trivial singular stateE0, boundary singular state E1, and the coexisting singular state E∗ of (4), in the r0−β parameter-space where r0∈0,0.1 and β∈0,1.0. Other parameter values are δ1=0.0125, γ=0.01, α=0.8, θ=0.4, and δ2=0.05. In the figure, (i) within the blue colored region, the equilibrium point E0 is locally asymptotically stable, but boundary singular point E1 and the interior equilibrium E∗ are unstable; (ii) within the green colored region, the equilibrium point E1 is locally asymptotically stable, but trivial equilibria E0 and interior equilibrium E∗ are unstable; and (iii) within the red colored region, the interior equilibrium point E∗ is locally asymptotically stable, but the trivial equilibrium E0 and boundary equilibrium E1 are unstable.


## 3.4. Comparative Study with Different Kinds of Functional Responses

Mainly, three types of response functions (Holling types-I, II, and III) and Beddington–DeAngelis are widely investigated by the researchers in the predator-prey system [29, 37, 38]. Here, we also investigate the role of different kinds of functional responses on the dynamics for the predator-prey system (2). To exemplify this matter, we performed a relative study for (2) with reference to three kinds of response functions, namely, Beddington–DeAngelis, Holling type-I, and Holling type-II (already used in the model (4)). Holling-type response functions are the function of prey population only; that is, the response functions are independent of predator population. But, in the nature, response functions are also the function of predator population. Widely studied Beddington–DeAngelis [13, 14] response function considered the density of predators also. In our study, we have also studied stability of the system by considering Beddington–DeAngelis functional responses and two types of response function. The Beddington–DeAngelis response function is of the following form:(18)gu,v=αu1+βu+ξv,where the parameters α,β, and ξ are nonnegative and have their usual biological meanings. The Beddington–DeAngelis response function is a similar kind of Holling type-II response function but contains an additional term ξv at the denominator. The term ξv indicates the mutual interaction between the predator population. Thus, when the prey density ut is constant and predator density vt is minimum, then the value of the response function gu,v is maximum; when predator density vt is maximum, then the value of the response function gu,v is minimum. This is more realistic because when the number of predator species are less than their prey, consumption due to predators per unit time will be more, but when predator density is more, then the prey consumption due to predators per unit time will be less due to mutual interaction of predator population. But, when ξ=0, then the Beddington–DeAngelis functional response is the same as Holling type-II functional response. Again, for β=0andξ=0, then Beddington–DeAngelis response function is the same as Holling type-I response function. Hence, we may conclude that Beddington–DeAngelis response function is a more generalized form of Holling type-II and type-I functional responses. By considering Beddington–DeAngelis response function, the system (2) becomes(19)dudt=r0u−δ1u−γu2−αuv1+βu+ξv,dvdt=αθuv1+βu+ξv−δ2v.The predator-prey system (19) with Beddington–DeAngelis functional response has three ecologically relevant singular points, namely, trivial singular point E00,0, boundary singular point E1u^,0, and an interior steady state E∗u^∗,v^∗. Here, u^=r0−δ1/γ, and u^∗ is the positive root for the following second-degree equation:(20)θγξu^∗2−r0−δ1θξ−αθ−βδ2u^∗−δ2=0,v^∗=1ξδ2αθ−βδ2u^∗−δ2.We have made a relative investigation of the predator-prey model (2), by considering three different functional responses, which are type-I (model (19) with β=0,ξ=0), Holling type-II (already used in model (4) with ξ=0), and Beddington–DeAngelis. We are not showing detailed analysis of the model (2) for the different functional responses; rather, we have presented the stability scenarios for these three types of response function in Table 1. The parameters r0, θ, and γ are varied to observed the stability behavior of the system. The detailed numerical study and graphical illustrations have been shown in numerical Section 5.Table 1
Stability conditions for different equilibrium points of system (19).
EquilibriaHolling type-I (β=0,ξ=0)Holling type-II (ξ=0)Beddington–DeAngelisE00,0r0<δ1r0<δ1r0<δ1E1u^,0r0>δ1 and αθr0−δ1<δ2r0>δ1 and αθγr0−δ1/γ+βr0−δ12<δ2r0>δ1 and αθγr0−δ1/γ+βr0−δ12<δ2E∗u^∗,v^∗M1+M4<0 and M1M4−M2M3>0 with β=0andξ=0M1+M4<0 and M1M4−M2M3>0 with ξ=0M1+M4<0 and M1M4−M2M3>0M1,M2,M3, and M4 are the elements of the variational matrix around interior singular point E∗u^∗,v^∗, where M1=−γu^∗−αβu^∗v^∗/1+βu^∗+ξv^∗2, M2=−αu^∗1−ξv^∗/1+βu^∗+ξv^∗2, M3=αθv^∗1−βu^∗/1+βu^∗+ξv^∗2, and M4=αθu^∗1−ξv^∗/1+βu^∗+ξv^∗2−δ2.


## 3.5. Global Stability Analysis aroundE∗u∗,v∗

In this subsection, we investigate the sufficient condition for global asymptotic stability of the positive coexisting singular pointE∗u∗,v∗ by establishing appropriate Lyapunov function.Theorem 6.
The positive interior steady stateE∗u∗,v∗ of the system (4) is globally asymptotically stable if γ/αβ>v∗.Proof.
The interior steady stateE∗u∗,v∗ of the prey-predator model (4) fulfills the following conditions:(21)r0−δ1=γu∗+αv∗1+βu∗,δ2=αθu∗1+βu∗.
With the help of these above two equations, the model (4) leads to(22)dudt=u−γu−u∗+αv∗1+βu∗−αv1+βu,dvdt=vαθu1+βu−αθu∗1+βu∗.
We define the Lyapunov functionalVu,v:ℜ2+⟶ℜ such that(23)Vu,v=V1u,v+V2u,v,where V1u,v=u−u∗−u∗lnu/u∗ and V2u,v=Kv−v∗−v∗lnv/v∗. Here, K is a positive constant and is defined as K=1/θ1+βu∗. This particular type of Lyapunov function to study the stability of the interior equilibrium has been used widely [29, 30, 39, 40]. The Lyapunov function Vu,v is continuous on Intℜ2+ and vanishes at E∗u∗,v∗. Also, Vu,v is positive definite in Intℜ2+. Furthermore, ∂V1/∂u>0 when u>u∗, ∂V1/∂u<0 when u<u∗, ∂V2/∂u>0 when v>v∗, and ∂V2/∂v<0 when v<v∗. Hence, E∗u∗,v∗ is a global minimum of Vu,v. The time derivative of V1u,v and V2u,v along the solution for (4) is(24)dV1dt=u−u∗ududt=u−u∗−γu−u∗+αv∗1+βu∗−αv1+βu,dV2dt=Kv−v∗vdvdt=Kαθv−v∗u1+βu−u∗1+βu∗.
Now, differentiating (23) with respect to time t and using the values of dV1/dt and dV2/dt, we obtain(25)dVdt=dV1dt+dV2dt=−γu−u∗2−αu−u∗v−v∗+βu∗v−uv∗1+βu1+βu∗+Kαθu−u∗v−v∗1+βu1+βu∗.
Using the relationu∗v−uv∗=u∗v−v∗−v∗u−u∗, we obtain(26)dVdt=−γu−u∗2+αβv∗u−u∗21+βu1+βu∗−α1+βu∗−Kθ1+βu1+βu∗u−u∗v−v∗.
Putting the value ofK, we obtain(27)dVdt=−γu−u∗2+αβv∗u−u∗21+βu1+βu∗≤−γ−αβv∗u−u∗2.
Hence,dV/dt<0 along the trajectories in ℜ2+ excluding at E∗u∗,v∗ if γ−αβv∗>0, that is, if γ/αβ>v∗. Thus, if γ/αβ>v∗, then the coexisting singular point E∗u∗,v∗ is globally asymptotically stable.


## 3.6. Hopf Bifurcation Analysis

In our study, the growth rater0 of prey population is one of the most important parameters. In this section, we investigate how predator-prey model behavior alters with respect to r0 by using the analysis of Hopf bifurcation. We start with a theorem to study the analysis of Hopf bifurcation.Theorem 7.
When the prey birth rater0 crosses a threshold value r0=r0∗=δ1+2γu∗1+1/βu∗−1/βu∗δ1+γu∗, the predator-prey model (4) experiences Hopf bifurcation around the interior singular point E∗. Necessary and sufficient conditions of Hopf bifurcation to occur are that there exists r0=r0∗ such that the transversality condition holds:(28)ddr0Reλr0r0=r0∗=βu∗1+βu∗≠0.Proof.
A necessary condition for the change in dynamics of the coexisting singular pointE∗u∗,v∗ of the prey-predator model (4) is that the characteristic polynomial (17) should have purely complex eigenvalues. The characteristic roots of the equation (17) will have purely complex roots if and only if K11=0 and K22>0. For K11=0 gives, r0=r0∗=δ1+2γu∗1+1/βu∗−1/βu∗δ1+γu∗, where u∗=δ2/αθ−βδ2. Due to Theorem 5, we can conclude that the positive coexisting singular point E∗u∗,v∗ of the predator-prey model (4) will be stable asymptotically if r0−δ1−2γu∗−αv∗/1+βu∗2<0 and the model (4) experiences Hopf bifurcation at r0=r0∗.
Let the roots of the characteristic (17) are of the form λ1,2=ϕr0±iψr0. Then, λ1+λ2=2ϕr0=K11; that is, Reλr0=ϕr0=1/2K11. To prove the transversality condition, we need to verify the following condition:(29)ddr0Reλr0r0=r0∗=βu∗1+βu∗≠0,ifβ≠0.
Hence, the transversality condition holds, and the system (4) experiences Hopf bifurcation at r0=r0∗.


## 4. Direction and Stability of Hopf Bifurcation

To study the stability of bifurcating limit cycle appearing through Hopf bifurcation, we use the theorem of center manifold. Center manifold theorem is an essential tool to reduce the dimension of a differential equation system in a neighborhood of a coexisting steady state (see the book by Guckenheimer and Holmes, 1983, [41]). For the predator-prey system (4), we obtain the variational matrix J∗ has a pair of complex conjugate eigenvalues for Hopf bifurcation [42]. Thus, we can investigate the model on a 2-dimensional center manifold [43].To determine center manifold and investigate the flow here, first interpret the origin for the coordinate system to an interior singular pointE∗u∗,v∗ by expressing u^=u−u∗ and v^=v−v∗. For simplicity u^ and v^ denoted by u and v, respectively, the system (4) leads to(30)dudt=r0−δ1u+u∗−γu+u∗2−αu+u∗v+v∗1+βu+u∗,dvdt=αθu+u∗v+v∗1+βu+u∗−δ2v+v∗.The system of (30) can be rewritten as(31)utvt=J∗r0uv+f1u,v,r0f2u,v,r0,where J∗r0 defined in (16), and(32)f1u,v,r0=a1u2+a2uv+a3v2+…,f2u,v,r0=−a1+γθu2−a2θuv+a3v2+…,where a1=−γ+αβv∗/1+βu∗3, a2=−α/1+βu∗2, and a3=0.Let the characteristic roots of the Jacobian matrixJ∗r0 are of the form λ1,2=ϕr0±iψr0, where ϕr0=1/2trJ∗ and ψr0=detJ∗−ϕ2r0. Here, the characteristics roots λ1,2 are complex conjugate if detJ∗−ϕ2r0>0 and the characteristic roots λ1,2 are purely imaginary; that is, λ1,2=±iψr0 when ϕr0=0; that is, when r0=r0∗. Now set the following matrix to get the normal form of the predator-prey system (4):(33)D=10PQ,where 1P−iQ is the eigenvector corresponding to λ=ϕr0+iψr0 with P=1+βu∗/αu∗αβu∗v∗/1+βu∗2−γu∗−ϕ and Q=1+βu∗ψ/αu∗. Then, we have(34)D−1=10−QP1Q.By using the transformation(35)uv=Dxy,System (30) becomes(36)xtyt=L0r0xy+Gx,y,r0Hx,y,r0,where(37)L0r0=ϕ−ψψϕ,with(38)Gx,y,r0=ϕ−βγu∗21+βu∗u∗x2−ϕ1+βu∗u∗xy+0×y2−βγu∗+ϕ1+βu∗2u∗x3+βϕ1+βu∗2u∗x2y+0×xy2+0×y3+…,Hx,y,r0=−PQGx,y,r0+1Qg2x,y,r0,g2x,y,r0=αθP1+βu∗2−βv∗1+βu∗3x2+0×y2+αθQ1+βu∗2xy−αβθP1+βu∗3−βv∗1+βu∗4x3−αβθQ1+βu∗3x2y+0×xy2+0×y3+….Now, rewriting the expression (36) in polar coordinate r,ζ and then expanding in the Taylor series expansion at r0=r0∗, we obtain (39)ζt=ψr0+cr0r2+…=ψr0∗+ψ′r0r0−r0∗+cr0∗r2+…,rt=ϕr0r+ar0r3+…=ϕ′r0∗r0−r0∗r+ar0∗r3+….To investigate the stability of Hopf bifurcating periodic solution, we must have to compute the sign of the coefficient ofar0, where ar0 is given by(40)ar0∗=116Gxxx+Gxyy+Hxxy+Hyyy0,0,r0∗+116ψr0∗GxyGxx+Gyy−HxyHxx+Hyy−GxxHxx+GyyHyy0,0,r0∗.The above expression ofar0∗ is obtained by well-known Hopf bifurcation formula to transform (36) into Jordan form [41]. Here, the subscript denotes partial derivative, and we have(41)Gxxx0,0,r0∗=−3βγu∗+ϕ1+βu∗2,Gxyy0,0,r0∗=Hyyy0,0,r0∗=Gyy0,0,r0∗=Hyy0,0,r0∗=0,Hxxy0,0,r0∗=−2PβϕQu∗1+βu∗2−2αβθ1+βu∗3,Gxy0,0,r0∗=−ϕ1+βu∗u∗,Gxx0,0,r0∗=2ϕ−βγu∗21+βu∗u∗,Hxy0,0,r0∗=PϕQ1+βu∗u∗+αθ1+βu∗2,Hxx0,0,r0∗=−2Pϕ−βγu∗2Qu∗1+βu∗+2αθQP1+βu∗2−βv∗1+βu∗3.Here,ϕr0=1/2trJ∗=1/2αβu∗v∗/1+βu∗2−γu∗ and dϕ/dr0=βu∗/21+βu∗>0.We obtain thatΛ=−ar0/ϕ′r0.From the above computation ofar0∗, the following theorem can be statedTheorem 8.
(i)
The direction of Hopf bifurcation is supercritical, and bifurcating periodic solutions are stable ifar0∗<0.(ii)
The direction of Hopf bifurcation is subcritical, and bifurcating periodic solutions are unstable ifar0∗>0.


## 5. Numerical Simulations

In this section, we compute some numerical simulations regarding stability and bifurcation for the predator-prey model (4) in the vicinity of the coexisting singular point E∗u∗,v∗. It is quite difficult to verify the mathematical model simulations with realistic parameter values. We have taken a hypothetical set of parameter values to illustrate our analytical findings. The set of model parameters is as follows: r0=0.01, δ1=0.02, γ=0.01, α=0.9, β=0.5, θ=0.4, and δ2=0.05. For this parameter, set the stability condition r0<δ1 for the trivial singular point E0 is satisfied, and hence, E00,0 is locally asymptotically stable. From the time series solution shown in Figure 3(a), it can be noticed that both prey and predator populations ut and vt initiate from the initial values 0.5,0.3 and goes to extinct; hence, the trivial steady state E0 is stable asymptotically. Again, we have chosen the parameter values as r0=0.03, δ1=0.01, γ=0.15, α=0.9, β=0.5, θ=0.4, and δ2=0.05. The condition of locally asymptotic stability δ1<r0<δ1+γδ2/αθ−βδ2 of the boundary equilibrium point E1u¯,0 is satisfied, and hence, E1u¯,0 is locally asymptotically stable. Time series solution demonstrated in Figure 3(b) represents that both predator and prey populations ut and vt initiate from an initial value 0.5,0.3 and go to their boundary equilibrium state, and hence, the boundary steady state E10.1333,0 is stable asymptotically. Also, we have chosen the parameter set as r0=0.03, δ1=0.02, γ=0.01, α=0.9, β=0.5, θ=0.4, and δ2=0.05. The condition of locally asymptotic stability r0−δ1−2γu∗−αv∗/1+βu∗2=−0.0009<0 of the coexisting singular point E∗u∗,v∗ is satisfied, and hence, E∗u∗,v∗ is locally asymptotically stable. From the time series analysis shown in Figure 3(c), it can be demonstrated that both prey and predators ut and vt initiate from the initial value 0.5,0.3 and go to their interior singular point; hence, the coexisting steady state E∗0.1493,0.0102 is locally asymptotically stable.Figure 3
The time series evolution of prey populationut and predator population vt for the model system (4) with parameters is defined in numerical section. The subplot (a) shows that trivial singular point E00,0 is locally asymptotically stable, (b) boundary singular point E10.1333,0 is locally asymptotically stable, and (c) the coexisting steady state E∗0.1493,0.0102 is stable asymptotically.
(a)(b)(c)The dynamical behavior around the coexisting steady stateE∗u∗,v∗ of the system (4) has been shown in the phase diagram (Figure 4). The system (4) is integrated by using fourth-order Runge–Kutta scheme with the parameter set δ1=0.01, γ=0.01, α=0.9, β=0.5, δ2=0.05, and θ=0.4 and varies the prey birth rate r0. For r0=0.018, the positive interior singular point is given by E∗0.1493,0.0078. To verify Theorem 5, we find K11=r0−δ1−2γu∗−αv∗/1+βu∗2=−0.1113<0 and K22=0.0003>0, which ensure that the asymptotic stability of coexisting state E∗u∗,v∗. Phase portrait diagram (see Figure 4(a)) demonstrates that the predator-prey model (4) is stable asymptotically at coexisting equilibrium state E∗u∗,v∗ as we have chosen an initial point E0.5,0.3 and both the species move around interior steady state E∗0.1493,0.0078 and ultimately go to E∗. The system (4) enters into the limit cycle behavior from stable one, for the value of prey birth rate r0 gradually increases, which means that, for increasing the size of prey population, the predator has taken more time to handle the prey species. Thus, the magnitude of the prey oscillation increases due to increasing size of the prey species. Figure 4(b) demonstrates the limit cycle oscillations around coexisting state E∗0.1493,0.0400 for r0=0.045, and the other parameter values are the same as in Figure 4(a). In this case, the values of the coefficients of the characteristic (17) are K11=0.000834>0 and K22=0.001559>0, which does not satisfy the conditions for locally asymptotically stability of E∗ in Theorem 5. As the values of K11 and K22 are both positive for r0=0.045, the solutions of the system (4) are unstable, and this ensures the limit cycle oscillation around the coexisting equilibrium E∗0.1493,0.0400.Figure 4
The phase portrait diagram of (4) around the coexisting singular point E∗(green filled circle) with initial values x0,y0=0.5,0.3. Black vertical line represents the predator isocline, and the red curve represents the prey isocline. The subplot (a) shows the phase portrait diagram is stable asymptotically at E∗0.1493,0.0078 for r0=0.018, and (b) the phase portrait diagram shows the limit cycle oscillation around E∗0.1493,0.0400 for r0=0.018; the other parameter values are δ1=0.01, γ=0.01, α=0.9, β=0.5, θ=0.4, and δ2=0.5.
(a)(b)Analytically, we have shown that the coexisting stateE∗u∗,v∗ for the system (4) experiences Hopf bifurcation with respect to the prey birth rate r0. By numerical simulations, Figure 5 exhibits that the predator-prey system (4) experiences Hopf bifurcation around the threshold values r0∗=0.0330. For the set of parameter values specified in this section, we obtain d/dr0Reλr0r0=r0∗=0.0694>0. This indicates that the transversality condition for Hopf bifurcation is verified. Hence, the coexistence steady state E∗u∗,v∗ of the predator-prey system (4) is asymptotically stable for r0<r0∗ and unstable for r0>r0∗. It can also be noticed that the system (4) exhibits oscillatory behavior from the stable one, when the value of prey birth rate r0 is increased gradually, and the high amplitude oscillation may lead to the crash of the species [44, 45]. Hence, the prey growth rate r0 has a critical role for the stability of the system (4).Figure 5
Bifurcation diagram with respect to the prey birth rater0 of the predator-prey system (4) with Holling type-II response function. Other parameter values are δ1=0.01, γ=0.01, α=0.9, β=0.5, δ2=0.05, and θ=0.4.
(a)(b)In our study, we have not shown analytically the bifurcation for the predator-prey system (4) with respect to the conversion coefficient θ and the intraspecific coefficient γ; here, we have numerically plotted the bifurcation figure for (4) with reference to the parameters θ and γ. The bifurcation figure for (4) with respect to θ is shown in Figure 6. Here, r0 is predetermined at 0.07, and the rest of parameters are the same as in Figure 5. The system (4) shows the limit cycle oscillation from equilibrium state as the bifurcation parameter θ increases. The system (4) experiences Hopf bifurcation around θ≈0.055. The bifurcation diagram for the system (4) with reference to the intraspecific competition coefficient γ is shown in Figure 7. Here, r0 is predetermined at 0.07, and the rest of parameters are the same as in Figure 5. This bifurcation diagram shows that the model alters the stability from limit cycle behavior to equilibrium state as bifurcation parameter γ increases.Figure 6
Bifurcation diagram with respect to conversion termθ for the predator-prey model (4) with Holling type-II response function. Rest of parameters: r0=0.07, δ1=0.01, γ=0.01, α=0.9, β=0.5, and δ2=0.05.
(a)(b)Figure 7
Bifurcation diagram with respect to intraspecific competition between prey speciesγ for predator-prey system (4) with Holling type-II response function. Rest of the parameters r0=0.07, δ1=0.01, α=0.9, β=0.5, θ=0.4, δ2=0.05.
(a)(b)Again, the dynamics around the interior singular point of the Beddington–DeAngelis-type predator-prey system (19) is shown in Figure 8. The system (19) is integrated using the Runge–Kutta scheme with parameter set δ1=0.01, γ=0.01, α=0.9, β=0.5, δ2=0.05, θ=0.4, ξ=0.3, and varies the prey birth rate r0. For r0=0.018, the predator-prey system (19) has one interior positive singular point E∗0.1496,0.0078. As shown in Table 1, we compute the value of M1+M4=−0.0056 and M1M4−M2M3=0.0002659; therefore, M1+M4<0 and M1M4−M2M3>0. This ensures the stability for the interior singular point for the predator-prey system (19). The phase portrait diagram for Beddington–DeAngelis-type prey-predator system (19) shown in Figure 8(a) demonstrates that the coexisting steady state is locally asymptotically stable. The Beddington–DeAngelis-type prey-predator system (19) enters into the limit cycle behavior from stable one for increasing value of r0. Also, the phase portrait (Figure 8(b)) shows the limit cycle oscillation at an interior singular state E∗0.1511,0.0405 for r0=0.045.Figure 8
Phase portrait diagram of the predator-prey system (19) with Beddington–DeAngelis-type response function at the coexisting singular point E∗ (green filled circle) with initial values x0,y0=0.5,0.3. Black line is the predator isocline, and the red curve is the prey isocline. (a)r0=0.018; (b)r0=0.045, and rest of the parameters are specified as δ1=0.01, γ=0.01, α=0.9, β=0.5, θ=0.4, δ2=0.05, and ξ=0.30.
(a)(b)The bifurcation diagram of the Beddington–DeAngelis-type predator-prey system (19) with respect to growth rate r0 for prey species, the conversion rate θ of prey biomass to predator biomass, and the intraspecific competition rate γ are plotted numerically in Figures 9–11, respectively. From the bifurcation diagram of the Holling Type-II predator-prey system (4) and Beddington–DeAngelis-type predator-prey system (19), it can be observed that we obtained the similar kind of dynamics of the prey and predator species except the changes in the bifurcation point. The bifurcation in Figure 9 represents that the predator-prey system (19) is asymptotically stable for lower critical values of prey birth rate r0, and if r0 crosses threshold value r0c≈0.035, the Beddington–DeAngelis-type predator-prey system (19) shows limit cycle behavior or oscillating behavior. The bifurcation diagram for the Beddington–DeAngelis-type predator-prey system (19) is shown in Figure 10 with respect to conversion parameter θ. Figure 10 exhibits that the system alters the stability from equilibrium state to limit cycle behavior as the bifurcation parameter θ increases. The bifurcation diagram shown in Figure 11 of the Beddington–DeAngelis-type predator-prey system (19) with respect to the intraspecific competition rate γ shows that the model alters the stability from limit cycle behavior to equilibrium state for increasing bifurcation term γ.Figure 9
Bifurcation diagram with respect to birth rate of prey speciesr0 of the predator-prey system (19) with Beddington–DeAngelis-type response function. Other parameter values are specified as δ1=0.01, γ=0.01, α=0.9, β=0.5, δ2=0.05, θ=0.4, and ξ=0.30.
(a)(b)Figure 10
Bifurcation diagram with respect to conversion termθ for predator-prey system (19) with Beddington–DeAngelis-type response function. Rest of parameters are r0=0.07, δ1=0.01,γ=0.01, α=0.9, β=0.5, δ2=0.05, and ξ=0.30.
(a)(b)Figure 11
Bifurcation diagram with respect to intraspecific competition between prey speciesγ for predator-prey system (19) with Beddington–DeAngelis response function. Rest of parameters are specified as r0=0.07, δ1=0.01, α=0.9, β=0.5, θ=0.4, δ2=0.05, and ξ=0.30.
(a)(b)Our predator-prey system (4) with Holling type-II response function undergoes transcritical bifurcation with respect to the prey growth rate r0, which is plotted in the bifurcation diagram in Figure 12, and the other parameter values are δ1=0.01, γ=0.01, α=0.9, β=0.5, θ=0.4, and δ2=0.05. The bifurcating parameter r0 varied from 0 to 0.04 and is plotted r0 along the x−axis and equilibrium prey (left panel in Figure 12) and predator (right panel in Figure 12) density along the y−axis . Black line represents the stable branch of the trivial equilibrium point E00,0, blue-dotted line represents the stable branch of the boundary equilibrium point E1u¯,0, and the green line represents the unstable branch of the boundary equilibrium point E1u¯,0. The red horizontal line (left panel) and the red oblique line (right panel) indicate the stable branch of the interior equilibrium point E∗u∗,v∗, and the magenta line indicates the unstable branch of the interior steady state E∗u∗,v∗. There is no stable branch of the boundary steady state E1u¯,0 in the plot of predator density (right panel in Figure 12) as there is no predator density in the boundary equilibrium E1u¯,0. Both the prey and predator populations oscillate in the blue-shaded region.Figure 12
The transcritical bifurcation diagram with respect to the prey birth rater0 for the predator-prey system (4) with Holling type-II response function. Other parameter values are specified as δ1=0.01, γ=0.01, α=0.9, β=0.5, θ=0.4, and δ2=0.05.
(a)(b)Our predator-prey system (19) with Beddington–DeAngelis response function also undergoes transcritical bifurcation with respect to the prey growth rate r0, which is plotted in the bifurcation diagram in Figure 13, and the other parameter values are δ1=0.01, γ=0.01, α=0.9, β=0.5, θ=0.4, δ2=0.05, and ξ=0.30. The bifurcating parameter r0 varied from 0 to 0.04 and is plotted r0 along the x−axis and equilibrium prey (left panel in Figure 13) and predator (right panel in Figure 13) density along the y−axis. Black line represents the stable branch of the trivial equilibrium point E00,0, blue-dotted curve represents the stable branch of the boundary equilibrium point E1u^,0, and the green line represents the unstable branch of the boundary equilibrium point E1u^,0. Red curve indicates the stable branch of the interior equilibrium point E∗u^∗,v^∗, and the magenta line indicates the unstable branch of the interior steady state E∗u^∗,v^∗. There is no stable branch of the boundary steady state E1u^,0 in the plot of the predator density (right panel in Figure 13) as there is no predator density in the boundary equilibrium E1u^,0. Both the prey and predator populations oscillate in the blue-shaded region. From the transcritical bifurcation for Holling type-II and Beddington–DeAngelis-type response function, we can conclude that density of prey species increases (see the read curve of Figures 12 and 13) for Beddington–DeAngelis-type response function. The red curve (left panel in Figure 13) is increasing as the value of r0 is increasing; that is, the density of prey species is increasing as the prey growth rate r0 is increasing. To better understand the increasing pattern of prey species for Beddington–DeAngelis-type response function, we put a portion of the red curve in the dash-dotted box, which has been zoomed and plotted in the inset (left panel of Figure 13).Figure 13
Bifurcation diagram with respect to the growth rater0 of the prey species for predator-prey model system (19) with Beddington–DeAngelis response function. Other parameter values are specified as δ1=0.01, γ=0.01, α=0.9, β=0.5, θ=0.4, δ2=0.05, and ξ=0.30. The dash-dotted box is zoomed in the inset, which shows the increasing trend of the stable branch of the interior steady state E∗u^∗,v^∗.
(a)(b)By constructing suitable Lyapunov function, we have analytically shown that (see the Section3.5) the coexisting steady state E∗u∗,v∗ for the predator-prey system (4) is globally asymptotically stable for γ/αβ>v∗. For the following parameter values r0=0.045, δ1=0.01, γ=0.15, α=0.9, β=0.5, θ=0.4, and δ2=0.05, we obtain the interior equilibrium is E∗0.1493,0.0151. Here, γ/αβ=0.3333 and v∗=0.0151, so it can be verified that γ/αβ>v∗ for the assumed set of parameter value. Therefore, the interior steady state E∗0.1493,0.0151 is globally asymptotically stable. Now, we have numerically integrated the predator-prey system (4) using the Runge–Kutta scheme with the assumed parameter values and considered different initial values 0.4,0.1, 0.2,0.2, 0.4,0.1, and 0.1,0.3. The phase diagram of this solution is plotted in Figure 14. All the solution trajectories for different initial values converge to an interior singular point E∗0.1493,0.0151. This numerically ensures that the interior steady state E∗ for the model (4) is globally asymptotically stable.Figure 14
Interior singular pointE∗0.1493,0.0151 for predator-prey system (4) is globally asymptotically stable with the parameters specified as r0=0.045, δ1=0.01, γ=0.15, α=0.9, β=0.5, θ=0.4, and δ2=0.05. The initial values are shown within the bracket.


## 6. Discussion

Mathematical modeling is an important powerful tool to better understand the interplays between various species that allow us to continue biodiversity and ecosystem in essence. In literature, a series of mathematical models has been explored to understand the dynamics for prey-predator model with different kinds of functional responses [12]. In this paper, we investigate a very simple mathematical model for predator-prey system with rich dynamics. Here, we consider a prey-predator system with Holling type-II functional response. The present investigation has been done in three folds. First, we study a prey-predator model with Holling type-II response function, and then, we perform a comparative study of different kinds of functional responses including Beddington–DeAngelis-type response function. Then, we confirm our analytical calculations with numerical illustrations for the hypothetical parameter set.The solutions of the appraised predator-prey model are showed to be bounded, which indicate none of the species would go at its infinite. To study the kinetics of (4), we discussed the existence of biologically meaningful singular points and their stability including global asymptotic stability. The singular points are identified using prey and predator isocline and found three biologically meaningful steady states of the system (4). The boundary steady state E1u¯,0 shows oscillatory dynamics to the present system (4) as the stability for the model around E1 indicates the infeasibility for the trivial equilibrium point E0. The boundary steady state E1u¯,0 is feasible if r0>δ1; that is, growth rate of the prey species is higher than the decay rate for the prey population. We find the stability region for all the three singular points, namely, trivial, boundary, and interior steady states. In the r0−β plane, the trivial equilibrium is stable in the blue-shaded region, boundary equilibrium is stable in the green-shaded region, and the interior equilibrium is stable in the red-shaded region (see Figure 2). The equilibrium points of the system (4) are stable in mutually disjoint regions.A comparative study has been made in this study by considering Holling type-I, Holling type-II, and Beddington–DeAngelis-type response functions of the predator-prey system (4). The stability conditions for different biologically feasible equilibrium points by considering three different types of response functions are presented in Table 1. Our results demonstrate that birth suppression of prey population by predator population lowers the density of equilibrium point of predator population. Moreover, we notice that, above the critical value of prey birth rate r0∗ of the interplay between prey and predator populations, the system produces limit cycle oscillations. Thus, the prey growth rate plays a critical role to destabilize the predator-prey dynamics through Hopf bifurcation (see Figure 5). Analytically, we have computed the global asymptotic stability condition for the interior equilibrium point of the prey-predator system (4) by formulating an appropriate Lyapunov functional. The coexisting singular point is globally asymptotically stable if γ/αβ>v∗; that is, the predator species is lower than the ratio of intraspecific competition between prey population and the product of attack rate and handling time to capture the prey population by predator population. For the hypothetical parameter set, we have shown numerically the condition for global asymptotic stability of the coexisting singular point E∗. Thus, Theorem 6 is verified numerically.The phase portrait diagram (Figure4) for the system (4) indicates that the dynamical behavior of (4) changes as the birth rate r0 of prey species increases. For increasing the birth rate r0 for prey species, the amplitude of limit cycle increases, which mean that the predator population will take more time to capture the prey population. We have explored the bifurcation scenario of the systems (4) and (19) by varying the birth rate r0 of prey species, the conversion coefficient θ of prey biomass to predator biomass, and the intraspecific competition γ between prey species. Analytically, we have shown that the system (4) experiences a Hopf bifurcation around E∗u∗,v∗ with respect to the prey birth rate at r0=r0∗. Bifurcation diagram with respect to another two parameters θ and γ have been plotted by numerical simulations and found that both the systems (4) and (19) undergo from limit cycle oscillations to equilibrium state situation as the parameter γ increases. Again both the systems (4) and (19) undergo equilibrium state to limit cycle oscillation as the conversion coefficient θ increases. Thus, the dynamical behavior of (19) remains the same as system (4) as only the position of the bifurcation points has been changed.We performed direction and stability for Hopf bifurcation of the predator-prey system (4). The main result of the direction of Hopf bifurcation is presented in Theorem 8. The sign of ar0∗ governs the stability of Hopf bifurcating periodic solution. When the bifurcating periodic solutions are stable, the direction for Hopf bifurcation becomes supercritical, and when the bifurcating periodic solutions are unstable, the direction for Hopf bifurcation becomes subcritical. Theoretical and numerical illustrations for our present paper describe that a simple predator-prey system can have rich dynamics which can be observed in the ecological system. Numeric simulations are employed for the complexity of theoretical findings, where achievable. For the aptitude suggestions of our investigation, there is a permanent need for observational work to test these forecasts. It can also be observed that the illustrations represented in our manuscript should be examined from qualitative, preferably, than quantitative view point. Our study is based on hypothetical set of parameter values and not based on any particular case as we have no experimental data in our hand, but this investigation might be useful for ecologists who are performing their observation in similar field on the basis of experimental data.



---

*Source: 4285294-2020-10-28.xml*