ON A PREDATOR PREY MODEL WITH NONLINEAR HARVESTING AND DISTRIBUTED DELAY

. A predator prey model with nonlinear harvesting (Holling type-II) with both constant and distributed delay is considered. The boundeness of solutions is proved and some suﬃcient conditions ensuring the persistence of the two populations are established. Also, a detailed study of the bifurcation of positive equilibria is provided. All the results are illustrated by some numerical simulations.

1. Introduction. In [5] the authors propose a predator prey model to study the impact of harvesting on a two species community. In particular the attention is focused on predator harvesting due to its importance in controlling the predator population and to prevent the extinction of the prey species. The published literature suggests that a nonlinear harvesting (see for example [2,4,7,10,14]) is capable describing complex behaviours. Moreover, the choice of nonlinear harvesting (with a functional response known as Holling type-II) is also motivated by the fact that nonlinear harvesting function exhibits saturation effects with respect to both the stock abundance and the effort-level of harvesting.
This preserves positivity of solutions and prevents blow up phenomena, that is, the model becomes more realistic. In fact, in [5] a detailed stability and bifurcation analysis of the following system is carried out: . (1) In [8] the authors consider an extended version of the model taking into account diffusion in order to describe individual motions within their habitats. This is motivated by an increasing interest in the literature (for example [12,13,15,16]). They consider the following model with and without time delay: Introducing time delay in the model is fully justified as can be seen in the previous literature (see for example [1]). Moreover, it is now admitted that diffusion stabilises pattern formation while time delay may destabilize. The model (2) has been extensively studied in [8] where some conditions ensuring Turing and Hopf bifurcations are provided. Moreover, numerical simulations show the formation of many different spatial patterns (such as spots, strips, mixture of spots and strips, spiral, patchy structure, chaos...) which are affected and transformed in presence of time delay. Motivated by both works [5] and [8], in this paper we come back to the finite dimensional model and consider a delayed version of the original model in Section 2. In Section 3 we introduce the model with distributed delay and in Section 4 we study persistence of the system and boundedness of solutions. Section 5 is devoted to the analysis of stability for positive equilibria and the problem of bifurcations, while, in Section 6, some remarks for future investigation are included.
2. The model. In this section we consider the case of constant delay; in other words, we analyse the following model: where u and v denote population biomass of prey and predator, respectively, τ is the time delay, and all parameters ρ, α, η, ε are positive. In the previous model, a Holling type II response is considered for the population of predators.
The positive equilibria of model (3) are the same as those of the model without delay. They are the solutions of system where u > 0 and v > 0. The analysis of the number of positive equilibria, when τ = 0, has been carried on in [5], where the authors showed that the number of these equilibria depends on the expression of η + αε − ε. In fact, they proved the following result: Theorem 1. The sufficient conditions ensuring the existence of positive steady states, for the zero-delay case, can be classified into four cases: 2. Let η + αε > ε, with α + ε < 1 and (α + ε − 1) 2 = 4(η + αε − ε). Then, system (3) has a unique positive equilibrium (ū,v), which is a saddle point, Let η + αε = ε, with α + ε < 1. Then, system (3) possesses a unique positive equilibrium (u * + , v * + ) which is also a saddle point. 4. Let η + αε < ε. Then, system (3) has a unique positive equilibrium (u * + , v * + ), which is locally asymptotically stable whenever ρ < . Now we analyse the system with a constant delay. Let (u * , v * ) be an equilibrium of (3). which is locally asymptotically stable in the zero-delay case. Translating this point to the origin, the linearization of the resulting system of (3) at this point possesses a characteristic equation given by where By Rouche's Theorem [3] and the continuity with respect to τ , the existence of roots of Eq. (4) with positive real parts guarantees the existence of purely imaginary roots and viceversa. From this, we shall be able to find conditions for all eigenvalues to have negative real parts. Let λ(τ ) = β(τ ) + iω(τ ), with β and ω real. As the equilibrium (u * , v * ) is stable, we have β(0) < 0. By continuity, if τ > 0 is sufficiently small, we still have β(τ ) < 0 and (u * , v * ) is still stable. The change of stability will occur at some values of τ for which β(τ ) = 0 and ω(τ ) = 0, i.e. λ(τ ) is purely imaginary. Let λ = iω (ω > 0) be a root of (4). Then Equating the real and imaginary parts of both sides, we have Squaring and adding Eqs. (5) we obtain where Then, it is not difficult to prove the following result: Lemma 2. The number of positive roots of (6) depends on the relationship between A and B as follows.

1)
The fixed point (u * , v * ) is locally asymptotically stable for all τ ≥ 0 if one of the following conditions is fulfilled: , its stability change a finite number of times for τ > τ (+) 0 ) and eventually it becomes unstable for τ > τ for which a stability switch occurs.
In the following lines we illustrate the case 3), that is, if In the first simulation τ = 2 and the fixed point is locally asymptotically stable (see figure 1), in the second one there is a stability switch, the fixed point is unstable and a stable limit cycle appears (see figure 2). 3. Model with distributed delays: introduction. In the present section, we consider a different version of model (3) obtained by introducing distributed delays where g(·) is a gamma distribution, i.e.
with m a positive integer, that determines the shape of the weight function, and T ≥ 0 a parameter associated with the mean time delay of the distribution. We will consider only the cases of the so-called weak kernel function (m = 1), i.e. the importance of events in the past simply decreases exponentially the further one looks into the past, and strong generic function (m = 2), i.e. a particular time in the past is more important than any other. To this end, we define the new variables then using the linear chain trick technique [9], system (11) can be transformed into the following equivalent system without delay and .
In the next sections we investigate the persistence together with the boundness of solutions and stability of positive equilibria of systems (13) and (14). 4. Model with distributed delays: persistence. For the sake of clarity in our presentation, we analyse in details only the case m = 1, since the same arguments work for the case m = 2.
4.1. Case m = 1. We start by studying the dynamics on the boundary of R 3 + . The plane u = 0 is invariant and the dynamics is described by We have two fixed points: (0,0) and (0, −ε − η/α), which is negative and, as a consequence, is outside ∂R 3 + . We observe thatv < −αρv.
Thus all solution on the plane u = 0 satisfy The plane v = 0 is invariant and the system restricted to it becomes We have two fixed points (0, 0) and (1,1), and linearisation provides the following Jacobian matrix , whose determinant and trace are respectively We conclude that (0, 0) is unstable with its stable manifold being the x axis, while the point (1, 1) is stable.
If T < 1, by the Bendixon Criterion we exclude the existence of periodic orbits in the whole positive quadrant of v = 0. In particular, if T < 1 a periodic orbit may lie across the line x = 1 − 1/T , however we can exclude the existence of a limit cycle in a box with center the fixed point (1,1) with side 2/T .
If T ≥ 1 we can construct the simple Dulac fuction In fact we have ∂ ∂u Then, by Dulac's theorem, there are no periodic orbits in the same set.
In the following lines we will prove that solutions are bounded in this set. By the equation it is easy to exclude that one or both u and x monotonically diverges, then we will show that orbits cannot spiral away around point (1,1) (see figure 3 below).
(1) Suppose that where D is a positive constant. Then there exists a sequence t n → ∞ such that From the second equation of the system we havė then by properties of t n there exists N > 0 such thaṫ and this is a contradiction.
Then there exists a sequence t n → ∞ such that Then, from the second equation we havė from whichẋ(t n ) → ∞, which is a contradiction.

2711
(3) The remaining case, that is lim sup x(t) = lim sup u(t) = +∞, can be excluded using the same reasonings. if the solutions start on the v axis (u = x = 0) then v → 0, since the v axis is invariant.
Now we pass to study the fixed points of system (13): where the last fixed point is solution of the system The functional Jacobian of the system is For the origin we have The eigenvalues are then the fixed point is unstable with its stable manifold entirely contained in the invariant plane u = 0. For the other boundary fixed point we have the fixed point (1, 1, 0) is unstable with stable manifold contained in the invariant plane v = 0. If condition (15) is fulfilled, then all the boundary fixed points are unstable. Moreover, the only invariant sets contained in the boundary of R 3 + are those fixed points.
Since no cycles are possible connecting the invariant sets of ∂R 3 + , and all the stable manifolds of the invariant sets are entirely contained in ∂R 3 + , we obtain the following result (see [11]) where we establish a sufficient condition ensuring the existence of at least one positive equilibrium: Theorem 5. If condition (15) is fulfilled then the system is persistent.
In fact, this sufficient condition implies more as will be proved in the next result. Proposition 6. If condition (15) is fulfilled then there exists a unique interior fixed point.
Then, by (15), these points are always real, v * 1 is positive, while v * 2 is negative. For simplicity we set v * = v * 1 . The other coordinate satisfies A PREDATOR PREY MODEL ...

2713
Then u * is positive if Observe that h can be rewritten in the following way Then, u * is positive and, as a consequence, x * = u * is also positive.
In the next section we will study in detail the stability of this fixed point while we end this section by proving the boundedness of solutions. Then there exists a sequence t n → +∞ such that u(t n ) is increasing witḣ In other words at each t n the function u(t) attains a local maximum. Then, by the first equation of the system and as a consequence x(t n ), v(t n ) ≤ 1, for all n. From the second equation of the systeṁ since u(t n ) is increasing, there exists N ∈ N such that u(t n ) > 1, ∀n ≥ N Moreover, passing to the limit whence x(t n ) → ∞, and this is a contradiction. Then u is bounded. We set From the second equation of the system we havė from which Then there exists a sequence t n → ∞ such that v(t n ) is increasing and We observe that,by the third equation of the system, we havė .
We can exclude that u(t n ) → 0: for if u(t n ) → 0 there exists N ∈ N such that u(t n ) < α, for n > N, and we havev(t n ) < 0 for n > N , that is not compatible with the hypothesis that v → ∞.
From the first equation of the system we can writė Then, passing to the limits, from which there exists N ∈ N such that u(t n ) < α, ∀n > N. The quantity v M satisfiesũ whereũ is the value attained by v when it attains the value v M . As a consequence, we obtain the lower bound Note that, from the third equation of the system it follows thaṫ v ≤ ρv u M − α − η ε + v .
Since the system is persistent we deduce the following lower bound for u M u M ≥ α + η ε .

Case m=2.
We have previously noted that the same arguments of the previous subsection work for this case. We have again only two boundary fixed points O = (0, 0, 0, 0), P = (1, 1, 1, 0).
The analysis of the stability of O is the same as in the case m = 1, in particular it is unstable and its three dimensional stable manifold is contained in R 4 + . For P = (1, 1, 1, 0) the Jacobian matrix is Then, an eigenvalue is given by Thus we must have again (as for m = 1) that This is a sufficient condition for persistence as in the case m = 1. The remaining eigenvalues are the roots of the following polynomial where we have set c = 2/T . We observe that This means that P (λ) is increasing for λ > 0 and therefore it has no positive root and at least one negative root. We observe that if the polynomial P (λ) has three negative roots and, as a consequence the above condition, (16) becomes also necessary for the persistence of the system. When c < 27/4 we have a pair of complex eigenvalues and we must look at the sign of the real part. If λ 2 = h is the real negative root of P (λ), then the real part of the remaining eigenvalues is −(2c + h). Then condition (16) is necessary if which means c > 1/2. Thus, when c < 1/2, condition (16) is not necessary for the instability of P , however, as in the case m = 1, it ensures the existence of one positive equilibrium (u * , z * , y * , v * ).

5.
Model with distributed delays: stability and bifurcation of the positive equilibrium. In this section we analyse the stability of the positive equilibrium for the model with distributed delay with m = 1 and m = 2. In particular we prove the existence of Hopf Bifurcation and study its stability.

5.1.
Case m = 1. For simplicity we perform a change of variable in order that the positive equilibrium (u * , x * , v * ) of (13) is shifted to the origin. The characteristic equation of the linearised system of (13) (after the change of variable) at the origin is λ 3 + a 1 λ 2 + a 2 λ + a 3 = 0, (17) where By the Routh-Hurwitz criterion, the necessary and sufficient conditions for Eq. (17) to have negative real parts are a 1 > 0, a 3 > 0 and a 1 a 2 > a 3 .
Thus, we have a pair of purely imaginary roots λ 1,2 = ±iω * , where ω * = (u * − b)/T * + ρu * v * > 0, and a real root λ 3 = −(1−bT * )/T * < 0. According to the Hopf bifurcation theorem, we need to verify that λ = iω * is a simple root of (17) and the transversality condition holds. Differentiating (17) with respect to T , we obtain where . If λ = iω * is a multiple root of (17), then one has −a 1 (T * )ω 2 * + ia 2 (T * )ω * + a 3 (T * ) = 0, leading to the contradiction ω * = 0. Recalling that ω * = a 2 (T * ), after some calculations, it follows from (24) that A positive sign in the previous expression corresponds to crossings of the imaginary axis from right to left, and a negative sign implies crossings from left to right. Notice that the sign is positive if u * − b ≥ 0, while it can assume any value otherwise.

T. CARABALLO, R. COLUCCI AND L. GUERRINI
We consider two cases: T = 40 and T = 40.5. In the first case (see figure 4) the fixed point is locally asymptotically stable and the solutions converge to it. In figure 5 we represent the second case, where a stable limit cycle appears and the fixed point becomes unstable. where By the Routh-Hurwitz criterion we have that the positive equilibrium (u * , y * , z * , v * ) of (14) is locally asymptotically stable if a 1 > 0, a 3 > 0, a 4 > 0 and a 1 a 2 a 3 > a 2 3 + a 2 1 a 4 . In particular, it follows a 2 > 0. A direct calculation shows that it must hold and The ϕ(T ) = 0 locus partitions the plane into a stable region and unstable region, and it is called the partition curve. Assume that there exists T = T * such that ϕ(T * ) = 0, that is a 1 a 2 a 3 − a 2 3 − a 2 1 a 4 = 0, when T = T * . On the partition curve, the characteristic equation (25) can be factored as follows a 1 λ 2 + a 3 a 1 λ 2 + a 2 1 λ + a 1 a 2 − a 3 = 0. Its solutions are which are clearly purely imaginary, and whose real part is different from zero. Next, we select the delay T as the bifurcation parameter and consider the roots of the characteristic equation (25) as continuous functions of T. Plugging λ = λ(T ) into (25), and differentiating it with respect to T , we obtain i.e. where . At λ = iω * , ω * = ω(T * ) > 0, the real part of the derivative is obtained as where ϕ (T * ) = a 1 a 2 a 3 + a 1 a 2 a 3 + a 1 a 2 a 3 − 2a 3 a 3 − 2a 1 a 1 a 4 − a 2 1 a 4 , with all the a j and a j (j = 1, 2, 3, 4) evaluated at T = T * . Since we can conclude that there is a crossing of the imaginary axis at T = T * from the left half plane to the right half plane if ϕ (T * ) < 0, and from right to left as T increases if ϕ (T * ) > 0. Summarizing all the previous analysis, we have the following result.
The function ϕ has only one positive and one negative zero at 2.13505 and -1.32525 respectively. Then we have T * = 2.13505. In the following simulation we start with T = 1.5 < T * , for which the fixed point (u * , z * , y * , v * ) is locally asymptotically stable (see figure 6) and let us increase it to values greater than the critical one. We observe that for T = 2, a stable limit cycle appears (se figure 7) and then for the values T = 2.5, 3.132, 3.2 we observe an increasing period for the limit cycle (see figures 7, 9 and 10) until reaching a possible chaotic attractor for T = 4 (see figure 11).  Figure 6. For T = 1.5 < T * the fixed point (u * , z * , y * , v * ) is locally asymptotically stable.
6. Concluding remarks. We have observed an interesting phenomenon: the choice of the parameter m seems to be important for the dynamical behaviour of the system. For m = 1 the numerical simulations suggest two possibilities (in the case of persistence), that is, attracting fixed point or attracting limit cycle (periodic orbit) generated by Hopf-Bifurcation. On the other side, for m = 2, we have observed a richer dynamics. In fact, cycles with increasing periods have been detected by the simulations with the possibility that a chaotic attractor may exist. We consider interesting for further research to investigate the dependence on m and to carry out an exhaustive analysis of all possible bifurcations which may arise for m ≥ 2. Another interesting target for research is to study the existence and dimension of the (possible) chaotic attractor observed in the simulations. . For T = 2.5 > T * the fixed point (u * , z * , y * , v * ) is unstable and a stable limit cycle appears. In the figures it is represented the limit cycle together with the time series of u and v respectively.   Figure 9. For T = 3.132 > T * we observe a limit cycle with three periods. In the figures it is represented the limit cycle together with the time series of u and v respectively.   Figure 11. For T = 4 > T * we observe a possible chaotic attractor which is represented together with the time series of u and v respectively.