# Dynamics of a Diffusive Predator-Prey Model with Allee Effect on Predator

Authors: Xiaoqin Wang; Yongli Cai; Huihai Ma
Journal: Discrete Dynamics in Nature and Society
 (2013)

---

## Abstract

The reaction-diffusion Holling-Tanner prey-predator model considering the Allee effect on predator, under zero-flux boundary conditions, is discussed. Some properties of the solutions, such as dissipation and persistence, are obtained. Local and global stability of the positive equilibrium and Turing instability are studied. With the help of the numerical simulations, the rich Turing patterns, including holes, stripes, and spots patterns, are obtained.


---

## Body

## 1. Introduction

The Holling-Tanner prey-predator model is an important and interesting predator-prey model in both biological and mathematical sense [1–4]. The reaction-diffusion Holling-Tanner prey-predator model takes the following form:
(1)∂u∂t=ru(1-uK)-m1uvu+a1+D1Δu,x=(ξ,η)∈Ω,t>0,∂v∂t=s1v(1-hvu)+D2Δv,x=(ξ,η)∈Ω,t>0,
where u and v represent population density of prey and predator at time t, respectively. The parameters r,  K,  m1,  a1,  s1, and  h are all positive. r stands for the intrinsic growth rate of prey, K is the prey carrying capacity, m1 is the maximum predation rate, a1 is the self-saturation prey density, s1 is predator intrinsic growth rate, and h is conversion rate of prey into predator biomass. D1 and D2 are the diffusion coefficients of u and v, respectively, and we always assume that D1>0,  D2>0. Δ=∂2/∂x2=∂2/∂ξ2+∂2/∂η2 is the usual Laplacian operator in 2-dimensional space. Ω⊂ℝ2 is a bounded domain with smooth boundary ∂Ω.Set(2)(u,v,t)=(Ku~,Kv~,t~r).
For the sake of convenience, we still use variables u, v instead of u~, v~. Thus, model (1) is converted into
(3)∂u∂t=u(1-u)-muva+u+d1Δu,x=(ξ,η)∈Ω,t>0,∂v∂t=sv(1-hvu)+d2Δv,x=(ξ,η)∈Ω,t>0,
where the new parameters are
(4)m=Km1r,a=a1r,s=s1r,d1=D1r,d2=D2r.The dynamics of the reaction-diffusion Holling-Tanner prey-predator model has proven quite interesting and received intensive study by both ecologists and mathematicians in many articles, see, for example, [5–10] and the references therein. Peng and Wang [5, 6] analyzed the global stability of the unique positive constant steady state and established the results for the existence and nonexistence of positive nonconstant steady states; Wang et al. [7] studied positive steady-state solutions and investigated the appearance of sharp spatial patterns arising from the model. Shi and coworkers [8] 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. [9] considered the Turing and Hopf bifurcations of the equilibrium solutions; Liu and Xue [10] found the model exhibits the spotted, black-eye, and labyrinthine patterns.On the other hand, in population dynamics, Allee effect [11] is an ecological phenomenon caused by any mechanism leading to a positive relationship between a component of individual fitness and either the number or density of conspecific [12–15]. In ecological studies, the understanding of the influence of Allee effect plays a central role since the Allee effect can greatly increase the likelihood of local and global extinctions [16, 17]. As a result both ecologists and mathematicians are interested in Allee effect in the predator-prey model, and much progress has been seen in the study of Allee effect, see [14, 18–27], and many more investigations were done in recent years. But there have been few papers discussing the impact of the Allee effect on predator in the predator-prey models.Based on the previous discussion, in the present paper we adopt the reaction-diffusion Holling-Tanner prey-predator model with Allee effect on predator.If we assume that the predator population is subject to an Allee effect, taking into account zero-flux boundary conditions, model (3) can be rewritten as
(5)∂u∂t=u(1-u)-muva+u+d1Δu,x=(ξ,η)∈Ω,t>0,∂v∂t=sv(vv+b-hvu)+d2Δv,x=(ξ,η)∈Ω,t>0,∂u∂ν=∂v∂ν=0,x=(ξ,η)∈∂Ω,t>0,u(x,0)=u0(x)≥0,v(x,0)=v0(x)≥0,x=(ξ,η)∈Ω,
where v/(v+b) is the term for the Allee effect, and b can be defined as the Allee effect constant. The per capita growth rate of the predator is reduced from s to sv/(v+b) due to the Allee effect [17, 28]. ν is the outward unit normal vector on ∂Ω, and the zero-flux boundary conditions mean that model (5) is self-contained and has no population flux across the boundary ∂Ω [29, 30]. The initial data u0(x) and v0(x) are continuous functions on Ω¯.The corresponding kinetic equation to model (5) is
(6)u˙=u(1-u)-muva+u≜f(u,v),v˙=sv(vv+b-hvu)≜g(u,v).The plan of the paper is as follows. Section2 is dedicated to furnish some properties of the solutions concerning the mathematical model used. In Section 3, the local and global stability of the positive equilibrium of the model is considered. Section 4 is devoted to the diffusion-driven instability (Turing effect) and illustrates the different Turing patterns by using the numerical simulations. Finally, in Section 5, some conclusions and discussions are given.


## 2. Large Time Behavior of Solution to Model (5)

In this section, we give some properties of the solutions, and these results will be often used later.


### 2.1. Dissipation

Theorem 1.
For any solution(u(x,t),v(x,t)) of model (5) with nonnegative initial conditions, then
(7)limsupt→∞maxΩ-u(x,t)≤1,limsupt→∞maxΩ-v(x,t)≤1h.
Hence, for any ε>0, the rectangle [0,1+ε]×[0,1/h+ε] is a global attractor of model (5) in ℝ+.Proof.
u satisfies
(8)∂u∂t-d1Δu≤u(1-u),x∈Ω,t>0,∂u∂ν=0,x∈∂Ω,t>0,u(x,0)=u0(x),x∈Ω.
Letz(t) be a solution of the ordinary differential equation,
(9)z˙(t)=z(1-u),t≥0,z(0)=maxΩ-u(x,0)>0.
Then, limt→∞z(t)=1. From the comparison principle, one can get u(x,t)≤z(t); hence,
(10)limsupt→∞maxΩ-u(x,t)≤1.
As a result, for anyε>0, there exists a T>0, such that u(x,t)≤1+ε for all x∈Ω- and t≥T. Similarly, v satisfies
(11)∂v∂t-d2Δv≤sv(1-hv1+ε),x∈Ω,t>0,∂v∂ν=0,x∈∂Ω,t>0,v(x,0)=v0(x),x∈Ω.
Thus, limsupt→∞maxΩ-v(x,t)≤(1+ε)/h. From the arbitrariness of ε>0, we can get that
(12)limsupt→∞maxΩ-v(x,t)≤1h.
This ends the proof.


### 2.2. Persistence

Definition 2 (see [31]).
Model (5) is said to have the persistence property if for any nonnegative initial data (u0(x),v0(x)), there exists a positive constant ε=ε(u0,v0), such that the corresponding solution (u(x,t),v(x,t)) of model (5) satisfies
(13)liminft→∞minΩ-u(x,t)≥ε,liminft→∞minΩ-v(x,t)≥ε.
In the following, we will show that model (5) is persistent. From the viewpoint of biology, this implies that the two species of prey and predator will always coexist at any time and any location of the inhabit domain, no matter what their diffusion coefficients are, under certain conditions on parameters.Theorem 3.
Ifm<ah and 2bh<1-a+(1-a)2+4(a-mh-1), then model (5) has the persistence property.Proof.
The proof is based on comparison principles. From (12), for 0<ε≪1, it is clear that there exists a t≫1, such that v(x,t)<1/h+ε for all x∈Ω- and t≥t0. Hence, u(x,t) is an upper solution of the following problem:
(14)∂z∂t-d1Δz=z-z2+(1-a)z+a-m(h-1+ε)z+a,x∈Ω,t>t0,∂z∂ν=0,x∈∂Ω,t>t0,z(x,t0)=u0(x,t0)≥0,x∈Ω-.
Letω(t) be the unique positive solution to the following problem:
(15)dwdt=w-w2+(1-a)w+a-m(h-1+ε)w+a,t>t0,w(t0)=minΩ-u0(x,t0)≥0.
Sincem<ah, then limt→∞w(t)=(1-a+(1-a)2+4(a-mh-1))/2. By comparison, it follows that limt→∞z(x,t)=(1-a+(1-a)2+4(a-mh-1))/2. This implies that
(16)liminft→∞minΩ-u(x,t)≥1-a+(1-a)2+4(a-mh-1)2≜α.
Hence, u(x,t)>α-ε for t>t0 and x∈Ω-.
Similarly, by the second equation in model (5), we have that v(x,t) is an upper solution of problem
(17)∂z∂t-d2Δz=sz-hz2+(α-ε-bh)z(α-ε)(v+b),x∈Ω,t>t0,∂z∂ν=0,x∈∂Ω,t>t0,z(x,t0)=I0(x)≥0,x∈Ω-.
Letv(t) be the unique positive solution to the following problem:
(18)dwdt=sw-hw2+(α-ε-bh)w(α-ε)(w+b),t>t0,w(t0)=minΩ-v0(x,t0)≥0.
Since 2bh<1-a+(1-a)2+4(a-mh-1), there exists a ε>0, such that bh<α-ε. Hence, we have limt→∞w(t)=(1-a+(1-a)2+4(a-mh-1)-2bh)/2h for the arbitrariness of ε, and an application of the comparison principle gives
(19)liminft→∞minΩ-v(x,t)≥1-a+(1-a)2+4(a-mh-1)-2bh2h≜β.
The proof is complete.


## 2.1. Dissipation

Theorem 1.
For any solution(u(x,t),v(x,t)) of model (5) with nonnegative initial conditions, then
(7)limsupt→∞maxΩ-u(x,t)≤1,limsupt→∞maxΩ-v(x,t)≤1h.
Hence, for any ε>0, the rectangle [0,1+ε]×[0,1/h+ε] is a global attractor of model (5) in ℝ+.Proof.
u satisfies
(8)∂u∂t-d1Δu≤u(1-u),x∈Ω,t>0,∂u∂ν=0,x∈∂Ω,t>0,u(x,0)=u0(x),x∈Ω.
Letz(t) be a solution of the ordinary differential equation,
(9)z˙(t)=z(1-u),t≥0,z(0)=maxΩ-u(x,0)>0.
Then, limt→∞z(t)=1. From the comparison principle, one can get u(x,t)≤z(t); hence,
(10)limsupt→∞maxΩ-u(x,t)≤1.
As a result, for anyε>0, there exists a T>0, such that u(x,t)≤1+ε for all x∈Ω- and t≥T. Similarly, v satisfies
(11)∂v∂t-d2Δv≤sv(1-hv1+ε),x∈Ω,t>0,∂v∂ν=0,x∈∂Ω,t>0,v(x,0)=v0(x),x∈Ω.
Thus, limsupt→∞maxΩ-v(x,t)≤(1+ε)/h. From the arbitrariness of ε>0, we can get that
(12)limsupt→∞maxΩ-v(x,t)≤1h.
This ends the proof.


## 2.2. Persistence

Definition 2 (see [31]).
Model (5) is said to have the persistence property if for any nonnegative initial data (u0(x),v0(x)), there exists a positive constant ε=ε(u0,v0), such that the corresponding solution (u(x,t),v(x,t)) of model (5) satisfies
(13)liminft→∞minΩ-u(x,t)≥ε,liminft→∞minΩ-v(x,t)≥ε.
In the following, we will show that model (5) is persistent. From the viewpoint of biology, this implies that the two species of prey and predator will always coexist at any time and any location of the inhabit domain, no matter what their diffusion coefficients are, under certain conditions on parameters.Theorem 3.
Ifm<ah and 2bh<1-a+(1-a)2+4(a-mh-1), then model (5) has the persistence property.Proof.
The proof is based on comparison principles. From (12), for 0<ε≪1, it is clear that there exists a t≫1, such that v(x,t)<1/h+ε for all x∈Ω- and t≥t0. Hence, u(x,t) is an upper solution of the following problem:
(14)∂z∂t-d1Δz=z-z2+(1-a)z+a-m(h-1+ε)z+a,x∈Ω,t>t0,∂z∂ν=0,x∈∂Ω,t>t0,z(x,t0)=u0(x,t0)≥0,x∈Ω-.
Letω(t) be the unique positive solution to the following problem:
(15)dwdt=w-w2+(1-a)w+a-m(h-1+ε)w+a,t>t0,w(t0)=minΩ-u0(x,t0)≥0.
Sincem<ah, then limt→∞w(t)=(1-a+(1-a)2+4(a-mh-1))/2. By comparison, it follows that limt→∞z(x,t)=(1-a+(1-a)2+4(a-mh-1))/2. This implies that
(16)liminft→∞minΩ-u(x,t)≥1-a+(1-a)2+4(a-mh-1)2≜α.
Hence, u(x,t)>α-ε for t>t0 and x∈Ω-.
Similarly, by the second equation in model (5), we have that v(x,t) is an upper solution of problem
(17)∂z∂t-d2Δz=sz-hz2+(α-ε-bh)z(α-ε)(v+b),x∈Ω,t>t0,∂z∂ν=0,x∈∂Ω,t>t0,z(x,t0)=I0(x)≥0,x∈Ω-.
Letv(t) be the unique positive solution to the following problem:
(18)dwdt=sw-hw2+(α-ε-bh)w(α-ε)(w+b),t>t0,w(t0)=minΩ-v0(x,t0)≥0.
Since 2bh<1-a+(1-a)2+4(a-mh-1), there exists a ε>0, such that bh<α-ε. Hence, we have limt→∞w(t)=(1-a+(1-a)2+4(a-mh-1)-2bh)/2h for the arbitrariness of ε, and an application of the comparison principle gives
(19)liminft→∞minΩ-v(x,t)≥1-a+(1-a)2+4(a-mh-1)-2bh2h≜β.
The proof is complete.


## 3. Stability

In this section, we will devote consideration to the stability of the positive equilibrium for model (5).Clearly, model (5) has a unique positive equilibrium E*=(u*,v*), where u*=h(b+v*) and
(20)v*=h-m-ah-2bh22h2+(h-m-ah-2bh2)2+4h2(a+bh)(1-bh)2h2,
with bh<1.For the sake of simplicity, we rewrite model (5) as the vectorial form
(21)wt=DΔw+ℋ(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)ℋ(w)=(u(1-u)-muvu+asv(vv+b-hvu)).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)}  is an orthonormal basis of E(μi), and Xij={cϕij∣c∈ℝ2}, then
(24)X=⨁i=1∞Xi,
where Xi=⨁j=1dimE(μi)Xij.The linearization of model (5) at the positive equilibrium E*=(u*,v*) can be expressed by
(25)wt=£(w)=DΔw+Jw,
where
(26)J=(-u*+mu*v*(a+u*)2-mu*a+u*shv*2u*2-sv*2(v*+b)2)≜(J11J12J21J22).From [32], 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 part, then the stability of E*=(u*,v*) cannot be determined by the linearization.For eachi≥0, Xi is invariant under the operator £ and λ is an eigenvalue of £ if and only if λ is an eigenvalue of the matrix Ai=-μiD+J(u,v) for some i≥0.So, the local stability of the positive equilibriumE*=(u*,v*) can be analyzed as follows.Theorem 4.
Assume thats>u*(v*+b)2(mv*-(a+u*)2)/v*2(a+u*)2 and the first eigenvalue μ1 subject to the zero-flux boundary conditions satisfies
(27)μ1>max{0,u*(mv*-(a+u*)2)d1(a+u*)2-sv*2d2(v*+b)2}.
Then, the positive equilibrium E*=(u*,v*) is uniformly asymptotically stable.Proof.
The stability of the positive equilibriumE*=(u*,v*) is reduced to consider the characteristic equation
(28)det(λI-Ai)=λ2-tr(Ai)λ+det(Ai),
with
(29)tr(Ai)=-μi(d1+d2)+tr(J),det(Ai)=d1d2μi2-(J11d2+J22d1)μi+det(J).
In view ofs>u*(v*+b)2(mv*-(a+u*)2)/v*(a+u*)2, it follows that
(30)tr(J)=mu*v*(a+u*)2-u*-sv*2(v*+b)2<0.
Remark that for any i≥0, we have tr(Ai)<0.
In view of the relationu*=h(v*+b), one can calculate that
(31)det(J)=-sv*2(v*+b)2(-u*+mu*v*(a+u*)2)+mu*a+u*shv*2u*2=(sv*2(h3v*2+(2ah2+2bh3)v*+b2h3+2abh2+a2h+bhm+am))×((v*+b)(a+hv*+h)2)-1>0.
Recall that μ1>max{0,u*(mv*-(a+u*)2)/d1(a+u*)2-sv*2/d2(v*+b)2}, we conclude that
(32)det(Ai)=μi(d1d2μi-(J11d2+J22d1))+det(J)>μi(d1d2μ1-(J11d2+J22d1))+det(J)>0
for all i≥0.
Therefore, the eigenvalues of the matrix-μiD+J have negative real parts. It thus follows from the Routh-Hurwitz criterion that, for each i≥0, the two roots λi1 and λi2 of det(λI-Ai)=0 all have negative real parts.
In the following, we prove that there existsδ>0 such that
(33)Re{λi1}≤-δ,Re{λi2}≤-δ.
Letλ=μiξ, then
(34)φ~i(λ)≜μi2ξ2-tr(Ai)μiξ+det(Ai).
Since μi→∞ as i→∞, it follows that
(35)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~=min{d1,d2}, 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,  for all i≥i0. In turn, Re{λi1},Re{λi2}≤-μid~/2≤-d~/2, for all i≥i0. Let
(36)-δ~=max1≤i≤i0{Re{λi1},Re{λi2}}.
Then δ~>0 and (33) holds for δ=min{δ~,d~/2}.
Consequently, the spectrum of£, which consists of eigenvalues, lies in {Reλ≤-δ}. In the sense of [32], we obtain that the positive constant steady-state solution E*=(u*,u*+k2) of model (5) is uniformly asymptotically stable. This ends the proof.In the following, we shall prove that the positive equilibriumE*=(u*,v*) of model (5) is globally asymptotically stable.Theorem 5.
Suppose thatbh<1, m<ah, and 2bh<1-a+(1-a)2+4(a-mh-1). The positive equilibrium E*=(u*,v*) of model (5) is globally asymptotically stable, if(a1)
mh-1<(a+u*)(a+α),(a2)
m2+1/α2<(4u*/hv*)(h(a+u*)-b(a+u*)/(v*+b)(β+b)-m/α(a+α)),where
(37)α=1-a+(1-a)2+4(a-mh-1)2,β=1-a+(1-a)2+4(a-mph-1)-2bh2h.Proof.
We adopt the Lyapunov function(38)V(t)=∫Ω[V1(u(x,t))+V2(v(x,t))]dx,
where V1(u)=(u*+a)∫u*u‍((ξ-u*)/ξ)dξ, V2(v)=(u*/hsv*)∫v*v‍((η-v*)/η)dη. It can be easily verified that the function V(t) is zero at the positive equilibrium E*=(u*,v*) and is positive for all other positive values of u and v.
Then,(39)dVdt=∫Ω((u*+a)(u-u*)u∂u∂t+u*(v-v*)hsv*v∂v∂t)dx=∫Ω(∑(u*+a)(u-u*)(1-u-mva+u)+u*(v-v*)hv*(vv+b-hvu)∑)dx+∫Ω(d1(u*+a)(u-u*)uΔu+d2u*(v-v*)hsv*vΔv)dx=-∫Ω(a+u*-mva+u)(u-u*)2dx-∫Ωu*hv*(hu-b(v*+b)(v+b))(v-v*)2dx+∫Ω(1u-m)(u-u*)(v-v*)dx-∫Ω(d1u*(a+u*)u2|∇u|2+d2u*hsv2|∇v|2)dx=M-∫Ω(d1u*(a+u*)u2|∇u|2+d2u*hsv2|∇v|2)dx,
where
(40)M(u,v)=-∫Ω(a+u*-mva+u)(u-u*)2dx-∫Ωu*hv*(hu-b(v*+b)(v+b))(v-v*)2dx+∫Ω(1u-m)(u-u*)(v-v*)dx.
It is obvious thatdV/dt<0 if M(u,v) is negative definite.
M
(
u
,
v
) can be expressed in a quadratic form -XBXT, where
(41)X=(u-u*,v-v*),B=(a+u*-mva+um2-12um2-12uu*hv*(hu-b(v*+b)(v+b))).
M
(
u
,
v
) is negative definite if the symmetric matrices B is positive. It can be easily shown that the symmetric matrix B is positive definite if the following conditions are true:(i)
a+u*-mv/(a+u)>0,(ii)
Φ(u,v)≜(u*/hv*)(a+u*-mv/(a+u))(h/u-b/(v*+b)(v+b))-(1/4)(m-1/u)2>0.
Proof of (i). Applying Theorems 1 and 3, we get
(42)a+u*-mva+u>a+u*-mh-1a+α=(a+u*)(a+α)-mh-1a+α.
Therefore, if (a1) holds, then a+u*-mv/(a+u)>0 for all t≥0.
Proof of (ii).
Consider(43)Φ(u,v)=u*hv*(a+u*-mva+u)(hu-b(v*+b)(v+b))-14(m-1u)2>u*hv*(h(a+u*)u-b(a+u*)(v*+b)(v+b)-hmvu(a+u))-m24-14u2>u*hv*(h(a+u*)-b(a+u*)(v*+b)(β+b)-mα(a+α))-m24-14α2.
Consequently, if (a2) holds, Φ(u,v)>0.
Hence,E* is globally asymptotically stable for model (5) following the well-known theorem of Lyapunov stability.


## 4. Diffusion-Driven Instability: Turing Effect

In this section, we will investigate Turing instability and bifurcation for our model problem. We will also study pattern formation of the predator-prey solutions.


### 4.1. Turing Instability

Mathematically speaking, an equilibrium is Turing instability (diffusion-driven instability) means that it is an asymptotically stable equilibrium of model (6) but is unstable with respect to the solutions of diffusion model (5). In this subsection, we mainly focus on the emergency of the Turing instability of the positive equilibrium E*=(u*,v*).Now, the conditions for the positive equilibrium to be stable for the ODE are given by(44)det(J)=J11J22-J12J21>0,tr(J)=J11+J22<0.
Hence, Ai (the matrix Ai=-μiD+J) 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. A necessary condition for the Turing instability of model (5) is
(45)d2J11+d1J22>0,
Otherwise, det(Ai)>0 for all eigenvalues μi of the operator Δ since det(J)>0. For the Turing instability, we must have det(Ai)<0 for some μi. And we notice that det(Ai) achieves its minimum
(46)minμidet(Ai)=4d1d2det(J)-(d2J11+d1J22)24d1d2
at the critical value μ*>0 when
(47)μ*=d2J11+d1J222d1d2.However, the inequalityminμidet(Ai)<0 is necessary but not sufficient for the Turing instability in the bounded domain Ω. The possible eigenvalues μi are discrete. In this case, det(Ai)=0 has two positive roots k1 and k2(48)k1=d2J11+d1J22-(d2J11+d1J22)2-4d1d2det(J)2d1d2,k2=d2J11+d1J22+(d2J11+d1J22)2-4d1d2det(J)2d1d2,
so if we can find some μi such that k1<μi<k2, then det(Ai)<0, and the positive equilibrium E*=(u*,v*) of model (5) is unstable.Summarizing the previous analysis and calculations, we have the following results.Theorem 6.
Assume that the positive equilibriumE*=(u*,v*) exists. If the following conditions are true:(i)
J11+J22<0, that is,
(49)s>u*(v*+b)2(mv*-(a+u*)2)v*2(a+u*)2,(ii)
d2J11+d1J22>0, that is,
(50)s<d2u*(v*+b)2(mv*-(a+u*)2)d1v*2(a+u*)2,(iii)
d2J11+d1J22>2det(D)det(J), that is,
(51)d2u*(mv*-(a+u*)2)(a+u*)2-d1sv*2(v*+b)2>2d1d2det(J),then the positive equilibriumE* of model (5) is Turing unstable if 0<k1<μi<k2 for some μi.In Figure1, we show the Turing bifurcation diagram for model (5) with parameters m=0.1,  a=0.003,  h=0.066,  d1=0.015, and  d2=1 in b-s parameters plane. The Turing bifurcation breaks spatial symmetry, leading to the formation of patterns that are stationary in time and oscillatory in space [33]. Below the Turing bifurcation curve, the solution of the model is unstable, and Turing instability emerges, that is, Turing patterns emerge. This domain is called the “Turing space.” We will focus on the Turing pattern formation in this domain.Figure 1
Turing bifurcation diagram for model (5) using b and s as parameters. Other parameters are taken as: m=0.1,  a=0.003,  h=0.066,  d1=0.015,  d2=1.


### 4.2. Pattern Formation

In this section, we perform extensive numerical simulations of the spatially extended model (5) 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. Other parameters are set as a=0.03,  m=0.1,  h=0.066,  d1=0.015, and  d1=1.The numerical integration of model (5) is performed by using a finite difference approximation for the spatial derivatives and an explicit Euler method for the time integration [34, 35] with a time stepsize of 1/1000 and the space stepsize of h=1/10. The initial condition is always a small amplitude random perturbation around the positive equilibrium E*=(u*,v*). After the initial period during which the perturbation spreads, either the model goes into a time dependent state, or to an essentially steady-state solution (time independent).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.Figure2 shows the evolution process of the holes pattern of prey for the parameters (b,s)=(3.6,0.9) at 0,1×105,  2×105 and 3×105 iterations. In this case, one can see that for model (5), the random perturbations lead to the formation of stripes holes (cf. Figure 2(b)), and the later random perturbations make these stripes decay ending with the holes pattern (cf. Figure 2(d))—the prey are isolated zones with low population density.Holes pattern formation for model (5) by taking (b,s)=(3.6,0.9). Iterations: (a) 0; (b) 1×105; (c) 2×105; (d) 3×105.
(a)
(b)
(c)
(d)Figure3 shows the process of spatial pattern formation of prey for the parameters (b,s)=(0.8,2) at 0,0.5×105,  1.5×105, and 3×105 iterations. The random perturbations lead to the formation of stripe-holes patterns (cf. Figure 3(b)), and the later random perturbations make these holes decay, ending with a time-independent stripe pattern (cf. Figure 3(d))—the prey are interlaced stripes of high and low population densities.Stripes pattern formation for model (5) by taking (b,s)=(0.8,2). Iterations: (a) 0; (b) 0.5×105; (c) 1.5×105; (d) 3×105.
(a)
(b)
(c)
(d)Figure4 shows the process of spot pattern formation of prey for (b,s)=(0.8,10) at 0,2×105, 5×105, and 1×106 iterations. There is a competition exhibited between stripes and spots. The pattern takes a long time to settle down, starting with a homogeneous state E*(cf. Figure 4(a)), and the random perturbations lead to the formation of stripes and spots (cf. Figure 4(b)), ending with spots only (cf. Figure 4(d))—the prey are isolated zones with high population density.Spots pattern formation for model (5) by taking (b,s)=(0.8,10). Iterations: (a) 0; (b) 2×105; (c) 5×105; (d) 1×106.
(a)
(b)
(c)
(d)Ecologically speaking, spots pattern shows that the prey population is driven by predators to a high level in those regions, while holes pattern shows that the prey population is driven by predators to a very low level in those regions.


## 4.1. Turing Instability

Mathematically speaking, an equilibrium is Turing instability (diffusion-driven instability) means that it is an asymptotically stable equilibrium of model (6) but is unstable with respect to the solutions of diffusion model (5). In this subsection, we mainly focus on the emergency of the Turing instability of the positive equilibrium E*=(u*,v*).Now, the conditions for the positive equilibrium to be stable for the ODE are given by(44)det(J)=J11J22-J12J21>0,tr(J)=J11+J22<0.
Hence, Ai (the matrix Ai=-μiD+J) 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. A necessary condition for the Turing instability of model (5) is
(45)d2J11+d1J22>0,
Otherwise, det(Ai)>0 for all eigenvalues μi of the operator Δ since det(J)>0. For the Turing instability, we must have det(Ai)<0 for some μi. And we notice that det(Ai) achieves its minimum
(46)minμidet(Ai)=4d1d2det(J)-(d2J11+d1J22)24d1d2
at the critical value μ*>0 when
(47)μ*=d2J11+d1J222d1d2.However, the inequalityminμidet(Ai)<0 is necessary but not sufficient for the Turing instability in the bounded domain Ω. The possible eigenvalues μi are discrete. In this case, det(Ai)=0 has two positive roots k1 and k2(48)k1=d2J11+d1J22-(d2J11+d1J22)2-4d1d2det(J)2d1d2,k2=d2J11+d1J22+(d2J11+d1J22)2-4d1d2det(J)2d1d2,
so if we can find some μi such that k1<μi<k2, then det(Ai)<0, and the positive equilibrium E*=(u*,v*) of model (5) is unstable.Summarizing the previous analysis and calculations, we have the following results.Theorem 6.
Assume that the positive equilibriumE*=(u*,v*) exists. If the following conditions are true:(i)
J11+J22<0, that is,
(49)s>u*(v*+b)2(mv*-(a+u*)2)v*2(a+u*)2,(ii)
d2J11+d1J22>0, that is,
(50)s<d2u*(v*+b)2(mv*-(a+u*)2)d1v*2(a+u*)2,(iii)
d2J11+d1J22>2det(D)det(J), that is,
(51)d2u*(mv*-(a+u*)2)(a+u*)2-d1sv*2(v*+b)2>2d1d2det(J),then the positive equilibriumE* of model (5) is Turing unstable if 0<k1<μi<k2 for some μi.In Figure1, we show the Turing bifurcation diagram for model (5) with parameters m=0.1,  a=0.003,  h=0.066,  d1=0.015, and  d2=1 in b-s parameters plane. The Turing bifurcation breaks spatial symmetry, leading to the formation of patterns that are stationary in time and oscillatory in space [33]. Below the Turing bifurcation curve, the solution of the model is unstable, and Turing instability emerges, that is, Turing patterns emerge. This domain is called the “Turing space.” We will focus on the Turing pattern formation in this domain.Figure 1
Turing bifurcation diagram for model (5) using b and s as parameters. Other parameters are taken as: m=0.1,  a=0.003,  h=0.066,  d1=0.015,  d2=1.


## 4.2. Pattern Formation

In this section, we perform extensive numerical simulations of the spatially extended model (5) 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. Other parameters are set as a=0.03,  m=0.1,  h=0.066,  d1=0.015, and  d1=1.The numerical integration of model (5) is performed by using a finite difference approximation for the spatial derivatives and an explicit Euler method for the time integration [34, 35] with a time stepsize of 1/1000 and the space stepsize of h=1/10. The initial condition is always a small amplitude random perturbation around the positive equilibrium E*=(u*,v*). After the initial period during which the perturbation spreads, either the model goes into a time dependent state, or to an essentially steady-state solution (time independent).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.Figure2 shows the evolution process of the holes pattern of prey for the parameters (b,s)=(3.6,0.9) at 0,1×105,  2×105 and 3×105 iterations. In this case, one can see that for model (5), the random perturbations lead to the formation of stripes holes (cf. Figure 2(b)), and the later random perturbations make these stripes decay ending with the holes pattern (cf. Figure 2(d))—the prey are isolated zones with low population density.Holes pattern formation for model (5) by taking (b,s)=(3.6,0.9). Iterations: (a) 0; (b) 1×105; (c) 2×105; (d) 3×105.
(a)
(b)
(c)
(d)Figure3 shows the process of spatial pattern formation of prey for the parameters (b,s)=(0.8,2) at 0,0.5×105,  1.5×105, and 3×105 iterations. The random perturbations lead to the formation of stripe-holes patterns (cf. Figure 3(b)), and the later random perturbations make these holes decay, ending with a time-independent stripe pattern (cf. Figure 3(d))—the prey are interlaced stripes of high and low population densities.Stripes pattern formation for model (5) by taking (b,s)=(0.8,2). Iterations: (a) 0; (b) 0.5×105; (c) 1.5×105; (d) 3×105.
(a)
(b)
(c)
(d)Figure4 shows the process of spot pattern formation of prey for (b,s)=(0.8,10) at 0,2×105, 5×105, and 1×106 iterations. There is a competition exhibited between stripes and spots. The pattern takes a long time to settle down, starting with a homogeneous state E*(cf. Figure 4(a)), and the random perturbations lead to the formation of stripes and spots (cf. Figure 4(b)), ending with spots only (cf. Figure 4(d))—the prey are isolated zones with high population density.Spots pattern formation for model (5) by taking (b,s)=(0.8,10). Iterations: (a) 0; (b) 2×105; (c) 5×105; (d) 1×106.
(a)
(b)
(c)
(d)Ecologically speaking, spots pattern shows that the prey population is driven by predators to a high level in those regions, while holes pattern shows that the prey population is driven by predators to a very low level in those regions.


## 5.  Conclusions and Remarks

In this paper, we have studied the dynamics of a reaction-diffusion Holling-Tanner prey-predator model where the predator population is subject to Allee effect under the zero-flux boundary conditions. The value of this study lies in threefolds. First, it investigates qualitative properties of solutions to this reaction-diffusion model. Second, it gives local and global stability of the positive equilibrium of the model. Third, it rigorously proves the Turing instability and illustrates three categories of Turing patterns close to the onset Turing bifurcation, which shows that the model dynamics exhibits complex pattern replication.It is seen that if Allee effect constantb is low, then the persistence of the model is guaranteed. It is interesting to notice that, from the result of Theorem 5, the condition for global stability of E*=(u*,v*) is independent of the diffusion coefficients di (i=1,2). So, it can be said that when the conditions on the parameters are satisfied, E*=(u*,v*) is stabilized under arbitrary spatially inhomogeneous perturbation. A very interesting observation can be made from the result of the numerical simulations. It indicates that the spatial model dynamics exhibits a diffusion-controlled formation growth not only to holes (cf. Figure 1) and stripes (cf. Figure 2) but also to spots replication (cf. Figure 3).Comparing Figures1 and 3, we can conclude that the model exhibits holes and spots Turing patterns due to Allee effect constant b. It is believed that the observations made in this investigation related to Allee effect on predator population remind us of the importance of the Allee effect. Therefore, if the prey or the predator to be protected is subject to an Allee effect, the measures taken should take this into account.



---

*Source: 984960-2013-02-17.xml*