Stability and bifurcation of a delayed diffusive predator-prey system with food-limited and nonlinear harvesting

Based on ecological significance, a delayed diffusive predator-prey system with foodlimited and nonlinear harvesting subject to the Neumann boundary conditions is investigated in this paper. Firstly, the sufficient conditions of the stability of nonnegative constant steady state solutions of system are derived. The existence of Hopf bifurcation is obtained by analyzing the associated characteristic equation and the conditions of Turing instability are derived when the system has no delay. Furthermore, the occurrence conditions the Hopf bifurcation are discussed by regarding delay expressing the gestation time of the predator as the bifurcation parameter. Secondly, by using upper-lower solution method, the global asymptotical stability of a unique positive constant steady state solution of system is investigated. Moreover, we also give the detailed formulas to determine the direction, stability of Hopf bifurcation by applying the normal form theory and center manifold reduction. Finally, numerical simulations are carried out to demonstrate our theoretical results.


Introduction
The ordinary differential equation is well known as the logistic equation in both ecology and mathematical biology, where r and K are positive constants that stand for the intrinsic growth rate and the carrying capacity, respectively. It follows from Eq (1.1) that the U (t)/U(t) is a linear growth function of the density U(t). However, Smith [1] concluded that the Eq (1.1) does not have practical significance for a food-limited population under the influence of environmental toxicants. Based on the above fact, Smith established a new growth function in [1]. Besides, Smith also pointed out that a food-limited population requires food for both preservation and growth in its growth. For another thing, when the specie is mature, food is needed for preservation only. Therefore, the modified system is as follows where r, K and α are positive constants, and r α is the replacement of mass in the population at K. In [2], Wan et al. considered a single population food-limited system with time delay: , τ > 0, (1.3) they studied existence of Hopf bifurcations and global existence of the periodic solutions at the positive equilibrium. Eq (1.3) has been discussed in the literature by numerous scholars, they mainly concluded global attractivity of positive constant equilibrium and oscillatory behaviour of solutions of Eq (1.3) [3,4]. Su et al. [5] considered the conditions of steady state bifurcation and existence of Hopf bifurcation of food-limited population system under the dirichlet boundary condition for Eq (1.1). In addition, Gourley et al. [6] investigated the global stability, boundedness and bifurcations phenomenon of Eq (1.2) with nonlocal delay. About more interesting conclusions for food-limited system, we refer to the literature [7][8][9][10][11]. However, in real world, the interaction of prey and predator is one of the basic relations in biology and ecology. The dynamical analysis of the predator-prey system is a hot issue in biomathematics all the time, an important reason is that compared with single population system, multi-population system can exhibit complex dynamical behavior. The well-known predator-prey system has been widely studied by many ecologists and mathematicians. The authors in [12,13] developed food-limited population system to prey-predator system of functional response. Moreover, the reproduction of predator after preying upon prey is not instantaneous, but is mediated by some reaction time delay τ for gestation. Compared with the predator-prey system without delay, the time-delay predator system is more ecological significance. The delay has an effect on population dynamics and induces very rich dynamical phenomenon, see [14][15][16][17][18][19][20]. Here, we will take ratio-dependent functional response into consideration, i.e., the characteristic of consumption of prey is mUV βU+γV . The predation and reproduction of predator are not simultaneous. Taking into account the delay, the reproduction of predator from consuming the prey is nU(t−τ)V βU(t−τ)+γV . Hence, the prey-predator system with ratio-dependent and food-limited as follows where the variables U(t) and V(t) denote the densities of the prey and predator at time t, respectively. r, K, m, n and e are all positive constants that stand for the intrinsic growth rate of the prey, the carrying capacity of the prey species, predate rate of prey, converation rate from prey, and mortality rate of predator, positive constants β and γ are half saturation constant, τ(> 0) is a time delay which occurs in the predator response term and represents a gestation time of the predators.
In fact, biological resources in the predator-prey system are most likely to be harvested for economic benefit, human need to develop biological resources and capture some biological species, such as in fishery, forestry and wildlife management [21]. Hence, the demand of sustainable development for suitable resources is felt in different region of human activities to maintenance the stability of the ecosystem. We know that harvesting in population has a significant impact on the dynamic behavior of species, because of the reduction of food in the space. There are some basic types of harvesting being considering in the literature, see [21][22][23][24][25][26]. Nevertheless, a great number of mathematicians have a strong interest in nonlinear harvesting, because Michaelis-Menten type harvesting is more realistic in biology and ecology [22]. Inspired by the above discussion, system (1.4) with nonlinear prey harvesting transforms into the following system where q, E, m 1 , m 2 are also positive constants. q is the catchability coefficient, E is the effort applied to harvest prey species, m 1 and m 2 are suitable constants. Fang gave some sufficient conditions of the existence of positive periodic solutions for a food-limited predator-prey system with harvesting effect in [24,25] . In nature, the populations all require food, space and spouse and so on for survival and reproduction. When the larger the density of population is, the higher require on the environment. The shortage of food and the change in space always limit its survival and development. Naturally, the populations change position or have a diffusion to search better environment. We assume that the populations are in an isolate patch and ignore the impact of migration, including immigration and emigration. Individuals tend to migrate towards regions with lower population densities for each population. To take spatial effects into consideration, reaction diffusion system become more and more important, see [12,13,23,[27][28][29][30][31][32][33][34][35]. In this paper, we study the following reaction diffusion system: where Ω ⊂ R n (n ≥ 1) is a bounded region and it has smooth boundary ∂Ω. D 1 and D 2 , respectively, denote the diffusion coefficients of the prey and predator, and they are positive constants. ∆ denotes the Laplacian operator in R n , ϑ is outer normal vector of a boundary ∂Ω. To simplify the system (1.6), we use the following nondimensionalization: Then rewrite the system (1.6) as follows: In this paper, the domain of system is confined to Ω = [0, lπ], we define a real-valued Hilbert space The corresponding complexification is The rest of the paper is arranged as follows. In section 2, the existence and priori bound of solution of the system (1.7) are considered. In section 3, the existence of nonnegative constant steady state solutions is investigated. In section 4, the stability of the nonnegative constant steady state solutions of system (1.7) and the conditions of Hopf bifurcation and Turing instability are discussed by stability analysis and bifurcation theory. In section 5, we give the detailed formulas to determine the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions by the normal form theory and center manifold theorem for PFDEs. In section 6, some numerical simulations are carried out to illustrate the correctness of the theoretical results. Finally, some conclusions and discussions are given.

Existence of solution and priori estimate
In this section, we establish the existence of solution of system (1.7) and a priori estimate of the solution. Firstly, we have Theorem 2.1. The following statements are true for system (1.7).
(1) Obviously, F(u, v) and G(u, v) are mixed quasi-monotone in set R 2 According to the definition of upper and lower solutions in [16], denoted (u, v) = (0, 0) and (u, v) = ( u, v), where ( u, v) is the unique solution of the following system, and 0 ≤ u 0 (t, x) ≤ u 0 , 0 ≤ v 0 (t, x) ≤ v 0 . Then (u, v) and (u, v) are the upper-solution and lower-solution of system (1.7). So we have that any solution of system (1.7) is nonnegative and exists on [0, ∞), it exhibits that system (1.7) has a global solution u(t, x), v(t, x) satisfying Then u(t, x) > 0, v(t, x) > 0 by the strong maximum principle for all t > 0 and x ∈ Ω.
(3) Consider any of conditions (3a) − (3b), then it follows from the first equation of system (1.7) that ∂u ∂t There exists a t 1 , such that u(t, x) ≤ m + ε for t ≥ t 1 and x ∈ Ω, where ε is a arbitrarily small positive constant.

The nonnegative constant steady states of system (1.7)
In reality, we are interested in all the nonnegative steady state solution. Next, we give the conditions of existence of nonnegative constant steady state solutions for system (1.7).
under c − aδ > 0 and u * satisfies the following quadratic equation For the distribution of roots of Eq (3.1), we are able to get the following results about the existence of a positive constant steady state solution.
Lemma 3.1. Suppose c − aδ > 0 holds, then the following statements are true.
has a unique positive root u * = u + * . By Lemma 3.1, the following Proposition is existing.
Proposition 3.2. Suppose c − aδ > 0 holds, then the following statements are true.

Stability and bifurcation
In this section, we consider the stability of nonnegative steady state and the conditions of Hopf bifurcation and Turing bifurcation. In [37], we know that Laplacian operator −∆ exists the eigenvalue , and where I stands for 2 × 2 identity matrix and D n = − n 2 l 2 diag(d 1 , d 2 ), n ∈ N 0 . Then we obtain n 2 l 2 + a 11 a 22 , Q = −a 12 a 21 . Firstly, we discuss the stability of the singularity E 0 = (0, 0) and the semi-trivial steady state solutions E i0 (i = 1, 2, 3).
For n ≥ 0, the corresponding characteristic equation of system (1.7) at E 10 is we have N n > 0 and M n > 0 for n ≥ 0, which implies E 10 is locally asymptotically stable. In a similar method, by c − aδ > 0 or d(ρ + h) > 1, it is obvious that J(n, E 10 ) has at least one eigenvalue with a positive real part for n = 0, which implies E 10 is unstable.
We can have J(n, and d > 1, it is obvious that J(n, E 30 ) has at least one eigenvalue with a positive real part for n = 0, which implies E 30 is unstable.
(4) The proof of (4) is similar to (3), so we omit it. This completes the proof.
. Then we can have that the trivial steady state E 0 = (0, 0) is globally asymptotically stable for system (1.7) with τ = 0. This completes the proof.
Next, we investigate the stability of positive constant steady state E * . In the following discussion, we always assume For the convenience of discussion, we make the following hypothesis: Then the unique positive constant steady state E * = (u * , v * ) of system (1.7) with τ = 0 is locally asymptotically stable, that is, system (1.7) has no stationary pattern under these hypothesis.
Proof. If (H 1 ) holds, the system (1.7) has a unique positive constant steady state E * . When τ = 0, the corresponding characteristic equation of system (1.7) at E * is Obviously, n 2 l 2 + a 11 a 22 − a 12 a 21 . It follows easily from (H 2 ) ∼ (H 4 ) that all roots of Eq (4.2) have negative real parts. Hence, by Routh−Hurwitz stability criterion, the unique positive constant steady state E * is locally asymptotically stable for τ = 0 when hypothesis (H 2 ) ∼ (H 4 ) hold.
This completes the proof.
According to the work by Turing [38], positive constant steady state E * is Turing instability, implying that E * is asymptotically stable for non-spatial system (1.7) but is unstable for spatial system (1.7) with τ = 0. So we make the following hypothesis:  Proof. We know that the positive equilibrium E * = (u * , v * ) is stable for the non-spatial system (1.7), and is unstable with respect to the constant steady state solution of the spatial system (1.7) with τ = 0. The stability of non-spatial system (1.7) is guaranteed if hypothesis (H 2 ) ∼ (H 3 ) hold. When τ = 0, the corresponding characteristic equation of system (1.7) at E * is Eq (4.2). Obviously, for spatial system (1.7), it follows from (H 5 ) ∼ (H 6 ) that if there is a n ∈ N 0 such that N n + Q < 0 for 0 < k 1 < n < k 2 , which implies that Eq.(4.2) has a eigenvalue with positive real part, it is shown that the positive constant steady state E * is unstable for spatial system (1.7), that is, the diffusion-driven instability occurs. This completes the proof.
Now, by regarding τ as the bifurcation parameter, we investigate the stability and the Hopf bifurcation near the unique positive constant steady state E * . Assume that iω(ω > 0) is a root of Eq (4.1), ω satisfies the following equation By Eq (4.4), adding the squared terms for both equations yields ω 4 + P n ω 2 + Q n = 0, (4.5) where P n = M 2 n − 2N n = a 11 + d 1 Let S = ω 2 , we have S 2 + P n S + Q n = 0.   .7) is locally asymptotically stable for all τ ≥ 0.
Proof. From Eq (4.6), we have P n > 0. By (H 1 ) ∼ (H 4 ) and Theorem 4.3, we have N n + Q > 0. If (H 7 ) holds, then n 2 l 2 + a 11 a 22 + a 12 a 21 > 0 for all τ ≥ 0, which implies that Eq (4.8) has no positive roots, according to [19,Lemma 2.3], therefore the characteristic Eq (4.1) has no purely imaginary roots. Combined with Theorem 4.3, we are able to obtain that all roots of Eq (4.1) have negative real parts for any τ ≥ 0.
This completes the proof.
Then from the second equation of system (1.7), we have By Lemma 4.1 again and the any ε, we obtain that On the other hand, from the first equation of system (1.7), we have It follows from Lemma 4.1 that lim inf for t > t 2 and any ε > 0, there exists t 3 > t 2 such that for any t > t 3 , u(t, x) ≥ u 1 − ε.
Then from the second equation of system (1.7), we have By Lemma 4.1 again and the any ε, we obtain for t > t 3 + τ and any ε > 0, there exists t 4 > t 3 such that for any Meanwhile, the first equation of system (1.7) can be written as Similarly, by d > ρ, ρ − dh − dρs b > 0, we know that B(u, v 1 − ε) = 0 has a unique positive root u(v 1 − ε) and lim sup for t > t 4 and any ε > 0, there exists t 5 > t 4 such that for any t > t 5 , u(t, x) ≤ u 2 + ε. It is easily known that and u 1 , u 1 , v 1 , v 1 satisfy the following inequalities The Eq (4.9) reveals that (u 1 , v 1 ) and (u 1 , v 1 ) are coupled upper and lower solutions of system (1.7) by Definition 2.1 in [16].
Moreover, we derive the following inequality: It is easy to display that there exists a positive constant K i (i = 1, 2) such that the following Lipschitz condition holds So we can define sequences (u n , v n ) and (u n , v n ) as follows where n = 1, 2, · · ·, (u 0 , v 0 ) = (u 1 , v 1 ) and (u 0 , v 0 ) = (u 1 , v 1 ). It is easy to deduce that the sequences (u n , v n ) and (u n , v n ) satisfy the following a series of inequalities and that the limits lim Therefore, we have the following Lemma.
Clearly, we know that τ j+1 n > τ j n , therefore the following Lemma exhibits that we get a complete ordering of the Hopf bifurcation parameters τ j n .
for j ∈ N 0 .
Proof. From Eq (4.11), we have Obviously, by Eq (4.6) and Eq (4.7), we know that M 2 n − 2N n is increasing in n and Q 2 − N 2 n is decreasing in n. Hence we have that And notice that N n is strictly increasing in n ∈ [0, N * ], then we deduce that This completes the proof.
It follows from Eq (4.12) we obtain Let λ(τ) = p(τ) ± iq(τ) be the pair of root of Eq (4.1) near τ = τ j n satisfies p(τ j n ) = 0 and q(τ j n ) = ω n . Then we have the following transversality condition. Proof. Differentiating two sides of Eq (4.1) with respect to τ, we get 2ω n cos(ω n τ j n ) + M n sin(ω n τ j n )

Thus, by Eq (4.3) and Eq (4.4), we have
This completes the proof.

Stability and direction of Hopf bifurcation
In this section, we investigate the stability of the bifurcating periodic solution and direction of Hopf bifurcation by applying center manifold theorem and normal form theory of PFDEs [40]. For convenience, for fixed j ∈ N 0 , 0 ≤ n ≤ N * , we denote τ * = τ j n . Firstly, we let u(t, x) = u(τt, x) − u * , v(t, x) = v(τt, x) − v * and drop the tilde. Then system (1.7) can be transformed into: T . After the system (5.1) can be rewritten in an abstract form in the phase space C = C [−1, 0], X aṡ Linearizing Eq (5.2) at (0, 0), we can obtain the following equation About the discussion of characteristic roots in section 4, we get that the characteristic equation of Eq (5.3) has a pair of simple purely imaginary eigenvalues Λ n = {iω n τ * , −iω n τ * } and consider the following functional differential equation According to the Riesz representation theorem, there exists a 2 × 2 matrix function η(θ, τ, n)(−1 ≤ θ ≤ 0), whose elements are of bounded variation functions such that for φ ∈ C. Actually, we choose Let us define C * = C([0, 1], R 2 * ), where R 2 * is the two-dimensional vector space of row vectors, A(τ * ) denotes the infinitesimal generator of the strongly continuous semigroup induced by the solution of Eq (5.4) and A * with domain dense in C * and is the formal adjoint of A * under the bilinear form for φ ∈ C, ψ ∈ C * . Let P and P * be the center subspace, that is, the generalized eigenspace of A(τ * ) and A * associated with Λ n . A(τ * ) has a pair of simple purely imaginary eigenvalues ±iω n τ * , and A * has also a pair of simple purely imaginary eigenvalues ±iω n τ * .
According to [40], the center subspace of linear Eq (5.4) is given by P CN C, where P CN C = Φ(Ψ, φ, f n ) · f n , φ ∈ C, and C = P CN C P s C, P s C denotes the complement subspace of P CN C in C.
By the definition of A(τ * ) and Eq (5.15), we havė That is Utilizing the definition of A(τ * ) and combining Eq (5.15) and the above discussions, it follows that Thus, we can compute the following formulas which determine the direction and stability of bifurcating periodic orbits: In fact, µ 2 determines the directions of the Hopf bifurcation, β 2 determines the stability of the bifurcating periodic solutions, T 2 determines the period of bifurcating periodic solutions. Moreover, by [41], we have the following results.

Numerical simulations
In this section, we perform some results of the numerical simulations to illustrate our mathematical findings in the previous sections. All of our numerical simulations employ the Neumann boundary conditions.

Global stability
To portray the global stability of trivial steady state E 0 and the positive constant steady state E * , let a = 1.35, b = 0.01, c = 1.353, d = 0.5676, δ = 1, h = 0.28, ρ = 0.045, s = 0.001, by simple calculation, we are able to obtain (ρ + h − 1) 2 < 4(dh − ρ). It follows from Theorem 4.2 that E 0 is globally asymptotically stable. It can be seen Figure 1a   for the PDE system, the positive steady state E * becomes unstable, and system (1.7) has a stable nonconstant steady state solution, which means that a Turing instability occurs. This is shown in Figure 2. It portrays that the population is irregularly distributed in space. We also observe that the system (1.7) has a stationary spatial pattern, that is shown in Figure 3. Moreover, choose d 1 = 0.01, we are able to observe that under the effect of diffusion coefficients of d 2 , system (1.7) portrays different spatial patterns that can be periodic, stationary, as shown in      It follows from Theorem 4.7 that homogeneous Hopf bifurcations occur at τ j 0 ≈ 6.1868 + 96.3679 j for j ∈ N 0 . We use the forward Euler method to find numerical solutions to system (1.7) with τ = 0, 5.85, 6.20, respectively. From Theorem 4.3, the unique positive constant steady state E * = (0.1293, 0.4311) of system (1.7) with τ = 0 is locally asymptotically stable, it can be seen from Figure 6.
Hence, it follows from Theorem 5.1 that the direction of the bifurcation is forward, and the bifurcating period solutions are locally asymptotically stable. Moreover, the period of bifurcating periodic solutions increases. This is shown in Figure 8. If τ is increasing crosses the critical value τ 0 1 ≈ 10.0462, a spatially inhomogeneous periodic solution occurs near the positive equilibrium E * = (0.1293, 0.4311). However, the bifurcating periodic solution bifurcating from the critical value τ 0 1 must be unstable on the whole phase space since the characteristic equation always has roots with positive real parts for τ > τ 0 0 ≈ 6.1868. This is shown in Figure 9. Besides, as τ increases further, with the same initial values, the solution of system (1.7) converges to (0, 0), which implies that the increasing delay may cause the prey and predator to go extinct. As shown in Figure 10.     Figure 11. Stability region exploring the dynamics of the system (1.7) in the (h, τ) parameter space.

The effect of nonlinear harvesting
Considering the effect of nonlinear harvesting, we denote the parameters be the same as section 6.3 and h vary in [0.05, 0.12]. The stability and instability regions of the unique positive constant steady state E * for system (1.7) is depicted by mapping the nonlinear harvesting to the critical value of the delay τ in Figure 11. We observe that the nonlinear harvesting effect parameter h increases from 0.05 and 0.12, the Hopf bifurcation is occurred for lower critical values of τ. Meanwhile, we observe that the harvesting effect parameter h has a significant effect on the stability of ecosystem.

Conclusions and discussions
In this paper, we focused on a delayed diffusive predator-prey system with food-limited and nonlinear harvesting effect. Firstly, in order to prove the global stability of the solution, we gave the existence of solution and priori bound. Meanwhile, we also derived the conditions of stability of the nonnegative constant steady state solution. Moreover, the global stability of the trivial and positive constant steady state is investigated. The conditions of Turing instability and Hopf bifurcation are obtained for system (1.7), respecitvely. For the positive constant steady state solution, it is demonstrated that Hopf bifurcation can be occurred when bifurcation parameter τ increases beyond a critical value. We derived the detailed formulas to determine the properties of Hopf bifurcation.
For the system (1.7), it follows from Theorem 4.2 and Theorem 4.6 that the trivial steady state E 0 and the positive constant steady state E * are globally asymptotically stable under the certain conditions (Figure 1), respectively. From an ecological point of view, due to the food-limited in the natural environment, intraspecific competition in prey population increases, so nonlinear prey harvesting is taken into consideration. Increasing the harvesting effect parameter h properly can decrease the density of prey population so that relieve the pressure of intraspecific competition and the system becomes stable under the certain conditions. However, if h crosses a certain value, the density of prey population begins decreasing and may even become extinct (Figures 10 and 11). The diffusion induces the occurrence of Turing instability for the positive steady state E * , which means that system (1.7) has a stable nonconstant steady state solution (Figures 2 and 3). Our results also reveal the fact that delay can induce very complex dynamics phenomenon, and a positive constant steady state E * is depicted to be locally asymptotically stable when the parameter τ is less than the critical value τ * (Figure 7), and a stable periodic solutions will bifurcate from the constant steady state E * (Figure 8), when the delay τ increase and it crosses through the critical value τ * , which means that a stable and spatially homogeneous periodic solutions will occur at the critical value of delay τ * , when the delay τ continues to increase to a certain value, spatially inhomogeneous periodic solution bifurcates from the positive constant steady state E * for system (1.7) (Figure 9). When the delay τ is larger, the solution of system (1.7) tends to (0, 0), that is, the population becomes extinct eventually ( Figure 10).
From the biological point of view, it is clear that delay, nonlinear harvesting and diffusion effect have a significant impact on the stability for ecosystem.
However, there exists many problems will need to be investigated for system (1.7). Firstly, the selection results of Turing patterns are obtained by weakly nonlinear analysis [29,30]. Secondly, the diffusion of the population is random and homogeneous in this paper, actually, individuals often perform a nonlocal diffusion or heterogeneous diffusion in the real world. Finally, spatial memory and cognition of animals has drawn much attention in the mathematical modeling of animal movement