# Pattern Formation in a Predator-Prey Model with Both Cross Diffusion and Time Delay

Authors: Boli Xie; Zhijun Wang; Yakui Xue
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A predator-prey model with both cross diffusion and time delay is considered. We give the conditions for emerging Turing instability in detail. Furthermore, we illustrate the spatial patterns via numerical simulations, which show that the model dynamics exhibits a delay and diffusion controlled formation growth not only of spots and stripe-like patterns, but also of the two coexist. The obtained results show that this system has rich dynamics; these patterns show that it is useful for the diffusive predation model with a delay effect to reveal the spatial dynamics in the real model.


---

## Body

## 1. Introduction

In recent years, the Lotka-Volterra model has been one of the important predator-prey models. However, this model has the unavoidable limitations to describe many realistic phenomena in biology. In order to well describe the real ecological interactions between the predator-prey species, the following predator-prey model has been proposed and studied [1]:
(1)
d
U
d
t
~
=
R
U
(
1
-
U
K
)
-
A
U
V
U
+
C
,
d
V
d
t
~
=
V
[
θ
(
1
-
B
V
U
)
]
,
where U and V stand for prey and predator density, respectively. The first equation states that in absence of predation the prey grow logistically with carrying capacity K and intrinsic growth rate R. The saturating predator functional response A
U
V
/
(
U
+
C
) used in (1) is of Michaelis-Menten type in enzyme-substrate kinetics. The parameter A is the maximum specific rate of product formation and C (the half-saturation constant) is the substrate density at which the rate of product formation is half-maximal. The second equation shows that predators grow logistically with intrinsic growth rate θ. The parameter B is the number of prey required to support one predator at equilibrium when V equals U
/
B [2–4].Following Hsu and Huang [4], with the next scaling
(2)
u
=
U
K
,
v
=
A
V
R
K
,
t
~
=
R
t
,
d
=
C
K
,
s
=
θ
R
,
γ
=
B
R
A
,
we arrive at the following equations containing dimensionless quantities:
(3)
d
u
d
t
=
u
(
1
-
u
)
-
u
v
u
+
d
,
d
v
d
t
=
v
[
s
(
1
-
γ
v
u
)
]
.Spatial patterns are ubiquitous in nature; these patterns modify the temporal dynamics and stability properties of the population densities in a range of spatial scales. Their effects must be incorporated in temporal ecological models that do not represent space explicitly. When combined with spatial factor and diffusion terms, the original spatially extended model is written as the following system:(4)
d
u
d
t
=
u
(
1
-
u
)
-
u
v
u
+
d
+
D
u
∇
2
u
=
f
(
u
,
v
)
+
D
u
∇
2
u
d
u
d
t
=
u
(
1
-
u
)
-
h
(
u
,
v
)
+
D
u
∇
2
u
,
d
v
d
t
=
v
[
s
(
1
-
γ
v
u
)
]
+
D
v
∇
2
v
=
g
(
u
,
v
)
+
D
v
∇
2
v
d
u
d
t
=
v
[
s
(
1
-
q
(
u
,
v
)
)
]
+
D
v
∇
2
v
,
where ∇
2
=
∂
2
/
∂
x
2 or ∇
2
=
∂
2
/
∂
x
2
+
∂
2
/
∂
y
2 is the usual Laplacian operator in the one- or two-dimensional space. The diffusion coefficients are denoted by by D
u and D
v, respectively.On the other hand, time delay plays an important role in many biological dynamical systems, being particularly relevant in ecology, where time delays have been recognized to contribute critically to the outcome for prey densities under predation being stable or unstable [5]. Time delay due to gestation is included in some predator-prey models, because generally a duration of τ time units elapses between the time when an individual prey is killed and the moment when a corresponding increase in the predator population is realized [6]. The effect of this kind of delay on the dynamical behavior of populations has been studied by a number of papers [5–8].However, to the best of our knowledge, there is little work on the dynamical behavior of both time delay and diffusion in the predator-prey model. As a result, in the present paper, we aim to study the predator-prey model with both cross diffusion and time delay. More specifically, the present paper is mainly to investigate the spatial patterns. And the model is given by(5)
∂
u
∂
t
=
u
(
1
-
u
)
-
u
v
(
t
-
τ
)
u
+
d
+
D
u
∇
2
u
,
∂
v
∂
t
=
v
[
s
(
1
-
γ
v
(
t
-
τ
)
u
)
]
+
D
v
∇
2
v
,
where τ
>
0 is a constant delay due to gestation.Model (5) needs to be analyzed with the initial populations
(6)
u
(
0
)
>
0
,
v
(
0
)
>
0
.We also assume that no external input is imposed from outside. Hence, the boundary conditions are taken as(7)
∂
u
∂
n
|
(
x
,
y
)
=
∂
u
∂
n
|
(
x
,
y
)
=
0
,
where (
x
,
y
)
∈
∂
Ω and Ω is the spatial domain.This paper is organized as follows. In Section2, by using the method of linear stability analysis, we deduce the conditions under which instability might occur. In Section 3, we perform a series of numerical simulations to show the evolution process of prey u. Finally, in Section 4, we give some concluding remarks.


## 2. Analysis for the Model

In this section, we will discuss the stability of model (5). It is easy to see that model (5) has the same equilibria as model (3). However, we find that model (3) exhibits two equilibria:(i)
E
0
=
(
1,0
), which is corresponding to extinction of the predator;
(ii)
interior equilibrium pointE
*
(
u
*
,
v
*
), which is corresponding to coexistence of prey and predator and
(8)
u
*
=
γ
-
1
-
γ
d
+
(
γ
-
1
-
γ
d
)
2
+
4
d
γ
2
2
γ
,
v
*
=
u
*
γ
.There has been some works on the stability analysis of model (3) [4]. However, the main purpose of the present paper is to investigate the effect of both cross diffusion and time delay on the spatial pattern. Following [9, 10], assume that τ is small enough; we replace v
(
x
,
y
,
t
-
τ
) as follows:
(9)
v
(
x
,
y
,
t
-
τ
)
=
v
(
x
,
y
,
t
)
-
τ
∂
v
(
x
,
y
,
t
)
∂
t
.Substituting (9) into model (5), we obtain
(10)
∂
u
∂
t
=
u
(
1
-
u
)
-
h
(
u
(
x
,
y
,
t
)
,
v
(
x
,
y
,
t
)
-
τ
∂
v
(
x
,
y
,
t
)
∂
t
)
∂
u
∂
t
=
+
D
u
∇
2
u
,
∂
v
∂
t
=
v
[
s
(
1
-
q
(
u
(
x
,
y
,
t
)
,
v
(
x
,
y
,
t
)
-
τ
∂
v
(
x
,
y
,
t
)
∂
t
)
)
]
∂
u
∂
t
=
+
D
v
∇
2
v
.Expanding (10) in Taylor series and neglecting the higher order nonlinearities, then (10) becomes
(11)
∂
u
∂
t
=
u
(
1
-
u
)
-
h
(
u
,
v
)
+
τ
h
v
(
u
,
v
)
∂
v
(
x
,
y
,
t
)
∂
t
∂
u
∂
t
=
+
D
u
∇
2
u
,
∂
v
∂
t
=
v
[
s
(
1
-
q
(
u
,
v
)
+
τ
q
v
(
u
,
v
)
∂
v
(
x
,
y
,
t
)
∂
t
)
]
∂
u
∂
t
=
+
D
v
∇
2
v
,
where h
v
(
u
,
v
)
=
∂
h
(
u
,
v
)
/
∂
v, q
v
(
u
,
v
)
=
∂
q
(
u
,
v
)
/
∂
v.From (11), we finally obtain
(12)
∂
u
∂
t
=
f
(
u
,
v
)
+
τ
h
v
(
u
,
v
)
1
1
-
v
τ
s
q
v
(
u
,
v
)
(
g
(
u
,
v
)
+
D
v
∇
2
v
)
∂
u
∂
t
=
+
D
u
∇
2
u
,
∂
v
∂
t
=
1
1
-
v
τ
s
q
v
(
u
,
v
)
(
g
(
u
,
v
)
+
D
v
∇
2
v
)
.To see how the system responds when the steady stateE
* is perturbed, we consider small spatiotemporal perturbations δ
u
(
x
,
y
,
t
) and δ
v
(
x
,
y
,
t
) around the steady state E
*
(
u
*
,
v
*
) as follows:
(13)
u
(
x
,
y
,
t
)
=
u
*
+
δ
u
(
x
,
y
,
t
)
,
v
(
x
,
y
,
t
)
=
v
*
+
δ
v
(
x
,
y
,
t
)
.Linearizing model (13) around E
*
(
u
*
,
v
*
), we obtain
(14)
∂
δ
u
∂
t
=
(
δ
u
)
f
u
+
(
δ
v
)
f
v
+
τ
h
v
1
1
-
v
*
τ
s
q
v
(
u
,
v
)
∂
δ
u
∂
t
=
×
[
(
δ
u
)
g
u
+
(
δ
v
)
g
v
+
D
v
∇
2
(
δ
v
)
]
+
D
u
∇
2
(
δ
u
)
,
∂
δ
v
∂
t
=
1
1
-
v
*
τ
s
q
v
(
u
,
v
)
∂
δ
u
∂
t
=
×
[
(
δ
u
)
g
u
+
(
δ
v
)
g
v
+
D
v
∇
2
(
δ
v
)
]
.From (14), we obtain
(15)
∂
(
δ
u
)
∂
t
=
(
f
u
+
τ
h
v
g
u
1
-
v
*
τ
s
q
v
)
(
δ
u
)
+
(
f
v
+
τ
h
v
g
v
1
-
v
*
τ
s
q
v
)
(
δ
v
)
∂
(
δ
u
)
∂
t
=
+
τ
h
v
1
-
v
*
τ
s
q
v
D
v
∇
2
(
δ
v
)
+
D
u
∇
2
(
δ
u
)
∂
(
δ
u
)
∂
t
≜
(
f
u
+
χ
τ
h
v
g
u
)
(
δ
u
)
+
(
f
v
+
χ
τ
h
v
g
v
)
(
δ
v
)
∂
(
δ
u
)
∂
t
=
+
χ
τ
h
v
D
v
∇
2
(
δ
v
)
+
D
u
∇
2
(
δ
u
)
,
∂
(
δ
v
)
∂
t
=
1
1
-
v
*
τ
s
q
v
(
u
,
v
)
∂
(
δ
v
)
∂
t
=
×
[
g
u
(
δ
u
)
+
g
v
(
δ
v
)
+
D
v
∇
2
(
δ
v
)
]
∂
(
δ
v
)
∂
t
≜
χ
g
u
(
δ
u
)
+
χ
g
v
(
δ
v
)
+
χ
D
v
∇
2
(
δ
v
)
,
where h
v
=
(
∂
h
/
∂
v
)
|
(
u
*
,
v
*
), q
v
=
(
∂
q
/
∂
v
)
|
(
u
*
,
v
*
), and (
1
/
(
1
-
v
*
τ
s
q
v
)
)
|
(
u
*
,
v
*
)
=
1
/
(
1
-
s
τ
)
≜
χ. Equation (12) can be used to analyze the dynamic behavior of model (5) when τ is small, so we only consider the case of τ
<
1
/
s  (i.e., χ
>
0), in this paper.Assume that the solution of (15) takes the form
(16)
δ
u
(
x
,
y
,
t
)
=
δ
u
*
e
λ
t
cos
⁡
⁡
(
k
x
x
)
cos
⁡
⁡
(
k
y
y
)
,
δ
v
(
x
,
y
,
t
)
=
δ
v
*
e
λ
t
cos
⁡
⁡
(
k
x
x
)
cos
⁡
⁡
(
k
y
y
)
,
where λ is the growth rate of the perturbation in time t, δ
u
* and δ
v
* represent the amplitudes, and k
x and k
y are the wavenumbers of the solutions. And upon inserting them in (15), we obtain the characteristic equation at E
* of model (5):
(17)
det
⁡
⁡
(
λ
I
-
J
k
)
=
λ
2
-
tr
⁡
(
J
k
)
λ
+
det
⁡
⁡
(
J
k
)
=
0
,
where
(18)
J
k
=
(
f
u
+
χ
τ
h
v
g
u
-
D
u
k
2
f
v
+
χ
τ
h
v
g
v
-
χ
τ
h
v
D
v
k
2
χ
g
u
χ
g
v
-
χ
D
v
k
2
)
,
(19)
tr
⁡
(
J
k
)
=
f
u
+
(
g
v
+
τ
h
v
g
u
)
χ
-
(
D
u
+
χ
D
v
)
k
2
,
(20)
det
⁡
⁡
(
J
k
)
=
χ
D
u
D
v
k
4
-
χ
(
f
u
D
v
+
g
v
D
u
)
k
2
det
⁡
⁡
(
J
k
)
=
+
χ
(
f
u
g
v
-
f
v
g
u
)
.The roots of (17) can be obtained by the following form:
(21)
λ
1,2
(
k
)
=
tr
⁡
(
J
k
)
±
tr
⁡
(
J
k
)
2
-
4
det
⁡
⁡
(
J
k
)
2
.Turing instability means that it is stable for nonspatial model (3) but is unstable with respect to the solutions of the spatial model (5). The stability of nonspatial model (3) is guaranteed if
(22)
tr
⁡
(
A
)
=
f
u
+
g
v
<
0
,
det
⁡
⁡
(
A
)
=
f
u
g
v
-
f
v
g
u
>
0
,
where
(23)
A
=
(
∂
f
∂
u
∂
f
∂
v
∂
g
∂
u
∂
g
∂
v
)
(
u
*
,
v
*
)
A
≜
(
f
u
f
v
g
u
g
v
)
A
=
(
1
-
2
u
*
-
v
*
d
(
u
*
+
d
)
2
-
u
*
u
*
+
d
s
γ
-
s
)
and the Turing instability sets in when at least one of tr
⁡
(
J
k
)
<
0 and det
⁡
(
J
k
)
>
0 is violated. Thus, we consider the emergence of the instability in the following two cases: (i)
det
⁡
(
J
k
)
>
0 is violated;
(ii)
tr
⁡
(
J
k
)
<
0 is violated.First, we considerdet
⁡
(
J
k
)
>
0 is violated.From(24)
det
⁡
⁡
(
J
k
)
=
χ
(
D
u
D
v
k
4
-
(
f
u
D
v
+
g
v
D
u
)
k
2
+
f
u
g
v
-
f
v
g
u
)
det
⁡
⁡
(
J
k
)
<
0
,
simple algebraic computation leads to
(25)
D
v
f
u
+
D
u
g
v
>
0
,
(
D
v
f
u
+
D
u
g
v
)
2
>
4
D
u
D
v
(
f
u
g
v
-
f
v
g
u
)
.Sincetr
⁡
(
J
k
)
<
0 has to be negative for some values of k, from (19), we notice that the following conditions must be satisfied:
(26)
f
u
+
(
g
v
+
τ
h
v
g
u
)
χ
<
0
,
which is equivalent to
(27)
0
≤
τ
<
-
f
u
+
g
v
h
v
g
u
-
f
u
s
≜
τ
c
.The above uses the fact that(28)
h
v
g
u
-
f
u
s
=
(
-
f
v
)
g
u
-
f
u
(
-
g
v
)
=
f
u
g
v
-
f
v
g
u
>
0
.Hence, in this case,τ must satisfy
(29)
τ
<
min
⁡
⁡
(
τ
c
,
1
s
)
.If the following conditions (22), (25), and (29) hold, the positive equilibrium (
u
*
,
v
*
) of model (5) is unstable.To well see the effect of cross diffusion and time delay, we plot the dispersion relation keeping the parameter values fixed in Figure1. It can be seen from Figure 1 that Turing modes Re
(
λ
)
>
0 can be available.Figure 1
An illustration of the dispersion relation from (21). We set the parameter values such that (a) γ
=
1, d
=
0.08, s
=
1.55, τ
=
0.1, D
u
=
0.1, and D
v
=
3.68; (b) γ
=
1, d
=
0.08, s
=
1.31, τ
=
0.2, D
u
=
0.1, and D
v
=
3.68; (c) γ
=
1, d
=
0.05, s
=
0.5, τ
=
0.3, D
u
=
0.25, and D
v
=
3.Next, we consider the second case, wheretr
⁡
(
J
k
)
<
0 is violated. On the basis of the same discussions as the above, it is well known that tr
⁡
(
J
k
)
<
0 is violated when the following inequality is satisfied:
(30)
f
u
+
(
g
v
+
τ
h
v
g
u
)
χ
>
0
;
that is,
(31)
τ
>
-
f
u
+
g
v
h
v
g
u
-
f
u
s
≜
τ
c
,
(32)
τ
c
<
τ
<
1
s
.It is to be noted that if(33)
D
v
f
u
+
D
u
g
v
<
0
,det
⁡
(
J
k
)
>
0 is valid for all k.If the following conditions (22), (32), and (33) hold, the positive equilibrium (
u
*
,
v
*
) of model (5) is unstable.Furthermore, the second case is similar to that in the first case. In the following section, we only discuss the first case, namely,det
⁡
(
J
k
)
<
0.


## 3. Pattern Structures

In practice, the continuous problem defined by the reaction-diffusion system in two-dimensional space is solved in a discrete domain withM
×
N lattice sites (i.e., abscissa axis and ordinate axis, resp., M and N). The spacing between the lattice points is defined by the lattice constant Δ
h. For Δ
h
→
0 the differences approach the derivatives. The time evolution is also discrete; that is, the time goes in step of Δ
t. In the present paper, we set Δ
h
=
1, Δ
t
=
0.01, and M
=
N
=
200. Note that when Δ
h
,
Δ
t are further decreased, the dynamics does not change any more.We run the simulations until they reach a stationary state or until they show a behavior that does not seem to change its characteristics anymore. In the simulations different types of dynamics are observed and we have found that the distributions ofu and v are always of the same type. As a result, we can restrict our analysis of pattern formation to one distribution (in this paper, we show the distribution of u, for instance).In Figure2, we set γ
=
1, d
=
0.08, s
=
1.55, τ
=
0.1, D
u
=
0.1, and D
v
=
3.68 and the steady state solution is (
u
*
,
v
*
)
=
(
0.2456571372, 0.2456571372
). After irregular transient pattern, we can see that the regular spotted patterns with the same radius prevail over the whole domain finally, and the dynamics of the system do not undergo any further changes.Snapshots of the time evolution of the prey at different instants withγ
=
1, d
=
0.08, s
=
1.55, τ
=
0.1, D
u
=
0.1, and D
v
=
3.68, which are in the Turing space. (a) 0 iteration; (b) 10000 iterations; (c) 50000 iterations; and (d) 300000 iterations.
(a)
(b)
(c)
(d)In Figure3, we set γ
=
1, d
=
0.08, s
=
1.31, τ
=
0.2, D
u
=
0.1, and D
v
=
3.68 and the steady state solution is (
u
*
,
v
*
)
=
(
0.2456571372, 0.2456571372
). After the irregular pattern forms, stripe like and spotted patterns emerge mixed in the distribution of the infected population density, and the dynamics of the system do not undergo any further changes.Snapshots of the time evolution of the prey at different instants withγ
=
1, d
=
0.08, s
=
1.31, τ
=
0.2, D
u
=
0.1, and D
v
=
3.68, which are in the Turing space. (a) 0 iteration; (b) 10000 iterations; (c) 50000 iterations; and (d) 300000 iterations.
(a)
(b)
(c)
(d)In Figure4, we set γ
=
1, d
=
0.05, s
=
0.5, τ
=
0.3, D
u
=
0.25, and D
v
=
3 and the steady state solution is (
u
*
,
v
*
)
=
(
0.2, 0.2
). We can see that the regular stripe patterns prevail over the whole domain at last, and the dynamics of the system do not undergo any further changes.Snapshots of the time evolution of the prey at different instants withγ
=
1, d
=
0.05, s
=
0.5, τ
=
0.3, D
u
=
0.25, and D
v
=
3, which are in the Turing space. (a) 0 iteration; (b) 10000 iterations; (c) 50000 iterations; and (d) 300000 iterations.
(a)
(b)
(c)
(d)


## 4. Discussions

In this paper, we analyze the spatiotemporal dynamics of a spatial predator-prey model with both time delay and cross diffusion. A series of numerical simulations reveal that the typical dynamics of population density variation is the formation of isolated groups, that is, spotted or stripe-like or coexistence of both. We have already presented three kinds of figures showing different patterns when different delays are used. That is to say, the interaction of delay and diffusion can create stationary patterns.Although more work is needed, in principle, it seems that delay and diffusion are able to generate many different kinds of spatiotemporal patterns. For such reasons, we can predict that delay and diffusion can be considered as an important mechanism for the appearance of complex spatiotemporal dynamics in ecology models.



---

*Source: 392435-2014-05-12.xml*