# Spatiotemporal Pattern in a Self- and Cross-Diffusive Predation Model with the Allee Effect

Authors: Feng Rao
Journal: Discrete Dynamics in Nature and Society
 (2013)

---

## Abstract

This paper proposes and analyzes a mathematical model for a predator-prey interaction with the Allee effect on prey species and with self- and cross-diffusion. The effect of diffusion which can drive the model with zero-flux boundary conditions to Turing instability is investigated. We present numerical evidence of time evolution of patterns controlled by self- and cross-diffusion in the model and find that the model dynamics exhibits a cross-diffusion controlled formation growth to spotted and striped-like coexisting and spotted pattern replication. Moreover, we discuss the effect of cross-diffusivity on the stability of the nontrivial equilibrium of the model, which depends upon the magnitudes of the self- and cross-diffusion coefficients. The obtained results show that cross-diffusion plays an important role in the pattern formation of the predator-prey model. It is also useful to apply the reaction-diffusion model to reveal the spatial predation in the real world.


---

## Body

## 1. Introduction

It is well known that predator-prey interaction is one of the basic interspecies relations for ecological and social models. The mathematical model of predator-prey type has played a major role in the studies of biological invasion of foreign species, epidemics spreading, and extinction/spread of flame balls in combustion or autocatalytic chemical reaction [1]. The first differential equation model of predator-prey type Lotka-Volterra equation was formulated by Lotka and Volterra in 1920s, when attempts were first made to find ecological laws of nature [2, 3]. In recent years, one of the important predator-prey models is the Holling type II model, originally due to Holling [4, 5], which has been studied in many articles; see, for example, [6–12].In population dynamics, any mechanism that can lead to the positive relationship between population density and the per capita growth rate can be termed an Allee effect [13–15], starting with the pioneer work of Allee [16]. The Allee effect can be caused by a number of sources such as difficulties in finding mates at small densities, reproductive facilitation, predation, environment conditioning, and inbreeding depression. This effect usually saturates or disappears as populations get larger [17, 18]. The Allee effect’s strong potential impact on population dynamics has been attracting much attention recently. Detailed investigations relating to the Allee effect may be found in [19–33].On the other hand, in nature, there is a tendency that the preys would keep away from predators and the escape velocity of the preys may be taken as proportional to the dispersive velocity of the predators. In the same manner, there is a tendency that the predators would get closer to the preys and the chase velocity of predators may be considered to be proportional to the dispersive velocity of the preys [34, 35]. In view of these, cross-diffusion arises, which was proposed first by Kerner [36] and first applied in competitive population system by Shigesada et al. [37] to describe the population pressure by other species. Cross-diffusion expresses the population fluxes of one species due to the presence of the other species. In the realistic situation, self-diffusion implies passive diffusion where the diffusing species moves along its concentration gradient, while cross-diffusion implies countertransport. The cross-diffusion coefficient may be positive or negative. When it is positive, one species moves along the concentration gradient of another species, and when it is negative, one species moves against the concentration gradient of another species [38, 39].Among theoretical and mathematical biologists, there has been considerable interest to investigate the dynamical behavior of a predation system of interacting populations by taking into account the effect of self- and cross-diffusion [40–48]. But in the studies on spatiotemporal dynamics of a Holling type II predator-prey model with an Allee effect on prey, little attention has been paid to study the effect of self- and cross-diffusion. Based on the previous discussions, the main purpose of this paper is to examine out the stability behavior of the Holling-type-II predator-prey model with Allee effect on prey by taking into account self- and cross-diffusion.The organization of this paper is as follows. In Section2, we present a Holling-type-II predation model in the presence of Allee effect with self- and cross-diffusion. In Section 3, we analyze the model and derive the condition of Turing instability with respect to parameters. Furthermore, we use numerical simulations to reveal the emergence of different patterns and the influence of cross-diffusion on the dynamical behavior of the model. Finally, a brief discussion is given in Section 4.


## 2. Model

In this paper, we consider the following model of two-dimensional (i.e., 2D) population dynamics:
(1)N˙=F(N)-f(N,P),P˙=μf(N,P)-MP,
where N and P are the densities of prey and predator at time T>0, respectively. The function F(N) describes prey multiplication, f(N,P) describes predation, and the term MP stands for predator mortality. μ is the food utilization coefficient. We assume that the predator response is of Holling-type-II which is usually parameterized as follows:
(2)f(N,P)=ANPN+B,
where A is the predation rate and B is the half-saturation density.Referring to [18], we assume that Allee dynamics for the prey population, its growth rate can be parameterized as follows:
(3)F(N)=4ω(K-N0)2N(N-N0)(K-N),
where K is the prey carrying capacity, ω is the maximum per capita growth rate, and N0 is the “threshold” (N0<K), so that for N<N0, the growth rate becomes negative. The value of N0 can be considered as a measure of the intensity of the Allee effect: the less the value of N0 is, the less prominent is the Allee effect [21]. In other words, the Allee effect is called “strong” if 0<N0<K (when the growth rate becomes negative for N<N0) and “weak” if -K<N0≤0. For N0≤-K, the Allee effect is absent [49].Based on the Holling-type-II predator-prey model with the Allee effect for the prey, in the following, we investigate the dynamical complexity induced by self- and cross-diffusion in the model. By choosing appropriate scales for the variables of (1), (2), and (3), the number of parameters can be lessened. Considering dimensionless variables as in [21] with the following scaling:
(4)u=NK,v=PμK,t=μAKTB,
we obtain and employ the corresponding diffusive model as follows:
(5)∂u∂t=γu(u-β)(1-u)-uv1+αu+D11∇2u+D12∇2v,∂v∂t=uv1+αu-δv+D21∇2u+D22∇2v.
Here,
(6)α=KB,β=N0K,γ=4ωBKμA(K-N0)2,δ=MBμAK.In model (5), the nonnegative constants D11 and D22 are the self-diffusion coefficients of u and v, respectively, which imply the movement of individuals from a higher to lower concentration region. The cross-diffusion coefficients D12 and D21 express the respective population fluxes of the preys and predators resulting from the presence of the other species, respectively. That is to say, D12 in (5) represents the tendency of the prey u to keep away from its predators v, and D21 represents the tendency of the predator to chase its prey. The cross-diffusion coefficients, D12 and D21, may be positive or negative. Positive cross-diffusion coefficient denotes that one species tends to move in the direction of lower concentration of another species, and negative cross-diffusion expresses the population fluxes of one species in the direction of higher concentration of another species [39, 50]. Note that ∇2=∂2/∂x2+∂2/∂y2 is the usual Laplacian operator in a two-dimensional space.Model (5) is to be analyzed with the following nonzero initial conditions:
(7)u(x,y,t)>0,v(x,y,t)>0,(x,y)∈Ω=[0,Lx]×[0,Ly]
and the following zero-flux boundary conditions:
(8)∂u∂n=∂v∂n=0,(x,y)∈∂Ω.In the previous part,Lx and Ly give the size of the model in the directions of x and y, respectively. n is the outward unit normal vector of the boundary ∂Ω which we will assume to be smooth. Zero-flux boundary conditions imply that no external input is imposed from outside [1, 51].For model (5), without diffusion terms, there are four stationary states: E0=(0,0), E1=(1,0), E2=(β,0), and the coexistence state E*=(u*,v*), where
(9)u*=δ1-αδ,v*=γ(u*-β)(1-u*)(1+αu*).Using the standard linear stability analysis, it can be readily seen that(i)
E0=(0,0) is always a stable node and corresponds to extinction of both prey and predator;(ii)
two uniform steady statesE1=(1,0) and E2=(β,0) are either nodes or saddle points, depending on the relation between the parameters, and correspond to extinction of the predator;(iii)
the coexistence stateE*=(u*,v*) lies in the biologically meaningful domain u>0, v>0 under condition β<u*<1 (for these parameter values E1=(1,0) and E2=(β,0) are saddle points) and it can be either a stable or unstable focus or a node.For more details about model (5) without diffusions, we can refer to reference [29].When we take the parameters asα=1.11, β=0.1, γ=11.25, and δ=0.36 for (5), E0=(0,0) is a stable node, E1=(1,0) and E2=(0.1,0) are saddle points, and the positive equilibrium E*=(0.6,3.75) exists, which is unstable and shown in Figure 1.Figure 1
Phase portrait of model (5). The parameters are taken as α=1.11, β=0.1, γ=11.25, and δ=0.36. The horizontal axis is the prey population u, and the vertical axis is the predator population v. E0=(0,0) is a stable node, E1=(1,0) and E2=(0.1,0) are saddle points, and E*=(0.6,3.75) is unstable.


## 3. Pattern Formation

To study patterns of the self- and cross-diffusion model given by (5), we must consider a spatially homogeneous system.In this section, we perform extensive numerical simulations of the spatially extended model (5) in a two-dimensional space, and the qualitative results are shown here. All our numerical simulations employ the zero-flux boundary conditions (8) with a system size of Lx×Ly discretized through x→(x0,x1,x2,…,xn) and y→(y0,y1,y2,…,yn) with n=200, and the time step of t=1/100.We use the standard five-point approximation for the 2D Laplacian with the zero-flux boundary conditions [52, 53]. That is, the concentrations (ui,jn+1,vi,jn+1) at the time (n+1)t at the mesh position (xi,yj) are given by
(10)ui,jn+1=ui,jn+tD11Δhui,jn+tD12Δhvi,jn+tf(ui,jn,vi,jn),vi,jn+1=vi,jn+tD21Δhui,jn+tD22Δhvi,jn+tg(ui,jn,vi,jn),
with the Laplacian defined by
(11)Δhui,jn=ui+1,jn+ui-1,jn+ui,j+1n+ui,j-1n-4ui,jnh2,
where the space stepsize h=1/3.We are mainly interested in investigating the behavior of model (5) around the interior equilibrium point; so we will put emphasis on the positive equilibrium point E*=(u*,v*). The entire system is initially placed in the stationary state (u*,v*), and the propagation velocity of the initial perturbation is thus of the order of 5×10-4 space units per time unit. After the initial period during which the perturbation spreads, either the system goes into a time-dependent state or goes into a time-independent steady state [54].To study the effect of diffusion on model (5), we consider the linearized form of the model about E*=(u*,v*) as follows:
(12)∂u^∂t=J11u^+J12v^+D11∇2u^+D12∇2v^,∂v^∂t=J21u^+J22v^+D21∇2u^+D22∇2v^,
where u=u*+u^, v=v*+v^ and
(13)J11=-(d((b-a-ab)n2+a2d(2n+d)+(ab+a+b)d2-2bdn))×(bn(n-d)2)-1,J12=-dmn,J21=-(ad+d-n)(ad+bd-bn)bm(n-d),J22=0.
We note that (u^,v^) are small perturbations of (u,v) about the equilibrium point (u*,v*) and the stationary state can become unstable to a nonuniform perturbation
(14)u^(r,t)~exp(λt+ik→·r),v^(r,t)~exp(λt+ik→·r),r=(x,y).
Any solution of model (12) can be expanded into a Fourier series so that
(15)u^(r,t)=∑i,j=0∞‍uij(r,t)=∑i,j=0∞‍αij(t)sin(k→r),v^(r,t)=∑i,j=0∞‍vij(r,t)=∑i,j=0∞‍βij(t)sin(k→r),
where 0<x<Lx, 0<y<Ly. Note that k→=(ki,kj), ki=iπ/Lx, and kj=jπ/Ly are the corresponding wavenumbers.Substitutinguij and vij into (12), we obtain
(16)dαijdt=(J11-D11k2)αij+(J12-D12k2)βij,dβijdt=(J21-D21k2)αij-D22k2βij,
where k2=ki2+kj2. The characteristic equation of the linearized system is given by
(17)λ2-tr(𝒥)λ+det(𝒥)=0,
where
(18)𝒥=(J11-D11k2J12-D12k2J21-D21k2-D22k2),tr(𝒥)=J11-(D11+D22)k2,det(𝒥)=det(D)k4+(D12J21+D21J12-D22J11)k2-J12J21.
Here, the diffusive matrix is
(19)D=(D11D12D21D22),
and the determinant of D is det(D)=D11D22-D12D21. Then the roots of (17) yield the dispersion relation
(20)λ1,2=12(tr(𝒥)±(tr(𝒥))2-4det(𝒥)).In order to have diffusion-induced instability (i.e., Turing instability), which means that it is an asymptotically stable equilibrium of the diffusion-less model (5) but is unstable with respect to solutions of the spatial model (5), at least one of the conditions tr(𝒥)<0 and det(𝒥)>0 must be violated by Routh-Hurwitz criterion. It is evident that the first condition tr(𝒥)<0 is not violated when the condition J11<0 is met. Hence, in order to obtain diffusion instability, we need det(𝒥)<0. Thus, the condition is given by
(21)G(k2)=det(D)k4+(D12J21+D21J12-D22J11)k2-J12J21<0.
If the minimum of G(k2) is negative, then a range of modes k-2<k2<k+2 will grow, where
(22)k-2=12det(D)×[-2D22J11(D21J12+D12J21))1/2D22J11-D21J12-D12J21-(D222J112+(D21J12+D12J21)2+4J12J21det(D)-2D22J11(D21J12+D12J21)D222J112)1/2],k+2=12det(D)×[-2D22J11(D21J12+D12J21))1/2D22J11-D21J12-D12J21+(D222J112+(D21J12+D12J21)2+4J12J21det(D)-2D22J11(D21J12+D12J21)D222J112)1/2].In addition, when(23)k2=kc2=D22J11-D12J21-D21J122det(D)>0,
we obtain
(24)Gmin=-J12J21-(D22J11-D12J21-D21J12)24det(D).
Then we have
(25)G(k2)<0⟺k-2<k2<k+2.
If both k-2 and k+2 exist and have positive values, they limit the range of instability for a locally stable equilibrium (see Figure 2).Figure 2
Model (5) with wavenumbers lying between the zeros of G(k2), k-2, and k+2 growing in the Turing instability. Parameters are taken as α=1.11, β=0.1, γ=11.25, δ=0.36, D11=0.2, D22=2, D12=-0.5, and D21=-0.4; hence, k-2=0.427433403, k+2=2.807586678, and Gmin=-0.283261256.Summarizing the above discussions, we can obtain the following theorem.Theorem 1.
The criterion for Turing instability for model (5) satisfies the following condition:
(26)D22J11-D12J21-D21J12>2-J12J21det(D)>0.In Figure3, we plot the dispersion relation corresponding to several values of parameter D12 while keeping the others fixed as
(27)α=1.11,β=0.1,γ=11.25,δ=0.36,D11=0.2,D22=2,D21=-0.4.Figure 3
Variation of dispersion relation of model (5) for different D12: (A) D12=-1.6<D12T, the Turing instability occurs; (B) D12=D12T=-0.31, the critical Turing value; and (C) D12=0.4>D12T, the Turing instability decays. Other parameters values are fixed as (27) in the text.Curve (B) in Figure3 corresponds to the critical Turing value D12T=-0.31. When the parameter D12=-1.6 below the D12T (i.e., curve (A) in Figure 3), the Turing instability occurs, while D12=0.4 above the D12T (i.e., curve (C) in Figure 3), the Turing instability decays. That is to say, when D12>D12T, the steady state is the only stable solution of model (5). The spatial patterns are generated when D12 passes through the critical Turing bifurcation point D12T.Figure4 shows the spatial patterns of prey u and predator v and the corresponding phase portraits of model (5). Figures 4(a1) and 4(a2) show stationary spotted and striped-like patterns emerge mixed in the distribution of prey u and predator v population density, respectively. In Figures 4(a1) and 4(a2), when D12=-0.365<D12T=-0.31, starting with a homogeneous state (0.6,3.75), the random perturbations make the formations end with spotted and striped-like coexistence in the spatial domain at t=10000 (i.e., 1000000 iterations), and the mixed patterns prevail over the whole domain at last. In Figure 4(a3), one can see that as t→∞, u and v coordinates fix to (0.745,3.428) which is the steady state.Figure 4
Dynamical behavior of model (5). D12=-0.365 (the first row); D12=-0.5 (the second row); and other parameters are fixed as (27) in the text. (a1), (b1) Two categories of Turing patterns of prey u at t=10000; (a2), (b2) two categories of Turing patterns of predator v at t=10000; and (a3), (b3) phase portraits.Figures4(b1) and 4(b2) show spotted pattern formations of prey u and predator v, respectively. In Figures 4(b1) and 4(b2), D12=-0.5<D12T=-0.31, starting with a homogeneous state (0.6,3.75), the random perturbations make the formations end with a time-independent spotted-patterns at t=10000 (i.e., 1000000 iterations). From Figures 4(b1) and 4(b2), we can see that the spotted patterns prevail over the whole domain and the dynamics of model (5) do not undergo any further changes. In Figure 4(b3), which is a phase portrait, the steady state (0.725,3.061) exists and is asymptotically stable. From the figure, we can see that (u,v)→(0.725,3.061) with time.In Figure4, two different snapshots during the temporal evolution of model (5) are presented in a two-dimensional space under different D12 values. These results indicate that the effect of cross-diffusion D12 for pattern formation is tremendous, and model (5) is conditionally stable at unique positive equilibrium.In order to investigate quantitatively the evolution of model (5), we depict the typical emergence of Turing pattern in the spacetime plots of Figure 5. The method of spacetime plots is to let y be a constant, choose the line y=100 from each pattern snapshot, and pile these lines in time order [55–59]. Figures 5(a) and 5(b) show the evolution process of prey u throughout time t and space x with different D12, respectively, where time increases from bottom to top, and the horizontal axis represents the spatial location. From Figures 5(a) and 5(b), the single parameter of model (5) namely D12, which is the diffusion coefficient of the prey, can lead to dramatic changes in the qualitative dynamics of solutions. From a biological perspective, the cross term plays a constructive role in the pattern formation of the prey, that is, changes the Turing pattern into regular spatial pattern with time.Spacetime plots of preyu of model (5). The time interval shown is 1500. (a) D12=-0.365 and (b) D12=-0.5. Other parameters are fixed as (27) in the text.
(a)
(b)According to the works of many researchers [53, 60–62], the choice of initial conditions can affect the spatiotemporal dynamics of a reaction-diffusion model in ecosystems. In Figure 6, when D12=-0.5 and other parameters are fixed as (27), we choose the special initial condition of model (5) to investigate the evolution of the spatial pattern of prey u at a different evolution time. The initial condition is introduced as the following form:
(28)u(x,y,0)=u*,v(x,y,0)=v*,90<x,y<110,
which is a square. Initially, the entire system is placed in the stationary state (u*,v*)=(0.6,3.75), and the propagation velocity of the initial random perturbations is thus of the order of 5×10-4 space units per time unit. Also the system is then integrated for 105 time steps, and some images are saved.The process of pattern formation of preyu in model (5) with D12=-0.5; other parameters are fixed as (27) in the text. Time steps: (a) t=10; (b) t=105; (c) t=151; (d) t=179; (e) t=300; and (f) t=1000.
(a)
(b)
(c)
(d)
(e)
(f)From Figure6(a), we can see that a target wave pattern emerges after perturbation of the stationary solutions u* and v* of model (5). Some iterations later, it breaks towards the core, and a spot pattern (interior) with target wave pattern (exterior) occurs (see Figure 6(b)). Then one can see that after the symmetrical target-spotted mixture patterns (see Figures 6(c) and 6(d)), it grows slightly and the spots increase with time (see Figures 6(e) and 6(f)). From Figure 6, we observe that when the iteration time is large enough, the time-independent striped-like and target-spotted mixed patterns take a long time to settle down and stay stable. It finally prevails over the whole domain at t=1000.Figure6 illustrates the formation of the spatial patterns behavior of prey u in two dimensions with self- and cross-diffusion, where the initial condition (at t=0) adds random and nonuniform small perturbations to the equilibrium values.


## 4. Conclusions and Remarks

In this paper, we investigate the pattern formation of a Holling-type-II predator-prey interacting model with an Allee effect in prey population and with self- and cross-diffusion under the zero-flux boundary conditions. We show that cross-diffusion in both predator and prey species tends to stabilize the predation model. That is to say, the unstable equilibrium(u*,v*) of model (5) becomes stable under certain conditions. Moreover, the numerical results show that model (5) dynamics exhibits a cross-diffusion controlled formation growth not only to spotted-striped, but also to spotted pattern replication.By a series of analysis and numerical simulations, we find that the evolutionary process depends on several parameters. In this paper, we fix the values ofα, β, γ, δ, D11, D22, and D21 as (27) and perform a large number of computer simulations corresponding to different D12. The results show that model (5) has rich spatiotemporal patterns and different patterns appear on account of the change of cross-diffusion coefficient D12. Biological speaking, cross-diffusion may have effect on the distribution of the species. It is also shown that if the equilibrium state (u*,v*) of model (5) with no diffusion is unstable, then the corresponding uniform steady state of the model with self- and cross-diffusion can be made stable by changing cross-diffusion coefficients. Therefore, we hope that the results presented here will be useful for studying the dynamical complexity of ecosystems or physical systems.



---

*Source: 681641-2013-07-01.xml*