Stability of Prey-Predator Model with Holling type Response Function and Selective Harvesting

The use of mathematical models in prey predator interplay is common to solve the interdisciplinary natural problems. This paper reports analytical advancement of measuring selective harvesting activity of prey proportional to their population size and studied the stability of the model using Holling type functional response. In this paper, we analysed four prey-predatory model and considered prey and predator as a X and Y axis respectively followed by applied variational matrix and Holling I and II type response function for equilibrium and local stability measurement. Simulation experiments were carried out. Further, numerical analysis was done with help of MATLAB packages at MS window 7. Analysis of result showed prey and predator population converges asymptotically to their equilibrium values when t (time) tends to infinity and corresponding spiral phase portraits obtained. Interestingly analysis of result showed the behaviour of prey and predator with respect to time and phase portrait of the system near the equilibrium point. Above analysis indicated that application of vibrational matrix and holing type response function give better understand ability of prey predator interplay of biological forces. Citation: Jha PK, Ghorai S (2017) Stability of Prey-Predator Model with Holling type Response Function and Selective Harvesting. J Appl Computat Math 6: 358. doi: 10.4172/2168-9679.1000358


Introduction
The prey-predator models were first modelled by Lotka-Volterra [1][2][3][4][5]. They used simple response function proportional to the number of predators. In prey-predator models, species normally follow different growth functions and among these [6][7][8], Logistic growth function is important one, which was first used by Verhulst [2] for human growth. Later Feller [9] assumed that almost every population that increases asymptotically will fit to the Logistic growth law to some degree. There also exist some other growth functions suggested by Gompertz [10], May [11]. For prey-predator system, response functions in between prey and predator play an important role for the long term existence of the ecosystem. There are several types of response functions such as ratio dependent, Holling types response functions [12], Michaelis-Menten type, Beddington-DeAngelis [13,14] response function, etc.
The stability of ecological systems and the persistence of species within them are fundamental concerns in ecology. Mathematical models of ecological systems, reflecting these concerns, have been sued to investigate the stability of a variety of systems. For example, see [15][16][17][18][19][20][21][22][23][24]. The dynamic relationship between predator and their prey has long been and will continue to be one of the dominant themes in both ecology and mathematical ecology due to its universal existence and importance [17]. Many excellent works have been done for the Lotka-Volterra type predator-prey system but considerable uncertainties and errors in interpretation aroused out this procedure. In ref. [18], Holling proposed that there exist three functional response of the predator which usually called Holling type I, Holling type II and Holling type III [25]. He proposed the form as a Holling type II response function, it usually describes the uptake of substrate by the microorganisms in microbial dynamics kinetics [26,27]. If the predator is the invertebrate, it always the case. Freedmaan et al. emphasizes the application of holling type response function in preypredator interplay but also Chen and Kaung specifically mentioned the need of intercomparative analysis of holling response function for the better understanding of prey-predator interplay. Existing work does not cover comparative analysis of holling response function with vibrational matrix. For that in this paper, we have consider four preypredator models with Holling type-I and II predation function and selective harvesting of prey or predator. Keeping in view the practical implication of this study in understanding biological interplay. We elucidated the theoretical basis, equilibrium and local stability followed by verification and computation simulation for derived holling functions. This paper presents computation approach for real application of these mathematical models in biological analysis.

Theoretical Basis for Prey-Predator Interplay and Equilibrium Estimation
Standard theoretical model for prey-predator Theoretical model I: Consider prey-predator system with following prey-predator model with logistic growth of prey and predation function is proportional to prey density dy y xy dt Theoretical model II: Consider the following pre-predator system with Holling-II type predation function Theoretical model III: Consider the following pre-predator system with Holling-II type predation function and constant harvesting of prey Theoretical model IV: Consider the following pre-predator system with Holling-II type predation function and constant harvesting of predator Where, α is the intrinsic growth rate of prey, ϒ is the death rate of predator, β is the predation rate of prey, a is the half saturation constant, h 1 is the constant rate of harvesting rate of prey, h 2 is the constant rate of harvesting rate of predator and δ is the conversion of predator and also α, ϒ, β, a, h 1 , h 2 and δ are all positive.

Existence and determination of equilibrium between preypredator
As prey predator stability is based on existence of equilibrium so we determine the existence of equilibrium for above four models.
Where x * = δ Υ and y * = For theoretical model II: Equation (2.3) and (2.4) of a particular system was used to determine the equilibrium point. We have * For theoretical model III: Equation (2.5) and (2.6) of a particular system was used to determine the equilibrium point. We have Solving eqns. (3.6) and (3.7) we obtain the equilibrium points are (x´, 0), (x*, y*).
Where, x´= x*, y* are positive as δ > ϒ and ( ) For theoretical model IV: Equation (2.7) and (2.8) of a particular system was used to determine the equilibrium point. We have Then from eqn. (3.10) we get There are two changes of sign. When we replace x by -x, there are no changes of sign.
So, by Descarte's rule of sign there are exactly two positive roots. Say, one root be x*.

Experimental Setting and Verification for Local Stability
Local stability of these four models is discussed with variational matrix.
For theoretical model I: The variational matrix of the system eqns.
Now discuss the stability near the (0, 0) At the point (0, 0), The characteristic equation of the corresponding variational matrix is Eigen values are real distinct and opposite sign so the equilibrium pint is a saddle point and therefore the system is unstable.
At the point (k, 0), Therefore, the roots are real distinct and negative therefore the equilibrium point is a node. So the system is asymptotically stable if The characteristic equation corresponding varitional matrix is There are two changes of sign. So, by Descarte's rule of sign there are two negative roots. Therefore, the system have stable node.

For theoretical model II:
The variational matrix of the system (2.3) and (2.4) is At the point (k, 0), The characteristic equation is Therefore, the roots are real distinct and negative therefore the equilibrium point is a node. So the system is asymptotically stable if a K δ ϒ < − ϒ . Now, at the point (x*, y*), The characteristic equation is,  Hence, therefore, the roots are real distinct and negative or complex with negative real part therefore the equilibrium point is a node. So the system is asymptotically stable when The characteristic equation of the corresponding variational matrix Therefore, either λ = 0 or λ = ' 1 ' x h k x α − So, the system is not stable. Now, at the point (x*, y*), The characteristic equation of the corresponding variational matrix is  If we replace λ by -λ, then there are two changes of sign when ( ) So, by Descarte's rule of sign there are two negative roots when Hence, therefore, the roots are real distinct and negative or complex with negative real part therefore the equilibrium point is a node.
The characteristic equation of the corresponding variational matrix is Hence, therefore, the roots are real distinct and negative or complex with negative real part therefore the equilibrium point is a node. So the system is asymptotically stable when

Results and Discussion
It was observed that the bifurcation of the theoretical model- Numerical simulation of Model-4 is similar to Model-3 (Table 1).

Discussion, Conclusion and Future Extension: Four mathematical prey-predatory population
In this paper, four mathematical prey-predator population models are considered and analysed. In these models we consider the Logistic law of growth of prey. We consider Holling-I type predator response function in model-1 and in models-2,3 and 4 Holling -II type predator response function. Also consider selective harvesting of prey and predator. In model 3, consider constant prey harvesting only and constant predator harvesting in modes l-4. The conditions of existence of several equilibrium points have been examined. Local stability of the equilibrium points are discussed by variational matrix and also derive the conditions of asymptotical stability of equilibrium points. Discuss the bifurcation of the model-2 with respect to the parameter a (half saturation constant) and also derive the condition of existence of bifurcation. Numerical simulation is done by MATLAB software. The behaviour of prey and predator population with respect to time and phase portrait of the system near the equilibrium point is presented. It is observed that with the given values of parameters the prey and predator population converges asymptotically to their equilibrium values when t (time) tends to infinity and corresponding spiral phase portraits are obtained which are presented. Bifurcation of the model-2 is discussed for two values of a (half saturation constant) (i.e., a=9 and a=8). For a=9, the condition of stability is satisfied i.e. 9>8.81 and for a=8, the condition of stability is not satisfied i.e., 8<8.81 and oscillatory behaviour of the prey and predator population and periodic phase portrait are presented.
These models can be further extended by other type of response functions such as Beddington-De Angelis sigmoid response function etc., Selective harvesting by combined harvesting and also constant rate of harvesting by variable rate of harvesting. The prey-predator model can be extended by three species food chain model i.e. prey, predator and super predator.