- Email: [email protected]
Cross-diffusion induced spatiotemporal patterns in a predator–prey model with herd behavior
Nonlinear Analysis: Real World Applications 24 (2015) 36–49
Contents lists available at ScienceDirect
Nonlinear Analysis: Real World Applications journal homepage: www.elsevier.com/locate/nonrwa
Cross-diffusion induced spatiotemporal patterns in a predator–prey model with herd behavior Xiaosong Tang a,b , Yongli Song a,∗ a
Department of Mathematics, Tongji University, Shanghai 200092, China
b
College of Mathematical and Physics, Jinggangshan University, Ji’an 343009, China
article
info
Article history: Received 15 August 2014 Received in revised form 25 December 2014 Accepted 31 December 2014
Keywords: Predator–prey model Herd behavior Cross-diffusion Turing instability Amplitude equations Spatiotemporal pattern
abstract In this paper, we investigate the phenomena of pattern formation for a predator–prey model with herd behavior and cross-diffusion. We give the conditions for cross-diffusion induced instability in detail, and have derived amplitude equations for the excited modes. Furthermore, we illustrate the spatial patterns via numerical simulations, which show that the model dynamics exhibits complex pattern replication: spotted pattern, stripe pattern, and coexistence of the two by increasing the control parameter β and cross-diffusion coefficient d21 respectively. Meanwhile, we also obtain the other patterns: ‘‘labyrinth’’ patterns and ‘‘black-eye’’ patterns. These results indicate that the cross-diffusion plays an important role in the pattern selection. The spatiotemporal patterns induced by cross-diffusion may enrich the pattern dynamics for diffusive predator–prey model and be helpful for the biologists to understand the phenomenon in which a flux of one species is induced by a gradient of another species. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction In ecological systems, the interactions of different species within a fluctuating natural environment are the most distinctive features. Studies in mathematical models are informative in understanding the interactions of predator and prey in these systems. As a result, predator–prey models have been important in ecological science since the pioneering works of Lotka and Volterra [1,2]. However, these models have the unavoidable limitations to describe many realistic phenomena in biology. Recently, a predator–prey model is considered in which the prey exhibits herd behavior, so that the predator interacts with the prey along the outer corridor of the herd of prey, which is more appropriate to model the response functions of prey that exhibit herd behavior in terms of the square root of the prey population, due originally to [3,4]:
√ dX = rX 1 − X − α X Y√ , dt K √ 1 + th α X (1.1) cα X Y dY = −sY + √ , dt 1 + th α X where X (t ) and Y (t ) stand for the prey and predator densities at time t, respectively. The parameter r is the growth rate of the prey, K is its carrying capacity, and s is the death rate of the predator in the absence of prey. The parameter α is the search efficiency of Y for X , c is the conversion or consumption rate of prey to predator, and th is average handling time. ∗
Corresponding author. E-mail address: [email protected] (Y. Song).
http://dx.doi.org/10.1016/j.nonrwa.2014.12.006 1468-1218/© 2015 Elsevier Ltd. All rights reserved.
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
37
For simplicity, we nondimensionalize model (1.1) with the following scaling: u=
X K
αY v= √ ,
,
r
K
tnew = rtold ,
√
rL1 ,
r1 =
√ r2 =
rL2 ,
and the other parameters are made dimensionless as follows:
√ snew =
K
α
sold ,
cnew =
√ α K r
cold ,
a = th α
√
K.
Then model (1.1) converts into
√ uv du √ , dt = u(1 − u) − 1+a u √ dv c uv = −sv + √ . dt 1+a u
(1.2)
On the other hand, the issue of spatial and temporal pattern formation in biological communities is probably one of the most exciting problems in modern biology and ecology. Spatial pattern formation arose from the observation in chemistry by Turing [5] that diffusion can also destabilize equilibrium solution, a scenario well known as Turing instability. Although this is counter-intuitive, the interest in the diffusion-driven instability has long expanded from chemical system to biological system [6–12]. Thus, the following reaction–diffusion system is obtained as a model for predator–prey with herd behavior:
√ ∂u uv = u ( 1 − u ) − √ + d11 ∇ 2 u, ∂t 1+a u √ ∂v c uv = −sv + √ + d22 ∇ 2 v, ∂t 1+a u
(1.3)
where the nonnegative constants d11 and d22 , are the diffusion coefficients of prey and predator respectively. Nevertheless, from [4], we know that the phenomenon of spatial pattern formation in (1.3) cannot occur under all possible diffusion rates. So, the authors in [4,11] change linear mortality sv into quadratic mortality sv 2 in (1.3) and have investigated Turing patterns, stability, Turing instability and Hopf bifurcations. However, the prey may be able to recognize the predator, and thus the prey would keep away from predators in order to not to be caught. In nature, the escape velocity of the prey might be proportional to the diffusive velocity of the predators. On the other hand, the tendency of predators would get closer to the prey, and hence the chase velocity of predators may be considered to be proportional to the diffusive velocity of the prey [13]. This phenomenon is called cross diffusion. Therefore the problem of cross-diffusion arises, which was proposed first by Kerner [14] and first applied to competitive population systems by Shigesada et al. [15]. Recently, there has been increasing interest in the stability behavior of predator–prey systems, in which both self-diffusion and cross-diffusion effects are considered [16–23], and references therein. The cross-diffusion describes the population flux of one species due to the presence of other species. Unfortunately, most studies on the spatial and temporal predator–prey systems pay little attention to the study on the effect of cross-diffusion and the selection of Turing patterns. They just focus on the bifurcation phenomena by varying the controlled parameter or only give some numerical results of the patterns [24–32]. Based on the previous work of many scientists, we can now focus on the following predator–prey model with herd behavior and cross-diffusion:
√ uv ∂u √ + d11 ∇ 2 u + d12 ∇ 2 v, ∂ t = u(1 − u) − 1+a u √ c uv ∂v = −sv + √ + d21 ∇ 2 u + d22 ∇ 2 v. ∂t 1+a u
(1.4)
In [3], the author has used the simplifying assumption that a = 0, which implies that the average handling time is zero. For the convenience of discussion and in line with the work of [3], we also assume that a = 0 and let γ = c , β = s/γ . Then Eq. (1.4) is simplified as
√ ∂u = u(1 − u) − uv + d11 ∇ 2 u + d12 ∇ 2 v, ∂t ∂v = γ v(−β + √u) + d21 ∇ 2 u + d22 ∇ 2 v, ∂t
(1.5)
where d12 and d21 , called cross-diffusion coefficients, describe the respective population fluxes of preys and predators resulting from the presence of the other species, respectively. Model (1.5) is to be analysed under the following non-zero initial conditions: u(x, y, 0) > 0,
v(x, y, 0) > 0,
(x, y) ∈ Ω = [0, L] × [0, L],
38
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
and zero-flux boundary conditions ∂u ∂v = = 0, ∂n ∂n where L denotes the size of the spatial variables x and y, n is the outward unit normal vector of the boundary ∂ Ω . The goal of this paper is to explore whether cross-diffusion can drive the emergence of Turing pattern and how different cross-diffusion coefficient effects the spatially inhomogeneous distribution of population density. Our analyses reveal that the cross-diffusion may induce spatial patterns and different cross-diffusion coefficient may lead to pattern transition between spots and stripes. The paper is organized as follows. In Section 2, we show the spatial dynamics of model (1.5) under the fixed crossdiffusion and varied cross-diffusion. In Section 3, we present and discuss the results of Turing pattern formation through extensive numerical simulations. Finally, conclusions and discussions are presented in Section 4. 2. Spatial dynamics of model (1.5) 2.1. Linear stability analysis In the absence of diffusion, simple calculations show that model (1.5) has three equilibrium points which consist of two boundary equilibriums (0, 0), (1, 0) and a positive equilibrium (u∗ , v ∗ ) = (β 2 , β(1 − β 2 )), β ∈ (0, 1). From the biological perspective, we are interested in studying the stability behavior of the positive equilibrium point. The Jacobian matrix corresponding to (u∗ , v ∗ ) is as follows:
J :=
a12 a22
a11 a21
1
2 2 (1 − 3β ) = 1 γ (1 − β 2 )
−β .
(2.1)
0
2 Following the standard linear analysis of the cross-diffusion system [33,34], we address the temporal stability of the uniform state to nonuniform perturbations
∗
u
v
u
=
v∗
+ε
uk
vk
eλt +ik·r + O(ε 2 ) + c .c .,
(2.2)
where k · k = k2 , k and λ are the wave number and frequency respectively; i is the imaginary unit and i2 = −1; r = (u, v) is the spatial vector in two-dimensional space and c .c. stands for the complex conjugate terms. The linear instability(ε ≪ 1) of the uniform state is deduced from the dispersion relations. We obtain that the eigenvalue is the root of the following equation:
λ2 − Tk λ + Dk = 0, (2.3) where Tk = a11 + a22 − (d1 + d2 )k2 , Dk = a11 a22 − a12 a21 − (a11 d22 + a22 d11 − a12 d21 − a21 d12 )k2 + (d11 d22 − d12 d21 )k4 . The roots of Eq. (2.3) can be obtained by the following form:
λk =
Tk ±
Tk2 − 4Dk
.
2 We choose β as the bifurcation parameter. The Hopf bifurcation occurs when Im(λk ) ̸= 0,
Re(λk ) = 0,
at k = 0,
where λk is the root of Eq. (2.3). Then we can get the critical value of the Hopf bifurcation parameter β :
√
βH =
3
.
(2.4)
3 The Turing bifurcation occurs when Im(λk ) = 0,
Re(λk ) = 0,
at k = kT ̸= 0,
and the wave number kT satisfies k2T =
βT :
a11 a22 −a12 a21 . d11 d22 −d12 d21
We can obtain the critical value of the Turing bifurcation parameter
(1 − 3β 2 )d22 + 2β d21 − (1 − β 2 )γ d12 − 2 2(d11 d22 − d12 d21 )βγ (1 − β 2 ) = 0.
(2.5)
Remark 2.1. When d12 = d21 = 0 in (1.5), a general linear analysis shows that the necessary conditions for yielding Turing √ patterns are given by T0 < 0, D0 > 0 and d22 a11 + d11 a22 − 2 d11 d22 D0 > 0. However, these conditions cannot hold simultaneously. √ In fact, T0 < 0 implies that a11 < 0 since a22 = 0. In addition, noticing that d11 > 0 and d22 > 0, we have d22 a11 + d11 a22 − 2 d11 d22 D0 < 0. So, the phenomenon of Turing pattern formation in (1.5) cannot occur when the cross-diffusion coefficients d12 and d21 are equal to zero. This implies that the cross-diffusion coefficients are crucial for the Turing instability.
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
39
Fig. 1. Bifurcation diagram of model (1.5). We set the parameter values as d11 = d12 = 1, d21 = 11, d22 = 15. The figure shows the Turing space D11 which is located below the positive equilibrium existence(the black curve), i.e. the region bounded by the Turing bifurcation curve (the blue (upper) curve) and the Hopf bifurcation curve (the red (lower) curve), D12 denotes Hopf space and D13 denotes Turing–Hopf space. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 2. (Color online) Bifurcation diagram of model (1.5). We set the parameter values as d11 = d12 = 1, d22 = 15. Turing and Hopf instability boundaries varying d21 . The instabilities stay below the lines.
Now, let us discuss the bifurcations represented by these formulas in the parameter region spanned by the parameters β and γ that can be seen from Fig. 1. The whole class of spatial model is included in this parameter region, where the domain D11 is called as Turing space, in which the Turing instability occurs, the domain D12 is called as Hopf space, in which the Hopf instability occurs and the region D13 is called Turing–Hopf space, where a competition between the two instabilities occurs. In D13 which one would develop, depends on the locations of the respective instability boundaries: as β decreases, if βT > βH , Turing instability occurs prior to the oscillatory instability. In Fig. 2, we report the Hopf and Turing instability boundaries for different values of d21 , and observe that, while the Hopf instability region does not depend on d21 , the Turing pattern region decreases by decreasing d21 due to the fact that βT decreases with d21 . Then we plot, in Fig. 3, the dispersion relations corresponding to several values of bifurcation parameter β , while keeping other parameters fixed. In Fig. 3, curve (2) corresponds to the critical Turing value β2 = 0.7691, curve (4) corresponds to the critical Hopf value β4 = 0.5773. On the other hand, we focus on how cross diffusion effects the dynamics for fixed deterministic parameters. To this end, we fix the deterministic model parameter values d11 = d12 = 1, d22 = 5, γ = 0.4, and vary d21 as a function of β , which is the coefficient of the cross diffusion of the prey. We show the bifurcation diagram of linear stability analysis in Fig. 4, where the domain D21 is called as Turing space, the domain D22 is called as Hopf space and the region D23 is called Turing–Hopf space. In D23 which one would develop, depends on the locations of the respective instability boundaries: as β increases, if βT > βH , Turing instability occurs prior to the oscillatory instability. In Fig. 5, we report the Hopf and Turing instability boundaries for different values of γ , and observe that, while the Hopf instability region does not depend on γ , the Turing pattern region increases by decreasing γ due to the fact that βT increases with γ . To see the Turing space properly, we plot the dispersion relation corresponding to several values of parameter d21 by keeping the β = 0.678. It can be seen from Fig. 6 that when d21 is decreased, the available Turing modes [Re(λ) > 0] are decreased, and all available modes are weakened.
40
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
Fig. 3. The relation between Re(λ) (the real part of the eigenvalue λ) and k with d11 = d12 = 1, d21 = 11, d22 = 15, γ = 2 and different β . Curve (1): β1 = 0.7820; Curve (2): β2 = 0.7691; Curve (3): β3 = 0.6250; Curve (4): β4 = 0.5773; Curve (5): β5 = 0.5560.
Fig. 4. Bifurcation diagram of model (1.5). We set the parameter values as d11 = d12 = 1, d22 = 5, γ = 0.4. The figure shows the Turing space D21 which is located below the positive equilibrium existence (the black curve), i.e. the region bounded by the Turing bifurcation curve (the blue (upper) curve) and the Hopf bifurcation curve (the red (lower) curve), D22 denotes Hopf space and D23 denotes Turing–Hopf space. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
From the above analysis, we can obtain the conditions of Turing instability for model (1.5), that is determining the conditions of the emergency of Turing patterns. But we cannot determine the selection and competition of Turing patterns. In the following section, we shall consider the selection of Turing patterns close to the onset β = βT = β2 .
2.2. Amplitude equations In this subsection, we can deduce the well-known amplitude equations by using the standard multiple-scale analysis [33–42], and references therein. Since the methods used are standard, we omit the detailed process and only give the main results. The readers can refer to [33] for more details on the derivation process. In what follows, we will briefly divide into three steps to obtain the amplitude equations. Step 1: we write out the linearized form of model (1.5) at the positive equilibrium point (u∗ , v ∗ ) as follows:
1 − 9β 2 2 1 1 − β2 3 1 2 ∂u 3 2 2 ∂ t = a11 u + a12 v + 8β 2 u − 2β uv − 16β 4 u + 8β 3 u v + o(ρ ) + d11 ∇ u + d12 ∇ v, γ γ (1 − β 2 ) 3 γ 2 ∂v γ (1 − β 2 ) 2 = a21 u + a22 v − u + uv + u − u v + o(ρ 3 ) + d21 ∇ 2 u + d22 ∇ 2 v. 2 4 ∂t 8β 2β 16β 8β 3
(2.6)
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
41
Fig. 5. (Color online) Bifurcation diagram of model (1.5). We set the parameter values as d11 = d12 = 1, d22 = 5. Turing and Hopf instability boundaries varying γ . The instabilities stay below the lines.
Fig. 6. The relation between Re(λ) (the real part of the eigenvalue λ) and k with d11 = d12 = 1, d22 = 5, β = 0.6780, γ = 0.4 and different d12 . Curve (1 )
(2)
(3)
(4 )
(5)
(6)
(1): d21 = 2.678; Curve (2): d21 = 2.850; Curve (3): d21 = 3.358; Curve (4): d21 = 4.000; Curve (5): d21 = 4.508; Curve (6): d21 = 4.952.
and setting U = (u, v)T , model (1.5) can be converted to the following system:
∂ U = L U + N, ∂t
(2.7)
where
L = LT + (βT − β)M =
1 − 9β 2
8β 2
N=
u2 −
aT11 + d11 ∇ 2
aT12 + d12 ∇ 2
aT21 + d21 ∇ 2
aT22 + d22 ∇ 2
1 2β
uv −
1 − β2 16β 4
u3 +
1 8β
+ (βT − β)
u2 v + o(ρ 3 ) 3
m11 m21
m12 m22
,
, γ γ 2 γ (1 − β ) 2 γ (1 − β ) 3 3 u + u − u v + o (ρ ) − u v + 8β 2 2β 16β 4 8β 3 2
2
where
−
aTij
= aij (βT )(i, j = 1, 2), m11 = 3β , m12 = 1, m21 = γ β , m22 = 0, ( ) = T
aT12 −(k2T )∗ d12 aT11 −(k2T )∗ d11
aT22 −(k2T )∗ d22 aT21 −(k2T )∗ d21
,g = − 1 L+ =0 T
=− f LT = 0, 1
aT11 −(k2T )∗ d11 aT21 −(k2T )∗ d21
T
=−
g
+
respectively, where LT is the adjoint operator of LT .
aT12 −k2T d12 aT22 −(k2T )∗ d22
k2T ∗
k2T
(βT ) =
satisfy the following equalities
aT11 aT22 −aT12 aT21 d11 d22 −d12 d21
, and f =
42
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
Step 2: Close to onset β = βT , we compute the solutions of model (1.5) that can be expanded as in the following asymptotic expansion:
u1
U=ε
v1
+ ε2
u2
v2
+ ε3
u3
v3
+ o(ε 3 ),
where
u1
=
v1
3 f 1
Wj exp(ikj · r) + c .c . ,
j =1
and
u2
U0 V0
=
v2
+
3 U j
j =1
+
U12 V12
Vj
exp(ikj · r) +
3 U jj
Vjj
j=1
exp(i(k1 − k2 ) · r) +
U23 V23
exp(i2kj · r)
exp(i(k2 − k3 ) · r) +
U31 V31
exp(i(k3 − k1 ) · r) + c .c .
where
U0 V0
u00
(|W1 |2 + |W2 |2 + |W3 |2 ), Uj = fVj , Ujj u11 Ujl u 2 = Wj , = ∗ Wj W l , Vjj v11 Vjl v∗ 2 γ γ (1−β 2 ) 2 1−9β aT22 β1 f − 4β 2 f 2 − aT12 f − f 2 β 4β T T T T a a − a a u00 11 22 12 21 = , v00 aT γ (1−β 2 ) f 2 − γ f − aT 1 f − 1−9β 2 f 2 11 21 β β 4β 2 4β 2 =
v00
aT11 aT22 − aT12 aT21
u11
v11
(aT22 − 4d22 (k2T )∗ )
1 f 2β
−
1−9β 2 2 f 8β 2
− (aT12 − 4d12 (k2T )∗ )
γ (1−β 2 ) 2 f 8β 2
−
γ
2β
f
(aT − 4d (k2 )∗ )(aT − 4d (k2 )∗ ) − (aT − 4d (k2 )∗ )(aT − 4d (k2 )∗ ) 11 T 22 T 12 T 21 T 11 22 12 21 = , T γ γ (1−β 2 ) 2 1−9β 2 2 1 T 2 ∗ 2 ∗ f − 2β f − (a12 − 4d12 (kT ) ) 2β f − 8β 2 f (a11 − 4d11 (kT ) ) 8β 2
(aT11 − 4d11 (k2T )∗ )(aT22 − 4d22 (k2T )∗ ) − (aT12 − 4d12 (k2T )∗ )(aT21 − 4d21 (k2T )∗ ) 2) 2 γ 2 (aT22 − 3d22 (k2T )∗ ) β1 f − 1−4β92β f 2 − (aT12 − 3d12 (k2T )∗ ) γ (14−β f − f 2 β β (aT − 3d (k2 )∗ )(aT − 3d (k2 )∗ ) − (aT − 3d (k2 )∗ )(aT − 3d (k2 )∗ ) 11 T 22 T 12 T 21 T u∗ 11 22 12 21 = . v∗ T γ (1−β 2 ) 2 γ 1−9β 2 2 1 2 ∗ T 2 ∗ f − β f − (a12 − 3d12 (kT ) ) β f − 4β 2 f (a11 − 3d11 (kT ) ) 4β 2
(aT11 − 3d11 (k2T )∗ )(aT22 − 3d22 (k2T )∗ ) − (aT12 − 3d12 (k2T )∗ )(aT21 − 3d21 (k2T )∗ ) Step 3: By the above expressions, together with Aj = ε Wj + ε 2 Vj + o(ε 3 ), we can get the amplitude equation corresponding to A1 , A2 , A3 as follows:
∂A 1 τ0 = µA1 + hA2 A3 − [g1 |A1 |2 + g2 (|A2 |2 + A3 |2 )]A1 , ∂t ∂ A2 τ0 = µA2 + hA1 A3 − [g1 |A2 |2 + g2 (|A1 |2 + A3 |2 )]A2 , ∂ t ∂ A3 τ = µA3 + hA1 A2 − [g1 |A3 |2 + g2 (|A1 |2 + A2 |2 )]A3 , 0 ∂t
(2.8)
where
τ0 = g1 =
f +g , βT [fm11 + m12 + g (fm21 + m22 )] G1
βT [fm11 + m12 + g (fm21 + m22 )]
,
µ= g2 =
βT − β , βT
h=
H
βT [fm11 + m12 + g (fm21 + m22 )]
G2
βT [fm11 + m12 + g (fm21 + m22 )]
,
,
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
H =
1 − 9β 2 4β
2
f2 −
1 − 9β 2
1
β
f +g 1
−
γ (1 − β 2 ) 2 γ f + f 4β 2 β
43
,
1
3
(u00 + u11 ) − f (v00 + v11 ) + f2 4β 2β 2β 8β 3 3(1 − β 2 ) 3 γ 3γ 2 γ (1 − β 2 ) γ 3γ (1 − β 2 ) 3 − ( u + u ) + f (v + v ) − f + g − f + f + f , 00 11 00 11 16β 4 4β 2 2β 2β 8β 3 16β 4 1 3(1 − β 2 ) 3 1 3 2 1 − 9β 2 f − (u00 + u∗ ) − f (v00 + v∗ ) + f − f −G2 = 2 3 4β 2β 2β 4β 8β 4 γ (1 − β 2 ) γ γ 3γ 2 3γ (1 − β 2 ) 3 +g − f + ( u + u ) + f (v + v ) − f + f . 00 ∗ 00 ∗ 4β 2 2β 2β 4β 3 8β 4
−G1 =
2
f −
2.2.1. Amplitude instability Each amplitude in Eq. (2.8) can be decomposed to mode ρj = |Aj | and a corresponding phase angle ϕj . Then, substituting Aj = ρj exp(iϕj ) into Eq. (2.8) and separating the real and imaginary parts, we can get four differential equations of the real variables as follows:
∂ϕ ρ 2 ρ 2 + ρ12 ρ32 + ρ22 ρ32 = − h 1 2 sin ϕ, τ 0 ∂t ρ1 ρ2 ρ3 ∂ρ τ0 1 = µρ1 + hρ2 ρ3 cos ϕ − g1 ρ13 − g2 (ρ22 + ρ32 )ρ1 , ∂t ∂ρ 2 = µρ2 + hρ1 ρ3 cos ϕ − g1 ρ23 − g2 (ρ12 + ρ31 )ρ2 , τ0 ∂t ∂ρ3 = µρ3 + hρ1 ρ2 cos ϕ − g1 ρ33 − g2 (ρ12 + ρ22 )ρ3 , τ0 ∂t where ϕ = ϕ1 + ϕ2 + ϕ3 .
(2.9)
The dynamical system (2.9) possesses five kinds of solutions in [33,34]. (1) The stationary state (O), given by
ρ1 = ρ2 = ρ3 = 0, is stable for µ < µ2 = 0 and unstable for µ > µ2 . (2) Stripe patterns (S ), given by
ρ1 =
µ g1
̸= 0,
ρ2 = ρ3 = 0, h2 g
is stable for µ > µ3 = (g −g1 )2 and unstable for µ < µ3 . 2 1 (3) Hexagon patterns (H0 , Hπ ) are given by
ρ1 = ρ2 = ρ3 =
| h| ±
h2 + 4(g1 + 2g2 )µ
2(g1 + 2g2 )
, √
2
with ϕ = 0 or π , and exist when µ > µ1 = 4(g−+h2g ) . The solution ρ + = 1 2 √ |h|−
h2 +4(g1 +2g2 )µ 2(g1 +2g2 )
and ρ − = (4) The mixed states are given by 2g1 +g2 h2 , (g2 −g1 )2
| h| ρ1 = , g2 − g1
ρ2 = ρ3 =
|h|+
h2 +4(g1 +2g2 )µ 2(g1 +2g2 )
is stable only for µ < µ4 =
is always unstable.
µ − g1 ρ12 , g1 + g2
with g2 > g1 , µ > g1 ρ12 and are always unstable. With fixed parameters d11 = d12 = 1, d21 = 11, d22 = 15, and γ = 2, we can find that h = 0.7958, g1 = 8.8043, g2 = 3.7966, µ1 = −0.0097, µ2 = 0, µ3 = 0.2223, µ4 = 0.5405. Summarizing the above analysis, we can show our results in Fig. 7. The system has a bistable region for µ ∈ (µ1 , µ2 ). In other words, when the control parameter µ lies in this region, the H0 patterns and the stationary state are all stable. The Hπ patterns are always unstable when µ > µ2 . When µ lies in region (µ2 , µ3 ), the stripe patterns are unstable, and the H0 patterns are stable. In region (µ3 , µ4 ), the system exists in another bistable state (i.e., the bistable state between the hexagon patterns and the stripe patterns). When the parameter µ is more than µ4 , stripe patterns emerge in the two-dimensional space.
44
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
Fig. 7. Bifurcation diagram of model (1.5) with d11 = d12 = 1, d21 = 11, d22 = 15, and γ = 2. H0 : hexagonal patterns with ϕ = 0; Hπ : hexagonal patterns with ϕ = π ; S: stripe patterns. Solid lines: stable states; dashed lines: unstable states. µ1 = −0.0097, µ2 = 0, µ3 = 0.2223, µ4 = 0.5405.
3. Pattern formation of model (1.5) In this section, we perform numerical simulations of model (1.5) in two dimensional space and the simulated results are shown here. To solve differential equations by computer, one has to discretize the space and time of the problem. The continuous problem defined by the reaction–diffusion system in two-dimensional space is solved in a discrete domain with M × N lattice sites. The spacing between the lattice points is defined by the lattice constant ∆h. In the discrete system the Laplacian describing diffusion is calculated using finite differences, i.e., the derivatives are approximated by differences over ∆h. The time evolution is also discrete, i.e., the time goes in steps of ∆t. The time evolution can be solved by using the Euler method. In the present paper, we set ∆h = 0.5, different ∆t values: 0.001 or 0.01, and M = N = 200. All of our numerical simulations employ the Neumann boundary conditions. And the initial data is all taken as a uniformly distributed random perturbation around the equilibrium (u∗ , v ∗ ). Through the Introduction section, we know that β is proportional to the death rate s, which plays a vital role in the spatial predator–prey model. So, from the motion point of view, we need to consider the effect of β on pattern dynamics of model (1.5), and also consider the effect of d21 on pattern dynamics of model (1.5). In what follows, we preform the numerical simulations from the two aspects of the varied β and varied cross-diffusion d21 , respectively. 3.1. Pattern formation on the effect of the varied β We first keep d11 = d12 = 1, d21 = 11, d22 = 15, and γ = 2. Then we can obtain g1 , g2 , h, µ1 , µ2 , µ3 , µ4 . In the numerical simulations, different types of patterns are observed and the distributions of predator and prey are always of the similar type. As a result, we only show our results of pattern formation about one species distribution (in this paper, we show the distribution of prey, for instance). For the different values of β located in the ‘‘Turing space’’ (the domain D11 in Fig. 1), we show and obtain three categories of Turing patterns for the distribution of prey u of model (1.5). In every pattern, the red (blue) represents the high (low) density of prey u. In Fig. 8, for prey u, Hπ hexagon patterns prevail over the whole domain eventually. Parameter values (β = 0.712) set in this figure satisfy µ = 0.07424262 ∈ (µ2 , µ3 ). According to the analysis above, there should be only Hπ hexagon patterns under this circumstance. In other words, the numerical simulation is corresponding to the theoretical analysis. We should also pay attention to the situation that µ is very close to the critical value µ2 . In this case, the Hπ hexagon patterns come into being very slowly. This is the universal phenomenon of critical slowing down. On the other hand, Fig. 8(f) consists of blue (minimum density of u) hexagons on a red (maximum density of u) background, i.e., the prey is isolated with low population density. We call this spot pattern. It exhibits a competition between stripes and spots. The pattern takes a long time to settle down, starting with a homogeneous state (u∗ , v ∗ ) = (0.5069, 0.3511) (cf. Fig. 8(a)), and the random perturbations lead to the formation of stripes and spots (cf. Fig. 8(b)–(e)), and ending with spots only (cf. Fig. 8(f)). In Fig. 9, parameter values (β = 0.619) set in this figure satisfy µ = 0.19516318 ∈ (µ2 , µ3 ). A few stripes emerge, and the remaining part of the spots pattern remains time independent, i.e., spots–stripes pattern. It exhibits a competition between stripes and spots. Starting with a homogeneous state (u∗ , v ∗ ) = (0.3832, 0.3818) (cf. Fig. 9(a)), the random perturbations lead to the formation of spots–stripes (cf. Fig. 9(b)–(e)), ending with the time-independent spots–stripes pattern (cf. Fig. 9(f)). In this case, the numerical simulations cannot correspond to the theoretical analysis. This phenomenon cannot be explained by the amplitude equations, because the control parameter β is far away from the critical value β = 0.7691. In Fig. 10, Parameter values (β = 0.5778) set in this figure satisfy µ = 0.2458718 ∈ (µ3 , µ4 ). Model dynamics exhibits a transition from spots–stripes growth to stripes replication, i.e., spots decay and the stripes pattern emerges. Starting with a homogeneous state (u∗ , v ∗ ) = (0.3339, 0.3849) (cf. Fig. 10(a)), the random perturbations lead to the formation
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
45
Fig. 8. Snapshots of contour pictures of the time evolution of the prey at different instants with d11 = d12 = 1, d21 = 11, d22 = 15, γ = 2, β = 0.712 and ∆t = 0.001. Moments: (a) t = 0; (b) t = 5100; (c) t = 7500; (d) t = 10 500; (e) t = 15 000; (f) t = 25 000.
Fig. 9. Snapshots of contour pictures of the time evolution of the prey at different instants with d11 = d12 = 1, d21 = 11, d22 = 15, γ = 2, β = 0.619 and ∆t = 0.001. Moments: (a) t = 0; (b) t = 2500; (c) t = 4000; (d) t = 10 000; (e) t = 15 000; (f) t = 25 000.
46
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
Fig. 10. Snapshots of contour pictures of the time evolution of the prey at different instants with d11 = d12 = 1, d21 = 11, d22 = 15, γ = 2, β = 0.5778 and ∆t = 0.001. Moments: (a) t = 0; (b) t = 2300; (c) t = 5000; (d) t = 9000; (e) t = 15 000; (f) t = 25 000.
Fig. 11. Snapshots of contour pictures of the time evolution of the prey at final instants. (a): β = 0.5713, (b): β = 0.5210, (c): β = 0.4712 with d11 = d12 = 1, d21 = 11, d22 = 15, γ = 2 and ∆t = 0.001. Moments: (a) t = 25 000; (b) t = 25 000; (c) t = 9900.
of spots–stripes (cf. Fig. 10(b)–(e)), and the later random perturbations make these spots decay, ending with the timeindependent stripes pattern (cf. Fig. 10(f)). In this case, the numerical simulations cannot also correspond to the theoretical analysis. This phenomenon cannot also be explained by the amplitude equations, because the control parameter β is far away from the critical value β = 0.7691. Finally, fixing d11 = d12 = 1, d21 = 11, d22 = 15, we choose some points in the Hopf and Turing instability region D13 in Fig. 1 to simulate, that is choosing four types of points: (I): β = 0.5713, 0.5210, 0.4712, γ = 2; (II): β = 0.528, 0.485, 0.380, γ = 7; (III): two points next to Turing–Hopf point: β = 0.572, γ = 5.75 and β = 0.535, γ = 6.25; (IV): one point far from Turing–Hopf point and between type (I) and type (II) points: β = 0.435, γ = 4.68. And we obtain some interesting simulation results, (a)–(c) in Fig. 11, and (a) and (b) in Fig. 13 are called as ‘‘labyrinth’’ patterns and the stripe-like patterns prevail over the whole domain (cf. Fig. 12(a)), and that the stripe-spotted patterns coexist in the space (cf. Fig. 12(b) and Fig. 13(c)), spotted patterns in the space (cf. Fig. 12(c)) which are different from the above results. Of course, we have to mention the fact that, due to the Hopf and Turing instability, the resulting pattern could exhibit oscillations and they could be so strong that our numerical scheme would be not convergent.
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
47
Fig. 12. Snapshots of contour pictures of the time evolution of the prey at final instants. (a): β = 0.528, (b): β = 0.485, (c): β = 0.380 with d11 = d12 = 1, d21 = 11, d22 = 15, γ = 7 and ∆t = 0.001. Moments: (a) t = 8800, (b) t = 10 000; (c) t = 9000.
Fig. 13. Snapshots of contour pictures of the time evolution of the prey at final instants. (a): β = 0.572, γ = 5.75, (b): β = 0.535, γ = 6.25, (c): β = 0.435, γ = 4.68 with d11 = d12 = 1, d21 = 11, d22 = 15 and ∆t = 0.001. Moments: (a) t = 13 000, (b) t = 8300; (c) t = 10 500.
3.2. Pattern formation on the effect of varied d21 Next, we continue to keep d11 = d12 = 1, d22 = 5, γ = 0.4, β = 0.678, and let d21 vary. In the 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 still restrict our analysis of pattern formation to the distribution of prey. In Fig. 14, we show the snapshots of prey spatial pattern at different values of d21 with d11 = d12 = 1, d22 = 15, γ = 0.4, β = 0.678 and 25,000 iterations. With small random perturbation of the stationary solution (u∗ , v ∗ ) = (0.4597, 0.3663) of the spatially homogeneous systems with d21 = 4.952, one can see that the stripe-like patterns prevail over the whole domain (cf. Fig. 14(e)), and that the stripe-like and spotted patterns coexist in the space (cf. Fig. 14(d)) when d21 is being decreased to 4.508. However, as d21 is being decreased to 4.000, 3.358 and 2.850, we find that the spotted patterns prevail over the whole domain, the color varies from crimson (blue), red (blue) to green (blue) (cf. Fig. 14(c)–(a)). If we continue to decrease to d21 = 2.678, the system (1.5) cannot induce spatial pattern. Therefore, we find that when the cross diffusion d21 decreases, the modes of the patterns convert from the stripes to the spots. Finally, we choose some points in the Hopf and Turing instability regions D23 in Fig. 4 to simulate, that is choosing d21 = 3.0, 3.8, 4.5, β = 0.2, γ = 0.4, d11 = d12 = 1, d22 = 5. And we obtain some interesting simulation results, (a): large yellow spots, (b): large dark blue dots, (c): small dark blue dots, in Fig. 15, which are different from the above results. Of course, we have to mention the fact that, owing to the interaction of Hopf and Turing instabilities, the resulting pattern could exhibit oscillations and they could be so strong that our numerical scheme would be not convergent and there exist oscillatory patterns. Through the numerical simulations of u above, we can know that the distributions of prey and predators are always of the same type. This means that there is a tendency that the prey stays away from the predators and the predators will get closer to the prey, which reflects the effect of the cross-diffusion. What is more, the local density of prey is higher than that of predators. This is the precondition of the coexistence of prey and predators. 4. Discussion and conclusion In this paper, we deduce the amplitude equations and present the Turing pattern selection of a predator–prey model with cross diffusion. By using mathematical analysis and numerical simulations, we found that this model can give rise to all three categories: spotted, spots–stripes mixtures, and stripe-like patterns and other patterns: ‘‘labyrinth’’ patterns and
48
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
Fig. 14. Snapshots of contour pictures of the evolution of the prey at different values of d21 with d11 = d12 = 1, d22 = 5, γ = 0.4, β = 0.678. (a): d21 = 2.850; (b): d21 = 3.358; (c): d21 = 4.000; (d): d21 = 4.508; (e): d21 = 4.952 and ∆t = 0.01. Moments: (a) t = 60 000; (b) t = 7500; (c) t = 8000; (d) t = 3500; (e) t = 2500.
Fig. 15. Snapshots of contour pictures of the time evolution of the prey at final instants. (a): d21 = 3.0, (b): d21 = 3.8, (c): d21 = 4.5 with β = 0.2,
γ = 0.4, d11 = d12 = 1, d22 = 5 and ∆t = 0.01. (a) t = 2000; (b) t = 10 500; (c) t = 10 000.
‘‘black-eye’’ patterns. That is to say, cross diffusion can create stationary patterns. It should be noted that, as cross-diffusion coefficients d12 = d21 = 0 and the predator mortality described by the linear form, the spatial predator–prey model cannot give rise to Turing structures that has been pointed out in [4], which reflects the effect of the cross-diffusion. From biological point of view, our results have more biological meanings. In a word, our results are different from [4] and well enrich the findings in predator–prey models with herd behavior. By the above analysis, we can find that the qualitative dynamics of the model (1.5) are fundamentally different when parameter β and cross-diffusion coefficient d21 in model (1.5) slightly change. The parameter β is proportional to the (non-dimensional) death rate of the predator. From the biological point of view, our results show that the death rate of predator may play a vital role in the spatial predator–prey model. By varying the value of β , we obtain three different typical types of pattern: spot pattern (Fig. 8(f)), spot–stripe pattern (Fig. 9(f)), and stripe pattern (Fig. 10(f)). From the view of population dynamics, one can see that there exists the spot pattern replication—the prey u is the isolated zone with low density, and the remaining region is high density, which means the prey may break out in the area in Fig. 8(f). In other words, the prey in this area is safe. The biological significance of the other cases can be determined in the same way as the above case. On the other hand, we only let d21 vary, that is decreasing the value of d21 from 4.952 to 2.850, and find that the system (1.5) gives rise to all three categories: stripe pattern (Fig. 14(e)), spots–stripes mixtures (Fig. 14(d)), spot pattern (Fig. 14(c–a)). But when the value of d21 decreases to 2.678, the system (1.5) cannot induce spatial pattern. In Figs. 11, 13 and 15, we also obtain other patterns: ‘‘labyrinth’’ patterns and ‘‘black-eye’’ patterns, respectively.
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
49
Of course, the methods and the results in the present paper can be applied to the pattern formation in the diffusive predator–prey model or other reaction–diffusion systems. We hope that our work could be instructive to study the effect of cross-diffusion on the pattern formation in population dynamics. Acknowledgments The authors would like to thank the referee for his/her careful reading and some comments on improving the presentation of this paper. The work was supported by the Shanghai Committee of Science and Technology (No. 13ZR1444500) and the Program for New Century Excellent Talents in University (NCET-11-0385). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42]
V. Volterra, Sui tentutive di applicazione delle mathematiche alle seienze biologiche e sociali, Ann. Radioelectr. Univ. Romandes (1901-2) 3–28. V. Volterra, Variazione e fluttuazini del numero d’individui in specie animali conviventi, Mem. R. Accad. Naz. dei Lincei 2 (1926) 31–113. P.A. Braza, Predator–prey dynamics with square root functional responses, Nonlinear Anal. RWA 13 (2012) 1837–1843. S. Yuan, C. Xu, T. Zhang, Spatial dynamics in a predator–prey model with herd behavior, Chaos 23 (2013) 0331023. A.M. Turing, The chemical basis of morphogenesis, Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. 237 (1952) 37–72. J.D. Murray, Mathematical Biology, Springer-Verlag, Berlin, 1993. S. Aly, I. Kim, D. Sheen, Turing instability for a ratio-dependent predator–prey model with diffusion, Appl. Math. Comput. 217 (2011) 7265–7281. L.N. Guin, P.K. Mandal, Spatial pattern in a diffusive predator–prey model with sigmoid ratio-dependent functional response, Int. J. Biomath. 7 (2014) 1–26. 1450047. Y. Song, X. Zou, Spatiotemporal dynamics in a diffusive ratio-dependent predator–prey model near a Hopf–Turing bifurcation point, Comput. Math. Appl. 67 (2014) 1978–1997. Y. Song, X. Zou, Bifurcation analysis of a diffusive ratio-dependent predator–prey model, Nonlinear Dynam. 78 (2014) 49–70. Z. Xu, Y. Song, Bifurcation analysis of a diffusive predator–prey system with a herd behavior and quadratic mortality, Math. Meth. Appl. Sci. (2014). http://dx.doi.org/10.1002/mma.3275. J. Shukla, S. Verma, Effects of convective and dispersive interactions on the stability of 2 species, Bull. Math. Biol. 43 (1981) 593–610. E.H. Kerner, A statistical mechanics of interacting biological species, Bull. Math. Biol. 19 (1957) 121–146. N. Shigesada, K. Kawasaki, E. Teramoto, Spatial segregation of interacting species, J. Theoret. Biol. 79 (1979) 83–99. W. Wang, Q. Liu, Z. Jin, Spatiotemporal complexity of a ratio-dependent predator–prey system, Phys. Rev. E 75 (2007) 051913. K. Kuto, Stability of steady-state solutions to a prey–predator system with cross-diffusion, J. Differential Equations 197 (2004) 293–314. K. Kuto, Y. Yamada, Multiple coexistence states for a prey–predator system with cross-diffusion, J. Differential Equations 197 (2004) 315–348. V.K. Vanag, I.R. Epstein, Cross-diffusion and pattern formation in reaction–diffusion systems, Phys. Chem. Chem. Phys. 11 (2009) 897–912. X. Zhang, G. Sun, Z. Jin, Spatial dynamics in a predator–prey model with Beddington–DeAngelis functional response, Phys. Rev. E 85 (2012) 021924. L. Xue, Pattern formation in a predator–prey model with spatial effect, Physica A 391 (2012) 5987–5996. G. Gambino, M.C. Lombardo, M. Sammartino, Pattern formation driven by cross-diffusion in a 2D domain, Nonlinear Anal. RWA 14 (2013) 1755–1779. L.N. Guin, Existence of spatial patterns in a predator–prey model with self- and cross-diffusion, Appl. Math. Comput. 226 (2014) 320–335. G. Galiano, M.L. Garzon, A. Jüngel, Semi-discretization in time and numerical convergence of solutions of a nonlinear cross-diffusion population model, Numer. Math. 93 (2003) 655–673. L.N. Guin, P.K. Mandal, Spatiotemporal dynamics of reaction–diffusion models of interacting populations, Appl. Math. Model. 38 (2014) 4417–4427. C. Tian, Z. Lin, M. Pedersen, Instability induced by cross-diffusion in reaction–diffusion systems, Nonlinear Anal. RWA 11 (2010) 1036–1045. C. Tian, Z. Ling, Z. Lin, Turing pattern formation in a predator–prey-mutualist system, Nonlinear Anal. RWA 12 (2011) 3224–3237. C. Tan, L. Zhang, Z. Lin, Pattern formation for a model of plankton allelopathy with cross-diffusion, J. Franklin Inst. 348 (2011) 1947–1964. R.R. Baier, C. Tian, Mathematical analysis and numerical simulation of pattern formation under cross-diffusion, Nonlinear Anal. RWA 14 (2013) 601–612. X. Wang, Y. Cai, Cross-diffusion-driven instability in a reaction–diffusion harrison predator–prey model, Abst. Appl. Anal. 2013 (2013) 306467. W. Wang, Z. Guo, R.K. Upadhyay, Y. Lin, Pattern formation in a cross-diffusive Holling–Tanner model, Discrete Dyn. Nat. Soc. 2012 (2012) 828219. M.C. Cross, P.C. Hohenberg, Pattern formation outside of equilibrium, Rev. Modern Phys. 65 (1993) 851–1112. G. Gambino, M.C. Lombardo, M. Sammartino, Turing instability and pattern formation for the Lengyel–Epstein system with nonlinear diffusion, Acta Appl. Math. 132 (2014) 283–294. Q. Ouyang, Nonlinear Science and the Pattern Dynamics Introduction, Peking University Press, Beijing, 2010. Y. Kuramoto, T. Tsuzuki, On the formation of dissipative structures in reaction–diffusion systems, Progr. Theoret. Phys. 54 (1975) 687–699. G.H. Gunaratne, Q. Ouyang, H.L. Swinney, Pattern formation in the presence of symmetries, Phys. Rev. E 50 (1994) 2802–2820. B. Peña, C. Pérez-García, Stability of Turing patterns in the Brusselator model, Phys. Rev. E 64 (2001) 056213. W. Wang, Y. Lin, L. Zhang, F. Rao, Y. Tan, Complex patterns in a predator–prey model with self and cross-diffusion, Commun. Nonlinear Sci. Numer. Simul. 16 (2011) 2006–2015. G. Sun, Z. Jin, Q. Liu, L. Li, Pattern formation induced by cross-diffusion in a predator–prey system, Chin. Phys. B 17 (2008) 3936–3941. W. Wang, Y. Lin, F. Rao, L. Zhang, Y. Tan, Pattern selection in a ratio-dependent predator–prey model, J. Stat. Mech. 11 (2010) P11036. F. Rossi, V.K. Vanag, I.R. Epstein, Pentanary cross-diffusion in water-in-oil microemulsions loaded with two components of the Belousov–Zhabotinsky reaction, Chem. Eur. J. 17 (2011) 2138–2145. G. Gambino, M.C. Lombardo, M. Sammartino, V. Sciacca, Turing pattern formation in the Brusselator system with nonlinear diffusion, Phys. Rev. E 88 (4) (2013) 042925. Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence, Springer-Verlag, Berlin, 1984.
Contents lists available at ScienceDirect
Nonlinear Analysis: Real World Applications journal homepage: www.elsevier.com/locate/nonrwa
Cross-diffusion induced spatiotemporal patterns in a predator–prey model with herd behavior Xiaosong Tang a,b , Yongli Song a,∗ a
Department of Mathematics, Tongji University, Shanghai 200092, China
b
College of Mathematical and Physics, Jinggangshan University, Ji’an 343009, China
article
info
Article history: Received 15 August 2014 Received in revised form 25 December 2014 Accepted 31 December 2014
Keywords: Predator–prey model Herd behavior Cross-diffusion Turing instability Amplitude equations Spatiotemporal pattern
abstract In this paper, we investigate the phenomena of pattern formation for a predator–prey model with herd behavior and cross-diffusion. We give the conditions for cross-diffusion induced instability in detail, and have derived amplitude equations for the excited modes. Furthermore, we illustrate the spatial patterns via numerical simulations, which show that the model dynamics exhibits complex pattern replication: spotted pattern, stripe pattern, and coexistence of the two by increasing the control parameter β and cross-diffusion coefficient d21 respectively. Meanwhile, we also obtain the other patterns: ‘‘labyrinth’’ patterns and ‘‘black-eye’’ patterns. These results indicate that the cross-diffusion plays an important role in the pattern selection. The spatiotemporal patterns induced by cross-diffusion may enrich the pattern dynamics for diffusive predator–prey model and be helpful for the biologists to understand the phenomenon in which a flux of one species is induced by a gradient of another species. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction In ecological systems, the interactions of different species within a fluctuating natural environment are the most distinctive features. Studies in mathematical models are informative in understanding the interactions of predator and prey in these systems. As a result, predator–prey models have been important in ecological science since the pioneering works of Lotka and Volterra [1,2]. However, these models have the unavoidable limitations to describe many realistic phenomena in biology. Recently, a predator–prey model is considered in which the prey exhibits herd behavior, so that the predator interacts with the prey along the outer corridor of the herd of prey, which is more appropriate to model the response functions of prey that exhibit herd behavior in terms of the square root of the prey population, due originally to [3,4]:
√ dX = rX 1 − X − α X Y√ , dt K √ 1 + th α X (1.1) cα X Y dY = −sY + √ , dt 1 + th α X where X (t ) and Y (t ) stand for the prey and predator densities at time t, respectively. The parameter r is the growth rate of the prey, K is its carrying capacity, and s is the death rate of the predator in the absence of prey. The parameter α is the search efficiency of Y for X , c is the conversion or consumption rate of prey to predator, and th is average handling time. ∗
Corresponding author. E-mail address: [email protected] (Y. Song).
http://dx.doi.org/10.1016/j.nonrwa.2014.12.006 1468-1218/© 2015 Elsevier Ltd. All rights reserved.
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
37
For simplicity, we nondimensionalize model (1.1) with the following scaling: u=
X K
αY v= √ ,
,
r
K
tnew = rtold ,
√
rL1 ,
r1 =
√ r2 =
rL2 ,
and the other parameters are made dimensionless as follows:
√ snew =
K
α
sold ,
cnew =
√ α K r
cold ,
a = th α
√
K.
Then model (1.1) converts into
√ uv du √ , dt = u(1 − u) − 1+a u √ dv c uv = −sv + √ . dt 1+a u
(1.2)
On the other hand, the issue of spatial and temporal pattern formation in biological communities is probably one of the most exciting problems in modern biology and ecology. Spatial pattern formation arose from the observation in chemistry by Turing [5] that diffusion can also destabilize equilibrium solution, a scenario well known as Turing instability. Although this is counter-intuitive, the interest in the diffusion-driven instability has long expanded from chemical system to biological system [6–12]. Thus, the following reaction–diffusion system is obtained as a model for predator–prey with herd behavior:
√ ∂u uv = u ( 1 − u ) − √ + d11 ∇ 2 u, ∂t 1+a u √ ∂v c uv = −sv + √ + d22 ∇ 2 v, ∂t 1+a u
(1.3)
where the nonnegative constants d11 and d22 , are the diffusion coefficients of prey and predator respectively. Nevertheless, from [4], we know that the phenomenon of spatial pattern formation in (1.3) cannot occur under all possible diffusion rates. So, the authors in [4,11] change linear mortality sv into quadratic mortality sv 2 in (1.3) and have investigated Turing patterns, stability, Turing instability and Hopf bifurcations. However, the prey may be able to recognize the predator, and thus the prey would keep away from predators in order to not to be caught. In nature, the escape velocity of the prey might be proportional to the diffusive velocity of the predators. On the other hand, the tendency of predators would get closer to the prey, and hence the chase velocity of predators may be considered to be proportional to the diffusive velocity of the prey [13]. This phenomenon is called cross diffusion. Therefore the problem of cross-diffusion arises, which was proposed first by Kerner [14] and first applied to competitive population systems by Shigesada et al. [15]. Recently, there has been increasing interest in the stability behavior of predator–prey systems, in which both self-diffusion and cross-diffusion effects are considered [16–23], and references therein. The cross-diffusion describes the population flux of one species due to the presence of other species. Unfortunately, most studies on the spatial and temporal predator–prey systems pay little attention to the study on the effect of cross-diffusion and the selection of Turing patterns. They just focus on the bifurcation phenomena by varying the controlled parameter or only give some numerical results of the patterns [24–32]. Based on the previous work of many scientists, we can now focus on the following predator–prey model with herd behavior and cross-diffusion:
√ uv ∂u √ + d11 ∇ 2 u + d12 ∇ 2 v, ∂ t = u(1 − u) − 1+a u √ c uv ∂v = −sv + √ + d21 ∇ 2 u + d22 ∇ 2 v. ∂t 1+a u
(1.4)
In [3], the author has used the simplifying assumption that a = 0, which implies that the average handling time is zero. For the convenience of discussion and in line with the work of [3], we also assume that a = 0 and let γ = c , β = s/γ . Then Eq. (1.4) is simplified as
√ ∂u = u(1 − u) − uv + d11 ∇ 2 u + d12 ∇ 2 v, ∂t ∂v = γ v(−β + √u) + d21 ∇ 2 u + d22 ∇ 2 v, ∂t
(1.5)
where d12 and d21 , called cross-diffusion coefficients, describe the respective population fluxes of preys and predators resulting from the presence of the other species, respectively. Model (1.5) is to be analysed under the following non-zero initial conditions: u(x, y, 0) > 0,
v(x, y, 0) > 0,
(x, y) ∈ Ω = [0, L] × [0, L],
38
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
and zero-flux boundary conditions ∂u ∂v = = 0, ∂n ∂n where L denotes the size of the spatial variables x and y, n is the outward unit normal vector of the boundary ∂ Ω . The goal of this paper is to explore whether cross-diffusion can drive the emergence of Turing pattern and how different cross-diffusion coefficient effects the spatially inhomogeneous distribution of population density. Our analyses reveal that the cross-diffusion may induce spatial patterns and different cross-diffusion coefficient may lead to pattern transition between spots and stripes. The paper is organized as follows. In Section 2, we show the spatial dynamics of model (1.5) under the fixed crossdiffusion and varied cross-diffusion. In Section 3, we present and discuss the results of Turing pattern formation through extensive numerical simulations. Finally, conclusions and discussions are presented in Section 4. 2. Spatial dynamics of model (1.5) 2.1. Linear stability analysis In the absence of diffusion, simple calculations show that model (1.5) has three equilibrium points which consist of two boundary equilibriums (0, 0), (1, 0) and a positive equilibrium (u∗ , v ∗ ) = (β 2 , β(1 − β 2 )), β ∈ (0, 1). From the biological perspective, we are interested in studying the stability behavior of the positive equilibrium point. The Jacobian matrix corresponding to (u∗ , v ∗ ) is as follows:
J :=
a12 a22
a11 a21
1
2 2 (1 − 3β ) = 1 γ (1 − β 2 )
−β .
(2.1)
0
2 Following the standard linear analysis of the cross-diffusion system [33,34], we address the temporal stability of the uniform state to nonuniform perturbations
∗
u
v
u
=
v∗
+ε
uk
vk
eλt +ik·r + O(ε 2 ) + c .c .,
(2.2)
where k · k = k2 , k and λ are the wave number and frequency respectively; i is the imaginary unit and i2 = −1; r = (u, v) is the spatial vector in two-dimensional space and c .c. stands for the complex conjugate terms. The linear instability(ε ≪ 1) of the uniform state is deduced from the dispersion relations. We obtain that the eigenvalue is the root of the following equation:
λ2 − Tk λ + Dk = 0, (2.3) where Tk = a11 + a22 − (d1 + d2 )k2 , Dk = a11 a22 − a12 a21 − (a11 d22 + a22 d11 − a12 d21 − a21 d12 )k2 + (d11 d22 − d12 d21 )k4 . The roots of Eq. (2.3) can be obtained by the following form:
λk =
Tk ±
Tk2 − 4Dk
.
2 We choose β as the bifurcation parameter. The Hopf bifurcation occurs when Im(λk ) ̸= 0,
Re(λk ) = 0,
at k = 0,
where λk is the root of Eq. (2.3). Then we can get the critical value of the Hopf bifurcation parameter β :
√
βH =
3
.
(2.4)
3 The Turing bifurcation occurs when Im(λk ) = 0,
Re(λk ) = 0,
at k = kT ̸= 0,
and the wave number kT satisfies k2T =
βT :
a11 a22 −a12 a21 . d11 d22 −d12 d21
We can obtain the critical value of the Turing bifurcation parameter
(1 − 3β 2 )d22 + 2β d21 − (1 − β 2 )γ d12 − 2 2(d11 d22 − d12 d21 )βγ (1 − β 2 ) = 0.
(2.5)
Remark 2.1. When d12 = d21 = 0 in (1.5), a general linear analysis shows that the necessary conditions for yielding Turing √ patterns are given by T0 < 0, D0 > 0 and d22 a11 + d11 a22 − 2 d11 d22 D0 > 0. However, these conditions cannot hold simultaneously. √ In fact, T0 < 0 implies that a11 < 0 since a22 = 0. In addition, noticing that d11 > 0 and d22 > 0, we have d22 a11 + d11 a22 − 2 d11 d22 D0 < 0. So, the phenomenon of Turing pattern formation in (1.5) cannot occur when the cross-diffusion coefficients d12 and d21 are equal to zero. This implies that the cross-diffusion coefficients are crucial for the Turing instability.
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
39
Fig. 1. Bifurcation diagram of model (1.5). We set the parameter values as d11 = d12 = 1, d21 = 11, d22 = 15. The figure shows the Turing space D11 which is located below the positive equilibrium existence(the black curve), i.e. the region bounded by the Turing bifurcation curve (the blue (upper) curve) and the Hopf bifurcation curve (the red (lower) curve), D12 denotes Hopf space and D13 denotes Turing–Hopf space. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 2. (Color online) Bifurcation diagram of model (1.5). We set the parameter values as d11 = d12 = 1, d22 = 15. Turing and Hopf instability boundaries varying d21 . The instabilities stay below the lines.
Now, let us discuss the bifurcations represented by these formulas in the parameter region spanned by the parameters β and γ that can be seen from Fig. 1. The whole class of spatial model is included in this parameter region, where the domain D11 is called as Turing space, in which the Turing instability occurs, the domain D12 is called as Hopf space, in which the Hopf instability occurs and the region D13 is called Turing–Hopf space, where a competition between the two instabilities occurs. In D13 which one would develop, depends on the locations of the respective instability boundaries: as β decreases, if βT > βH , Turing instability occurs prior to the oscillatory instability. In Fig. 2, we report the Hopf and Turing instability boundaries for different values of d21 , and observe that, while the Hopf instability region does not depend on d21 , the Turing pattern region decreases by decreasing d21 due to the fact that βT decreases with d21 . Then we plot, in Fig. 3, the dispersion relations corresponding to several values of bifurcation parameter β , while keeping other parameters fixed. In Fig. 3, curve (2) corresponds to the critical Turing value β2 = 0.7691, curve (4) corresponds to the critical Hopf value β4 = 0.5773. On the other hand, we focus on how cross diffusion effects the dynamics for fixed deterministic parameters. To this end, we fix the deterministic model parameter values d11 = d12 = 1, d22 = 5, γ = 0.4, and vary d21 as a function of β , which is the coefficient of the cross diffusion of the prey. We show the bifurcation diagram of linear stability analysis in Fig. 4, where the domain D21 is called as Turing space, the domain D22 is called as Hopf space and the region D23 is called Turing–Hopf space. In D23 which one would develop, depends on the locations of the respective instability boundaries: as β increases, if βT > βH , Turing instability occurs prior to the oscillatory instability. In Fig. 5, we report the Hopf and Turing instability boundaries for different values of γ , and observe that, while the Hopf instability region does not depend on γ , the Turing pattern region increases by decreasing γ due to the fact that βT increases with γ . To see the Turing space properly, we plot the dispersion relation corresponding to several values of parameter d21 by keeping the β = 0.678. It can be seen from Fig. 6 that when d21 is decreased, the available Turing modes [Re(λ) > 0] are decreased, and all available modes are weakened.
40
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
Fig. 3. The relation between Re(λ) (the real part of the eigenvalue λ) and k with d11 = d12 = 1, d21 = 11, d22 = 15, γ = 2 and different β . Curve (1): β1 = 0.7820; Curve (2): β2 = 0.7691; Curve (3): β3 = 0.6250; Curve (4): β4 = 0.5773; Curve (5): β5 = 0.5560.
Fig. 4. Bifurcation diagram of model (1.5). We set the parameter values as d11 = d12 = 1, d22 = 5, γ = 0.4. The figure shows the Turing space D21 which is located below the positive equilibrium existence (the black curve), i.e. the region bounded by the Turing bifurcation curve (the blue (upper) curve) and the Hopf bifurcation curve (the red (lower) curve), D22 denotes Hopf space and D23 denotes Turing–Hopf space. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
From the above analysis, we can obtain the conditions of Turing instability for model (1.5), that is determining the conditions of the emergency of Turing patterns. But we cannot determine the selection and competition of Turing patterns. In the following section, we shall consider the selection of Turing patterns close to the onset β = βT = β2 .
2.2. Amplitude equations In this subsection, we can deduce the well-known amplitude equations by using the standard multiple-scale analysis [33–42], and references therein. Since the methods used are standard, we omit the detailed process and only give the main results. The readers can refer to [33] for more details on the derivation process. In what follows, we will briefly divide into three steps to obtain the amplitude equations. Step 1: we write out the linearized form of model (1.5) at the positive equilibrium point (u∗ , v ∗ ) as follows:
1 − 9β 2 2 1 1 − β2 3 1 2 ∂u 3 2 2 ∂ t = a11 u + a12 v + 8β 2 u − 2β uv − 16β 4 u + 8β 3 u v + o(ρ ) + d11 ∇ u + d12 ∇ v, γ γ (1 − β 2 ) 3 γ 2 ∂v γ (1 − β 2 ) 2 = a21 u + a22 v − u + uv + u − u v + o(ρ 3 ) + d21 ∇ 2 u + d22 ∇ 2 v. 2 4 ∂t 8β 2β 16β 8β 3
(2.6)
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
41
Fig. 5. (Color online) Bifurcation diagram of model (1.5). We set the parameter values as d11 = d12 = 1, d22 = 5. Turing and Hopf instability boundaries varying γ . The instabilities stay below the lines.
Fig. 6. The relation between Re(λ) (the real part of the eigenvalue λ) and k with d11 = d12 = 1, d22 = 5, β = 0.6780, γ = 0.4 and different d12 . Curve (1 )
(2)
(3)
(4 )
(5)
(6)
(1): d21 = 2.678; Curve (2): d21 = 2.850; Curve (3): d21 = 3.358; Curve (4): d21 = 4.000; Curve (5): d21 = 4.508; Curve (6): d21 = 4.952.
and setting U = (u, v)T , model (1.5) can be converted to the following system:
∂ U = L U + N, ∂t
(2.7)
where
L = LT + (βT − β)M =
1 − 9β 2
8β 2
N=
u2 −
aT11 + d11 ∇ 2
aT12 + d12 ∇ 2
aT21 + d21 ∇ 2
aT22 + d22 ∇ 2
1 2β
uv −
1 − β2 16β 4
u3 +
1 8β
+ (βT − β)
u2 v + o(ρ 3 ) 3
m11 m21
m12 m22
,
, γ γ 2 γ (1 − β ) 2 γ (1 − β ) 3 3 u + u − u v + o (ρ ) − u v + 8β 2 2β 16β 4 8β 3 2
2
where
−
aTij
= aij (βT )(i, j = 1, 2), m11 = 3β , m12 = 1, m21 = γ β , m22 = 0, ( ) = T
aT12 −(k2T )∗ d12 aT11 −(k2T )∗ d11
aT22 −(k2T )∗ d22 aT21 −(k2T )∗ d21
,g = − 1 L+ =0 T
=− f LT = 0, 1
aT11 −(k2T )∗ d11 aT21 −(k2T )∗ d21
T
=−
g
+
respectively, where LT is the adjoint operator of LT .
aT12 −k2T d12 aT22 −(k2T )∗ d22
k2T ∗
k2T
(βT ) =
satisfy the following equalities
aT11 aT22 −aT12 aT21 d11 d22 −d12 d21
, and f =
42
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
Step 2: Close to onset β = βT , we compute the solutions of model (1.5) that can be expanded as in the following asymptotic expansion:
u1
U=ε
v1
+ ε2
u2
v2
+ ε3
u3
v3
+ o(ε 3 ),
where
u1
=
v1
3 f 1
Wj exp(ikj · r) + c .c . ,
j =1
and
u2
U0 V0
=
v2
+
3 U j
j =1
+
U12 V12
Vj
exp(ikj · r) +
3 U jj
Vjj
j=1
exp(i(k1 − k2 ) · r) +
U23 V23
exp(i2kj · r)
exp(i(k2 − k3 ) · r) +
U31 V31
exp(i(k3 − k1 ) · r) + c .c .
where
U0 V0
u00
(|W1 |2 + |W2 |2 + |W3 |2 ), Uj = fVj , Ujj u11 Ujl u 2 = Wj , = ∗ Wj W l , Vjj v11 Vjl v∗ 2 γ γ (1−β 2 ) 2 1−9β aT22 β1 f − 4β 2 f 2 − aT12 f − f 2 β 4β T T T T a a − a a u00 11 22 12 21 = , v00 aT γ (1−β 2 ) f 2 − γ f − aT 1 f − 1−9β 2 f 2 11 21 β β 4β 2 4β 2 =
v00
aT11 aT22 − aT12 aT21
u11
v11
(aT22 − 4d22 (k2T )∗ )
1 f 2β
−
1−9β 2 2 f 8β 2
− (aT12 − 4d12 (k2T )∗ )
γ (1−β 2 ) 2 f 8β 2
−
γ
2β
f
(aT − 4d (k2 )∗ )(aT − 4d (k2 )∗ ) − (aT − 4d (k2 )∗ )(aT − 4d (k2 )∗ ) 11 T 22 T 12 T 21 T 11 22 12 21 = , T γ γ (1−β 2 ) 2 1−9β 2 2 1 T 2 ∗ 2 ∗ f − 2β f − (a12 − 4d12 (kT ) ) 2β f − 8β 2 f (a11 − 4d11 (kT ) ) 8β 2
(aT11 − 4d11 (k2T )∗ )(aT22 − 4d22 (k2T )∗ ) − (aT12 − 4d12 (k2T )∗ )(aT21 − 4d21 (k2T )∗ ) 2) 2 γ 2 (aT22 − 3d22 (k2T )∗ ) β1 f − 1−4β92β f 2 − (aT12 − 3d12 (k2T )∗ ) γ (14−β f − f 2 β β (aT − 3d (k2 )∗ )(aT − 3d (k2 )∗ ) − (aT − 3d (k2 )∗ )(aT − 3d (k2 )∗ ) 11 T 22 T 12 T 21 T u∗ 11 22 12 21 = . v∗ T γ (1−β 2 ) 2 γ 1−9β 2 2 1 2 ∗ T 2 ∗ f − β f − (a12 − 3d12 (kT ) ) β f − 4β 2 f (a11 − 3d11 (kT ) ) 4β 2
(aT11 − 3d11 (k2T )∗ )(aT22 − 3d22 (k2T )∗ ) − (aT12 − 3d12 (k2T )∗ )(aT21 − 3d21 (k2T )∗ ) Step 3: By the above expressions, together with Aj = ε Wj + ε 2 Vj + o(ε 3 ), we can get the amplitude equation corresponding to A1 , A2 , A3 as follows:
∂A 1 τ0 = µA1 + hA2 A3 − [g1 |A1 |2 + g2 (|A2 |2 + A3 |2 )]A1 , ∂t ∂ A2 τ0 = µA2 + hA1 A3 − [g1 |A2 |2 + g2 (|A1 |2 + A3 |2 )]A2 , ∂ t ∂ A3 τ = µA3 + hA1 A2 − [g1 |A3 |2 + g2 (|A1 |2 + A2 |2 )]A3 , 0 ∂t
(2.8)
where
τ0 = g1 =
f +g , βT [fm11 + m12 + g (fm21 + m22 )] G1
βT [fm11 + m12 + g (fm21 + m22 )]
,
µ= g2 =
βT − β , βT
h=
H
βT [fm11 + m12 + g (fm21 + m22 )]
G2
βT [fm11 + m12 + g (fm21 + m22 )]
,
,
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
H =
1 − 9β 2 4β
2
f2 −
1 − 9β 2
1
β
f +g 1
−
γ (1 − β 2 ) 2 γ f + f 4β 2 β
43
,
1
3
(u00 + u11 ) − f (v00 + v11 ) + f2 4β 2β 2β 8β 3 3(1 − β 2 ) 3 γ 3γ 2 γ (1 − β 2 ) γ 3γ (1 − β 2 ) 3 − ( u + u ) + f (v + v ) − f + g − f + f + f , 00 11 00 11 16β 4 4β 2 2β 2β 8β 3 16β 4 1 3(1 − β 2 ) 3 1 3 2 1 − 9β 2 f − (u00 + u∗ ) − f (v00 + v∗ ) + f − f −G2 = 2 3 4β 2β 2β 4β 8β 4 γ (1 − β 2 ) γ γ 3γ 2 3γ (1 − β 2 ) 3 +g − f + ( u + u ) + f (v + v ) − f + f . 00 ∗ 00 ∗ 4β 2 2β 2β 4β 3 8β 4
−G1 =
2
f −
2.2.1. Amplitude instability Each amplitude in Eq. (2.8) can be decomposed to mode ρj = |Aj | and a corresponding phase angle ϕj . Then, substituting Aj = ρj exp(iϕj ) into Eq. (2.8) and separating the real and imaginary parts, we can get four differential equations of the real variables as follows:
∂ϕ ρ 2 ρ 2 + ρ12 ρ32 + ρ22 ρ32 = − h 1 2 sin ϕ, τ 0 ∂t ρ1 ρ2 ρ3 ∂ρ τ0 1 = µρ1 + hρ2 ρ3 cos ϕ − g1 ρ13 − g2 (ρ22 + ρ32 )ρ1 , ∂t ∂ρ 2 = µρ2 + hρ1 ρ3 cos ϕ − g1 ρ23 − g2 (ρ12 + ρ31 )ρ2 , τ0 ∂t ∂ρ3 = µρ3 + hρ1 ρ2 cos ϕ − g1 ρ33 − g2 (ρ12 + ρ22 )ρ3 , τ0 ∂t where ϕ = ϕ1 + ϕ2 + ϕ3 .
(2.9)
The dynamical system (2.9) possesses five kinds of solutions in [33,34]. (1) The stationary state (O), given by
ρ1 = ρ2 = ρ3 = 0, is stable for µ < µ2 = 0 and unstable for µ > µ2 . (2) Stripe patterns (S ), given by
ρ1 =
µ g1
̸= 0,
ρ2 = ρ3 = 0, h2 g
is stable for µ > µ3 = (g −g1 )2 and unstable for µ < µ3 . 2 1 (3) Hexagon patterns (H0 , Hπ ) are given by
ρ1 = ρ2 = ρ3 =
| h| ±
h2 + 4(g1 + 2g2 )µ
2(g1 + 2g2 )
, √
2
with ϕ = 0 or π , and exist when µ > µ1 = 4(g−+h2g ) . The solution ρ + = 1 2 √ |h|−
h2 +4(g1 +2g2 )µ 2(g1 +2g2 )
and ρ − = (4) The mixed states are given by 2g1 +g2 h2 , (g2 −g1 )2
| h| ρ1 = , g2 − g1
ρ2 = ρ3 =
|h|+
h2 +4(g1 +2g2 )µ 2(g1 +2g2 )
is stable only for µ < µ4 =
is always unstable.
µ − g1 ρ12 , g1 + g2
with g2 > g1 , µ > g1 ρ12 and are always unstable. With fixed parameters d11 = d12 = 1, d21 = 11, d22 = 15, and γ = 2, we can find that h = 0.7958, g1 = 8.8043, g2 = 3.7966, µ1 = −0.0097, µ2 = 0, µ3 = 0.2223, µ4 = 0.5405. Summarizing the above analysis, we can show our results in Fig. 7. The system has a bistable region for µ ∈ (µ1 , µ2 ). In other words, when the control parameter µ lies in this region, the H0 patterns and the stationary state are all stable. The Hπ patterns are always unstable when µ > µ2 . When µ lies in region (µ2 , µ3 ), the stripe patterns are unstable, and the H0 patterns are stable. In region (µ3 , µ4 ), the system exists in another bistable state (i.e., the bistable state between the hexagon patterns and the stripe patterns). When the parameter µ is more than µ4 , stripe patterns emerge in the two-dimensional space.
44
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
Fig. 7. Bifurcation diagram of model (1.5) with d11 = d12 = 1, d21 = 11, d22 = 15, and γ = 2. H0 : hexagonal patterns with ϕ = 0; Hπ : hexagonal patterns with ϕ = π ; S: stripe patterns. Solid lines: stable states; dashed lines: unstable states. µ1 = −0.0097, µ2 = 0, µ3 = 0.2223, µ4 = 0.5405.
3. Pattern formation of model (1.5) In this section, we perform numerical simulations of model (1.5) in two dimensional space and the simulated results are shown here. To solve differential equations by computer, one has to discretize the space and time of the problem. The continuous problem defined by the reaction–diffusion system in two-dimensional space is solved in a discrete domain with M × N lattice sites. The spacing between the lattice points is defined by the lattice constant ∆h. In the discrete system the Laplacian describing diffusion is calculated using finite differences, i.e., the derivatives are approximated by differences over ∆h. The time evolution is also discrete, i.e., the time goes in steps of ∆t. The time evolution can be solved by using the Euler method. In the present paper, we set ∆h = 0.5, different ∆t values: 0.001 or 0.01, and M = N = 200. All of our numerical simulations employ the Neumann boundary conditions. And the initial data is all taken as a uniformly distributed random perturbation around the equilibrium (u∗ , v ∗ ). Through the Introduction section, we know that β is proportional to the death rate s, which plays a vital role in the spatial predator–prey model. So, from the motion point of view, we need to consider the effect of β on pattern dynamics of model (1.5), and also consider the effect of d21 on pattern dynamics of model (1.5). In what follows, we preform the numerical simulations from the two aspects of the varied β and varied cross-diffusion d21 , respectively. 3.1. Pattern formation on the effect of the varied β We first keep d11 = d12 = 1, d21 = 11, d22 = 15, and γ = 2. Then we can obtain g1 , g2 , h, µ1 , µ2 , µ3 , µ4 . In the numerical simulations, different types of patterns are observed and the distributions of predator and prey are always of the similar type. As a result, we only show our results of pattern formation about one species distribution (in this paper, we show the distribution of prey, for instance). For the different values of β located in the ‘‘Turing space’’ (the domain D11 in Fig. 1), we show and obtain three categories of Turing patterns for the distribution of prey u of model (1.5). In every pattern, the red (blue) represents the high (low) density of prey u. In Fig. 8, for prey u, Hπ hexagon patterns prevail over the whole domain eventually. Parameter values (β = 0.712) set in this figure satisfy µ = 0.07424262 ∈ (µ2 , µ3 ). According to the analysis above, there should be only Hπ hexagon patterns under this circumstance. In other words, the numerical simulation is corresponding to the theoretical analysis. We should also pay attention to the situation that µ is very close to the critical value µ2 . In this case, the Hπ hexagon patterns come into being very slowly. This is the universal phenomenon of critical slowing down. On the other hand, Fig. 8(f) consists of blue (minimum density of u) hexagons on a red (maximum density of u) background, i.e., the prey is isolated with low population density. We call this spot pattern. It exhibits a competition between stripes and spots. The pattern takes a long time to settle down, starting with a homogeneous state (u∗ , v ∗ ) = (0.5069, 0.3511) (cf. Fig. 8(a)), and the random perturbations lead to the formation of stripes and spots (cf. Fig. 8(b)–(e)), and ending with spots only (cf. Fig. 8(f)). In Fig. 9, parameter values (β = 0.619) set in this figure satisfy µ = 0.19516318 ∈ (µ2 , µ3 ). A few stripes emerge, and the remaining part of the spots pattern remains time independent, i.e., spots–stripes pattern. It exhibits a competition between stripes and spots. Starting with a homogeneous state (u∗ , v ∗ ) = (0.3832, 0.3818) (cf. Fig. 9(a)), the random perturbations lead to the formation of spots–stripes (cf. Fig. 9(b)–(e)), ending with the time-independent spots–stripes pattern (cf. Fig. 9(f)). In this case, the numerical simulations cannot correspond to the theoretical analysis. This phenomenon cannot be explained by the amplitude equations, because the control parameter β is far away from the critical value β = 0.7691. In Fig. 10, Parameter values (β = 0.5778) set in this figure satisfy µ = 0.2458718 ∈ (µ3 , µ4 ). Model dynamics exhibits a transition from spots–stripes growth to stripes replication, i.e., spots decay and the stripes pattern emerges. Starting with a homogeneous state (u∗ , v ∗ ) = (0.3339, 0.3849) (cf. Fig. 10(a)), the random perturbations lead to the formation
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
45
Fig. 8. Snapshots of contour pictures of the time evolution of the prey at different instants with d11 = d12 = 1, d21 = 11, d22 = 15, γ = 2, β = 0.712 and ∆t = 0.001. Moments: (a) t = 0; (b) t = 5100; (c) t = 7500; (d) t = 10 500; (e) t = 15 000; (f) t = 25 000.
Fig. 9. Snapshots of contour pictures of the time evolution of the prey at different instants with d11 = d12 = 1, d21 = 11, d22 = 15, γ = 2, β = 0.619 and ∆t = 0.001. Moments: (a) t = 0; (b) t = 2500; (c) t = 4000; (d) t = 10 000; (e) t = 15 000; (f) t = 25 000.
46
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
Fig. 10. Snapshots of contour pictures of the time evolution of the prey at different instants with d11 = d12 = 1, d21 = 11, d22 = 15, γ = 2, β = 0.5778 and ∆t = 0.001. Moments: (a) t = 0; (b) t = 2300; (c) t = 5000; (d) t = 9000; (e) t = 15 000; (f) t = 25 000.
Fig. 11. Snapshots of contour pictures of the time evolution of the prey at final instants. (a): β = 0.5713, (b): β = 0.5210, (c): β = 0.4712 with d11 = d12 = 1, d21 = 11, d22 = 15, γ = 2 and ∆t = 0.001. Moments: (a) t = 25 000; (b) t = 25 000; (c) t = 9900.
of spots–stripes (cf. Fig. 10(b)–(e)), and the later random perturbations make these spots decay, ending with the timeindependent stripes pattern (cf. Fig. 10(f)). In this case, the numerical simulations cannot also correspond to the theoretical analysis. This phenomenon cannot also be explained by the amplitude equations, because the control parameter β is far away from the critical value β = 0.7691. Finally, fixing d11 = d12 = 1, d21 = 11, d22 = 15, we choose some points in the Hopf and Turing instability region D13 in Fig. 1 to simulate, that is choosing four types of points: (I): β = 0.5713, 0.5210, 0.4712, γ = 2; (II): β = 0.528, 0.485, 0.380, γ = 7; (III): two points next to Turing–Hopf point: β = 0.572, γ = 5.75 and β = 0.535, γ = 6.25; (IV): one point far from Turing–Hopf point and between type (I) and type (II) points: β = 0.435, γ = 4.68. And we obtain some interesting simulation results, (a)–(c) in Fig. 11, and (a) and (b) in Fig. 13 are called as ‘‘labyrinth’’ patterns and the stripe-like patterns prevail over the whole domain (cf. Fig. 12(a)), and that the stripe-spotted patterns coexist in the space (cf. Fig. 12(b) and Fig. 13(c)), spotted patterns in the space (cf. Fig. 12(c)) which are different from the above results. Of course, we have to mention the fact that, due to the Hopf and Turing instability, the resulting pattern could exhibit oscillations and they could be so strong that our numerical scheme would be not convergent.
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
47
Fig. 12. Snapshots of contour pictures of the time evolution of the prey at final instants. (a): β = 0.528, (b): β = 0.485, (c): β = 0.380 with d11 = d12 = 1, d21 = 11, d22 = 15, γ = 7 and ∆t = 0.001. Moments: (a) t = 8800, (b) t = 10 000; (c) t = 9000.
Fig. 13. Snapshots of contour pictures of the time evolution of the prey at final instants. (a): β = 0.572, γ = 5.75, (b): β = 0.535, γ = 6.25, (c): β = 0.435, γ = 4.68 with d11 = d12 = 1, d21 = 11, d22 = 15 and ∆t = 0.001. Moments: (a) t = 13 000, (b) t = 8300; (c) t = 10 500.
3.2. Pattern formation on the effect of varied d21 Next, we continue to keep d11 = d12 = 1, d22 = 5, γ = 0.4, β = 0.678, and let d21 vary. In the 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 still restrict our analysis of pattern formation to the distribution of prey. In Fig. 14, we show the snapshots of prey spatial pattern at different values of d21 with d11 = d12 = 1, d22 = 15, γ = 0.4, β = 0.678 and 25,000 iterations. With small random perturbation of the stationary solution (u∗ , v ∗ ) = (0.4597, 0.3663) of the spatially homogeneous systems with d21 = 4.952, one can see that the stripe-like patterns prevail over the whole domain (cf. Fig. 14(e)), and that the stripe-like and spotted patterns coexist in the space (cf. Fig. 14(d)) when d21 is being decreased to 4.508. However, as d21 is being decreased to 4.000, 3.358 and 2.850, we find that the spotted patterns prevail over the whole domain, the color varies from crimson (blue), red (blue) to green (blue) (cf. Fig. 14(c)–(a)). If we continue to decrease to d21 = 2.678, the system (1.5) cannot induce spatial pattern. Therefore, we find that when the cross diffusion d21 decreases, the modes of the patterns convert from the stripes to the spots. Finally, we choose some points in the Hopf and Turing instability regions D23 in Fig. 4 to simulate, that is choosing d21 = 3.0, 3.8, 4.5, β = 0.2, γ = 0.4, d11 = d12 = 1, d22 = 5. And we obtain some interesting simulation results, (a): large yellow spots, (b): large dark blue dots, (c): small dark blue dots, in Fig. 15, which are different from the above results. Of course, we have to mention the fact that, owing to the interaction of Hopf and Turing instabilities, the resulting pattern could exhibit oscillations and they could be so strong that our numerical scheme would be not convergent and there exist oscillatory patterns. Through the numerical simulations of u above, we can know that the distributions of prey and predators are always of the same type. This means that there is a tendency that the prey stays away from the predators and the predators will get closer to the prey, which reflects the effect of the cross-diffusion. What is more, the local density of prey is higher than that of predators. This is the precondition of the coexistence of prey and predators. 4. Discussion and conclusion In this paper, we deduce the amplitude equations and present the Turing pattern selection of a predator–prey model with cross diffusion. By using mathematical analysis and numerical simulations, we found that this model can give rise to all three categories: spotted, spots–stripes mixtures, and stripe-like patterns and other patterns: ‘‘labyrinth’’ patterns and
48
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
Fig. 14. Snapshots of contour pictures of the evolution of the prey at different values of d21 with d11 = d12 = 1, d22 = 5, γ = 0.4, β = 0.678. (a): d21 = 2.850; (b): d21 = 3.358; (c): d21 = 4.000; (d): d21 = 4.508; (e): d21 = 4.952 and ∆t = 0.01. Moments: (a) t = 60 000; (b) t = 7500; (c) t = 8000; (d) t = 3500; (e) t = 2500.
Fig. 15. Snapshots of contour pictures of the time evolution of the prey at final instants. (a): d21 = 3.0, (b): d21 = 3.8, (c): d21 = 4.5 with β = 0.2,
γ = 0.4, d11 = d12 = 1, d22 = 5 and ∆t = 0.01. (a) t = 2000; (b) t = 10 500; (c) t = 10 000.
‘‘black-eye’’ patterns. That is to say, cross diffusion can create stationary patterns. It should be noted that, as cross-diffusion coefficients d12 = d21 = 0 and the predator mortality described by the linear form, the spatial predator–prey model cannot give rise to Turing structures that has been pointed out in [4], which reflects the effect of the cross-diffusion. From biological point of view, our results have more biological meanings. In a word, our results are different from [4] and well enrich the findings in predator–prey models with herd behavior. By the above analysis, we can find that the qualitative dynamics of the model (1.5) are fundamentally different when parameter β and cross-diffusion coefficient d21 in model (1.5) slightly change. The parameter β is proportional to the (non-dimensional) death rate of the predator. From the biological point of view, our results show that the death rate of predator may play a vital role in the spatial predator–prey model. By varying the value of β , we obtain three different typical types of pattern: spot pattern (Fig. 8(f)), spot–stripe pattern (Fig. 9(f)), and stripe pattern (Fig. 10(f)). From the view of population dynamics, one can see that there exists the spot pattern replication—the prey u is the isolated zone with low density, and the remaining region is high density, which means the prey may break out in the area in Fig. 8(f). In other words, the prey in this area is safe. The biological significance of the other cases can be determined in the same way as the above case. On the other hand, we only let d21 vary, that is decreasing the value of d21 from 4.952 to 2.850, and find that the system (1.5) gives rise to all three categories: stripe pattern (Fig. 14(e)), spots–stripes mixtures (Fig. 14(d)), spot pattern (Fig. 14(c–a)). But when the value of d21 decreases to 2.678, the system (1.5) cannot induce spatial pattern. In Figs. 11, 13 and 15, we also obtain other patterns: ‘‘labyrinth’’ patterns and ‘‘black-eye’’ patterns, respectively.
X. Tang, Y. Song / Nonlinear Analysis: Real World Applications 24 (2015) 36–49
49
Of course, the methods and the results in the present paper can be applied to the pattern formation in the diffusive predator–prey model or other reaction–diffusion systems. We hope that our work could be instructive to study the effect of cross-diffusion on the pattern formation in population dynamics. Acknowledgments The authors would like to thank the referee for his/her careful reading and some comments on improving the presentation of this paper. The work was supported by the Shanghai Committee of Science and Technology (No. 13ZR1444500) and the Program for New Century Excellent Talents in University (NCET-11-0385). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42]
V. Volterra, Sui tentutive di applicazione delle mathematiche alle seienze biologiche e sociali, Ann. Radioelectr. Univ. Romandes (1901-2) 3–28. V. Volterra, Variazione e fluttuazini del numero d’individui in specie animali conviventi, Mem. R. Accad. Naz. dei Lincei 2 (1926) 31–113. P.A. Braza, Predator–prey dynamics with square root functional responses, Nonlinear Anal. RWA 13 (2012) 1837–1843. S. Yuan, C. Xu, T. Zhang, Spatial dynamics in a predator–prey model with herd behavior, Chaos 23 (2013) 0331023. A.M. Turing, The chemical basis of morphogenesis, Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. 237 (1952) 37–72. J.D. Murray, Mathematical Biology, Springer-Verlag, Berlin, 1993. S. Aly, I. Kim, D. Sheen, Turing instability for a ratio-dependent predator–prey model with diffusion, Appl. Math. Comput. 217 (2011) 7265–7281. L.N. Guin, P.K. Mandal, Spatial pattern in a diffusive predator–prey model with sigmoid ratio-dependent functional response, Int. J. Biomath. 7 (2014) 1–26. 1450047. Y. Song, X. Zou, Spatiotemporal dynamics in a diffusive ratio-dependent predator–prey model near a Hopf–Turing bifurcation point, Comput. Math. Appl. 67 (2014) 1978–1997. Y. Song, X. Zou, Bifurcation analysis of a diffusive ratio-dependent predator–prey model, Nonlinear Dynam. 78 (2014) 49–70. Z. Xu, Y. Song, Bifurcation analysis of a diffusive predator–prey system with a herd behavior and quadratic mortality, Math. Meth. Appl. Sci. (2014). http://dx.doi.org/10.1002/mma.3275. J. Shukla, S. Verma, Effects of convective and dispersive interactions on the stability of 2 species, Bull. Math. Biol. 43 (1981) 593–610. E.H. Kerner, A statistical mechanics of interacting biological species, Bull. Math. Biol. 19 (1957) 121–146. N. Shigesada, K. Kawasaki, E. Teramoto, Spatial segregation of interacting species, J. Theoret. Biol. 79 (1979) 83–99. W. Wang, Q. Liu, Z. Jin, Spatiotemporal complexity of a ratio-dependent predator–prey system, Phys. Rev. E 75 (2007) 051913. K. Kuto, Stability of steady-state solutions to a prey–predator system with cross-diffusion, J. Differential Equations 197 (2004) 293–314. K. Kuto, Y. Yamada, Multiple coexistence states for a prey–predator system with cross-diffusion, J. Differential Equations 197 (2004) 315–348. V.K. Vanag, I.R. Epstein, Cross-diffusion and pattern formation in reaction–diffusion systems, Phys. Chem. Chem. Phys. 11 (2009) 897–912. X. Zhang, G. Sun, Z. Jin, Spatial dynamics in a predator–prey model with Beddington–DeAngelis functional response, Phys. Rev. E 85 (2012) 021924. L. Xue, Pattern formation in a predator–prey model with spatial effect, Physica A 391 (2012) 5987–5996. G. Gambino, M.C. Lombardo, M. Sammartino, Pattern formation driven by cross-diffusion in a 2D domain, Nonlinear Anal. RWA 14 (2013) 1755–1779. L.N. Guin, Existence of spatial patterns in a predator–prey model with self- and cross-diffusion, Appl. Math. Comput. 226 (2014) 320–335. G. Galiano, M.L. Garzon, A. Jüngel, Semi-discretization in time and numerical convergence of solutions of a nonlinear cross-diffusion population model, Numer. Math. 93 (2003) 655–673. L.N. Guin, P.K. Mandal, Spatiotemporal dynamics of reaction–diffusion models of interacting populations, Appl. Math. Model. 38 (2014) 4417–4427. C. Tian, Z. Lin, M. Pedersen, Instability induced by cross-diffusion in reaction–diffusion systems, Nonlinear Anal. RWA 11 (2010) 1036–1045. C. Tian, Z. Ling, Z. Lin, Turing pattern formation in a predator–prey-mutualist system, Nonlinear Anal. RWA 12 (2011) 3224–3237. C. Tan, L. Zhang, Z. Lin, Pattern formation for a model of plankton allelopathy with cross-diffusion, J. Franklin Inst. 348 (2011) 1947–1964. R.R. Baier, C. Tian, Mathematical analysis and numerical simulation of pattern formation under cross-diffusion, Nonlinear Anal. RWA 14 (2013) 601–612. X. Wang, Y. Cai, Cross-diffusion-driven instability in a reaction–diffusion harrison predator–prey model, Abst. Appl. Anal. 2013 (2013) 306467. W. Wang, Z. Guo, R.K. Upadhyay, Y. Lin, Pattern formation in a cross-diffusive Holling–Tanner model, Discrete Dyn. Nat. Soc. 2012 (2012) 828219. M.C. Cross, P.C. Hohenberg, Pattern formation outside of equilibrium, Rev. Modern Phys. 65 (1993) 851–1112. G. Gambino, M.C. Lombardo, M. Sammartino, Turing instability and pattern formation for the Lengyel–Epstein system with nonlinear diffusion, Acta Appl. Math. 132 (2014) 283–294. Q. Ouyang, Nonlinear Science and the Pattern Dynamics Introduction, Peking University Press, Beijing, 2010. Y. Kuramoto, T. Tsuzuki, On the formation of dissipative structures in reaction–diffusion systems, Progr. Theoret. Phys. 54 (1975) 687–699. G.H. Gunaratne, Q. Ouyang, H.L. Swinney, Pattern formation in the presence of symmetries, Phys. Rev. E 50 (1994) 2802–2820. B. Peña, C. Pérez-García, Stability of Turing patterns in the Brusselator model, Phys. Rev. E 64 (2001) 056213. W. Wang, Y. Lin, L. Zhang, F. Rao, Y. Tan, Complex patterns in a predator–prey model with self and cross-diffusion, Commun. Nonlinear Sci. Numer. Simul. 16 (2011) 2006–2015. G. Sun, Z. Jin, Q. Liu, L. Li, Pattern formation induced by cross-diffusion in a predator–prey system, Chin. Phys. B 17 (2008) 3936–3941. W. Wang, Y. Lin, F. Rao, L. Zhang, Y. Tan, Pattern selection in a ratio-dependent predator–prey model, J. Stat. Mech. 11 (2010) P11036. F. Rossi, V.K. Vanag, I.R. Epstein, Pentanary cross-diffusion in water-in-oil microemulsions loaded with two components of the Belousov–Zhabotinsky reaction, Chem. Eur. J. 17 (2011) 2138–2145. G. Gambino, M.C. Lombardo, M. Sammartino, V. Sciacca, Turing pattern formation in the Brusselator system with nonlinear diffusion, Phys. Rev. E 88 (4) (2013) 042925. Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence, Springer-Verlag, Berlin, 1984.