Global Asymptotic Stability of 3-Species Mutualism Models with Diffusion and Delay Effects

Abstract

In this paper, the Lotka-Volterra 3-species mutualism models with diffusion and delay effects is investigated. A simple and easily verifiable condition is given to ensure the global asymptotic stability of the unique positive steady-state solution of the corresponding steady-state problem in a bounded domain with Neumann boundary condition. Our approach to the problem is based on inequality skill and the method of the upper and lower solutions for a more general reaction—diffusion system. Finally, some numerical simulations are given to illustrate our results.

1. Introduction

Modeling and analysis of the dynamics of biological populations by means of differential equations is one of the primary concerns in population growth problems. A well-known and extensively studied class of models in population dynamics is the Lotka-Volterra system which models the interaction among various species. In the earlier literature most of the discussions are devoted to coupled systems of two equations (cf. [1–11]). In the recent years, attention has been given to reaction-diffusion systems with three population species, and the main concerns in these works are the prey-predator and competition models with or without time delays (cf. [12–24]). The purpose of this paper is to investigate the asymptotic behavior of the time-dependent solution for a 3-species mutualism model, where the effect of diffusion and time delays is both taken into consideration by obtaining existence of positive solution for the corresponding inequalities. The system of equations under consideration is the Lotka-Volterra 3-species mutualism model, where the population densities do not move across the boundary of a habitat, and time delays may appear in the opposing species; the densities of three populations 𝑢,𝑣,and 𝑤 are governed by the following coupled equations:

𝜕𝑢𝜕𝑡−𝑑₁𝐿₁𝑢=𝑢(𝑎₁−𝑐₁𝑢+𝑒₁𝑣+𝛽₁𝐽₂∗𝑣),𝜕𝑣𝜕𝑡−𝑑₂𝐿₂𝑣=𝑣(𝑎₂+𝑏₂𝑢−𝑐₂𝑣+𝑒₂𝑤+𝛼₂𝐽₁∗𝑢+𝛽₂𝐽₃∗𝑤),𝜕𝑤𝜕𝑡−𝑑₃𝐿₃𝑤=𝑤(𝑎₃+𝑒₃𝑣−𝑐₃𝑤+𝛽₃𝐽₂∗𝑣)(𝑥∈Ω,𝑡>0), (1.1) with the boundary and initial conditions

𝜕𝑢(𝑥,𝑡)=𝜕𝑛𝜕𝑣(𝑥,𝑡)=𝜕𝑛𝜕𝑤(𝑥,𝑡)𝜕𝑛=0(𝑥∈𝜕Ω,𝑡>0),𝑢(𝑥,𝑡)=𝜂₁(𝑥,𝑡)(𝑡∈𝐼₁),𝑣(𝑥,𝑡)=𝜂₂(𝑥,𝑡)(𝑡∈𝐼₂),𝑤(𝑥,𝑡)=𝜂₃(𝑥,𝑡)(𝑡∈𝐼₃,𝑥∈Ω). (1.2) Here Ω is a bounded domain in 𝑅^𝑛 with boundary 𝜕Ω, 𝜕/𝜕𝑛 denotes the outward normal derivative on 𝜕Ω, and 𝜂_𝑖(𝑥,𝑡) is a smooth function. The above problem (1.1)-(1.2) arises in a simple food chain describing three interacting species in a spatial habitat Ω. For each 𝑖=1,2,3, 𝑎_𝑖,𝑐_𝑖 are positive constants. 𝑏_𝑖,𝑒_𝑖,𝛼_𝑖, and 𝛽_𝑖 are nonnegative constants with 𝑒^∗_𝑖≡𝑒_𝑖+𝛽_𝑖>0, 𝑏^∗₂≡𝑏₂+𝛼₂>0, and its respective diffusion rate is denoted by 𝑑_𝑖. The real number 𝑎_𝑖 is the net birth rate of the 𝑖th species, and 𝑐_𝑖 is its respective intraspecific competition. The parameters 𝑏_𝑖,𝑒_𝑖,𝛼_𝑖 and 𝛽_𝑖 are interspecific cooperation. 𝐿_𝑖 is a uniformly elliptic operator in the form

𝐿_𝑖𝑢=−^𝑁_𝑗,𝑘=1𝑎^(𝑖)_𝑗,𝑘𝜕(𝑥,𝑡)²𝑢𝜕𝑥_𝑗𝜕𝑥_𝑘+^𝑁_𝑗=1𝑏_𝑗^(𝑖)(𝑥,𝑡)𝜕𝑢𝜕𝑥_𝑗(𝑖=1,2,3) (1.3) (cf. [1]). The functions 𝐽_𝑖∗𝑢_𝑖,𝑖=1,2,3, with (𝑢₁,𝑢₂,𝑢₃)=(𝑢,𝑣,𝑤) are given either by the discrete time delay

𝐽_𝑖∗𝑢_𝑖=𝑢_𝑖(𝑡−𝜏_𝑖,𝑥), (1.4) or by the continuous time delay 𝐽_𝑖∗𝑢_𝑖=^𝑡_−𝜏𝑖𝐽_𝑖(𝑡−𝑠)𝑢_𝑖(𝑠,𝑥)𝑑𝑠, (1.5) and the interval 𝐼_𝑖 is given by 𝐼_𝑖=[−𝜏_𝑖,0], where 𝜏_𝑖>0 is a constant representing the time delay. It is allowed that the type of time delays and the values of 𝜏_𝑖 may be different for different 𝑢_𝑖. This consideration includes various combination of discrete and continuous time delays for the species 𝑢,𝑣, and 𝑤.

Throughout the paper we assumed that the function 𝐽_𝑖(𝑡) is piecewise continuous in 𝑅₊≡[0,+∞) and possesses the property

𝐽_𝑖(𝑡)≥0for𝑡≥0,𝐽_𝑖(𝑡)≡0for𝑡>𝜏_𝑖,_𝜏𝑖0𝐽_𝑖(𝑡)𝑑𝑡=1(𝑖=1,2,3). (1.6) The above property implies that for any constant function 𝑢_𝑖(𝑡,𝑥)=𝜌_𝑖, we have

^𝑡_−𝜏𝑖𝐽_𝑖(𝑡−𝑠)𝑢_𝑖(𝑠,𝑥)𝑑𝑠=𝜌_𝑖(𝑡≥0). (1.7) The same is obviously true if 𝐽_𝑖∗𝑢_𝑖 is given by (1.4). It is also assumed that the domain Ω is smooth and the coefficients 𝑎^(𝑖)_𝑗,𝑘,𝑏_𝑗^(𝑖) of 𝐿_𝑖 are smooth functions in Ω (cf. [12]). In the special case of the diffusion operator 𝐿_𝑖𝑢_𝑖=𝐷_𝑖(𝑥)𝐷²𝑢_𝑖, it suffices to assume that 𝐷_𝑖(𝑥) is strictly positive on Ω and in 𝐶^1+𝛼(Ω) for some 𝛼∈(0,1). We allow 𝐿_𝑖=0 (and without the corresponding boundary condition) for some or all 𝑖. In particular, if 𝐿_𝑖𝑢_𝑖≡0 for all 𝑖 then the equations in (1.1)-(1.2) are reduced to the ordinary differential system (with time delays)

𝑑𝑢𝑑𝑡=𝑢(𝑎₁−𝑐₁𝑢+𝑒₁𝑣+𝛽₁𝐽₂∗𝑣)(𝑡>0),𝑑𝑣𝑑𝑡=𝑣(𝑎₂+𝑏₂𝑢−𝑐₂𝑣+𝑒₂𝑤+𝛼₂𝐽₁∗𝑢+𝛽₂𝐽₃∗𝑤)(𝑡>0),𝑑𝑤𝑑𝑡=𝑤(𝑎₃+𝑒₃𝑣−𝑐₃𝑤+𝛽₃𝐽₂∗𝑣)(𝑡>0),𝑢(𝑡)=𝜂₁(𝑡)(𝑡∈𝐼₁),𝑣(𝑡)=𝜂₂(𝑡)(𝑡∈𝐼₂),𝑤(𝑥,𝑡)=𝜂₃(𝑡)(𝑡∈𝐼₃). (1.8) Problem (1.8) and various similar problems have been investigated by many investigators in the framework of ordinary differential systems (cf. [25–28] and references therein). It is to be noted that if Δ𝑢_𝑖≠0 and 𝛼_𝑖=𝛽_𝑖=0 for every 𝑖 then problem (1.1) is reduced to the following 3-species mutualism models without time delays: 𝑢_1𝑡−𝑑₁Δ𝑢₁=𝑢₁(𝑎₁−𝑐₁𝑢₁+𝑒₁𝑢₂),inΩ_𝑇,𝑢_2𝑡−𝑑₂Δ𝑢₂=𝑢₂(𝑎₂+𝑏₂𝑢₁−𝑐₂𝑢₂+𝑒₂𝑢₃),inΩ_𝑇,𝑢_3𝑡−𝑑₃Δ𝑢₃=𝑢₃(𝑎₃+𝑏₃𝑢₂−𝑐₃𝑢₃),inΩ_𝑇,𝑢₁(𝑥,𝑡)=𝑢₂(𝑥,𝑡)=𝑢₃(𝑥,𝑡)=0,on𝜕Ω_𝑇,𝑢_𝑖(𝑥,0)=𝑢_𝑖,0(𝑥),for𝑖=1,2,3,inΩ, (1.9) where the domains Ω_𝑇 and 𝜕Ω_𝑇 are defined as, respectively, Ω×(0,𝑇) and 𝜕Ω×(0,𝑇). Here Ω is a bounded domain in 𝑅^𝑁 with smooth boundary 𝜕Ω. 𝑇 is the maximal existence time of the solution. 𝑢_𝑖,0 is a smooth function satisfying the compatibility condition 𝑢_𝑖,0(𝑥)=0 for 𝑥∈𝜕Ω. The meaning of the parameters of systems (1.9) is as same as those of (1.1)-(1.2). Kim and Lin [22] proved that global solutions of (1.9) exist if the intraspecific competitions are strong; whereas blow-up solutions exist under certain conditions of that the intraspecific competitions are weak, and [23] obtained the upper bound of blow-up rate for any 𝑁 and the lower bound of blow-up rate for 𝑁=1. On the other hand, if 𝑒₂=𝛽₂=0 then it reduces to the 2-species cooperating models, which have been extensively investigated in the current literature (cf. [2, 10, 11] and references therein). Our conclusion about the global asymptotic stability of a positive steady-state solution to system (1.1)-(1.2) is directly applicable to the above special cases (see Corollary 3.2).

This work is motivated from the following two-prey one-predator model: 𝜕𝑢𝜕𝑡−𝐿₁𝑢=𝑎₁(𝑥)𝑢(1−𝑢−𝑏₁𝑣−𝑐₁𝑤−𝛽₁𝐽₂∗𝑣−𝛾₁𝐽₃∗𝑤),𝜕𝑣𝜕𝑡−𝐿₂𝑣=𝑎₂(𝑥)𝑣(1−𝑣−𝑏₂𝑢−𝑐₂𝑤−𝛽₂𝐽₁∗𝑢−𝛾₂𝐽₃∗𝑤),𝜕𝑤𝜕𝑡−𝐿₃𝑤=𝑎₃(𝑥)𝑤(1−𝑤+𝑏₃𝑢+𝑐₃𝑣+𝛽₃𝐽₁∗𝑢+𝛾₃𝐽₂∗𝑣)(𝑥∈Ω,𝑡>0),𝜕𝑢=𝜕𝑛𝜕𝑣(𝑥,𝑡)=𝜕𝑛𝜕𝑤(𝑥,𝑡)𝜕𝑛=0(𝑥∈𝜕Ω,𝑡>0),𝑢(𝑥,𝑡)=𝜂₁(𝑥,𝑡)(𝑡∈𝐼₁),𝑣(𝑥,𝑡)=𝜂₂(𝑥,𝑡)(𝑡∈𝐼₂),𝑤(𝑥,𝑡)=𝜂₃(𝑥,𝑡)(𝑡∈𝐼₃,𝑥∈Ω), (1.10) and one-prey two-predator model 𝜕𝑢𝜕𝑡−𝐿₁𝑢=𝑎₁(𝑥)𝑢(1−𝑢−𝑏₁𝑣−𝑐₁𝑤−𝛽₁𝐽₂∗𝑣−𝛾₁𝐽₃∗𝑤),𝜕𝑣𝜕𝑡−𝐿₂𝑣=𝑎₂(𝑥)𝑣(1−𝑣+𝑏₂𝑢−𝑐₂𝑤+𝛽₂𝐽₁∗𝑢−𝛾₂𝐽₃∗𝑤),𝜕𝑤𝜕𝑡−𝐿₃𝑤=𝑎₃(𝑥)𝑤(1−𝑤+𝑏₃𝑢−𝑐₃𝑣+𝛽₃𝐽₁∗𝑢−𝛾₃𝐽₂∗𝑣)(𝑥∈Ω,𝑡>0),𝜕𝑢=𝜕𝑛𝜕𝑣(𝑥,𝑡)=𝜕𝑛𝜕𝑤(𝑥,𝑡)𝜕𝑛=0(𝑥∈𝜕Ω,𝑡>0),𝑢(𝑥,𝑡)=𝜂₁(𝑥,𝑡)(𝑡∈𝐼₁),𝑣(𝑥,𝑡)=𝜂₂(𝑥,𝑡)(𝑡∈𝐼₂),𝑤(𝑥,𝑡)=𝜂₃(𝑥,𝑡)(𝑡∈𝐼₃,𝑥∈Ω) (1.11) as well as three-species food-chain model

𝜕𝑢𝜕𝑡−𝐿₁𝑢=𝑎₁(𝑥)𝑢(1−𝑢−𝑏₁𝑣−𝑐₁𝑤−𝛽₁𝐽₂∗𝑣−𝛾₁𝐽₃∗𝑤),𝜕𝑣𝜕𝑡−𝐿₂𝑣=𝑎₂(𝑥)𝑣(1−𝑣+𝑏₂𝑢−𝑐₂𝑤+𝛽₂𝐽₁∗𝑢−𝛾₂𝐽₃∗𝑤),𝜕𝑤𝜕𝑡−𝐿₃𝑤=𝑎₃(𝑥)𝑤(1−𝑤+𝑏₃𝑢+𝑐₃𝑣+𝛽₃𝐽₁∗𝑢+𝛾₃𝐽₂∗𝑣)(𝑥∈Ω,𝑡>0),𝜕𝑢=𝜕𝑛𝜕𝑣(𝑥,𝑡)=𝜕𝑛𝜕𝑤(𝑥,𝑡)𝜕𝑛=0(𝑥∈𝜕Ω,𝑡>0),𝑢(𝑥,𝑡)=𝜂₁(𝑥,𝑡)(𝑡∈𝐼₁),𝑣(𝑥,𝑡)=𝜂₂(𝑥,𝑡)(𝑡∈𝐼₂),𝑤(𝑥,𝑡)=𝜂₃(𝑥,𝑡)(𝑡∈𝐼₃,𝑥∈Ω). (1.12) Pao [29] obtained some simple and easily verifiable conditions for the existence and global asymptotic stability of a positive steady-steady solution for each of the above three-model problems.

Based on the above results, we are mainly interested in studying the asymptotic behavior of the solution of (1.1)-(1.2). The corresponding steady-state problem is given by

−𝑑₁𝐿₁𝑢=𝑢(𝑎₁−𝑐₁𝑢+𝑒₁𝑣+𝛽₁𝑣)(𝑥∈Ω),−𝑑₂𝐿₂𝑣=𝑣(𝑎₂+𝑏₂𝑢−𝑐₂𝑣+𝑒₂𝑤+𝛼₂𝑢+𝛽₂𝑤)(𝑥∈Ω),−𝑑₃𝐿₃𝑤=𝑤(𝑎₃+𝑒₃𝑣−𝑐₃𝑤+𝛽₃𝑣)(𝑥∈Ω),𝜕𝑢=𝜕𝑛𝜕𝑣=𝜕𝑛𝜕𝑤𝜕𝑛=0(𝑥∈𝜕Ω). (1.13) It is clear that the above system has always the trivial solution (0, 0, 0) and various forms of constant semitrivial solutions, that is, constant nonnegative solutions with at least one component zero and one component positive. Our aim is to prove the existence and global asymptotic stability of a positive constant steady-state solution in the mutualism model (with respect to nonnegative initial perturbations). This global asymptotic stability result implies that system (1.1)-(1.2) is permanent, the trivial and all forms of semitrivial solutions are unstable, and the nonuniform steady-state solution does not exist.

This paper is arranged as follows. In Section 2, we give some preliminary results for a more general time-delayed reaction-diffusion system. The main results and their proofs are given in Section 3. Some numerical simulations are shown in Section 4 to illustrate our theoretical analysis.

2. Some Preliminary Results

To prove the main results in this paper we first give some preliminary results for a more general time-delayed parabolic system in the form

𝜕𝑢_𝑖𝜕𝑡−𝐿_𝑖𝑢_𝑖=𝑓_𝑖(𝑥,𝑈,𝐽∗𝑈)(𝑥∈Ω,𝑡>0),𝜕𝑢_𝑖𝑢𝜕𝑛=0(𝑥∈𝜕Ω,𝑡>0),_𝑖(𝑥,𝑡)=𝜂_𝑖(𝑥,𝑡)(𝑥∈Ω,𝑡∈𝐼_𝑖),𝑖=1,2,…,𝑁, (2.1) where 𝑈=(𝑢₁,𝑢₂,…,𝑢_𝑁),𝐽∗𝑈=(𝐽₁∗𝑢₁,𝐽₂∗𝑢₂,…,𝐽_𝑁∗𝑢_𝑁), and for each 𝑖=1,2,…,𝑁,𝑓_𝑖(⋅,𝑈,𝐽∗𝑈) is, in general, a nonlinear function of 𝑈 and 𝐽∗𝑈, and 𝐿_𝑖 is a uniformly elliptic operator in the form of (1.3). The components 𝐽_𝑖∗𝑢_𝑖 of 𝐽∗𝑈 are given by (1.4) for 𝑖=1,2,…,𝑛₀ and by (1.5) for 𝑖=𝑛₀+1,…,𝑁, where 𝑛₀ is a nonnegative integer. (The case 𝑛₀=0 corresponds to (1.5) for every 𝑖=1,2,…,𝑁.)

By writing 𝑈 and 𝑉 in the split forms 𝑈≡(𝑢_𝑖,[𝑢]_𝑎𝑖,[𝑢]_𝑏𝑖) and 𝑉≡([𝑣]_𝑐𝑖,[𝑣]_𝑑𝑖), respectively, where 𝑎_𝑖,𝑏_𝑖,𝑐_𝑖, and 𝑑_𝑖 are nonnegative integers satisfying 𝑎_𝑖+𝑏_𝑖=𝑁−1,𝑐_𝑖+𝑑_𝑖=𝑁(𝑖=1,2,…,𝑁), (2.2) and [𝑤]_𝛿𝑖 denotes a vector with 𝛿_𝑖 components of 𝑤, we write 𝑓_𝑖(𝑥,𝑈,𝑉) in the form

𝑓_𝑖(𝑥,𝑈,𝑉)=𝑓_𝑖(𝑥,𝑢_𝑖,[𝑢]_𝑎𝑖,[𝑢]_𝑏𝑖,[𝑣]_𝑐𝑖,[𝑣]_𝑑𝑖)forevery𝑖=1,2,…,𝑁. (2.3)

Definition 2.1. We say that the vector function 𝐹(⋅,𝑈,𝑉)≡(𝑓₁(⋅,𝑈,𝑉),𝑓₂(⋅,𝑈,𝑉),…,𝑓_𝑁(⋅,𝑈,𝑉)) possesses a mixed quasimonotone property in a subset 𝑆 of 𝑅^𝑁, if for each 𝑖 there exist nonnegative integers 𝑎_𝑖,𝑏_𝑖,𝑐_𝑖, and 𝑑_𝑖 satisfying (2.2) such that 𝑓_𝑖(⋅,𝑢_𝑖,[𝑢]_𝑎𝑖,[𝑢]_𝑏𝑖,[𝑣]_𝑐𝑖,[𝑣]_𝑑𝑖) is nondecreasing in [𝑢]_𝑎𝑖,[𝑣]_𝑐𝑖 and is nonincreasing in [𝑢]_𝑏𝑖,[𝑣]_𝑑𝑖 for every 𝑈,𝑉 in 𝑆 (cf. [1, 12]). In particular, if 𝑏_𝑖=𝑑_𝑖=0 for all 𝑖 then 𝐹(⋅,𝑈,𝑉) is said to be quasimonotone nondecreasing in 𝑆. It is easily seen that the reaction functions in the models (1.1) are quasimonotone nondecreasing in 𝑅^𝑁₊.

Definition 2.2. For quasimonotone nondecreasing functions 𝐹(⋅,𝑈,𝑉) we call a pair of smooth functions 𝑈≡(𝑢₁,𝑢₂,…,𝑢_𝑁), 𝑈≡(𝑢₁,𝑢₂,…,𝑢_𝑁) coupled upper and lower solutions of (2.1) if 𝑈>𝑈 and if 𝜕𝑢_𝑖𝜕𝑡−𝐿_𝑖𝑢_𝑖≥𝑓_𝑖𝜕𝑢(𝑥,𝑈,𝐽∗𝑈)(𝑥∈Ω,𝑡>0),_𝑖𝜕𝑡−𝐿_𝑖𝑢_𝑖≤𝑓_𝑖𝜕𝑢(𝑥,𝑈,𝐽∗𝑈)(𝑥∈Ω,𝑡>0),_𝑖𝜕𝑢𝜕𝑛≥0≥^𝑖𝑢𝜕𝑛(𝑥∈𝜕Ω,𝑡>0),_𝑖(𝑥,𝑡)≥𝜂_𝑖(𝑥,𝑡)≥𝑢_𝑖(𝑥,𝑡)(𝑥∈Ω,𝑡∈𝐼_𝑖,𝑖=1,2,…,𝑁), (2.4) where inequality between vectors is in the componentwise sense. It is obvious that if 𝑈≡𝑀≡(𝑀₁,𝑀₂,…,𝑀_𝑁)  and 𝑈≡𝛿≡(𝛿₁,𝛿₂,…,𝛿_𝑁) are constant vectors, then the above inequalities become 𝑓_𝑖𝑓(𝑥,𝑀,𝑀)≤0(𝑥∈Ω),_𝑖𝑀(𝑥,𝛿,𝛿)≥0(𝑥∈Ω),_𝑖≥𝜂_𝑖(𝑥,𝑡)≥𝛿_𝑖(𝑥∈Ω,𝑡∈𝐼_𝑖),𝑖=1,2,…,𝑁. (2.5) Notice that the boundary inequalities in (2.4) are trivially satisfied. For a given pair of constant vectors 𝑀, 𝜹 satisfying (2.5) we set 𝑆≡{𝑈∈𝐶(Ω)∶𝛿≤𝑈≤𝑀}, (2.6) where 𝐶(Ω) denotes the set of continuous functions on Ω. In the following we give our basic hypotheses on 𝐹(⋅,𝑈,𝑉) with 𝑉≡𝐽∗𝑈.

(H) The function 𝐹(⋅,𝑈,𝑉) is quasimonotone nondecreasing for 𝑈,𝑉∈𝑆, and for each 𝑖, 𝑓_𝑖(⋅,𝑈,𝑉)  satisfies the Lipschitz condition

|𝑓_𝑖(𝑥,𝑈,𝑉)−𝑓_𝑖(𝑥,𝑈^′,𝑉^′)|≤𝐾_𝑖(‖𝑈−𝑈^′‖+‖𝑉−𝑉^′‖)(𝑈,𝑈^′,𝑉,𝑉^′∈𝑆), (2.7) where 𝐾_𝑖 is a positive constant, and ‖𝑤‖=|𝑤₁|+|𝑤₂|+⋯+|𝑤_𝑁| for 𝑤∈𝑅^𝑁(|⋅| is a maximum norm).

It is clear that if there exists a positive constant vector 𝑀=(𝑀₁,𝑀₂,…,𝑀_𝑁) and a zero vector 𝑜=(0,0,…,0) such that

𝑓_𝑖𝑓(𝑥,𝑀,𝑀)≤0(𝑥∈Ω),_𝑖(𝑥,𝑜,𝑜)≥0(𝑥∈Ω),𝑖=1,2,…,𝑁, (2.8) then the pairs𝑈=𝑀,𝑈=𝑜 are coupled upper and lower solutions of (2.1) whenever 𝑀_𝑖≥𝜂_𝑖(𝑥,𝑡)≥0. By an application of [12,Theorem  2.1] we have the following global existence result.

Lemma 2.3. Let 𝑈≡𝑀≡(𝑀₁,𝑀₂,…,𝑀_𝑁), 𝑈≡𝛿≡(𝛿₁,𝛿₂,…,𝛿_𝑁)  be a pair of constant vectors satisfying (2.5) and 𝑀≥𝛿, and let hypothesis (H) hold. Then for any 𝜂(𝑥,𝑡) in 𝑆, problem (2.1) has a unique global solution 𝑈(𝑥,𝑡)∈𝑆 for all 𝑡≥0. In particular, if condition (2.8) holds for every 𝑀≥𝑜 then problem (2.1) has a unique bounded nonnegative global solution whenever 𝜂(𝑥,𝑡)≥𝑜.

In system (2.1) if 𝑓_𝑖 is given in the form

𝑓_𝑖(𝑥,𝑈,𝑉)=𝑢_𝑖𝑓^∗_𝑖(𝑥,𝑈,𝑉),𝑖=1,2,…,𝑁, (2.9) then the pairs 𝑈=𝑀,𝑈=𝑜 satisfy (2.8) whenever

𝑓^∗_𝑖(𝑥,𝑀,𝑀)≤0(𝑥∈Ω),𝑖=1,2,…,𝑁. (2.10) Moreover, by the positivity lemma for the linear scalar parabolic problem

𝜕𝑊_𝑖𝜕𝑡−𝐿_𝑖𝑊_𝑖=𝑊_𝑖𝑓^∗_𝑖(𝑥,𝑈,𝑉)(𝑥∈Ω,𝑡>0),𝜕𝑊_𝑖𝑊𝜕𝑛=0(𝑥∈𝜕Ω,𝑡>0),_𝑖(𝑥,0)=𝜂_𝑖(𝑥,0)𝑥∈Ω,𝑖=1,2,…,𝑁, (2.11) where 𝑈≡𝑈(𝑥,𝑡) is the solution of (2.1), the solution 𝑊_𝑖≡𝑊_𝑖(𝑥,𝑡) is positive in (0,+∞)×Ω whenever 𝜂_𝑖(𝑥,0)≥0 and not identical 0 (cf. [1]). Since for each 𝑖 the solution component 𝑢_𝑖≡𝑢_𝑖(𝑥,𝑡) is also a solution of (2.11) the uniqueness property of 𝑊_𝑖 ensures that 𝑢_𝑖(𝑥,𝑡)=𝑊_𝑖(𝑥,𝑡)>0 in (0,+∞)×Ω. This observation leads to the following.

Lemma 2.4. Let 𝑓_𝑖(𝑥,𝑈,𝑉) be given by (2.9) and satisfy hypothesis  (H),   and let 𝑀≡(𝑀₁,𝑀₂,…,𝑀_𝑁) be a constant vector satisfying (2.10). Then for any nontrivial nonnegative 𝜂_𝑖(𝑥,𝑡)≤𝑀_𝑖 problem (2.1) has a unique bounded positive solution 𝑢(𝑥,𝑡) in (0,+∞)×Ω.

To investigate the asymptotic behavior of the solution of (2.1), we consider the special case 𝑓_𝑖(𝑥,𝑈,𝑉)=𝑔_𝑖(𝑈,𝑉),𝑖=1,2,…,𝑁, (2.12) where 𝑔_𝑖(𝑈,𝑉) is not explicitly dependent on 𝑥. The corresponding steady-state problem of (2.1) is given by

−𝐿_𝑖𝑢_𝑖=𝑔_𝑖(𝑈,𝑉)(𝑥∈Ω),𝜕𝑢_𝑖𝜕𝑛=0(𝑥∈𝜕Ω,),𝑖=1,2,…,𝑁, (2.13) Assume that 𝐺(𝑈,𝑉)≡(𝑔₁(𝑈,𝑉),…,𝑔_𝑁(𝑈,𝑉)) is quasimonotone nondecreasing in 𝑆. Then by using 𝑈⁽⁰⁾=𝑀,𝑈⁽⁰⁾=𝛿 as a pair of coupled initial iterations we construct two sequences {𝑈^(𝑚)}≡{𝑢₁,…,𝑢_𝑁},{𝑈^(𝑚)}≡{𝑢₁,…,𝑢_𝑁} from the linear iteration process

−𝐿_𝑖𝑢_𝑖^(𝑚)+𝐾_𝑖𝑢_𝑖^(𝑚)=𝐾_𝑖𝑢_𝑖^(𝑚−1)+𝑔_𝑖(𝑢₁^(𝑚−1),…,𝑢_𝑁^(𝑚−1),𝐽∗𝑢₁^(𝑚−1),…,𝐽∗𝑢_𝑁^(𝑚−1)),−𝐿_𝑖𝑢_𝑖^(𝑚)+𝐾_𝑖𝑢_𝑖^(𝑚)=𝐾_𝑖𝑢_𝑖^(𝑚−1)+𝑔_𝑖(𝑢₁^(𝑚−1),…,𝑢_𝑁^(𝑚−1),𝐽∗𝑢₁^(𝑚−1),…,𝐽∗𝑢_𝑁^(𝑚−1)𝜕),𝑢_𝑖^(𝑚)=𝜕𝑛𝜕𝑢_𝑖^(𝑚)𝜕𝑛=0(𝑖=1,2,…,𝑁), (2.14) where 𝐾_𝑖 is the Lipschitz constant in (2.7). It is clear that the sequences {𝑈^(𝑚)},{𝑈^(𝑚)} are well defined. Since the initial iterations in (2.14) are the constant vectors 𝑀,𝛿 and 𝑔_𝑖(𝑈,𝑉) is independent of 𝑥 whenever 𝑈 and 𝑉 are constants we conclude from the uniqueness property of solution for linear boundary-value problems that the solutions 𝑢_𝑖^(𝑚),𝑢_𝑖^(𝑚) of (2.14) are constants and are given by 𝑢_𝑖^(𝑚)=𝑢_𝑖^(𝑚−1)+1𝐾_𝑖𝑔_𝑖(𝑢₁^(𝑚−1),…,𝑢_𝑁^(𝑚−1),𝑢₁^(𝑚−1),…,𝑢_𝑁^(𝑚−1)𝑢),_𝑖(𝑚)=𝑢_𝑖^(𝑚−1)+1𝐾_𝑖𝑔_𝑖(𝑢₁^(𝑚−1),…,𝑢_𝑁^(𝑚−1),𝑢₁^(𝑚−1),…,𝑢_𝑁^(𝑚−1))(𝑖=1,2,…,𝑁). (2.15) It is easy to show from the quasimonotone nondecreasing property of 𝐺(𝑈,𝑉) that the sequences 𝑢_𝑖^(𝑚),𝑢_𝑖^(𝑚) possess the monotone property

𝛿≤𝑈^(𝑚)≤𝑈^(𝑚+1)≤𝑈^(𝑚+1)≤𝑈^(𝑚)≤𝑀(𝑚=1,2,…) (2.16) (cf. [12]). The above property implies that the constant limits lim_𝑚→∞𝑈^(𝑚)=𝜌,lim_𝑚→∞𝑈^(𝑚)=𝜌 (2.17) exist and satisfy the relation 𝛿≤𝜌≤𝜌≤𝑀. Letting 𝑚→∞ in (2.15) shows that 𝜌 and 𝜌 satisfy the equations

𝑔_𝑖(𝜌,𝜌)=0,𝑔_𝑖(𝜌,𝜌)=0(𝑖=1,2,…,𝑁). (2.18) It is easy to verify that the limits 𝜌 and 𝜌 are maximal and minimal solutions of systems (2.13) in 𝑆, respectively, furthermore, if 𝜌=𝜌(≡𝜌^∗) then 𝜌^∗ is the unique solution in 𝑆 (cf. [1, 12]). In the latter case, the time-dependent solution of (2.1) (with 𝑓_𝑖(𝑥,𝑈,𝑉) given by (2.12) converges to 𝜌^∗ as 𝑡→∞ (cf. [12]). To summarize the above conclusions we have the following.

Lemma 2.5. Let the conditions in Lemma 2.3 be satisfied with respect to the function given by (2.12).Then (i)the sequences {𝑈^(𝑚)},{𝑈^(𝑚)} given by (2.14) with 𝑈⁽⁰⁾=𝑀,𝑈⁽⁰⁾=𝛿 are constant functions and converge monotonically to their respective constant limits𝜌 and 𝜌 that satisfy (2.18);(ii)for any 𝑚∈𝑁, if 𝑈^(𝑚)≤𝜂(𝑥,𝑡)≤𝑈^(𝑚), then the constant pairs 𝑈=𝑈^(𝑚),𝑈=𝑈^(𝑚) are also coupled upper and lower solutions of (2.1), and the problem (2.1) has a unique positive solution 𝑈(𝑥,𝑡) such that 𝑈^(𝑚)≤𝑈(𝑥,𝑡)≤𝑈^(𝑚) and 𝜌≤𝑢(𝑥,𝑡)≤𝜌,as𝑡→∞(𝑥∈Ω);(iii)if 𝜌=𝜌(≡𝜌^∗) then 𝜌^∗ is the unique solution of (2.13) in 𝑆, and for any initial function 𝜂(𝑥,𝑡)∈𝑆 the corresponding solution 𝑈(𝑥,𝑡) of (2.1) converges to 𝜌^∗as 𝑡→∞.

In Lemma 2.5 the convergence of the time-dependent solution 𝑈(𝑥,𝑡) to 𝜌^∗ is for the class of initial functions in 𝑆. For arbitrary nontrivial nonnegative initial functions we have the following result from [12].

Lemma 2.6. Let the conditions in Lemma 2.5 be satisfied, and let 𝑈≡(𝑢₁,𝑢₂,…,𝑢_𝑁) be the solution of (2.1) corresponding to an arbitrary nontrivial nonnegative 𝜂(𝑥,𝑡). Assume that 𝜌=𝜌(≡𝜌^∗) and there exists𝑡^∗>0 such that 𝛿_𝑖≤𝑢_𝑖(𝑥,𝑡)≤𝑀_𝑖for𝑡^∗−𝜏_𝑖≤𝑡≤𝑡^∗,𝑥∈Ω,(𝑖=1,2,…,𝑁). (2.19) Then 𝑈(𝑥,𝑡)converges to  𝜌^∗as 𝑡→∞.

To prove the existence and uniqueness of constant positive solution of the steady-state problem (1.13) we also need the following result.

Lemma 2.7. Let 𝑎_𝑖>0,𝑏_𝑖>0,𝑐_𝑖>0(𝑖=1,2,3). If 𝑎₁𝑏₂𝑐₃−𝑎₂𝑏₁𝑐₃−𝑎₁𝑏₃𝑐₂>0 then the inequalities 𝑎₁𝑥−𝑏₁𝑦>0,−𝑎₂𝑥+𝑏₂𝑦−𝑐₂𝑧>0,−𝑏₃𝑦+𝑐₃𝑧>0 (2.20) have a positive solution(𝜌₁,𝜌₂,𝜌₃).

Proof. With the inequalities (2.20) exists a positive solution if and only if there exist positive constants 𝜆₁,𝜆₂,𝜆₃ such that the following equations: 𝑎₁𝑥−𝑏₁𝑦=𝜆₁,−𝑎₂𝑥+𝑏₂𝑦−𝑐₂𝑧=𝜆₂,−𝑏₃𝑦+𝑐₃𝑧=𝜆₃ (2.21) have a positive solution. Let 𝐴=(^{𝑎1−𝑏10−𝑎}₂^𝑏₂^−𝑐_{20−𝑏3𝑐3}), we have |𝐴|>0. From the Cramer theorem, the solutions of (2.21) are as follows: 𝜆𝑥=¹|^{𝑏2−𝑐2}_{−𝑏3𝑐3}|−𝜆₂|^−𝑏10_{−𝑏3𝑐3}|+𝜆₃|^−𝑏10_{𝑏2−𝑐2}|=𝜆|𝐴|₁(𝑏₂𝑐₃−𝑏₃𝑐₂)+𝜆₂𝑏₁𝑐₃+𝜆₃𝑏₁𝑐₃,|𝐴|𝑦=−𝜆₁|^{−𝑎2−𝑐2}_0𝑐3|+𝜆₂|^𝑎10_0𝑐3|−𝜆₃|^𝑎10_{−𝑎2−𝑐2}|=𝜆|𝐴|₁𝑎₂𝑐₃+𝜆₂𝑎₁𝑐₃+𝜆₃𝑎₁𝑐₂,𝜆|𝐴|𝑧=¹|^{−𝑎2𝑏2}_{0−𝑏𝑐3}|−𝜆₂|^{𝑎1−𝑏1}_0−𝑏3|+𝜆₃|^{𝑎1−𝑏1}_0−𝑏3|=𝜆|𝐴|₁𝑎₂𝑏₃+𝜆₂𝑎₁𝑏₃+𝜆₃(𝑎₁𝑏₂−𝑎₂𝑏₁).|𝐴| (2.22) Let 𝜆₁=𝜆₃=1, 𝜆₂=max{(1−𝑏₁𝑐₃−𝑏₂𝑐₃+𝑏₃𝑐₂)/𝑏₁𝑐₃,(1−𝑎₁𝑏₂+𝑎₂𝑏₁−𝑎₂𝑏₃)/𝑎₁𝑏₃,1}, then with the equations (2.21) exists a positive solution (𝜌₁,𝜌₂,𝜌₃). Therefore the inequalities (2.20) have a positive solution.

3. The Main Theorems and Proof

To prove the main results in this paper, we apply Lemmas 2.3–2.7 for the models in (1.1)-(1.2) with 𝑈=(𝑢,𝑣,𝑤) by constructing a suitable pair of constant upper and lower solutions 𝑈=𝑀 and 𝑈=𝛿. This is equivalent to show that 𝑀 and 𝜹 satisfy condition (2.5) for the reaction functions 𝑓_𝑖(𝑥,𝑈,𝑉).

Theorem 3.1. Let 𝑒^∗_𝑖≡𝑒_𝑖+𝛽_𝑖>0(𝑖=1,2,3) and 𝑏^∗₂≡𝑏₂+𝛼₂>0. If 𝑐₁𝑐₂𝑐₃−𝑐₁𝑒^∗₂𝑒^∗₃−𝑒^∗₁𝑏^∗₂𝑐₃>0, then the following statements hold true. (i)The steady-state problem (1.9) has a unique constant positive solution.(ii) For any nontrivial nonnegative (𝜂₁,𝜂₂,𝜂₃), problems (1.1)-(1.2) have a unique nonnegative global solution (𝑢,𝑣,𝑤), and the solution is uniformly boundary in Ω×[0,∞).(iii)The positive constant solution 𝜌^∗of the steady-state problem (1.9) is globally asymptotic stable.

Proof. (i)Since the boundary conditions are   Neumann type, we seek a pair of positive constant vectors 𝑀=(𝑀₁,𝑀₂,𝑀₃),  𝛿=(𝛿₁,𝛿₂,𝛿₃)satisfying 𝑀₁(𝑎₁−𝑐₁𝑀₁+𝑒^∗₁𝑀₂)≤0≤𝛿₁(𝑎₁−𝑐₁𝛿₁+𝑒^∗₁𝛿₂𝑀),₂(𝑎₂+𝑏^∗₂𝑀₁−𝑐₂𝑀₂+𝑒^∗₂𝑀₃)≤0≤𝛿₂(𝑎₂+𝑏^∗₂𝛿₁−𝑐₂𝛿₂+𝑒^∗₂𝛿₃𝑀),₃(𝑎₃+𝑒^∗₃𝑀₂−𝑐₃𝑀₃)≤0≤𝛿₃(𝑎₃+𝑒^∗₃𝛿₂−𝑐₃𝛿₃),𝜕𝛿_𝑖𝜕𝑛≤0≤𝜕𝑀_𝑖𝜕𝑛(𝑖=1,2,3). (3.1) It is clear that the boundary conditions are satisfied, and from the positive property of 𝑀_𝑖and 𝛿_𝑖 that condition (3.1) is satisfied if (𝑎₁−𝑐₁𝑀₁+𝑒^∗₁𝑀₂)≤0≤(𝑎₁−𝑐₁𝛿₁+𝑒^∗₁𝛿₂),(𝑎₂+𝑏^∗₂𝑀₁−𝑐₂𝑀₂+𝑒^∗₂𝑀₃)≤0≤(𝑎₂+𝑏^∗₂𝛿₁−𝑐₂𝛿₂+𝑒^∗₂𝛿₃),(𝑎₃+𝑒^∗₃𝑀₂−𝑐₃𝑀₃)≤0≤(𝑎₃+𝑒^∗₃𝛿₂−𝑐₃𝛿₃). (3.2) We construct the following inequalities: 𝑐₁𝑥−𝑒^∗₁𝑦>0,−𝑏^∗₂𝑥+𝑐₂𝑦−𝑒^∗₂𝑧>0,−𝑒^∗₃𝑦+𝑐₃𝑧>0. (3.3) Let 𝐴=(^{𝑐1−𝑒∗10−𝑏∗}₂^𝑐₂^−𝑒∗_{20−𝑒∗3𝑐3}), we have |𝐴|>0. Since 𝑒^∗_𝑖≡𝑒_𝑖+𝛽_𝑖>0(𝑖=1,2,3) and 𝑏^∗₂≡𝑏₂+𝛼₂>0, from the Lemma 2.7 the inequalities (3.3) have a positive solution (𝜆₁,𝜆₂,𝜆₃). By choosing 𝑀=(𝑀₁,𝑀₂,𝑀₃)=(𝑃𝜆₁,𝑃𝜆₂,𝑃𝜆₃), 𝛿=(𝛿₁,𝛿₂,𝛿₃)=(𝛾𝜆₁,𝛾𝜆₂,𝛾𝜆₃), where 𝑃 and 𝛾 are a sufficiently large and small positive constant, respectively, we see that all the inequalities in (3.2) are satisfied. That is, 𝑀 and 𝛿 is a pair of positive constant upper-lower solution of (1.1)-(1.2). By Lemma 2.5 the sequences governed by (2.14) with the reaction function (𝑓₁,𝑓₂,𝑓₃) given by (1.1) converge monotonically, respectively, to some constant limits 𝜌≡(𝜌₁,𝜌₂,𝜌₃),𝜌≡(𝜌₁,𝜌₂,𝜌₃) that satisfy 𝜌₁(𝑎₁−𝑐₁𝜌₁+𝑒^∗₁𝜌₂)=0,𝜌₁(𝑎₁−𝑐₁𝜌₁+𝑒^∗₁𝜌₂)=0,𝜌₂(𝑎₂+𝑏^∗₂𝜌₁−𝑐₂𝜌₂+𝑒^∗₂𝜌₃)=0,𝜌₂(𝑎₂+𝑏^∗₂𝜌₁−𝑐₂𝜌₂+𝑒^∗₂𝜌₃)=0,𝜌₃(𝑎₃+𝑒^∗₃𝜌₂−𝑐₃𝜌₃)=0,𝜌₃(𝑎₃+𝑒^∗₃𝜌₂−𝑐₃𝜌₃)=0. (3.4) Let 𝜌_𝑖=𝜌_𝑖−𝜌_𝑖(𝑖=1,2,3). Then by the positivity of 𝜌_𝑖 and 𝜌_𝑖, a subtraction of the above pairs of equations, leads to 𝑐₁𝜌₁−𝑒^∗₁𝜌₂=0,−𝑏^∗₂𝜌₁+𝑐₂𝜌₂−𝑒^∗₂𝜌₃=0,−𝑒^∗₃𝜌₂+𝑐₃𝜌₃=0. (3.5) Since these equations are equivalent to 𝐴𝜌=0, where 𝐴 is coefficient matrix and 𝜌 is the column vector (𝜌₁,𝜌₂,𝜌₃)^𝑇, we conclude from the nonsingular property of 𝐴 that 𝜌=𝑜. This proves 𝜌=𝜌(≡𝜌^∗) and 𝜌^∗ is the unique positive solution of (1.13) in 𝑆.
(ii)It is obvious that the pairs 𝑀=(𝑀₁,𝑀₂,𝑀₃)=(𝑃𝜆₁,𝑃𝜆₂,𝑃𝜆₃),𝛿=(𝛿₁,𝛿₂,𝛿₃)=(0,0,0) satisfy all the equalities in (3.1). By Lemma 2.3 and the arbitrariness of P, problems (1.1)-(1.2) have a unique bounded nonnegative global solution (𝑢,𝑣,𝑤) for any nontrivial nonnegative (η_1,η_2,η_3).
(iii)Lemma 2.5 ensures that for any initial function 𝜂(𝑥,𝑡)≡(𝜂₁(𝑥,𝑡),𝜂₂(𝑥,𝑡),𝜂₃(𝑥,𝑡)) with 𝛾(𝜆₁,𝜆₂,𝜆₃)≡𝛿≤𝜂(𝑥,𝑡)≤𝑀≡𝑃(𝜆₁,𝜆₂,𝜆₃), where 𝑃 and 𝛾are a sufficiently large and small positiveconstants, respectively, the corresponding solution 𝑈(𝑥,𝑡)≡(𝑢(𝑥,𝑡),𝑣(𝑥,𝑡),𝑤(𝑥,𝑡)) of (1.1)-(1.2) converges to 𝜌^∗ as 𝑡→∞. To show this convergence property for an arbitrary nontrivial nonnegative 𝜂(𝑥,𝑡), we observe from Lemma 2.4, the arbitrary smallness of 𝛾 and the arbitrary largeness of 𝑃 that there exists 𝑡₀>0 such that 𝛿≤𝑈(𝑥,𝑡)≤𝑀 on [𝑡₀,∞)×Ω. The upper and lower bounds of 𝑈(𝑥,𝑡)show that there exists 𝑡^∗>0 such that condition (2.19) is satisfied for 𝑈≡(𝑢₁,𝑢₂,𝑢₃)=(𝑢,𝑣,𝑤). The conclusion of 𝑈(𝑥,𝑡)→𝜌^∗ as 𝑡→∞ follows from Lemma 2.6. We next consider the stability of 𝑈(𝑥,𝑡). We observe from Lemma 2.5 that for any 𝑈(𝜌^∗,𝜀) there exists 𝛿>0 such that 𝑈(𝑥,𝑡)∈𝑈^∗≡𝑈(𝜌^∗,𝛿)⊆(𝑈^(𝑚),𝑈^(𝑚))⊂𝑈(𝜌^∗,𝜀) as initial function 𝜂(𝑥,𝑡)∈𝑈^∗. This completes the proof.

Corollary 3.2. All the results of Theorem 3.1 hold true if 𝑑_𝑖=0 for some or all 𝑖. In particular, these results hold true for the ordinary differential system (1.8). They are also true if the reaction functions involve no time delays.

Proof. This follows from the argument in the proof of Theorem 3.1 by letting 𝑑_𝑖=0 (and without the corresponding boundary condition) and by letting 𝜏_𝑖=0(𝑖=1,2,3), respectively.

Remark 3.3. Our model and result are different from the existence ones such as those of Pao [29] and Kim and Lin [22, 23]. In some sense, we enrich the results of the 3-specics Lotka-Volterra reaction-diffusion systems.

4. Numerical Simulations

In this section, we give numerical simulations supporting our theoretical analysis. As an example, we consider system (1.1) with different diffusion rates 𝑑_𝑖, birth rate 𝑎_𝑖, and time delays 𝜏_𝑖, that is, the following systems: 𝜕𝑢𝜕𝑡−𝑢_𝑥𝑥1=𝑢(𝑡,𝑥)[1−𝑢(𝑡,𝑥)+21𝑣(𝑡,𝑥)+2𝑣(𝑡−0.1,𝑥)](2𝜋>𝑥>0,𝑡>0),𝜕𝑣𝜕𝑡−𝑣_𝑥𝑥1=𝑣(𝑡,𝑥)[1+21𝑢(𝑡,𝑥)−3𝑣(𝑡,𝑥)+2+1𝑤(𝑡,𝑥),21𝑢(𝑡−0.1,𝑥)+2𝑤(𝑡−0.1,𝑥)](2𝜋>𝑥>0,𝑡>0),𝜕𝑤𝜕𝑡−𝑤_𝑥𝑥1=𝑤(𝑡,𝑥)[1+21𝑣(𝑡,𝑥)−𝑤+(2𝑣(𝑡−0.1,𝑥))](2𝜋>𝑥>0,𝑡>0),𝜕𝑢=𝜕𝑛𝜕𝑣=𝜕𝑛𝜕𝑤𝜋𝜕𝑛=0(𝑡>0,𝑥=0,2𝜋),𝑢(𝑡,𝑥)=𝑣(𝑡,𝑥)=𝑤(𝑡,𝑥)=(1+10𝑡)[1−sin(𝑥+2)](𝑡,𝑥)∈[−0.1,0]×[0,2𝜋].𝜕𝑢𝜕𝑡−𝑢_𝑥𝑥1=𝑢(𝑡,𝑥)[1−𝑢(𝑡,𝑥)+21𝑣(𝑡,𝑥)+2𝑣(𝑡−0.01,𝑥)](2𝜋>𝑥>0,𝑡>0),𝜕𝑣𝜕𝑡−𝑣_𝑥𝑥1=𝑣(𝑡,𝑥)[1+21𝑢(𝑡,𝑥)−3𝑣(𝑡,𝑥)+2+1𝑤(𝑡,𝑥),21𝑢(𝑡−0.01,𝑥)+2𝑤(𝑡−0.01,𝑥)](2𝜋>𝑥>0,𝑡>0),𝜕𝑤𝜕𝑡−𝑤_𝑥𝑥1=𝑤(𝑡,𝑥)[1+21𝑣(𝑡,𝑥)−𝑤+2𝑣(𝑡−0.01,𝑥)](2𝜋>𝑥>0,𝑡>0),𝜕𝑢=𝜕𝑛𝜕𝑣=𝜕𝑛𝜕𝑤𝜋𝜕𝑛=0(𝑡>0,𝑥=0,2𝜋),𝑢(𝑡,𝑥)=𝑣(𝑡,𝑥)=𝑤(𝑡,𝑥)=(1+𝑡)[1−sin(𝑥+2)](𝑡,𝑥)∈[−0.01,0]×[0,2𝜋].𝜕𝑢𝜕𝑡−𝑢_𝑥𝑥1=𝑢(𝑡,𝑥)[1−𝑢(𝑡,𝑥)+21𝑣(𝑡,𝑥)+2𝑣(𝑡−0.1,𝑥)](2𝜋>𝑥>0,𝑡>0),𝜕𝑣𝜕𝑡−𝑣_𝑥𝑥1=𝑣(𝑡,𝑥)[2+21𝑢(𝑡,𝑥)−3𝑣(𝑡,𝑥)+2+1𝑤(𝑡,𝑥),21𝑢(𝑡−0.1,𝑥)+2𝑤(𝑡−0.1,𝑥)](2𝜋>𝑥>0,𝑡>0),𝜕𝑤𝜕𝑡−𝑤_𝑥𝑥1=𝑤(𝑡,𝑥)[3+21𝑣(𝑡,𝑥)−𝑤+2𝑣(𝑡−0.1,𝑥)](2𝜋>𝑥>0,𝑡>0),𝜕𝑢=𝜕𝑛𝜕𝑣=𝜕𝑛𝜕𝑤𝜋𝜕𝑛=0(𝑡>0,𝑥=0,2𝜋),𝑢(𝑡,𝑥)=𝑣(𝑡,𝑥)=𝑤(𝑡,𝑥)=(1+10𝑡)[1−sin(𝑥+2)](𝑡,𝑥)∈[−0.1,0]×[0,2𝜋].𝜕𝑢𝜕𝑡−𝑢_𝑥𝑥1=𝑢(𝑡,𝑥)[1−𝑢(𝑡,𝑥)+21𝑣(𝑡,𝑥)+2𝑣(𝑡−0.1,𝑥)](2𝜋>𝑥>0,𝑡>0),𝜕𝑣𝜕𝑡−2𝑣_𝑥𝑥1=𝑣(𝑡,𝑥)[1+21𝑢(𝑡,𝑥)−3𝑣(𝑡,𝑥)+2+1𝑤(𝑡,𝑥),21𝑢(𝑡−0.1,𝑥)+2𝑤(𝑡−0.1,𝑥)](2𝜋>𝑥>0,𝑡>0),𝜕𝑤𝜕𝑡−3𝑤_𝑥𝑥1=𝑤(𝑡,𝑥)[1+21𝑣(𝑡,𝑥)−𝑤+2𝑣(𝑡−0.1,𝑥)](2𝜋>𝑥>0,𝑡>0),𝜕𝑢=𝜕𝑛𝜕𝑣=𝜕𝑛𝜕𝑤𝜋𝜕𝑛=0(𝑡>0,𝑥=0,2𝜋),𝑢(𝑡,𝑥)=𝑣(𝑡,𝑥)=𝑤(𝑡,𝑥)=(1+10𝑡)[1−sin(𝑥+2)](𝑡,𝑥)∈[−0.1,0]×[0,2𝜋]. (4.1)(1)(2)(3)(4)(5)(4.2)(4.3)(4.4) By using the classical implicit format solving the partial differential equations and the method of steps for differential difference equations and employing the software package MATLAB7.0, we can solve the numerical solutions of systems (4.1), (4.2), (4.3), and (4.4) which are shown respectively in Figures 1, 2, 3, and 4.

(a)

(b)

(c)

Figure 1 

The solution of system (4.1).

(a)

(b)

(c)

Figure 2 

The solution of system (4.2).

(a)

(b)

(c)

Figure 3 

The solution of system (4.3).

(a)

(b)

(c)

Figure 4 

The solution of system (4.4).

5. Conclusions

It is well known that the analysis of stability for a system of delay reaction-diffusion multispecies model is quite difficult since the reaction among multispecies is more complex. Therefore, the works on this subject are very rare. A detailed analysis on the stability for a two-prey one-predator model, one-prey two-predator model, and three-species food-chain model with delay and diffusion was given by Pao [29], and he obtained some simple and easily verifiable conditions for the existence and global asymptotic stability of a positive steady-state solution for each of the three model problems.

In this paper, based on the ideas of Pao [12], we have considered a delay cooperative three-species system with Neumann boundary condition. It is shown that the system has a positive equilibrium under some certain conditions. We have obtained the similar conclusions to those of Pao [29]. More precisely, we have obtained the following results.

(a) The steady-state problem (1.9) has a unique constant positive solution if 𝑐₁𝑐₂𝑐₃−𝑐₁𝑒^∗₂𝑒^∗₃−𝑒^∗₁𝑏^∗₂𝑐₃>0.

(b) For any nontrivial nonnegativ (𝜂₁,𝜂₂,𝜂₃), problems (1.1)-(1.2) have a unique nonnegative global solution (𝑢,𝑣,𝑤), and the solution is uniformly boundary in Ω×[0,∞).

(c) The positive constant solution 𝜌^∗ of the steady-state problem (1.9) is globally asymptotic stable.

The condition 𝑐₁𝑐₂𝑐₃−𝑐₁𝑒^∗₂𝑒^∗₃−𝑒^∗₁𝑏^∗₂𝑐₃>0 involves only the reaction rate constants, which shows that the diffusion rates 𝑑_𝑖, the birth rate 𝑎_𝑖, and time delays 𝜏_𝑖 do not bring effect on permanence of one species as well as contribution to its extinction. The result of global asymptotic stability implies that the three-species model system coexists, is permanent, and the trivial and all semitrivial solutions are unstable.

Acknowledgments

The authors are grateful to the referee for her/his comments. This work is supported by Science and Technology Study Project of Chongqing Municipal Education Commission (Grant no. KJ 080511) of China, Natural Science Foundation Project of CQ CSTC (Grant no. 2008BB7415 and 2007BB2450) of China, Foundation Project of Doctor Graduate Student Innovation of Beijing University of Technology of China, the NSFC (Grant nos.10471009), and BSFC (Grant no. 1052001) of China.