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

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.

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) where 𝑢 and 𝑣 are prey and predator densities, respectively. 𝑟, 𝑎, 𝑑, 𝑏_𝑖, and 𝑐_𝑖(𝑖=1,2) are positive constants. 𝑟 is the growth rate of preys, 𝑑 is the growth rate of predators, 𝑎 is the strength of competition among individuals of species 𝑢, 𝑏₁ (resp., 𝑏₂) measures the extent to which environment provides protection to prey 𝑢 (resp., to the predator 𝑣), and 𝑐₁ is the maximum value of the per capita reduction of 𝑢 due to 𝑣. 𝑐₂ has a similar meaning to 𝑐₁. 𝑑₁ and 𝑑₂ are the diffusion coefficients that imply the movement of individuals from a higher to a lower concentration region. ∇²=(𝜕²/𝜕𝑥²)+(𝜕²/𝜕𝑦²) 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 𝑚𝑢 of the prey in the case of 𝑐₁=c₂≠0 and 𝑏₁=𝑏₂≠0, where 𝑚∈[0,1) is a constant. This leaves (1−𝑚)𝑢 of the prey available to the predator, and model (1) becomes 𝜕𝑢()−()𝜕𝑡=𝑢𝑟−𝑎𝑢1−𝑚𝑢𝑣𝑏+(1−𝑚)𝑢+𝑑₁₁∇²𝑢≜𝑓(𝑢,𝑣)+𝑑₁₁∇²𝑢,𝜕𝑣𝜕𝑡=𝑣𝑑−𝑐𝑣𝑏+(1−𝑚)𝑢+𝑑₂₂∇²𝑣≜𝑔(𝑢,𝑣)+𝑑₂₂∇²𝑣≜𝑣(𝑑−ℎ(𝑢,𝑣))+𝑑₂₂∇²𝑣. (2) 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: 𝜕𝑢()−()𝜕𝑡=𝑢𝑟−𝑎𝑢1−𝑚𝑢𝑣𝑏+(1−𝑚)𝑢+𝑑₁₁∇²𝑢≜𝑓(𝑢,𝑣)+𝑑₁₁∇²𝑢,𝜕𝑣()𝜕𝑡=𝑣𝑑−𝑐𝑣𝑡−𝜏𝑏+(1−𝑚)𝑢(𝑡−𝜏)+𝑑₂₂∇²𝑣≜𝑣{𝑑−ℎ[𝑢(𝑡−𝜏),𝑣(𝑡−𝜏)]}+𝑑₂₂∇²𝑣, (3) where 𝜏>0 is a constant delay due to gestation and 𝑢(𝑥,𝑦,𝑡−𝜏) and 𝑣(𝑥,𝑦,𝑡−𝜏) are densities of prey and predator, respectively, at moment 𝑡−𝜏 and position (𝑥,𝑦). For simplicity, set 𝑢(𝑥,𝑦,𝑡−𝜏)∶=𝑢(𝑡−𝜏),𝑣(𝑥,𝑦,𝑡−𝜏)∶=𝑣(𝑡−𝜏). 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) where the cross-diffusion coefficients 𝑑₁₂ and 𝑑₂₁ 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 𝑑𝐷=¹¹𝑑₁₂d₂₁𝑑₂₂ (5) the diffusive matrix and assume that 𝑑₁₁>0, 𝑑₂₂>0, and det(𝐷)=𝑑₁₁𝑑₂₂−𝑑₁₂𝑑₂₁>0.

Following [44], we explain the meaning and units of each variable and constant in model (4).(i)The prey 𝑢 and predator 𝑣 are usually measured in mg of dry weight per liter [mg·dw/l].(ii)Time 𝑡 and time delay 𝜏 are measured in days [𝑑].(iii)The length of 𝑥,𝑦∈(0,𝐿) is measured in meter [m].(iv)𝑟 is the growth rate of prey 𝑢 and 𝑑 is the growth rate of predator 𝑣, measured in [d⁻¹].(v)𝑐 is the maximum value of the per capita reduction of 𝑣 due to 𝑢, which is not available in abundance and is measured in [d⁻¹].(vi)𝑎 describes the strength of competition among individuals of species 𝑢 and is measured in [(mg⋅dw/l)⁻¹d⁻¹].(vii)𝑏 describes the extent to which environment provides protection to prey 𝑢 and predator 𝑣 and is usually measured in mg of dry weight per liter [mg·dw/l].(viii)𝑑_𝑖𝑗>0(𝑖,𝑗=1,2) are the diffusion coefficients corresponding to the prey and predator, respectively, and are measured in [m²d⁻¹].

Model (4) is to be analyzed under the following nonzero initial conditions: 𝑢𝑥,𝑦,𝑡₀>0,𝑣𝑥,𝑦,𝑡₀()×(),𝑡>0,𝑥,𝑦∈Ω=0,𝐿0,𝐿₀∈[−𝜏,0], (6) and zero-flux boundary conditions: 𝜕𝑢𝜕𝑛_𝜕Ω=𝜕𝑣𝜕𝑛_𝜕Ω=0. (7) In the above, 𝐿 denotes the size of the system in square domain, and 𝑛 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 Section 3, we perform a series of numerical simulation to show evolution process of prey 𝑢. 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 𝑢(𝑥,𝑦,𝑡−𝜏) and 𝑣(𝑥,𝑦,𝑡−𝜏) as 𝑢𝑥,𝑦,𝑡−𝜏=𝑢𝑥,𝑦,𝑡−𝜏𝜕𝑢𝑥,𝑦,𝑡,𝑣𝜕𝑡𝑥,𝑦,𝑡−𝜏=𝑣𝑥,𝑦,𝑡−𝜏𝜕𝑣𝑥,𝑦,𝑡.𝜕𝑡 (8) Substituting (8) into (4), we obtain 𝜕𝑢()𝜕𝑡=𝑓𝑢,𝑣+𝑑₁₁∇²𝑢+𝑑₁₂∇²𝑣,𝜕𝑣𝑢𝜕𝑡=𝑣𝑑−ℎ𝑥,𝑦,𝑡−𝜏𝜕𝑢𝑥,𝑦,𝑡,𝑣𝜕𝑡𝑥,𝑦,𝑡−𝜏𝜕𝑣𝑥,𝑦,𝑡𝜕𝑡+𝑑₂₁∇²𝑢+𝑑₂₂∇²𝑣. (9) Expanding (9) in Taylor series and neglecting the higher-order nonlinearities, then (9) becomes 𝜕𝑢()𝜕𝑡=𝑓𝑢,𝑣+𝑑₁₁∇²𝑢+𝑑₁₂∇²𝑣,𝜕𝑣()𝜕𝑡=𝑣𝑑−ℎ𝑢,𝑣+𝜏ℎ_𝑢(𝑢,𝑣)𝜕𝑢𝑥,𝑦,𝑡𝜕𝑡+𝜏ℎ_𝑣(𝑢,𝑣)𝜕𝑣𝑥,𝑦,𝑡𝜕𝑡+𝑑₂₁∇²𝑢+𝑑₂₂∇²𝑣, (10) where ℎ_𝑢(𝑢,𝑣)=(𝜕ℎ(𝑢,𝑣)/𝜕𝑢) and ℎ_𝑣(𝑢,𝑣)=(𝜕ℎ(𝑢,𝑣)/𝜕𝑣).

From (10), we finally obtain that 𝜕𝑢()𝜕𝑡=𝑓𝑢,𝑣+𝑑₁₁∇²𝑢+𝑑₁₂∇²𝑣,𝜕𝑣=1𝜕𝑡1−𝑣𝜏ℎ_𝑣()×𝑔()𝑢,𝑣𝑢,𝑣+𝑣𝜏ℎ_𝑢()𝑓()+𝑢,𝑣𝑢,𝑣𝑣𝜏ℎ_𝑢(𝑢,𝑣)𝑑₁₁+𝑑₂₁∇²𝑢+𝑣𝜏ℎ_𝑢(𝑢,𝑣)𝑑₁₂+𝑑₂₂∇²𝑣. (11)(12) 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 𝐸₁=(0,0), 𝐸₂=((𝑎/𝑟),0), and 𝐸₃=(0,(𝑏𝑑/𝑟)) and a unique positive equilibrium 𝐸^∗=(𝑢^∗,𝑣^∗)=(((𝑑𝑚−𝑑+𝑐𝑟)/𝑎𝑐),(((𝑐𝑎𝑏+(𝑑𝑚−𝑑+𝑐𝑟)(1−𝑚))𝑑)/𝑎𝑐²)).

Remark 1. In the case without time delay and cross-diffusion, that is, 𝜏=0 and 𝑑₁₂=𝑑₂₁=0, model (11) becomes (2). It is easy to see that 𝐸₁,𝐸₂,𝐸₃, and 𝐸^∗ 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).

Let 0=𝜇₀<𝜇₁<𝜇₂<⋯ be the eigenvalues of the operator −Δ on Ω with the zero-flux boundary conditions. And setting 𝐗=(𝑢,𝑣)∈𝐻²(Ω)²∣𝜕_𝑛𝑢=𝜕_𝑛,𝐸𝜇=𝑣=0on𝜕Ω𝜙∣−Δ𝜙=𝜇𝜙inΩ,𝜕_𝜈𝜙=0on𝜕Ω,with𝜇∈ℝ¹, (13){𝜙_𝑖𝑗∣𝑗=1,…,dim𝐸(𝜇_𝑖)} will be an orthonormal basis of 𝐸(𝜇_𝑖), and 𝐗_𝑖𝑗={𝐜𝜙_𝑖𝑗∣𝐜∈ℝ²}; then 𝐗=⨁^∞_𝑖=1𝐗_𝑖, where 𝐗_𝑖=⨁^{dim𝐸(𝜇𝑖)}_𝑗=1𝐗_𝑖𝑗.

The linearization of model (11) around 𝐸^∗=(𝑢^∗,𝑣^∗) can be expressed by 𝜕𝑢𝑣𝜕𝑡=𝐽_𝐸^∗𝑢𝑣+𝐷∇²𝑢𝑣, (14) with domain 𝐗. Here, 𝐽_𝐸^∗=𝑓_𝑢𝑓_𝑣𝒳𝑔_𝑢+𝑣^∗𝜏ℎ_𝑢𝑓_𝑢𝒳𝑔_𝑣+𝑣^∗𝜏ℎ_𝑢𝑓_𝑣≜𝐽₁₁𝐽₁₂𝐽₂₁𝐽₂₂,𝑑𝐷=¹¹𝑑₁₂𝒳𝑣^∗𝜏ℎ_𝑢𝑑₁₁+𝑑₂₁𝒳𝑣^∗𝜏ℎ_𝑢𝑑₁₂+𝑑₂₂, (15) with 𝑓_𝑢=(𝜕𝑓/𝜕𝑢)∣_{(𝑢∗,𝑣∗)}, 𝑓_𝑣=(𝜕𝑓/𝜕𝑣)∣_{(𝑢∗,𝑣∗)}, 𝑔_𝑢=(𝜕𝑔/𝜕𝑢)∣_{(𝑢∗,𝑣∗)}, 𝑔_𝑣=(𝜕𝑔/𝜕𝑣)∣_{(𝑢∗,𝑣∗)}, ℎ_𝑢=(𝜕ℎ/𝜕𝑢)∣_{(𝑢∗,𝑣∗)}, ℎ_𝑣=(𝜕ℎ/𝜕𝑣)∣_{(𝑢∗,𝑣∗)}, and (1/(1−𝑣𝜏ℎ_𝑣)|_{(𝑢∗,𝑣∗)})=(1/(1−𝜏𝑑))≜𝒳. Since 𝜏 is small, we only consider this case of 𝜏<(1/𝑑) (i.e. 𝒳>0) in this paper.

For each 𝑗(𝑗=0,1,2,…), 𝐗_𝑗 is invariant under the operator 𝐷∇²+𝐿_𝐸^∗, and 𝜆 is an eigenvalue of 𝐷∇²+𝐿_𝐸 on 𝐗_𝑗 if and only if 𝜆 is an eigenvalue of the matrix 𝑀_𝑗=−𝐷𝜇_𝑗+𝐽_𝐸^∗ for some 𝑗≥0.

So the stability is reduced to consider the characteristic equation det𝜆𝐼−𝑀_𝑗=𝜆²−trace𝑀_𝑗𝜆+det𝑀_𝑗=0, (16) where trace𝑀_𝑗=trace𝐷𝜇_𝑗−trace𝐽_𝐸^∗,det𝑀_𝑗=𝒳det(𝐷)𝜇²_𝑗−𝒳𝑑₁₁𝑔_𝑣−𝑑₁₂𝑔_𝑢−𝑑₂₁𝑓_𝑣+𝑑_22𝑓𝑢𝜇_𝑗+det𝐽_𝐸^∗. (17)

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 𝐸^∗, the matrix 𝑀_𝑗=−𝐷𝜇_𝑗+𝐽_𝐸^∗ must have an eigenvalues on the imaginary axis [52]; that is, trace(𝑀_𝑗)≡0. Therefore, the only possible critical values of 𝜏 are 𝜏(𝑖) such that, for 𝑖≥0, trace𝐷𝜇_𝑖−trace𝐽_𝐸^∗=trace𝐷𝜇_𝑖−𝑓_𝑢+𝑔_𝑣+𝑣^∗𝜏(𝑖)ℎ_𝑢𝑓_𝑣𝒳=0. (18) At 𝜏=𝜏(𝑖), trace(𝑀_𝑖)=trace(𝐷(𝜇_𝑖−𝜇_𝑗)), and therefore, if 𝑖≥1, for all 0≤𝑗<𝑖, trace(𝑀_𝑖)>0 and the matrix 𝑀_𝑗=−𝐷𝜇_𝑗+𝐽_𝐸^∗ has at least 2𝑖 eigenvalues with positive real parts. Therefore, the only value of 𝜏 at which Hopf bifurcation hypotheses may be satisfied is 𝜏_𝑐(0)𝑓=𝜏∶=^𝑢+𝑔_𝑣𝑣^∗𝑓_𝑢ℎ_𝑣−ℎ_𝑢𝑓_𝑣. (19) Near 𝜏_𝑐, the complex conjugate pair 𝜅(𝜏)±𝑖𝜔(𝜏) is given by 𝜅(𝜏)=12trace𝐽_𝐸^∗,𝜔²(𝜏)=det𝐽_𝐸^∗−𝜅²(𝜏). (20)

Since 𝜅(𝜏_𝑐)=0 and det(𝐽_𝐸^∗)=(((𝑑𝑚+𝑟𝑐−𝑑)𝑑)/𝑐)𝒳>0, there exists an interval 𝐼 containing 𝜏_𝑐 such that(i)𝜔²(𝜏)≠0 for all 𝜏∈𝐼;(ii)trace(𝑀_𝑗)<0 and det(𝑀_𝑗)>0 with 𝑑₁₁𝑔_𝑣−𝑑₁₂𝑔_𝑢−𝑑₂₁𝑓_𝑣+𝑑_22𝑓𝑢<0 uniformly for 𝑗≥1 and 𝜏∈𝐼;(iii)(𝑑/𝑑𝜏)𝜅(𝜏)∣_{𝜏=𝜏𝑐}=𝑣^∗𝑓_𝑣ℎ_𝑢+(𝑑𝑔_𝑣/(1−𝑑𝜏_𝑐)²)>0.

Thus, we have the following theorem.

Theorem 2. If 𝑑₁₁𝑔_𝑣−𝑑₁₂𝑔_𝑢−𝑑₂₁𝑓_𝑣+𝑑_22𝑓𝑢<0, then model (11) undergoes a Hopf bifurcation that occurs around the positive equilibrium 𝐸^∗ at 𝜏=𝜏_𝑐.

2.2. Turing Instability

In this subsection, we will state the Turing instability for the positive equilibrium 𝐸^∗ of model (11).

Mathematically speaking, the positive equilibrium 𝐸^∗ 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(𝑀_𝑗)<0 or the condition det(𝑀_𝑗)>0 is violated.

Since the positive equilibrium 𝐸^∗ is stable without diffusion which means that trace(𝐽_𝐸^∗)<0 and det(𝐽_𝐸^∗)>0 hold, then trace(𝑀_𝑗)<0 is always true. Therefore, for the emergency of the diffusion-driven instability in model (11), it is needed that det(𝑀_𝑗)<0 for some 𝜇_𝑗>0. A necessary condition is 𝑑₁₁𝑔_𝑣−𝑑₁₂𝑔_𝑢−𝑑₂₁𝑓_𝑣+𝑑_22𝑓𝑢>0; (21) otherwise det(𝑀_𝑗)>0 for all 𝑗>0 since 𝒳det(𝐷)>0 and det(𝐽_𝐸^∗)>0. And we notice that det(𝑀_𝑗) achieves its minimum: min_𝑗∈𝐑+det𝑀_𝑗=4𝒳det(𝐷)det𝐽_𝐸^∗−𝒳𝑑₁₁𝑔_𝑣−𝑑₁₂𝑔_𝑢−𝑑₂₁𝑓_𝑣+𝑑_22𝑓𝑢²×4det𝐷⁻¹ (22) at the critical value 𝜇²_𝑐>0 where 𝜇²_𝑐=𝑑₁₁𝑔_𝑣−𝑑₁₂𝑔_𝑢−𝑑₂₁𝑓_𝑣+𝑑_22𝑓𝑢(𝐷).2det (23) Summarizing the above calculation, we conclude the following.

Theorem 3. If 1𝜏<min𝑑,𝜏_𝑐,4𝒳det(𝐷)det𝐽_𝐸^∗−𝒳𝑑₁₁𝑔_𝑣−𝑑₁₂𝑔_𝑢−𝑑₂₁𝑓_𝑣+𝑑_22𝑓𝑢²𝑑<0,₁₁𝑔_𝑣−𝑑₁₂𝑔_𝑢−𝑑₂₁𝑓_𝑣+𝑑_22𝑓𝑢>0, (24) then the positive equilibrium 𝐸^∗ of model (11) is Turing unstable.

In Figure 1, we show the bifurcation diagram for model (11). The Turing instability curve and Hopf bifurcation curve separate the 𝑑−𝑑₂₁ 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 𝑑−𝑑21 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 𝑑𝑏=0.22,𝑐=0.4,𝑚=0.17,₁₂=0.05,𝜏=0.46. (25) The numerical integration of model (4) was performed by means of forward Euler integration, with a time step of Δ𝑡=0.001 and a space step of Δ𝑥=Δ𝑦=(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 Figure 1, which is pure Turing unstable domain.

Initially, the entire system is placed in the steady state (𝑢^∗,𝑣^∗), and the propagation velocity of the initial perturbation is thus on the order of 5×10⁻⁴ space units per time unit. The initial data of prey 𝑢 is a 101×101 matrix, such as 𝑢𝑥,𝑦,𝑡₀=𝑢^∗𝑢^∗⋯𝑢^∗𝑢^∗𝑢^∗⋯𝑢^∗𝑢⋮⋮⋯⋮_∗𝑢^∗⋯𝑢^∗+(rand(101,101)−0.5)×10⁻³,𝑡₀∈[−𝜏,0], (26) 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 prey 𝑢.

In Figure 2, we show three typical Turing patterns and their three-dimension pattern for model (4) at 𝑡=2000. In the case of (𝑑,𝑑₂₁,𝑑₂₂)=(0.08,0.89,0.37), the equilibrium is (𝑢^∗,𝑣^∗)=(2.8697,0.5203) and the holes patterns are isolated zones with low prey density (cf. Figure 2(a)). When increasing 𝑑 to 0.4, (𝑑,𝑑₂₁,𝑑₂₂)=(0.47,0.89,0.37), (𝑢^∗,𝑣^∗)=(0.9877,1.2218), the stripes patterns appear. While with (𝑑,𝑑₂₁,𝑑₂₂)=(0.61,0.39,0.75), (𝑢^∗,𝑣^∗)=(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).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2 

Stationary Turing patterns of 𝑢 in model (4) at time 𝑡=2000 with parameters 𝑎=0.43, 𝑏=0.22, 𝑐=0.4, 𝑟=1.4, 𝑚=0.17, 𝑑11=0.4, 𝑑12=0.05, 𝑑22=0.37, 𝜏=0.46, and (a) (d) (𝑑,𝑑21,𝑑22)=(0.08,0.89,0.37), (b) (e) (𝑑,𝑑21,𝑑22)=(0.4,0.89,0.37), and (c) (f) (𝑑,𝑑21,𝑑22)=(0.61,0.39,0.75), and the other parameters are fixed as (25). Times: 𝑡=2000.

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(𝑀_𝑗)<0 and det(𝑀_𝑗)>0 for all 𝑗≥0 if the cross-diffusion coefficients 𝑑₁₂=0 and 𝑑₂₁=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 Figure 1, the domain of pure Hopf-Turing instability.

Figure 3 shows the evolution process of the spatial pattern of 𝑢 at time 𝑡=0,100,1000,2000 with parameter set (25) and 𝑎=0.1, 𝑏=0.71, 𝑟=1.54, 𝑑=0.47, 𝑑₁₁=0.2, 𝑑₂₁=0.44, 𝑑₂₂=0.88, and  𝜏=0.46. The pattern takes a long time to settle down; starting with a homogeneous state 𝐸^∗=(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 𝑢 and 𝑣 strongly oscillates over time 𝑡. 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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

Chaotic patterns of 𝑢 in model (4) with parameters 𝑎=0.01, 𝑏=0.22, 𝑐=0.4, 𝑟=1.54, 𝑚=0.17, 𝑑11=0.2, 𝑑12=0.05, 𝑑21=0.44, 𝑑22=0.88, and  𝜏=0.46. Times: (a) 𝑡=0, (b) 𝑡=100, (c) 𝑡=1000, (d) 𝑡=2000, (e) time-series plots, and (f) phase portraits.

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 Figure 1, 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 Δ𝑥=Δ𝑦=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 𝑢𝑥,𝑦,𝑡₀=𝑢^∗−𝜀₁(),𝑣𝑥−100𝑥,𝑦,𝑡₀=𝑣^∗−𝜀₂,𝑡𝑦−100₀∈[−𝜏,0], (27) where 𝜀₁=2×10⁻⁴ and 𝜀₂=3×10⁻⁴. For more clarity, we put the initial data of prey 𝑢 as follow 𝑢𝑥,𝑦,𝑡₀=𝑢^∗+𝜀₁(100−0)𝑢^∗+𝜀₁(100−0)⋯𝑢^∗+𝜀₁()𝑢100−0_∗+𝜀₁(100−1)𝑢^∗+𝜀₁(100−1)⋯𝑢^∗+𝜀₁()𝑢100−1⋮⋮⋯⋮_∗+𝜀₁(100−200)𝑢^∗+𝜀₁(100−200)⋯𝑢^∗+𝜀₁(),𝑡100−200₀∈[−𝜏,0]. (28)

Figure 4 shows the evolution process of the spatial pattern of prey 𝑢 at 𝑡=0,500,1500,3000 with (𝑑,𝑑₂₁,𝑑₂₂)=(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 (𝑥_cr,𝑦_cr)=(100,100), where 𝑢(𝑥_cr,𝑦_cr)=𝑢^∗,𝑣(𝑥_cr,𝑦_cr)=𝑣^∗. 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.

(a)

(b)

(c)

(d)

Figure 4 

The process of formation of spiral pattern of 𝑢 in model (4) with parameters 𝑎=0.43, 𝑏=0.22, 𝑐=0.4, 𝑟=1.4, 𝑚=0.17, 𝑑=0.5, 𝑑11=0.4, 𝑑12=0.05, 𝑑21=0.15, 𝑑22=0.37, and 𝜏=0.46. Times: (a) 𝑡=0, (b) 𝑡=500, (c) 𝑡=1500, and (d) 𝑡=3000.

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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors would like to thank the anonymous referee for the very helpful suggestions and comments which led to the improvement of our original paper. This research was supported by the National Science Foundation of China (61373005), Zhejiang Provincial Natural Science Foundation (LY12A01014), and the Fund Project of Zhejiang Provincial Education Department (Y201223449).