Persistence and Stability for a Generalized Leslie-Gower Model with Stage Structure and Dispersal

Abstract

A generalized version of the Leslie-Gower predator-prey model that incorporates the prey structure and predator dispersal in two-patch environments is introduced. The focus is on the study of the boundedness of solution, permanence, and extinction of the model. Sufficient conditions for global asymptotic stability of the positive equilibrium are derived by constructing a Lyapunov functional. Numerical simulations are also presented to illustrate our main results.

1. Introduction

Lotka-Volterra predator-prey models have been extensively and deeply investigated (see monographs [1–5]). If we let 𝑥(𝑡) denote the density of prey and let 𝑦(𝑡) be the density of predator, then the classical Lotka-Volterra predator-prey model is given by the following system: 𝑑𝑥𝑑𝑡=(𝑟₁−𝑐₁𝑦−𝑏₁𝑥)𝑥,𝑑𝑦𝑑𝑡=(−𝜀₂+𝜌₂𝑥)𝑦. (1.1) The equations in system (1.1) set no upper limit on the per-capita growth rate of the predator (the second term of model (1.1) which is unrealistic. For example, for mammals, such a limit will be determined in part by physiological factors (length of the gestation period, the shortest interval between litters, the maximum average number of daughters per-litter, the age at which breeding first starts, and so on [6]). Leslie modeled the effect of such limitations via a predator-prey model where the “carrying capacity’’ of the predator's environment was assumed to be proportional to the number of prey. Hence, if 𝑥(𝑡) and 𝑦(𝑡) denote the prey and predator density, respectively, then Leslie's model is given by the following system: 𝑑𝑥𝑑𝑡=(𝑟₁−𝑐₁𝑦−𝑏₁𝑥)𝑥,𝑑𝑦𝑑𝑡=(𝑟₂−𝑐₂𝑦𝑥)𝑦, (1.2) where 𝑟_𝑖,𝑐_𝑖,𝑖=1,2, and 𝑏₁ are positive constants. The first equation of system (1.2) is standard but the second is not because it contains the “so-called’’ Leslie-Gower term, namely, 𝑐₂(𝑦/𝑥). The rationale behind this term is based on the view that as the prey becomes numerous (𝑥→∞) then the per-capita growth rate of the predator achieves its maximum ((1/𝑦)(𝑑𝑦/𝑑𝑡)→𝑟₂). Conversely as the prey becomes scarce 𝑥→0, we have that (1/𝑦)(𝑑𝑦/𝑑𝑡)→−∞. That is, the predator must go extinct. Recently, the use of a Holling-type II functional for the prey has led various researchers [7, 8] to the consideration of the following model (a modification of system (1.2): 𝑑𝑥𝑑𝑡=(𝑟₁−𝑐₁𝑦𝑥+𝑘₁−𝑏₁𝑥)𝑥,𝑑𝑦𝑑𝑡=(𝑟₂−𝑐₂𝑦𝑥+𝑘₂)𝑦, (1.3) where 𝑟₁ is the per-capita growth rate of the prey 𝑥; 𝑏₁ is a measure of the strength of prey (on prey) interference competition; 𝑐₁ is the maximum value of the per-capita reduction rate of 𝑥 due to 𝑦; 𝑘₁ measures the extent to which the environment provides protection to prey 𝑥 (𝑘₂ for 𝑦); 𝑟₂ gives the maximal per-capita growth rate of 𝑦; 𝑐₂ has a similar meaning to that of 𝑐₁.

In [9], the global stability of the unique coexisting interior equilibrium of system (1.2) is established. In [7], the existence and boundedness of solutions (including that of an attracting set) are established as well as the global stability of the coexisting interior equilibrium for model (1.3). There have been additional extensions, for example, in [10, 11] a Leslie-Gower type model with impulse was introduced and investigated.

The study of the role of dispersal in continuous-time metapopulation models is extensive (see [12–16] and the references cited therein). They show that dispersal can have a stabilizing influence on the system (see [12, 13]) and also can have a destabilizing influence on the system (see [14, 15]).

On the other hand, most prey species have a life history that includes multiple stages (juvenile and adults or immature and mature). In Aiello and Freedman [17], the population dynamics of a single species with two identifiable stages was modeled by the following system: 𝑥^′₁(𝑡)=𝛼𝑥₂(𝑡)−𝛾𝑥₁(𝑡)−𝛼𝑒^−𝛾𝜏𝑥₂𝑥(𝑡−𝜏),_′2(𝑡)=𝛼𝑒^−𝛾𝜏𝑥₂(𝑡−𝜏)−𝛽𝑥²₂(𝑡), (1.4) where 𝑥₁(𝑡),𝑥₂(𝑡) denote the immature and mature population densities, respectively. Here, 𝛼>0 represents the per-capita birth rate; 𝛾>0 is the per-capita immature death rate; 𝛽>0 is the death rate due to overcrowding, and 𝜏 is the “fixed’’ time to maturity; the term 𝛼𝑒^−𝛾𝜏𝑥₂(𝑡−𝜏) models the immature individuals who were born at time 𝑡−𝜏 (i.e., 𝛼𝑥₂(𝑡−𝜏)) and survive and mature at time 𝑡. The derivation and analysis of system (1.4) can be found in [17]. More and More researchers (see [16–22] and the references cited therein) have investigated many kinds of predator-prey model under various stage-structure assumptions. In Xu et al. [16], they discussed a Lotka-Volterra-type predator-prey model with stage structure for predator and prey dispersal in two-patch environments. They obtained sufficient conditions of permanence and impermanence and global asymptotic stability of the positive equilibrium; they also discussed the local stability of the positive equilibrium. In [22], they studied a generalized version of the Leslie-Gower predator-prey model that incorporates the prey structure and obtained sufficient conditions of permanence and stability of the nonnegative equilibrium.

Motivated by the above works, in this paper we study the effects of stage structure for prey and predator dispersal on the global dynamics of modified version of the Leslie-Gower and Holling-type II predator-prey system. Following [16, 23], we assumethe following.

(A 1) The prey population : the prey only lives in patch 1. For immature prey, 𝛼 is birth rate, 𝑟₁ is death rate, and the term 𝛼𝑒^−𝑟1𝜏𝑥₂(𝑡−𝜏) represents the number of immature prey that was born at time 𝑡−𝜏, which still survive at time 𝑡 and are transferred from the immature stage to the mature stage at time 𝑡. For mature prey, 𝑟₂ is death rate, 𝑟₃ is the intraspecific competition rate of mature prey, 𝑎₁ is the maximum value of the per-capita reduction rate of 𝑥₂ due to 𝑦₁, and 𝑘₁ (resp., 𝑘₂) measures the extent to which environment provides protection to prey 𝑥₂ (resp., to the predator 𝑦₁).

(A 2) The predator population : 𝛽_𝑖 are the birth rate of predator in patch 𝑖, 𝑖=1,2; 𝐷_𝑖 is the dispersion rate of predator between two patches; 𝑟₄ is death rate of predator in patch 2; 𝑎₂ has a similar meaning to 𝑎₁. It is assumed that predators in patch 1 do not capture immature prey, then we have the following delayed differential system: 𝑥̇₁(𝑡)=𝛼𝑥₂(𝑡)−𝑟₁𝑥₁(𝑡)−𝛼𝑒^−𝑟1𝜏𝑥₂𝑥̇(𝑡−𝜏),₂(𝑡)=𝛼𝑒^−𝑟1𝜏𝑥₂(𝑡−𝜏)−𝑟₂𝑥₂(𝑡)−𝑟₃𝑥²₂𝑎(𝑡)−¹𝑦₁(𝑡)𝑥₂(𝑡)𝑥₂(𝑡)+𝑘₁,𝑦̇₁(𝑡)=(𝛽₁−𝑎₂𝑦₁(𝑡)𝑥₂(𝑡)+𝑘₂)𝑦₁(𝑡)+𝐷₁(𝑦₂(𝑡)−𝑦₁𝑦̇(𝑡)),₂(𝑡)=(𝛽₂−𝑟₄𝑦₂(𝑡))𝑦₂(𝑡)+𝐷₂(𝑦₁(𝑡)−𝑦₂(𝑡)), (1.5) where 𝑥₁(𝑡) and 𝑥₂(𝑡) represent the densities of immature and mature individual prey in patch 1 at time 𝑡, 𝑦_𝑖(𝑡) denote the density of predator species in patch 𝑖,𝑖=1,2 at time 𝑡, all parameters of (1.5) are positive constants.

The initial conditions for system (1.5) take the form of 𝑥_𝑖(𝜃)=Φ_𝑖(𝜃),𝑦_𝑖(𝜃)=Ψ_𝑖(𝜃),𝑥_𝑖(0)>0,𝑦_𝑖(0)>0,𝑖=1,2, (1.6) where (Φ₁(𝜃),Φ₂(𝜃),Ψ₁(𝜃),Ψ₂(𝜃))∈𝐶([−𝜏,0],𝑅⁴₊₀), the Banach space of continuous function mapping the interval [−𝜏,0] into 𝑅⁴₊₀, where 𝑅⁴₊₀={(𝑥₁,𝑥₂,𝑥₃,𝑥₄)∶𝑥_𝑖≥0,𝑖=1,2,3,4}.

For continuity of the initial conditions, we further require 𝑥₁(0)=⁰_−𝜏𝛼𝑒^𝑟1𝑠Φ₂(𝑠)𝑑𝑠. (1.7)

The paper is organized as follows. In Section 2, we will discuss the uniform persistence of system (1.5). In Section 3, we are concerned with the global stability of a positive equilibrium of system (1.5) by constructing Lyapunov functional and also present two numerical simulations to illustrate our main results.

2. Uniform Persistence

In this section, we will discuss the uniform persistence of system (1.5) with initial conditions (1.6) and (1.7).

Definition 2.1. System (1.5) is said to be uniformly persistent if there exists a compact region 𝐷⊂Int𝑅⁴₊₀ such that every solution of system (1.5) with initial conditions (1.6) and (1.7) eventually enters and remains in the region 𝐷.

Lemma 2.2. Solutions of system (1.5) with initial conditions (1.6) and (1.7) are positive for all 𝑡≥0.

Proof. Let (𝑥₁(𝑡),𝑥₂(𝑡),𝑦₁(𝑡),𝑦₂(𝑡)) be a solution of system (1.5) with initial conditions (1.6) and (1.7); we first consider 𝑦₁(𝑡) and 𝑦₂(𝑡) for 𝑡∈[0,𝜏],𝑦̇₁(𝑡)|_𝑦1=0=𝐷₁𝑦₂(𝑡)>0for𝑦₂𝑦̇>0,₂(𝑡)|_𝑦2=0=𝐷₂𝑦₁(𝑡)>0for𝑦₁>0. (2.1) Thus, it follows that 𝑦₁(𝑡)>0,𝑦₂(𝑡)>0 for 𝑡∈[0,𝜏].
For 𝑡∈[0,𝜏], it follows from the second equation of system (1.5) that 𝑥̇₂(𝑡)≥[−𝑟₂−𝑟₃𝑥₂𝑎(𝑡)−¹𝑦₁(𝑡)𝑥₂(𝑡)+𝑘₁]𝑥₂(𝑡). (2.2)
Consider the following auxiliary equation: 𝑢̇(𝑡)≥[−𝑟₂−𝑟₃𝑎𝑢(𝑡)−¹𝑦₁(𝑡)𝑢(𝑡)+𝑘₁]𝑢(𝑡),𝑢(𝑡)=𝑢(0)exp(−^𝑡₀(𝑟₂+𝑟₃𝑎𝑢(𝑠)+¹𝑦₁(𝑠)𝑢(𝑠)+𝑘₁)𝑑_𝑠)>0. (2.3) For 𝑡∈[0,𝜏],𝑢(0)=𝑥₂(0)>0; thus, 𝑥₂(𝑡)≥𝑢(𝑡)>0.
In a similar way, we consider the intervals [𝜏,2𝜏]⋅⋅⋅[𝑛𝜏,(𝑛+1)𝜏],𝑛∈𝑁. Thus, 𝑥₂(𝑡)>0,𝑦₁(𝑡)>0,𝑦₂(𝑡)>0 for all 𝑡≥0.
By (1.7) and the first equation of (1.5) we can obtain that 𝑥₁(𝑡)=𝛼^𝑡_𝑡−𝜏𝑒^{−𝑟1(𝑡−𝑠)}𝑥₂(𝑠)𝑑_𝑠. (2.4) Therefore the positivity of 𝑥₁(𝑡) for 𝑡≥0 follows, this completes the proof.

In order to discuss the uniform persistence, we need the following result from [24].

Lemma 2.3. Consider the following equation: 𝑥̇(𝑡)=𝑎𝑥(𝑡−𝜏)−𝑏𝑥(𝑡)−𝑐𝑥²(𝑡), (2.5) where 𝑎,𝑏,𝑐, and 𝜏 are positive constants, 𝑥(𝑡)>0 for 𝑡∈[−𝜏,0]. We have the following:(i)if 𝑎>𝑏, then lim_𝑡→+∞𝑥(𝑡)=(𝑎−𝑏)/𝑐;(ii)if 𝑎<𝑏, then lim_𝑡→+∞𝑥(𝑡)=0.

Lemma 2.4. Let (𝑥₁(𝑡),𝑥₂(𝑡),𝑦₁(𝑡),𝑦₂(𝑡)) be a solution of system (1.5) with initial conditions (1.6) and (1.7). Then there exists a 𝑇₃>0 such that 𝑥_𝑖(𝑡)≤𝑁,𝑦_𝑖(𝑡)≤𝑁,(𝑖=1,2)for𝑡≥𝑇₃, (2.6) where 𝑁 is a constant and 𝑁>max{𝑁₁,𝑁₂,𝑁^∗𝑁},₁=𝛼𝑁₂𝑟₁(1−𝑒^−𝑟1𝜏𝑁),₂=𝛼𝑒^−𝑟1𝜏𝑟₃𝑁+𝜀,_∗=𝛼²4𝐴𝑟₃+(𝐴+𝐷₂+𝛽₁)²(𝑁₂+𝑘₂)4𝐴𝑎₂+(𝐴+𝐷₁+𝛽₂)²4𝐴𝑟₁,𝐴=min{𝑟_1,𝑟₂}. (2.7)

Proof. Suppose 𝑋(𝑡)=(𝑥₁(𝑡),𝑥₂(𝑡),𝑦₁(𝑡),𝑦₂(𝑡))to be any positive solution of system (1.5) with initial conditions (1.6) and (1.7). It follows from the second equation of system (1.5) that for 𝑡≥𝜏,𝑥̇₂(𝑡)≤𝛼𝑒^−𝑟1𝜏𝑥₂(𝑡−𝜏)−𝑟₃𝑥²₂(𝑡). (2.8)
Consider the following auxiliary equation: 𝑢̇(𝑡)=𝛼𝑒^−𝑟1𝜏𝑢(𝑡−𝜏)−𝑟₃𝑢²(𝑡). (2.9) By Lemma 2.3 we obtain that lim_𝑡→+∞𝑢(𝑡)=𝛼𝑒^−𝑟1𝜏𝑟₃. (2.10) Using comparison principle, it follows that lim_𝑡→+∞sup𝑥₂(𝑡)≤𝛼𝑒^−𝑟1𝜏𝑟₃. (2.11) Therefore, for sufficiently small 𝜀>0 there is a 𝑇₁>𝜏 such that if 𝑡≥𝑇₁,𝑥₂(𝑡)≤𝛼𝑒^−𝑟1𝜏𝑟₃+𝜀∶=𝑁₂. (2.12) Setting 𝑇₂=𝑇₁+𝜏, it then follows (2.4) and (2.12) that, for 𝑡≥𝑇₂,𝑥₁(𝑡)≤𝛼𝑁₂𝑟₁(1−𝑒^−𝑟1𝜏)∶=𝑁₁. (2.13) We define 𝜌(𝑡)=𝑥₁(𝑡)+𝑥₂(𝑡)+𝑦₁(𝑡)+𝑦₂̇𝛼(𝑡),𝜌(𝑡)≤−𝐴𝜌(𝑡)+²4𝑟₃+(𝐴+𝐷₂+𝛽₁)²(𝑁₂+𝑘₂)4𝑎₂+(𝐴+𝐷₁+𝛽₂)²4𝑟₁, (2.14) where 𝐴=min{𝑟_1,𝑟₂}.
It follows from (2.14) that lim_𝑡→+∞𝛼sup𝜌(𝑡)≤²4𝐴𝑟₃+(𝐴+𝐷₂+𝛽₁)²(𝑁₂+𝑘₂)4𝐴𝑎₂+(𝐴+𝐷₁+𝛽₂)²4𝐴𝑟₁∶=𝑁^∗. (2.15) Therefore, there exists a 𝑇₃=𝑇₂+𝜏 and 𝑁>max{𝑁₁,𝑁₂,𝑁^∗}. (2.16) Such that if 𝑡≥𝑇₃,𝑥_𝑖(𝑡)≤𝑁,𝑦_𝑖(𝑡)≤𝑁(𝑖=1,2). This completes the proof.

Theorem 2.5. System (1.5) with initial conditions (1.6) and (1.7) is uniformly persistent provided that(𝐻1)𝛼𝑒^−𝑟1𝜏>𝑟₂+𝑎₁𝑁/𝑘₁, where 𝑁 is defined by (2.7).

Proof. Suppose 𝑋(𝑡)=(𝑥₁(𝑡),𝑥₂(𝑡),𝑦₁(𝑡),𝑦₂(𝑡))to be any positive solution of system (1.5) with initial conditions (1.6) and (1.7). It follows from the second equation of system (1.5) that for 𝑡≥𝑇₃+𝜏,𝑥̇₂(𝑡)≥𝛼𝑒^−𝑟1𝜏𝑥₂(𝑡−𝜏)−(𝑟₂+𝑎₁𝑁𝑘₁)𝑥₂(𝑡)−𝑟₃𝑥²₂(𝑡). (2.17) Consider the following auxiliary equation: 𝑢̇(𝑡)=𝛼𝑒^−𝑟1𝜏𝑢(𝑡−𝜏)−(𝑟₂+𝑎₁𝑁𝑘₁)𝑢(𝑡)−𝑟₃𝑢²(𝑡). (2.18) By Lemma 2.3, we obtain that lim_𝑡→+∞𝑢(𝑡)=𝛼𝑒^−𝑟1𝜏−𝑟₂−𝑎₁𝑁/𝑘₁𝑟₃. (2.19) According to comparison principle it follows that lim_𝑡→+∞inf𝑥₂(𝑡)≥𝛼𝑒^−𝑟1𝜏−𝑟₂−𝑎₁𝑁/𝑘₁𝑟₃. (2.20) Therefore, for sufficiently small 𝜀>0 there is a 𝑇₄=𝑇₃+𝜏 such that if 𝑡≥𝑇₄,𝑥₂(𝑡)≥𝛼𝑒^−𝑟1𝜏−𝑟₂−𝑎₁𝑁/𝑘₁𝑟₃−𝜀∶=𝑛₂. (2.21) By the third and forth equation of system (1.5), we have 𝑦̇₁(𝑡)≥(𝛽₁−𝑎₂𝑦₁(𝑡)𝑛₂+𝑘₂)𝑦₁(𝑡)+𝐷₁(𝑦₂(𝑡)−𝑦₁𝑦̇(𝑡)),₂(𝑡)=(𝛽₂−𝑟₄𝑦₂(𝑡))𝑦₂(𝑡)+𝐷₂(𝑦₁(𝑡)−𝑦₂(𝑡)),𝑡≥𝑇₄+𝜏. (2.22) Consider the following auxiliary equation: 𝑢̇₁(𝑡)=(𝛽₁−𝑎₂𝑢₁(𝑡)𝑛₂+𝑘₂)𝑢₁(𝑡)+𝐷₁(𝑢₂(𝑡)−𝑢₁𝑢̇(𝑡)),₂(𝑡)=(𝛽₂−𝑟₄𝑢₂(𝑡))𝑢₂(𝑡)+𝐷₂(𝑢₁(𝑡)−𝑢₂(𝑡)). (2.23) Define 𝑉₁₁(𝑡)=min{𝑢₁(𝑡),𝑢₂(𝑡)}. (2.24) Using a similar argument in the proof of [25, Lemma??2.1] we obtain lim_𝑡→+∞inf𝑉₁₁𝛽(𝑡)≥min{¹(𝑛₂+𝑘₂)𝑎₂,𝛽₂𝑟₄}∶=𝑛^∗₃. (2.25) Therefore, for sufficiently small 𝜀>0 there is a 𝑇₅=𝑇₄+𝜏 such that if 𝑡≥𝑇₅,𝑦₁(𝑡)≥𝑛^∗₃−𝜀∶=𝑛₃,𝑦₂(𝑡)≥𝑛^∗₃−𝜀∶=𝑛₃. (2.26) Setting 𝑇₆=𝑇₅+𝜏, then by (2.4), we have 𝑥₁(𝑡)≥𝛼𝑛₂(1−𝑒^−𝑟1𝜏)𝑟₁∶=𝑛₁,𝑡≥𝑇₆. (2.27) This completes the proof.

We now state a result on the extinction of the mature and immature prey.

Theorem 2.6. The mature and immature prey population will go to extinction if (𝐻2) holds(𝐻2)𝛼𝑒^−𝑟1𝜏<𝑟₂.

Remark 2.7. From the (𝐻2), we know that if the death rate of mature prey 𝑟₂ is more than the product of birth rate of immature prey 𝛼 and the surviving probability of each immature prey becomes mature 𝑒^−𝑟1𝜏, then the mature and immature prey population will go to extinction.

Proof. Suppose 𝑋(𝑡)=(𝑥₁(𝑡),𝑥₂(𝑡),𝑦₁(𝑡),𝑦₂(𝑡))to be any positive solution of system (1.5) with initial conditions (1.6) and (1.7). It follows from the second equation of system (1.5) that there is a 𝑇₁₁>0,𝑥̇₂(𝑡)≤𝛼𝑒^−𝑟1𝜏𝑥₂(𝑡−𝜏)−𝑟₂𝑥₂(𝑡)−𝑟₃𝑥²₂(𝑡). (2.28) Consider the following auxiliary equation: 𝑢̇(𝑡)=𝛼𝑒^−𝑟1𝜏𝑢(𝑡−𝜏)−𝑟₂𝑢(𝑡)−𝑟₃𝑢²(𝑡). (2.29) By Lemma 2.3, we derived from (2.29) and (𝐻2) that lim_𝑡→+∞𝑢(𝑡)=0. (2.30) A standard comparison argument shows that lim_𝑡→+∞𝑥₂(𝑡)=0. (2.31) Therefore, ∀𝜀>0, there is a 𝑇₇>𝑇₆ such that if 𝑡≥𝑇₇,0<𝑥₂(𝑡)<𝑟₁𝜀/2𝛼(1−𝑒^−𝑟1𝜏). Thus, we derive from (2.4) that for 𝑡≥𝑇₇+𝜏,𝑥₁(𝑡)≤𝛼^𝑡_𝑡−𝜏𝑒^{−𝑟1(𝑡−𝑠)}𝑟₁𝜀2𝛼(1−𝑒^−𝑟1𝜏)𝑑_𝑠<𝜀. (2.32) We therefore obtain that lim_𝑡→+∞𝑥₁(𝑡)=0. (2.33) This completes the proof.

3. Global Stability

In this section, we study the global asymptotic stability of a positive equilibrium of system (1.5). By Theorem 2.5 we see that if (𝐻1) satisfies, system (1.5) is uniformly persistent, which implies that system (1.5) must have at least one positive equilibrium. So in the following we assume that a positive equilibrium exists and denote it by 𝐸^∗(𝑥^∗₁,𝑥^∗₂,𝑦^∗₁,𝑦^∗₂).

Theorem 3.1. Let (𝐻1) hold. Assume further that(𝐻3)𝐴_𝑖>0,𝑖=1,3, where 𝐴₁=𝑎₂𝑁+𝑘₂+𝑎₁𝑥^∗₂4(𝑥^∗₂+𝑘₁)−𝑎₂𝑦^∗₁4(𝑛₂+𝑘₂)(𝑥^∗₂+𝑘₂),𝐴₃=𝑟₃𝑛₂+𝑎₁𝑘₁𝑛₃(𝑁+𝑘₁)(𝑥^∗₂+𝑘₁)+𝑎₁(𝑥^∗₂−𝑦^∗₁)𝑥^∗₂+𝑘₁−𝑎₂𝑦^∗₁(𝑛₂+𝑘₂)(𝑥^∗₂+𝑘₂), (3.1)where 𝑛₂=((𝛼𝑒^−𝑟1𝜏−𝑟₂)−𝑎₁𝑁/𝑘₁)/𝑟₃−𝜀,𝑛₃=min{𝛽₁(𝑛₂+𝑘₂)/𝑎₂,𝛽₂/𝑟₄}−𝜀,𝜀>0 is a sufficient small constant, and 𝑁 is defined by (2.7).
Then the positive equilibrium 𝐸^∗(𝑥^∗₁,𝑥^∗₂,𝑦^∗₁,𝑦^∗₂) of system (1.5) is globally asymptotically stable.

Remark 3.2. Theorem 3.1 shows that if the time delay due to maturity is sufficiently small, the positive equilibrium of system (1.5) is globally asymptotically stable.

Proof. We first consider the following subsystem: 𝑥̇₂(𝑡)=𝛼𝑒^−𝑟1𝜏𝑥₂(𝑡−𝜏)−𝑟₂𝑥₂(𝑡)−𝑟₃𝑥²₂𝑎(𝑡)−¹𝑦₁(𝑡)𝑥₂(𝑡)𝑥₂(𝑡)+𝑘₁,𝑦̇₁(𝑡)=(𝛽₁−𝑎₂𝑦₁(𝑡)𝑥₂(𝑡)+𝑘₂)𝑦₁(𝑡)+𝐷₁(𝑦₂(𝑡)−𝑦₁𝑦̇(𝑡)),₂(𝑡)=(𝛽₂−𝑟₄𝑦₂(𝑡))𝑦₂(𝑡)+𝐷₂(𝑦₁(𝑡)−𝑦₂(𝑡)). (3.2) Noting that 𝐸^∗(𝑥^∗₂,𝑦^∗₁,𝑦^∗₂) is a positive equilibrium of system (3.2), we can rewrite system (3.2) as 𝑥̇₂(𝑡)=𝛼𝑒^−𝑟1𝜏(𝑥₂(𝑡−𝜏)−𝑥^∗₂)−𝑟₂(𝑥₂(𝑡)−𝑥^∗₂)−𝑟₃(𝑥₂(𝑡)+𝑥^∗₂)(𝑥₂(𝑡)−𝑥^∗₂)−𝑎₁𝑘₁𝑦₁(𝑡)(𝑥₂(𝑡)−𝑥^∗₂)(𝑥₂(𝑡)+𝑘₁)(𝑥^∗₂+𝑘₁)−𝑎₁𝑥^∗₂(𝑦₁(𝑡)−𝑦^∗₁)𝑥^∗₂+𝑘₁,𝑦̇₁𝑎(𝑡)=(−²(𝑦₁(𝑡)−𝑦^∗₁)𝑥₂(𝑡)+𝑘₂+𝑎₂(𝑦₁(𝑡)−𝑦^∗₁)(𝑥₂(𝑡)+𝑘₂)(𝑥^∗₂+𝑘₂))𝑦₁−𝐷(𝑡)₁𝑦^∗₁𝑦₂(𝑡)(𝑦₁(𝑡)−𝑦^∗₁𝐷)+¹𝑦^∗₁𝑦₁(𝑡)(𝑦₂(𝑡)−𝑦^∗₂𝑦̇),₂(𝑡)=(−𝑟₄(𝑦₂(𝑡)−𝑦^∗₂))𝑦₂𝐷(𝑡)−²𝑦^∗₂𝑦₁(𝑡)(𝑦₂(𝑡)−𝑦^∗₂𝐷)+²𝑦^∗₂𝑦₂(𝑡)(𝑦₁(𝑡)−𝑦^∗₁). (3.3) Define 𝑉₁(𝑡)=²_𝑖=1𝑐_𝑖(𝑦_𝑖(𝑡)−𝑦^∗_𝑖−𝑦^∗_𝑖𝑦ln^𝑖(𝑡)𝑦^∗_𝑖1)+2𝑐₃(𝑥₂(𝑡)−𝑥^∗₂)². (3.4) Calculating the derivative of 𝑉₁(𝑡) along solution of system (1.5), we have 𝑑𝑉₁(𝑡)𝑑_𝑡=²_𝑖=1𝑐_𝑖(𝑦_𝑖(𝑡)−𝑦^∗_𝑖)𝑦̇_𝑖(𝑡)𝑦_𝑖(𝑡)+𝑐₃(𝑥₂(𝑡)−𝑥^∗₂)𝑥̇₂(𝑡)=−𝑐₁𝑎₂(𝑦₁(𝑡)−𝑦^∗₁)²𝑥₂(𝑡)+𝑘₂+𝑐₁𝑎₂𝑦^∗₁(𝑥₂(𝑡)−𝑥^∗₂)(𝑦₁(𝑡)−𝑦^∗₁)(𝑥₂(𝑡)+𝑘₂)(𝑥^∗₂+𝑘₂)−𝑐₁𝐷₁𝑦^∗₁𝑦₁𝑦(𝑡)²(𝑡)(𝑦₁(𝑡)−𝑦^∗₁)²+𝑐₁𝐷₁𝑦^∗₁(𝑦₁(𝑡)−𝑦^∗₁)(𝑦₂(𝑡)−𝑦^∗₂)−𝑐₂𝑟₄(𝑦₂(𝑡)−𝑦^∗₂)²−𝑐₂𝐷₂𝑦^∗₂𝑦₂𝑦(𝑡)¹(𝑡)(𝑦₂(𝑡)−𝑦^∗₂)²+𝑐₂𝐷₂𝑦^∗₂(𝑦₁(𝑡)−𝑦^∗₁)(𝑦₂(𝑡)−𝑦^∗₂)+𝑐₃𝛼𝑒^−𝑟1𝜏(𝑥₂(𝑡−𝜏)−𝑥^∗₂)(𝑥₂(𝑡)−𝑥^∗₂)−𝑐₃𝑟₂(𝑥₂(𝑡)−𝑥^∗₂)²−𝑐₃𝑟₃(𝑥₂(𝑡)+𝑥^∗₂)(𝑥₂(𝑡)−𝑥^∗₂)²−𝑐₃𝑎₁𝑘₁𝑦₁(𝑡)(𝑥₂(𝑡)−𝑥^∗₂)²(𝑥₂(𝑡)+𝑘₁)(𝑥^∗₂+𝑘₁)−𝑐₃𝑎₁𝑥^∗₂(𝑥₂(𝑡)−𝑥^∗₂)(𝑦₁(𝑡)−𝑦^∗₁)𝑥^∗₂+𝑘₁. (3.5) Setting 𝑐₁=1,𝑐₂=𝐷₁𝑦^∗₂/𝐷₂𝑦^∗₁. By (3.5) we obtain 𝑑𝑉₁(𝑡)𝑑_𝑡𝑎=−²(𝑦₁(𝑡)−𝑦^∗₁)²𝑥₂(𝑡)+𝑘₂−𝐷₁𝑦^∗₂𝐷₂𝑦^∗₁𝑟₄(𝑦₂(𝑡)−𝑦^∗₂)²−𝐷₁𝑦^∗₁[𝑦₂(𝑡)𝑦₁(𝑡)(𝑦₁(𝑡)−𝑦^∗₁)−𝑦₁(𝑡)𝑦₂(𝑡)(𝑦₂(𝑡)−𝑦^∗₂)]²+𝑎₂𝑦^∗₁(𝑥₂(𝑡)−𝑥^∗₂)(𝑦₁(𝑡)−𝑦^∗₁)(𝑥₂(𝑡)+𝑘₂)(𝑥^∗₂+𝑘₂)−𝑐₃𝑎₁𝑥^∗₂(𝑥₂(𝑡)−𝑥^∗₂)(𝑦₁(𝑡)−𝑦^∗₁)𝑥^∗₂+𝑘₁+𝑐₃𝛼𝑒^−𝑟1𝜏(𝑥₂(𝑡−𝜏)−𝑥^∗₂)(𝑥₂(𝑡)−𝑥^∗₂)−𝑐₃𝑟₂(𝑥₂(𝑡)−𝑥^∗₂)²−𝑐₃𝑟₃(𝑥₂(𝑡)+𝑥^∗₂)(𝑥₂(𝑡)−𝑥^∗₂)²−𝑐₃𝑎₁𝑘₁𝑦₁(𝑡)(𝑥₂(𝑡)−𝑥^∗₂)²(𝑥₂(𝑡)+𝑘₁)(𝑥^∗₂+𝑘₁). (3.6) Using the inequality 𝑎𝑏≤(1/2)𝑘𝑎²+(1/2𝑘)𝑏², it follows from (3.6) that 𝑑𝑉₁(𝑡)𝑑_𝑡𝑎≤−²𝑥₂(𝑡)+𝑘₂(𝑦₁(𝑡)−𝑦^∗₁)²−𝐷₁𝑦^∗₂𝐷₂𝑦^∗₁𝑟₄(𝑦₂(𝑡)−𝑦^∗₂)²𝑎+(²𝑦^∗₁(𝑥₂(𝑡)+𝑘₂)(𝑥^∗₂+𝑘₂)−𝑐₃𝑎₁𝑥^∗₂𝑥^∗₂+𝑘₁)(𝐴(𝑥₂(𝑡)−𝑥^∗₂)²2+(𝑦₁(𝑡)−𝑦^∗₁)²)2𝐴+𝑐₃𝛼𝑒^−𝑟1𝜏(𝐵(𝑥₂(𝑡)−𝑥^∗₂)²2+(𝑥₂(𝑡−𝜏)−𝑥^∗₂)²2𝐵)−𝑐₃𝑟₂(𝑥₂(𝑡)−𝑥^∗₂)²−𝑐₃𝑟₃(𝑥₂(𝑡)+𝑥^∗₂)(𝑥₂(𝑡)−𝑥^∗₂)²−𝑐₃𝑎₁𝑘₁𝑦₁(𝑡)(𝑥₂(𝑡)+𝑘₁)(𝑥^∗₂+𝑘₁)(𝑥₂(𝑡)−𝑥^∗₂)², (3.7) where parameters 𝐴,𝐵 are positive constants to be determined.
Define 𝑉(𝑡)=𝑉₁1(𝑡)+𝑐2𝐵³𝛼𝑒^{−𝑟1𝜏𝑡}_𝑡−𝜏(𝑥₂(𝑠)−𝑥^∗₂)²𝑑_𝑠. (3.8) Setting 𝐴=2,𝐵=1,𝑐₃=1, then it follows from (3.7) and (3.8) that 𝑑𝑉₁(𝑡)𝑑_𝑡𝑎≤−{²𝑥₂(𝑡)+𝑘₂+𝑎₁𝑥^∗₂4(𝑥^∗₂+𝑘₁)−𝑎₂𝑦^∗₁4(𝑥₂(𝑡)+𝑘₂)(𝑥^∗₂+𝑘₂)}(𝑦₁(𝑡)−𝑦^∗₁)²−𝐷₁𝑦^∗₂𝐷₂𝑦^∗₁𝑟₄(𝑦₂(𝑡)−𝑦^∗₂)²−{𝑟₃𝑥₂𝑎(𝑡)+¹𝑘₁𝑦₁(𝑡)(𝑥₂(𝑡)+𝑘₁)(𝑥^∗₂+𝑘₁)+𝑎₁(𝑥^∗₂−𝑦^∗₁)𝑥^∗₂+𝑘₁−𝑎₂𝑦^∗₁(𝑥₂(𝑡)+𝑘₂)(𝑥^∗₂+𝑘₂)}(𝑥₂(𝑡)−𝑥^∗₂)²𝑎≤−{²𝑁+𝑘₂+𝑎₁𝑥^∗₂4(𝑥^∗₂+𝑘₁)−𝑎₂𝑦^∗₁4(𝑛₂+𝑘₂)(𝑥^∗₂+𝑘₂)}(𝑦₁(𝑡)−𝑦^∗₁)²−𝐷₁𝑦^∗₂𝐷₂𝑦^∗₁𝑟₄(𝑦₂(𝑡)−𝑦^∗₂)²−{𝑟₃𝑛₂+𝑎₁𝑘₁𝑛₃(𝑁+𝑘₁)(𝑥^∗₂+𝑘₁)+𝑎₁(𝑥^∗₂−𝑦^∗₁)𝑥^∗₂+𝑘₁−𝑎₂𝑦^∗₁(𝑛₂+𝑘₂)(𝑥^∗₂+𝑘₂)}(𝑥₂(𝑡)−𝑥^∗₂)²∶=−𝐴₁(𝑦₁(𝑡)−𝑦^∗₁)²−𝐴₂(𝑦₂(𝑡)−𝑦^∗₂)²−𝐴₃(𝑥₂(𝑡)−𝑥^∗₂)², (3.9) where 𝐴₁=𝑎₂𝑁+𝑘₂+𝑎₁𝑥^∗₂4(𝑥^∗₂+𝑘₁)−𝑎₂𝑦^∗₁4(𝑛₂+𝑘₂)(𝑥^∗₂+𝑘₂),𝐴₂=𝐷₁𝑦^∗₂𝐷₂𝑦^∗₁𝑟₄,𝐴₃=𝑟₃𝑛₂+𝑎₁𝑘₁𝑛₃(𝑁+𝑘₁)(𝑥^∗₂+𝑘₁)+𝑎₁(𝑥^∗₂−𝑦^∗₁)𝑥^∗₂+𝑘₁−𝑎₂𝑦^∗₁(𝑛₂+𝑘₂)(𝑥^∗₂+𝑘₂). (3.10)𝑁,𝑛₂, and 𝑛₃ are defined in (2.16), (2.21), and (2.26), respectively.
If (𝐻1) and (𝐻3) hold and 𝜀>0 is sufficiently small, we have 𝐴_𝑖>0,𝑖=1,3. In view of Lyapunov theorem [26], we conclude that the positive equilibrium 𝐸^∗(𝑥^∗₂,𝑦^∗₁,𝑦^∗₂) of system (3.2) is globally asymptotically stable. Thus, we have lim_𝑡→+∞𝑥₂(𝑡)=𝑥^∗₂,lim_𝑡→+∞𝑦₁(𝑡)=𝑦^∗₁,lim_𝑡→+∞𝑦₂(𝑡)=𝑦^∗₂. (3.11) Using L'Hospital's rule, it follows from (2.4) and (3.11) that lim_𝑡→+∞𝑥₁(𝑡)=lim_𝑡→+∞𝛼^𝑡_𝑡−𝜏𝑒^{−𝑟1(𝑡−𝑠)}𝑥₂(𝑠)𝑑_𝑠=lim_𝑡→+∞𝛼𝑟₁𝑥₂(𝑡)−𝑒^−𝑟1𝜏𝑥₂=(𝑡−𝜏)𝛼𝑥^∗₂𝑟₁(1−𝑒^−𝑟1𝜏)=𝑥^∗₁. (3.12) This completes the proof.

It is interesting to discuss the local stability of the positive equilibrium 𝐸^∗(𝑥^∗₁,𝑥^∗₂,𝑦^∗₁,𝑦^∗₂) of system (1.5).

The characteristic equation of the positive equilibrium 𝐸^∗ of system (1.5) is of the form (𝜆+𝑟₁)[𝑃(𝜆)+𝑄(𝜆)𝑒^−𝜆𝜏]=0, (3.13) where 𝑃(𝜆)=𝜆³+𝑎₂𝜆²+𝑎₁𝜆+𝑎₀,𝑄(𝜆)=𝑏₂𝜆²+𝑏₁𝜆+𝑏₀, (3.14) here 𝑎₀=(𝑟₂+2𝑟₃+𝑎₁𝑘₁𝑦^∗₁(𝑥^∗₂+𝑘₁)²𝑎)[(²𝑦^∗₁𝑥^∗₂+𝑘₂+𝐷₁𝑦^∗₂𝑦^∗₁)(𝑟₄𝑦^∗₂+𝐷₂𝑦^∗₁𝑦^∗₂)−𝐷₁𝐷₂]+𝑎₁𝑎₂𝑥^∗₂𝑦₁^∗2(𝑥^∗₂+𝑘₁)(𝑥^∗₂+𝑘₂)²(𝑟₄𝑦^∗₂+𝐷₂𝑦^∗₁𝑦^∗₂𝑎),₁=(𝑟₂+2𝑟₃+𝑎₁𝑘₁𝑦^∗₁(𝑥^∗₂+𝑘₁)²𝑎)(²𝑦^∗₁𝑥^∗₂+𝑘₂+𝐷₁𝑦^∗₂𝑦^∗₁+𝑟₄𝑦^∗₂+𝐷₂𝑦^∗₁𝑦^∗₂)𝑎+[(²𝑦^∗₁𝑥^∗₂+𝑘₂+𝐷₁𝑦^∗₂𝑦^∗₁)(𝑟₄𝑦^∗₂+𝐷₂𝑦^∗₁𝑦^∗₂)−𝐷₁𝐷₂𝑎]+¹𝑎₂𝑥^∗₂𝑦₁^∗2(𝑥^∗₂+𝑘₁)(𝑥^∗₂+𝑘₂)²,𝑎₂=𝑟₂+2𝑟₃+𝑎₁𝑘₁𝑦^∗₁(𝑥^∗₂+𝑘₁)²+𝑎₂𝑦^∗₁𝑥^∗₂+𝑘₂+𝐷₁𝑦^∗₂𝑦^∗₁+𝑟₄𝑦^∗₂+𝐷₂𝑦^∗₁𝑦^∗₂,𝑏₀=−𝛼𝑒^−𝑟1𝜏𝑎[(²𝑦^∗₁𝑥^∗₂+𝑘₂+𝐷₁𝑦^∗₂𝑦^∗₁)(𝑟₄𝑦^∗₂+𝐷₂𝑦^∗₁𝑦^∗₂)−𝐷₁𝐷₂𝑏],₁=−𝛼𝑒^−𝑟1𝜏(𝑎₂𝑦^∗₁𝑥^∗₂+𝑘₂+𝐷₁𝑦^∗₂𝑦^∗₁+𝑟₄𝑦^∗₂+𝐷₂𝑦^∗₁𝑦^∗₂𝑏),₂=−𝛼𝑒^−𝑟1𝜏. (3.15) Clearly, 𝜆=−𝑟₁ is a negative eigenvalue. If 𝑟₃𝑥^∗₂−𝑎₁𝑥^∗₂𝑦^∗₁/(𝑥^∗₂+𝑘₁)²>0, which implies that 𝑎_𝑖+𝑏_𝑖>0(𝑖=1,2,3),and (𝑎₁+𝑏₁)(𝑎₂+𝑏₂)−(𝑎₀+𝑏₀)>0, then by Routh-Hurwitz Theorem the positive equilibrium 𝐸^∗ of system (1.5) is locally asymptotically stable when 𝜏=0.

Let 𝐹(𝑦)=|𝑃(𝑖𝑦)|²−|𝑄(𝑖𝑦)|²=𝑦⁶+𝑙𝑦⁴+𝑚𝑦²+𝑛=0, (3.16) where 𝑙=𝑎²₂−2𝑎₁−𝑏²₂,𝑚=𝑎²₁−2𝑎₀𝑎₂+2𝑏₀𝑏₂−𝑏²₁,𝑛=𝑎²₀−𝑏²₀. (3.17) Let 𝑧=𝑦², and then (3.16) becomes 𝑧³+𝑙𝑧²+𝑚𝑧+𝑛𝑧=0. (3.18) By applying the results on the distribution of roots of (3.16) and (3.18) in [27] and [26, Theorem??4.1, page 83], we therefore derive the following results on the stability of the positive equilibrium 𝐸^∗.

Theorem 3.3. Suppose that system (1.5) admits a positive equilibrium 𝐸^∗(𝑥^∗₁,𝑥^∗₂,𝑦^∗₁,𝑦^∗₂) and 𝑟₃𝑥^∗₂−𝑎₁𝑥^∗₂𝑦^∗₁/(𝑥^∗₂+𝑘₁)²>0.(1)If Δ=𝑙²−3𝑚≤0, then the positive equilibrium 𝐸^∗ of system (1.5) is locally asymptotically stable.(2)If Δ=𝑙²−3𝑚>0 and ℎ(𝑧^∗₁)≤0, then there exists a positive number 𝜏₀ such that the positive equilibrium 𝐸^∗ of system (1.5) is locally asymptotically stable if 0<𝜏<𝜏₀ and is locally unstable if 𝜏>𝜏₀; further, as 𝜏 increases through 𝜏₀, 𝐸^∗ bifurcates into small amplitude periodic solutions, here, 𝑧^∗₁=(−𝑙+√Δ)/3,ℎ(𝑧)=𝑧³+𝑙𝑧²+𝑚𝑧+𝑛𝑧.

4. Two Examples

In this section, we give two examples to illustrate our main results.

Example 4.1. Consider the following system: 𝑥̇₁(𝑡)=5𝑥₂(𝑡)−𝑥₁(𝑡)−5𝑒^−𝜏𝑥₂𝑥̇(𝑡−𝜏),₂(𝑡)=5𝑒^−𝜏𝑥₂(𝑡−𝜏)−1.5𝑥₂(𝑡)−3𝑥²₂(𝑡)−0.8𝑦₁(𝑡)𝑥₂(𝑡)𝑥₂,𝑦̇(𝑡)+8₁(𝑡)=(0.2−1.5𝑦₁(𝑡)𝑥₂(𝑡)+1.5)𝑦₁(𝑡)+0.5(𝑦₂(𝑡)−𝑦₁𝑦̇(𝑡)),₂(𝑡)=(1.5−𝑦₂(𝑡))𝑦₂(𝑡)+0.5(𝑦₁(𝑡)−𝑦₂(𝑡)), (4.1) where the parameter 𝜏 is a positive constant.

System (4.1) has a unique positive equilibrium 𝐸^∗(0.9589,0.4874,0.7466,1.2895). It is easy to show that if 𝜏<0.8973, then (𝐻1) and (𝐻3) hold for system (4.1). By Theorem 2.5 we see that system (4.1) is uniformly persistent when 𝜏<0.8973. By Theorem 3.1 we see that the positive equilibrium of system (4.1) is globally asymptotically stable when 𝜏=0.5. Numerical integration can be carried out using standard MATLAB algorithm. Numerical simulation also confirms the fact (see Figure 1).

Figure 1 

The temporal solution found by numerical integration of system (4.1) with 𝜏=0.5 and (Φ1(𝜃),Φ2(𝜃),Ψ1(𝜃),Ψ2(𝜃))=(5(1−𝑒−0.5),1,0.6,0.6).

Example 4.2. Consider the following system: 𝑥̇₁(𝑡)=5𝑥₂(𝑡)−𝑥₁(𝑡)−5𝑒⁻¹𝑥₂𝑥̇(𝑡−1),₂(𝑡)=5𝑒⁻¹𝑥₂(𝑡−1)−2𝑥₂(𝑡)−3𝑥²₂(𝑡)−2𝑦₁(𝑡)𝑥₂(𝑡)𝑥₂,𝑦̇(𝑡)+8₁(𝑡)=(1−2𝑦₁(𝑡)𝑥₂(𝑡)+2)𝑦₁(𝑡)+𝑦₂(𝑡)−𝑦₁𝑦̇(𝑡),₂(𝑡)=(1−𝑦₂(𝑡))𝑦₂(𝑡)+𝑦₁(𝑡)−𝑦₂(𝑡). (4.2)
System (4.2) has a unique boundary equilibrium 𝐸^∗(0,0,1,1). It is easy to show that (𝐻2)holds for system (4.2). By Theorem 2.6 we see that mature and immature prey population goes to extinction. Numerical integration can be carried out using standard MATLAB algorithm. Numerical simulation also confirms the fact (see Figure 2).

Figure 2 

The temporal solution found by numerical integration of system (4.2) with (Φ1(𝜃),Φ2(𝜃),Ψ1(𝜃),Ψ2(𝜃))=(5(1−𝑒−1),2,2,2).

5. Discussion

In this paper, we discussed a generalized Leslie-Gower-type predator-prey model with stage structure for prey and predator dispersal in two-patch environments. By using comparison arguments we established sufficient conditions for system (1.5) to be permanent. By constructing Lyapunov functionals, sufficient conditions are derived for the global asymptotic stability of the positive equilibrium of system (1.5). By Theorem 3.1 we see that if the birth rate of immature prey and the extent to which environment provides protection to mature prey and predator in patch 1, respectively, are high and the maximum value of the per-capita reduction rate of mature prey due to predator in patch 1 is low satisfying (H 1) and (H 3), the positive equilibrium of system (1.5) is globally asymptotically stable. By Theorem 2.6 we see that if the death rate of mature prey is more than the transformation rate of immatures to matures satisfying (H2), the immature and mature prey population will go to extinction.

Acknowledgments

The first author was partially supported by the Key Project of Chinese Ministry of Education (209131), the Project Sponsored by the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry, the NSF of Bureau of Education of Gansu Province of China for Postgraduate Tutors (0803-01), the Development Program for Outstanding Young Teachers in Lanzhou University of Technology (Q200703), and the Doctor's Foundation of Lanzhou University of Technology. The second author was partially supported by the Young Teacher's Foundation of Dali University (2008X34).