Limit cycles in the Holling-Tanner model

Publicacions Matem`atiques, Vol 41 (1997), 149–167. LIMIT CYCLES IN THE HOLLING-TANNER MODEL Armengol Gasull, Robert E. Kooij and Joan Torregrosa Abstract This paper deals with the following question: does the asymptotic stability of the positive equilibrium of the Holling-Tanner model imply it is also globally stable? We will show that the answer to this question is negative. The main tool we use is the computation of Poincar´e-Lyapunov constants in case a weak focus occurs. In this way we are able to construct an example with two limit cycles. 1. Introduction The two main types of interaction between any pair of biological species, which are of interest to the ecologist, are either when they are competing together for some common source of food supply, or when one of the species preys upon the other. In this paper we will restrict our attention to the latter case. The existence and the number of isolated periodic solutions (limit cycles) is one of the most delicate problems connected with two-dimensional predator-prey models. One of the first examples of a biological system modelling the interaction between prey and predators was formulated by Lotka in 1925 [11] and Volterra in 1927 [17]: (1.1) dx dt =αx −βxy, dy dt =−δy +γxy. In system (1.1) x(t) and y(t) denote prey and predator densities respectively, as functions of time. Furthermore, all constants are assumed to be positive. Obviously, the attention is restricted to x≥0, y≥0. 150 A. Gasull, R. E. Kooij, J. Torregrosa It is a well-known fact that (1.1) has a family of periodic orbits, but no limit cycles. Due to the fact that (1.1) has a center system (1.1) lacks structural stability. That is, the phase portrait of (1.1) can be changed if we take into account arbitrarily small additional effects. If, for example, (1.1) is modified to include the effect of competition among the prey (by adding −x2to the first equation of (1.1)) then the resulting system no longer has a center and the population oscillations decay. A first generalization of system (1.1) was suggested by Gause in 1934 [6]: (1.2) dx dt =αx −p(x)y, dy dt =−δy +γp(x)y. Here α>0 is the growth rate of the prey in absence of the predator; δ>0 is the death rate of the predator in absence of the prey; γ>0is the rate of conversion of consumed prey to predator. Finally, p(x) is the capture rate of prey per predator or functional response of a predator. For most examples that appear in the literature (see the bibliography in [5]) it is assumed that p(0) = 0 and p(x)>0 for all x>0. The generalized Gause model for the interaction of the two species is (see [5]): (1.3) dx dt =xg(x)−p(x)y, dy dt =−δy +h(x)y. System (1.3) incorporates density-dependent prey growth in absence of the predator. This is introduced in the model because it is quite unrealistic to assume that the prey will grow to infinity in absence of predators, as will happen for (1.1) and (1.2). The growth rate g(x) satisfies g(0) >0, g(x)<0 for all x>0 and there exists a K>0 such that g(K)=0. Kis called the carrying capacity of the prey. A growth rate of this type is thought to model the situation where the food supply for the prey is limited. For high densities of prey they will compete for the resources. A famous example that belongs to systems of type (1.3) is a system first mentioned by Rosenzweig and McArthur in 1963 [14]: (1.4) dx dt =rx1−x K−mx A+xy, dy dt =−δy +γmx A+xy, Limit cycles in the Holling-Tanner model 151 where r, K, m, A, δ and γare positive constants. System (1.4) is important because it is a structurally stable model which can exhibit persistent oscillations. The uniqueness of the limit cycle of system (1.4) was first proved by Cheng [3]. However, his arguments are rather tedious; a simplified proof was given by Kuang and Freedman [10]. The functional response that occurs in system (1.4), p(x)= mx A+x, was suggested by the biologist Holling in 1959 [7]. In fact, based on actual field data, he argued that the functional response should not only be a monotonically increasing, but also a bounded function. Holling suggested p(x)= mx A+x, referred to as a functional response of Holling type II, to represent invertebrate predators. In this function p(x), mis the maximum rate of predation while Acan be shown to be proportional to the time required for the predator to search and find a prey. A slightly different model was suggested by Tanner [16]: (1.5) dx dt =rx1−x K−mx A+xy, dy dt =sy 1−hy x. In the literature model (1.5) is referred to as the Holling-Tanner model. In (1.5) the predator grows logistically with intrinsic growth rate sand carrying capacity proportional to the size of the prey. The parameter his the number of prey required to support one predator at equilibrium. Clearly this model incorporates intraspecific competition among the predators. A study of several pairs of interacting species in [16] shows that the theoretical predictions of (1.5) based on estimated parameter values are broadly in line with practical reality. The local stability of the unique positive equilibrium of (1.5) was investigated by May [12] and Murray [13]. Recently Hsu and Huang [8] obtained some results on the global stability of the positive equilibrium. To be more precise, they obtained conditions under which local stability of the positive equilibrium implies its global stability. In fact, it is sometimes conjectured that for predator-prey systems with a unique positive equilibrium local and global stability are equivalent. For instance, in Arrowsmith and Place [2], it is remarked that this conjecture is true for system (1.5). However, we will show in this paper that this conjecture is not true for system (1.5). We do this by constructing an example where the positive equilibrium is stable while being surrounded by two limit cycles, the innermost being unstable and the outermost being stable. 152 A. Gasull, R. E. Kooij, J. Torregrosa This implies the interesting phenomenon of coexistence of a stable equilibrium and persistent oscillations. Our results confirm the numerical results obtained by Wollkind et al. [18], who used system (1.5) to model temperature-dependent mite interaction on fruit trees. The rest of the paper is organized as follows. In section 2 we discuss some general properties of (1.5). To this end we will transform (1.5) to a generalized Li´enard equation. This makes it possible to simplify the proofs of the results obtained in [8]. In section 3 we discuss the case where the linearization matrix at the positive equilibrium has a pair of pure imaginary roots. For this case nonlinear terms have to be taken into account in order to determine the local stability of the positive equilibrium. For this case the stability depends on the sign of the first non-vanishing Poincar´e-Lyapunov constant. For the effective computation of these constants and for the study of their sign the use of a system for mathematical computation like Maple is inevitable. In section 4 we discuss the phase portraits of system (1.5), specially when one or two limit cycles occur. In section 5 the analytic results are illustrated by numerical experiments. Finally we state some open problems for system (1.5). 2. General properties of system (1.5) After applying the rescaling t=rt,x=x K,y=my rK system (1.5) becomes (after omitting the bars): (2.1) dx dt =x(1 −x)−x a+xy, dy dt =yδ−βy x, where a, β, δ > 0. To study the phase portrait we will first investigate the singular points of (2.1). It is easy to see that S(1,0) is a saddle and the character of O(0,0) can be obtained after rescaling the time in (2.1) by dt dτ =x: (2.2) dx dτ =x2(1 −x)−x2 a+xy, dy dτ =y(δx −βy). For the system (2.2) the character of the origin can be obtained by applying one blow up. It appears that the character of this singularity depends on the sign of δ−1, see Figure 2.1. Limit cycles in the Holling-Tanner model 153 (a) 0 <δ<1(b) δ>1 > > > > > > > >> Figure 2.1. Phase portrait of system (2.2) near the origin. The behaviour of the trajectories near the origin for system (2.1) for x>0 is the same as in Figure 2.1. We will see later that we are only interested in the case 0 <δ<1 because for δ≥1 system (2.1) has no limit cycles. Furthermore there is a positive equilibrium E(x∗,y∗), where 0<x∗<1 is the unique positive zero of (2.3) δx −β(1 −x)(a+x)=0. It is easy to show that Eis an antisaddle i.e. the product of the eigenvalues of the linearization matrix at Eis positive. The following lemma, which is proved in [8, Lemma 2.1], will appear to be useful in the sequel. Lemma 2.1. The solutions of (2.1) are positive and bounded, and furthermore there exists T≥0such that x(t)<1,y(t)<δ βfor t≥T. The following lemma shows that system (2.1) can be transformed to a generalized Li´enard system. This is useful because there are many criteria available which guarantee nonexistence or uniqueness of limit cycles for (generalized) Li´enard equations, see for instance Ye et al. [19]. Lemma 2.2. System (2.1) can be transformed to a generalized Li´enard system of the form (2.4) dx dτ =ϕ(z)−x 1 f(τ)dτ, dz dτ =−g(x), 154 A. Gasull, R. E. Kooij, J. Torregrosa where ϕ(z)=ez, f(x)=f(x)x−2−βexp aβ x, g(x)=g(x)x−2−βexp aβ x, with, f(x)=(1−a)x2−2x3−β(1 −x)(a+x)2, g(x)=(δx −β(1 −x)(a+x))(a+x). Proof: First rewrite system (2.1) as (2.5) dx dt =f0(x)−f1(x)y, dy dt =g0(x)+g1(x)y+g2(x)y2, where, f0(x)=x(1 −x),g 0(x)=0, f1(x)= x a+x,g 1(x)=δ, g2(x)=−β x. It is well-known that systems of the form (2.5) can be transformed toaLi´enard equation, see for instance [19]. However, if we know a solution y=ψ(x) of (2.5) then we can transform (2.5) to a simpler generalized Li´enard equation, see Zegeling and Kooij [20]. Because y≡0 is a solution of (2.1) we can apply the transformation y=ω(x)ez,dt dτ =−1 f1(x)ω(x), where ω(x) satisfies ω(x) ω(x)=−g2(x) f1(x)and hence ω(x)=xβexp −aβ x, see Lemma 2.2 in [20], to obtain the generalized Li´enard equation (2.4). Note that it follows from this transformation that F(x)=x 1 f(τ)dτ =f0 f1ω(x)=(1−x)(a+x)x−βexp aβ x and from this we can obtain f(x)= d dx (F(x)). This completes the proof. Limit cycles in the Holling-Tanner model 155 The equilibrium E(x∗,y∗) of system (2.1) is mapped unto the equilibrium D(x∗,z∗) of system (2.4), where z∗satisfies y∗=ω(x∗)ez∗. Its stability depends on the sign of −f(x∗) which equals the sign of −f(x∗). If we solve βfrom δx∗−β(1 −x∗)(a+x∗) = 0 and substitute this into f(x) then we obtain f(x∗)=x∗(−2x∗2+(1−a−δ)x∗−aδ)=−x∗P(x∗). This defines a function (2.6) P(x)=2x2+(a+δ−1)x+aδ, which is also derived in [8, formula (3.3)]. If P(x) has a fixed sign (i.e. ≥0) then there are no limit cycles. This is shown in [8, Theorem 3.2 (i)], by using a Dulac function x a+x−1y−2 for system (2.1). It also can be deduced by using a Dulac function e−y for system (2.4) and the fact that f(x)−g(x)=−xP(x). Hence it follows directly that for δ≥1 system (2.1) has no limit cycles. Therefore the only case that interests us is where P(x) has two positive roots: 0 <α 1<α 2<1. In the aδ-plane this means we restrict our attention to the region D, defined by (2.7) 0<δ<1+3a−22a+2a2, 0<a<1, see Figure 2.2. 1 0.8 0.6 0.4 0.2 0D a 00.20.40.60.81 δ Figure 2.2. The region D. 156 A. Gasull, R. E. Kooij, J. Torregrosa For this case the function F(x) has at most two extrema. We will assume that F(x) has exactly two extrema, β1and β2, where 0 <β 1< β2<1, see Figure 2.3, because if F(x) is monotone then f(x) has a fixed sign and then the nonexistence of limit cycles follows by Bendixson’s criterion. β1 ˜xβ21x F(x) Figure 2.3. The function F(x). As mentioned before the local stability of E(x∗,y∗) depends on the sign of −f(x∗). Therefore we obtain immediately, see also [8, Lemma 3.1]: Lemma 2.3. (a) If 0<x ∗<α 1or α2<x ∗<1, then E(x∗,y∗)is asymptotically stable; (b) If α1<x ∗<α 2then E(x∗,y∗)is unstable. In [8, Theorem 3.2], is proved: Lemma 2.4. (a) If α2<x ∗<1then E(x∗,y∗)is globally stable; (b) If x∗is sufficiently small (to be precise, if x∗<˜x, with F(˜x)= F(β2), see Figure 2.3) then E(x∗,y∗)is globally stable. These results can be obtained by showing that under the conditions (a) and (b) no limit cycles surround E(x∗,y∗). Then the result follows from Lemma 2.1. The proof given by [8], who use different transformations for case (a) and (b), can be unified by studying the generalized Li´enard equation (2.4). We leave this as an exercise to the reader. Because for α1<x ∗<α 2,E(x∗,y∗) is unstable it follows by Lemma 2.1 that we can apply the Poincar´e-Bendixson theorem to deduce the existence of at least one stable limit cycle for this case. Limit cycles in the Holling-Tanner model 157 For x∗=α1or x∗=α2the eigenvalues of the linearization matrix of D(x∗,z∗) (and hence at E(x∗,y∗)) are purely imaginary so we have to take into account higher order terms in order to determine whether this equilibrium is a center or a weak focus. Recall that a weak focus is a singularity that is a center for the linearized system but not necessarily for the nonlinear system. If the origin is a weak focus then the canonical form of the system reads: (2.8) dx dt =−y+F2(x, y), dy dt =x+G2(x, y), where F2and G2denote terms of at least order two. It is known (see [1]) that if Π denotes the Poincar´e return map defined by the flow of (2.8) on the positive x-axis at a neigbourhood of zero, then either Π(x)≡x(the origin of (2.8) is a center) or there exists k∈N and Vk= 0 such that Π(x)=x+Vkx2k+1 +O(x2k+2). The value Vk is called the Poincar´e-Lyapunov constant of the origin of (2.8) and the stability of this point is determined by its sign. In this situation we will say that the order of the weak focus is k. Under perturbation of the coefficients of (2.8) at most klimit cycles can bifurcate out of a weak focus of order k. Such limit cycles are said to be of small amplitude. For an ample discussion on the computation of Poincar´e-Lyapunov constants we refer to Andronov et al. [1]. In the next section we will compute the Poincar´e-Lyapunov constants for system (2.1) when x∗=α1or x∗=α2. However, we will first state a result on the stability of E(x∗,y∗) in the case that x∗=α2, without making an actual computation. Theorem 2.1. If x∗=α2then E(x∗,y∗)is asymptotically stable. Proof: First suppose that for x∗=α2,E(x∗,y∗) is locally unstable i.e. the first non-vanishing Poincar´e-Lyapunov constant is positive. Then by Lemma 2.1 and the Poincar´e-Bendixson theorem E(x∗,y∗)is surrounded by at least one stable limit cycle. Then, by continuity, system (2.1) still would have a limit cycle if 0 <x ∗−α21, but this contradicts Lemma 2.4 (a). We can exclude the possibility that all Poincar´eLyapunov constants are zero because this would mean that E(x∗,y∗)isa center i.e. a singularity surrounded by a family of closed orbits. However, by Lemma 2.1 it would follow that this family has an outermost closed orbit without singularities on it. But this is known to be impossible for 164 A. Gasull, R. E. Kooij, J. Torregrosa Problem 5.2. If 0<x ∗≤α1and V1≤0then E(x∗,y∗)is globally stable. Note that this problem is already proved for 0 <x ∗≤˜x, where F(˜x)= F(β2), see Lemma 2.4 (b). In Figure 5.3 we have obtained some trajectories of system (2.1) numerically for the parameter values a=.2, β=.0682, δ=.2 for which ˜x<x ∗≤α1and V1<0. 00.20.40.60.81 0 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 Holling-Tanner model: a=.2, β=.682, δ=.2 Figure 5.3. System (2.1) with a=.2, β=.0682, δ=.2. The graph of the function f(x) g(x)for this case is depicted in Figure 5.4. 0.20.1 0 −1 x y 1 00.30.40.5 −2 −3 2 3 Figure 5.4. f(x) g(x)with a=.2, β=.0682, δ=.2. Limit cycles in the Holling-Tanner model 165 Because for every c∈Rthe equation f(x) g(x)=chas a solution, it follows that it is not possible to use a Dulac function of the form e−cy for the generalized Li´enard equation (2.4) to prove nonexistence of limit cycles in this case. Problem 5.3. If 0<x ∗<α 1and V1>0then E(x∗,y∗)is surrounded by at most two limit cycles. This is by far the most difficult of the three open problems because there are no theorems available which guarantee the existence of at most two limit cycles for a general system. Maybe one can obtain some results in case it is assumed that some of the parameters in the system are small. This would lead to the study of the zeros of Abelian integrals. This approach is followed by Rothe and Shafer [15]. A numerical example with two limit cycles is displayed in Figure 5.5. The parameter values are a=.1, β=.005, δ=.1, such that ˜x<x ∗<α 1 and V1>0. 00.20.40.60.81 0 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 Holling-Tanner model: a=.1, β=.005, δ=.1 Figure 5.5. System (2.1) with a=.1, β=.005, δ=.1. References 1. A. A. Andronov et al.,“Theory of bifurcation of dynamic systems on a plane,” translated by John Wiley and Sons, New York, 1973. 166 A. Gasull, R. E. Kooij, J. Torregrosa 2. D. K. Arrowsmith and C. M. Place,“Dynamical Systems,” Chapman and Hall, London, 1992. 3. Cheng Kuosheng, Uniqueness of a limit cycle for a predator-prey system, SIAM J. Math. Anal. 12 (1981), 541–548. 4. A. Cima, A. Gasull, V. Ma˜ nosa and F. Ma˜ nosas, Algebraic properties of the Lyapunov and period constants, Rocky Mountain J. Math. 27 (1997) (to appear). 5. H. I. Freedman,“Deterministic mathematical models in population ecology,” Marcel Dekker, New York, 1980. 6. G. F. Gause,“The struggle for existence,” Williams and Wilkins, Baltimore, 1934. 7. C. S. Holling, The components of predation as revealed by a study of small-mammal predation of the European pine sawfly, Can. Entomol. 91 (1959), 293–320. 8. Sze-bi Hsu and Tzy-wei Huang, Global stability for a class of predator-prey systems, SIAM J. Appl. Math. 55(3) (1995), 763–783. 9. R. E. Kooij and A. Zegeling, Qualitative properties of two-dimensional predator-prey systems, J. Nonlinear Anal. TMA (1996) (to appear). 10. Kuang Yang and H. I. Freedman, Uniqueness of limit cycles in Gause type models of predator-prey systems, Math. Biosci. 88 (1988), 67–84. 11. A. J. Lotka,“Elements of physical biology,” Williams and Wilkins, Baltimore, 1925. 12. R. M. May,“Stability and complexity in model ecosystems,” Princeton University Press, Princeton, N.J., 1974. 13. J. D. Murray,“Mathematical Biology,” Springer-Verlag, Berlin, 1989. 14. M. L. Rosenzweig and R. H. McArthur, Graphical representation and stability conditions of predator-prey interactions, Am. Nat. 47 (1963), 209–223. 15. F. Rothe and D. S. Shafer, Multiple bifurcation in a predator-prey system with non-monotonic predator response, Proc. Roy. Soc. Edinburgh 120A (1992), 313–347. 16. J. T. Tanner, The stability and the intrinsic growth rates of prey and predator populations, Ecology 56 (1975), 855–867. 17. V. Volterra, Variazioni e fluttuazioni del numero d’individui in specie animali conviventi, in Italian, Mem. R. Com. Tolassogr. Ital. 131 (1927), 1–142. Limit cycles in the Holling-Tanner model 167 18. D. J. Wollkind, J. B. Collings and J. A. Logan, Metastability in a temperature-dependent model system for predator-prey mite outbreak interactions on fruit trees, Bull. Math. Biol. 50(4) (1988), 379–409. 19. Ye Yanqian et al.,“Theory of limit cycles,” Transl. of Math. Monographs 66, AMS, 1986. 20. A. Zegeling and R. E. Kooij, Uniqueness of limit cycles in polynomial systems with algebraic invariants, Bull. Austral. Math. Soc. 49 (1994), 7–20. 21. Zhang Zhifen, On the uniqueness of limit cycles of certain equations of nonlinear oscillations, in Russian, Dokl. Akad. Nauk SSSR 119 (1958), 659–662. 22. Zhang Zhifen, Proof of the uniqueness theorem of limit cycles of generalized Li´enard equations, Appl. Anal. 23 (1986), 63–76. Armengol Gasull: Departament de Matem`atiques Universitat Aut`onoma de Barcelona 08193 Bellaterra (Barcelona) SPAIN e-mail: [email protected] Robert E. Kooij: University of Technology Delft Fac. of Tech. Mathematics & Informatics Mekelweg 4 2628 CD Delft THE NETHERLANDS e-mail: ko[email protected]wi.tudelft.nl Joan Torregrosa: Departament de Matem`atiques Universitat Aut`onoma de Barcelona 08193 Bellaterra (Barcelona) SPAIN e-mail: [email protected] Rebut el 30 de Novembre de 1996