Relaxation oscillator profile of limit cycle in predator–prey system, Discrete

It is known that some predator-prey system can possess a unique 
limit cycle which is globally asymptotically stable. For a 
prototypical predator-prey system, we show that the solution curve 
of the limit cycle exhibits temporal patterns of a relaxation 
oscillator, or a Heaviside function, when certain parameter is 
small.


1.
Introduction. For a class of conventional predator-prey interaction models, it is known that a stable limit cycle exists for a range of parameters [19]. A typical model is where the prey U satisfies a logistic growth pattern; γ > 0 represents the intrinsic growth rate of the prey; K > 0 is the carrying capacity of the prey; D > 0 is the death rate of the predator; C > 0 measures the relative loss of the prey; the function φ(U ) is the functional response of the predator, which corresponds to saturation of their appetites and reproductive capacity, and like effects [7,19]. A functional response (of Type II [7]) usually satisfies φ(0) = 0, φ(U ) is increasing and concave, and φ(U ) → M > 0 for some M > 0 as U → ∞. Examples include φ(U ) = M U/(A + U ) (Holling) and φ(U ) = 1 − e −AU (Ivlev).
In this paper, we consider the Holling type II functional response. The system considered is

SZE-BI HSU AND JUNPING SHI
We introduce a change of variables: then we obtain a dimensionless equation: where m = M γ , d = D γ , and a = A K .
Here a, m, d > 0 are dimensionless parameters. Phase portrait analysis can show for certain parameters, a prey-only or coexistence equilibrium is globally stable (see Section 2 or [8]); and for other parameters, a periodic solution exists. It has been shown that for a class of systems including (4), the periodic solution is unique thus a globally stable limit cycle. The first such uniqueness result was proved by Cheng [2], and more general uniqueness results for limit cycle in predator-prey systems have been proved later in [13,14,16,27,30,31]. A main idea of later result is to transform (4) or a more general predator-prey system into a Liénard equation.
Our interest in this article is on the asymptotic behavior of the limit cycle of (4) when the predator death rate d tends to zero. A bifurcation point of view could ease the understanding of our result. If we fix other parameters in the system (4) so that 0 < a < 1, and take d as a bifurcation parameter, then the behavior of the system changes as the v-isocline mu a + u = d slides when d changes. It is more convenient to solve this v-isocline as u = λ ≡ ad m − d . When λ ≥ 1, the semi-trivial steady state (1, 0) is globally stable; when (1 − a)/2 < λ < 1 (u = λ intersects with the falling part of u-isocline), then the coexistence steady state is globally stable; and when 0 < λ < (1 − a)/2 (u = λ intersects with the rising part of u-isocline), then the limit cycle is globally stable. Notice that λ = (1 − a)/2 is the Hopf bifurcation point, where a subcritical Hopf bifurcation occurs, and a small amplitude periodic solutions emerges for λ < (1 − a)/2. Our main result in the article is on the limiting behavior of the unique limit cycle Σ λ when the death rate of predator d tends to zero (or equivalently λ tends to zero). When d is not very small, the periodic functions u(t) and v(t) are still sinusoidal-like (see Figure 1). Some sharp patterns emerge as d → 0 (or λ → 0). For small λ, we show that the period of Σ λ is in an order of O(λ −1 ); the prey population u(t) is low in order of O(λ) for a time scale of O(λ −1 ), then it has a spike to reach the maximum but only for a time scale of O(| ln λ|), hence the graph of u(t) is a periodic pulse; the predator population v(t) reaches the maximum value from the minimum value in a time scale of O(| ln λ|), then it slowly decays to the minimum value in a time scale of O(λ −1 ) and the decay is exponentially slow (see Figure 2). See Theorem 3.5 for a more mathematical description.
The phenomenon which we describe above makes the limit cycle of predator-prey system (4) behave similar to a nonlinear relaxation oscillator. Well-known examples of nonlinear relaxation oscillators are Van der Pol oscillator in electrical circuits employing vacuum tubes, Fitzhugh-Nagumo oscillator in action potentials of neurons (see [6,10,26,28,29].) The existence of relaxation dynamics in predator-prey model (4) seems to be first discovered in this article. For a two competing predators and one prey model considered in [11,12], it is known that stable relaxation oscillations exist for some parameter ranges by using singular perturbation methods [18,21]. In these work, it is assumed that the prey population has fast dynamics, i.e. the prey population grows much faster than those of the predators. In the current article, we assume that the predator has small death rate, and our method is totally different. Notice that from (5), fast prey growth rate (large γ) implies small m and d, and we only assume small d and fix m. Yet another example of singular perturbation in predator-prey system can be found in Deng et. al. [4].
In comparison we also consider the case as the parameter a tends to zero (λ also tends to zero in this case). The total period of the limit cycle also tends to infinity as a → 0 + . But the asymptotic profile of the limit cycle is quite different. In this limit, the predator population v(t) shows a spiky pulse shape, and the temporal length of the pulse is in a scale of O(| ln λ|); on the other hand, the prey population u(t) shows a profile of Heaviside function, with slow time scales O(λ −1 ) when u(t) ≈ 0 or u(t) ≈ 1 (the carrying capacity), connected by fast time scales of O(| ln λ|) between them (see Figure 3). See Theorem 4.2 for a more mathematical description.
The latter result provides further answer to an old question in ecological studies. In [24], Rosenzweig argues that enrichment of the environment (larger carrying capacity K in (2)) leads to destabilizing of the coexistence equilibrium, which is socalled paradox of enrichment. From (5), when other parameters are fixed, increasing K is necessarily equivalent to decreasing a. Our result shows that the time interval when the prey is population near zero is extremely long when the carrying capacity K is extremely large. That could make the prey population even more vulnerable to catastrophe perturbation with long time with very low population density.
Our result is rigorously proved by using basic differential and integral calculus, a Lyapunov function, and phase plane analysis. It is noteworthy that the orbit of the limit cycle in (4) does not follow the slow manifold as other nonlinear relaxation oscillators. It is well known that (4) can be converted to a generalized Liénard equation with a nontrivial transformation (see [16]), but the relaxation oscillation found here does not follow from known results for Liénard equations or Van der Pol equations. In fact, by using the change of variables: we can convert (4) to where a i (i = 1, 2, 3) are positive constants defined by The system (6) has a unique coexistence equilibrium point (x, y) = (1, y 0 ≡ a 1 + a 2 − a 3 ) and A further change of variables transform (6) into a generalized Liénard equation:

SZE-BI HSU AND JUNPING SHI
where φ(v) = y 0 (e v/y0 − 1), F (u) = a 2 − a 3 − a 2 e u + a 3 e 2u , and h(u) = y 0 (e u − 1). Using this form and a uniqueness result of limit cycle of Liénard equation by Zhang [31], one can prove the uniqueness of limit cycle of (4) (see [16]). But when d → 0, we have y 0 → ∞ and the profile of the limit cycle does not follow from any existing results. We point out that the relaxation oscillation property of Van der Pol system when ε → 0 has been studied by Liénard [17], Ponzo and Wax [22,23], Grasman [6]. More delicate limiting behavior of the limit cycle of the special system has been recently obtained by Dumortier and Roussarie [5], Krupa and Szmolyan [15] and others. We recall some well-known results regarding the dynamics of system (4) in Section 2, and we prove our main results in Section 3 and 4 for the case d → 0 and a → 0 respectively. We will use δ i and C i , (i ∈ N), to denote various positive constants. These constants are independent of d in Section 3, and are independent of a in Section 4.

SZE-BI HSU AND JUNPING SHI
The dynamics of (4) under (13) is completely understood. The local stability of (λ, v λ ) can be determined from the linearization at the equilibrium. We use λ as the bifurcation parameter. The Jacobian at (λ, v λ ) is Then is locally asymptotically stable. Indeed the local stability indeed implies the global asymptotical stability of (λ, v λ ) from the Poincaré-Bendixon theory and Dulac Theorem [11], and the global stability of (λ, v λ ) can also be proved through a mixed type Lyapunov function (see [1,3,9]). Finally when 0 < λ < 1 − a 2 , (λ, v λ ) is locally unstable, and (4) possesses a unique limit cycle which is globally asymptotically orbital stable (see [2,16]).
3. Asymptotic behavior of the limit cycle for d small. In the equation, we define In the first part we construct an invariant region where the limit cycle is located. For this part, we always assume that m, d > 0, 0 < a < 1, and λ = ad/(m − d) satisfies 0 < λ < (1 − a)/2. We first give an estimate of the unstable manifold U = {(u 1 (t), v 1 (t)) : t ∈ R} at the saddle point (1, 0). From the phase portrait, it satisfies 0 < u 1 (t) < 1 for all Proof. From the equation (4), we have Since the unstable manifold satisfies 0 < u 1 (t) < 1 for all t ∈ R, then along U , we Integrating along the portion of U from u = 1 to some u < λ, we obtain For the upper bound, we notice that the tangent line of the unstable manifold is Hence we only need to show that the vector field (f (u, v), g(u, v)) points towards the region below the That proves the upper bound v 1 (u) ≤ v 2 (u).
From Lemma 3.1, the unstable manifold reaches its maximum v-value when u = λ, and the maximum value v * can be estimated as From the phase portrait of the system, the limit cycle is below the unstable manifold U , then we also have the following upper bound for the location of limit cycle.
Then the orbit of the limit cycle In constructing a more precise region R 2 ⊂ R 1 containing Σ, we prove that for a sub-region R 3 containing (λ, v λ ), Σ ∩ R 3 = ∅. Define where W (u, v) is the function defined in (9). We notice that (1 − a − λ, v λ ) is the reflection of (λ, v λ ) with respect to the line u = (1 − a)/2. Such reflection technique is a key in proving the uniqueness of the limit cycle of (4) ( [2]). Lemma 3.3. Let R 3 be defined as in (20). Then R 3 is a bounded convex subset of It is easy to see that W 1 (u) is strictly decreasing in [0, λ) and is strictly increasing in (λ, ∞); and W 2 (v) is strictly decreasing in [0, v λ ) and is strictly increasing in (v λ , ∞  are both convex one-variable functions. For R 3 defined in (20), (1 − a − λ, v λ ) is the right-most point of R 3 . Thus for any solution orbit (u(t), v(t)) passing through points outwards. Hence from the properties of periodic orbit, Σ ∩ R 3 = ∅.
From Lemmas 3.2 and 3.3, we have obtained an invariant region R 2 where the limit cycle is located. Next we give some estimates for the extremal points on the orbit of limit cycle as d → 0 + . The other two parameters 0 < a < 1 and m > 0 are fixed, while λ = ad/(m − d) → 0 as d → 0 + . Hence d and λ are two equivalent parameters which tend to zero. Define u λ,− = min{u(t) : (u(t), v(t)) ∈ Σ}, u λ,+ = max{u(t) : (u(t), v(t)) ∈ Σ}, v λ,− = min{v(t) : (u(t), v(t)) ∈ Σ}, v λ,+ = max{v(t) : (u(t), v(t)) ∈ Σ}. (21) Notice that the both the upper and lower portions of the limit cycle are monotone functions, thus we define such that v − (λ, u) < v 0 (u) < v + (λ, u) for u λ,− < u < u λ,+ . That is, {(u, v + (λ, u))} is the upper portion of the limit cycle Σ, and {(u, v − (λ, u))} is the lower portion. From the equations, it is easy to see that u λ,− and u λ,+ are achieved when Σ intersects with the isocline v = v 0 (u), and v λ,− and v λ,+ are achieved when Σ intersects with the line u = λ. Our estimates are mainly based on the inner boundary of the region R 2 , i.e. the level curve and where The function h(x, b) satisfies and for fixed b > 0, h(·, b) achieves its global minimum 0 at x = b, and lim x→0 + h(x, b) = lim x→∞ h(x, b) = ∞. Thus for any b > 0, h(x, b) = c has exactly two roots for any c > 0.
To obtain the global asymptotical behavior of the limit cycle Σ, we divide the orbit with four reference points (see Figure 4): Let T = T (λ) be the period of Σ.
. We also assume that u(0) = λ and v(0) = v λ,+ , i.e. the orbit starts from the highest point of v(t). Our main result in this section is , v(t) : t ∈ R} be the orbit of the unique periodic solution of (4) when 0 < λ < (1 − a)/2. Assume that 0 < a < 1 and m > 0 are fixed, the extremal points of Σ are defined as in (21), and O i , T i (i = 1, 2, 3, 4) and the period T are defined as above. When λ > 0 is sufficiently small (or equivalently d > 0 is small), then there exist constants C 4 , C 5 > 0 independent of λ, such that C 5 λ −1 ≥ T ≥ C 4 λ −1 . Moreover, for λ > 0 sufficiently small, there exists some C 6 > 0, such that as λ → 0 + .

RELAXATION OSCILLATION IN PREDATOR-PREY SYSTEM 905
Proof. We prove the theorem in several steps.
Step 3: We show that For this portion, u(t) ≥ (1 − a)/2. From the equation of v, we have Hence which implies (42).
Step 4: We show there exist constants δ 5 , δ 6 , C 2 > 0 such that when 0 < λ < δ 6 , This is similar to Step 2. Now we have (44) Here the first inequality is from Lemma 3.4, and the second inequality is from that fact that v 0 (u) is increasing while v 5 (u) is decreasing in [λ, (1 − a)/2), and v 0 (u) < v 5 (u). Similar to Step 2, we obtain that when 0 < λ < δ 6 , The remaining part is same as Step 2.
From Lemma 3.1 and (18), we obtain the estimate of upper bound of v λ,+ by letting v 3 = (m + a + 1)/m. For the estimate of v 4 , we notice that any solution orbit satisfies Recall that O 1 = (λ, v λ,+ ) and O 2 = (λ, v λ,− ) are the highest and lowest points on the orbit of the limit cycle Σ. Let the leftmost point on Σ be O 5 = (u λ,− , v * ). Then from (48), we obtain that Since v 2 − δ 0 < v λ,+ < v 3 for small λ, then the right hand side of (50) is bounded. On the other hand, for the first integral in (49), Thus − ln v λ,− is bounded from above from (49), (50) and (51), and consequently v λ,− is bounded from below by some v 4 > 0 for all small λ > 0.
Step 7: The completion of the proof. From Lemma 3.4 and Step 6, when λ > 0 is small, v 4 < v λ,− < v 1 + δ 0 and v 2 − δ 0 < v λ,+ < v 3 where v 1 and v 2 are the two roots of h(v, a/m) = 1 − a such that v 1 < a/m < v 2 , and also lim λ→0 + λ −1 u λ,− = 0. Also from the definition of λ, d −1 = a + λ λm > a m λ −1 , and d −1 < a m (1 + δ 10 )λ −1 for any small δ 10 > 0 and we assume λ small. Thus from Step 1 and Step 5, for any 0 < δ 11 < 1, as long as λ > 0 is sufficiently small, Hence we obtain the estimate for T 1 in the theorem, since all constants except λ are independent of λ. The estimate for T 3 can also be obtained from Step 3 and  4. Asymptotic behavior of the limit cycle for a small. In this section, we assume that d and m are fixed so that m > d > 0, and a > 0 is small (thus a < 1). We will use a lot of estimates established in Section 3, and we will also use the same notations in Section 3 as well.