# Stability in a Simple Food Chain System with Michaelis-Menten Functional Response and Nonlocal Delays

Authors: Wenzhen Gan; Canrong Tian; Qunying Zhang; Zhigui Lin
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

This paper is concerned with the asymptotical behavior of solutions to the reaction-diffusion system under homogeneous Neumann boundary condition. By taking food ingestion and species' moving into account, the model is further coupled with Michaelis-Menten type functional response and nonlocal delay. Sufficient conditions are derived for the global stability of the positive steady state and the semitrivial steady state of the proposed problem by using the Lyapunov functional. Our results show that intraspecific competition benefits the coexistence of prey and predator. Furthermore, the introduction of Michaelis-Menten type functional response positively affects the coexistence of prey and predator, and the nonlocal delay is harmless for stabilities of all nonnegative steady states of the system. Numerical simulations are carried out to illustrate the main results.


---

## Body

## 1. Introduction

The overall behavior of ecological systems continues to be of great interest to both applied mathematicians and ecologists. Two species predator-prey models have been extensively investigated in the literature. But recently more and more attention has been focused on systems with three or more trophic levels. For example, the predator-prey system for three species with Michaelis-Menten type functional response was studied by many authors [1–4]. However, the systems in [1–4] are either with discrete delay or without delay or without diffusion. In view of individuals taking time to move, spatial dispersal was dealt with by introducing diffusion term to corresponding delayed ODE model in previous literatures, namely, adding a Laplacian term to the ODE model. In recent years, it has been recognized that there are modelling difficulties with this approach. The difficulty is that diffusion and time delay are independent of each other, since individuals have not been at the same point at previous times. Britton [5] made a first comprehensive attempt to address this difficulty by introducing a nonlocal delay; that is, the delay term involves a weighted-temporal average over the whole of the infinite domain and the whole of the previous times.There are many results for reaction-diffusion equations with nonlocal delays [5–18]. The existence and stability of traveling wave fronts were studied in reaction-diffusion equations with nonlocal delay [5–9]. The stability of impulsive cellular neural networks with time varying was discussed in [10] by means of new Poincare integral inequality. The asymptotic behavior of solutions of the reaction-diffusion equations with nonlocal delay was investigated in [11, 12] by using an iterative technique and in [13–15] by the Lyapunov functional. The stability and Hopf bifurcation were discussed in [16] for a diffusive logistic population model with nonlocal delay effect.Motivated by the work above, we are concerned with the following food chain model with Michaelis-Menten type functional response:(1)∂u1∂t-d1Δu1=u1(a1-a11u1-a12u2m1+u1),∂u2∂t-d2Δu2=u2(-a2+a21∫Ω‍∫-∞t‍u1(s,y)m1+u1(s,y)11111111111111111111×K1(x,y,t-s)dsdy11111111111111111111-a22u2-a23u3m2+u2(-a2+a21∫Ω‍∫-∞t‍u1(s,y)m1+u1(s,y)),∂u3∂t-d3Δu3=u3(-a3+a32∫Ω‍∫-∞t‍u2(s,y)m2+u2(s,y)11111111111111111111×K2(x,y,t-s)dsdy11111111111111111111-a33u3(-a2+a21∫Ω‍∫-∞t‍u1(s,y)m1+u1(s,y)),
for t>0, x∈Ω with homogeneous Neumann boundary conditions
(2)∂u1∂ν=∂u2∂ν=∂u3∂ν=0,t>0,x∈∂Ω,
and initial conditions
(3)ui(t,x)=ϕi(t,x)≥0(i=1,2),(t,x)∈(-∞,0]×Ω,u3(0,x)=ϕ3(x)≥0,x∈Ω,
where ϕi is bounded, Hölder continuous function and satisfies ∂ϕi/∂ν=0 (i=1,2,3) on (-∞,0]×∂Ω. Here, Ω is a bounded domain in ℝn with smooth boundary ∂Ω and ∂/∂ν is the outward normal derivative on ∂Ω. ui(t,x) represents the density of the ith species (prey, predator, and top predator resp.) at time t and location x and thus only nonnegative ui(t,x) is of interest. The parameter a1 is the intrinsic growth rate of the prey, and a2 and a3 are the death rates of the predator and top-predator. aii is the intraspecific competitive rate of the ith species. a12 is the maximum predation rate. a23 and a32 are the efficiencies of food utilization of the predator and top predator, respectively. We assume the predator and top predator show the Michaelis-Menten (or Holling type II) functional response with u1/(m1+u1) and u2/(m2+u2), respectively, where m1 and m2 are half-saturation constants. For a through biological background of similar models, see [18, 19]. As our most knowledge, the tritrophic food chain model has been found to have many interesting biological properties, such as the coexistence and the Hopf bifurcation. However, the effect of nonlocal time delays on the coexistence has not been reported. Our paper mainly concerns this perspective.Additionally,∫Ω‍∫-∞t‍(ui(s,y)/(mi+ui(s,y)))Ki(x,y,t-s)dsdy (i=1,2) represents the nonlocal delay due to the ingestion of predator; that is, mature adult predator can only contribute to the production of predator biomass. The boundary condition in (2) implies that there is no migration across the boundary of Ω.The main purpose of this paper is to study the global asymptotic behavior of the solution of system (1)–(3). The preliminary results are presented in Section 2. Section 3 contains sufficient conditions for the global asymptotic behaviors of the equilibria of system (1)–(3) by means of the Lyapunov functional. Numerical simulations are carried out to show the feasibility of the conditions in Theorems 8–10 in Section 4. Finally, a brief discussion is given to conclude this work.


## 2. Preliminary Results

In this section, we present several preliminary results that will be employed in the sequel.Lemma 1 (see [3]).
Leta and b be positive constants. Assume that ϕ,φ∈C1(a,+∞),φ(t)≥0, and ϕ is bounded from below. If ϕ′(t)≤-bφ(t) and φ′(t)≤K in [a,+∞) for some positive constant K, then limt→+∞φ(t)=0.The following lemma is the Positivity Lemma in [20].Lemma 2.
Letui∈C([0,T]×Ω¯)⋂‍C1,2((0,T]×Ω)(i=1,2,3) and satisfy
(4)uit-diΔui≥∑j=13‍bij(t,x)uj(t,x)1111111111+∑j=13‍cijuj(t-τj,x)in(0,T]×Ω,∂ui(t,x)∂ν≥0on(0,T]×∂Ω,ui(t,x)≥0in[-τi,0]×Ω,
where bij,cij∈C([0,T]×Ω¯),uit=(ui)t. If bij≥0 for j≠i and cij≥0 for all i,j=1,2,3. Then ui(t,x)≥0 on [0,T]×Ω¯. Moreover, if the initial function is nontrivial, then ui>0 in (0,T]×Ω¯.Lemma 3 (see [20]).
Letc^ and c~ be a pair of constant vector satisfying c~≥c^ and let the reaction functions satisfy local Lipschitz condition with Λ=〈c^,c~〉. Then system (1)–(3) admits a unique global solution u(t,x) such that
(5)c^≤u(t,x)≤c~,∀t>0,x∈Ω¯,
whenever c^≤ϕ(t,x)≤c~, (t,x)∈(-∞,0]×Ω.The following result was obtained by the method of upper and lower solutions and the associated iterations in [21].Lemma 4 (see [21]).
Letu(t,x)∈C([0,∞)×Ω¯)∩C2,1((0,∞)×Ω) be the nontrivial positive solution of the system
(6)ut-dΔu=Bu(t-τ,x)±A1u(t,x)200000000-A2u2(t,x),(t,x)∈(0,∞)×Ω,∂u∂ν=0,(t,x)∈(0,∞)×∂Ω,u(t,x)=ϕ(t,x)≥0,(t,x)∈(-τ,0)×Ω,
where A1≥0,B,A2,τ>0. The following results hold:(1)
ifB±A1>0, then u→(B±A1)/A2 uniformly on Ω¯ as t→∞;(2)
ifB±A1≤0, then u→0 uniformly on Ω¯ as t→∞.Throughout this paper, we assume that(7)Ki(x,y,t)=Gi(x,y,t)ki(t),x,y∈Ω,ki(t)≥0;∫Ω‍Gi(x,y,t)dx=∫Ω‍Gi(x,y,t)dy=1,t≥0;∫0+∞‍ki(t)dt=1,tki(t)∈L1((0,+∞);R),
where Gi(x,y,t) is a nonnegative function, which is continuous in (x,y)∈Ω¯×Ω¯ for each t∈[0,+∞) and measurable in t∈[0,+∞) for each pair (x,y)∈Ω¯×Ω¯.Now we prove the following propositions which will be used in the sequel.Proposition 5.
For any nonnegative initial function, the corresponding solution of system (1)–(3) is nonnegative.Proof.
Suppose that(u1,u2,u3) is a solution and satisfies ui∈C([0,T)×Ω¯)∩C1,2((0,T)×Ω) with T≤+∞. Choose 0<τ<T. Then
(8)u1t-d1Δu1=A1u1,0<t<τ,x∈Ω,u2t-d2Δu2=A2u2,0<t<τ,x∈Ω,u3t-d3Δu3=A3u3,0<t<τ,x∈Ω,∂ui∂ν=0,0<t<τ,x∈∂Ω,ui(t,x)≥0(i=1,2,3),t≤0,x∈Ω¯,
where
(9)A1=a1-a11u1-a12u2m1+u1,A2=-a2+a21∫Ω‍∫-∞t‍u1(s,y)m1+u1(s,y)111111111111111111×K1(x,y,t-s)dsdy111111111111111111-a22u2-a23u3m2+u3,A3=-a3+a32∫Ω‍∫-∞t‍u2(s,y)m2+u2(s,y)111111111111111111×K2(x,y,t-s)dsdy-a33u3.
It follows from Lemma 2 that ui≥0,(t,x)∈[0,τ]×Ω¯. Due to the arbitrariness of τ, we have ui≥0 for (t,x)∈[0,T)×Ω¯.Proposition 6.
System (1)–(3) with initial functions ϕi admits a unique global solution (u1,u2,u3) satisfying 0≤ui(t,x)≤Mi for i=1,2,3, where Mi is defined as
(10)M1=max{a1a11,∥ϕ1(t,x)∥L∞((-∞,0]×Ω¯)},M2=max{(a21-a2)M1-m1a2a22(m1+M1),∥ϕ2(t,x)∥L∞((-∞,0]×Ω¯)},M3=max{(a32-a3)M2-m2a3a22(m2+M2),∥ϕ3(x)∥L∞(Ω¯)}.Proof.
It follows from standard PDE theory that there exists aT>0 such that problem (1)–(3) admits a unique solution in [0,T)×Ω¯. From Lemma 2 we know that ui(t,x)≥0 for (t,x)∈[0,T)×Ω¯. Considering the first equation in (1), we have
(11)∂u1∂t-d1Δu1≤u1(a1-a11u1).
Using the maximum principle gives that u1≤M1. In a similar way, we have ui≤Mi for i=2,3. It is easy to see that (0,0,0) and (M1,M2,M3) are a pair of coupled upper and lower solutions to system (1)–(3) from the direct computation. In virtue of Lemma 3, system (1)–(3) admits a unique global solution (u1,u2,u3) satisfying 0≤ui(t,x)≤Mi for i=1,2,3. In addition, if ϕ1(t,x),ϕ2(t,x),ϕ3(x), (t,x)∈(-∞,0]×Ω is nontrivial, it follows from the strong maximum principle that ui(t,x)>0(i=1,2,3) for all t>0,x∈Ω¯.


## 3. Global Stability

In this section, we study the asymptotic behavior of the equilibrium of system (1)–(3). In the beginning, we show the existence and uniqueness of positive steady state.Let us consider the following equations:(12)a1-a11v1-a12v2m1+v1=0,-a2+a21v1m1+v1-a22v2=0.
A direct computation shows that the above equations have only one positive solution (v1*,v2*) if and only if
(13)H1:a1(a21-a2)>m1a2a11,
is satisfied. Taking H1 into account, we consider the equation a3=a32v2*/(m2+v2*), which corresponds to L3 in Figure 1.Figure 1
Graphs of the equations in (12). L3 is the boundary condition (a3=a32v2/(m2+v2)) and L2 corresponds to -a2+a21(v1/(m1+v1))-a22v2=0, while L1 corresponds to a1-a11v1-a12v2/(m1+v1)=0. Intersection point between L1 and L2 gives the positive solution (v1*,v2*) to (12). Parameter values are listed in the example in Section 4.It is clear that ifa3<a32v2*/(m2+v2*), the intersection point (v1*,v2*) always lies in L3. Let
(14)H2:a3m2<(a32-a3)v2*.
Suppose that H1 and H2 hold. We take u3 as a parameter and consider the following system:
(15)a1-a11u1-a12u2m1+u1=0,-a2+a21u1m1+u1-a22u2-a23u3m2+u2=0,-a3+a32u2m2+u2-a33u3=0.
When the parameter u3 is sufficiently small, the first two equations in (15) can be approximated by (12). Moreover, by continuously increasing the value of u3,L3 goes up, and meanwhile the intersection point between L1 and L2 also goes up. However, L3 goes up faster than L2, while L1 keeps still. In other words, there exists a critical value u3 such that the intersection point lies in L3, which implies that there is a unique positive solution E*(u1*,u2*,u3*) to (15), or equivalently (1)–(3).Lemma 7.
Assume thatH2 and G2 hold, and then the positive steady state E* of system (1)–(3) is locally asymptotically stable.Proof.
We get the local stability ofE* by performing a linearization and analyzing the corresponding characteristic equation. Similarly as in [22], let 0<μ1<μ2<μ3<… be the eigenvalues of -Δ on Ω with the homogeneous Neumann boundary condition. Let E(μi) be the eigenspace corresponding to μi in C1(Ω¯), for i=1,2,3,…. Let
(16)X={u=(u1,u2,u3)∈[C1(Ω¯)]3∂ηu=0,x∈∂Ω},{φij,j=1,…,dimE(μi)} be an orthonormal basis of E(μi), and Xij={c·φij∣c∈R3}. Then
(17)Xi=⨁j=1dimE(μi)Xij,X=⨁i=1∞Xi.
Let D=diag(D1,D2,D3) and v=(v1,v2,v3),vi=ui-ui*,(i=1,2,3). Then the linearization of (1) is vt=Lv=DΔv+Fv(E*)v,Fv(E*)v={cij}, where cii=biivi,i=1,2,3,c12=b12v2,c23=b23v3,c21=b21∫Ω‍∫-∞t‍v1(s,y)K1(x,y,t-s)dsdy,andc32=b32∫Ω‍∫-∞t‍v2(s,y)K2(x,y,t-s)dsdy. The coefficients bij are defined as follows:
(18)b11=-a11u1*+a12u1*u2*(m1+u1*)2,b12=-a12u1*m1+u1*,b13=0,b21=a21m1u2*(m1+u1*)2,b22=-a22u2*+a23u2*u3*(m2+u2*)2,b23=-a23u2*m2+u2*,b31=0,b32=a32m2u3*(m2+u2*)2,b33=-a33u3*.
Since Xi is invariant under the operator L for each i≥1, then the operator L on Xi is vt=Lv=Dμiv+Fv(E*)v. Let vi=cieλt(i=1,2,3) and we can get the characteristic equation φi(λ)=λ3+Aiλ2+Biλ+Ci=0, where Ai, Bi, and Ci are defined as follows:
(19)Ai=(d1+d2+d3)μi-b11-b22-b33,Bi=(d1d2+d2d3+d3d1)μi2111+[+d2(-b11-b33)]d3(-b11-b22)+d1(-b33-b22)111111+d2(-b11-b33)]μi111+b11b22+b11b33+b22b33111-b12b21∫0∞‍k1(s)e-λsds-b23b32∫0∞‍k2(s)e-λsds,Ci=d1d2d3μi3+[-b33d1d2-b11d2d3-b22d1d3]μi2111+[∫0∞‍b11b22d3+b11b33d2+b22b33d1+b22b33d1111111-d3b12b21∫0∞‍k1(s)e-λsds111111-d1b23b32∫0∞‍k2(s)e-λsds]μi111+[-b11b22b33+b12b21b33∫0∞‍k1(s)e-λsds111111+b11b23b32∫0∞‍k2(s)e-λsds].
It is easy to see that b12,b23,b33<0 and b21,b32>0. It follows from assumption G1 and G2 that a11<0,a22<0. So Ai>0,Bi>0,Ci>0 and AiBi-Ci>0 for i≥1 from the direct calculation. According to the Routh-Hurwitz criterion, the three roots λi,1,λi,2,λi,3 of φi(λ)=0 all have negative real parts.
By continuity of the roots with respect toμi and Routh-Hurwitz criterion, we can conclude that there exists a positive constant ε such that
(20)Re{λi,1},Re{λi,2},Re{λi,3}≤-ε,i≥1.
Consequently, the spectrum of L, consisting only of eigenvalues, lies in {Reλ≤-ε}. It is easy to see that E* is locally asymptotically stable and follows from Theorem 5.1.1 of [23].Theorem 8.
Assume that(21)H3:a1(a21-a2)m1a2>a11>a12u2*m12,H2 and G2 hold, and the positive steady state E* of system (1)–(3) with nontrivial initial function is globally asymptotically stable.Proof.
It is easy to see that the equations in (1) can be rewritten as
(22)∂u1∂t-d1Δu1=u1[a12(u2-u2*)m1+u1-a11(u1-u1*)2000000000000311-a12(u2-u2*)m1+u1+a12u2*(u1-u1*)(m1+u1)(m1+u1*)],∂u2∂t-d2Δu2=u2[+a23u3*(u2-u2*)(m2+u2)(m2+u2*)]∫Ω‍∫-∞t‍K1(x,y,t-s)2000000000000311×a21m1(u1(s,y)-u1*)(m1+u1(s,y))(m1+u1*)dsdy2000000000000311-a22(u2-u2*)-a23(u3-u3*)m2+u22000000000000311+a23u3*(u2-u2*)(m2+u2)(m2+u2*)],∂u3∂t-d3Δu3=u3[∫Ω‍∫-∞t‍K2(x,y,t-s)2000000000000311×a32m2(u2(s,y)-u2*)(m2+u2(s,y))(m2+u2*)dsdy2000000000000311-a33(u3-u3*)[∫Ω‍∫-∞t‍K2(x,y,t-s)].
Define
(23)V(t)=α∫Ω‍(u1-u1*-u1*lnu1u1*)dx+∫Ω‍(u2-u2*-u2*lnu2u2*)dx+β∫Ω‍(u2-u3*-u3*lnu3u3*)dx,
where α and β are positive constants to be determined. Calculating the derivatives V(t) along the positive solution to the system (1)–(3) yields
(24)V′(t)=Φ1(t)+Φ2(t),
where
(25)Φ1(t)=-∫Ω‍(αd1u1*u12|∇u1|2+d2u2*u22|∇u2|2111111111111+βd3u3*u32|∇u3|2)dx,Φ2(t)=-∫Ω‍α[a11-a12u2*(m1+u1)(m1+u1*)](u1-u1*)2dx111111,-∫Ω‍[a22-a23u3*(m2+u2)(m2+u2*)](u2-u2*)2dx111111,-∫Ω‍a33β(u3-u3*)2dx111111,-∫Ω‍a23(u3-u3*)(u2-u2*)m2+u2dx111111,-∫Ω‍αa12(u2-u2*)(u1-u1*)m1+u1dx111111,+∬Ω∫-∞t‍K1(x,y,t-s)111111,×a21m1(u1(s,y)-u1*)(u2(t,x)-u2*)(m1+u1(s,y))(m1+u1*)dsdydx111111,+β∬Ω∫-∞t‍K2(x,y,t-s)111111,×a32m2(u2(s,y)-u2*)(u3(t,x)-u3*)(m2+u2(s,y))(m2+u2*)dsdydx.
Applying the inequality ab≤ϵa2+(1/4ϵ)b2, we derive from (25) that
(26)Φ2(t)≤-∫Ω‍α[a11-a12u2*(m1+u1)(m1+u1*)](u1-u1*)2dx111111,-∫Ω‍a33β(u3-u3*)2dx111111,-∫Ω‍[a22-a23u3*(m2+u2)(m2+u2*)](u2-u2*)2dx111111,+∫Ω‍αa12m1+u1[ϵ1(u1-u1*)2+14ϵ1(u2-u2*)2]dx111111,+∫Ω‍a23m2+u2[14ϵ2(u2-u2*)2+ϵ2(u3-u3*)2]dx111111,+∬Ω∫-∞t‍K1(x,y,t-s)111111,×a21m1(m1+u1(s,y))(m1+u1*)111111,×[14ϵ3ϵ3(u1(s,y)-u1*)2111111111+14ϵ3(u2(t,x)-u2*)2]dsdydx111111,+∬Ω∫-∞t‍K2(x,y,t-s)111111,×a32βm2(m2+u2(s,y))(m2+u2*)111111,×[14ϵ4(u2(s,y)-u2*)2111111111+ϵ4(u3(t,x)-u3*)214ϵ3]dsdydx.
According to the property of the Kernel functions Ki(x,y,t), (i=1,2), we know that
(27)Φ2(t)≤-∫Ω‍[αa11-αa12u2*m12-αa12ϵ1m1](u1-u1*)2dx-∫Ω‍[a22-a23u3*m22-αa124m1ϵ1-a234m2ϵ220031000-a214m1ϵ3](u2-u2*)2dx-∫Ω‍[a33β-a23ϵ2m2-a32βϵ4m2](u3-u3*)2dx+a21ϵ3m1∬Ω∫-∞t‍K1(x,y,t-s)111111111111111111×(u1(s,y)-u1*)2dsdydx+a32β4m2ϵ4∬Ω∫-∞t‍K2(x,y,t-s)111111111111111111×(u2(s,y)-u2*)2dsdydx.
Define a new Lyapunov functional as follows:
(28)E(t)=V(t)+a21ϵ3m1∬Ω∫0+∞‍∫t-rt‍K1(x,y,r)1111111111111111111×(u1(s,y)-u1*)2dsdrdydx+a32β4m2ϵ4∬Ω∫0+∞‍∫t-rt‍K2(x,y,r)1111111111111111111×(u2(s,y)-u2*)2dsdrdydx.
It is derived from (27) and (28) that
(29)E′(t)≤-∫Ω‍(αd1u1*u12|∇u1|2+d2u2*u22|∇u2|211111111+βd3u3*u32|∇u3|2)dx-[αa11-αa12u2*m12-αa12ϵ1m111111-a21ϵ3m1]∫Ω‍(u1-u1*)2dx-[a22-a23u3*m22-αa124m1ϵ1-a234m2ϵ2-a214m1ϵ3111111-a32β4m2ϵ4]∫Ω‍(u2-u2*)2dx-[a33β-a23ϵ2m2-βa32ϵ4m2]∫Ω‍(u3-u3*)2dx.
Since (u1,u2,u3) is the unique positive solution of system (1). Using Proposition 5, there exists a constant C which does not depend on x∈Ω¯ or t≥0 such that ∥ui(·,t)∥∞≤C (i=1,2,3) for t≥0. By Theorem A2 in [24], we have
(30)∥ui(·,t)∥C2,α(Ω¯)≤C,∀t≥1.
Assume that
(31)G1:a11m12>a12u2*,G2:a22m22>a23u3*+2a23a32a33+2a12a21m22a11m12-a12u2*,
and denote
(32)l1=αa11-αa12u2*m12-αa12ϵ1m1-a21ϵ3m1,l2=a22-a23u3*m22-αa124m1ϵ1-a234m2ϵ2-a214m1ϵ3-a32β4m2ϵ4,l3=a33β-a23ϵ2m2-βa32ϵ4m2.
Then (29) is transformed into
(33)E′(t)≤-∫Ω‍(αd1u1*u12|∇u1|2+d2u2*u22|∇u2|211111111+βd3u3*u32|∇u3|2)dx-∑i=13‍li∫Ω‍(ui-ui*)2dx.
If we choose
(34)α=a21a12,β=a23a32,ϵ1=ϵ3=a11m12-a12u2*4m1a12,ϵ2=ϵ4=m2a334a32,
then li>0(i=1,2,3). Therefore, we have
(35)E′(t)≤-∑i=13‍li∫Ω‍(ui-ui*)2dx.
From Proposition6 we can see that the solution of system (1) and (3) is bounded, and so are the derivatives of (ui-ui*) (i=1,2,3) by the equations in (1). Applying Lemma 1, we obtain that
(36)limt→∞∫Ω‍(ui-ui*)dx=0,(i=1,2,3).
Recomputing E′(t) gives
(37)E′(t)≤-∫Ω‍(αd1u1*u12|∇u1|2+d2u2*u22|∇u2|211111111+βd3u3*u32|∇u3|2)dx≤-c∫Ω‍(|∇u1|2+|∇u2|2+|∇u3|2)dx=-g(t),
where c=min{αd1u1*/M12,d2u2*/M22,βd3u3*/M32}. Using (30) and (1), we obtain that the derivative of g(t) is bounded in [1,+∞). From Lemma 1, we conclude that g(t)→0 as t→∞. Therefore
(38)limt→∞∫Ω‍(|∇u1|2+|∇u2|2+|∇u3|2)dx=0.
Applying the Poincaré inequality
(39)∫Ω‍λ|u-u¯|2dx≤∫Ω‍|∇u|2dx,
leads to
(40)limt→∞∫Ω‍(ui-u¯i)2dx=0,(i=1,2,3),
where u¯i=(1/|Ω|)∫Ω‍uidx and λ is the smallest positive eigenvalue of -Δ with the homogeneous Neumann condition. Therefore,
(41)|Ω|(u¯1(t)-u1*)2=∫Ω‍(u¯1(t)-u1*)2dx=∫Ω‍(u¯1(t)-u1(t,x)+u1(t,x)-u1*)2dx≤2∫Ω‍(u¯1(t)-u1(t,x))2dx+2∫Ω‍(u1(t,x)-u1*)2dx.
So we have u¯1(t)→u1* as t→∞. Similarly, u¯2(t)→u2* and u¯3(t)→u3* as t→∞. According to (30), there exists a subsequence tm, and non-negative functions wi∈C2(Ω¯), such that
(42)limm→∞∥ui(·,tm)-wi(·)∥C2(Ω¯)=0(i=1,2,3).
Applying (40) and noting that u¯i(t)→ui*, we then have wi=ui*, (i=1,2,3). That is,
(43)limm→∞∥ui(·,tm)-ui*∥C2(Ω¯)=0(i=1,2,3).
Furthermore, the local stability ofE* combining with (43) gives the following global stability.Theorem 9.
Assume that(44)H5:a1(a21-a2)m1a2>a11>a12v2*m12,H4 and G4 hold, and then the semi-trivial steady state E2 of system (1)–(3) with non-trivial initial functions is globally asymptotically stable.Proof.
It is obvious that system (1)–(3) always has two non-negative equilibria (u1,u2,u3) as follows: E0=(0,0,0) and E1=(a1/a11,0,0). If H1 and H4 are satisfied, system (1) has the other semitrivial solution denoted by E2(v1*,v2*,0), where
(45)H4:a3m2>(a32-a3)v2*.
We consider the stability ofE2 under condition H1 and H4. Equation (1) can be rewritten as
(46)∂u1∂t-d1Δu1=u1[-a11(u1-v1*)-a12(u2-v2*)m1+u11111111111111111+a12v2*(u1-v1*)(m1+u1)(m1+v1*)],∂u2∂t-d2Δu2=u2[∫Ω‍∫-∞t‍K1(x,y,t-s)11111111111111×a21m1(u1(s,y)-v1*)(m1+u1(s,y))(m1+v1*)dsdy11111111111111-a22(u2-v2*)-a23u3m2+u2[∫Ω‍∫-∞t‍K1(x,y,t-s)],∂u3∂t-d3Δu3=u3[∫Ω‍∫-∞t‍K2(x,y,t-s)11111111111111×a32m2(u2(s,y)-v2*)(m2+u2(s,y))(m2+v2*)dsdy11111111111111-a33u3[∫Ω‍∫-∞t‍K2(x,y,t-s)].
Similar to the argument of Theorem 8, we have u1(x,t)→v1* and u2(x,t)→v2* as t→∞ uniformly on Ω¯ provided that the following additional condition holds:
(47)G3:a11m12>a12v2*,G4:a22m22>a23a32a33+2a12a21m22a11m12-a12v2*.
Next, we consider the asymptotic behavior of u3(t,x). Let
(48)θ=a3-(a32-a3)v2*2δ,δ=|a32-a3|.
Then there exists t1>0 such that
(49)v2*-θ<u2(t,x)<v2*+θ,∀t>t1,x∈Ω¯.
Consider the following two systems:
(50)w3t-d3Δw3=w3[-a3+a32(v2*+θ)m2+v2*+θ-a33w3],111111111111111111111,,(t,x)∈[t1,+∞)×Ω,∂w3∂ν=0,(t,x)∈[t1,+∞)×∂Ω,w3=u3,(t,x)∈(-∞,t1]×Ω,W3t-d3ΔW3=W3[-a3+a32(v2*-θ)m2+v2*-θ-a33v3],1111111111111,,111111111(t,x)∈[t1,+∞)×Ω,∂W3∂ν=0,(t,x)∈[t1,+∞)×∂Ω,W3=u3,(t,x)∈(-∞,t1]×Ω.
Combining comparison principle with (50), we obtain that
(51)W3(t,x)≤u3(t,x)≤w3(t,x),∀t>t1,x∈Ω¯.
By Lemma 4, we obtain
(52)0=limt→∞W3(t,x)≤infu3(t,x)≤supu3(t,x)≤limt→∞w3(t,x)=0,∀x∈Ω¯,
which implies that limt→∞u3(t,x)=0 uniformly on Ω¯.Theorem 10.
Suppose thatG5 and G6 hold, and then the semi-trivial steady state E1 of system (1)–(3) with non-trivial initial functions is globally asymptotically stable.Proof.
We study the stability of the semi-trivial solutionE1=(u~1,0,0). Similarly, the equations in (1) can be written as
(53)∂u1∂t-d1Δu1=u1[-a11(u1-u~1)-a12u2m1+u1],∂u2∂t-d2Δu2=u2[-a2+∫Ω‍∫-∞t‍K1(x,y,t-s)1111111111111111111×a21m1(u1(s,y)-u~1)(m1+u1(s,y))(m1+u~1)dsdy1111111111111111111+a21u~1m1+u~1-a22u2-a23u3m2+u2],∂u3∂t-d3Δu3=u3[-a3+∫Ω‍∫-∞t‍K2(x,y,t-s)1111111111111111111×a32u2(s,y)m2+u2(s,y)dsdy-a33u3].
Define
(54)V(t)=α∫Ω‍[u1-u~1-u~1logu1u~1]dx+∫Ω‍u2dx+β∫Ω‍u3dx.
Calculating the derivative of V(t) along E1, we get from (54) that
(55)V′(t)≤-∫Ω‍[α(a11-a12ϵ1m1)](u1-u~1)2dx+∫Ω‍[a22-αa124m1ϵ1-a234ϵ2m2-a214m1ϵ3]u22dx+∫Ω‍[βa33-a23ϵ2m2]u32dx-∫Ω‍[(a2-a21u~1m1+u~1)u2+a3u3]dx+a21ϵ3m1∫Ω‍∫-∞t‍K1(x,y,t-s)11111111111111×(u1(s,y)-u~1)2dsdydx+a32β4ϵ4m2∬Ω∫-∞t‍K2(x,y,t-s)111111111111111111×u22(s,y)dsdydx.
Define
(56)E(t)=V(t)+a21ϵ3m1∬Ω∫0+∞‍∫t-rt‍K1(x,y,r)11111111111111111111×(u1(s,y)-u1*)2dsdrdydx+a32β4ϵ4m2∬Ω∫0+∞‍∫t-rt‍K2(x,y,r)11111111111111111111×(u2(s,y)-u2*)2dsdrdydx.
It is easy to see that
(57)E′(t)≤-αd1u~1∫Ω‍|∇u1|2u12dx-∫Ω‍[α(a11-a12ϵ1m1)-a21ϵ3m1](u1-u~1)2dx+∫Ω‍[a22-αa124m1ϵ1-a234m2ϵ2-a214m1ϵ311111111-a324m2ϵ4]u22dx+∫Ω‍[βa33-a23ϵ2m2-a32βϵ4m2]u32dx-∫Ω‍[(a2-a21u~1m1+u~1)u2-a3βu3]dx.
Assume that
(58)G5:a1(a21-a2)<m1a2a11,G6:a22m22>a23a32a33+2a12a21m22a11m12.
Let
(59)l1=αa11-αa12ϵ1m1-a21ϵ3m1,l2=a22-αa124m1ϵ1-a234m2ϵ2-a214m1ϵ3-a32β4m2ϵ4,l3=a33β-a23ϵ2m2-βa32ϵ4m2.
Choose
(60)α=a21a12,β=a23a32,ϵ1=ϵ3=a11m14a12,ϵ2=ϵ4=m2a334a32.
Then we get li>0 (i=1,2,3). Therefore we have
(61)limt→∞u1(t,x)=u~1,
uniformly on Ω¯.Next, we consider the asymptotic behavior ofu2(t,x) and u3(t,x). For any T>0, integrating (57) over [0,T] yields
(62)E(T)+αd1∥|∇u1u1|∥22+l1∥|u1-u~1|∥22+l2∥|u2|∥22+l3∥|u3|∥22≤E(0),
where ∥|ui|∥22=∫0T‍∫Ω‍ui2dxdt. It implies that ∥|ui|∥2≤Ci(i=2,3) for the constant Ci which is independent of T. Now we consider the boundedness of ∥|∇u2|∥2 and ∥|∇u3|∥2. From the Green's identity, we obtain
(63)d2∫0T‍∫Ω‍|∇u2(t,x)|2dxdt=-d2∫0T‍∫Ω‍u2(t,x)Δu2(t,x)dxdt=-∫0T‍∫Ω‍u2∂u2∂tdxdt-a2∫0T‍∫Ω‍u22dxdt-a22∫0T‍∫Ω‍u23dxdt-a23∫0T‍∫Ω‍u22u3m2+u3dxdt+a21∫0T‍∬Ω∫-∞t‍u1(s,y)u22(t,x)m1+u1(s,y)1111111111111111111×K1(x,y,t-s)dsdydxdt.
Note that
(64)∫0T‍∫Ω‍u2(t,x)∂u2(t,x)∂tdxdt=∫0T‍∫Ω‍12∂u22∂tdtdx=12∫Ω‍u22(T,x)dx-12∫Ω‍u22(0,x)dx≤M22|Ω|,∫0T‍∫Ω‍u22u3dxdt≤M3∫0T‍∫Ω‍u22dxdt,∫0T‍∫Ω‍u23dxdt≤M2∫0T‍∫Ω‍u22dxdt,∫0T‍∬Ω∫-∞t‍u1(s,y)u22(t,x)m1+u1(s,y)K1(x,y,t-s)dsdydxdt≤∫0T‍∫Ω‍u22(t,x)dxdt.
Thus, we get ∥|∇u2|∥2≤C4. In a similar way, we have ∥|∇u3|∥2≤C5. Here C4 and C5 are independent of T.It is easy to see thatu2(t,x),u3(t,x)∈L2((0,∞);W1,2(Ω)). These imply that
(65)limt→∞∥ui(·,t)∥W1,2(Ω)=0,i=2,3.
From the Sobolev compact embedding theorem, we know
(66)limt→∞∥ui(·,t)∥C(Ω¯)=0,i=2,3.In the end, we show that the trivial solutionE0 is an unstable equilibrium. Similarly to the local stability to E*, we can get the characteristic equation of E0 as
(67)(λ+μiD1-a1)(λ+μiD2+a2)(λ+μiD3+a3)=0.
If i=1, then μ1=0. It is easy to see that this equation admits a positive solution λ=a1. According to Theorem 5.1 in [23], we have the following result.Theorem 11.
The trivial equilibriumE0 is an unstable equilibrium of system (1)–(3).


## 4. Numerical Illustrations

In this section, we perform numerical simulations to illustrate the theoretical results given in Section3.In the following, we always takeΩ=[0,π],Ki(x,y,t)=Gi(x,y,t)ki(t), where ki(t)=(1/τ*)e-t/τ*(i=1,2) and
(68)G1(x,y,t)=1π+2π∑n=1∞‍e-d2n2tcosnxcosny,G2(x,y,t)=1π+2π∑n=1∞‍e-d3n2tcosnxcosny.
However, it is difficult for us to simulate our results directly because of the nonlocal term. Similar to [25], the equations in (1) can be rewritten as follows:
(69)∂u1∂t-d1Δu1=u1(a1-a11u1-a12u2m1+u1),∂u2∂t-d2Δu2=u2(-a2+a21v1-a22u2-a23u3m2+u3),∂u3∂t-d3Δu3=u3(-a3+a32v2-a33u3),∂v1∂t-d2Δv1=1τ*(u1m1+u1-v1),∂v2∂t-d3Δv2=1τ*(u2m2+u2-v2),
where
(70)vi=∫0π‍∫-∞t‍Gi(x,y,t-s)1τ*e-(t-s)/τ*ui(s,y)mi+ui(s,y)dyds,11111111i=1,2.
Each component is considered with homogeneous Neumann boundary conditions, and the initial condition of vi is
(71)vi(0,x)=∫0π‍∫-∞0‍Gi(x,y,-s)1τ*es/τ*ui(s,y)mi+ui(s,y)dyds,111111111i=1,2.
In the following examples, we fix some coefficients and assume that d1=d2=d3=1,a1=3,a12=1,m1=m2=1,a2=1,a22=12,a23=1,a33=12, and τ*=1. The asymptotic behaviors of system (1)–(3) are shown by choosing different coefficients a11,a21,a3, and a32.Example 12.
Leta11=53/9,a21=81/13,a3=1/13 and a32=14. Then it is easy to see that the system admits a unique positive equilibrium E*(1/2,1/12,1/12). By Theorem 8, we see that the positive solution (u1(t,x),u2(t,x),u3(t,x)) of system (1)–(3) converges to E* as t→∞. See Figure 2.Global stability of the positive equilibriumE* with an initial condition (ϕ1,ϕ2,ϕ3)=(2.7+0.5sinx,2.7+0.5cosx,2.7+0.5sinx).
(a)
(b)
(c)Example 13.
Leta11=6,a21=12,a3=1/4 and a32=1. Clearly, H2 does not hold. Hence, the positive steady state is not feasible. System (1) admits two semi-trivial steady state E2(0.4732,0.2372,0) and E1(1/2,0,0). According to Theorem 9, we know that the positive solution (u1(t,x),u2(t,x),u3(t,x)) of system (1)–(3) converges to E2 as t→∞. See Figure 3.Global stability of the positive equilibriumE2 with an initial condition (ϕ1,ϕ2,ϕ3)=(2.7+0.5sinx,2.7+0.5cosx,2.7+0.5sinx).
(a)
(b)
(c)Example 14.
Leta11=6,a21=2,a3=1/4 and a32=1. Clearly, H1 does not hold. Hence, E* and E2 are not feasible. System (1) has a unique semi-trivial steady state E1(1/2,0,0). According to Theorem 10 we know that the positive solution (u1(t,x),u2(t,x),u3(t,x)) of system (1)–(3) converges to E1 as t→∞. See Figure 4.Global stability of the positive equilibriumE1 with an initial condition (ϕ1,ϕ2,ϕ3)=(2.7+0.5sinx,2.7+0.5cosx,2.7+0.5sinx).
(a)
(b)
(c)


## 5. Discussion

In this paper, we incorporate nonlocal delay into a three-species food chain model with Michaelis-Menten functional response to represent a delay due to the gestation of the predator. The conditions, under which the spatial homogeneous equilibria are asymptotically stable, are given by using the Lyapunov functional.We now summarize the ecological meanings of our theoretical results. Firstly, the positive equilibriumE* of system (1)–(3) exists under the high birth rate of the prey (a1) and low death rates (a2 and a3) of predator and top predator. E* is globally stable if the intraspecific competition a11 is neither too big nor too small and the maximum harvest rates a12, a23 are small enough. Secondly, the semi-trivial equilibrium E2 of system (1)–(3) exists if the birth rate of the prey (a1) is high, death rate of predator a2 is low, and the death rate (a3) exceeds the conversion rate from predator to top predator (a32). E2 is globally stable if the maximum harvest rates a12, a23 are small and the intra-specific competition a11 is neither too big nor too small. Thirdly, system (1)–(3) has only one semi-trivial equilibrium E1 when the death rate (a2) exceeds the conversion rate from prey to predator (a21). E1 is globally stable if intra-specific competitions (a11, a22, and a33) are strong. Finally, E0 is unstable and the non-stability of trivial equilibrium tells us that not all of the populations go to extinction. Furthermore, our main results imply that the nonlocal delay is harmless for stabilities of all non-negative steady states of system (1)–(3).There are still many interesting and challenging problems with respect to system (1)–(3), for example, the permanence and stability of periodic solution or almost periodic solution. These problems are clearly worthy for further investigations.



---

*Source: 936952-2013-09-05.xml*