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 [312]. 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] where and are prey and predator densities, respectively. , , , , and 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 [1820] 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 and , where is a constant. This leaves of the prey available to the predator, and model (1) becomes 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 [2229]. 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 [3133]. 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 [3444]. In our previous work [43], we studied the spatiotemporal complexity of a Leslie-Gower model incorporating a prey refuge with time delay effect: where 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 [4549]. 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: 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 the diffusive matrix and assume that , , and .

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 is measured in meter [m].(iv) is the growth rate of prey and is the growth rate of predator , measured in [].(v) is the maximum value of the per capita reduction of due to , which is not available in abundance and is measured in [].(vi) describes the strength of competition among individuals of species and is measured in .(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) are the diffusion coefficients corresponding to the prey and predator, respectively, and are measured in .

Model (4) is to be analyzed under the following nonzero initial conditions: and zero-flux boundary conditions: 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 Substituting (8) into (4), we obtain Expanding (9) in Taylor series and neglecting the higher-order nonlinearities, then (9) becomes where and .

From (10), we finally obtain that 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 , , and and a unique positive equilibrium .

Remark 1. In the case without time delay and cross-diffusion, that is, and , 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 be the eigenvalues of the operator − on with the zero-flux boundary conditions. And setting will be an orthonormal basis of , and ; then , where .

The linearization of model (11) around can be expressed by with domain . Here, with , , , , , , and . Since is small, we only consider this case of (i.e. ) in this paper.

For each , is invariant under the operator , and is an eigenvalue of on if and only if is an eigenvalue of the matrix for some .

So the stability is reduced to consider the characteristic equation where

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, . Therefore, the only possible critical values of are such that, for , At , , and therefore, if , for all , and the matrix has at least eigenvalues with positive real parts. Therefore, the only value of at which Hopf bifurcation hypotheses may be satisfied is Near , the complex conjugate pair is given by

Since and , there exists an interval containing such that(i) for all ;(ii) and with uniformly for and ;(iii).

Thus, we have the following theorem.

Theorem 2. If , 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 or the condition is violated.

Since the positive equilibrium is stable without diffusion which means that and hold, then is always true. Therefore, for the emergency of the diffusion-driven instability in model (11), it is needed that for some . A necessary condition is otherwise for all since and . And we notice that achieves its minimum: at the critical value where Summarizing the above calculation, we conclude the following.

Theorem 3. If 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.

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 space units. And some parameters are taken as The numerical integration of model (4) was performed by means of forward Euler integration, with a time step of and a space step of , 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 space units per time unit. The initial data of prey is a matrix, such as where is MATLAB procedure, which returns a matrix containing pseudorandom values drawn from the standard uniform distribution on the open interval . 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 . In the case of , the equilibrium is and the holes patterns are isolated zones with low prey density (cf. Figure 2(a)). When increasing to 0.4, , , the stripes patterns appear. While with , , 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).

On the other hand, it is easy to check that holes, stripes, and spots also appear in model (4) without time delay for and other the parameters are taken the same way as in Figures 2(a)2(c), respectively. However, in parameters of Figure 2, and for all if the cross-diffusion coefficients and . 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 with parameter set (25) and , , , , , , , and  . The pattern takes a long time to settle down; starting with a homogeneous state (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.

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 and the system size is .

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 where and . For more clarity, we put the initial data of prey as follow

Figure 4 shows the evolution process of the spatial pattern of prey at with ; 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 , where . 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.

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.


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).