Dynamics of a harvested prey-predator model with prey refuge depended on both species

TL;DR

This paper investigates a prey-predator model with prey refuge for both species and independent harvesting, showing how refuge influences system stability and establishing optimal harvesting policies through analytical and numerical methods.

Contribution

It introduces a prey-predator model incorporating refuge for both species and analyzes how refuge affects stability and harvesting strategies.

Findings

Refuge coefficient significantly influences system stability.

Refuge can break limit cycles and achieve equilibrium.

Numerical results support theoretical analysis.

Abstract

The present paper deals with a prey-predator model with prey refuge proportion to both species and independent harvesting of each species. Our study shows that using refuge as control, it can break the limit circle of the system and reach the required state of equilibrium level. It is established the optimal harvesting policy. The boundedness, feasibility of interior equilibria, bionomic equilibrium have been determined. The main observation is that coefficient of refuge plays an important role in regulating the dynamics of the present system. Moreover the variation of the coefficient of refuge change the system from stable to unstable and vice-versa. Some numerical illustration are given in order to support of our analytical and theoretical findings.

Tables

Table 1. Table 1: Schematic representation of the feasibility and stability conditions of the equilibria of the proposed model ( 1.2 ): LAS = Locally Asymptotically Stable; GAS = Globally Asymptotically Stable

Equilibria	Feasibility conditions	Stability/Persistent    conditions	Nature
$E^0$	No condition		LAS
$E^1$			LAS
$E^{\ast }$	,	See the conditions (4.1)-(4.2) of Section 4	LAS
$E^{\ast }$	,	$4rme(1+ak)(1+a{x}^{\ast })>p{m}^2{y}^{\ast }{}^2$	GAS

Table 2. Table 2: Schematic representation of existence/feasibility conditions of the bionomic equilibria of the system ( 1.2 )

Bionomic     Equilibrium	Existence/Feasibility conditions
$(y_{\infty },E_1^{\infty })$	provided prey fishing is in operational.
$(x_{\infty },E_2^{\infty })$	$E_2^{\infty }\ge \frac{1}{q_2}\left(d-\frac{p_2{q}_{2}erk}{4c_2}\right),$ provided predator fishing is in operational.
$(x_{\infty },y_{\infty },E_1^{\infty },E_2^{\infty })$	$\frac{r}{q_1}\left(1-\frac{c_1}{p_1{q}_{1}k}\right)>\frac{p\left(p_2{q}_2-m{c}_2\right)c_2}{p_2{q}_2\left(p_1{q}_1{p}_2{q}_2+a{c}_1\left(p_2{q}_2-m{c}_2\right)\right)},$ ,  provided  both species fishing are in operational.

Table 3. Table 3: The values of the components of equilibrium point x ∗ superscript 𝑥 x^{*} and y ∗ superscript 𝑦 y^{*} corresponding to the value of m 𝑚 m .

No.	Fixed Parameters	$m$	$x^{\ast },y^{\ast }$	Figure 5
1	$r_1=3,a=0.008,d=0.04,p=0.2,q_1=0.2,q_2=0.6,E_1=2,E_2=2,k=500,e=0.15.$	$0.010$	$77.14,19.94$	blue
2	$\dots \dots \dots \dots \dots \dots \dots$	$0.015$	$94.99,23.33$	black
3	$\dots \dots \dots \dots \dots \dots \dots$	$0.045$	$363.40,18.45$	yellow
4	$\dots \dots \dots \dots \dots \dots \dots$	$0.060$	$383.05,13.98$	green
5	$\dots \dots \dots \dots \dots \dots \dots$	$0.075$	$394.02,11.24$	red
6	$\dots \dots \dots \dots \dots \dots \dots$	$0.500$	$427.82,1.71$	cyan
7	$\dots \dots \dots \dots \dots \dots \dots$	$0.800$	$429.90,1.07$	magenta

Equations

$$\frac{dx}{dt}$$

$$\frac{dy}{dt}$$

$$\frac{dx}{dt}$$

$$\frac{dy}{dt}$$

$$\frac{d\xi }{dt}+\zeta \xi$$

$$\Phi (x,y)$$

$$\Phi (x_1,0)$$

$$r(1-\frac{x^{\ast }}{k})-\frac{p(1-m{y}^{\ast })y^{\ast }}{1+a{x}^{\ast }(1-m{y}^{\ast })}-q_1{E}_1=0,\frac{ep(1-m{y}^{\ast })}{1+a{x}^{\ast }(1-m{y}^{\ast })x^{\ast }}-d-q_2{E}_2=0.$$

$$J=\left[\begin{array}[]{cc}r-2\,\frac{rx}{k}-q_{1}E_{1}-\frac{py(1-my)}{(1+ax(1-my))^{2}}&-\frac{px(ax(1-my)^{2}+1-2my}{(1+ax(1-my))^{2}})\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\frac{epy(1-my)}{(1+ax(1-my))^{2}}&-d-q_{2}E_{2}+\frac{epx(ax(1-my)^{2}+1-2my}{(1+ax(1-my))^{2}})\end{array}\right].$$

$$J_{0}=\left[\begin{array}[]{cc}r-\,q_{1}E_{1}&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&-d-q_{2}E_{2}\end{array}\right].$$

$$J_{1}=\left[\begin{array}[]{cc}r-2\,\frac{rx_{1}}{k}-q_{1}E_{1}&-\frac{px_{1}}{(1+ax_{1})}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&-d-q_{2}E_{2}+\frac{epx_{1}}{(1+ax_{1})}\end{array}\right].$$

$$J_{*}=\left[\begin{array}[]{cc}r-2\,\frac{rx^{*}}{k}-q_{1}E_{1}-\frac{py^{*}(1-my^{*})}{(1+ax^{*}(1-my^{*}))^{2}}&-\frac{px^{*}(ax^{*}(1-my^{*})^{2}+1-2my^{*}}{(1+ax^{*}(1-my^{*}))^{2}})\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\frac{epy^{*}(1-my^{*})}{(1+ax^{*}(1-my^{*}))^{2}}&-d-q_{2}E_{2}+\frac{epx^{*}(ax^{*}(1-my^{*})^{2}+1-2my^{*}}{(1+ax^{*}(1-my^{*}))^{2}})\end{array}\right].$$

$$ea{x}^{\ast }(1-m{y}^{\ast }{)}^2(p+y^{\ast }-p{x}^{\ast })+e{y}^{\ast }$$

$$p{y}^{\ast }(1-m{y}^{\ast })$$

$$\frac{dV}{dt}$$

$$=$$

$$4rme(1+ak)(1+a{x}^{\ast })>p{m}^2{y}^{\ast }{}^2.$$

$$r(1-\frac{x}{k})-\frac{p(1-my)y}{1+ax(1-my)}-q_1{E}_1=0,$$

$$\frac{ep(1-my)x}{1+ax(1-my)}-d-q_2{E}_2=0,$$

$$(p_1{q}_{1}x-c_1)E_1+(p_2{q}_{2}y-c_2)E_2=0.$$

$$J^1={\int }_0^{\infty }{e}^{-\delta t}[(p_1{q}_{1}x-c_1)E_1(t)+(p_2{q}_{2}y-c_2)E_2(t)]dt,$$

$$H$$

$$\frac{\partial H}{\partial E_1}=0,⟹{\lambda }_1(t)=e^{-\delta t}(p_1-\frac{c_1}{q_{1}x});\frac{\partial H}{\partial E_2}=0,⟹{\lambda }_2(t)=e^{-\delta t}(p_2-\frac{c_2}{q_{2}y}).$$

$$\frac{d{\lambda }_1}{dt}$$

$$\displaystyle\delta\bigl{(}p_{1}-\frac{c_{1}}{q_{1}x}\bigr{)}$$

$$\frac{d{\lambda }_2}{dt}$$

$$\delta (p_2-\frac{c_2}{q_{2}y})$$

$$E_1^{\ast }$$

$$E_2^{\ast }$$

$$\dot{x_1}$$

$$\dot{y_1}$$

$$\sigma$$

Full text

Mathematics Subject Classification: 92D25, 92D30, 92D40.

Keywords: Harvesting, Refuge, Bionomic equilibrium, Limit cycle, Hopf bifurcation, Numerical simulations

1 Introduction

The growing needs for more food and more energy have led to increase exploitation of natural resources. The problem related to multi-species fisheries have been drawing attention of researcher in the recent years (cf. Clark [2, 3, 4, 5], Mesterton-Gibbons [6]). The economic progress and ecological balance always have a conflicting interests. Therefore, concerning the conservation for the long term benefits of humanity, there is a wide range of interest in the use of bionomic modeling to gain greater insight in the scientific management of renewable resources like fisheries and forestry. An excellent introduction of optimal management of renewable resources has been presented by (cf. Clark [2]). The techniques and issues associated with the dynamic economic models of natural resources exploitation are further developed by Clark (cf. Clark [3, 4, 7]). Harvesting is one of the most important issue to control extinction and minimization of exploitation of renewable resources. Renewable harvesting policy is indisputably one of the major and interesting problems from ecological and economical point of view. There is also a realistic phenomena for surviving himself is refuge. The use of spatial refuges by the prey is one of the more relevant behavioral traits that affect the dynamics of predator prey systems. By using refuges some fraction of prey population are partially protected against predators. The existence of refuges has significant influence on the coexistence of the predator-prey systems. It has a strong impact for stabilizing and destabilizing the dynamical nature of the system. The study of consequence of prey refuge on the dynamic of prey-predator interaction can be recognized as a major effect in applied mathematics and theoretical ecology (cf. Hassell and May [8], Hassell [9], Holling [10], Hoy [11], Smith [12]). Some of the empirical and theoretical works have investigated the effect of refuge and gave a decision that the refuge used by prey have a stabilizing effect on prey predator interaction and prey species can be prevented from extinction by using this policy (cf. Collings [13], Freedman [14], Gonzalez-Olivares and Ramos-Jiliberto [15], Hochberg and Holt [16], Huang et al. [17], Krivan [18], McNair [19], Ruxton [20], Sih [21], Taylor [22] ).

Also, Brauer and Soudak (cf. Brauer and Soudak [23, 24, 25, 26]) studied a class of predator prey models under constant rate of harvesting and under constant percentage of harvesting of both species simultaneously. They have shown that how to classify the possibilities of the quantitative behavior of solution to locate the set of initial values in which the trajectories of the solution approaches to either an asymptotic stable equilibrium or an asymptotic stable limit cycle. Dai and Tang [27] studied the following predator-prey model in which two ecological interacting species are harvested independently with constant rates of harvesting:

[TABLE]

where $x(t)$, $y(t)$ represents prey and predator species respectively. $r$, $k$, $d$, $a$, $c$, $\phi(x)$ are the intrinsic growth rate, carrying capacity of prey, death rate of predator, maximum per capita consumption rate of predator, conversion rate and general predator response function on prey respectively; $\mu$, $h$ are constant harvesting rates. The most crucial element in an interacting population model is the “functional response” the expression that describes the rate at which the number of prey is consumed by a predator. They have shown the complicated dynamics of the system. Kar (cf. Kar [28]) studied the prey-predator model with prey refuge and individuals are subject to proportional harvesting rates. As far as our knowledge goes there is no work in predator-prey model system with independent harvesting and prey refuge proportional to both the species have done even consideration of prey refuge proportional to both species is one step closure to reality. Keeping this in mind an attempt is made in the present investigation to study the effect of refuge as well harvesting on a Holling Type-II prey-predator model. In this paper we take the following model:

[TABLE]

where $x(t)$, $y(t)$ represents the prey and predator populations respectively at any time t. $r$, $k$, $p$, $m$, $d$, $q_{1}$, $q_{2}$, are all positive constants and have their biological meanings, accordingly $r$, represent the intrinsic growth rate of prey species, $k$ is the carrying capacity of the prey in absence of predator and harvesting, $p>0$ is the maximum per capita consumption rate of predator, $e$ $(0<e<1)$ is the efficiency by which predator converts consumed prey into new predator, $d>0$ is the death rate of predator. $E_{1}\geq 0$, $E_{2}\geq 0$, denotes the harvesting effort for the prey and predator respectively. $q_{1}E_{1}x$ and $q_{2}E_{2}y,$ represent the catch of the respective species, where $q_{1}$, $q_{2},$ are catchability coefficients of the prey and predator species respectively. The present model incorporates a refuge proportional to both the prey-predator determines. i.e., $mxy$ from the predator species, where $m\in[0,\,1]$ is a constant. Incorporation of prey refuge leaves $(1-my)x$ of the prey available to be hunted by the predator.

2 Some preliminaries

2.1 Existence and positive invariance

Letting $X=(x,y)^{t}$, $F:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$, $F=(F_{1},F_{2})^{t}$, the system (1.2) can be rewritten as $\frac{d{X}}{dt}=F({X})$. Here $F_{i}\in C^{\infty}(\mathbb{R})$,  $i=1,2.$ Since the vector function $F$ is a smooth function of the variables $(x,y)$ in the positive quadrant $\Omega=\{(x,y);x>0,\,y>0\}\subset\mathbb{R}^{2}_{+}$, local existence and uniqueness of the solution set hold.

2.2 Boundedness

Proposition 1

All the solutions of the model (1.2) are bounded uniformly.

Proof. Let us consider a function $\xi=x+\frac{y}{e}.$ For any $\zeta=d+q_{2}E_{2}>0$,

[TABLE]

By applying the theory of differential inequality (cf. Brickhoff and Rota [29]), we have the following inequality $0<\xi(x,y)<\frac{\kappa}{\zeta}(1-e^{-\zeta})+\xi(0)e^{-\zeta t}<\max{\{\frac{\kappa}{\zeta},\xi(0)\}}$. Therefore, $\underset{t\rightarrow+\infty}{\limsup}\,\xi(t)\leq\frac{\kappa}{\zeta}$ with last bound independent of initial conditions. Hence, all the solutions of the system (1.2) starting from $\mathbb{R}^{2}_{+}$ evolve with respect to time and remain in the compact region $\mathbb{R}_{xy}={\{{(x,y)\in\mathbb{R}^{2}_{+}}:\xi(x,y)\leq\frac{\kappa}{\zeta}+\epsilon\}}$ for any $\epsilon>0.$

2.3 Persistence

Persistence of a predator-prey model system plays an important role in mathematical ecology since the criteria of persistence for ecological systems is a condition that ensuring the long-term survival of all the species. Here we have shown the persistence by using average Lyapunov function (cf. Gard and Hallam [30]), together with its boundedness.

Considering the average Lyapunov function $P(x,\,y)=x^{\rho_{1}}y^{\rho_{2}},$ where $\rho_{1}$ and $\rho_{2}$ are undetermined positive constant. Let us define the function $\phi$ as follows:

[TABLE]

Now, $\Phi(0,0)=\rho_{1}r-\rho_{1}q_{1}E_{1}-\rho_{2}d-\rho_{2}q_{2}E_{2}>0$, if $\rho_{1}r>\rho_{1}q_{1}E_{1}+\rho_{2}d+\rho_{2}q_{2}E_{2}.$

[TABLE]

Hence, the solution of the system (1.2) is permanent, if the conditions $\rho_{1}r>\rho_{1}q_{1}E_{1}+\rho_{2}d+\rho_{2}q_{2}E_{2}$ and p>\frac{d+q_{2}E_{2}}{e}\bigl{(}a+\frac{1}{k(1-\frac{q_{1}E_{1}}{r})}\bigr{)} are satisfied.

3 Analysis of equilibria

The equilibria of the system $(\ref{eqn1})$ are (i) the trivial equilibrium $E^{0}(0,\,0)$;  (ii) the axial equilibrium E^{1}(x_{1},\,0)=\bigl{(}k(1-\frac{q_{1}E_{1}}{r}),\,0\bigr{)},   and the interior equilibrium point (iii) $E^{*}(x^{*},\,y^{*})$.

We are interested on the interior equilibrium point $E^{*}$, where $x^{*}$ and $y^{*}$ are satisfying the following system of equations:

[TABLE]

One can easily found that the interior equilibrium $E^{*}$ is feasible if the conditions $E_{1}<\frac{r}{q_{1}}(1-\frac{x^{*}}{k})$ and $E_{2}<\frac{1}{aq_{2}}(ep-ad)$ are satisfied.

4 Local stability

We now investigate the local asymptotically stability (LAS) of $(\ref{eqn1})$ around the feasible equilibrium points. The Jacobian matrix at an arbitrary point $(x,\,y)$ is

[TABLE]

Thus, the Jacobian matrix of the system $(\ref{eqn1})$ around the trivial equilibrium point $E^{0}(0,\,0)$ is

[TABLE]

The eigenvalues for the steady state $E^{0}(0,\,0)$ are $(r-\,q_{1}E_{1})$ and $-(d+q_{2}E_{2}).$ The eigenvalue $r-\,q_{1}E_{1}$ is positive or negative according to $E_{1}<\frac{r}{q_{1}}$ or $E_{1}>\frac{r}{q_{1}},$ i.e., the equilibrium point $E^{0}(0,\,0)$ is locally asymptotically stable (cf. Fig 1(a)) if the BTP (Biotechnical Productivity ) of the prey species $(\frac{r}{q_{1}})$ is less than the effort $(E_{1})$ for the prey species. If $BTP>E_{1}$, then equilibrium point $E^{0}(0,0)$ will be a saddle for the system (1.2).

The Jacobian matrix at the point $E^{1}(x_{1},\,0)$ of the system $(\ref{eqn1})$ is

[TABLE]

The eigenvalues of $J_{1}$ are $\lambda_{1}=r-2\,\frac{rx_{1}}{k}-q_{1}E_{1}$ and $\lambda_{2}=-d-q_{2}E_{2}+\frac{epx_{1}}{(1+ax_{1})}.$ The eigenvalues are negative if the condition \bigl{(}1-\frac{d+q_{2}E_{2}}{epk-ak(d+q_{2}E_{2})}\bigr{)}<\frac{q_{1}E_{1}}{r}<1 is satisfied and hence $E^{1}(x_{1},\,0)$ is locally asymptotically stable (cf. Fig 1(b)).

The Jacobian matrix around the interior equilibrium point $E^{*}(x^{*},\,y^{*})$ of the system (1.2) is

[TABLE]

Both the eigenvalues of $J_{*}$ will be negative if the following conditions hold:

[TABLE]

Hence, under this parametric conditions (4.1)-(4.2) the interior equilibrium point $E^{*}(x^{*},\,y^{*})$ will be locally asymptotically stable.

5 Global stability around the interior equilibrium $E^{*}$

Theorem 5.1

The interior equilibrium $E^{*}$ is globally asymptotically stable if the condition $4rme(1+ak)(1+ax^{*})>pm^{2}y^{*}{{}^{2}}$ is satisfied.

Proof. Let us consider the suitable Lyapunov function as followers:

$V=(x-x^{*}-x^{*}ln\frac{x}{x^{*}})+(y-y^{*}-y^{*}ln\frac{y}{y^{*}})$. Clearly V is positive definite for all $(x,\,y)\in\mathbf{R^{2}_{+}}\setminus(x^{*},y^{*})$. Taking the time derivative along solution of the system (1.2), we have

[TABLE]

where $\alpha=\frac{r}{k}+\frac{pm(y+y^{*})}{A}-\frac{py(1+myy^{*})}{A},$ $\beta=\frac{mxep}{A},$ $\gamma=\frac{px(1+myy^{*})+p(1-m(y+y^{*}))-ep-epy^{*}-p(y+y^{*})}{A},$ $A=(1+ax(1-my))(1+ax^{*}(1-my^{*})).$

Now, $4\alpha\beta-\gamma^{2}>0,$ if

[TABLE]

Thus, the quadratic form (5.2) is positive definite hence $\frac{dV}{dt}<0$ along all the trajectories in the first quadrant except $(x^{*},\,y^{*}).$ Also $\frac{dV}{dt}|_{E^{*}}=0.$ The proof follows from the suitable Lyapunov function $V$ and Lyapunov-LaSalle’s invariant principle (cf. Hale [31]). Hence the equilibrium point is globally asymptotically stable if the condition (5.3) is satisfied.

6 Bionomic equilibrium

Let $c_{1}$ and $c_{2}$ are the fishing cost of prey and predator species per unit effort; $p_{1}$ and $p_{2}$ are the prices of prey and predator per unit biomass.

Therefore, the economic rent (net revenue) at any time can be taken as $\pi=(p_{1}q_{1}x-c_{1})E_{1}+(p_{2}q_{2}x-c_{2})E_{2}$ $=\pi_{x}+\pi_{y}$ , where $\pi_{x}=(p_{1}q_{1}x-c_{1})E_{1}$ and $\pi_{y}=(p_{2}q_{2}y-c_{2})E_{2}$ are the net revenue for the prey and predator species respectively. The bionomic equilibrium is defined as a point where the biological and the economical equilibrium takes places.

The bionomic equilibrium $(x_{\infty},\,y_{\infty},\,E_{1}^{\infty},\,E_{2}^{\infty})$ is given by the solutions of following simultaneous equations:

[TABLE]

The bionomic equilibria are determined in different case as follows:

$\bf CaseI$: If $c_{2}>p_{2}q_{2}y$, i.e., the fishing cost is greater than the revenue for the predator species, then the predator fishing will be stopped (i.e., $E_{2}$=0). Only the prey fishing will be in operational (i.e., $E_{1}<p_{1}q_{1}x$ ).

We then have $x_{\infty}=\frac{c_{1}}{p_{1}q_{2}}$. Therefore, putting this value in $(\ref{eqn3})$ and $(\ref{eqn4})$ we have the bionomic equilibrium point $(y_{\infty},E_{1}^{\infty})$ in the y-$E_{1}$ plane if the parametric condition $E_{1}^{\infty}<\frac{r}{q_{1}}(1-\frac{c_{1}}{p_{1}q_{1}k})$ holds.

$\bf CaseII$: If $c_{1}>p_{1}q_{1}x,$ i.e., the fishing cost is greater than the revenue in the prey fishing then the prey fishery will be closed (i.e., $E_{1}=0$) . Only the predator fishing will be in operational (i.e., $c_{2}<p_{2}q_{2}y$). In this case the bionomic equilibrium $(x_{\infty},E_{2}^{\infty})$ will be in the first quadrant of x-$E_{2}$ plane if $E_{2}^{\infty}\geq\frac{1}{q_{2}}(d-\frac{p_{2}q_{2}erk}{4c_{2}}).$

$\bf CaseIII$: If $c_{1}>p_{1}q_{1}x$ and $c_{2}>p_{2}q_{2}y$, then the fishing cost is greater than the revenue for both the species and the whole fishery will be closed.

$\bf CaseIV$: If $c_{1}<p_{1}q_{1}x$ and $c_{2}<p_{2}q_{2}y$ holds simultaneously, then the revenue for both the species being positive then the whole fishery will be in operational, in this case $x_{\infty}=\frac{c_{1}}{p_{1}q_{1}}$, $y_{\infty}=\frac{c_{2}}{p_{2}q_{2}}$ , therefore using this values in $(\ref{eqn3})$ and $(\ref{eqn4})$ we have the non trivial bionomic equilibrium point $(x_{\infty},\,y_{\infty},\,E_{1}^{\infty},\,E_{2}^{\infty})$ exists if the parametric conditions (i) \frac{r}{q_{1}}(1-\frac{c_{1}}{p_{1}q_{1}k})>\frac{p(p_{2}q_{2}-mc_{2})c_{2}}{p_{2}q_{2}\bigl{(}p_{1}q_{1}p_{2}q_{2}+ac_{1}(p_{2}q_{2}-mc_{2})\bigr{)}}  and (ii) $d<\frac{ep(p_{2}q_{2}-mc_{2})c_{1}}{p_{1}q_{1}p_{2}q_{2}+ac_{1}(p_{2}q_{2}-mc_{2})}$ are satisfied.

7 Optimal harvesting policy

The fundamental problem in commercial exploitation of renewable resources is to determine the optimal trade-off between current and future harvests. As observed by Clark (cf. Clark [2]). This problem which is very assured of resources conservation is an exceedingly different one not from the Mathematical view point perhaps, but certainly from a political and philosophical view point. However, we look at the problem from the economic view point only and we have to use the time discounting policy to handle the question of inter temporal benefits. This discounting is a normal practice in business management (cf. Solow [32]). For determination of an optimal harvesting policy, we consider the present value $J^{1}$ of continuous stream of revenue as follows:

[TABLE]

where $\delta$ denote the instantaneous annual rate of discount. We are to optimize the equation $(\ref{eqn6})$ subject to the state equation $(\ref{eqn1})$ by using Pontryagian’s maximum principle (cf. Pontryagin [33]).

Let us consider the Hamiltonian function $H$ as follows:

[TABLE]

where $\lambda_{1}(t)$ and $\lambda_{2}(t)$ are adjoint variables. $E_{1}(t),E_{2}(t)$ are the control variables subject to the constraint $0\leq E_{i}(t)\leq(E_{i})_{max}$, $(i=1,2)$. The control variable $E_{1}(t),E_{2}(t)$ appear linearly in the Hamiltonian function $H$. Assuming that the control constraint are not binding, i.e., optimal solution does not occur at $(E_{i})_{max}$ or $(E_{i})_{min}$, we have singular control (cf. Clark [2]) given by $\frac{\partial H}{\partial E_{i}}=0.$

Therefore,

[TABLE]

Thus, the shadow prices $e^{\delta t}\lambda_{i}(t)$, $(i=1,2)$ do not vary with time in optimal equilibrium. Hence they satisfy the transversally condition at $t\rightarrow+\infty$, i.e., they remain bounded as $t\rightarrow\infty.$ Again, $\frac{\partial H}{\partial E_{1}}=0$ $\implies\lambda_{1}q_{1}xe^{\delta t}=\frac{\partial\pi_{x}}{\partial E_{1}}$;  $\frac{\partial H}{\partial E_{2}}=0$ $\implies\lambda_{2}q_{2}ye^{\delta t}=\frac{\partial\pi_{y}}{\partial E_{2}}$.

From these relations one can say that for each species the user cost of harvesting per unit effort must be equal to the discounted value of the future marginal profit of effort at the steady state level.

Now we are to find out the optimal solution of the problem as follows:

[TABLE]

Substituting the values of $\lambda_{1}$ and $\lambda_{2}$ from $(\ref{eqn7})$, we have the relation between $x$ and $y$ as follows:

[TABLE]

Again

[TABLE]

Substituting the values of $\lambda_{1}$ and $\lambda_{2}$ in $(\ref{eqn10})$ we have the following relation in $x$ and $y$

[TABLE]

Solving equation $(\ref{eqn9})$ and $(\ref{eqn11})$ for $x$ and $y$ we have the optimal equilibrium $(x^{*},y^{*})$ and the optimal harvesting efforts $E_{1}^{*}$ and $E_{2}^{*}$ can be determined by the following equations:

[TABLE]

8 Bifurcation analysis

The subject of bifurcation is the study of structurally unstable systems. This topic is a branch of mathematics concerned with dynamical systems which suffer sudden qualitative changes in parameters. A small change in parameter causes a topological change. The important question in the field of local bifurcation theory is that a system depends on a control parameters, as parameter changes, what happens to the non hyperbolic equilibria. Structurally unstable dynamic systems can be classified according to the number of parameters that appears in the differential equations describing the dynamics of the system. For good introduction on the basics of local bifurcation analysis the interested readers are referred to check out the books by Guckenheimer [34], Wiggins [35], and Kuznetsov [36]. In our study of bifurcation the parameter $r$ has been chosen for the transcritical bifurcation and $m$ has been taken as bifurcation parameter for the analysis of Hopf- bifurcation for the present system (1.2).

8.1 Existence of transcritical bifurcation around $E_{1}$

The the system (1.2) experiences a transcritical bifurcation around the axial equilibrium $E_{1}$ as the parameter $r$ crosses its critical value $r=q_{1}E_{1}=r_{tc}$. One of the eigenvalues of $J_{1}$ will be zero iff $\text{Det}(J_{1})=0$, which implies either $r-E_{1}q_{1}=0$, or $-d-q_{2}E_{2}+\frac{epx_{1}}{1+ax_{1}}=0$. Let $v$, $w$ are the eigenvector corresponding to the eigenvalue $\lambda=0$ for the matrices $J_{1}$ and $J_{1}^{T}$. The eigenvectors $v$, $w$ are found as $v=(v_{1},v_{2})^{T}=(1,0)$ and w=(w_{1},w_{2})^{T}=\bigl{(}1,\,\frac{px_{1}}{epx_{1}-(1+ax_{1})(d+E_{2}q_{2})}\bigr{)}^{T}. With these eigenvectors, it is found that (i) $w^{T}F_{r}(E_{1},\,r_{tc})=0$, (ii) $w^{T}[DF_{r}(E_{1},\,r_{tc})v]=(1-\frac{2x_{1}}{k})\neq 0,$ (iii) $w^{T}[D^{2}F(E_{1},\,r_{tc})(v,v)]=-\frac{2r_{tc}}{k}\neq 0.$ Hence, due to Satomayor theorem (cf. Sotomayar [37]), the system experience transcritical bifurcation around the axial equilibrium $E^{1}$ at $r=r_{tc}$.

8.2 Hopf bifurcation

It can be easily conclude that the equilibrium point $E^{*}$ may loss its stability through Hopf bifurcation under certain parametric condition. Considering $m$ as a bifurcation parameter one can detect the threshold value of $m=m_{h}$, which satisfy $\text{Tr}(J_{E^{*}})|_{m=m_{h}}=0$ and $\text{Det}(J_{E^{*}})|_{m=m_{h}}>0$.

The transversality condition for the Hopf bifurcation (cf. Carr [38], Hassard and Kazarinoff [39], Perko [40]) is $\frac{d}{dm}(\text{Tr}(J_{E^{*}}))|_{m=m_{h}}=\frac{py_{m}(2aex_{m}^{2}-amx_{m}y_{m}^{2}+ax_{m}y_{m}+2ex_{m}-y_{m})}{(1+ax_{m}(-my_{m}+1))^{3}}|_{m=m_{h}}\neq 0$, where $x_{m}$ and $y_{m}$ indicate their functionality with respect to the parameter $m$.

The interior equilibrium point $E^{*}$ loss its stability through the non-degenerate Hopf-bifurcation when the parametric restriction $\text{Tr}(J_{E^{*}})|_{m=m_{h}}=0$ and the transversality conditions are satisfied simultaneously.

Now we calculate the Lyapunov number to determine the nature of Hopf-bifurcating periodic solutions. Introducing perturbations $x=x_{1}+x_{m}|_{m=m_{h}},$ $y=y_{1}+y_{m}|_{m=m_{h}}$ in (1.2) and then expanding in Taylor series, we have

[TABLE]

where $a_{10}$, $a_{01}$, $b_{10},$ $b_{01}$ are the elements of the Jacobian matrix evaluated at the equilibrium point $E^{*}$ with $m=m_{h}$, hence $a_{10}+b_{01}=0$ and $\Delta=a_{10}b_{01}-a_{01}b_{10}>0.$ The expression of the coefficients $a_{ij}$ and $b_{ij}$ are given bellow:

$a_{10}=\frac{\partial F_{1}}{\partial x}|_{(E^{*},\,m_{h})}$,  $a_{01}=\frac{\partial F_{1}}{\partial y}|_{(E^{*},\,m_{h})}$, $a_{12}=\frac{1}{2}\frac{\partial^{3}F_{1}}{\partial x\partial y^{2}}|_{(E^{*},\,m_{h})}$,  $a_{21}=\frac{1}{2}\frac{\partial^{3}F_{1}}{\partial x^{2}\partial y}|_{(E^{*},\,m_{h})}$,

$a_{20}=\frac{1}{2}\frac{\partial^{2}F_{1}}{\partial x^{2}}|_{(E^{*},\,m_{h})}$,  $a_{11}=\frac{\partial^{2}F_{1}}{\partial x\partial y}|_{(E^{*},\,m_{h})}$, $a_{30}=\frac{1}{6}\frac{\partial^{3}F_{1}}{\partial x^{3}}|_{(E^{*},\,m_{h})}$;

$b_{10}=\frac{\partial F_{2}}{\partial x}|_{(E^{*},\,m_{h})}$,  $b_{01}=\frac{\partial F_{2}}{\partial y}|_{(E^{*},\,m_{h})}$, $b_{12}=\frac{1}{2}\frac{\partial^{3}F_{2}}{\partial x\partial y^{2}}|_{(E^{*},\,m_{h})}$,  $b_{21}=\frac{1}{2}\frac{\partial^{3}F_{3}}{\partial x^{2}\partial y}|_{(E^{*},\,m_{h})}$,

$b_{20}=\frac{1}{2}\frac{\partial^{2}F_{2}}{\partial x^{2}}|_{(E^{*},\,m_{h})}$,  $b_{11}=\frac{\partial^{2}F_{2}}{\partial x\partial y}|_{(E^{*},\,m_{h})}$, $b_{30}=\frac{1}{6}\frac{\partial^{3}F_{2}}{\partial x^{3}}|_{(E^{*},\,m_{h})}$.

The value of first Lyapunov number (cf. Perko [40]), which helps to determine the nature of the stability of limit cycle arising through Hopf-bifurcation is given by

[TABLE]

where the values of $a_{ij}$ and $b_{ij}$, $i,\,j=0,1,2,3$ are included in Appendix A.

If the $\sigma<0$, the equilibrium point $E^{*}$ destabilized through a supercritical Hopf-bifurcation, and if $\sigma>0$ then the Hopf bifurcation is subcritical.

9 Numerical simulations

In this section we perform numerical simulation to validate our analytical findings of the previous Section by making uses of the computing Software Packages MATLAB-R2015a and Maple-18. The analytical findings of the present study are summarized and presented Schematically in Tables 1, Table 2 and Table 3. It is very difficult to validate the model results with realistic data so far proportional refuse and harvesting are considered in natural field. These results are verified by means of numerical illustrations of which some chosen ones are shown in figures. We took a hypothetical set of parameter values to illustrate our results. In Section $(7)$ optimal harvesting policy and corresponding effort are determined. For numerical justification we have taken a set of parameter values: $r=3$, $a=0.008,$ $d=0.04,$ $m=0.02,$ $p=0.2,$ $q_{1}=0.2,$ $q_{2}=0.6$, $k=500$, $e=0.15$, $p_{1}=2$, $p_{2}=3$, $c_{1}=1$, $c_{2}=2$, $\delta=.004$ for this set of parameter values the optimal equilibrium is found at $(188.5858,\,30.6567)$ and the corresponding harvesting efforts are $E_{1}^{*}=1.8534$ and $E_{2}^{*}=5.8875.$ In Section $(8)$, it is found that the present system (1.2) experiences Hopf-bifurcation for the bifurcation parameter $m$ and we find out the threshold value of $m$ and Lyapunov number $\sigma$ to determine the nature of Hopf-bifurcation (cf. Sen et al. [41]). For numerical validation we take the fixed set of parameter values: $r=3,$ $a=0.008,$ $d=0.04,$ $p=0.2$, $q_{1}=0.2$, $q_{2}=0.6,$ $E_{1}=2,$ $E_{1}=2,$ $k=500,$ $e=0.15$. It is found that the threshold value of the bifurcation parameter $m=m_{h}=0.010695$ at which the system (1.2) experiences Hopf bifurcation around the interior equilibrium $E^{*}$. For this choice of parameter values the first Lypunov number is $\sigma=-000143<0$. Hence the Hopf-bifurcation is supercritical. It is observed that the trivial and axial equilibria are locally asymptotically stable (LAS) by starting the solution plots from different initial conditions in the neighborhood of $E^{0}$ and $E^{1}$, all the solution plots eventually converge to the equilibria respectively (cf Fig. 1: (a)-(b)). Fig. 2: (a)-(b) shows the limit cycle behavior of the system (1.2) around the equilibrium point $E^{*}=(67.86,\,18.00)$ for the set of parameter values: $r=3,$ $a=0.008,$ $d=0.04,$ $m=0.005$, $p=0.2$ $q_{1}=0.2$, $q_{2}=0.6,$ $E_{1}=2,$ $E_{1}=2,$ $k=500,$ $e=0.15$. In Fig. 3 shows that the system experiences Hopf-bifurcation for bifurcation parameter $m=0.01$ and the other parameters are same as Fig. 2. Fig. 4: (a)-(b) shows that the system (1.2) is globally asymptotically stable and conversing to the point $E^{*}(94.99,\,23.33)$ for $m=0.015$ and the other parameters are same as Fig. 2. Fig. 5 Shows that gradual increase of the coefficient of refuge gives more protection to the prey species and reduce the rate of predation of predator, as a result the volume of prey species became larger and larger and the predator became smaller and smaller until it goes to extinct. This fact is presented in tabular form (cf. Table 3).

10 Conclusion

In this paper we consider a prey-predator harvesting model with Holling type -II response function incorporating prey refuge proportional to both the species. The novelty of our work lies in taking such kind of refuge function, which is more realistic phenomenon in ecosystem. We emphasize, on refuge coefficient that how the refuge function changes the system dynamics. We also observed from the both mathematical and empirical point of view that traits of behavioral policy of prey refuge has a stabilizing effect on a predator prey dynamics and this policy can help the prey species from extinction (cf. Anderson [42], Cressman and Garay [43], Magalhaes et al. [44], Rudolf et al. [45], Sarwardi et at. [46], Sarwardi et al. [47], Mukherjee [48], Ma et al. [49]). The key question is now the refuge alters the evolutionary dynamics of the system. In this paper we have shown that the dynamics of our system depends upon the functionality of refuge construction. In the present paper, it is assumed that functionality of prey refuge depends not only on prey size but on both the species. From the analytical view point the study of influence of prey refuge on the dynamics of interacting population is at present recognized as a significant and challenging issues (cf. Collings [13], Huang [50]). We have shown that depending upon the bifurcation parameter $m$ the system exhibits stability as well as bifurcation around some of the equilibrium points. Also we have found that the optimal harvesting policy and the corresponding optimal effort using $E_{1}$ and $E_{2}$ as control parameter.

Before ending our conclusion we must say that, there are still some options to improve our model system to have much richer dynamics that what we have found in the present study. Here are given some rooms for our future studies. Firstly, it would be more logical and proven reality that the term representing time delay used in digestion or gestation period for the predator species to produce new born have not been taken into account. Secondly, the harvesting efforts $E_{1}\text{and}E_{2}$ can be taken as time dependent functions. Thirdly, the harvesting effort can be taken as non-linear functions of both species. Fourthly, on the basis of the fact that more species could give more stable ecosystem, one can consider one more prey or predator or a pair of predator and prey into the exited system to have more stable system from the biological point of view. The incorporation of all such relaxations existed in the present system in the future updated model would certainly be of some help to empirical researcher to predict their findings one step closer to the real situation from the ecological point of view. In our next paper we will study the dynamics of the system taking non-linear harvesting effort incorporating the taxation to sustain the renewable resources.

Acknowledgements: Authors are thankful to the Department of Mathematics, Aliah University for providing opportunities to perform the present work. The corresponding author Dr. S. Sarwardi is thankful to his Ph.D. supervisor Prof. Prashanta Kumar Mandal, Department of Mathematics, Visva-Bharati (a Central University) for his generous help and continuous encouragement while preparing this manuscript.

11 Appendix

The expressions $a_{ij}$ and $b_{ij}$, appearing in (8.3) are defined by $a_{ij}=\frac{1}{i!j!}\frac{\partial^{i+j}(F_{1})}{\partial x^{i}y^{j}}|_{(E^{*},\,m_{h})}$,   $b_{ij}=\frac{1}{i!j!}\frac{\partial^{i+j}(F_{2})}{\partial x^{i}y^{j}}|_{(E^{*},\,m_{h})}$ and their explicit values are as follows:

$a_{10}=r-q_{1}E_{1}-\frac{2x}{k}-\frac{p(1-my)y}{(1+ax(1-my))^{2}}|_{E^{*}},$ a_{01}=-\frac{px\bigl{(}ax(1-my)^{2}+1-2my\bigr{)}}{(1+ax(-my+1))^{2}}|_{E^{*}}, $b_{10}=\frac{ep(1-my)y}{(1+ax(1-my))^{2}}|_{E^{*}}$,

b_{01}=-d-q_{2}E_{2}+\frac{epx\bigl{(}ax(1-my)^{2}+1-2my\bigr{)}}{(1+ax(1-my))^{2}}|_{E^{*}}, $a_{11}=\frac{-p(amxy-ax+2my-1)}{(1+ax(-my+1))^{3}}|_{E^{*}},$ $a_{20}=-\frac{r}{k}+\frac{p(my-1)^{2}ya}{(1+ax(1-my))^{3}}|_{E^{*}}$, $a_{02}=\frac{-2pmx(ax+1)}{(1+ax(1-my))^{3}}|_{E^{*}},$ $a_{30}=\frac{p(my-1)^{3}ya^{2}}{(1+ax(1-my))^{4}}|_{E^{*}},$ $b_{20}=\frac{ep(my-1)^{2}ya}{(1+ax(1-my))^{3}}|_{E^{*}},$

$a_{21}=\frac{pa(my-1)(amxy-ax+3my-1)}{(1+ax(1-my))^{4}}|_{E^{*}},$ a_{12}=\frac{pm(a^{2}mx^{2}y-a^{2}x^{2}+2amxy+1)}{(1+ax(1-my))^{4}}|_{E^{*}},$$b_{30}=\frac{-ep(my-1)^{3}ya^{2}}{(1+ax(1-my))^{4}}|_{E^{*}}, $b_{11}=\frac{ep(amxy-ax+2my-1)}{(1+ax(1-my))^{3}}|_{E^{*}},$ $b_{21}=\frac{-ep(my-1)a(amxy-ax+3my-1)}{(1+ax(1-my))^{4}}|_{E^{*}},$

$b_{12}=\frac{-epm(a^{2}mx^{2}y-a^{2}x^{2}+2amxy+1)}{(1+ax(1-my))^{4}}|_{E^{*}}$.