On a diffusive predator-prey model with nonlinear harvesting

In this paper, we study the dynamics of a diffusive Leslie-Gower model with a nonlinear harvesting term on the prey. We analyze the existence of positive equilibria and their dynamical behaviors. In particular, we consider the model with a weak harvesting term and find the conditions for the local and global asymptotic stability of the interior equilibrium. The global stability is established by considering a proper Lyapunov function. In contrast, the model with strong harvesting term has two interior equilibria and bi-stability may occur for this system. We also give the conditions of Turing instability and perform a series of numerical simulations and find that the model exhibits complex patterns.


Introduction
Reaction diffusion systems have long been of the interest to both mathematicians and ecologists to understand the dynamics of biological populations.Predator-prey models are arguably the building blocks of ecosystems due to its universal existence and importance.One of the earliest and also the best known predator and prey model is Lotka-Volterra model.However, there are no upper limits to the rates of increase of both prey H and predator P in this model.Leslie where H and P are the density of prey species and the predator species, respectively.We observe that in (1), the carrying capacity of the predator's environment is proportional to the number of prey.This model does not admit limit cycles.Only until relatively recently has it been shown by A. Korobeninikov [8] that the positive equilibrium is globally asymptotically stable by constructing a Lyapunov function.
An important aspect of predator-prey models are the functional responses of the predator to the prey density, which refers to the change of the density of prey attached per time unit and predator unit as the prey density changes.In general, these functional responses, denoted by p(H), is monotone and continuously differentiable on [0, ∞).The following functional responses are called Holling-type I, II, III and IV, respectively where a,b,r are positive constants.
The dynamics of Holling-type predator-prey models are very interesting and have been extensively studied.For example, in [12], the author studied the dynamic behavior of the following predator-prey system incorporating Holling-type II response where the constants r 1 and r 2 are the birth rate of H and P respectively, k is the prey environment carrying capacity, c measures the conversion rate from prey into the predator births, a is the maximum number of prey that can be consumed by the predator per time unit and b measures the extent to which environment provides protection to prey H.Many interesting behaviors, such as stable limit cycles, bifurcations and global stability of constant equilibrium solutions, have been studied.We refer the readers to [1,3,14].
In [2], the authors proposed the following modified Leslie-Gower model with Holling-type II response scheme: where n 1 and n 2 measures the extent to which the environment protects the prey and predator respectively.Many aspects of this model, including permanence, boundedness and global stability of solutions, have already been investigated in [2,4,18].We also refer the readers to [13] for the model incorporating time delay.
The harvesting of the prey species and(or) the predator species is another important issue from the ecological perspective and the economic perspective.In many cases, the goal is to achieve maximum and sustainable yield of the prey species.Sometimes the harvesting of the predator species is introduced in order to control the populations of the species.These are commonly practiced in wildlife management.Predator-prey systems incorporating harvesting terms have also been studied by many authors.We refer the readers to [9,10,15] and the references within.Recently, Gupta [7] proposed the following modified version of Leslie-Gower model with Michaelis-Menten type prey harvesting: Note that they made the assumption that the environment provides the same protection to both the prey and the predator.By considering the following non-dimensional scheme: We can further reduce the system to the following simpler system: The goal of this paper to study the associated reaction diffusion system: The main purpose is to study the dynamical behaviors of this model.We focus on the stability of the equilibria and pattern formation analysis of the model.In particular, we study the stability of the trivial and axial equilibria, as well as the local and global asymptotic stability of interior equilibrium.Another issue we will address is the nonexistence of non-constant positive equilibrium solutions.The rest of this article is organized as follows: In Section 2, We investigate the persistence property of the system.In Section 3, we study the asymptotic behavior of the non-interior equilibria.In Section 4, we consider the local and global stability of the interior equilibrium for the weak harvesting case.In this case, the interior equilibrium is unique.In Section 5, we give the lower bounds of the the diffusion coefficients for which the model has no nonconstant equilibrium solutions.In Section 5, we carry out some numerical simulations to illustrate the diffusion driven instability for our model.Finally in Section 6, we end our investigation with some concluding remarks.

Persistence property
In this section, we show that the solutions to our model are bounded.We also show that our system is persistent under certain condition.Lemma 1. (See [16]) Assume that u(x, t) is defined by then lim t→∞ u(x, t) = K.
Theorem 1.All solutions (u, v) of ( 6) are nonnegative and defined for t > 0. Furthermore, the nonnegative solution (u, v) satisfies Proof.Since (u, v) = (u(x, t), v(x, t)) is a solution of (6), u satisfies From the standard comparison principle of parabolic equations and Lemma 1, it follows that for any arbitrary > 0, there exists t 1 such that for any It then follows that v satisfies Similarly, by comparison principle, there exists t 2 > t 1 > 0 such that for t > t 2 , we have v(x, t) ≤ a + 1 + 2 .
Definition 1. System ( 6) is said to be persistent if for any nonnegative initial data (u 0 (x), v 0 (x)) with u 0 (x) and v 0 (x) are both not identically equal to 0, there exists positive constants 1 , 2 such that In our next result we give the condition under which the system is persistent.Note that this result implies that prey and predator will coexist, no matter what their diffusion coefficients are.
Proof.From the first equation of ( 6), we have We may now apply this lower bound of u to the second equation of system (6) and we conclude there exists t 4 > t 3 such that for any t > t 4 , Therefore we conclude that lim inf

Local stability property of equilibria
We can easily obtain the following equilibria (1) To investigate the local stability of the equilibrium points, we consider the linearized system at E i .The linearized system takes the following form and Proof.The linearized system at E 0 takes the following form Consider the following associated eigenvalue problem To prove that E 0 is unstable under the hypothesis, we need to show that the largest eigenvalue of this system is positive.Let η be an eigenvalue of this system with eigenfunction (U 1 , U 2 ).If U 1 = 0, then η is an eigenvalue of ) with homogeneous Neumann boundary condition.Therefore, η must be real.Similarly, if U 2 = 0, η is also real.Hence all eigenvalues of (15) are real.Since 1 − h c > 0, the principle eigenvalue λ 1 of is positive and the associated eigenvector Ũ1 > 0. It follows that (U 1 , U 2 ) = ( Ũ1 , 0) solves ( 15) with η = λ 1 , i.e., λ 1 is an eigenvalue of (15).Hence the largest eigenvalue of ( 15) is positive and E 0 is unstable.
Proof.At E 2 = (0, a), the linearized system is As in previous theorem, we obtain the following eigenvalue problem is positive and the corresponding eigenfunction Ũ1 > 0.
If Ũ1 ≡ 0, then η 1 is an eigenvalue of The largest eigenvalue of this problem is obviously −δ which is negative.Therefore, we also have η 1 < 0. Hence in this case, E 2 is stable.
Theorem 5. E 3 and E 4 are both unstable.
Proof.At E 3 and E 4 , the Jacobian matrix takes the following form We shall consider the following eigenvalue problem: Since δ > 0, let λ 1 > 0 be the principle eigenvalue of with the associated eigenfunction Ũ2 > 0.

Local and global stability of interior equilibrium
In this section, we will investigate the dynamics of the spatial model in the case of weak harvesting h < c(1 − k).Theorem 6.If 1 − 2u * − ak a+u * − hc (c+u * ) 2 < 0, then the interior equilibrium (u * , v * ) is locally asymptotically stable.
Proof.The linearized problem around the interior equilibrium is where We then expand the solution U of (26) via where z j ∈ R 2 and φ j (x) is the jth eigenfunction of −∆ on Ω with Neumann boundary condition.Substituting (28) into (26) and equating the coefficients of each φ j , we have where A j = −λ j D + J(E 5,6 ) and λ j is the j − th eigenvalue satisfying 0 = λ 0 < λ 1 < λ 2 < . . . .The eigenvalues η 1,2 of A j are determined by where To guarantee that each A j has two eigenvalues with negative real parts, we shall need tr(A j ) < 0 and det(A j ) > 0.
It is easy to see that A < 0 guarantees both conditions hold.
Theorem 7. The positive equilibrium is globally asymptotical stable if Proof.We consider the following Lyapunov function: where Note that and The above equation can be written as where It is obvious that dt < 0 if the matrix above is positive definite.Since k > 0, we shall need the following conditions: Note that implies Φ(u, v) > 0 and the global stability follows.

Nonexistence of nonconstant positive solutions
In this section we shall derive the lower bounds of the diffusion rates under which the system (6) has no nonconstant positive steady-state solutions, that is the nonexistence of nonconstant positive classical solutions of the following elliptic system For any φ ∈ L 1 (Ω), we denote by φ the average value of φ over Ω, that is, Our main result is stated in the following theorem.
Theorem 8.The system (32) does not admit any nonconstant positive solutions if where µ 1 is the second eigenvalue of the Laplacian in Ω with Neumann boundary conditions.
The proof is based on Poincaré inequality and Cauchy inequality [5].

Proof.
In what follows, we shall apply the Cauchy inequality Let (u, v) be a nonconstant positive solution of (32).Then we have u ≤ 1 and v ≤ a + 1. Multiplying the two equations in (32) by ξ := u − ū and η := v − v, and then integrating over Ω and applying the Cauchy inequality, we obtain Applying Poincaré inequality, we have Hence if then system (32) does not admit any nonconstant positive solutions.

Turing instability and pattern formation
As usual, we introduce small perturbations We consider the following linearized system about (u * , v * ) as follows: where J 11 , J 12 , J 21 and J 22 are defined as in (13).
We know any solution of (35) can be expanded into the following Fourier series so that where x = (ξ, η), and 0 < ξ < Lx, 0 < η < Ly. k = (k n , k m ) = (nπ/Lx, mπ/Ly) are the corresponding wavenumbers.Substituting u nm and v nm into (35), we obtain: where k 2 = k 2 n + k 2 m .A general solution of (38) has the form C 1 e λ 1 t +C 2 e λ 2 t , where the constant C 1 , C 2 are determined by the initial conditions and the exponents λ 1 , λ 2 are the eigenvalues of the following matrix Correspondingly, λ 1 , λ 2 are the solutions of the following characteristic equation: where The roots yield the dispersion relation Diffusive instability occurs when at least one of the following conditions is violated: Clearly ρ 1 < 0 is not violated when J 11 + J 22 < 0. Hence, only violation of the condition ρ 2 < 0 will give rise to diffusion instability.Then reversal of the inequality gives the following where Hence the sufficient condition for instability is Hence we arrive at the conclusion that the equilibrium solution is diffusively unstable if Summarizing the above discussions, we obtain the following theorem.

Numerical simulations
In this section, we numerically solve system (6) in two dimensional space.All our numerical simulations employ the zero-flux boundary conditions with a system size 200 × 200 discretized through x → (x 0 , x 1 , . . ., x n ) and y → (y 0 , y 1 , . . ., y n ), with n = 200.The numerical integration is done by using a forward Euler integration with a rather small time step τ = 0.01 × h 2 where h = 1/4 is the space time step.Then we use the standard five-point approximation for the 2D Laplacian with the zero-flux boundary conditions.The concentration at mesh point (x i , y j ) at the moment (n + 1)τ is denoted by (u n+1 i,j , v n+1 i,j ) is given by ).We perform simulations for the following two sets of parameters: In the numerical simulations, we can observe different types of dynamics and we take snapshots of the distribution of the prey u at different time.
Fig 1 shows the evolution of patterns of prey started with small amplitude random perturbation around the stationary solutions (u * , v * ) = (3.01174,3.11174).In this case, one can see that the random initial distribution leads to irregular patterns in the domain as shown at t = 10 and t = 30.Eventually, the spotted pattern forms and the dynamics of the model does not undergo any further changes.We remark that the changes at the very beginning is rather rapid.

Conclusions
In this paper, we investigated the complex dynamics in a diffusive Leslie-Gower Model with a nonlinear harvesting term on the prey.We gave the conditions under which the system is persistent.We also analyzed the stability of the trivial equilibrium and the axial equilibria.We note that the model behaves quite differently with a strong or weak harvesting term.In the case of weak harvesting term, we showed the local and global stability of the unique interior equilibrium under certain conditions.We further analyzed the conditions under which there is no nonconstant positive steady-state solutions.
Furthermore, we also studied the conditions under which the model undergoes diffusion driven instability.We derived the conditions of Turing instability in terms of our parameters analytically.Numerical simulations show that our model can exhibit different patterns such as stripes and spots patterns.

2 PH
addressed this issue and proposed the now well-known Leslie-Gower model dH dt = (r 1 − a 1 P − b 1 H)H, dP dt = (r 2 − a )P,