Medical, Pharma, Engineering, Science, Technology and Business

^{1}School of Medical Science and Technology, Indian Institute of Technology, Kharagpur, West Bengal, India

^{2}National Institute of Animal Welfare, Ministry of Environment Forest and Climate Change, Government of India, Faridabad, New Delhi, India

- *Corresponding Author:
- Pradeep K Jha

Senior Research Scientist

School of Medical Science and Technology

Indian Institute of Technology, Kharagpur

West Bengal, India

**Tel:**03222 255 221

**E-mail:**[email protected]

**Received Date:** July 07, 2017; **Accepted Date:** July 31, 2017; **Published Date:** August 07, 2017

**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

**Copyright:** © 2017 Jha PK, et al. This is an open-access article distributed under
the terms of the Creative Commons Attribution License, which permits unrestricted
use, distribution, and reproduction in any medium, provided the original author and
source are credited.

**Visit for more related articles at** Journal of Applied & Computational Mathematics

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

Holling type-I and II predation function; Selective harvesting; Prey-predator; Variational matrix

The prey-predator models were first modelled by Lotka-Volterra [1-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-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-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.

**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

(2.1)

(2.2)

**Theoretical model II:** Consider the following pre-predator system
with Holling-II type predation function

(2.3)

(2.4)

**Theoretical model III:** Consider the following pre-predator system
with Holling-II type predation function and constant harvesting of prey

(2.5)

(2.6)

**Theoretical model IV:** Consider the following pre-predator system
with Holling-II type predation function and constant harvesting of
predator

(2.7)

(2.8)

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,
h1 is the constant rate of harvesting rate of prey, h2 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 prey-predator**

As prey predator stability is based on existence of equilibrium so we determine the existence of equilibrium for above four models.

**For theoretical model I:** To find the equilibrium points of the
system eqns. (2.1) and (2.2) we’ve,

(3.1)

(3.2)

Solving eqns. (3.1) and (3.2) we get the equilibrium points are (0, 0), (k, 0), (x*,y*).

**For theoretical model II:** Equation (2.3) and (2.4) of a particular
system was used to determine the equilibrium point. We have

(3.3)

(3.4)

Solving eqns. (3.3) and (3.4) we get the equilibrium points are (0, 0), (k, 0), (x*, y*).

Where and and x* is
positive when *δ* > ϒ (3.5)

**For theoretical model III:** Equation (2.5) and (2.6) of a particular
system was used to determine the equilibrium point. We have

(3.6)

(3.7)

Solving eqns. (3.6) and (3.7) we obtain the equilibrium points are (x´, 0), (x*, y*).

Where,

and and

x*, y* are positive as *δ* > ϒ and (3.8)

**For theoretical model IV:** Equation (2.7) and (2.8) of a particular
system was used to determine the equilibrium point. We have

(3.9)

(3.10)

From eqn. (3.9) we get, either *x* = 0 or (3.11)

From eqn. (3.11) we get, .

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*.

Therefore, the equilibrium point is (x*, y*).

Local stability of these four models is discussed with variational matrix.

**For theoretical model I:** The variational matrix of the system eqns.
(2.1) and (2.2) is

Where,

g(x, y) = −Υ*y* +*δxy*.

Therefore,

Now discuss the stability near the (0, 0)

At the point (0, 0),

The characteristic equation of the corresponding variational matrix is

Therefore, *λ* =*α* > 0

*λ* = −*γ*

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),

The characteristic equation corresponding variational matrix is

Therefore, *λ* = −*α* < 0

*λ* = −(*γ* −*δ K*) < 0 when and *λ* = −(*γ* −*δ K*) > 0 when .

Therefore, the roots are real distinct and negative therefore the equilibrium point is a node. So the system is asymptotically stable if .

Now at the point (x*, y*),

The characteristic equation corresponding varitional matrix is

Or,

Now

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

Where

Therefore,

At the point (0, 0),

So, the characteristic equation is

Therefore, *λ* = *α* > 0 and *λ* = − ϒ < 0

Hence, the system is unstable.

At the point (k, 0),

The characteristic equation is

Therefore, *λ* = −*α* < 0 and when.

Therefore, the roots are real distinct and negative therefore the equilibrium point is a node. So the system is asymptotically stable if .

Now, at the point (x*, y*),

The characteristic equation 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. So the system is asymptotically stable when . (A)

**For theoretical model III:** The vibrational matrix of the system
eqns. (2.5) and (2.6) is

Where ,

Therefore,

At the point (x´, 0),

The characteristic equation of the corresponding variational matrix is

Or,

Therefore, either λ = 0 or

So, the system is not stable.

Now, at the point (x*, y*),

The characteristic equation of the corresponding variational matrix is

Or,

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.

So the system is asymptotically stable if

For Theoretical model III: The variational matrix of the system (2.7) and (2.8) is

Where

Therefore,

Now, at the point (x*, y*),

The characteristic equation of the corresponding variational matrix is

Or,

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. So the system is asymptotically stable when

It was observed that the bifurcation of the theoretical model- II with respect to the parameter a half saturation constant. The roots of characteristic equation are real distinct negative or complex with negative real part if (From A). So the system is asymptotically stable i.e..

For theoretical model-III the system is asymptotically stable if

For theoretical model-IV the system is asymptotically stable if .

Numerical example: In this section numerical simulation is present to illustrate the behaviour of prey and predator with respect to time and phase portrait of the system near the equilibrium point which results obtained in previous sections. Now choose the values of parameters for different model.

For Model-I: α= 1.8;K=20; β=0.03; δ=0.05; γ=0.3; (**Figure 1A**).

For Model-2: α=1.8; K=10; β=0.6; δ=0.9;γ=0.1; a=9; for these values and the condition of stability with respect to the
parameter a is satisfied as 9>8.81 (**Figure 1B**).

For Model-IV α=1.8; K=10; β=0.6; δ=0.9; γ =0.1; a= 8; for these
values and the condition of stability with respect to
the parameter a is not satisfied as 8 < 8.81( 8 is not greater than 8.86)
(**Figure 2A and 2B**).

**For Model-IV:**

α=1.8;K=10; β=0.5; δ=0.6; γ =0.1; α = 10 h1=10; (**Figures 3 and 4**).

Numerical simulation of Model-4 is similar to Model-3 (**Table 1**).

Model | at (0,0) | At (x,0) | at (x,y) | Phase portrait |
---|---|---|---|---|

Model-1 | Unstable | Asymptotically stable | Stable | Figure 1A and 1B |

Model-2 | Unstable | Asymptotically stable | Asymptotically stable if | Figures 2A and 2B |

Model-3 | - | Unstable | Asymptotically stable if | Figures 3 and 4 |

Model-4 | - | - | Asymptotically stable if | Same as model-3 |

**Table 1:** Nature of equilibria of four models.

**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.

- Malthus TR (1798) An Essay on the principles of populations. St. Paul’s London.
- Verhulst PF (1838) Notice sur la loi que la population suit dans son accoroissemnt.Correspondence Mathematique et Physique 10: 113-121.
- Brauer F, Sanchez DA (1975) Constant rate population harvesting: equilibrium and stability. Theoritical Population Biology 8: 12-30.
- Lotka AJ (1925) Elements of Physical Biology. Baltimore: Williams and Willkins.
- Volterra V (1926) Variazioni e fluttazionidelnumero di individui in species animaliconviventi. MemorieAccademiaNazionaledeiLincei 2: 31-313.
- Gause GF (1934) The Struggle for Existence. Dover Phoenix Editions. New York: Hafner.
- Holling CS (1959) Some characteristic of simple types of predation and parasitism. Canadian Entomologist 91: 385-398.
- Xiao D, Ruan S (2001) Global dynamics of a ratio-dependent predator-prey system. Journal of Mathematical Biology 43: 268-290.
- Feller W (1940) On the logistic law of growth and its empirical verification biology. Acta Biotheoretica 5: 51-66.
- Gompertz B (1825) On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. Philosophical Transactions of the Royal Society 115: 513-585.
- May RM (1973) Stability and Complexity in Model Ecosystem. Monograph in Population Biology VI, Princeton: Princeton University Press.
- Holling CS (1959) Some characteristics of simple types of predation and parasitism. Canadian Entomologist 91: 385-398.
- Beddington JR, May RM (1975) Time delays are not necessarily destabilizing. Mathematical Biosciences 27: 109-118.
- DeAngelis DL, Goldstein RA, Neill RVO (1975) A model for trophic interaction. Ecology 56: 881-892.
- Reeve JD (1988) Environmental variability, migration and persistence in host-parasitoid systems. Am. Nat. 132: 810-836.
- Murdoch WW, Briggs CJ, Nisbet RM, Gurney WSC, Stewart-Oaten A (1992) Aggregation and stability in metapopulation models. Am Nat 140: 41-58.
- Berryman AA (1992) The origins and evolutions of predator-prey theory. Ecology 73: 1530-1535.
- Holling CS (1965) The functional response of predators to prey density and its role in mimicry and population regulations. Mem. Entomol. Soc Can 97: 5-60.
- Maynard SJ (1974) Models in ecology. Cambridge University Press.
- Kar TK (2005) Stability analysis of a prey-predator model incorporation a prey refuge, Commun. Nonlinear Sci. Numer. Simul 10: 681-691.
- Myerscough MR, Darwen MJ, Hogarth WL (1996) Stability, persistence and structural stability in a classical predator-prey model. Ecol. Model 89: 31-42.
- Chen FD, Chen XX (2003) The n-competing Volterra-Lotka almost periodic systems with grazing rate. J Biomath 18: 411-416.
- Chen FD (2005) Positive periodic solutions of neutral Lotka-Volterra system with feedback control. Appl. Math. Comput 162: 1279-1302.
- Peng QL, Chen LS (1994) Asymptotic behavior of the nonautonomous two-species Lotka-Volterra competition models. Computers Math. Appl 27: 53-60.
- Chen JP, Hong-de Z (1986) The qualitative analysis of two species predator-prey model with Holling type III functional response. Appl Math Mech 7: 73-80.
- Freedman HI (1980) Deterministic Mathematical Models in Population Ecology. Marcel-Dekker. New York.
- Kuang Y, Freedman HI (1988) Uniqueness of limit cycle in Gause-type predator-prey systems. Math. Biosci 88: 67-84.

Select your language of interest to view the total content in your interested language

- Adomian Decomposition Method
- Algebraic Geometry
- Analytical Geometry
- Applied Mathematics
- Axioms
- Balance Law
- Behaviometrics
- Big Data Analytics
- Binary and Non-normal Continuous Data
- Binomial Regression
- Biometrics
- Biostatistics methods
- Clinical Trail
- Complex Analysis
- Computational Model
- Convection Diffusion Equations
- Cross-Covariance and Cross-Correlation
- Differential Equations
- Differential Transform Method
- Fourier Analysis
- Fuzzy Boundary Value
- Fuzzy Environments
- Fuzzy Quasi-Metric Space
- Genetic Linkage
- Hamilton Mechanics
- Hypothesis Testing
- Integrated Analysis
- Integration
- Large-scale Survey Data
- Matrix
- Microarray Studies
- Mixed Initial-boundary Value
- Molecular Modelling
- Multivariate-Normal Model
- Noether's theorem
- Non rigid Image Registration
- Nonlinear Differential Equations
- Number Theory
- Numerical Solutions
- Physical Mathematics
- Quantum Mechanics
- Quantum electrodynamics
- Quasilinear Hyperbolic Systems
- Regressions
- Relativity
- Riemannian Geometry
- Robust Method
- Semi Analytical-Solution
- Sensitivity Analysis
- Smooth Complexities
- Soft biometrics
- Spatial Gaussian Markov Random Fields
- Statistical Methods
- Theoretical Physics
- Theory of Mathematical Modeling
- Three Dimensional Steady State
- Topology
- mirror symmetry
- vector bundle

- Total views:
**130** - [From(publication date):

September-2017 - Oct 22, 2017] - Breakdown by view type
- HTML page views :
**107** - PDF downloads :
**23**

Peer Reviewed Journals

International Conferences 2017-18