Bifurcation Analysis of a Delayed Predator-Prey Model with Holling Type III Functional Response and Predator Harvesting
This paper tries to highlight a delayed prey-predator model with Holling type III functional response and harvesting to predator species. In this context, we have discussed local stability of the equilibria, and the occurrence of Hopf bifurcation of the system is examined by considering the harvesting effort as bifurcation parameter along with the influences of harvesting effort of the system when time delay is zero. Direction of Hopf bifurcation and the stability of bifurcating periodic solutions are also studied by applying the normal form theory and the center manifold theorem. Lastly some numerical simulations are carried out to draw for the validity of the theoretical results.
1. Introduction and Model Description
Differential equation models for interactions between species are one of the classical applications of mathematics to biology, dating back to the first half of this century. The development and use of analytical techniques and the growth of computer power have progressively improved our understanding of these types of models.
The study of population dynamics with harvesting is a subject of mathematical bioeconomics, which in turn related to the optimal management of renewable resources, Clark . Generally the concept of optimal resource management is based on the standard cost benefit criterion which maximizes present values of net economic revenues. This criterion is relevant to both private and public management decisions, although the specification of costs and benefits are not necessarily the same in both cases.
Regulation of exploitation of biological resources has become a problem of major concern nowadays in view of the dwindling resource stocks and the deteriorating environment. Exploitation reduces the biomass of the concerned species, exhibits oscillation, and even causes extinction of some other species. Rosenweig-MacAurtho model experiences oscillation under selective effort of Kar and Ghosh . Legovic et al.  have concluded that harvesting the prey species at maximum level causes the extinction of the predator species in traditional prey-predator system. Kar and Ghosh  and Ghosh and Kar  show that harvesting the prey species at maximum sustainable yield (MSY) level never causes the extinction of the predator species in both ratio-dependent and Holling-Tanner prey-predator systems. It is shown that harvesting the prey species at MSY level may or may not drive the predator population to extinction if intraspecific competition is present among the predator species, Kar and Ghosh . More recently, Ghosh and Kar  somehow managed an example to show that prey harvesting at MSY level may be sustainable in tritrophic food chain model. In conclusion, it is reasonable that harvesting the prey species to achieve MSY will cause the extinction of the predator if predation process follows the Holling type functional response. As MSY policy has a severe impact on other tropic levels in ecosystem, suitable management tools are required for smooth running of sustainable fishing. Creating marine protected areas, maximizing the economic rent, and taxation are some good examples to protect fish population from extinction, over exploitation, as well as to achieve desired stock levels for future generation. Kar and Matsuda  and Kar and Ghosh  have shown that marine protection can have positive impact from biological stock assessment but virtually reduces the economic rent of the fishermen. Recently, Kar and Das  have analyzed a prey-predator model with harvesting where they introduce a tax to regulate the fishery.
It has been noticed over the years by many authors that time delays have come to play an important role in almost all branches of science, for example, ecology and biology. The time delay is considered in the population dynamics when the rate of change of the population is not only a function of the present population but also depends on the past population. Delay is frequently used in a predator-prey model to represent the biological process more accurately.
Wangersky and Cunningham  have considered a predator-prey system with time delay and studied the effect of time delay on the system. Martin and Ruan  have studied the combined effects of prey harvesting and delay on the dynamics of predator-prey systems and focused on three very well-studied delay predator-prey models. For first and third models they were shown that the time delay could cause not only instability and oscillations but also the switching of stabilities, while the prey harvesting changes only the equilibrium values but not the properties of solutions. But in the second model, the time delay induces instability and bifurcation but there is no switching of stabilities; however, increasing the prey harvesting level will help the system to regain its stability. This indicates that the prey harvesting has a stabilizing effect on the dynamics of the model. Chen and Changming  have considered a delayed model with stage structured for prey and investigated nonnegative equilibria, existence of Hopf bifurcation, the stability and direction of Hopf bifurcation by applying the normal from theory, and the centre manifold argument. Zhang et al.  have considered a delayed predator-prey model with constant rate prey harvesting and Holling type III functional response. They studied local stability of equilibria and existence of Hopf bifurcation at interior equilibrium. Direction of Hopf bifurcation and the stability of bifurcating periodic solutions are also studied by applying the normal form theory and the center manifold theorem. Also Kar , Kar and Matsuda , Kar and Pahari [16, 17], Kar and Ghorai , Kuang , and some other authors have discussed delayed predator-prey system.
The purpose of the work is to illustrate the combined effects of harvesting and delay on the dynamics of predator-prey system.
In this paper, we consider the following prey-predator system: with initial conditions , , , , , , and . Here and denote the prey and predator populations, respectively, at time . and , respectively, represent the intrinsic growth rate and carrying capacity of the prey. is the maximum capture rate of the prey, is the death rate of the predator, and is the conversion factor. The term denotes the functional response of the predator, which is termed as Holling type III response function (see ). We assume that the predator population is harvested according to catch-per-unit effort (CPUE) hypothesis  which describes that catch-per-unit effort is proportional to the stock level. Therefore, the harvest at time is , where is the catchability coefficient and is the harvesting effort. The time delay is a constant based on the assumption that the change rate of the predators depends on the number of prey and of predators present at some previous time. , , , , , , , , and are positive constants.
For , system (1) becomes
Our paper is organized in the following way: existence of equilibria of the system (2), the local stability criterion, and uniqueness of limit cycles at interior equilibrium are discussed in Section 2. Influence of harvesting effort is discussed in Section 3. In presence of delay we studied stability and Hopf bifurcation in Section 4. Direction of Hopf bifurcation and the stability of bifurcating periodic solutions are also studied by applying the normal form theory and the center manifold theorem in Section 5. Numerical simulations are given in Section 6. A brief concluding remark is given in Section 7.
2. The Steady States and Their Stability
We now study the existence and nature of the steady states. Particularly we are interested in the interior equilibrium of the system. To begin with, we list all possible steady states of the system (2) as follows: , , and the interior equilibrium , where If holds, then the interior equilibrium exists.
For the stability analysis of the equilibria, let the variational matrix of the system (2) at an arbitrary point be .
At , the eigenvalues of the variational matrix are and . Therefore, is unstable.
At , the eigenvalues of the variational matrix are and . Therefore, is asymptotically stable.
The characteristic equation of the variational matrix is given by where Routh-Hurwitz criterion states that all roots of the characteristic equation (4) have negative real parts if .
Therefore, interior equilibrium will be asymptotically stable if and unstable if .
Theorem 1. If exists with , then is locally asymptotically stable.
2.1. Uniqueness of Limit Cycles
It is known that, for prey-predator systems, existence and stability of a limit cycle are related to the existence and stability of a positive equilibrium. If the limit cycles do not exist, in this case the equilibrium is globally asymptotically stable. On the other hand if the positive equilibrium exists and is unstable, there must occur at least one limit cycle.
Let us consider system (2) in the form where , , and .
Now we consider the following Theorem 4.2 of Kuang and Freedman  regarding uniqueness of limit cycles of the above system.
Theorem 2. Suppose for system (6) in and . Then the above system has exactly one limit cycles which is globally asymptotically stable with respect to the set
Following Theorem 2, we may state that when , then system (6) has unique globally stable limit cycle. Thus we see that when the system is unstable, there exists unique globally stable limit cycle.
2.2. Bifurcation Analysis for the Parameter
The characteristic equation (4) has two purely imaginary roots if and only if and for some values of , say . Therefore, there is only one value of at which we have a Hopf bifurcation. Since and contain , it is not easy to find explicitly. So we are trying to find numerically by solving the system of equations with suitable values of the parameters: Therefore the system undergoes a Hopf bifurcation at . The system is stable for and unstable .
3. Influence of the Harvesting Effort
Let the interior equilibrium of the system (2) be , where Now, Therefore is strictly increasing function of .
Also Since , therefore Thus as the harvesting effort increases, the predator population increases when is smaller than the threshold value. But if the harvesting effort gradually increases above the threshold value, that is, as the harvesting effort becomes large enough, then the increasing amount of the harvesting effort can decreases predator population.
To construct Figure 1, we take , , , , , , and in appropriate units. From the above analysis we see that, as the harvesting effort increases, the prey and predator species increase. But if the harvesting effort increases larger than a threshold value, that is, as the harvesting effort becomes large enough, then increasing amount of the harvesting effort can increase the prey species but decrease predator species and go to extinction of the predator species when is large.
4. Stability and Hopf Bifurcation in the Presence of Delay
We will investigate the dynamics of delay system (1).
Let be the only interior equilibrium of the system (1), and let , be the perturbed variables. After removing the nonlinear terms, we obtain the linear variational system, by using equilibrium conditions as where The characteristic equation of the linearized system (14) is the transcendental equation of the following form: where , , , and .
It is known that the steady state is asymptotically stable if all roots of the characteristic equation (16) have negative real parts. For a nonlinear delay equation, there are two types of stabilities: absolute stability (independent of the delay) and conditional stability (depending of the delay).
Case 1 (). Equation (16) becomes Assume that all roots of (16) have negative real parts, which is true if ,. Therefore, the interior equilibrium is locally asymptotically stable if and hold.
Case 2 (). We want to determine if the real part of some root increases to reach zero and eventually becomes positive as varies. If is a root of (16), then by separating the real and imaginary parts, we get
From (19), we obtain the fourth-order equation for as
From (19) we also have and .
From (20), it follows that, if and ,then (20) does not have positive roots. Therefore, characteristic equation (16) does not have purely imaginary roots. Since and ensure that all roots of (17) have negative real parts by Rouche’s theorem, it follows that the roots of (20) have negative real part too. This is summarized as follows.
Theorem 3. If hold, then all roots of (16) have negative real parts for all .
On the other hand, if,then (20) has a unique positive root .
Substituting into (19) and solving for , we get If , , and , then (20) has two positive roots and . Substituting into (19), we get Differentiating (16) with respect to , we obtain by using .
Thus by using .
Theorem 4. If , , and hold, then equilibrium is asymptotically stable for and unstable . Further, as increases through , bifurcates into small amplitude periodic solutions, where
Proof. For , is asymptotically stable if conditions and hold. Hence, by Butler’s lemma, remain stable for . We now have to show that
This will signify that there exists at least one eigenvalue with positive real part for . Moreover, the conditions of Hopf bifurcation (see ) are then satisfied yielding the required periodic solution.
From (24), we have Therefore,
Therefore, the transversality condition holds, and hence, Hopf bifurcation occurs at , . This completes the proof.
Theorem 5. Let be defined in (22). If , , and hold, then there is a positive integer such that there are switches from stability to instability and to stability. In other words, when , the equilibrium is stable, and when , the equilibrium is unstable. Therefore, there are bifurcations at when , .
Proof. If conditions , , and hold, then to prove the theorem, we need only to verify the transversality conditions by Cushing  and Cushing and Saleem  From (24), we have Therefore, Again, Therefore, Hence, the transversality conditions hold. This completes the proof.
5. Stability and Direction of Hopf Bifurcation
Without loss of generality, denote the critical values by , and set . Then is a Hopf bifurcation value of the system (1). Thus, we can work in the phase space . Let , .
Then system (1) is transformed into where We rewrite (34) as where , is defined by and , are given by respectively. By the Riesz representation theorem, there exists a two-dimensional matrix , , whose elements are functions of bounded variation such that In fact, we can choose where is the Dirac delta function.
For , define the operator as Then system (36) is equivalent to the operator equation of the form where and for .
For , define and define the bilinear form where . Then and are adjoint operators.
By the discussion in Section 4, we know that are eigenvalues of , and therefore they are also eigenvalues of . We first need to compute the eigenvector of and corresponding to and , respectively.
By direct computation, we can obtain that and is the eigenvector of corresponding to the eigenvalue , where Also we have , and is the eigenvector of corresponding to the eigenvalue , where Then we choose such that and .
Using the same algorithms and similar computation process as in Hassard et al. , we obtain the following quantities: where where and are constant vectors, where Therefore, we can compute the following values: Here determines the direction of the Hopf bifurcation, determines the stability of the bifurcating periodic solutions, and determines the periods of the bifurcating periodic solutions.
Based on the discussion above, we can obtain the following results.
Theorem 6. The sign of determines the direction of Hopf bifurcation: if , then the Hopf bifurcation is supercritical (subcritical); determines the stability of the bifurcating periodic solution; the bifurcating periodic solution is stable (unstable) if ; determines the period of the bifurcating periodic solution; the period increases (decreases) if .
6. Numerical Simulation
As the problem is not a case study, the real-world data are not available for this model. We, therefore, take here some hypothetical data with the sole purpose of illustrating the results that we have established in the previous sections.
(i) Let us consider the parameters of the system (2) as , , , , , , and in appropriate units. For this value of parameters we get the critical value of as . Thus it is easy to verify that for this set of parameters the system (2) is locally asymptotically stable around its interior equilibrium for and is unstable for . Thus for , the system (2) undergoes a Hopf bifurcation. Now for , we have interior equilibrium which is asymptotically stable (see Figures 2 and 3), but for , the interior equilibrium is unstable (see Figures 4 and 5). Thus taking as control parameter, it is possible to drive the prey-predator system to require equilibrium and to prevent the cycle behaviour of the system.
(ii) For the values , , , , , , , and , then the system (1) has a positive equilibrium . The conditions and are satisfied, so the positive equilibrium is asymptotically stable for (see Figure 6). The conditions , , and are satisfied, so we compute that and . Therefore, by Theorem 4, we see that the positive equilibrium is stable when (see Figures 7 and 8) and unstable when (see Figures 9 and 10). Thus the system undergoes a Hopf bifurcation at . Thus when , , we compute , , , , and . Since , , and , therefore, the Hopf bifurcation is subcritical, the bifurcating periodic solutions at are stable, and period of bifurcating periodic solutions increases.
7. Concluding Remarks
This paper deals with a delayed prey-predator model with Holling type III functional response and harvesting of predator species. Oscillatory behavior and existence of limit cycles in harvested predator-prey system are common in nature. In the absence of delay, we have proved that exactly one stable limit cycle occurs when positive equilibrium is unstable. Since harvesting is associated with economic interests, our analysis shows that the harvesting of predator species plays an important role in the shaping of the dynamical behavior of the system. From Figure 1, we see that harvesting effort affects the interior equilibrium of system (2), and, as the harvesting effort increases, the interior equilibrium may be vanished, which means that the population of predator may become extinct. It is also observed from the obtained results that harvesting effort can cause an unstable equilibrium to become stable and even a simple Hopf bifurcation occurred when the parameter passes through its critical value.
Taking delay as bifurcation parameters, we see that a Hopf bifurcation occurs, whenever delay increases a critical value. Furthermore, the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions are determined by the normal form theory and center manifold theorem.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
C. W. Clark, Mathematical Bioeconomics: The Optimal Management of Renewable Resource, John Wiley & Sons, New York, NY, USA, 2nd edition, 1990.
T. K. Kar and B. Ghosh, “Sustainability and optimal control of an exploited prey predator system through provision of alternative food to predator,” BioSystems, vol. 109, no. 2, pp. 220–232, 2012.View at: Publisher Site | Google Scholar
T. Legovic, J. Klanjšček, and S. Geček, “Maximum sustainable yield and species extinction in ecosystems,” Ecological Modelling, vol. 221, no. 12, pp. 1569–1574, 2010.View at: Publisher Site | Google Scholar
T. K. Kar and B. Ghosh, “Impacts of maximum sustainable yield policy to prey-predator systems,” Ecological Modelling, vol. 250, pp. 134–142, 2013.View at: Google Scholar
B. Ghosh and T. K. Kar, “Possible ecosystem impacts of applying maximum sustainable yield policy in food chain models,” Journal of Theoretical Biology, vol. 329, pp. 6–14, 2013.View at: Google Scholar
B. Ghosh and T. K. Kar, “Maximum sustainable yield and species extinction in a prey-predator system: some new results,” Journal of Biological Physics, vol. 39, no. 3, pp. 453–467, 2013.View at: Publisher Site | Google Scholar
T. K. Kar and H. Matsuda, “A bioeconomic model of a single-species fishery with a marine reserve,” Journal of Environmental Management, vol. 86, no. 1, pp. 171–180, 2008.View at: Publisher Site | Google Scholar
T. K. Kar and B. Ghosh, “Sustainability and economic consequences of creating marine protected areas in multispecies multi activity context,” Journal of Theoretical Biology, vol. 318, pp. 81–90, 2013.View at: Google Scholar
T. K. Kar and U. Das, “Regulation of an exploited prey-predator system: a dynamic reaction model,” International Journal of Ecological Economics and Statistics, vol. 31, no. 4, pp. 102–121, 2013.View at: Google Scholar
P. Wangersky and W. Cunningham, “Time lag in prey-predator population models,” Ecology, vol. 38, pp. 136–139, 1957.View at: Google Scholar
A. Martin and S. Ruan, “Predator-prey models with delay and prey harvesting,” Journal of Mathematical Biology, vol. 43, no. 3, pp. 247–267, 2001.View at: Google Scholar
Y. Chen and S. Changming, “Stability and Hopf bifurcation analysis in a prey-predator system with stage-structure for prey and time delay,” Chaos, Solitons & Fractals, vol. 38, no. 4, pp. 1104–1114, 2008.View at: Publisher Site | Google Scholar
X. Zhang, R. Xu, and Q. Gan, “Periodic solution in a delayed predator prey model with Holling type III functional response and harvesting term,” World Journal of Modelling and Simulation, vol. 7, no. 1, pp. 70–80, 2011.View at: Google Scholar
T. K. Kar, “Selective harvesting in a prey-predator fishery with time delay,” Mathematical and Computer Modelling, vol. 38, no. 3-4, pp. 449–458, 2003.View at: Google Scholar
T. K. Kar and H. Matsuda, “Controllability of a harvested prey-predator system with time delay,” Journal of Biological Systems, vol. 14, no. 2, pp. 243–254, 2006.View at: Publisher Site | Google Scholar
T. K. Kar and U. K. Pahari, “Non-selective harvesting in prey-predator models with delay,” Communications in Nonlinear Science and Numerical Simulation, vol. 11, no. 4, pp. 499–509, 2006.View at: Publisher Site | Google Scholar
T. K. Kar and U. K. Pahari, “Modelling and analysis of a prey-predator system with stage-structure and harvesting,” Nonlinear Analysis: Real World Applications, vol. 8, no. 2, pp. 601–609, 2007.View at: Publisher Site | Google Scholar
T. K. Kar and A. Ghorai, “Dynamic behaviour of a delayed predator-prey model with harvesting,” Applied Mathematics and Computation, vol. 217, no. 22, pp. 9085–9104, 2011.View at: Publisher Site | Google Scholar
Y. Kuang, Delay Differential Equations: With Applications in Population Dynamics, Academic Press, Boston, Mass, USA, 1993.
C. S. Holling, “The functional response of predator to prey density and its role mimicry and population regulation,” Memoirs of the Entomological Society of Canada, vol. 45, pp. 3–60, 1965.View at: Google Scholar
Y. Kuang and H. I. Freedman, “Uniqueness of limit cycles in Gause-type models of predator-prey systems,” Mathematical Biosciences, vol. 88, no. 1, pp. 67–84, 1988.View at: Google Scholar
J. Hale, Ordinary Differential Equations, John Wiley & Sons, New York, NY, USA, 1969.
J. M. Cushing, Integro-Differential Equations and Delay Models in Population Dynamics, Springer, Heidelberg, Germany, 1977.
J. M. Cushing and M. Saleem, “A predator prey model with age structure,” Journal of Mathematical Biology, vol. 14, no. 2, pp. 231–250, 1982.View at: Google Scholar
B. Hassard, N. Kazarinoff, and Y. Wan, Theory and Application of Hopf Bifurcation, Cambridge University Press, Cambridge, Mass, USA, 1981.