# Pattern Formation in Predator-Prey Model with Delay and Cross Diffusion

Authors: Xinze Lian; Shuling Yan; Hailing Wang
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We consider the effect of time delay and cross diffusion on the dynamics of a modified Leslie-Gower predator-prey model incorporating a prey refuge. Based on the stability analysis, we demonstrate that delayed feedback may generate Hopf and Turing instability under some conditions, resulting in spatial patterns. One of the most interesting findings is that the model exhibits complex pattern replication: the model dynamics exhibits a delay and diffusion controlled formation growth not only to spots, stripes, and holes, but also to spiral pattern self-replication. The results indicate that time delay and cross diffusion play important roles in pattern formation.


---

## Body

## 1. Introduction

The spatial component of ecological interaction has been identified as an important factor in how ecological communities are shaped. Understanding the role of space is challenging both theoretics and empirism [1]. Empirical evidence suggests that the spatial scale and structure of environment can influence population interactions [2]. And pattern formation in nonlinear complex systems is one of the central problems of the natural, social, and technological sciences [3–12]. The classical approach to the origin of spatial structures was developed by Turing [13]. Understanding of patterns and mechanisms of species spatial dispersal is an issue of significant current interest in conservation biology and ecology. It arises from many ecological applications; particularly, it plays a major role in connection to biological invasion and disease spreading [14, 15].In recent years, one of the important predator-prey models is Leslie-Gower model. A modified version of reaction-diffusion Leslie-Gower predator-prey model with Holling-type II functional response takes the form [16, 17]
(1)∂u∂t=u(r-au)-c1uvb1+u+d1∇2u,∂v∂t=v(d-c2vb2+u)+d2∇2v,
where u and v are prey and predator densities, respectively. r, a, d, bi, and ci(i=1,2) are positive constants. r is the growth rate of preys, d is the growth rate of predators, a is the strength of competition among individuals of species u, b1 (resp., b2) measures the extent to which environment provides protection to prey u (resp., to the predator v), and c1 is the maximum value of the per capita reduction of u due to v. c2 has a similar meaning to c1. d1 and d2 are the diffusion coefficients that imply the movement of individuals from a higher to a lower concentration region. ∇2=(∂2/∂x2)+(∂2/∂y2) is the usual Laplacian operator in two-dimensional space.For model (1), Camara and Aziz-Alaoui [18–20] paid more attention to pattern formation in the spatial domain and observed the labyrinth, chaos, and spiral wave patterns. In addition, Guan et al. [21] studied model (1) incorporating a refuge protecting mu of the prey in the case of c1=c2≠0 and b1=b2≠0, where m∈[0,1) is a constant. This leaves (1-m)u of the prey available to the predator, and model (1) becomes
(2)∂u∂t=u(r-au)-(1-m)uvb+(1-m)u+d11∇2u≜f(u,v)+d11∇2u,∂v∂t=v(d-cvb+(1-m)u)+d22∇2v≜g(u,v)+d22∇2v≜v(d-h(u,v))+d22∇2v.
In [21], the authors proved the global stability of the positive equilibria, determined the Turing space in the spatial domain, and performed a series of numerical simulations and found that model (2) exhibits complex Turing pattern replication: spots, stripes, and holes patterns.On the other hand, a time delay occurs naturally in just about every interaction of the real world. The original motivations for studying delay differential equations mainly come from their applications in feedback control theory [22–29]. And time delay also is an ubiquitous phenomenon in biological systems and has an important role in affecting population dynamics. In a system described by ordinary differential equations, time delay could change qualitatively the nature of the equilibrium from a stable equilibrium to an unstable one and thus induces bifurcation [30]. For a predator-prey system, time delay was often realized as feedback. The effect of this kind of delay on the dynamical behavior of populations has been studied by a number of papers [31–33]. Because of the importance of time delay to control dynamics of a system, a number of works have been done on effect of delay on the spatial-temporal system [34–44]. In our previous work [43], we studied the spatiotemporal complexity of a Leslie-Gower model incorporating a prey refuge with time delay effect:
(3)∂u∂t=u(r-au)-(1-m)uvb+(1-m)u+d11∇2u≜f(u,v)+d11∇2u,∂v∂t=v(d-cv(t-τ)b+(1-m)u(t-τ))+d22∇2v≜v{d-h[u(t-τ),v(t-τ)]}+d22∇2v,
where τ>0 is a constant delay due to gestation and u(x,y,t-τ) and v(x,y,t-τ) are densities of prey and predator, respectively, at moment t-τ and position (x,y). For simplicity, set u(x,y,t-τ)∶=u(t-τ),v(x,y,t-τ)∶=v(t-τ). In [43], we analyzed the effects of delay and diffusion on the spatiotemporal dynamics of model (3) and gave the pattern formation via numerical simulation.Furthermore, 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 [45–49]. 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. Keeping these in view, cross-diffusion arises, which was proposed first by Kerner [50] and first applied in competitive population system by Shigesada et al. [51]. And there comes a question: how does cross-diffusion affect pattern formation in a delayed Leslie-Gower predator-prey model incorporating a prey refuge?In this paper, we mainly focus on the effect of cross-diffusion on the pattern formation in the following Leslie-Gower predator-prey model incorporating a prey refuge:(4)∂u∂t=f(u,v)+d11∇2u+d12∇2v,∂v∂t=v{d-h[u(t-τ),v(t-τ)]}+d21∇2u+d22∇2v,
where 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. Other parameters are the same definitions as model (3). And we call
(5)D=(d11d12d21d22)
the diffusive matrix and assume that d11>0, d22>0, and det(D)=d11d22-d12d21>0.Following [44], we explain the meaning and units of each variable and constant in model (4).(i)
The preyu and predator v are usually measured in mg of dry weight per liter [mg·dw/l].(ii)
Timet and time delay τ are measured in days [d].(iii)
The length ofx,y∈(0,L) is measured in meter [m].(iv)
r is the growth rate of prey u and d is the growth rate of predator v, measured in [d-1].(v)
c is the maximum value of the per capita reduction of v due to u, which is not available in abundance and is measured in [d-1].(vi)
a describes the strength of competition among individuals of species u and is measured in [(mg·dw/l)-1d-1].(vii)
b describes the extent to which environment provides protection to prey u and predator v and is usually measured in mg of dry weight per liter [mg·dw/l].(viii)
dij>0(i,j=1,2) are the diffusion coefficients corresponding to the prey and predator, respectively, and are measured in [m2d-1].Model (4) is to be analyzed under the following nonzero initial conditions:
(6)u(x,y,t0)>0,v(x,y,t0)>0,(x,y)∈Ω=(0,L)×(0,L),t0∈[-τ,0],
and zero-flux boundary conditions:
(7)∂u∂n|∂Ω=∂v∂n|∂Ω=0.
In the above, L denotes the size of the system in square domain, 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 zero-flux conditions imply no external input [5].This paper is organized as follows. In the next section, we mainly present the Hopf and Turing bifurcation analysis. In Section3, we perform a series of numerical simulation to show evolution process of prey u. And in the last section, we give some concluding remarks.


## 2. Bifurcation Analysis

In this section, we will focus on the conditions of Hopf and Turing bifurcation. Following [40, 43], if time delay τ is small, one can expand u(x,y,t-τ) and v(x,y,t-τ) as
(8)u(x,y,t-τ)=u(x,y,t)-τ∂u(x,y,t)∂t,v(x,y,t-τ)=v(x,y,t)-τ∂v(x,y,t)∂t.
Substituting (8) into (4), we obtain
(9)∂u∂t=f(u,v)+d11∇2u+d12∇2v,∂v∂t=v[d-h(u(x,y,t)-τ∂u(x,y,t)∂t,v(x,y,t)-τ∂v(x,y,t)∂t)]+d21∇2u+d22∇2v.
Expanding (9) in Taylor series and neglecting the higher-order nonlinearities, then (9) becomes
(10)∂u∂t=f(u,v)+d11∇2u+d12∇2v,∂v∂t=v(d-h(u,v)+τhu(u,v)∂u(x,y,t)∂t+τhv(u,v)∂v(x,y,t)∂t)+d21∇2u+d22∇2v,
where hu(u,v)=(∂h(u,v)/∂u) and hv(u,v)=(∂h(u,v)/∂v).From (10), we finally obtain that
(11)∂u∂t=f(u,v)+d11∇2u+d12∇2v,(12)∂v∂t=11-vτhv(u,v)×([vτhu(u,v)d12+d22]∇2g(u,v)+vτhu(u,v)f(u,v)+[vτhu(u,v)d11+d21]∇2u+[vτhu(u,v)d12+d22]∇2v).
Equation (11) can be used to analyze the bifurcation of model (4) when τ is small.It is easy to obtain that model (11) has three boundary equilibria E1=(0,0), E2=((a/r),0), and E3=(0,(bd/r)) and a unique positive equilibrium E*=(u*,v*)=(((dm-d+cr)/ac),(((cab+(dm-d+cr)(1-m))d)/ac2)).Remark 1.
In the case without time delay and cross-diffusion, that is,τ=0 and d12=d21=0, model (11) becomes (2). It is easy to see that E1,E2,E3, and E* are also equilibria of model (2). The stabilities of these equilibria of model (2) can be seen in [21].In the following, we are mainly concerned with the properties of Hopf and Turing bifurcations for model (11).Let0=μ0<μ1<μ2<⋯ be the eigenvalues of the operator −Δ on Ω with the zero-flux boundary conditions. And setting
(13)X={(u,v)∈[H2(Ω)]2∣∂nu=∂nv=0on∂Ω},E(μ)={ϕ∣-Δϕ=μϕinΩ,∂νϕ=0on∂Ω},withμ∈ℝ1,{ϕij∣j=1,…,dimE(μi)} will be an orthonormal basis of E(μi), and Xij={cϕij∣c∈ℝ2}; then X=⨁i=1∞Xi, where Xi=⨁j=1dimE(μi)Xij.The linearization of model (11) around E*=(u*,v*) can be expressed by
(14)∂∂t(uv)=JE*(uv)+D~∇2(uv),
with domain X. Here,
(15)JE*=(fufv𝒳(gu+v*τhufu)𝒳(gv+v*τhufv))≜(J11J12J21J22),D~=(d11d12𝒳(v*τhud11+d21)𝒳(v*τhud12+d22)),
with fu=(∂f/∂u)∣(u*,v*), fv=(∂f/∂v)∣(u*,v*), gu=(∂g/∂u)∣(u*,v*), gv=(∂g/∂v)∣(u*,v*), hu=(∂h/∂u)∣(u*,v*), hv=(∂h/∂v)∣(u*,v*), and (1/(1-vτhv)|(u*,v*))=(1/(1-τd))≜𝒳. Since τ is small, we only consider this case of τ<(1/d) (i.e. 𝒳>0) in this paper.For eachj(j=0,1,2,…), Xj is invariant under the operator D~∇2+LE*, and λ is an eigenvalue of D∇2+LE on Xj if and only if λ is an eigenvalue of the matrix Mj=-D~μj+JE* for some j≥0.So the stability is reduced to consider the characteristic equation(16)det(λI-Mj)=λ2-trace(Mj)λ+det(Mj)=0,
where
(17)trace(Mj)=trace(D~)μj-trace(JE*),det(Mj)=𝒳det(D)μj2-𝒳(d11gv-d12gu-d21fv+d22fu)μj+det(JE*).


### 2.1. Hopf Bifurcation

In this subsection, we analyze the properties of Hopf bifurcation for model (11). For Hopf bifurcation to occur at the positive equilibrium E*, the matrix Mj=-D~μj+JE* must have an eigenvalues on the imaginary axis [52]; that is, trace(Mj)≡0. Therefore, the only possible critical values of τ are τ(i) such that, for i≥0,
(18)trace(D~)μi-trace(JE*)=trace(D~)μi-fu+(gv+v*τ(i)hufv)𝒳=0.
At τ=τ(i), trace(Mi)=trace(D~(μi-μj)), and therefore, if i≥1, for all 0≤j<i, trace(Mi)>0 and the matrix Mj=-D~μj+JE* has at least 2i eigenvalues with positive real parts. Therefore, the only value of τ at which Hopf bifurcation hypotheses may be satisfied is
(19)τc=τ(0)∶=fu+gvv*(fuhv-hufv).
Near τc, the complex conjugate pair κ(τ)±iω(τ) is given by
(20)κ(τ)=12trace(JE*),ω2(τ)=det(JE*)-κ2(τ).Sinceκ(τc)=0 and det(JE*)=(((dm+rc-d)d)/c)𝒳>0, there exists an interval I containing τc such that(i)
ω2(τ)≠0 for all τ∈I;(ii)
trace(Mj)<0 and det(Mj)>0 with d11gv-d12gu-d21fv+d22fu<0 uniformly for j≥1 and τ∈I;(iii)
(d/dτ)κ(τ)∣τ=τc=v*fvhu+(dgv/(1-dτc)2)>0.Thus, we have the following theorem.Theorem 2.
Ifd11gv-d12gu-d21fv+d22fu<0, then model (11) undergoes a Hopf bifurcation that occurs around the positive equilibrium E* at τ=τc.


### 2.2. Turing Instability

In this subsection, we will state the Turing instability for the positive equilibriumE* of model (11).Mathematically speaking, the positive equilibriumE* is Turing instability, which was emphasized by Turing in his pioneering work in 1952 [13], which means that it is asymptotically stable in the absence of diffusion in model (11) but unstable with respect to solutions in model (11). Hence, Turing instability occurs when either the condition trace(Mj)<0 or the condition det(Mj)>0 is violated.Since the positive equilibriumE* is stable without diffusion which means that trace(JE*)<0 and det(JE*)>0 hold, then trace(Mj)<0 is always true. Therefore, for the emergency of the diffusion-driven instability in model (11), it is needed that det(Mj)<0 for some μj>0. A necessary condition is
(21)d11gv-d12gu-d21fv+d22fu>0;
otherwise det(Mj)>0 for all j>0 since 𝒳det(D)>0 and det(JE*)>0. And we notice that det(Mj) achieves its minimum:
(22)minj∈R+det(Mj)=(((4det(D~))-14𝒳det(D)det(JE*)-𝒳(d11gv-d12gu-d21fv+d22fu)2)×(4det(D~))-1)
at the critical value μc2>0 where
(23)μc2=d11gv-d12gu-d21fv+d22fu2det(D).
Summarizing the above calculation, we conclude the following.Theorem 3.
If(24)τ<min{1d,τc},4𝒳det(D)det(JE*)-𝒳(d11gv-d12gu-d21fv+d22fu)2<0,d11gv-d12gu-d21fv+d22fu>0,
then the positive equilibrium E* of model (11) is Turing unstable.In Figure1, we show the bifurcation diagram for model (11). The Turing instability curve and Hopf bifurcation curve separate the d-d21 parameter space into four domains. The domain on the left of the Hopf bifurcation curve is stable; the domain above Turing bifurcation curve is unstable. Hence, among these domains, only the domain I satisfies conditions of Theorem 3, and we call domain I as Turing space, where the Turing instability occurs and the Turing patterns may be undergone. Domain IV gives birth to spiral patterns.Figure 1
Bifurcation diagram for model (11) in the d-d21 parameter space. The red line is Turing bifurcation curve, and the blue line denotes Hopf bifurcation. Domain I is the region of pure Turing instability. In domain II, Hopf-Turing instability occurs. And domain III is the stable region. Domain IV is the region of pure Hopf instability.


## 2.1. Hopf Bifurcation

In this subsection, we analyze the properties of Hopf bifurcation for model (11). For Hopf bifurcation to occur at the positive equilibrium E*, the matrix Mj=-D~μj+JE* must have an eigenvalues on the imaginary axis [52]; that is, trace(Mj)≡0. Therefore, the only possible critical values of τ are τ(i) such that, for i≥0,
(18)trace(D~)μi-trace(JE*)=trace(D~)μi-fu+(gv+v*τ(i)hufv)𝒳=0.
At τ=τ(i), trace(Mi)=trace(D~(μi-μj)), and therefore, if i≥1, for all 0≤j<i, trace(Mi)>0 and the matrix Mj=-D~μj+JE* has at least 2i eigenvalues with positive real parts. Therefore, the only value of τ at which Hopf bifurcation hypotheses may be satisfied is
(19)τc=τ(0)∶=fu+gvv*(fuhv-hufv).
Near τc, the complex conjugate pair κ(τ)±iω(τ) is given by
(20)κ(τ)=12trace(JE*),ω2(τ)=det(JE*)-κ2(τ).Sinceκ(τc)=0 and det(JE*)=(((dm+rc-d)d)/c)𝒳>0, there exists an interval I containing τc such that(i)
ω2(τ)≠0 for all τ∈I;(ii)
trace(Mj)<0 and det(Mj)>0 with d11gv-d12gu-d21fv+d22fu<0 uniformly for j≥1 and τ∈I;(iii)
(d/dτ)κ(τ)∣τ=τc=v*fvhu+(dgv/(1-dτc)2)>0.Thus, we have the following theorem.Theorem 2.
Ifd11gv-d12gu-d21fv+d22fu<0, then model (11) undergoes a Hopf bifurcation that occurs around the positive equilibrium E* at τ=τc.


## 2.2. Turing Instability

In this subsection, we will state the Turing instability for the positive equilibriumE* of model (11).Mathematically speaking, the positive equilibriumE* is Turing instability, which was emphasized by Turing in his pioneering work in 1952 [13], which means that it is asymptotically stable in the absence of diffusion in model (11) but unstable with respect to solutions in model (11). Hence, Turing instability occurs when either the condition trace(Mj)<0 or the condition det(Mj)>0 is violated.Since the positive equilibriumE* is stable without diffusion which means that trace(JE*)<0 and det(JE*)>0 hold, then trace(Mj)<0 is always true. Therefore, for the emergency of the diffusion-driven instability in model (11), it is needed that det(Mj)<0 for some μj>0. A necessary condition is
(21)d11gv-d12gu-d21fv+d22fu>0;
otherwise det(Mj)>0 for all j>0 since 𝒳det(D)>0 and det(JE*)>0. And we notice that det(Mj) achieves its minimum:
(22)minj∈R+det(Mj)=(((4det(D~))-14𝒳det(D)det(JE*)-𝒳(d11gv-d12gu-d21fv+d22fu)2)×(4det(D~))-1)
at the critical value μc2>0 where
(23)μc2=d11gv-d12gu-d21fv+d22fu2det(D).
Summarizing the above calculation, we conclude the following.Theorem 3.
If(24)τ<min{1d,τc},4𝒳det(D)det(JE*)-𝒳(d11gv-d12gu-d21fv+d22fu)2<0,d11gv-d12gu-d21fv+d22fu>0,
then the positive equilibrium E* of model (11) is Turing unstable.In Figure1, we show the bifurcation diagram for model (11). The Turing instability curve and Hopf bifurcation curve separate the d-d21 parameter space into four domains. The domain on the left of the Hopf bifurcation curve is stable; the domain above Turing bifurcation curve is unstable. Hence, among these domains, only the domain I satisfies conditions of Theorem 3, and we call domain I as Turing space, where the Turing instability occurs and the Turing patterns may be undergone. Domain IV gives birth to spiral patterns.Figure 1
Bifurcation diagram for model (11) in the d-d21 parameter space. The red line is Turing bifurcation curve, and the blue line denotes Hopf bifurcation. Domain I is the region of pure Turing instability. In domain II, Hopf-Turing instability occurs. And domain III is the stable region. Domain IV is the region of pure Hopf instability.


## 3. Pattern Formation via Numerical Simulations

In this section, we perform extensive numerical simulations of the spatially extended model (4) in two-dimension spaces, and the qualitative results are shown here. Our numerical simulations employ a system size of 100×100 space units. And some parameters are taken as
(25)b=0.22,c=0.4,m=0.17,d12=0.05,τ=0.46.
The numerical integration of model (4) was performed by means of forward Euler integration, with a time step of Δt=0.001 and a space step of Δx=Δy=(1/4), using the standard five-point approximation for the 2D Laplacian with the zero-flux boundary conditions [53].


### 3.1. Turing Pattern Formation in Domain I

In this subsection, we first consider the pattern formation when parameters are located in domain I of Figure1, which is pure Turing unstable domain.Initially, the entire system is placed in the steady state (u*,v*), and the propagation velocity of the initial perturbation is thus on the order of 5×10-4 space units per time unit. The initial data of prey u is a 101×101 matrix, such as
(26)u(x,y,t0)=(u*u*⋯u*u*u*⋯u*⋮⋮⋯⋮u*u*⋯u*)+(rand(101,101)-0.5)×10-3,t0∈[-τ,0],
where rand(101,101) is MATLAB procedure, which returns a 101×101 matrix containing pseudorandom values drawn from the standard uniform distribution on the open interval (0,1). And the system is then integrated for 2,000,000 time steps, and the last images are saved. After the initial period during which the perturbation spreads, the system goes into an essentially steady state (time independent).In the numerical simulations, different types of dynamics are observed, and it is found that the distributions of predator and prey are always of the same type. Consequently, we can restrict our analysis of pattern formation to one distribution. In the following, we show the distribution of preyu.In Figure2, we show three typical Turing patterns and their three-dimension pattern for model (4) at t=2000. In the case of (d,d21,d22)=(0.08,0.89,0.37), the equilibrium is (u*,v*)=(2.8697,0.5203) and the holes patterns are isolated zones with low prey density (cf. Figure 2(a)). When increasing d to 0.4, (d,d21,d22)=(0.47,0.89,0.37), (u*,v*)=(0.9877,1.2218), the stripes patterns appear. While with (d,d21,d22)=(0.61,0.39,0.75), (u*,v*)=(0.3122,0.7306), spots patterns are isolated zones with high prey density (cf. Figure 2(c)). And Figures 2(d)-2(e) are three-dimension patterns corresponding to Figures 2(a)–2(c).Stationary Turing patterns ofu in model (4) at time t=2000 with parameters a=0.43, b=0.22, c=0.4, r=1.4, m=0.17, d11=0.4, d12=0.05, d22=0.37, τ=0.46, and (a) (d) (d,d21,d22)=(0.08,0.89,0.37), (b) (e) (d,d21,d22)=(0.4,0.89,0.37), and (c) (f) (d,d21,d22)=(0.61,0.39,0.75), and the other parameters are fixed as (25). Times: t=2000.
(a)
(b)
(c)
(d)
(e)
(f)On the other hand, it is easy to check that holes, stripes, and spots also appear in model (4) without time delay for τ=0 and other the parameters are taken the same way as in Figures 2(a)–2(c), respectively. However, in parameters of Figure 2, trace(Mj)<0 and det(Mj)>0 for all j≥0 if the cross-diffusion coefficients d12=0 and d21=0. That is to say, model (4) is stable in absence of cross-diffusion, and this instability is induced by cross-diffusion.


### 3.2. Chaotic Pattern Formation Domain II

In this subsection, we will focus on the pattern formation when parameters are located in domain II of Figure1, the domain of pure Hopf-Turing instability.Figure3 shows the evolution process of the spatial pattern of u at time t=0,100,1000,2000 with parameter set (25) and a=0.1, b=0.71, r=1.54, d=0.47, d11=0.2, d21=0.44, d22=0.88, and  τ=0.46. The pattern takes a long time to settle down; starting with a homogeneous state E*=(56.4750,55.3357) (cf. Figure 3(a)), the random initial distribution leads to the formation of a strongly irregular transient pattern in the domain. After the irregular pattern form (cf. Figure 3(b)), it grows slightly and “jumps” alternately for a certain time, and finally the chaos patterns prevail over the whole domain (cf. Figures 3(c)-3(d)), which is time-dependent. For the sake of learning the spatiotemporal dynamics of model (4) further, we illustrate the time-series plots and phase portraits in Figures 3(e)-3(f). From the time-series plots in Figure 3(e), one can see that u and v strongly oscillates over time t. And in Figure 3(f), there exhibits a “local” phase plane of the system invaded by the irregular spatiotemporal oscillations; the trajectory now fills nearly the whole domain inside the quasi-limit cycle. This regime of the system dynamics corresponds to spatiotemporal chaos. And the spatial symmetry of model (4) is broken, and the patterns are oscillatory in space.Chaotic patterns ofu in model (4) with parameters a=0.01, b=0.22, c=0.4, r=1.54, m=0.17, d11=0.2, d12=0.05, d21=0.44, d22=0.88, and  τ=0.46. Times: (a) t=0, (b) t=100, (c) t=1000, (d) t=2000, (e) time-series plots, and (f) phase portraits.
(a)
(b)
(c)
(d)
(e)
(f)


### 3.3. Spiral Wave Pattern Formation Domain IV

In this subsection, we will focus on the pattern formation when parameters are located in domain IV of Figure1, the domain of pure Hopf instability.In order to make any influence of the corners of the domain more visible, and thanks to the insightful work of Medvinsky et al. [4] and Upadhyay et al. [54], we have studied the spiral wave pattern for the following initial condition (27). In this case, we employ Δx=Δy=1 and the system size is 200×200.In order to make any influence of the corner of the domain more visible, following [43, 55], the initial conditions are deliberately chosen to be unsymmetrical as
(27)u(x,y,t0)=u*-ε1(x-100),v(x,y,t0)=v*-ε2(y-100),t0∈[-τ,0],
where ε1=2×10-4 and ε2=3×10-4. For more clarity, we put the initial data of prey u as follow
(28)u(x,y,t0)=(u*+ε1(100-0)u*+ε1(100-0)⋯u*+ε1(100-0)u*+ε1(100-1)u*+ε1(100-1)⋯u*+ε1(100-1)⋮⋮⋯⋮u*+ε1(100-200)u*+ε1(100-200)⋯u*+ε1(100-200)),t0∈[-τ,0].Figure4 shows the evolution process of the spatial pattern of prey u at t=0,500,1500,3000 with (d,d21,d22)=(0.5,0.15,0.37); the other parameters are fixed as in (25). Figure 4(a) shows that, for the model (4) with initial conditions (27), the formation of the irregular patchy structure can be preceded by the evolution of a regular spiral spatial pattern. Note that the appearance of the spirals is not induced by the initial conditions. The center of each spiral is situated in a critical point (xcr,ycr)=(100,100), where u(xcr,ycr)=u*,v(xcr,ycr)=v*. After the spirals form (Figure 4(b)), they grow slightly for a certain time and their spatial structure becomes more distinct (Figures 4(c) and 4(d)). Obviously, spiral pattern arises from Hopf instability.The process of formation of spiral pattern ofu in model (4) with parameters a=0.43, b=0.22, c=0.4, r=1.4, m=0.17, d=0.5, d11=0.4, d12=0.05, d21=0.15, d22=0.37, and τ=0.46. Times: (a) t=0, (b) t=500, (c) t=1500, and (d) t=3000.
(a)
(b)
(c)
(d)


## 3.1. Turing Pattern Formation in Domain I

In this subsection, we first consider the pattern formation when parameters are located in domain I of Figure1, which is pure Turing unstable domain.Initially, the entire system is placed in the steady state (u*,v*), and the propagation velocity of the initial perturbation is thus on the order of 5×10-4 space units per time unit. The initial data of prey u is a 101×101 matrix, such as
(26)u(x,y,t0)=(u*u*⋯u*u*u*⋯u*⋮⋮⋯⋮u*u*⋯u*)+(rand(101,101)-0.5)×10-3,t0∈[-τ,0],
where rand(101,101) is MATLAB procedure, which returns a 101×101 matrix containing pseudorandom values drawn from the standard uniform distribution on the open interval (0,1). And the system is then integrated for 2,000,000 time steps, and the last images are saved. After the initial period during which the perturbation spreads, the system goes into an essentially steady state (time independent).In the numerical simulations, different types of dynamics are observed, and it is found that the distributions of predator and prey are always of the same type. Consequently, we can restrict our analysis of pattern formation to one distribution. In the following, we show the distribution of preyu.In Figure2, we show three typical Turing patterns and their three-dimension pattern for model (4) at t=2000. In the case of (d,d21,d22)=(0.08,0.89,0.37), the equilibrium is (u*,v*)=(2.8697,0.5203) and the holes patterns are isolated zones with low prey density (cf. Figure 2(a)). When increasing d to 0.4, (d,d21,d22)=(0.47,0.89,0.37), (u*,v*)=(0.9877,1.2218), the stripes patterns appear. While with (d,d21,d22)=(0.61,0.39,0.75), (u*,v*)=(0.3122,0.7306), spots patterns are isolated zones with high prey density (cf. Figure 2(c)). And Figures 2(d)-2(e) are three-dimension patterns corresponding to Figures 2(a)–2(c).Stationary Turing patterns ofu in model (4) at time t=2000 with parameters a=0.43, b=0.22, c=0.4, r=1.4, m=0.17, d11=0.4, d12=0.05, d22=0.37, τ=0.46, and (a) (d) (d,d21,d22)=(0.08,0.89,0.37), (b) (e) (d,d21,d22)=(0.4,0.89,0.37), and (c) (f) (d,d21,d22)=(0.61,0.39,0.75), and the other parameters are fixed as (25). Times: t=2000.
(a)
(b)
(c)
(d)
(e)
(f)On the other hand, it is easy to check that holes, stripes, and spots also appear in model (4) without time delay for τ=0 and other the parameters are taken the same way as in Figures 2(a)–2(c), respectively. However, in parameters of Figure 2, trace(Mj)<0 and det(Mj)>0 for all j≥0 if the cross-diffusion coefficients d12=0 and d21=0. That is to say, model (4) is stable in absence of cross-diffusion, and this instability is induced by cross-diffusion.


## 3.2. Chaotic Pattern Formation Domain II

In this subsection, we will focus on the pattern formation when parameters are located in domain II of Figure1, the domain of pure Hopf-Turing instability.Figure3 shows the evolution process of the spatial pattern of u at time t=0,100,1000,2000 with parameter set (25) and a=0.1, b=0.71, r=1.54, d=0.47, d11=0.2, d21=0.44, d22=0.88, and  τ=0.46. The pattern takes a long time to settle down; starting with a homogeneous state E*=(56.4750,55.3357) (cf. Figure 3(a)), the random initial distribution leads to the formation of a strongly irregular transient pattern in the domain. After the irregular pattern form (cf. Figure 3(b)), it grows slightly and “jumps” alternately for a certain time, and finally the chaos patterns prevail over the whole domain (cf. Figures 3(c)-3(d)), which is time-dependent. For the sake of learning the spatiotemporal dynamics of model (4) further, we illustrate the time-series plots and phase portraits in Figures 3(e)-3(f). From the time-series plots in Figure 3(e), one can see that u and v strongly oscillates over time t. And in Figure 3(f), there exhibits a “local” phase plane of the system invaded by the irregular spatiotemporal oscillations; the trajectory now fills nearly the whole domain inside the quasi-limit cycle. This regime of the system dynamics corresponds to spatiotemporal chaos. And the spatial symmetry of model (4) is broken, and the patterns are oscillatory in space.Chaotic patterns ofu in model (4) with parameters a=0.01, b=0.22, c=0.4, r=1.54, m=0.17, d11=0.2, d12=0.05, d21=0.44, d22=0.88, and  τ=0.46. Times: (a) t=0, (b) t=100, (c) t=1000, (d) t=2000, (e) time-series plots, and (f) phase portraits.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.3. Spiral Wave Pattern Formation Domain IV

In this subsection, we will focus on the pattern formation when parameters are located in domain IV of Figure1, the domain of pure Hopf instability.In order to make any influence of the corners of the domain more visible, and thanks to the insightful work of Medvinsky et al. [4] and Upadhyay et al. [54], we have studied the spiral wave pattern for the following initial condition (27). In this case, we employ Δx=Δy=1 and the system size is 200×200.In order to make any influence of the corner of the domain more visible, following [43, 55], the initial conditions are deliberately chosen to be unsymmetrical as
(27)u(x,y,t0)=u*-ε1(x-100),v(x,y,t0)=v*-ε2(y-100),t0∈[-τ,0],
where ε1=2×10-4 and ε2=3×10-4. For more clarity, we put the initial data of prey u as follow
(28)u(x,y,t0)=(u*+ε1(100-0)u*+ε1(100-0)⋯u*+ε1(100-0)u*+ε1(100-1)u*+ε1(100-1)⋯u*+ε1(100-1)⋮⋮⋯⋮u*+ε1(100-200)u*+ε1(100-200)⋯u*+ε1(100-200)),t0∈[-τ,0].Figure4 shows the evolution process of the spatial pattern of prey u at t=0,500,1500,3000 with (d,d21,d22)=(0.5,0.15,0.37); the other parameters are fixed as in (25). Figure 4(a) shows that, for the model (4) with initial conditions (27), the formation of the irregular patchy structure can be preceded by the evolution of a regular spiral spatial pattern. Note that the appearance of the spirals is not induced by the initial conditions. The center of each spiral is situated in a critical point (xcr,ycr)=(100,100), where u(xcr,ycr)=u*,v(xcr,ycr)=v*. After the spirals form (Figure 4(b)), they grow slightly for a certain time and their spatial structure becomes more distinct (Figures 4(c) and 4(d)). Obviously, spiral pattern arises from Hopf instability.The process of formation of spiral pattern ofu in model (4) with parameters a=0.43, b=0.22, c=0.4, r=1.4, m=0.17, d=0.5, d11=0.4, d12=0.05, d21=0.15, d22=0.37, and τ=0.46. Times: (a) t=0, (b) t=500, (c) t=1500, and (d) t=3000.
(a)
(b)
(c)
(d)


## 4. Concluding Remarks

In this paper, we analyze the spatiotemporal dynamics of a spatial Leslie-Gower predator-prey model with time delay and diffusion under the zero-flux boundary conditions. The value of this study is twofold. First, it gives the conditions of Hopf and Turing instability induced by the effects of time delay and cross-diffusion. Second, it gives three domains for three kinds of patterns of model (4) and finds that time delay and cross-diffusion may induce Turing instability, resulting in stationary Turing patterns, Hopf instability, resulting in spiral wave patterns, and Turing-Hopf instability, resulting in chaotic wave patterns.Ifτ is small enough, it is worthy to note that, for the delayed model (4), two typical patterns are similar to those in our previous work [43] (a special case of model (4) with self diffusion): pure Turing instability gives birth to the holes, stripes, and spots patterns (cf. Figure 2) and pure Hopf instability to spiral wave pattern (cf. Figure 4). But this is not forever. In the present paper, we find that the delayed reaction-diffusion model (4) exhibits more complex pattern formation, for there exists the third pattern, chaotic wave pattern, which is induced by Turing-Hopf instability. The complexity of the pattern formation is induced by the cross-diffusion effect.The results of the formation of these patterns may indicate the vital role of phase transience regimes in the spatiotemporal organization of the prey and predator in the real world situation. It is also important to distinguish between “intrinsic” patterns, that is, Turing patterns arising due to predation, and “forced” patterns induced by the heterogeneity of the environment, that is, spiral wave pattern. Because of the environmental heterogeneity, the dispersion of varying quantities and typical times and lengths can be essentially different in different cases and thus can induce different spatial patterns.



---

*Source: 147232-2013-11-07.xml*