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

In this study, we investigate a mathematical model that describes the interactive dynamics of a predator-prey system with diﬀerent 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 diﬀerent 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 ﬁxed 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 ﬁndings are also discussed in this study.


Introduction
e interaction between prey species and predator species constitutes biological processes of major ecological importance that are ubiquitous in nature. e 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][3][4][5][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. e 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. e dynamical behavior for the prey-predator system takes part in a significant role in theoretical ecology. ere 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. erefore, 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][6][7][8][9][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. e response functions are the rate of prey consumption by the predator per unit time. e functional responses proposed by Holling [12] are the functions of prey density only. e 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. e 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][17][18][19][20][21][22][23][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. ey 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. ey also defined the basic reproduction number R 0 and show that it is a threshold parameter determining whether the predators can coexist with the prey species. Kohnke et al. [18] studied a predatorprey 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. eir 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. e 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. eir model demonstrates nonnegative equilibria as the global attractor. Kon et al. [27] established the criteria for the permanence of stage structure prey-predator interactions.
ese are some existing articles which focused on Beddington-DeAngelis-type functional responses [28][29][30][31]. Beddington-DeAngelis-type response function confesses plentiful and biologically relevant prey-predator dynamics [29], which attract researchers to further study the predatorprey 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. Predatorprey 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][33][34][35].
e 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 rate r 0 , which gives more complicated dynamics of the system. e predator-prey system also experiences transcritical bifurcation with respect to the prey birth rate r 0 . e biological implications of the theoretical and numerical results are also investigated.
is article is organized in the following way: we present a predator-prey model with different kinds of functional responses in the next section. In Section 3, we discuss the qualitative behavior for the model, and positivity for the model, boundedness for the solutions, and the position of all 2 Complexity possible feasible equilibria and the existence criteria of the equilibria are presented. e 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.

The Model
e interaction between prey and predator species in simplest form is represented by the following system of coupled differential equations [2,36]: where v(t) and u(t) represent the predator and prey populations, respectively, at any time t and f(u) 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 g(u, 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 time t can be represented by u(t) and v(t), respectively. Both the prey and predator populations are continuous in time t. e predator-prey interaction with response function g(u, v) can be expressed by the system of coupled differential equations: With initial values u(0) � u 0 ≥ 0 and v(0) � v 0 ≥ 0, all the model parameters r 0 , δ 1 , c, δ 2 , and θ are positive constants and can be interpreted as follows: r 0 represents the birth rate of prey species in absence of predators, δ 1 represents the death rate for prey population, c 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. e dimension of prey and predator biomass is Mass. Also, the dimensions of r 0 , δ 1 , c, 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 (r 0 − δ 1 ) > 0.

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: where α represents the attack rate and β is the handling time to capture the prey population; both α and β are positive constants. e dimensions of α and β are Mass − 1 Time − 1 and Mass − 1 , respectively. e predation rate increases as prey density increases, but after a certain stage, rate of predation becomes constant although prey density increases. e functional response depends on prey density but not on predator density. By using the functional form of g(u, v) given by (3), system (2) becomes

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 From the above expressions of u(t) and v(t), it is clear that both u(t) and v(t) remain nonnegative for all finite time, that is, is positively invariant of system (4). In order to prove the boundedness of the prey-predator model (4), we consider the function Complexity Differentiation with respect to t gives Now, for each ϕ(t) > 0, we get where ρ � min δ 1 , (δ 2 /θ) . e maximum value of r 0 u − cu 2 is (r 2 0 /4c). erefore, Due to the theory of differential inequality, we have for t ⟶ ∞, and we have 0 ≤ ϕ(t) ≤ (m/ρ). erefore, ϕ(t) is bounded. us, all the solution of (4) are confined in is implies that the solution of system (4) is bounded. us, we have the following theorem. Theorem 1. For any initial value (u 0 , v 0 ) ∈ R + 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. e 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.  Proof. From the proof of eorem 1, we have lim sup t⟶∞ u(t) ≤ (m/ρ) and limsup t⟶∞ v(t) ≤ (mθ/ρ). We choose τ � max (m/ρ), (mθ/ρ) . Hence, there exists τ > 0 such that max limsup t⟶∞ u(t), limsup t⟶∞ v(t) ≤ τ. Now, from the first equation of system (4), we have where v is the upper limit of v in Ω. us, we can conclude that lim inf t⟶∞ u(t) ≥ ((r 0 − δ 1 − αv)/c) with r 0 > δ 1 + αv. Again, from the second equation of (4), we have erefore, system (4) will be uniformly persistent if r 0 > δ 1 + αv.

Equilibria and eir Existence.
e singular points for the model (4) are nonnegative solutions for F 1 (u, v) � 0 and F 2 (u, v) � 0 in R 2 + ; that is, the equilibrium points of system (4) are the positive point of intersection for zero growth of prey and predator isoclines. e prey isoclines are given by u � 0, and v � ((r 0 − δ 1 − cu)(1 + β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 state E 0 (0, 0). (ii) Boundary singular point E 1 (u, 0), where u � (r 0 − δ 1 )/c, which is feasible only when r 0 > δ 1 ; that is, growth of prey population is higher than the death rate for prey population. (iii) e positive interior equilibrium point e shape and position of the curve v � (r 0 − δ 1 − cu)(1 + βu)/α depends upon multiple parameters. For different values of β, the zero growth prey isoclines are convex curve joining the boundary equilibrium E 1 (u, 0) and (0, ((r 0 − δ 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 E 1 (u, 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 r 0 , the zero growth prey isoclines are convex curve, but the zero growth predator isoclines do not change position. As r 0 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)).

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. e variational matrix of the model (4) where Theorem 3. If r 0 < δ 1 , then the trivial singular point E 0 (0, 0) of (4) is locally asymptotically stable and otherwise unstable.
Proof. e eigenvalues of the Jacobian matrix (14) of the system (4) at the trivial equilibrium point E 0 (0, 0) are r 0 − δ 1 and − δ 2 . Here, δ 2 > 0; therefore, the trivial equilibrium point E 0 (0, 0) of the model (4) will be locally asymptotically stable if r 0 < δ 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 E 0 (0, 0) is asymptotically stable, then the boundary singular point E 1 (u, 0) does not exist. For r 0 � δ 1 , the trivial equilibrium point E 0 (0, 0) has eigenvalues 0 and − δ 2 . So E 0 (0, 0) is a nonhyperbolic equilibrium point [6]. (4) is locally asymptotically stable and otherwise unstable. Proof.
Proof. e variational matrix around the coexisting singular point E * (u * , v * ) for the system (4) is given by e characteristic equation of the above matrix J * can be expressed as where Due to well-known Routh-Hurwitz criterion, the characteristic (17) will have negative real roots if K 11 < 0 and It is obvious that v * > 0 because E * (u * , v * ) will be positive interior singular point of (4). erefore, the positive interior steady state E * (u * , v * ) of (4) will be asymptotically stable if and otherwise unstable. Figure 2 shows the stability region of different singular points of (4) in r 0 − β parameter-space. e growth rate of prey species r 0 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, c � 0.01, α � 0.8, θ � 0.4, and δ 2 � 0.05. For a very small value of r 0 , that is, when r 0 < δ 1 , the trivial equilibrium E 0 is locally asymptotically stable. e stability region of E 0 is shown by blue-shaded region in the r 0 − β plane. In this blue-shaded region, the boundary equilibrium E 1 and the interior equilibrium E * are unstable. e boundary equilibrium E 1 is locally asymptotically stable in the green-shaded region. But the trivial equilibrium E 0 and the interior equilibrium E * are unstable in this green-shaded region. e 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 r 0 (such as r 0 > 0.06) and for higher value of β (such as β > 0.5), the coexisting equilibrium E * loses stability. e trivial equilibrium E 0 and the boundary equilibrium E 1 are unstable in the red-shaded region. All the equilibria lose stability in the white-shaded region.

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. e Beddington-DeAngelis response function is of the following form: where the parameters α, β, and ξ are nonnegative and have their usual biological meanings. e Beddington-DeAngelis response function is a similar kind of Holling type-II response function but contains an additional term ξv at the denominator. e term ξv indicates the mutual interaction between the predator population. us, when the prey density u(t) is constant and predator density v(t) is minimum, then the value of the response function g(u, v) is maximum; when predator density v(t) is maximum, then the value of the response function g(u, v) is minimum. is 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 β � 0 and ξ � 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 6 Complexity responses. By considering Beddington-DeAngelis response function, the system (2) becomes e predator-prey system (19) with Beddington-DeAngelis functional response has three ecologically relevant singular points, namely, trivial singular point E 0 (0, 0), boundary singular point E 1 (u, 0), and an interior steady state E * (u * , v * ). Here, u � (r 0 − δ 1 /c), and u * is the positive root for the following second-degree equation: 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. e parameters r 0 , θ, and c are varied to observed the stability behavior of the system. e detailed numerical study and graphical illustrations have been shown in numerical Section 5.

Global Stability
Analysis around E * (u * , v * ). In this subsection, we investigate the sufficient condition for global asymptotic stability of the positive coexisting singular point E * (u * , v * ) by establishing appropriate Lyapunov function.
e interior steady state E * (u * , v * ) of the preypredator model (4) fulfills the following conditions: With the help of these above two equations, the model (4) leads to We define the Lyapunov functional V(u, v): where Here, K is a positive constant and is defined as K � (1/θ)(1 + βu * ). is particular type of Lyapunov function to study the stability of the interior equilibrium has been used widely [29,30,39,40]. (ii) within the green colored region, the equilibrium point E 1 is locally asymptotically stable, but trivial equilibria E 0 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 E 0 and boundary equilibrium E 1 are unstable.
Now, differentiating (23) with respect to time t and using the values of (dV 1 /dt) and (dV 2 /dt), we obtain Using the relation Putting the value of K, we obtain Hence, (dV/dt) < 0 along the trajectories in R + 2 ex- us, if (c/αβ) > v * , then the coexisting singular point E * (u * , v * ) is globally asymptotically stable. □ 3.6. Hopf Bifurcation Analysis. In our study, the growth rate r 0 of prey population is one of the most important parameters. In this section, we investigate how predatorprey model behavior alters with respect to r 0 by using the analysis of Hopf bifurcation. We start with a theorem to study the analysis of Hopf bifurcation.

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 predatorprey system (4), we obtain the variational matrix J * has a pair of complex conjugate eigenvalues for Hopf bifurcation Table 1: Stability conditions for different equilibrium points of system (19). Equilibria Holling type-I (β � 0, ξ � 0) Holling type-II (ξ � 0) Beddington-DeAngelis E 0 (0, 0) 8 Complexity [42]. us, we can investigate the model on a 2-dimensional center manifold [43].

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. e set of model parameters is as follows: r 0 � 0.01, δ 1 � 0.02, c � 0.01, α � 0.9, β � 0.5, θ � 0.4, and δ 2 � 0.05. For this parameter, set the stability condition r 0 < δ 1 for the trivial singular point E 0 is satisfied, and hence, E 0 (0, 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 u(t) and v(t) initiate from the initial values (0.5, 0.3) and goes to extinct; hence, the trivial steady state E 0 is stable asymptotically. Again, we have chosen the parameter values as r 0 � 0.03, δ 1 � 0.01, c � 0.15, α � 0.9, β � 0.5, θ � 0.4, and δ 2 � 0.05. e condition of locally asymptotic stability δ 1 < r 0 < δ 1 + (cδ 2 /(αθ − βδ 2 )) of the boundary equilibrium point E 1 (u, 0) is satisfied, and hence, E 1 (u, 0) is locally asymptotically stable. Time series solution demonstrated in Figure 3(b) represents that both predator and prey populations u(t) and v(t) initiate from an initial value (0.5, 0.3) and go to their boundary equilibrium state, and hence, the boundary steady state E 1 (0.1333, 0) is stable asymptotically. Also, we have chosen the parameter set as r 0 � 0.03, δ 1 � 0.02, c � 0.01, α � 0.9, β � 0.5, θ � 0.4, and δ 2 � 0.05. e condition of locally asymptotic stability r 0 − δ 1 − 2cu * − (α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 u(t) and v(t) 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. e dynamical behavior around the coexisting steady state E * (u * , v * ) of the system (4) has been shown in the phase diagram (Figure 4). e system (4) is integrated by using fourth-order Runge-Kutta scheme with the parameter set δ 1 � 0.01, c � 0.01, α � 0.9, β � 0.5, δ 2 � 0.05, and θ � 0.4 and varies the prey birth rate r 0 . For r 0 � 0.018, the positive interior singular point is given by E * (0.1493, 0.0078). To verify eorem 5, we find K 11 � r 0 − δ 1 − 2cu * − (αv * /(1 + βu * ) 2 ) � − 0.1113 < 0 and K 22 � 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 E(0.5, 0.3) and both the species move around interior steady state E * (0.1493, 0.0078) and ultimately go to E * . e system (4) enters into the limit cycle behavior from stable one, for the value of prey birth rate r 0 gradually increases, which means that, for increasing the size of prey population, the predator has taken more time to handle the prey species. us, 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 r 0 � 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 K 11 � 0.000834 > 0 and K 22 � 0.001559 > 0, which does not satisfy the conditions for locally asymptotically stability of E * in eorem 5. As the values of K 11 and K 22 are both positive for r 0 � 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).
Analytically, we have shown that the coexisting state E * (u * , v * ) for the system (4) experiences Hopf bifurcation with respect to the prey birth rate r 0 . By numerical simulations, Figure 5 exhibits that the predator-prey system (4) experiences Hopf bifurcation around the threshold values r * 0 � 0.0330. For the set of parameter values specified in this section, we obtain (d/dr 0 )(Re(λ(r 0 )))| r 0 �r * 0 � 0.0694 > 0. is 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 r 0 < r * 0 and unstable for r 0 > r * 0 . It can also be noticed that the system (4) exhibits oscillatory behavior from the stable one, when the value of prey birth rate r 0 is increased gradually, and the high amplitude oscillation may lead to the crash of the species [44,45]. Hence, the prey growth rate (r 0 ) has a critical role for the stability of the system (4).
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 c; here, we have numerically plotted the bifurcation figure for (4) with reference to the parameters θ and c. e bifurcation figure for (4) with respect to θ is shown in Figure 6. Here, r 0 is predetermined at 0.07, and the rest of parameters are the same as in Figure 5. e system (4) shows the limit cycle oscillation from equilibrium state as the bifurcation parameter θ increases. e system (4) experiences Hopf bifurcation around θ ≈ 0.055. e bifurcation diagram for the system (4) with reference to the intraspecific competition coefficient c is shown in Figure 7. Here, r 0 is predetermined at 0.07, and the rest of parameters are the same as in Figure 5.
is bifurcation diagram shows that the model alters the stability from limit cycle behavior to equilibrium state as bifurcation parameter c increases.
Again, the dynamics around the interior singular point of the Beddington-DeAngelis-type predator-prey system (19) is shown in Figure 8. e system (19) is integrated using the Runge-Kutta scheme with parameter set δ 1 � 0.01, c � 0.01, α � 0.9, β � 0.5, δ 2 � 0.05, θ � 0.4, ξ � 0.3, and varies the prey birth rate r 0 . For r 0 � 0.018, the predatorprey system (19) has one interior positive singular point E * (0.1496, 0.0078). As shown in Table 1 is ensures the stability for the interior singular point for the predator-prey system (19). e 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.
e Beddington-DeAngelis-type prey-predator system (19) enters into the limit cycle behavior from stable one for increasing value of r 0 . Also, the phase portrait (Figure 8(b)) shows the limit cycle oscillation at an interior singular state E * (0.1511, 0.0405) for r 0 � 0.045. e bifurcation diagram of the Beddington-DeAngelistype predator-prey system (19) with respect to growth rate r 0 for prey species, the conversion rate θ of prey biomass to predator biomass, and the intraspecific competition rate c 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. e bifurcation in Figure 9 represents that the predator-prey system (19) is asymptotically stable for lower critical values of prey birth rate r 0 , and if r 0 crosses threshold value r c 0 ≈ 0.035, the Beddington-DeAngelis-type predator-prey system (19) shows limit cycle behavior or oscillating behavior. e 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. e bifurcation diagram shown in Figure 11 of the Beddington-DeAngelis-type predator-prey system (19) with respect to the intraspecific competition rate c shows that the model alters the stability from limit cycle behavior to equilibrium state for increasing bifurcation term c.
Our predator-prey system (4) with Holling type-II response function undergoes transcritical bifurcation with respect to the prey growth rate r 0 , which is plotted in the bifurcation diagram in Figure 12, and the other parameter values are δ 1 � 0.01, c � 0.01, α � 0.9, β � 0.5, θ � 0.4, and δ 2 � 0.05. e bifurcating parameter r 0 varied from 0 to 0.04 and is plotted r 0 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 E 0 (0, 0), blue-dotted line represents the stable branch of the boundary equilibrium point E 1 (u, 0), and the green line represents the unstable branch of the boundary equilibrium point E 1 (u, 0). e 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 * ).
ere is no stable branch of the boundary steady state E 1 (u, 0) in the plot of predator density (right panel in Figure 12) as there is no predator density in the boundary equilibrium E 1 (u, 0). Both the prey and predator populations oscillate in the blue-shaded region.
Our predator-prey system (19) with Beddington-DeAngelis response function also undergoes transcritical bifurcation with respect to the prey growth rate r 0 , which is plotted in the bifurcation diagram in Figure 13, and  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 E 0 (0, 0), blue-dotted curve represents the stable branch of the boundary equilibrium point E 1 (u, 0), and the green line represents the unstable branch of the boundary equilibrium point E 1 (u, 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 * ). ere is no stable branch of the boundary steady state E 1 (u, 0) in the plot of the predator density (right panel in Figure 13) as there is no predator density in the boundary equilibrium E 1 (u, 0). Both the prey and predator populations oscillate in the blue-shaded region. From the transcritical bifurcation for Holling type-II and Beddington-DeAngelistype 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. e red curve (left panel in Figure 13) is increasing as the value of r 0 is increasing; that is, the density of prey species is increasing as the prey growth rate r 0 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). 14 Complexity By constructing suitable Lyapunov function, we have analytically shown that (see the Section 3.5) the coexisting steady state E * (u * , v * ) for the predator-prey system (4) is globally asymptotically stable for (c/αβ) > v * . For the following parameter values r 0 � 0.045, δ 1 � 0.01, c � 0.15, α � 0.9, β � 0.5, θ � 0.4, and δ 2 � 0.05, we obtain the interior equilibrium is E * (0.1493, 0.0151).
Here, (c/αβ) � 0.3333 and v * � 0.0151, so it can be verified that (c/αβ) > v * for the assumed set of parameter value. erefore, 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 is numerically ensures that the interior steady state E * for the model (4) is globally asymptotically stable.

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. e 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. en, we confirm our analytical calculations with numerical illustrations for the hypothetical parameter set. e 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. e singular points are identified using prey and predator isocline and found three biologically meaningful steady states of the system (4). e boundary steady state E 1 (u, 0) shows oscillatory dynamics to the present system (4) as the stability for the model around E 1 indicates the infeasibility for the trivial equilibrium point E 0 . e boundary steady state E 1 (u, 0) is feasible if r 0 > δ 1 ; that is, growth rate of the prey species is higher than the decay rate Periodic solution Periodic solution (b) Figure 12: e transcritical bifurcation diagram with respect to the prey birth rate r 0 for the predator-prey system (4)  for the prey population. We find the stability region for all the three singular points, namely, trivial, boundary, and interior steady states. In the r 0 − β 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). e 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). e 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 r * 0 of the interplay between prey and predator populations, the system produces limit cycle oscillations. us, 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. e coexisting singular point is globally asymptotically stable if (c/αβ) > 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 * . us, eorem 6 is verified numerically. e phase portrait diagram (Figure 4) for the system (4) indicates that the dynamical behavior of (4) changes as the birth rate r 0 of prey species increases. For increasing the birth rate r 0 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 r 0 of prey species, the conversion coefficient θ of prey biomass to predator biomass, and the intraspecific competition c 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 r 0 � r 0 * . Bifurcation diagram with respect to another two parameters θ and c 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 c increases. Again both the systems (4) and (19) undergo equilibrium state to limit cycle oscillation as the conversion coefficient (θ) increases. us, 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). e main result of the direction of Hopf bifurcation is presented in eorem 8. e sign of a(r * 0 ) 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.
eoretical 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.

Data Availability
All relevant data are included within the paper.

Conflicts of Interest
e authors declare that they have no conflicts of interest. Acknowledgments e work of J. J. Nieto has been partially supported by the Agencia Estatal de Investigación (AEI) of Spain, cofinanced by the European Fund for Regional Development (FEDER) corresponding to the 2014-2020 multiyear financial framework (project MTM2016-75140-P) and Xunta de Galicia (grant ED431C 2019-02). e work of Subhas Khajanchi was supported by the Science and Engineering