# Global Dynamics of an Exploited Prey-Predator Model with Constant Prey Refuge

Authors: Uttam Das; T. K. Kar; U. K. Pahari
Journal: ISRN Biomathematics
 (2013)

---

## Abstract

This paper describes a prey-predator model with Holling type II functional response incorporating constant prey refuge and harvesting to both prey and predator species. We have analyzed the boundedness of the system and existence of all possible feasible equilibria and discussed local as well as global stabilities at interior equilibrium of the system. The occurrence of Hopf bifurcation of the system is examined, and it was observed that the bifurcation is either supercritical or subcritical. Influences of prey refuge and harvesting efforts are also discussed. Some numerical simulations are carried out for the validity of theoretical results.


---

## Body

## 1. Introduction

Prey-predator models are of great interest to researchers in mathematics and ecology because they deal with environmental problems such as community’s morbidity and how to control it and optimal harvest policy to sustain a community. In the physical sciences, generic models can be constructed to explain a variety of phenomena. However, in the life sciences, a model only describes a particular situation. So a variety of models are needed due to the complexity of the ecosystem. Theoretical and numerical studies of these models are able to give us an understanding of the interactions that are taking place.While investigating biological phenomena, there are many factors which affect dynamical properties of biological and mathematical models. One of the familiar nonlinear factors is the functional response. In population dynamics, a functional response of the predator to the prey density refers to the change in the density of prey attached per unit time per predator as the prey density changes [1]. Holling [2] suggested three different kinds of functional response for different kinds of species to model the phenomena of predation, which made the classical Lotka-Volterra system more realistic. The other factor which affects dynamical properties is harvesting. The subject of harvesting in predator-prey systems has been of interest to economists, ecologists, and natural resource managers for some time now. There are basically three types of harvesting reported in the literature: (i) constant harvesting where a constant number of individuals are harvested per unit time, (ii) proportional harvesting where the number of individuals harvested per unit time is proportional to the current population, and (iii) nonlinear harvesting. Most research has focused attention on optimal exploitation, guided entirely by profits from harvesting. First of all, depending on the nature of the applied harvesting strategy, the long-run stationary density of population may be significantly smaller than the long-run stationary density of a population in the absence of harvesting. Therefore, while a population can, in the absence of harvesting, be free of extinction risk, harvesting can lead to the incorporation of a positive extinction probability and therefore, to potential extinction in finite time. Secondly, if a population is subjected to a positive extinction rate, then harvesting can drive the population density to a dangerously low level at which extinction becomes a sure thing no matter how the harvester affects the population afterwards.Effects of harvesting on various types of prey-predator models have been considered by many researchers. Mathematical modeling with harvesting renewable resources started with the studies of Clark [3, 4]. Xiao et al. [5] have investigated the dynamics of a system with constant harvesting on the predator, whereas Leard et al. [6] have studied the dynamics of ratio-dependent models that include nonconstant harvesting on the prey. Xia et al. [7] studied the combined effects of harvesting and time delay. Recently, Pahari and Kar [8] analyzed a prey-predator fishery model with harvesting where they applied a tax to regulate the fishery.The existence of prey refuges can clearly have important effects on the coexistence of predators and their prey. According to Taylor [9], the different kinds of refuges can be arranged into three types: (a) those which provide permanent spatial protection for a small subset of the prey population, (b) those which provide temporary spatial protection, and (c) those which provide a temporal refuge in numbers, that is, decrease the risk of predation by increasing the abundance of vulnerable prey. It is also observed that refuge has a stabilizing effect on the equilibrium for a simple Lotka-Volterra model. The problem of predator-prey interactions under a prey refuge has been studied by some authors. McNair [10] showed that several kinds of refuges could exert a locally destabilizing effort and create stable, large-amplitude oscillations which would damp out if no refuge was present. McNair [11] obtained that a prey refuge with legitimate entry existing dynamics was quite capable of amplifying rather than damping predator-prey oscillation. Kar [12] proposed a predator-prey model incorporating a prey refuge and independent harvesting on either species. He showed that, using the harvesting efforts as control, it is possible to break the cyclic behavior of the system. Wang et al. [13] considered a prey-predator system where individuals from a prey fish population could hide in holes where predators were unable to enter. Ji and Wu [14] considered a predator-prey model incorporating a constant prey refuge and a constant-rate prey harvesting and showed that the system is controlled by using constant harvesting or prey refuge. Prepredator interactions with prey refuge also may be found in the work of Huang et al. [15], Wang and Pan [16], Kar et al. [17], and so forth.The objective of this paper is to study the combined effects of harvesting and constant prey refuge on the dynamics of predator-prey model. We also examined the existence of a Hopf bifurcation and determined the conditions for which the bifurcation is either supercritical or subcritical.González-Olivares and Ramos-Jiliberto [18] studied the following prey-predator model with constant number of prey using refuges as shown in the Figure 1:
(1)dxdt=αx(1-xk)-β(x-m)y1+a(x-m),dydt=-dy+cβ(x-m)y1+a(x-m),
with x(0)>m, y(0)>0, where x and y denote the prey and predator populations, respectively, at any time t and α and k, respectively, represent the intrinsic growth rate and carrying capacity of the prey.  m is a constant, which represents number of prey which seeks refuge from predation. d is the death rate of the predator and c is the conversion factor. The term βx∣(1+ax) denotes the functional response of the predator, which is termed Holling type II response function (see [2]). They examined the local stability of equilibria and existence of limit cycle. Chen et al. [19] considered the same model and examined also global stability and uniqueness of limit cycle. Also Ji and Wu [14] considered the prey-predator model with constant prey refuge and a constant-rate prey harvesting h>0,
(2)dxdt=αx(1-xk)-β(x-m)y1+a(x-m)-h,dydt=-dy+cβ(x-m)y1+a(x-m),
and studied instability and global stability of the equilibria and uniqueness of limit cycle and showed also the influence of constant prey refuge and constant-rate prey harvesting.Figure 1
Conceptual diagram of the Rosenzweig-MacArthur model with prey refuge. Small quadrilateral and oval shapes represent the prey and predator, respectively. The vertical dashed line represents the boundary of a refuge, where prey species are invulnerable to predators. The refuge could protect either a constant number or a constant fraction of the total prey population.Motivated by the paper of González-Olivares and Ramos-Jiliberto [18], we considered the following prey-predator system:
(3)dxdt=αx(1-xk)-β(x-m)y1+a(x-m)-q1E1x,dydt=-dy+cβ(x-m)y1+a(x-m)-q2E2y,
where E1≥0 and E2≥0 denote the harvesting efforts for the prey and predator, respectively. q1E1x and q2E2y represent the catch of the prey and predator population, where q1 and q2 represent the catchability coefficients, respectively. α, β, k, d, a, c, m, q1, and q2 are positive constants.Our paper is organized in the following way. In Section2, we have discussed the boundedness of the solutions of system (3). All possible equilibria of the system (8) and the stability criterion at those equilibria and uniqueness of limit cycles at interior equilibrium are discussed in Section 3. In Section 4, we have discussed Hopf bifurcation at interior equilibrium. Influence of refuge parameter m and harvesting efforts E1 and E2 are discussed in Section 5. Numerical simulations are given in Section 6. A brief concluding remark is given in Section 7.


## 2. Boundedness of the System

Boundedness of a model guarantees its validity. The following theorem establishes the uniform boundedness of the model system (3).Theorem 1.
All solutions of the system (3) which start in R+2 are uniformly bounded.Proof.
Let(x(t),y(t)) be any solution of the system with positive initial conditions x(0)>m and y(0)>0.
Now we define the functionW=x+y/c.
Therefore, time derivative gives(4)dWdt=dxdt+1cdydt=αx-αx2k-q1E1x-dyc-q2E2yc.
Now for each v>0, we have
(5)dWdt+vW=αx-αx2k-q1E1x+vx+vyc-dyc-q2E2yc.
We have dW/dt+vW≤(k/4α)(α+v-q1E1)2-(1/c)(d+q2E2-v).
Let us choosev>d+q2E2; then the right-hand side is positive. As we assume that both E1 and E2 are bounded, the right-hand side is bounded for all (x,y)∈R+2.
Thus we choose au>0 such that dW/dt+vW<u.
Applying the theory of differential inequality [20], we obtain
(6)0<W(x,y)<uv(1-e-vt)+W(x(0),y(0))evt,fort⟶∞,0<W<uv.
Thus, all solutions of the system (3) that start in R+2 are confined to the region B, where
(7)B={(x,y)∈R+2:W=uv+ε,foranyε>0}.
This completes the theorem.Now for simplicity, let us introduceX=x-m; then the system (3) of equations changes to (still denote X=x)
(8)dxdt=α(x+m)(1-x+mk)-βxy1+ax-q1E1(x+m),dydt=-dy+cβxy1+ax-q2E2y.


## 3. 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 (8) as follows:(i)
equilibrium in the absence of predatorP1(x1,0), where x1=(kα-mα-q1E1k)/α exists provided kα>mα+q1E1k,(ii)
the interior equilibriumP*(x*,y*), where
(9)x*=d+q2E2βc-ad-aq2E2,y*=1+ax*βkx*(α(x*+m)(k-m-x*)-q1E1k(x*+m)).If βc>ad+aq2E2 holds and any of the following conditions hold:(H1)
A>0, B>0, A2≥4B, and (A-A2-4B)/2<m<(A+A2-4B)/2,(H2)
A>0, B<0, and 0<m<(A+A2-4B)/2,(H3)
A<0, B<0, and 0<m<(A+A2-4B)/2,where A=k-2x*-q1E1k/α and B=(x*)2-kx*+q1E1kx*/α, then the interior equilibrium P*(x*,y*) exists.For the stability analysis of the equilibria, let the variational matrix of the system (8) at an arbitrary point (x,y) be M(x,y).AtP1(x1,0), the eigenvalues of the variational matrix M(x1,0) are E1q1-α and cβx1/(1+ax1)-d-q2E2.Now if(α/q1){(k-m)/k-(d+q2E2)/k(cβ-ad-aq2E2)}<E1<α/q1 and βc>ad+aq2E2, then both eigenvalues are negative and hence P1(x1,0) is locally asymptotically stable.To determine the local stability character of the interior equilibriumP*(x*,y*), we compute the variational matrix M(x*,y*) at (x*,y*).The characteristic equation of the variational matrixM(x*,y*) is given by
(10)γ2+γ(2mαk-M1)+M2=0,
where
(11)M1=α-q1E1-βy*(1+ax)2-2x*αk,M2=cβ2x*y*(1+ax*)3>0.
Routh-Hurwitz criterion states that all roots of the characteristic equation (10) have negative real parts if (2mα/k-M1)>0.Therefore, interior equilibriumP*(x*,y*) will be asymptotically stable if m>M1k/2α, and unstable if m<M1k/2α.Theorem 2.
IfP*(x*,y*) exists with m>M1k/2α then P*(x*,y*) is locally asymptotically stable.


### 3.1. Global Stability for the Interior Equilibrium

Theorem 3.
IfR(x)>0, then the system (8) will be globally asymptotically stable around the interior equilibrium P*(x*,y*).Proof.
To show the global stability of system (8), we define a Lyapunov function as follows:
(12)V(x,y)=A1(x-x*-x*lnxx*)+B1(y-y*-y*lnyy*),
where A1 and B1 are positive constants to be determined in the subsequent steps.
The time derivative along the trajectories of (8) is
(13)dVdt=A1x-x*xdxdt+B1y-y*ydydt=A1(x-x*)[βy1+axq1E1mxα(1+mx)(1-x+mk)-βy1+ax-q1E1-q1E1mx(1-x+mk)]+B1(y-y*)(-d-q2E2+cβx1+ax)=-A1(x-x*)2x*[αx*k+αmx(1-mk)-q1E1mx-aβx*y*(1+ax)(1+ax*)]+β[cB1-A1(1+ax*)](x-x*)(y-y*)(1+ax)(1+ax*)=-A1(x-x*)2x*R(x)+β[cB1-A1(1+ax*)](x-x*)(y-y*)(1+ax)(1+ax*),
where R(x)=[αx*/k+(αm/x)(1-m/k)-q1E1m/x-aβx*y*/(1+ax)(1+ax*)].
ChoosingA1=1 and B1=(1+ax*)/c, we have
(14)dVdt=-(x-x*)2x*R(x).
Thus if R(x)>0, then dV/dt<0. This completes the proof.In Figure2, we show that R(x) may be positive for some positive value of x by using the numerical value of the parameters α=8, a=0.5, β=4, c=0.8, d=2, k=50, q1=0.6, q2=0.2, E1=0.4, E2=0.2, and m=0.75.Figure 2
Graph ofR(x).


### 3.2. 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 (8) in the form
(15)dxdt=xg(x)-yh(x),x(0)>0,dydt=y(-d-q2E2+q(x)),y(0)>0,
where g(x)=α(1-x/k)-2αm/k-q1E1+(1/x)(αm-αm2/k-q1E1m),h(x)=βx/(1+ax), and q(x)=cβx/(1+ax).Now we consider the following theorem (see [21]) regarding uniqueness of limit cycles of the previous system.Theorem 4.
Suppose for system (15) that
(16)ddx(xg′(x)+g(x)-xg(x)(h′(x)/h(x))-d-q2E2+q(x))≤0
in 0≤x<x* and x*<x≤k. Then the previous system has exactly one limit cycle which is globally asymptotically stable with respect to the set
(17){(x,y):x>0,y>0}∖{P*(x*,y*)}.Following Theorem4, we may state that when m<M1k/2α, system (8) has unique globally stable limit cycle. Thus we see that when the system is unstable, there exists a unique globally stable limit cycle.


## 3.1. Global Stability for the Interior Equilibrium

Theorem 3.
IfR(x)>0, then the system (8) will be globally asymptotically stable around the interior equilibrium P*(x*,y*).Proof.
To show the global stability of system (8), we define a Lyapunov function as follows:
(12)V(x,y)=A1(x-x*-x*lnxx*)+B1(y-y*-y*lnyy*),
where A1 and B1 are positive constants to be determined in the subsequent steps.
The time derivative along the trajectories of (8) is
(13)dVdt=A1x-x*xdxdt+B1y-y*ydydt=A1(x-x*)[βy1+axq1E1mxα(1+mx)(1-x+mk)-βy1+ax-q1E1-q1E1mx(1-x+mk)]+B1(y-y*)(-d-q2E2+cβx1+ax)=-A1(x-x*)2x*[αx*k+αmx(1-mk)-q1E1mx-aβx*y*(1+ax)(1+ax*)]+β[cB1-A1(1+ax*)](x-x*)(y-y*)(1+ax)(1+ax*)=-A1(x-x*)2x*R(x)+β[cB1-A1(1+ax*)](x-x*)(y-y*)(1+ax)(1+ax*),
where R(x)=[αx*/k+(αm/x)(1-m/k)-q1E1m/x-aβx*y*/(1+ax)(1+ax*)].
ChoosingA1=1 and B1=(1+ax*)/c, we have
(14)dVdt=-(x-x*)2x*R(x).
Thus if R(x)>0, then dV/dt<0. This completes the proof.In Figure2, we show that R(x) may be positive for some positive value of x by using the numerical value of the parameters α=8, a=0.5, β=4, c=0.8, d=2, k=50, q1=0.6, q2=0.2, E1=0.4, E2=0.2, and m=0.75.Figure 2
Graph ofR(x).


## 3.2. 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 (8) in the form
(15)dxdt=xg(x)-yh(x),x(0)>0,dydt=y(-d-q2E2+q(x)),y(0)>0,
where g(x)=α(1-x/k)-2αm/k-q1E1+(1/x)(αm-αm2/k-q1E1m),h(x)=βx/(1+ax), and q(x)=cβx/(1+ax).Now we consider the following theorem (see [21]) regarding uniqueness of limit cycles of the previous system.Theorem 4.
Suppose for system (15) that
(16)ddx(xg′(x)+g(x)-xg(x)(h′(x)/h(x))-d-q2E2+q(x))≤0
in 0≤x<x* and x*<x≤k. Then the previous system has exactly one limit cycle which is globally asymptotically stable with respect to the set
(17){(x,y):x>0,y>0}∖{P*(x*,y*)}.Following Theorem4, we may state that when m<M1k/2α, system (8) has unique globally stable limit cycle. Thus we see that when the system is unstable, there exists a unique globally stable limit cycle.


## 4. Hopf Bifurcation

To discuss Hopf bifurcation of the system (8) we take the help of the paper [22].Let us now consider system (8) in the form
(18)dxdt=h(x)(f(x)-y),dydt=cy(h(x)-d′),
where f(x)=(α/kβ)(k-x)(1+ax)-(2αm/kβ)(1+ax)-(q1E1/β)(1+ax)+(1/xβ)(αm-αm2/k-q1E1m)(1+ax),h(x)=βx/(1+ax), and d′=(d+q2E2)/c.Let(δ,yδ) be the interior equilibrium of the system (18), where
(19)δ=d+q2E2βc-ad-aq2E2,yδ=1+aδβkδ(α(δ+m)(k-m-δ)-q1E1k(δ+m)).
In this section, we examined the Hopf bifurcation occurring at (δ-,yδ-) taking δ as a bifurcation parameter.Variational matrix of the system (18) at (δ,yδ) is
(20)J=(W(δ)Y(δ)Z(δ)0),
where W(δ)=h(δ)f′(δ), Y(δ)=-h(δ), and Z(δ)=cf(δ)h′(δ).The characteristic equation of the variational matrixJ is given by μ2-μT+D=0, where T=trJ=W(δ) and D=detJ=-Y(δ)Z(δ).Now ifT<0 and D>0, then (δ,yδ) is locally asymptotically stable. On the other hand, if T>0 and D>0 then (δ,yδ), is unstable. Also T(δ-)=0 and the Jacobian matrix J(δ-) has a pair of imaginary eigenvalues μ=±i-Y(δ-)Z(δ-). Let μ=β(δ)±iω(δ) be the roots of μ2-μT+D=0 when δ is near δ-; then
(21)β(δ)=A(δ)2,ω(δ)=-4Y(δ)Z(δ)-{W(δ)}22,β′(δ)={h(δ)f′(δ)}′2={h′(δ)f′(δ)+h(δ)f′′(δ)}2.δ- is a root of f′(δ)=0  (as h(δ)=d′≠0); that is δ=δ-, satisfies the equation
(22)2αaδ3+(2αam+q1E1ak+α-aαk)δ2+(αmk-αm2-q1E1km)=0,
and also f″(δ-)≠0.Therefore,β′(δ)|δ=δ-=h(δ-)f′′(δ-)/2≠0 as h(δ-)≠0.The system (18) undergoes a Hopf bifurcation at (δ-,yδ-) when δ=δ-. By further analysis, we determined that the bifurcation is either supercritical or subcritical by the first Lyapunov coefficient (see [23–25])
(23)a(δ-)=116h′(δ-)[f′′′(δ-)h(δ-)h′(δ-)+2{h′(δ-)}2f′′(δ-)-f′′(δ-)h′′(δ-)h(δ-){h′(δ-)}2].
The computation of a(δ-) is technical, and the detailed calculations are given in the appendix and only the results are stated in the following theorem.Theorem 5.
The system (18) undergoes a Hopf bifurcation at (δ-,yδ-); the Hopf bifurcation is supercritical and backward (subcritical and forward, resp.) if a(δ-)<0(a(δ-)>0), where a(δ-) is defined in (23).


## 5. Influence of the Parameterm and Harvesting Efforts E1 and E2

Interior equilibrium of the system (8) is P*(x*,y*), where
(24)x*=d+q2E2βc-ad-aq2E2,y*=1+ax*βkx*(α(x*+m)(k-m-x*)-q1E1k(x*+m)).
Notice that the scaling from system (3) to system (8) is X=x-m. Let the interior equilibrium of the system (3) be P′′(x′′,y′′), where
(25)x′′=d+q2E2βc-ad-aq2E2+m,y′′=1+a(x′′-m)βk(x′′-m)(αx′′(k-x′′)-q1E1kx′′).
The model (3) is a one-prey-one-predator model with constant prey refuge. Now we discuss the influence of the prey refuge on the model dynamics.Nowdx′′/dm=1>0.Therefore,x′′ is a strictly increasing function of m.That is, increasing the amount of prey refuge can increase prey population.To check the effect ofm on predator, differentiating y′′ with respect to m, we get dy′′/dm=-(2α(1+ax*)/βx*k)(m+x*-k/2+q1E1k/2α), where x*=(d+q2E2)/(βc-ad-aq2E2).Case 1.
Ifm>k/2-x*-q1E1k/2α, then y′′ is a monotonically decreasing function of m.
That is, if the value of prey refugem gradually increases above threshold value k/2-x*-q1E1k/2α, then the predator population gradually decreases.Case 2.
Ifm<k/2-x*-q1E1k/2α, then y′′ is a monotonically increasing function of m.
That is, if the value of prey refugem gradually decreases below threshold value k/2-x*-q1E1k/2α, then the predator population gradually increases.Case 3.
If=k/2-x*-q1E1k/2α, theny′′ has a maximum value.To construct Table1 and Figure 3, we take α=8, a=0.5, β=4, c=0.8, d=2, k=50, q1=0.6, q2=0.2, E1=0.4, and E2=0.2 in appropriate units. From the analysis shown we see that increasing the amount of prey refuge can increase prey population and that increasing the amount of prey refuge can increase the density of predator species, and this happened due to predator species still having enough food for predation with m being small, but if the prey refuge is larger than a threshold, that is, as the prey refuge becomes large enough, then the increasing amount of prey refuge can decrease predator species, and this happened due to the loss of food for predator species.Table 1
m
x
′
′
y
′
′
0.1
1.03578
3.0847
0.2
1.13578
3.37539
0.398969
1.33578
3.95004
0.6
1.53578
4.5256
0.75
1.68578
4.95175
1
1.93578
5.65572
3
3.93578
11.0052
10
10.93578
25.7753
20
20.93578
36.2088
23.31422*
24.25002
36.8980
24
24.93578
36.8685
25
25.93578
36.7197
30
30.93578
34.0934
35
35.93578
28.3297
40
40.93578
19.4288
47
47.93578
7.39071
47.5
48.43578
0.195172
*At m = 23.31422, the value of y′′ is maximum.Figure 3
Change ofx′′ and y′′ with the change of the refuge parameter m.Influence of the Harvesting Efforts E
1
and E
2
. Interior equilibrium of the system (8) is P*(x*,y*), where
(26)x*=d+q2E2βc-ad-aq2E2,y*=1+ax*βkx*(α(x*+m)(k-m-x*)-q1E1k(x*+m)).
At first let m and E2 be fixed.It is observed thatx* is independent of the parameter E1, whereas the value of y* depends on E1. Therefore E1 has no effect on the interior equilibrium level on x.To check the effect ofE1 on predator, differentiating y* with respect to E1, we get
(27)dy*dE1=-q1(1+ax*)(m+x*)βx*<0.This shows thaty* is a strictly decreasing function of E1.To construct Table2 and Figure 4, we take α=8, a=0.5, β=4, c=0.8, d=2, k=50, q1=0.6, q2=0.2, m=0.35, and E2=0.2 in appropriate units. From the analysis shown we see that as the harvesting effort E1 increases, the prey species remain unchanged but predator species decrease and this happens due to loss of food for predator species and goes to extinction when E1 is large.Table 2
E
1
x
*
y
*
1
0.93578
3.62755
3
0.93578
3.02248
5
0.93578
2.4174
6
0.93578
2.11487
7.62878
0.93578
1.6221
9.6
0.93578
1.02574
11
0.93578
0.602186
12
0.93578
0.29965
12.5
0.93578
0.148382
12.9
0.93578
0.0273671Figure 4
Change ofy* with the change of the harvesting effort E1.Now letm and E1 be fixed; then we have
(28)dx*dE2=βcq2(βc-ad-aE2q2)2>0.
Therefore x* is a strictly increasing function of E2.Also(29)dy*dE2=[αm2(x*)2+q1E1km(x*)2+2amα+akα-α-2αax*-q1E1ka-kαm(x*)2αm2(x*)2q1E1km(x*)2]dx*dE2.
Since dx*/dE2>0, therefore
(30)dy*dE2>0ifαm2(x*)2+q1E1km(x*)2+2amα+akα>α+2aαx*+q1E1ka+kαm(x*)2,dy*dE2<0ifαm2(x*)2+q1E1km(x*)2+2amα+akα<α+2aαx*+q1E1ka+kαm(x*)2.
Thus as the harvesting effort E2 increases, the predator population increases when E2 is smaller than the threshold value. But if the harvesting effort E2 gradually increases above the threshold value, that is, as the harvesting effort E2 becomes large enough, then the increasing amount of the harvesting effort E2 can decrease predator population.To construct Table3 and Figure 5, we take α=8, a=0.5, β=4, c=0.8, d=2, k=50, q1=0.6, q2=0.2, E1=0.5, and m=0.35 in appropriate units. From the analysis shown we see that as the harvesting effort E2 increases, the prey and predator species increase. But if the harvesting effort E2 increases larger than a threshold value, that is, as the harvesting effort E2 becomes large enough, then increasing amount of the harvesting effort E2 can increase the prey species but will decrease predator species and go to extinction of the predator species when E2 is large.Table 3
E
2
x
*
y
*
1
1.04762
3.79968
5
1.76471
4.15139
10
3.33333
5.23819
15
7.14286
7.79394
17
10.8
9.77236
19.4068
22.6799
12.578
19.6
24.6667
12.4993
20
30.00
11.5087
20.5
40.6667
6.11798
20.7
47.2308
0.539829Figure 5
Variation ofx* and y* with the change of the harvesting effort E2.


## 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. Let us consider the parameters of the system asα=8, a=0.5, β=4, c=0.8, d=2, k=50, q1=0.6, q2=0.2, E1=0.4, and E2=0.2 in appropriate units. For these value of parameters, we get the critical value of m as m*=0.398969. Thus it is easy to verify that for this set of parameters, the system (8) is locally asymptotically stable around its interior equilibrium P*(x*,y*) for m>m* and is unstable for m<m*. Thus for m=m*=0.398969 the system (8) undergoes a Hopf bifurcation. Now for m=0.75, we have interior equilibrium (0.93578,4.95175) which is asymptotically stable (see Figure 6), but m=0.2, and the interior equilibrium (0.93578,3.37539) is unstable (see Figure 7). Thus taking m as a control parameter, it is possible to drive the prey-predator system to require equilibrium and to prevent the cycle behaviour of the system. From Figures 8, 9, 10, and 11, we see that E1 and E2 may also be used as controls for the system (8); Hopf bifurcation occurs when E1=E1*=7.62878 (here E2=0.2, m=0.35) and E2=E2*=19.4068 (here E1=0.5, m=0.35). For the previous values of parameters and m=0.398969, we obtained one value of δ say δ-=0.93578 and a(δ-)=-0.332712<0. Thus we may conclude that the Hopf bifurcation around the interior equilibrium is supercritical and backward.Figure 6
Phase space trajectories corresponding to different initial levels, which shows that(0.93578,4.95175) is a global attractor.Figure 7
There is a stable limit cycle surrounding(0.93578,3.3754) with m=0.2.Figure 8
Phase space trajectories corresponding to different initial levels. HereE1=9.6, E2=0.2, and m=0.35.Figure 9
There is a stable limit cycle surrounding(0.93578,2.11487) with E1=6, E2=0.2, and m=0.35.Figure 10
Phase space trajectories corresponding to different initial levels withE1=0.5, E2=19.6, and m=0.35.Figure 11
There exist Hopf-bifurcating small amplitude periodic solutions withE1=0.5, E2=17, and m=0.35.


## 7. Concluding Remarks

This paper deals with a prey-predator model with Holling type II functional response incorporating a constant prey refuge and independent harvesting in either species. Oscillatory behavior and existence of limit cycles in harvested prey-predator system are common in nature. It is noted that constant prey refuge plays an important role in the dynamics of the proposed model system. It is also observed from the obtained results that constant prey refuge can cause an unstable equilibrium to become stable, and even a simple Hopf bifurcation occurred when the parameterm passes through its critical value. There exists a threshold value of m, such that for the prey refuge smaller than this threshold, increasing the amount of prey refuge can increase the predator population, and if the prey refuge is larger than the threshold, increasing the amount of prey refuge can decrease the predator population. We have proved that exactly one stable limit cycle occurs when the positive equilibrium is unstable. We also determined the critical value of δ at which Hopf bifurcation occurs and observed that the bifurcation is supercritical and backward. It was also found that it is possible to control the system in such a way that the system approaches a required state, using the efforts E1 and E2 as controls.Our analytical results and numerical simulation also indicate that dynamic behavior of the model not only depends on the prey refuge parameterm but also depends on harvesting efforts E1 and E2. Hence it is possible to control the system in such a way that the system approaches a required state, using the harvesting efforts E1 and E2 or prey refuge m as controls.In our model, we have considered the catch-rate function based on catch-per-unit-effort hypothesis; that is,h=qEu (u and E denote the prey or predator population and harvesting effort, resp.). But this type of catch-rate function embodies some defects in that (i) it assumes random search for fish, (ii) it assumes equal likelihood of being captured for every fish, (iii) there is unbounded linear increase of h with respect to E for fixed u, and (iv) there is unbounded linear increase h with respect to u for a fixed E. These unrealistic features can largely be removed by adopting the alternative functional form h=qEu/(n1E+n2u), where n1 and n2 are positive constants, but we leave it for our future research work. The entire study of the paper is mainly based on the deterministic framework. On the other hand, it will be more realistic if it is possible to consider the model system in the stochastic environment due to some ecological fluctuations and other factors. Thus, a future research problem would be considered in stochastic environment.



---

*Source: 637640-2013-09-15.xml*