Turing Patterns in a Predator-Prey System with Self-Diffusion

Abstract

For a predator-prey system, cross-diffusion has been confirmed to emerge Turing patterns. However, in the real world, the tendency for prey and predators moving along the direction of lower density of their own species, called self-diffusion, should be considered. For this, we investigate Turing instability for a predator-prey system with nonlinear diffusion terms including the normal diffusion, cross-diffusion, and self-diffusion. A sufficient condition of Turing instability for this system is obtained by analyzing the linear stability of spatial homogeneous equilibrium state of this model. A series of numerical simulations reveal Turing parameter regions of the interaction of diffusion parameters. According to these regions, we further demonstrate dispersion relations and spatial patterns. Our results indicate that self-diffusion plays an important role in the spatial patterns.

1. Introduction

In ecological systems, the interactions of different species indicate abundant dynamical features. It is informative to use mathematical model to study the interactions of species in these systems. Among these models, predator-prey systems, which were based on the pioneering works of Volterra [1], have been important in ecological problems. However, since we live in a spatial world, the predator-prey systems should include spatial factors. Thus, these systems should be described by using reaction-diffusion equations. As a result, it is an open problem to understand spatiotemporal behaviors of the temporal-spatial predator-prey systems [2–11]. Thereinto, the formation of spatial patterns of predator-prey systems is a very active research area [12–15], which is based on the pioneering work of Turing [16] in 1952.

In recent years, there are a lot of bodies of literature to study the predator-prey system by taking into account the normal diffusion as well as cross-diffusion [17–20]. Normal diffusion is a natural phenomenon of the movement of the prey or the predators from higher-density regions to lower-density ones. Cross-diffusion of the prey expresses a flux of the prey because of the presence of the predators and vice versa. Furthermore, in predator-prey systems, cross-diffusion can induce Turing instability to produce spatial patterns even though spatial homogeneous equilibrium states for the corresponding system in the absence of cross-diffusion are stable [8, 11, 21–28].

Besides the normal diffusion and the cross-diffusion of the predator and prey in ecological systems, there exists, in fact, another diffusion form—self-diffusion for the pressure of their own species. It can describe the tendency to move along the direction of lower density of the predator's and prey's own species [29]. Unfortunately, most of the studies mainly focused on well-posedness of solutions for predator-prey systems with self-diffusion [30, 31]. Little attention was paid to examine Turing patterns of these systems. Based on the above discussion, in this paper we mainly concentrate on Turing instability of a predator-prey system that includes a normal diffusion, cross-diffusion, and self-diffusion terms. To this end, we find a sufficient condition to generate Turing patterns. By using numerical simulation, for this system we examine parameter regions of forming patterns and show snapshots of spatial patterns.

The paper is organized as follows. In Section 2, we build the predator-prey model with nonlinear diffusion terms including normal diffusion, cross-diffusion, and self-diffusion terms and the biological meaning of these parameters are interpreted. In Section 3, we find the sufficient condition of Turing instability. By performing a series of numerical simulations, we locate the Turing parameter spaces when parts of the parameters are fixed in Section 4. In Section 5, by choosing values of some parameters from the Turing parameter spaces, we illustrate Turing patterns. Finally, some conclusions and discussions are given.

2. A Predator-Prey Model with Nonlinear Diffusion

In this paper, we are interested in the spatiotemporal patterns of the following predator-prey system with nonlinear-diffusion terms. Mathematical properties for this system have been investigated in [30, 31]: 𝜕𝑢₁𝑑𝜕𝑡−Δ₁+𝑎₁₁𝑢₁+𝑎₁₂𝑢₂𝑢₁=𝑢₁𝑎₁−𝑏₁𝑢₁−𝑐₁𝑢₂,inΩ×[0,∞),𝜕𝑢₂𝑑𝜕𝑡−Δ₂+𝑎₂₁𝑢₁+𝑎₂₂𝑢₂𝑢₂=𝑢₂𝑎₂+𝑏₂𝑢₁−𝑐₂𝑢₂,inΩ×[0,∞),𝜕𝑢₁=𝜕𝑛𝜕𝑢₂[),𝑢𝜕𝑛=0,on𝜕Ω×0,∞₁(𝑥,0)=𝑢⁰₁(𝑥)𝑢≥0,₂(𝑥,0)=𝑢⁰₂(𝑥)≥0,inΩ, (1) where Ω is a bounded domain in ℝ² with smooth boundary 𝜕Ω and represents the domain that these two species inhabit. The vector 𝑛 is the outward unit normal vector of Ω. In this model, 𝑢₁ and 𝑢₂ are the prey and the predator densities, respectively; 𝑎_𝑖, 𝑏_𝑖, and 𝑐_𝑖 (𝑖=1,2) are positive constants. 𝑎₁,𝑎₂ are rates of the prey and predator proliferations for food source; 𝑎₁/𝑏₁ and 𝑎₂/𝑐₂ are environmental carrying capacities for prey and predator, respectively. 𝑐₁ is a consumption rate; 𝑏₂ is a conversion rate. In the diffusion terms, the constant 𝑑_𝑖(𝑖=1,2), which is usually termed as a normal-diffusion coefficient, represents the natural dispersive force of movement of a species. The positive constants 𝑎₁₂ and 𝑎₂₁ are referred to as cross-diffusion coefficients, which describe that the prey tends to avoid higher density of the predators and vice versa by diffusing away. In addition, for the predators and the prey, the positive constants 𝑎₁₁ and 𝑎₂₂ are self-diffusion rates due to pressure within their own species.

Next, we want to look for the condition on the parameter values such that a positive homogeneous equilibrium state is linearly stable in the absence of the cross-diffusion and the self-diffusion (i.e., a normal reaction-diffusion system) but unstable in the present of the cross-diffusion and the self-diffusion.

For simplicity, we set up the following notation as in [22, 24].

Notation 1. Let 0=𝜇₁<𝜇₂<⋯→∞ be the eigenvalues of −Δ on Ω under no-flux boundary condition, and let 𝐸(𝜇_𝑖) be the space of eigenfunction corresponding to 𝜇_𝑖. We define the following space decomposition:(i)𝐗_𝑖𝑗∶={𝜑_𝑖𝑗𝐜∶𝐜∈ℝ²}, where 𝜑_𝑖𝑗 are orthonormal bases of 𝐸(𝜇_𝑖) for 𝑗=1,…,dim𝐸(𝜇_𝑖);(ii)𝐗∶={𝐮∈[𝒞¹(Ω)]²∶𝜕𝑢₁/𝜕𝑛=𝜕𝑢₂/𝜕𝑛=0on𝜕Ω}, and thus 𝐗=⊕^∞_𝑖=1𝐗_𝑖, where 𝐗_𝑖=⊕^{dim𝐸(𝜇𝑖)}_𝑗=1𝐗_𝑖𝑗.

Notation 2. For the sake of simplicity, we denote reaction terms for systems (1) by 𝐆(𝐮)=𝐺₁(𝐮)𝐺₂(𝐮)=𝑢₁𝑔₁(𝐮)∶=𝑢₁𝑎₁−𝑏₁𝑢₁−𝑐₁𝑢₂𝑢₂𝑔₂(𝐮)∶=𝑢₂𝑎₂+𝑏₂𝑢₁−𝑐₂𝑢₂, (2) and its Jacobian matrix at the point 𝐮^∗ is 𝐆_𝐮𝐮^∗=𝐺₁₁𝐺₁₂𝐺₂₁𝐺₂₂∶=−𝑏₁𝑢^∗₁−𝑐₁𝑢^∗₁𝑏₂𝑢^∗₂−𝑐₂𝑢^∗₂. (3) The diffusion term is denoted as Φ(𝐮)=Φ₁(𝐮)Φ₂(𝐮)𝑑∶=¹+𝑎₁₁𝑢₁+𝑎₁₂𝑢₂𝑢₁𝑑₂+𝑎₂₁𝑢₁+𝑎₂₂𝑢₂𝑢₂, (4) and its corresponding Jacobian matrix at the point 𝐮^∗ is Φ_𝐮𝐮^∗=Φ₁₁Φ₁₂Φ₂₁Φ₂₂𝑑∶=¹+2𝑎₁₁𝑢^∗₁+𝑎₁₂𝑢^∗₂𝑎₁₂𝑢^∗₁𝑎₂₁𝑢^∗₂𝑑₂+𝑎₂₁𝑢^∗₁+2𝑎₂₂𝑢^∗₂. (5)

3. Linear Stability Analysis of System (1) with a Normal Diffusion

Let 𝑎_𝑖𝑗=0 for 𝑖,𝑗=1,2; then the system (1) degenerates into standard reaction-diffusion equations: 𝜕𝑢₁𝜕𝑡−𝑑₁Δ𝑢₁=𝑢₁𝑎₁−𝑏₁𝑢₁−𝑐₁𝑢₂,inΩ×[0,∞),𝜕𝑢₂𝜕𝑡−𝑑₂Δ𝑢₂=𝑢₂𝑎₂+𝑏₂𝑢₁−𝑐₂𝑢₂,inΩ×[0,∞),𝜕𝑢₁=𝜕𝑛𝜕𝑢₂[),𝑢𝜕𝑛=0,on𝜕Ω×0,∞₁(𝑥,0)=𝑢⁰₁(𝑥)𝑢≥0,₂(𝑥,0)=𝑢⁰₂(𝑥)≥0,inΩ. (6) Through this paper, we assume the following conditions: []𝑐𝐶1²𝑐₁>𝑎₂𝑎₁,[]𝑏𝐶2²𝑏₁>𝑐₂𝑐₁−2𝑎₂𝑎₁. (7)

Then, there exists a unique positive equilibrium state for system (6), denoted by 𝐮^∗=(𝑢^∗₁,𝑢^∗₂), where 𝑢^∗₁=𝑎₁𝑐₂−𝑎₂𝑐₁𝑏₁𝑐₂+𝑏₂𝑐₁,𝑢^∗₂=𝑎₁𝑏₂+𝑎₂𝑏₁𝑏₁𝑐₂+𝑏₂𝑐₁. (8)

Theorem 1. If there are no cross-diffusion and self-diffusion, the positive equilibrium state 𝐮^∗ of the system (6) is locally asymptotically stable when condition [𝐶1] holds.

Proof. The linearization of (6) around the steady state 𝐮^∗ can be therefore expressed by 𝜕𝐮=𝜕𝑡𝐷Δ+𝐆_𝐮𝐮^∗𝐮, (9) where 𝐷=diag(𝑑₁,𝑑₂). Obviously, the operator 𝐷Δ+𝐆_𝐮(𝐮^∗) is invariant in the subspace 𝐗_𝑖, and 𝜆_𝑖 is an eigenvalue of this operator on 𝐗_𝑖, if and only if it is an eigenvalue of −𝜇_𝑖𝐷+𝐆_𝐮(𝐮^∗). Direct calculation shows the characteristic equation 𝜆²_𝑖−𝑚_𝑖𝜆_𝑖+𝑛_𝑖=0, (10) where 𝑚_𝑖=−𝑢^∗₁𝑏₁−𝑢^∗₂𝑐₂−(𝑑₁+𝑑₂)𝜇_𝑖 and 𝑛_𝑖=(𝑢^∗₁𝑏₁+𝑑₁𝜇_𝑖)(𝑢^∗₂𝑐₂+𝑑₂𝜇_𝑖)+𝑢^∗₁𝑢^∗₂𝑏₂𝑐₁. It is easy to be verified that 𝑚_𝑖 is negative and 𝑛_𝑖 positive. Thus, for each 𝑖>1, the two roots of (10) have negative real parts. Consequently, we complete this proof.

Remark 2. According to the proof of Theorem 1, we can also calculate out the eigenvector 𝐜 corresponding to the eigenvalue 𝜆_𝑖 for the operator −𝜇_𝑖𝐷+𝐆_𝐮(𝐮^∗), which satisfies 𝜆_𝑖+𝜇_𝑖𝐷−𝐆_𝐮𝐮^∗𝐜=0. (11) Furthermore, we can yield corresponding eigenspace {𝜑_𝑖𝑗𝐜} in 𝐗_𝑖.

By Theorem 1, under condition [𝐶1] system (6) cannot destabilize 𝐮^∗. Next, the cross-diffusion and self-diffusion are taken into account.

Theorem 3. Assume that conditions [𝐶1] and [𝐶2] hold; then there exists a sufficiently large positive constant 𝑎^∗₂₁ such that the homogeneous state 𝐮^∗ of the system (1) is unstable provided that 𝑎₂₁>𝑎^∗₂₁.

Proof. We first linearize the system (1) around (𝑢₁,𝑢₂)=(𝑢^∗₁,𝑢^∗₂) and arrive at 𝜕𝐮=Φ𝜕𝑡^𝐮𝐮^∗Δ+𝐆_𝐮𝐮^∗𝐮. (12) By the same method as the proof of Theorem 1, we consider the operator Φ_𝐮(𝐮^∗)Δ+𝐆_𝐮(𝐮^∗) on the subspace 𝐗_𝑖. The eigenvalue of this operator on 𝐗_𝑖 is denoted as 𝜆, and then it is also an eigenvalue of the matrix −𝜇_𝑖Φ_𝐮(𝐮^∗)+𝐆_𝐮(𝐮^∗). Thus, 𝜆 satisfies the following equation: 𝜆²−𝐺₁₁+𝐺₂₂−𝜇_𝑖Φ₁₁+Φ₂₂𝜆+det𝐆_𝐮𝐮^∗−𝜇_𝑖Φ_𝐮=0. (13) According to Notation 2, we can calculate and get 𝐺₁₁+𝐺₂₂−𝜇_𝑖Φ₁₁+Φ₂₂=−𝜇_𝑖𝑑₁+𝑑₂+2𝑎₁₁+𝑎₂₁𝑢^∗₁+𝑎₁₂+2𝑎₁₂𝑢^∗₂−𝑏₁𝑢^∗₁−𝑐₂𝑢^∗₂𝐆<0,det_𝐮𝐮^∗−𝜇_𝑖Φ_𝐮=𝑎𝜇²_𝑖+𝑏𝜇_𝑖+𝑐, (14)(15) where 𝑎=𝑑₁𝑑₂+𝑑₁𝑎₂₁𝑢^∗₁+2𝑎₂₂𝑢^∗₂+𝑑₂2𝑎₁₁𝑢^∗₁+𝑎₁₂𝑢^∗₂+2𝑎₁₁𝑎₂₁𝑢^∗₁²+4𝑎₁₁𝑎₂₂𝑢^∗₁𝑢^∗₂+2𝑎₁₂𝑎₂₂𝑢^∗₂²>0,𝑏=𝑐₂𝑢^∗₂𝑑₁+2𝑎₁₁𝑢^∗₁+𝑎₁₂𝑢^∗₂+𝑏₁𝑢^∗₁𝑑₂+𝑎₂₁𝑢^∗₁+2𝑎₂₂𝑢^∗₂+𝑎₁₂𝑏₂𝑢^∗₁𝑢^∗₂−𝑎₂₁𝑐₁𝑢^∗₁𝑢^∗₂,𝑐=𝑏₁𝑐₂𝑢^∗₁𝑢^∗₂+𝑏₂𝑐₁𝑢^∗₁𝑢^∗₂>0. (16) Let 𝜆₁(𝜇_𝑖) and 𝜆₂(𝜇_𝑖) be the two roots of (13); then we have 𝜆₁𝜇_𝑖+𝜆₂𝜇_𝑖=𝐺₁₁+𝐺₂₂−𝜇_𝑖Φ₁₁+Φ₂₂,𝜆₁𝜇_𝑖𝜆₂𝜇_𝑖=det𝐆_𝐮𝐮^∗−𝜇_𝑖Φ_𝐮. (17) In order to ensure that Re𝜆₁(𝜇_𝑖)<0 and Re𝜆₂(𝜇_𝑖)>0, a sufficient condition of Turing instability about homogeneous equilibrium 𝐮^∗ is det𝐆_𝐮𝐮^∗−𝜇_𝑖Φ_𝐮<0. (18)
Next, we look for the diffusion conditions such that det(𝐆_𝐮(𝐮^∗)−𝜇_𝑖Φ_𝐮)<0 holds. Furthermore, we have lim_𝑎21→∞𝑎𝑎₂₁=𝑑₁+2𝑎₁₁𝑢^∗₁𝑢^∗₁,lim_𝑎21→∞𝑏𝑎₂₁=𝑏₁𝑢^∗₁−𝑐₁𝑢^∗₂𝑢^∗₁,lim_𝑎21→∞𝑐𝑎₂₁=0. (19) Taking condition [𝐶2] in consideration, we can obtain that lim_𝑎21→∞det𝐆_𝐮𝐮^∗−𝜇_𝑖Φ_𝐮𝑎₂₁=𝑑₁+2𝑎₁₁𝑢^∗₁𝑢^∗₁𝜇²_𝑖+𝑏₁𝑢^∗₁−𝑐₁𝑢^∗₂𝑢^∗₁𝜇_𝑖∶=𝑃𝜇_𝑖. (20) Then, 𝑃(𝜇_𝑖)=0 has two real roots, one being 0 and the other being positive. A continuity argument shows that there exists a positive constant 𝑎^∗₂₁ such that when 𝑎₂₁>𝑎^∗₂₁, (15) and 𝜇_𝑖 axis intersect on two real and positive points, denoted by 𝜇_𝑖1, 𝜇_𝑖2(𝜇_𝑖1<𝜇_𝑖2). Hence, there exists a 𝜇_𝑖∈(𝜇_𝑖1,𝜇_𝑖2) such that det(𝐆_𝐮(𝐮^∗)−𝜇_𝑖Φ_𝐮)<0.

Remark 4. Theorem 3 is available for a case of the system (1) equipped with cross-diffusion and self-diffusion; that is, 𝑎₁₁≠0, 𝑎₁₂≠0 and 𝑎₂₂≠0. When 𝑎₁₁=𝑎₁₂=𝑎₂₂=0 and 𝑎₂₁≠0, the system (1) possesses a diffusion term the same as in [22].

Corollary 5. If 𝑎₂₁=0, then the homogeneous steady state 𝐮^∗ of the system (1) is always stable.

4. Turing Parameter Space

In this section, we will find some parameter regions of nonlinear diffusion coefficients where the equilibrium state 𝐮^∗ is unstable. For this, according to (15) the sufficient condition of Turing instability is 𝑏<0 and 𝑏²−4𝑎𝑐>0 besides [𝐶1] and [𝐶2]. In this paper, the parameter values satisfying conditions [𝐶1] and [𝐶2] are taken as follows: 𝑎₁=1,𝑎₂=2,𝑏₁𝑏=0.08,₂=0.05,𝑐₁=0.05,𝑐₂𝑑=0.15,₁=𝑑₂=0.1. (21) Then, for these fixed parameters the homogeneous steady state 𝐮^∗ is given by (𝑢^∗₁,𝑢^∗₂)=(3.4483,14.4828). In Figure 1, we examine the parameter regions where the homogeneous steady state 𝐮^∗ is expected to be unstable. These charts are obtained by fixed parameters in (21) as well as 𝑎₁₂=𝑎₂₂=0.01 for Figure 1(a), 𝑎₁₁=𝑎₂₂=0.01 for Figure 1(b), and 𝑎₁₁=𝑎₁₂=0.01 for Figure 1(c).

(a)

(b)

(c)

Figure 1 

Turing parameter spaces for the system (1). Inside the gray zones, 𝐮∗ is unstable, while, for parameters outside these zones, 𝐮∗ remains stable. These graphics are obtained by fixing parameters in (21); besides 𝑎12=𝑎22=0.01 (a), 𝑎11=𝑎22=0.01 (b), and 𝑎11=𝑎12=0.01 (c).

From the mathematical viewpoint, the Turing bifurcation occurs when for the characteristic root of (13), Im(𝜆(𝜇))=0 and Re(𝜆(𝜇))=0 at 𝜇=𝜇_𝑐≠0. Next, we will look for the critical wave of spatial patterns and note the relationship of 𝜇 and the wave number 𝑘; that is, 𝜇=𝑘² [26]. Thus, we only need to confirm that min_𝑘det𝐆_𝐮𝐮^∗−𝑘²Φ_𝐮=0. (22)

Then, the Turing bifurcation thresholds of parameters satisfy the following equation: 4det𝐆_𝐮𝐮^∗detΦ_𝐮𝐮^∗=Φ₁₁𝐆₂₂+Φ₂₂𝐆₁₁−Φ₂₁𝐆₁₂−Φ₁₂𝐆₂₁², (23) and the critical wavenumber 𝑘_𝑐 satisfies 𝑘²_𝑐=Φ₁₁𝐆₂₂+Φ₂₂𝐆₁₁−Φ₂₁𝐆₁₂−Φ₁₂𝐆₂₁2detΦ_𝐮(𝐮^∗). (24) To well see the effect of the nonlinear diffusion, according to Turing parameter regions in Figure 1 we plot the dispersion relations in Figure 2. The critical parameter values in Figure 2(a) correspond to 𝑎₁₁=0.201, 𝑎₁₂=0.353 in Figure 2(b), and 𝑎₂₂=0.558 in Figure 2(c). In addition, we find that the lowest limit of wavenumber 𝑘 corresponding to the available Turing modes Re(𝜆)>0 turns small with 𝑎₁₁ increasing in Figure 2(a), with 𝑎₁₂ increasing in Figure 2(b), and with 𝑎₂₂ increasing in Figure 2(c).

(a)

(b)

(c)

Figure 2 

Dispersion relations for different parameters: in (a), curves of the real part of the eigenvalue are derived by taking 𝑎21=10, 𝑎12=𝑎22=0.01, and 𝑎11=0.05,0.1,0.15,0.18,0.201,0.25 from top to bottom; in (b) 𝑎21=10, 𝑎11=𝑎22=0.01, and 𝑎12=0.15,0.2,0.25,0.3,0.353,0.4 from top to bottom; and in (c) 𝑎21=5.06, 𝑎11=𝑎12=0.01, and 𝑎22=0.35,0.4,0.45,0.508,0.558,0.608 from top to bottom. Other parameters are fixed as in (21).

5. Pattern Formation

In this section, using numerical methods, we perform numerical simulations of the system (1) in a two-dimensional space and illustrate that cross-diffusion and self-diffusion induce spatial patterns. Throughout this section, we assume that the region of the system (1) is [0,𝐿]×[0,𝐿]⊂ℝ². Hence, according to the definition of the eigenvalue 𝜇_𝑖 in Notation 2, we can obtain 𝜇_𝑖=((𝑖−1)𝜋/𝐿)². To numerically solve partial differential equations, we first have to discretize the space-time of the system (1). The region of [0,𝐿]×[0,𝐿] is solved in a discrete domain with 𝑀×𝑁 lattice sites. The length of the lattices is defined by a constant ℎ. The time is also discrete by a constant step 𝜏. All our numerical simulations employ the Neumann boundary condition. Here, we use the standard five-point approximation for the two-dimensional Laplacian derivative and the time evolution is solved by using the Euler method. More precisely, the value (𝑢^𝑛+1_1,𝑖,𝑗,𝑢^𝑛+1_2,𝑖,𝑗) at the time (𝑛+1)𝜏 at the mesh position (𝑥_𝑖,𝑦_𝑖) is obtained by 𝑢^𝑛+1_1,𝑖,𝑗=𝑢^𝑛_1,𝑖,𝑗+𝜏ΔΦ₁+𝜏𝐺₁𝑢^𝑛_1,𝑖,𝑗,𝑢^𝑛_2,𝑖,𝑗,𝑢^𝑛+1_2,𝑖,𝑗=𝑢^𝑛_2,𝑖,𝑗+𝜏ΔΦ₂+𝜏𝐺₂𝑢^𝑛_1,𝑖,𝑗,𝑢^𝑛_2,𝑖,𝑗, (25) with Laplacian defined by ΔΦ₁=Φ₁𝑢_{1,𝑖+1,𝑗},𝑢_{2,𝑖+1,𝑗}+Φ₁𝑢_{1,𝑖−1,𝑗},𝑢_{2,𝑖−1,𝑗}ℎ²+Φ₁𝑢_{1,𝑖,𝑗+1},𝑢_{2,𝑖,𝑗+1}+Φ₁𝑢_{1,𝑖,𝑗−1},𝑢_{2,𝑖,𝑗−1}ℎ²−4Φ₁𝑢_1,𝑖,𝑗,𝑢_2,𝑖,𝑗ℎ². (26) In this paper, we set ℎ=1, 𝜏=0.001, and 𝑀=𝑁=100. The initial data of the system (1) is taken as a uniformly distributed random perturbation in order of 1×10⁻⁴ around the homogeneous equilibrium state 𝐮^∗. More precisely, 𝑢₁(𝑥,0)=𝑢^∗₁+𝜂₁(𝑥),𝑢₂(𝑥,0)=𝑢^∗₂+𝜂₂(𝑥), (27) where 𝜂₁,𝜂₂∈[−1×10⁻⁴,1×10⁻⁴]. We simulate different patterns according to the dispersion relations in Figure 2.

In Figure 3, we show the evolution of the spatial patterns of the prey and the predators at 2×10³, 4×10⁴, and 1×10⁶ iterations for 𝑎₁₁=0.05,0.1 when we set 𝑎₂₁=10, 𝑎₁₂=0.01, and 𝑎₂₂=0.01. One can see that the patterns arise from random initial conditions. After the cold spot patterns for the prey and the hot spot patterns for the predator arise, they turn steadily with time until these patterns are temporally independent. In addition, for 𝑎₁₁=0.1 the cold spot patterns for the prey and the hot spot patterns for the predator are looser compared with those for 𝑎₁₁=0.05.

Figure 3 

Snapshots of contour pictures of time evolution of 𝑢1 and 𝑢2. From top to bottom: patterns of 𝑢1,𝑢2 for 𝑎11=0.05 and of 𝑢1,𝑢2 for 𝑎11=0.1. From left to right: 2×103 iterations, 4×104 iterations, and 1×106 iterations. Other diffusion parameters: 𝑎12=𝑎22=0.01 and 𝑎21=10.

In Figure 4, we fix 𝑎₁₁=0.01, 𝑎₂₂=0.01, and 𝑎₂₁=10 and obtain the spatial patterns of species 𝑢₁,𝑢₂ of time evolution for 𝑎₁₂=0.15 and 𝑎₁₂=0.3. For the case of 𝑎₁₂=0.15, the random initial distribution leads to the formation of irregular patterns. After a long time evolution, we find that the cold spot-strip patterns emerge for the prey 𝑢₁ and that the hot spot-strip patterns for the predator 𝑢₂ arise. However, in the case of 𝑎₁₂=0.3, the steady patterns of the prey 𝑢₁ consist of hot spots in a bigger size, while the steady patterns of the predator 𝑢₂ are in the formation of bigger cold spots.

Figure 4 

Snapshots of contour pictures of time evolution of 𝑢1 and 𝑢2. From top to bottom: patterns of 𝑢1,𝑢2 for 𝑎12=0.15 and of 𝑢1,𝑢2 for 𝑎12=0.3. From left to right: 2×103 iterations, 4×104 iterations, and 1×106 iterations. Other diffusion parameters: 𝑎11=𝑎22=0.01 and 𝑎21=10.

In Figure 5, diffusion parameters are set as 𝑎₁₁=𝑎₁₂=0 and 𝑎₂₁=5.06. We plot the patterns for 𝑎₂₂=0.35 and 𝑎₂₂=0.508, respectively. For both cases, one can see that as time goes on, the cold spot patterns of the prey 𝑢₁ and the hot spot patterns of the predator 𝑢₂ ultimately form.

Figure 5 

Snapshots of contour pictures of time evolution of 𝑢1 and 𝑢2. From top to bottom: patterns of 𝑢1,𝑢2 for 𝑎22=0.35 and of 𝑢1,𝑢2 for 𝑎22=0.508. From left to right: 2×103 iterations, 4×104 iterations, and 1×106 iterations. Other diffusion parameters: 𝑎11=𝑎12=0.01 and 𝑎21=5.06.

6. Conclusion and Discussion

In this paper, we have studied the prey-predator model with the nonlinear diffusions including normal diffusion, cross-diffusion, and self-diffusion. By applying the mathematical analysis and suitable numerical simulations, we obtain the sufficient conditions of the formation of Turing patterns for this nonlinear diffusion and illustrate Turing parameter regions and Turing patterns when some parameters in system (1) are set.

In our results, we have provided Theorems 1 and 3 to demonstrate that for the nonlinear diffusion including self-diffusion and the cross-diffusion of the predator, the parameter 𝑎₂₁ plays an important role to induce Turing instability. Furthermore, if 𝑎₂₁=0, then the homogeneous equilibrium state 𝐮^∗ is always stable; that is, the system (1) has no Turing patterns. By performing numerical simulations, we find the Turing parameter regions of the interaction between cross-diffusion 𝑎₂₁ and other diffusion terms including cross-diffusion of the other species and self-diffusions. Besides, according to these parameter regions, we show the corresponding dispersion relations and the corresponding patterns. These results indicate that Turing patterns can emerge through the interaction between the cross-diffusion 𝑎₂₁ and self-diffusion as well as other cross-diffusions in the system (1).

It is well known that for a prey-predator system, the formation of patterns can occur by introducing the cross-diffusion. However, our results further show that under condition (see Theorem 3), self-diffusion can produce Turing patterns.

Acknowledgments

The authors would like to thank the referees for their careful reading and constructive comments to improve their paper and writing. This work is supported by Doctoral Initial Foundation of Nanchang University, China (no. 06301004), and the Youth Natural Science Foundation of China (no. 61304161).