Dynamics of a prey-predator system with modified Leslie-Gower and Holling type II schemes incorporating a prey refuge

We study a modiﬁed version of a prey-predator system with mod-iﬁed Leslie-Gower and Holling type II functional responses studied by M.A. Aziz-Alaoui and M. Daher-Okiye. The modiﬁcation consists in incorporating a refuge for preys, and substantially complicates the dynamics of the system. We study the local and global dynamics and the existence of cycles. We also investigate conditions for extinction or existence of a stationary distribution, in the case of a stochastic perturbation of the system.


Introduction
We study a two-dimensional prey-predator system with modified Leslie-Gower and Holling type II functional responses.This system is a generalization of the system investigated in the papers by M.A. Aziz-Alaoui and M. Daher-Okiye [3,9].
A novelty of the present paper is that we add a refuge in a way which is different from [7], since the density of prey in our refuge is not proportional to the total density of prey.This kind of refuge entails a qualitatively different behavior of the solutions, even for a small refuge, contrarily to the type of refuge investigated in [7].Let us emphasize that, even in the case without refuge, our study provides new results.
In the first and main part of the paper (Section 2), we study the system of 2 [3,9] with refuge, but without stochastic perturbation: (1.1) In this system, • x ≥ 0 is the density of prey, • y ≥ 0 is the density of predator, • µ ≥ 0 models a refuge for the prey, i.e, the quantity (x−µ) + := max(0, x− µ) is the density of prey which is accessible to the predator, • ρ 1 > 0 (resp.ρ 2 > 0) is the growth rate of prey (resp. of predator), • β > 0 measures the strength of competition among individuals of the prey species, • α 1 > 0 (resp.α 2 > 0) is the rate of reduction of preys (resp. of predators) • κ 1 > 0 (resp.κ 2 > 0) measures the extent to which the environment provides protection to the prey (resp.to the predator).
When the predator is absent, the density of prey x satisfies a logistic equation and converges to ρ1 β , so we assume that 0 ≤ µ < ρ 1 β .
The last term in the right hand side of the first equation of (1.1), which expresses the loss of prey population due to the predation, is a modified Holling type II functional response, where the modification consists in the introduction of the refuge µ.The predation rate of the predators decreases when they are driven to satiety, so that the consumption rate of preys decreases when the density of prey increases.Similarly, if its favorite prey is absent (or hidden in the refuge), the predator has a logistic dynamic, which means that it survives with other prey species, but with limited growth.The last term in the right hand side of the second equation, of (1.1) is a modified Leslie-Gower functional response, see [20,30].Here, the modification lies in the addition of the constant κ 2 , as in [3,9], as well as in the introduction of the refuge µ.It models the loss of predator population when the prey becomes less available, due its rarity and the refuge.
Setting, for i = 1, 2, we get the simpler equivalent system (1.2) where 0 ≤ m < 1, all other parameters are positive, and (x, y) takes its values in the quadrant R + × R + .
In this first part, we study the dynamics of Equation (1.2), which is complicated by the refuge parameter m.However, even in the case when m = 0, we provide some new results.We first show the persistence and the existence of a compact attracting set.Then, we study in detail the equilibrium points (there can be 3 distinct non trivial such points when m > 0) and their local stability.We also give sufficient conditions for the existence of a globally asymptotically stable equilibrium, and we give some sufficient conditions for the absence of periodic orbits.A stable limit cycle may surround several limit points, as we show numerically.
In a second part (Section 3), we study the stochastically perturbed system where w = (w 1 , w 2 ) is a standard Brownian motion defined on the filtered probability space (Ω, F, (F t ), P), and σ 1 and σ 2 are constant real numbers.This perturbation represents the environmental fluctuations.There are many ways to model the randomness of the environment, for example using random parameters in Equation (1.2).Since the right hand side of Equation (1.2) depends nonlinearly on many parameters, the approach using Itô stochastic differential equations with Gaussian centered noise models in a simpler way the fuzzyness of the solutions.The choice of a multiplicative noise in this context is classical, see [28], and it has the great advantage over additive noise that solutions starting in the quadrant [0, +∞[×[0, +∞[ remain in it.Furthermore, the independence of the Brownian motions w 1 and w 2 reflects the independence of the parameters in both equations of (1.2).Another possible choice of stochastic perturbation would be to center the noise on an equilibrium point of the deterministic system, as in [4].But we shall see in Theorem 2.3 that Equation (1.2) may have three distinct equilibrium points.Furthermore, as in the case of additive noise, this type of noise would allow the solutions to have excursions outside the quadrant [0, +∞[×[0, +∞[, which of course would be unrealistic.
We show in Section 3 the existence and uniqueness of the global positive solution with any initial positive value of the stochastic system (1.3), and we show that, when the diffusion coefficients σ 1 > 0 and σ 2 > 0 are small, the solutions to (1.3) converge to a unique ergodic stationary distribution, whereas, when they are large, the system (1.3) goes asymptotically to extinction.Small values of σ 1 and σ 2 are more interesting for ecological modeling, because they make solutions of (1.3) closer to the prey-predator dynamics.The effect of such a small or moderate perturbation is the disparition of all equilibrium points of the open quadrant ]0, +∞[×]0, +∞[, replaced by a unique equilibrium, the stationary ergodic distribution, which is an attractor.
The last part of the paper is Section 4, where we make numerical simulation to illustrate our results.

Dynamics of the deterministic system
In this section, we study the dynamics of (1.2).
Throughout, we denote by v the vector field associated with (1.2), and The right hand side of (1.2) is locally Lipschitz, thus, for any initial condition, (1.2) has a unique solution defined on a maximal time interval.

Persistence and compact attracting set
The next result shows that there is no explosion of the system (1.2).It also shows a qualitative difference brought by the refuge: when m = 0, the density of prey may converge to 0, whereas, when m > 0, the system (1.2) is always uniformly persistent.
In particular, the system (1.2) is uniformly persistent.
In the case when m > 0, we need to prove that lim inf x(t) ≥ m, provided that x(0) > 0. Actually we have a better result, since, if x(0) ≤ m, then x coincides with the solution to the logistic equation ẋ = x(1 − x) as long as x does not reach the value m, that is, .
In particular, we have This implies that, for any > 0, and for t large enough (depending on x(0)), we have x(t) ≤ 1 + .We deduce that, for any > 0, and for t large enough, we have which implies that, for t large enough, say, t ≥ t 0 , (2.10) .
(b) We have already seen that x(t) ≥ m for t large enough, let us now check that x(t) ≤ 1 for t large enough.Since A is invariant, we only need to prove this for x(0) > 1.Let > 0 such that k 2 − > 0. Let δ > 0 such that δ + m < 1 and such that From the first inequality in (2.12), we have y(t) ≥ k 2 − for t large enough, say t ≥ t 0 .From (2.8), we can take t 0 large enough such that, for t ≥ t 0 , we have also x(t) ≤ 1 + .Using (2.13), we deduce, for t ≥ t 0 and x(t) ≥ 1 − δ, Thus x decreases with speed less than − aδ(k2− ) 1+ −m < 0. Thus x(t) ≤ 1 − δ for t large enough.
We can now repeat the reasoning of (2.9) and (2.10), replacing by −δ, which yields that lim sup y(t) ≤ L − δ.In particular, y(t) < L for t large enough.
To prove that y(t) > k 2 for t large enough, let us first sharpen the result of (2.5).This is where we use that m > 0. Let δ > 0, with m + δ < 1.If |x − m| < δ, we have From (2.12), we deduce that, for any > 0, and t large enough, depending on , we have y(t) ≤ L + and x(t) ≥ m, (we do not write t here for the sake of simplicity).For δ small enough, we have D > 0. Thus, if m > 0, we can find δ > 0 small enough (depending on m), such that, when .
This proves that lim inf and that y > k 2 for t large enough.
(c) Assume now that m = 0. Since the first part of the proof of (b) is valid for all m ≥ 0, we have already proved that x(t) < 1 and y(t) < L for t large enough.

k1
. By the second inequality in (2.12), we have, for t large enough Thus lim inf x(t) ≥ K .As is arbitrary, this proves (2.2).From (2.2) and the first inequality in (2.12), we deduce that (1.2) is uniformly persistent.
(cii) Assume now that ak 2 < k 1 ≤ aL.Observe first that, if lim sup x(t) < l for some l > 0, then, for t large enough, we have x(t) < l, thus ẏ(t) < by(1 − y/(k 2 + l)).We deduce that Let us now rewrite the first equation of (1.2) as where U (x) = (−1/a)(x − 1)(x + k 1 ).Since ak 2 < k 1 , the point E 2 lies below the parabola y = U (x), thus in the neighborhood of E 2 , for x > 0, we have ẋ > 0. By (2.16), if lim sup x(t) < l for some l > 0, then for t large enough, the point (x(t), y(t)) remains in the rectangle R = [0, l] × [0, k 2 + l].But if, furthermore, l is small enough such that R lies entirely below the parabola y = U (x), then, when (x(t), y(t)) ∈ R, we have ẋ(t) > 0, which entails that x(t) is eventually greater than l, a contradiction.This shows that, for l > 0 small enough, we have necessarily lim sup t→∞ x(t) ≥ l.
Let us now calculate the largest value of l such that (x, y) ∈ R implies y < U (x), that is, the largest l such that min From the concavity of U , the minimum of U on the interval [0, l] is attained at 0 or l.Thus the optimal value of l is the minimum of This proves (2.3).
• Assume that 1 Then, the parabola ỹ = V (x) is above the line ỹ = x for all x in the interval ]0, l[, where l is the non-zero solution to V (x) = x, that is, Let us show that lim sup x(t) ≥ l.Assume the contrary, that is, lim sup x(t) < δ for some δ < l.For t large enough, say, t ≥ t δ , we have x(t) < δ.Let us first prove that |ỹ(t)| < δ for t large enough.If ỹ(t δ ) < δ, we have, for all t ≥ t δ , as long as ỹ(t) < δ, which proves that ỹ(t) < δ after a finite time.
We have proved that, for t large enough, ( which shows that x(t) > δ for t large enough, a contradiction.This proves (2.4).
Then, the portion of the parabola ỹ = V (x) which lies in ]0, +∞[×]−k 2 , +∞[, is below the line ỹ = x.This means that, for any > 0 such that k 2 − > 0, the system (1.2) has no other equilibrium point than E 2 in the invariant attracting compact set we can use exactly the same arguments as in the case when

Trivial critical points
The right hand side of (1.2) has continuous partial derivatives in the first quadrant R + × R + , except on the line x = m if m > 0. The Jacobian matrix of the right hand side of (1.2 We start with a result on the obvious critical points of (1.2) which lie on the axes.
Proposition 1.The system (1.2) has three trivial critical points on the axes: • E 0 = (0, 0), which is an hyperbolic unstable node, • E 1 = (1, 0), which is an hyperbolic saddle point whose stable manifold is the x axis, and with an unstable manifold which is tangent to the line an hyperbolic saddle point whose stable manifold is the y axis, with an unstable manifold which is tangent to the line bx In this case, the y axis is the stable manifold, and there is a center manifold which is tangent to the line y − k 2 = x.
(Compare with the case (c) of Theorem 2.1).
Proof.The nature of E 0 , E 1 , and E 2 , is obvious since The results on stable and unstable manifolds of hyperbolic saddles are straightforward.In the case when E 2 is semi-hyperbolic, since it is either a topological node or a topological saddle (see [11,Theorem 2.19]), the nature of E 2 follows from Part (ciii) of Theorem 2.1.In the topological saddle case, that is, when m = 0 with ak 2 = k 1 and 1 − k 1 − a > 0, the eigen values of J (0, k 2 ) are −b and 1, with corresponding eigenvectors (0, 1) and (1, 1).Clearly, the y axis is the stable manifold.The change of variables We can thus write where A and B are analytic and their jacobian matrix at (0, 0) is 0. In the neighborhood of (0, 0), the equation 0 = −Y b+B(X, Y ) has the unique solution Y = f (X), where and From [11, Theorem 2.19], we deduce that there exists an unstable center manifold which is infinitely tangent to the line Y = 0.

Counting and localizing equilibrium points
Let us now look for critical points outside the axes, i.e., critical points E = (x, y) with x > 0 and y > 0. From the results of Section 2.1, such points are necessarily in A, in particular they satisfy x ≥ m.We have, obviously: Lemma 2.2.The set of equilibrium points of (1.2) which lie in the open quadrant ]0, +∞[×]0, +∞[ consists of the intersection points of the curves Furthermore, these points lie in A.
We shall see that, when m > 0, the system (1.2) has always at least one equilibrium point in ]0, +∞[×]0, +∞[, whereas, for m = 0, some condition is necessary for the existence of such a point.
• When m > 0, the solutions to (2.20) lie at the abscissa of the intersection of the parabola z 0 and, for x > 1, we have P (x) < 0 and Q(x) > 0, thus P (x) − Q(x) > 0. This implies that the curves of P and Q have at least one intersection whose abscissa is greater than m, and that the abscissa of any such intersection lies necessarily in the interval ]m, 1[.The change of variable By Routh's scheme (see [14]), the number p of roots of (2.22) with positive real part, counted with multiplicities, is equal to the number of changes of sign of the sequence (2.24) provided that all terms of V are non zero.Thus p = 3 when and, in all other cases, p = 1.When p = 1, we know that the number n of real positive roots of R is exactly 1.When p = 3, we have either n = 1 if R has two complex conjugate roots, or n = 3.So, we need to examine when all roots of R are real numbers.A very simple method to do that for cubic polynomials is described by Tong [32]: a necessary and sufficient condition for R to have three distinct real roots is that R has a local maximum and a local minimum, and that these extrema have opposite signs.The abscissa of these extrema are the roots of the derivative R (X) = 3X 2 + 2α 2 X + α 1 , thus R has three distinct real roots if, and only if, the following conditions are simultaneously satisfied: (ii) R(x)R(x) < 0, where x and x are the distinct roots of R .
If R(x)R(x) = 0 with ∆ R > 0, the polynomial R still has three real roots, two of which coincide and differ from the third one.If R(x)R(x) = 0 with ∆ R = 0, it has a real root with multiplicity 3, which is x = x, and if ∆ R = 0 with R(x)R(x) = 0, it has only one real root.Fortunately, all radicals disappear in the calculation of R(x)R(x): In particular, Conditions (i) and (ii) can be summarized as Let us now examine what happens when one term of the sequence V in (2.24) is zero.We skip temporarily the case α 0 = 0, which is equivalent to m = 0.
and α 2 and α 1 have opposite signs, because α 0 < 0. Thus, in that case, R has a unique positive root, which is and increasing in [ √ −α 1 , +∞[.Since R(0) < 0, R has only one positive root.Thus, in that case too, R has a unique positive root.
From the preceding discussion, we deduce the following theorem: Theorem 2.3.Assume that m > 0. With the notations of (2.23), the number n of distinct equilibrium points of the system Remark 2. Numerical computations show that all cases considered in Theorem 2.3 are nonempty.See Figure 1 for an example of positive numbers (a, k 1 , k 2 , m) satisfying (2.25) and (2.26).
• When m = 0, the system (1.2) is exactly the system studied by M.A. Aziz-Alaoui and M. Daher-Okiye [3,9].As x is assumed to be positive, (2.20) is equivalent to the quadratic equation which can be written where α 2 = a+k 1 −1 and α 1 = ak 2 −k 1 as in (2.23).The associated discriminant is thus a sufficient and necessary condition for the existence of solutions to (2.27) in R is ∆ ≥ 0, i.e., k 2 must not be too large: Since the sum of the solutions to (2.27) is −α 2 and their product is α 1 , we deduce the following result: Theorem 2.4.Assume that m = 0.With the notations of (2.23), the number n of distinct equilibrium points of the system Remark 3. If m = 0 and n = 0, the point E 2 is the only equilibrium point in the compact invariant attracting set ).This gives a more general condition of global attractivity of E 2 than the result given in Parts (ciii) and (civ) of Theorem 2.1.
Remark 4. Since the roots of the polynomial R defined by (2.22) depend continuously on its coefficients, Theorem 2.4 expresses the limiting localization of the equilibrium points of (1.2) when m goes to 0. In particular, the case (a) of Theorem 2.4 is the limiting case of (a) in Theorem 2.3.Indeed, it is easy to check that Condition (2.30), with m = 0, is a limit case of (2.25) and (2.26).This means that, in the case (a) of Theorem 2.
thus it has at most one positive solution.In that case, the coordinates of the unique non trivial equilibrium point E * can be explicited in a simple way, and we have

Local stability
Let E * = (x * , y * ) be an equilibrium point of ( The characteristic polynomial of J (x * , y * ) is where (2.33) The roots of χ are real if, and only if, ∆ χ ≥ 0, where The point E * is non-hyperbolic if one of the roots of χ is zero (that is, if p = 0), or if χ has two conjugate purely imaginary roots (that is, if s = 0 with p > 0).If only one root of χ is zero, that is, if p = 0 with s = 0, the point E * is semi-hyperbolic.
Remark 6.An obvious sufficient condition for any equilibrium point E * ∈ A to be stable hyperbolic is m ≥ 1/2, since x * > m.This condition can be slightly improved, as we shall see in the study of global stability (see Theorem 2.11).
Application of the Poincaré index theorem When E * is an hyperbolic equilibrium, its index is either 1 (if it is a node or focus) or −1 (if it is a saddle).Let n be the number of distinct equilibrium points, which we denote by E * 1 , ..., E * n , and let I 1 , ..., I n their respective indices.As we shall see in the proof of the next theorem, by a generalized version of the Poincaré index theorem, we have I 1 + ... + I n = 1.When all equilibrium points are hyperbolic, this allows us to count the number of nodes or foci and of saddles.• If n = 1, the unique equilibrium point in the interior of A is a node or a focus.
• If n = 3, the system (1.2) has one saddle point and two nodes or foci in the interior of A.
• If n = 2, one equilibrium point is a node or focus, and the other is a saddle.
• If n = 1, the unique equilibrium point in the interior of A is a node or a focus.
Proof.Let N (respectively S) denote the number of nodes or foci (respectively of saddles) among the hyperbolic singular points which lie in A. 2. Assume now that m = 0. We use the same reasoning as for m > 0, but with a different domain.Instead of A, we consider the domain for a small > 0. Thus B contains E 2 .
• With the notations of (2.17), if ak 2 > k 1 , we have y > U (x) for x = 0 and for all y ∈ [k 2 , L].We have By continuity of v, we can choose > 0, with < k 1 , such that the inequality y > U (x) remains true on the rectangle • If ak 2 < k 1 , E 2 is a saddle point, thus, constructing B and B as precedingly, we have now S = S − 1 and N = N .Furthermore, the vector field v is no more outward directed along the whole boundary of B .We use Pugh's algorithm [31] to compute N − S : taking small enough such that the vector field v does not vanish on ∂B , we have − ), where χ denotes the Euler characteristic, R 1 − is the part of the boundary of B where v is directed outward, and R 2 − is the part of ∂R 1 − where v points to the exterior of R 1 − .Since k 2 < k 1 /a, we see that the parabola y = U (x) crosses the line {x = − ; y > k 2 } at some point (− , r), so that the part of the boundary of B where v points outward is the segment {− } × [k 2 − , min(r, L)].Thus, for small , R 1 − is an arc whose extremities are tangency points.Observe also that, since v 1 < 0 for x < 0 and v 2 < 0 for y > k 2 + x > 0, the field v points toward the interior of R 1 − at those tangency points, thus R 2 − is empty.Formula (2.34) becomes Since, by Theorem 2.4, we have N + S = 1, we deduce that N = 1 and S = 0.
For p = 0, the Jacobian matrix J (x * , y * ) is The change of variables The coordinates of v are, in the basis We can thus write where A and B are analytic and their jacobian matrix at (0, 0) is 0 and λ > 0. It is not easy to determine v = f (u) the solution to the equation λv + B(u, v) = 0 in a neighborhood of the point (0, 0), for that we use implicit function theorem.We find: Case 1: If κ 3 − ky * a + ak 1 κ = 0, we have and g(u) = A(u, f (u)) has the form We apply [11,Theorem 2.19] to System (2.35).Since the power of u in f (u) is even, we deduce from Part (iii) of [11,Theorem 2.19]: Lemma 2.6.If E * is a semi-hyperbolic equilibrium of (1.2) in the positive quadrant ]0, +∞[×]0, +∞[, and if κ 3 − ky * a + ak 1 κ = 0, then E * is a saddlenode, that is, its phase portrait is the union of one parabolic and two hyperbolic sectors.In this case, the index of E * is 0.
Case 2: if κ 3 − ky * a + ak 1 κ = 0, we have and g(u) = A(u, f (u)) has the form Again, we apply [11,Theorem 2.19] to System (2.35).Since the power of u in f (u) is odd, we look at the cofficient of u 3 and we have two possibilities: , then E * is a saddle.In this case, the index of E * is -1.
Remark 7. From Theorem 2.5, when the system (1.2) has one equilibrium point, this point cannot be a saddle.
Hopf bifurcation When ∆ χ < 0, the roots of χ are . The values of x * , y * and p do not depend on the parameter b, whereas s is an affine function of b, so that the eigenvalues of χ cross the imaginary axis at speed −1/2 when b passes through the value Let us check the genericity condition for Hopf bifurcations.We use the condition of Guckenheimer and Holmes [16,Formula (3.4.11)].Let us denote and V 1uv = ∂V1 ∂u∂v , etc.We have 1 κcy * 2 + 2k 2 1 y * 3 a 2 c.If λ < 0, then the periodic solutions are stable limit cycles, while if λ > 0, the periodic solutions are repelling.See Figure 3 for a numerical exemple.

c-Non-elementary equilibria Let us rewrite the vector field
For simplification, we denote Using the equality .
Since y * = x * + k 2 − m, we have also This shows in particular that the linear part of v is never zero.Thus the only non-hyperbolic cases are the nilpotent case and the case when E * is a center for the linear part of v. Let us now investigate these cases: This is when p = 0 = s.From the discussion at the beginning of Case b, it is clear that this case is nonempty.
In this case, the Jacobian matrix J (x * , y * ) is With the preceding notations, we thus have .

The case of a center of the linearized vector field
The point E * is a center of the linear part of v if the Jacobian J (x * , y * ) has purely imaginary eigenvalues ±i √ p, that is, when p > 0 and s = 0. Again, this case is nonempty.Let us denote With the notations of (2.36), we have p > 0 and s = 0 if, and only if, Note that x * , y * , as well as b 0 , a, ρ, and the sign of p do not depend on the parameter b, and that s = b − b 0 .Let us fix all parameters except b, and assume that ∆ χ < 0, that is, the eigenvalues of J (x * , y * ) are These eigenvalues cross the imaginary axis at speed −1/2 when b passes through the value b 0 .Let us denote c = aρ.By (2.37) and (2.38), we have Let us denote by (i, j) the standard basis of R 2 .In this basis, the matrix of the linear part ϕ of (X, The matrix of ϕ in the basis (u, v) is In particular, for b = b 0 ,

Existence of a globally asymptotically stable equilibrium point
When m = 0, in the case (c) of Theorem 2.
Proof.Let E * = (x * , y * ) ∈ A be an equilibrium point in the interior of A. Let us denote and let us set Then, using (2.20) and (2.21), we have For x ≥ m, a sufficient condition for V to be negative when (x, y) = (x * , y * ) is that g be nonincreasing.Let us make the change of variable X = x − m.We have which leads to Thus, if 2m + k 1 ≥ 1, g (X) remains negative for X > 0, i.e., for x > m.Thus, for x > m, under the assumption (2.42), V is negative.
We have seen that the first part of (2.42) implies that the equilibrium point E * , if it exists, is globally asymptotically stable.Note that Condition (2.42) is independent of the coordinates of E * , and the global stability implies that the equilibrium point E * , if it exists, is unique.
The second part of (2.42) is a necessary and sufficient condition for the existence of such an equilibrium point.
When m > 0, we already know that there exists at least one equilibrium point in A. Actually, Condition (2.42) implies that the coefficient α 2 = a + k 1 − 1 + 2m of (2.23) is positive.Thus, when m > 0, (2.42) is a particular case of (c) in Theorem 2.3.
When m = 0, by Theorem 2.4-(c), since α 2 > 0, there exists an equilibrium point in the interior of A if, and only if, (2.29) is satisfied.

Cycles
Let us investigate the existence of periodic orbits of (1.2).By Theorem 2.1 such orbits can take place only in A.
Proof.In the case (c), the only equilibrium points of (1.2) in R + × R + are the trivial points E 0 , E 1 , and E 2 , on the axes.Thus (1.2) has no cycle, because the compact set delimited by a cycle would contain a critical point, see [5,Theorem V.3.8].
In the case (a), if there was a cycle inside A, we could apply the Poincaré-Hopf Index Theorem to the compact manifold whose boundary is delineated by this cycle (see [26] for a version of this theorem when the vector field is tangent to the boundary).Denoting N the number of nodes or foci and S the number of saddles in the open quadrant ]0, +∞[×]0, +∞[, we would have N − S = 1.But Theorem 2.5 shows that N − S = 0, a contradiction.
In the case (b), if s < 0 and p > 0, the system (1.2) has an unstable equilibrium point.From Theorem 2.1 and Poincaré-Bendixson Theorem, there exists at least one limit cycle around this equilibrium.
Note that the conditions of Lemma 2.12 do not involve the value of b.Using Bendixson-Dulac criterion, M.A. Aziz-Alaoui and M. Daher-Okiye obtain another criterion: Lemma 2.13.[9, Theorem 7] if b + k 1 ≥ 1, then the system (1.2) has no limit cycle.

Case with refuge (m > 0)
By Theorem 2.11, if Condition (2.42) is satisfied, there can be no periodic orbits.
Let us now give some sufficient conditions for the absence of periodic orbits, using Bendixson-Dulac criterion.Let us denote by v 1 (x, y) and v 2 (x, y) the coordinates of the vector field in (1.2).For a Dulac function, we choose Let us look for conditions that ensure that ∂(v1D) ∂x + ∂(v2D) ∂y < 0 in A. We have For (x, y) ∈ A, we have Since the maximum of −3m 2 + m is 1/12 and the maximum of −m 2 + m is 1/4, we deduce: In particular, a condition that ensures that ∂(v1D ).
On the other hand, for (x, y) ∈ A, ∂(v2D) ∂y (x, y) has the same sign as x + k 2 − m − 2y, and we have The same technique does not provide any sufficient condition for Now, we consider the existence of limit cycles which are not occuring from a Hopf bifurcation.The special configuration of the existence of a limit cycle enclosing three equilibrium points is numerically investigated.In particular, when the system parameters satisfy a = 0.5, k1 = 0.08, k2 = 0.2, b = 0.1, m = 0.0025, then three hyperbolic equilibrium points exist, namely, E * 1 = (0.0222589; 0.2197589), E * 2 = (0.0299525; 0.2274525), E * 3 = (0.3702886; 0.5677886).They define respectively a stable focus, a saddle point and an unstable focus.Accordingly to the Poincaré index theorem, the sum of the corresponding indexes is equal to 1.
The numerical simulations show that there exists a limit cycle, which is hyperbolic and stable, see Figure 1.

Stochastic model
We now study the dynamics of the system (1.3), with initial conditions x 0 > 0 and y 0 > 0. In the case when m = 0 and k 1 = k 2 , the persistence and boundedness of solutions have been investigated in by Ji, Jiang and Shi in [17].A similar model has been studied by Fu, Jiang, Shi, Hayat and Alsaedi in [13].
Proof.Since the coefficients of (1.3) are locally Lipschitz, uniqueness of the solution until explosion time is guaranteed for any initial condition.
Let us now prove global existence of the solution.
The case when (x 0 , y 0 ) ∈ R + × {0} ∪ {0} × R + is trivial because both equations in (1.3) become independent, for example if y 0 = 0 with x 0 = 0, we have y(t) = 0 for all t ≥ 0, and x is a solution to the stochastic logistic equation which is well known (see Section 3.2), thus x(t) is defined for every t ≥ 0. Assume now that x 0 > 0 and y 0 > 0. Since the coordinate axes are stable by (1.3), we deduce, applying locally the comparaison theorem for SDEs (see [12, Theorem 1], this theorem is given for globally Lipschitz coefficients), that the solution to (1.3) remains in ]0, +∞[×]0, +∞[ until its explosion time.

Comparison results
In this section, we compare the dynamics of (1.3) with some simpler models, in view of applications to the long time behaviour of the solutions to (1.3).Applying locally the comparaison theorem for SDEs (see [12, Theorem 1], this theorem is given for globally Lipschitz coefficients), we have, for every t ≥ 0, (3.3) 0 ≤ x(t) ≤ u(t) a.s.
The process u is well known and can be written explicitely, see [19, page 125]:  s+σ2w2(s) ds .

Theorem 2 . 5 .
Assume that all equilibrium points of the system (1.2) which lie in the open quadrant ]0, +∞[×]0, +∞[ (equivalently, in the interior of A) are hyperbolic, and let n be their number.1.Assume that m > 0. Then n is equal to 3 or 1.

1 . 1 ∂∂x + v 2 ∂
Assume that m > 0. By Theorem 2.1, the vector field v = v ∂y generated by (1.2) is directed inward along the boundary of A. By continuity of v, we can round the corners of A and define a compact domain A ⊂ A with smooth boundary which contains all critical points of A, and such that v is directed inward along the boundary of A .Applying a generalized version of the Poincaré index theorem (see e.g.[23,15,31]) to v in A , we get N − S = 1.Since 1 ≤ N + S ≤ 3, the only possibilities are (N = 1 and S = 0) or (N = 2 and S = 1).

and v 2 <
0 for y = L, the field v is directed inward along the boundary of B. Again, by rounding the corners, we can modify B into a a compact domain B with smooth boundary which contains the same critical points as B and such that v is directed inward along the boundary of B .By the Poincaré Index Theorem, we have N − S = 1, where N (respectively S ) is the number of nodes or foci (respectively of saddles) in the interior of B .If we have chosen small enough, the singularities of v in B are those which are in the interior of A, with the addition of the point E 2 , which is a node by Proposition 1. Thus N = N − 1 and S = S which entails N − S = 0. Thus, taking into account Theorem 2.4, we have N = S = 1 (if n = 2), or N = S = 0 (if n = 0).

Figure 1 :
Figure 1: A phase portrait of (1.2) with three equilibrium points and a cycle in the interior of A. The dashed lines are isoclines y = x(1−x)(k1+x−m) a(x−m) and y = k 2 + x − m.The grey region is the invariant attracting domain A. m = 0.0025, a = 0.5, k 1 = 0.08, k 2 = 0.2, b = 0.1.
Remark 1.A more general sufficient condition of global attractivity of E 2 is provided by Theorem 2.4 (see Remark 3).