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

Authors: Hongwei Yin; Xiaoyong Xiao; Xiaoqing Wen
Journal: Abstract and Applied Analysis
 (2013)

---

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


---

## Body

## 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 Section2, 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]:
(1)∂u1∂t-Δ[(d1+a11u1+a12u2)u1]=u1(a1-b1u1-c1u2),inΩ×[0,∞),∂u2∂t-Δ[(d2+a21u1+a22u2)u2]=u2(a2+b2u1-c2u2),inΩ×[0,∞),∂u1∂n=∂u2∂n=0,on∂Ω×[0,∞),u1(x,0)=u10(x)≥0,u2(x,0)=u20(x)≥0,inΩ,
where Ω is a bounded domain in ℝ2 with smooth boundary ∂Ω and represents the domain that these two species inhabit. The vector n is the outward unit normal vector of Ω. In this model, u1 and u2 are the prey and the predator densities, respectively; ai, bi, and ci (i=1,2) are positive constants. a1,a2 are rates of the prey and predator proliferations for food source; a1/b1 and a2/c2 are environmental carrying capacities for prey and predator, respectively. c1 is a consumption rate; b2 is a conversion rate. In the diffusion terms, the constant di(i=1,2), which is usually termed as a normal-diffusion coefficient, represents the natural dispersive force of movement of a species. The positive constants a12 and a21 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 a11 and a22 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.
Let0=μ1<μ2<⋯→∞ be the eigenvalues of -Δ on Ω under no-flux boundary condition, and let E(μi) be the space of eigenfunction corresponding to μi. We define the following space decomposition:(i)
Xij≔{φijc:c∈ℝ2}, where φij are orthonormal bases of E(μi) for j=1,…,dimE(μi);(ii)
X≔{u∈[𝒞1(Ω¯)]2:∂u1/∂n=∂u2/∂n=0on∂Ω}, and thus X=⊕i=1∞Xi, where Xi=⊕j=1dimE(μi)Xij.Notation 2.
For the sake of simplicity, we denote reaction terms for systems (1) by
(2)G(u)=(G1(u)G2(u))=(u1g1(u)≔u1(a1-b1u1-c1u2)u2g2(u)≔u2(a2+b2u1-c2u2)),
and its Jacobian matrix at  the point u*  is
(3)Gu(u*)=(G11G12G21G22)≔(-b1u1*-c1u1*b2u2*-c2u2*).
The diffusion term is denoted as
(4)Φ(u)=(Φ1(u)Φ2(u))≔((d1+a11u1+a12u2)u1(d2+a21u1+a22u2)u2),
and its corresponding Jacobian matrix at the point u* is
(5)Φu(u*)=(Φ11Φ12Φ21Φ22)≔(d1+2a11u1*+a12u2*a12u1*a21u2*d2+a21u1*+2a22u2*).


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

Letaij=0 for i,j=1,2; then the system (1) degenerates into standard reaction-diffusion equations:
(6)∂u1∂t-d1Δu1=u1(a1-b1u1-c1u2),inΩ×[0,∞),∂u2∂t-d2Δu2=u2(a2+b2u1-c2u2),inΩ×[0,∞),∂u1∂n=∂u2∂n=0,on∂Ω×[0,∞),u1(x,0)=u10(x)≥0,u2(x,0)=u20(x)≥0,inΩ.
Through this paper, we assume the following conditions:
(7)[C1]c2c1>a2a1,[C2]b2b1>c2c1-2a2a1.Then, there exists a unique positive equilibrium state for system (6), denoted by u*=(u1*,u2*), where
(8)u1*=a1c2-a2c1b1c2+b2c1,u2*=a1b2+a2b1b1c2+b2c1.Theorem 1.
If there are no cross-diffusion and self-diffusion, the positive equilibrium stateu* of the system (6) is locally asymptotically stable when condition [C1] holds.Proof.
The  linearization of (6) around the steady state u* can be therefore expressed by
(9)∂u∂t=(DΔ+Gu(u*))u,
where D=diag(d1,d2). Obviously, the operator DΔ+Gu(u*) is invariant in the subspace Xi, and λi is an eigenvalue of this operator on Xi, if and only if it is an eigenvalue of -μiD+Gu(u*). Direct calculation shows the characteristic equation
(10)λi2-miλi+ni=0,
where mi=-u1*b1-u2*c2-(d1+d2)μi and ni=(u1*b1+d1μi)(u2*c2+d2μi)+u1*u2*b2c1. It is easy to be verified that mi is negative and ni positive. Thus, for each i>1, the two roots of (10) have negative real parts. Consequently, we complete this proof.Remark 2.
According to the proof of Theorem1, we can also calculate out the eigenvector c corresponding to the eigenvalue λi for the operator -μiD+Gu(u*), which satisfies
(11)(λi+μiD-Gu(u*))c=0.
Furthermore, we can yield corresponding eigenspace {φijc} in Xi.By Theorem1, under condition [C1] system (6) cannot destabilize u*. Next, the cross-diffusion and self-diffusion are taken into account.Theorem 3.
Assume that conditions[C1] and [C2] hold; then there exists a  sufficiently large positive constant a21* such that the homogeneous state u* of the system (1) is unstable provided that a21>a21*.Proof.
We first linearize the system (1) around (u1,u2)=(u1*,u2*) and arrive at
(12)∂u∂t=(Φu(u*)Δ+Gu(u*))u.
By the same method as the proof of Theorem 1, we consider the operator Φu(u*)Δ+Gu(u*) on the subspace Xi. The eigenvalue of this operator on Xi is denoted as λ, and then it is also an eigenvalue of the matrix -μiΦu(u*)+Gu(u*). Thus, λ satisfies the following equation:
(13)λ2-[G11+G22-μi(Φ11+Φ22)]λ+det(Gu(u*)-μiΦu)=0.
According to Notation 2, we can calculate and get
(14)G11+G22-μi(Φ11+Φ22)=-μi[d1+d2+(2a11+a21)u1*+(a12+2a12)u2*]-b1u1*-c2u2*<0,(15)det(Gu(u*)-μiΦu)=aμi2+bμi+c,
where
(16)a=d1d2+d1(a21u1*+2a22u2*)+d2(2a11u1*+a12u2*)+2a11a21(u1*)2+4a11a22u1*u2*+2a12a22(u2*)2>0,b=c2u2*(d1+2a11u1*+a12u2*)+b1u1*(d2+a21u1*+2a22u2*)+a12b2u1*u2*-a21c1u1*u2*,c=b1c2u1*u2*+b2c1u1*u2*>0.
Let λ1(μi) and λ2(μi) be the two roots of (13); then we have
(17)λ1(μi)+λ2(μi)=G11+G22-μi(Φ11+Φ22),λ1(μi)λ2(μi)=det(Gu(u*)-μiΦu).
In order to ensure that Reλ1(μi)<0 and Reλ2(μi)>0, a sufficient condition of Turing instability about homogeneous equilibrium u* is
(18)det(Gu(u*)-μiΦu)<0.
Next, we look for the diffusion conditions such thatdet(Gu(u*)-μiΦu)<0 holds. Furthermore, we have
(19)lima21→∞aa21=(d1+2a11u1*)u1*,lima21→∞ba21=(b1u1*-c1u2*)u1*,lima21→∞ca21=0.
Taking condition [C2] in consideration, we can obtain that
(20)lima21→∞det(Gu(u*)-μiΦu)a21=(d1+2a11u1*)u1*μi2+(b1u1*-c1u2*)u1*μi≔P(μi).
Then, P(μi)=0 has two real roots, one being 0 and the other being positive. A continuity argument shows that there exists a positive constant a21* such that when a21>a21*, (15) and μi axis intersect on two real and positive points, denoted by μi1, μi2(μi1<μi2). Hence, there exists a μi∈(μi1,μi2) such that det(Gu(u*)-μiΦu)<0.Remark 4.
Theorem3 is available for a case of the system (1) equipped with cross-diffusion and self-diffusion; that is, a11≠0, a12≠0 and a22≠0. When a11=a12=a22=0 and a21≠0, the system (1) possesses a diffusion term the same as in [22].Corollary 5.
Ifa21=0, then the homogeneous steady state u* 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 stateu* is unstable. For this, according to (15) the sufficient condition of Turing instability is b<0 and b2-4ac>0 besides [C1] and [C2]. In this paper, the parameter values satisfying conditions [C1] and [C2] are taken as follows:
(21)a1=1,a2=2,b1=0.08,b2=0.05,c1=0.05,c2=0.15,d1=d2=0.1.
Then, for these fixed parameters the homogeneous steady state u* is given by (u1*,u2*)=(3.4483,14.4828). In Figure 1, we examine the parameter regions where the homogeneous steady state u* is expected to be  unstable. These charts are obtained by fixed parameters in (21) as well as a12=a22=0.01 for Figure 1(a), a11=a22=0.01 for Figure 1(b), and a11=a12=0.01 for Figure 1(c).Turing parameter spaces for the system (1). Inside the gray zones, u* is  unstable, while, for parameters outside these zones, u* remains stable. These graphics are obtained by fixing parameters in (21); besides a12=a22=0.01 (a), a11=a22=0.01 (b), and a11=a12=0.01 (c).
(a)
(b)
(c)From the mathematical viewpoint, the Turing bifurcation occurs when for the characteristic root of (13), Im(λ(μ))=0 and Re(λ(μ))=0 at μ=μc≠0. Next, we will look for the critical wave of spatial patterns and note the relationship of μ and the wave number k; that is, μ=k2 [26]. Thus, we only need to confirm that
(22)minkdet(Gu(u*)-k2Φu)=0.Then, the Turing bifurcation thresholds of parameters satisfy the following equation:(23)4det(Gu(u*))det(Φu(u*))=(Φ11G22+Φ22G11-Φ21G12-Φ12G21)2,
and the critical wavenumber kc satisfies
(24)kc2=(Φ11G22+Φ22G11-Φ21G12-Φ12G21)(2det(Φu(u*))).
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 a11=0.201, a12=0.353 in Figure 2(b), and a22=0.558 in Figure 2(c). In addition, we find that the lowest limit of wavenumber k corresponding to the available Turing modes Re(λ)>0 turns small with a11 increasing in Figure 2(a), with a12 increasing in Figure 2(b), and with a22 increasing in Figure 2(c).Dispersion relations for different parameters: in (a), curves of the real part of the eigenvalue are derived by takinga21=10, a12=a22=0.01, and a11=0.05,0.1,0.15,0.18,0.201,0.25 from top to bottom; in (b)  a21=10, a11=a22=0.01, and a12=0.15,0.2,0.25,0.3,0.353,0.4 from top to bottom; and in (c)  a21=5.06, a11=a12=0.01, and a22=0.35,0.4,0.45,0.508,0.558,0.608 from top to bottom. Other parameters are fixed as in (21).
(a)
(b)
(c)


## 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,L]×[0,L]⊂ℝ2. Hence, according to the definition of the eigenvalue μi in Notation 2, we can obtain μi=((i-1)π/L)2. To numerically solve partial differential equations,  we first have to  discretize the space-time of the system (1). The region of [0,L]×[0,L] is solved in a discrete domain with M×N lattice sites. The length of the lattices is defined by a constant h. 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 (u1,i,jn+1,u2,i,jn+1) at the time (n+1)τ at the mesh position (xi,yi) is obtained by
(25)u1,i,jn+1=u1,i,jn+τΔΦ1+τG1(u1,i,jn,u2,i,jn),u2,i,jn+1=u2,i,jn+τΔΦ2+τG2(u1,i,jn,u2,i,jn),
with Laplacian defined by
(26)ΔΦ1=Φ1(u1,i+1,j,u2,i+1,j)+Φ1(u1,i-1,j,u2,i-1,j)h2+Φ1(u1,i,j+1,u2,i,j+1)+Φ1(u1,i,j-1,u2,i,j-1)h2-4Φ1(u1,i,j,u2,i,j)h2.
In this paper, we set h=1, τ=0.001, and M=N=100. The initial data of the system (1) is taken as a uniformly distributed random perturbation in order of 1×10-4 around the homogeneous equilibrium state u*. More precisely,
(27)u1(x,0)=u1*+η1(x),u2(x,0)=u2*+η2(x),
where η1,η2∈[-1×10-4,1×10-4]. We simulate different patterns according to the dispersion relations in Figure 2.In Figure3, we show the evolution of the spatial patterns of the prey and the predators at 2×103, 4×104, and 1×106 iterations for a11=0.05,0.1 when we set a21=10, a12=0.01, and a22=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 a11=0.1 the cold spot patterns for the prey and the hot spot patterns for the predator are looser compared with those for a11=0.05.Figure 3
Snapshots of contour pictures of time evolution ofu1 and u2. From top to bottom: patterns of u1,u2 for a11=0.05 and of u1,u2 for a11=0.1. From left to right: 2×103 iterations, 4×104 iterations, and 1×106 iterations. Other diffusion parameters: a12=a22=0.01 and a21=10.In Figure4, we fix a11=0.01, a22=0.01, and a21=10 and obtain the spatial patterns of species u1,u2 of time evolution for a12=0.15 and a12=0.3. For the case of a12=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 u1 and that the hot spot-strip patterns for the predator u2 arise. However, in the case of a12=0.3, the steady patterns of the prey u1 consist of hot spots in a bigger size, while the steady patterns of the predator u2 are in the formation of bigger cold spots.Figure 4
Snapshots of contour pictures of time evolution ofu1 and u2. From top to bottom: patterns of u1,u2 for a12=0.15 and of u1,u2 for a12=0.3. From left to right: 2×103 iterations, 4×104 iterations, and 1×106 iterations. Other diffusion parameters: a11=a22=0.01 and a21=10.In Figure5, diffusion parameters are set as a11=a12=0 and a21=5.06. We plot the patterns for a22=0.35 and a22=0.508, respectively. For both cases, one can see that as time goes on, the cold spot patterns of the prey u1 and the hot spot patterns of the predator u2 ultimately form.Figure 5
Snapshots of contour pictures of time evolution ofu1 and u2. From top to bottom: patterns of u1,u2 for a22=0.35 and of u1,u2 for a22=0.508. From left to right: 2×103 iterations, 4×104 iterations, and 1×106 iterations. Other diffusion parameters: a11=a12=0.01 and a21=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 Theorems1 and 3 to demonstrate that for the nonlinear diffusion including self-diffusion and the cross-diffusion of the predator, the parameter a21 plays an important role to induce Turing instability. Furthermore, if a21=0, then the homogeneous equilibrium state u* 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 a21 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 a21 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 Theorem3), self-diffusion can produce Turing patterns.



---

*Source: 891738-2013-11-12.xml*