Prey-predator model in drainage system with migration and harvesting


 In this paper, we consider a prey-predator model with a reserve region of predator where generalist predator cannot enter. Based on the intake capacity of food and other factors, we introduce the predator population which consumes the prey population with Holling type-II functional response; and generalist predator population consumes the predator population with Beddington-DeAngelis functional response. The density-dependent mortality rate for prey and generalist predator are considered. The equilibria of proposed system are determined. Local stability for the system are discussed. The environmental carrying capacity is considered as a bifurcation parameter to evaluate Hopf bifurcation in the neighbourhood at an interior equilibrium point. Here the fishing effort is used as a control parameter to harvest the generalist predator population of the system. With the help of this control parameter, a dynamic framework is developed to investigate the optimal utilization of resources, sustainability properties of the stock and the resource rent. Finally, we present a numerical simulation to verify the analytical results, and the system is analyzed through graphical illustrations. The main findings with future research directions are described at last.


Introduction
Predation is an important factor in a food chain. In predation, predators may or may not kill their prey prior to feeding on them, but the act of predation often results in the death of its prey through consumption. Moreover, predation is an interaction between the species in which one species uses another species as food. Generally, successful predation leads to increase in a population size of predator and to decrease in a population size of prey ( [1], [2], [3], [4]). In this way, predation controls a food chain as well as an ecological system. To describe this phenomenon, many researchers ( [5], [6], [7]) have shown their interest in this eld. Holling-Tanner prey-predator model with Beddington-DeAngelis functional response including delay has been discussed by Jana and Roy [8]. Predator-prey system with prey-taxis has been considered by Qiu et al. [9]. Tripathi et al. [10] analyzed global dynamics and parameter identi ability in a predator-prey interaction model. Banerjee et al. [11] discussed a prey-predator model with Holling Type-III functional response. Chaotic behavior of predator-prey model with group defense and non-linear harvesting in prey has been discussed by Kumar and Kharbanda [12]. Abid et al. [13] discussed about optimality in selective harvesting predator-prey model with modi ed Leslie-Gower and Holling Type II. Harvesting is the process of catching or killing those substances which are used to prepare foods. To ll up the demand of human needs and balance the species for future harvesting, harvesting e ort is very much essential to optimize economically. Clark [14] considered the harvesting of one species in a shery consisting of two competing species. Mathematical modelling in harvesting of such sheries were studied by Chaudhuri et al. [15], Roy and Roy [16], Kar and Ghosh [17]. A harvesting strategy for a prey-predator dynamical system was studied by Sen et al. [18].
A nature reserve is a protected area which is very important for an ecological system. It may be designated by government institutions in some countries or by private land owners. In natural reserve region, limited number of species lived, all species cannot enter. In last few decades, researchers ( [19], [20]) have shown their interest to study prey-predator system in a nature reserve region. They considered only one species i.e., prey lives in a reserve region in their works. But, many cases, there are more than one species live in a reserve region and they are related with prey-predator relationship. We consider here a real-life example on our proposed problem for showing the feasibility and e ectiveness of this work. In geomorphology, drainage systems are the patterns formed by the streams, rivers, and lakes in a particular drainage basin. Mosquito larvae grow in drainage system from where di erent types of mosquitoes are born. They are the carrier of di erent types of harmful diseases like Dengue, Malaria, Meningitis etc. Nowadays, harmful diseases are attacking to human population. A large number of peoples is a ected by these harmful diseases. Controlling of mosquito population is an essential measure to prevent the harmful diseases. Mosquito shes (Gambezi) live and also cultivate in a drainage system to prevent the mosquito populations. Here, gambezi is considered as predator and mosquito larva is treated as prey. A drainage basin is an area of land where all surface water from drainage system usually fall into the basin, as well as the water connects another body of water, such as a river, lake, reservoir, estuary, wetland, sea or ocean. Fishes like cat sh lives in the other body of water like river, lake, reservoir etc. Cat sh is considered here as generalist predator, which consumes mosquito shes for their diet but cat sh does not enter in drainage system. For this reason, drainage system is here the reserve region for the mosquito sh. A discrete model for releases of sterile mosquitoes with Beverton-Holt-type of survivability has been described by Li and Li [21]. Many researchers ( [19], [20]) have paid their attention on di erent regions such as refuge region and predatory region and shown their impacts on prey-predator in these regions. To the best of our knowledge, for the rst time, we consider the interaction among the generalist predator, predator and prey in reserve region. Here are two regions are considered, one is drainage system and other is form by lake, stream and river etc. where the drainage opens. The mosquito larvae generally grow in drainage system. To control their populations we have to cultivate gambezi in drainage system. But the drainage system is not a closed system, it opens at lake, river, stream etc. So gambezi comes from the drainage system. When they come out, we treat them as prey for the generalized predator like cat sh. This is the new consideration in this paper, which is shown mathematically and analyzed the obtained results. In addition to the above, we choose the generalist predator for harvesting e ort in this paper which is more realistic to analyze the whole system. Motivated by this idea, in this study, a mathematical model is designed and then analyzed. The rest of the paper is sorted out as follows: mathematical model is presented in Section 2. In Sections 3 and 4, we describe the equilibria and their existence criteria and stability, respectively of the designed model. Optimal control is introduced to analyze the proposed model in Section 5. In order to verify the model, numerical simulation is presented in Section 6. Finally, conclusion with the future studies are described in Section 7.

Mathematical structure
To analyze the aforementioned situation practically, we choose here drainage basin in which one part is drainage system and the other part is formed by streams, rivers and lakes as shown in Figure 1. The drainage system is opened in a certain place like as river and ocean. There are so many species which do not enter in drainage system like cat sh (generalist predator) which consumes mosquito sh (Gambezi). Assuming that, at any time t, x and x represent the mosquito larva populations (preys) and gambezi populations (predators), respectively. Here, for mosquito sh, drainage system is a reserve region. Considering that, at any time t, x and x denote the gambezi population and cat sh population in other region, respectively. Generally, the growth rate of a predator in the reserve region is di erent from other area. Due to this reason, di erent growth rates are considered for the predator in two di erent regions. Since the mosquito larva is considered as prey and gambezi as predator, the predator population consumes the prey population with Holling type-II functional response or Michaelis-Menten functional response which is chosen as x a +x . The generalist predator population i.e., cat sh consumes the gambezi, so predator population with Beddington-DeAngelis functional response is treated as β x a+bx +mx , where β, a, b and m are positive integers. In a reserved region, we select two regions which are connected to each other. So, migration as well as emigration should occur between the considered regions. Considering that, at any time t, m is per unit emigration and m denotes per unit migration of the predator population in the reserve region. The predator population can interact with prey population in all regions but the generalist predator population cannot interact with predator population in the reserve region. The generalist predators i.e., cat shes are harvested for human needs. Without loss of generality, at any time t, the harvest rate is denoted by h(t) and is de ned by h(t) = qEx , where q denotes as catchability co-e cient and E represents the shing e ort for harvesting the generalist predator population. So, the system of equations is described as follows: with initial conditions x ( ) ≥ , x ( ) ≥ , x ( ) ≥ and x ( ) ≥ . Since x and x are both predator populations, one is in the reserve region and another is in other region, respectively. To simplicity of the calculation and the predators x and x are the same species, the term of functional response of predator on prey in reserve region and other region is treated as same. So, the system of equations (1) is described as follows: with initial conditions x ( ) ≥ , x ( ) ≥ , x ( ) ≥ , x ( ) ≥ and α = γ. Description of the parameters which are mentioned in mathematical structure of the model is given in Table .

Equilibria and their existence criteria
To analyze the system (2) at its interior equilibrium point we rst try to nd all possible non-negative equilibria of the system, though we particularly, interest about the analysis around the interior equilibrium point P(x ,x ,x ,x ) to the system (2) for its natural importance.
So, the su cient conditions for the system at the equilibrium point (x * , x * , , ) are s > m and r − i.e., logistic growth of predator in reserve region is greater than the migration of predator and growth rate of prey is greater than death rate of prey. (iv) The steady state of system of equations (2) is calculated by solving the equations. Let the interior equilibrium point be denoted by P(x ,x ,x ,x ) wherex ,x ,x andx are the positive roots of the equationṡ x =ẋ =ẋ =ẋ = , whereẋ i = dx i dt , i = , , and .
On solving the equations we have the following: andx is the positive root of the following equation where So, the su cient conditions for the system with a positive interior equilibrium point are r > d , β > b(d +qE) and for all B i (i = , , , , , ) are positive. i.e., condition for a positive equilibrium point, the intrinsic growth rate of prey is greater than the natural death rate of prey, and growth rate of generalist predator is greater than the sum of natural death rate of generalist predator and harvesting rate of generalist predator.

Local stability
The stability criterion of the system is analyzed at the interior equilibrium point. The Jacobian matrix of the system (2) at the interior equilibrium point P is denoted by V and is de ned as follows: where I represents an identity matrix of order 4 and , .
ββ x x (a + bx )(a + mx ) (a + bx + mx ) + r r +m m r αx x (a +x ) + αm a s x x r (a +x ) , In that case, we assume that K is the common positive root of φ(K ) = and ψ(K ) = . Now using Routh-Hurwitz criteria around at the interior equilibrium point, we can now state and prove two theorems for analyzing the local stability of the system (2). Theorem 1: Considering that all Q j (j = , , . . . , ), h h − h h and ψ(K ) are positive. Then at the equilibrium point P(x ,x ,x ,x ) of the system (2) is locally asymptotically stable, i.e., after a resilience the system comes to the equilibrium state. Proof: Using Routh-Hurwitz criteria, the conclusion becomes all eigenvalues of the system (2) around its interior equilibrium point P(x ,x ,x ,x ) have negative real parts. Consequently the system will be locally asymptotically stable. This completes the proof of the theorem. Lemma 1: In Theorem 1, it is proved that for a known parameter set, we nd all Q j in terms of K . In this case ifK is the only common positive root. Then for K >K , h h − h h and ψ(K ) are positive. If h h − h h and ψ(K ) are negative then the system (2) will be unstable aroundP(x ,x ,x ,x ). Also if ψ(K ) = , then the system (2) undergoes through a bifurcation. Again K = K , thus we have, for K <K, h h − h h and ψ(K ) are positive. If h h − h h and ψ(K ) are negative then the system (2) must be unstable aroundP(x ,x ,x ,x ). Also if ψ(K ) = , then the system (2) undergoes through a bifurcation. After a certain time by cultivating gambezi in a drainage system, mosquito larvae will be controlled automatically for K >K . Theorem 2: System (2) undergoes through Hopf bifurcation aroundP(x ,x ,x ,x ) for K =K . Proof: For K =K , we have ψ(K ) = and then the eigenvalues of the system atP(x ,x ,x ,x ) which can be expressed as λ , = ϕ (K ) ± iϕ (K ) and λ , = ϕ (K ) ± iϕ (K ). Now it is easy to show that dϕ dK is non-zero at the point K =K . Also we see ϕ(K ) = and ∂ ∂K tr(V|P) K =K = −rx < . Hence, using the given conditions in Venkatsubramanian et al. [22], it is clear to show that our system (2) undergoes a Hopf bifurcation at its interior equilibrium point for the critical value of K , i.e., for K =K . This completes the proof of the theorem. Now we present the condition that the system will be globally stable at the equilibrium point P(x ,x ,x ,x ).
x −x x dx where p i (i = , , ) are suitable positive constants to be determined in the following subsequent steps. Time derivative of the equation along the solutions of the system (2) is given by Then the above expression reduces as follows: (a +x )(a +x ) , then from the above expression we conclude that dW dt ≤ . Hence the theorem is proved.

Optimal control
For a better economical point of the shery as well as the economic development of human life, we optimize the pro t in bionomic equilibrium state. The bionomic equilibrium is a combination of economic equilibrium as well as biological equilibrium. Assuming that p is a constant price per unit biomass, c is a constant cost of harvesting e ort. Then the net economic revenue calculated from the shery is P (x , x , x , x , E, t) = the total revenue obtained by selling the harvested biomass − the total cost for e ort devoted to harvesting = pqEx − cE, then we consider an integral J which is a continuous time-stream of revenues as follows: where δ denotes the instantaneous annual rate of discount [2]. Our problem is to maximize J subject to the state equations (2) using Pontryagin's maximum principle [23]. The control variable E(t) is subject to the constraint set, ≤ E ≤ Emax. At rst, we construct the corresponding Hamiltonian function of this optimal control problem as stated below: where λ i (i = , , , ) are called the adjoint variables. By Pontryagin's maximum principle, the adjoint equations are described as follows: Now we derive an optimal equilibrium solution of the problem at the interior equilibrium point P(x ,x ,x ,x ). Then from equations (9) to (12), we get where Using the values of λ , λ and λ in equations (13) to (16), we get Hamiltonian function by which we can calculate the optimality using Pontryagin's maximum principle. The numerical illustration of the system (2) is described in the next Section i.e., Section 6.    Describing the discussed phenomenon with their relative factors, we justify the analytical results of the proposed system with the hypothetical data set. The results of the simulation experiment are taken using MATLAB 7.10. Choosing the parameter set arbitrarily as P = {r = . , α = . , a = , d = . , m = . , m = . , s = . , s = . , β = . , β = . , d = . , a = , b = . , m = . , q = . , E = . }. Now, K can be taken as (i) K(= ) <K and (ii) K(= ) >K. We consider two cases as follows: Case (i), we draw Figures 2, 3, 4 and 5 by choosing the parameter set P . From Figure 2, we see that the proposed system is stable after certain time. Figure 3 shows the phase space diagram of the system with  respect to x , x and x . Similarly, Figures 4 and 5 show the phase diagram of the system (2) with respect to x , x , x and x , x , x , respectively. Case (ii), we depict Figures 6, 7, 8 and 9 by taking the parameter set P . From Figure 6, we observe that the proposed system is unstable. Figure 7 shows the phase space diagram of the system with respect to x , x and x . Similarly, Figures 8 and 9 show the phase diagram of the system (2) with respect to x , x , x and x , x , x , respectively. Figure 10 shows the bifurcation diagram of the system (2) with the parameter set P . For existence of interior equilibrium point P(x ,x ,x ,x ), the conditions r > d , β > b(d + qE) will hold. The parameter set P and K = satisfy these conditions. By choosing the parameter set P , Figure 2 is drawn which shows stable equilibrium of the system. Now we consider r(= . ) < d (= ) and then draw Figure 11, which shows unstable equilibrium. Again, considering P = {r = . , α = . , a = , d = , m = . , m = . , s = . , s = . , β = . , β = . , d = . , a = , b = , m = . , q = . , E = . } as the set of parameters, where β < b(d + qE), we depict Figure 12. From Figures 11 and 12, it is easily seen that some populations become endanger species.

Conclusion
In this paper, we have formulated the mathematical model on prey-predator relationship in a drainage basin. To the best of our knowledge, the presented work is the rst attempt to consider reserve region for the predator population. Based on the discussed phenomena, we have considered a food chain with a reserve region of predator. A study of potential e ects of generalist predator in a predator and predator in a prey has been introduced with Beddington-DeAngelis functional response and Holling type II functional response, respectively. The density-dependent mortality rates for prey and generalist predator have been assumed. Depending on di erent environments, we have considered di erent consumption rates of predator populations. The equilibria of the proposed system have been determined with the discussions on local and global stability for the   system. We optimize the pro t in bionomic equilibrium state. With the help of this control parameter, a dynamic framework is developed to investigate the optimal utilization of resources, sustainability properties of the stock and the resource rent. Existence criteria at interior equilibrium point and stability of the proposed model has been discussed in numerical simulation and analyzed through the geometrical gures.
The main theoretical results have been demonstrated by numerical simulation with hypothetical parameter values. Someone can verify the validity of the system by considering real data sets. Mosquito-borne diseases control model can be developed for reducing mosquito larva. Hence, future research work is necessary to achieve the needs in this direction.