# Dynamic Analysis of a Heterogeneous Diffusive Prey-Predator System in Time-Periodic Environment

Authors: Chuanjun Dai; He Liu; Zhan Jin; Qing Guo; Yi Wang; Hengguo Yu; Qi Wang; Zengling Ma; Min Zhao
Journal: Complexity
 (2020)

---

## Abstract

In this paper, a heterogeneous diffusive prey-predator system is first proposed and then studied analytically and numerically. Some sufficient conditions are derived, including permanence and extinction of system and the boundedness of the solution. The existence of periodic solution and its stability are discussed as well. Furthermore, numerical results indicate that both the spatial heterogeneity and the time-periodic environment can influence the permanence and extinction of the system directly. Our numerical results are consistent with the analytical analysis.


---

## Body

## 1. Introduction

Due to the complexity of ecosystems, prey-predator dynamics have always drawn interest among mathematical ecologists, as well as experimental ecologists [1–3]. The significance of studying prey-predator dynamics is to gain insights into the complex ecological processes. Prey-predator models, as the base of researching prey-predator dynamics, have attracted increasing attention [4–7]. Since Holling [8] introduced the concept of the functional response, a lot of studies have been devoted to the understanding of the effect of functional response on prey-predator dynamics [9]. Usually, the functional response is assumed to be either prey dependent or ratio dependent in prey-predator models [10, 11].A classical general prey-predator system can be written as follows [12]:(1a)dNdt=fNN−gN,PP,(1b)dPdt=hgN,P,PP,where N and P denote the prey and predator densities, respectively, fN is the prey growth rate, gN,P is the functional response, and hgN,P,P is the per capita growth rate of predators. Let hgN,P,P=egN,P−mP, then equation (1b) can be rewritten as follows:(2)dPdt=egN,P−mPP,where e is the conversion efficiency and mP is the specific mortality of predators in absence of prey. For the function mP, the most widely accepted assumption [13] is mP=μ, where μ is a constant describing the death rate of the predator. However, Cavani and Farkas [14] introduced another function for mP:(3)mP=γ+δP1+P,where γ is the mortality at low density and δ is the limiting, maximal mortality (obviously, γ<δ). The specific mortality (3) depends on the quantity of predators, which suggests that the predator mortality is neither a constant nor an unbounded function, and increasing with quantity. Obviously, when γ=δ, equation (3) can be simplified to a constant death rate type. Prey-predator systems with this nonconstant death rate have been studied by some researchers [15–17].Additionally, in order to understand patterns and the mechanisms of spatial distribution of interacting species, the dispersal process is taken into consideration [18–20]. Thus, the spatiotemporal dynamics of a prey-predator system can be presented by a couple of reaction-diffusion equations based on equations (1a) and (2) [10, 21, 22]:(4a)∂Ndt=fNN−gN,PP+DNΔN,(4b)∂Pdt=egN,P−mPP+DPΔP,where DN and DP are the prey and predator diffusion coefficients, respectively, and the Laplace operator Δ describes the spatial dispersal.Because of the emergence of Lotka–Volterra models [23, 24], a logistic type growth fN is usually assumed for the prey species in the models. Some functional response gN,P are taken into account in many works, such as Holling type [25], Michaelis–Menten type [26, 27], and Beddington–DeAngelis type [28, 29]. Especially, many biologists argued that the ratio-dependent theory is more suitable for describing prey-predator systems in many situations [13, 30–32]. Since Ardini and Ginzburg proposed the ratio-dependent prey-predator system, the prey-predator systems with ratio-dependent functional response are widely studied [13, 33–36], and many interesting results are obtained.Based on model (4a) and (4b), in this paper, we employ the ratio-dependent functional response and the nonconstant death rate (i.e., equation (3)) and assume that the growth rate of prey population follows the logistic growth type. Moreover, let u and v be the prey density and the predator density, respectively. Then, the resulting system is(5a)∂u∂t=ru1−uK−auvbv+u+μ1Δu,x∈Ω,t>0,(5b)∂v∂t=euvbv+u−c+δv1+vv+μ2Δv,x∈Ω,t>0,(5c)∂u∂n=∂v∂n=0,x∈∂Ω,t>0,where Ω∈ℝn is a bounded domain with smooth boundary ∂Ω.In system (5a)–(5c), when c=δ=μ, the system without diffusion is so-called the Michaelis–Menten ratio-dependent predator-prey system, which has been studied by many researchers. Kuang and Beretta [37] systematically studied the global behaviors of solutions and obtained some new and significant results, but many important open questions remain to be unsolved. For these open questions, Hsu et al. [38] resolved the global stability of all equilibria in various cases and the uniqueness of limit cycles by transforming the Michaelis–Menten-type ratio-dependent model. Xiao and Ruan [39] investigated the qualitative behavior of the Michaelis–Menten-type ratio-dependent model at the origin in the interior of the first quadrant and confirmed that the origin is indeed a critical point inducing rich and complicated dynamics. Additionally, when the diffusion process is considered, the Michaelis–Menten ratio-dependent predator-prey system with diffusion can produce rich spatial patterns, which makes it a widely studied system for pattern formation [10, 40–43].Whilec<δ, Kovács et al. [44] incorporated delays into system (5a)–(5c) and studied the qualitative behaviour of the system without diffusion. Yun et al. [45] presented an efficient and accurate numerical method for solving system (5a)–(5c) with a Turing instability and studied the existence of nonconstant stationary solutions. Aly et al. [46] studied Turing instability for system (5a)–(5c) and showed that diffusion-driven instability occurs at a certain critical value analytically. In these works, parameters in system (5a)–(5c) are always considered as constants.However, it seems that there is no research for considering spatial heterogeneity and time-periodic environment in system (5a)–(5c). It is well known that spatial heterogeneity occurs at all scales of the environment [47]. Additionally, interactive populations often live in a fluctuating environment [48], where some environmental conditions such as temperature, light, availability of food, and other resources usually vary in time. Specially, some data depending on season in systems may be periodic functions of time. Thus, more realistic models to describe ecosystem should be nonautonomous systems with spatial heterogeneity. With this mind, we propose the following system to study effects of spatial heterogeneity and time-periodic environment on prey-predator dynamics:(6a)∂ut,x∂t=rt,xut,x1−ut,xKt,x−at,xut,xvt,xbt,xvt,x+ut,x+μ1Δut,x,(6b)∂vt,x∂t=et,xut,xvt,xbt,xvt,x+ut,x−ct,x+δt,xvt,x1+vt,xvt,x+μ2Δvt,x,(6c)∂ut,x∂n=∂vt,x∂n=0,x∈∂Ω,t>0,where ut,x and vt,x represent the densities of the prey and predator, respectively, at a space point x and time t; for simplification, ut,x and vt,x are rewritten as u and v in the rest of this paper, respectively; rt,x is the intrinsic growth rate of prey population; Kt,x denotes the environmental carrying capacity of prey population; at,x is the capturing rate of the predator; bt,x is the half saturation; and et,x denotes the conversion rate. The term ct,x+δt,xvt,x/1+vt,x describes the specific mortality of predators in absence of prey population, where ct,x is the mortality at low density and δt,x is the limiting, maximal mortality. The terms μ1Δut,x and μ2Δvt,x with positive diffusion coefficients μ1 and μ2 represent the nonhomogeneous dispersion of the prey and the predator, respectively. Neumann boundary conditions (see equation (6c)) are employed, which characterize the absence of migration. Here, we assume that prey and predator populations are confined to a fixed bounded space domain Ω∈Rn with smooth boundary ∂Ω and Ω¯=Ω∪∂Ω.The rest of the paper is organized as follows. In Section2, some conditions and definitions are given. In Section 3, dynamics of system (6a)–(6c) are studied, including boundedness, permanence, extinction, and periodic solution. Moreover, a series of numerical simulations are carried out for further study of the dynamics of system (6a)–(6c) in Section 4. Finally, the paper ends with conclusion in Section 5.


## 2. Preliminaries

Letℝ,ℤ, and ℕ be the sets of all real numbers, integers, and positive integers, repectively, and ℝ+=0,+∞. We assume that the following condition holds throughout the paper:(H) The functions rt,x,Kt,x,at,x,bt,x,ct,x,et,x,δt,x are bounded positive-valued functions on ℝ×Ω¯, continuously differentiable in t and x, and are periodic in t with period τ>0.Moreover, for a continuous functionϕt,x, we denote ϕL=inft,xϕt,x and ϕM=supt,xϕt,x.Definition 1.
Solutions of system (6a)–(6c) are ultimately bounded if there exist positive constants N1andN2 such that for every solution ut,x,u0,v0,vt,x,u0,v0, there exists a moment of time T=Tu0,v0>0 such that(7)ut,x,u0,v0≤N1,vt,x,u0,v0≤N2,for all x∈Ω¯ and t≥T.Definition 2.
System (6a)–(6c) is permanent if there exist positive constants ζ and η such that for every solution with nonnegative initial functions u0x≢0 and v0x≢0, there exists a moment of time t^=t^u0,v0 such that(8)ζ≤ut,x,u0,v0≤η,ζ≤vt,x,u0,v0≤η,for all x∈Ω¯ and t≥t^.
Consider the following equations:(9a)∂u∂t−dΔu+ft,x,u=0,t,x∈0,T×Ω,(9b)∂u∂n=0,t,x∈0,T×∂Ω.
Then, we have the following definition.Definition 3.
A functionu^:0,T×Ω⟶ℝ is called a lower solution of equations (9a) and 9b if it satisfies(10a)∂u^∂t−dΔu^+ft,x,u^≤0,t,x∈0,T×Ω,(10b)∂u^∂n≤0,t,x∈0,T×∂Ω.
To analyze dynamics of system (6a)–(6c), the following results will be needed.Theorem 1 (Walter [49]).
Suppose that vector-functionsvt,x=v1t,x,…,vmt,x and wt,x=w1t,x,…,wmt,x,m≥1, satisfy the following conditions:(i)
They are of classC2 in x, x∈Ω and of class C1 in t,x∈a,b×Ω¯, where Ω∈ℝn is a bounded domain with a smooth boundary;(ii)
vt−μΔv−gt,x,v≤wt−μΔw−gt,x,w, where t,x∈a,b×Ω,μ=μ1,…,μm>0 (inequalities between vectors are satisfied coordinate-wise), and vector function gt,x,u=g1t,x,u,…,gmt,x,u is continuously differentiable and quasimonotonically increasing with respect to u=u1,…,um:(11)∂git,x,u1,…,um∂uj≥0,i,j=1,…,m,i≠j;(iii)
∂v/∂n=∂w/∂n=0,t,x∈a,b×∂Ω.
Then,vt,x≤wt,x for t,x∈a,b×Ω¯.Theorem 2 (Smith [50]).
Assume thatT and d are positive real numbers, a function ut,x is continuous on 0,T×Ω¯, continuously differentiable in x∈Ω¯, with continuous derivatives ∂2u/∂xi∂xj and ∂u/∂t on 0,T×Ω, and ut,x satisfies the following inequalities:(12)∂u∂t−dΔu+ct,xu≥0,t,x∈0,T×Ω,∂u∂n≥0,t,x∈0,T×∂Ω,u0,x≥0,x∈Ω,where ct,x is bounded on 0,T×Ω. Then, ut,x≥0 on 0,T×Ω¯.Moreover,ut,x is strictly positive on 0,T×Ω¯ if ut,x is not identically zero.


## 3. Main Results

### 3.1. Boundedness

From the biological and ecological viewpoint, we are always interested in the nonnegative solutions. Thus, the following theorem is given first in system (6a)–(6c).Theorem 3.
Suppose that the condition (H) holds, then nonnegative and positive quadrants of ℝ2 are positively invariant for system (6a)–(6c).Proof.
Letut,x,u0,v0,vt,x,u0,v0 be a solution of system (6a)–(6c) with initial condition u0x≥0≢0,v0x≥0≢0. Additionally, u^ is a solution of the following system:(13)∂u^∂t−μ1Δu^−u^rL−aMbL−rMKLu^=0,u^0,x=u0x.
From system (6a), we can obtain(14)∂u∂t−μ1Δu−rt,xu1−uKt,x+at,xuvbt,xv+u≤∂u∂t−μ1Δu−urL−aMbL−rMKLu,which implies u^t,x is a lower solution of system (6a). According to Theorem 2, it is obvious that u^t,x≥0 for all x∈Ω¯ and t>0. Furthermore, due to u0x≥0≢0, u^t,x>0 holds for all x∈Ω¯ and t>0. Thus, ut,x>0 holds because ut,x is bounded from below by positive function u^t,x.
For system (6b), it can be simply verified that v^t,x is a lower solution of system (6b), where v^t,x satisfies(15)∂v^∂t−μ2Δv^+δMv^=0,v^0,x=v0x.
By the similar argument tout,x, we can prove the positiveness of vt,x.
This completes the proof.
Based on Theorem3, we will discuss ultimate boundedness of solutions in system (6a)–(6c), and then the following theorem can be obtained.Theorem 4.
If the condition (H) holds, then all solutions of system (6a)–(6c) with nonnegative initial conditions are ultimately bounded.Proof.
From system (6a), it can be found that the following inequality holds:(16)0=∂u∂t−μ1Δu−rt,xu1−uKt,x+at,xuvbt,xv+u≥∂u∂t−μ1Δu−urM−rLKMu.
Letu¯t,x,u0 be a solution of(17)∂u¯∂t−μ1Δu¯−u¯rM−rLKMu¯=0,then(18)∂u¯∂t−μ1Δu¯−u¯rM−rLKMu¯=0≥∂u∂t−μ1Δu−urM−rLKMu.
According to Theorem1, we can get ut,x,u0,v0≤u¯t,Mu, where Mu satisfies u0xC=maxx∈Ω¯u0x≤Mu. By the uniqueness theorem, it is obvious that the solution u¯t,Mu with initial conditions independent of x does not depend on x for t>0. Therefore, u¯t,Mu is the solution of the following ordinary differential equation:(19)du¯dt=u¯rM−rLKMu¯.u¯0,Mu=Mu.
Hence, we have(20)ut,x,u0,v0≤u¯t,Mu⟶rMKMrL,ast⟶∞.
Thus, there exists a positive constantM1 in system (6a)–(6c) such that ut,x≤M1, starting with some moment of time.
For predator populationv, by system (6b), we have(21)0=∂v∂t−μ2Δv−vet,xubt,xv+u−ct,x+δt,xv1+v=∂v∂t−μ2Δv−vet,xubt,xv+u−δt,x−ct,x−δt,x1+v≥∂v∂t−μ2Δv+cLv−eMM1bL,which implies that vt,x,u0,v0≤v¯t,Mv, where v¯t,Mv is a solution of the following initial value problem:(22)dv¯dt=−cLv¯+eMM1bL,v¯0,Mv=Mv,and Mv satisfies v0xC=maxx∈Ω¯v0x≤Mv. Obviously, we can obtain that(23)v¯t,Mv=Mve−cLt+eMM1bLcL⟶eMM1bLcL,ast⟶∞.
Therefore,vt,x,u0,v0 is also ultimately bounded.
This completes the proof.


### 3.2. Permanence

Theorem 5.
Under the condition (H), if the following inequalities(24a)rL−aMbL>0,(24b)eL−δM>0,hold, then system (6a)–(6c) is permanent, i.e., there exist positive constants mi and Mi (i = 1, 2) such that any solution of system (6a)–(6c) with nonnegative initial functions u0x≢0 and v0x≢0 satisfies ut,x,vt,x∈S=u,v:m1≤ut,x≤M1,m2≤vt,x≤M2, starting with a certain time.Proof.
Under the condition (H), we can know from Theorem 4 that there exists Mii=1,2 such that ut,x≤M1,vt,x≤M2, starting with some moment of time. By comparison principle, if u0x≥0≢0 and v0x≥0≢0, then ut,x,u0,v0>0 and vt,x,u0,v0>0 for all x∈Ω¯ and t>0.
Thus, for some smallε>0, we can get initial conditions uε,x,u0,v0,vε,x,u0,v0 separated from zero by the solution on the interval t≥ε. Without loss of generality, we assume that minx∈Ω¯u0x=mu,minx∈Ω¯v0x=mv. Then, the following inequality holds:(25)0=∂u∂t−μ1Δu−urt,x−rt,xKt,xu+at,xuvbt,xv+u≤∂u∂t−μ1Δu−urL−aMbL−rMKLu.
Obviously, we can get(26)0=∂u^∂t−μ1Δu^−u^rL−aMbL−rMKLu^≤∂u∂t−μ1Δu−urL−aMbL−rMKLu.
Consequently, fort≥0, we have(27)ut,x,u0,v0≥u^t,mu.
Thus, the solutionut,x,u0,v0 is bounded from below by a solution of the following logistic equation:(28)du^dt=u^rL−aMbL−rMKLu^,u^0=mu.
Thus, by Theorem1 and condition (24a) and (24b), we have(29)ut,x,u0,v0≥u^t,x⟶KLrL−aM/bLrM,ast⟶∞
Therefore, there exists a positive constantm1 such that ut,x,u0,v0≥m1 for t large enough.
By system (6b), the following inequality holds:(30)∂v∂t−μ2Δv−vet,xubt,xv+u−ct,x+δt,xv1+v≤∂v∂t−μ2Δv+δM−eLv+bMeLm1v2.
By a similar analysis tou, we have vt,x,u0,v0≥v^t,mv, where v^t,mv is a solution of the following system:(31)∂v^∂t−μ2Δv^+δM−eLv^+bMeLm1v^2,v^0=mv.
According to condition (24b), we can obtain that there exists a positive m2 such that vt,x,u0,v0≥m2 for t large enough. Thus, system (6a)–(6c) is permanent, starting with a certain time.
This completes the proof.


### 3.3. Extinction

In this section, we will discuss the extinction of predator species, and then the following theorem arrives in system (6a)–(6c).Theorem 6.
If the condition (H) holds, and(32)eM−cL<0,then, vt,x⟶0 as t⟶∞.Proof.
SupposeMv is a fixed positive constant guaranteeing Mv≤v0x, and v¯t,Mv is the solution of the following initial value problem:(33)∂v¯∂t=v¯eM−cL,v¯0,Mv=Mv.
By system (6b), we have(34)0=∂v∂t−μ2Δv+v−et,xubt,xv+u+ct,x+δt,xv1+v≥∂v∂t−μ2Δv+cL−eMv.
Thus, according to Theorem1, we can deduce that vt,x,u0,v0≤v¯t,Mv⟶0 as t⟶∞ if inequality (32) holds.
This completes the proof.


### 3.4. Periodic Solution

In this section, we will study the periodic solutions in system (6a)–(6c) by constructing a proper Lyapunov function.Theorem 7.
Under the condition (H), assume that system (6a)–(6c) is permanent, that is, there exist positive constants N and M such that an arbitrary solution of system (6a)–(6c) with nonnegative initial functions not identically equal to zero satisfies the condition:(35)ut,x,vt,x∈E=u,v:N≤ut,x≤M,N≤vt,x≤M,starting with a certain moment of time. If(36)λMW<0,where λM is the maximal eigenvalue of the following matrix:(37)E11E12E21E22,where(38)E11=2rM−rLKMN−aLbLN2bMM+M2,E22=2−δL+δM−cL1+N2+eMM2bLN+N2,E12=E21=aM+eMbL.Then system (6a)–(6c) has a unique and strictly positive τ-periodic solution, which is globally asymptotically stable.Proof.
Letut,x,vt,x and u¯t,x,v¯t,x be two solutions of system (6a)–(6c) bounded by constants N and M from below and above, respectively. Consider the following function:(39)Lt=∫Ωut,x−u¯t,x2+vt,x−v¯t,x2dx.
By system (6a)–(6c), we can get its derivative:(40)dLtdt=2∫Ωu−u¯∂u∂t−∂u¯∂tdx+2∫Ωv−v¯∂v∂t−∂v¯∂tdx=2μ1∫Ωu−u¯Δu−u¯dx+2μ2∫Ωv−v¯Δv−v¯dx+2∫Ωu−u¯ur−rKu−auvbv+u−u¯r−rKu¯−au¯v¯bv¯+u¯dx+2∫Ωv−v¯euvbv+u−c+δv1+vv−eu¯v¯bv¯+u¯−c+δv¯1+v¯v¯dx=I1+I2+I3+I4.
Then, from the boundary condition (6c),(41)I1+I2=−2μ1∫Ω∇2u−u¯dx−2μ2∫Ω∇2v−v¯dx≤−2μ1∫Ω∇u−u¯2dx−2μ2∫Ω∇v−v¯2dx≤0.
For other termsI3andI4,(42)I3+I4=2∫Ωu−u¯ru−u¯−rKu−u¯u+u¯+au¯v¯bv¯+u¯−auvbv+udx+2∫Ωv−v¯−δv−v¯+c−δ1+v¯v¯−c−δ1+vv+euvbv+u−eu¯v¯bv¯+u¯dx=2∫Ωu−u¯2r−rKu−u¯−abvv¯bv+ubv¯+u¯dx+2∫Ωv−v¯2euu¯bv+ubv¯+u¯−c−δ1+v1+v¯−δ+2∫Ωu−u¯v−v¯ebvv¯bv+ubv¯+u¯−auu¯bv+ubv¯+u¯dx≤2∫Ωu−u¯2rM−rLKMN−aLbLN2bMM+M2dx+2∫Ωv−v¯2eMM2bLN+N2+δM−cL1+N2−δL+2∫Ωu−u¯v−v¯aM+eMbLdx≤λM∫Ωu−u¯2+v−v¯2dx.
By condition (36), we have(43)Lt≤L0eλMt⟶0,ast⟶∞,which implies that ut,x−u¯t,x⟶0 and vt,x−v¯t,x⟶0 as t⟶∞, where ⋅ is the norm of the space L2Ω. Additionally, by condition (35), solutions of system (6a)–(6c) are bounded in the space C1+vΩ¯,ℝ2, where 0<v<2l−n/p and 1/2+n/2p<l<1.
Therefore,(44)limt⟶∞supx∈Ω¯ut,x−u¯t,x=0,limt⟶∞supx∈Ω¯vt,x−v¯t,x=0.
Consider the sequenceukτ,x,u0,v0,vkτ,x,u0,v0=Wkτ,W0,k∈N. Then, Wkτ,W0,k∈ℕ is compact in the space CΩ¯×CΩ¯. Let W¯ be a limit of this sequence, then Wτ,W¯=W¯.
Actually, becauseWτ,Wknτ,W0=Wknτ,Wτ,W0 and Wknτ,Wτ,W0−Wknτ,W0⟶0 as kn⟶∞, we have(45)Wτ,W¯−W¯C≤Wτ,W¯−Wτ,Wknτ,W0C+Wτ,Wknτ,W0−Wknτ,W0C+Wknτ,W0−W¯C⟶0asn⟶∞.
Thus, the sequenceWkτ,W0,k∈ℕ has a unique limit point. Otherwise, suppose that the sequence has two limit points W¯=limn⟶∞Wknτ,W0 and W˜=limn⟶∞Wknτ,W0, then we can get the following result from (45) and W˜=Wknτ,W˜:(46)W¯−W˜C≤W¯−Wknτ,W0C+Wknτ,W0−W˜C⟶0,n⟶∞.
Hence,W¯=W˜. The solution ut,x,u¯,v¯,vt,x,u¯,v¯ is the unique periodic solution of system (6a)–(6c), and it is asymptotically stable using equation (44).
This completes the proof.


## 3.1. Boundedness

From the biological and ecological viewpoint, we are always interested in the nonnegative solutions. Thus, the following theorem is given first in system (6a)–(6c).Theorem 3.
Suppose that the condition (H) holds, then nonnegative and positive quadrants of ℝ2 are positively invariant for system (6a)–(6c).Proof.
Letut,x,u0,v0,vt,x,u0,v0 be a solution of system (6a)–(6c) with initial condition u0x≥0≢0,v0x≥0≢0. Additionally, u^ is a solution of the following system:(13)∂u^∂t−μ1Δu^−u^rL−aMbL−rMKLu^=0,u^0,x=u0x.
From system (6a), we can obtain(14)∂u∂t−μ1Δu−rt,xu1−uKt,x+at,xuvbt,xv+u≤∂u∂t−μ1Δu−urL−aMbL−rMKLu,which implies u^t,x is a lower solution of system (6a). According to Theorem 2, it is obvious that u^t,x≥0 for all x∈Ω¯ and t>0. Furthermore, due to u0x≥0≢0, u^t,x>0 holds for all x∈Ω¯ and t>0. Thus, ut,x>0 holds because ut,x is bounded from below by positive function u^t,x.
For system (6b), it can be simply verified that v^t,x is a lower solution of system (6b), where v^t,x satisfies(15)∂v^∂t−μ2Δv^+δMv^=0,v^0,x=v0x.
By the similar argument tout,x, we can prove the positiveness of vt,x.
This completes the proof.
Based on Theorem3, we will discuss ultimate boundedness of solutions in system (6a)–(6c), and then the following theorem can be obtained.Theorem 4.
If the condition (H) holds, then all solutions of system (6a)–(6c) with nonnegative initial conditions are ultimately bounded.Proof.
From system (6a), it can be found that the following inequality holds:(16)0=∂u∂t−μ1Δu−rt,xu1−uKt,x+at,xuvbt,xv+u≥∂u∂t−μ1Δu−urM−rLKMu.
Letu¯t,x,u0 be a solution of(17)∂u¯∂t−μ1Δu¯−u¯rM−rLKMu¯=0,then(18)∂u¯∂t−μ1Δu¯−u¯rM−rLKMu¯=0≥∂u∂t−μ1Δu−urM−rLKMu.
According to Theorem1, we can get ut,x,u0,v0≤u¯t,Mu, where Mu satisfies u0xC=maxx∈Ω¯u0x≤Mu. By the uniqueness theorem, it is obvious that the solution u¯t,Mu with initial conditions independent of x does not depend on x for t>0. Therefore, u¯t,Mu is the solution of the following ordinary differential equation:(19)du¯dt=u¯rM−rLKMu¯.u¯0,Mu=Mu.
Hence, we have(20)ut,x,u0,v0≤u¯t,Mu⟶rMKMrL,ast⟶∞.
Thus, there exists a positive constantM1 in system (6a)–(6c) such that ut,x≤M1, starting with some moment of time.
For predator populationv, by system (6b), we have(21)0=∂v∂t−μ2Δv−vet,xubt,xv+u−ct,x+δt,xv1+v=∂v∂t−μ2Δv−vet,xubt,xv+u−δt,x−ct,x−δt,x1+v≥∂v∂t−μ2Δv+cLv−eMM1bL,which implies that vt,x,u0,v0≤v¯t,Mv, where v¯t,Mv is a solution of the following initial value problem:(22)dv¯dt=−cLv¯+eMM1bL,v¯0,Mv=Mv,and Mv satisfies v0xC=maxx∈Ω¯v0x≤Mv. Obviously, we can obtain that(23)v¯t,Mv=Mve−cLt+eMM1bLcL⟶eMM1bLcL,ast⟶∞.
Therefore,vt,x,u0,v0 is also ultimately bounded.
This completes the proof.


## 3.2. Permanence

Theorem 5.
Under the condition (H), if the following inequalities(24a)rL−aMbL>0,(24b)eL−δM>0,hold, then system (6a)–(6c) is permanent, i.e., there exist positive constants mi and Mi (i = 1, 2) such that any solution of system (6a)–(6c) with nonnegative initial functions u0x≢0 and v0x≢0 satisfies ut,x,vt,x∈S=u,v:m1≤ut,x≤M1,m2≤vt,x≤M2, starting with a certain time.Proof.
Under the condition (H), we can know from Theorem 4 that there exists Mii=1,2 such that ut,x≤M1,vt,x≤M2, starting with some moment of time. By comparison principle, if u0x≥0≢0 and v0x≥0≢0, then ut,x,u0,v0>0 and vt,x,u0,v0>0 for all x∈Ω¯ and t>0.
Thus, for some smallε>0, we can get initial conditions uε,x,u0,v0,vε,x,u0,v0 separated from zero by the solution on the interval t≥ε. Without loss of generality, we assume that minx∈Ω¯u0x=mu,minx∈Ω¯v0x=mv. Then, the following inequality holds:(25)0=∂u∂t−μ1Δu−urt,x−rt,xKt,xu+at,xuvbt,xv+u≤∂u∂t−μ1Δu−urL−aMbL−rMKLu.
Obviously, we can get(26)0=∂u^∂t−μ1Δu^−u^rL−aMbL−rMKLu^≤∂u∂t−μ1Δu−urL−aMbL−rMKLu.
Consequently, fort≥0, we have(27)ut,x,u0,v0≥u^t,mu.
Thus, the solutionut,x,u0,v0 is bounded from below by a solution of the following logistic equation:(28)du^dt=u^rL−aMbL−rMKLu^,u^0=mu.
Thus, by Theorem1 and condition (24a) and (24b), we have(29)ut,x,u0,v0≥u^t,x⟶KLrL−aM/bLrM,ast⟶∞
Therefore, there exists a positive constantm1 such that ut,x,u0,v0≥m1 for t large enough.
By system (6b), the following inequality holds:(30)∂v∂t−μ2Δv−vet,xubt,xv+u−ct,x+δt,xv1+v≤∂v∂t−μ2Δv+δM−eLv+bMeLm1v2.
By a similar analysis tou, we have vt,x,u0,v0≥v^t,mv, where v^t,mv is a solution of the following system:(31)∂v^∂t−μ2Δv^+δM−eLv^+bMeLm1v^2,v^0=mv.
According to condition (24b), we can obtain that there exists a positive m2 such that vt,x,u0,v0≥m2 for t large enough. Thus, system (6a)–(6c) is permanent, starting with a certain time.
This completes the proof.


## 3.3. Extinction

In this section, we will discuss the extinction of predator species, and then the following theorem arrives in system (6a)–(6c).Theorem 6.
If the condition (H) holds, and(32)eM−cL<0,then, vt,x⟶0 as t⟶∞.Proof.
SupposeMv is a fixed positive constant guaranteeing Mv≤v0x, and v¯t,Mv is the solution of the following initial value problem:(33)∂v¯∂t=v¯eM−cL,v¯0,Mv=Mv.
By system (6b), we have(34)0=∂v∂t−μ2Δv+v−et,xubt,xv+u+ct,x+δt,xv1+v≥∂v∂t−μ2Δv+cL−eMv.
Thus, according to Theorem1, we can deduce that vt,x,u0,v0≤v¯t,Mv⟶0 as t⟶∞ if inequality (32) holds.
This completes the proof.


## 3.4. Periodic Solution

In this section, we will study the periodic solutions in system (6a)–(6c) by constructing a proper Lyapunov function.Theorem 7.
Under the condition (H), assume that system (6a)–(6c) is permanent, that is, there exist positive constants N and M such that an arbitrary solution of system (6a)–(6c) with nonnegative initial functions not identically equal to zero satisfies the condition:(35)ut,x,vt,x∈E=u,v:N≤ut,x≤M,N≤vt,x≤M,starting with a certain moment of time. If(36)λMW<0,where λM is the maximal eigenvalue of the following matrix:(37)E11E12E21E22,where(38)E11=2rM−rLKMN−aLbLN2bMM+M2,E22=2−δL+δM−cL1+N2+eMM2bLN+N2,E12=E21=aM+eMbL.Then system (6a)–(6c) has a unique and strictly positive τ-periodic solution, which is globally asymptotically stable.Proof.
Letut,x,vt,x and u¯t,x,v¯t,x be two solutions of system (6a)–(6c) bounded by constants N and M from below and above, respectively. Consider the following function:(39)Lt=∫Ωut,x−u¯t,x2+vt,x−v¯t,x2dx.
By system (6a)–(6c), we can get its derivative:(40)dLtdt=2∫Ωu−u¯∂u∂t−∂u¯∂tdx+2∫Ωv−v¯∂v∂t−∂v¯∂tdx=2μ1∫Ωu−u¯Δu−u¯dx+2μ2∫Ωv−v¯Δv−v¯dx+2∫Ωu−u¯ur−rKu−auvbv+u−u¯r−rKu¯−au¯v¯bv¯+u¯dx+2∫Ωv−v¯euvbv+u−c+δv1+vv−eu¯v¯bv¯+u¯−c+δv¯1+v¯v¯dx=I1+I2+I3+I4.
Then, from the boundary condition (6c),(41)I1+I2=−2μ1∫Ω∇2u−u¯dx−2μ2∫Ω∇2v−v¯dx≤−2μ1∫Ω∇u−u¯2dx−2μ2∫Ω∇v−v¯2dx≤0.
For other termsI3andI4,(42)I3+I4=2∫Ωu−u¯ru−u¯−rKu−u¯u+u¯+au¯v¯bv¯+u¯−auvbv+udx+2∫Ωv−v¯−δv−v¯+c−δ1+v¯v¯−c−δ1+vv+euvbv+u−eu¯v¯bv¯+u¯dx=2∫Ωu−u¯2r−rKu−u¯−abvv¯bv+ubv¯+u¯dx+2∫Ωv−v¯2euu¯bv+ubv¯+u¯−c−δ1+v1+v¯−δ+2∫Ωu−u¯v−v¯ebvv¯bv+ubv¯+u¯−auu¯bv+ubv¯+u¯dx≤2∫Ωu−u¯2rM−rLKMN−aLbLN2bMM+M2dx+2∫Ωv−v¯2eMM2bLN+N2+δM−cL1+N2−δL+2∫Ωu−u¯v−v¯aM+eMbLdx≤λM∫Ωu−u¯2+v−v¯2dx.
By condition (36), we have(43)Lt≤L0eλMt⟶0,ast⟶∞,which implies that ut,x−u¯t,x⟶0 and vt,x−v¯t,x⟶0 as t⟶∞, where ⋅ is the norm of the space L2Ω. Additionally, by condition (35), solutions of system (6a)–(6c) are bounded in the space C1+vΩ¯,ℝ2, where 0<v<2l−n/p and 1/2+n/2p<l<1.
Therefore,(44)limt⟶∞supx∈Ω¯ut,x−u¯t,x=0,limt⟶∞supx∈Ω¯vt,x−v¯t,x=0.
Consider the sequenceukτ,x,u0,v0,vkτ,x,u0,v0=Wkτ,W0,k∈N. Then, Wkτ,W0,k∈ℕ is compact in the space CΩ¯×CΩ¯. Let W¯ be a limit of this sequence, then Wτ,W¯=W¯.
Actually, becauseWτ,Wknτ,W0=Wknτ,Wτ,W0 and Wknτ,Wτ,W0−Wknτ,W0⟶0 as kn⟶∞, we have(45)Wτ,W¯−W¯C≤Wτ,W¯−Wτ,Wknτ,W0C+Wτ,Wknτ,W0−Wknτ,W0C+Wknτ,W0−W¯C⟶0asn⟶∞.
Thus, the sequenceWkτ,W0,k∈ℕ has a unique limit point. Otherwise, suppose that the sequence has two limit points W¯=limn⟶∞Wknτ,W0 and W˜=limn⟶∞Wknτ,W0, then we can get the following result from (45) and W˜=Wknτ,W˜:(46)W¯−W˜C≤W¯−Wknτ,W0C+Wknτ,W0−W˜C⟶0,n⟶∞.
Hence,W¯=W˜. The solution ut,x,u¯,v¯,vt,x,u¯,v¯ is the unique periodic solution of system (6a)–(6c), and it is asymptotically stable using equation (44).
This completes the proof.


## 4. Numerical Results

In the previous section, we have obtained some interesting results of system (6a)–(6c). However, due to the complexity of system (6a)–(6c), it becomes much more difficult to provide in-depth analysis. Thus, here, we perform some numerical simulations to investigate prey-predator dynamics further.According to Theorem5, when rL−aM/bL>0 and eL−δM>0 holds, system (6a)–(6c) is permanent under condition (H). Figure 1 shows that system (6a)–(6c) is permanent, where rL−aM/bL≈0.1>0 and eL−δM=0.528>0. When e=0.005−0.002sinπ∗t/10, other parameters are the same as the ones in Figure 1, and we can get a numerical solution of system (6a)–(6c) (see Figure 2). It is obvious that predator population v is extinct ultimately, which is consistent with Theorem 6 because eM−cL=−0.008<0.Figure 1
Numerical solutions of system (6a)–(6c): (a) prey population u and (b) predator population v, where rt,x=1.1+0.1sinπ∗t/10, Kt,x=2+0.5∗cosπ∗t/10, at,x=0.8+0.005∗cosπ∗t/10, bt,x=0.9+0.005∗cosπ∗t/10, et,x=0.7−0.002∗cosπ∗t/10, ct,x=0.02+0.005∗cosπ∗t/10, δt,x=0.12+0.05∗cosπ∗t/10, μ1=1, and μ2=1.
(a)
(b)Figure 2
Numerical solutions of system (6a)–(6c): (a) prey population u and (b) predator population v.
(a)
(b)In section3.4, the existence of periodic solution was discussed, and its stability and uniqueness were analyzed as well. In fact, Figure 1 has shown the existence of a periodic solution. Yet, we here take another set of function corresponding to the parameters of system (6a)–(6c), which is only the periodic function of time t with period 200. The corresponding numerical solutions are shown in Figure 3. Clearly, the numerical solution is periodic in t with the period of 200 (see Figures 3(c) and 3(d)), but it is homogeneous in space (see Figures 3(a) and 3(b)). Compared to Figure 3, we consider another situation that the parameters of system (6a)–(6c) are functions with respect to both time t and space x. We find the solution is still periodic, but it is heterogeneous in space (see Figure 4). It is evident that the spatial heterogeneity is the reason giving rise to the oscillation of the solution in space.Figure 3
Numerical solutions of system (6a)–(6c): (a) prey population u, (b) predator population v, (c) profile of ut at x=50, and (d) profile of vt at x=50. Here, rt,x=1.2+0.1sinπ∗t/100, Kt,x=2+0.5cosπ∗t/100, at,x=0.8+0.005cosπ∗t/100, bt,x=0.9+0.005cosπ∗t/100, et,x=0.7−0.005sinπ∗t/100, ct,x=0.02+0.01cosπ∗t/100, δt,x=0.12+0.05cosπ∗t/100, μ1=1, and μ2=1.
(a)
(b)
(c)
(d)Figure 4
Numerical solutions of system (6a)–(6c): (a) prey population u, (b) predator population v, (c) profile of ut at x=50, and (d) profile of ux at t=1200. Here, rt,x=1.2+0.1sinπ∗t/100+0.01cos4∗π∗x/100∗sin4∗π∗x/100, Kt,x=2+0.5cosπ∗t/100+cos4∗π∗x/100, at,x=0.8+0.005cosπ∗t/100+0.1cos4∗π∗x/100, bt,x=0.9+0.005cosπ∗t/100+0.1cos4∗π∗x/100, et,x=0.7−0.005sinπ∗t/100+0.1cos4∗π∗x/100, ct,x=0.02+0.01cosπ∗t/100+0.001cos4∗π∗x/100, δt,x=0.12+0.05cosπ∗t/100+0.01cos4∗π∗x/100, μ1=1, and μ2=1.
(a)
(b)
(c)
(d)Additionally, we find that the spatial heterogeneity can promote the permanence of the system. Figure5 indicates that the extinction occurs in system (6a)–(6c). However, when we set rt,x=0.9+0.001sin5∗π∗t/100+0.1cos5∗π∗x/100 (other parameters are the same as the ones in Figure 5), we get a very interesting result, that is, system (6a)–(6c) becomes permanent (see Figure 6). Likewise, other parameters are explored by repeating the same procedure, and similar results are obtained, which are omitted here. Obviously, the spatial heterogeneity plays an important role in dynamics of system (6a)–(6c).Figure 5
Numerical solutions of system (6a)–(6c): (a) prey population u and (b) predator population v, where rt,x=0.9+0.001sin5∗π∗t/100, Kt,x=10+0.5cos5∗π∗t/100, at,x=0.8+0.005cos5∗π∗t/100, bt,x=0.78+0.005cos5∗π∗t/100, et,x=0.7−0.005sin5∗π∗t/100, ct,x=0.02+0.001cos5∗π∗t/100, δt,x=0.12+0.05cos5∗π∗t/100, μ1=1, and μ2=1.
(a)
(b)Figure 6
Numerical solutions of system (6a)–(6c): (a) prey population u and (b) predator population v.
(a)
(b)Let all the parameters be constant, then there exists a nonconstant stationary solution in system (6a)–(6c), as shown in Figure 7. Furthermore, we consider the parameters depending on time t based on Figure 7, but the result shows that the extinction of both prey and predator occurs (see Figure 8). However, when the parameter r depends on space x except for time t, it can be found from Figure 9 that system (6a)–(6c) is still permanent. Although both Figures 7 and 9 show the spatial heterogeneity of population distribution, their natures are different. The spatial heterogeneity in Figure 7 is induced by the diffusion, while in Figure 9 it depends on the space variation of parameters. For other parameters of system (6a)–(6c), we can obtain similar results, which are omitted here.Figure 7
Numerical solutions of system (6a)–(6c): (a) prey population u and (b) predator population v, where rt,x=0.9, Kt,x=10, at,x=0.8, bt,x=0.8, et,x=0.7, ct,x=0.02, δt,x=0.12, μ1=0.01, and μ2=20.
(a)
(b)Figure 8
Numerical solutions of system (6a)–(6c): (a) prey population u and (b) predator population v, where rt,x=0.9+0.001sin5∗π∗t/100, Kt,x=10+0.5sin5∗π∗t/100, at,x=0.8+0.005sin5∗π∗t/100, bt,x=0.8+0.005sin5∗π∗t/100, et,x=0.7−0.005sin5∗π∗t/100, ct,x=0.02+0.001sin5∗π∗t/100, δt,x=0.12+0.05sin5∗π∗t/100, μ1=0.01, and μ2=20.
(a)
(b)Figure 9
Numerical solutions of system (6a)–(6c): (a) prey population u and (b) predator population v, where rt,x=0.9+0.001sin5∗π∗t/100+0.1cost∗π∗x/100 and other parameters are the same as the ones in Figure 8.
(a)
(b)


## 5. Conclusion

In this paper, we first propose a reaction-diffusion system (6a)–(6c) to describe the interaction between the prey and the predator, where the spatial heterogeneity and the time-periodic environment are considered. In order to study the boundedness of solution, the positive invariance of system (6a)–(6c) is discussed, and the results demonstrate that nonnegative and positive quadrants of ℝ2 are always positively invariant for system (6a)–(6c) when the condition (H) holds. Based on this, we find that all solutions of system (6a)–(6c) are ultimately bounded as long as the initial conditions are nonnegative. Also, we discuss the permanence of system (6a)–(6c) and obtain the sufficient conditions. Moreover, we derive the sufficient conditions for the extinction of predator population. Obviously, these conditions are very significant for studies of permanence and extinction of the system. When system (6a)–(6c) is permanent, we discuss the existence of a periodic solution, which suggests that a unique and strictly positive periodic solution with fixed period exists under certain conditions.According to theoretical analysis, some numerical results are given, which show further dynamics in system (6a)–(6c). Results from literature [45, 46] indicate that Turing patterns can exist in system (5a)–(5c) (i.e., system (6a)–(6c) without spatial heterogeneity and time-periodic environment). After taking time-periodic environment into account, we find that both prey population and predator population are extinct. However, when the combination of spatial heterogeneity and time-periodic environment is considered, it is demonstrated that both prey population and predator population are permanent, which means that spatial heterogeneity tends to enhance the persistence of prey and predator population. Additionally, when prey population and predator population are permanent, our results show that solutions of system (6a)–(6c) seem to be periodic because of the time-periodic environment and the spatial heterogeneity. Thus, we want to emphasize that spatial heterogeneity and time-periodic environment indeed play a significant role in prey-predator dynamics.



---

*Source: 7134869-2020-04-14.xml*