# Dynamical Analysis of a Beddington–DeAngelis Interacting Species System with Prey Harvesting

Authors: Lakshmi Narayan Guin; Reeta Murmu; Hunki Baek; Kyoung-Hwan Kim
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

A reaction–diffusion interacting species system with Beddington–DeAngelis functional response that has been proposed in the environment of mathematical ecology, which provides the rise to spatial pattern formation, is investigated and associated with the models of deterministic dynamics. The dynamical behaviour of a generalist predator–prey system with linear harvesting of each species and predator-dependent functional response is fully analyzed. Conditions of stability behaviour of the interior equilibrium point are established properly. Furthermore, we have recognized that the unique positive equilibrium point of the system is globally stable via appropriate Lyapunov function structure, which signifies that appropriate harvesting has no impact on the persistence property of the harvesting system. Also, we establish the conditions for the existence of bifurcation phenomena including a saddle-node bifurcation and a Hopf bifurcation. Subsequently, complete analysis regarding the impact of harvesting is carried out, and an interesting decision is that under some appropriate constraints, harvesting has immense impact on the final size of the interacting species. In addition, in accordance with Turing’s ideas on morphogenesis , our analysis shows that harvesting effort in a reaction–diffusion predator–prey system plays a vital function for geological conservation of interacting species. Finally, we discuss sufficient conditions for the existence of bionomic equilibrium point and the optimal harvesting policy attained by using the Pontryagin maximal principle.


---

## Body

## 1. Introduction

There has been great attention during the last few decades in dynamical features of population models and along with these models, predator–prey systems [2–6] play an imperative function in population dynamics. The dynamical correlation between predators and preys has been a vital topic in theoretical ecology since the well-known Lotka–Volterra interacting species system. The vital elements in these forms of models are the predator functional response on prey species per unit time for given sizes (densities) of prey and predator. Different structures of functional responses have become the focus of considerable attention from time to time in ecological literature, for instance, prey-dependent functional response, predator-dependent functional response, ratio-dependent functional response, etc. To originate a more rational ecological system, an appropriate and correct variety of functional response is important.In the current literature, usually, there are two main types of predators, namely, generalist and specialist. Specialist predators feed more or less entirely on one prey species. However, the generalist predator is rather special from specialist predator. Generalist predators feed on various categories of prey species. Once a focal prey species is exposed to extinction, the predator species is capable of changing its food habit to a different species and may prolong to persist; accordingly, their dynamics is not controlled by the activities of a particular prey species [7–10].Mathematical ecology aims at the mathematical depiction, behaviour, and modelling of complex ecosystems using a choice of applied mathematical techniques and tools to get a better understanding of the system [3, 9, 11–13]. Predator–prey relationship plays a significant role in regulating the numbers of prey and predators in the environment of complex ecosystem [14, 15]. The harvesting used by both the species has a strong impact on the system dynamics of the predator–prey interaction.In ecology, to maintain the persistence of the species or the coexistence of the species, there are different feasible techniques such as appropriate amount of harvesting policy, constructing reserved zones or refuges, etc., which should be taken in order to help the species grow in the regions without annoyances. The functions of such kind of measures were investigated by different eminent scholars such as suitable harvesting policy in predator–prey dynamics which attracts the major attention of the researchers. By constructing a suitable Lyapunov function, some scholars proposed that the coexistence equilibrium point of the system is globally stable which reflects on the event that suitable harvesting has no influence on the persistence property of the harvesting system, whereas under some appropriate parametric restrictions, harvesting has no influence on the final size of the prey species but the size of the predator species is strictly decreasing function of the harvesting efforts [10, 16–18].Some eminent evolutionary biologists and theoretical ecologists have recognized that the interactions between predator and their preys should not be simply explained by direct predation through prey-dependent functional responses [16, 19–21] and that the importance of predator-dependent functional responses should be considered [18, 22, 23]; so far our knowledge goes, no mathematical models have been projected to investigate whether or the extent to which harvested generalist predator–prey model with Beddington–DeAngelis type functional response can affect interacting populations. In this investigation, the interactions of two populations have been proposed and analyzed incorporating harvesting phenomena on both species to explore the effect that harvesting can have on reaction–diffusion population dynamics in generalist predator–prey systems. Our goal in this research is to investigate the spatiotemporal dynamics of a harvested reaction–diffusion system. We develop the existence of limit cycle through Hopf bifurcation theorem and spatial patterns through diffusion-driven instability, by using the harvesting coefficient as a parameter of bifurcation.In the present investigation, we have categorized the research work as follows. In Section2, we propose the basic assumptions and then accordingly formulated the mathematical model. With the nondimensional model, we checked for system’s (2) existence, boundedness, dissipativity, and permanence. In Section 3, we look for equilibria, and their dynamical behaviours are well studied herein. Conditions for stability are achieved and phase portraits are plotted numerically. System (2) is found to have Hopf bifurcation and saddle-node bifurcation about coexistence equilibrium point. Following in Section 4, the effect of diffusion is reflected on the system and the related pattern formation is investigated. In Discussion section, we have talked about bionomic equilibrium and optimal harvesting policy. Finally, the paper ends with a brief conclusion and significance of the obtained results in Section 6.


## 2. Mathematical Models

Assuming the importance of harvesting on a generalist predator–prey model with Beddington–DeAngelis functional response, an effort is prepared to study the influence of linear harvesting of each species in a two-dimensional2D Beddington–DeAngelis predator–prey model which is a modification of the following model investigated extensively by Haque [24] and Sarwadi et al. [25].(1a)dudt=r1u1−uK1−auvb+u+cv−H1u,(1b)dvdt=r2v1−vK2+eauvb+u+cv−H2v,(1c)u0≥0,v0≥0,where ut and vt represents the population sizes of prey and generalist predator, respectively, at any instant of time t, and other system parameters, viz., r1,r2,K1,K2,e,a,b,c,H1,H2 are all positive. Here, r1 and r2 are the intrinsic growth rates of the prey and predator, respectively; K1 and K2 are the environmental carrying capacities of the prey and predator, respectively; e is the conversion factor of prey to predator; a is the maximum consumption rate and the parameter b is a normalization coefficient that relates the sizes of the predator and prey to the environment in which they interact; c scales the impact of the predator interference; c<0 is the case where predators benefits from cofeeding; H1 and H2 are the linear harvesting rate of prey and predator species, respectively, which means that the number of individuals harvested per unit of time is proportional to the current population 0≤H1<r1,0≤H2<r2. The predator species consumes the prey with functional response of Beddington–DeAngelis type auv/b+u+cv and contributes to its growth with eauv/b+u+av [26, 27]. Throughout the investigation, the environment is believed to be uniform; that is, the ecological parameters do not depend on space and time.Using nondimensional variables such asu˜=u/K1, v˜=v/K2, and t˜=ta, the corresponding model can be written in dimensionless form as (dropping tildes)(2a)dudt=u1−u−εuvβ+u+αv−δu≡F1u,v,(2b)dvdt=γv1−v+λuvβ+u+αv−σv≡F2u,v,(2c)u0≥0,v0≥0,where α=cK2/K1, β=b/K1, γ=r2/r1, ε=aK2/r1K1, δ=H1/r1, λ=ea/r1, and σ=H2/r1.Now, we study the boundedness of the solutions and dissipative behaviour of the system (2). Permanence is an important behaviour of the system in the sense that it describes long-term behaviour of the system. Analytically, a system is said to be permanent if it is dissipative; i.e., all trajectories of the system (2) initiating from the point in IR+2 ultimately lie in the fixed bounded region for sufficiently large t.Lemma 1.
Every solution of the system (2) exists in the interval 0,∞ and ut≥0, vt≥0 for all t≥0.Proof.
SinceF1u,v and F2u,v are completely continuous and locally Lipschitzian on C1IR+2, the solution ut,vt of (2) exits and is unique on 0,ξ, where 0<ξ<+∞. From the system (2), we have(3)ut=u0exp∫0t1−us−δ−εvsβ+us+αvsds≥0,vt=v0exp∫0tγ1−vs−σ+λusβ+us+αvsds≥0,which completes the proof.Proposition 1.
All the solutions of the system (2) in IR+2 are uniformly bounded.Proof.
Define a functionηt=ut+ε/λvt. Differentiating ηt with respect to t, we get(4)dηdt=dudt+ελdvdt=u1−u−δu+γελv1−v−σελv.
For eachξ>0, we have(5)dηdt+ξη=u1−u−δu+γελv1−v−σελv+ξu+ξελv=u1−u−δ+ξ+ελvγ1−v−σ+ξ=u1−δ+ξ−u+ελvγ−σ+ξ−γv=1−δ−ξ24+ε4λγγ−σ+ξ2.
Hence, one can findϕ>0, such that dη/dt+εη≤ϕ. From the above, we find(6)ηt≤e−ξtη0+ϕξ1−e−ξt≤maxη0,ϕξ.
Sinceu0≥0,v0≥0,η0=u0+ε/λv0 and letting t⟶0, we have limsupηt≤ϕ/ξ, which is independent of the initial condition. We observe that all the solutions of the system initiating in IR+2 eventually lie in the region S defined by(7)S=u,v∈IR+2:η=ϕξ+κ,for someκ>0.
Hence, the proof is complete.Proposition 2.
The system (2) is dissipative if γ+λ>σ.Proof.
By usual straightforward arguments, one can show that the solution of the system (2) always exists and keeps up positive. Since all the interacting populations and parameters associated with the system (2) are positive, for the first equation of system (2), we have du/dt≤u1−u. Consequently, by solving, one can observe that(8)ut≤11/u0−1e−t+1.
Thus, any solution of the system (2) must satisfy limsupt⟶∞ut≤1.
Now, from the second equation of the system (2), we examine that(9)dvdt≤vγ1−v−σ+λ.
Therefore, we havelimsupt⟶∞vt≤γ+λ−σ/γ=v¯, where we denote by v¯ a unique value of vt which is positive if γ+λ>σ. Hence, the proof is complete.Proposition 3.
The system (2) is permanent if ε/α+δ<1 and γ+λ>σ.Proof.
From the prey equation of system (2), we have(10)dudt=u1−u−εuvβ+u+αv−δu≥u1−u−δ−εα.
Ifu¯ is the root of 1−u¯−δ−ε/α=0, then again by standard comparison argument, we observe that the above inequality implies liminft⟶∞ut≤α1−δ−ε/α=u¯.
Hence, for larget, we have ut≥u¯.
Again, from the predator system of (2), we have(11)dvdt=γv1−v+λuvβ+u+αv−σv≥vγ1−v−σ+λ.
Ifv¯ is the positive root of the quadratic equation γ1−v−σ+λ=0, then(12)v¯=−ασ−γ+γu+β+ασ−γ+γu+β2−4γαβσ−γ+uσ−γ−λ2γα.
Therefore, we obtainliminft⟶∞vt≥v¯.
Note thatu¯>0 if α>0 and ε/α+δ<1; again v¯>0 if γ+λ>σ. Therefore, combining these along with Proposition 2 yields the conclusion.
Hence, the proof is complete.


## 3. Equilibria and Their Stability

This section compacts with stability and bifurcation analysis of the proposed system. First, we analyze the behaviour of the proposed system through the existence of equilibrium of the system (2). The corresponding system (2) has at most four feasible equilibria in the first quadrant IR+2 as follows:(13)ie0u0,v0,iie1u1,v1,iiie2u2,v2,ive3u3,v3,where u0=0,v0=0,u1=1−δ, v1=0,u2=0,v2=1−σ/γ, u3=β+αv3σ−γ1−v3/λ−σ+γ1−v3; 0<σ−γ1−v3<λ and v3 are the roots of the following cubic equation:(14)a0x3+a1x2+a2x+a3=0,a0≠0,where(15)a0=α2λγ+εγ2>0,a1=2εγσ−λ−α2λγ−σ−2εγ2+γαλ2β−δ+1,a2=2εγλ−αλ2+λσα−δλγβ+ελ2+εσ2+β2λγ−λγα+λγβ−2εσλ+λ2δα−2εγσ+δλγα+εγ2−2βλγα+2βλσα−δλσα,a3=β2λσ−γ+δβλ2+γλ−λσ+λβσ−γ.Using Descartes’ rule of signs, one may identify that the cubic equation has one and only one positive root if any one of the following inequalities holds:(16)aa1<0,a2<0,a3<0,ba1>0,a2<0,a3<0,ca1>0,a2>0,a3<0.As a result of the ecological significance, the main concentration has been focused on investigating the qualitative dynamics around the unique and feasible coexistence equilibrium pointe3u3,v3.Example 1.
In Figure1, we have shown the coexistence equilibrium point marked by the green circle which is the point of intersection of the prey zero-growth isoclines and the predator zero-growth isoclines for the system (2) related to the parameter values ε=1.0, β=0.0041, α=0.91881, δ=0.2, γ=0.0327, λ=0.15, and σ=0.1. The prey isoclines attains the axis u=0 and the curve 1−u−εv/β+u+αv−δ=0 and the predator isoclines attains v=0 and γ1−v+λu/β+u+αv−σ=0.Figure 1
Isocline diagram of the system (2) around the positive interior point e3 corresponding to the different parameters ε=1.0, β=0.0041, α=0.91881, δ=0.2, γ=0.0327, λ=0.15, and σ=0.1.


### 3.1. Linearization and Stability Analysis

Let us commence the following notations. The system (2) can be written in the form(17)X˙=FX,ε=F1u,v,F2u,vT,where X=u,vT. Here, XT and X˙, respectively, represent the transpose and the derivative with respect to t of X. The Jacobian matrix J≡DFX,ε of the system (2) at any arbitrary point u,v is given by(18)J=1−2u−δ−εvβ+αvβ+u+αv2−εuβ+uβ+u+αv2λvβ+αvβ+u+αv2γ1−v−γv−σ+λuβ+uβ+u+αv2=εvuβ+u+αv2−u−εuβ+uβ+u+αv2λvβ+αvβ+u+αv2−γv−λuαvβ+u+αv2=Γ11Γ12Γ21Γ22.We denoteJk=J, the Jacobian evaluated at ek,i=1,2;j=1,2;k=0,1,2,3 and the determinant Jk=detJk, trace Jk=trJk.The dynamical behaviours of the system (2) around the equilibrium points e0,e1 and e2 are summarized in the following theorems.Theorem 1.
The trivial equilibrium pointe00,0 is(i)
saddle point if1−δ>0,γ−σ<0 or 1−δ<0,γ−σ>0,(ii)
unstable point if1−δ>0,γ−σ>0,(iii)
stable point if1<δ and γ<σ.Proof.
The Jacobian matrixJ0 of the system (2) around the boundary equilibrium point e00,0 is given by(19)J0=1−δ00γ−σ.
The eigenvalues corresponding toJ0 are given by 1−δ and γ−σ. Thus, one can easily identify the dynamical behaviour around the equilibrium point e00,0, via eigenvalues analysis, and accordingly the equilibrium point e00,0 is(i)
saddle point if1−δ>0,γ−σ<0 or 1−δ<0,γ−σ>0(ii)
unstable point if1−δ>0,γ−σ>0(iii)
stable point if1<δ and γ<σExample 2.
If one takes the parameters in the system (2) as ε=1.0,β=0.0041,α=0.91881,δ=1.2,γ=0.0327,λ=0.15, and σ=0.1, one can have 1<δ and γ<σ. Thus, it follows from Theorem 1 that e00,0 is locally stable. Figure 2 shows such phenomena.Figure 2
Local asymptotic stability of the system (2) around the boundary equilibrium point e0u0,v0 corresponding to the parameters ε=1.0,β=0.0041,α=0.91881,δ=1.2,γ=0.0327,λ=0.15,σ=0.1 and varying the initial conditions for prey species (u).Theorem 2.
The boundary equilibrium pointe11−δ,0 is(1)
saddle point if(i)
δ<1 with λ+σ<γ,β+1<δ or λ+γ>σ,β+1>δ(ii)
δ>1 with σ>γ+λ,β+1<δ or γ<σ+λ,β+1>δ(2)
unstable point ifδ>1,σ<γ+λ,β+1<δ or δ>1,λ+σ<γ,β+1>δ(3)
stable point ifδ<1,γ<σ−λ,β+1>δ or δ>1,λ+σ>γ,β+1<δProof.
The Jacobian matrixJ1 of the system (2) around the boundary equilibrium point e11−δ,0 is given by(20)J1=δ−1−ε1−δβ+1−δ0γ+λ1−δβ+1−δ−σ.
The eigenvalues corresponding toJ1 are given by δ−1 and γ+λ1−δ/β+1−δ−σ. Straightforwardly, one can observe that at the equilibrium points e11−δ,0, the system is(1)
saddle if(i)
δ<1 with λ+σ<γ,β+1<δ or λ+γ>σ,β+1>δ(ii)
δ>1 with σ>γ+λ,β+1<δ or γ<σ+λ,β+1>δ(2)
unstable ifδ>1,σ<γ+λ,β+1<δ or δ>1,λ+σ<γ,β+1>δ(3)
stable ifδ<1,γ<σ−λ,β+1>δ or δ>1,λ+σ>γ,β+1<δExample 3.
The system (2) with ε=1.0,β=1.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15, and σ=0.2 satisfies the conditions δ<1,γ<σ−λ and β+1>δ. Thus, equilibrium point e11−δ,0 is locally stable. We illustrate this phenomenon in Figure 3.Figure 3
Local asymptotic stability of the system (2) around the boundary equilibrium point e1u1,v1 corresponding to the parameters ε=1.0,β=1.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15,σ=0.2 satisfying the stability condition δ<1,γ<σ−λ,β+1>δ.Theorem 3.
The axial equilibrium pointe20,1−σ/γ is(1)
saddle point if(i)
σ>γ with γβ>α,γβ+ε<α+δ or γβ<α,γβ+ε>α−δ(ii)
σ<γ with γβ>ασ−γ,γβ−ε>α+δ or γβ<ασ−γ,γβ+α<ε−δ(2)
unstable point ifσ>γ,γβ>α,γβ+ε>α+δ or σ>γ,γβ<α,γβ+ε<α−δ(3)
stable point ifσ<γ,γβ>ασ−γ,γβ+α<ε+δ or σ<γ,γβ<ασ−γ,γβ+α<ε−δProof.
The Jacobian matrixJ2 of the system (2) around the axial equilibrium point e20,1−σ/γ is given by(21)J2=1−ε1−σ/γβ+α1−σ/γ−δ0λ1−σ/γβ+α1−σ/γ−γ1−σ/γ.
The eigenvalues corresponding toJ2 are given by γβ+αγ−ασ+εσ−γ−δγβ+αγ−ασ/γβ+αγ−ασ and σ−γ.
Similarly, using the same argument, one can find out that the equilibrium pointe20,1−σ/γ is(1)
saddle point if(i)
σ>γ with γβ>α,γβ+ε<α+δ or γβ<α,γβ+ε>α−δ(ii)
σ<γ with γβ>ασ−γ,γβ−ε>α+δ or γβ<ασ−γ,γβ+α<ε−δ(2)
unstable point ifσ>γ,γβ>α,γβ+ε>α+δ or σ>γ,γβ<α,γβ+ε<α−δ(3)
stable point ifσ<γ,γβ>ασ−γ,γβ+α<ε+δ or σ<γ,γβ<ασ−γ,γβ+α<ε−δExample 4.
If one chooses the parameters in the system (2) as ε=1.0,β=0.0041,α=0.91881,δ=0.2,γ=0.1327,λ=0.15, and σ=0.1, it is easy to see that the conditions δ<1,γ<σ−λ and β+1>δ are satisfied. Thus, from Theorem 10, one can know that the equilibrium point e20,1−σ/γ is locally stable as shown in Figure 4.
Now we mention the dynamical behaviours of the system (2) around the interior equilibrium point e3u3,v3 as follows.Figure 4
Local asymptotic stability of the system (2) around the boundary equilibrium point e2u2,v2 corresponding to the parameters ε=1.0,β=0.0041,α=0.91881,δ=0.2,γ=0.1327,λ=0.15,σ=0.1 satisfying the stability condition σ<γ,γβ>ασ−γ,γβ+α<ε+δ.Theorem 4.
The coexistence equilibrium pointe3u3,v3 is(i)
locally asymptotically stable ifεv3<β+u3+αv32(ii)
If1−ε/β2v3>0, then the local stability of e3u3,v3 ensures its global stabilityProof.
(i)
In order to prove that the system is locally asymptotically stable around the equilibrium pointe3u3,v3, firstly, we need to observe the Jacobian matrix J3 of the system around the point, which is given by(22)J3=1−2u3−δ−εv3β+αv3β+u3+αv32−εu3β+u3β+u3+αv32λv3β+αv3β+u3+αv32γ1−v3−γv3+λu3β+u3β+u3+αv32−σ=εu3v3β+u3+αv32−u3−εu3β+u3β+u3+αv32λv3β+αv3β+u3+αv32−λu3αv3β+u3+αv32−γv3=Γ11Γ12Γ21Γ22.From the Routh–Hurwitz criterion, given as trJ3<0 and det J3>0, we can get Γ11<0, which gives εv3<β+u3+αv32.(ii)
To show the global stability of the system (2), we construct a Lyapunov function as follows:(23)Lu,v=L1∫u3uξ1−u3ξ1dξ1+L2∫v3vξ2−v3ξ2dξ2=L1u−u3−u3loguu3+L2v−v3−v3logvv3,where L1 and L2 are positive constants to be determined. Clearly, Lu,v=0 at e3u3,v3 and is positive for all other positive values of u and v. The time derivative of Lu,v along the solution of the system (2) is(24)dLdt=L1dudt−u3ududt+L2dvdt−v3vdvdt=L1u−u3ududt+L2v−v3vdvdt=L1u−u3uu1−u−δu−εuvβ+u+αv+L2v−v3vγv1−v−σv+λuvβ+u+αv=L1u−u31−u−δ−εvβ+u+αv+L2v−v3γ1−v−σ+λuβ+u+αv=L1u−u3u3+εv3β+u3+αv3−u−εvβ+u+αv+L2v−v3γv3−λu3β+u3+αv3−γv+λuβ+u+αv=L1u−u3−u−u3+εβ+u+αvβ+u3+αv3βv3+uv3−vβ−vu3+L2v−v3−γv−v3+λβ+u+αvβ+u3+αv3βu+αuv3−u3β−αvu3≤L1u−u3−u−u3+εβ2βv3−v+uv3−vu3+L2v−v3−γv−v3+λβ2βu−u3+αuv3−vu3=−L1u−u32−L1u−u3εβv−v3+L1u−u3εβ2uv3−vu3−γL2v−v32+L2v−v3λβu−u3+L2v−v3λαβ2uv3−vu3=−L1u−u32−L1εβu−u3v−v3+L1εβ2v3u−u32−L1εβ2u3u−u3v−v3−γL2v−v32+L2λβu−u3v−v3+L2λαβ2v3u−u3v−v3−L2λαβ2u3v−v32=−L1u−u321−εβ2v3−L2v−v32γ+λαβ2u3−u−u3v−v3L1εβ1+u3β−L2λβ1+αv3β2,
ChoosingL1=1 and L2=εβ+u3/λβ+αv3, one can write(25)dLdt≤−u−u321−εβ2v3+v−v32γ+λαβ2u3εβ+u3λβ+αv3.
Thus, if1−ε/β2v3>0, then dL/dt<0 and dL/dt=0 if u,v=u3,v3. Hence, the proof is complete.Example 5.
Numerically, we display the solution of the system (2) with ε=1.0,β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15, and σ=0.1 in Figure 5. It follows from Theorem 4 that the equilibrium point e30.2753798784,0.2830660001 is globally stable due to 1−ε/β2v3>0.Figure 5
Global stability of the system (2) around e30.2753798784,0.2830660001 corresponding to ε=1.0,β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15,σ=0.1.


### 3.2. Bifurcation Analysis of the System

Now, we investigate bifurcation phenomena of the system according to the parameter valueε.Theorem 5.
The system (2) undergoes a saddle-node bifurcation around e3u3,v3 at ε=εSN, where(26)εSN=γv33α3+3α2v33γβ+u3+α2λu3v3+3αv3γβ2+βu3+u32+αu3λu3+β+βγβ2+3βu3+3u32+u33γαγv32+u3+βv3γ−λβ.Proof.
For a saddle-node bifurcation, one of the eigenvalues ofJ3 has to be zero which is possible iff det J3=Γ11Γ22−Γ12Γ21=0, which gives ε=εSN. The other eigenvalue is obtained from the relation tr J3=Γ11+Γ22 evaluated at ε=εSN and should be negative [28].
Let℧ and Υ be the eigenvectors corresponding to eigenvalue 0 of the matrix J3 and its transpose, respectively. We obtain that ℧=θ1,θ2T,Υ=ℏ1,ℏ2T where(27)θ1=−Γ12θ2Γ11=−Γ22θ2Γ21,ℏ1=−Γ21ℏ2Γ11=−Γ22ℏ2Γ12.
Andθ2,ℏ2 are any two nonzero real numbers.
Now, for the system (2), we can get(28)ΥTFεX3,εSN=ℏ1ℏ2∂F1∂ε∂F2∂ε=ℏ1ℏ2−uvβ+u+αv0=−ℏ1uvβ+u+αv≠0ΥTD2FX3,εSN℧,℧=ℏ1ℏ2∂2F1∂u2θ12+∂2F1∂u∂vθ1θ2+∂2F1∂v∂uθ2θ1+∂2F1∂v2θ22∂2F2∂u2θ12+∂2F2∂u∂vθ1θ2+∂2F2∂v∂uθ2θ1+∂2F2∂v2θ22=ℏ1∂2F1∂u2θ12+∂2F1∂u∂vθ1θ2+∂2F1∂v∂uθ2θ1+∂2F1∂v2θ22+ℏ2∂2F2∂u2θ12+∂2F2∂u∂vθ1θ2+∂2F2∂v∂uθ2θ1+∂2F2∂v2θ22=ℏ1∂∂u1−2u−δ−εvβ+αvβ+u+αv2θ12−∂∂uεuβ+uβ+u+αv2θ1θ2+∂∂v1−2u−δ−εvβ+αvβ+u+αv2θ2θ1−∂∂vεuβ+uβ+u+αv2θ22+ℏ2∂∂uλvβ+αvβ+u+αv2θ12+∂∂uγ1−v−γv+λuβ+uβ+u+αv2−σθ1θ2+∂∂vλvβ+αvβ+u+αv2θ2θ1+∂∂vγ1−v−γv+λuβ+uβ+u+αv2−σθ22=ℏ1−2+2εvβ+αvβ+u+αv3θ12−2εβ2+βu+βαv+2uαvβ+u+αv3θ1θ2+2εαuβ+uβ+u+αv3θ22−ℏ22λvβ+αvβ+u+αv3θ12+λβ2+βu+βαv+2uvαβ+u+αv3θ1θ2−2γ+2uλαβ+uβ+u+αv3θ22=−2θ12ℏ1+2γθ22ℏ2+1β+u+αv3ℏ12εvβ+vαθ12−2εβ2+uβ+vβα+2uvαθ1θ2+2uεαβ+uθ22−ℏ22vλβ+vαθ22−2λβ2+uβ+vβα+2uvαθ1θ2−2uλαβ+uθ22≠0.
Here, we writeu,v instead of u3,v3 for algebraic convenience.
So, the system experiences a saddle-node bifurcation [28] around the positive interior equilibrium e3u3,v3 at ε=εSN. Again, since ΥTFεX2,εSN≠0, therefore, the system does not attain neither transcritical bifurcation nor pitch-fork bifurcation.Theorem 6.
The system (2) admits a Hopf bifurcation around e3u3,v3 at ε=εHB, where ε=εHB=λα+u3+γv3β+u3+αv32/u3v3.Proof.
In a system, a Hopf bifurcation typically occurs when a complex conjugate pair of eigenvalues of the linearized flow at the fixed pointe3u3,v3 becomes purely imaginary, i.e., if det J3>0 and if tr J3=0 for some value of ε say εHB.
The characteristic equation is given by(29)μ2−Γ11+Γ22μ+Γ11Γ22−Γ12Γ21=0,where u−u3≈eμt and v−v3≈eμt.
If trJ3=Γ11+Γ22=0 at ε=εHB, then both the eigenvalues will be purely imaginary if det J3=Γ11Γ22−Γ12Γ21 is positive.
Replacingμ=μ1+iμ2 into the corresponding characteristic equation and separating real and imaginary parts, we get(30)μ12−μ22−trJ3μ1+detJ3=0,(31)2μ1μ2−trJ3μ2=0.
Elementary differential of (30) with respect to ε and considering μ1=0, we get(32)dμ1dεε=εHB=v3u32β+u3+αv32≠0,
Hence, the system goes through a Hopf bifurcation arounde3u3,v3 at ε=εHB, where(33)εHB=λα+u3+γv3β+u3+αv32u3v3.Example 6.
It follows from Theorem6 that the system could have a limit cycle on account of Hopf bifurcation. In fact, Figure 6 shows the existence of a limit cycle numerically when the parameters are taken as ε=1.0596222,β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15, and σ=0.1 in the system (2).Figure 6
A limit cycle of the system (2) around e30.2371368969,0.2503112874 corresponding to ε=1.0596222,β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15,σ=0.1.Remark 1.
For better understanding of qualitative change of the temporal dynamics of the system (2), one may observe the bifurcation diagrams of both prey and predator species with respect to the system parameter ϵ (cf. Figures 7 and 8). In support of this technique of bifurcation offered in Figure 7, the successive maxima of u and v in the ranges 0.01,0.8 and 0.01,0.33, respectively, as a function of ϵ in the range 0.95≤ϵ≤1.2 and the other system parameters are captioned in the figure. The other bifurcation diagram is generated for the successive maxima of population sizes of u and v in the ranges 0.01,0.9 and 0.0,0.7, respectively, as a function of ϵ in the range 0.8≤ϵ≤1.5 and the other system parameters are captioned in Figure 8. It is interesting to note from Figures 7 and 8 that maximum concentration of predator species reduces drastically in the presence of prey harvesting factor δ at positive level. Due to the certain amount of prey harvesting from the interacting species system, predator species cannot consume prey species in earlier manner and accordingly this scenario supports the natural phenomena very well.Figure 7
The bifurcation of the system (2) around e3u3,v3 corresponding to the bifurcating parameter ε, where other parameter values are β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15,σ=0.1.Figure 8
The bifurcation of the system (2) around e3u3,v3 corresponding to the bifurcating parameter ε, where other parameter values are β=0.031,α=0.871881,δ=0.0,γ=0.0327,λ=0.15,σ=0.05.


## 3.1. Linearization and Stability Analysis

Let us commence the following notations. The system (2) can be written in the form(17)X˙=FX,ε=F1u,v,F2u,vT,where X=u,vT. Here, XT and X˙, respectively, represent the transpose and the derivative with respect to t of X. The Jacobian matrix J≡DFX,ε of the system (2) at any arbitrary point u,v is given by(18)J=1−2u−δ−εvβ+αvβ+u+αv2−εuβ+uβ+u+αv2λvβ+αvβ+u+αv2γ1−v−γv−σ+λuβ+uβ+u+αv2=εvuβ+u+αv2−u−εuβ+uβ+u+αv2λvβ+αvβ+u+αv2−γv−λuαvβ+u+αv2=Γ11Γ12Γ21Γ22.We denoteJk=J, the Jacobian evaluated at ek,i=1,2;j=1,2;k=0,1,2,3 and the determinant Jk=detJk, trace Jk=trJk.The dynamical behaviours of the system (2) around the equilibrium points e0,e1 and e2 are summarized in the following theorems.Theorem 1.
The trivial equilibrium pointe00,0 is(i)
saddle point if1−δ>0,γ−σ<0 or 1−δ<0,γ−σ>0,(ii)
unstable point if1−δ>0,γ−σ>0,(iii)
stable point if1<δ and γ<σ.Proof.
The Jacobian matrixJ0 of the system (2) around the boundary equilibrium point e00,0 is given by(19)J0=1−δ00γ−σ.
The eigenvalues corresponding toJ0 are given by 1−δ and γ−σ. Thus, one can easily identify the dynamical behaviour around the equilibrium point e00,0, via eigenvalues analysis, and accordingly the equilibrium point e00,0 is(i)
saddle point if1−δ>0,γ−σ<0 or 1−δ<0,γ−σ>0(ii)
unstable point if1−δ>0,γ−σ>0(iii)
stable point if1<δ and γ<σExample 2.
If one takes the parameters in the system (2) as ε=1.0,β=0.0041,α=0.91881,δ=1.2,γ=0.0327,λ=0.15, and σ=0.1, one can have 1<δ and γ<σ. Thus, it follows from Theorem 1 that e00,0 is locally stable. Figure 2 shows such phenomena.Figure 2
Local asymptotic stability of the system (2) around the boundary equilibrium point e0u0,v0 corresponding to the parameters ε=1.0,β=0.0041,α=0.91881,δ=1.2,γ=0.0327,λ=0.15,σ=0.1 and varying the initial conditions for prey species (u).Theorem 2.
The boundary equilibrium pointe11−δ,0 is(1)
saddle point if(i)
δ<1 with λ+σ<γ,β+1<δ or λ+γ>σ,β+1>δ(ii)
δ>1 with σ>γ+λ,β+1<δ or γ<σ+λ,β+1>δ(2)
unstable point ifδ>1,σ<γ+λ,β+1<δ or δ>1,λ+σ<γ,β+1>δ(3)
stable point ifδ<1,γ<σ−λ,β+1>δ or δ>1,λ+σ>γ,β+1<δProof.
The Jacobian matrixJ1 of the system (2) around the boundary equilibrium point e11−δ,0 is given by(20)J1=δ−1−ε1−δβ+1−δ0γ+λ1−δβ+1−δ−σ.
The eigenvalues corresponding toJ1 are given by δ−1 and γ+λ1−δ/β+1−δ−σ. Straightforwardly, one can observe that at the equilibrium points e11−δ,0, the system is(1)
saddle if(i)
δ<1 with λ+σ<γ,β+1<δ or λ+γ>σ,β+1>δ(ii)
δ>1 with σ>γ+λ,β+1<δ or γ<σ+λ,β+1>δ(2)
unstable ifδ>1,σ<γ+λ,β+1<δ or δ>1,λ+σ<γ,β+1>δ(3)
stable ifδ<1,γ<σ−λ,β+1>δ or δ>1,λ+σ>γ,β+1<δExample 3.
The system (2) with ε=1.0,β=1.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15, and σ=0.2 satisfies the conditions δ<1,γ<σ−λ and β+1>δ. Thus, equilibrium point e11−δ,0 is locally stable. We illustrate this phenomenon in Figure 3.Figure 3
Local asymptotic stability of the system (2) around the boundary equilibrium point e1u1,v1 corresponding to the parameters ε=1.0,β=1.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15,σ=0.2 satisfying the stability condition δ<1,γ<σ−λ,β+1>δ.Theorem 3.
The axial equilibrium pointe20,1−σ/γ is(1)
saddle point if(i)
σ>γ with γβ>α,γβ+ε<α+δ or γβ<α,γβ+ε>α−δ(ii)
σ<γ with γβ>ασ−γ,γβ−ε>α+δ or γβ<ασ−γ,γβ+α<ε−δ(2)
unstable point ifσ>γ,γβ>α,γβ+ε>α+δ or σ>γ,γβ<α,γβ+ε<α−δ(3)
stable point ifσ<γ,γβ>ασ−γ,γβ+α<ε+δ or σ<γ,γβ<ασ−γ,γβ+α<ε−δProof.
The Jacobian matrixJ2 of the system (2) around the axial equilibrium point e20,1−σ/γ is given by(21)J2=1−ε1−σ/γβ+α1−σ/γ−δ0λ1−σ/γβ+α1−σ/γ−γ1−σ/γ.
The eigenvalues corresponding toJ2 are given by γβ+αγ−ασ+εσ−γ−δγβ+αγ−ασ/γβ+αγ−ασ and σ−γ.
Similarly, using the same argument, one can find out that the equilibrium pointe20,1−σ/γ is(1)
saddle point if(i)
σ>γ with γβ>α,γβ+ε<α+δ or γβ<α,γβ+ε>α−δ(ii)
σ<γ with γβ>ασ−γ,γβ−ε>α+δ or γβ<ασ−γ,γβ+α<ε−δ(2)
unstable point ifσ>γ,γβ>α,γβ+ε>α+δ or σ>γ,γβ<α,γβ+ε<α−δ(3)
stable point ifσ<γ,γβ>ασ−γ,γβ+α<ε+δ or σ<γ,γβ<ασ−γ,γβ+α<ε−δExample 4.
If one chooses the parameters in the system (2) as ε=1.0,β=0.0041,α=0.91881,δ=0.2,γ=0.1327,λ=0.15, and σ=0.1, it is easy to see that the conditions δ<1,γ<σ−λ and β+1>δ are satisfied. Thus, from Theorem 10, one can know that the equilibrium point e20,1−σ/γ is locally stable as shown in Figure 4.
Now we mention the dynamical behaviours of the system (2) around the interior equilibrium point e3u3,v3 as follows.Figure 4
Local asymptotic stability of the system (2) around the boundary equilibrium point e2u2,v2 corresponding to the parameters ε=1.0,β=0.0041,α=0.91881,δ=0.2,γ=0.1327,λ=0.15,σ=0.1 satisfying the stability condition σ<γ,γβ>ασ−γ,γβ+α<ε+δ.Theorem 4.
The coexistence equilibrium pointe3u3,v3 is(i)
locally asymptotically stable ifεv3<β+u3+αv32(ii)
If1−ε/β2v3>0, then the local stability of e3u3,v3 ensures its global stabilityProof.
(i)
In order to prove that the system is locally asymptotically stable around the equilibrium pointe3u3,v3, firstly, we need to observe the Jacobian matrix J3 of the system around the point, which is given by(22)J3=1−2u3−δ−εv3β+αv3β+u3+αv32−εu3β+u3β+u3+αv32λv3β+αv3β+u3+αv32γ1−v3−γv3+λu3β+u3β+u3+αv32−σ=εu3v3β+u3+αv32−u3−εu3β+u3β+u3+αv32λv3β+αv3β+u3+αv32−λu3αv3β+u3+αv32−γv3=Γ11Γ12Γ21Γ22.From the Routh–Hurwitz criterion, given as trJ3<0 and det J3>0, we can get Γ11<0, which gives εv3<β+u3+αv32.(ii)
To show the global stability of the system (2), we construct a Lyapunov function as follows:(23)Lu,v=L1∫u3uξ1−u3ξ1dξ1+L2∫v3vξ2−v3ξ2dξ2=L1u−u3−u3loguu3+L2v−v3−v3logvv3,where L1 and L2 are positive constants to be determined. Clearly, Lu,v=0 at e3u3,v3 and is positive for all other positive values of u and v. The time derivative of Lu,v along the solution of the system (2) is(24)dLdt=L1dudt−u3ududt+L2dvdt−v3vdvdt=L1u−u3ududt+L2v−v3vdvdt=L1u−u3uu1−u−δu−εuvβ+u+αv+L2v−v3vγv1−v−σv+λuvβ+u+αv=L1u−u31−u−δ−εvβ+u+αv+L2v−v3γ1−v−σ+λuβ+u+αv=L1u−u3u3+εv3β+u3+αv3−u−εvβ+u+αv+L2v−v3γv3−λu3β+u3+αv3−γv+λuβ+u+αv=L1u−u3−u−u3+εβ+u+αvβ+u3+αv3βv3+uv3−vβ−vu3+L2v−v3−γv−v3+λβ+u+αvβ+u3+αv3βu+αuv3−u3β−αvu3≤L1u−u3−u−u3+εβ2βv3−v+uv3−vu3+L2v−v3−γv−v3+λβ2βu−u3+αuv3−vu3=−L1u−u32−L1u−u3εβv−v3+L1u−u3εβ2uv3−vu3−γL2v−v32+L2v−v3λβu−u3+L2v−v3λαβ2uv3−vu3=−L1u−u32−L1εβu−u3v−v3+L1εβ2v3u−u32−L1εβ2u3u−u3v−v3−γL2v−v32+L2λβu−u3v−v3+L2λαβ2v3u−u3v−v3−L2λαβ2u3v−v32=−L1u−u321−εβ2v3−L2v−v32γ+λαβ2u3−u−u3v−v3L1εβ1+u3β−L2λβ1+αv3β2,
ChoosingL1=1 and L2=εβ+u3/λβ+αv3, one can write(25)dLdt≤−u−u321−εβ2v3+v−v32γ+λαβ2u3εβ+u3λβ+αv3.
Thus, if1−ε/β2v3>0, then dL/dt<0 and dL/dt=0 if u,v=u3,v3. Hence, the proof is complete.Example 5.
Numerically, we display the solution of the system (2) with ε=1.0,β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15, and σ=0.1 in Figure 5. It follows from Theorem 4 that the equilibrium point e30.2753798784,0.2830660001 is globally stable due to 1−ε/β2v3>0.Figure 5
Global stability of the system (2) around e30.2753798784,0.2830660001 corresponding to ε=1.0,β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15,σ=0.1.


## 3.2. Bifurcation Analysis of the System

Now, we investigate bifurcation phenomena of the system according to the parameter valueε.Theorem 5.
The system (2) undergoes a saddle-node bifurcation around e3u3,v3 at ε=εSN, where(26)εSN=γv33α3+3α2v33γβ+u3+α2λu3v3+3αv3γβ2+βu3+u32+αu3λu3+β+βγβ2+3βu3+3u32+u33γαγv32+u3+βv3γ−λβ.Proof.
For a saddle-node bifurcation, one of the eigenvalues ofJ3 has to be zero which is possible iff det J3=Γ11Γ22−Γ12Γ21=0, which gives ε=εSN. The other eigenvalue is obtained from the relation tr J3=Γ11+Γ22 evaluated at ε=εSN and should be negative [28].
Let℧ and Υ be the eigenvectors corresponding to eigenvalue 0 of the matrix J3 and its transpose, respectively. We obtain that ℧=θ1,θ2T,Υ=ℏ1,ℏ2T where(27)θ1=−Γ12θ2Γ11=−Γ22θ2Γ21,ℏ1=−Γ21ℏ2Γ11=−Γ22ℏ2Γ12.
Andθ2,ℏ2 are any two nonzero real numbers.
Now, for the system (2), we can get(28)ΥTFεX3,εSN=ℏ1ℏ2∂F1∂ε∂F2∂ε=ℏ1ℏ2−uvβ+u+αv0=−ℏ1uvβ+u+αv≠0ΥTD2FX3,εSN℧,℧=ℏ1ℏ2∂2F1∂u2θ12+∂2F1∂u∂vθ1θ2+∂2F1∂v∂uθ2θ1+∂2F1∂v2θ22∂2F2∂u2θ12+∂2F2∂u∂vθ1θ2+∂2F2∂v∂uθ2θ1+∂2F2∂v2θ22=ℏ1∂2F1∂u2θ12+∂2F1∂u∂vθ1θ2+∂2F1∂v∂uθ2θ1+∂2F1∂v2θ22+ℏ2∂2F2∂u2θ12+∂2F2∂u∂vθ1θ2+∂2F2∂v∂uθ2θ1+∂2F2∂v2θ22=ℏ1∂∂u1−2u−δ−εvβ+αvβ+u+αv2θ12−∂∂uεuβ+uβ+u+αv2θ1θ2+∂∂v1−2u−δ−εvβ+αvβ+u+αv2θ2θ1−∂∂vεuβ+uβ+u+αv2θ22+ℏ2∂∂uλvβ+αvβ+u+αv2θ12+∂∂uγ1−v−γv+λuβ+uβ+u+αv2−σθ1θ2+∂∂vλvβ+αvβ+u+αv2θ2θ1+∂∂vγ1−v−γv+λuβ+uβ+u+αv2−σθ22=ℏ1−2+2εvβ+αvβ+u+αv3θ12−2εβ2+βu+βαv+2uαvβ+u+αv3θ1θ2+2εαuβ+uβ+u+αv3θ22−ℏ22λvβ+αvβ+u+αv3θ12+λβ2+βu+βαv+2uvαβ+u+αv3θ1θ2−2γ+2uλαβ+uβ+u+αv3θ22=−2θ12ℏ1+2γθ22ℏ2+1β+u+αv3ℏ12εvβ+vαθ12−2εβ2+uβ+vβα+2uvαθ1θ2+2uεαβ+uθ22−ℏ22vλβ+vαθ22−2λβ2+uβ+vβα+2uvαθ1θ2−2uλαβ+uθ22≠0.
Here, we writeu,v instead of u3,v3 for algebraic convenience.
So, the system experiences a saddle-node bifurcation [28] around the positive interior equilibrium e3u3,v3 at ε=εSN. Again, since ΥTFεX2,εSN≠0, therefore, the system does not attain neither transcritical bifurcation nor pitch-fork bifurcation.Theorem 6.
The system (2) admits a Hopf bifurcation around e3u3,v3 at ε=εHB, where ε=εHB=λα+u3+γv3β+u3+αv32/u3v3.Proof.
In a system, a Hopf bifurcation typically occurs when a complex conjugate pair of eigenvalues of the linearized flow at the fixed pointe3u3,v3 becomes purely imaginary, i.e., if det J3>0 and if tr J3=0 for some value of ε say εHB.
The characteristic equation is given by(29)μ2−Γ11+Γ22μ+Γ11Γ22−Γ12Γ21=0,where u−u3≈eμt and v−v3≈eμt.
If trJ3=Γ11+Γ22=0 at ε=εHB, then both the eigenvalues will be purely imaginary if det J3=Γ11Γ22−Γ12Γ21 is positive.
Replacingμ=μ1+iμ2 into the corresponding characteristic equation and separating real and imaginary parts, we get(30)μ12−μ22−trJ3μ1+detJ3=0,(31)2μ1μ2−trJ3μ2=0.
Elementary differential of (30) with respect to ε and considering μ1=0, we get(32)dμ1dεε=εHB=v3u32β+u3+αv32≠0,
Hence, the system goes through a Hopf bifurcation arounde3u3,v3 at ε=εHB, where(33)εHB=λα+u3+γv3β+u3+αv32u3v3.Example 6.
It follows from Theorem6 that the system could have a limit cycle on account of Hopf bifurcation. In fact, Figure 6 shows the existence of a limit cycle numerically when the parameters are taken as ε=1.0596222,β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15, and σ=0.1 in the system (2).Figure 6
A limit cycle of the system (2) around e30.2371368969,0.2503112874 corresponding to ε=1.0596222,β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15,σ=0.1.Remark 1.
For better understanding of qualitative change of the temporal dynamics of the system (2), one may observe the bifurcation diagrams of both prey and predator species with respect to the system parameter ϵ (cf. Figures 7 and 8). In support of this technique of bifurcation offered in Figure 7, the successive maxima of u and v in the ranges 0.01,0.8 and 0.01,0.33, respectively, as a function of ϵ in the range 0.95≤ϵ≤1.2 and the other system parameters are captioned in the figure. The other bifurcation diagram is generated for the successive maxima of population sizes of u and v in the ranges 0.01,0.9 and 0.0,0.7, respectively, as a function of ϵ in the range 0.8≤ϵ≤1.5 and the other system parameters are captioned in Figure 8. It is interesting to note from Figures 7 and 8 that maximum concentration of predator species reduces drastically in the presence of prey harvesting factor δ at positive level. Due to the certain amount of prey harvesting from the interacting species system, predator species cannot consume prey species in earlier manner and accordingly this scenario supports the natural phenomena very well.Figure 7
The bifurcation of the system (2) around e3u3,v3 corresponding to the bifurcating parameter ε, where other parameter values are β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15,σ=0.1.Figure 8
The bifurcation of the system (2) around e3u3,v3 corresponding to the bifurcating parameter ε, where other parameter values are β=0.031,α=0.871881,δ=0.0,γ=0.0327,λ=0.15,σ=0.05.


## 4. Pattern Formation Analysis of the System

In 1952, Alan Turing recommended in his seminal work on “The Chemical Basis of Morphogenesis” that the formation of biological patterns can be understood through reaction–diffusion framework, that is, in a form of coupled system of nonlinear parabolic equations [1]. The novelty of Turing’s 1952 paper has influenced many evolutionary biologists, theoretical ecologists, and mathematicians. Under this development, chemicals produced by cells interact like activators or inhibitors and diffuse through the medium at diverse rates. This can generate a symmetry breaking of the uniform concentrations, a mechanism that is frequently called a “diffusion-driven instability” [29]. Turing recommended an ingenious mathematical theory to produce spatial patterns in a chemical system. For this, he reflected on a system of morphogens reacting and diffusing in such an approach that, without diffusion, they revealed a spatially uniform steady state which would be stable. Then, he explained that the presence of diffusion could lead to diffusion-driven instability resulting in a spatially heterogeneous pattern of chemical concentrations. A reaction–diffusion interacting species model equation with harvesting of such dynamics has been offered herein to study biological pattern formation. Consequently, the system (2) develops into a reaction–diffusion predator–prey system:(34a)∂u∂t=u1−u−εuvβ+u+αv−δu+d1∇2u,(34b)∂v∂t=γv1−v+λuvβ+u+αv−σv+d2∇2v,(34c)u0,x,y≥0,v0,x,y≥0,where u0,x,y and v0,x,y stand for the prey and predator species sizes at any time t and at the spatial location x,y in 2D spatial domain; ∇2≡∂2/∂x2+∂2/∂y2 stands for the common Laplacian operator in 2D spatial domain Ω; and the remaining system parameters have the same meaning as those mentioned in case of system (2). One may believe that the environment is uniform; that is, the system parameters are independent of space or time.Mathematically, the offered reaction–diffusion model in this section is coupled nonlinear partial differential equations (PDEs). The parametersd2 and d1 are the positive self-diffusion coefficients of the v and u, respectively. As a general rule, self-diffusion is assumed as a spatial transmission way, which shifts from higher concentration on the way to lower concentration. Note that since the reaction–diffusion model (34a)–(34c) takes into account possible existence of alternative food sources, generalist predator can survive without prey and extinction of prey does not lead to extinction of predator species. Also, one may believe the subsequent zero-flux boundary conditions:(35)∂u∂n=∂v∂n=0,x,y∈∂Ω,where n is the outward unit normal vector of the smooth boundary ∂Ω. The zero-flux boundary conditions signify that the system is self-contained with zero-species flux across the boundary.Previously, we have observed that the unique and positive coexistence equilibrium pointe3u3,v3 exists under some parametric restrictions. Moreover, one may notice that the coexistence equilibrium point e3u3,v3 for the nonspatial model of (34a)–(34c) is a spatially homogeneous steady state for the reaction–diffusion system (34a)–(34c). With definite parametric restrictions, such a steady state could be linearly stable without the presence of diffusion factor but linearly unstable with diffusion factor; this is the well-known recognizable information of diffusion-driven instability [1].Linear stability analysis about the homogeneous steady statee3u3,v3 for the reaction–diffusion system (34a)–(34c) and the essential conditions for the onset of diffusion-driven instability are well known. On a square 2D domain with zero-flux boundary conditions, diffusion-driven instability needs the following conditions (see [3, 9, 10, 13] for detailed information):(36)iΓ11+Γ22<0,iiΓ11Γ22−Γ12Γ21>0,iiid1Γ22+d2Γ11>0,ivd1Γ22+d2Γ1124d1d2>Γ11Γ22−Γ12Γ21,⟹d1Γ22+d2Γ11>2d1d2Γ11Γ22−Γ12Γ21,where the Jacobian matrix J of the nonspatial system of (34a)–(34c) is given by(37)J=1−2u3−δ−εv3β+αv3β+u3+αv32−εu3β+u3β+u3+αv32λv3β+αv3β+u3+αv32γ1−v3−γv3−σ+λu3β+u3β+u3+αv32=Γij2x2i,j=1,2,Out of four conditions, the first two conditions of diffusion-driven instability are well-known Routh–Hurwitz criterion for which the coexistence equilibrium pointe3 related to the nonspatial system of (34a)–(34c) is linearly stable. One may interestingly observe that the last condition iv summarizes condition iii as well. Moreover, one may examine that the three restrictions i,ii,iv are essential conditions for spatial pattern formation given the proper perturbation and system length scale [30, 31]. One may define “Turing’s space” as the parametric space which satisfies the four conditions i−iv of diffusion-driven instability.Figure9 shows the diffusion-driven instability region of the diffusive system (34a)–(34c) in d1d2-plane related to α=0.91881, β=0.01, γ=0.0327, δ=0.05, ε=0.9, λ=0.5, and σ=0.05, and we have captured contour pictures of Turing’s patterns in this region.Figure 9
Turing’s space of the reaction–diffusion model (34a)–(34c) in d1,d2 parametric space related to the parameter values α=0.91881, β=0.01, γ=0.0327, δ=0.05, ε=0.9, λ=0.5, and σ=0.05.


### 4.1. The Spatial Patterns through Diffusion-Driven Instability of the System (34a)–(34c)

This section deals with the presentation of extensive numerical simulations of the reaction–diffusion predator–prey system (34a)–(34c) with the effect of harvesting and shows qualitative consequences in two-dimensional space. The numerical integration of the system (34a)–(34c) is executed through explicit Euler method with time step size 0.01; space step size 1.0; system size 200×200, and the standard five-point approximation for the two-dimensional Laplacian with zero-flux boundary conditions. The initial conditions are always a small amplitude random perturbation about e3.Different kinds of dynamics are examined through numerical simulation and we have figured out that the distributions of predator and prey are always of the same type. For this reason, one may confine their interest to the pattern formation of prey species or predator species only where the parameter values are in the region of Turing’s space. But we have incorporated spatial patterns for predator species to understand the fluctuation of predator’s concentration properly.Figure10 demonstrates the development of spatial pattern formation for the reaction–diffusion system (34a)–(34c) corresponding to the system parameter values, captioned in the connected figure. In this case, the pattern takes a long time to settle down, starting with a homogeneous state e3u3,v3 and the random perturbation leads to the formation of cold spots and typical circular patterns mixture (cf. Figure 10(c)) from cold spots patterns (cf. Figure 10(a)) at stationary level. Also, from Figure 11, one may identify predator spatial patterns sequence simultaneously. From the colour bar of each panel, one may notice that due to the increasing value of prey harvesting parameter δ, the size/concentration of prey species is also increasing but in the other way generalist predator’s size/concentration is subsequently decreasing. Biologically, one may conclude that optimal harvesting policy is essential for both species to survive in the system.Figure 10
Numerical integration of the reaction–diffusion system (34a)–(34c) with Neumann boundaries; to generate stationary Turing’s patterns by prey species u at (a) δ=0.0; (b) δ=0.01; and (c) δ=0.023. Other parameter values are α=0.91881, β=0.005, γ=0.0327, ε=1.0, λ=0.15, σ=0.05, d1=0.1, and d2=20.0.
(a)(b)(c)Figure 11
Numerical integration of the reaction–diffusion system (34a)–(34c) with Neumann boundaries; to generate stationary Turing’s patterns by predator species v at (a) δ=0.0; (b) δ=0.01; and (c) δ=0.023. Other parameter values are α=0.91881, β=0.005, γ=0.0327, ε=1.0, λ=0.15, σ=0.05, d1=0.1, and d2=20.0.
(a)(b)(c)Whenδ=0.184, we illustrate the cold spots patterns of the prey species in the spatial domain at t=20000 in Figure 12(a) corresponding to α=0.91881, β=0.0041, γ=0.0327, ε=1.0, λ=0.164, σ=0.1, d1=0.1, and d2=20.0. Although the qualitative feature of all the figures (12(a) for δ=0.184, 12(b) for δ=0.192, and 12(c) for δ=0.2) are similar, there is a fundamental difference for the spatially extended mode. At last, we find out that the cold spots patterns with smaller radius prevail over the whole spatial domain and the system dynamics does not undergo any further changes. From Figure 13, one may recognize predator spatial patterns sequence, in a similar manner.Figure 12
Numerical integration of the reaction–diffusion system (34a)–(34c) with Neumann boundaries; to generate stationary Turing’s patterns by prey species u at (a) δ=0.184; (b) δ=0.192; and (c) δ=0.2. Other parameter values are α=0.91881, β=0.0041, γ=0.0327, ε=1.0, λ=0.164, σ=0.1, d1=0.1, and d2=20.0.
(a)(b)(c)Figure 13
Numerical integration of the reaction–diffusion system (34a)–(34c) with Neumann boundaries; to generate stationary Turing’s patterns by predator species v at (a) δ=0.184; (b) δ=0.192; and (c) δ=0.2. Other parameter values are α=0.91881, β=0.0041, γ=0.0327, ε=1.0, λ=0.164, σ=0.1, d1=0.1, and d2=20.0.
(a)(b)(c)On increasing the value of significant parameterδ, one may observe the following prey spatial patterns sequence from Figure 14: “cold spots (cf. Figure 14(a)) ⟶ labyrinthine (cf. Figure 14(b)) ⟶ stripes (cf. Figure 14(c)) ⟶ hot spots-stripes mixtures (cf. Figure 14(d)) ⟶ hot spots-short stripes mixtures (cf. Figure 14(e)) ⟶ hot spots (cf. Figure 14(f)).” Moreover, from Figure 15, one may detect predator spatial patterns sequence simultaneously. It is interesting to note from the colour bar of each panel that due to the increasing value of prey harvesting parameter δ from 5% to 9%, the size/concentration of prey species is also increasing but in the other approach generalist predator’s size/concentration is subsequently decreasing. Biologically, one may conclude that optimal harvesting policy is essential for both species to survive in the system.Figure 14
Numerical integration of the reaction–diffusion system (34a)–(34c) with Neumann boundaries; to generate stationary Turing’s patterns by prey species u at (a) δ=0.05; (b) δ=0.053; (c) δ=0.055; (d) δ=0.057; (e) δ=0.06; and (f) δ=0.09. Other parameter values are α=0.91881, β=0.01, γ=0.0327, ε=0.9, λ=0.5, σ=0.05, d1=0.1, and d2=20.0.
(a)(b)(c)(d)(e)(f)Figure 15
Numerical integration of the reaction–diffusion system (34a)–(34c) with Neumann boundaries; to generate stationary Turing’s patterns by predator species v at (a) δ=0.05; (b) δ=0.053; (c) δ=0.055; (d) δ=0.057; (e) δ=0.06; and (f) δ=0.09. Other parameter values are α=0.91881, β=0.01, γ=0.0327, ε=0.9, λ=0.5, σ=0.05, d1=0.1, and d2=20.0.
(a)(b)(c)(d)(e)(f)However, if we fix the harvesting coefficient for prey species (δ) at some positive level, then for increasing values of harvesting coefficient for predator species σ at positive level, the numerical simulation develops to a mixed-mode steady state whose graphical depiction is not presented here for the sake of brevity.


## 4.1. The Spatial Patterns through Diffusion-Driven Instability of the System (34a)–(34c)

This section deals with the presentation of extensive numerical simulations of the reaction–diffusion predator–prey system (34a)–(34c) with the effect of harvesting and shows qualitative consequences in two-dimensional space. The numerical integration of the system (34a)–(34c) is executed through explicit Euler method with time step size 0.01; space step size 1.0; system size 200×200, and the standard five-point approximation for the two-dimensional Laplacian with zero-flux boundary conditions. The initial conditions are always a small amplitude random perturbation about e3.Different kinds of dynamics are examined through numerical simulation and we have figured out that the distributions of predator and prey are always of the same type. For this reason, one may confine their interest to the pattern formation of prey species or predator species only where the parameter values are in the region of Turing’s space. But we have incorporated spatial patterns for predator species to understand the fluctuation of predator’s concentration properly.Figure10 demonstrates the development of spatial pattern formation for the reaction–diffusion system (34a)–(34c) corresponding to the system parameter values, captioned in the connected figure. In this case, the pattern takes a long time to settle down, starting with a homogeneous state e3u3,v3 and the random perturbation leads to the formation of cold spots and typical circular patterns mixture (cf. Figure 10(c)) from cold spots patterns (cf. Figure 10(a)) at stationary level. Also, from Figure 11, one may identify predator spatial patterns sequence simultaneously. From the colour bar of each panel, one may notice that due to the increasing value of prey harvesting parameter δ, the size/concentration of prey species is also increasing but in the other way generalist predator’s size/concentration is subsequently decreasing. Biologically, one may conclude that optimal harvesting policy is essential for both species to survive in the system.Figure 10
Numerical integration of the reaction–diffusion system (34a)–(34c) with Neumann boundaries; to generate stationary Turing’s patterns by prey species u at (a) δ=0.0; (b) δ=0.01; and (c) δ=0.023. Other parameter values are α=0.91881, β=0.005, γ=0.0327, ε=1.0, λ=0.15, σ=0.05, d1=0.1, and d2=20.0.
(a)(b)(c)Figure 11
Numerical integration of the reaction–diffusion system (34a)–(34c) with Neumann boundaries; to generate stationary Turing’s patterns by predator species v at (a) δ=0.0; (b) δ=0.01; and (c) δ=0.023. Other parameter values are α=0.91881, β=0.005, γ=0.0327, ε=1.0, λ=0.15, σ=0.05, d1=0.1, and d2=20.0.
(a)(b)(c)Whenδ=0.184, we illustrate the cold spots patterns of the prey species in the spatial domain at t=20000 in Figure 12(a) corresponding to α=0.91881, β=0.0041, γ=0.0327, ε=1.0, λ=0.164, σ=0.1, d1=0.1, and d2=20.0. Although the qualitative feature of all the figures (12(a) for δ=0.184, 12(b) for δ=0.192, and 12(c) for δ=0.2) are similar, there is a fundamental difference for the spatially extended mode. At last, we find out that the cold spots patterns with smaller radius prevail over the whole spatial domain and the system dynamics does not undergo any further changes. From Figure 13, one may recognize predator spatial patterns sequence, in a similar manner.Figure 12
Numerical integration of the reaction–diffusion system (34a)–(34c) with Neumann boundaries; to generate stationary Turing’s patterns by prey species u at (a) δ=0.184; (b) δ=0.192; and (c) δ=0.2. Other parameter values are α=0.91881, β=0.0041, γ=0.0327, ε=1.0, λ=0.164, σ=0.1, d1=0.1, and d2=20.0.
(a)(b)(c)Figure 13
Numerical integration of the reaction–diffusion system (34a)–(34c) with Neumann boundaries; to generate stationary Turing’s patterns by predator species v at (a) δ=0.184; (b) δ=0.192; and (c) δ=0.2. Other parameter values are α=0.91881, β=0.0041, γ=0.0327, ε=1.0, λ=0.164, σ=0.1, d1=0.1, and d2=20.0.
(a)(b)(c)On increasing the value of significant parameterδ, one may observe the following prey spatial patterns sequence from Figure 14: “cold spots (cf. Figure 14(a)) ⟶ labyrinthine (cf. Figure 14(b)) ⟶ stripes (cf. Figure 14(c)) ⟶ hot spots-stripes mixtures (cf. Figure 14(d)) ⟶ hot spots-short stripes mixtures (cf. Figure 14(e)) ⟶ hot spots (cf. Figure 14(f)).” Moreover, from Figure 15, one may detect predator spatial patterns sequence simultaneously. It is interesting to note from the colour bar of each panel that due to the increasing value of prey harvesting parameter δ from 5% to 9%, the size/concentration of prey species is also increasing but in the other approach generalist predator’s size/concentration is subsequently decreasing. Biologically, one may conclude that optimal harvesting policy is essential for both species to survive in the system.Figure 14
Numerical integration of the reaction–diffusion system (34a)–(34c) with Neumann boundaries; to generate stationary Turing’s patterns by prey species u at (a) δ=0.05; (b) δ=0.053; (c) δ=0.055; (d) δ=0.057; (e) δ=0.06; and (f) δ=0.09. Other parameter values are α=0.91881, β=0.01, γ=0.0327, ε=0.9, λ=0.5, σ=0.05, d1=0.1, and d2=20.0.
(a)(b)(c)(d)(e)(f)Figure 15
Numerical integration of the reaction–diffusion system (34a)–(34c) with Neumann boundaries; to generate stationary Turing’s patterns by predator species v at (a) δ=0.05; (b) δ=0.053; (c) δ=0.055; (d) δ=0.057; (e) δ=0.06; and (f) δ=0.09. Other parameter values are α=0.91881, β=0.01, γ=0.0327, ε=0.9, λ=0.5, σ=0.05, d1=0.1, and d2=20.0.
(a)(b)(c)(d)(e)(f)However, if we fix the harvesting coefficient for prey species (δ) at some positive level, then for increasing values of harvesting coefficient for predator species σ at positive level, the numerical simulation develops to a mixed-mode steady state whose graphical depiction is not presented here for the sake of brevity.


## 5. Discussion

In this section, we discuss optimal harvesting problems and bionomic equilibrium of the system (2). In fact, until now, although there are many researches dealt with predator–prey system with harvesting (cf. [10, 21, 32–37]), it is not prevalent to study optimal harvesting strategy and economic cost theoretically.


### 5.1. Bionomic Equilibrium

The term is derived from the concept concerning with both biological and economic equilibrium and is given bydu/dt=dv/dt=0. The economic equilibrium is said to be achieved when the TR (total revenue collected by selling the harvested species) equals TC (the total cost for the effort devoted for harvesting).The revenue at any timet (economic rent) is given by(38)N=TR−TC=p1u−q1δ+p2v−q2σ,where p1  = price per unit biomass of the prey species u, p2  = price per unit biomass of the predator species v, q1  = fishing cost per unit effort of the prey species u, and q2  = fishing cost per unit effort of the predator species v.So, the bionomic equilibrium is given by the simultaneous equations:(39)u˙=u1−u−εuvβ+u+αv−δu,(40)v˙=γv1−v+λuvβ+u+αv−σv,(41)N=p1u−q1δ+p2v−q2σ=0.Since the price and the cost of the predator are not sure, we will consider the following cases in order to determine the bionomic equilibrium:Case 1If(42)q1p1>u,i.e.,(43)p1u−q1<0.holds, i.e., to say that the total cost exceeds the total revenue for the harvesting of prey, obviously the prey harvesting will be stopped (i.e.,δ=0) and the predator harvesting remains operational if p2v−q2>0.Then, from equation (40), we have(44)v1∞=q2p2.Substituting it in equation (39), we get(45)u1−u−εvβ+u+αv=0,⟹u=1−εvβ+u+αv,⟹u=121−β−αv±121+2β+2αv+β+αv2−4εv,⟹u1∞=121−β−αv1∞±121+2β+2αv1∞+β+αv1∞2−4εv1∞.Again, substituting the value ofu1∞ and v1∞ in equation (40), we get(46)γ1−v+λuβ+u+αv−σ=0,⟹σ=γ1−v+λuβ+u+αv,⟹σ1∞=γ1−v1∞+λu1∞β+u1∞+αv1∞.Hence, the following conditions need to be satisfied in order to have the bionomic equilibriumu1∞,v1∞,0,σ1∞:(47)ip2>q2,iiγ>γq2p2.Case 2If(48)q2p2<v,q1p1<u,i.e.,(49)p2v−q2>0,p1u−q1>0,holds; in this case, the total revenue exceeds the total cost for two populations and the harvesting is in operational because it can bring profit for capturing species.From (41), we have(50)u1∞=q1p1,v1∞=q2p2.Substituting the values from (49) in (39) and (40), it is easy to obtain(51)δ1∞=1−u1∞+εq2p2β+u1∞+αv1∞,σ1∞=γ1−v1∞+λq1p1β+u1∞+αv1∞.So, to have the bionomic equilibrium, i.e.,u,v,δ,σ⟶u1∞,v1∞,δ1∞,σ1∞, the following conditions must be satisfied:(52)i1>u1∞+εq2p2β+u1∞+αv1∞,ii1>v1∞,iiiv1∞<γ+λq1p1β+u1∞+αv1∞.It is very much clear that when the intrinsic growth rates of the two species exceed some value, then only the existence of bionomic equilibrium is established.Example 7.
For Case 1, we display trajectories of the solutions of the system (2) with respect to the values σ and δ in Figures 16 and 17, respectively, when ε=1.0,β=0.031,δ=0.0,α=0.871881,γ=0.0327,λ=0.15 and for Case 2, in Figures 18 and 19, respectively, when ε=1.0,β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15.Figure 16
Time series plot of the system (2) for Case 1 with different values of σ corresponding to the parameter values ε=1.0,β=0.031,δ=0.0,α=0.871881,γ=0.0327,λ=0.15, and u0,v0=0.266,0.278.Figure 17
Time series plot of the system (2) for Case 1 with different values of δ corresponding to the parameter values ε=1.0,β=0.031,α=0.871881,σ=0.05,γ=0.0327,λ=0.15, and u0,v0=0.266,0.278.Figure 18
Time series plot of the system (2) for Case 2 with different values of σ corresponding to the parameter values ε=1.0,β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15, and u0,v0=0.266,0.278.Figure 19
Time series plot of the system (2) for Case 2 with different values of δ corresponding to the parameter values ε=1.0,β=0.0041,α=0.91881,γ=0.0327,λ=0.15,σ=0.1, andu0,v0=0.266,0.278.


### 5.2. Optimal Harvesting

The one major problem in commercial exploitation of renewable resources is to determine the optical trade-off between present and future harvests. In order to calculate the optimal harvesting policy, we consider the present valueJ of a continuous time-stream of revenues given by(53)J=∫0∞e−ωtp1u−q1δt+p2v−q2σtdt,where ω represents the instantaneous annual rate of discount and δ and σ are the control variables. We intend to maximize J subject to the state equations (2a)–(2c) by introducing Pontryagin’s maximal principle. The Hamiltonian for the problem (52) is given by(54)H¯=e−ωtp1u−q1δt+p2v−q2σt+Ω1u1−u−εuvβ+u+αv−δu+Ω2γv1−v+λuvβ+u+αv−σv,where Ω1t and Ω2t are the adjoint variables.For optimality, the necessary conditions are(55)∂H¯∂δ=0,∂H¯∂σ=0.Then, we have(56)∂H¯∂δ=0,e−ωtp1u−q1−Ω1u=0,Ω1=e−ωtp1−q1u.Similarly,(57)Ω2=e−ωtp2−q2v.Equations (55) and (56) can also be written as(58)Ω1eωt=p1−q1u,Ω2eωt=p2−q2v.Therefore, the shadow priceseωtΩit for i=1,2 do not vary with time in the optimal equilibrium. Hence, they remain bounded as t⟶∞ (i.e., they satisfy the transversality condition at ∞).We intend to derive here an optimal equilibrium solution of the problem. Since we are considering an equilibrium solution,u,v,δ, and σ are to be treated as constants in the subsequent steps.Now,(59)dΩ1dt=−∂H¯∂u⟹−ωe−ωtp1−q1u=−e−ωtp1δ+p1−q1u1−2u−εvβ+αvβ+u+αv2−δ+p2−q2vλvβ+αvβ+u+αv2.Similarly,(60)dΩ2dt=−∂H¯∂v,⟹−ωe−ωtp2−q2v=−e−ωtp2σ−p1−q1uεuβ+uβ+u+αv2+p2−q2vγ1−v−γv+λuβ+uβ+u+αv2−σ,From equations (58) and (59), we may find the positive values of u,v=uω,vω.Substitutinguω,vω into (2), we get the equations as follows:(61)δ1ω=1−u−εvβ+u+αv,σ1ω=γ1−v+λuβ+u+αv.Then, the optimal equilibrium efforts levelsδ1ω and σ1ω can be obtained if the following conditions are satisfied:(62)i1>u+εvβ+u+αv,ii1>v,withγv<γ+λuβ+u+αv.Example 8.
If we consider the set of values whereε=1.0,β=0.041,α=0.91881,γ=0.0327,λ=0.15,p1=0.5,p2=1.0,δ=0.083,σ=0.5 and ω=0.8, respectively, we can get the optimal equilibrium uω,vω=0.1848452592,0.5314707140 and the corresponding optimal harvesting efforts are δ1ω=0.0709708763 and σ1ω=0.05414492630.


## 5.1. Bionomic Equilibrium

The term is derived from the concept concerning with both biological and economic equilibrium and is given bydu/dt=dv/dt=0. The economic equilibrium is said to be achieved when the TR (total revenue collected by selling the harvested species) equals TC (the total cost for the effort devoted for harvesting).The revenue at any timet (economic rent) is given by(38)N=TR−TC=p1u−q1δ+p2v−q2σ,where p1  = price per unit biomass of the prey species u, p2  = price per unit biomass of the predator species v, q1  = fishing cost per unit effort of the prey species u, and q2  = fishing cost per unit effort of the predator species v.So, the bionomic equilibrium is given by the simultaneous equations:(39)u˙=u1−u−εuvβ+u+αv−δu,(40)v˙=γv1−v+λuvβ+u+αv−σv,(41)N=p1u−q1δ+p2v−q2σ=0.Since the price and the cost of the predator are not sure, we will consider the following cases in order to determine the bionomic equilibrium:Case 1If(42)q1p1>u,i.e.,(43)p1u−q1<0.holds, i.e., to say that the total cost exceeds the total revenue for the harvesting of prey, obviously the prey harvesting will be stopped (i.e.,δ=0) and the predator harvesting remains operational if p2v−q2>0.Then, from equation (40), we have(44)v1∞=q2p2.Substituting it in equation (39), we get(45)u1−u−εvβ+u+αv=0,⟹u=1−εvβ+u+αv,⟹u=121−β−αv±121+2β+2αv+β+αv2−4εv,⟹u1∞=121−β−αv1∞±121+2β+2αv1∞+β+αv1∞2−4εv1∞.Again, substituting the value ofu1∞ and v1∞ in equation (40), we get(46)γ1−v+λuβ+u+αv−σ=0,⟹σ=γ1−v+λuβ+u+αv,⟹σ1∞=γ1−v1∞+λu1∞β+u1∞+αv1∞.Hence, the following conditions need to be satisfied in order to have the bionomic equilibriumu1∞,v1∞,0,σ1∞:(47)ip2>q2,iiγ>γq2p2.Case 2If(48)q2p2<v,q1p1<u,i.e.,(49)p2v−q2>0,p1u−q1>0,holds; in this case, the total revenue exceeds the total cost for two populations and the harvesting is in operational because it can bring profit for capturing species.From (41), we have(50)u1∞=q1p1,v1∞=q2p2.Substituting the values from (49) in (39) and (40), it is easy to obtain(51)δ1∞=1−u1∞+εq2p2β+u1∞+αv1∞,σ1∞=γ1−v1∞+λq1p1β+u1∞+αv1∞.So, to have the bionomic equilibrium, i.e.,u,v,δ,σ⟶u1∞,v1∞,δ1∞,σ1∞, the following conditions must be satisfied:(52)i1>u1∞+εq2p2β+u1∞+αv1∞,ii1>v1∞,iiiv1∞<γ+λq1p1β+u1∞+αv1∞.It is very much clear that when the intrinsic growth rates of the two species exceed some value, then only the existence of bionomic equilibrium is established.Example 7.
For Case 1, we display trajectories of the solutions of the system (2) with respect to the values σ and δ in Figures 16 and 17, respectively, when ε=1.0,β=0.031,δ=0.0,α=0.871881,γ=0.0327,λ=0.15 and for Case 2, in Figures 18 and 19, respectively, when ε=1.0,β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15.Figure 16
Time series plot of the system (2) for Case 1 with different values of σ corresponding to the parameter values ε=1.0,β=0.031,δ=0.0,α=0.871881,γ=0.0327,λ=0.15, and u0,v0=0.266,0.278.Figure 17
Time series plot of the system (2) for Case 1 with different values of δ corresponding to the parameter values ε=1.0,β=0.031,α=0.871881,σ=0.05,γ=0.0327,λ=0.15, and u0,v0=0.266,0.278.Figure 18
Time series plot of the system (2) for Case 2 with different values of σ corresponding to the parameter values ε=1.0,β=0.0041,α=0.91881,δ=0.2,γ=0.0327,λ=0.15, and u0,v0=0.266,0.278.Figure 19
Time series plot of the system (2) for Case 2 with different values of δ corresponding to the parameter values ε=1.0,β=0.0041,α=0.91881,γ=0.0327,λ=0.15,σ=0.1, andu0,v0=0.266,0.278.


## 5.2. Optimal Harvesting

The one major problem in commercial exploitation of renewable resources is to determine the optical trade-off between present and future harvests. In order to calculate the optimal harvesting policy, we consider the present valueJ of a continuous time-stream of revenues given by(53)J=∫0∞e−ωtp1u−q1δt+p2v−q2σtdt,where ω represents the instantaneous annual rate of discount and δ and σ are the control variables. We intend to maximize J subject to the state equations (2a)–(2c) by introducing Pontryagin’s maximal principle. The Hamiltonian for the problem (52) is given by(54)H¯=e−ωtp1u−q1δt+p2v−q2σt+Ω1u1−u−εuvβ+u+αv−δu+Ω2γv1−v+λuvβ+u+αv−σv,where Ω1t and Ω2t are the adjoint variables.For optimality, the necessary conditions are(55)∂H¯∂δ=0,∂H¯∂σ=0.Then, we have(56)∂H¯∂δ=0,e−ωtp1u−q1−Ω1u=0,Ω1=e−ωtp1−q1u.Similarly,(57)Ω2=e−ωtp2−q2v.Equations (55) and (56) can also be written as(58)Ω1eωt=p1−q1u,Ω2eωt=p2−q2v.Therefore, the shadow priceseωtΩit for i=1,2 do not vary with time in the optimal equilibrium. Hence, they remain bounded as t⟶∞ (i.e., they satisfy the transversality condition at ∞).We intend to derive here an optimal equilibrium solution of the problem. Since we are considering an equilibrium solution,u,v,δ, and σ are to be treated as constants in the subsequent steps.Now,(59)dΩ1dt=−∂H¯∂u⟹−ωe−ωtp1−q1u=−e−ωtp1δ+p1−q1u1−2u−εvβ+αvβ+u+αv2−δ+p2−q2vλvβ+αvβ+u+αv2.Similarly,(60)dΩ2dt=−∂H¯∂v,⟹−ωe−ωtp2−q2v=−e−ωtp2σ−p1−q1uεuβ+uβ+u+αv2+p2−q2vγ1−v−γv+λuβ+uβ+u+αv2−σ,From equations (58) and (59), we may find the positive values of u,v=uω,vω.Substitutinguω,vω into (2), we get the equations as follows:(61)δ1ω=1−u−εvβ+u+αv,σ1ω=γ1−v+λuβ+u+αv.Then, the optimal equilibrium efforts levelsδ1ω and σ1ω can be obtained if the following conditions are satisfied:(62)i1>u+εvβ+u+αv,ii1>v,withγv<γ+λuβ+u+αv.Example 8.
If we consider the set of values whereε=1.0,β=0.041,α=0.91881,γ=0.0327,λ=0.15,p1=0.5,p2=1.0,δ=0.083,σ=0.5 and ω=0.8, respectively, we can get the optimal equilibrium uω,vω=0.1848452592,0.5314707140 and the corresponding optimal harvesting efforts are δ1ω=0.0709708763 and σ1ω=0.05414492630.


## 6. Conclusion

A large number of interacting species models have been suggested and were explored for a decade to recognize the underlying mechanisms that are happening in real-world ecological systems. In this paper, the dynamical behaviour of a generalist predator–prey model with linear harvesting of each species and Beddington–DeAngelis functional response is totally analyzed. At first, we nondimensionalize the main model, so that it reduces the number of ecological parameters and then the new model (2) experiences for its behavioural studies. In the course of stability theory of ordinary differential equations (ODEs), it has been verified that the interior equilibrium exists under certain conditions and it is globally asymptotically stable via suitable Lyapunov function structure, which signifies that appropriate harvesting has no impact on the persistence property of the harvesting system.The system under consideration (2) undergoes Hopf and saddle-node bifurcations, with ε being the bifurcating parameter for both the bifurcations. Subsequently, complete analysis regarding the impact of harvesting is carried out and a conclusion is made that, under some appropriate constraints, harvesting has enormous impact on the final size of the interacting species. Also, the predator harvesting should be greater than prey harvesting (i.e., σ>δ); if this condition is not fulfilled, then at 0,0 the prey species u will meet extinction. Then, for the practical significance, we consider the economic profit of the harvesting. The bionomic equilibrium and optimal harvesting policy are studied. Using Pontryagin’s maximum principle, the optimal harvesting policy is discussed. The results show that the optimal harvesting policy may exist and an example is given to show that the optimal harvesting policy of the system (2) is realizable. Finally, our analysis shows that harvesting effort in a generalist predator–prey model plays an imperative role for ecological conservation of interacting species.A variety of spatiotemporal patterns produced by a Beddington–DeAngelis predator–prey system with reaction–diffusion platform are well studied theoretically and numerically (cf. Figures10–15). The spatial parametric space is subdivided into stable and unstable regions particularly through Turing’s bifurcation curve, as revealed in Figure 9. To generate spatiotemporal patterns, many researchers pointed out that spatiotemporal dynamics of a reaction–diffusion system depend on the choice of initial conditions. Throughout the manuscript, we use a particular type of initial condition for discussing the evolution process of the spatial pattern formation of system (34a)–(34c) [38, 39]. The current investigation points out that the variations of harvesting can additionally show a variety of spatiotemporal patterns for certain range of harvesting of both species. Based on this information, we have concluded that two interacting species reaction–diffusion system with the effect of harvesting can have complex dynamical phenomena containing Turing’s patterns scenario.



---

*Source: 7596394-2020-10-28.xml*