# Cross-Diffusion-Driven Instability in a Reaction-Diffusion Harrison Predator-Prey Model

Authors: Xiaoqin Wang; Yongli Cai
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We present a theoretical analysis of processes of pattern formation that involves organisms distribution and their interaction of spatially distributed population with cross-diffusion in a Harrison-type predator-prey model. We analyze the global behaviour of the model by establishing a Lyapunov function. We carry out the analytical study in detail and find out the certain conditions for Turing’s instability induced by cross-diffusion. And the numerical results reveal that, on increasing the value of the half capturing saturation constant, the sequences “spots → spot-stripe mixtures → stripes → hole-stripe mixtures → holes” are observed. The results show that the model dynamics exhibits complex pattern replication controlled by the cross-diffusion.


---

## Body

## 1. Introduction

Understanding of spatial and temporal behaviors of interacting species in ecological systems is one of the central scientific problems in population ecology [1–15], since the pioneering work of Turing [16]. Throughout the history of theoretical ecology, reaction-diffusion equations have been intensively used to describe spatiotemporal dynamics [15, 17].In recent years, the effect of cross-diffusion in reaction-diffusion systems has received much attention by both ecologists and mathematicians, for example, see [18–25] and the references therein. Kerner [18] was the first to examine that cross-diffusion can induce pattern forming instability in an ecological situation. Cross-diffusion expresses the population fluxes of one species due to the presence of the other species. And Gurtin [19] developed some mathematical models for population dynamics with the inclusion of cross-diffusion as well as self-diffusion and showed that the effect of cross-diffusion may give rise to the segregation of two species.In this paper, we are attempting to study the effect of cross-diffusion in a predator-prey model with Harrison-type functional response [26]. The model can be written as
(1)∂u∂t=ru(1-uK)-c1uvm1v+1+D11Δu+D12Δv,x=(ξ,η)∈Ω,t>0,∂v∂t=v(-h1+b1um1v+1)+D21Δu+D22Δv,x=(ξ,η)∈Ω,t>0,u(x,0)=u0(x),v(x,0)=v0(x),x=(ξ,η)∈Ω,
where u and v represent population density of prey and predator at time t, respectively. r is the intrinsic growth rate of prey, K is the prey carrying capacity, c1 is the capture rate, m1 the half capturing saturation constant, h1 is the death rate of predator, and b1 is conversion rate. D11 and D22 are the self-diffusion coefficients of u and v, respectively, D12 and D21 are the cross-diffusion coefficients of u and v, respectively. We always assume that D11>0, D22>0 and D11D22-D12D21>0. The value of the cross-diffusion coefficient may be positive, negative, or zero. Positive cross-diffusion coefficient denotes, that one species tends to move in the direction of lower concentration of another species, while negative cross-diffusion expresses the population fluxes of one species in the direction of higher concentration of the other species [27]. Δ=∂2/∂x2=∂2/∂ξ2+∂2/∂η2 is the usual Laplacian operator in 2-dimensional space. Ω⊂ℝ2 is a bounded domain with smooth boundary ∂Ω. The initial data u0(x) and v0(x) are continuous functions on Ω¯.We make a change of variables:(2)(u,v,t)=(Ku~,Kv~,t~r).
For the sake of convenience, we still use variables u, v instead of u~, v~. Thus, considering zero-flux boundary conditions, model (1) is converted into
(3)∂u∂t=u(1-u)-cuvmv+1+d11Δu+d12Δv,x=(ξ,η)∈Ω,t>0,∂v∂t=v(-h+bumv+1)+d21Δu+d22Δv,x=(ξ,η)∈Ω,t>0,∂u∂ν=∂v∂ν=0,x=(ξ,η)∈∂Ω,t>0,u(x,0)=u0(x)≥0,v(x,0)=v0(x)≥0,x=(ξ,η)∈Ω,
where the new parameters are
(4)c=Kc1r,m=Km1,h=h1r,b=b1Kr,d11=D11r,d12=D12r,d21=D21r,d22=D22r.D=(d11d12d21d22) is the diffusion matrix, d11>0, d22>0, and det(D)=d11d22-d12d21>0. ν is the outward unit normal vector on ∂Ω and the zero-flux boundary conditions mean that model (3) is self-contained and has no population flux across the boundary ∂Ω [28, 29].In particular, whend12=d21=0, that is, the cross-diffusion coefficients are equal 0, we can obtain the following model:
(5)∂u∂t=u(1-u)-cuvmv+1+d11Δu,x=(ξ,η)∈Ω,t>0,∂v∂t=v(-h+bumv+1)+d22Δv,x=(ξ,η)∈Ω,t>0,∂u∂ν=∂v∂ν=0,x=(ξ,η)∈∂Ω,t>0,u(x,0)=u0(x)≥0,v(x,0)=v0(x)≥0,x=(ξ,η)∈Ω.
We call model (5) as self-diffusion model, while we call model (3) cross-diffusion model.The corresponding kinetic equation to models (3) and (5) is:
(6)u˙=u(1-u)-cuvmv+1≜f(u,v),v˙=v(-h+bumv+1)≜g(u,v).In recent years, there has been considerable interest to investigate the stability behavior of a predator-prey system by taking into account the effect of self- as well as cross-diffusion [3, 6, 8–10, 12–15]. But in the studies on the spatiotemporal dynamics of predator-prey system with functional response, little attention has been paid to study on the effect of cross-diffusion.Mathematically speaking, an equilibrium in Turing’s instability (diffusion-driven instability) means that it is an asymptotically stable equilibriumE* of model (6) but is unstable with respect to the solutions of reaction-diffusion model (3) or (5). Especially, if E* is also stable with respect to the solutions of the self-diffusion model (3), that is, d12=d21=0 in the cross-diffusion model (5), then there is nonexistence of Turing’s instability in this situation.And there comes a question: if there is nonexistence of Turing’s instability in the case of self-diffusion (i.e.,d12=d21=0), does model (5) exhibit Turing’s instability induced by cross-diffusion?The main purpose of this paper is to focus on the effect of cross-diffusion on the spatiotemporal dynamics of the reaction-diffusion predation model. The paper is organized as follows. In Section2, we give some properties of the solutions of the model. In Section 3, we give the linearized stability analysis to show (i.e., no cross-diffusion), deduce the conditions of Turing’s instability induced by cross-diffusion, and illustrate the different Turing patterns by using the numerical simulations. Finally, in Section 4, some conclusions and discussions are given.


## 2. Dynamics Analysis

In this section, we present some preliminary results, including dissipativeness, boundedness, permanence of the solutions, and the equilibria stability analysis of the models.


### 2.1. Dissipativeness

Theorem 1.
For any solution(u,v) of model (6),
(7)limsupt→∞u(t)≤1,limsupt→∞v(t)≤max{0,b-hhm}.
Hence, model (6) is dissipative.Proof.
From the first equation of model (6), it can be easily shown that
(8)limsupt→∞u(t)≤1.
If b>h, from the second equation of model (6), one has:
(9)v˙≤v(b-h-hmv).
A standard comparison argument shows that
(10)limsupt→∞v(t)≤b-hhm≜α.
If b≤h, we have the following differential inequality:
(11)v˙≤-hmv2mv+1,
and the same argument above yields
(12)limsupt→∞v(t)≤0.
In either case, the second inequality of (7) holds.


### 2.2. Boundedness

Theorem 2.
All the solutions of model (6) which initiate in ℝ2+ are uniformly bounded within the region Γ, where
(13)Γ={(u,v):0≤u+cbv≤1+14h}.Proof.
Let us define the function:(14)w(t)=u(t)+cbv(t).
Calculating the time derivative of w(t) along the trajectories of model (6), we get
(15)w˙(t)=u˙+v˙=u(1-u)-chbv.
Then,
(16)w˙(t)+hw(t)=u(1-u)+hu<14+h.
Using the theory of differential inequality, for all t≥T≥0, we have
(17)0≤w(t)≤1+14h-(1+14h-w(T))e-(t-T).
Hence, we have
(18)limsupt→∞w(t)≤1+14h.
Hence, all the solutions of model (6) that initiate in ℝ2+ are confined in the region Γ.


### 2.3. Permanence

Theorem 3.
Ifc<m and h<(b/2)(1-c/m), then model (6) has the permanence property.Proof.
By the first equation of model (6), we have
(19)u˙=u(1-u)-cuvmv+1=u(1-u)-cu(1/m)(mv+1)-(1/m)mv+1≥u(1-c/m-u).
Since c<m, by the famous comparison theorem, we have
(20)liminft→∞u(t)≥1-cm>0.
Hence, for large t, u(t)>(1/2)(1-c/m)≜η.
As a result, for larget, v satisfies
(21)v˙≥v(-h+bηmv+1)=v(bη-h-hmv)mv+1.
Since h<(b/2)(1-c/m), by the famous comparison theorem, we can get
(22)liminft→∞v(t)≥b(m-c)-2hm2hm2≜β>0.
The proof is complete.


### 2.4. Stability Analysis of the Equilibria

The nonspatial model (6) has three equilibria, which correspond to spatially homogeneous equilibria of model (3) and model (5), in the positive quadrant:(i)
E0=(0,0) (total extinct) is a saddle point;(ii)
E1=(1,0) (extinct of the predator, or prey only) is a saddle when b>h, or stable node when b<h;(iii)
E3=(u*,v*) (coexistence of prey and predator), where u*=h(mv*+1)/b, and v* satisfies
(23)hm2v2-(m(b-h)-hm-bc)v+h-b=0.It is easy to verify that it has a unique positive equilibrium if b>h.The Jacobian matrix for the positive equilibriumE3=(u*,v*) is given by
(24)J=(-h(mv*+1)b-chb(mv*+1)bv*mv*+1-hmv*mv*+1)≜(J11J12J21J22).
Obviously,
(25)det(J)=hv*(hm3v*2+2hm2v*+hm+bc)b(mv*+1)2>0,tr(J)=-h(m2v*2+m(2+b)v*+1)b(mv*+1)<0.Therefore, we can obtain the following.Theorem 4.
Assume that the positive equilibriumE3=(u*,v*) exists, then E3=(u*,v*) is locally stable for model (6).In the following, we shall prove that the positive equilibriumE3=(u*,v*) of model (3) is globally asymptotically stable.Theorem 5.
Suppose thath<b, c<m, and h<(b/2)(1-c/m). The positive equilibrium E3=(u*,v*) of model (3) is globally asymptotically stable, if,(a1)
h(mv*+1)/b-c/m+βc/(mv*+1)>0;(a2)
h(h(mv*+1)/b+βc/(mv*+1)-1)(1-1/(βm+1))-h2/4b2-αb2/4(mv*+1)2>0;(a3)
4d11d22>(d12+d21)2,
whereα=(b-h)/hm>0, β=(b(m-c)-2hm)/2hm2.Proof.
We adopt the Lyapunov function:(26)V(t)=∫Ω[V1(u(x,t))+V2(v(x,t))]dx,
where V1(u)=(1/2)(u-u*)2, V2(v)=(1/2)(v-v*)2.
Then,(27)dVdt=∫Ω(∂u∂t(u-u*)+∂v∂t(v-v*))dx=∫Ω((u-u*)(u-u2-cuvmv+1)+(v-v*)(-hv+buvmv+1))dx+∫Ω((u-u*)(d11Δu+d12Δv)+(v-v*)(d21Δu+d22Δv))dx=I1+I2,
where
(28)I1=∫Ω(u-u*)×(u-u2-cuvmv+1-u*+u*2+cu*v*mv*+1)dx+∫Ω(v-v*)(-hv+buvmv+1+hv*-bu*v*mv*+1)dx,I2=∫Ω(u-u*)(d11Δu+d12Δv-d11Δu*-d12Δv*)dx+∫Ω(v-v*)(d21Δu+d22Δv-d21Δu*-d22Δv*)dx.
By some computational analysis, we obtain(29)I1=-∫Ω(u-u*)2(u+u*+cvmv*+1-1)dx-∫Ω(v-v*)2(h-h(mv+1))dx-∫Ω(u-u*)(v-v*)u*-bv(mv+1)(mv*+1)(mv+1)dx.
Considering the zero-flux boundary conditions, we have
(30)I2=-d11∫Ω|∇u|2dx-d12∫Ω∇u∇vdx+2d11∫Ω∇u∇u*dx+d12∫Ω∇u∇v*dx+d12∫Ω∇v∇v*dx-d11∫Ω|∇u*|2dx-d12∫Ω∇u*∇v*dx-d22∫Ω|∇v|2dx-d21∫Ω∇u∇vdx+d21∫Ω∇v∇u*dx+2d22∫Ω∇v∇v*dx+d21∫Ω∇u∇v*dx-d21∫Ω∇u*∇v*dx-d22∫Ω|∇v*|2dx=-d11∫Ω|∇u|2dx-(d12+d21)∫Ω∇u∇vdx+2d11∫Ω∇u∇u*dx-d11∫Ω|∇u*|2dx+(d12+d21)∫Ω∇u∇v*dx-(d12+d21)∫Ω∇u*∇v*dx+d12∫Ω∇v∇v*dx+2d22∫Ω∇v∇v*dx-d22∫Ω|∇v|2dx-d22∫Ω|∇v*|2dx+d21∫Ω∇v∇u*dx.
I
1 and I2 in the curly brackets can be expressed in the form -XAXT and -YBYT, respectively, where
(31)X=(u-u*,v-v*),Y=(∇u,∇v,∇u*,∇v*),A=(u+u*+cvmv*+1-1u*-bv(mv+1)2(mv*+1)(mv+1)u*-bv(mv+1)2(mv*+1)(mv+1)h-hmv+1)≜(φ1φ2φ3φ4),B=(d11d12+d212-d11d12+d212d12+d212d22d12+d212-d22-d11d12+d212d11d12+d212d12+d212-d22d12+d212d22).dV/dt is negative definite if the symmetric matrices A and B are positive. It can be easily shown that the symmetric matrix B is positive definite if
(32)4d11d22>(d12+d21)2.
The symmetric matrixA is positive definite if
(33)φ1>0,φ4>0,Φ(u,v)=φ1φ4-φ2φ3>0.
Since
(34)φ1=u+u*+cvmv*+1-1>h(mv*+1)b-cm+βcmv*+1,
due to (a1), φ1>0 is true. It is easy to verify that φ4=h-(h/(mv+1))>0.
Since(35)Φ(u,v)=h(u+u*+cvmv*+1-1)×(1-1mv+1)-14(u*-bv(mv+1)(mv*+1)(mv+1))2=h(u+h(mv*+1)b+cvmv*+1-1)×(1-1mv+1)-14(hb(mv+1)-bvmv*+1)2,
then
(36)∂Φ(u,v)∂u=h(1-1mv+1)>0.
Hence, Φ(u,v) is strictly increasing in ℝ+, with respect to u, and
(37)Φ(0,v)=h(h(mv*+1)b+cvmv*+1-1)(1-1mv+1)-14(hb(mv+1)-bvmv*+1)2≥h(h(mv*+1)b+βcmv*+1-1)(1-1βm+1)-h24b2-αb24(mv*+1)2.
Consequently, if (a2) holds, Φ(0,v)>0. As a result, Φ(u,v)>Φ(0,v)>0.
Hence,V is a Lyapunov function and the positive equilibrium E3 of model (3) is globally asymptotically stable. This completes the proof.Remark 6.
Whend12=d21=0, Theorem 5 is true, too. That is, the positive equilibrium E3 of the self-diffusion model (5) is globally asymptotically stable.


## 2.1. Dissipativeness

Theorem 1.
For any solution(u,v) of model (6),
(7)limsupt→∞u(t)≤1,limsupt→∞v(t)≤max{0,b-hhm}.
Hence, model (6) is dissipative.Proof.
From the first equation of model (6), it can be easily shown that
(8)limsupt→∞u(t)≤1.
If b>h, from the second equation of model (6), one has:
(9)v˙≤v(b-h-hmv).
A standard comparison argument shows that
(10)limsupt→∞v(t)≤b-hhm≜α.
If b≤h, we have the following differential inequality:
(11)v˙≤-hmv2mv+1,
and the same argument above yields
(12)limsupt→∞v(t)≤0.
In either case, the second inequality of (7) holds.


## 2.2. Boundedness

Theorem 2.
All the solutions of model (6) which initiate in ℝ2+ are uniformly bounded within the region Γ, where
(13)Γ={(u,v):0≤u+cbv≤1+14h}.Proof.
Let us define the function:(14)w(t)=u(t)+cbv(t).
Calculating the time derivative of w(t) along the trajectories of model (6), we get
(15)w˙(t)=u˙+v˙=u(1-u)-chbv.
Then,
(16)w˙(t)+hw(t)=u(1-u)+hu<14+h.
Using the theory of differential inequality, for all t≥T≥0, we have
(17)0≤w(t)≤1+14h-(1+14h-w(T))e-(t-T).
Hence, we have
(18)limsupt→∞w(t)≤1+14h.
Hence, all the solutions of model (6) that initiate in ℝ2+ are confined in the region Γ.


## 2.3. Permanence

Theorem 3.
Ifc<m and h<(b/2)(1-c/m), then model (6) has the permanence property.Proof.
By the first equation of model (6), we have
(19)u˙=u(1-u)-cuvmv+1=u(1-u)-cu(1/m)(mv+1)-(1/m)mv+1≥u(1-c/m-u).
Since c<m, by the famous comparison theorem, we have
(20)liminft→∞u(t)≥1-cm>0.
Hence, for large t, u(t)>(1/2)(1-c/m)≜η.
As a result, for larget, v satisfies
(21)v˙≥v(-h+bηmv+1)=v(bη-h-hmv)mv+1.
Since h<(b/2)(1-c/m), by the famous comparison theorem, we can get
(22)liminft→∞v(t)≥b(m-c)-2hm2hm2≜β>0.
The proof is complete.


## 2.4. Stability Analysis of the Equilibria

The nonspatial model (6) has three equilibria, which correspond to spatially homogeneous equilibria of model (3) and model (5), in the positive quadrant:(i)
E0=(0,0) (total extinct) is a saddle point;(ii)
E1=(1,0) (extinct of the predator, or prey only) is a saddle when b>h, or stable node when b<h;(iii)
E3=(u*,v*) (coexistence of prey and predator), where u*=h(mv*+1)/b, and v* satisfies
(23)hm2v2-(m(b-h)-hm-bc)v+h-b=0.It is easy to verify that it has a unique positive equilibrium if b>h.The Jacobian matrix for the positive equilibriumE3=(u*,v*) is given by
(24)J=(-h(mv*+1)b-chb(mv*+1)bv*mv*+1-hmv*mv*+1)≜(J11J12J21J22).
Obviously,
(25)det(J)=hv*(hm3v*2+2hm2v*+hm+bc)b(mv*+1)2>0,tr(J)=-h(m2v*2+m(2+b)v*+1)b(mv*+1)<0.Therefore, we can obtain the following.Theorem 4.
Assume that the positive equilibriumE3=(u*,v*) exists, then E3=(u*,v*) is locally stable for model (6).In the following, we shall prove that the positive equilibriumE3=(u*,v*) of model (3) is globally asymptotically stable.Theorem 5.
Suppose thath<b, c<m, and h<(b/2)(1-c/m). The positive equilibrium E3=(u*,v*) of model (3) is globally asymptotically stable, if,(a1)
h(mv*+1)/b-c/m+βc/(mv*+1)>0;(a2)
h(h(mv*+1)/b+βc/(mv*+1)-1)(1-1/(βm+1))-h2/4b2-αb2/4(mv*+1)2>0;(a3)
4d11d22>(d12+d21)2,
whereα=(b-h)/hm>0, β=(b(m-c)-2hm)/2hm2.Proof.
We adopt the Lyapunov function:(26)V(t)=∫Ω[V1(u(x,t))+V2(v(x,t))]dx,
where V1(u)=(1/2)(u-u*)2, V2(v)=(1/2)(v-v*)2.
Then,(27)dVdt=∫Ω(∂u∂t(u-u*)+∂v∂t(v-v*))dx=∫Ω((u-u*)(u-u2-cuvmv+1)+(v-v*)(-hv+buvmv+1))dx+∫Ω((u-u*)(d11Δu+d12Δv)+(v-v*)(d21Δu+d22Δv))dx=I1+I2,
where
(28)I1=∫Ω(u-u*)×(u-u2-cuvmv+1-u*+u*2+cu*v*mv*+1)dx+∫Ω(v-v*)(-hv+buvmv+1+hv*-bu*v*mv*+1)dx,I2=∫Ω(u-u*)(d11Δu+d12Δv-d11Δu*-d12Δv*)dx+∫Ω(v-v*)(d21Δu+d22Δv-d21Δu*-d22Δv*)dx.
By some computational analysis, we obtain(29)I1=-∫Ω(u-u*)2(u+u*+cvmv*+1-1)dx-∫Ω(v-v*)2(h-h(mv+1))dx-∫Ω(u-u*)(v-v*)u*-bv(mv+1)(mv*+1)(mv+1)dx.
Considering the zero-flux boundary conditions, we have
(30)I2=-d11∫Ω|∇u|2dx-d12∫Ω∇u∇vdx+2d11∫Ω∇u∇u*dx+d12∫Ω∇u∇v*dx+d12∫Ω∇v∇v*dx-d11∫Ω|∇u*|2dx-d12∫Ω∇u*∇v*dx-d22∫Ω|∇v|2dx-d21∫Ω∇u∇vdx+d21∫Ω∇v∇u*dx+2d22∫Ω∇v∇v*dx+d21∫Ω∇u∇v*dx-d21∫Ω∇u*∇v*dx-d22∫Ω|∇v*|2dx=-d11∫Ω|∇u|2dx-(d12+d21)∫Ω∇u∇vdx+2d11∫Ω∇u∇u*dx-d11∫Ω|∇u*|2dx+(d12+d21)∫Ω∇u∇v*dx-(d12+d21)∫Ω∇u*∇v*dx+d12∫Ω∇v∇v*dx+2d22∫Ω∇v∇v*dx-d22∫Ω|∇v|2dx-d22∫Ω|∇v*|2dx+d21∫Ω∇v∇u*dx.
I
1 and I2 in the curly brackets can be expressed in the form -XAXT and -YBYT, respectively, where
(31)X=(u-u*,v-v*),Y=(∇u,∇v,∇u*,∇v*),A=(u+u*+cvmv*+1-1u*-bv(mv+1)2(mv*+1)(mv+1)u*-bv(mv+1)2(mv*+1)(mv+1)h-hmv+1)≜(φ1φ2φ3φ4),B=(d11d12+d212-d11d12+d212d12+d212d22d12+d212-d22-d11d12+d212d11d12+d212d12+d212-d22d12+d212d22).dV/dt is negative definite if the symmetric matrices A and B are positive. It can be easily shown that the symmetric matrix B is positive definite if
(32)4d11d22>(d12+d21)2.
The symmetric matrixA is positive definite if
(33)φ1>0,φ4>0,Φ(u,v)=φ1φ4-φ2φ3>0.
Since
(34)φ1=u+u*+cvmv*+1-1>h(mv*+1)b-cm+βcmv*+1,
due to (a1), φ1>0 is true. It is easy to verify that φ4=h-(h/(mv+1))>0.
Since(35)Φ(u,v)=h(u+u*+cvmv*+1-1)×(1-1mv+1)-14(u*-bv(mv+1)(mv*+1)(mv+1))2=h(u+h(mv*+1)b+cvmv*+1-1)×(1-1mv+1)-14(hb(mv+1)-bvmv*+1)2,
then
(36)∂Φ(u,v)∂u=h(1-1mv+1)>0.
Hence, Φ(u,v) is strictly increasing in ℝ+, with respect to u, and
(37)Φ(0,v)=h(h(mv*+1)b+cvmv*+1-1)(1-1mv+1)-14(hb(mv+1)-bvmv*+1)2≥h(h(mv*+1)b+βcmv*+1-1)(1-1βm+1)-h24b2-αb24(mv*+1)2.
Consequently, if (a2) holds, Φ(0,v)>0. As a result, Φ(u,v)>Φ(0,v)>0.
Hence,V is a Lyapunov function and the positive equilibrium E3 of model (3) is globally asymptotically stable. This completes the proof.Remark 6.
Whend12=d21=0, Theorem 5 is true, too. That is, the positive equilibrium E3 of the self-diffusion model (5) is globally asymptotically stable.


## 3. Turing’s Instability and Pattern Formation

### 3.1. Nonexistence of Turing’s Instability in the Self-Diffusion Model (5)

And in the presence of diffusion, we will introduce small perturbationsU1=u-u*, U2=v-v*, where |U1|,|U2|≪1. To study the effect of self-diffusion on model (5), we consider the linearized form of system about E*=(u*,v*) as follows:
(38)∂U1∂t=J11U1+J12U2+d11ΔU1,∂U2∂t=J21U1+J22U2+d22ΔU2,
where J11,J12,J21,andJ22 are defined as (24).Following Malchow et al. [11], we can know any solution of model (38) can be expanded into a Fourier series so that
(39)U1(x,t)=∑n,m=0∞‍unm(x,t)=∑n,m=0∞‍αnm(t)sinkx,U2(x,t)=∑n,m=0∞‍vnm(x,t)=∑n,m=0∞‍βnm(t)sinkx,
where x=(ξ,η), and 0<ξ<Lx, 0<η<Ly. k=(kn,km) and kn=nπ/Lx, km=mπ/Ly are the corresponding wavenumbers.Substitutingunm and vnm into (38), we obtain
(40)dαnmdt=(J11-d11k2)αnm+J12βnm,dβnmdt=J21αnm+(J22-d22k2)βnm,
where k2=kn2+km2.A general solution of (40) has the form C1exp(λ1t)+C2exp(λ2t), where the constants C1 and C2 are determined by the initial conditions and the exponents λ1,λ2 are the eigenvalues of the following matrix:
(41)D~=(J11-d11k2J12J21J22-d22k2).Correspondingly,λ1, λ2 are the solution of the following characteristic equation:
(42)λ2-ρ1λ+ρ2=0,
where
(43)ρ1=-k2(d11+d22)+tr(J),ρ2=d11d22k4-(d22J11+d11J22)k2+det(J).From (24), one can easily obtain
(44)ρ1<0,ρ2>0.
Then, we can conclude that the equilibrium E*=(u*,v*) is also stable for self-diffusion model (5). That is to say, there is nonexistence of Turing’s instability in model (5).


### 3.2. Turing’s Instability in the Cross-Diffusion Model (3)

The linearized form of the cross-diffusion model (3) about E*=(u*,v*) is as follows:
(45)∂U1∂t=J11U1+J12U2+d11ΔU1+d12ΔU2,∂U2∂t=J21U1+J22U2+d21ΔU1+d22ΔU2.The characteristic equation of the linearized model (3) is:
(46)λ2-σ1λ+σ2=0,
where
(47)σ1=-k2(d11+d22)+tr(J),σ2=det(D)k4-(d22J11+d11J22-d21J12-d12J21)k2+det(J).Diffusive instability occurs when at least one of the following conditions is violated [2]:
(48)σ1<0orσ2<0.It is evident that the conditionσ1>0 is not violated when the requirement J11+J22<0 is met because we assume d11>0 and d22>0. Hence, only violation of the condition σ2>0 will give rise to diffusion instability, that is, Turing’s instability. Then the condition for diffusive instability is given by
(49)d22J11+d11J22-d21J12-d12J21>0,
otherwise σ2>0 for all k>0 since det(D)>0 and det(J)>0.For Turing’s instability, we must haveσ2<0 for some k. And we notice that σ2 achieves its minimum:
(50)minμiσ2=4det(D)det(J)-(d22J11+d11J22-d21J12-d12J21)24det(D)
at the critical value kc2>0 when
(51)kc2=d22J11+d11J22-d21J12-d12J212det(D).
As a consequence, if d22J11+d11J22-d21J12-d12J21>0 and σ2<0 hold, then E3=(u*,v*) is an unstable equilibrium with respect to model (3). In this case, σ2=0 has two positive roots k12 and k22 which satisfy
(52)k1,22=d22J11+d11J22-d21J12-d12J21±Λ2det(D),
where Λ=(d22J11+d11J22-d21J12)2-4det(D)det(J). Therefore, if we can find some k2 such that k12<k2<k22, then σ2<0.Summarizing the above calculation, we obtain the following.Theorem 7.
Assume that the positive equilibriumE3=(u*,v*) exists. If the following conditions are true:(i)
d22J11+d11J22>d21J12+d12J21, that is,
(53)chd21>d22hv*2m2+hv*(bd11+2d22)m+b2d12v*+hd22;(ii)
d22J11+d11J22-d21J12-d12J21>2det(D)det(J), that is,(54)
c
h
d
21
-
(
d
22
h
v
*
2
m
2
+
h
v
*
(
b
d
11
+
2
d
22
)
m
+
b
2
d
12
v
*
+
h
d
22
)
>
2
b
h
v
*
(
d
11
d
22
-
d
12
d
21
)
(
h
m
3
v
*
2
+
2
h
m
2
v
*
+
h
m
+
b
c
)
,then the positive equilibriumE3 of model (3) is Turing unstable if  0<k12<k2<k22  for some k.


### 3.3. Pattern Formation

In this section, we perform extensive numerical simulations of the spatially extended model (3) in two-dimensional space, and the qualitative results are shown here. All our numerical simulations employ the zero-flux boundary conditions with a system size of 100×100. Other parameters are set as b=9, c=0.5, d=0.45, d11=0.01, d22=1, d12=-0.025, and d21=0.01.The numerical integration of model (3) is performed by using a finite difference approximation for the spatial derivatives and an explicit Euler method for the time integration [30, 31] with a time stepsize of 1/1000 and the space stepsize h=1/10. The initial condition is always a small amplitude random perturbation around the positive equilibrium E3=(u*,v*). After the initial period during which the perturbation spread, either the model goes into a time-dependent state or to an essentially steady-state solution (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 this section, we show the distribution of preyu, for instance. We have taken some snapshots with red (blue) corresponding to the high (low) value of prey u.Figure1 shows five typical Turing’s patterns of prey u in model (3) arising from random initial conditions for several values of the control parameter m.Five typical Turing’s patterns ofu in model (3) with fixed parameters b=9, c=0.5, d=0.45, d11=0.01, d22=1, d12=-0.025, and d21=0.01. (a) Spots pattern, m=0.05; (b) spot-stripe mixtures pattern, m=0.125; (c) stripes pattern, m=0.25; (d) hole-stripe mixtures, m=0.45; (e) holes pattern, m=0.5. Iterations: pattern (a): 5×105, pattern (c): 1×105, and others: 3×105.
(a)
(b)
(c)
(d)
(e)In Figure1(a)—H0-pattern, m=0.05, consists of red (maximum density of u) hexagons on a blue (minimum density of u) background, that is, isolated zones with high population densities. In this paper, we call this pattern as “spots.”In Figure1(b), when increasing m to 0.125, a few of stripes emerge, and the remainder of the spots pattern remains time independent. Pattern (b) is called H0-hexagon-stripe mixtures pattern.While increasingm to 0.25, model dynamics exhibits a transition from stripes-spots growth to stripes replication, that is, spots decay and the stripes pattern emerges (cf. Figure 1(c)).In Figure1(d), m=0.45, on increasing of m, a few of blue hexagons (i.e., holes, named by Von Hardenberg et al. [32], associated with low population densities) fill in the stripes, that is, the stripes-holes pattern emerges. Pattern (d) is called Hπ-hexagon-stripe mixtures pattern.When increasingm to 0.5, model dynamics exhibits a transition from stripe-holes growth to spots replication, that is, stripes decay and the holes pattern (Hπ-pattern) emerges (cf. Figure 1(e)).From Figure1, one can see that, on increasing the control parameter m, the sequences “spots → spot-stripe mixtures → stripes → hole-stripe mixtures → holes” are observed. Ecologically speaking, spots pattern shows that the prey population are driven by predators to a high level in those regions, while holes pattern shows that the prey population are driven by predators to a very low level in those regions. The final result is the formation of patches of high prey density surrounded by areas of low prey densities [3].


## 3.1. Nonexistence of Turing’s Instability in the Self-Diffusion Model (5)

And in the presence of diffusion, we will introduce small perturbationsU1=u-u*, U2=v-v*, where |U1|,|U2|≪1. To study the effect of self-diffusion on model (5), we consider the linearized form of system about E*=(u*,v*) as follows:
(38)∂U1∂t=J11U1+J12U2+d11ΔU1,∂U2∂t=J21U1+J22U2+d22ΔU2,
where J11,J12,J21,andJ22 are defined as (24).Following Malchow et al. [11], we can know any solution of model (38) can be expanded into a Fourier series so that
(39)U1(x,t)=∑n,m=0∞‍unm(x,t)=∑n,m=0∞‍αnm(t)sinkx,U2(x,t)=∑n,m=0∞‍vnm(x,t)=∑n,m=0∞‍βnm(t)sinkx,
where x=(ξ,η), and 0<ξ<Lx, 0<η<Ly. k=(kn,km) and kn=nπ/Lx, km=mπ/Ly are the corresponding wavenumbers.Substitutingunm and vnm into (38), we obtain
(40)dαnmdt=(J11-d11k2)αnm+J12βnm,dβnmdt=J21αnm+(J22-d22k2)βnm,
where k2=kn2+km2.A general solution of (40) has the form C1exp(λ1t)+C2exp(λ2t), where the constants C1 and C2 are determined by the initial conditions and the exponents λ1,λ2 are the eigenvalues of the following matrix:
(41)D~=(J11-d11k2J12J21J22-d22k2).Correspondingly,λ1, λ2 are the solution of the following characteristic equation:
(42)λ2-ρ1λ+ρ2=0,
where
(43)ρ1=-k2(d11+d22)+tr(J),ρ2=d11d22k4-(d22J11+d11J22)k2+det(J).From (24), one can easily obtain
(44)ρ1<0,ρ2>0.
Then, we can conclude that the equilibrium E*=(u*,v*) is also stable for self-diffusion model (5). That is to say, there is nonexistence of Turing’s instability in model (5).


## 3.2. Turing’s Instability in the Cross-Diffusion Model (3)

The linearized form of the cross-diffusion model (3) about E*=(u*,v*) is as follows:
(45)∂U1∂t=J11U1+J12U2+d11ΔU1+d12ΔU2,∂U2∂t=J21U1+J22U2+d21ΔU1+d22ΔU2.The characteristic equation of the linearized model (3) is:
(46)λ2-σ1λ+σ2=0,
where
(47)σ1=-k2(d11+d22)+tr(J),σ2=det(D)k4-(d22J11+d11J22-d21J12-d12J21)k2+det(J).Diffusive instability occurs when at least one of the following conditions is violated [2]:
(48)σ1<0orσ2<0.It is evident that the conditionσ1>0 is not violated when the requirement J11+J22<0 is met because we assume d11>0 and d22>0. Hence, only violation of the condition σ2>0 will give rise to diffusion instability, that is, Turing’s instability. Then the condition for diffusive instability is given by
(49)d22J11+d11J22-d21J12-d12J21>0,
otherwise σ2>0 for all k>0 since det(D)>0 and det(J)>0.For Turing’s instability, we must haveσ2<0 for some k. And we notice that σ2 achieves its minimum:
(50)minμiσ2=4det(D)det(J)-(d22J11+d11J22-d21J12-d12J21)24det(D)
at the critical value kc2>0 when
(51)kc2=d22J11+d11J22-d21J12-d12J212det(D).
As a consequence, if d22J11+d11J22-d21J12-d12J21>0 and σ2<0 hold, then E3=(u*,v*) is an unstable equilibrium with respect to model (3). In this case, σ2=0 has two positive roots k12 and k22 which satisfy
(52)k1,22=d22J11+d11J22-d21J12-d12J21±Λ2det(D),
where Λ=(d22J11+d11J22-d21J12)2-4det(D)det(J). Therefore, if we can find some k2 such that k12<k2<k22, then σ2<0.Summarizing the above calculation, we obtain the following.Theorem 7.
Assume that the positive equilibriumE3=(u*,v*) exists. If the following conditions are true:(i)
d22J11+d11J22>d21J12+d12J21, that is,
(53)chd21>d22hv*2m2+hv*(bd11+2d22)m+b2d12v*+hd22;(ii)
d22J11+d11J22-d21J12-d12J21>2det(D)det(J), that is,(54)
c
h
d
21
-
(
d
22
h
v
*
2
m
2
+
h
v
*
(
b
d
11
+
2
d
22
)
m
+
b
2
d
12
v
*
+
h
d
22
)
>
2
b
h
v
*
(
d
11
d
22
-
d
12
d
21
)
(
h
m
3
v
*
2
+
2
h
m
2
v
*
+
h
m
+
b
c
)
,then the positive equilibriumE3 of model (3) is Turing unstable if  0<k12<k2<k22  for some k.


## 3.3. Pattern Formation

In this section, we perform extensive numerical simulations of the spatially extended model (3) in two-dimensional space, and the qualitative results are shown here. All our numerical simulations employ the zero-flux boundary conditions with a system size of 100×100. Other parameters are set as b=9, c=0.5, d=0.45, d11=0.01, d22=1, d12=-0.025, and d21=0.01.The numerical integration of model (3) is performed by using a finite difference approximation for the spatial derivatives and an explicit Euler method for the time integration [30, 31] with a time stepsize of 1/1000 and the space stepsize h=1/10. The initial condition is always a small amplitude random perturbation around the positive equilibrium E3=(u*,v*). After the initial period during which the perturbation spread, either the model goes into a time-dependent state or to an essentially steady-state solution (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 this section, we show the distribution of preyu, for instance. We have taken some snapshots with red (blue) corresponding to the high (low) value of prey u.Figure1 shows five typical Turing’s patterns of prey u in model (3) arising from random initial conditions for several values of the control parameter m.Five typical Turing’s patterns ofu in model (3) with fixed parameters b=9, c=0.5, d=0.45, d11=0.01, d22=1, d12=-0.025, and d21=0.01. (a) Spots pattern, m=0.05; (b) spot-stripe mixtures pattern, m=0.125; (c) stripes pattern, m=0.25; (d) hole-stripe mixtures, m=0.45; (e) holes pattern, m=0.5. Iterations: pattern (a): 5×105, pattern (c): 1×105, and others: 3×105.
(a)
(b)
(c)
(d)
(e)In Figure1(a)—H0-pattern, m=0.05, consists of red (maximum density of u) hexagons on a blue (minimum density of u) background, that is, isolated zones with high population densities. In this paper, we call this pattern as “spots.”In Figure1(b), when increasing m to 0.125, a few of stripes emerge, and the remainder of the spots pattern remains time independent. Pattern (b) is called H0-hexagon-stripe mixtures pattern.While increasingm to 0.25, model dynamics exhibits a transition from stripes-spots growth to stripes replication, that is, spots decay and the stripes pattern emerges (cf. Figure 1(c)).In Figure1(d), m=0.45, on increasing of m, a few of blue hexagons (i.e., holes, named by Von Hardenberg et al. [32], associated with low population densities) fill in the stripes, that is, the stripes-holes pattern emerges. Pattern (d) is called Hπ-hexagon-stripe mixtures pattern.When increasingm to 0.5, model dynamics exhibits a transition from stripe-holes growth to spots replication, that is, stripes decay and the holes pattern (Hπ-pattern) emerges (cf. Figure 1(e)).From Figure1, one can see that, on increasing the control parameter m, the sequences “spots → spot-stripe mixtures → stripes → hole-stripe mixtures → holes” are observed. Ecologically speaking, spots pattern shows that the prey population are driven by predators to a high level in those regions, while holes pattern shows that the prey population are driven by predators to a very low level in those regions. The final result is the formation of patches of high prey density surrounded by areas of low prey densities [3].


## 4. Conclusions and Remarks

In this paper, we study the spatiotemporal dynamics of a Harrison predator-prey model with self- and cross-diffusions under the zero-flux boundary conditions. The value of this study lies in twofold. First, it gives the global stability of the positive equilibrium of the model by establishing a Lyapunov function. Second, it rigorously proves that the Turing instability can be induced by cross-diffusion, which shows that the model dynamics exhibits complex pattern replication controlled by the cross-diffusion.The most important observation in this paper is that the cross-diffusion terms are necessary for the emergence of Turing’s instability and pattern formation in the model. More precisely, with the help of the numerical simulations, the sequences “spots → spot-stripe mixtures → stripes → hole-stripe mixtures → holes” can be observed.On the other hand, population dynamics in the real world is inevitably affected by environmental noise which is an important component in an ecosystem. The deterministic models, such as model (3) or (5), assume that parameters in the systems are all deterministic irrespective environmental fluctuations. It is well known that the fact that due to environmental noise, the birth rate, carrying capacity, competition coefficient, and other parameters involved in the system exhibit random fluctuation to a greater or lesser extent [33]. We think that there may be exist other noise-controlled self-replicating patterns in models (3) and (5). This is desirable in future studies.It is believed that our results related to cross-diffusion in predator-prey interactions model would certainly be of some help to theoretical mathematicians and ecologists who are engaged in performing experimental work.



---

*Source: 306467-2013-02-02.xml*