# Complex Dynamics of a Diffusive Holling-Tanner Predator-Prey Model with the Allee Effect

Authors: Zongmin Yue; Xiaoqin Wang; Haifeng Liu
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We investigate the complex dynamics of a diffusive Holling-Tanner predation model with the Allee effect on prey analytically and numerically. We examine the existence of the positive equilibria and the related dynamical behaviors of the model and find that when the model is with weak Allee effect, the solutions are local and global stability for some conditions around the positive equilibrium. In contrast, when the model is with strong Allee effect, this may lead to the phenomenon of bistability; that is to say, there is a separatrix curve that separates the behavior of trajectories of the system, implying that the model is highly sensitive to the initial conditions. Furthermore, we give the conditions of Turing instability and determine the Turing space in the parameters space. Based on these results, we perform a series of numerical simulations and find that the model exhibits complex pattern replication: spots, spots-stripes mixtures, and stripes patterns. The results show that the impact of the Allee effect essentially increases the models spatiotemporal complexity.


---

## Body

## 1. Introduction

Recently, there has been a great interest in studying nonlinear difference/differential equations and systems [1–6]. One of the reasons for this is a necessity for some techniques which can be used in investigating equations arising in mathematical models describing real-life situations in population biology, economy, probability theory, genetics, psychology, sociology, and so forth. And the bases for analyzing the dynamics of complex ecological systems are the interactions between two species, particularly the dynamical relationship between predators and their preys [7]. From the Lotka-Volterra model [8, 9], several alternatives for modeling continuous time consumer-resource interactions have been proposed. In recent years, one of the important predator-prey models is Holling-Tanner model, which was described by May [10]. This model reads as follows:
(1)dHdt=r1H(1-HK)-c1HPk1+H,dPdt=s1P(1-PδH),
where H and P represent prey and predator population densities at time t, respectively. r, K, c1, k1, s1, and δ are positive constants. r1 and s1 are the intrinsic growth rate of prey and predator, respectively. K is the carrying capacity of the prey, and δ takes on the role of the prey-dependent carrying capacity for the predator. The rate at which the predator consumes the prey, c1HP/(k1+H), is known as the Holling type-II functional response [11]. The parameter c1 is the maximum number of the prey that can be eaten per predator per time, and k1 is the saturation value that corresponds to the number of the preys necessary to achieve one half of the maximum rate c1.The dynamics of model (1) has been considered in many articles. For example, Hsu and Huang [12] obtained some results on the global stability of the positive equilibrium. More precise, under the conditions which local stability of the positive equilibrium implies its global stability. Gasull and coworkers [13] investigated the conditions of the asymptotic stability of the positive equilibrium which does not imply global stability. Sáez and González-Olivares [14] showed the asymptotic stability of a positive equilibrium and gave a qualitative description of the bifurcation curve.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 conspecifics constitutes what is usually called an Allee effect [15–24], starting with the pioneer work of Allee [25]. The outflux of prey to constant rate can be considered as Allee effect because a change on interaction dynamics is provoked, for instance, due to difficulty of encountering mates [17]. Nowadays, it is widely accepted that the Allee effect greatly increases the likelihood of local and global extinction [18] and can lead to a rich variety of dynamical effects.From an ecological point of view, the Allee effect has been denominated in different ways [19–22] and modeled into strong and weak ones [15, 16, 19], depending on the degree of positive density dependence. Mathematically speaking, if H=H(t) indicates the population size, we assume that the growth function G(H) satisfies the following:(i)
ifG(0)=0, G′(0)>0, G(H) is called weak Allee effect;(ii)
ifG(0)=0, G′(0)<0, G(H) is called strong Allee effect.The most common mathematical form describing this phenomenon for a single species is given by(2)G(H)=H(1-H)(H-m),
where 0<m<1 or -1<m≤0, which is named the multiplicative Allee effect; here, a threshold value m is incorporated such that population growth is negative below m. When m<H<1, the per capita growth rate is positive.Furthermore, Boukal et al. [22] proposed that the prey exhibits a demographic Allee effect at low population densities due to reasons other than predation by the focal predator as follows:
(3)G(H)=AH(1-H)(1-b+cH+c),
where b is the Allee threshold, and c is an auxiliary parameter (c>0 and b≥-c). The auxiliary parameter c affects the overall shape of the per capita growth curve of the prey. When c is fixed, the unit growth rate of the species is only in connection with the Allee threshold.For model (1), we make a change of variables as follows:
(4)(H,P,t)=(KH~,KP~,t~r1).
For the sake of convenience, we still use variables H and P instead of H~ and P~.(H1) The basic model is a Holling-Tanner type as the form
(5)dHdt=H(1-H)-aHP1+H,dPdt=rP(1-PδH),
where a=c1/K, k1=K, and r=s1/r.(H2)  Following Boukal et al. [22], in Allee effect equation (3), we choose the auxiliary parameter c=1, and b+1≡m is the Allee threshold. That is, prey H has the population growth function
(6)G(H)=H(1-H)(1-mH+1).Obviously, we have the following:(i)
if0<m≤1, G(0)=0, G′(0)>0, the Allee effect (6) is the weak one;(ii)
ifm>1, G(0)=0, G′(0)<0, the Allee effect (6) is the strong one;(iii)
ifm=0, the Allee effect will disappear.And we can get the following model with the Allee effect on prey:(7)dHdt=H(1-H)(1-mH+1)-aHP1+H,dPdt=rP(1-PδH).(H3) Assume that the individuals in populations H and P move randomly described as Brownian random motion [26]. We can get a simple spatial model corresponding to model (7) as follows:
(8)∂H∂t=H(1-H)(1-mH+1)-aHP1+H+d1∇2H,∂P∂t=rP(1-PδH)+d2∇2P,H(x,y,0)=H0>0,P(x,y,0)=P0>0,(x,y)∈Ω=(0,L)×(0,L).
Here, the nonnegative constants d1 and d2 are the diffusion coefficients of H(t) and P(t), respectively. ∇2=∂2/∂x2+∂2/∂y2 is the Laplacian operator in two-dimensional space, which describes the random moving. The initial distribution of species N0 and P0 are continuous functions. And the boundary condition is assumed to be zero-flux one as follows:
(9)∂H∂n=∂P∂n=0,(x,y)∈∂Ω.L indicates the size of the model in the directions of x and y, respectively, and n is the outward unit normal vector of the boundary ∂Ω. The main reason for choosing such boundary conditions is that we are interested in the self-organization of pattern, and the zero-flux boundary conditions imply that no external input is imposed form exterior [27].There are some excellent works on a Holling-Tanner model considering the diffusion [28–33] and the references therein. In [28], Guan and co-workers studied the spatiotemporal dynamics of a modified version of the Leslie-Gower predator-prey model incorporating a prey refuge and showed that the model dynamics exhibits complex Turing pattern replication: stripes, cold/hot spots-stripes coexistence, and cold/hot spots patterns. Without the Allee effect, Peng and Wang [29, 30] analyzed the global stability of the unique positive constant steady state and established some results for the existence and nonexistence of positive nonconstant steady states. Wang et al. [31] considered the Turing and Hopf bifurcations of the equilibrium solutions. Liu and Xue [32] investigated the pattern formation and found that spots, black-eye, and labyrinthine patterns can be observed in the model. Chen and Shi [33] proved global stability of the unique constant equilibrium.However, the research about the influence of Allee effect on pattern formation of diffusive Holling-Tanner model seems rare. The main purpose of this paper is to study dynamical behaviors of a Holling-Tanner predator-prey model with the Allee effect. We will determine how the Allee effect affects the dynamics of the model and focus on the stability of the positive steady state and bifurcation mechanism and patterns formation analysis of the model.The rest of the paper is organized as follows. In Sections2 and 3, we present our main results about the stability and bifurcation analysis of the nonspatial model (7) and the spatial model (8), respectively. Especially, in regards to the spatial model (8) in Section 3, we will give the conditions of the Turing instability and determine the Turing space, and by performing a series of numerical simulations, we illustrate the emergence of different patterns. Finally, in Section 4, some conclusions and remarks are given.


## 2. Dynamics Analysis of the Nonspatial Model (7)

### 2.1. Boundedness

Now, we prove that all solutions are eventually bounded.Theorem 1.
All the solutions of model (7) which are initiated in ℝ+2 are uniformly bounded.Proof.
LetH(t) and P(t) be any solution of model (7) with initial conditions (H(0),P(0))=(H0,P0) such that H0>0, P0>0. From the first equation of model (7), we have
(10)dHdt≤H(1-H);
a standard comparison theorem shows that
(11)limsupt→∞H(t)≤1.
Then, from the second equation of model (7), we get dP/dt≤rP(1-(P/δ)), which implies that
(12)limsupt→∞P(t)≤δ.
Define the functionW(t)=H(t)+P(t), differentiating both sides with respect to t; we get
(13)dWdt=dHdt+dPdt≤H(1-H)+rP(1-Pδ)≤14+rP(1-Pδ).
Then,
(14)dWdt+W≤14+rP(1-Pδ)+H+P≤54+P(1+r-rPδ)≤54+δ(r+1)24r≜M.
Using the theory of differential inequality, for allt≥T≥0, we have
(15)0≤W(t)≤M-(M-W(T))e-(t-T).
Hence, limsupt→∞(H(t)+P(t))≤M. This completes the proof.Remark 2.
In fact, ifm≥2, dH/dt<0 always holds, which means that the prey and predator will extinct. Hence, we will later only focus on the case of 0≤m<2.
Next, we will investigate the existence of equilibria and their local and global stability with respect to model (7).


### 2.2. Equilibria Analysis in the Case of the Strong Allee Effect (i.e.,1<m<2)

In this subsection, we consider the existence and stability of the equilibrium of model (7) with strong Allee effect; that is, 1<m<2.We note that model (7) is not defined at the P-axis, particularly at the point (0,0), but both isoclines pass through this point, and in this case, it is a point of particular interest [34]. The character of (0,0) can be obtained after rescaling the time in model (7) by t=τH(1+H) as follows:
(16)dHdτ=H2(1-H)(1+H-m)-aH2P,dPdτ=r(1+H)P(H-Pδ).Lemma 3.
The point(0,0) of model (16) has a hyperbolic and a parabolic sector [20, 35] determined for the line P=(δ(m-1+r)/r)H. That is, there exists a separatrix curve in the phase plane that divides the behavior of trajectories; the point (0,0) is then an attractor point for certain trajectories and a saddle point for others.Proof.
As the Jacobian matrix of the point(0,0) for model (16) is the zero matrix, we follow the methodology used in [20, 35] given by the function φ(u,v)=(uv,v)=(H,P). Then, we have that
(17)dudτ=1v(dHdτ-udPdτ),dvdτ=dPdτ,
and rescaling the time by T=vτ, it becomes
(18)dudT=u((u-1δ)((1-uv)(uv+1-m)-av)u-r(1+uv)(u-1δ)),dvdT=r(1+uv)v(u-1δ).
Clearly, ifv=0, then dv/dT=0. Moreover, du/dT=u((1-m)u-r(u-(1/δ))).
The singularities of model (18) are (0,0) and (r/(δ(m-1+r)),0); that is, a separatrix straight exists in the phase plane uv, given by u=r/(δ(m-1+r)). The Jacobian matrixes of (0,0) and (r/(δ(m-1+r)),0) for model (18) are
(19)J(0,0)=(rδ00-rδ),J(r/(δ(m-1+r)),0)=(-rδ-r2(-2rm+aδ(m-1+r)+r)δ3(m-1+r)30-r(m-1)δ(m-1+r)).
Then,(0,0) is a hyperbolic saddle point, and (r/(δ(m-1+r)),0) is an attractor point. Using the blowing down, the point (0,0) is a saddle node in model (16), and the line P=((δ(m-1+r))/r)H divides the behavior of trajectories on the phase plane. The proof is completed.Moreover, it is easy to verify that model (7) always has two boundary equilibria E0=(m-1,0) and E1=(1,0). And the behavior of model (7) around E0 and E1 is found as follows.The Jacobian matrix of model (7) at the equilibrium E0=(m-1,0) takes the form
(20)JE1=((m-1)(2-m)ma(m-1)m0r).
Hence, the equilibrium E0=(m-1,0) is an unstable node point (nodal source).The Jacobian matrix of model (7) at the equilibrium E0=(1,0) takes the form
(21)JE0=(m2-1-a20r).
Hence, the equilibrium E1=(1,0) is a saddle point.And model (7) has a positive equilibrium E=(H,δH), where H satisfies
(22)H2-(m-aδ)H-(1-m)=0.
For simplicity, we consider A=m-aδ and B=(m-aδ)2+4(1-m); then, the two roots of (22) are given by
(23)H+=12(A+B),H-=12(A-B).Lemma 4.
(i)
Suppose thatm-aδ>0 and 1<m<2.
(a)
IfB2>0 holds, model (7) has two positive equilibria E+=(H+,δH+) and E-=(H-,δH-).
(b)
IfB2=0 holds, model (7) has a unique positive equilibrium Ee=(He,δHe). Note that in this case He=H+=H-=A/2=m-1.
(c)
IfB2<0, model (7) has no positive equilibrium.
(ii)
Ifm-aδ≤0, model (7) has no positive equilibrium.LetE=(N,P) be an arbitrary positive equilibrium. The Jacobian matrix of model (7) at the positive equilibrium E=(H,δH) takes the form
(24)JE=(H(aδH+2m-(1+H)2)(1+H)2-aH1+Hrδ-r).
Then, we can get
(25)det(JE)=rH(H2+2H+1+aδ-2m)(H+1)2,tr(JE)=H(aδH+2m-(1+H)2)(1+H)2-r.We can see that the sign ofdet(JE) is determined by
(26)F(H)≜H2+2H+1+aδ-2m=H2+2H+1-A-m.Thus, we can obtain(27)F(H+)=(A+B)24+B-m+1=B(1+A2)+12B2>0;F(H-)=(A-B)24-B-m+1=-12B(A-B)-B<0;F(He)=A4+1-m=0.Hence, we obtaindet(JE+)>0, det(JE-)<0, and det(JEe)=0. And the positive equilibrium E-=(H-,δH-) is a saddle point. The nature of the equilibrium point E+ is dependent on the sign of the trace of the Jacobian matrix evaluated in this point. Whether E+ is a node or a focus depends on the sign of (tr(JE+))2-4det(JE+).In the following results, we study the stability of the positive equilibriumE+ and the unique positive equilibria Ee.Theorem 5.
Define(28)r+=H+(1+H+)2(aδH++2m-(1+H+)2).(a)
Ifr>r+, the positive equilibrium E+=(H+,δH+) is a locally asymptotically stable point;
(a1)
if(r+-r)2<4det(JE+), then E+ is a stable focus,(a2)
if(r+-r)2>4det(JE+), then E+ is a stable node point.(b)
Ifr<r+, the positive equilibrium E+=(H+,δH+) is an unstable point;
(b1)
if(r+-r)2<4det(JE+), then E+ is an unstable focus surrounded by a stable limit cycle,(b2)
if(r+-r)2>4det(JE+), then E+ is an unstable node and the limit cycle disappears.(c)
A Hopf bifurcation occurs atr=r+ around the positive equilibrium E+=(H+,δH+). That is to say, model (7) has at least one positive periodic orbit.Proof.
Here, we only give the proof of the existence of Hopf bifurcation. It is easy to see that(i)
tr(J(E+))|r=r+=0 holds,(ii)
the characteristic equation isλ2+det(J(E+))|r=r+=0, whose roots are purely imaginary,(iii)
(d/dr)[tr(J(E+))]r=r+=-1≠0.
From the Poincaré-Andronov-Hopf Bifurcation Theorem [36], we know that model (7) undergoes a Hopf bifurcation at E+ as r passes through the value r+. The proof is completed.Figure1 illustrates the local stability of the positive equilibrium E+ and the separatrix curve generated by the stable manifold of the positive equilibrium E-. The orbits initiating the right of the separatrix curve tend to E+, while the orbits initiating the left of the separatrix curve tend to (0,0) that represents the extinction of the population. Figure 2 illustrates a Hopf bifurcation situation of the model around E+. The parameter values are given in the figures.Figure 1
The phase portrait of model (7) with the strong Allee effect. The parameters are taken as a=0.25, δ=0.8, m=1.3 and r=0.15. In this case, E0=(0.3,0) is an unstable node point, E1=(1,0) and E-=(0.5,0.4) are saddle points; the positive equilibrium E+=(0.6,0.48) is local asymptotically stable. There exists a separatrix curve determined by the stable manifold of the equilibrium point E-. The dotted curves are the nullclines.Figure 2
The phase portrait of model (7) with the strong Allee effect. The parameters are taken as a=0.25, δ=0.8, m=1.3, and r=0.0375. The model enters into a Hopf bifurcation around E+=(0.6,0.48) at r=r+.Theorem 6.
The unique equilibrium pointEe=(m-1,(1/δ)m-1) is(i)
a nonhyperbolic attractor node, if and only ifr>m-1(1-m-1)2/(1+m-1);(ii)
a nonhyperbolic repellor node, if and only ifr>m-1(1-m-1)2/(1+m-1);(iii)
a cusp point, if and only ifr=m-1(1-m-1)2/(1+m-1), and in this case, there exists a unique trajectory which attains the point Ee. And in this case, the point (0,0) is a global attractor.Proof.
We have(29)tr(JEe)=m-1(1-m-1)21+m-1-r.
Hence (i) and (ii) hold.
Moreover,tr(JEe)=0, if and only if r=m-1(1-m-1)2/(1-m+1). In this case, we obtain the Jacobian matrix of (16) as follows:
(30)J~(m-1,1δm-1)=((m-1)(1-m-1)2-(m-1)(1-m-1)2(m-1)(1-m-1)2-(m-1)(1-m-1)2)=(m-1)(1-m-1)2(1-11-1);
and the associate Jordan matrix is
(31)(0-(m-1)(1-m-1)200).
Then, the singularity (m-1,(1/δ)m-1) is a cusp point, since it is a point of codimension 2, and we have a Bogdanov-Takens bifurcation [37].The cusp point is shown in Figure3.Figure 3
The phase portrait of model (7) with the strong Allee effect. The parameters are taken as a=0.2556936062, δ=0.8, m=1.3, and r=0.07238979895. E0=(0.3,0) is an unstable node point; E1=(1,0) is saddle point; the positive equilibrium E2=(0.2,0.4) is a cusp point. In this case, the point (0,0) is globally asymptotically stable. The dotted curves are the nullclines.


### 2.3. Equilibria Analysis in the Case of the Weak Allee Effect (i.e.,0<m≤1)

In this subsection, we consider the stability of the equilibrium of model (7) with weak Allee effect (0<m≤1).It is easy to verify that model (7) always has one boundary equilibrium E0=(1,0) which is a saddle point and a positive equilibrium E*=(H*,δH*), where
(32)H*=m-aδ+(m-aδ)2+4(1-m)2.From (24), we have
(33)det(JE*)=rH*(H*2+2H*+1+aδ-2m)(H*+1)2>0,tr(JE*)=H*(aδH*+2m-(1+H*)2)(1+H*)2-r.Hence, we have the following results on the stability of the positive equilibriumE*=(H*,δH*).Theorem 7.
Define(34)r*=H*(aδH*+2m-(1+H*)2)(1+H*)2.(a)
Ifr>r*, the positive equilibrium E*=(H*,δH*) is a locally asymptotically stable point, and
(a1)
if(r*-r)2<4det(JE*), then E* is a stable focus;(a2)
if(r*-r)2>4det(JE*), then E* is a stable node point.(b)
Ifr<r*, the positive equilibrium E*=(H*,δH*) is an unstable point, and
(b1)
if(r*-r)2<4det(JE*), then E* is an unstable focus surrounded by a stable limit cycle;(b2)
if(r*-r)2>4det(JE*), then E* is an unstable node and the limit cycle disappears.(c)
A Hopf bifurcation occurs atr=r* around the positive equilibrium E*=(H*,δH*). That is to say, model (7) has at least one positive periodic orbit.In the following theorem, we study the global behavior of the positive equilibriumE*.Theorem 8.
If0<m<1/(1+aδ), the positive equilibrium E*=(H*,δH*) is globally asymptotically stable.Proof.
Construct the following Lyapunov function:(35)V(H,P)=∫H*H‍ξ-H+ξϕ(ξ)dξ+1r∫δH*P‍η-δH*ηdη,
where
(36)ϕ(H)=aHH+1.
Then,
(37)dVdt=H-H*Hϕ(H)dHdt+P-δH*rPdPdt.
Substituting the value ofdH/dt and dP/dt from the model of (7), we obtained
(38)dVdt=H-H*aH[(1-H)(1+H-m)-aδH*]-δH(P-δH*)2.
Note thataδH*=H*(1-H*)(1+H*-m); we obtain
(39)dVdt=-(H-H*)2aH(H+H*-m)-δH(P-δH*)2.
Hence, if 0<m<1/(1+aδ), H*-m>0 which is equivalent to dV/dt<0.
Hence, the positive equilibriumE*=(H*,δH*) is globally asymptotically stable. This completes the proof.Figure4 demonstrates the global stability situation of model (7) around E*. Figure 5 illustrates a Hopf bifurcation situation of the model around E*. Figure 6 shows a stable limit cycle around E* which is an unstable focus. The parameter values are given in the figures.Figure 4
The phase portrait of model (7) with the weak Allee effect. The parameters are taken as a=0.25, δ=0.8, m=1.3, and r=0.1. In this case, E0=(1,0) is saddle point; the positive equilibrium E*=(0.3232928050,0.4849392075) is globally asymptotically stable. The dotted curves are the nullclines.Figure 5
The phase portrait of model (7) with the weak Allee effect. The parameters are taken as a=0.25, δ=0.8, m=1.3, and r=0.025265. The model enters into a Hopf bifurcation around E*=(0.3232928050,0.4849392075) at r=r*.Figure 6
The phase portrait of model (7) with the weak Allee effect. The parameters are taken as a=0.8, δ=1.5, m=0.75, and r=0.025. The positive equilibrium E*=(0.3232928050,0.4849392075) is an unstable focus surrounded by a stable limit cycle.


## 2.1. Boundedness

Now, we prove that all solutions are eventually bounded.Theorem 1.
All the solutions of model (7) which are initiated in ℝ+2 are uniformly bounded.Proof.
LetH(t) and P(t) be any solution of model (7) with initial conditions (H(0),P(0))=(H0,P0) such that H0>0, P0>0. From the first equation of model (7), we have
(10)dHdt≤H(1-H);
a standard comparison theorem shows that
(11)limsupt→∞H(t)≤1.
Then, from the second equation of model (7), we get dP/dt≤rP(1-(P/δ)), which implies that
(12)limsupt→∞P(t)≤δ.
Define the functionW(t)=H(t)+P(t), differentiating both sides with respect to t; we get
(13)dWdt=dHdt+dPdt≤H(1-H)+rP(1-Pδ)≤14+rP(1-Pδ).
Then,
(14)dWdt+W≤14+rP(1-Pδ)+H+P≤54+P(1+r-rPδ)≤54+δ(r+1)24r≜M.
Using the theory of differential inequality, for allt≥T≥0, we have
(15)0≤W(t)≤M-(M-W(T))e-(t-T).
Hence, limsupt→∞(H(t)+P(t))≤M. This completes the proof.Remark 2.
In fact, ifm≥2, dH/dt<0 always holds, which means that the prey and predator will extinct. Hence, we will later only focus on the case of 0≤m<2.
Next, we will investigate the existence of equilibria and their local and global stability with respect to model (7).


## 2.2. Equilibria Analysis in the Case of the Strong Allee Effect (i.e.,1<m<2)

In this subsection, we consider the existence and stability of the equilibrium of model (7) with strong Allee effect; that is, 1<m<2.We note that model (7) is not defined at the P-axis, particularly at the point (0,0), but both isoclines pass through this point, and in this case, it is a point of particular interest [34]. The character of (0,0) can be obtained after rescaling the time in model (7) by t=τH(1+H) as follows:
(16)dHdτ=H2(1-H)(1+H-m)-aH2P,dPdτ=r(1+H)P(H-Pδ).Lemma 3.
The point(0,0) of model (16) has a hyperbolic and a parabolic sector [20, 35] determined for the line P=(δ(m-1+r)/r)H. That is, there exists a separatrix curve in the phase plane that divides the behavior of trajectories; the point (0,0) is then an attractor point for certain trajectories and a saddle point for others.Proof.
As the Jacobian matrix of the point(0,0) for model (16) is the zero matrix, we follow the methodology used in [20, 35] given by the function φ(u,v)=(uv,v)=(H,P). Then, we have that
(17)dudτ=1v(dHdτ-udPdτ),dvdτ=dPdτ,
and rescaling the time by T=vτ, it becomes
(18)dudT=u((u-1δ)((1-uv)(uv+1-m)-av)u-r(1+uv)(u-1δ)),dvdT=r(1+uv)v(u-1δ).
Clearly, ifv=0, then dv/dT=0. Moreover, du/dT=u((1-m)u-r(u-(1/δ))).
The singularities of model (18) are (0,0) and (r/(δ(m-1+r)),0); that is, a separatrix straight exists in the phase plane uv, given by u=r/(δ(m-1+r)). The Jacobian matrixes of (0,0) and (r/(δ(m-1+r)),0) for model (18) are
(19)J(0,0)=(rδ00-rδ),J(r/(δ(m-1+r)),0)=(-rδ-r2(-2rm+aδ(m-1+r)+r)δ3(m-1+r)30-r(m-1)δ(m-1+r)).
Then,(0,0) is a hyperbolic saddle point, and (r/(δ(m-1+r)),0) is an attractor point. Using the blowing down, the point (0,0) is a saddle node in model (16), and the line P=((δ(m-1+r))/r)H divides the behavior of trajectories on the phase plane. The proof is completed.Moreover, it is easy to verify that model (7) always has two boundary equilibria E0=(m-1,0) and E1=(1,0). And the behavior of model (7) around E0 and E1 is found as follows.The Jacobian matrix of model (7) at the equilibrium E0=(m-1,0) takes the form
(20)JE1=((m-1)(2-m)ma(m-1)m0r).
Hence, the equilibrium E0=(m-1,0) is an unstable node point (nodal source).The Jacobian matrix of model (7) at the equilibrium E0=(1,0) takes the form
(21)JE0=(m2-1-a20r).
Hence, the equilibrium E1=(1,0) is a saddle point.And model (7) has a positive equilibrium E=(H,δH), where H satisfies
(22)H2-(m-aδ)H-(1-m)=0.
For simplicity, we consider A=m-aδ and B=(m-aδ)2+4(1-m); then, the two roots of (22) are given by
(23)H+=12(A+B),H-=12(A-B).Lemma 4.
(i)
Suppose thatm-aδ>0 and 1<m<2.
(a)
IfB2>0 holds, model (7) has two positive equilibria E+=(H+,δH+) and E-=(H-,δH-).
(b)
IfB2=0 holds, model (7) has a unique positive equilibrium Ee=(He,δHe). Note that in this case He=H+=H-=A/2=m-1.
(c)
IfB2<0, model (7) has no positive equilibrium.
(ii)
Ifm-aδ≤0, model (7) has no positive equilibrium.LetE=(N,P) be an arbitrary positive equilibrium. The Jacobian matrix of model (7) at the positive equilibrium E=(H,δH) takes the form
(24)JE=(H(aδH+2m-(1+H)2)(1+H)2-aH1+Hrδ-r).
Then, we can get
(25)det(JE)=rH(H2+2H+1+aδ-2m)(H+1)2,tr(JE)=H(aδH+2m-(1+H)2)(1+H)2-r.We can see that the sign ofdet(JE) is determined by
(26)F(H)≜H2+2H+1+aδ-2m=H2+2H+1-A-m.Thus, we can obtain(27)F(H+)=(A+B)24+B-m+1=B(1+A2)+12B2>0;F(H-)=(A-B)24-B-m+1=-12B(A-B)-B<0;F(He)=A4+1-m=0.Hence, we obtaindet(JE+)>0, det(JE-)<0, and det(JEe)=0. And the positive equilibrium E-=(H-,δH-) is a saddle point. The nature of the equilibrium point E+ is dependent on the sign of the trace of the Jacobian matrix evaluated in this point. Whether E+ is a node or a focus depends on the sign of (tr(JE+))2-4det(JE+).In the following results, we study the stability of the positive equilibriumE+ and the unique positive equilibria Ee.Theorem 5.
Define(28)r+=H+(1+H+)2(aδH++2m-(1+H+)2).(a)
Ifr>r+, the positive equilibrium E+=(H+,δH+) is a locally asymptotically stable point;
(a1)
if(r+-r)2<4det(JE+), then E+ is a stable focus,(a2)
if(r+-r)2>4det(JE+), then E+ is a stable node point.(b)
Ifr<r+, the positive equilibrium E+=(H+,δH+) is an unstable point;
(b1)
if(r+-r)2<4det(JE+), then E+ is an unstable focus surrounded by a stable limit cycle,(b2)
if(r+-r)2>4det(JE+), then E+ is an unstable node and the limit cycle disappears.(c)
A Hopf bifurcation occurs atr=r+ around the positive equilibrium E+=(H+,δH+). That is to say, model (7) has at least one positive periodic orbit.Proof.
Here, we only give the proof of the existence of Hopf bifurcation. It is easy to see that(i)
tr(J(E+))|r=r+=0 holds,(ii)
the characteristic equation isλ2+det(J(E+))|r=r+=0, whose roots are purely imaginary,(iii)
(d/dr)[tr(J(E+))]r=r+=-1≠0.
From the Poincaré-Andronov-Hopf Bifurcation Theorem [36], we know that model (7) undergoes a Hopf bifurcation at E+ as r passes through the value r+. The proof is completed.Figure1 illustrates the local stability of the positive equilibrium E+ and the separatrix curve generated by the stable manifold of the positive equilibrium E-. The orbits initiating the right of the separatrix curve tend to E+, while the orbits initiating the left of the separatrix curve tend to (0,0) that represents the extinction of the population. Figure 2 illustrates a Hopf bifurcation situation of the model around E+. The parameter values are given in the figures.Figure 1
The phase portrait of model (7) with the strong Allee effect. The parameters are taken as a=0.25, δ=0.8, m=1.3 and r=0.15. In this case, E0=(0.3,0) is an unstable node point, E1=(1,0) and E-=(0.5,0.4) are saddle points; the positive equilibrium E+=(0.6,0.48) is local asymptotically stable. There exists a separatrix curve determined by the stable manifold of the equilibrium point E-. The dotted curves are the nullclines.Figure 2
The phase portrait of model (7) with the strong Allee effect. The parameters are taken as a=0.25, δ=0.8, m=1.3, and r=0.0375. The model enters into a Hopf bifurcation around E+=(0.6,0.48) at r=r+.Theorem 6.
The unique equilibrium pointEe=(m-1,(1/δ)m-1) is(i)
a nonhyperbolic attractor node, if and only ifr>m-1(1-m-1)2/(1+m-1);(ii)
a nonhyperbolic repellor node, if and only ifr>m-1(1-m-1)2/(1+m-1);(iii)
a cusp point, if and only ifr=m-1(1-m-1)2/(1+m-1), and in this case, there exists a unique trajectory which attains the point Ee. And in this case, the point (0,0) is a global attractor.Proof.
We have(29)tr(JEe)=m-1(1-m-1)21+m-1-r.
Hence (i) and (ii) hold.
Moreover,tr(JEe)=0, if and only if r=m-1(1-m-1)2/(1-m+1). In this case, we obtain the Jacobian matrix of (16) as follows:
(30)J~(m-1,1δm-1)=((m-1)(1-m-1)2-(m-1)(1-m-1)2(m-1)(1-m-1)2-(m-1)(1-m-1)2)=(m-1)(1-m-1)2(1-11-1);
and the associate Jordan matrix is
(31)(0-(m-1)(1-m-1)200).
Then, the singularity (m-1,(1/δ)m-1) is a cusp point, since it is a point of codimension 2, and we have a Bogdanov-Takens bifurcation [37].The cusp point is shown in Figure3.Figure 3
The phase portrait of model (7) with the strong Allee effect. The parameters are taken as a=0.2556936062, δ=0.8, m=1.3, and r=0.07238979895. E0=(0.3,0) is an unstable node point; E1=(1,0) is saddle point; the positive equilibrium E2=(0.2,0.4) is a cusp point. In this case, the point (0,0) is globally asymptotically stable. The dotted curves are the nullclines.


## 2.3. Equilibria Analysis in the Case of the Weak Allee Effect (i.e.,0<m≤1)

In this subsection, we consider the stability of the equilibrium of model (7) with weak Allee effect (0<m≤1).It is easy to verify that model (7) always has one boundary equilibrium E0=(1,0) which is a saddle point and a positive equilibrium E*=(H*,δH*), where
(32)H*=m-aδ+(m-aδ)2+4(1-m)2.From (24), we have
(33)det(JE*)=rH*(H*2+2H*+1+aδ-2m)(H*+1)2>0,tr(JE*)=H*(aδH*+2m-(1+H*)2)(1+H*)2-r.Hence, we have the following results on the stability of the positive equilibriumE*=(H*,δH*).Theorem 7.
Define(34)r*=H*(aδH*+2m-(1+H*)2)(1+H*)2.(a)
Ifr>r*, the positive equilibrium E*=(H*,δH*) is a locally asymptotically stable point, and
(a1)
if(r*-r)2<4det(JE*), then E* is a stable focus;(a2)
if(r*-r)2>4det(JE*), then E* is a stable node point.(b)
Ifr<r*, the positive equilibrium E*=(H*,δH*) is an unstable point, and
(b1)
if(r*-r)2<4det(JE*), then E* is an unstable focus surrounded by a stable limit cycle;(b2)
if(r*-r)2>4det(JE*), then E* is an unstable node and the limit cycle disappears.(c)
A Hopf bifurcation occurs atr=r* around the positive equilibrium E*=(H*,δH*). That is to say, model (7) has at least one positive periodic orbit.In the following theorem, we study the global behavior of the positive equilibriumE*.Theorem 8.
If0<m<1/(1+aδ), the positive equilibrium E*=(H*,δH*) is globally asymptotically stable.Proof.
Construct the following Lyapunov function:(35)V(H,P)=∫H*H‍ξ-H+ξϕ(ξ)dξ+1r∫δH*P‍η-δH*ηdη,
where
(36)ϕ(H)=aHH+1.
Then,
(37)dVdt=H-H*Hϕ(H)dHdt+P-δH*rPdPdt.
Substituting the value ofdH/dt and dP/dt from the model of (7), we obtained
(38)dVdt=H-H*aH[(1-H)(1+H-m)-aδH*]-δH(P-δH*)2.
Note thataδH*=H*(1-H*)(1+H*-m); we obtain
(39)dVdt=-(H-H*)2aH(H+H*-m)-δH(P-δH*)2.
Hence, if 0<m<1/(1+aδ), H*-m>0 which is equivalent to dV/dt<0.
Hence, the positive equilibriumE*=(H*,δH*) is globally asymptotically stable. This completes the proof.Figure4 demonstrates the global stability situation of model (7) around E*. Figure 5 illustrates a Hopf bifurcation situation of the model around E*. Figure 6 shows a stable limit cycle around E* which is an unstable focus. The parameter values are given in the figures.Figure 4
The phase portrait of model (7) with the weak Allee effect. The parameters are taken as a=0.25, δ=0.8, m=1.3, and r=0.1. In this case, E0=(1,0) is saddle point; the positive equilibrium E*=(0.3232928050,0.4849392075) is globally asymptotically stable. The dotted curves are the nullclines.Figure 5
The phase portrait of model (7) with the weak Allee effect. The parameters are taken as a=0.25, δ=0.8, m=1.3, and r=0.025265. The model enters into a Hopf bifurcation around E*=(0.3232928050,0.4849392075) at r=r*.Figure 6
The phase portrait of model (7) with the weak Allee effect. The parameters are taken as a=0.8, δ=1.5, m=0.75, and r=0.025. The positive equilibrium E*=(0.3232928050,0.4849392075) is an unstable focus surrounded by a stable limit cycle.


## 3. Dynamics of the Spatial Model (8)

In this section, we will investigate the dynamics of the spatial model (8). As an example, we only focus on the positive equilibrium point E*=(H*,δH*) in the case of weak Allee effect (0<m<1).


### 3.1. Turing Instability

Mathematically speaking, an equilibrium is Turing instability (diffusion-driven instability) means that it is an asymptotically stable equilibriumE* of model (7) but is unstable with respect to the solutions of reaction-diffusion model (8).In the presence of diffusion, we will introduce small perturbationsU1=H-H*, U2=P-δH*, where |U1|,|U2|≪1. To study the effect of diffusion on the model, we have considered the linearized form of model as follows:
(40)∂U1∂t=r*U1-aH*1+H*U2+d1∇2U1,∂U2∂t=rδU1-rU2+d2∇2U2,
where r* is defined as  (34).Following Malchow et al. [38], we can know that any solution of model (40) can be expanded into a Fourier series so that
(41)U1(r,t)=∑n,m=0∞‍unm(r,t)=∑n,m=0∞‍αnm(t)sinkr,U2(r,t)=∑n,m=0∞‍vnm(r,t)=∑n,m=0∞‍βnm(t)sinkr,
where r=(x,y), and 0<x<L and 0<y<L. k=(kn,km), and kn=nπ/L and km=mπ/L are the corresponding wavenumbers.Having substitutedunm and vnm into (40), we obtain
(42)dαnmdt=(r*-d1k2)αnm-aH*1+H*βnm,dβnmdt=rδαnm-(r+d2k2)βnm,
where k2=kn2+km2.A general solution of (42) has the form C1exp(λ1t)+C2exp(λ2t), where the constants C1 and C2 are determined by the initial conditions (3), and the exponents λ1 and λ2 are the eigenvalues of the following matrix:
(43)D~=(r*-d1k2-aH*1+H*rδ-r-d2k2).Correspondingly,λ1 and λ2 are the solutions of the following equation:
(44)λ2-tr(D~)λ+det(D~)=0,
where
(45)tr(D~)=r*-r-(d1+d2)k2,det(D~)=d1d2k4+(rd1-r*d2)k2+det(J(E*)).Summarizing the previous discussions, we can get the following theorem immediately.Theorem 9.
(i) The positive equilibriumE* of model (8) is locally asymptotically stable if r>max{r*,r*d2/d1} holds.
(ii) If the positive equilibriumE* of model (7) is globally asymptotically stable, then the corresponding steady state E* of model (8) is also globally asymptotically stable.Proof.
(i) Using Routh-Hurwitz criteria, we can know that the positive equilibriumE* is locally asymptotically stable, if and only if tr(D~)<0 and det(D~)>0. So, we obtain r>max{r*,r*d2/d1}.
(ii) We select the Lyapunov function for model (8) as follows:
(46)V2(t)=∬Ω‍V(H,P)dxdy,
where V(H,P) is the same as defined in (35). So,
(47)dV2dt=∬ΩdVdtdxdy+∬Ω{∂V∂Hd1∇2H+∂V∂Pd2∇2P}dxdy.
Using Green’s first identity in the plane,(48)∬ΩF∇2Gdxdy=∫∂Ω‍F∂G∂nds-∬Ω‍(∇F·∇G)dxdy.
Considering the zero-flux boundary conditions, one can show that
(49)∬Ω‍∂V∂Hd1∇2Hdxdy=-d1∬Ω‍∂2V∂H2[(∂H∂x)2+(∂H∂y)2]dxdy≤0,∬Ω‍∂V∂Pd2∇2Pdxdy=-d2∬Ω‍∂2V∂P2[(∂P∂x)2+(∂P∂y)2]dxdy≤0.
From the previous analysis, we note thatdV2/dt<0 is valid if dV/dt<0 is true. This implies that the equilibrium E* of both model (7) and model (8) is globally asymptotically stable if 0<m<1/(1+aδ) holds. This ends the proof.On the other hand, Turing instability sets in when at least one of the conditions is eithertr(D~)<0 or det(D~)>0. It is evident that the condition tr(D~)<0 is not violated when the requirement r*-r<0 is met [39]. Hence, only violation of condition det(D~)>0 gives rise to diffusion instability. Then, the condition for diffusive instability is given by
(50)G(k2)≡d1d2k4+(rd1-r*d2)k2+det(J(E*))<0.G
(
·
) is quadratic in k2, and the graph of G(k2)=0 is a parabola. The minimum of G(k2) occurs at k2=km2, where
(51)km2=r*d2-rd12d1d2>0.
The critical wave number kc of the first perturbations to grow is found by evaluating km from (51).Thus, a sufficient condition for Turing instability is thatG(km2) is negative. Therefore,
(52)G(km2)=det(J(E*))-(rd1-r*d2)24d1d2<0.
Combination of (51) and (52) leads to the following final criterion for diffusive instability:
(53)(rd1-r*d2)2>4d1d2det(J(E*)).Summarizing the previous discussions, we can obtain the following theorem.Theorem 10.
Ifrd1/d2<r*<r and r*d2-rd1>2d1d2det(J(E*)) hold, the criterion for Turing instability for model (8) emerges, and the critical wave number kc=(r*d2-rd1)/2d1d2.The Turing instability (or bifurcation) breaks spatial symmetry, leading to the formation of patterns that are stationary in time and oscillatory in space [40, 41]. We adopt the intrinsic growth rates of predator r as the bifurcation parameter, and the linear stability analysis yields the bifurcation diagram shown in Figure 7. The Turing bifurcation curve separates the parametric space into two domains. Above the curve, the solutions of model (8) are stable for all pairs of (m,r); that is, there is no Turing instability. While below the curve, the solutions of model (8) are unstable for (m,r) and diffusive instability emerges; that is, Turing patterns emerge. This domain is called the Turing space.Figure 7
Turing bifurcation diagram for model (8) using m and r as parameters. Other parameters are taken as a=0.8, δ=1.75, d1=0.01,  and d2=1. Above the curve, the positive equilibrium E* is the only stable solution of model (3). Below the curve, the positive equilibrium E* loses its stability with respect to model (3), and Turing instability occurs; this domain is called the Turing space.


### 3.2. Pattern Formation

In this subsection, we performed extensive numerical simulations of the spatially extended model (8) in two-dimension spaces, and the qualitative results are shown here. All our numerical simulations employ the zero-flux boundary conditions with a model size of L×L, with L=100 discretized through x→(x0,x1,x2,…,xn) and y→(y0,y1,y2,…,yn), with n=200. Other parameters are fixed as a=0.8, δ=1.75, r=0.15,  d1=0.01,  and d2=1.The numerical integration of model (8) is performed by using a finite difference approximation for the spatial derivatives and an explicit Euler method for the time integration [42] with a time stepsize of τ=1/100. The initial condition is always a small amplitude random perturbation around the positive constant steady state solution E*. After the initial period during which the perturbation spread, either the model goes into a time-dependent state or to an essentially steady state solution (time independent).More precisely, the concentrations(Hi,jn+1,Pi,jn+1) at the moment (n+1)τ at the mesh position (i,j) are given by
(54)Hi,jn+1=Hi,jn+τd1ΔhHi,jn+τf(Hi,jn,Pi,jn),Pi,jn+1=Pi,jn+τd2ΔhPi,jn+τg(Hi,jn,Pi,jn),
with the Laplacian defined by
(55)ΔhHi,jn=Hi+1,jn+Hi-1,jn+Hi,j+1n+Hi,j-1n-4Hi,jnh2,
where f(H,P)=H(1-H)(1-m/(H+1))-aHP/(1+H), g(H,P)=rP(1-P/δH), 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 preyH, for instance.Figure8 shows the evolution of the spatial pattern of prey H at t=0, 500, 1000, 2000, with small random perturbation of the stationary solution E* of the spatially homogeneous systems when m is located in “Turing space.” In this case, one can see that for model (8), the random initial distribution leads to the formation of a strongly irregular transient pattern in the domain. After the irregular pattern is formed (c.f., Figures 8(b) and 8(c)), it grows slightly and jumps alternately for a certain time, and finally spots patterns, which are isolated zones with low prey densities, prevail over the whole domain, and the dynamics of the model does not undergo any further changes (c.f., Figure 8(d)).Spots pattern ofH in model (8) for m=0.75. Times: (a) 0, (b) 500, (c) 1000, and (d) 2000.
(a)
(b)
(c)
(d)Figure9 shows stripe patterns are interlaced stripes of high and low population densities of prey H for the parameter m=0.78 at t=1000. In Figure 10, with the parameter m=0.85, the spot-stripe mixtures pattern is time independent with low prey densities.Figure 9
Stripes pattern ofH in model (8) for m=0.78 at t=1000.Figure 10
Spot-stripe mixtures pattern ofH in model (8) for m=0.85 at t=1000.


## 3.1. Turing Instability

Mathematically speaking, an equilibrium is Turing instability (diffusion-driven instability) means that it is an asymptotically stable equilibriumE* of model (7) but is unstable with respect to the solutions of reaction-diffusion model (8).In the presence of diffusion, we will introduce small perturbationsU1=H-H*, U2=P-δH*, where |U1|,|U2|≪1. To study the effect of diffusion on the model, we have considered the linearized form of model as follows:
(40)∂U1∂t=r*U1-aH*1+H*U2+d1∇2U1,∂U2∂t=rδU1-rU2+d2∇2U2,
where r* is defined as  (34).Following Malchow et al. [38], we can know that any solution of model (40) can be expanded into a Fourier series so that
(41)U1(r,t)=∑n,m=0∞‍unm(r,t)=∑n,m=0∞‍αnm(t)sinkr,U2(r,t)=∑n,m=0∞‍vnm(r,t)=∑n,m=0∞‍βnm(t)sinkr,
where r=(x,y), and 0<x<L and 0<y<L. k=(kn,km), and kn=nπ/L and km=mπ/L are the corresponding wavenumbers.Having substitutedunm and vnm into (40), we obtain
(42)dαnmdt=(r*-d1k2)αnm-aH*1+H*βnm,dβnmdt=rδαnm-(r+d2k2)βnm,
where k2=kn2+km2.A general solution of (42) has the form C1exp(λ1t)+C2exp(λ2t), where the constants C1 and C2 are determined by the initial conditions (3), and the exponents λ1 and λ2 are the eigenvalues of the following matrix:
(43)D~=(r*-d1k2-aH*1+H*rδ-r-d2k2).Correspondingly,λ1 and λ2 are the solutions of the following equation:
(44)λ2-tr(D~)λ+det(D~)=0,
where
(45)tr(D~)=r*-r-(d1+d2)k2,det(D~)=d1d2k4+(rd1-r*d2)k2+det(J(E*)).Summarizing the previous discussions, we can get the following theorem immediately.Theorem 9.
(i) The positive equilibriumE* of model (8) is locally asymptotically stable if r>max{r*,r*d2/d1} holds.
(ii) If the positive equilibriumE* of model (7) is globally asymptotically stable, then the corresponding steady state E* of model (8) is also globally asymptotically stable.Proof.
(i) Using Routh-Hurwitz criteria, we can know that the positive equilibriumE* is locally asymptotically stable, if and only if tr(D~)<0 and det(D~)>0. So, we obtain r>max{r*,r*d2/d1}.
(ii) We select the Lyapunov function for model (8) as follows:
(46)V2(t)=∬Ω‍V(H,P)dxdy,
where V(H,P) is the same as defined in (35). So,
(47)dV2dt=∬ΩdVdtdxdy+∬Ω{∂V∂Hd1∇2H+∂V∂Pd2∇2P}dxdy.
Using Green’s first identity in the plane,(48)∬ΩF∇2Gdxdy=∫∂Ω‍F∂G∂nds-∬Ω‍(∇F·∇G)dxdy.
Considering the zero-flux boundary conditions, one can show that
(49)∬Ω‍∂V∂Hd1∇2Hdxdy=-d1∬Ω‍∂2V∂H2[(∂H∂x)2+(∂H∂y)2]dxdy≤0,∬Ω‍∂V∂Pd2∇2Pdxdy=-d2∬Ω‍∂2V∂P2[(∂P∂x)2+(∂P∂y)2]dxdy≤0.
From the previous analysis, we note thatdV2/dt<0 is valid if dV/dt<0 is true. This implies that the equilibrium E* of both model (7) and model (8) is globally asymptotically stable if 0<m<1/(1+aδ) holds. This ends the proof.On the other hand, Turing instability sets in when at least one of the conditions is eithertr(D~)<0 or det(D~)>0. It is evident that the condition tr(D~)<0 is not violated when the requirement r*-r<0 is met [39]. Hence, only violation of condition det(D~)>0 gives rise to diffusion instability. Then, the condition for diffusive instability is given by
(50)G(k2)≡d1d2k4+(rd1-r*d2)k2+det(J(E*))<0.G
(
·
) is quadratic in k2, and the graph of G(k2)=0 is a parabola. The minimum of G(k2) occurs at k2=km2, where
(51)km2=r*d2-rd12d1d2>0.
The critical wave number kc of the first perturbations to grow is found by evaluating km from (51).Thus, a sufficient condition for Turing instability is thatG(km2) is negative. Therefore,
(52)G(km2)=det(J(E*))-(rd1-r*d2)24d1d2<0.
Combination of (51) and (52) leads to the following final criterion for diffusive instability:
(53)(rd1-r*d2)2>4d1d2det(J(E*)).Summarizing the previous discussions, we can obtain the following theorem.Theorem 10.
Ifrd1/d2<r*<r and r*d2-rd1>2d1d2det(J(E*)) hold, the criterion for Turing instability for model (8) emerges, and the critical wave number kc=(r*d2-rd1)/2d1d2.The Turing instability (or bifurcation) breaks spatial symmetry, leading to the formation of patterns that are stationary in time and oscillatory in space [40, 41]. We adopt the intrinsic growth rates of predator r as the bifurcation parameter, and the linear stability analysis yields the bifurcation diagram shown in Figure 7. The Turing bifurcation curve separates the parametric space into two domains. Above the curve, the solutions of model (8) are stable for all pairs of (m,r); that is, there is no Turing instability. While below the curve, the solutions of model (8) are unstable for (m,r) and diffusive instability emerges; that is, Turing patterns emerge. This domain is called the Turing space.Figure 7
Turing bifurcation diagram for model (8) using m and r as parameters. Other parameters are taken as a=0.8, δ=1.75, d1=0.01,  and d2=1. Above the curve, the positive equilibrium E* is the only stable solution of model (3). Below the curve, the positive equilibrium E* loses its stability with respect to model (3), and Turing instability occurs; this domain is called the Turing space.


## 3.2. Pattern Formation

In this subsection, we performed extensive numerical simulations of the spatially extended model (8) in two-dimension spaces, and the qualitative results are shown here. All our numerical simulations employ the zero-flux boundary conditions with a model size of L×L, with L=100 discretized through x→(x0,x1,x2,…,xn) and y→(y0,y1,y2,…,yn), with n=200. Other parameters are fixed as a=0.8, δ=1.75, r=0.15,  d1=0.01,  and d2=1.The numerical integration of model (8) is performed by using a finite difference approximation for the spatial derivatives and an explicit Euler method for the time integration [42] with a time stepsize of τ=1/100. The initial condition is always a small amplitude random perturbation around the positive constant steady state solution E*. After the initial period during which the perturbation spread, either the model goes into a time-dependent state or to an essentially steady state solution (time independent).More precisely, the concentrations(Hi,jn+1,Pi,jn+1) at the moment (n+1)τ at the mesh position (i,j) are given by
(54)Hi,jn+1=Hi,jn+τd1ΔhHi,jn+τf(Hi,jn,Pi,jn),Pi,jn+1=Pi,jn+τd2ΔhPi,jn+τg(Hi,jn,Pi,jn),
with the Laplacian defined by
(55)ΔhHi,jn=Hi+1,jn+Hi-1,jn+Hi,j+1n+Hi,j-1n-4Hi,jnh2,
where f(H,P)=H(1-H)(1-m/(H+1))-aHP/(1+H), g(H,P)=rP(1-P/δH), 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 preyH, for instance.Figure8 shows the evolution of the spatial pattern of prey H at t=0, 500, 1000, 2000, with small random perturbation of the stationary solution E* of the spatially homogeneous systems when m is located in “Turing space.” In this case, one can see that for model (8), the random initial distribution leads to the formation of a strongly irregular transient pattern in the domain. After the irregular pattern is formed (c.f., Figures 8(b) and 8(c)), it grows slightly and jumps alternately for a certain time, and finally spots patterns, which are isolated zones with low prey densities, prevail over the whole domain, and the dynamics of the model does not undergo any further changes (c.f., Figure 8(d)).Spots pattern ofH in model (8) for m=0.75. Times: (a) 0, (b) 500, (c) 1000, and (d) 2000.
(a)
(b)
(c)
(d)Figure9 shows stripe patterns are interlaced stripes of high and low population densities of prey H for the parameter m=0.78 at t=1000. In Figure 10, with the parameter m=0.85, the spot-stripe mixtures pattern is time independent with low prey densities.Figure 9
Stripes pattern ofH in model (8) for m=0.78 at t=1000.Figure 10
Spot-stripe mixtures pattern ofH in model (8) for m=0.85 at t=1000.


## 4. Concluding Remarks

In this paper, we are concerned with the complex dynamics in a diffusive Holling-Tanner predator-prey model with the Allee effect on prey. The value of this study lies in two folds. First, the local asymptotic stability conditions for coexisting equilibrium and conditions for Hopf bifurcation are described briefly for the model with the weak and strong Allee effects. Second, it gives the analysis of Turing instability which determines the Turing space in the spatial domain and meanwhile illustrates the Turing pattern formation close to the onset Turing bifurcation via numerical simulations, which shows that the model dynamics exhibits complex pattern replication.We note that in the analyzed models, a big difference between the dynamics of model with strong or weak Allee effect exists. In the case of strong Allee effect, two positive equilibria can coexist for a subset of parameters with a varied dynamics but different to other Holling-Tanner models analyzed earlier [12–14]. We have shown that one of these equilibria is always a saddle point and proved the existence of a separatrix curve. And there is no global asymptotically stable positive equilibrium. In this case, the point (0,0) is an attractor in addition to locally stable positive equilibrium E+ for determined parameter values, which leads to the existence of a bistability phenomenon. The dynamics of the model is determined by the initial conditions; the predator and prey may be extinction or coexistence. This means that the strong Allee effect could easily lead to the risk of population extinction. Nevertheless, in the case of weak Allee effect, model (7) can only have one unique positive equilibrium, which is globally asymptotically stable under some conditions. Therefore, the predators and preys can coexist in stable conditions.Furthermore, we have investigated the conditions for the predator-prey model which experiences spatial patterns through diffusion-driven instability. We have derived the conditions of Turing instability in terms of our model parameters analytically. In addition, to get a deeper insight into the model’s dynamics behaviour, we select the different values of parameterm. An increase of m, from the numerical results, one can see that our model has rich and complex spatiotemporal behavior. We find three typical Turing patterns, that is, spots pattern, stripes pattern, and spots-stripes mixtures pattern. To the best of our knowledge, the Turing pattern we illustrate here is the first reported case to our model. And our complete analysis of the spatial model will give new suggestion to the models with the Allee effect.



---

*Source: 270191-2013-02-11.xml*