# Spatiotemporal Complexity of a Leslie-Gower Predator-Prey Model with the Weak Allee Effect

Authors: Yongli Cai; Caidi Zhao; Weiming Wang
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

We investigate a diffusive Leslie-Gower predator-prey model with the additive Allee effect on prey subject to the zero-flux boundary conditions. Some results of solutions to this model and its corresponding steady-state problem are shown. More precisely, we give the
stability of the positive constant steady-state solution, the refineda priori estimates of positive solution, and the nonexistence and existence of the positive nonconstant solutions. We carry out the analytical study for two-dimensional system in detail and find out the certain conditions for Turing instability. Furthermore, we perform numerical simulations and show that the model exhibits a transition from stripe-spot mixtures growth to isolated spots and
also to stripes. These results show that the impact of the Allee effect essentially increases the model spatiotemporal complexity.


---

## Body

## 1. Introduction

The dynamics of a predator-prey model in a homogeneous environment can be described by the following reaction-diffusion equations:(1)∂u(x,t)∂t=uF(u)-G(u)v+D1Δu,∂v(x,t)∂t=vQ(u,v)+D2Δv,
where u(x,t) and v(x,t) are the densities of prey and predator at time t and position x∈Ω⊂ℝm, respectively. The Laplace operator Δ describes the spatial dispersal with passive diffusion; D1>0 and D2>0 are the diffusion coefficients corresponding to species u and v. F(u) describes the per-capita growth rate of the prey; G(u) is the functional response of the predator, which corresponds to the saturation of their appetites and reproductive capacity; Q(u,v), the so-called numerical response, is the per-capita growth rate of the predator [1–4].FunctionsF(u), G(u), and Q(u,v) can be formulated in various specific situations. In general, F(u) is of the standard logistic growth:
(2)F(u)=r1(1-uK),
which was first created by Verhulst [5]. Here K>0 is the prey carrying capacity and r1>0 is the intrinsic growth rate of prey.Some conventional functional response functionsG(u) include Holling types I, II, and III (see [6–10]). Among many possible choices of G(u), the Holling type-II functional response is most commonly used in the ecological literature, which is defined by [11]:
(3)G(u)=c1uu+K1,
where c1>0 describes the maximum predation rate and K1≥0 measures the extent to which environment provides protection to prey u. The Leslie-Gower type numerical response Q(u,v) is given by
(4)Q(u,v)=r2-c2vδu,
which was first proposed by Leslie [12], and has been discussed by Leslie and Gower [13] and Pielou [14]. A modified version of Leslie-Gower functional response
(5)Q(u,v)=r2-c2vu+K2
is given by Aziz-Alaoui et al. [15, 16]. Here, r2>0 describes the growth rate of the predator v; c2>0 has a similar meaning to c1; δ>0 takes on the role of the prey-dependent carrying capacity for the predator; K2≥0 is the extent to which environment provides protection to predator v. Hence, we can rewrite model (1) as follows:
(6)∂u(x,t)∂t=r1u(1-uK)-c1uvu+K1+D1Δu,∂v(x,t)∂t=v(r2-c2vu+K2)+D2Δv.
The biological significance of all parameters in model (6) is as above.For model (6), in the case of K1=K2=0, Du et al. [17, 18] mainly focused attention on the steady-state problem and observed some quite interesting phenomena of pattern formation. In the case of K1>0,K2=0, the so-called Holling-Tanner model, Peng and Wang [19, 20] analyzed the global stability of the unique positive constant steady-state and established the results for the existence and nonexistence of positive nonconstant steadystates; Shi and coworkers [21] studied the global attractor and persistence property, local and global asymptotic stability of the unique positive constant equilibrium, and the existence and nonexistence of nonconstant positive steady-states; Li et al. [22] considered the Turing and Hopf bifurcations of the equilibrium solutions; Liu and Xue [23] found the model exhibits the spotted, black-eye, and labyrinthine patterns. For model (6), that is, K1≠K2≠0, Camara and Aziz-Alaoui [24–26] paid more attention to pattern formation in the spatial domain and observed the labyrinth, chaos, and spiral wave patterns.On the other hand, in population dynamics, any mechanism that can lead to a positive relationship between a component of individual fitness and either the number or density of conspecific can be termed an Allee effect [27–30], starting with the pioneer work of ecologist Allee [31]. In particular, theoretical investigations have shown that an Allee effect can greatly increase the likelihood of local and global extinction [32] and can lead to a rich variety of dynamical effects. As a consequence, it is necessary to explore the influence of Allee effect in the growth of a population.The Allee effect has been modeled in different ways [33–37]. From an ecological viewpoint, the Allee effect has been modeled into strong and weak ones [33, 38–42]. In a recent analytic approach by Wang and Kot [38], the Allee effect is “strong” if the sign of the growth function uF(u) in the limit of law density is negative; that is,
(7)d(uF(u))dt|u=0<0.
It is “weak” if the sign of the growth function uF(u) in the limit of law density is positive; that is,
(8)d(uF(u))dt|u=0>0.
The strong Allee effect introduces a population threshold, and the population must surpass this threshold to grow. In contrast, the weak case has not any threshold [10, 35, 38, 42].In particular, the growth function considering Allee effect is expressed by the equation:(9)uF(u)=r1u(1-uK-m1u+b1),
having an additive Allee effect, which was first deduced in [43] and applied in [34–36]. Where m1/(u+b1) is the term of additive Allee effect, m1 and b1 are the Allee effect constants. It should be noted that, if m1<b1, the Allee effect in (9) is weak, while if m1>b1, the Allee effect in (9) is strong.Based on the above discussions, in this paper, we rigorously consider the spatiotemporal dynamics of the following modified Leslie-Gower predation model with the additive Allee effect on prey:(10)∂u(x,t)∂t=r1u(1-uK-m1u+b1)-c1uvu+K1+D1Δu,∂v(x,t)∂t=v(r2-c2vu+K2)+D2Δv.We make a change of variables:(11)(u,v,t)=(Ku~,Kv~,t~r1).
And for the sake of convenience, we still use variables u, v instead of u~, v~. Thus, the model to be studied is as follows:
(12)∂u∂t=d1Δu+u(1-u-mu+b-cvu+k1),x∈Ω,t>0,∂v∂t=d2Δv+sv(1-vu+k2),x∈Ω,t>0,∂u∂ν=∂v∂ν=0,x∈∂Ω,t>0,u(x,0)=u0(x)≥0,v(x,0)=v0(x)≥0,x∈Ω,
where
(13)m=m1K,b=b1K,c=c1r1,k1=K1K,s=r2r1,c2=r2,k2=K2K,d1=D1r1,d2=D2r1.
Here, s is the growth rate of the predator v. mu/(u+b) is the term of additive Allee effect, and m and b are the Allee effect constants. Ω⊂ℝ2 is a bounded domain with smooth boundary ∂Ω, and ν is the outward unit normal vector on ∂Ω. The initial data u0(x) and v0(x) are continuous functions on Ω¯, and the zero-flux boundary conditions mean that model (12) is self-contained and has no population flux across the boundary ∂Ω [44, 45].By the standard theory for semilinear parabolic systems (see, [46]), we have model (12) that admits a unique classical solution (u(x,t),v(x,t)) for all time.The stationary problem of model (12), which may display the dynamical behavior of solutions to model (12) as time goes to infinity, satisfies the following elliptic system:
(14)-d1Δu=u(1-u-mu+b-cvu+k1),x∈Ω,-d2Δv=sv(1-vu+k2),x∈Ω,∂u∂ν=∂v∂ν=0,x∈∂Ω.
Unless otherwise specified, in this paper, we always assume that m<b; that is, we only focus on the case of weak Allee effect.The rest of the paper is organized as follows. In Section2, we investigate the stability of nonnegative constant steady-state solutions. In Section 3, we mainly give a priori upper and lower bounds for positive solutions of model (14). In Section 4 we discuss existence and nonexistence of nonconstant positive solutions, which might give us some suggestions on the conditions under which the patterns may or may not occur. In Section 5 we first use the method of linearized stability analysis to deduce the conditions under which the Turing instability might occur, and next we perform a series of numerical simulations to show the occurrence of different patterns. Finally, in the last section we make a summary to our results and give some concluding remarks.


## 2. Dynamics Analysis of Model (12)

### 2.1. The Existence of the Constant Steady-State Solution

It is easy to verify that model (12) has the following nonnegative constant steady-state solutions:(i)
the trivial constant solutionE0=(0,0) (extinction of two species);(ii)
the semitrivial constant solutionE1=(0,k2) (extinction of the prey);(iii)
the semitrivial constant solutionE2=((1-b+(1-b)2-4(m-b))/2,0) (extinction of the predator);(iv)
the unique positive constant solutionE3=(u*,v*) (coexistence of two species), where v*=u*+k2, and u* is a real positive root of the cubic
(15)ψ(w)=w3+3η1w2+3η2w+η3=0,
where 3η1=b+c+k1-1,3η2=bc+bk1+ck2+m-b-k1,η3=bck2+k1(m-b).By the transformationz=w+η1, (15) is reduced to
(16)h(z)=z3+3pz+q=0,
where
(17)p=η2-η12,(18)q=η3-3η1η2+2η13.It is worthy to note that ifq=η3-3η1η2+2η13<0 holds, (16) has a real positive root. Considering (17) and (18), from q=0, one can determine m=m*, where
(19)m*≜2b-c-k1+19(b+c-2k1-1)×(b2+b-bc-bk1-2c2+4c-4ck1-2k12-5k1+9ck2-2).Hence, we have the following lemma regarding the existence of the positive constant steady-state solution of model (12).Lemma 1.
If either of the following inequalities holds:(20)2k1+1-c-b>0,0<m<min{b,m*};2k1+1-c-b<0,0<m*<m<b,
model (12) has a unique positive constant steady-state solution E3=(u*,v*).


### 2.2. Stability of the Constant Steady-State Solution

In this subsection, we will analyze the asymptotical stability of the nonnegative constant solutions for the corresponding reaction-diffusion dynamics (12).For sake of simplicity, we rewrite model (12) in the vectorial form:
(21)wt=DΔw+H(w)x∈Ω,t>0∂w∂ν=0,x∈∂Ω,t>0,w(x,0)=(u0(x),v0(x)T),x∈Ω,
where w=(u,v)T, D=diag(d1,d2), and
(22)H(w)=(u(1-u-mu+b-cvu+k1)sv(1-vu+k2)).Let0=μ0<μ1<μ2<⋯ be the eigenvalues of the operator-Δ on Ω with the zero-flux boundary conditions. And set
(23)X={w∈[H2(Ω)]2∣∂νw=0on∂Ω},E(μ)={ϕ∣-Δϕ=μϕinΩ,∂νϕ=0on∂Ω},withμ∈ℝ1,{ϕij∣j=1,…,dimE(μi)} being an orthonormal basis of E(μi), and Xij={cϕij∣c∈ℝ2}; then
(24)X=⨁i=1∞Xi,
where Xi=⨁j=1dimE(μi)Xij.LetE=(u,v) be any arbitrary constant steady-state solution of model (12). And the linearization of model (12) at the constant steady-state solution can be expressed by
(25)wt=£(w)=DΔw+J(u,v)w,
where J(u,v)=(ξ11ξ12ξ21ξ22) and
(26)ξ11=1-u-mu+b+u(-1+m(u+b)2)-cvk1+u+cuv(k1+u)2,ξ12=-cuk1+u,ξ21=sv2(u+k2)2,ξ22=s(1-vu+k2)-svu+k2.From [46], it is known that if all the eigenvalues of the operator £ have negative real parts, then E=(u,v) is asymptotically stable; if there is an eigenvalue with positive real part, then E=(u,v) is unstable; if all the eigenvalues have nonpositive real parts while some eigenvalues have zero real parts, then the stability of E=(u,v) cannot be determined by the linearization [10].For eachi≥0, Xi is invariant under the operator £, and λ is an eigenvalue of £ if and only if λ is an eigenvalue of the matrix Mi=-μiD+J(u,v) for some i≥0.In the following, we denoteξjj[k]=ξjj evaluated at Ek, j=1,2 and k=0,1,2,3. So, the local stability of the constant steady-state solution can be analyzed as follows.Theorem 2.
For any positive parameters,(a)
the trivial constant solutionE0=(0,0) is unstable;(b)
the semitrivial constant solutionE1=(0,k2) is
(b1)
locally asymptotically stable ifc>(k1/k2) or (k1/k2)(1-m/b)<c≤(k1/k2) holds(b2)
unstable ifc<(k1/k2)(1-m/b) holds;(c)
the semitrivial constant solutionE2=((1-b+(1-b)2-4(m-b))/2,0) is unstable.Proof.
The stability of the constant steady-state solution is reduced to consider the characteristic equation:(27)det(λI-Mi)=λ2-tr(Mi)λ+det(Mi),
with
(28)tr(Mi)=-μi(d1+d2)+ξ11+ξ22,det(Mi)=d1d2μi2-(ξ11d2+ξ22d1)+detJ(u,v).(a)
JE0=(1-(b/m)00s), for i=0; the eigenvalues are 1-(b/m)>0 and s>0, so E0=(0,0) is unstable.(b)
We can obtainJE1=(1-(b/m)-(ck2/k1)0s-s).
(b1)
Ifk1<ck2 or k1≥ck2andb(1-(ck2/k1))<m<b hold, then 1-(b/m)-(ck2/k1)<0, so for i≥0,
(29)tr(Mi)=-μi(d1+d2)+1-bm-ck2k1-s<0,det(Mi)=d1d2μi2-(1-bm-ck2k1d2-sd1)-s(1-bm-ck2k1)>0.Hence,E1 is locally asymptotically stable.(b2)
Whenk1≥ck2and0<m<b(1-(ck2/k1)), then 1-(b/m)-(ck2/k1)>0. For i≥0, det(Mi)=-s(1-(b/m)-(ck2/k1))<0, which implies that (27) has at least one root with positive real part. Hence, E1 is an unstable steady-state solution of model (12).(c)
Sinceξ11[2]=-((2(1+b)2-4m((1+b)2-4m+1+b-2m))/(b+1+(1+b)2-4m)2),ξ22[2]=s, for i=0, one of the eigenvalues is s>0, so E2 is unstable.
The proof is complete.Straightforward calculations show that(30)JE3=(ξ11[3]-cu*u*+k1s-s),
where ξ11[3]=((cu*(u*+k2))/(u*+k1)2)+(mu*/(u*+b)2)-u*.The determinant ofJE3 is given by(31)det(JE*)=su*((u*+k1)2((u*+b)2-m)+c(u*+b)2(k1-k2))(u*+k1)2(u*+b)2;then, the sign of det(JE3) depends on the factor ϕ(u*):
(32)ϕ(u*)=(u*+k1)2((u*+b)2-m)+c(u*+b)2(k1-k2)=(u*+b+k1+1-c)ψ(u*)+η1u*2+η2u*+η3=ρ1u*2+ρ2u*+ρ3,
where ψ(·) is the same definition as (15), and
(33)ρ1=1+b+k1+c2+ck1+bk1-2m-cb-2ck2-2c,ρ2=-ck2+b-m+2bk1-2cb+b2-ck1+c2k2+c2b,+mc-bm-cb2-4mk1+k12-4cbk2-ck2k1+b2k1+k1+bk12+2bk1c,ρ3=b2k12-bk1c-mbk1+mck1+b2k1c+c2k2b-bk1ck2-cbk2+bk1-mk1+b2k1+bk12-2mk12-2ck2b2.Therefore, we have the following.Theorem 3.
Assume thatξ11[3]<s, ρ1u*2+ρ2u*+ρ3>0, and the first eigenvalue μ1 subject to the zero-flux boundary conditions satisfies
(34)μ1>max{0,ξ11[3]d1-sd2}.
Then the positive constant steady-state solution E3=(u*,u*+k2) of model (15) is uniformly asymptotically stable.Proof.
Whenξ11[3]<s, ρ1u*2+ρ2u*+ρ3>0, then tr(JE3)<0 and det(JE3)>0. So, for i≥0,
(35)tr(Mi)=-(d1+d2)μi+ξ11[3]-sdet(Mi)=μi(d1d2μi-ξ11[3]d2+sd1)+det(JE3).
Note that for any i≥0, we have det(Ai)>0>tr(Ai). Therefore, the eigenvalues of the matrix -μiD+JE3 have negative real parts. It thus follows from the Routh-Hurwitz criterion that, for each i≥0, the two roots λi1 and λi2 of φi(λ)=0 all have negative real parts.
In the following, we prove that there existsδ>0 such that
(36)Re{λi1}≤-δ,Re{λi2}≤-δ.
Letλ=μiξ; then
(37)φ~i(λ)≜μi2ξ2-tr(Mi)μiξ+det(Mi).
Since μi→∞ as i→∞, it follows that
(38)limi→∞φ~i(λ)μi2=ξ2+(d1+d2)ξ+d1d2.
By the Routh-Hurwitz criterion, it follows that the two rootsξ1, ξ2 of φ~i(λ)=0 all have negative real parts. Thus, there exists a positive constant d~, such that Re{ξ1},Re{ξ2}≤-d~. By continuity, we see that there exists i0 such that the two roots ξi1, ξi2 of φ~i(λ)=0 satisfy Re{ξi1}≤-d~/2,Re{ξi2}≤-d~/2, ∀i≥i0. In turn, Re{λi1},Re{λi2}≤-μid~/2≤-d~/2, ∀i≥i0. Let
(39)-δ~=max1≤i≤i0{Re{λi1},Re{λi2}}.
Then δ~>0 and (36) holds for δ=min{δ~,(d~/2)}.
Consequently, the spectrum of£, which consists of eigenvalues, lies in {Reλ≤-δ}. In the sense of [46], we obtain that the positive constant steady-state solution E3=(u*,u*+k2) of model (12) is uniformly asymptotically stable. This ends the proof.


## 2.1. The Existence of the Constant Steady-State Solution

It is easy to verify that model (12) has the following nonnegative constant steady-state solutions:(i)
the trivial constant solutionE0=(0,0) (extinction of two species);(ii)
the semitrivial constant solutionE1=(0,k2) (extinction of the prey);(iii)
the semitrivial constant solutionE2=((1-b+(1-b)2-4(m-b))/2,0) (extinction of the predator);(iv)
the unique positive constant solutionE3=(u*,v*) (coexistence of two species), where v*=u*+k2, and u* is a real positive root of the cubic
(15)ψ(w)=w3+3η1w2+3η2w+η3=0,
where 3η1=b+c+k1-1,3η2=bc+bk1+ck2+m-b-k1,η3=bck2+k1(m-b).By the transformationz=w+η1, (15) is reduced to
(16)h(z)=z3+3pz+q=0,
where
(17)p=η2-η12,(18)q=η3-3η1η2+2η13.It is worthy to note that ifq=η3-3η1η2+2η13<0 holds, (16) has a real positive root. Considering (17) and (18), from q=0, one can determine m=m*, where
(19)m*≜2b-c-k1+19(b+c-2k1-1)×(b2+b-bc-bk1-2c2+4c-4ck1-2k12-5k1+9ck2-2).Hence, we have the following lemma regarding the existence of the positive constant steady-state solution of model (12).Lemma 1.
If either of the following inequalities holds:(20)2k1+1-c-b>0,0<m<min{b,m*};2k1+1-c-b<0,0<m*<m<b,
model (12) has a unique positive constant steady-state solution E3=(u*,v*).


## 2.2. Stability of the Constant Steady-State Solution

In this subsection, we will analyze the asymptotical stability of the nonnegative constant solutions for the corresponding reaction-diffusion dynamics (12).For sake of simplicity, we rewrite model (12) in the vectorial form:
(21)wt=DΔw+H(w)x∈Ω,t>0∂w∂ν=0,x∈∂Ω,t>0,w(x,0)=(u0(x),v0(x)T),x∈Ω,
where w=(u,v)T, D=diag(d1,d2), and
(22)H(w)=(u(1-u-mu+b-cvu+k1)sv(1-vu+k2)).Let0=μ0<μ1<μ2<⋯ be the eigenvalues of the operator-Δ on Ω with the zero-flux boundary conditions. And set
(23)X={w∈[H2(Ω)]2∣∂νw=0on∂Ω},E(μ)={ϕ∣-Δϕ=μϕinΩ,∂νϕ=0on∂Ω},withμ∈ℝ1,{ϕij∣j=1,…,dimE(μi)} being an orthonormal basis of E(μi), and Xij={cϕij∣c∈ℝ2}; then
(24)X=⨁i=1∞Xi,
where Xi=⨁j=1dimE(μi)Xij.LetE=(u,v) be any arbitrary constant steady-state solution of model (12). And the linearization of model (12) at the constant steady-state solution can be expressed by
(25)wt=£(w)=DΔw+J(u,v)w,
where J(u,v)=(ξ11ξ12ξ21ξ22) and
(26)ξ11=1-u-mu+b+u(-1+m(u+b)2)-cvk1+u+cuv(k1+u)2,ξ12=-cuk1+u,ξ21=sv2(u+k2)2,ξ22=s(1-vu+k2)-svu+k2.From [46], it is known that if all the eigenvalues of the operator £ have negative real parts, then E=(u,v) is asymptotically stable; if there is an eigenvalue with positive real part, then E=(u,v) is unstable; if all the eigenvalues have nonpositive real parts while some eigenvalues have zero real parts, then the stability of E=(u,v) cannot be determined by the linearization [10].For eachi≥0, Xi is invariant under the operator £, and λ is an eigenvalue of £ if and only if λ is an eigenvalue of the matrix Mi=-μiD+J(u,v) for some i≥0.In the following, we denoteξjj[k]=ξjj evaluated at Ek, j=1,2 and k=0,1,2,3. So, the local stability of the constant steady-state solution can be analyzed as follows.Theorem 2.
For any positive parameters,(a)
the trivial constant solutionE0=(0,0) is unstable;(b)
the semitrivial constant solutionE1=(0,k2) is
(b1)
locally asymptotically stable ifc>(k1/k2) or (k1/k2)(1-m/b)<c≤(k1/k2) holds(b2)
unstable ifc<(k1/k2)(1-m/b) holds;(c)
the semitrivial constant solutionE2=((1-b+(1-b)2-4(m-b))/2,0) is unstable.Proof.
The stability of the constant steady-state solution is reduced to consider the characteristic equation:(27)det(λI-Mi)=λ2-tr(Mi)λ+det(Mi),
with
(28)tr(Mi)=-μi(d1+d2)+ξ11+ξ22,det(Mi)=d1d2μi2-(ξ11d2+ξ22d1)+detJ(u,v).(a)
JE0=(1-(b/m)00s), for i=0; the eigenvalues are 1-(b/m)>0 and s>0, so E0=(0,0) is unstable.(b)
We can obtainJE1=(1-(b/m)-(ck2/k1)0s-s).
(b1)
Ifk1<ck2 or k1≥ck2andb(1-(ck2/k1))<m<b hold, then 1-(b/m)-(ck2/k1)<0, so for i≥0,
(29)tr(Mi)=-μi(d1+d2)+1-bm-ck2k1-s<0,det(Mi)=d1d2μi2-(1-bm-ck2k1d2-sd1)-s(1-bm-ck2k1)>0.Hence,E1 is locally asymptotically stable.(b2)
Whenk1≥ck2and0<m<b(1-(ck2/k1)), then 1-(b/m)-(ck2/k1)>0. For i≥0, det(Mi)=-s(1-(b/m)-(ck2/k1))<0, which implies that (27) has at least one root with positive real part. Hence, E1 is an unstable steady-state solution of model (12).(c)
Sinceξ11[2]=-((2(1+b)2-4m((1+b)2-4m+1+b-2m))/(b+1+(1+b)2-4m)2),ξ22[2]=s, for i=0, one of the eigenvalues is s>0, so E2 is unstable.
The proof is complete.Straightforward calculations show that(30)JE3=(ξ11[3]-cu*u*+k1s-s),
where ξ11[3]=((cu*(u*+k2))/(u*+k1)2)+(mu*/(u*+b)2)-u*.The determinant ofJE3 is given by(31)det(JE*)=su*((u*+k1)2((u*+b)2-m)+c(u*+b)2(k1-k2))(u*+k1)2(u*+b)2;then, the sign of det(JE3) depends on the factor ϕ(u*):
(32)ϕ(u*)=(u*+k1)2((u*+b)2-m)+c(u*+b)2(k1-k2)=(u*+b+k1+1-c)ψ(u*)+η1u*2+η2u*+η3=ρ1u*2+ρ2u*+ρ3,
where ψ(·) is the same definition as (15), and
(33)ρ1=1+b+k1+c2+ck1+bk1-2m-cb-2ck2-2c,ρ2=-ck2+b-m+2bk1-2cb+b2-ck1+c2k2+c2b,+mc-bm-cb2-4mk1+k12-4cbk2-ck2k1+b2k1+k1+bk12+2bk1c,ρ3=b2k12-bk1c-mbk1+mck1+b2k1c+c2k2b-bk1ck2-cbk2+bk1-mk1+b2k1+bk12-2mk12-2ck2b2.Therefore, we have the following.Theorem 3.
Assume thatξ11[3]<s, ρ1u*2+ρ2u*+ρ3>0, and the first eigenvalue μ1 subject to the zero-flux boundary conditions satisfies
(34)μ1>max{0,ξ11[3]d1-sd2}.
Then the positive constant steady-state solution E3=(u*,u*+k2) of model (15) is uniformly asymptotically stable.Proof.
Whenξ11[3]<s, ρ1u*2+ρ2u*+ρ3>0, then tr(JE3)<0 and det(JE3)>0. So, for i≥0,
(35)tr(Mi)=-(d1+d2)μi+ξ11[3]-sdet(Mi)=μi(d1d2μi-ξ11[3]d2+sd1)+det(JE3).
Note that for any i≥0, we have det(Ai)>0>tr(Ai). Therefore, the eigenvalues of the matrix -μiD+JE3 have negative real parts. It thus follows from the Routh-Hurwitz criterion that, for each i≥0, the two roots λi1 and λi2 of φi(λ)=0 all have negative real parts.
In the following, we prove that there existsδ>0 such that
(36)Re{λi1}≤-δ,Re{λi2}≤-δ.
Letλ=μiξ; then
(37)φ~i(λ)≜μi2ξ2-tr(Mi)μiξ+det(Mi).
Since μi→∞ as i→∞, it follows that
(38)limi→∞φ~i(λ)μi2=ξ2+(d1+d2)ξ+d1d2.
By the Routh-Hurwitz criterion, it follows that the two rootsξ1, ξ2 of φ~i(λ)=0 all have negative real parts. Thus, there exists a positive constant d~, such that Re{ξ1},Re{ξ2}≤-d~. By continuity, we see that there exists i0 such that the two roots ξi1, ξi2 of φ~i(λ)=0 satisfy Re{ξi1}≤-d~/2,Re{ξi2}≤-d~/2, ∀i≥i0. In turn, Re{λi1},Re{λi2}≤-μid~/2≤-d~/2, ∀i≥i0. Let
(39)-δ~=max1≤i≤i0{Re{λi1},Re{λi2}}.
Then δ~>0 and (36) holds for δ=min{δ~,(d~/2)}.
Consequently, the spectrum of£, which consists of eigenvalues, lies in {Reλ≤-δ}. In the sense of [46], we obtain that the positive constant steady-state solution E3=(u*,u*+k2) of model (12) is uniformly asymptotically stable. This ends the proof.


## 3.A Priori Estimates

In this section, we givea priori estimates for the steady-state solutions of model (14). To prove that we recall the following maximum principle [47].Lemma 4 (see [47, maximum principle]).
LetΩ be a bounded Lipschitz domain in ℝn and g∈C(Ω¯×ℝ).(a)
Assume thatw∈C2(Ω)∩C1(Ω¯) and satisfies
(40)Δw(x)+g(x,w(x))≥0inΩ,∂w∂ν≤0on∂Ω.
If w(x0)=maxΩ¯w(x), then g(x0,w(x0))≥0.(b)
Assume thatw∈C2(Ω)∩C1(Ω¯) and satisfies
(41)Δw(x)+g(x,w(x))≤0inΩ,∂w∂ν≥0on∂Ω.
If w(x0)=minΩ¯w(x), then g(x0,w(x0))≤0.Theorem 5.
Suppose thatc(1+k2)<k1(1-(m/b)). Any positive solution (u(x),v(x)) of model (14) satisfies
(42)B<u(x)<AinΩ,(1+k2)B<v(x)<(1+k2)AinΩ,
where
(43)A=12((1-m1+b-k1)2-4(c(B+k2)-k1(1-m1+b))1-m1+b-k1+(1-m1+b-k1)2-4(c(B+k2)-k1(1-m1+b))),B=12((1-mb-k1)2-4(c(1+k2)-k1(1-mb))1-mb-k1+(1-mb-k1)2-4(c(1+k2)-k1(1-mb))).Proof.
Assume that(u(x),v(x)) is a positive solution of (14). Set
(44)u(x10)=maxΩ¯u(x),v(x20)=maxΩ¯v(x),u(x30)=minΩ¯u(x),v(x40)=minΩ¯v(x).
Applying Lemma 4 to model (14), we obtain that
(45)1-u(x10)-mu(x10)+b-cv(x10)u(x10)+k1≥0,hhlh1-v(x20)u(x20)+k2≥0,1-u(x30)-mu(x30)+b-cv(x30)u(x30)+k1≤0,hhlh1-v(x40)u(x40)+k2≤0.
By virtue of the definitions of xi0 (i=1,2,3,4), it follows from (45) that u(x10)<1 and v(x20)<1+k2, and
(46)1-u(x10)-m1+b-cv(x10)u(x10)+k1≥0,1-u(x30)-mb-cv(x30)u(x30)+k1≤0,v(x40)≥u(x40)+k2.
So, we have
(47)u2(x10)-(1-m1+b-k1)u(x10)+cu(x30)+ck2-k1(1-m1+b)≤0,(48)u2(x30)-(1-mb-k1)u(x30)+c(1+k2)-k1(1-mb)≥0.
If c(1+k2)<k1(1-(m/b)) or c(1+k2)=k1(1-(m/b)) and c(1+k2)>k22 hold, from (48), we get that u(x30)>B.
Ifck2<k1(1-(m/(1+b))) or ck2=k1(1-(m/(1+b))) and ck2>k22 hold, from (47), we get that
(49)u(x10)≤12((1-m1+b-k1)2-4(c(u(x3)+k2)-k1(1-m1+b))1-m1+b-k1hh+(1-m1+b-k1)2-4(c(u(x3)+k2)-k1(1-m1+b)))<A.
By simple computations, c(1+k2)<k1(1-(m/b)) indicates that ck2<k1(1-(m/(1+b))). So, if c(1+k2)<k1(1-(m/b)) holds, we can obtain B<u(x)<A and (1+k2)B<v(x)<(1+k2)A. The proof is complete.In order to obtain the desired bounds, we need to use the following Harnack inequality due to [48].Lemma 6 (see [48, Harnack inequality]).
Letw∈C2(Ω)∩C1(Ω¯) be a positive solution to Δw(x)+c(x)w(x)=0, where c∈C(Ω¯), satisfying zero-flux boundary conditions. Then there exists a positive constant C*=C*(∥c∥∞,Ω), such that
(50)maxΩ¯w≤C*minΩ¯w.Theorem 7.
Letd* be an arbitrary fixed positive number. Then there exist positive C_=C_(m,b,c,k1,k2,Ω), such that if d1≥d* and d2>0, any positive solution (u,v) of model (14) satisfies
(51)C_<u(x)<1,C_<v(x)<1+k2.Proof.
By Theorem5, we note that u(x)<1 and v(x)<1+k2. And it suffices to verify the lower bounds of (u,v). We will verify the conclusion by a contradiction argument.
On the contrary, suppose that the conclusion is not true; then there exist sequences{d1,i}i=1∞ and {d2,i}i=1∞ with d1,i≥d,d2,i>0 and the positive solution (ui,vi) of model (14) corresponding to (d1,d2)=(d1,i,d2,i), such that
(52)minΩ¯ui(x)⟶0orminΩ¯vi(x)⟶0asi⟶∞,
and (ui,vi) satisfies
(53)-d1,iΔui=ui(1-ui-mui+b-cvui+k1),(x,y)∈Ω,-d2,iΔvi=svi(1-viui+k2),(x,y)∈Ω,∂ui∂ν=∂vi∂ν=0,(x,y)∈∂Ω.
We observe that Lemma 4 guarantees
(54)k2minΩ¯ui(x)≤minΩ¯vi(x)≤maxΩ¯vi(x)≤k2+maxΩ¯ui(x)
for all i≥1 by use of the second equation of (53). On the other hand, by the Harnack inequality (Lemma 6), we know that there is a positive constant C independent of i, such that maxΩ¯ui(x)≤CminΩ¯ui(x) for all i≥1. Consequently,
(55)maxΩ¯ui(x)⟶0asi⟶∞,
which contradicts Theorem 5. We finish the proof.


## 4. Nonexistence and Existence of the Nonconstant Solutions

In this section, we apply the energy method and the topological degree theory [49] to establish some results on the nonexistence and existence of the positive nonconstant solutions of model (14), respectively.


### 4.1. Nonexistence of the Nonconstant Solutions

#### 4.1.1. Nonexistence of Nonconstant Semitrivial Steady-State Solution

In the case that model (14) has the semitrivial steady-state solutions (u(x),0) and (0,v(x)), such (u(x),0) and (0,v(x)) should satisfy, respectively,
(56)d1Δu+u(1-u-mu+b)=0,x∈Ω,∂u∂ν=0,x∈∂Ω,(57)d2Δv+sv(1-vk2)=0,x∈Ω,∂v∂ν=0,x∈∂Ω.Note thatμ1 is the smallest positive eigenvalue of the operator-Δ in Ω subject to the zero-flux boundary conditions. Now, using the energy estimates, for the semitrivial solution, we can claim the following.Theorem 8.
(i) Ifd1>1/μ1, the only nonnegative solutions of model (56) are u=0 and u=((1-b+(1-b)2-4(m-b))/2).
(ii) Ifd2>s/μ1, the only nonnegative solution of model (57) is v=k2.Proof.
Letu(x),v(x) be a nonnegative solution of models (56) and (57), respectively. Denote u¯=|Ω|-1∫Ω‍udx and v¯=|Ω|-1∫Ω‍vdx. Then ∫Ω‍(u-u¯)dx=∫Ω‍(v-v¯)dx=0.
Multiplying the equation in (56) by u-u¯, we get
(58)d1∫Ω‍|∇(u-u¯)|2dx=∫Ω‍(u-u¯)u(1-u-mu+b)dx=∫Ω‍(u-u¯)2(1-(u+u¯)-m(u¯+b)(u+b))dx≤∫Ω‍(u-u¯)2dx.
Then with the Poincaré inequality
(59)∫Ω‍(u-u¯)2dx≤1μ1∫Ω‍|∇(u-u¯)|2dx,
we find that
(60)d1∫Ω‍|∇(u-u¯)|2dx≤1μ1∫Ω‍|∇(u-u¯)|2dx,
which implies that d1<1/μ1 unless u-u¯≡0, and u must be a constant solution.
In a similar manner, we multiply the equation in (57) by (v-v¯) to have
(61)d2∫Ω‍|∇(v-v¯)|2dx=∫Ω‍s(1-vk2)v(v-v¯)dx=∫Ω‍s(1-1k2(v+v¯))(v-v¯)2dx≤s∫Ω‍(v-v¯)2dx≤sμ1∫Ω‍|∇(v-v¯)|2dx.
In view of d2>s/μ1, we have v-v¯≡0 and v must be a constant solution. This ends the proof.


#### 4.1.2. Nonexistence of the Positive Nonconstant Solutions

This subsection is devoted to the consideration of the nonexistence for the positive nonconstant solutions of model (12), and in the below results, the diffusion coefficients do play a significant role.Theorem 9.
There exists a positive constantd*=d*(b,m,c,s,k1,k2,Ω) such that model (14) has no positive nonconstant solution provided that min{d1,d2}>d*.Proof.
Let(u,v) be any positive solution of model (14). Then, multiplying the first equation of model (14) by (u-u¯), integrating over Ω and using Theorem 7, we have that
(62)d1∫Ω‍|∇(u-u¯)|2dx=∫Ω‍(u-u¯)u(1-u-mu+b)dx-∫Ω‍cuv(u-u¯)u+k1dx=∫Ω‍(u-u¯)u(1-u-mu+b)dx-∫Ω‍cv(u-u¯)2u+k1dx-∫Ω‍cu¯v(u-u¯)u+k1dx≤∫Ω‍(u-u¯)2dx+∫Ω‍-cu¯v(u-u¯)u+k1dx.
In a similar manner, we multiply the second equation in model (14) by (v-v¯) to have
(63)d2∫Ω‍|∇(v-v¯)|2dx=∫Ω‍s(1-vu+k2)v(v-v¯)dx=∫Ω‍s(1-vu+k2)(v-v¯)2dx+∫Ω‍s(1-vu+k2)v¯(v-v¯)dx≤∫Ω‍s(v-v¯)2dx+∫Ω‍s(1-vu+k2)v¯(v-v¯)dx.
Using Theorem 7 again, we have that
(64)∫Ω‍-cu¯v(u-u¯)u+k1dx=∫Ω‍cu¯(v¯u¯+k1-vu+k1)(u-u¯)dx=∫Ω‍cu¯(u-u¯)(u¯+k1)(u+k1)(v¯(u-u¯)-(u¯+k1)(v-v¯))dx≤c(1+k2)k1∫Ω‍(u-u¯)2dx+ck1∫Ω‍|u-u¯||v-u¯|dx≤c(3+2k2)2k1∫Ω‍(u-u¯)2dx+c2k1∫Ω‍(v-v¯)2dx.
Similarly,
(65)∫Ω‍s(1-vu+k2)v¯(v-v¯)dx=∫Ω‍sv¯(v¯u¯+k2-vu+k2)(v-v¯)dx=∫Ω‍sv¯(v-v¯)(u¯+k2)(u+k2)(v¯(u-u¯)-(u¯+k2)(v-v¯))dx≤sv¯2(u¯+k1)(u+k1)∫Ω‍|u-u¯||v-u¯|dx≤s(1+k2)22k12∫Ω‍(u-u¯)2dx+s(1+k2)22k12∫Ω‍(v-v¯)2dx.
From the well-known Poincaré inequality (62), (63), (64), and (65), we obtain that
(66)d1∫Ω‍|∇(u-u¯)|2dx+d2∫Ω‍|∇(v-v¯)|2dx≤1μ1(α1∫Ω‍|∇(u-u¯)|2dx+α2∫Ω‍|∇(v-v¯)|2dx),
where α1=1+(c(3+2k2)/2k1)+(s(1+k2)2/2k12), α2=s+(c/2k1)+(s(1+k2)2/2k12).
This shows that if(67)min{d1,d2}>1μ1max{α1,α2},
then ∇(u-u¯)=∇(v-v¯)=0, which asserts our results.


### 4.2. Existence of the Nonconstant Positive Solutions

In this subsection, we will discuss the existence of the positive nonconstant solution of model (14).Unless otherwise specified, in this subsection, we always require that Lemma1 holds, which guarantees that model (14) has the unique positive constant solution E3=(u*,u*+k2). From now on, we denote w0=E3.LetX be the space defined in (23) and let
(68)X+={(u,v)∈Xu,v>0onΩ¯},B(C)={(u,v)∈XC-1<u,v<ConΩ¯},C>0.
We write model (14) in the form:
(69)-Δw=𝒢(w),w∈X+,∂νw=0on∂Ω,
where
(70)𝒢(w)=(ud1(1-u-mu+b-cvu+k1)svd2(1-vu+k2)).
Then w is a positive solution of model (69) if and only if  w satisfies
(71)ℱ(w)=w-(I-Δ)-1{𝒢(w)+w}=0,inX+,
where (I-Δ)-1 is the inverse operator of I-Δ subject to the zero-flux boundary condition. Then
(72)∇ℱ(w0)=I-(I-Δ)-1(I+𝒜),
where
(73)𝒜≜∇𝒢(w0)=(ξ11[3]d1-cu*d1(u*+k1)sd2-sd2).If∇ℱ(w0) is invertible, by Theorem 2.8.1 of [50] the index of ℱ at w0 is given by
(74)index(ℱ,w0)=(-1)γ,
where γ is the multiplicity of negative eigenvalues of ∇ℱ(w0).On the other hand, using the decomposition (24), we have that Xi is an invariant space under ∇ℱ(w0) and ξ∈ℝ is an eigenvalue of ∇ℱ(w0) in Xi, if and only if ξ is an eigenvalue of (μi+1)-1(μiI-𝒜). Therefore, ∇ℱ(w0) is invertible, if and only if for any i≥0 the matrix μiI-𝒜 is invertible.Letm(ui) be the multiplicity of μi. For the sake of convenience, we denote
(75)H(d2,d2,μ)=det(μI-𝒜).
Then, if μiI-𝒜 is invertible for any i≥0, with the same arguments as in [51], we can assert the following conclusion.Lemma 10.
Assume that, for alli≥0, the matrix μiI-𝒜 is nonsingular; then
(76)index(ℱ,w0)=(-1)γ,whereγ=∑i≥0,H(ui)<0‍m(ui).To computeindex(ℱ,w0), we have to consider the sign of H(μ). A straightforward computation yields
(77)H(d2,d2,μ)=μ2-θ1μ+θ2,
where θ1=(ξ11[3]/d1)-(s/d2),θ2=(s/d1d2)((cu*/(u*+k1))-ξ11[3]).Ifθ12-4θ2>0, then H(d2,d2,μ)=0 has two positive solutions μ± given by
(78)μ±=12(θ1±θ12-4θ2).Theorem 11.
Assume thatθ12-4θ2>0. If μ-∈(μi,μi+1) and μ+∈(μj,μj+1) for some 0≤i<j, and ∑k=i+1j‍m(uk) is odd, then model (14) has at least one nonconstant solution.Proof.
By Theorem9, we can fixed d¯1>d1 and d¯2>d2 such that(i)
model (12) with diffusion coefficients d¯1 and d¯2 has no nonconstant solutions;(ii)
H(d¯1,d¯2,μ)>0 for all μ≥0.
By virtue of Theorem7, there exists a positive constant M=M(m,b,c,k1,k2) such that M-1<u,v<M.
Set(79)ℳ={(u,v)∈C(Ω¯)×C(Ω¯):M-1<u,v<MinΩ¯},
and define
(80)Ψ:ℳ×[0,1]⟶C(Ω¯)×C(Ω¯)
by
(81)Ψ(w,t)=(I-Δ)-1{G(w,t)+w},
where
(82)G(w,t)=((td1+(1-t)d¯1)-1(u(1-u-mu+b-cvu+k1))(td2+(1-t)d¯2)-1(sv(1-vu+k2))).
It is clear that solving model (14) is equivalent to finding a fixed point of Ψ(w,1) in ℳ. Further, by virtue of the definition of ℳ, we have that Ψ(w,1) has no fixed point in ∂ℳ for all 0≤t≤1.
SinceΨ(w,t) is compact, the Leray-Schauder topological degree deg(I-Ψ(w,t),ℳ,0) is well defined. From the invariance of Leray-Schauder degree at the homotopy, we deduce
(83)deg(I-Ψ(w,1),ℳ,0)=deg(I-Ψ(w,0),ℳ,0).
Clearly,I-Ψ(w,1)=ℱ. Thus, if model (14) has no other solutions except the constant one w0, then Lemma 10 shows that
(84)deg(I-Ψ(w,1),ℳ,0)=index(ℱ,w0)=(-1)∑k=i+1jm(uk)=-1.
On the contrary, by the choice ofd¯1 and d¯2, we have that w0 is the only solution of Ψ(w,0)=w. Furthermore, by (ii) above, we have
(85)deg(I-Ψ(w,0),ℳ,0)=index(I-Ψ(w,0),w0)=1.
From (83)–(85), we get a contradiction, and the proof is completed.Corollary 12.
Letd2 be fixed. If ξ11[3]>0 and all the eigenvalues μi have odd multiplicity, then, there exists a sequence of intervals {(kn,Kn)} with 0<kn<Kn<kn+1→0 (as n→∞) such that the steady-state model (14) has at least one nonconstant solution for all d1∈⋃n≥1‍(kn,Kn).Corollary 13.
Letd1 be fixed. If ξ11[3]>0 and ∑i≥0,0<μi<ξ11[3]‍m(μi) is odd, then there exists D*>0 such that the steady-state model (14) has at least one nonconstant solution for any d2>D*.We omit the proofs of Corollaries12 and 13 here and refer the reader to more detailed proofs in [52].Remark 14.
Our results suggest that if the parameters are properly chosen, both the general stationary pattern and more interesting Turing pattern can arise as a result of diffusion.


## 4.1. Nonexistence of the Nonconstant Solutions

### 4.1.1. Nonexistence of Nonconstant Semitrivial Steady-State Solution

In the case that model (14) has the semitrivial steady-state solutions (u(x),0) and (0,v(x)), such (u(x),0) and (0,v(x)) should satisfy, respectively,
(56)d1Δu+u(1-u-mu+b)=0,x∈Ω,∂u∂ν=0,x∈∂Ω,(57)d2Δv+sv(1-vk2)=0,x∈Ω,∂v∂ν=0,x∈∂Ω.Note thatμ1 is the smallest positive eigenvalue of the operator-Δ in Ω subject to the zero-flux boundary conditions. Now, using the energy estimates, for the semitrivial solution, we can claim the following.Theorem 8.
(i) Ifd1>1/μ1, the only nonnegative solutions of model (56) are u=0 and u=((1-b+(1-b)2-4(m-b))/2).
(ii) Ifd2>s/μ1, the only nonnegative solution of model (57) is v=k2.Proof.
Letu(x),v(x) be a nonnegative solution of models (56) and (57), respectively. Denote u¯=|Ω|-1∫Ω‍udx and v¯=|Ω|-1∫Ω‍vdx. Then ∫Ω‍(u-u¯)dx=∫Ω‍(v-v¯)dx=0.
Multiplying the equation in (56) by u-u¯, we get
(58)d1∫Ω‍|∇(u-u¯)|2dx=∫Ω‍(u-u¯)u(1-u-mu+b)dx=∫Ω‍(u-u¯)2(1-(u+u¯)-m(u¯+b)(u+b))dx≤∫Ω‍(u-u¯)2dx.
Then with the Poincaré inequality
(59)∫Ω‍(u-u¯)2dx≤1μ1∫Ω‍|∇(u-u¯)|2dx,
we find that
(60)d1∫Ω‍|∇(u-u¯)|2dx≤1μ1∫Ω‍|∇(u-u¯)|2dx,
which implies that d1<1/μ1 unless u-u¯≡0, and u must be a constant solution.
In a similar manner, we multiply the equation in (57) by (v-v¯) to have
(61)d2∫Ω‍|∇(v-v¯)|2dx=∫Ω‍s(1-vk2)v(v-v¯)dx=∫Ω‍s(1-1k2(v+v¯))(v-v¯)2dx≤s∫Ω‍(v-v¯)2dx≤sμ1∫Ω‍|∇(v-v¯)|2dx.
In view of d2>s/μ1, we have v-v¯≡0 and v must be a constant solution. This ends the proof.


### 4.1.2. Nonexistence of the Positive Nonconstant Solutions

This subsection is devoted to the consideration of the nonexistence for the positive nonconstant solutions of model (12), and in the below results, the diffusion coefficients do play a significant role.Theorem 9.
There exists a positive constantd*=d*(b,m,c,s,k1,k2,Ω) such that model (14) has no positive nonconstant solution provided that min{d1,d2}>d*.Proof.
Let(u,v) be any positive solution of model (14). Then, multiplying the first equation of model (14) by (u-u¯), integrating over Ω and using Theorem 7, we have that
(62)d1∫Ω‍|∇(u-u¯)|2dx=∫Ω‍(u-u¯)u(1-u-mu+b)dx-∫Ω‍cuv(u-u¯)u+k1dx=∫Ω‍(u-u¯)u(1-u-mu+b)dx-∫Ω‍cv(u-u¯)2u+k1dx-∫Ω‍cu¯v(u-u¯)u+k1dx≤∫Ω‍(u-u¯)2dx+∫Ω‍-cu¯v(u-u¯)u+k1dx.
In a similar manner, we multiply the second equation in model (14) by (v-v¯) to have
(63)d2∫Ω‍|∇(v-v¯)|2dx=∫Ω‍s(1-vu+k2)v(v-v¯)dx=∫Ω‍s(1-vu+k2)(v-v¯)2dx+∫Ω‍s(1-vu+k2)v¯(v-v¯)dx≤∫Ω‍s(v-v¯)2dx+∫Ω‍s(1-vu+k2)v¯(v-v¯)dx.
Using Theorem 7 again, we have that
(64)∫Ω‍-cu¯v(u-u¯)u+k1dx=∫Ω‍cu¯(v¯u¯+k1-vu+k1)(u-u¯)dx=∫Ω‍cu¯(u-u¯)(u¯+k1)(u+k1)(v¯(u-u¯)-(u¯+k1)(v-v¯))dx≤c(1+k2)k1∫Ω‍(u-u¯)2dx+ck1∫Ω‍|u-u¯||v-u¯|dx≤c(3+2k2)2k1∫Ω‍(u-u¯)2dx+c2k1∫Ω‍(v-v¯)2dx.
Similarly,
(65)∫Ω‍s(1-vu+k2)v¯(v-v¯)dx=∫Ω‍sv¯(v¯u¯+k2-vu+k2)(v-v¯)dx=∫Ω‍sv¯(v-v¯)(u¯+k2)(u+k2)(v¯(u-u¯)-(u¯+k2)(v-v¯))dx≤sv¯2(u¯+k1)(u+k1)∫Ω‍|u-u¯||v-u¯|dx≤s(1+k2)22k12∫Ω‍(u-u¯)2dx+s(1+k2)22k12∫Ω‍(v-v¯)2dx.
From the well-known Poincaré inequality (62), (63), (64), and (65), we obtain that
(66)d1∫Ω‍|∇(u-u¯)|2dx+d2∫Ω‍|∇(v-v¯)|2dx≤1μ1(α1∫Ω‍|∇(u-u¯)|2dx+α2∫Ω‍|∇(v-v¯)|2dx),
where α1=1+(c(3+2k2)/2k1)+(s(1+k2)2/2k12), α2=s+(c/2k1)+(s(1+k2)2/2k12).
This shows that if(67)min{d1,d2}>1μ1max{α1,α2},
then ∇(u-u¯)=∇(v-v¯)=0, which asserts our results.


## 4.1.1. Nonexistence of Nonconstant Semitrivial Steady-State Solution

In the case that model (14) has the semitrivial steady-state solutions (u(x),0) and (0,v(x)), such (u(x),0) and (0,v(x)) should satisfy, respectively,
(56)d1Δu+u(1-u-mu+b)=0,x∈Ω,∂u∂ν=0,x∈∂Ω,(57)d2Δv+sv(1-vk2)=0,x∈Ω,∂v∂ν=0,x∈∂Ω.Note thatμ1 is the smallest positive eigenvalue of the operator-Δ in Ω subject to the zero-flux boundary conditions. Now, using the energy estimates, for the semitrivial solution, we can claim the following.Theorem 8.
(i) Ifd1>1/μ1, the only nonnegative solutions of model (56) are u=0 and u=((1-b+(1-b)2-4(m-b))/2).
(ii) Ifd2>s/μ1, the only nonnegative solution of model (57) is v=k2.Proof.
Letu(x),v(x) be a nonnegative solution of models (56) and (57), respectively. Denote u¯=|Ω|-1∫Ω‍udx and v¯=|Ω|-1∫Ω‍vdx. Then ∫Ω‍(u-u¯)dx=∫Ω‍(v-v¯)dx=0.
Multiplying the equation in (56) by u-u¯, we get
(58)d1∫Ω‍|∇(u-u¯)|2dx=∫Ω‍(u-u¯)u(1-u-mu+b)dx=∫Ω‍(u-u¯)2(1-(u+u¯)-m(u¯+b)(u+b))dx≤∫Ω‍(u-u¯)2dx.
Then with the Poincaré inequality
(59)∫Ω‍(u-u¯)2dx≤1μ1∫Ω‍|∇(u-u¯)|2dx,
we find that
(60)d1∫Ω‍|∇(u-u¯)|2dx≤1μ1∫Ω‍|∇(u-u¯)|2dx,
which implies that d1<1/μ1 unless u-u¯≡0, and u must be a constant solution.
In a similar manner, we multiply the equation in (57) by (v-v¯) to have
(61)d2∫Ω‍|∇(v-v¯)|2dx=∫Ω‍s(1-vk2)v(v-v¯)dx=∫Ω‍s(1-1k2(v+v¯))(v-v¯)2dx≤s∫Ω‍(v-v¯)2dx≤sμ1∫Ω‍|∇(v-v¯)|2dx.
In view of d2>s/μ1, we have v-v¯≡0 and v must be a constant solution. This ends the proof.


## 4.1.2. Nonexistence of the Positive Nonconstant Solutions

This subsection is devoted to the consideration of the nonexistence for the positive nonconstant solutions of model (12), and in the below results, the diffusion coefficients do play a significant role.Theorem 9.
There exists a positive constantd*=d*(b,m,c,s,k1,k2,Ω) such that model (14) has no positive nonconstant solution provided that min{d1,d2}>d*.Proof.
Let(u,v) be any positive solution of model (14). Then, multiplying the first equation of model (14) by (u-u¯), integrating over Ω and using Theorem 7, we have that
(62)d1∫Ω‍|∇(u-u¯)|2dx=∫Ω‍(u-u¯)u(1-u-mu+b)dx-∫Ω‍cuv(u-u¯)u+k1dx=∫Ω‍(u-u¯)u(1-u-mu+b)dx-∫Ω‍cv(u-u¯)2u+k1dx-∫Ω‍cu¯v(u-u¯)u+k1dx≤∫Ω‍(u-u¯)2dx+∫Ω‍-cu¯v(u-u¯)u+k1dx.
In a similar manner, we multiply the second equation in model (14) by (v-v¯) to have
(63)d2∫Ω‍|∇(v-v¯)|2dx=∫Ω‍s(1-vu+k2)v(v-v¯)dx=∫Ω‍s(1-vu+k2)(v-v¯)2dx+∫Ω‍s(1-vu+k2)v¯(v-v¯)dx≤∫Ω‍s(v-v¯)2dx+∫Ω‍s(1-vu+k2)v¯(v-v¯)dx.
Using Theorem 7 again, we have that
(64)∫Ω‍-cu¯v(u-u¯)u+k1dx=∫Ω‍cu¯(v¯u¯+k1-vu+k1)(u-u¯)dx=∫Ω‍cu¯(u-u¯)(u¯+k1)(u+k1)(v¯(u-u¯)-(u¯+k1)(v-v¯))dx≤c(1+k2)k1∫Ω‍(u-u¯)2dx+ck1∫Ω‍|u-u¯||v-u¯|dx≤c(3+2k2)2k1∫Ω‍(u-u¯)2dx+c2k1∫Ω‍(v-v¯)2dx.
Similarly,
(65)∫Ω‍s(1-vu+k2)v¯(v-v¯)dx=∫Ω‍sv¯(v¯u¯+k2-vu+k2)(v-v¯)dx=∫Ω‍sv¯(v-v¯)(u¯+k2)(u+k2)(v¯(u-u¯)-(u¯+k2)(v-v¯))dx≤sv¯2(u¯+k1)(u+k1)∫Ω‍|u-u¯||v-u¯|dx≤s(1+k2)22k12∫Ω‍(u-u¯)2dx+s(1+k2)22k12∫Ω‍(v-v¯)2dx.
From the well-known Poincaré inequality (62), (63), (64), and (65), we obtain that
(66)d1∫Ω‍|∇(u-u¯)|2dx+d2∫Ω‍|∇(v-v¯)|2dx≤1μ1(α1∫Ω‍|∇(u-u¯)|2dx+α2∫Ω‍|∇(v-v¯)|2dx),
where α1=1+(c(3+2k2)/2k1)+(s(1+k2)2/2k12), α2=s+(c/2k1)+(s(1+k2)2/2k12).
This shows that if(67)min{d1,d2}>1μ1max{α1,α2},
then ∇(u-u¯)=∇(v-v¯)=0, which asserts our results.


## 4.2. Existence of the Nonconstant Positive Solutions

In this subsection, we will discuss the existence of the positive nonconstant solution of model (14).Unless otherwise specified, in this subsection, we always require that Lemma1 holds, which guarantees that model (14) has the unique positive constant solution E3=(u*,u*+k2). From now on, we denote w0=E3.LetX be the space defined in (23) and let
(68)X+={(u,v)∈Xu,v>0onΩ¯},B(C)={(u,v)∈XC-1<u,v<ConΩ¯},C>0.
We write model (14) in the form:
(69)-Δw=𝒢(w),w∈X+,∂νw=0on∂Ω,
where
(70)𝒢(w)=(ud1(1-u-mu+b-cvu+k1)svd2(1-vu+k2)).
Then w is a positive solution of model (69) if and only if  w satisfies
(71)ℱ(w)=w-(I-Δ)-1{𝒢(w)+w}=0,inX+,
where (I-Δ)-1 is the inverse operator of I-Δ subject to the zero-flux boundary condition. Then
(72)∇ℱ(w0)=I-(I-Δ)-1(I+𝒜),
where
(73)𝒜≜∇𝒢(w0)=(ξ11[3]d1-cu*d1(u*+k1)sd2-sd2).If∇ℱ(w0) is invertible, by Theorem 2.8.1 of [50] the index of ℱ at w0 is given by
(74)index(ℱ,w0)=(-1)γ,
where γ is the multiplicity of negative eigenvalues of ∇ℱ(w0).On the other hand, using the decomposition (24), we have that Xi is an invariant space under ∇ℱ(w0) and ξ∈ℝ is an eigenvalue of ∇ℱ(w0) in Xi, if and only if ξ is an eigenvalue of (μi+1)-1(μiI-𝒜). Therefore, ∇ℱ(w0) is invertible, if and only if for any i≥0 the matrix μiI-𝒜 is invertible.Letm(ui) be the multiplicity of μi. For the sake of convenience, we denote
(75)H(d2,d2,μ)=det(μI-𝒜).
Then, if μiI-𝒜 is invertible for any i≥0, with the same arguments as in [51], we can assert the following conclusion.Lemma 10.
Assume that, for alli≥0, the matrix μiI-𝒜 is nonsingular; then
(76)index(ℱ,w0)=(-1)γ,whereγ=∑i≥0,H(ui)<0‍m(ui).To computeindex(ℱ,w0), we have to consider the sign of H(μ). A straightforward computation yields
(77)H(d2,d2,μ)=μ2-θ1μ+θ2,
where θ1=(ξ11[3]/d1)-(s/d2),θ2=(s/d1d2)((cu*/(u*+k1))-ξ11[3]).Ifθ12-4θ2>0, then H(d2,d2,μ)=0 has two positive solutions μ± given by
(78)μ±=12(θ1±θ12-4θ2).Theorem 11.
Assume thatθ12-4θ2>0. If μ-∈(μi,μi+1) and μ+∈(μj,μj+1) for some 0≤i<j, and ∑k=i+1j‍m(uk) is odd, then model (14) has at least one nonconstant solution.Proof.
By Theorem9, we can fixed d¯1>d1 and d¯2>d2 such that(i)
model (12) with diffusion coefficients d¯1 and d¯2 has no nonconstant solutions;(ii)
H(d¯1,d¯2,μ)>0 for all μ≥0.
By virtue of Theorem7, there exists a positive constant M=M(m,b,c,k1,k2) such that M-1<u,v<M.
Set(79)ℳ={(u,v)∈C(Ω¯)×C(Ω¯):M-1<u,v<MinΩ¯},
and define
(80)Ψ:ℳ×[0,1]⟶C(Ω¯)×C(Ω¯)
by
(81)Ψ(w,t)=(I-Δ)-1{G(w,t)+w},
where
(82)G(w,t)=((td1+(1-t)d¯1)-1(u(1-u-mu+b-cvu+k1))(td2+(1-t)d¯2)-1(sv(1-vu+k2))).
It is clear that solving model (14) is equivalent to finding a fixed point of Ψ(w,1) in ℳ. Further, by virtue of the definition of ℳ, we have that Ψ(w,1) has no fixed point in ∂ℳ for all 0≤t≤1.
SinceΨ(w,t) is compact, the Leray-Schauder topological degree deg(I-Ψ(w,t),ℳ,0) is well defined. From the invariance of Leray-Schauder degree at the homotopy, we deduce
(83)deg(I-Ψ(w,1),ℳ,0)=deg(I-Ψ(w,0),ℳ,0).
Clearly,I-Ψ(w,1)=ℱ. Thus, if model (14) has no other solutions except the constant one w0, then Lemma 10 shows that
(84)deg(I-Ψ(w,1),ℳ,0)=index(ℱ,w0)=(-1)∑k=i+1jm(uk)=-1.
On the contrary, by the choice ofd¯1 and d¯2, we have that w0 is the only solution of Ψ(w,0)=w. Furthermore, by (ii) above, we have
(85)deg(I-Ψ(w,0),ℳ,0)=index(I-Ψ(w,0),w0)=1.
From (83)–(85), we get a contradiction, and the proof is completed.Corollary 12.
Letd2 be fixed. If ξ11[3]>0 and all the eigenvalues μi have odd multiplicity, then, there exists a sequence of intervals {(kn,Kn)} with 0<kn<Kn<kn+1→0 (as n→∞) such that the steady-state model (14) has at least one nonconstant solution for all d1∈⋃n≥1‍(kn,Kn).Corollary 13.
Letd1 be fixed. If ξ11[3]>0 and ∑i≥0,0<μi<ξ11[3]‍m(μi) is odd, then there exists D*>0 such that the steady-state model (14) has at least one nonconstant solution for any d2>D*.We omit the proofs of Corollaries12 and 13 here and refer the reader to more detailed proofs in [52].Remark 14.
Our results suggest that if the parameters are properly chosen, both the general stationary pattern and more interesting Turing pattern can arise as a result of diffusion.


## 5. Turing Instability and Pattern Formation

### 5.1. Turing Instability

In this subsection, we mainly focus on the effects of reaction-diffusion on Turing instability for model (12).Let us consider the spatially homogeneous system corresponding to model (12):
(86)dudt=u(1-u-mu+b-cvu+k1),t>0,dvdt=sv(1-vu+k2),t>0.Mathematically speaking, a positive constant steady-state solutionE3 is Turing instability, which was emphasized by Turing in his pioneering work in 1952 [53], meaning that it is an asymptotically stable steady-state solution of model (86) but is unstable with respect to the solutions of spatial model (12). Since ξ11[33]-s<0, then tr(Ai)<0 (the matrix Ai=-μiD+JE3) is always true. Hence Ai has an eigenvalue with a positive real part; then it must be a real value and the other eigenvalue must be a negative real one. As a consequence, we can obtain the following results.Theorem 15.
Assume that the following conditions are true:(a)
ξ11[3]<s;(b)
η1u*2+η2u*+μ3>0;(c)
d2ξ11[3]-d1s>2d1d2detJE3;
then the positive constant steady-state solutionE3=(u*,u*+k2) of model (12) is Turing instability if 0<k1<μi<k2 for some μi, where
(87)k1=d2ξ11[3]-d1s-(d2ξ11[3]-d1s)2-4d1d2detJE32d1d2,k2=d2ξ11[3]-d1s+(d2ξ11[3]-d1s)2-4d1d2detJE32d1d2.Proof.
Using the same approach as in Theorem3, we have that the positive constant steady state solution E3 of model (86) is asymptotically stable provided (a) and (b).
For the Turing instability, we must havedet(Mi)<0 for some μi. And we notice that det(Mi) achieves its minimum:
(88)minμidet(Mi)=4d1d2det(JE3)-(d2ξ11[3]-d1s)24d1d2
at the critical value μ*>0 when μ*=((d2ξ11[3]-d1s)/2d1d2).
We now discuss the conditions for Turing instability in the further calculation. The conditionminμidet(Mi)<0 is equivalent to 4d1d2det(JEe2)-(d2ξ11[3]-d1s)2<0. In this case, det(Mi)=0 has two positive roots k1 and k2 which satisfy
(89)k1=d2ξ11[3]-d1s-(d2ξ11[3]-d1s)2-4d1d2detJE32d1d2,k2=d2ξ11[3]-d1s+(d2ξ11[3]-d1s)2-4d1d2detJE32d1d2.
Therefore, if we can find some μi such that k1<μi<k2, then det(Ai)<0. By [46] it follows that E3 is uniformly asymptotically instable. This finishes the proof.


### 5.2. Pattern Formation

In this section, we perform extensive numerical simulations of the spatially extended model (12) in two-dimensional space, and the qualitative results are shown here. All our numerical simulations employ the zero-flux boundary conditions with a system size of 100×100.The numerical integration of model (12) is performed by using a finite difference approximation for the spatial derivatives and an explicit Euler method for the time integration [54] with a time stepsize of 1/100. The initial condition is always a small amplitude random perturbation around the positive constant steady-state solution E3=(u*,u*+k2). After the initial period during which the perturbation spread, the model goes either into a time-dependent state, or to an essentially steady-state solution (time independent).We use the standard five-point approximation [55] for the 2D Laplacian with the zero-flux boundary conditions. More precisely, the concentrations (ui,jn+1,vi,jn+1) at the moment (n+1)τ at the mesh position (xi,yj) are given by
(90)ui,jn+1=ui,jn+τd1Δhui,jn+τf(ui,jn,vi,jn),vi,jn+1=vi,jn+τd2Δhvi,jn+τg(ui,jn,vi,jn),
with the Laplacian defined by
(91)Δhui,jn=ui+1,jn+ui-1,jn+ui,j+1n+ui,j-1n-4ui,jnh2,
where f(u,v)=u(1-u-(m/(u+b)))-(cuv/(u+k1)) and g(u,v)=sv(1-(v/(u+k2))), and the space stepsize h=1/3.In the numerical simulations, different types of dynamics are observed and it is found that the distributions of predator and prey are always of the same type. Consequently, we can restrict our analysis of pattern formation to one distribution. In this section, we show the distribution of preyu, for instance. We have taken some snapshots with red (blue) corresponding to the high (low) value of prey u.Now, we show the Turing pattern for the different values of the parameters. Via numerical simulation, one can see that the model dynamics exhibits spatiotemporal complexity of pattern formation, including spots, stripes, and holes Turing patterns.In Figure1, we show the time process of spots pattern formation of the prey u at t=0,300,1000,3000 for b=0.495, c=0.5, k1=0.3, k2=0.2, m=0.32, s=0.25, d1=0.1, and d2=2. In this case, there appears a competition between stripes and spots. The pattern takes a long time to settle down; starting with a homogeneous state E3=(0.066566,0.266566) (cf., Figure 1(a)), the random perturbations lead to the formation of stripes and spots (cf., Figure 1(b)), and the later random perturbations make these stripes decay and end with the time-independent regular spots pattern (cf., Figure 1(d)), which is isolated zones with low prey densities. Ecologically, spots pattern shows that the prey population is driven by predators to a very high level in those regions. The final result is the formation of patches of high prey density surrounded by areas of high prey densities; that is to say, under the control of these parameters, the prey is predominant in the area.The formation of spots patterns of preyu with parameters b=0.495,c=0.5,k1=0.3,k2=0.2,m=0.32,s=0.25,d1=0.1, and d2=2. Time: (a) t=0, (b) t=300, (c) t=1000, and (c) t=3000.
(a)
(b)
(c)
(d)In Figure2, we show the time process of stripes pattern formation of the prey u at t=0,300,1000,3000 for b=0.495,c=0.5,k1=0.3,k2=0.2,m=0.3,s=0.25,d1=0.1,d2=2. In this case, starting with a homogeneous state E3=(0.142318,0.342318) (cf., Figure 2(a)), the random perturbations lead to the formation of stripes and spots (cf., Figure 2(b)), and the later random perturbations make these stripes decay and end with the time-independent regular spots pattern (cf., Figure 2(d)).The formation of stripes patterns of preyu with parameters b=0.495,c=0.5,k1=0.3,k2=0.2,m=0.3,s=0.25,d1=0.1, and d2=2. Time: (a) t=0, (b) t=300, (c) t=1000, and (c) t=3000.
(a)
(b)
(c)
(d)In Figure3, we show the time process of holes pattern formation of the prey u at t=0,300,1000,3000 for b=0.4,c=0.5,k1=0.3,k2=0.2,m=0.24,s=0.25,d1=0.1, and d2=2. In this case, starting with a homogeneous state E3=(0.2,0.4) (cf., Figure 3(a)), the random perturbations lead to the formation of independent regular holes pattern (cf., Figure 1(d)), which is isolated zones with low prey densities.The formation of holes patterns of preyu with parameters b=0.4,c=0.5,k1=0.3,k2=0.2,m=0.24,s=0.25,d1=0.1, and d2=2. Time: (a) t=0, (b) t=300, (c) t=1000, and (c) t=3000.
(a)
(b)
(c)
(d)Comparing Figure1(d) with Figure 3(d), we find that they share similarities. Figure 1(d) consists of red (maximum density of u) spots on a blue (minimum density of u) background; that is, the preys are isolated zones with high population density. Figure 3(d) consists of blue (minimum density of u) spots on a red (maximum density of u) background; that is, the preys are isolated zones with low population density. For the sake of learning the pattern dynamics of model (12) further, we illustrate the three-dimensional patterns in the space (u,x,y). From Figure 4, one can realize the relations of the patterns (e.g., Figures 4(a) and 4(c)) with their corresponding numerical solutions (e.g., Figures 4(b) and 4(d)). In fact, the patterns (e.g., Figures 1–3) are the projections in xy plane of the numerical solutions u(x,y,t) to the model (12). In the software MATLAB, imagesc(u) displays u(x,y;t) as an image; each numerical solution of u(x,y;t) corresponds to a rectangular area in the image; that is, the values of the numerical solutions u(x,y;t) to the model (12) are indices into the current colormap that determine the color of each patch.Pattern dynamics of model (12) in two- and three-dimensional space. (a) spots pattern; (b) spots pattern in three-dimensional space; (c) holes pattern; (d) holes pattern in three-dimensional space.
(a)
(b)
(c)
(d)


## 5.1. Turing Instability

In this subsection, we mainly focus on the effects of reaction-diffusion on Turing instability for model (12).Let us consider the spatially homogeneous system corresponding to model (12):
(86)dudt=u(1-u-mu+b-cvu+k1),t>0,dvdt=sv(1-vu+k2),t>0.Mathematically speaking, a positive constant steady-state solutionE3 is Turing instability, which was emphasized by Turing in his pioneering work in 1952 [53], meaning that it is an asymptotically stable steady-state solution of model (86) but is unstable with respect to the solutions of spatial model (12). Since ξ11[33]-s<0, then tr(Ai)<0 (the matrix Ai=-μiD+JE3) is always true. Hence Ai has an eigenvalue with a positive real part; then it must be a real value and the other eigenvalue must be a negative real one. As a consequence, we can obtain the following results.Theorem 15.
Assume that the following conditions are true:(a)
ξ11[3]<s;(b)
η1u*2+η2u*+μ3>0;(c)
d2ξ11[3]-d1s>2d1d2detJE3;
then the positive constant steady-state solutionE3=(u*,u*+k2) of model (12) is Turing instability if 0<k1<μi<k2 for some μi, where
(87)k1=d2ξ11[3]-d1s-(d2ξ11[3]-d1s)2-4d1d2detJE32d1d2,k2=d2ξ11[3]-d1s+(d2ξ11[3]-d1s)2-4d1d2detJE32d1d2.Proof.
Using the same approach as in Theorem3, we have that the positive constant steady state solution E3 of model (86) is asymptotically stable provided (a) and (b).
For the Turing instability, we must havedet(Mi)<0 for some μi. And we notice that det(Mi) achieves its minimum:
(88)minμidet(Mi)=4d1d2det(JE3)-(d2ξ11[3]-d1s)24d1d2
at the critical value μ*>0 when μ*=((d2ξ11[3]-d1s)/2d1d2).
We now discuss the conditions for Turing instability in the further calculation. The conditionminμidet(Mi)<0 is equivalent to 4d1d2det(JEe2)-(d2ξ11[3]-d1s)2<0. In this case, det(Mi)=0 has two positive roots k1 and k2 which satisfy
(89)k1=d2ξ11[3]-d1s-(d2ξ11[3]-d1s)2-4d1d2detJE32d1d2,k2=d2ξ11[3]-d1s+(d2ξ11[3]-d1s)2-4d1d2detJE32d1d2.
Therefore, if we can find some μi such that k1<μi<k2, then det(Ai)<0. By [46] it follows that E3 is uniformly asymptotically instable. This finishes the proof.


## 5.2. Pattern Formation

In this section, we perform extensive numerical simulations of the spatially extended model (12) in two-dimensional space, and the qualitative results are shown here. All our numerical simulations employ the zero-flux boundary conditions with a system size of 100×100.The numerical integration of model (12) is performed by using a finite difference approximation for the spatial derivatives and an explicit Euler method for the time integration [54] with a time stepsize of 1/100. The initial condition is always a small amplitude random perturbation around the positive constant steady-state solution E3=(u*,u*+k2). After the initial period during which the perturbation spread, the model goes either into a time-dependent state, or to an essentially steady-state solution (time independent).We use the standard five-point approximation [55] for the 2D Laplacian with the zero-flux boundary conditions. More precisely, the concentrations (ui,jn+1,vi,jn+1) at the moment (n+1)τ at the mesh position (xi,yj) are given by
(90)ui,jn+1=ui,jn+τd1Δhui,jn+τf(ui,jn,vi,jn),vi,jn+1=vi,jn+τd2Δhvi,jn+τg(ui,jn,vi,jn),
with the Laplacian defined by
(91)Δhui,jn=ui+1,jn+ui-1,jn+ui,j+1n+ui,j-1n-4ui,jnh2,
where f(u,v)=u(1-u-(m/(u+b)))-(cuv/(u+k1)) and g(u,v)=sv(1-(v/(u+k2))), and the space stepsize h=1/3.In the numerical simulations, different types of dynamics are observed and it is found that the distributions of predator and prey are always of the same type. Consequently, we can restrict our analysis of pattern formation to one distribution. In this section, we show the distribution of preyu, for instance. We have taken some snapshots with red (blue) corresponding to the high (low) value of prey u.Now, we show the Turing pattern for the different values of the parameters. Via numerical simulation, one can see that the model dynamics exhibits spatiotemporal complexity of pattern formation, including spots, stripes, and holes Turing patterns.In Figure1, we show the time process of spots pattern formation of the prey u at t=0,300,1000,3000 for b=0.495, c=0.5, k1=0.3, k2=0.2, m=0.32, s=0.25, d1=0.1, and d2=2. In this case, there appears a competition between stripes and spots. The pattern takes a long time to settle down; starting with a homogeneous state E3=(0.066566,0.266566) (cf., Figure 1(a)), the random perturbations lead to the formation of stripes and spots (cf., Figure 1(b)), and the later random perturbations make these stripes decay and end with the time-independent regular spots pattern (cf., Figure 1(d)), which is isolated zones with low prey densities. Ecologically, spots pattern shows that the prey population is driven by predators to a very high level in those regions. The final result is the formation of patches of high prey density surrounded by areas of high prey densities; that is to say, under the control of these parameters, the prey is predominant in the area.The formation of spots patterns of preyu with parameters b=0.495,c=0.5,k1=0.3,k2=0.2,m=0.32,s=0.25,d1=0.1, and d2=2. Time: (a) t=0, (b) t=300, (c) t=1000, and (c) t=3000.
(a)
(b)
(c)
(d)In Figure2, we show the time process of stripes pattern formation of the prey u at t=0,300,1000,3000 for b=0.495,c=0.5,k1=0.3,k2=0.2,m=0.3,s=0.25,d1=0.1,d2=2. In this case, starting with a homogeneous state E3=(0.142318,0.342318) (cf., Figure 2(a)), the random perturbations lead to the formation of stripes and spots (cf., Figure 2(b)), and the later random perturbations make these stripes decay and end with the time-independent regular spots pattern (cf., Figure 2(d)).The formation of stripes patterns of preyu with parameters b=0.495,c=0.5,k1=0.3,k2=0.2,m=0.3,s=0.25,d1=0.1, and d2=2. Time: (a) t=0, (b) t=300, (c) t=1000, and (c) t=3000.
(a)
(b)
(c)
(d)In Figure3, we show the time process of holes pattern formation of the prey u at t=0,300,1000,3000 for b=0.4,c=0.5,k1=0.3,k2=0.2,m=0.24,s=0.25,d1=0.1, and d2=2. In this case, starting with a homogeneous state E3=(0.2,0.4) (cf., Figure 3(a)), the random perturbations lead to the formation of independent regular holes pattern (cf., Figure 1(d)), which is isolated zones with low prey densities.The formation of holes patterns of preyu with parameters b=0.4,c=0.5,k1=0.3,k2=0.2,m=0.24,s=0.25,d1=0.1, and d2=2. Time: (a) t=0, (b) t=300, (c) t=1000, and (c) t=3000.
(a)
(b)
(c)
(d)Comparing Figure1(d) with Figure 3(d), we find that they share similarities. Figure 1(d) consists of red (maximum density of u) spots on a blue (minimum density of u) background; that is, the preys are isolated zones with high population density. Figure 3(d) consists of blue (minimum density of u) spots on a red (maximum density of u) background; that is, the preys are isolated zones with low population density. For the sake of learning the pattern dynamics of model (12) further, we illustrate the three-dimensional patterns in the space (u,x,y). From Figure 4, one can realize the relations of the patterns (e.g., Figures 4(a) and 4(c)) with their corresponding numerical solutions (e.g., Figures 4(b) and 4(d)). In fact, the patterns (e.g., Figures 1–3) are the projections in xy plane of the numerical solutions u(x,y,t) to the model (12). In the software MATLAB, imagesc(u) displays u(x,y;t) as an image; each numerical solution of u(x,y;t) corresponds to a rectangular area in the image; that is, the values of the numerical solutions u(x,y;t) to the model (12) are indices into the current colormap that determine the color of each patch.Pattern dynamics of model (12) in two- and three-dimensional space. (a) spots pattern; (b) spots pattern in three-dimensional space; (c) holes pattern; (d) holes pattern in three-dimensional space.
(a)
(b)
(c)
(d)


## 6. Concluding Remarks

In the current investigation, we propose and analyze the dynamics of a reaction-diffusive Leslie-Gower predator-prey model with the additive Allee effect on prey. We are mainly concerned with the coexistence of the predator and prey and focus on the case of weak Allee effect (i.e.,m<b). The value of this study is twofold. First, it shows the nonexistence and existence of the nonconstant positive steady-states, which guarantees the existence of Turing patterns. Second, it rigorously proves Turing instability by linear stability analysis and illustrates all three categories of Turing patterns close to the onset of Turing bifurcation via numerical simulations which indicates that the model dynamics exhibits complex pattern replication.We summarize our findings as well as their related biological implications as follows.(1)
Theorems2 and 3 provide us with a full picture on the dynamics of the model with weak Allee effect. The dynamics of the model introduced can be very complicated due to Allee effects. In Theorem 2(a), we find that the trivial constant solution E0=(0,0) of the model which is subject to an Allee effect is unstable, which is exactly consistent with the model without Allee effect [24–26]. This demonstrates that there is no extinction of a species in the present of Allee effects. By comparing them to their corresponding models without Allee effect, we can conclude that Allee effect can make the extinction of the prey although the maximum predation rate c is small (see Theorem 2(b)). In this case, the predator species intends to change its food habits as predator has sufficient resources for alternative foods [56]. Furthermore, from Theorem 3, one can obtain that E3=(u*,u*+k2) is locally uniformly asymptotically stable, which means that nonconstant positive solution (stationary pattern) of model (12) unlikely exists.(2)
Theorems9 and 11 indicate the existence and nonexistence of nonconstant steady-states with respect to various parameters. Roughly speaking, we can state that there is no pattern if the diffusion coefficients are suitably chosen. While pattern occurs provided that θ11-4θ2>0 and all the eigenvalues μi have odd multiplicity (see Theorem 11). Hence, we can conclude that the multiplicity of patterns seems very interesting from the viewpoint of mathematics.(3)
Theorem15 and numerical simulations give us the existence of conditions of Turing instability and the types of Turing pattern. From the numerical results, one can see that our model has rich and complex spatiotemporal behavior. We find three typical Turing patterns, that is, stripes pattern, spot-stripe mixtures pattern, and spots pattern. That is to say, the effect of the Allee effect for pattern formation is tremendous. Therefore, the results of the present paper and [24, 25] show that the types of Turing pattern in the biological models depend on the effect of the Allee effect. In other words, the Allee effect may be one of the determining factors in producing spots and spot-stripe mixtures Turing patterns.It is believed that our results made in this investigation related to predator-prey interactions due to the effect of Allee effect would certainly be of some help to theoretical mathematicians and the ecologists who are engaged in performing experimental work. Further studies are necessary to analyze the behaviour of a reaction-diffusion predator-prey model with the strong Allee effect (i.e.,m>b).



---

*Source: 535746-2013-12-30.xml*