- Email: [email protected]

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Pattern dynamics in a vegetation model with time delay✩ Bo-li Xie Department of Mathematics, North University of China, Taiyuan, Shan’xi, 030051, People’s Republic of China

highlights • A vegetation model with both cross diffusion and time delay is considered. • There exists a threshold of time delay in predator–prey interactions. • The spatial patterns via numerical simulations are illustrated.

article

abstract

info

Article history: Received 3 April 2015 Available online 9 October 2015

A vegetation model with both cross diffusion and time delay is considered. Based upon a stability analysis, we demonstrate that the delay affects the stability and spatial patterns under some conditions. In addition, through numerical simulations, we obtain different spatial patterns, including spot patterns and stripe-like patterns. These results may help us better understand the dynamics of a vegetation model with time delay. © 2015 Elsevier B.V. All rights reserved.

Keywords: Vegetation Pattern formation Diffusion Time delay

1. Introduction Vegetation plays a critical role in life. On the one hand, vegetation absorbs carbon dioxide and releases oxygen, which maintains the balance of the ecosystem. On the other hand, vegetation evaporates soil moisture into the atmosphere, which promotes the global water circulation cycle. Klausmeier proposed a very simple partial differential equation model with equations for surface water R(X , Y , T ) and for vegetation U (X , Y , T ) [1]: dR dT dU dT

= d − eR − fU 2 R, (1.1)

= aU 2 R − bU ,

where a is the plant growth rate, b is the plant loss rate, d is the amount of rainfall, e is the water evaporation rate, and f is the uptake rate of the plant. According to Sherratt [2,3], through a non-dimensional transformation a

γ = R√ , ef

u=U

f e

,

t = Te,

✩ This work is supported by the National Sciences Foundation of China (10471040) and the National Sciences Foundation of Shanxi Province (2009011005-1). E-mail address: [email protected]

http://dx.doi.org/10.1016/j.physa.2015.09.075 0378-4371/© 2015 Elsevier B.V. All rights reserved.

B.-l. Xie / Physica A 443 (2016) 460–466

461

we arrive at the following equations: dγ dt du dt

= S − γ − u2 γ , (1.2)

= u2 γ − Bu,

where S = ade3/2

f,

B = b/e.

In most practically relevant cases, the state of the ecosystem can exhibit temporal long-range correlations and should first be affected by its immediate past, with additional correction arising from time delays. In practice, time delays always exist and play a significant role in the dynamics; for example, time delays can be used to introduce oscillations [4–8]. When combined with the time delay factor, the model is given by

∂γ = S − γ − [u(t − τ )]2 γ , ∂t ∂u = [u(t − τ )]2 γ − Bu(t − τ ), ∂t where τ > 0 is a constant delay due to gestation.

(1.3)

Spatial patterns are ubiquitous in nature. These patterns modify the temporal dynamics and stability properties of population densities over 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

∂γ = S − γ − [u(t − τ )]2 γ + d1 ∇ 2 γ ∂t ∂u = [u(t − τ )]2 γ − Bu(t − τ ) + d2 ∇ 2 u, ∂t

(1.4)

where ∇ 2 = ∂∂x2 or ∇ 2 = ∂∂x2 + ∂∂y2 is the typical Laplacian operator in one- or two-dimensional space. d12 and d21 are the cross-diffusion coefficients of surface water and vegetation, respectively. However, to the best of our knowledge, little work has been conducted on the dynamical behavior of both time delay and diffusion in a vegetation model. Therefore, in the present study, our objective is to investigate a vegetation model with both cross diffusion and time delay. More specifically, the primary objective of the present study is to investigate the spatial patterns. 2

2

2

2. Analysis In this section, we consider the dynamics of model (1.4) with discrete time delay. We determined that model (1.4) and model (1.2) possessed the same positive equilibrium point E ∗ (γ ∗ , u∗ ), where

√

γ = ∗

u∗ =

S−

S 2 − 4B2

2 2B

√ S−

S 2 − 4B2

, .

It is easy to obtain that the condition for ensuring that γ ∗ and u∗ are positive is that S > 2B. We always assume that (γ ∗ , u∗ ) is linearly stable, thus, the eigenvalues of the Jacobian

∂f ∂γ A= ∂g ∂γ

∂f ∂u ∂g ∂u

,

a11 a21

(γ ∗ ,u∗ )

a12 a22

=

2B

−

√ S−

(S −

S 2 − 4B2 4B2

√

S 2 − 4B2 )2

−2B B

at (γ ∗ , u∗ ) must have negative real parts, which is equivalent to tr (A) = a11 + a22 < 0,

(2.1)

det (A) = a11 a22 − a12 a21 > 0.

(2.2)

and

462

B.-l. Xie / Physica A 443 (2016) 460–466

Next, we consider small spatiotemporal perturbations h and p on a homogeneous steady state E ∗ (γ ∗ , u∗ ). Let h = γ − γ ∗ and p = u − u∗ ; then, we derive that

∂h = a11 h + a12 p(t − τ ) + d1 ∇ 2 h, ∂t ∂p = a21 h + a22 p(t − τ ) + d2 ∇ 2 p. ∂t

(2.3)

The characteristic equation for the linear system (2.3) is given by

λ2 − (a11 − d1 k2 − d2 k2 )λ − (a11 − d1 k2 )d2 k2 +[a22 (a11 − d1 k2 ) − a12 a21 ]e−λτ − a22 λe−λτ = 0.

(2.4)

Spatial patterns form if Eq. (2.4) has root λ = iω, which are called delay-driven spatial patterns. Moreover, the critical value of the delay τ is called the Turing bifurcation. If iω is a root of Eq. (2.4), then we have

ω2 + (a11 − d1 k2 )d2 k2 = [a22 (a11 − d1 k2 ) − a12 a21 ] cos ωτ − ωa22 sin ωτ , −[(a11 − d1 k2 − d2 k2 )]ω = [a22 (a11 − d1 k2 ) − a12 a21 ] sin ωτ + ωa22 cos ωτ

(2.5)

which leads to

ω4 + αk ω2 + βk = 0,

(2.6)

where

αk = (a11 − d1 k2 − d2 k2 )2 − a222 , βk = [(a11 − d1 k2 )d2 k2 ]2 − [a22 (a11 − d1 k2 ) − a12 a21 ]2 . Using mathematical calculations, the Hopf bifurcation occurs when Im(λ) ̸= 0,

Re(λ) = 0

at k = 0.

(2.7)

From Eq. (2.4), we obtain

ω4 + (a211 − a222 )ω2 − (a11 a22 − a12 a21 )2 = 0.

(2.8)

Then, Eq. (2.8) has the solution

−(a211 − a222 ) +

ωc2 =

(a211 − a222 )2 + 4(a211 − a222 )(a11 a22 − a12 a21 )2 2

.

(2.9)

From Eq. (2.5), we can obtain

τc =

1

ωc

arccos

−a12 a21 ω2 . (a22 a11 − a12 a21 )2 + ω2 a222

(2.10)

dRe(λ)

Now, we investigate the sign of dτ |τ =τc . When k = 0, from (2.4), we can see that

λ2 − a11 λ + (a22 a11 − a12 a21 )e−λτ − a22 λe−λτ = 0.

(2.11)

Both sides of the equation derivation of τ , we can get 2λ

dλ dτ

dλ

dλ

+ (a22 a11 − a12 a21 ) −τ e − λe dτ dλ dλ −a22 e−λτ − a22 λ −τ e−λτ − λe−λτ = 0. dτ dτ

− a11

−λτ

−λτ

dτ

(2.12)

From Eq. (2.12), we obtain

dλ dτ

−1 =

2λ − a11 a22 τ − − . λe−λτ [(a22 a11 − a12 a21 ) − λa22 ] λ[(a22 a11 − a12 a21 ) − λa22 ] λ

(2.13)

B.-l. Xie / Physica A 443 (2016) 460–466

463

Fig. 1. (A) Time series of the induced defense and herbivore; (B) Phase diagram of induced defense and herbivore. Parameter values are used as: S = 10, B = 2, τ = 0.45 < τc .

By substituting λ = iω into model (2.13), we derive that:

dλ

2iω − a11

−1 =

dτ

iω(cos ωτ − i sin ωτ )[(a22 a11 − a12 a21 ) − iωa22 ]

−

a22 iω[(a22 a11 − a12 a21 ) − iωa22 ]

−

τ iω

.

(2.14)

Thus we obtain

Re

dλ

−1

dτ

|τ =τc =

1

ω[(a22 a11 − a12 a21 )2 + ω2 a222 ]

(−[a11 (a22 a11 − a12 a21 ) + 2ω2 a22 ] sin ωτ

+ [2ω(a22 a11 − a12 a21 ) − ωa11 a22 ] cos ωτ − ωa222 ) =

1

[(a22 a11 − a12 a21 )2 + ω2 a222 ]2

((a22 a11 − a12 a21 )2 (a211 + 2ω2 )

+ ω2 a211 a222 + 2ω4 a222 − a222 (a22 a11 − a12 a21 )2 − a422 ω2 ) =

(a22 a11 − a12 a21 )2 (a211 + 2ω2 ) + ω4 a222 > 0. [(a22 a11 − a12 a21 )2 + ω2 a222 ]2

(2.15)

So we get the following conclusions: If the delay τ is satisfied τ ̸= 0, then the system (1.4) exist a Hopf bifurcation critical τ = τc . When τ > τc , the positive equilibrium E ∗ is unstable. We take the following values: S = 10, B = 2. Through calculations, we obtain the critical value τc = 0.4692, then E ∗ = (0.41742, 4.791337). The initial value is (0.4, 4). We adopt τ = 0.45 < τc . From Fig. 1, we can see that the positive equilibrium E ∗ is stable. We adopt τ = 0.5 > τc . From Fig. 2, we can see that the positive equilibrium E ∗ is unstable. Using mathematical calculations, the Turing bifurcation occurs when Im(λ) = 0,

Re(λ) = 0

at k = kT ̸= 0.

(2.16)

By substituting λ = 0 into the equation (2.4). By setting βkmin = d1 d2 k4min − (a11 d2 + a22 d1 )k2min + a11 a22 − a12 a21 = 0 and √ the wave number k2min = (a11 a22 − a11 a21 )/d1 d2 we can obtain the Turing bifurcation parameter ST , which is equal to

ST =

B d2 a12 d1 (8d22 + 3d21 a212 + 7d2 a12 d1 − 2d1 a12 P − 4d2 P ) d2 (Bd1 + d2 )

where P =

2d21 B2 − 2d2 Bd1 .

,

(2.17)

464

B.-l. Xie / Physica A 443 (2016) 460–466

Fig. 2. (A) Time series of the induced defense and herbivore; (B) Phase diagram of induced defense and herbivore. Parameter values are used as: S = 10, B = 2, τ = 0.5 > τc .

Fig. 3. Bifurcation diagram for the system (1.4). We set the parameter values as d1 = 1000, d2 = 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

To clearly show the bifurcation condition, we let d1 = 1000, d2 = 1. The red line of Γ1 is the Hopf bifurcation critical line that corresponds to τ = 0, and the green line of Γ4 is the Turing bifurcation critical line. T is located Turing space between Γ1 and Γ4 . When τ = 0.02, the Hopf bifurcation line Γ2 (the blue line) is upon Γ1 . When τ = 0.1, the Hopf bifurcation line Γ3 (the yellow line) is also upon Γ1 (see Fig. 3). 3. Pattern structures In the following, we will expand the space of the two-dimensional model (1.4) to perform a series of numerical simulations using the aforementioned flow during simulations under zero boundary conditions and discrete areas of a M × N lattice. For model (1.4), space and time were approximated using the finite difference method and Euler’s method, taking the time step as △t = 0.0001, the space step as △h = 1, and M = N = 200. The results indicated that △h and △t are reduced and do not lead to considerable changes in the results. To conduct the numerical simulations, we take the following values: S = 10,

B = 2,

d1 = 1000,

d2 = 1 .

B.-l. Xie / Physica A 443 (2016) 460–466

465

Fig. 4. (Color online) Snapshots of the time evolution of the surface water at different instants with S = 10, B = 2, d1 = 1000, d2 = 1, τ = 0.1, which are in the Turing space. (A): 0 iterations; (B): 10 000 iterations; (C): 50 000 iterations; (D): 500 000 iterations.

Through calculations, we obtain the critical value τc = 0.5868, and we adopt τ

(0.41742, 4.791337).

= 0.1 < τc ; then, (u∗ , v ∗ ) =

In Fig. 4, we set S = 10, B = 2, d1 = 1000, d2 = 1, τ = 0.1. As time passes, regular spotted patterns appear in space, and the dynamics of the system do not undergo any further changes. We take the following values: S = 10,

B = 4,

d1 = 1000,

d2 = 1.

Through calculations, we obtain the critical value τc = 0.11898, and we adopt τ = 0.02 < τc ; then, (u∗ , v ∗ ) = (2, 2). In Fig. 5, we set S = 10, B = 4, d1 = 1000, d2 = 1, τ = 0.02. In this case, the infected populations exhibit stripe-like patterns. 4. Discussions In this paper, a vegetation model with a time delay was considered. First, we showed that there exists a Hopf bifurcation threshold τc of time delay; when τ < τc , the positive equilibrium E ∗ of system (1.4) is stable. However, when τ > τc , the positive equilibrium E ∗ of system (1.4) is unstable and period solution will emerge. The time evolution and solution of system are plotted by simulations in the end. Second, the spatial patterns via numerical simulations are illustrated, which show that the model dynamics exhibit rich parameter-space Turing structures, including spot patterns and stripe-like patterns. The obtained results show that this system has rich dynamics, and these patterns show that this system is useful for a predation model both with predator cannibalism and a delay effect to reveal the spatial dynamics in the real model. Although further work is needed, in principle, it appears that delay and diffusion are able to generate many different types of spatiotemporal patterns. For such reasons, we hypothesize that delay and diffusion can be considered important mechanisms for the appearance of complex spatiotemporal dynamics in other models, such as predator–prey models and mutualistic models.

466

B.-l. Xie / Physica A 443 (2016) 460–466

Fig. 5. (Color online) Snapshots of the time evolution of the surface water at different instants with S = 10, B = 4, d1 = 1000, d2 = 1, τ = 0.02, which are in the Turing space. (A): 0 iterations; (B): 10 000 iterations; (C): 50 000 iterations; (D): 500 000 iterations.

References [1] C.A. Klausmeier, Regular and irregular patterns in semiarid vegetation, Science 284 (1999) 1826–1828. [2] J.A. Sherratt, An analysis of vegetative stripe formation in semi-arid landscape, J. Math. Biol. 51 (2005) 183–197. [3] J.A. Sherratt, G.J. Lord, Nonlinear dynamics and pattern bifurcations in a model for vegetation stripes in semi-arid environments, Theor. Popul. Biol. 71 (2007) 1–11. [4] N.B. Janson, A.G. Balanov, E. Schöll, Delayed feedback as a means of control of noise-induced motion, Phys. Rev. Lett. 93 (1) (2004) 010601. [5] D. Bratsun, D. Volfson, L.S. Tsimring, J. Hasty, Delay-induced stochastic oscillations in gene regulation, Proc. Natl. Acad. Sci. USA 102 (41) (2005) 14593–14598. [6] B.D. Aguda, Y. Kim, M.G. Piper-Hunter, A. Friedman, C.B. Marsh, MicroRNA regulation of a cancer network: Consequences of the feedback loops involving miR-17-92, E2F, and Myc, Proc. Natl. Acad. Sci. USA 105 (50) (2008) 19678–19683. [7] L.F. Lafuerza, R. Toral, Exact solution of a stochastic protein dynamics model with delayed degradation, Phys. Rev. E 84 (2011) 051121. [8] J. Miekisz, J. Poleszczuk, M. Bodnar, U. Foryś, Stochastic models of gene expression with delayed degradation, Bull. Math. Biol. 73 (9) (2011) 2231–2247.