STABILITY AND BIFURCATION IN A DELAYED RATIO-DEPENDENT PREDATOR-PREY SYSTEM

Dongmei Xiao & Wenxia Li ^{1} Department of Mathematics

Central China Normal University, Wuhan 430079, P.R. China E-mail: dmxiao@ccnu-1.ccnu.edu.cn; wxli@mail.ccnu.edu.cn

Abstract. Recently, ratio-dependent predator-prey systems have been regarded by some researchers to be more appropriate for predator-prey interactions where predation involves serious searching processes. However, such models with delay have set up a challenging issue regarding their dynamics. In this paper, the qualitative behavior of a class of ratio-dependent predator-prey system with delay at the interior equilibrium in the interior of the first quadrant is studied. It is shown that the interior equilibrium cannot be absolutely stable and there exist nontrivial periodic solutions for the model. Moreover, by choosing delay τ as bifurcation parameter we study the Hopf bifurcation and the stability of the periodic solutions.

AMS Subject Classifications: Primary 34C25, 92D25; Secondary 58F14.

Keywords: Stability, bifurcation, Predator-prey system, ratio-dependent response, time delay.

1Research supported by the National Natural Science Foundations of China (No.10071027) and the Laboratory of Nonlinear Analysis .

2 STABILITY AND BIFURCATION IN A PREDATOR-PREY SYSTEM

1. Introduction

In population dynamics, the models most often used by ecologists to describe interactions between predator and prey populations are variations of the standard Lotka-Volterra type models. The standard Lotka-Volterra model is built by assuming that the per capita rate predation depends on the prey numbers only. Recently the traditional prey-dependent predator-prey models have been challenged by several biologists (see Arditi and Ginzburg [2], Arditi , Ginzburg and Akcakaya [3], Akcakaya [1], Gutierrez [12], and the references therein) based on the fact that functional and numerical responses over typical ecological timescales ought to depend on the densities of both prey and predators, especially when predators have to share or compete for food. Such a functional response is called a ratio-dependent response function, which can be roughly stated as that the per capita predator growth rate should be a function of the ratio of prey to predator abundance. Thus an alternative assumption is that, as the numbers of predators change slowly (relative to prey change), there is often competition among the predators, and the capita rate of predation depends on the numbers of both prey and predator, most likely and simply on their ratio. Based on the Michaelis-Menten or Holling type II function, Arditi and Ginzburg [2] proposed a ratio-dependent function of the form

p (\frac{x}{y}) = \frac{c\frac{x}{y}}{m + \frac{x}{y}} = \frac{cx}{my + x}

and the following ratio-dependent predator-prey model

\begin{aligned}ẋ & = x(a − bx) − \frac{cxy}{my + x}, \\ ẏ & = y (−d + \frac{fx}{my + x}) . \end{aligned} ; (1.1)

Here x(t) and y(t) represent population densities of prey and predator at time t, respectively; \frac{a}{b} > 0 is the carrying capacity of the prey, d > 0 is the death rate of the predator, and a, c, m and f are positive constants that stand for prey intrinsic growth rate, capturing rate, half saturation constant and conversion rate, respectively.

As it is typical for the predator-prey systems, the x−axis, y−axis and the interior of the first quadrant are all invariant under system (1.1), and solutions with positive initial values are positive and bounded. However differing from the prey-dependent

STABILITY AND BIFURCATION IN A PREDATOR-PREY SYSTEM 3

predator-prey models, the ratio-dependent predator-prey systems have two principal predictions: (a) equilibrium abundances are positively correlated along a gradient of enrichment (see Arditi and Gurzburg [2]) and (b) the “paradox of enrichment” (see Rosenzweig [18]) either completely disappears or enrichment is linked to stability in a more complex way. The ratio-dependent predator-prey model (1.1) has been studied by many researchers recently and very rich dynamics have been observed (see Jost, Arino and Aridti [12], Kuang [14], Kuang and Beretta [15], Xiao and Ruan [19] etc.). On the other hand, the effect of time delays on stability of the ratio-dependent predator-prey system has been investigated by Beretta and Kuang . In [6] Beretta and Kuang proposed a ratio-dependent model with a single discrete positive delay τ

\begin{aligned}ẋ(t) & = x(t)(a − bx(t)) − \frac{cx(t)y(t)}{my(t) + x(t)}, \\ ẏ(t) & = y(t) (−d + \frac{fx(t − τ)}{my(t − τ) + x(t − τ)}) . \end{aligned} ; (1.2)

System (1.2) satisfies the theorem of existence and uniqueness of solutions under the initial value condition

(1.3) x_{0}(θ) = ϕ_{1}(θ) ≥ 0, y_{0}(θ) = ϕ_{2}(θ) ≥ 0, θ ∈ [−τ, 0], x_{0}(0) > 0, y_{0}(0) > 0,

where ϕ = (ϕ_{1}, ϕ_{2}) ∈ C([−τ, 0], R_{+}^{2}), R_{+}^{2} = {(x, y) : x ≥ 0, y ≥ 0} and ||ϕ|| = max{|ϕ(θ)| : θ ∈ [−τ, 0]} with |ϕ| any norm in R^{2}.

Beretta and Kuang in [6] made use of a rather novel and nontrivial way of constructing proper Liapunov functions to obtain some new and significant global stability or convergence results for system (1.2). However it is far from completely answering some subtle mathematical questions on the behavior of solutions of the model, for example, the existence of periodic solutions for the model.

In this paper, we pay our attention on the interior equilibrium of system (1.2). We first study the effect of time delay on local stability of the interior equilibrium, then find out conditions on the delay and parameters of system (1.2) so that the interior equilibrium is conditional stable or unstable. And it will be shown that the interior equilibrium of system (1.2) cannot be absolutely stable for any parameters.

It is well-known that the oscillatory phenomenon of the population densities is a rather characteristic behavior of predator-prey dynamics. We then apply the

4 STABILITY AND BIFURCATION IN A PREDATOR-PREY SYSTEM

method in [10] to analyze Hopf bifurcation of system (1.2) by choosing delay τ as a bifurcation parameter, and show the existence and stability of periodic solutions to system (1.2).

This paper is organized as follows. In the next section we show that delay does cause instability, present results on non-occurrence of absolute stability and on conditional stability of the interior equilibrium. In section 3, the analysis of singularities is performed by using the normal form theory for functional differential equations in [10] and choosing delay τ as bifurcation parameter. We show that when the delay takes some critical values the interior equilibrium becomes unstable and Hopf bifurcation occurs, i.e., periodic solutions bifurcate from the interior equilibrium as delay τ passes through the critical values. We give the conditions which guarantee the stability of the periodic solution and use an example to illustrate the result. The paper ends with a brief discussion.

2. Local stability analysis

In this section, we focus on investigating the local stability of interior equilibrium of system (1.2). As shown in [16], system (1.2)(or (1.1)) has a unique interior equilibrium if and only if any one of the following two conditions is true:

(C1) d < f and c ≤ ma;

(C2) d < f < \frac{cd}{c−am} and am < c. We denote this unique interior equilibrium by E^{∗} = (x^{∗}, y^{∗}) and use this notation throughout this paper, here

x^{∗} = \frac{a}{b} − \frac{c(f − d)}{bmf}, y^{∗} = \frac{f − d}{dm}x^{∗}.

To study the local stability of the interior equilibrium E^{∗}, we consider linearization of system (1.2) at E^{∗}(x^{∗}, y^{∗})

\begin{aligned}Ẋ(t) & = α_{1}X(t) + α_{2}Y (t), \\ Ẏ (t) & = β_{1}X(t − τ) + β_{2}Y (t − τ), \end{aligned} ; (2.1)

where

α_{1} = −bx^{∗} + \frac{cx^{∗}y^{∗}}{(my^{∗} + x^{∗})^{2}} , α_{2} = −\frac{cx^{∗}^{2}}{(my^{∗} + x^{∗})^{2}} ,

β_{1} = \frac{fmy^{∗}^{2}}{(my^{∗} + x^{∗})^{2}} , β_{2} = −\frac{fmx^{∗}y^{∗}}{(my^{∗} + x^{∗})^{2}} .

STABILITY AND BIFURCATION IN A PREDATOR-PREY SYSTEM 5

The characteristic equation for (2.1) takes the form

(2.2)  λ^{2} − α_{1}λ − β_{2}λe^{−λτ} + (α_{1}β_{2} − α_{2}β_{1})e^{−λτ} = 0,

where α_{1}β_{2} − α_{2}β_{1} > 0.

We recall some definitions (see [5], [8] and [14]). The equilibrium (x^{∗}, y^{∗}) is called asymptotically stable if there exists a δ > 0 such that

sup_{−τ≤θ≤0}[|ϕ_{1}(θ) − x^{∗}| + |ϕ_{2}(θ) − y^{∗}|] < δ

implies that

lim_{t→∞}(x(t), y(t)) = (x^{∗}, y^{∗}),

where (x(t), y(t)) is the solution of system (1.2) with initial values (1.3). (x^{∗}, y^{∗}) is absolutely stable if it is asymptotically stable for all delay τ ≥ 0 and is conditionally stable if it is asymptotically stable for τ in some finite intervals. (x^{∗}, y^{∗}) is absolutely unstable if it is unstable for all delay τ ≥ 0.

It is known that the equilibrium is asymptotically stable if all roots of the corresponding characteristic equation have negative real parts. For the equilibrium E^{∗} of system (1.2) there are two types of stability: absolute stability (independent of the delay) and conditional stability (depending on the delay). The equilibrium E^{∗} of system (1.2) is absolutely stable if and only if this equilibrium of corresponding ODE system (i.e. system (1.1)) is asymptotically stable and the characteristic equation (2.2) has no purely imaginary roots for any τ > 0. Thus if all roots of characteristic equation (2.2) have negative real parts as τ = 0 and there exist some positive values τ such that the characteristic equation (2.2) has a pair of purely imaginary roots, say ±iω_{0}, then the equilibrium E^{∗} of system (1.2) is not absolutely stable but could be conditionally stable. Suppose that ω_{0} is obtained when the delay reaches a positive value τ_{0}, here τ_{0} satisfies that when 0 < τ < τ_{0} the real parts of all roots of the characteristic equation still remain negative, and when τ = τ_{0}, the characteristic equation (2.2) has a pair of purely imaginary roots ±iω_{0}. If moreover the transversal condition holds at τ = τ_{0}, then when τ > τ_{0} the characteristic equation (2.2) will have at least one root with positive real part and the equilibrium E^{∗} of system (1.2) becomes unstable by Rouché’s theorem (Dieudonné [9]). Hence, the equilibrium E^{∗} of system (1.2) is conditionally stable.

6 STABILITY AND BIFURCATION IN A PREDATOR-PREY SYSTEM

In order to study local stability of system (1.2), we have to recall some results of the corresponding ODE system (i.e. system (1.1)) on stability of the interior equilibrium as the following lemma (the proof can be found in [19]).

Lemma 2.1. Suppose that system (1.1) has a unique interior equilibrium (x^{∗}, y^{∗}), i.e., any one of conditions (C1) and (C2) holds. Denote Δ = m^{2}d^{4} +4cd^{2}(c−am− dm). Then there are only three possibilities.

(1) the equilibrium (x^{∗}, y^{∗}) of system (1.1) is locally asymptotically stable if any one of the following conditions holds:

(i) d < f and c ≤ ma; (ii) d < f < \frac{cd}{c−am} and am < c ≤ am + dm; (iii) d < f < \frac{−md^{2}+\sqrt{Δ}}{2(c−am−dm)} and 0 < c − am − dm.

(2) (x^{∗}, y^{∗}) is unstable if

(iv) \frac{−md^{2}+\sqrt{Δ}}{2(c−am−dm)} < f < \frac{cd}{c−am} and 0 < c − am − dm.

(3) (x^{∗}, y^{∗}) is nonhyperbolic if

(v) f = \frac{−md^{2}+\sqrt{Δ}}{2(c−am−dm)} and 0 < c − am − dm.

Remark 2.1: According to the results in [16]: the limit cycle of system (1.1) is unique and stable if it exists, we obtain that the equilibrium (x^{∗}, y^{∗}) in case (3) is stable if a nontrivial periodic orbit can be bifurcated by (x^{∗}, y^{∗}).

By applying Theorem in Kuang [14], we see that if (x^{∗}, y^{∗}) is unstable for system (1.1), then it will remain so for all τ > 0. Hence, from Lemma 2.1 we have

Theorem 2.2. When \frac{−md^{2}+\sqrt{Δ}}{2(c−am−dm)} < f < \frac{cd}{c−am} and 0 < c − am − dm with Δ = m^{2}d^{4} + 4cd^{2}(c − am − dm), the interior equilibrium (x^{∗}, y^{∗}) of system (1.2) is absolutely unstable.

Now we turn to investigate the roots of the characteristic equation (2.2) in λ. Obviously these roots dependent on τ. In the following we explore when Eq.(2.2) in λ has a pair of purely imaginary roots. Assume that for some τ ≥ 0, iω with ω > 0 is a root of Eq.(2.2) in λ. Then we have

−ω^{2} − α_{1}ωi − β_{2}ωie^{−ωτi} + (α_{1}β_{2} − α_{2}β_{1})e^{−ωτi} = 0.

STABILITY AND BIFURCATION IN A PREDATOR-PREY SYSTEM 7

Separating the real and imaginary parts, we obtain

−ω^{2} − β_{2}ω sin ωτ + (α_{1}β_{2} − α_{2}β_{1}) cos ωτ = 0,
−α_{1}ω − β_{2}ω cos ωτ − (α_{1}β_{2} − α_{2}β_{1}) sin ωτ = 0.

(2.3)

Eqs.(2.3) is equivalent to

\begin{aligned}cos ωτ & = −\frac{α_{2}β_{1}ω^{2}}{β_{2}^{2}ω^{2} + (α_{1}β_{2} − α_{2}β_{1})^{2}} , \\ sin ωτ & = −\frac{β_{2}ω^{3} + α_{1}ω(α_{1}β_{2} − α_{2}β_{1})}{β_{2}^{2}ω^{2} + (α_{1}β_{2} − α_{2}β_{1})^{2}}. \end{aligned} ; (2.4)

By Eqs.(2.4) we have

(2.5) (−\frac{α_{2}β_{1}ω^{2}}{β_{2}^{2}ω^{2} + (α_{1}β_{2} − α_{2}β_{1})^{2}} )^{2} + (−\frac{β_{2}ω^{3} + α_{1}ω(α_{1}β_{2} − α_{2}β_{1})}{β_{2}^{2}ω^{2} + (α_{1}β_{2} − α_{2}β_{1})^{2}})^{2} = 1.

Thus if Eq.(2.5) in ω has positive real solutions, then we can get the corresponding τ ≥ 0 from Eqs.(2.4) such that the characteristic equation (2.2) in λ has purely imaginary roots. It is not difficult to check Eq.(2.5) does have positive real roots by letting s = ω^{2} to change Eq.(2.5) into a polynomial equation in s. Furthermore, in order to determine the exact value of ω, by (2.3) we have

ω^{4} + (α_{1}^{2} − β_{2}^{2})ω^{2} − (α_{1}β_{2} − α_{2}β_{1})^{2} = 0,

which yields that ω takes the unique positive real number

(2.6)  ω_{0} = \sqrt{\frac{−(α_{1}^{2} − β_{2}^{2}) + \sqrt{(α_{1}^{2} − β_{2}^{2})^{2} + 4(α_{1}β_{2} − α_{2}β_{1})^{2}}}{2}} .

Therefore we have the following lemma. Lemma 2.3. Let Δ = m^{2}d^{4} + 4cd^{2}(c − am − dm).

(1) suppose that any one of the following conditions holds (i) d < f and c ≤ ma; (ii) d < f < \frac{cd}{c−am} and am < c ≤ am + dm; (iii) d < f < \frac{−md^{2}+\sqrt{Δ}}{2(c−am−dm)} and 0 < c − am − dm.

Then, only when τ = τ_{k} = ω_{0}^{−1} arcsin (−\frac{β_{2}ω_{0}^{3}+α_{1}ω_{0}(α_{1}β_{2}−α_{2}β_{1})}{β_{2}^{2}ω_{0}^{2}+(α_{1}β_{2}−α_{2}β_{1})^{2}})+2kπ for some k = 0, 1, 2 ⋯ (or some k = 1, 2, ⋯) if τ_{0} > 0 (τ_{0} < 0, respectively), the characteristic equation (2.2) in λ has a unique pair of purely imaginary roots ±iω_{0} where ω_{0} is given by (2.6).

8 STABILITY AND BIFURCATION IN A PREDATOR-PREY SYSTEM

(2) If f = \frac{−md^{2}+\sqrt{Δ}}{2(c−am−dm)} and 0 < c − am − dm, then only when τ = τ_{k}^{∗} = \frac{2kπ}{ω_{0}} for some k = 0, 1, 2, ⋯, the characteristic equation (2.2) in λ has a unique pair of purely imaginary roots ±iω_{0} where ω_{0} = \sqrt{α_{1}β_{2} − α_{2}β_{1}}. Proof (1) From above analysis, we can see that the characteristic equation (2.2) in λ has a unique pair of purely imaginary roots ±iω_{0} only when the nonnegative real number τ satisfies Eqs.(2.4). Noting that −α_{2}β_{1} > 0, from Eqs.(2.4) we get

τ = τ_{k} = ω_{0}^{−1} arcsin (−\frac{β_{2}ω_{0}^{3} + α_{1}ω_{0}(α_{1}β_{2} − α_{2}β_{1})}{β_{2}^{2}ω_{0}^{2} + (α_{1}β_{2} − α_{2}β_{1})^{2}}) + 2kπ.

Here we would like to point out that τ_{0} ≠ 0 when the one of the conditions (i)-(iii) holds. Hence, to guarantee τ_{k} > 0 we take k = 0, 1, 2 ⋯ (or k = 1, 2, ⋯) if τ_{0} > 0 ( τ_{0} < 0, respectively).

(2) Note that f = \frac{−md^{2}+\sqrt{Δ}}{2(c−am−dm)} implies α_{1} +β_{2} = 0. So we have ω_{0} = \sqrt{α_{1}β_{2} − α_{2}β_{1}} by (2.6). Correspondingly from Eqs.(2.4) we get τ = τ_{k}^{∗} = \frac{2kπ}{ω_{0}} , k = 0, 1, 2, · · ·. QED

Now we investigate how the real part of the roots of Eq.(2.2) in λ varies as τ varies in the small neighborhood of τ_{k}. Let u + iω be a root of Eq.(2.2) in λ. Separating the real and imaginary parts of Eq.(2.2), we then have

\begin{aligned}H_{1}(u, ω, τ) & = 0, \\ H_{2}(u, ω, τ) & = 0, \end{aligned} ; (2.7)

where H_{1}(u, ω, τ) = u^{2} − ω^{2} − α_{1}u − e^{−uτ} β_{2}u cos ωτ − e^{−uτ} β_{2}ω sin ωτ + (α_{1}β_{2} − α_{2}β_{1})e^{−uτ} cos ωτ and H_{2}(u, ω, τ) = 2uω −α_{1}ω −e^{−uτ} β_{2}ω cos ωτ +e^{−uτ} β_{2}u sin ωτ − (α_{1}β_{2} −α_{2}β_{1})e^{−uτ} sin ωτ. By Lemma 2.3 we have H_{1}(0, ω_{0}, τ_{k}) = H_{2}(0, ω_{0}, τ_{k}) = 0.

It is easy to check the Jacobian matrix J = ( \begin{matrix}\frac{∂H_{1}}{∂u} & \frac{∂H_{1}}{∂ω} \\ \frac{∂H_{2}}{∂u} & \frac{∂H_{2}}{∂ω} \end{matrix} ) satisfies |J|_{(0,ω_{0},τ_{k})} >

0. By means of implicit function theorem we deduce that Eqs.(2.7) defines u, ω as functions of τ in a neighborhood of (0, ω_{0}, τ_{k}) such that u(τ_{k}) = 0 and ω(τ_{k}) = ω_{0}. Moreover, \frac{du(τ_{k})}{dt} > 0. Thus this root of Eq.(2.2) crosses the imaginary axis from the left to the right as τ varies from the left to the right of τ_{k} by Rouché’s theorem and Lemma 2.1. The following theorem is an addition to Theorem 2.2, which can be proven by the above arguments and Lemma 2.3. Theorem 2.4. Let Δ = m^{2}d^{4} + 4cd^{2}(c − am − dm).

STABILITY AND BIFURCATION IN A PREDATOR-PREY SYSTEM 9

(1) When any one of the following conditions holds (i) d < f and c ≤ ma; (ii) d < f < \frac{cd}{c−am} and am < c ≤ am + dm; (iii) d < f < \frac{−md^{2}+\sqrt{Δ}}{2(c−am−dm)} and 0 < c − am − dm,

the interior equilibrium (x^{∗}, y^{∗}) of system (1.2) is conditionally stable. Thus small variation of delay doesn’t change the stability of the interior equilibrium.

(2) When f = \frac{−md^{2}+\sqrt{Δ}}{2(c−am−dm)} and 0 < c − am − dm, the interior equilibrium (x^{∗}, y^{∗}) of system (1.2) is unstable for all τ > 0.

Remark 2.2: From Remark 2.1 we can see that any small delay can change the stability of the interior equilibrium in the case (2) of Theorem 2.4. Remark 2.3: Theorem 2.2 and Theorem 2.4 imply the non-occurrence of absolute stability of the interior equilibrium (x^{∗}, y^{∗}) of system (1.2).

3. The Hopf Bifurcation

In this section, we will discuss the Hopf bifurcation of system (1.2) by choosing delay τ as a bifurcation parameter. According to the analysis in section 2, one can see that when the interior equilibrium (x^{∗}, y^{∗}) of system (1.1) is locally asymptotically stable, the characteristic equation (2.2) of system (1.2) has pure imaginary roots ±iω_{0} if and only if τ = τ_{k}, and the transversal condition holds at τ = τ_{k}. Thus the following theorem is implied by Theorem 2.4 (see [14]). Theorem 3.1. If any one of the following conditions holds

(i) d < f and c ≤ ma; (ii) d < f < \frac{cd}{c−am} and am < c ≤ am + dm; (iii) d < f < \frac{−md^{2}+\sqrt{Δ}}{2(c−am−dm)} and 0 < c − am − dm, here Δ = m^{2}d^{4} + 4cd^{2}(c −

am − dm), then a Hopf bifurcation occurs on a 2-dimensional center manifold for system (1.2) at the interior equilibrium (x^{∗}, y^{∗}) for τ = τ_{k} (here k ∈ N_{0} = {0, 1, 2, ⋯ } or k ∈ N = {1, 2, ⋯ } as stated in Lemma 2.3) which is associated with the simple eigenvalues ±iω_{0}. This center manifold is locally stable for τ = τ_{0} (or τ = τ_{1}, respectively) and unstable for τ = τ_{k}, k = 1, 2, ⋯ (or k = 2, 3, ⋯, respectively).

10 STABILITY AND BIFURCATION IN A PREDATOR-PREY SYSTEM

To determine the direction of the bifurcation and stability of the periodic solution created by Hopf bifurcation, the normal form theory for FDEs developed by Faria and Magalhães [10] will be employed for discussing the explicit expressions of the normal form of system (1.2) in terms of the original parameters in the small neighborhood of τ_{k}. By time scaling t → \frac{t}{τ}, system (1.2) is transformed into

\begin{aligned}ẋ(t) & = τ (x(t)(a − bx(t)) − \frac{cx(t)y(t)}{my(t) + x(t)}) =^{Δ} τf^{(1)}(x(t), y(t)), \\ ẏ(t) & = τy(t) (−d + \frac{fx(t − 1)}{my(t − 1) + x(t − 1)}) =^{Δ} τf^{(2)}(y(t), x(t − 1), y(t − 1)) \end{aligned} ; (3.1)

in a fixed phase space C_{1} = C([−1, 0]; R_{+}^{2}).

Even though the normal form procedure for FDEs is given in [10], it is very difficult to completely treat for all permissible parameters of system (3.1) since it involves hard computation.

To simplify computation, we only consider the Hopf bifurcation of system (3.1) in the case that f = \sqrt{\frac{cd^{2}}{c−am}}. In this case, we have α_{1} = 0 in system (2.1).

Theorem 3.2. Suppose am < c < am + dm and f = \sqrt{\frac{cd^{2}}{c−am}}. Then system (3.1) has a unique interior equilibrium (x_{0}, y_{0}) where x_{0} = \frac{a}{b} − \frac{c−\sqrt{c^{2}−amc}}{bm}, y_{0} = \frac{\sqrt{c}−\sqrt{c−am}}{m\sqrt{c−am}}x_{0}. Furthermore there exists a τ^{(0)} such that (x_{0}, y_{0}) is locally asymptotically stable (unstable) whenever 0 < τ < τ^{(0)} (τ > τ^{(0)}, respectively). When τ = τ^{(k)}, there exists a unique pair of simple purely imaginary roots ±iω_{k} for the characteristic equation of the linearized equation of system (3.1) around (x_{0}, y_{0}), where τ^{(k)} = \sqrt{\frac{ω_{k}^{2}(\sqrt{b_{2}^{4}+4a_{2}^{2}b_{1}^{2}}−b_{2}^{2})}{2a_{2}^{2}b_{1}^{2}}}, ω_{k} = arccos ^{\frac{b_{2}^{2}−\sqrt{b_{2}^{4}+4a_{2}^{2}b_{1}^{2}}}{2a_{2}b_{1}}} + 2kπ, k = 0, 1, 2, ⋯, a_{2} = −\frac{c(x_{0})^{2}}{(my_{0}+x_{0})^{2}} , b_{1} = \frac{m(y_{0})^{2}}{(my_{0}+x_{0})^{2}} \sqrt{\frac{cd^{2}}{c−am}} and b_{2} = −\frac{mx_{0}y_{0}}{(my_{0}+x_{0})^{2}} \sqrt{\frac{cd^{2}}{c−am}}.

Proof First of all, we translate the equilibrium (x_{0}, y_{0}) of system (3.1) to the origin by letting x_{1} = x − x_{0}, x_{2} = y − y_{0}. By noting that f = \sqrt{\frac{cd^{2}}{c−am}} and separating the linear terms from the nonlinear terms, system (3.1) can be changed into

\begin{aligned}ẋ_{1}(t) & = τ(a_{2}x_{2}(t) + F_{1}(x_{1}(t), x_{2}(t))), \\ ẋ_{2}(t) & = τ(b_{1}x_{1}(t − 1) + b_{2}x_{2}(t − 1) + F_{2}(x_{2}(t), x_{1}(t − 1), x_{2}(t − 1))), \end{aligned} ; (3.2)

STABILITY AND BIFURCATION IN A PREDATOR-PREY SYSTEM 11

here a_{2} = −\frac{c(x_{0})^{2}}{(my_{0}+x_{0})^{2}} , b_{1} = \frac{m(y_{0})^{2}}{(my_{0}+x_{0})^{2}} \sqrt{\frac{cd^{2}}{c−am}}, b_{2} = −\frac{mx_{0}y_{0}}{(my_{0}+x_{0})^{2}} \sqrt{\frac{cd^{2}}{c−am}} and

F_{1}(x_{1}(t), x_{2}(t)) = f^{(1)}(x_{1}(t) + x_{0}, x_{2}(t) + y_{0}) − \frac{∂f^{(1)}(x_{0},y_{0})}{∂x}x_{1}(t) − \frac{∂f^{(1)}(x_{0},y_{0})}{∂y}x_{2}(t),

F_{2}(x_{2}(t), x_{1}(t − 1), x_{2}(t − 1)) = f^{(2)}(x_{2}(t) + y_{0}, x_{1}(t − 1) + x_{0}, x_{2}(t − 1) + y_{0})

−\frac{∂f^{(2)}(y_{0},x_{0},y_{0})}{∂y}x_{2}(t) − \frac{∂f^{(2)}(y_{0},x_{0},y_{0})}{∂x(t−1)}x_{1}(t − 1) − \frac{∂f^{(2)}(y_{0},x_{0},y_{0})}{∂y(t−1)}x_{2}(t − 1).

To study the local stability of the zero solution of system (3.2), we consider the characteristic equation of the linearized equation

\begin{aligned}ẋ_{1}(t) & = τa_{2}x_{2}(t), \\ ẋ_{2}(t) & = τ(b_{1}x_{1}(t − 1) + b_{2}x_{2}(t − 1)), \end{aligned} ; (3.3)

around the equilibrium (0, 0), which is given by

(3.4)  Δ(λ, τ) = λ^{2} − b_{2}τλe^{−λ} − a_{2}b_{1}τ^{2}e^{−λ} = 0,

Assume that for some τ > 0 equation (3.4) in λ has a root iω, ω > 0. Then

−ω^{2} − b_{2}τω sin ω − a_{2}b_{1}τ^{2} cos ω = 0,

−b_{2}τω cos ω + a_{2}b_{1}τ^{2} sin ω = 0.

(3.5)

It is easy to see that Eqs.(3.5) is equivalent to

\begin{matrix}−ω^{2} cos ω − a_{2}b_{1}τ^{2} = 0, \\ −b_{2}τ − ω sin ω = 0. \end{matrix} ; (3.6)

Solving Eqs.(3.6) yields the following solutions

ω_{k} = arccos \frac{b_{2}^{2} − \sqrt{b_{2}^{4} + 4a_{2}^{2}b_{1}^{2}}}{2a_{2}b_{1}} + 2kπ, τ^{(k)} = \sqrt{\frac{ω_{k}^{2}(\sqrt{b_{2}^{4} + 4a_{2}^{2}b_{1}^{2}} − b_{2}^{2})}{2a_{2}^{2}b_{1}^{2}}}, k = 0, 1, · · · .

As a result, Eq.(3.4) in λ has the root iω_{k} if and only if τ = τ^{(k)}. Now let H_{1}(u, ω, τ) = ReΔ(u+iω, τ), H_{2}(u, ω, τ) = ImΔ(u+iω, τ) and H(u, ω, τ) = (H_{1}(u, ω, τ), H_{2}(u, ω, τ)). By the implicit function theorem, in a small neighborhood of τ^{(k)} the equation H(u, ω, τ) = 0 defines u, ω as functions of τ such that u(τ^{(k)}) = 0, ω(τ^{(k)}) = ω_{k} and u′(τ^{(k)}) > 0. Thus the root ω_{k}i of equation Δ(λ, τ^{(k)}) = 0 is simple. Using the similar method, we can prove that the root −iω_{k} of equation Δ(λ, τ^{(k)}) = 0 is simple too. This completes the proof of the theorem. QED

We now introduce some notations. Considering the phase space C_{1}, system (3.2) can be written as an equation in C_{1}

(3.7)  Ẋ(t) = τL(X_{t}) + τF (X_{t}),

where X(t) is a vector (x_{1}(t), x_{2}(t)), X_{t} ∈ C_{1}, X_{t}(θ) = X(t + θ), −1 ≤ θ ≤ 0 and L : C_{1} → R_{+}^{2}, F : C_{1} → R_{+}^{2} are given by

L(φ) = ( \begin{matrix}a_{2}φ_{2}(0) \\ b_{1}φ_{1}(−1) + b_{2}φ_{2}(−1) \end{matrix} ) , F (φ) = ( \begin{matrix}F_{1}(φ_{1}(0), φ_{2}(0)) \\ F_{2}(φ_{2}(0), φ_{1}(−1), φ_{2}(−1)) \end{matrix} ) ,

for φ = (φ_{1}, φ_{2}) ∈ C_{1}.

Let A be the generator of the linear semigroup corresponding to (3.3). When τ = τ^{(k)}, A has a pair of simple purely imaginary characteristic roots ±iω_{k} and no other characteristic roots with zero real part. Hence, Hopf bifurcation may occur at τ = τ^{(k)}.

For a fixed k ∈ N_{0}, define Λ = {−iω_{k}, iω_{k}} and introduce a new parameter ν = τ − τ^{(k)}. System (3.7) can be written as

(3.8)  Ẋ(t) = τ^{(k)}L(X_{t}) + F_{0}(X_{t}, ν),

where F_{0}(φ, ν) = νL(φ) + (τ^{(k)} + ν)F (φ).

To describe the Hopf bifurcation occurring at τ^{(k)} and the stability of the periodic solutions by Hopf bifurcation, we will apply the normal form theory for FDES in [10] to system (3.8). Let the phase space C_{1} be decomposed by Λ as C_{1} = P ⊕Q, where P is the generalized eigenspace associated with Λ. Consider the adjoint bilinear form (·, ·) on C_{1}^{∗} × C_{1} associated with the linear equation Ẋ(t) = τ_{k}L(X_{t}),

(ψ, ϕ) = ψ(0)ϕ(0) − ∫_{−r}^{0} ∫_{0}^{θ} ψ(ξ − θ)dη(θ)ϕ(ξ)dξ.

Assume that Φ = (ϕ_{1}, ϕ_{2}) is a basis for P and Ψ = col(ψ_{1}, ψ_{2}) is a basis for the dual space P^{∗} in C_{1}^{∗} associated with the eigenvalues ±iω_{k} of the adjoint equations. Then it can be normalized so that (Φ, Ψ) = I. Here it is convenient to combine one complex coordinate and two complex conjugate basis vectors to describe a twodimensional real subspace P. Consider system (3.8) in C([−1, 0]; ℂ), still denoted by C_{1}. Note that Φ̇ = ΦB, where B = diag(iω_{k}, −iω_{k}) is a diagonal matrix. Therefore Φ and Ψ are 2 × 2 matrices of the form

Φ(θ) = [ϕ_{1}(θ), ϕ_{2}(θ)], ϕ_{1}(θ) = e^{iω_{k}θ}v, ϕ_{2}(θ) = \overline{ϕ_{1}(θ)}, −1 ≤ θ ≤ 0,

Ψ(s) = ( \begin{matrix}ψ_{1}(s) \\ ψ_{2}(s) \end{matrix} ) , ψ_{1}(s) = e^{−iω_{k}s}u^{T} , ψ_{2}(s) = \overline{ψ_{1}(s)}, 0 ≤ s ≤ 1,

where the bar means complex conjugation, u^{T} is the transpose of u,

u = ( \begin{matrix}u_{1} \\ u_{2} \end{matrix} ) = ( u_{1}_{\frac{iω_{k}e^{iω_{k}} u_{1}}{b_{1}τ^{(k)}}} ) ∈ ℂ^{2}, v = ( \begin{matrix}v_{1} \\ v_{2} \end{matrix} ) = ( 1_{\frac{iω_{k}}{a_{2}τ^{(k)}}} ) ∈ ℂ^{2},

with u_{1} = \frac{a_{2}b_{1}τ^{(k)}}{(2a_{2}b_{1}τ^{(k)}−b_{2}ω_{k}^{2})+(a_{2}b_{1}ω_{k}τ^{(k)}+b_{2}ω_{k})i}.

Enlarge the phase space C_{1} by considering the space BC and using the decomposition X_{t} = Φz(t) + y_{t}, z ∈ ℂ^{2}, y_{t} ∈ Q′. We decompose system (3.8) as follows

\begin{aligned}ż & = Bz + Ψ(0)F_{0}(Φz + y, ν), \\ ẏ & = A_{Q′}y + (I − π)X_{0}F_{0}(Φz + y, ν). \end{aligned} ; (3.9)

Here and forward we refer to [10] for results and explanations of some notations involved. Following the procedure of reducing normal form in [10], we consider the Taylor formula

Ψ(0)F_{0}(Φz + y, ν) = \frac{1}{2}g_{2}^{(1)}(z, y, ν) + \frac{1}{3!}g_{3}^{(1)}(z, y, ν) + h.o.t.,

(I − π)X_{0}F_{0}(Φz + y, ν) = \frac{1}{2}g_{2}^{(2)}(z, y, ν) + \frac{1}{3!}g_{3}^{(2)}(z, y, ν) + h.o.t.,

where g_{j}^{(1)}(z, y, ν), g_{j}^{(2)}(z, y, ν)(j = 2, 3) are homogeneous polynomials in (z, y, ν) of degree j with coefficients in ℂ^{2} and Kerπ respectively, and h.o.t. stands for higher order terms. Thus in a finite dimensional locally invariant manifold tangent to the invariant subspace P of (3.9) at x = 0, ν = 0, the normal form of (3.9) is given by

(3.10)  ż = Bz + \frac{1}{2}ḡ_{2}^{(1)}(z, 0, ν) + \frac{1}{3!}ḡ_{3}^{(1)}(z, 0, ν) + h.o.t.,

where ḡ_{2}^{(1)}, ḡ_{3}^{(1)} are the second and third order terms in (z, ν) respectively. Using the notations in [10], we have

ḡ_{2}^{(1)}(z, 0, ν) = Proj_{Ker(M′_{2})}g_{2}^{(1)}(z, 0, ν),

where

Ker(M′_{2}) = span \{( \begin{matrix}z_{1}ν \\ 0 \end{matrix} ) , ( \begin{matrix}0 \\ z_{2}ν \end{matrix} )\} ,

g_{2}^{(1)}(z, 0, ν) = ( \begin{matrix}2i\frac{ω_{k}}{τ^{(k)}} u^{T} vz_{1}ν + a_{20}z_{1}^{2} + a_{11}z_{1}z_{2} + a_{02}z_{2}^{2} \\ −2i\frac{ω_{k}}{τ^{(k)}} ū^{T} v̄z_{2}ν + ā_{02}z_{1}^{2} + ā_{11}z_{1}z_{2} + ā_{20}z_{2}^{2} \end{matrix} ) .

We would like to point out that the coefficients a_{20}, a_{11} and a_{02} can be expressed in terms of the parameters in original equations (3.8) by the algorithm of normal form. However the explicit formulas for a_{20}, a_{11} and a_{02} are tediously long since both Φ and Ψ have quite complicated expressions as shown above. Therefore, for the sake of convenience we prefer to use the notations a_{20}, a_{11} and a_{02} instead of their exact expressions in the following context.

To eliminate these nonresonant terms in the quadratic terms g_{2}^{(1)}(z, 0, ν) we have to make a series of transformations of variables, which may change the coefficients of the cubic terms of g_{3}^{(1)}(z, 0, ν). From the canonical basis for Ker(M′_{2}), we have

\frac{1}{2}ḡ_{2}^{(1)}(z, 0, ν) = ( \begin{matrix}i\frac{ω_{k}}{τ^{(k)}} u^{T} vz_{1}ν \\ −i\frac{ω_{k}}{τ^{(k)}} ū^{T} v̄z_{2}ν \end{matrix} ) .

Furthermore, recalling the operator M′_{3} we know

Ker(M′_{3}) = span \{( \begin{matrix}z_{1}^{2}z_{2} \\ 0 \end{matrix} ) , ( \begin{matrix}z_{1}ν^{2} \\ 0 \end{matrix} ) , ( \begin{matrix}0 \\ z_{1}z_{2}^{2} \end{matrix} ) , ( \begin{matrix}0 \\ z_{2}ν^{2} \end{matrix} )\} .

However the terms O(|z|ν^{2}) are irrelevant to determining the generic Hopf bifurcation. Hence we only need to compute the coefficient of z_{1}^{2}z_{2}. After some computations we find that the coefficient of z_{1}^{2}z_{2} is

c = \frac{i}{2ω_{k}} (a_{20}a_{11} − 2|a_{11}|^{2} − \frac{1}{3}|a_{02}|^{2}) + \frac{1}{2}a_{21},

where \frac{1}{2}a_{21} is the coefficient of the term z_{1}^{2}z_{2} in g_{3}^{(1)}(z, 0, ν). Thus

\frac{1}{3!}ḡ_{3}^{(1)}(z, 0, ν) = ( \begin{matrix}cz_{1}^{2}z_{2} \\ c̄z_{1}z_{2}^{2} \end{matrix} ) + O(|z|ν^{2}).

The normal form (3.10) relative to P can be written in real coordinates (x, y), through the change of variables z_{1} = x − iy, z_{2} = x + iy. Followed by the use of polar coordinates (r, θ), x = r cos θ, y = r sin θ, this normal form becomes

\begin{aligned}ṙ & = c_{1}νr + c_{2}r^{3} + O(ν^{2}r + |(r, ν)|^{4}), \\ θ̇ & = −ω_{k} + O(|(r, ν)|), \end{aligned}

where c_{1} = Re(i\frac{ω_{k}}{τ^{(k)}} u^{T} v), c_{2} = Re(c). Therefore these arguments lead to the following theorem.

STABILITY AND BIFURCATION IN A PREDATOR-PREY SYSTEM 15

Theorem 3.3. If c_{2} ≠ 0 and τ^{(k)} > 0, then system (3.8) exhibits a generic Hopf bifurcation. The periodic orbits of system (3.8) bifurcating from the origin and ν = 0 satisfy

r(t, ν) = \sqrt{−\frac{c_{1}ν}{c_{2}}} + O(ν), θ(t, ν) = −2ω_{k}t + O(|ν|^{\frac{1}{2}} )

so that

(i) if c_{1}c_{2} < 0 (c_{1}c_{2} > 0 respectively), there exists a unique nontrivial periodic orbit in the neighborhood of r = 0 for ν > 0 (ν < 0 respectively) and no nontrivial periodic orbits for ν < 0 (ν > 0 respectively);

(ii) the nontrivial periodic solutions in the center manifold are stable if c_{2} < 0

and unstable if c_{2} > 0. As an example of the preceding results, we now consider

\begin{aligned}ẋ(t) & = τ (x(t)(1 − x(t)) − \frac{2x(t)y(t)}{y(t) + x(t)}) , \\ ẏ(t) & = τy(t) (−\sqrt{2} + \frac{2x(t − 1)}{y(t − 1) + x(t − 1)}) \end{aligned} ; (3.11)

in a fixed phase space C_{1} = C([−1, 0]; R_{+}^{2}).

System (3.11) has a unique interior equilibrium E^{∗} = (\sqrt{2} − 1, (\sqrt{2} − 1)^{2}). The purpose here is to consider the dynamics of (3.11) in the neighborhood of E^{∗} at a bifurcation point of the parameter τ. In the following, we shall apply the Theorem 3.2 and Theorem 3.3 to system (3.11) and obtain Theorem 3.4. Let ω_{1} ∈ (2π, 3π) and τ_{1} be defined by

ω_{1} = arccos \frac{\sqrt{5} − 1}{2} + 2π, τ_{1} = ω_{1} \frac{\sqrt{\sqrt{5} − 1}}{2 − \sqrt{2}}.

Then for τ = τ_{1}, ±iω_{1} are simple roots of the characteristic equation of system (3.11) at E^{∗} and the remaining roots have nonzero real parts. Furthermore, for system (3.11) there is a subcritical Hopf bifurcation at τ = τ_{1}, with the nontrivial periodic orbit being unstable in the center manifold, and the center manifold is unstable. Proof According to Theorem 3.2, we get immediately when τ = τ_{k} = ω_{k} ^{\frac{\sqrt{\sqrt{5}−1}}{2−\sqrt{2}}}, ±iω_{k} are simple roots of the characteristic equation of system (3.11) at E^{∗} and the remaining roots have nonzero real parts, here ω_{k} = arccos \frac{\sqrt{5}−1}{2} + 2kπ, k = 0, 1, 2, · · ·. More precisely, as τ ∈ (0, τ_{0}), all roots of the characteristic equation

16 STABILITY AND BIFURCATION IN A PREDATOR-PREY SYSTEM

of system (3.11) at E^{∗} have negative real parts, system (3.11) then has a locally asymptotically stable equilibrium E^{∗}; as τ = τ_{0}, ±iω_{0} are simple roots of the characteristic equation of system (3.11) at E^{∗} and the remaining roots have negative real parts, system (3.11) then has a stable center manifold; as τ = τ_{k}, ±iω_{k} are simple roots of the characteristic equation of system (3.11) and the remaining roots have nonzero real parts, system (3.11) then has a unstable center manifold for k = 1, 2, · · ·.

We are interested in studying the behavior of periodic solutions bifurcating from τ = τ_{1}. Taking ν = τ − τ_{1}, the system (3.11), in the notations used in the inference of Theorem 3.3, becomes

(3.12)  Ẋ(t) = τ_{1}L(X_{t}) + F_{0}(X_{t}, ν),

where F_{0}(φ, ν) = νL(φ) + (τ_{1} + ν)F (φ) for φ = (φ_{1}, φ_{2}) ∈ C_{1},

L(φ) = ( \begin{matrix}−φ_{2}(0) \\ (\sqrt{2} − 1)^{2}φ_{1}(−1) + (−\sqrt{2} + 1)φ_{2}(−1) \end{matrix} ) ,

F (φ) = ( \begin{matrix}F_{1}(φ_{1}(0), φ_{2}(0)) \\ F_{2}(φ_{2}(0), φ_{1}(−1), φ_{2}(−1)) \end{matrix} ) .

We can compute the vectors u and v according to their expressions given in the inference of Theorem 3.3. Omitting their complicated expressions, we obtain directly the numerical results by means of the software Mathematica

u = ( \begin{matrix}u_{1} \\ u_{2} \end{matrix} ) = ( \begin{matrix}−0.05827 − 0.06901i \\ 0.27164 + 0.05601i \end{matrix} ) , v = ( \begin{matrix}v_{1} \\ v_{2} \end{matrix} ) = ( \begin{matrix}1 \\ −0.52689i \end{matrix} ) .

Following the procedure of reducing normal form, we can get the coefficients as ν = 0, given in the inference of Theorem 3.3, by means of the software Mathematica

a_{20} = 0.765562 − 0.380496i, a_{11} = 0.345173 − 0.291445i,
a_{02} = 0.184747 − 1.41851i, a_{21} = 0.28828 + 0.85934i.

By the expressions of c_{1} and c_{2} in Theorem 3.3, we have

c_{1} = 0.111771, c_{2} = 0.168797.

STABILITY AND BIFURCATION IN A PREDATOR-PREY SYSTEM 17

Therefore, Theorem 3.3 implies the occurrence of a subcritical Hopf bifurcation at the point τ = τ_{1} from equilibrium E^{∗} with the associated periodic solution being unstable in the center manifold.

4. Discussion

It is interesting to compare our local stability results with the global results in [6]. Though we don’t make global qualitative analysis of the model, we make the local stable analysis of the interior equilibrium depending on all parameters and bifurcation analysis by choosing time delay τ as parameter parameter. Our analysis supports their results (i.e. the global stability of the interior equilibrium must require conditions on τ) and also complements their results in the sense that we give many conditions depending on all parameters in system (1.2) such that the interior equilibrium is not global stable, and the existence and stability of the nontrivial periodic solution of system (1.2). Our analysis indicates that the dynamics of the ratio-dependent predator-prey systems with time delay can be much more complicated than we may have expected. It is still interesting and challenging to systematically describe global dynamics of the model for all parameters by means of the local properties of the interior equilibrium.

Acknowledgments. The authors are grateful to the anonymous referee for his/her helpful comments and valuable suggestions.

References

[1] H. R. Akcakaya, Population cycles of mammals: evidence for a ratio-dependent predation hypothesis, Ecol. Monogr. 62(1992), 119-142.

[2] R. Arditi and L. R. Ginzburg, Coupling in predator-prey dynamics: ratio-dependence, J. Theoret. Biol. 139 (1989),311–326.

[3] R. Arditi, L. R. Ginzburg and H. R. Akcakaya, Variation in plankton densities among lakes: a case of ratio-dependent models, American Naturalist 138 (1991), 1287-1296.

[4] F. G. Boese, Stability with respect to the delay: on a paper of K.L.Cooke and P. Van den Driessche, J. Math. Analysis and Appl. 228(1998), 293-321.

[5] A. A. Berryman, The origins and evolution of predator-prey theory, Ecology 75(1992), 1530-1535.

[6] E. Beretta and Y. Kuang, Global analyses in some delayed ratio-dependent predator-prey systems, Nonlinear Analysis, Theory, Methods & Applications 32 (1998), 381–408.

[7] F. Brauer, Absolute stability in delay equations, J. Differential Equations 69(1987), 185-191.

[8] K. L. Cooke and Z. Grossman, Discrete delay, distributed delay and stability switches, J. Math.Anal.Appl. 86(1982), 592-627.

[9] J. Dieudonné, Foundations of Modern Analysis, Academic Press, New York,1960.

[10] T.Faria and L.T.Magalhães , Normal forms for retarded functional differential equations with parameters and applications to Hopf bifurcation, J. Differential Equations 122 (1995), 181–200.

[11] T.Faria , Stability and bifurcation for a delayed predator-prey model and the effect of diffusion, preprint.

[12] A. P. Gutierrez, The physiological basis of ratio-dependent predator-prey theory: a metabolic pool model of Nicholson’s blowflies as an example, Ecology 73 (1992), 1552-1563.

[13] C. Jost, O. Arino and R. Arditi, About deterministic extinction in ratio-dependent predatorprey models, Bull. Math. Biol. 61 (1999), 19–32.

[14] Y. Kuang, Delay Differential Equations with Applications in Population Dynamics. Academic Press, New York, 1993.

[15] Y. Kuang, Rich dynamics of Gause-type ratio-dependent predator-prey system, Fields Institute Communications 21 (1999), 325–337.

[16] Y. Kuang and E. Beretta, Global qualitative analysis of a ratio-dependent predator-prey system, J. Math. Biol. 36 (1998), 389–406.

[17] P. Lundberg and J. M. Fryxell, Expected population density versus productivity in ratiodependent and prey-dependent models, American Naturalist 146 (1995), 153-161.

[18] M. L. Rosenzweig, Paradox of enrichment: destabilization of exploitation ecosystems in ecological time, Science 171 (1971), 385-387.

[19] D. Xiao and S. Ruan, Global dynamics of a ratio-dependent predator-prey system, J. Math. Biol. 43 (2001), 268-290.