On the Dynamics of Two-consumers-one-resource
Competing Systems with Beddington-DeAngelis
Functional Response
Sze-Bi Hsu
∗1, Shigui Ruan
†2, and Ting-Hui Yang
‡31Department 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.
1
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
(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.
2
Local Analysis
2.1
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.
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)
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.
2.2
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).
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, m2e2−d2
de(me−d)+mre2}
λ < K < K∗
(¯x, ¯y) is globally asymptotically stable
and b < min{1e, m2e2−d2
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)
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 ,
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 +(a1+¯xm11+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 ea22m+¯2xx¯11 − d2 < 0
if and only if ¯x1 < λ2 if the assumption (H2) holds. There are four cases for the
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 +(a2+¯xm22+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
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)
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 = e1md11b−1d1 and M2 = e2md22b−2d2 for simplifying.
Assume
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
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 .
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.
2.3
Hopf Bifurcation
In this section, we will verify that the Hopf bifurcation indeed occurs. It is obvious
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 .
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.
3
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.
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) .
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.
4
Competition of Two Identical Species with
Dif-ferent Interference Effects
In this section we consider two identical predators competing for a shared prey
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→∞
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.
5
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
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.
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.
6
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
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
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.
7
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.
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.
[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,
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.
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
