# On the Dynamics of Two-consumers-one-resource Competing Systems with Beddington-DeAngelis Functional Response

(1)

## Functional Response

∗1

†2

### , and Ting-Hui Yang

‡3

1Department of Mathematics and The National Center for Theoretical Science, National Tsing-Hua University, Hsinchu 30013, Taiwan, Republic of China 2Department of Mathematics, University of Miami, Coral Gables, FL 33124-4250, USA

3Department of Mathematics, Tamkang University, 151 Ying-chuan Road, New Taipei City 25137, Taiwan, Republic of China

### October 1, 2012

Abstract

In this paper we study a two-consumers-one-resource competing system with Beddington-DeAngelis functional response. The two consumers com-peting for a renewable resource have intraspecific competition among their own populations. We investigate the extinction and uniform persistence of predators, Hopf bifurcation of positive equilibrium. We compare the system to the system without interference effects. Analytically we study the com-petition of two identically species with different interference effects. Also we study the relaxation oscillation in the case of interference effects. In the dis-cussion section, we present extensive numerical simulations to understand the interference effects on the competition outcomes.

Research was partially supported by National Council of Science, Republic of China.Research was partially supported by National Science Foundation (DMS-1022728).Research was partially supported by National Council of Science, Republic of China.

(2)

### Introduction

In this paper we shall study a two-consumers-one-resource system with Beddington-DeAngelis functional response. The two consumers (predators) competing for a renewable resource (prey) have interference competition among their own popu-lation. The mathematical model takes the following system of three nonlinear ordinary differential equations [2, 7, 13]:

dx dt = rx(1 − x K) − m1x a1+ x + b1y1 y1− m2x a2+ x + b2y2 y2 dy1 dt = ( e1m1x a1+ x + b1y1 − d1)y1 (1.1) dy2 dt = ( e2m2x a2+ x + b2y2 − d2)y2

with initial values x(0) = x0 > 0, y1(0) = y10 > 0, y2(0) = y20 > 0.

In (1.1) x(t), y1(t), and y2(t) represent the population density of prey and

two predators respectively at time t. In the absence of predation, the prey grows logistically with intrinsic growth rate r and carrying capacity K. The i-th preda-tor consumes the prey according to the Beddington-DeAngelis functional response

mixyi

ai+x+biyi and its growth rate is

eimixyi

ai+x+biyi where ei is the conversion efficiency

coef-ficient ; mi is the maximal consumption rate; ai is the half-satuation constant and

bi measures the intraspecific interference among the population of i-th predator;

di is the death rate.

Note that if b1 = b2 = 0 the system (1.1) is reduced to a system with Holling

type II functional responses: dx dt = rx(1 − x K) − m1x a1+ x y1− m2x a2+ x y2 dy1 dt = ( e1m1x a1+ x − d1)y1 (1.2) dy2 dt = ( e2m2x a2+ x − d2)y2

Hsu, Hubbel and Waltman [12, 11], Butler and Waltman [5]. Keener [16], Muratri and Rinaldi [18], Smith [19], Liu, Xiao and Yi [17], among others, have showed that the system (1.2) exhibits coexistence in the sense of Armstrong and McGehee [1], that is, for appropriate parameter values and suitable initial population densities

(3)

(x(0), y(0), z(0)), the model does predicts coexistence of the two predators via a locally attracting periodic orbit.

This paper is organized as follows. In Section 2, we study existence and stability of equilibria in system (1.1), including the semi-trivial equilibria( i.e., with survival of only one predator species ) and the positive equilibrium ( with the coexistence of both competing predators). Sufficient conditions for the uniform persistence of the system are obtained. In Section 3, we construct a Lyapunov function to establish the global stability of the positive equilibrium. In Section 4, we consider the competition of two identical predators with different interference effects. In Section

5, we study relaxation oscillations to system (1.1) with r ≫ 1 and bi = O(ε1+µi)

where ε = 1/r and µi > 0, i = 1, 2. Numerical simulations are presented to explain

the obtained results.

### Subsystems

Consider the following predator-prey system which is a subsystem of (1.1):              dx dt = rx(1 − x K) − mx a + x + byy dy dt = ( emx a + x + by − d)y x(0) > 0, y(0) > 0. (2.1)

With the scaling:

t → rt, x → x/K, y → by/K (2.2)

the system (2.1) becomes        dx dt = x(1 − x) − sxy x + y + A dy dt = δy(−D + x x + y + A) (2.3) where s = m br, δ = me r , D = d me, A = a K.

(4)

From the analysis in [6, 14, 15], we have the following results about the asymp-totic behavior of the solutions of (2.3). First results is about the extinction of predator.

Proposition 2.1. If em ≤ d or K ≤ λ = a

(me d )−1

, then the equilibrium (1, 0) of

system (2.3) is globally asymptotically stable, or equivalently the equilibrium (K, 0)

of system (2.1) is globally asymptotically stable.

Now we assume (H1) K > λ > 0.

Under the assumption (H1), there exists three equilibria (0, 0), (1, 0) and (x∗, y∗)

where x∗ and y∗ are positive and satisfy

     1 − x∗− sy∗ x∗+ y∗+ A = 0 x∗ x∗ + y∗+ A = D. (2.4) Obviously, we have s > sy∗ x∗+ y∗+ A = 1 − x∗

and from (2.4) it follows that      y∗ = (1 − x∗)(x∗+ A) x∗+ s − 1 , x2∗+ (s − 1 − Ds)x∗− DAs = 0. (2.5)

From the second equation of (2.5), we have

x∗+ s − 1 > x∗ + s − 1 − Ds =

DAs x∗

, y∗ > 0 . The variational matrix of system (2.3) is given by

J(x, y) = "

1 − 2x − x+y+Asy + (x+y+A)sxy 2

−sx x+y+A + sxy (x+y+A)2 δy(y+A) (x+y+A)2 δx x+y+A− δxy (x+y+A)2 − Dδ # . (2.6)

(5)

Proposition 2.2. Let the assumption (H1) hold.

(i) If tr J(x∗, y∗) < 0 then the equilibrium (x∗, y∗) of system (2.3) is globally

asymptotically stable.

(ii) If tr J(x∗, y∗) > 0 then there exists a unique limit cycle for system (2.3).

Furthermore, (1) If s ≤ max{δ, Dδ 1+D+ 1 1−D2} or equivalently b ≥ min{1 e, m2e2− d2 de(me − d) + mre2} (2.7) then tr J(x∗, y∗) ≤ 0. (2) If s > max{δ, Dδ 1+D + 1 1−D2} or equivalently 0 ≤ b < min{1 e, m2e2 − d2 de(me − d) + mre2} (2.8)

then there exists 0 < A∗ < 1−DD such that tr J(x∗, y∗) < 0 (> 0) if and only

of A > A∗ (A < A∗).

Remark 2.1. In above (ii), if we set K∗ = a/A∗, then the prey and predator coexist

in equilibrium if the carrying capacity K satisfies λ < K < K∗ and the prey and

predator populations exhibit periodic oscillation if K > K∗.

Let ¯x = Kx∗, ¯y =

K

by∗. From (2.3), (¯x, ¯y) is the positive equilibrium of the

system (2.1). We summarize the results for the system (2.1) in the Table I.

### Equilibria Analysis and Uniform Persistence

In this section, we shall find all equilibria of system (1.1) and do their stability analysis. Consider dx dt = rx(1 − x K) − m1x a1+ x + b1y1 y1− m2x a2+ x + b2y2 y2 := f (x, y1, y2), dy1 dt = ( e1m1x a1+ x + b1y1 − d1)y1 := g(x, y1), dy2 dt = ( e2m2x a2+ x + b2y2 − d2)y2 := h(x, y2).

(6)

Table I: Stability of equilibria for system (2.1)

Condition

em ≤ d or K ≤ λ (K, 0) is globally asymptotically stable

K > λ

(¯x, ¯y) is globally asymptotically stable

and

b ≥ min{1e, m2e2d2

de(me−d)+mre2}

λ < K < K∗

(¯x, ¯y) is globally asymptotically stable

and b < min{1e, m2e2d2

de(me−d)+mre2}

K > K∗ > λ

(¯x, ¯y) is an unstable focus and there exists a unique limit cycle

and b < min{1e, m2e2−d2

de(me−d)+mre2}

Then the Jacobian matrix of the system (1.1) takes the form

J(x, y1, y2) =    fx fy1 fy2 gx gy1 0 hx 0 hy2    (2.9)

(7)

where fx = r(1 − x K) − m1y1 a1+ x + b1y1 − m2y2 a2+ x + b2y2 + x(−r K + m1y1 (a1+ x + b1y1)2 + m2y2 (a2+ x + b2y2)2 ) fy1 = − m1x(a1+ x) (a1 + x + b1y1)2 fy2 = − m2x(a2+ x) (a2 + x + b2y2)2 gx = e1m1y1(a1+ b1y1) (a1+ x + b1y1)2 gy1 = e1m1x a1+ x + b1y1 − d1− b1e1m1xy1 (a1+ x + b1y1)2 = e1m1x(a1+ x) (a1+ x + b1y1)2 − d1 hx = e2m2y2(a2+ b2y2) (a2+ x + b2y2)2 hy2 = e2m2x a2+ x + b2y2 − d2− b2e2m2xy2 (a2+ x + b2y2)2 = e2m2x(a2+ x) (a2+ x + b2y2)2 − d2. We now consider the equilibria and periodic solutions on the boundary.

(a) E0 = (0, 0, 0). The trivial equilibrium E0always exists and is a saddle with

two-dimensional stable manifold {(x, y, z) : x = 0, y1 > 0, y2 > 0} and one-dimensional

unstable manifold {(x, y, z) : y1 = 0, y2 = 0}.

(b) EK = (K, 0, 0). The semi-trivial equilibrium EK always exists. The Jacobian

matrix at EK is given by J(EK) =    −r ∗ ∗ 0 e1m1K a1+K − d1 0 0 0 e2m2K a2+K − d2   .

Then EK is asymptotically stable if

e1m1K a1+ K − d1 < 0 and e2m2K a2+ K − d2 < 0 . We note that eimiK ai+K − di < 0 if and only if eimi ≤ di or K < λi = ai (eimi di ) − 1 ,

(8)

where λi is the break-even density for the i-th predator where there is no

intraspe-cific competition within the population of the i-th predator. If K > λ1and K > λ2

then EK is a saddle with one-dimensional stable manifold {(x, y1, y2) : x > 0, y1=

y2 = 0}.

Actually, we can verify the global asymptotical stability of EK under a weaker

condition in the following lemma.

Lemma 2.3. If eimi ≤ di then lim supt→∞yi(t) = 0 for i = 1, 2.

Proof. We only prove the case of i = 1. By the first equation of (1.1), we know

that lim supt→∞x(t) ≤ K. So we assume x(t) ≤ K for t large enough. It is easy

to see that

e1m1K ≤ d1K < d1(a1+ K).

Let µ = d1−ea11m+K1K > 0. According to the monotonicity of the function

e1m1x a+x with respect to x, we have ˙y1 y1 = e1m1x a1+ x + by1 − d1 < e1m1x a1+ x − d1 ≤ e1m1K a1+ K − d1 = −µ.

This implies that lim supt→∞y1(t) = 0. We complete the proof.

From now on we always assume that (H2) e1m1 > d1 and e2m2 > d2. Hence eimiK

ai+K − di < 0 if and only if K < λi if (H2) holds.

(c) E1 = (¯x1, ¯y1, 0). The semi-trivial equilibrium E1 is a boundary equilibrium on

the (x, y1)-plane, where ¯x1, ¯y1 are obtained by restricting to the system of the first

predator y1 and the prey x. The Jacobian matrix at E1 is given by

J(E1) =     ¯ x1(−Kr +(a1xm11+by¯11y¯1)2) − m1x¯1(a1+¯x1) (a1+¯x1+b1y¯1)2 − m2x¯1 a2+¯x1 e1m1y¯1(a1+b1¯y1) (a1+¯x1+b1y¯1)2 − b1e1m1x¯1y¯1 (a1+¯x1+b1y¯1)2 0 0 0 e2m2x¯1 a2+¯x1 − d2     .

We note that the top left 2 × 2 submatrix is exactly the Jacobian matrix J in (2.6)

for the subsystem (2.1) at the equilibrium (x∗, y∗) where a, b, m, d are replaced

by a1, b1, m1, d1 (The conditions are presented in Table I). And ea22m2x11 − d2 < 0

if and only if ¯x1 < λ2 if the assumption (H2) holds. There are four cases for the

(9)

Case A1: The equilibrium E1 is asymptotically stable in R3 if (¯x1, ¯y1) is an asymptotically stable equilibrium for the system (2.1) where a, b, m, d

replaced by a1, b1, m1, d1 (The conditions are presented in Table I) and

e2m2x¯1

a2+¯x1 − d2 < 0.

Case A2: If (¯x1, ¯y1) is an asymptotically stable equilibrium for the system (2.1)

and ¯x1 > λ2, then E1 is a saddle with 1-dimensional unstable manifold W1u

and 2-dimensional stable manifold (x, y1) plane

Case A3: If (¯x1, ¯y1) is an unstable focus for the system (2.1) and ¯x1 < λ2, then

E1 is a saddle with 1-dimensional stable manifold W1s and a unique limit

cycle Γ1 on the (x, y1) plane.

Case A4: If (¯x1, ¯y1) is an unstable focus for the system (2.1) and ¯x1 > λ2, then

E1 is a repeller.

(d) E2 = (¯x2, 0, ¯y2). Similar to the above case (c), the Jacobian matrix at E2 is

given by J(E1) =     ¯ x2(−Kr +(a2xm22+by¯22y¯2)2) − m1¯x2 a1+¯x2 − m2x¯2(a2+¯x2) (a2+¯x2+b2y¯2)2 0 e1m1¯x2 a1+¯x2 − d1 0 e2m2y¯2(a2+b2¯y2) (a2+¯x2+b2y¯2)2 0 − b2e2m2x¯2y¯2 (a2+¯x2+b2y¯2)2     .

We note that the 2 × 2 submatrix gotten by deleting the second row and second column of above matrix is exactly the Jacobian matrix J in (2.6) for the subsystem

(2.1) at the equilibrium (x∗, y∗) where a, b, m, d are replaced by a2, b2, m2, d2.

We have four cases:

Case B1: The equilibrium E2 is asymptotically stable in R3 if (¯x2, ¯y2) is an

asymptotically stable equilibrium for the system (2.1) where a, b, m, d re-placed by a2, b2, m2, d2 and ¯x2 < λ1.

Case B2: If (¯x2, ¯y2) is an asymptotically stable equilibrium for the system (2.1)

and ¯x2 > λ1, then E2 is a saddle with 1-dimensional unstable manifold W2u

(10)

Case B3: If (¯x2, ¯y2) is an unstable focus for the system (2.1) and ¯x2 < λ1, then

E2 is a saddle with 1-dimensional stable manifold W2s and a unique limit

cycle Γ2 on the (x, y2) plane.

Case B4: If (¯x2, ¯y2) is an unstable focus for the system (2.1) and ¯x2 > λ1, then

E2 is a repeller.

(e) EΓ1 = (φ1, ψ1, 0). If the condition in Proposition 2.2 (ii) is satisfied, then the

equilibrium ¯E = (¯x1, ¯y1) on the (x, y1) plane is unstable and there is a unique

stable limit cycle Γ1 on the (x, y1) plane, denoted by (φ1(t), ψ1(t)). Consequently,

EΓ1 = (φ1, ψ1, 0) is a boundary periodic solution for the system (1.1). Since EΓ1 is

stable restricted to the (x, y1) plane, we only need to discuss its stability in y2-axis

direction.

The stability of EΓ1 is determined by the Floquet multipliers of the variational

system

˙Φ(t) = J(φ1, ψ1, 0)Φ(t), Φ(0) = I (2.10)

where J(x, y1, y2) is defined in (2.9) and I is the 3×3 identity matrix. Let ω1be the

period of the periodic solution (φ1, ψ1). Then the Floquet multiplier corresponding

to the y2-direction is given by

exp[ 1 ω1 Z ω1 0 (m2e2φ1(t) a2+ φ1(t) − d2) dt]. Thus EΓ1 is stable if d2 > Z ω1 0 m2e2φ1(t) a2+ φ1(t) dt (2.11) and unstable if d2 < Z ω1 0 m2e2φ1(t) a2+ φ1(t) dt . (2.12)

(f) Similarly, if the boundary periodic solution EΓ2 = (φ2(t), 0, ψ2(t)) with period

ω2 exists then it is stable if

d1 > Z ω2 0 m1e1φ2(t) a1+ φ2(t) dt (2.13)

(11)

and unstable if d1 < Z ω2 0 m1e1φ2(t) a1+ φ2(t) dt . (2.14)

We now have the following results on the uniform persistence of system (1.1). (Bulter et. al [4], Butler and Waltman [3], Freedman et. al [8], Smith and Thieme [20])

Theorem 2.4. Assume one of the following cases holds:

(i) Let Case A2 and Case B2 holds, i.e., E1 and E2 are unstable in y2-axis

and y1-axis direction, respectively.

(ii) Let Case A2, Case B4 and (2.14) hold, i.e., E1 and EΓ2 are unstable in

y2-axis and y1-axis direction, respectively.

(iii) Let Case B2, Case A4 and (2.12) hold, i.e., E2 and EΓ1 are unstable in

y1-axis and y2-axis direction, respectively.

(iv) Let Case A4, (2.12), Case B4 and (2.14) hold, i.e., EΓ1 and EΓ2 are

un-stable in y2-axis and y1-axis direction, respectively.

Then system (1.1) is uniformly persistent.

(g) Ec = (xc, y1c, y2c). From the 2nd and 3rd equations of (1.1), xc, y1c, y2c satisfy

eimix ai+ x + biyi

= di (2.15)

for i = 1, 2 or

y1c= M1(xc− λ1) > 0, y2c= M2(xc− λ2) > 0 (2.16)

where we use the notations M1 = e1md11b1d1 and M2 = e2md22b2d2 for simplifying.

Assume

(12)

From the first equation of (1.1), xc satisfies the equation rx(1 − x K) − d1 e1 M1(x − λ1) − d2 e2 M2(x − λ2) = 0. Let F (x) = rx(1 − x K) − d1 e1 M1(x − λ1) − d2 e2 M2(x − λ2). Then F (K) < 0, F (0) > 0, F (λ1) > 0, and F (λ2) = rλ2(1 − λ2 K) − d1 e1 M1(λ2− λ1). Hence if F (λ2) > 0 (2.17)

then Ec = (xc, y1c, y2c) exists and is unique. If

F (λ2) < 0 (2.18)

then Ec does not exist. Rewrite

F (x) = (− r K)x 2 + x r − d1 e1 M1− d2 e2 M2 + d1 e1 M1λ1+ d2 e2 M2λ2.

Then xc is the unique positive root of F (x) = 0,

xc = K(B +pB2 + 4rC/K) 2r (2.19) where B = r −d1 e1M1− d2 e2M2 and C = d1 e1M1λ1+ d2

e2M2λ2. The condition (2.17) for

the existence of Ec is equivalent to

K > λ2 1 − d1 re1λ2 M1(λ2− λ1) −1 = ˜K > 0 (2.20)

or xc > λ2 . We note that in (2.20) we need

r > d1 e1 M1(1 − λ1 λ2 ) (2.21)

The Jacobian matrix of the system (1.1) at Ec takes the form

J(Ec) =    f∗ x f ∗ y1 f ∗ y2 g∗ x g ∗ y1 0 h∗ x 0 h ∗ y2   

(13)

where f∗ x = xc − r K + m1y1c (a1+ xc+ b1y1c)2 + m2y2c (a2 + xc+ b2y2c)2  f∗ y1 = − m1xc(a1+ xc) (a1+ xc+ b1y1c)2 < 0 f∗ y2 = − m2xc(a2+ xc) (a2+ xc+ b2y2c)2 < 0 g∗ x = e1m1y1c(a1+ b1y1c) (a1+ xc+ b1y1c)2 > 0 (2.22) g∗ y1 = − b1e1m1xcy1c (a1+ xc+ b1y1c)2 < 0 h∗ x = e2m2y2c(a2+ b2y2c) (a2+ xc+ b2y2c)2 > 0 h∗ y2 = − b2e2m2xcy2c (a2+ xc+ b2y2c)2 < 0 .

The characteristic polynomial of J(Ec) is given by

λ3+ α1λ2 + α2λ + α3 = 0 where α1 = −(f ∗ x + g ∗ y1+ h ∗ y2), α2 = f ∗ xg ∗ y1 + f ∗ xh ∗ y2 + g ∗ y1h ∗ y2 − f ∗ y2h ∗ x− f ∗ y1g ∗ x , α3 = fy∗1g ∗ xh ∗ y2 + f ∗ y2g ∗ y1h ∗ x− f ∗ xg ∗ y1h ∗ y2 .

By Routh-Hurwitz criterion we have the following result on the local stability of Ec.

Proposition 2.5. Assume that

α1 > 0, α3 > 0, and α1α2 > α3

then Ec is locally asymptotically stable.

Remark 2.2. If f∗

x < 0, then α1 > 0 and α2 > 0. From equations (2.22), (2.15),

and (2.16), f∗ x < 0 if and only if r Kxc > m1xcy1c (a1+ xc+ b1y1c)2 + m2xcy2c (a2+ xc+ b2y2c)2 = ( d1 e1m1 ) m1y1c a1+ xc+ b1y1c + ( d2 e2m2 ) m2y2c a2+ xc + b2y2c .

(14)

Then ( d1 e1m1 ) m1y1c a1+ xc + b1y1c + ( d2 e2m2 ) m2y2c a2 + xc+ b2y2c ≤ max{ d1 e1m1 , d2 e2m2 } m1y1c a1+ xc+ b1y1c + m2y2c a2+ xc + b2y2c  = max{ d1 e1m1 , d2 e2m2 }r(1 −xc K). If r Kxc > max{ d1 e1m1 , d2 e2m2 }r(1 − xc K) or equivalent ¯ M 1 + ¯MK < xc < K, where ¯M = max{ d1 e1m1, d2 e2m2}, then f ∗ x < 0.

### Hopf Bifurcation

In this section, we will verify that the Hopf bifurcation indeed occurs. It is obvious

(15)

expressions of α1 and α3 α1 = −(f ∗ x + g ∗ y1 + h ∗ y2), = rxc K − m1xcy1c (b1y1c+ xc+ a1)2 − m2xcy2c (b2y2c+ xc+ a2)2 + b1e1m1xcy1c (b1y1c+ xc+ a1)2 + b2e2m2xcy2c (b2y2c+ xc+ a2)2 , α3 = f ∗ y1g ∗ xh ∗ y2 + f ∗ y2g ∗ y1h ∗ x− f ∗ xg ∗ y1h ∗ y2 = b2e1e2m1 2m 2xc2 (xc+ a1) y1c (b1y1c+ a1) y2c (b1y1c+ xc+ a1)4(b2y2c+ xc+ a2)2 + b1e1e2m1m22xc2 (xc+ a2) y1cy2c (b2y2c+ a2) (b1y1c+ xc + a1)2(b2y2c+ xc + a2)4 + b1b2e1e2m1m2xc2y1cy2c  r xc K − m2xcy2c (b2y2c+xc+a2)2 − m1xcy1c (b1y1c+xc+a1)2  (b1y1c+ xc + a1)2(b2y2c+ xc+ a2)2 = a1b2e1e2m1 2m 2xc2y1cy2c (b1y1c+ xc+ a1)3(b2y2c+ xc+ a2)2 + a2b1e1e2m1m2 2x c2y1cy2c (b1y1c+ xc+ a1)2(b2y2c+ xc+ a2)3 + rb1b2e1e2m1m2xc3y1cy2c K (b1y1c+ xc + a1)2 (b2y2c+ xc+ a2)2 > 0 .

Hence, by Proposition 2.5, the positive equilibrium Ec will lose its stability if

α1α2− α3 ≤ 0. We take K as the bifurcation parameter. It is easy to see that xc,

y1c, and y2c are functions of K by the equations (2.19) and (2.16). The expression

of α1α2− α3 has the form,

α1α2− α3 = −(f ∗ x + g ∗ y1 + h ∗ y2)(f ∗ xg ∗ y1 + f ∗ xh ∗ y2 + g ∗ y1h ∗ y2 − f ∗ y2h ∗ x− f ∗ y1g ∗ x)− (f∗ y1g ∗ xh ∗ y2 + f ∗ y2g ∗ y1h ∗ x− f ∗ xg ∗ y1h ∗ y2) = −(f∗ x) 2g y1 − (f ∗ x) 2h y2 − (g ∗ y1) 2h y2 + f ∗ y1g ∗ xg ∗ y1 − g ∗ y1(h ∗ y2) 2+ f∗ y2h ∗ xh ∗ y2+ f∗ x f ∗ y2h ∗ x+ f ∗ y1g ∗ x− (g ∗ y1) 2 − 2g∗ y1h ∗ y2 − (h ∗ y2) 2 .

In the last formula, we have two classes

−(f∗ x)2gy1 − (f ∗ x)2hy2 − (g ∗ y1) 2h y2 + f ∗ y1g ∗ xg ∗ y1 − g ∗ y1(h ∗ y2) 2+ f∗ y2h ∗ xh ∗ y2 and f∗ f∗ h∗ + f∗ g∗ − (g∗ )2− 2g∗ h∗ − (h∗ )2 .

(16)

All terms of the first class are positive and all term of another one are negative

except for the function f∗

x. So we should clarify the behavior of f

x as a function

of K.

By the representation of xc, (2.19), it is easy to see that if B = r − de11M1 −

d2

e2M2 > 0 then limK→0

+xc(K) = 0, limK→0+xc(K)/K > 0, limK→∞xc(K) = ∞,

and limK→∞xc(K)/K = B/r > 0. These implies limK→0+fx∗(K) < 0. But

the restriction of K, (2.20), it is required that f∗

c( ˜K) < 0. It is easy to see that

d1 e1M1+ d2 e2M2 > d1 e1M1(1− λ1

λ2) which is the restriction of r to guarantee the existence

of Ec in (2.21), so we assume r > de11M1 + de22M2. A necessary condition for the

occurrence of Hopf bifurcation is limK→∞f

x(K) > 0. Easy computation shows

that lim k→∞f ∗ c(K) = −r + d1M1 e1 + d2M2 e2 + m1M1 (1 + b1M1)2 + m2M2 (1 + b2M2)2 . Hence we assume (H4) 0 < r − d1 e1 M1− d2 e2 M2 < m1M1 (1 + b1M1)2 + m2M2 (1 + b2M2)2 . Proposition 2.6. Assume the assumption (H4) holds and

(i) b1e1 ≥ 1 and b2e2 ≥ 1,

(ii) there is a K∗

> 0 such that α1(K∗)α2(K∗) = α3(K∗) and

d dK K=K∗α1(K)α2(K) < d dK K=K∗α3(K),

then the positive equilibriumEc is locally stable whenK < K

and loses its stability

when K = K∗

. When K > K∗

, Ec becomes unstable and a family of periodic

solutions bifurcates from Ec.

### Global Stability of Coexistence State

Using the Lyapunov function constructed in [9, 10] we give sufficient conditions

for the global stability of the positive equilibrium Ec.

(17)

Lemma 3.1. The solutions of (1.1) are positive and bounded for t ≥ 0. The proof of Lemma 3.1 is similar to that in [12].

Theorem 3.2. Let the assumption (H3) hold. Assume Ec exists, i.e., (2.20) and

(2.21) hold. If

K < 1

max{1/a1, 1/a2}

+ xc (3.1)

then the positive equilibrium Ec is globally stable.

Proof. Choose a Lyapunov function as follows

V (x, y1, y2) = Z x xc ξ − xc ξ dξ + α Z y1 y1c ξ − y1c ξ dξ + β Z y2 y2c ξ − y2c ξ dξ ,

where α and β are positive constants to be determined. Along the trajectories of the system (1.1) we have

dV dt = (x − xc)  r(1 − x K) − m1y1 a1+ x + b1y1 − m2y2 a2+ x + b2y2  + α(y1− y1c) m1e1x a1+ x + b1y1 − d1 + β(y2− y2c) m2e2x a2+ x + b2y2 − d2  = (x − xc) n − r K(x − xc) − m1y1 a1+ x + b1y1 − m1y1c a1+ xc+ b1y1c − m2y2 a2+ x + b2y2 − m2y2c a2+ xc+ b2y2c o + α(y1− y1c) m1e1x a1+ x + b1y1 − m1e1xc a1 + xc+ b1y1c  + β(y2− y2c) m2e2x a2+ x + b2y2 − m2e2xc a2+ xc+ b2y2c  = (x − xc) n − r K(x − xc) − m1 (a1+ xc)(y1− y1c) − y1c(x − xc)  (a1+ x + b1y1)(a1+ xc + b1y1c) − m2 (a2+ xc)(y2− y2c) − y2c(x − xc)  (a2+ x + b2y2)(a2+ xc + b2y2c) o + α(y1− y1c) m1e1 (a1+ b1y1c)(x − xc) − b1xc(y1− y1c)  (a1+ x + b1y1)(a1+ xc+ b1y1c) + β(y2− y2c) m2e2 (a2+ b2y2c)(x − xc) − b2xc(y2− y2c)  (a2+ x + b2y2)(a2+ xc + b2y2c) .

(18)

Choose α = a1+xc e1(a1+b1y1c) and β = a2+xc e2(a2+b2y2c). Therefore, dV dt = (x − xc) 2 − r K + m1y1c (a1+ x + b1y1)(a1+ xc+ b1y1c) + m2y2c (a2+ x + b2y1)(a2 + xc+ b2y2c) − αb1xc(y1− y1c) 2 (a1+ x + b1y1)(a1+ xc + b1y1c) − βb2xc(y2− y2c) 2 (a2+ x + b2y2)(a2+ xc + b2y2c) .

The coefficients of (y1−y1c)2and (y2−y2c)2are negative. The coefficient of (x−xc)2

is − r K + m1y1c (a1+ x + b1y1)(a1+ xc+ b1y1c) + m2y2c (a2+ x + b2y1)(a2+ xc+ b2y2c) ≤ − r K + m1y1c a1(a1 + xc+ b1y1c) + m2y2c a2(a2+ xc+ b2y2c) ≤ − r K + max{ 1 a1 , 1 a2 }r(1 −xc K) = − r K 1 − max{ 1 a1 , 1 a2 }(K − xc).

If (3.1) is satisfied, then dV /dt ≤ 0 and dV /dt = 0 if and only if x = xc, y1 = y1c,

and y2 = y2c. The largest invariant set of {dV /dt = 0} is {(xc, y1c, y2c)}. Therefore,

from Lemma 3.1 and LaSalle’s Invariant Principle implies that Ec = (xc, y1c, y2c)

is globally stable. Thus we complete the proof.

Remark 3.1. Under the assumption (H2) and (2.20), (2.21), Ec exists and xc > λ2.

Let ˜K = λ2(1 − reλ12me1b1−d1(λ2− λ1))−1. If r is sufficient large then ˜ K < 1 max{1/a1, 1/a2} + λ2 < 1 max{1/a1, 1/a2} + xc .

Thus the condition (3.1) is feasible when r is sufficiently large.

### Dif-ferent Interference Effects

In this section we consider two identical predators competing for a shared prey

(19)

following: x′ = rx(1 − x K) − mxy1 a + x + b1y1 − mxy2 a + x + b2y2 , y1 ′ = ( emx a + x + b1y1 − d)y1 , (4.1) y2 ′ = ( emx a + x + b2y2 − d)y2 ,

with initial conditions x(0) > 0, y1(0) > 0, y2(0) > 0. Let

K > λ1 = λ2 = a/( em d − 1). (4.2) Assume b2 > b1. Then y1′ = ( emx a + x + b1y1 − d)y1 > ( emx a + x + b2y1 − d)y1 .

Thus, if y1(0) ≥ y2(0) then y1(t) > y2(t) for all t ≥ 0. If y1(0) < y2(0) then either

there exists t0 > 0 such that y1(t0) = y2(t0) or y1(t) < y2(t) for all t ≥ 0. If

y1(t0) = y2(t0) then

y1(t) > y2(t) for all t ≥ t0. (4.3)

If y1(t) < y2(t) for all t ≥ 0 then y1′ y1 = emx a + x + b1y1 − d > emx a + x + b2y2 − d = y2 ′ y2 . Thus we have y1(t) y1(0) > y2(t) y2(0) . (4.4)

Thus, we have either y1(t0) > y2(t0) for some t0 > 0 or y2(0)y1(t) > y1(0)y2(t) for

all t ≥ 0. If y1(t) → 0 as t → ∞ then y2(t) → 0 as t → ∞. Hence we obtain a

contradiction to the assumption (4.2). Hence lim sup

t→∞

(20)

On the other hand, assume y2(t) → 0 as t → ∞. Let Case A1 hold. Then x(t) → ¯x1 and y1(t) → ¯y1 as t → ∞ and a+¯em¯x1+bx11y¯1 = d. Thus

em¯x1 a + ¯x1

− d > 0 . (4.6)

Let Case A3 hold. Then (x(t), y1(t)) → (φ1(t), ψ1(t)) as t → ∞ and

Z ω1 0 emφ1(t) a + φ1(t) + b1ψ1(t) − d dt = 0. Hence Z ω1 0 emφ1(t) a + φ1(t) − d dt > 0 . (4.7)

However (4.6) and (4.7) imply that E1and EΓ1 are unstable in the y2-axis direction

respectively. Thus the assumption y2(t) → 0 as t → ∞ leads to a contradiction.

Hence we have the following results.

Theorem 4.1. For system (4.1), if (4.2) holds then lim supt→∞y1(t) > 0 and

lim supt→∞y2(t) > 0.

### Relaxation Oscillations

Consider system (1.1) with large prey intrinsic growth rate, i.e., r ≫ 1. Let ε = 1/r. Then 0 < ε ≪ 1, With the scaling:

x → x/K, a1 → a1/K, a2 → a2/K, y1 = y1/(Kr), y2 = y2/(Kr), system (1.1) becomes εx′ = x(1 − x) − m1xy1 a1+ x + (bε1)y1 − m2xy2 a2 + x + (bε2)y2 y1 ′ = ( e1m1x a1+ x + (bε1)y1 − d1)y1 (5.1) y2 ′ = ( e2m2x a2+ x + (bε2)y2 − d2)y2

Assume b1 = b1(ε), b2 = b2(ε) such that

(21)

for some µi > 0, i = 1, 2. Under the assumption (5.2) we apply the geometric singular perturbation method as in [17] to prove the existence of periodic solutions.

Setting ε = 0 in (5.1) results in the so-called limiting slow system

xF (x, y1, y2) = x 1 − x − m1y1 a1+ x − m2y2 a2+ x , y1′ = ( e1m1x a1+ x − d1)y1, (5.3) y2 ′ = (e2m2x a2+ x − d2)y2,

which is generally defined on the slow manifold S0 = {(x, y1, y2) : xF (x, y1, y2) =

0, x ≥ 0, y1 ≥ 0, y2 ≥ 0}. Orbits or parts of orbits of the system (5.3) on S0

are called the slow orbits of system (5.1) and the variables y1, y2 are called slow

variables. For system (5.3), the slow manifold S0 consists of two portions S1 and

S2, where S1 = {(x, y, z) ∈ S0 : x = 0}, S2 = {(x, y1, y2) : F (x, y1, y2) = 0}.

In term of the fast time scale τ = t/ε, system (5.1) becomes dy1 dτ = εy1( e1m1x a1+ x + (bε1)y1 − d1), dy2 dτ = εy2( e2m2x a2+ x + (bε2)y2 − d2), (5.4) dx dτ = x 1 − x − m1y1 a1+ x + (bε1)y1 − m2y2 a2+ x + (bε2)y2 .

The system (5.5) is referred to as the fast system. Its limit, the limiting fast system, is obtained by setting ε = 0:

dy1 dτ = 0, dy2 dτ = 0, dx dτ = xF (x, y1, y2). (5.5)

The orbits of system (5.5) are parallel to the x-axis and their directions are

char-acterized by the sign of xF (x, y1, y2). We refer to these orbits as fast orbits of

system (5.1) and the variable x is the fast variable.

A continuous and piecewise smooth curve is said to be a limiting orbit of system (5.1) if it is the union of a finitely many fast and slow orbits with compatible orientations. A limiting orbit is called a limiting periodic orbit if it is a simple closed curve and contains no equilibrium of system (5.1). A periodic orbit of system (5.1) is called a relaxation oscillation if its limiting as ε → 0 is a limiting periodic orbit consisting of both fast and slow orbits.

(22)

Table II: Parameter Values in the General Case.

r = 2.0 a1 = 3 b1 = 0.6 d1= 0.4 e1 = 0.6 m1 = 1.5

K = ∗ a2 = 6 b2 = 2.0 d2 = 0.45 e2 = 0.7 m2 = 1.5

In the following theorem, we first prove that under the assumption (5.2) there is no positive equilibrium for system (5.1). Then following the methods in [17] we construct a limiting periodic orbit consisting of both fast and slow orbits. By the theorem of geometric singular perturbation method, there exists a stable relaxation oscillation.

Theorem 5.1. Let (H3) and (5.2) hold. Assume that the relaxation cycle Γε

1 on

the (x, y1)-plane is unstable in the y2-axis direction and the relaxation cycle Γε2 on

the (x, y2)-plane is unstable in the y1-axis direction. Then there is at least one

stable relaxation oscillation in the positive octant of R3.

Proof. If Eε

c = (xεc, y1cε , y2cε ) exists then from (2.16) and (5.2), y1cε → ∞ as ε → 0.

Thus the equilibrium Eε

c is not on the surface S0 and the limiting periodic orbit

does not contain Eε

c. From Theorem 3.4 in [17], there exists a stable relaxation

oscillation in the positive octant of R3. We complete the proof.

### Numerical Simulations

Choose the values of parameters as in Table II and calculate the values λ1 =

a1d1

e1m1−d1 = 2.4 and λ2 =

a2d2

e2m2−d2 = 4.5. Now, using K (the carrying capacity of

the resource) as a bifurcation parameter, increase K from 4.5 to 80 and calculate f∗

x(K) as a function of K in (2.22). We can see that f

x is monotonically increasing

from negative to positive (see the first graph of Fig. I). The values of functions

α1, α3, and α1α2 − α3 are also calculated (see the 2nd - 4th graphs of Fig. I).

The dynamics of solutions with respect to the capacity K are illustrated in Figure II(a)-(e).

(i) 0 < K = 2 < λ1. The semi-trivial equilibrium EK is globally asymptotically

(23)

0 10 20 30 40 50 60 70 80 −1.5 −1.0 −0.5 0.0 0.5 fx 0 10 20 30 40 50 60 70 80 0.0 0.2 0.4 0.6 0.8 1.0 a1 0 10 20 30 40 50 60 70 80 −0.02 −0.010.00 0.01 0.02 0.03 0.04 a3 0 10 20 30 40 50 60 70 80 K −0.050.00 0.05 0.10 0.15 0.20 0.25 0.30 a1 a2 -a 3

Figure I: The graphs of f∗

x(K), α1(K), α3(K) and α1(K)α2(K) − α3(K) in terms

of K as K increases from 4.5 to 80.

(ii) λ1 < K = 3 < λ2. The semi-trivial equilibrium E1 is globally asymptotically

stable, (see Fig. II:(b))

(iii) λ2 < K = 10. The solution converges to the positive equilibrium Ec as

t → ∞. We can see that the positive equilibrium is asymptotically stable, (see Fig. II:(c))

(iv) K = 75. The positive equilibrium Ec loses its stability and a periodic solution

bifurcates from it. (see Fig. II:(d))

Next, we do some numerical simulations of system (1.1) with interference

ef-fects, i.e., b1 6= 0 and b2 6= 0. In order to compare the differences of solutions of

system (1.1) with or without interference effects, we choose the same parameters as those in Fig. 3 of [11] in Table III. We plot limit cycles of population of predator

(24)

Table III: Parameter Values for the Case with Interference.

r = 20 · ln 2 a1 = 200 d1 = ln 2/2 e1 = 0.1 m1 = 10 · ln 2

K = 1100 a2 = 500 d2 = ln 2 e2 = 1.4 m2 = 2 · ln 2

.

b1 = 0, b2 = 1, and (c) is for b1 = 1, b2 = 0. All above three limit cycles are

plotted in a graph showed in (d). With the same parameters, we compute the

numerical solutions of (1.1) with various parameters b1 and b2. In Fig III (e) and

(f) show the numerical results where b1, b2 are varied from 0 to 10 with step-size

0.1 in (e) and b1, b2 are varied from 0 to 1 with step-size 0.01 in (e). The white

region represents that the solutions are periodic and the black region means that the solutions approach a positive equilibrium.

### Discussion

In systems where there is no intraspecific or interspecific interference but only ex-ploitative competition for resources between consumers, the conventional wisdom is that the number of consumer species which can coexist is less than or equal to the number of distinct resources. This has been shown in a number of mod-els including chemostat modmod-els (Smith and Waltman [21]). In the case of two consumers, our results indicate that sufficiently strong intraspecific and mutual interference can ensure the long term survival of both competing consumers not only in terms of uniform persistence but also in the sense of global stability. A possible explanation of this phenomenon is that intraspecific feeding interference in one consumer reduces its equilibrium density and allows the other consumer to have better access to the resource. It is also interesting to observe that when the carrying capacity K of the resource is increased, the interior equilibrium loses stability and a three dimensional positive periodic solution arises via Hopf bifur-cation. This indicates that the paradox of enrichment phenomenon may occur for the two-predators-one-prey model as well.

(25)

### References

[1] R. A. Armstrong and R. McGehee. Competitive exclusion. American Natu-ralist, 115(2):151–170, 1980.

[2] J. R. Beddington. Mutual Interference Between Parasites or Predators and its Effect on Searching Efficiency. The Journal of Animal Ecology, 44(1):331, Feb. 1975.

[3] G. Butler. Persistence in dynamical systems. J Differential Equations, 1986. [4] G. Butler, H. I. Freedman, and P. Waltman. Uniformly persistent systems.

Proceedings of the American Mathematical Society, 96(3):425–430, 1986.

[5] G. J. Butler and P. Waltman. Bifurcation from a limit cycle in a two predator-one prey ecosystem modeled on a chemostat. Journal of Mathematical Biology, 12(3):295–310, 1981.

[6] R. S. Cantrell and C. Cosner. On the dynamics of predator-prey models with the Beddington-DeAngelis functional response. Journal of Mathematical Analysis and Applications, 257(1):206–222, 2001.

[7] D. L. DeAngelis, R. A. Goldstein, and R. V. O’neill. A model for tropic interaction. Ecology, 56(4):881–892, 1975.

[8] H. I. Freedman, S. Ruan, and M. Tang. Uniform Persistence and Flows Near a Closed Positively Invariant Set. Journal of Dynamics and Differential Equations, 6(4):583–600, 1994.

[9] S. B. Hsu. Limiting Behavior for Competing Species. SIAM Journal on Applied Mathematics, 34(4):760–763, June 1978.

[10] S. B. Hsu. A survey of constructing Lyapunov functions for mathematical models in population biology. Taiwanese Journal of Mathematics, 9(2):151– 173, 2005.

[11] S. B. Hsu, S. P. Hubbell, and P. Waltman. A contribution to the theory of competing predators. Ecological Monographs, 48(3):337–349, 1978.

(26)

[12] S. B. Hsu, S. P. Hubbell, and P. Waltman. Competing predators. SIAM Journal on Applied Mathematics, 35(4):617–625, 1978.

[13] G. Huisman and R. J. DeBoer. A formal derivation of the ”Beddington” functional response. Journal of Theoretical Biology, 185(3):389–400, 1997. [14] T.-W. Hwang. Global analysis of the predator-prey system with

Beddington-DeAngelis functional response. Journal of Mathematical Analysis and

Appli-cations, 281(1):395–401, 2003.

[15] T.-W. Hwang. Uniqueness of limit cycles of the predator-prey system with Beddington-DeAngelis functional response. Journal of Mathematical Analysis

and Applications, 290(1):113–122, 2004.

[16] J. P. Keener. Oscillatory coexistence in the chemostat: a codimension two unfolding. SIAM Journal on Applied Mathematics, 43(5):1005–1018, 1983. [17] W. Liu, D. Xiao, and Y. Yi. Relaxation oscillations in a class of predator-prey

systems. Journal of Differential Equations, 188(1):306–331, 2003.

[18] S. Muratori and S. Rinaldi. Remarks on competitive coexistence. SIAM Journal on Applied Mathematics, 49(5):1462–1472, 1989.

[19] H. L. Smith. The interaction of steady state and Hopf bifurcations in a two-predator-one-prey competition model. SIAM Journal on Applied Mathe-matics, 42(1):27–43, 1982.

[20] H. L. Smith and H. R. Thieme. Dynamical systems and population persistence, volume 118 of Graduate Studies in Mathematics. American Mathematical Society, Providence, RI, 2011.

[21] H. L. Smith and P. Waltman. The Theory of the Chemostat, volume 13 of

Cambridge Studies in Mathematical Biology. Cambridge University Press,

(27)

0 50 100 150 200 time (t) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 x y1 y2 (a) 0 50 100 150 200 time (t) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 x y1 y2 (b) 0 50 time (t)100 150 200 0 1 2 3 4 5 6 7 8 x y1 y2 (c) 0 50 time (t)100 150 200 0 10 20 30 40 50 60 70 x y1 y2 (d)

Figure II: The parameters are given in Table II. In Fig II (a), K = 2, EK =

(K, 0, 0) is globally asymptotically stable. In Fig II (b), K = 3, E1 = (¯x1, ¯y1, 0) is

globally asymptotically stable. In Fig II (c), K = 10, Ec = (xc, y1c, y2c) is globally

asymptotically stable. In Fig II (d), K = 75, the periodic solution exist. Hopf bifurcation occurs between K = 70 and K = 75.

(28)

180 200 220 240 260 280 300 320 340 360 500 1000 1500 2000 2500 3000 3500 4000 4500 b_1=b_2=0 y1 y2 (a) 100 200 300 400 500 600 700 800 900 0 50 100 150 200 250 300 b_1=0, b_2=1 y1 y2 (b) 62 64 66 68 70 4500 5000 5500 6000 b_1=1, b_2=0 y1 y2 (c) 0 100 200 300 400 500 600 700 800 900 0 1000 2000 3000 4000 5000 6000 b_1=0, b_2=0 b_1=1, b_2=0 b_1=0, b_2=1 y1 y2 (d) 0.0 0.2 0.4 0.6 0.8 1.0 b_2 0.0 0.2 0.4 0.6 0.8 1.0 b_1 (e)

Figure III: The parameters are given in Table III. The graphs in Fig III (a), (b),

(c) are the limit cycle solutions of system (1.1) projected in (y1, y2)-plane with

b1 = b2 = 0 in Fig III(a), b1 = 0, b2 = 1 in Fig III (b), b1 = 1, b2 = 0 in Fig III

(c). We put Fig III (a), (b), (c) in the same graph in Fig III (d). In Fig III (e),

in the b1-B2 parameter region, 0 ≤ b1, b2 ≤ 1, the white region represents that the

numerical solutions are periodic and the black region represents that the numerical 28