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.
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 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 rate , which gives more complicated dynamics of the system. The predator-prey system also experiences transcritical bifurcation with respect to the prey birth rate . 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 Section 3, 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]:where and represent the predator and prey populations, respectively, at any time and 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, denotes the natural death rate for predator species and stands the conversion coefficient of prey’s biomass to predator’s biomass, and 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 can be represented by and , respectively. Both the prey and predator populations are continuous in time . The predator-prey interaction with response function can be expressed by the system of coupled differential equations:
With initial values and , all the model parameters , , , , and are positive constants and can be interpreted as follows: represents the birth rate of prey species in absence of predators, represents the death rate for prey population, denotes the intraspecific competition of prey species, 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 . Also, the dimensions of , , , and are , , , and , respectively, and the parameter is dimensionless. In exclusion of predator species, the prey population obeys the logistic growth with an intrinsic growth rate .
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: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 and , 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 given by (3), system (2) becomes
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 . Hence, after integration on both sides of each equation of system (4), we get
From the above expressions of and , it is clear that both and remain nonnegative for all finite time, that is, . Hence, is positively invariant of system (4). In order to prove the boundedness of the prey-predator model (4), we consider the function
Differentiation with respect to gives
Now, for each , we getwhere . The maximum value of is . Therefore,
Due to the theory of differential inequality, we havefor , and we have . Therefore, is bounded. Thus, all the solution of (4) are confined in
This implies that the solution of system (4) is bounded. Thus, we have the following theorem.
Theorem 1. For any initial value , 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 such that and , where and are independent of the initial conditions.
Theorem 2. System (4) is uniformly persistent if , where is the upper bound of .
Proof. From the proof of Theorem 1, we have and . We choose . Hence, there exists such that .
Now, from the first equation of system (4), we havewhere is the upper limit of in .
Thus, we can conclude that with . Again, from the second equation of (4), we haveHence, we have . Thus, if we choose with , then . Therefore, system (4) will be uniformly persistent if .
3.2. Equilibria and Their Existence
The singular points for the model (4) are nonnegative solutions for and in ; 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 , and with predator isoclines is stated as and the vertical straight line . System (4) possesses three biologically meaningful singular points, as follows:(i)Trivial steady state .(ii)Boundary singular point , where , which is feasible only when ; that is, growth of prey population is higher than the death rate for prey population.(iii)The positive interior equilibrium point , where and .
The shape and position of the curve depends upon multiple parameters. For different values of , the zero growth prey isoclines are convex curve joining the boundary equilibrium and (see Figure 1(a)). As increases, the vertical predator isocline shifted towards increasing . For different values of , zero growth prey isocline becomes a convex curve, which meets at boundary equilibrium (see Figure 1(b)). As increases, the vertical predator isocline shifted towards decreasing and the prey isocline shifted towards decreasing . Also, for different values of , the zero growth prey isoclines are convex curve, but the zero growth predator isoclines do not change position. As increases, the zero growth prey isocline shifted towards increasing. In all these situations, the interior equilibrium point exists only if (see Figure 1(c)).
(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 iswhere
Theorem 3. If , then the trivial singular point 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 are and . Here, ; therefore, the trivial equilibrium point of the model (4) will be locally asymptotically stable if ; 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 is asymptotically stable, then the boundary singular point does not exist. For , the trivial equilibrium point has eigenvalues 0 and . So is a nonhyperbolic equilibrium point [6].
Theorem 4. If , then the boundary equilibrium point 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 are and . The boundary equilibrium point will be locally asymptotically stable if and , and this condition implies . Biologically, implies that the rate of growth of prey population is higher than the natural rate of death of the prey. Also, implies that the prey growth rate cannot be unbounded. It is to be noted that if the boundary equilibrium state is locally asymptotically stable, then the trivial equilibrium state is unstable.
Theorem 5. If , then the positive interior singular point of (4) will be asymptotically stable and otherwise unstable.
Proof. The variational matrix around the coexisting singular point for the system (4) is given byThe characteristic equation of the above matrix can be expressed aswhere and .
Due to well-known Routh-Hurwitz criterion, the characteristic (17) will have negative real roots if and . Hence, and since . It is obvious that because will be positive interior singular point of (4). Therefore, the positive interior steady state of (4) will be asymptotically stable if and otherwise unstable.
Figure 2 shows the stability region of different singular points of (4) in parameter-space. The growth rate of prey species varies from 0 to 0.1, and the half saturation constant varies from 0 to 1; the other parameters are fixed as , , , , and . For a very small value of , that is, when , the trivial equilibrium is locally asymptotically stable. The stability region of is shown by blue-shaded region in the plane. In this blue-shaded region, the boundary equilibrium and the interior equilibrium are unstable. The boundary equilibrium is locally asymptotically stable in the green-shaded region. But the trivial equilibrium and the interior equilibrium are unstable in this green-shaded region. The red-shaded region of Figure 2 shows the locally asymptotically stable behavior of the interior equilibrium ; that is, both the prey and predator populations can exist together in the red-shaded region. For higher value of (such as ) and for higher value of (such as ), the coexisting equilibrium loses stability. The trivial equilibrium and the boundary equilibrium are unstable in the red-shaded region. All the equilibria lose stability in the white-shaded region.
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: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 at the denominator. The term indicates the mutual interaction between the predator population. Thus, when the prey density is constant and predator density is minimum, then the value of the response function is maximum; when predator density is maximum, then the value of the response function 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 , then the Beddington–DeAngelis functional response is the same as Holling type-II functional response. Again, for , 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
The predator-prey system (19) with Beddington–DeAngelis functional response has three ecologically relevant singular points, namely, trivial singular point , boundary singular point , and an interior steady state . Here, , and 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 ), Holling type-II (already used in model (4) with ), 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 , , 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.
3.5. Global Stability Analysis around
In this subsection, we investigate the sufficient condition for global asymptotic stability of the positive coexisting singular point by establishing appropriate Lyapunov function.
Theorem 6. The positive interior steady state of the system (4) is globally asymptotically stable if .
Proof. The interior steady state of the prey-predator model (4) fulfills the following conditions:With the help of these above two equations, the model (4) leads toWe define the Lyapunov functional such thatwhere and . Here, is a positive constant and is defined as . 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 is continuous on and vanishes at . Also, is positive definite in . Furthermore, when , when , when , and when . Hence, is a global minimum of . The time derivative of and along the solution for (4) isNow, differentiating (23) with respect to time and using the values of and , we obtainUsing the relation , we obtainPutting the value of , we obtainHence, along the trajectories in excluding at if , that is, if . Thus, if , then the coexisting singular point is globally asymptotically stable.
3.6. Hopf Bifurcation Analysis
In our study, the growth rate of prey population is one of the most important parameters. In this section, we investigate how predator-prey model behavior alters with respect to 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 rate crosses a threshold value , the predator-prey model (4) experiences Hopf bifurcation around the interior singular point . Necessary and sufficient conditions of Hopf bifurcation to occur are that there exists such that the transversality condition holds:
Proof. A necessary condition for the change in dynamics of the coexisting singular point 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 and . For gives, , where . Due to Theorem 5, we can conclude that the positive coexisting singular point of the predator-prey model (4) will be stable asymptotically if and the model (4) experiences Hopf bifurcation at .
Let the roots of the characteristic (17) are of the form . Then, ; that is, . To prove the transversality condition, we need to verify the following condition:Hence, the transversality condition holds, and the system (4) experiences Hopf bifurcation at .
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 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 point by expressing and . For simplicity and denoted by and , respectively, the system (4) leads to
The system of (30) can be rewritten aswhere defined in (16), andwhere , , and .
Let the characteristic roots of the Jacobian matrix are of the form , where and . Here, the characteristics roots are complex conjugate if and the characteristic roots are purely imaginary; that is, when ; that is, when . Now set the following matrix to get the normal form of the predator-prey system (4):where is the eigenvector corresponding to with and . Then, we have
By using the transformation
System (30) becomeswherewith
Now, rewriting the expression (36) in polar coordinate and then expanding in the Taylor series expansion at , we obtain
To investigate the stability of Hopf bifurcating periodic solution, we must have to compute the sign of the coefficient of , where is given by
The above expression of is obtained by well-known Hopf bifurcation formula to transform (36) into Jordan form [41]. Here, the subscript denotes partial derivative, and we have
Here, and .
We obtain that .
From the above computation of , the following theorem can be stated
Theorem 8. (i)The direction of Hopf bifurcation is supercritical, and bifurcating periodic solutions are stable if .(ii)The direction of Hopf bifurcation is subcritical, and bifurcating periodic solutions are unstable if .
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 . 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: , , , , , , and . For this parameter, set the stability condition for the trivial singular point is satisfied, and hence, is locally asymptotically stable. From the time series solution shown in Figure 3(a), it can be noticed that both prey and predator populations and initiate from the initial values and goes to extinct; hence, the trivial steady state is stable asymptotically. Again, we have chosen the parameter values as , , , , , , and . The condition of locally asymptotic stability of the boundary equilibrium point is satisfied, and hence, is locally asymptotically stable. Time series solution demonstrated in Figure 3(b) represents that both predator and prey populations and initiate from an initial value and go to their boundary equilibrium state, and hence, the boundary steady state is stable asymptotically. Also, we have chosen the parameter set as , , , , , , and . The condition of locally asymptotic stability of the coexisting singular point is satisfied, and hence, is locally asymptotically stable. From the time series analysis shown in Figure 3(c), it can be demonstrated that both prey and predators and initiate from the initial value and go to their interior singular point; hence, the coexisting steady state is locally asymptotically stable.
(a)
(b)
(c)
The dynamical behavior around the coexisting steady state 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 , , , , , and and varies the prey birth rate . For , the positive interior singular point is given by . To verify Theorem 5, we find and , which ensure that the asymptotic stability of coexisting state . Phase portrait diagram (see Figure 4(a)) demonstrates that the predator-prey model (4) is stable asymptotically at coexisting equilibrium state as we have chosen an initial point and both the species move around interior steady state and ultimately go to . The system (4) enters into the limit cycle behavior from stable one, for the value of prey birth rate 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 for , 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 and , which does not satisfy the conditions for locally asymptotically stability of in Theorem 5. As the values of and are both positive for , the solutions of the system (4) are unstable, and this ensures the limit cycle oscillation around the coexisting equilibrium .
(a)
(b)
Analytically, we have shown that the coexisting state for the system (4) experiences Hopf bifurcation with respect to the prey birth rate . By numerical simulations, Figure 5 exhibits that the predator-prey system (4) experiences Hopf bifurcation around the threshold values . For the set of parameter values specified in this section, we obtain . This indicates that the transversality condition for Hopf bifurcation is verified. Hence, the coexistence steady state of the predator-prey system (4) is asymptotically stable for and unstable for . It can also be noticed that the system (4) exhibits oscillatory behavior from the stable one, when the value of prey birth rate is increased gradually, and the high amplitude oscillation may lead to the crash of the species [44, 45]. Hence, the prey growth rate has a critical role for the stability of the system (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, 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 . The bifurcation diagram for the system (4) with reference to the intraspecific competition coefficient is shown in Figure 7. Here, 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.
(a)
(b)
(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 , , , , , , , and varies the prey birth rate . For , the predator-prey system (19) has one interior positive singular point . As shown in Table 1, we compute the value of and ; therefore, and . 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 . Also, the phase portrait (Figure 8(b)) shows the limit cycle oscillation at an interior singular state for .
(a)
(b)
The bifurcation diagram of the Beddington–DeAngelis-type predator-prey system (19) with respect to growth rate 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 , and if crosses threshold value , 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 .
(a)
(b)
(a)
(b)
(a)
(b)
Our predator-prey system (4) with Holling type-II response function undergoes transcritical bifurcation with respect to the prey growth rate , which is plotted in the bifurcation diagram in Figure 12, and the other parameter values are , , , , , and . The bifurcating parameter varied from 0 to 0.04 and is plotted along the and equilibrium prey (left panel in Figure 12) and predator (right panel in Figure 12) density along the . Black line represents the stable branch of the trivial equilibrium point , blue-dotted line represents the stable branch of the boundary equilibrium point , and the green line represents the unstable branch of the boundary equilibrium point . The red horizontal line (left panel) and the red oblique line (right panel) indicate the stable branch of the interior equilibrium point , and the magenta line indicates the unstable branch of the interior steady state . There is no stable branch of the boundary steady state in the plot of predator density (right panel in Figure 12) as there is no predator density in the boundary equilibrium . Both the prey and predator populations oscillate in the blue-shaded region.
(a)
(b)
Our predator-prey system (19) with Beddington–DeAngelis response function also undergoes transcritical bifurcation with respect to the prey growth rate , which is plotted in the bifurcation diagram in Figure 13, and the other parameter values are , , , , , , and . The bifurcating parameter varied from 0 to 0.04 and is plotted along the and equilibrium prey (left panel in Figure 13) and predator (right panel in Figure 13) density along the . Black line represents the stable branch of the trivial equilibrium point , blue-dotted curve represents the stable branch of the boundary equilibrium point , and the green line represents the unstable branch of the boundary equilibrium point . Red curve indicates the stable branch of the interior equilibrium point , and the magenta line indicates the unstable branch of the interior steady state . There is no stable branch of the boundary steady state in the plot of the predator density (right panel in Figure 13) as there is no predator density in the boundary equilibrium . 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 is increasing; that is, the density of prey species is increasing as the prey growth rate 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).
(a)
(b)
By constructing suitable Lyapunov function, we have analytically shown that (see the Section 3.5) the coexisting steady state for the predator-prey system (4) is globally asymptotically stable for . For the following parameter values , , , , , , and , we obtain the interior equilibrium is . Here, and , so it can be verified that for the assumed set of parameter value. Therefore, the interior steady state 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 , , , and . 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 . This numerically ensures that the interior steady state for the model (4) is globally asymptotically stable.
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 shows oscillatory dynamics to the present system (4) as the stability for the model around indicates the infeasibility for the trivial equilibrium point . The boundary steady state is feasible if ; 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 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 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 ; 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 . Thus, Theorem 6 is verified numerically.
The phase portrait diagram (Figure 4) for the system (4) indicates that the dynamical behavior of (4) changes as the birth rate of prey species increases. For increasing the birth rate 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 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 with respect to the prey birth rate at . 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 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.
Data Availability
All relevant data are included within the paper.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The 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). The work of Subhas Khajanchi was supported by the Science and Engineering Research Board (SERB) (file no. ECR/2017/000234), Department of Science and Technology, Government of India.