# Boundedness and Global Stability for a Predator-Prey System with the Beddington-DeAngelis Functional Response

Authors: Wahiba Khellaf; Nasreddine Hamri
Journal: Differential Equations and Nonlinear Mechanics
 (2010)

---

## Abstract

We study the qualitative behavior of a class of predator-prey models with Beddington-DeAngelis-type functional response, primarily from the viewpoint of permanence (uniform persistence). The Beddington-DeAngelis functional response is similar to the Holling type-II functional response but contains a term describing mutual interference by predators. We establish criteria under which we have boundedness of solutions, existence of an attracting set, and global stability of the coexisting interior equilibrium via Lyapunov function.


---

## Body

## 1. Introduction and Mathematical Model

Standard Lotka-Volterra systems are also known as the predator-prey systems, on which a large body of existing predator-prey theory is built by assuming that the per capita rate of predation depends on the prey numbers only [1]. Recently, the traditional prey-dependent predator-prey models have been challenged by several biologists based on the fact that functional and numerical response over typical ecological time scale sought to depend on the densities of both prey and predator, especially when predators have to search, share, or compete for food. A more suitable general predator-prey model should be based on the ratio-dependent theory [2–4]. This roughly states that the per capita predator growth rate should be a function of the ratio of prey to predator abundance. Moreover, as the number of predators often changes slowly (relative to prey number), there is often a competition among the predators, and the per capita rate of predation should therefore depend on the numbers of both prey and predator, most probably and simply on their ratios. These hypotheses are strongly supported by numerous field and laboratory experiments and observations [5–9].The general model describing the dynamics of prey-predator populations in continuous time can be written as(1.1)dXdτ=f(X)X-g(X,Y)Y,dYdτ=h(X,Y)Y,
where X and Y are the densities (or biomasses) of prey and of predators at time T, respectively. f(X) is the per capita net prey production in the absence of predation, whereas g(X,Y) is the functional response of predators (the number of preys eaten per predator per unit time). Natural mortality of prey is considered to be negligible compared to mortality due to predation.The function h(X,Y) represents the numerical response of predators (measures the growth rate of predators). The function f will be taken either as the Malthusian growth f(X)=rX or as the logistic model f(X)=rX(1-X/K). The key role in prey-predator models is played by the functional response g (Solomon 1949). Traditionally, it is assumed that the functional response is a function of prey density only prey-dependent feeding,g=g(X), without any dependence on predator density [9, 10]. The hypothesis is based on an analogy with the law of mass action in chemistry assuming that prey and predator individuals encounter each other randomly in space and time [11]. Therefore, the prey-dependent model can be applied to systems which are spatially homogeneous and in which the time scale of prey removal by predators is of the same order of magnitude as that of population reproduction [2].Many questions in predator-prey theory, including the question of interference between predators, revolve around the expression that is used for the functional responseg=g(X,Y). Arditi and Ginzburg [2] have argued that, in many cases, this predator dependence could be simplified as a ratio-dependent model g=g(X/Y) instead of modeling explicitly all conceivable interference mechanisms (and thus adding parameters to the model).The Beddington-DeAngelis-type functional response performed even better. Although the predator-dependent models that they considered fit those data reasonably well, no single functional response best describes all the data sets. The Beddington-DeAngelis response can be generated by a number of natural mechanisms [5, 12], and because it admits rich but biologically reasonable dynamics [6], it is worthy for us to further study the Beddington-DeAngelis model.Therefore, it is interesting and important to study the following autonomous predator-prey model with the Beddington-DeAngelis functional response:(1.2)dXdτ=(a1-b1X-m1Yα1X+β1Y+γ1)X,dYdτ=(a2-m2YX+k1)Y,
with the initial values X(0)>0 and Y(0)>0. The constants a1,a2,b1,m1,m2,α1,β1,γ1, and k1 are the parameters of model and are assumed to be nonnegative with β1 nontrivial (if β1=0, then the model (1.2) is the same as that in [13]).These parameters are defined as follows:a1 (resp., a2) describes the growth rate of prey (resp., of predator), b1 measures the strength of competition among individuals of prey's species, m1 is the maximum value which per capita reduction rate of prey can attain, γ1 (resp., k1) measures the extent to which environment provides protection to prey (resp., to predator), and m2 has a similar meaning to m1. The functional response in (1.2) was introduced by Beddington [5] and DeAnglis et al. (1975) in [7]. It is similar to the well-known Holling type-II functional response but has an extra term β1Y in the first right term equation modeling mutual interference among predators. Hence this kind of type functional response given in (1.2) is affected by both predator and prey, that is, the so-called predator dependence by Arditi and Ginzburg [2]. Dynamics for the Holling type-II model have been much studied (see, e.g., [13–15]). Then how the mutual interference term affects the dynamic of the whole system is an interesting problem.Introducing the following scaling (see [16]), t=a1τ,x(t)=(b1/a1)X(τ), and y(t)=(m2b1/a1a2)Y(τ), then Beddington-DeAngelis predator-prey model (1.2) should take the following nondimensional form:
(1.3)dxdt=x(1-x)-axyαx+βy+γ,dydt=b(1-yx+k)y,
where a=(a2/a1)(m1/m2),b=a2/a1,α=α1,β=β1(a2/m2),γ=γ1(b1/a1), and k=k1(b1/a1).


## 2. Boundedness of the Model and Existence of a Positively Invariant Attracting Set

We denote byℝ+2 the nonnegative quadrant, and by Int(ℝ+2) the positive quadrant.Lemma 2.1.
Positive quadrantInt(ℝ+2) is invariant for system (1.3).Proof.
From system (1.3), we observe that the boundaries of the nonnegative quadrant ℝ+2, which are the positives x-axis and y-axis are invariant; this is immediately obvious from the system (1.3). Therefore, densities x(t) and y(t) are positive for all t≥0 if x(0)>0 and y(0)>0. Theorem of existence and uniqueness ensures that the positive solutions of the autonomous system (1.3) and the axis cannot intersect.Next, we will show that, under some assumptions, the solutions of system (1.3) which start in ℝ+2 are ultimately bounded. First, let us give the following comparison lemma.Lemma 2.2.
Letϕ be an absolutely-continuous function satisfying the differential inequality
(2.1)ddtϕ(t)+α1ϕ(t)≤α2,t≥0,
where (α1,α2)∈ℝ2; then
(2.2)∀t≥T¯:ϕ(t)≤α1α2-(α1α2-ϕ(T¯))e-α1(t-T¯).We define the functionψ(t) to be the solution of differential equation
(2.3)ddtϕ(t)+α1ϕ(t)=α2,t≥0,
and we have the Gronwall's lemma.Lemma 2.3.
Letx(t) satisfy for t≥t0 the linear, scalar equation
(2.4)dxdt=a(t)x+b(t),x(t0)=x0,
with a(t) and b(t) being continuous functions. If y(t) satisfies for t≥t0 the inequalities
(2.5)dydt≤a(t)x+b(t),y(t0)≤x0,
then
(2.6)y(t)≤x(t),t≥t0.Definition 2.4.
A solutionϕ(t,t0,x0,y0) of system (1.3) is said to be ultimately bounded with respect to ℝ+2 if there exists a compact region A∈ℝ+2 and a finite time T(T=T(t0,x0,y0)) such that, for any (t0,x0,y0)∈ℝ×ℝ+2,(2.7)ϕ(t,t0,x0,y0)∈A,∀t≥T.Theorem 2.5.
LetA be the set defined by
(2.8)A={(x,y)∈ℝ+2:0≤x≤1,0≤x+y≤L1},
where
(2.9)L1=14b{5b+(1+b)2(1+k)}.
Then(1)
A is positively invariant; and(2)
all solutions of (1.3) initiating in ℝ+2 are ultimately bounded with respect to ℝ+2 and eventually enter the attracting set A.Proof.
Let(x(0),y(0))∈A, we will show that (x(t),y(t))∈A for all t≥0. Obviously, from Lemma 2.1, as (x(0),y(0))∈A,(x(t),y(t)) is in Int(ℝ+2). Then, we have to show that for all t≥0,0≤x(t)≤1, and 0≤x(t)+y(t)≤L1.(1)
(a)
First, we prove that for allt≥0,0≤x(t)≤1. We have x>0 and y>0 in Int(ℝ+2); then every solution ϕ(t)=(x(t),y(t)) of system (1.3), which starts in Int(ℝ+2), satisfies the differential inequality dx/dt≤(1-x(t))x(t). This is obvious by considering the first equation of (1.3). Thus, x(t) may be compared with solutions of du/dt=(1-u(t))u(t) and u(0)=x(0)>0 which is a Bernoulli's equation; then the solution is u(t)=1/(1+ce-t), with c=1/x(0)-1, where 0≤x(0)≤1, which implies that c=1/x(0)-1≥1.It follows that every nonnegative solution ϕ(t) satisfies
(2.10)x(t)≤1∀t≥0.(b)
We prove now that, for allt≥0,0≤x+y≤L1.We define the function σ(t)=x(t)+y(t); the time derivative of this function is
(2.11)dσdt=dxdt+dydt=(1-x-ayαx+βy+γ)x+b(1-yx+k)y.
Since all parameters are positive and solutions initiating in ℝ+2 remain in the nonnegative quadrant, then
(2.12)dσdt≤(1-x)x+b(1-yx+k)y
holds for all x and y being nonnegative. Thus, as
(2.13)maxℝ+(1-x)x=14,
we have
(2.14)dσdt≤14+b(1-yx+k)y;
then
(2.15)dσdt≤14+b(1-yx+k)y+σ(t)-σ(t),
which implies that
(2.16)dσdt+σ(t)≤14+x+(b+1-byx+k)y.
Since in A, x(t)≤1 for all t≥0, we obtain
(2.17)dσdt+σ(t)≤54+(b+1-by1+k)y.
Moreover, it can be easily verified that
(2.18)maxℝ+[(b+1-by1+k)y]=14b(1+b)2(1+k).
Consequently
(2.19)dσdt+σ(t)≤L1,
with
(2.20)L1=14b5b+(1+b)2(1+k).
Using Lemma 2.2, with
(2.21)α1=1,α2=L1,
then we get
(2.22)∀t≥T¯≥0:σ(t)≤L1-(L1-σ(T¯))e-(t-T¯).
Then, if T¯=0,(2.23)σ(t)≤L1-(L1-σ(0))e-t.
Hence, for (x(0),y(0))∈A,(2.24)σ(t)=x(t)+y(t)≤L1.
Then
(2.25)(x(t),y(t))∈A,∀t≥0.(2)
We have to prove that, for(x(0),y(0))∈ℝ+2,(x(t),y(t))→A when t→+∞. We will show that lim¯t→+∞x(t)≤1 and lim¯t→+∞(x(t)+y(t))≤L1.
(a)
The first result,lim¯t→+∞x(t)≤1, is obtained directly from (1a) and Lemma 2.2, since solutions of the initial value problem
(2.26)dxdt=x(t)(1-x(t)),x(0)≥0,
verify
(2.27)lim¯t→+∞x(t)≤1.(b)
For the second result, letε>0 be given, and T1>0 exists, such that
(2.28)x(t)≤1+ε2∀t≥T1.
From (2.22) with T¯=T1, we get, for all t≥T1,(2.29)σ(t)=x(t)+y(t)≤L1-(L1-σ(T1))e-(t-T1)≤L1-{L1eT1-(x(T1)+y(T1))eT1}e-t≤L1-{L1-(x(T1)+y(T1))eT1}e-t.
Then
(2.30)σ(t)=x(t)+y(t)≤(L1+ε2)-{(L1+ε2)-(x(T1)+y(T1))eT1}e-t,∀t≥T1≥0.
Let T2≥T1 be such that
(2.31)(L1+ε2)-(x(T1)+y(T1))eT1e-t≤ε2,∀t≥T2.
Then
(2.32)x(t)+y(t)≤L1+ε∀t≥T2.
Hence,
(2.33)lim¯t→+∞(x(t)+y(t))≤L1.
This completes the proof; then we conclude that system (1.3) is dissipative in ℝ+2.


## 3. Linear Stability

First of all, it is easy to verify that this system has three trivial equilibria, (belonging to the boundary ofℝ+2, i.e., at which one or more of populations has zero density or is extinct) (3.1)P1(0,0),P2(0,k),P3(1,0).The other equilibria are defined by the system(3.2)ayαx+βy+γ=1-x,y=x+k.Proposition 3.1.
The system (1.3) has a unique interior equilibria P*(x*,y*) (i.e., x*>0 and y*>0) if the following condition is verified:
(3.3)k(a-β)≤γ.Proof.
We introduce the second equation of (3.2) in the first one; then
(3.4)a(x+k)=(1-x)((α+β)x+(βk+γ)),
and we obtain
(3.5)(α+β)x2+(βk+γ+a-α-β)x+ak-βk-γ=0.The discriminant of this equation(3.6)Δ=(βk+γ+a-α-β)2-4(ak-βk-γ)(α+β)=((βk+γ+a)-(α+β))2+4((βk+γ)-ak)(α+β)=(βk+γ+a)2+(α+β)2-2(βk+γ+a)(α+β)+4(βk+γ+a)(α+β)-4a(k+1)(α+β)=((βk+γ+a)+(α+β))2-4a(1+k)(α+β).
Therefore, if (3.3) holds, then
(3.7)Δ=((βk+γ)+a+(α+β))2-4a(1+k)(α+β)≥(a(k+1)+(α+β))2-4a(1+k)(α+β)≥(a(k+1)-(α+β))2≥0.
Consequently, Δ is positive, and the system (1.3) has two other equilibriums P1*(x1,y1) and P2*(x2,y2), where
(3.8)x1,2=((βk+γ)-(α+β)+a)±Δ2(α+β),y1,2=(x1,2+k).
Now, we show, under the condition (3.3), that one of this equilibriums is not in ℝ+2; let
(3.9)x2=((βk+γ+a)-(α+β))-Δ2(α+β);
then
(3.10)2x2=1-(βk+γ+a)+Δ(α+β),
from (3.7)
(3.11)Δ≥|a(k+1)-(α+β)|,
and due to (3.3)
(3.12)2x2≤1-(a(k+1)+|a(k+1)-(α+β)|)(α+β),
which implies that(1)
ifa(k+1)<(α+β)(3.13)2x2≤1-a(k+1)-a(k+1)+(α+β)(α+β)≤0;(2)
ifa(k+1)>(α+β),
(3.14)2x2≤1-a(k+1)+a(k+1)-(α+β)(α+β)≤2-2a(k+1)α+β≤0.it results that P2*(x2,y2) is not in ℝ+2, such that
(3.15)x1x2=ak-βk-γ≤0;
then the first point P1*(x1,y1) is in ℝ+2.The Jacobian matrix is given by(3.16)J(Xi)=(∂f1∂x(xi,yi)∂f1∂y(xi,yi)∂f2∂x(xi,yi)∂f2∂y(xi,yi)),
where (3.17)∂f1∂x(xi,yi)=1-2x-ay(βy+γ)(αx+βy+γ)2,∂f1∂y(xi,yi)=-ax(αx+γ)(αx+βy+γ)2,∂f2∂x(xi,yi)=b(yx+k)2,∂f2∂y(xi,yi)=b-2byx+k.(1)
AtP0(0,0),(3.18)J(P0)=(100b).
The eigenvalues of this matrix are
(3.19)λ1=1,λ2=b.
Hence, all parameters are positive; then P0(0,0) is an unstable node.(2)
AtP1(0,k),(3.20)J(P1)=(1-akβk+γ0b-b).
The eigenvalues are
(3.21)λ1=1-akβk+γ,λ2=-b<0.
then, we have
(a)
ifk(a-β)>γ, P1(0,k) is stable node;(b)
ifk(a-β)≤γ,P1(0,k) is unstable with the positivey-axis as its stable manifold.(3)
AtP2(1,0),(3.22)J(P2)=(-1-aα+γ0b).
The eigenvalues are
(3.23)λ1=-1<0,λ2=b>0.
Then the equilibrium P2(1,0) is a saddle point with the stable manifold being the x-axis.AroundP*(x*,y*), the Jacobian matrix takes the form(3.24)J(P*)=(1-2x*-ay*(βy*+γ)(αx*+βy*+γ)2-ax*(αx*+γ)(αx*+βy*+γ)2b-b).The characteristic equation is(3.25)λ2-trJ(P*)λ+detJ(P*)=0,
where (3.26)detJ(x*,y*)=-b(1-2x*-ay*(βy*+γ)(αx*+βy*+γ)2)+abx*(αx*+γ)(αx*+βy*+γ)2=b(αx*+βy*+γ)2{ay*(βy*+γ)+(2x*-1)(αx*+βy*+γ)2+ax*(αx*+γ)(αx*+βy*+γ)2}=b(αx*+βy*+γ)2{(x*-1)(αx*+βy*+γ)2+x*(αx*+βy*+γ)2+ay*(βy*+γ)+ax*(αx*+γ)(αx*+βy*+γ)2}.From (3.2), we get
(3.27)detJ(P*)=b(αx*+βy*+γ)2{-ay*(αx*+βy*+γ)+x*(αx*+βy*+γ)2+ay*(βy*+γ)+ax*(αx*+γ)(αx*+βy*+γ)2}=b(αx*+βy*+γ)2{-aαx*y*+x*(αx*+βy*+γ)2+ax*(αx*+γ)}=b(αx*+βy*+γ)2{x*(αx*+βy*+γ)2+ax*(αx*-αy*+γ)}=b(αx*+βy*+γ)2{x*(αx*+βy*+γ)2-ax*(αk-γ)};
then(3.28)detJ(P*)=bx*(αx*+βy*+γ)2{(αx*+βy*+γ)2-a(αk-γ)}.We observe thatdetJ(P*) is positive if
(3.29){(αx*+βy*+γ)2-a(αk-γ)}>0.
To simplify, we developed detJ(P*) respecting one variable, from (3.2); then
(3.30)(αx*+βy*+γ)2-a(αk-γ)=(a(x*+k)1-x*)2-a(αk-γ)=a2(x*+k)2(1-x*)2-a(αk-γ),
which implies that detJ(P*) has the same sign of
(3.31)a(x*+k)2+(1-x*)2(γ-αk).
We rewrite
(3.32)a(x*)2+2kax*+ak2+(γ-αk)((x*)2-2x*+1).
Let (3.33)f(x)=(a+γ-αk)x2+2(ak+αk-γ)x+(γ-αk+ak2).The discriminant is(3.34)Δ′=(ak+αk-γ)2-(a+γ-αk)(γ-αk+ak2)=(a+α)2k2-2γk(a+α)+γ2-(a+γ)(ak2+γ)+(ak2+a+2γ)αk-(αk)2=(a2+2aα+α2)k2-2γk(a+α)+γ2-(a2k2+a(1+k2)γ+γ2)+aαk3+γ+(a+2γ)αk-α2k2=2aαk2-2aγk-aγk2-aγ+aαk3+aαk=a{αk3+(2α-γ)k2+(α-γ)k-γ}=a{αk(k2+2k+1)-γ(k2+2k+1)}=a(αk-γ)(k+1)2.We get three cases.(1)
Ifαk<γ,Δ′ is negative, f(x) has the same sign of (a+γ-αk), and we have (a+γ-αk)>a>0.Then,detJ(P*) is positive.(2)
Ifαk>γ,Δ′ is positive and f(x) has at list two solutions x1 and x2, then
(a)
ifx∈]-∞,x1]∪[x2,+∞[, f(x) has the same sign of (a+γ-αk);(b)
ifx∈]x1,x2[ has sign of -(a+γ-αk).(3)
Ifαk=γ, then
(3.35)f(x)=ax2+2akx+ak2=a(x+k)2>0.
Then, detJ(P*) is positive.Remark 3.2.
From the expression (3.28), we find that det(J(P*)) is positive, if αk≤γ, hence the eigenvalues associated to P* have the same sign.To determine the sign of these eigenvalues, we have(3.36)trJ(P*)=1-2x*-ay*(βy*+γ)(αx*+βy*+γ)2-b=1(αx*+βy*+γ)2{(1-2x*-b)(αx*+βy*+γ)2-ay*(βy*+γ)}=1(αx*+βy*+γ)2{(1-x*)(αx*+βy*+γ)2-(x*+b)(αx*+βy*+γ)2-ay*(βy*+γ)(αx*+βy*+γ)2}.
From (3.2), we get
(3.37)trJ(P*)=1(αx*+βy*+γ)2{ay*(αx*+βy*+γ)-(x*+b)(αx*+βy*+γ)2-ay*(βy*+γ)(αx*+βy*+γ)2}=1(αx*+βy*+γ)2{aαx*y*-(x*+b)(αx*+βy*+γ)2}=αx*ay*(αx*+βy*+γ)2-(x*+b)=αx*ay*(1-x*)2a2(y*)2-(x*+b)=αx*(1-x*)2ay*-(x*+b)=1a(x*+k){αx*(1-x*)2-a(x*+b)(x*+k)}=1a(x*+k){α(x*)3-(a+2α)(x*)2-(a(b+k)-α)x*-abk}.
Let
(3.38)P3(x)=αx3-(a+2α)x2-(a(b+k)-α)x-abk.
We obtain the following lemma.Lemma 3.3.
Ifαk<γ is verified, the interior equilibrium P*(x*,y*) is locally asymptotically stable if P3(x*)>0 and it is unstable if P3(x*)<0.We used theCardan's method to solve the cubic equation P3(x)=0. Then we consider the equation
(3.39)a3x3+a2x2+a1x+a0=0,
with a3=α,a2=-(a+2α),a1=-(a(b+k)-α), and a0=abk. Making the substitution y=a3x+a2/3 reduces the equation to the standard form
(3.40)y3-py-q=0,
where p and q depend on a3,a2,a1, and a0(3.41)p=-a1a3+a223,q=-a0a32-2(a23)3+a3a2a13.
Let
(3.42)y=u+v.
Then
(3.43)y3-py-q=(u+v)3-3uv(u+v)-(u3+v3)=0.
We get after the identification of the coefficients
(3.44)3uv=p,u3+v3=q.
Then
(3.45)u3v3=(p3)3,u3+v3=q,
and we obtain that u3 and v3 are solutions of the quadratic equation
(3.46)z2-qz+(p3)3=0.
Then we constitute three cases.(1)
if27q2-4p3>0, then
(3.47)u3=q2+(q2)2-(p3)3,v3=q2-(q2)2-(p3)3.
We have
(3.48)y=q2-1227q2-4p3273+q2+1227q2-4p3273;
then (3.39) has at list one real root
(3.49)r0=1a3(q2-1227q2-4p3273+q2+1227q2-4p3273-a23).
So, we have P3(x)<0 if 0<x<r0, and P3(x)>0 if r0<x.(2)
if27q2-4p3<0, then there are three real roots r1,r2, and r3, and therefore, r1r2r3=abk>0 then one of them is positive, then
(a)
ifr2<r3<0<r1, then P3(x)>0 if 0<r1<x, and P3(x)<0 if x<r1;(b)
if0<r1<r2<r3, then P3(x)>0 if 0<r1<x, and P3(x)<0 if x<r1.(3)
if27q2=4p3, then there are one real root positive r0=2q/23 and a double root r1,2=-q/23; we also have P3(x)<0 if 0<x<r0, and P3(x)>0 if r0<x.


## 4. Uniform Permanence

In this section we shall prove the permanence [8, 17–19], that is, the uniform persistence and dissipativity, of system (1.3).The principal notion of persistence theory is uniform persistence or permanence. Before the study of the permanence of system (1.3), we introduce some necessary definitions. Consider an ODE model for n interacting biological species(4.1)dxidt=fi(x1,x2,…,xn),i=1,2,…,n,
where xi(t) denotes the density of the ith species. Let (x1(t),x2(t),…,xn(t)) denote the solution of (4.1) with componentwise positive initial values. The system (4.1) is said to be weakly persistent if(4.2)lim supt→+∞xi(t)>0,i=1,2,…,n,
persistent if(4.3)lim inft→+∞xi(t)>0,i=1,2,…,n,
and uniformly persistent if there is an ε0>0 such that(4.4)lim inft→+∞xi(t)≥ε0,i=1,2,…,n.
The system (4.1) is said to be permanent if for each i=1,2,…,n there are constants ε0 and Mi such that(4.5)0<ε0≤lim inft→+∞xi(t)≤lim supt→+∞xi(t)≤Mi.
Clearly, a permanent system is uniformly persistent which in turn is persistent, and persistence implies weak persistence; a dissipative uniformly persistent system is permanent. For further discussion about various definitions of persistence and permanence and their connections, see [18].Suppose thatY is a complete metric space with Y=Y0∪∂Y0 for an open set Y0. We will choose Y0 to be the positive cone in ℝ2. For the following definitions and theorems, one can see [6], and for the proof of the theorem, see [8].Definition 4.1.
A flow or semiflow onY under which Y0 and ∂Y0 are forward invariant is said to be permanent if it is dissipative and if there is a number ε>0 such that any trajectory starting in Y0 will be at least a distance ε from ∂Y0 for all sufficiently large t.Letω(∂Y0)∈∂Y0 denote the union of the sets ω(u) over u∈∂Y0.Definition 4.2.
Theω-limit set ω(∂Y0) is said to be isolated if it has a covering Ω=⋃k=1NΩk of pairwise disjoint sets Ωk which are isolated and invariant with respect to the flow or the semiflow both on ∂Y0 and on Y=Y0∪∂Y0, (M is called an isolated covering). The set ω(∂Y0) is said to be acyclic if there exists an isolated covering ⋃k=1NΩk such that no subset of Ωk is a cycle.Theorem 4.3 (Hale and Waltman 1989).
Suppose that a semiflow onY leaves both Y0 and ∂Y0 forward invariant, maps bounded sets in Y to precompact set for t>0, and it is dissipative. If in addition(1)
ω(∂Y0) is isolated and acyclic;(2)
Ws(Ωk)∩Y0=∅ for all K, where ⋃k=1NΩk is the isolated covering used in the definition of acyclicity of ∂Y0, and Ws denote the stable manifold.Then the semiflow is permanent.And, we have this theorem.Theorem 4.4.
Let us assume the following condition:(4.6)k(a-β)≤γ.
Then, system (1.3) is permanent.Proof.
We takeY the strictly positive quadrant of R2; then ω(∂Y0) consists of the equilibria P0(0,0),P1(0,k), and P2(1,0). P0(0,0) is an unstable node, P2(1,0) is saddle point, and its stable manifold is x-axis. If ak≤βk+γ,P1(0,k) is a saddle point stable along the y-axis and unstable along the x-axis.
Then, all trajectories on the axis(ox) other than P0(0,0) approach the point P2(1,0) and all trajectories on the axis (oy) other than P0(0,0) approach the point P1(0,k). It follows from these structural features that the flow in ∂Y0 is acyclic. So ω(∂Y0) is isolated and acyclic. The stable manifold of P2(1,0) is the x-axis and the stable manifold of P1(0,k) is the y-axis, and we know, from Theorem 2.5, that these stable manifolds cannot intersect the interior of Y0.
In this case, Theorem4.3 implies permanence of the flow defined by (1.3).


## 5. Global Stability

In this section, we shall prove the global stability of system (1.3) by constructing a suitable Lyapunov function. First, we have to show that there exists one interior equilibrium P*(x*,y*). The linear analysis shows that if αk≤γ and x*<x(0), then P*(x*,y*) is locally stable. We prove now that, under some assumptions, this steady state is globally asymptotically stable.Theorem 5.1.
The interior equilibriumP*(x*,y*) is globally asymptotically stable if
(5.1)β<α,(5.2)2aL1<γ,(5.3)(βk+γ)(αL1+γ)(1+k)<4βk2γ2,(5.4)α(1+2k)-β(1+k)<γ,(5.5)a-4βL1<4γ.Proof.
The proof is based on construction of a positive definite Lyapunov function. Let(5.6)V(x,y)=V1(x,y)+V2(x,y),
where
(5.7)V1(x,y)=(αx*+βy*+γ)(x-x*-x*ln(xx*)),V2(x,y)=ab(x*+k)(y-y*-y*ln(yy*)).
This function is defined and continuous on Int(ℝ+2).
We can easily verify that the functionV(x,y) is zero at the equilibrium (x*,y*) and is positive for all other positive values of x and y, and thus, P*(x*,y*) is the global minimum of V.
Since the solutions of (1.3) are bounded and ultimately enter the set A, we restrict the study for this set. The time derivative of V1 and V2 along the solutions of system (1.3) is
(5.8)dVidt=dVidxẋ+dVidyẏ,i=1,2.
Then
(5.9)dV1dt(x,y)=(αx*+βy*+γ)(1-x*x)ẋ,dV2dt(x,y)=ab(x*+k)(1-y*y)ẏ,
and using (3.2), we get after simplifing
(5.10)dV1dt(x,y)=-(αx*+βy*+γ)(x-x*)2-(x-x*)(y-y*)(aαx+aγαx+βy+γ)+(aαyαx+βy+γ)(x-x*)2,dV2dt(x,y)=-a(y-y*)2+ayx+k(x-x*)(y-y*).
Therefore, computing dV/dt via dV1/dt and dV2/dt yields
(5.11)dVdt=(-(αx*+βy*+γ)+(aαyαx+βy+γ))(x-x*)2+(x-x*)(y-y*)(-(aαx+aγαx+βy+γ)+ayx+k)-a(y-y*)2.
The above equation can be written as
(5.12)dVdt(x,y)=-(x-x*,y-y*)(-g(x,y)-h(x,y)-h(x,y)a)(x-x*y-y*),
where
(5.13)g(x,y)=-(αx*+βy*+γ)+aαyαx+βy+γ,h(x,y)=a2(-αx+γαx+βy+γ+yx+k).
From (5.11), it is obvious that dV/dt<0 if the matrix above is positive definite. This matrix is positive definite if only if all upper-left submatrices are positive (Sylvester's criteria), that is, since a>0, if only if(1)
g(x,y)<0;(2)
ϕ(x,y)=ag(x,y)+h2(x,y)<0.
Proof.
It is of (1)(5.14)g(x,y)=-(αx*+βy*+γ)+aαyαx+βy+γ<0.
So, asA is an attracting positively invariant set, where, all solutions satisfy 0≤x≤1 and 0≤x+y≤L1, then
(5.15)g(x,y)≤-αx*+aαyαx+βy+γ,g(x,y)≤α(-1+ay*αx*+βy*+γ+ayαx+βy+γ),g(x,y)≤α(-1+2aL1γ).
Therefore, if (5.2) holds, then
(5.16)g(x,y)<0,∀(x,y)∈A.Proof.
It is of (2)(5.17)ϕ(x,y)=-a(αx*+βy*+γ)+a2αyαx+βy+γ+a24(-αx+γαx+βy+γ+yx+k)2<0.
Since (for x fixed)
(5.18)∂ϕ∂y=a2α(αx+γ)(αx+βy+γ)2+a22(-αx+γαx+βy+γ+yx+k)(β(αx+γ)(αx+βy+γ)2+1x+k),
then
(5.19)∂2ϕ∂y2=-2a2αβ(αx+γ)(αx+βy+γ)3+a22(β(αx+γ)(αx+βy+γ)2+1x+k)2+a22(-2β2(αx+γ)(αx+βy+γ)3)(-αx+γαx+βy+γ+yx+k)=-2a2αβ(αx+γ)(αx+βy+γ)3+a2β2(αx+γ)2(αx+βy+γ)4+a22(x+k)2+a2β(αx+γ)(αx+βy+γ)2(x+k)-a2β2y(αx+γ)2(αx+βy+γ)3(x+k).
We have
(5.20)∂2ϕ∂y2≤-2a2αβ(αx+γ)(αx+βy+γ)3-a2β2y(αx+γ)2(αx+βy+γ)3(x+k)+a2β2(αx+γ)2+a22(x+k)2+a2β(αx+γ)(x+k)≤-a2β(αx+γ)2(αx+βy+γ)3(x+k)(4α(x+k)+βy)+a2β2(αx+γ)2+a22(x+k)2+a2β(αx+γ)(x+k).
We note, using (5.1), that for (x,y) in A(5.21)1a2∂2ϕ∂y2≤-2kβ2(αL1+γ)3(1+k)+β2γ2+12k2+βkγ.
If (5.3) holds, then
(5.22)∂2ϕ∂y2≤0.
Hence, ∂ϕ/∂y is strictly decreasing in ℝ+, with respect to y.
Now(5.23)∂ϕ∂y|y=0=a2ααx+γ-a22(βαx+γ+1x+k)=a22(αx+γ)(x+k)((α-β)x+k(2α-β)-γ).
In A, all solutions satisfy 0≤x≤1, and from (5.1)
(5.24)(α-β)x+k(2α-β)-γ≤(α-β)+k(2α-β)-γ,
then, if (5.4) holds, (∂ϕ/∂y)|y=0≤0 in ℝ+. Hence, ϕ(x,y) is strictly decreasing in ℝ+. This yields ϕ(x,y)<D(x,0) for (x,y)∈A; that is, using (5.1), we get
(5.25)ϕ(x,y)<-a(αx*+βy*+γ)+a24<-a(β(x*+y*)+γ-a4)<-a(βL1+γ-a4).
Consequently, due to (5.5),
(5.26)ϕ(x,y)<0,∀(x,y)∈A.It follows that if the hypotheses of Theorem 5.1 are satisfied, then dV/dt<0 along all trajectories in the first quadrant except (x*,y*); so P*(x*,y*) is globally asymptotically stable.


## 6. Conclusion

The Beddington-DeAngelis functional response admits a range of dynamics which include the possibilities of extinction, persistence, and stable or unstable equilibria. The criteria for persistence are the same as for systems with a Holling-type 2 response.The future research will complete the qualitative analysis by studying the limit cycles of the model. It will also contain the numerical simulations to justify the obtained results.



---

*Source: 813289-2010-08-24.xml*