THREE LIMIT CYCLES IN A LESLIEâGOWER PREDATOR-PREY ...
ÁREA INFORMACIÓN Y DOCUMENTACIÓN SERVICIO DE REFERENCIA ELECTRÓNICA SISTEMA DE BIBLIOTECAS USM CONTACTO: [email protected] SIAM J. APPL. MATH. Vol. 69, No. 5, pp. 1244–1262
c 2009 Society for Industrial and Applied Mathematics
THREE LIMIT CYCLES IN A LESLIE–GOWER PREDATOR-PREY MODEL WITH ADDITIVE ALLEE EFFECT∗ ‡ , AND EDUARDO SAEZ † ´ ´ PABLO AGUIRRE† , EDUARDO GONZALEZ-OLIVARES
Abstract. In this work, a bidimensional continuous-time differential equations system is analyzed which is derived of Leslie-type predator-prey schemes by considering a nonmonotonic functional response and Allee effect on population prey. For the system obtained we describe the bifurcation diagram of limit cycles that appears in the first quadrant, the only quadrant of interest for the sake of realism. We show that, under certain conditions over the parameters, the system allows the existence of three limit cycles: The first two cycles are infinitesimal ones generated by Hopf bifurcation; the third one arises from a homoclinic bifurcation. Furthermore, we give conditions over the parameters such that the model allows long-term extinction or survival of both populations. In particular, the presence of a weak Allee effect does not imply extinction of populations necessarily for our model. Key words. stability, limit cycles, homoclinic orbits, bifurcations, predator-prey models, Allee effect AMS subject classifications. 92D25, 34C, 58F14, 58F21 DOI. 10.1137/070705210
1. Introduction. This work deals with a continuous predator-prey model considering the following: (i) the Allee effect [6, 13, 16, 34] affecting the prey population, (ii) the functional response of predators of nonmonotonic type, and (iii) a predator growth function of logistic type. Other inherent assumptions for the model are that population size varies only in time and it is uniformly distributed in space, there is no division of ages or sex, and it is not affected by abiotic factors. It is known that the Allee effect refers to a positive density dependence in prey population growth at low prey densities [34], and it occurs whenever fitness of an individual in a small or sparse population decreases as the population size or density also declines [6]. To be able to recognize the consequences of the Allee effect on reproduction, conservation and behavior of species has become an important aim over the last years. The analysis and understanding of this phenomenon can bring important benefits not only for ecology but also for various applied engineering disciplines such as agropecuary, fishing, and forestal industries. Different mechanisms generating Allee effects have been suggested (see Table 1 in [6]), being largely studied as singular entities and usually describing a situation in which the population growth rate decreases under some minimum critical density [5] or when a limited population growth capacity is observed [15]. In some cases this growth rate might be even negative, causing an extinction threshold [5]. In other words, the Allee effect may be understood as the cause of the increase in extinction risk at low densities [16], introducing in some cases a population threshold that has to be exceeded by population to be able to grow. This effect is also named in population dynamics as the negative competition effect [38]; in fisheries sciences, it is called depensation ∗ Received by the editors October 12, 2007; accepted for publication (in revised form) October 6, 2008; published electronically February 4, 2009. This work was partially supported by USM grant 12.08.04. http://www.siam.org/journals/siap/69-5/70521.html † Departamento de Matem´ atica, Universidad T´ecnica Federico Santa Mar´ıa, Casilla 110-V, Valpara´ıso, Chile ([email protected], [email protected]). ‡ Instituto de Matem´ aticas, Pontificia Universidad Cat´ olica de Valpara´ıso, Casilla 4059, Valpara´ıso, Chile ([email protected]).
1244
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
THREE LIMIT CYCLES IN A LESLIE–GOWER
1245
[9, 16], and in epidemiology, its analogous is the eradication threshold, the population level of susceptible individuals below which an infectious illness is eliminated from a population [5]. More precisely, the (component) Allee effect is any positive relationship between any measurable component of individual fitness and population size or density [6, 13, 34], and it might be the strong Allee effect [6] or critical depensation [9] that implies the existence of a threshold population level [5, 9], whereas the weak Allee effect [34] or noncritical depensation [9] does not have it. When the species is submitted to a strong Allee effect, it may have a bigger tendency to be less able to overcome these additional mortality causes, to have a slower recovery, and to be more prone to extinction than other species [34]. In most predation models it has been considered that the Allee effect has influence only on the prey population, and this effect is independent of the functional response or consumption rate that reflects the change on predation due to the prey’s population size. Quantitatively, it is assumed that the functional response has influence on the extension of the bistability region [15]. The Allee effect has been modeled in different ways using various mathematical tools and, in a first approach, like a deterministic phenomenon frequently associated to population’s stochastic fluctuations [5]. For instance, if x = x (t) indicates the population size, the most usual continuous growth equation to express the Allee effect is given by dx x =r 1− (x − m) x, dt K featuring the multiplicative Allee effect. Clearly, if m = 0, we have the weak Allee effect and if m > 0, it has the strong Allee effect. Other mathematical forms have been proposed to describe this phenomenon [8]. In this work, we consider the natural growth function deduced in [16] and [34] given by the equation dx x m (1) =r 1− − x dt K x+b that we call additive Allee effect. Recent ecological research suggests the possibility that two or more Allee effects generate mechanisms acting simultaneously on a single population (see Table 2 in [6]), and the combined influence of some of these phenomena has been named the multiple Allee effects. In the interacting populations, the predation can be largely reduced due to better ability of prey to avoid predation when their population size is large enough [42, 43]; but at low population densities, there could be a low effectiveness of antipredator vigilance, which reflects an Allee effect. For some marine species it has been shown that the per capita growth decreases as the size population is reduced below some critical level, and two proposed causes of depensation or Allee effect are the following: reduced breedings success at low densities of the population and increased relative predation on small populations [21]. This situation has happened in many real fisheries as a result of overfishing when man acts as predator [9]. Considering these aspects, (1) can be rewritten as dx rx x = 1− (x − xm )x, (2) dt x+b xK
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1246
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
where xm = 12 (K−b− 1r (r(K + b)2 − 4mK)) and xK = 12 (K−b+ 1r (r(K + b)2 − 4mK)) for r(K + b)2 − 4mK > 0. Hence, (2) represents two types of Allee effect affecting the same population since xm expresses the minimum of viable population and the r x factor r (x) = x+b indicates the impact of an Allee effect due to other causes affecting the intrinsic growth rate, for example, the predation reducing breeding success at low densities [9, 21]. On the other hand, the functional response of predators or consumption rate function refers to the change in the density of prey attached per unit time per predator as the prey density changes [20, 42]. In most predator-prey models considered in the ecological literature, the predator response to prey density is assumed to be monotonic increasing, the inherent assumption being that as there is more prey in the environment, it is better for the predator [20]. However, there is evidence that indicates that this need not always be the case, for instance, when a type of antipredator behavior (APB) exists. Group defense is one of these, and the term is used to describe the phenomenon whereby predators decrease, or are even prevented altogether, due to the increased ability of the prey to better defend or disguise themselves when their number is large enough [20, 42, 39, 43], and in this case a nonmonotonic functional response is better. For example, lone musk ox can be successfully attacked by wolves; however, large herds of them can be attacked but with rare success. Another manifestation of an APB in which a nonmonotonic functional response (or Holling-type IV or Monod Haldane) can be used is the phenomenon of aggregation, a social behavior of prey in which prey congregate on a fine scale relative to the predator so that the predator’s hunting is not spatially homogeneous [36], such as succeeds with mile-long schools of certain class of fishes. In this case, a primary advantage of schooling seems to be confusion of the predator when it attacks. The more important benefit of aggregation is an increasing in wariness. Moreover, aggregation can both decrease the vulnerability to be attacked and increase the time that group members can devote to activities other than surveillance [36]. Other related examples of nonmonotone consumption occur at the microbial level where evidence indicates that when faced with an overabundance of nutrient the effectiveness of the consumer can begin to decline. This is often seen when microorganisms are used for waste decomposition or for water purification, a phenomenon that is called inhibition [20, 42, 43]. In these cases, the functional response curves have an upper bound on the rate of predation per individual predator at some prey density, in contrast to the old Lotka– Volterra model which assumed a linear relationship between prey density and the rate of predation over the entire range of prey densities [36]. In this work, we use the function h(x) = x2qx +a , also employed in [28, 42, 39, 43] and corresponding to the qx Holling-type IV functional response [36], which is generalized as h(x) = x2 +bx+a in [44, 40] for a Gause model. This generalized expression is derived by Collings [12], who affirms that this type of functional response seems a reasonable possibility if it is assumed that prey and webbing densities are directly related. In [40], for the corresponding Gause model the existence of two limit cycles is proven. Moreover, this Gause model exhibits bifurcation of cusp type with codimension two [4, 23] or Bogdanov–Takens bifurcation [44, 40]. We note that the phenomena of the Allee effect and aggregation described by a nonmonotonic functional response are quite compatible and justify our assumptions in the model studied.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
THREE LIMIT CYCLES IN A LESLIE–GOWER
1247
The last aspect in our equations is a feature of Leslie-type predator-prey models [37] or the Leslie–Gower model [27], in which the conventional environmental carrying capacity Ky is proportional to prey abundance x [31], that is, Ky = nx, as in the May–Holling–Tanner model [4] and other models recently analyzed [28, 45]. In [37] it is affirmed that Leslie models can lead to anomalies in their predictions since they predict that even at very low prey density, when the consumption rate by an individual predator is essentially zero, predator population can nevertheless increase if predator population size is even smaller than prey population size. However, Leslie models are recently employed to model vole-weasel dynamics with a parameter time dependent [25], and the autonomous model proposed in [24] is analyzed in [45]. This scheme of modeling differs from a more common Gause-type model [40] in which the predator equation is based on the mass action principle, since the numerical response is dependent on functional response. Although it may seem that the two aspects considered in the model contradict each other since the prey population exhibits the Allee effect for low densities, while a nonmonotone functional response is suggested for the aggregation (group defense) when the prey population size is large, it is known that predation induces an Allee effect. Strikingly, a wide range of predator-driven Allee effects have been reported (see Table 2 in [22]); in particular, there is the case of the Atlantic cod (Gadus morhua) that forms schools during the day, since commercial fishing (man as predator) provokes stock collapse because a higher proportion of this aggregative population is caught per unit effort when population declines [10]. Also, for obligately cooperative breeders as the African wild dog (Lycaon pictus) and meerkat (Suricata suricatta), there is a similar situation, because juvenile survival is lower in small groups than large groups in areas with high predator densities but lower in large groups than small groups in areas with low predator densities [14, 22]. From a mathematical point of view, simple models for the Allee effect may reveal a lot about its dynamics, and, reciprocally, different nonequivalent dynamics for the same model will have different biological interpretations. For instance, in our model we obtain the existence of a subset of parameter values for which three limit cycles appear but surrounding different equilibrium points. These limit cycles are not only periodic solutions of the system but also along the attractive behavior of the origin, allowing the phenomenon of multistability of our predator-prey system, that is, the existence of four ω-limit sets in the first quadrant. This paper is structured as follows: The model and the main results are presented in section 2. In section 3 we give the proofs of the main statements. Finally, an interpretation of the results is given in section 4, complemented with some numerical simulations. 2. The model. Let us consider the bidimensional system of ordinary differential equations
m x − xqxy x˙ = r 1 − xk − x+b 2 +a , Xμ : (3)
y y˙ = sy 1 − nx , where (x, y) ∈ A = {(x, y)| x > 0, y ≥ 0} and μ = (r, a, b, k, m, n, q, s) ∈ R8+ . In system (3), x(t) and y(t) denote prey and predator densities, respectively, as functions of time t. Furthermore, the parameters have the following meanings: (a) r and s are the intrinsic growth rates or biotic potential of the prey and predators, respectively.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1248
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
(b) q is the maximal predator per capita consumption rate, i.e., the maximum number of prey that can be eaten by a predator in each time unit. (c) a is the number of prey necessary to achieve one-half of the maximum rate q. (d) n is a measure of the food quality that the prey provides for conversion into predator births. (e) k is the prey environment carrying capacity [12]. (f) m and b are constants that indicate the severity of the Allee effect that has been modeled. In order to describe the dynamics of Xμ we will consider a C ∞ -equivalent polynomial extension by the following change of coordinates: u(u + b)(u2 + a) τ + 2 = (x, y, t), ζ : R2 ×R+ −→ R ×R such that ζ(u, v, τ ) = u, nv, 0 0 s (4) 2 +a) where detDζ(u, v, τ ) = nu(u+b)(u > 0. s Also, let us consider new parameters given by ϕ1 : R8+ −→ R × R7+ , m nq br − m , a, b, k, , n, , s := (L, a, b, k, M, n, c, s), ϕ1 (r, a, b, k, m, n, q, s) = s s s 3
s with Jacobian detDϕ1 (r, a, b, k, m, n, q, s) = bn > 0, which says that ϕ1 is invertible. For simplicity, let us rename the new coordinates u → x, v → y. Then in the new −1 coordinates we have the vector field Yη = ζ∗ Xμ = (Dζ) ◦ Xμ ◦ ζ, where ⎧
⎨x˙ = x (a + x2 ) (b+x)(L+M ) 1 − x − M x − cxy(b + x) , b k Yη : (5) ⎩ 2 y˙ = (b + x)(a + x )(x − y)y,
¯ = {(x, y)|x, y ≥ 0} and the new vector of parameters with (x, y) ∈ Ω η = (L, a, b, k, M, c) ∈ D0 is given by the natural projection, where D0 = η ∈ R × R5+ | M + L > 0 . (6) Notice that vector field (5) is a polynomial extension of the original system (3) ¯ including axis x = 0. Moreover, Yη and Xμ are clearly to the whole first quadrant Ω ∞ C -equivalent in Ω. First, we study the local behavior of singularities in the coordinate axis, that is, in absence of prey and predator, respectively. In system (5), it is immediate that ∂ Yη (0, y) = −aby 2 ∂y . Then the x = 0 axis in Ω is an invariant manifold of vector field (5) and the origin is its only singularity, being an attracting one. Moreover, if we ∂ ∂ ∂ denote Yη (x, y) = P (x, y) ∂x + Q(x, y) ∂y , we have Yη (x, 0) = P (x, 0) ∂x . Hence, the y = 0 axis in Ω is also an invariant manifold of (5). For vector field Yη let us define Sing(E) and Sing(Ω) as the sets of singularities of (5) in the y = 0 axis in Ω and in the open set Ω, respectively. In order to describe the dynamics of singularities in Sing(E), we consider the following regions in parameter space D0 : Λ1 Λ2 Λ4 Θ+ Γ+
= Δ−1 (−∞, 0) ∪ Δ−1 (0, ∞) ∩ Θ−1 (0, ∞) ∩ Υ−1 (−∞, 0) ; −1 −1 −1 = Δ (0, ∞) ∩ Θ (∞, 0) ∩ Υ (−∞, 0); = {η ∈ D0 | 0 < b ≤ L}; = Θ−1 [0, ∞) ∩ Υ−1 (0); = Γ−1 (0, ∞) ∩ Υ−1 (0);
Λ3 D Θ− Γ−
= {η ∈ D0 | 0 < L < b}; = Δ−1 (0) ∩ Θ−1 (0, ∞); = Θ−1 (−∞, 0) ∩ Υ−1 (0); = Γ−1 (−∞, 0] ∩ Υ−1 (0);
(7)
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
THREE LIMIT CYCLES IN A LESLIE–GOWER
b
1249
y
Λ1 y
• P0/+
• P0
D
L = −M
Λ2
y
x
b = L b = k •
y
• • P0 P±
Λ3
x
x y
y
• • P0/−P+
• • • x P0P− P+ L = 0
• • P0 P+
x
Λ4
x y
• • P0 P+
x
L
Fig. 1. Bifurcation diagram in plane Lb = {η ∈ D0 | a, k, M, c are constant} in parameter space for singularities in the y = 0 axis.
where Υ(η) = sign(L), Θ(η) = b − k, and Δ(η) = (M + L)(k − b)2 + 4bkL, √ √ √ Γ(η) = −2kcb + a M −a M + aM + 4k 2 c . Moreover, in the xy plane let us consider the points P0 = (0, 0), P+ = (x+ , 0), and P− = (x− , 0), with (M + L)(k − b) ± (M + L)Δ(η) (8) . x± = 2(M + L) Lemma 1. For vector field Yη the origin P0 is a non-hyperbolic singularity and it satisfies the following: (i) It has two hyperbolic sectors divided by a repulsing separatrix if η ∈ Γ+ ∩ Θ+ . (ii) It has only a hyperbolic sector if η ∈ Λ4 . (iii) It has a hyperbolic sector and a repulsing parabolic sector if η ∈ Λ3 ∪(Θ− ∩ Γ+ ) . (iv) It has a hyperbolic sector and an attracting parabolic sector if η ∈ Θ− ∩ Γ− . (v) It is a local attractor if η ∈ D ∪ Λ1 ∪ Λ2 ∪ (Θ+ ∩ Λ− ) . A qualitative diagram of these results in plane Lb = {η ∈ D0 | a, k, M, c are constant} is shown in Figure 1. Lemma 2. For vector field Yη the following statements hold: (i) P+ ∈ Sing(E) if η ∈ D ∪ Θ− ∪ Λ3 ∪ Λ4 . Moreover, P+ is a saddle node if η ∈ D; otherwise, it is a hyperbolic saddle. (ii) P− ∈ Sing(E) if η ∈ Λ2 . Moreover, P− is a hyperbolic repulsing node. A qualitative diagram of these results in plane Lb = {η ∈ D0 | a, k, M, c are constant} is shown in Figure 1. About the boundness of the system, in [1] it is proven that no trajectory of vector field Yη has the infinity as ω-limit. Since the coordinate axes are invariant, then our system is bounded. We now state the same result but with a different approach. Theorem 3. Consider Aw = {(x, y) ∈ R2 : 0 ≤ x ≤ w, 0 ≤ y}. For every η ∈ D0 , there is a w∗ ≥ 0 such that if w > w∗ , then Aw is a trapping domain for system (5), meaning that it is invariant for positive time evolution and also captures all trajectories starting in Ω.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1250
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
We take care now of the dynamics at the interior of the first quadrant. It is clear that every singularity of system (5) in the open quadrant Ω must be over the diagonal y = x. Furthermore, over the isocline y = x the vector field Yη has the form ∂ Yη (x, x) = x2 p4 (x) ∂x , where p4 (x) is a polynomial of degree 4 in x. Therefore, the existence of equilibrium points in the open first quadrant is given by the existence of real and positive roots of the polynomial p4 (x). In order to describe the bifurcation diagram for an equilibrium point of Yη in Ω, we change parameters η = (L, a, b, k, M, c) to get a new vector ξ = (L, a, b, k, V, q) by means of the transformation ϕ2 : R6 −→ R6 , (L, a, b, k, V, q) → (L, a, b, k, M (ξ), c(ξ)),
(9) with c(ξ) = M (ξ) =
aLq 4 (b−k+2q)−q 6 (2b2 +kL−b(2k+L−4q)−2(k+L)q+2q 2 ) kq 6 −bq 4 (b+q)6 (a+q 2 )2 (a+3q 2 )(aL+(L−2b)q 2 )+V . bq 6 (b+q)4 (a+q 2 )3 (a+3q 2 )
−
(b−k+q)V kq 6 b2 (b+q)5 (a+q 2 )2 (a+3q 2 ) ,
Let D1 = ϕ−1 2 (D0 ) be the new admissible region in parameter space, and let us consider the projection π : R6 −→ R3 such that
π(L, a, b, k, q, V ) = (L, b, k).
If we name
D2 = π D1 ∩ {(L, a, b, k, q, V )| a = 3b2 , q = b, V > 0} ∩ E −1 (0, ∞) ⊂ R × R2+ , with E(L, b, k) = −18432b19 L + V, it is easy to check that D2 = F −1 (0, ∞) ∩ E −1 (0, ∞), with F (L, b, k) = 12288b19 (b − 2L) + V. Changing the time t → bkt, let us call Yξ to the new vector field qualitatively equivalent to (5). Now, it is straightforward to check that there is a singularity in the isocline y = x at the point with coordinates (b, b). In order to describe the bifurcation diagram of this singularity of vector field Yξ in parameter space (L, b, k), we consider the projection (see Figure 2) 2 −→ R2 such that πR : SR
πR (L, b, k) = (L, b),
2 = {(L, b, k)| L2 + b2 + k 2 = R2 , b, k, R > 0}. Since function F (L, b, k) does where SR not depend on k, it is clear that in D2 ⊂ R × R2+ the surface F −1 (0) is a straight 2 half-cylinder with axis parallel to the k axis, so πR F −1 (0) ∩ SR is exactly the curve
2 , F = {(L, b)| F (L, b, k) = 0, k > 0} ∩ πR D2 ∩ SR
which is located in the interior of the first quadrant of plane Lb because equaV tion F (L, b, k) = 0 defines implicitly L = 2b + 24576b 19 > 0. Furthermore, (0, 0) ∈ −1 F (0, ∞), so the origin is in the admissible region for parameters in the upper
half2 . With plane Lb. The same arguments follow to sketch the curve πR E −1 (0) ∩ SR 2 this, the qualitative shape of domain πR (D2 ∩ SR ) is shown in Figure 3(a).
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1251
THREE LIMIT CYCLES IN A LESLIE–GOWER k 2 SR
L
b πR b = R k = 0
L = −R
L = R
2 to the plane Lb. Fig. 2. Projection πR from the sphere SR
2 πR D2 ∩ SR
b
Σ0 2
2 πR E −1 (0) ∩ SR
b p2 •
2 πR F −1 (0) ∩ SR
• p3 •
p1 •
p4
L
(a)
H
L
(b)
2 . (b) Qualitative shape of Fig. 3. (a) Qualitative shape of parameters domain πR D2 ∩ SR 2 curves H and Σ0 (R) in parameters domain.
2 in plane Lb, we define the sets Finally, in πR D2 ∩ SR
sg() 2 Σ1 (R) = πR {(L, b, k) ∈ D2 | sg(T (L, b, K)) = sg( )} ∩ SR ,
sg() 2 Σ2 (R) = πR {(L, b, k) ∈ D2 | sg(G(L, b, K)) = sg( )} ∩ SR , H = Σ01 (R),
(10)
with T (L, b, k) = −36864b21 + 55296b20 L − 18432b19 kL − 3bV + kV, G(L, b, k) = 1811939328b39 (b + L)(b + 2L) − 8bV 2 , and the qualitative shape of curves H and Σ02 (R), in terms of R, is as in Figure 3(b). Moreover, H ∩ Σ02 (R) consists of four points as in Figure 3(b), that is, (11)
H ∩ Σ02 (R) = {pi = (Li , bi ); i = 1, . . . , 4}.
Lemma 4. The singularity (b, b) of vector field Yξ is an order two weak focus if (L, b) ∈ H ∩ Σ02 (R) = {pi = (Li , bi ); i = 1, . . . , 4}. Moreover, the focus is repellor at p1 and p2 and attractor at p3 and p4 . Theorem 5. (i) The system Yξ has a bifurcation diagram as indicated in Figures 4(a)–(d) in a neighborhood of H ∩ Σ02 (R). (ii) The bifurcation curve S34 of semistable limit cycle of the system Yξ in the parameter space Σ− 1 (R) is as indicated in Figure 4(d). (iii) The bifurcation curve HC of a homoclinic loop surrounding the focus (b, b) of system Yξ is as indicated in Figure 4(a), in a tubular neighborhood of curve H in sector Σ− 2 (R) in parameter space. Lemma 6. In parameter space ξ = (L, a, b, k, q, V ) ∈ D1 of system Yξ , there exists −1 an open subset Υ ⊂ π −1 ◦ πR (p3 ) ∩{ξ | L < 0}, where the following statements hold:
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
1252
(a) S2 b
• p2 HC
• p3
H S1
S34 • p1
• p4 L
(b)
(c) H
•u0 •u1 S1
•u0
•u1
H
S2 • p2
• p1
•u0
•u0
•u2
•s0 •u2 •s1
•s1
•s0
•s1
(d) •s2
• p3
•u0 •u1 •s2 •s0
•s0
•s1
H S34
• p4
Fig. 4. (a) Bifurcation curves for focus (b, b) of vector field Yξ , where H is a Hopf bifurcation curve, S1 , S2 , and S34 are semistable limit cycle bifurcation curves, and HC is a Homoclinic bifurcation curve. (b) Bifurcation diagram for focus (b, b) in a neighborhood of p1 ∈ H. (c) Bifurcation diagram for focus (b, b) in a neighborhood of p2 ∈ H. (d) Bifurcation diagram for focus (b, b) in a tubular neighborhood of curve H between p3 and p4 . Notation •ij means a focus with local stability i ∈ {s, u} and weakness j ∈ {0, 1, 2}.
(i) Yξ has two foci (b, b), (xf , xf ) and two saddle points (xs , xs ), (xN , xN ) in the interior of the first quadrant Ω, where 0 < xs < b < xN < xf . (ii) There exists an open subset O ∈ R+ and a continuous function H : O −→ R+ such that if V = H(b), then there is a homoclinic loop L0 for vector field Yξ passing through saddle (xN , xN ) and surrounding focus (xf , xf ), with (xN , xN ) and (xf , xf ) as in part (i). Moreover, this homoclinic orbit is inner unstable. Lemma 7. Let Υ, H, and L0 be as in Lemma 6. There exist δ > 0 and ε > 0 such that, for every b, V with |H(b) − V | < ε, the vector field Yξ has, at most, one limit cycle L in the inner δ-neighborhood of L0 , and L is unstable. Theorem 8. There exists an open subset in parameter space such that vector field Yη in (5) has three hyperbolic limit cycles in the interior of the first quadrant Ω. Two of them are infinitesimal cycles surrounding a hyperbolic attracting focus (b, b); moreover, the outermost limit cycle is stable and the innermost unstable. The third limit cycle is unstable and surrounds a hyperbolic attracting focus (xf , xf ), with 0 < b < xf . 3. Proofs of the main results. The proofs of Lemmas 1 and 2 are straightforward and follow from the blowing-up method [17] for (0, 0) and the central manifold and Hartman’s theorems for P+ and P− ; for details, see [1].
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
THREE LIMIT CYCLES IN A LESLIE–GOWER •
1253
u
v
fw (x, y) = x − w = 0 Aw y x
• x = w
−1 (−∞, 0) ∩ Ω is an invariant set under Φ Fig. 5. The set Aw = fw Yη (t; · ), as stated in the proof of Theorem 3.
Proof of Theorem 3. We will prove that for every η ∈ D0 (see (6) for definition), there exists some w = w(η) > 0 such that the integral curves ΦYη (t; · ) of vector field Yη , for t > 0, enter an invariant subset Aw ⊂ Ω and do not leave again. Let us define the function fw : R2 −→ R given by fw (x, y) = x − w. For each w > 0 consider the line Lw = fw−1 (0) ∩ Ω. If we consider the Poincar´e compactification of system (5) given by + 2 Ψy : R2 × R+ 0 −→ R × R0 , Ψy (u, v, τ ) = (u/v, 1/v, τ ) = (x, y, t),
it is straightforward to check that the line Ψ−1 y (Lw ) contains the point (u, v) = (0, 0). Hence, we will show that for w sufficiently large, in the Poincar´e compactification of system (5), the compact set Aw = fw−1 (−∞, 0) ∩ Ω is an invariant set under ΦYη (t; · ) (see Figure 5). Recall that the axes x = 0 and y = 0 are invariant under ΦYη (t; · ), so it is sufficient to prove that the scalar product Yη · ∇fw < 0 in Lw for appropriate w > 0 as in Figure 5. We have that for (x, y) ∈ Lw
w2 a + w2 2 Yη · ∇fw |x=w = −cw (b + w)y + φ(w), bk where φ(w) = −(M + L)w2 + (k − b)(M + L)w + bkL. The discriminant of φ as a quadratic polynomial in w is (M + L)Δ(η) (see (7) for definition). Therefore, since M + L > 0 in D0 , if Δ < 0, then Yη · ∇fu |x=w < 0 for every w > 0. On the other side, if Δ ≥ 0, the roots of φ(w) are x− and x+ (see (8) for definition); hence, the statement holds for every w > x+ if x+ > 0 and for every w > 0 if x+ ≤ 0. This concludes the proof of Theorem 3. 2 Proof of Lemma 4. Let (L, b) ∈ πR (D2 ∩ SR ). In order to study the singularity (b, b) in vector field Yξ , we translate (b, b) to the origin and change the time by means of + 2 ζ : R2 × R + 0 −→ R × R0 , 8b3 E(L, b, k) t , (x, y, t) → x + b, y + b, 768b15
where E(L, b, k) = −18432b19 L + V.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1254
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
So if (L, b, k) ∈ E −1 (0, ∞), ζ is a C ∞ -equivalence. Let Zξ = ζ∗ Yξ ; then we have trDZξ (0, 0) =
8E(L, b, k)T (L, b, k) , 589824b25
with
T (L, b, k) = − 36864b21 − 55296b20 L + 18432b19 kL + 3bV − kV . 2 Moreover, if (L, b) ∈ H = πR (T −1 (0) ∩ SR ), we have detDZξ (0, 0) = V > 0, so the origin is a weak focus of order, at least, one for Zξ . In order to recognize the topological type of the origin as (L, b) ∈ H, we bring Zξ onto its Jordan canonical + 2 form at (0, 0). For that, consider the C ∞ -equivalence ϕ : R2 × R+ 0 −→ R × R0 , √ √ ϕ(u, v, τ ) = 192b10 u − V v, 192b10 u, τ / V = (x, y, t).
The qualitatively equivalent vector field in the new coordinates is ZξJ = ϕ∗ Zξ , and we have in a neighborhood of the origin ⎛ ⎛ ⎞ ⎞ 5 5 ∂ ∂ + ⎝u + , ZξJ (u, v) = ⎝−v + Ai,j ui v j + H.O.T.⎠ Bi,j ui v j + H.O.T.⎠ ∂u ∂v i,j=2 i,j=2 where H.O.T. denotes the higher order term and Ai,j = Ai,j (b, V, L), Bi,j = Bi,j (b, V, L). For j = 0, 1, 2, let lj be the first three Lyapunov quantities at the origin [2, 7, 33] of the vector field ZξJ . Since the trace of the linear part of the vector field at the origin is zero, we have that l0 = 0. Now, as l1 depends on the 3-jet of ZξJ , see [17], then l1 = (A02 A11 + A12 + A11 A20 + 3A30 + 2A02 B02 + 3B03 − B02 B11 − 2A20 B20 − B11 B20 + B21 )/8.
Using the Mathematica software [41] we have l1 =
36864b20 + V G(L, b, k), 8192b13 V 3/2
with G(L, b, k) = 1811939328b39 (b + L)(b + 2L) − 8bV 2 . Thereby, if (L, b) ∈ H ∩ Σ02 (R) (see (10) and (11)), the weakness of the origin depends only on l2 . On the other hand, it is known that l2 depends on the 5-jet of ZξJ . Then, again using the Mathematica software, as l0 = l1 = 0, we obtain l2± =
36864b20 + V N± , 98304b15 V 3/2
where
61 N± = − √−221259535220736b − 13589544960b41 V −√1963008b21 V 2 − 57bV 3 40 ± 3 30349983744b + 2211840b20 V + 43V 2 28311552b42 + b2 V 2
and the sign ± denotes the branch (12)
√ −27648b39 ± 3 g(b) L± = 36864b38
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1255
THREE LIMIT CYCLES IN A LESLIE–GOWER
0 -0.2 2 4 N+ -0.4 -0.6 6 0
0 0.2 b
0.2 V 0.4
0.4
0 -20 0 0 0.01∗N+ -40 -60 0 0.6 0.
20 40 V 60
0.8b 80 1001
(a)
(b)
Fig. 6. Graphics of the function N+ = N+ (b, V ) where we can see that N+ < 0: (a) graphic of N+ for (b, V ) ∈ [0, 0.5] × [0, 0.5]; (b) graphic of 10−2 × N+ for (b, V ) ∈ [0.6, 1] × [10, 100].
of curve Σ02 (R) given implicitly by equation G(L, b, k) = 0, where g(b) = 28311552b78 + b38 V 2 > 0. Since the branch L− intersects curve H at points p1 and p2 (see Figure 3(b)), then if (L, b) ∈ {p1 , p2 } ⊂ H ∩ Σ02 (R), l2 > 0 and the origin is a repellor weak focus of order two for ZξJ . On the other hand, the branch L+ intersects curve H at points p3 and p4 (see Figure 3(b)), and by looking at the graphic of N+ as a function of b and V it can be verified that N+ < 0. Figure 6(a) shows the graphic of N+ for the range (b, V ) ∈ [0, 0.5] × [0, 0.5], meanwhile in Figure 6(b) there is the graphic of 10−2 × N+ for (b, V ) ∈ [0.6, 1] × [10, 100], made with the software Mathematica [41]; here it can be seen that if b > 0, V > 0, as in our case, then by continuity N+ < 0. Therefore, if (L, b) ∈ {p3 , p4 } ⊂ H ∩ Σ02 (R), l2 < 0 and the origin is an attracting weak focus of order two for ZξJ . 2 Proof of Theorem 5(i). It is clear that T −1 (0) and SR are transversal manifolds in the parameter space and, by Lemma 4, there are four points
2 p1 , p2 , p3 , p4 ∈ πR T −1 (0) ∩ SR = H, where H is a Hopf bifurcation curve and where the singularity (b, b) is a repellor weak focus of order 2 of system Yξ at p1 and p2 and an attractor weak focus of order two at p3 and p4 ; i.e., the eigenvalues of DYξ (b, b) are on the imaginary axes and nonzero. Moreover, p1 , p2 , p3 , and p4 have codimension two, since the Lyapunov quantities l0 = l1 = 0 and l2 = 0. Thereby, let (L, b) ∈ H − {p1 , p2 , p3 , p4 }. Then we have the following (see Figure 7): (a) In a neighborhood of p4 , when (L, b) is below l1−1 (0), the singularity (b, b) of system Yξ is an order one attractor weak focus, because l1 < 0. Therefore, if > 0 is sufficiently small, the point (L + , b + ) is over the curve H; then the stability of system Yξ at the singularity (b, b) is reversed; hence, a unique hyperbolic stable limit cycle bifurcates (Hopf bifurcation). (b) When (L, b) is above l1−1 (0) in a neighborhood of p4 , the singularity (b, b) is an order one repellor weak focus surrounded by a hyperbolic stable limit cycle, because l1 > 0; then Hopf bifurcation from p4 occurs. Therefore, if > 0 is sufficiently small, the point (L + , b) is over the curve H; hence, the stability of system Yξ at the singularity (b, b) does not change and the limit cycle persists because it is hyperbolic. (c) When the point (L, b) is as in case (b), for > 0 sufficiently small, the point (L − , b) is under the curve H; then the stability of system Yξ at the singularity (b, b) is reversed; hence, a new Hopf bifurcation occurs and a
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1256
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ b
−1 Σ0 2 (R); l1 (0) p4 • −1 H; l0 (0) L
Fig. 7. Hopf bifurcation curve H near point p4 in parameter space, where the singularity (b, b) of vector field Yξ is an order two attracting weak focus.
second hyperbolic limit cycle bifurcates, which is surrounded by a hyperbolic stable limit cycle. By normal forms theory [17, 23], p4 has a 2-parameter versal unfolding. Under the same hypothesis that we have here, Takens [35] and Arrowsmith and Place [3] describe in detail the bifurcation diagram for codimension two singularity type. By using their results, we have that there exists a neighborhood Vδ (p4 ), with δ > 0, such that it has a diagram as in Figure 4, where S34 is a curve in which the unstable and stable limit cycles collapse (semistable limit cycle). The analysis near points p1 , p2 , and p3 follows in the same way. Proof of Theorem 5(ii). From part (i) of this Theorem, there is an open set C(R) ⊂ Σ− 1 (R) where system Yξ has two hyperbolic limit cycles surrounding the singularity (b, b). Moreover, C(R) is upperly bounded by the segment of curve H between p3 and p4 , because at any point in that segment focus (b, b) is surrounded by a unique infinitesimal limit cycle. On the other hand, at Σ− 1 (R)\ C(R), the singularity (b, b) has no infinitesimal limit cycle in its vicinity, so C(R) is bounded and the lower bound must correspond to a bifurcation curve where both limit cycles collapse into a semistable limit cycle; hence, by versality of the unfolding of the singularity (b, b) at p3 and p4 described in part (i) (see [35, 3]), it is clear that the semistable limit cycle bifurcation curve S34 must be as in Figure 4(d). Proof of Theorem 5(iii). By Hartman–Grobman’s theorem it is straightforward to see that vector field Yξ has a hyperbolic saddle singularity Ps = (s, s) with 0 < s < b if (L, b) ∈ H ∩ Σ− 2 (R). Let (L, b) ∈ H ∩Σ− 2 (R) in a neighborhood of p3 . By part (i), (b, b) is an attracting weak focus of order one with no limit cycle in its vicinity, and it can be checked that Ps ∈ α−lim(b, b). As we move along segment H ∩ Σ− 2 (R) towards p2 , the topological type of focus (b, b) remains the same according part (i); however, if (L, b) ∈ H ∩Σ− 2 (R) c in a neighborhood of p2 , Ps ∈ (α−lim(b, b)) and an infinitesimal limit cycle exists surrounding (b, b). Therefore, there exist δ > 0 sufficiently small and a differentiable function Ψ : Vδ (H) −→ R
such that (L, b) → Ψ(L, b),
defined in a tubular neighborhood Vδ (H) of radius δ of segment H ∩ Σ− 2 (R) such that vector field Yξ has a homoclinic loop passing through the saddle Ps and surrounding focus (b, b), if (L, b) ∈ HC = Ψ−1 (0). Remark. An extended and deeper proof of the existence of a homoclinic loop surrounding another singularity is given in Lemma 6. Nevertheless, the same arguments might have been used in this case as well.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1257
THREE LIMIT CYCLES IN A LESLIE–GOWER
Σ
Σ Φτ (ˆ r) p ˆ •
W s (N )
p ˆ •
ˆ • r • W s (N )
Φτ (ˆ r)
• p∗
•
• r ˆ •
W u (N )
• p∗
•
N
W u (N )
N
(a)
(b)
Fig. 8. Poincar´ e map for orbits near W u (N ) ∩ W s (xf ) = ∅.
W s (N )
in two cases: (a) W u (N ) ∩ W s (xf ) = ∅; (b)
Proof of Lemma 6(i). Let (L, b) = p3 ∈ H ∩ Σ02 . It is easy to check that Yξ has four singularities in isocline y = x at Ω, if L < 0. According to Lemma 4, one of these equilibria is an attracting weak focus of order two in (b, b). Due to Hartman–Grobman’s theorem, it is straightforward to see that the other singularities are a hyperbolic attracting focus in (xf , xf ) and two saddles (xs , xs ), (xN , xN ) with 0 < xs < b < xN < xf . The statement follows directly. Proof of Lemma 6(ii). Let (xN , xN ) and (xf , xf ) be as in the proof of part (i). If we name N = (xN , xN ), let us consider the function D(N ) = ∇ · Yξ (N ), and let Φt = ΦYξ (t; · ) be the flow of Yξ and p∗ = (x∗ , y ∗ ) be a point in W u (N ) with xN < x∗ . Let Σ be a one-dimensional local cross section to W s (N ) ∪ W u (N ) as in Figure 8. Let pˆ = (ˆ p1 , pˆ2 ) the unique point of Σ ∩ W s (N ) which never returns to Σ by the flow Φt with 0 < t < ∞, and let U ⊂ Σ be some neighborhood of pˆ. In Σ let us choose local coordinates given by the chart h : Σ −→ R, h(x, y) = y − pˆ2 . Let us consider the Poincar´e map P : U −→ Σ defined for a point r ∈ U by P(r) = Φτ (r), where τ = τ (r) is the time taken for the orbit Φτ (r) based at r to first return to Σ. Now consider the displacement function d : U −→ R given by r → d(r) = P(r) − r. Finally, let rˆ be a point in U located below pˆ as in Figure 8. Let b = 1 and V = 300. Then it can be seen by means of numerical simulations on Matlab [30] that D(N ) > 0, W u (N ) ∩ W s (xf ) = ∅, W u (N ) ∩ W s (N ) = ∅, and ω − lim Φt (p∗ ) = {(xf , xf )}. In the local coordinates given by h this implies that (13)
d(ˆ r) < 0.
A qualitatively equivalent situation is shown in Figure 8(a). Instead, if b = 1 and V = 400, it can be checked (numerically as well) that D(N ) > 0, W u (N ) ∩ W s (xf ) = ∅, and W u (N ) ∩ W s (N ) = ∅; there exists y † > xN such that the point (xN , y † ) ∈ Φt (p∗ ) ⊂ W u (N ). Therefore, we have (see Figure 8(b)) (14)
d(ˆ r) > 0.
Hence, by the continuous dependence of the orbits of the vector field on parameters, there exists V ∗ , 300 < V ∗ < 400, such that for b = 1 and V = V ∗ we have W u (N ) ∩ W s (N ) = ∅; that is, there is a homoclinic loop L0 ⊂ W u (N ) ∩ W s (N ).
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1258
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
L0 • • V
V ∗•
(1, V ∗ ) •
• • • 1
• F−1 (0) •
Bδ
• b
Fig. 9. Existence of a homoclinic orbit along the curve F−1 (0) in plane bV .
Therefore, there exists δ > 0 such that in a δ-ball Bδ (1, V ∗ ) ⊂ R2+ of the point (1, V ∗ ) in plane bV in parameter space, there is a locally differentiable function F : Bδ (1, V ∗ ) −→ R, (b, V ) → F(b, V ), such that if (b, V ) ∈ F−1 (0), there exists an inner unstable homoclinic loop L0 for vector field Yξ (see Figure 9). The existence of the function H : O −→ R+ in the statement follows from the implicit function theorem, due to the transversality shown in passing from (13) to (14). Proof of Lemma 7. Let the function d : U −→ R, the saddle N = (xN , xN ), and the point rˆ be as in the proof of Lemma 6 (see Figure 8). From the previous lemma, a homoclinic bifurcation (see [29]) occurs as V = H(b) since, by continuity, D(N ) = ∇ · Yξ (N ) > 0, N is a strong saddle [29], and L0 is inner unstable. Then, for U sufficiently small, by carrying out a perturbation on parameters, such that d(ˆ r) > 0, the homoclinic cycle breaks out and an unstable limit cycle appears surrounding the attracting focus (xf , xf ). Proof of Theorem 8. From Lemmas 4 and 7, in parameter space ξ = (L, a, b, k, q, V ) −1 ∈ D1 of system Yξ , there exists an open subset ΥL ⊂ π −1 ◦ πR (p3 ) ∩ {ξ | L < 0}, where system Yξ has four singularities in the interior of the first quadrant Ω: two saddles (xs , xs ), (xN , xN ), a two-order attracting weak focus (b, b), and a hyperbolic attracting focus (xf , xf ) surrounded by an unstable limit cycle L, with 0 < xs < b < xN < xf . Moreover, this limit cycle is hyperbolic, so it persists in an open subset of −1 π −1 (πR (C(R))) in parameter space, where C(R) is the bounded sector in plane Lb (see Figure 4(d)) where system Yξ also has two infinitesimal limit cycles surrounding focus (b, b), due to the proof of Theorem 5(ii). The statement of the theorem is immediate after recalling that vector field Yξ is induced by a vector field C ∞ -equivalent to Yη by means of the local difeomorphism ϕ2 in (9), since detDϕ2 (ξ) =
b3 q 12 (b
+
V > 0, + q 2 )4 (a + 3q 2 )3
q)10 (a
where ξ ∈ D1 ∩ π −1 (T −1 (0)). 4. Discussion. In this paper, a predator-prey model with the Allee effect [5, 16, 34, 38] and a Holling-type IV [36] functional response was considered, making a qualitative analysis of a bidimensional system of ordinary differential equations of polynomial type, which is topologically equivalent to the original one, with only six parameters. Employing a reparameterization and a time change, we showed that b, k, and L = br−m are the most significative parameters of this predator-prey model. s We note that this model predicts that, when L < 0, the system has three equilibria
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
THREE LIMIT CYCLES IN A LESLIE–GOWER
1259
P0 , P− , and P+ for which the predator density is null; moreover, at low population densities, both predator and prey populations may disappear, since the equilibrium P0 = (0, 0) is an attractor; furthermore, P− is a repelling node, and P+ is a saddle. Unlike other models, in this work we show that the weak Allee effect does not imply extinction of species necessarily for certain parameter values, which, as far as we know, is a completely new feature for the continuous model that has not been reported in the above works. As L = 0, P− collapses with the origin and a weak Allee effect is obtained if P+ remains an isolated singularity. In this case, the model predicts two phenomena depending on parameter values. In the first one, the origin has a separatrix dividing a repelling parabolic sector and a hyperbolic sector; see Lemma 1(iii). This means that every solution with initial conditions (x0 , y0 ) in the interior of the first quadrant near the origin ultimately tends to go away from the origin allowing both populations to not extinct. The second case is more usual, see Lemma 1(iv), and has been observed in previous works [19, 32]; that is, the origin has a parabolic attracting sector and a hyperbolic sector divided by an attracting separatrix, representing a critical line which has to be trespassed for populations in order to survive in time and not extinct. In this work, we also show in Theorem 8 that, for an open subset of the parameter values, three limit cycles can coexist in the open quadrant Ω. This result has been recently observed [26] but just in general cubic polynomial systems. In our case, only two of these limit cycles are concentric though, since infinitesimal ones are generated by Hopf bifurcation, and the third one is generated by a homoclinic bifurcation [23, 29] and surrounds another focus as stated in Lemma 6 and Theorem 8. Furthermore, the model can achieve the phenomenon of multistability by the existence of four ω-limit sets in the first quadrant as L < 0 and a strong Allee effect is present. A locally stable cycle surrounds a locally stable equilibrium point in the first quadrant as well as an unstable limit cycle, which serves as their separatrix in the phase plane; i.e., there is a range for population sizes for which there exists both autoregulation for the predator-prey system and prey and predator populations approach equilibrium, depending upon the population size. The third ω-limit set is the origin, stating that extinction is always possible in the presence of a strong Allee effect in our model. The fourth ω-limit set is a stable focus surrounded by the third limit cycle mentioned above. In order to illustrate the result stated in Theorem 8 and the multistability, Figure 10 shows a numerical simulation for system (5), which was affected with the software MATLAB [30]. Here a = 3, b = 1, k = 6.79211, M = 4.07697, c = 4.05116, and L = L+ + 0.000001 ≈ −0.498898 (see (12)). Initial conditions at (x(0), y(0)) are (0.8, 0.8), (1.01, 1.01), (2, 2), (2.2, 2.2), (1, 0.95), and (0.4, 0.1). In Figure 10 it is clear that a limit cycle lies between the orbits by the points (2, 2) and (2.2, 2.2), respectively, since the ω-limit of (2, 2) is the origin, but the ω-limit of (2.2, 2.2) is the stable focus located in the diagonal near x = 2.5. On the other hand, the two infinitesimal limit cycles are surrounding the singular point (1, 1), but they are too small for both the scale of this figure and these values of parameters. Hence, we present another numerical simulation in Figure 11 where both infinitesimal limit cycles are more visible. Here a = 3, b = 1, k = 6.56102, M = 4.16516, c = 4.16276, and L = L+ + 0.01 ≈ −0.480256 (see (12)). Initial conditions at (x(0), y(0)) are (0.6, 0.6), (0.8, 0.8), and (0.95, 0.95). Behaviors of the three positive oriented orbits in Figure 11 imply that, between two of them, attracting and repelling limit cycles exist. Furthermore, it can be easily verified that, for these values of parameters, the divergence of vector field (5) in singu-
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
1260
4 3.5 3
y
2.5 2 1.5 1 0.5 0 0
1
2
3
4
5
x Fig. 10. Numerical simulation of the results of Theorem 8.
1.5 1.4 1.3
y
1.2 1.1 1 0.9 0.8 0.7 0.6 0.4
0.6
0.8
1
1.2
1.4
1.6
x
Fig. 11. Numerical simulation of the infinitesimal limit cycles.
larity (1, 1) is negative, and, consequently, this equilibrium point is a local attractor and the innermost limit cycle is the unstable one. Several natural predator-prey communities have been studied, each one of them featuring an ecologically stable cycle, that is, a periodic orbit that must be somewhat insensitive to the perturbations outer to the interaction. In this work, we have answered partially one of the almost impossible projects proposed in [11] which is the following: Find a predator-prey or other interacting system in nature, or construct one in the laboratory, with at least two ecologically stable cycles. This shows the ecological relevance of the existence of multiple limit cycles in predation models and the importance of our result, which should serve for the outlined problem to be actually feasible in a biological lab with appropriate little creatures [11]. REFERENCES ´lez-Olivares, and E. Sa ´ez, Two limit cycles in a Leslie–Gower [1] P. Aguirre, E. Gonza predator-prey model with additive Allee effect, Nonlinear Anal. Real World Appl., to appear. [2] A. Andronov, E. Leontovich, I. Gordon, and A. Maier, Theory of Bifurcations of Dynamical Systems on a Plane, Israel Program for Scientific Translations, Jerusalem, 1971.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
THREE LIMIT CYCLES IN A LESLIE–GOWER
1261
[3] D. Arrowsmith and C. Place, An Introduction to Dynamical Systems, Cambridge University Press, New York, 1990. [4] D. K. Arrowsmith and C. M. Place, Dynamical Systems. Differential Equations, Maps and Chaotic Behaviour, Chapman and Hall, London, 1992. [5] J. Bascompte, Extinction thresholds: Insights from simple models, Ann. Zool. Fennici, 40 (2003), pp. 99–114. [6] L. Berec, E. Angulo, and F. Courchamp, Multiple Allee effects and population management, Trends Ecology Evol., 22 (2007), pp. 185–191. [7] T. R. Blows and N. G. Lloyd, The number of limit cycles of certain polynomial differential equations, Proc. Roy. Soc. Edinburgh Sect. A, 98 (1984), pp. 215–239. [8] D. S. Boukal and L. Berec, Single-species models and the Allee effect: Extinction boundaries, sex ratios and mate encounters, J. Theoret. Biol., 218 (2002), pp. 375–394. [9] C. W. Clark, Mathematical Bioeconomic. The Optimal Management of Renewable Resources, 2nd ed., John Wiley and Sons, New York, 1990. [10] C. W. Clark, The Worldwide Crisis in Fisheries: Economic Models and Human Behavior, Cambridge University Press, New York, 2006. [11] C. S. Coleman, Hilbert’s 16th problem: How many cycles? in Differential Equation Models, Modules Appl. Math. 1, M. Braun, C. S. Coleman, D. A. Drew, and W. F. Lucas, eds., Springe, New York, 1983. [12] J. B. Collings, Bifurcation and stability analysis of a temperature-dependent mite predatorprey interaction model incorporating a prey refuge, Bull. Math. Biol., 57 (1995), pp. 63–76. [13] F. Courchamp, T. Clutton-Brock, and B. Grenfell, Inverse dependence and the Allee effect, Trends Ecology Evol., 14 (1999), pp. 405–410. [14] F. Courchamp, B. T. Grenfell, and T. H. Clutton-Brock, Impact of natural enemies on obligately cooperative breeders, Oikos, 91 (2000), pp. 311–322. [15] A. M. De Roos and L. Persson, Size-dependent life-history traits promote catastrophic collapses of top predators, Proc. Natl. Acad. Sci. USA, 99 (2002), pp. 12907–12912. [16] B. Dennis, Allee effects: Population growth, critical density and the chance of extinction, Natur. Resource Modeling, 3 (1989), pp. 481–538. [17] F. Dumortier, Singularities of Vector Fields, Monogr. Mat. 32, IMPA, Rio de Janeiro, 1978. [18] R. Etienne, B. Wertheim, L. Hemerik, P. Schneider, and J. Powell, The interaction between dispersal, the Allee effect and scramble competition affects population dynamics, Ecol. Model., 148 (2002), pp. 153–168. ´lez-Yan ˜ez, and E. Gonza ´lez-Olivares, Conse[19] J. D. Flores, J. Mena-Lorca, B. Gonza quences of depensation in a Smith’s bioeconomic model for open-access fishery, in Proceedings of the 2006 International Symposium on Mathematical and Computational Biology BIOMAT 2006, Manaus, Brazil, 2007, pp. 219–232. [20] H. I. Freedman and G. S. K. Wolkowicz, Predator-prey systems with group defence: The paradox of enrichment revisted, Bull. Math. Biol., 8 (1986), pp. 493–508. [21] J. Gascoigne and R. N. Lipcius, Allee effects in marine systems, Marine Ecology Prog. Ser., 269 (2004), pp. 49–59. [22] J. Gascoigne and R. N. Lipcius, Allee effects driven by predation, J. Appl. Ecology, 41 (2004), pp. 801–810. [23] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, New York, 1983. [24] I. Hanski, L. Hansson, and H. Henttonen, Specialist predators, generalist predators and the microtine rodent cycle, J. Animal Ecology, 60 (1991), pp. 353–367. [25] I. Hanski, H. Hentonnen, E. Korpimaki, L. Oksanen, and P. Turchin, Small-rodent dynamics and predation, Ecology, 82 (2001), pp. 1505–1520. [26] X. Huang, Y. Wang, and L. Zhu, One and three limit cycles in a cubic predator-prey system, Math. Methods Appl. Sci., 30 (2007), pp. 501–511. [27] A. Korobeinikov, A Lyapunov function for Leslie-Gower predator-prey models, Appl. Math. Lett., 14 (2001), pp. 697–699. [28] Y. Li and D. Xiao, Bifurcations of a predator-prey system of Holling and Leslie types, Chaos Solitons Fractals, 34 (2007), pp. 606–620. [29] D. Luo, X. Wang, D. Zhu, and M. Han, Bifurcation Theory and Methods of Dynamical Systems, Adv. Ser. Dyn. Syst. 15, World Scientific, Singapore, 1997. [30] MATLAB: The Language of Technical Computing, Using MATLAB (Version 7.0), MathWorks, Natwick, MA, 2004. [31] R. M. May, Stability and Complexity in Model Ecosystems, 2nd ed., Princeton University Press, Princeton, NJ, 2001.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1262
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
´lez-Olivares, and B. Gonza ´lez-Yan ˜ez, The Leslie-Gower [32] J. Mena-Lorca, E. Gonza predator-prey model with Allee effect on prey: A simple model with a rich and interesting dynamics, in Proceedings of the 2006 International Symposium on Mathematical and Computational Biology BIOMAT 2006, Manaus, Brazil, 2007, pp. 105–132. ´ez and E. Gonza ´lez-Olivares, Dynamics of a predator-prey model, SIAM J. Appl. [33] E. Sa Math., 59 (1999), pp. 1867–1878. [34] P. A. Stephens and W. J. Sutherland, Consequences of the Allee effect for behaviour, ecology and conservation, Trends Ecology Evol., 14 (1999), pp. 401-405. [35] F. Takens, Unfoldings of certain singularities of vector fields: Generalized Hopf bifurcations, J. Differential Equations, 14 (1973), pp. 476–493. [36] R. J. Taylor, Predation, Chapman and Hall, New York, 1984. [37] P. Turchin, Complex Population Dynamics. A Theoretical/Empirical Synthesis, Monogr. Population Biol., 35, Princeton University Press, Princeton, NJ, 2003. [38] X. Wang, G. Liang, and F.-Z. Wang, The competitive dynamics of populations subject to an Allee effect, Ecol. Model., 124 (1999), pp. 130–168. [39] G. S. W. Wolkowicz, Bifurcation analysis of a predator-prey system involving group defense, SIAM J. Appl. Math., 48 (1988), pp. 592–606. [40] G. S. K. Wolkowicz, H. Zhu, and S. A. Campbell, Bifurcation analysis of a predator-prey system with nonmonotonic functional response, SIAM J. Appl. Math., 63 (2003), pp. 636– 682. [41] S. Wolfram, Mathematica: A System for Doing Mathematics by Computer, Addison–Wesley, Longman, Boston, 1988 [42] D. Xiao and S. Ruan, Global analysis in a predator-prey system with nonmonotonic functional response, SIAM J. Appl. Math., 61 (2001), pp. 1445–1472. [43] D. Xiao and S. Ruan, Bifurcations in a predator-prey system with group defense, Int. J. Bifurcation and Chaos, 11 (2001), pp. 2123–2131. [44] D. Xiao and H. Zhu, Multiple focus and Hopf bifurcations in a predator-prey system with nonmonotonic functional response, SIAM J. Appl. Math., 66 (2006), pp. 802–819. [45] D. Xiao and K. F. Zhang, Multiple bifurcations of a predator-prey system, Discrete Contin. Dyn. Syst. Ser. B, 8 (2007), pp. 417–433.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
c 2009 Society for Industrial and Applied Mathematics
THREE LIMIT CYCLES IN A LESLIE–GOWER PREDATOR-PREY MODEL WITH ADDITIVE ALLEE EFFECT∗ ‡ , AND EDUARDO SAEZ † ´ ´ PABLO AGUIRRE† , EDUARDO GONZALEZ-OLIVARES
Abstract. In this work, a bidimensional continuous-time differential equations system is analyzed which is derived of Leslie-type predator-prey schemes by considering a nonmonotonic functional response and Allee effect on population prey. For the system obtained we describe the bifurcation diagram of limit cycles that appears in the first quadrant, the only quadrant of interest for the sake of realism. We show that, under certain conditions over the parameters, the system allows the existence of three limit cycles: The first two cycles are infinitesimal ones generated by Hopf bifurcation; the third one arises from a homoclinic bifurcation. Furthermore, we give conditions over the parameters such that the model allows long-term extinction or survival of both populations. In particular, the presence of a weak Allee effect does not imply extinction of populations necessarily for our model. Key words. stability, limit cycles, homoclinic orbits, bifurcations, predator-prey models, Allee effect AMS subject classifications. 92D25, 34C, 58F14, 58F21 DOI. 10.1137/070705210
1. Introduction. This work deals with a continuous predator-prey model considering the following: (i) the Allee effect [6, 13, 16, 34] affecting the prey population, (ii) the functional response of predators of nonmonotonic type, and (iii) a predator growth function of logistic type. Other inherent assumptions for the model are that population size varies only in time and it is uniformly distributed in space, there is no division of ages or sex, and it is not affected by abiotic factors. It is known that the Allee effect refers to a positive density dependence in prey population growth at low prey densities [34], and it occurs whenever fitness of an individual in a small or sparse population decreases as the population size or density also declines [6]. To be able to recognize the consequences of the Allee effect on reproduction, conservation and behavior of species has become an important aim over the last years. The analysis and understanding of this phenomenon can bring important benefits not only for ecology but also for various applied engineering disciplines such as agropecuary, fishing, and forestal industries. Different mechanisms generating Allee effects have been suggested (see Table 1 in [6]), being largely studied as singular entities and usually describing a situation in which the population growth rate decreases under some minimum critical density [5] or when a limited population growth capacity is observed [15]. In some cases this growth rate might be even negative, causing an extinction threshold [5]. In other words, the Allee effect may be understood as the cause of the increase in extinction risk at low densities [16], introducing in some cases a population threshold that has to be exceeded by population to be able to grow. This effect is also named in population dynamics as the negative competition effect [38]; in fisheries sciences, it is called depensation ∗ Received by the editors October 12, 2007; accepted for publication (in revised form) October 6, 2008; published electronically February 4, 2009. This work was partially supported by USM grant 12.08.04. http://www.siam.org/journals/siap/69-5/70521.html † Departamento de Matem´ atica, Universidad T´ecnica Federico Santa Mar´ıa, Casilla 110-V, Valpara´ıso, Chile ([email protected], [email protected]). ‡ Instituto de Matem´ aticas, Pontificia Universidad Cat´ olica de Valpara´ıso, Casilla 4059, Valpara´ıso, Chile ([email protected]).
1244
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
THREE LIMIT CYCLES IN A LESLIE–GOWER
1245
[9, 16], and in epidemiology, its analogous is the eradication threshold, the population level of susceptible individuals below which an infectious illness is eliminated from a population [5]. More precisely, the (component) Allee effect is any positive relationship between any measurable component of individual fitness and population size or density [6, 13, 34], and it might be the strong Allee effect [6] or critical depensation [9] that implies the existence of a threshold population level [5, 9], whereas the weak Allee effect [34] or noncritical depensation [9] does not have it. When the species is submitted to a strong Allee effect, it may have a bigger tendency to be less able to overcome these additional mortality causes, to have a slower recovery, and to be more prone to extinction than other species [34]. In most predation models it has been considered that the Allee effect has influence only on the prey population, and this effect is independent of the functional response or consumption rate that reflects the change on predation due to the prey’s population size. Quantitatively, it is assumed that the functional response has influence on the extension of the bistability region [15]. The Allee effect has been modeled in different ways using various mathematical tools and, in a first approach, like a deterministic phenomenon frequently associated to population’s stochastic fluctuations [5]. For instance, if x = x (t) indicates the population size, the most usual continuous growth equation to express the Allee effect is given by dx x =r 1− (x − m) x, dt K featuring the multiplicative Allee effect. Clearly, if m = 0, we have the weak Allee effect and if m > 0, it has the strong Allee effect. Other mathematical forms have been proposed to describe this phenomenon [8]. In this work, we consider the natural growth function deduced in [16] and [34] given by the equation dx x m (1) =r 1− − x dt K x+b that we call additive Allee effect. Recent ecological research suggests the possibility that two or more Allee effects generate mechanisms acting simultaneously on a single population (see Table 2 in [6]), and the combined influence of some of these phenomena has been named the multiple Allee effects. In the interacting populations, the predation can be largely reduced due to better ability of prey to avoid predation when their population size is large enough [42, 43]; but at low population densities, there could be a low effectiveness of antipredator vigilance, which reflects an Allee effect. For some marine species it has been shown that the per capita growth decreases as the size population is reduced below some critical level, and two proposed causes of depensation or Allee effect are the following: reduced breedings success at low densities of the population and increased relative predation on small populations [21]. This situation has happened in many real fisheries as a result of overfishing when man acts as predator [9]. Considering these aspects, (1) can be rewritten as dx rx x = 1− (x − xm )x, (2) dt x+b xK
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1246
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
where xm = 12 (K−b− 1r (r(K + b)2 − 4mK)) and xK = 12 (K−b+ 1r (r(K + b)2 − 4mK)) for r(K + b)2 − 4mK > 0. Hence, (2) represents two types of Allee effect affecting the same population since xm expresses the minimum of viable population and the r x factor r (x) = x+b indicates the impact of an Allee effect due to other causes affecting the intrinsic growth rate, for example, the predation reducing breeding success at low densities [9, 21]. On the other hand, the functional response of predators or consumption rate function refers to the change in the density of prey attached per unit time per predator as the prey density changes [20, 42]. In most predator-prey models considered in the ecological literature, the predator response to prey density is assumed to be monotonic increasing, the inherent assumption being that as there is more prey in the environment, it is better for the predator [20]. However, there is evidence that indicates that this need not always be the case, for instance, when a type of antipredator behavior (APB) exists. Group defense is one of these, and the term is used to describe the phenomenon whereby predators decrease, or are even prevented altogether, due to the increased ability of the prey to better defend or disguise themselves when their number is large enough [20, 42, 39, 43], and in this case a nonmonotonic functional response is better. For example, lone musk ox can be successfully attacked by wolves; however, large herds of them can be attacked but with rare success. Another manifestation of an APB in which a nonmonotonic functional response (or Holling-type IV or Monod Haldane) can be used is the phenomenon of aggregation, a social behavior of prey in which prey congregate on a fine scale relative to the predator so that the predator’s hunting is not spatially homogeneous [36], such as succeeds with mile-long schools of certain class of fishes. In this case, a primary advantage of schooling seems to be confusion of the predator when it attacks. The more important benefit of aggregation is an increasing in wariness. Moreover, aggregation can both decrease the vulnerability to be attacked and increase the time that group members can devote to activities other than surveillance [36]. Other related examples of nonmonotone consumption occur at the microbial level where evidence indicates that when faced with an overabundance of nutrient the effectiveness of the consumer can begin to decline. This is often seen when microorganisms are used for waste decomposition or for water purification, a phenomenon that is called inhibition [20, 42, 43]. In these cases, the functional response curves have an upper bound on the rate of predation per individual predator at some prey density, in contrast to the old Lotka– Volterra model which assumed a linear relationship between prey density and the rate of predation over the entire range of prey densities [36]. In this work, we use the function h(x) = x2qx +a , also employed in [28, 42, 39, 43] and corresponding to the qx Holling-type IV functional response [36], which is generalized as h(x) = x2 +bx+a in [44, 40] for a Gause model. This generalized expression is derived by Collings [12], who affirms that this type of functional response seems a reasonable possibility if it is assumed that prey and webbing densities are directly related. In [40], for the corresponding Gause model the existence of two limit cycles is proven. Moreover, this Gause model exhibits bifurcation of cusp type with codimension two [4, 23] or Bogdanov–Takens bifurcation [44, 40]. We note that the phenomena of the Allee effect and aggregation described by a nonmonotonic functional response are quite compatible and justify our assumptions in the model studied.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
THREE LIMIT CYCLES IN A LESLIE–GOWER
1247
The last aspect in our equations is a feature of Leslie-type predator-prey models [37] or the Leslie–Gower model [27], in which the conventional environmental carrying capacity Ky is proportional to prey abundance x [31], that is, Ky = nx, as in the May–Holling–Tanner model [4] and other models recently analyzed [28, 45]. In [37] it is affirmed that Leslie models can lead to anomalies in their predictions since they predict that even at very low prey density, when the consumption rate by an individual predator is essentially zero, predator population can nevertheless increase if predator population size is even smaller than prey population size. However, Leslie models are recently employed to model vole-weasel dynamics with a parameter time dependent [25], and the autonomous model proposed in [24] is analyzed in [45]. This scheme of modeling differs from a more common Gause-type model [40] in which the predator equation is based on the mass action principle, since the numerical response is dependent on functional response. Although it may seem that the two aspects considered in the model contradict each other since the prey population exhibits the Allee effect for low densities, while a nonmonotone functional response is suggested for the aggregation (group defense) when the prey population size is large, it is known that predation induces an Allee effect. Strikingly, a wide range of predator-driven Allee effects have been reported (see Table 2 in [22]); in particular, there is the case of the Atlantic cod (Gadus morhua) that forms schools during the day, since commercial fishing (man as predator) provokes stock collapse because a higher proportion of this aggregative population is caught per unit effort when population declines [10]. Also, for obligately cooperative breeders as the African wild dog (Lycaon pictus) and meerkat (Suricata suricatta), there is a similar situation, because juvenile survival is lower in small groups than large groups in areas with high predator densities but lower in large groups than small groups in areas with low predator densities [14, 22]. From a mathematical point of view, simple models for the Allee effect may reveal a lot about its dynamics, and, reciprocally, different nonequivalent dynamics for the same model will have different biological interpretations. For instance, in our model we obtain the existence of a subset of parameter values for which three limit cycles appear but surrounding different equilibrium points. These limit cycles are not only periodic solutions of the system but also along the attractive behavior of the origin, allowing the phenomenon of multistability of our predator-prey system, that is, the existence of four ω-limit sets in the first quadrant. This paper is structured as follows: The model and the main results are presented in section 2. In section 3 we give the proofs of the main statements. Finally, an interpretation of the results is given in section 4, complemented with some numerical simulations. 2. The model. Let us consider the bidimensional system of ordinary differential equations
m x − xqxy x˙ = r 1 − xk − x+b 2 +a , Xμ : (3)
y y˙ = sy 1 − nx , where (x, y) ∈ A = {(x, y)| x > 0, y ≥ 0} and μ = (r, a, b, k, m, n, q, s) ∈ R8+ . In system (3), x(t) and y(t) denote prey and predator densities, respectively, as functions of time t. Furthermore, the parameters have the following meanings: (a) r and s are the intrinsic growth rates or biotic potential of the prey and predators, respectively.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1248
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
(b) q is the maximal predator per capita consumption rate, i.e., the maximum number of prey that can be eaten by a predator in each time unit. (c) a is the number of prey necessary to achieve one-half of the maximum rate q. (d) n is a measure of the food quality that the prey provides for conversion into predator births. (e) k is the prey environment carrying capacity [12]. (f) m and b are constants that indicate the severity of the Allee effect that has been modeled. In order to describe the dynamics of Xμ we will consider a C ∞ -equivalent polynomial extension by the following change of coordinates: u(u + b)(u2 + a) τ + 2 = (x, y, t), ζ : R2 ×R+ −→ R ×R such that ζ(u, v, τ ) = u, nv, 0 0 s (4) 2 +a) where detDζ(u, v, τ ) = nu(u+b)(u > 0. s Also, let us consider new parameters given by ϕ1 : R8+ −→ R × R7+ , m nq br − m , a, b, k, , n, , s := (L, a, b, k, M, n, c, s), ϕ1 (r, a, b, k, m, n, q, s) = s s s 3
s with Jacobian detDϕ1 (r, a, b, k, m, n, q, s) = bn > 0, which says that ϕ1 is invertible. For simplicity, let us rename the new coordinates u → x, v → y. Then in the new −1 coordinates we have the vector field Yη = ζ∗ Xμ = (Dζ) ◦ Xμ ◦ ζ, where ⎧
⎨x˙ = x (a + x2 ) (b+x)(L+M ) 1 − x − M x − cxy(b + x) , b k Yη : (5) ⎩ 2 y˙ = (b + x)(a + x )(x − y)y,
¯ = {(x, y)|x, y ≥ 0} and the new vector of parameters with (x, y) ∈ Ω η = (L, a, b, k, M, c) ∈ D0 is given by the natural projection, where D0 = η ∈ R × R5+ | M + L > 0 . (6) Notice that vector field (5) is a polynomial extension of the original system (3) ¯ including axis x = 0. Moreover, Yη and Xμ are clearly to the whole first quadrant Ω ∞ C -equivalent in Ω. First, we study the local behavior of singularities in the coordinate axis, that is, in absence of prey and predator, respectively. In system (5), it is immediate that ∂ Yη (0, y) = −aby 2 ∂y . Then the x = 0 axis in Ω is an invariant manifold of vector field (5) and the origin is its only singularity, being an attracting one. Moreover, if we ∂ ∂ ∂ denote Yη (x, y) = P (x, y) ∂x + Q(x, y) ∂y , we have Yη (x, 0) = P (x, 0) ∂x . Hence, the y = 0 axis in Ω is also an invariant manifold of (5). For vector field Yη let us define Sing(E) and Sing(Ω) as the sets of singularities of (5) in the y = 0 axis in Ω and in the open set Ω, respectively. In order to describe the dynamics of singularities in Sing(E), we consider the following regions in parameter space D0 : Λ1 Λ2 Λ4 Θ+ Γ+
= Δ−1 (−∞, 0) ∪ Δ−1 (0, ∞) ∩ Θ−1 (0, ∞) ∩ Υ−1 (−∞, 0) ; −1 −1 −1 = Δ (0, ∞) ∩ Θ (∞, 0) ∩ Υ (−∞, 0); = {η ∈ D0 | 0 < b ≤ L}; = Θ−1 [0, ∞) ∩ Υ−1 (0); = Γ−1 (0, ∞) ∩ Υ−1 (0);
Λ3 D Θ− Γ−
= {η ∈ D0 | 0 < L < b}; = Δ−1 (0) ∩ Θ−1 (0, ∞); = Θ−1 (−∞, 0) ∩ Υ−1 (0); = Γ−1 (−∞, 0] ∩ Υ−1 (0);
(7)
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
THREE LIMIT CYCLES IN A LESLIE–GOWER
b
1249
y
Λ1 y
• P0/+
• P0
D
L = −M
Λ2
y
x
b = L b = k •
y
• • P0 P±
Λ3
x
x y
y
• • P0/−P+
• • • x P0P− P+ L = 0
• • P0 P+
x
Λ4
x y
• • P0 P+
x
L
Fig. 1. Bifurcation diagram in plane Lb = {η ∈ D0 | a, k, M, c are constant} in parameter space for singularities in the y = 0 axis.
where Υ(η) = sign(L), Θ(η) = b − k, and Δ(η) = (M + L)(k − b)2 + 4bkL, √ √ √ Γ(η) = −2kcb + a M −a M + aM + 4k 2 c . Moreover, in the xy plane let us consider the points P0 = (0, 0), P+ = (x+ , 0), and P− = (x− , 0), with (M + L)(k − b) ± (M + L)Δ(η) (8) . x± = 2(M + L) Lemma 1. For vector field Yη the origin P0 is a non-hyperbolic singularity and it satisfies the following: (i) It has two hyperbolic sectors divided by a repulsing separatrix if η ∈ Γ+ ∩ Θ+ . (ii) It has only a hyperbolic sector if η ∈ Λ4 . (iii) It has a hyperbolic sector and a repulsing parabolic sector if η ∈ Λ3 ∪(Θ− ∩ Γ+ ) . (iv) It has a hyperbolic sector and an attracting parabolic sector if η ∈ Θ− ∩ Γ− . (v) It is a local attractor if η ∈ D ∪ Λ1 ∪ Λ2 ∪ (Θ+ ∩ Λ− ) . A qualitative diagram of these results in plane Lb = {η ∈ D0 | a, k, M, c are constant} is shown in Figure 1. Lemma 2. For vector field Yη the following statements hold: (i) P+ ∈ Sing(E) if η ∈ D ∪ Θ− ∪ Λ3 ∪ Λ4 . Moreover, P+ is a saddle node if η ∈ D; otherwise, it is a hyperbolic saddle. (ii) P− ∈ Sing(E) if η ∈ Λ2 . Moreover, P− is a hyperbolic repulsing node. A qualitative diagram of these results in plane Lb = {η ∈ D0 | a, k, M, c are constant} is shown in Figure 1. About the boundness of the system, in [1] it is proven that no trajectory of vector field Yη has the infinity as ω-limit. Since the coordinate axes are invariant, then our system is bounded. We now state the same result but with a different approach. Theorem 3. Consider Aw = {(x, y) ∈ R2 : 0 ≤ x ≤ w, 0 ≤ y}. For every η ∈ D0 , there is a w∗ ≥ 0 such that if w > w∗ , then Aw is a trapping domain for system (5), meaning that it is invariant for positive time evolution and also captures all trajectories starting in Ω.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1250
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
We take care now of the dynamics at the interior of the first quadrant. It is clear that every singularity of system (5) in the open quadrant Ω must be over the diagonal y = x. Furthermore, over the isocline y = x the vector field Yη has the form ∂ Yη (x, x) = x2 p4 (x) ∂x , where p4 (x) is a polynomial of degree 4 in x. Therefore, the existence of equilibrium points in the open first quadrant is given by the existence of real and positive roots of the polynomial p4 (x). In order to describe the bifurcation diagram for an equilibrium point of Yη in Ω, we change parameters η = (L, a, b, k, M, c) to get a new vector ξ = (L, a, b, k, V, q) by means of the transformation ϕ2 : R6 −→ R6 , (L, a, b, k, V, q) → (L, a, b, k, M (ξ), c(ξ)),
(9) with c(ξ) = M (ξ) =
aLq 4 (b−k+2q)−q 6 (2b2 +kL−b(2k+L−4q)−2(k+L)q+2q 2 ) kq 6 −bq 4 (b+q)6 (a+q 2 )2 (a+3q 2 )(aL+(L−2b)q 2 )+V . bq 6 (b+q)4 (a+q 2 )3 (a+3q 2 )
−
(b−k+q)V kq 6 b2 (b+q)5 (a+q 2 )2 (a+3q 2 ) ,
Let D1 = ϕ−1 2 (D0 ) be the new admissible region in parameter space, and let us consider the projection π : R6 −→ R3 such that
π(L, a, b, k, q, V ) = (L, b, k).
If we name
D2 = π D1 ∩ {(L, a, b, k, q, V )| a = 3b2 , q = b, V > 0} ∩ E −1 (0, ∞) ⊂ R × R2+ , with E(L, b, k) = −18432b19 L + V, it is easy to check that D2 = F −1 (0, ∞) ∩ E −1 (0, ∞), with F (L, b, k) = 12288b19 (b − 2L) + V. Changing the time t → bkt, let us call Yξ to the new vector field qualitatively equivalent to (5). Now, it is straightforward to check that there is a singularity in the isocline y = x at the point with coordinates (b, b). In order to describe the bifurcation diagram of this singularity of vector field Yξ in parameter space (L, b, k), we consider the projection (see Figure 2) 2 −→ R2 such that πR : SR
πR (L, b, k) = (L, b),
2 = {(L, b, k)| L2 + b2 + k 2 = R2 , b, k, R > 0}. Since function F (L, b, k) does where SR not depend on k, it is clear that in D2 ⊂ R × R2+ the surface F −1 (0) is a straight 2 half-cylinder with axis parallel to the k axis, so πR F −1 (0) ∩ SR is exactly the curve
2 , F = {(L, b)| F (L, b, k) = 0, k > 0} ∩ πR D2 ∩ SR
which is located in the interior of the first quadrant of plane Lb because equaV tion F (L, b, k) = 0 defines implicitly L = 2b + 24576b 19 > 0. Furthermore, (0, 0) ∈ −1 F (0, ∞), so the origin is in the admissible region for parameters in the upper
half2 . With plane Lb. The same arguments follow to sketch the curve πR E −1 (0) ∩ SR 2 this, the qualitative shape of domain πR (D2 ∩ SR ) is shown in Figure 3(a).
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1251
THREE LIMIT CYCLES IN A LESLIE–GOWER k 2 SR
L
b πR b = R k = 0
L = −R
L = R
2 to the plane Lb. Fig. 2. Projection πR from the sphere SR
2 πR D2 ∩ SR
b
Σ0 2
2 πR E −1 (0) ∩ SR
b p2 •
2 πR F −1 (0) ∩ SR
• p3 •
p1 •
p4
L
(a)
H
L
(b)
2 . (b) Qualitative shape of Fig. 3. (a) Qualitative shape of parameters domain πR D2 ∩ SR 2 curves H and Σ0 (R) in parameters domain.
2 in plane Lb, we define the sets Finally, in πR D2 ∩ SR
sg() 2 Σ1 (R) = πR {(L, b, k) ∈ D2 | sg(T (L, b, K)) = sg( )} ∩ SR ,
sg() 2 Σ2 (R) = πR {(L, b, k) ∈ D2 | sg(G(L, b, K)) = sg( )} ∩ SR , H = Σ01 (R),
(10)
with T (L, b, k) = −36864b21 + 55296b20 L − 18432b19 kL − 3bV + kV, G(L, b, k) = 1811939328b39 (b + L)(b + 2L) − 8bV 2 , and the qualitative shape of curves H and Σ02 (R), in terms of R, is as in Figure 3(b). Moreover, H ∩ Σ02 (R) consists of four points as in Figure 3(b), that is, (11)
H ∩ Σ02 (R) = {pi = (Li , bi ); i = 1, . . . , 4}.
Lemma 4. The singularity (b, b) of vector field Yξ is an order two weak focus if (L, b) ∈ H ∩ Σ02 (R) = {pi = (Li , bi ); i = 1, . . . , 4}. Moreover, the focus is repellor at p1 and p2 and attractor at p3 and p4 . Theorem 5. (i) The system Yξ has a bifurcation diagram as indicated in Figures 4(a)–(d) in a neighborhood of H ∩ Σ02 (R). (ii) The bifurcation curve S34 of semistable limit cycle of the system Yξ in the parameter space Σ− 1 (R) is as indicated in Figure 4(d). (iii) The bifurcation curve HC of a homoclinic loop surrounding the focus (b, b) of system Yξ is as indicated in Figure 4(a), in a tubular neighborhood of curve H in sector Σ− 2 (R) in parameter space. Lemma 6. In parameter space ξ = (L, a, b, k, q, V ) ∈ D1 of system Yξ , there exists −1 an open subset Υ ⊂ π −1 ◦ πR (p3 ) ∩{ξ | L < 0}, where the following statements hold:
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
1252
(a) S2 b
• p2 HC
• p3
H S1
S34 • p1
• p4 L
(b)
(c) H
•u0 •u1 S1
•u0
•u1
H
S2 • p2
• p1
•u0
•u0
•u2
•s0 •u2 •s1
•s1
•s0
•s1
(d) •s2
• p3
•u0 •u1 •s2 •s0
•s0
•s1
H S34
• p4
Fig. 4. (a) Bifurcation curves for focus (b, b) of vector field Yξ , where H is a Hopf bifurcation curve, S1 , S2 , and S34 are semistable limit cycle bifurcation curves, and HC is a Homoclinic bifurcation curve. (b) Bifurcation diagram for focus (b, b) in a neighborhood of p1 ∈ H. (c) Bifurcation diagram for focus (b, b) in a neighborhood of p2 ∈ H. (d) Bifurcation diagram for focus (b, b) in a tubular neighborhood of curve H between p3 and p4 . Notation •ij means a focus with local stability i ∈ {s, u} and weakness j ∈ {0, 1, 2}.
(i) Yξ has two foci (b, b), (xf , xf ) and two saddle points (xs , xs ), (xN , xN ) in the interior of the first quadrant Ω, where 0 < xs < b < xN < xf . (ii) There exists an open subset O ∈ R+ and a continuous function H : O −→ R+ such that if V = H(b), then there is a homoclinic loop L0 for vector field Yξ passing through saddle (xN , xN ) and surrounding focus (xf , xf ), with (xN , xN ) and (xf , xf ) as in part (i). Moreover, this homoclinic orbit is inner unstable. Lemma 7. Let Υ, H, and L0 be as in Lemma 6. There exist δ > 0 and ε > 0 such that, for every b, V with |H(b) − V | < ε, the vector field Yξ has, at most, one limit cycle L in the inner δ-neighborhood of L0 , and L is unstable. Theorem 8. There exists an open subset in parameter space such that vector field Yη in (5) has three hyperbolic limit cycles in the interior of the first quadrant Ω. Two of them are infinitesimal cycles surrounding a hyperbolic attracting focus (b, b); moreover, the outermost limit cycle is stable and the innermost unstable. The third limit cycle is unstable and surrounds a hyperbolic attracting focus (xf , xf ), with 0 < b < xf . 3. Proofs of the main results. The proofs of Lemmas 1 and 2 are straightforward and follow from the blowing-up method [17] for (0, 0) and the central manifold and Hartman’s theorems for P+ and P− ; for details, see [1].
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
THREE LIMIT CYCLES IN A LESLIE–GOWER •
1253
u
v
fw (x, y) = x − w = 0 Aw y x
• x = w
−1 (−∞, 0) ∩ Ω is an invariant set under Φ Fig. 5. The set Aw = fw Yη (t; · ), as stated in the proof of Theorem 3.
Proof of Theorem 3. We will prove that for every η ∈ D0 (see (6) for definition), there exists some w = w(η) > 0 such that the integral curves ΦYη (t; · ) of vector field Yη , for t > 0, enter an invariant subset Aw ⊂ Ω and do not leave again. Let us define the function fw : R2 −→ R given by fw (x, y) = x − w. For each w > 0 consider the line Lw = fw−1 (0) ∩ Ω. If we consider the Poincar´e compactification of system (5) given by + 2 Ψy : R2 × R+ 0 −→ R × R0 , Ψy (u, v, τ ) = (u/v, 1/v, τ ) = (x, y, t),
it is straightforward to check that the line Ψ−1 y (Lw ) contains the point (u, v) = (0, 0). Hence, we will show that for w sufficiently large, in the Poincar´e compactification of system (5), the compact set Aw = fw−1 (−∞, 0) ∩ Ω is an invariant set under ΦYη (t; · ) (see Figure 5). Recall that the axes x = 0 and y = 0 are invariant under ΦYη (t; · ), so it is sufficient to prove that the scalar product Yη · ∇fw < 0 in Lw for appropriate w > 0 as in Figure 5. We have that for (x, y) ∈ Lw
w2 a + w2 2 Yη · ∇fw |x=w = −cw (b + w)y + φ(w), bk where φ(w) = −(M + L)w2 + (k − b)(M + L)w + bkL. The discriminant of φ as a quadratic polynomial in w is (M + L)Δ(η) (see (7) for definition). Therefore, since M + L > 0 in D0 , if Δ < 0, then Yη · ∇fu |x=w < 0 for every w > 0. On the other side, if Δ ≥ 0, the roots of φ(w) are x− and x+ (see (8) for definition); hence, the statement holds for every w > x+ if x+ > 0 and for every w > 0 if x+ ≤ 0. This concludes the proof of Theorem 3. 2 Proof of Lemma 4. Let (L, b) ∈ πR (D2 ∩ SR ). In order to study the singularity (b, b) in vector field Yξ , we translate (b, b) to the origin and change the time by means of + 2 ζ : R2 × R + 0 −→ R × R0 , 8b3 E(L, b, k) t , (x, y, t) → x + b, y + b, 768b15
where E(L, b, k) = −18432b19 L + V.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1254
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
So if (L, b, k) ∈ E −1 (0, ∞), ζ is a C ∞ -equivalence. Let Zξ = ζ∗ Yξ ; then we have trDZξ (0, 0) =
8E(L, b, k)T (L, b, k) , 589824b25
with
T (L, b, k) = − 36864b21 − 55296b20 L + 18432b19 kL + 3bV − kV . 2 Moreover, if (L, b) ∈ H = πR (T −1 (0) ∩ SR ), we have detDZξ (0, 0) = V > 0, so the origin is a weak focus of order, at least, one for Zξ . In order to recognize the topological type of the origin as (L, b) ∈ H, we bring Zξ onto its Jordan canonical + 2 form at (0, 0). For that, consider the C ∞ -equivalence ϕ : R2 × R+ 0 −→ R × R0 , √ √ ϕ(u, v, τ ) = 192b10 u − V v, 192b10 u, τ / V = (x, y, t).
The qualitatively equivalent vector field in the new coordinates is ZξJ = ϕ∗ Zξ , and we have in a neighborhood of the origin ⎛ ⎛ ⎞ ⎞ 5 5 ∂ ∂ + ⎝u + , ZξJ (u, v) = ⎝−v + Ai,j ui v j + H.O.T.⎠ Bi,j ui v j + H.O.T.⎠ ∂u ∂v i,j=2 i,j=2 where H.O.T. denotes the higher order term and Ai,j = Ai,j (b, V, L), Bi,j = Bi,j (b, V, L). For j = 0, 1, 2, let lj be the first three Lyapunov quantities at the origin [2, 7, 33] of the vector field ZξJ . Since the trace of the linear part of the vector field at the origin is zero, we have that l0 = 0. Now, as l1 depends on the 3-jet of ZξJ , see [17], then l1 = (A02 A11 + A12 + A11 A20 + 3A30 + 2A02 B02 + 3B03 − B02 B11 − 2A20 B20 − B11 B20 + B21 )/8.
Using the Mathematica software [41] we have l1 =
36864b20 + V G(L, b, k), 8192b13 V 3/2
with G(L, b, k) = 1811939328b39 (b + L)(b + 2L) − 8bV 2 . Thereby, if (L, b) ∈ H ∩ Σ02 (R) (see (10) and (11)), the weakness of the origin depends only on l2 . On the other hand, it is known that l2 depends on the 5-jet of ZξJ . Then, again using the Mathematica software, as l0 = l1 = 0, we obtain l2± =
36864b20 + V N± , 98304b15 V 3/2
where
61 N± = − √−221259535220736b − 13589544960b41 V −√1963008b21 V 2 − 57bV 3 40 ± 3 30349983744b + 2211840b20 V + 43V 2 28311552b42 + b2 V 2
and the sign ± denotes the branch (12)
√ −27648b39 ± 3 g(b) L± = 36864b38
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1255
THREE LIMIT CYCLES IN A LESLIE–GOWER
0 -0.2 2 4 N+ -0.4 -0.6 6 0
0 0.2 b
0.2 V 0.4
0.4
0 -20 0 0 0.01∗N+ -40 -60 0 0.6 0.
20 40 V 60
0.8b 80 1001
(a)
(b)
Fig. 6. Graphics of the function N+ = N+ (b, V ) where we can see that N+ < 0: (a) graphic of N+ for (b, V ) ∈ [0, 0.5] × [0, 0.5]; (b) graphic of 10−2 × N+ for (b, V ) ∈ [0.6, 1] × [10, 100].
of curve Σ02 (R) given implicitly by equation G(L, b, k) = 0, where g(b) = 28311552b78 + b38 V 2 > 0. Since the branch L− intersects curve H at points p1 and p2 (see Figure 3(b)), then if (L, b) ∈ {p1 , p2 } ⊂ H ∩ Σ02 (R), l2 > 0 and the origin is a repellor weak focus of order two for ZξJ . On the other hand, the branch L+ intersects curve H at points p3 and p4 (see Figure 3(b)), and by looking at the graphic of N+ as a function of b and V it can be verified that N+ < 0. Figure 6(a) shows the graphic of N+ for the range (b, V ) ∈ [0, 0.5] × [0, 0.5], meanwhile in Figure 6(b) there is the graphic of 10−2 × N+ for (b, V ) ∈ [0.6, 1] × [10, 100], made with the software Mathematica [41]; here it can be seen that if b > 0, V > 0, as in our case, then by continuity N+ < 0. Therefore, if (L, b) ∈ {p3 , p4 } ⊂ H ∩ Σ02 (R), l2 < 0 and the origin is an attracting weak focus of order two for ZξJ . 2 Proof of Theorem 5(i). It is clear that T −1 (0) and SR are transversal manifolds in the parameter space and, by Lemma 4, there are four points
2 p1 , p2 , p3 , p4 ∈ πR T −1 (0) ∩ SR = H, where H is a Hopf bifurcation curve and where the singularity (b, b) is a repellor weak focus of order 2 of system Yξ at p1 and p2 and an attractor weak focus of order two at p3 and p4 ; i.e., the eigenvalues of DYξ (b, b) are on the imaginary axes and nonzero. Moreover, p1 , p2 , p3 , and p4 have codimension two, since the Lyapunov quantities l0 = l1 = 0 and l2 = 0. Thereby, let (L, b) ∈ H − {p1 , p2 , p3 , p4 }. Then we have the following (see Figure 7): (a) In a neighborhood of p4 , when (L, b) is below l1−1 (0), the singularity (b, b) of system Yξ is an order one attractor weak focus, because l1 < 0. Therefore, if > 0 is sufficiently small, the point (L + , b + ) is over the curve H; then the stability of system Yξ at the singularity (b, b) is reversed; hence, a unique hyperbolic stable limit cycle bifurcates (Hopf bifurcation). (b) When (L, b) is above l1−1 (0) in a neighborhood of p4 , the singularity (b, b) is an order one repellor weak focus surrounded by a hyperbolic stable limit cycle, because l1 > 0; then Hopf bifurcation from p4 occurs. Therefore, if > 0 is sufficiently small, the point (L + , b) is over the curve H; hence, the stability of system Yξ at the singularity (b, b) does not change and the limit cycle persists because it is hyperbolic. (c) When the point (L, b) is as in case (b), for > 0 sufficiently small, the point (L − , b) is under the curve H; then the stability of system Yξ at the singularity (b, b) is reversed; hence, a new Hopf bifurcation occurs and a
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1256
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ b
−1 Σ0 2 (R); l1 (0) p4 • −1 H; l0 (0) L
Fig. 7. Hopf bifurcation curve H near point p4 in parameter space, where the singularity (b, b) of vector field Yξ is an order two attracting weak focus.
second hyperbolic limit cycle bifurcates, which is surrounded by a hyperbolic stable limit cycle. By normal forms theory [17, 23], p4 has a 2-parameter versal unfolding. Under the same hypothesis that we have here, Takens [35] and Arrowsmith and Place [3] describe in detail the bifurcation diagram for codimension two singularity type. By using their results, we have that there exists a neighborhood Vδ (p4 ), with δ > 0, such that it has a diagram as in Figure 4, where S34 is a curve in which the unstable and stable limit cycles collapse (semistable limit cycle). The analysis near points p1 , p2 , and p3 follows in the same way. Proof of Theorem 5(ii). From part (i) of this Theorem, there is an open set C(R) ⊂ Σ− 1 (R) where system Yξ has two hyperbolic limit cycles surrounding the singularity (b, b). Moreover, C(R) is upperly bounded by the segment of curve H between p3 and p4 , because at any point in that segment focus (b, b) is surrounded by a unique infinitesimal limit cycle. On the other hand, at Σ− 1 (R)\ C(R), the singularity (b, b) has no infinitesimal limit cycle in its vicinity, so C(R) is bounded and the lower bound must correspond to a bifurcation curve where both limit cycles collapse into a semistable limit cycle; hence, by versality of the unfolding of the singularity (b, b) at p3 and p4 described in part (i) (see [35, 3]), it is clear that the semistable limit cycle bifurcation curve S34 must be as in Figure 4(d). Proof of Theorem 5(iii). By Hartman–Grobman’s theorem it is straightforward to see that vector field Yξ has a hyperbolic saddle singularity Ps = (s, s) with 0 < s < b if (L, b) ∈ H ∩ Σ− 2 (R). Let (L, b) ∈ H ∩Σ− 2 (R) in a neighborhood of p3 . By part (i), (b, b) is an attracting weak focus of order one with no limit cycle in its vicinity, and it can be checked that Ps ∈ α−lim(b, b). As we move along segment H ∩ Σ− 2 (R) towards p2 , the topological type of focus (b, b) remains the same according part (i); however, if (L, b) ∈ H ∩Σ− 2 (R) c in a neighborhood of p2 , Ps ∈ (α−lim(b, b)) and an infinitesimal limit cycle exists surrounding (b, b). Therefore, there exist δ > 0 sufficiently small and a differentiable function Ψ : Vδ (H) −→ R
such that (L, b) → Ψ(L, b),
defined in a tubular neighborhood Vδ (H) of radius δ of segment H ∩ Σ− 2 (R) such that vector field Yξ has a homoclinic loop passing through the saddle Ps and surrounding focus (b, b), if (L, b) ∈ HC = Ψ−1 (0). Remark. An extended and deeper proof of the existence of a homoclinic loop surrounding another singularity is given in Lemma 6. Nevertheless, the same arguments might have been used in this case as well.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1257
THREE LIMIT CYCLES IN A LESLIE–GOWER
Σ
Σ Φτ (ˆ r) p ˆ •
W s (N )
p ˆ •
ˆ • r • W s (N )
Φτ (ˆ r)
• p∗
•
• r ˆ •
W u (N )
• p∗
•
N
W u (N )
N
(a)
(b)
Fig. 8. Poincar´ e map for orbits near W u (N ) ∩ W s (xf ) = ∅.
W s (N )
in two cases: (a) W u (N ) ∩ W s (xf ) = ∅; (b)
Proof of Lemma 6(i). Let (L, b) = p3 ∈ H ∩ Σ02 . It is easy to check that Yξ has four singularities in isocline y = x at Ω, if L < 0. According to Lemma 4, one of these equilibria is an attracting weak focus of order two in (b, b). Due to Hartman–Grobman’s theorem, it is straightforward to see that the other singularities are a hyperbolic attracting focus in (xf , xf ) and two saddles (xs , xs ), (xN , xN ) with 0 < xs < b < xN < xf . The statement follows directly. Proof of Lemma 6(ii). Let (xN , xN ) and (xf , xf ) be as in the proof of part (i). If we name N = (xN , xN ), let us consider the function D(N ) = ∇ · Yξ (N ), and let Φt = ΦYξ (t; · ) be the flow of Yξ and p∗ = (x∗ , y ∗ ) be a point in W u (N ) with xN < x∗ . Let Σ be a one-dimensional local cross section to W s (N ) ∪ W u (N ) as in Figure 8. Let pˆ = (ˆ p1 , pˆ2 ) the unique point of Σ ∩ W s (N ) which never returns to Σ by the flow Φt with 0 < t < ∞, and let U ⊂ Σ be some neighborhood of pˆ. In Σ let us choose local coordinates given by the chart h : Σ −→ R, h(x, y) = y − pˆ2 . Let us consider the Poincar´e map P : U −→ Σ defined for a point r ∈ U by P(r) = Φτ (r), where τ = τ (r) is the time taken for the orbit Φτ (r) based at r to first return to Σ. Now consider the displacement function d : U −→ R given by r → d(r) = P(r) − r. Finally, let rˆ be a point in U located below pˆ as in Figure 8. Let b = 1 and V = 300. Then it can be seen by means of numerical simulations on Matlab [30] that D(N ) > 0, W u (N ) ∩ W s (xf ) = ∅, W u (N ) ∩ W s (N ) = ∅, and ω − lim Φt (p∗ ) = {(xf , xf )}. In the local coordinates given by h this implies that (13)
d(ˆ r) < 0.
A qualitatively equivalent situation is shown in Figure 8(a). Instead, if b = 1 and V = 400, it can be checked (numerically as well) that D(N ) > 0, W u (N ) ∩ W s (xf ) = ∅, and W u (N ) ∩ W s (N ) = ∅; there exists y † > xN such that the point (xN , y † ) ∈ Φt (p∗ ) ⊂ W u (N ). Therefore, we have (see Figure 8(b)) (14)
d(ˆ r) > 0.
Hence, by the continuous dependence of the orbits of the vector field on parameters, there exists V ∗ , 300 < V ∗ < 400, such that for b = 1 and V = V ∗ we have W u (N ) ∩ W s (N ) = ∅; that is, there is a homoclinic loop L0 ⊂ W u (N ) ∩ W s (N ).
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1258
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
L0 • • V
V ∗•
(1, V ∗ ) •
• • • 1
• F−1 (0) •
Bδ
• b
Fig. 9. Existence of a homoclinic orbit along the curve F−1 (0) in plane bV .
Therefore, there exists δ > 0 such that in a δ-ball Bδ (1, V ∗ ) ⊂ R2+ of the point (1, V ∗ ) in plane bV in parameter space, there is a locally differentiable function F : Bδ (1, V ∗ ) −→ R, (b, V ) → F(b, V ), such that if (b, V ) ∈ F−1 (0), there exists an inner unstable homoclinic loop L0 for vector field Yξ (see Figure 9). The existence of the function H : O −→ R+ in the statement follows from the implicit function theorem, due to the transversality shown in passing from (13) to (14). Proof of Lemma 7. Let the function d : U −→ R, the saddle N = (xN , xN ), and the point rˆ be as in the proof of Lemma 6 (see Figure 8). From the previous lemma, a homoclinic bifurcation (see [29]) occurs as V = H(b) since, by continuity, D(N ) = ∇ · Yξ (N ) > 0, N is a strong saddle [29], and L0 is inner unstable. Then, for U sufficiently small, by carrying out a perturbation on parameters, such that d(ˆ r) > 0, the homoclinic cycle breaks out and an unstable limit cycle appears surrounding the attracting focus (xf , xf ). Proof of Theorem 8. From Lemmas 4 and 7, in parameter space ξ = (L, a, b, k, q, V ) −1 ∈ D1 of system Yξ , there exists an open subset ΥL ⊂ π −1 ◦ πR (p3 ) ∩ {ξ | L < 0}, where system Yξ has four singularities in the interior of the first quadrant Ω: two saddles (xs , xs ), (xN , xN ), a two-order attracting weak focus (b, b), and a hyperbolic attracting focus (xf , xf ) surrounded by an unstable limit cycle L, with 0 < xs < b < xN < xf . Moreover, this limit cycle is hyperbolic, so it persists in an open subset of −1 π −1 (πR (C(R))) in parameter space, where C(R) is the bounded sector in plane Lb (see Figure 4(d)) where system Yξ also has two infinitesimal limit cycles surrounding focus (b, b), due to the proof of Theorem 5(ii). The statement of the theorem is immediate after recalling that vector field Yξ is induced by a vector field C ∞ -equivalent to Yη by means of the local difeomorphism ϕ2 in (9), since detDϕ2 (ξ) =
b3 q 12 (b
+
V > 0, + q 2 )4 (a + 3q 2 )3
q)10 (a
where ξ ∈ D1 ∩ π −1 (T −1 (0)). 4. Discussion. In this paper, a predator-prey model with the Allee effect [5, 16, 34, 38] and a Holling-type IV [36] functional response was considered, making a qualitative analysis of a bidimensional system of ordinary differential equations of polynomial type, which is topologically equivalent to the original one, with only six parameters. Employing a reparameterization and a time change, we showed that b, k, and L = br−m are the most significative parameters of this predator-prey model. s We note that this model predicts that, when L < 0, the system has three equilibria
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
THREE LIMIT CYCLES IN A LESLIE–GOWER
1259
P0 , P− , and P+ for which the predator density is null; moreover, at low population densities, both predator and prey populations may disappear, since the equilibrium P0 = (0, 0) is an attractor; furthermore, P− is a repelling node, and P+ is a saddle. Unlike other models, in this work we show that the weak Allee effect does not imply extinction of species necessarily for certain parameter values, which, as far as we know, is a completely new feature for the continuous model that has not been reported in the above works. As L = 0, P− collapses with the origin and a weak Allee effect is obtained if P+ remains an isolated singularity. In this case, the model predicts two phenomena depending on parameter values. In the first one, the origin has a separatrix dividing a repelling parabolic sector and a hyperbolic sector; see Lemma 1(iii). This means that every solution with initial conditions (x0 , y0 ) in the interior of the first quadrant near the origin ultimately tends to go away from the origin allowing both populations to not extinct. The second case is more usual, see Lemma 1(iv), and has been observed in previous works [19, 32]; that is, the origin has a parabolic attracting sector and a hyperbolic sector divided by an attracting separatrix, representing a critical line which has to be trespassed for populations in order to survive in time and not extinct. In this work, we also show in Theorem 8 that, for an open subset of the parameter values, three limit cycles can coexist in the open quadrant Ω. This result has been recently observed [26] but just in general cubic polynomial systems. In our case, only two of these limit cycles are concentric though, since infinitesimal ones are generated by Hopf bifurcation, and the third one is generated by a homoclinic bifurcation [23, 29] and surrounds another focus as stated in Lemma 6 and Theorem 8. Furthermore, the model can achieve the phenomenon of multistability by the existence of four ω-limit sets in the first quadrant as L < 0 and a strong Allee effect is present. A locally stable cycle surrounds a locally stable equilibrium point in the first quadrant as well as an unstable limit cycle, which serves as their separatrix in the phase plane; i.e., there is a range for population sizes for which there exists both autoregulation for the predator-prey system and prey and predator populations approach equilibrium, depending upon the population size. The third ω-limit set is the origin, stating that extinction is always possible in the presence of a strong Allee effect in our model. The fourth ω-limit set is a stable focus surrounded by the third limit cycle mentioned above. In order to illustrate the result stated in Theorem 8 and the multistability, Figure 10 shows a numerical simulation for system (5), which was affected with the software MATLAB [30]. Here a = 3, b = 1, k = 6.79211, M = 4.07697, c = 4.05116, and L = L+ + 0.000001 ≈ −0.498898 (see (12)). Initial conditions at (x(0), y(0)) are (0.8, 0.8), (1.01, 1.01), (2, 2), (2.2, 2.2), (1, 0.95), and (0.4, 0.1). In Figure 10 it is clear that a limit cycle lies between the orbits by the points (2, 2) and (2.2, 2.2), respectively, since the ω-limit of (2, 2) is the origin, but the ω-limit of (2.2, 2.2) is the stable focus located in the diagonal near x = 2.5. On the other hand, the two infinitesimal limit cycles are surrounding the singular point (1, 1), but they are too small for both the scale of this figure and these values of parameters. Hence, we present another numerical simulation in Figure 11 where both infinitesimal limit cycles are more visible. Here a = 3, b = 1, k = 6.56102, M = 4.16516, c = 4.16276, and L = L+ + 0.01 ≈ −0.480256 (see (12)). Initial conditions at (x(0), y(0)) are (0.6, 0.6), (0.8, 0.8), and (0.95, 0.95). Behaviors of the three positive oriented orbits in Figure 11 imply that, between two of them, attracting and repelling limit cycles exist. Furthermore, it can be easily verified that, for these values of parameters, the divergence of vector field (5) in singu-
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
1260
4 3.5 3
y
2.5 2 1.5 1 0.5 0 0
1
2
3
4
5
x Fig. 10. Numerical simulation of the results of Theorem 8.
1.5 1.4 1.3
y
1.2 1.1 1 0.9 0.8 0.7 0.6 0.4
0.6
0.8
1
1.2
1.4
1.6
x
Fig. 11. Numerical simulation of the infinitesimal limit cycles.
larity (1, 1) is negative, and, consequently, this equilibrium point is a local attractor and the innermost limit cycle is the unstable one. Several natural predator-prey communities have been studied, each one of them featuring an ecologically stable cycle, that is, a periodic orbit that must be somewhat insensitive to the perturbations outer to the interaction. In this work, we have answered partially one of the almost impossible projects proposed in [11] which is the following: Find a predator-prey or other interacting system in nature, or construct one in the laboratory, with at least two ecologically stable cycles. This shows the ecological relevance of the existence of multiple limit cycles in predation models and the importance of our result, which should serve for the outlined problem to be actually feasible in a biological lab with appropriate little creatures [11]. REFERENCES ´lez-Olivares, and E. Sa ´ez, Two limit cycles in a Leslie–Gower [1] P. Aguirre, E. Gonza predator-prey model with additive Allee effect, Nonlinear Anal. Real World Appl., to appear. [2] A. Andronov, E. Leontovich, I. Gordon, and A. Maier, Theory of Bifurcations of Dynamical Systems on a Plane, Israel Program for Scientific Translations, Jerusalem, 1971.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
THREE LIMIT CYCLES IN A LESLIE–GOWER
1261
[3] D. Arrowsmith and C. Place, An Introduction to Dynamical Systems, Cambridge University Press, New York, 1990. [4] D. K. Arrowsmith and C. M. Place, Dynamical Systems. Differential Equations, Maps and Chaotic Behaviour, Chapman and Hall, London, 1992. [5] J. Bascompte, Extinction thresholds: Insights from simple models, Ann. Zool. Fennici, 40 (2003), pp. 99–114. [6] L. Berec, E. Angulo, and F. Courchamp, Multiple Allee effects and population management, Trends Ecology Evol., 22 (2007), pp. 185–191. [7] T. R. Blows and N. G. Lloyd, The number of limit cycles of certain polynomial differential equations, Proc. Roy. Soc. Edinburgh Sect. A, 98 (1984), pp. 215–239. [8] D. S. Boukal and L. Berec, Single-species models and the Allee effect: Extinction boundaries, sex ratios and mate encounters, J. Theoret. Biol., 218 (2002), pp. 375–394. [9] C. W. Clark, Mathematical Bioeconomic. The Optimal Management of Renewable Resources, 2nd ed., John Wiley and Sons, New York, 1990. [10] C. W. Clark, The Worldwide Crisis in Fisheries: Economic Models and Human Behavior, Cambridge University Press, New York, 2006. [11] C. S. Coleman, Hilbert’s 16th problem: How many cycles? in Differential Equation Models, Modules Appl. Math. 1, M. Braun, C. S. Coleman, D. A. Drew, and W. F. Lucas, eds., Springe, New York, 1983. [12] J. B. Collings, Bifurcation and stability analysis of a temperature-dependent mite predatorprey interaction model incorporating a prey refuge, Bull. Math. Biol., 57 (1995), pp. 63–76. [13] F. Courchamp, T. Clutton-Brock, and B. Grenfell, Inverse dependence and the Allee effect, Trends Ecology Evol., 14 (1999), pp. 405–410. [14] F. Courchamp, B. T. Grenfell, and T. H. Clutton-Brock, Impact of natural enemies on obligately cooperative breeders, Oikos, 91 (2000), pp. 311–322. [15] A. M. De Roos and L. Persson, Size-dependent life-history traits promote catastrophic collapses of top predators, Proc. Natl. Acad. Sci. USA, 99 (2002), pp. 12907–12912. [16] B. Dennis, Allee effects: Population growth, critical density and the chance of extinction, Natur. Resource Modeling, 3 (1989), pp. 481–538. [17] F. Dumortier, Singularities of Vector Fields, Monogr. Mat. 32, IMPA, Rio de Janeiro, 1978. [18] R. Etienne, B. Wertheim, L. Hemerik, P. Schneider, and J. Powell, The interaction between dispersal, the Allee effect and scramble competition affects population dynamics, Ecol. Model., 148 (2002), pp. 153–168. ´lez-Yan ˜ez, and E. Gonza ´lez-Olivares, Conse[19] J. D. Flores, J. Mena-Lorca, B. Gonza quences of depensation in a Smith’s bioeconomic model for open-access fishery, in Proceedings of the 2006 International Symposium on Mathematical and Computational Biology BIOMAT 2006, Manaus, Brazil, 2007, pp. 219–232. [20] H. I. Freedman and G. S. K. Wolkowicz, Predator-prey systems with group defence: The paradox of enrichment revisted, Bull. Math. Biol., 8 (1986), pp. 493–508. [21] J. Gascoigne and R. N. Lipcius, Allee effects in marine systems, Marine Ecology Prog. Ser., 269 (2004), pp. 49–59. [22] J. Gascoigne and R. N. Lipcius, Allee effects driven by predation, J. Appl. Ecology, 41 (2004), pp. 801–810. [23] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, New York, 1983. [24] I. Hanski, L. Hansson, and H. Henttonen, Specialist predators, generalist predators and the microtine rodent cycle, J. Animal Ecology, 60 (1991), pp. 353–367. [25] I. Hanski, H. Hentonnen, E. Korpimaki, L. Oksanen, and P. Turchin, Small-rodent dynamics and predation, Ecology, 82 (2001), pp. 1505–1520. [26] X. Huang, Y. Wang, and L. Zhu, One and three limit cycles in a cubic predator-prey system, Math. Methods Appl. Sci., 30 (2007), pp. 501–511. [27] A. Korobeinikov, A Lyapunov function for Leslie-Gower predator-prey models, Appl. Math. Lett., 14 (2001), pp. 697–699. [28] Y. Li and D. Xiao, Bifurcations of a predator-prey system of Holling and Leslie types, Chaos Solitons Fractals, 34 (2007), pp. 606–620. [29] D. Luo, X. Wang, D. Zhu, and M. Han, Bifurcation Theory and Methods of Dynamical Systems, Adv. Ser. Dyn. Syst. 15, World Scientific, Singapore, 1997. [30] MATLAB: The Language of Technical Computing, Using MATLAB (Version 7.0), MathWorks, Natwick, MA, 2004. [31] R. M. May, Stability and Complexity in Model Ecosystems, 2nd ed., Princeton University Press, Princeton, NJ, 2001.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1262
´ ´ P. AGUIRRE, E. GONZALEZ-OLIVARES, AND E. SAEZ
´lez-Olivares, and B. Gonza ´lez-Yan ˜ez, The Leslie-Gower [32] J. Mena-Lorca, E. Gonza predator-prey model with Allee effect on prey: A simple model with a rich and interesting dynamics, in Proceedings of the 2006 International Symposium on Mathematical and Computational Biology BIOMAT 2006, Manaus, Brazil, 2007, pp. 105–132. ´ez and E. Gonza ´lez-Olivares, Dynamics of a predator-prey model, SIAM J. Appl. [33] E. Sa Math., 59 (1999), pp. 1867–1878. [34] P. A. Stephens and W. J. Sutherland, Consequences of the Allee effect for behaviour, ecology and conservation, Trends Ecology Evol., 14 (1999), pp. 401-405. [35] F. Takens, Unfoldings of certain singularities of vector fields: Generalized Hopf bifurcations, J. Differential Equations, 14 (1973), pp. 476–493. [36] R. J. Taylor, Predation, Chapman and Hall, New York, 1984. [37] P. Turchin, Complex Population Dynamics. A Theoretical/Empirical Synthesis, Monogr. Population Biol., 35, Princeton University Press, Princeton, NJ, 2003. [38] X. Wang, G. Liang, and F.-Z. Wang, The competitive dynamics of populations subject to an Allee effect, Ecol. Model., 124 (1999), pp. 130–168. [39] G. S. W. Wolkowicz, Bifurcation analysis of a predator-prey system involving group defense, SIAM J. Appl. Math., 48 (1988), pp. 592–606. [40] G. S. K. Wolkowicz, H. Zhu, and S. A. Campbell, Bifurcation analysis of a predator-prey system with nonmonotonic functional response, SIAM J. Appl. Math., 63 (2003), pp. 636– 682. [41] S. Wolfram, Mathematica: A System for Doing Mathematics by Computer, Addison–Wesley, Longman, Boston, 1988 [42] D. Xiao and S. Ruan, Global analysis in a predator-prey system with nonmonotonic functional response, SIAM J. Appl. Math., 61 (2001), pp. 1445–1472. [43] D. Xiao and S. Ruan, Bifurcations in a predator-prey system with group defense, Int. J. Bifurcation and Chaos, 11 (2001), pp. 2123–2131. [44] D. Xiao and H. Zhu, Multiple focus and Hopf bifurcations in a predator-prey system with nonmonotonic functional response, SIAM J. Appl. Math., 66 (2006), pp. 802–819. [45] D. Xiao and K. F. Zhang, Multiple bifurcations of a predator-prey system, Discrete Contin. Dyn. Syst. Ser. B, 8 (2007), pp. 417–433.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
