# Pattern Dynamics in a Spatial Predator-Prey System with Allee Effect

Authors: Gui-Quan Sun; Li Li; Zhen Jin; Zi-Ke Zhang; Tao Zhou
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We investigate the spatial dynamics of a predator-prey system with Allee effect. By using bifurcation analysis, the exact Turing domain is found in the parameters space. Furthermore, we obtain the amplitude equations and determine the stability of different patterns. In Turing
space, it is found that predator-prey systems with Allee effect have rich dynamics. Our results indicate that predator mortality plays an important role in the pattern formation of populations. More specifically, as predator mortality rate increases, coexistence of spotted and stripe patterns, stripe patterns, spotted patterns, and spiral wave emerge successively. The results enrich the finding in the spatial predator-prey systems well.


---

## Body

## 1. Introduction

The Allee effect, named after the ecologist Warder Clyde Allee, has been recognized as an important phenomenon of positive density dependence in low-density population [1–5]. Allee effect can occur whenever fitness of an individual in a small or sparse population decreases as the population size or density also declines [6, 7]. Since the outstanding work of Allee [1], the Allee effect has been regarded as one of the central and highly important issues in the population and community ecology. And its critical importance has widely been realized in the conservation biology that Allee effect is most likely to increase the extinction risk of low-density populations. As a result, studies on Allee effect have received more and more attention from both mathematicians and ecologists.Long time series of the density of both prey and predator is needed, so it is difficult to analyse their dynamics. As a result, it may provide useful information by constructing mathematical models to investigate the dynamical behaviors of predator-prey systems. There have been a large group of papers on predator-prey systems with Allee effect [8–13]. However, these previous works did not take into account the effect of space.There are also some works done on spatial predator-prey systems with Allee effect [14–16]. Petrovskii et al. found that the deterministic system with Allee effect can induce patch invasion [14]. Morozov et al. found that the temporal population oscillations can exhibit chaotic dynamics even when the distribution of the species in the space was regular [15]. Moreover, they found that the chaos accompanied with patch invasion even though the environments were heterogeneous [16]. However, their results were obtained by choosing particular initial conditions. Then, it is natural to ask what kind of patterns can be obtained in predator-prey systems with Allee effect by using other initial conditions. To understand that mechanism well, we will investigate a predator-prey system with Allee effect.Because of the insightful work of many scientists over recent years, we can make research on pattern selection by using the standard multiple scale analysis [17, 18], in which the control parameters and the derivatives are expanded in terms of a small enough parameter. In the neighborhood of the bifurcation points (Hopf and Turing bifurcation points), the critical amplitudes follow the normal forms, and thus their general forms can be obtained from the methods of symmetry-breaking bifurcations.The paper is organized as follows. In Section2, we present a predator-prey system with Allee effect and give Turing region in parameters space. In Section 3, by using multiple scale analysis, we obtain amplitude equations. In Section 4, we show the spatial patterns by a series of numerical simulations. Finally, conclusions and discussions are presented in Section 5.


## 2. A Predator-Prey System with Allee Effect

We consider the following model of two-dimensional spatiotemporal system [14–16, 19]:(1a)∂H∂T=F(H)-f(H,P)+D1ΔH,(1b)∂P∂T=κf(H,P)-D(P)+D2ΔP,whereH=H(X,Y,T)and  P=P(X,Y,T)  are densities of prey and predator, respectively, at time T and position  (X,Y). The function  F(H)  represents the intrinsic prey growth,  f(H,P)=f(H)P  represents predation term,κis the food utilization coefficient, D1 and D2 are diffusion coefficients, and D(P) describes predator mortality.It is assumed that the predation term is a bilinear form of prey and predator density and predator mortality is a nonlinear function of predator density. As a result, we choosef(H,P)=HP and D(P)=MP2 [20].When the prey population obeys Allee dynamics, its growth rate can be parameterized as follows [14, 15, 21]:
(2)F(H)=4ω(K-H0)2H(H-H0)(K-H),
whereKis the prey-carrying capacity,ωis the maximum per capita growth rate, andH0quantifies the intensity of the Allee effect. If  0<H0<K,  F(H)  is a strong Allee effect; if -K<H0<0, F(H) is a weak Allee effect; if H0≤-K, the Allee effect is absent.In order to minimize the number of parameters involved in the model system, it is extremely useful to write the system in a nondimensionalized form. Although there is no unique method of doing this, it is often a good idea to relate the variables to some key relevant parameters. Introducing dimensionless variables(3)u=HK,v=PκK,t=aT,X¯=XaD1,Y¯=YaD1,
we obtain the following equations:(4a)∂u∂t=γu(u-β)(1-u)-uv+Δu,(4b)∂v∂t=uv-δv2+εΔv,where(5)β=H0K,γ=4ωKAκ(K-H0)2,δ=Ma,ε=D2D1.
First of all, we need to investigate the dynamics of nonspatial  model of systems (4a) and (4b)(6a)dudt=γu(u-β)(1-u)-uv,(6b)dvdt=uv-δv2.Systems (6a) and (6b) have three boundary equilibrium named   E0=(0,0),  E1=(1,0), and E2=(β,0)and two interior equilibriums named   E3 and E4, where (7a)E3=(γδ+γβδ-1+Q2γδ,γδ+γβδ-1+Q2γδ2),(7b)E4=(γδ+γβδ-1-Q2γδ,γδ+γβδ-1-Q2γδ2),whereQ=(γδ)2-2(γδ)2β-2γβδ+1.From a biological point of view, we are concerned with the dynamics ofE3andE4. The Jacobian matrix corresponding to the equilibrium point is that
(8)J=(a11a12a21a22),
where
(9)a11=2γu*-γβ-3γ(u*)2+2γβu*,a12=-u*,a21=v*-δ,a22=u*.Diffusion-driven instability requires the stable, homogeneous steady state is driven unstable by the interaction of the dynamics and diffusion of the species; and therefore(10)a11+a22<0,a11a22-a12a21>0.It is found from direct calculations thatE3 is unstable and E4 is stable. Denote E4=(u*,v*).Following the standard linear analysis of the reaction-diffusion equation [22], we consider a perturbation near the steady state:
(11)u(r→,t)=u*+u¯(r,t),v(r→,t)=v*+v¯(r,t),
where u¯(r,t)≪u*,  v¯(r,t)≪v*, and  r=(X¯,Y¯). Assume that
(12)(u¯(r,t)v¯(r,t))=(α1α2)eλtei(κXX+κYY),
where  λis the growth rate of perturbation in time  t, α1 and α2 represent the amplitudes, and κXand κY are the wave number of the solutions.The characteristic equation of the systems (4a) and (4b) is
(13)(A-λI)(u¯v¯)=0,
where
(14)A=(a11-(κX2+κY2)a12a21a22-ε(κX2+κY2)).As a result, we have characteristic polynomial:(15)λ2-trκλ+Δκ=0,(16)trκ=a11+a22-κ2(1+ε)=ΔtrJ-κ2(1+ε),Δκ=a11a22-a12a21-κ2(a11ε+a22)+κ4ε=ΔΔJ-κ2(a11ε+a22)+κ4ε,
where  κ2=κX2+κY2.The roots of (15) can be obtained by the following form:
(17)λκ=12(trκ±trκ2-4Δκ).WhenIm(λκ)≠0  and  Re(λκ)=0, Hopf bifurcation will emerge. Then, we have that the critical value of Hopf bifurcation parameter-δ equals
(18)δH=γ(γ+β-1)γ2β2+γ2-2γ2β-1.Whenκ2=(κT)2=ΔJ/ε  and  Im(λκ)=0, Re(λκ)=0, Turing bifurcation will occur. Denote δT as the critical value ofδas Turing instability occurs. Since the expression is complicated, we omit it here.In Figure1, we show the two critical lines in the parameter space spanned by β andδ. The equilibria that can be found in the region, marked by T (Turing space), are stable with respect to the homogeneous perturbations, but they lose their stability with respect to the perturbations of specific wave numbers κ. In this region, stationary patterns can be observed. To see the effect of parameter  δ well, we plot in Figure 2 the dispersion relation corresponding to several values of δ  while keeping the other parameters fixed. We see that the available Turing modes shift to higher wave numbers when  β decreases.Figure 1
Bifurcation diagram for the systems (4a) and (4b). The green one is the Hopf bifurcation critical line and the red one, Turing bifurcation critical line. The figure shows the Turing space which is marked by T. Parameters values:  γ=1.5andε=0.15.Figure 2
Dispersion relation for differentδ. Parameters values:  β=0.02,  γ=1.5, and  ε=0.15. (a)  δ=1.08; (b)  δ=1.04; (c)  δ=1; (d)  δ=0.96; and (e)  δ=0.92.


## 3. Spatial Dynamics of Systems(4a) and  (4b)

In the following, we use multiple scale analysis to determine the amplitude equations when|κ|=κT. Denote δ as the controlled parameters. When the controlled parameter is larger than the critical value of Turing point, the solutions of systems (4a) and (4b) can be expanded as
(19)c=c0+∑i=1N‍(Aiexp(iκir→)+(A¯iexp(-iκir→)),
with|κ|=κT. Aj and the conjugate A-j are the amplitudes associated with the modes κj and -κj.Close to onsetβ=βT, one has that
(20)∂Ai∂t=siAi+Fi(Ai,Aj,…).Based on the center manifold near the Turing bifurcation point, it can be concluded that amplitudeAjsatisfies
(21)∂Ai∂t=Fi(Ai,A¯i,Aj,A¯j,…).In order to obtain the amplitude equations, we first need to investigate the linearized form of systems (4a) and (4b) at the equilibrium point E4. By setting u=u*+x  and  v=v*+y, we have the following equations:(22a)∂x∂t=[2γu*-3γ(u*)2+2γ(u*)2β-γβ-v*]x+(βγ-3γu*+γ)x2-γx3-xy+Δx,(22b)∂y∂t=v*x+u*y+xy-2δv*y-δy2+εΔy.Close to onsetδ=δT, the solutions of systems (4a) and (4b) can be expanded as series form:
(23)U=Us+∑j=13‍U0[Ajexp(iκjr→)+A¯jexp(-iκjr→)].System (19) can be expanded as
(24)U*=∑j=13‍U0[Ajexp(iκjr→)+A¯jexp(-iκjr→)],
where  U0=((a11*ε+a11*)/(2a21*),1)T  is the eigenvector of the linearized operator.From the standard multiple scale analysis, up to the third order in the perturbations, the spatiotemporal evolution of the amplitudes can be described as(25)τ∂Ak∂t=μAk+∑lm‍hlmAlAm+∑lmn‍glmnAlAmAn.Due to spatial translational symmetry, we have the following equation:(26)τ∂Ak∂texp(iκkr0)=μAkexp(iκkr0)+∑lm‍hlmAlAmexp[i(κl+κm)r0]+∑lmn‍glmnAlAmAnexp[i(κl+κm+κn)r0].Comparing (25) with (26), one can find that the two equations hold only if  κk=κl+⋯+κm. From the center manifold theory, we know that amplitude equations do not include the amplitude with unstable mode. As a result, we have the following equations:
(27)τ0∂A1dt=μA1+hA-2A-3-(g1|A1|2+g2(|A2|2+|A3|2))A1,τ0∂A2dt=μA2+hA-1A-3-(g1|A2|2+g2(|A1|2+|A3|2))A2,τ0∂A3dt=μA3+hA-1A-2-(g1|A3|2+g2(|A1|2+|A2|2))A3,
where μ=(δT-δ)/δT and  τ0  is a typical relaxation time.In the following part, we will give the expressions ofτ0,  h,  g1, and  g2. Let
(28)X=(xy),N=(N1N2).
Then systems (4a) and (4b) can be written as:
(29)∂X∂t=LX+N,
where
(30)L=(2γu*-3γ(u*)2+2γβ(u*)2-γβ-v*+Δ0v*u*-2δv*+εΔ),N=((βγ-3γu*+γ)-γx3-xyxy-δy2).We need to investigate the dynamical behavior whenδis close toδT, and thus we expandδas:
(31)δT-δ=ϵδ1+ϵ2δ2+ϵ3δ3+O(ϵ4),
where ϵ is a small enough parameter. We expand  X and  N as the series form of  ϵ:
(32)X=(xy)=ϵ(x1y1)+ϵ2(x2y2)+ϵ3(x3y3)+⋯,N=((βγ-3γu*+γ)(x12ϵ2+2x1x2ϵ3)-γx13ϵ3-x1y1ϵ2-(x2y1+x1y2)ϵ3+o(ϵ4)x1y1ϵ2+(x2y1+x1y2)ϵ3+o(ϵ4)).
Linear operator L can be expanded as
(33)L=LT+(δT-δ)M,
where
(34)LT=(a11*+Δa12*a21*a22*+εΔ),M=(b11b12b21b22).Let(35)T0=t,T1=ϵt,T2=ϵ2t,
and Ti is a dependent variable. For the derivation of time, we have that
(36)∂∂t=∂∂T0+ϵ∂∂T1+ϵ2∂∂T2+o(ϵ3).The solutions of systems (4a) and (4b) have the following form:
(37)X=(xy)=∑i=13‍(AixAiy)exp(iκir→)+⋯.
This expression implies that the bases of the solutions have nothing to do with time and the amplitude A is a variable that changes slowly. As a result, one has the following equation:
(38)∂A∂t=ϵ∂A∂T1+ϵ2∂A∂T2+o(ϵ3).Substituting the above equations into (29) and expanding (29) according to different orders of ϵ, we can obtain three equations as follows:
(39)ϵ:LT(x1y1)=0;ϵ2:LT(x2y2)=∂∂T1(x1y1)-δ1M(x1y1)-((βγ-3γu*+γ)x12-x1y1x1y1);ϵ3:LT(x3y3)=∂∂T1(x2y2)+∂∂T2(x1y1)-δ1M(x2y2)-δ2M(x1y1)-E,
where
(40)E=(2x1x2(βγ-3γu*+γ)-γx13-(x2y1+x1y2)x2y1+x1y2).We first consider the case of the first order ofε. Since  LT is the linear operator of the system close to the onset,  (x1,y1)T is the linear combination of the eigenvectors that corresponds to the eigenvalue zero. Since that
(41)(xy)=∑i=13‍(AixAiy)exp(iκir→)+c.c.,
we have that (42a)(a11*+Δ)x1+a12*y1=0,(42b)a21*x1+(a22*+εΔ)y1=0.Asεa12*=((a22*-εa11*)/2a21*)2, we can obtain that  x1=(a22*-εa11*)/(2a21*) by assuming  y1=1.LetR=(a11*ε-a22*)/2a21*  then
(43)(x1y1)=(R1)(W1exp(iκ1r→)+W2exp(iκ2r→)+W3exp(iκ3r→))+c.c.,
where  |κj|=κT*  and  Wj  is the amplitude of the mode exp(iκjr).Now, we consider the case of the second order ofε. Note that
(44)LT(x2y2)=∂∂T1(x1y1)-δT(b11x1+b12y1b21x1+b22y1)-((βγ-3γu*+γ)x12-x1y1x1y1)=(FxFy).According to the Fredholm solubility condition, the vector function of the right hand of the above equation must be orthogonal with the zero eigenvectors of operatorLc+. And the zero eigenvectors of operator Lc+ are
(45)(1-1εR)exp(iκjr→)+c.c(j=1,2,3).It can be found from the orthogonality condition that(46)(1-1εR)(FxiFyi)=0,
where  Fxi and  Fyi represent the coefficients corresponding to exp(iκjr) in Fx and  Fy.By investigatingexp(iκ1r→), one has
(47)(Fx1Fy1)=(R∂W1∂T1∂W1∂T1)-δ1(b11RW1+b12W1b21RW1+b22W1)-(2R2(βγ-3γu*+γ)W¯2W¯3+2RW¯2W¯32RW¯2W¯3).It can be obtained from the orthogonality condition that(48)ε-1εR∂W1∂T1=δ(Rb11+b22-Rε(Rb21+b22)W1)+2R2(βγ-3γu*+γ+1R-1ε)W¯2W¯3.By using the same methods, one has(49)(x2y2)=(X0Y0)+∑j=13‍(XjYj)exp(iκjr→)+∑j=13‍(XjjYjj)exp(2iκjr→)+Q+c.c.,
where
(50)Q=(X12Y12)exp(i(κ1-κ2)r→)+(X23Y23)exp(i(κ2-κ3)r→)+(X31Y31)exp(i(κ3-κ1)r→).By solving the sets of the linear equations aboutexp(0),  exp(iκjr→),  exp(2iκjr→), and exp(i(κj-κk)r→), we obtain that
(51)(X0Y0)=(a22*[-2R2(βγ-3γu*+γ)+2R]+2Ra12*a11*a22*-a12*a21*a21*[2R2(βγ-3γu*+γ)-2R]+2Ra11*a11*a22*-a12*a21*)×(|W1|2+|W2|2+|W3|2),(XjjYjj)=((a22*-4εκT2)[-R2(βγ-3γu*+γ)+R]+Ra12*(a11*-4κT2)(a22*-4εκT2)-a12*a21*a21*[R2(βγ-3γu*+γ)-2R]-R(a11*-4εκT2)(a11*-4κT2)(a22*-4εκT2)-a12*a21*)×Wj2,(XjkYjk)=((a22*-3εκT2)[-2R2(βγ-3γu*+γ)+2R]+2Ra12*(a11*-3κT2)(a22*-3εκT2)-a12*a21*a21*[2R2(βγ-3γu*+γ)-2R]-2R(a11*-3κT2)(a11*-3κT2)(a22*-3εκT2)-a12*a21*)×WjW¯k,
where κT2=(a11*a22*-a12*a21*)/ε.For the third order ofε, we have that
(52)LT(x3y3)=∂∂T1(x2y2)+∂∂T2(x1y1)-δ1M(x2y2)-δ2M(x1y1)-S,
where
(53)S=(2x1x2(βγ-3γu*)-γx13-(x2y1+x1y2)x2y1+x1y2).Using the Fredholm solubility condition, we can obtain(54)ε-1εR∂W1∂T2+ε-1εR∂Y1∂T1=δ2[Rb11+b12-1εR(Rb21+b22)]W1+δ1[Rb11+b12-1ε(Rb21+b22)]Y1+Z,
where
(55)Z=[2R2(βγ-3γu*+γ)-2R-21ε]×[W1Y0+W-2Y12+W3Y13+W¯1Y11+W¯2Y¯3+W¯3Y¯2]-(G1|W1|2+G2|W2|2+G3|W3|2)W1,G1=(1εR-1)[R(y11+y0)+x11+x0]-2R(x11+x0)(βγ-3γu*+γ)+3γR3,G2=(1εR-1)[R(y12+y0)+x12+x0]-2R(x12+x0)(βγ-3γu*+γ)+6γR3.By using the same methods, we can obtain the other two equations. The amplitudeAican be expanded as
(56)Ai=ϵWi+ϵ2Vi+o(ϵ3).As a result, we have(57)τ0∂A1∂t=μA1+hA-2A-3-(g1|A1|2+g2|A2|2+|A3|2)A1.The other two equations can be obtained through the transformation of the subscript ofA. By calculations, we obtain the expressions of the coefficients ofτ0,h,g1, andg2as follows:
(58)τ0=Rε-1δT[Rb11+b12-(R/ε)(Rb21+b22)],h=[2R2(βγ-3γu*+γ)-2R-2(R2/ε)]δT[Rb11+b12-(R/ε)(Rb21+b22)],g1=G1δT[Rb11+b12-(R/ε)(Rb21+b22)],g2=G2δT[Rb11+b12-(R/ε)(Rb21+b22)],
whereG1=((R/ε)-1)[R(y0+y11)+x0+x11]-2R(x0+x11)(βγ-3γu*+γ)+3γR3 and G2=((R/ε)-1)[R(y0+y12)+x0+x12]-2R(x0+x12)(βγ-3γu*+γ)+6γR3.By using substitutions, we have(59)τ0∂φdt=-hρ12ρ22+ρ12ρ32+ρ22ρ32ρ1ρ2ρ3sinφ,τ0∂ρ1dt=μρ1+hρ2ρ3cosφ-g1ρ13-g2(ρ22ρ32)ρ1,τ0∂ρ2dt=μρ2+hρ1ρ3cosφ-g1ρ23-g2(ρ12ρ32)ρ2,τ0∂ρ3dt=μρ3+hρ1ρ2cosφ-g1ρ33-g2(ρ12ρ22)ρ3,
where φ=φ1+φ2+φ3. In order to see the relationships between different parameters, we give the values of coefficients for different parameter sets in Table 1.Table 1
Coefficients for different parameter sets.
β
δ
h
g
1
g
2
μ
1
μ
2
μ
3
μ
4
0.02
0.92
−19.08604
7599.215
6906.578
0.0042531
0
5.770186
0.046076
0.02
0.96
2.1329690
−740.11
−1429.72
−0.000315
0
−0.00708
−0.00611
0.02
1
8.4304106
−207.521
−474.186
−0.015371
0
−0.20741
−0.01250
0.02
1.12
11.304093
−99.3194
−193.856
−0.0655924
0
−1.42005
−0.04391The dynamical systems (4a) and (4b) possess five kinds of solutions [23] as follows.(1)
The stationary state(O), given by
(60)ρ1=ρ2=ρ3=0,is stable forμ<μ2=0  and unstable for μ>μ2.(2)
Stripe patterns(S), given by
(61)ρ1=μg1≠0,ρ2=ρ3=0,are stable forμ>μ3=h2g1/(g2-g1)2, and unstable for μ<μ3.(3)
Hexagon patterns(H0,Hπ) are given by
(62)ρ1=ρ2=ρ3=|h|±h2+4(g1+2g2μ)2(g1+2g2),withφ=0 or π, and exist when
(63)μ>μ1=-h24(g1+2g2).The solutionρ+=|h|+h2+4(g1+2g2μ)/2(g1+2g2) is stable only for
(64)μ<μ4=2g1+g2(g2-g1)2h2,andρ-=(|h|-h2+4(g1+2g2μ))/2(g1+2g2) is always unstable.(4)
The mixed states are given by(65)ρ1=|h|g2-g1,ρ2=ρ3=μ-g1ρ12g1+g2,withg2>g1. They exist whenμ>μ3and are always unstable.


## 4. Spatial Pattern of Systems(4a) and (4b)

In this section, we perform extensive numerical simulations of the spatially extended systems (4a) and (4b) in two-dimensional spaces. All our numerical simulations employ the zero-flux boundary conditions with a system size of 200 × 200. The space step is ΔH=1, and the time step is Δt=0.00001.In Figure3, we show the spatial pattern of prey population at different time. In the parameter set, γ=1.5, ε=0.15, and δ=0.92, we find that μ∈(μ3,μ4), which means that there is coexistence of spotted and stripe patterns. As shown in this figure, our theoretical results are consistent with the numerical results.Spatial pattern of prey population at different time. Parameters set:γ=1.5,  ε=0.15, and  δ=0.92. (a)  t=0; (b)  t=100; (c)  t=200; (d)  t=500; (e)  t=1000; and (f)  t=2000.
(a)
(b)
(c)
(d)
(e)
(f)By settingγ=1.5, ε=0.15, and δ=0.96, one can obtain that μ>μ4. In Figure 4, we show the spatial pattern of prey population whentequals 0, 50, 100, 200, 500, and 1000. At the initial time, the prey population shows patched invasion. As time increases, stripe pattern appears and the structure does not change a lot. While keeping other parameters fixed and increasingδ, we find that stripe pattern will occupy the whole space. However, some stripe patterns connect with each other and cause the emergence of spotted patterns which are shown in Figure 5.Spatial pattern of prey population at different time. Parameters set:γ=1.5,  ε=0.15,  and  δ=0.96. (a)  t=0; (b)  t=50; (c)  t=100; (d)  t=200; (e)  t=500; and (f)  t=1000.
(a)
(b)
(c)
(d)
(e)
(f)Spatial pattern of prey population at different time. Parameters set:γ=1.5,ε=0.15, andδ=1. (a)t=0; (b)t=150; (c)t=300; (d)t=500; (e)t=600; and (f)t=1000.
(a)
(b)
(c)
(d)
(e)
(f)Figure6 shows the evolution of the spatial pattern of prey population at t=0, 100, 300, 500, 1000, and 2000 iterations, with small random perturbation of the stationary solution of the spatially homogeneous systems (4a) and (4b). The corresponding parameters values are γ=1.5, ε=0.15, and δ=1.04. By the amplitude equations, we can conclude that there are spotted patterns of prey population for this parameter set. In this case, one can see that for the systems (4a) and (4b), the random initial distribution leads to the formation of an irregular transient pattern in the domain. After these forms, it grows slightly and spotted patterns emerge. When the time is large enough, the spotted patterns prevail over the two-dimensional space. As time further increases, the pattern structures of the prey population do not undergo any further changes.Spatial pattern of prey population at different time. Parameters set:γ=1.5,  ε=0.15, and  δ=1.04. (a)  t=0; (b)  t=100; (c)  t=300; (d)  t=500; (e)  t=1000; and (f)  t=2000.
(a)
(b)
(c)
(d)
(e)
(f)


## 5. Conclusion and Discussion

Allee effect has been paid much attention due to its strong potential impact on population dynamics [24]. In this paper, we investigated the pattern dynamics of a spatial predator-prey systems with Allee effect. Based on the bifurcation analysis, exact Turing pattern region is obtained. By using amplitude equations, the Turing pattern selection of the predator-prey system is well presented. It is found that the predator-prey systems with Allee effect have rich spatial dynamics by performing a series of numerical simulations.It should be noted that our results were obtained under the assumption that predation is modeled by the bilinear function of the prey and predator densities. However, this function has limitations to describe many realistic phenomena in the biology. By numerical simulations, we find that the system exhibits similar behaviour when the functional response is of other types, such as Holling-II and Holling-III forms.To compare the spatial dynamics for different parameters, we give the spatial patterns of populationuwhen the parameter values are out of the domain of Turing space. For this parameter set, systems (4a) and (4b) have Hopf bifurcation, and spiral waves occupy the whole domain instead of stationary patterns, which is shown in Figure 7. The stability of spiral wave can be done by using the spectrum theory analysis [25, 26]. In the further study, we will use the spectrum theory to show the stability of spiral wave.Spatial pattern of prey population at different time. Parameters set:γ=1.5,  ε=0.15, and  δ=1.2. (a)  t=0; (b)  t=100; (c)  t=200; (d)  t=300; (e)  t=400; and (f)  t=500.
(a)
(b)
(c)
(d)
(e)
(f)In [15], they found that a spatial predator-prey model with Allee effect and linear death rate could increase the system’s complexity and enhance chaos in population dynamics. However, in this paper, we showed that a spatial population model with Allee effect and nonlinear death rate can induce stationary patterns, which is different from the previous results.From a biological point of view, our results show that predator mortality plays an important role in the spatial invasions of populations. More specifically, low predator mortality will induce stationary patterns (cf. Figures3–6), and high predator mortality corresponds to travelling patterns (cf. Figure 7). When the populations exhibit wave distribution in space, the dynamics of populations may be accompanied with chaotic properties [27, 28]. If the chaotic behavior occurs, it may lead to the extinction of the population, or the population may be out of control [29, 30]. In that case, we need to find out the best way to control the chaos or change the chaotic behavior.



---

*Source: 921879-2013-09-24.xml*