Delay-driven spatial patterns in a network-organized semiarid vegetation model

Applied Mathematics and Computation 367 (2020) 124778

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Delay-driven spatial patterns in a network-organized semiarid vegetation model Canrong Tian a, Zhi Ling b, Lai Zhang b,∗ a b

Department of Basic Sciences, Yancheng Institute of Technology, Yancheng 224003, China School of Mathematical Science, Yangzhou University, Yangzhou 225002, China

a r t i c l e

i n f o

Article history: Received 30 October 2018 Revised 23 July 2019 Accepted 23 September 2019

MSC: 35B36 35B40 92D25 Keywords: Spatial pattern Time delay Hopf bifurcation Complex network

a b s t r a c t Spatiotemporal dynamics of vegetation models are traditionally investigated on a spatially continuous domain. The increasingly fragmented agricultural landscape necessitates network-organized models. In this paper, we develop a semiarid vegetation model to describe the spatiotemporal dynamics between plant and water on a network accounting for fragmented habitats which are connected by dispersal of seeds. Time delay is introduced to account for time lag in water uptake. By linear stability analysis we show that the coexistence equilibrium is asymptotically stable in the absence of time delay, but loses its stability via Hopf bifurcation when time delay is beyond a threshold. Applying the center manifold theory, we derive the explicit formulas that determine the stability and direction of the Hopf bifurcation. Numerical simulations demonstrate the emergence of spatial patterns on a network. Comparing our network-organized model to other model variants, we find that increasing landscape fragmentation is more likely to generate the variation of plant density among different habitats. © 2019 Elsevier Inc. All rights reserved.

1. Introduction A variety of vegetation patterns have been observed particularly in arid regions, including spotted vegetation, labyrinths, gap patterns, and regular bands [1–5]. Apart from the well-recognized mechanisms (i.e. habitat heterogeneity) that underly the formation of vegetation patterns [6], great attention has been paid to the positive feedback between plant growth and water availability that triggers the emergence of self-organized pattern formation [7], and self-organized vegetation patterns have been widely investigated [8–15]. The widely employed modeling approach is reaction–diffusion systems. As an application, Klausmeier [14] developed a reaction–advection–diffusion system to describe the interaction between plant and water, where the motions of plant and water are, respectively, described by a reaction–diffusion equation and an advection equation. He found a close agreement between theoretical predictions and filed observations. Under the framework of reaction–diffusion systems, the classic mechanisms behind the formation of spatial pattern is Turing bifurcation, which means that a locally stable equilibrium becomes unstable due to certain spatially varying perturbations, leading to spatially heterogeneous steady state [16,17]. Turing–Hopf bifurcation, a time and space breaking symmetry, is a well-known mechanism for generating spatio-temporal patterns [18–20], which is also known as oscillatory Turing, or wave instability [21,22].

∗

Corresponding author. E-mail address: [email protected] (L. Zhang).

https://doi.org/10.1016/j.amc.2019.124778 0 096-30 03/© 2019 Elsevier Inc. All rights reserved.

2

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

However, most of the aforementioned studies were carried out on a spatially continuous domain, that is, the entire landscape is spatially continuous. The frequent anthropogenic activities have been breaking landscape into many fine-scale habitats. These fine-scale habitats form a network connected by dispersal of seeds. Recent studies of pattern formation on network-organized systems in physics show that Turing patterns can take place in large network [23]. A further recent study shows that Turing-like waves can emerge in a one-species delayed reaction–diffusion model on a network [24], differing from the convectional knowledge that Turing patterns can only occur in two-component reaction–diffusion systems [17,18]. Motivated by these observations, we are interested in understanding whether spatial pattern can emerge in semiarid vegetation models on networks standing for fragmented habitats, which remains unexplored. To this aim, we propose the following delayed reaction–diffusion vegetation system on a network

u˙ i (t ) = a − ui (t ) − ui (t )v2i (t ) + Du

n 

Li j u j (t ),

j=1

v˙ i (t ) = ui (t − τ )v2i (t ) − bvi (t ) + Dv

n 

Li j v j (t ),

∀i = 1, 2, . . . , n,

(1.1)

j=1

ui (ζ ) ≥ 0 is continuous on − τ ≤ ζ < 0, and ui (0 ), vi (0 ) > 0. Here u and v are densities of water and plant at time t, respectively. Parameters a and b are, respectively, the rates of rainfall and plant mortality. Plant takes up water at a rate ui (t − τ )v2i (t ). The time delay τ means that some time lag is required for plant to consume water. This semiarid vegetation model is defined on an undirected network with n nodes but no self-loops. G = (gi j )n×n is an adjacency matrix, whose entry gij equals 1, if there exists an edge from node j to node i, and 0, otherwise. It is obvious that gii = 0 for i = 1, 2, . . . , n. We assume that G is strongly connected, i.e., for any pair of distinct nodes, there exists a path from one node to the other. It is known that the adjacency matrix G is irreducible [25]. The Laplacian matrix L is defined as

⎛ n ⎜

L=⎝

−

k=1

g1k

−

g21 ··· gn1

gn12

k=1

··· gn2

g2k

··· ··· ··· ···

−

g1n g2n · · ·

n k=1

⎞

⎟ ⎠. gnk

Interaction between water and plant occurs in each node and their time-dependent concentrations in the ith node are denoted by ui (t) and vi (t), respectively. The water and plant move across the network links and the process is characterized by two diffusion coefficients Du > 0 and Dv > 0, and the links represent dispersal of plant seed and water flow from one habitat to another. Time delay was first considered in a prey–predator system by Volterra [26] and found to be able to drive the emergence of oscillatory behavior. Since then, many studies about delay-driven pattern formation followed (e.g., [27–33]). In addition, time delay was also introduced in experiments to induce spiral wave [34]. Time delay in a network-organized model with discrete domain was recently introduced in [24]. It was found that travelling waves can occur following a symmetry-breaking instability of homogeneous stationary stable solution. This Turing-like patterns emerge in a single species system. A followup work on a two-component delayed reaction–diffusion system on a complex network shows that the system can evolve towards a stationary Turing pattern or to a wave pattern associated with a diffusion induced symmetry-breaking instability [35]. The paper is structured as follows. In Section 2, we perform analytical analysis to determine the local stability of the coexistence equilibrium between plant and water, derive the conditions for Hopf bifurcation to occur, and finally apply the center manifold theory to determine the direction as well as stability of the Hopf bifurcation. In Section 3, we carry out numerical analysis to demonstrate the dynamic behavior of plant on a network domain. The paper is closed with a brief discussion. 2. Delay-driven spatial patterns In this section we first carry out a linear stability analysis around the positive equilibrium (Section 2.1), and (then explore the condition and stability of the emerging Hopf bifurcation by means of the normal form theory and the center manifold theory (Section 2.2). We find that the system (1.1) can generate spatial patterns on the discrete domain when time delay goes beyond a threshold value, which is, however, impossible in the absence of time delay. The system (1.1) can possibly have three spatially homogeneous equilibrium solutions E0 , E1 and E2 , where



E0 = (a, 0 ), E1 =

a+



 √ √ a2 − 4b2 2b a − a2 − 4b2 2b , , E2 = (uˆ, vˆ ) = , . √ √ 2 2 a + a2 − 4b2 a − a2 − 4b2

The boundary equilibrium E0 always exits, which stands for a desert state where there are no plants at all. Whenever the condition a > 2b is satisfied, the two interior equilibria E1 and E2 exist, implying a coexistence between water and plant. Routine linear stability analysis around E1 in the associated non-diffusion system shows that it is a saddle point. Although the boundary equilibrium and the unstable interior equilibrium can be of ecological relevance, we are more interested in

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

3

understanding how time delay affects the stable equilibrium E2 . However, the employed techniques here can be equally applied to the other two equilibria. 2.1. Linear stability analysis In this section, we conduct linear stability analysis in the system (1.1) around E2 . To this end, we denote by (ui , vi )1≤i≤n = (uˆ, vˆ ) the homogeneous equilibrium. In order to determine its stability we translate the equilibrium (uˆ, vˆ ) to the origin by δ ui = ui − uˆ and δvi = vi − vˆ . The linearized system of (1.1) around the equilibrium E2 can be expressed by







n  δ ui δ ui,τ δ u˙ i j=1 Li j δ u j = J1 + J2 + D n , δvi δvi,τ δv˙ i j=1 Li j δv j

where δ ui,τ = δ ui (t − τ ), δvi,τ = δvi (t − τ ), and



J1 =

−1 − vˆ 2

−2b

0

b



,

J2 =

0

0

vˆ 2

0

(2.1)

,

D=

Du

0

0

Dv

 .

(2.2)

We use the same notations as in [35] to deal with the Laplacian matrix in the network system (2.1). Let 0 = 1 > > · · · > n be the eigenvalues of the Laplacian matrix L, and the eigenbasis {α : α = 1, 2, . . . , n} be each eigenvector associated to a topological eigenvalue α . For any small perturbation (δ ui (t), δ vi (t))T (T denotes the transpose) around the equilibrium, we have the following basis decomposition

2



α n δ ui (t ) c1  = eλα t αi , (i = 1, 2, . . . , n ). δvi (t ) c2α

(2.3)

α =1

Inserting (2.3) into (2.1), noticing that each mode α £ (α = 1, . . . , n),

λα

α

α

cα

cα

c1

eλα t = J1

c1

eλα t + J2

α

2

2

n

c1

cα

α =

j=1 Li j  j

α αi and using the orthogonality of the eigenvectors, we get for

eλα (t−τ ) + α D

2

α c1

c2α

eλα t .

Hence we obtain the following characteristic equation

   (λα , τ ) = 0, where (λα , τ ) = det λα I − J1 + e−λα τ J2 + α D .

(2.4)

Substituting (2.2) into (2.4), the characteristic equation of (1.1) becomes

λα 2 + A1 λα + A2 + A3 e−λα τ = 0, vˆ 2

(2.5)

− Du α − Dv α , A2 = −b(1 + vˆ 2 ) + bDu α − Dv (1 + vˆ 2 )α + Du Dv α 2 , A3 = 2bvˆ 2 .

where A1 = 1 − b + By analyzing the characteristic Eq. (2.4), we have the following results on the local stability of E2 , the proof of which is deferred to Appendix A. Lemma 2.1. Suppose that





a > max 2b, b2 , 4bDu Dv (vˆ 2 − 1 ) > (bDu − Dv (vˆ 2 + 1 ))2 ,

(2.6)

then (i) If the delay is absent, that is, τ = 0, then all characteristic roots of (2.5) have a negative real part. (ii) If the delay is present, that is, τ > 0, then the characteristic equation of (2.5) has a pair of purely imaginary roots ± iω∗ at τ = τ j , where

ω∗2 = and

τj =





 1 2A2 − A21 + A41 − 4A21 A2 + 4A23 , 2

1

ω∗



2 jπ + arcsin

 A ω∗  1 A3

,

j = 0, 1, 2, . . . .

(2.7)

(2.8)

Here A1 = 1 − b + vˆ 2 − Du α − Dv α , A2 = −b(1 + vˆ 2 ) + bDu α − Dv (1 + vˆ 2 )α + Du Dv α 2 , A3 = 2bvˆ 2 . Theorem 2.1. Suppose that (2.6) holds true, then the system (1.1) is locally unstable around the equilibrium E2 for τ ∈ [τ 0 , ∞), where τ 0 is given in (2.8). Moreover, the system (1.1) undergoes a Hopf bifurcation at τ = τ0 . There exist a constant σ > 0 and a smooth curve λ(τ ) : (τ0 − σ , τ0 + σ ) → C such that λ(τ0 ) = iω∗ and (iω∗ , τ0 ) = 0 for all τ ∈ (τ0 − σ , τ0 + σ ).

4

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

Proof. In order to guarantee E2 being locally unstable for τ ∈ [τ 0 , ∞), we only have to show there exists at least one characteristic root with positive real part for τ ∈ [τ 0 , ∞). Actually, It is sufficient to show that the following transversality condition holds

d Reλ(τ )|τ =τ0 > 0. dτ

(2.9)

Substituting a complex eigenvalue λα = μ + iω into (2.5) and separating the real and the imaginary parts, we have

μ2 − ω2 + A1 μ + A2 + A3 e−μτ cos ωτ = 0, 2μω + A1 ω − A3 e−μτ sin ωτ = 0.

(2.10)

Differentiating (2.10) with respect to τ , we get

dμ dω − (2ω + A3 τ e−μτ sin ωτ ) = A3 e−μτ (μ cos ωτ + ω sin ωτ ), dτ dτ dω + (2μ + A1 − A3 τ e−μτ cos ωτ ) = A3 e−μτ (μ sin ωτ + ω cos ωτ ). dτ

(2μ + A1 − A3 τ e−μτ cos ωτ ) (2ω + A3 τ e−μτ sin ωτ )

dμ dτ

Eliminating ddω τ , we have



 dμ (2μ + A1 − A3 τ e−μτ cos ωτ )2 + (2ω + A3 τ e−μτ sin ωτ )2 dτ  = A3 e−μτ (μ cos ωτ + ω sin ωτ )(2μ + A1 − A3 τ e−μτ cos ωτ )  + (μ sin ωτ + ω cos ωτ )(2ω + A3 τ e−μτ sin ωτ ) .

(2.11)

In view of some positive terms of (2.11), to show ddμ τ > 0, it is sufficient to verify that

(μ cos ωτ + ω sin ωτ )(2μ + A1 − A3 τ e−μτ cos ωτ ) + (μ sin ωτ + ω cos ωτ )(2ω + A3 τ e−μτ sin ωτ ) > 0.

Now at τ = τ0 , μ = 0, ω =

−ω

∗2

ω∗ ,

(2.12)

Eq. (2.10) becomes

+ A2 + A3 cos ω∗ τ0 = 0, A1 ω∗ − A3 sin ω∗ τ0 = 0,

(2.13)

ω∗4 + (A21 − 2A2 )ω∗2 + A22 − A23 = 0. and (2.12) becomes

ω∗ sin ω∗ τ0 (A1 − A3 τ0 cos ω∗ τ0 ) + ω∗ cos ω∗ τ0 (2ω∗ + A3 τ0 sin ω∗ τ0 ) > 0. Simplifying the above inequality, we have to check

A1 ω∗ sin ω∗ τ0 + 2ω∗2 cos ω∗ τ0 > 0.

(2.14)

Substituting the first and second equations of (2.13) into (2.14), the inequality (2.14) becomes

1 1 2 ∗2 A ω − [2ω∗2 (A2 − 2ω∗2 )] > 0. A3 1 A3 Since A3 > 0, we only have to verify

2ω∗4 + (A21 − 2A2 )ω∗2 > 0.

(2.15)

Substituting the third equation of (2.13) into (2.15), we only need to verify

ω∗4 − A22 + A23 > 0.

(2.16)

The above inequality holds true because 0 < A2 < A3 . The conditions for Hopf bifurcation are then satisfied. Moreover, by applying the implicit function theorem to (λ, τ ), it is easy to see that there exist a constant σ > 0 and a smooth curve λ(τ ) : (τ0 − σ , τ0 + σ ) → C such that λ(τ0 ) = iω∗ and (iω∗ , τ0 ) = 0 for all τ ∈ (τ0 − σ , τ0 + σ ).  2.2. The direction and stability of the Hopf bifurcation By Theorem 2.1, we obtained conditions under which a family of periodic solutions bifurcate from the equilibrium E2 at the critical value of τ 0 . As pointed out in Hassard et al. [36], it is interesting and important to determine the direction, stability and period of these periodic solutions. To this aim, we shall derive explicit formulas that can be used to determine the properties of the emergent periodic solutions by means of the normal form theory and the center manifold theory, which have been applied in Hassard et al. [36]. In the following analysis, we always assume that the condition (2.6) holds, and ± iω are the purely imaginary roots, where ω = ω∗ . To simplify our analysis of the Hopf bifurcation around E2 , we place this equilibrium at the origin by performing a coordinate shift through a change of variables ui − uˆ = u¯ i , vi − vˆ = v¯ i , and τ − τ0 = ν . To minimize the number of notations,

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

5

we omit the bars of u¯ i and v¯ i , and rescale the time by t → τt to normalize the delay. Then the system (1.1) without the initial conditions can be reformulated as



u˙ i (t ) = (τ0 + ν ) −(1 + vˆ 2 )ui (t ) − 2bvi (t ) − 2vˆ ui (t )vi (t ) − uˆv2i (t ) − ui (t )v2i (t ) + Du

n 



Li j u j (t ) ,

j=1



v˙ i (t ) = (τ0 + ν ) vˆ ui (t − 1 ) + bvi (t ) + 2vˆ ui (t − 1 )vi (t ) + uˆv (t ) + ui (t − 1 )v (t ) + Dv 2

2 i

2 i

n 

 Li j v j (t ) .

j=1

Since ± iω are the only purely imaginary roots of the above linear operator, by using the orthogonality of the eigenvectors of Lij , the above equation can be rewritten by





u˙ i (t ) = (τ0 + ν ) −(1 + vˆ 2 )ui (t ) − 2bvi (t ) − 2vˆ ui (t )vi (t ) − uˆv2i (t ) − ui (t )v2i (t ) + Du α ui (t ) ,

  v˙ i (t ) = (τ0 + ν ) vˆ 2 ui (t − 1 ) + bvi (t ) + 2vˆ ui (t − 1 )vi (t ) + uˆv2i (t ) + ui (t − 1 )v2i (t ) + Dv α vi (t ) .

(2.17)

We denote the solution of (2.17) by X(t ) = (ui (t ), vi and set Xt (θ ) = X(t + θ ), θ ∈ [−1, 0], Xt ∈ C ([−1, 0],  C. It can be verified that the system (2.17) undergoes a Hopf bifurcation at the critical value τ = 0 to which the corresponding eigenvalues are ± iω∗ . We rewrite (2.17) as the following functional differential equation

(t ))T ,

X˙ t = Lν (Xt ) + Rν (Xt ), where Lν : C →

R2

(2.18)

and Rν : C →



R2 )

R2

are given by



 φ1 ( 0 ) φ1 (−1 ) + ( τ0 + ν ) M 1 , φ2 ( 0 ) φ2 (−1 )



 0 0 −1 − vˆ 2 + Du α −2b M0 = , M = , 1 0 b + Dv α vˆ 2 0

 −2vˆ φ1 (0 )φ2 (0 ) − uˆφ22 (0 ) − φ1 (0 )φ22 (0 ) R ν ( φ ) = ( τ0 + ν ) . 2vˆ φ1 (−1 )φ2 (0 ) + uˆφ22 (0 ) + φ1 (−1 )φ22 (0 ) L ν ( φ ) = ( τ0 + ν ) M 0

(2.19)

Here φ = (φ1 , φ2 )T ∈ C ([−1, 0], R2 )  C. We choose

η ( θ , ν ) = ( τ0 + ν ) M 0 δ ( θ ) + ( τ0 + ν ) M 1 δ ( θ + 1 ) ,

(2.20)

where δ (·) is the Dirac function. By the Riesz Representation Theorem, we have

Lν ( φ ) =



0 −1

φ (θ )dη (θ , ν ) for φ ∈ C,

(2.21)

where

d η ( θ , ν ) = ( τ0 + ν ) M 0 δ ( θ ) d θ + ( τ0 + ν ) M 1 δ ( θ + 1 ) d θ . Next, we define two operators Aν and Rν on C 1 ([−1, 0], R2 )  C 1

 



Aν φ ( θ ) =

⎧ ⎨ ⎩





Fν φ (θ ) =

dφ , dθ



0 −1

θ ∈ [−1, 0 ), φ ( ξ )d η ( ξ , ν ),

(2.22)

θ = 0;

0,

θ ∈ [−1, 0 ),

Rν ( φ ),

θ = 0.

(2.23)

Then the system (2.18) is transformed into

X˙ t = Aν (Xt ) + Fν (Xt ).

(2.24) A∗

Note that the adjoint operator ν of Aν is defined as

(Aν ψ )(θ ) = ∗

∗

⎧ dψ ⎨ − dθ ∗ , ⎩



0 −1

θ ∗ ∈ (0, 1],

ψ (−ξ )dηT (ξ , ν ),

θ ∗ = 0,

(2.25)

where ψ ∈ C 1 ([0, 1], R2 ). We take the complex space C2 instead of merely R2 as the solution space of (2.24), and define the following bilinear form T

ψ ( θ ∗ ), φ ( θ ) = ψ ( 0 )φ ( 0 ) −



0

 θ

θ =−1 ξ =0

T

ψ ( ξ − θ )d η ( θ )φ ( ξ )d ξ ,

(2.26)

6

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

for φ ∈ C ([−1, 0], C2 ), ψ ∈ C 1 ([0, 1], C2 ), where η (θ ) = η (θ , 0 ), to determine the coordinates of the center manifold near the origin of (2.18). Note that the overlines in Eq. (2.26) stand for complex conjugate. By the analysis in Lemma 2.1 (Appendix A), we know that ± iω∗ are eigenvalues of (2.5). Therefore, ± iω∗ are eigenvalues of A0 and A∗0 , respectively. We first need to compute the eigenvectors q(θ ) and q∗ (θ ∗ ) of A0 and A∗0 corresponding to iω∗ and −iω∗ , respectively, that is,

A0 q(θ ) = iω∗ q(θ ) and A∗0 q∗ (θ ∗ ) = −iω∗ q∗ (θ ∗ ),

(2.27)

which satisfy the normalized conditions q∗ , q = 1 and q∗ , q¯ = 0. Thus, we set



∗ q ( θ ) = q ( 0 ) e iω θ =

∗ ∗ q 1 iω ∗ θ e and q∗ (θ ∗ ) = q∗ (0 )e−iω θ = q2



q∗1 −iω∗ θ ∗ e , q∗2

(2.28)

for θ ∈ [−1, 0 ), θ ∗ ∈ (0, 1]. Noting that q(0) can also be solved by the characteristic Eq. (2.5), we have



q (θ ) =

2b(−1 − vˆ 2 + iω∗ ) iωθ (1 + vˆ 2 )2 + ω∗2 e . 1

Similarly substituting (2.28) into (2.27), we have



q∗ ( θ ∗ ) = ρ

2b(−b + iω∗ ) ∗ e−iωθ . b2 + ω∗2 1

In order to ensure q∗ , q = 1, we need to determine the value of ρ . From (2.26), we have

q∗ (s ), q(θ ) = q¯ ∗ (0 )q(0 ) −



0

 θ

−1

0

q¯ ∗ (ξ − θ )dη (θ )q(ξ )dξ 

0

 θ

−1

0

ρ¯ (q¯ ∗1 , q¯ ∗2 )e−iω (ξ −θ ) dη (θ )(q1 , q2 )T eiω ξ dξ     ∗ ∗ = ρ¯ q¯ ∗1 q1 + τ0 e−iω τ0 (0q1 + 0q2 ) + ρ¯ q¯ ∗2 q2 + τ0 e−iω τ0 (vˆ 2 q1 + 0q2 ) = ρ¯ (q¯ ∗1 , q¯ ∗2 )(q1 , q2 )T −

∗

∗

∗ = ρ¯ (q¯ ∗1 q1 + q¯ ∗2 q2 + q¯ ∗2 q1 vˆ 2 τ0 e−iω τ0 ).

Thus we can choose ρ as

ρ¯ = 1/(q¯ ∗1 q1 + q¯ ∗2 q2 + q¯ ∗2 q1 vˆ 2 τ0 e−iω τ0 ). ∗

Next, by using q and q∗ , we can construct a coordinate for the center manifold C0 at ν = 0 using the method of Hassard et al. [36]. Denote by xt = xt (θ ) the solution of (2.24) when ν = 0, and define

z(t )  q∗ , xt ,

(2.29)

W (t, θ )  xt (θ ) − 2Re{z(t )q(θ )}.

(2.30)

In terms of the center manifold reduction, we have W (t, θ ) = W (z(t ), z(t ), θ ) on C0 . According to the center eigenspace at the equilibrium, we further have

W (t, θ ) = W (z(t ), z(t ), θ )  W20 (θ )

2

z2 z z3 + W11 (θ )zz + W02 (θ ) + W30 (θ ) + · · · , 2 2 6

(2.31)

where z and z¯ are the local coordinates of the center manifold C0 in directions of q and q∗ , respectively. Note that W is real if X(t ) is real. We are only concerned with the real solutions. For solution X(t ) ∈ C0 of (2.24) with ν = 0, the Eq. (2.30) yields

z˙ (t ) = q∗ , x˙ t



= iω∗ z + q¯∗ (0 )F0 W (z, z¯, θ ) + 2Re{z(t )q(θ )}



 iω∗ z + q¯∗ (0 )F0 (z, z¯ ).

(2.32)

We rewrite the above equation as

z˙ (t ) = iω∗ z + g(z, z¯ ),

(2.33)

where

g(z, z¯ ) = q¯∗ (0 )F0 (z, z¯ ) g20 2 g02 2 g30 3 g21 2 g12 2 g03 3 g40 4  z + g11 zz¯ + z + z z¯ + zz¯ + z + ··· . z¯ + z¯ + 2 2 6 2 2 6 24

(2.34)

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

7

Fig. 1. Delay-driven instability in the parameter space spanned by parameter a and b. The right panel is a zoom of the left panel for b ∈ [0, 0.1]. Color bars indicate the critical value of τ 0 above which instability of the coexistence equilibrium E2 occurs. Attention is paid to the parameter space where plant can survive.

The explicit expressions of gij are derived in Appendix B. After a standard bifurcation stability argument as in Hassard et al. [36], we have





1 i 1 g20 g11 − 2|g11 |2 − |g02 |2 + g21 . 2ω 3 2 A(τ0 , ν ) = ω∗ τ0 + ν Im{λ (0 )}, C1 (0 ) =

1 Im{C1 (0 )}Re{λ (0 )}ν Re{C1 (0 )} − . 2 2 A ( τ0 , ν ) 2π T (τ0 , ν ) = ∗ (1 + N (τ0 )ν + o(|ν| ) ), B ( τ0 , ν ) =

ω τ0 Im{C1 (0 )}Re{λ (0 )} − Re{C1 (0 )}Im{λ (0 )} N ( τ0 ) = . ω∗ τ0 Re{C1 (0 )}

(2.35)

We now summarize the results about the stability and direction of Hopf bifurcation in Theorem 2.2. Since the proof is similar to that of Hassard et al. [36], we hence omit the details. Theorem 2.2. If Re{C1 (0)} = 0, then the system (1.1) has a branch of Hopf bifurcating solutions for τ = τ0 + ν with ν satisfying ν Re{λ (0)}B(τ 0 , ν ) < 0, where C1 (0) and B(τ 0 , ν ) will be defined in (2.35). Also the bifurcating periodic solutions have the following properties: (i) Re{C1 (0)} determines the stability of the bifurcating periodic solutions. The bifurcating periodic solutions are orbitally stable (resp., unstable) if Re{C1 (0)} < 0 (resp., Re{C1 (0)} > 0). Re{C (0 )} Re{C (0 )} Re{C (0 )} (ii) − Re{λ1 (0 )} determines the direction of the Hopf bifurcation. If − Re{λ1 (0 )} > 0 (resp., − Re{λ1 (0 )} < 0), the bifurcating periodic solutions exist for ν > 0 (resp., ν < 0), which is also called a supercritical (resp., subcritical) Hopf bifurcation. (iii) The period of the bifurcating periodic solution is ω2∗πτ as ν = 0, the period T(τ 0 , ν ) is increasing in parameter ν (resp., 0 decreasing) if N(τ 0 ) > 0 (resp., N(τ 0 ) < 0). Here T(τ 0 , ν ) and N(τ 0 ) are defined in (2.35). 3. Numerical simulations In this section, we carry out numerical simulations to demonstrate our analytical findings. To this end, we consider the following set of parameters based on the real data reported in [14]

a = 2, b = 0.45, Du = 100, Dv = 1.

(3.1)

Recall that the conditions (2.6) and (2.8) basically determine when Hopf bifurcation induced instability occurs around the positive equilibrium E2 . Fig. 1 shows how the critical time delay τ 0 varies with parameters a and b. As expected, large a and b promote the emergence of Hopf bifurcation since small time delay is required. Ecologically speaking, in case of large amount rainfall or high plant mortality, long time lag in water uptake (e.g., less frequent rainfall of approximately constant annual amount) can easily result in the instability of coexistence state between plant and ground water. Theorem 2.2 implies that spatial patterns can develop if time delay goes beyond a threshold. For selected parameters in (3.1), we find that the threshold is τ0 = 1.3496. For numerical simulation, we take a time delay τ = 1.4 and consider a

8

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

Fig. 2. Time series dynamics of plant at different nodes. Color bar indicates plant densities. Parameters are in (3.1) and τ = 0.

Fig. 3. Time series dynamics of plant at different nodes. Color bar indicates plant densities. Parameters are in (3.1) and τ = 1.4.

Watts–Strogatz network of 100 nodes with an average degree of 4 and a link rewiring probability of 0.15. The initial data is taken as a uniformly distributed random perturbation around the equilibrium E2 = (uˆ, vˆ ) = (0.107, 4.2067 ). The model is simulated by the standard Runge–Kutta scheme for the time integration. When τ = 0, the time evolutions of plant density at different habitats are shown in Fig. 2. It shows that the system (1.1) is asymptotically stable. When τ = 1.4, the time evolutions of plant density at different habitats are shown in Fig. 3 and density patterns on the network at different time instants are displayed in Fig. 4. Furthermore, we show, for a particular habitat, the time series dynamics of plant density and rainfall amount in Fig. 5. Our numerical simulations clearly demonstrate the occurrence of Hopf bifurcations and emergence of spatial patterns across the entire landscape. Moreover, the water and plant dynamics are in phase opposition. Larger time delay does not alter this observation but increases the magnitude of oscillation. It means that plants are more vulnerable to longer delay in water consumption, and more prone to go extinct, if stochastic effects are taken into account.

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

9

Fig. 4. Plant density (color bar) on the network at different time instants. The node size corresponds to node degree. Parameters are in (3.1) and τ = 1.4.

Fig. 5. Densities of plant and water (left panel) and the stable attractor (right panel) in one of the considered habitats. Parameters are in (3.1) and τ = 1.4.

To check the stability and direction of the Hopf bifurcations, we further compute the first Lyapunov number Re{C1 (0)}. Direct calculation shows that Re{C1 (0 )} = −20.4375 < 0, which means that the periodic solution bifurcating from the spatially uniform equilibrium is stable. 4. Discussion We have developed a theoretical framework for studying spatiotemporal dynamics of a semiarid vegetation model on a network domain. The network stands for fragmented habitats connected by dispersal of plant seed. Applying standard

10

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

Fig. 6. Turing instability region (1 ) for the reaction–diffusion system (4.2) with d1 = 242.5 and d2 = 1, (1 + 2 ) for the reaction–advection–diffusion system (4.1) with ν = 182.5, and Hopf bifurcation region (1 + 2 + 3 ) for the network-organized model (1.1).

stability analysis and the center manifold theory, we investigated the direction and stability of delay induced Hopf bifurcation, which leads to the emergence of oscillations that are periodic in time. An interesting finding is that small time delay can easily induce the onset of Hopf bifurcation when plant mortality is high, given the same rainfall rate (Fig. 6). Spatial dynamics of semiarid vegetation system has traditionally been investigated under the framework of reaction– advection–diffusion system such as [14,15]

∂u = a − u − uv2 + d1 ∇ u, ∂t ∂v = uv2 − bv + d2 ∇ 2 v. ∂t

(4.1)

Kealy and Wollkind [37] studied a model variant of above reaction–diffusion system

∂u = a − u − uv2 + d1 ∇ 2 u, ∂t ∂v = uv2 − bv + d2 ∇ 2 v. ∂t

(4.2)

In these models, water flow was described as either an advection or diffusion on a spatially continuous domain. In our model, water flow was described as a motion across different habitats. Interestingly, Lemma 2.1 shows that our system (1.1) never undergoes a transition toward a patterned state in the absence of time delay, which is contradictory to the reported results in [14,15]. This contradiction is due to the presence of advection terms for water in their models. Moreover, we compared the systems (4.1) and (4.2) with our network model in terms of equilibrium stability in the parameter space spanned by plant mortality and rainfall, and summarized the results in Fig. 6. The figure shows that the presence of advection term increases Turing instability region, and we did not find Turing instability in our network-organized model, however, we did find a substantial region of Hopf bifurcation, which highlights that time delay has a significantly destabilizing effect, leading to oscillatory vegetation dynamics that are periodic in time. Semi-arid ecosystems are an important composition to earth’s ecosystem, which together with rangelands make up around one quarter of the earth’s landscape. A typical example is Australia with almost 75% of the total land mainly being classified as arid and semi-arid lands. A pronounced scenario occurring on these lands is the increased land fragmentation due to the massive development of water sources, fences, roads, and agriculture, in addition to climate and natural resource bases. Representative areas are the Dalrymple Shire, South-west Queensland in Australia and the Northern Great Plains region in USA [38]. Fragmentation has substantially detrimental effects on animals which could previously exploit spatial heterogeneity to compensate for temporal heterogeneity but are now constrained to limited resource availability, and on ecosystems functions by altering species composition through alien invasion along roadside [39]. Besides these negative impacts, the frequent observation of spatial heterogeneity of vegetation density on networks (connected fragmented

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

11

parcels) in our work implies that resource availability on some parcels could be further suppressed, reducing local survival/reproduction rates, leading to long-term passive impacts on populations. These passive impacts could be exacerbated due to increased time delay arising from less frequent rainfall. However, our proposed vegetation model is based on an undirected network, which is actually a toy model that does not ground enough physical facts. In reality, water diffuses in a direction of downstream and plant seeds can diffuse via different approaches such as wind and animal paths. Thus a promising extension of the current work is to develop a network model with directional diffusion for water and multigraph approach for seed diffusion.

Acknowledgements The work is supported by the PRC Grant of NSFC (61877052, 11571301, and 11871065). Zhi Ling and Lai Zhang are additionally supported by the NSF of Jiangsu Province (BK20151305 and BK20181450). Lai Zhang further acknowledges the financial support by the Jiangsu Distinguished Professor Program, and the Yangzhou Talent Program of ‘LvYangJinFeng’.

Appendix A. Proof of Lemma 2.1 Proof. (i) We first show the case τ = 0. The characteristic Eq. (2.5) is reduced to

λα 2 + (1 − b + vˆ 2 − Du α − Dv α )λα − b(1 + vˆ 2 ) + bDu α − Dv (1 + vˆ 2 )α + Du Dv α 2 + 2bvˆ 2 = 0. In view of (2.6), it follows from vˆ = we have

√2b

a−

a2 −4b2

and a > 2b that vˆ > 1. By using 4bDu Dv (vˆ 2 − 1 ) > (bDu − Dv (vˆ 2 + 1 ))2 of (2.6),

−b(1 + vˆ 2 ) + bDu α − Dv (1 + vˆ 2 )α + Du Dv α + 2bvˆ 2 > 0. 2

Hence the roots of characteristic Eq. (2.4) have a negative real part if and only if

1 − b + vˆ 2 − Du α − Dv α > 0. Since α ≤ 0, the above inequality is guaranteed if vˆ 2 + 1 − b > 0. It follows from vˆ =

√2b

a−

a2 −4b2

that

√ 4b2 2a2 − 4b2 + 2a a2 − 4b2 2a2 − 4b2 a2 vˆ 2 = = > = − 1 > b − 1. √ 2 2 4b 4b 2b2 (a − a2 − 4b2 )2

The last inequality holds if a > 2b and a > b2 , which is guaranteed by the condition (2.6). (ii) Let ± iω be a pair of pure imaginary roots. Substituting ± iω into characteristic Eq. (2.5) and separating real and imaginary parts, we have

−ω2 + A2 + A3 cos ωτ = 0, A1 ω − A3 sin ωτ = 0,

(A.1)

which lead to

ω4 + (A21 − 2A2 )ω2 + A22 − A23 = 0. The above equation is essentially a quadratic equation with respect to ω2 and it has a unique positive real root ω∗2 if A22 − A23 < 0. In order to show A22 − A23 < 0, we choose α = 1 such that α = 0. Thus direct computation shows that if the condition (2.6) together with

| − b(1 + vˆ 2 ) + bDu α − Dv (1 + vˆ 2 )α + Du Dv α 2 | = b(1 + vˆ 2 ) < 2bvˆ 2 is satisfied, we have A22 − A23 < 0. Moreover, by solving the above equation, the expression of ω∗ is written in the form of (2.7). ω∗2 is always real and positive. Hence the corresponding τ j of (2.8) can be obtained by substituting (2.7) into (A.1).  Appendix B. Derivation of gij From (2.30) and (2.31), we obtain that

xt (θ ) = W (t, θ ) + 2Re{z(t )q(θ )} = W20 (θ )

z2 z¯2 + W11 (θ )zz¯ + W02 (θ ) + zq + z¯q¯ + · · · 2 2

= W20 (θ )

z2 z¯2 ∗ ∗ + W11 (θ )zz¯ + W02 (θ ) + (q1 , q2 )T eiω θ z + (q¯ 1 , q¯ 2 )T e−iω θ z¯ + · · · . 2 2

(B.1)

12

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

Substituting (B.1) into (2.34), we have



g(z, z¯ ) = q¯∗ (0 )F0 (z, z¯ ) = q¯∗ (0 )F0 xt (θ )





= ρ¯ (

q¯ ∗1 , q¯ ∗2



−2vˆ φ1 (0 )φ2 (0 ) − uˆφ22 (0 ) − φ1 (0 )φ22 (0 ) ) 2vˆ φ1 (−1 )φ2 (0 ) + uˆφ22 (0 ) + φ1 (−1 )φ22 (0 )

g20 2 g02 2 g30 3 g21 2 g12 2 g03 3 g40 4 z + g11 zz¯ + z + z z¯ + zz¯ + z + h.o.t., z¯ + z¯ + 2 2 6 2 2 6 24



where h.o.t. stands for higher order terms, and





∗ g20 = 2ρ¯ q¯ ∗1 (−uˆ − 2vˆ q1 ) + q¯ ∗2 (uˆ + 2vˆ q1 e−iω τ0 ) ,





∗ ∗ g11 = ρ¯ q¯ ∗1 (−2uˆ − 2vˆ q1 − 2vˆ q¯ 1 ) + q¯ ∗2 (2uˆ + 2vˆ q1 e−iω τ0 + 2vˆ q¯ 1 eiω τ0 ) ,





∗ g02 = 2ρ¯ q¯ ∗1 (−uˆ − 2vˆ q¯ 1 ) + q¯ ∗2 (uˆ + 2vˆ q¯ 1 eiω τ0 ) ,





∗ ∗ (1 ) (2 ) (1 ) (2 ) g30 = 6ρ¯ q¯ ∗1 (−q1 − vˆW20 (0 ) − (vˆ q1 + uˆ )W20 (0 )) + q¯ ∗2 (q1 e−iω τ0 + vˆW20 (−1 ) + (vˆ q1 e−iω τ0 + uˆ )W20 (0 )) ,

 

(1 ) (2 ) g21 = 2ρ¯ q¯ ∗1 − q¯ 1 − 2q1 − vˆW20 (0 ) − (uˆ + vˆ q¯ 1 )W20 (0 ) − 2vˆW11(1) (0 ) − (2uˆ + 2vˆ q1 )W11(2) (0 )





∗ ∗ (1 ) (2 ) + q¯ ∗2 q¯ 1 + 2q1 + vˆW20 (−1 ) + (uˆ + vˆ q¯ 1 eiω τ0 )W20 (0 ) + 2vˆW11(1) (−1 ) + (2uˆ + 2vˆ q1 e−iω τ0 )W11(2) (0 )





(1 ) (1 ) (2 ) g12 = 2ρ¯ q¯ ∗1 (−q1 − 2vˆW11 (0 ) − (2vˆ q¯ 1 + 2uˆ )W11(2) (0 ) − vˆW02 (0 ) − (uˆ + vˆ q1 )W02 (0 ))

,



∗ ∗ ∗ (1 ) (1 ) (2 ) + q¯ ∗2 (q1 e−iω τ0 + 2vˆW11 (−1 ) + (2vˆ q¯ 1 eiω τ0 + 2uˆ )W11(2) (0 ) + vˆW02 (−1 ) + (uˆ + vˆ q1e−iω τ0 )W02 (0 )) ,





∗ ∗ (1 ) (2 ) (2 ) (1 ) (2 ) g03 = 6ρ¯ q¯ ∗1 (−q¯ 1 − vˆW02 (0 ) − (vˆ q¯ 1 + uˆ )W20 (0 )) + q¯ ∗2 (q¯ 1 eiω τ0 + vˆW20 (−1 ) + vˆ q¯ 1 eiω τ0 W02 (0 ) + uˆW02 (0 )) ,



1 (1 ) uˆ (2) vˆ (1) (2 ) (2 ) g40 = 24ρ¯ q¯ ∗1 − W20 (0 ) − W20 (0 )W20 (0 ) − W20 (0 )W20 (0 ) 2 2 2



+ q¯ ∗2

1 (1 ) uˆ (2) vˆ (1) (2 ) (2 ) W (−1 ) + W20 (−1 )W20 (0 ) + W02 (0 )W02 (0 ) 2 20 2 2



! ,

(k ) (k ) where W20 and W11 are the kth components of W20 and W11 , respectively, for −1 ≤ θ ≤ 0. Since g21 depends on W20 (θ ) and W11 (θ ), we need to find the values of W20 (θ ) and W11 (θ ). From (2.24) and (2.30), we have

W˙ = x˙ t − z˙ q − z¯˙ q¯

"

=









A0W − 2Re q¯ ∗ (0 )F0 q(θ ) ,

−1 ≤ θ < 0,

A0W − 2Re q¯ ∗ (0 )F0 q(θ ) + F0 ,

θ = 0,

 A0W + H (z, z¯, θ ),

(B.2)

where

H (z, z¯, θ ) = H20

z2 z¯2 + H11 zz¯ + H02 + · · · . 2 2

(B.3)

From (2.31), we have

W˙ = W˙ z z˙ (t ) + W˙ z¯ z¯˙ (t ) = (W20 (θ )z + W11 (θ )z¯ + · · · )(iω∗ z(t ) + g(z, z¯ )) + (W11 (θ )z + W02 (θ )z¯ + · · · )(−iω∗ z¯ (t ) + g¯ (z, z¯ )). Substituting (2.31) to (B.2) yields that



W˙ = A0 W20 (θ )

z2 z¯2 + W11 (θ )zz¯ + W02 (θ ) + · · · 2 2

= (A0W20 (θ ) + H20 (θ ))

 + H20 (θ )

(B.4)

z2 z¯2 + H11 (θ )zz¯ + H02 (θ ) + · · · 2 2

z2 z¯2 + (A0W11 (θ ) + H11 (θ ))zz¯ + (A0W02 (θ ) + H02 (θ )) + · · · . 2 2

(B.5)

Comparing the coefficients of z2 and zz¯ between (B.4) and (B.5), we get

(A0 − 2iω∗ I )W20 (θ ) = −H20 (θ ), A0W11 (θ ) = −H11 (θ ). For θ ∈ [0, 1], it follows from (2.30), (2.31), (B.3) and (B.4) that

(B.6)

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

13

H (z, z¯, θ ) = −q¯ ∗ (0 )F0 q(θ ) − q∗ (0 )F¯0 q¯ (θ ) = −g(z, z¯ )q(θ ) − g¯ (z, z¯ )q¯ (θ ),





z2 z¯2 = − g20 + g11 zz¯ + g02 + · · · q(θ ) − 2 2





z2 z¯2 g¯ 20 + g¯ 11 zz¯ + g¯ 02 + · · · q¯ (θ ). 2 2

(B.7)

Comparing the coefficients of z2 and zz¯ between (B.3) and (B.7), we get

H20 (θ ) = −g20 q(θ ) − g¯ 20 q¯ (θ ), H11 (θ ) = −g11 q(θ ) − g¯ 11 q¯ (θ ).

(B.8)

From the definition of Aν (θ ) and (B.6) and (B.8), we have

W˙ 20 (θ ) = 2iω∗W20 (θ ) + g20 q(θ ) + g¯ 02 q¯ (θ ). Since q(θ ) = (q1 , q2 )T eiωθ , we obtain

ig20

W20 (θ ) =

ω

q(0 )eiωθ +

ig¯ 02 q¯ (0 )e−iωθ + G1 e2iωθ , 3ω

(B.9)

where G1 = (G1(1 ) , G1(2 ) )T is a constant vector. Similarly, we have

−ig11

W11 (θ ) =

ω

q(0 )eiωθ +

ig¯ 11

ω

q¯ (0 )e−iωθ + G2 ,

(B.10)

where G2 = (G2(1 ) , G2(2 ) )T is a constant vector. Now, we shall find the values of G1 and G2 . From the definition of A0 and (B.6), we have



0 −1

and



0 −1

dη (θ )W20 (θ ) = 2iω∗W20 (0 ) − H20 (0 ),

(B.11)

dη (θ )W11 (θ ) = −H11 (0 ),

(B.12)

where η (θ ) = η (θ , 0 ). In view of (B.2), we induce that when θ = 0,

H (z, z¯, 0 ) − 2Re{q¯ ∗ (0 )F0 q(0 )} + F0 = −g(z, z¯ )q(0 ) − g¯ (z, z¯ )q¯ (0 ) + F0 . Then we have

H20





z2 z2 z¯2 z¯2 + H11 zz¯ + H02 + · · · = − g20 + g11 zz¯ + g02 + · · · q(0 ) 2 2 2 2



−

g¯ 20

Comparing both sides of (B.13), we obtain





(B.14)



−2uˆ − 2vˆ q1 − 2vˆ q¯ 1 . ∗ 2uˆ + 2vˆ q1 e−iω + 2vˆ q¯ 1 eiω τ0

Since iω∗ is the eigenvalue of A0 and q(0) is the corresponding eigenvector, we get



iω I − ∗

and





0

e

iω ∗ θ

−1

−iω I − ∗



d η ( θ ) q ( 0 ) = 0,



0

e

−iω∗ θ

−1

dη (θ ) q¯ (0 ) = 0.

Therefore, substituting (B.9) and (B.11) into (B.14), we have

(2iω I − ∗





0

e −1

2 iω ∗ θ

dη (θ ))G1 =

(B.13)



−uˆ − 2vˆ q1 , ∗ uˆ + 2vˆ q1 e−iω τ0

H20 = −g20 q(0 ) − g¯ 02 q¯ (0 ) + H11 = −g11 q(0 ) − g¯ 11 q¯ (0 ) +



z2 z¯2 + g¯ 11 zz¯ + g¯ 02 + · · · q¯ (0 ) + F0 . 2 2

−uˆ − 2vˆ q1 ∗ uˆ + 2vˆ q1 e−iω τ0 ,

(B.15)

14

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

that is,



H ∗ G1 = where

H∗ =

−uˆ − 2vˆ q1 ∗ uˆ + 2vˆ q1 e−iω τ0 ,

2iω + 1 + vˆ 2 ∗ −ˆv2 e−iω τ0

Thus

2b 2iω − b .



−uˆ − 2vˆ q1 ∗ Det uˆ + 2vˆ q1 e−iω τ0 G1(1) =



,

Det(H ∗ )

2iω + 1 + vˆ 2 ∗ Det −ˆv2 e−iω τ0 G1(2) =

2b 2iω − b

−uˆ − 2vˆ q1 ∗ uˆ + 2vˆ q1 e−iω τ0

.

Det(H ∗ )

In a similar way, substituting (B.10) and (B.12) into (B.15), we have



P ∗ G2 = where

P∗ =

Thus



−2uˆ − 2vˆ q1 − 2vˆ q¯ 1 , ∗ ∗ 2uˆ + 2vˆ q1 e−iω τ0 + 2vˆ q¯ 1 eiω τ0

1 + vˆ 2 ∗ −ˆv2 e−iω τ0



2b . −b

Det G2(1) =

−2uˆ − 2vˆ q1 − 2vˆ q¯ 1 ∗ ∗ 2uˆ + 2vˆ q1 e−iω τ0 + 2vˆ q¯ 1 eiω τ0



1 + vˆ 2 ∗ Det −ˆv2 e−iω τ0 G2(2) =

Det(P ∗ )



2b −b

,

−2uˆ − 2vˆ q1 − 2vˆ q¯ 1 ∗ ∗ 2uˆ + 2vˆ q1 e−iω τ0 + 2vˆ q¯ 1 eiω τ0 Det(P ∗ )

.

Therefore, we can determine W20 (θ ) and W11 (θ ) from (B.9) and (B.10). References [1] J. Bromley, J. Brouwer, A.P. Barker, S.R. Gaze, C. Valentin, The role of surface water distribution in an area of patterned vegetation in a semiarid environment, south-west Niger, J. Hydrol. 198 (1997) 1–19. [2] M.R. Aguiar, O.E. Sala, Patch structure, dynamics and implications for the functioning of arid eco-systems, Trends Ecol. Evol. 14 (1999) 273–277. [3] J.C. Leprun, The influences of ecological factors on tiger bush and dotted bush patterns along a gradient from mali to northern Burkina Faso, Catena 37 (1999) 25–34. [4] P. Couteron, O. Lejeune, Periodic spotted patterns in semiarid vegetation explained by a propagation-inhibition model, J. Ecol. 89 (2001) 616–628. [5] J.V. Hardenberg, E. Meron, M. Shachak, Y. Zarmi, Diversity of vegetation patterns and desertification, Phys. Rev. Lett. 87 (2001) 198101. [6] J. Tews, U. Brose, V. Grimm, K. Tielbörger, M.C. Wichmann, M. Schwager, F. Jeltsch, Animal species diversity driven by habitat heterogeneity/diversity: the importance of key stone structures, J. Biogeogr. 31 (2004) 79–92. [7] M. Rietkerk, S.C. Dekker, P.C. deRuiter, J. vandeKoppel, Self-organized patchiness and catastrophic shifts in ecosystems, Science 305 (2004) 1926–1929. [8] M. Rietkerk, M.C. Boerlijst, F. vanLangevelde, R. HilleRisLambers, J. vandeKoppel, L. Kumar, H.H.T. Prins, A.M. DeRoos, Self-organization of vegetation in arid ecosystems, Am. Nat. 160 (2002) 524–530. [9] A.I. Borthagaray, M.A. Fuentes, P.A. Marquet, Vegetation pattern formation in a fog dependent ecosystem, J. Theor. Biol. 265 (2010) 18–26. [10] F. Borgogno, P. Dodorico, F. Laio, L. Ridolfi, Mathematical models of vegetation pattern formation in ecohydrology, Rev. Geophys. 47 (2009) 1–36. [11] E. Gilad, M. Shachak, E. Meron, Dynamics and spatial organization of plant communities in water limited systems, Theor. Popul. Biol. 72 (2007) 214–230. [12] P.M. Saco, G.R. Willgoose, G.R. Hancock, Ecogeomorphology and vegetation patterns in arid and semiarid regions, Hydrol. Earth Syst. Sci. Discuss. 3 (2006) 2559–2593. [13] R. Lefever, O. Lejeune, On the origin of tiger bush, Bull. Math. Biol. 59 (1997) 263–294. [14] C.A. Klausmeier, Regular and irregular patterns in semiarid vegetation, Science 284 (1999) 1826–1828. [15] J.A. Sherratt, An analysis of vegetative stripe formation in semi-arid landscape, J. Math. Biol. 51 (2005) 183–197. [16] A. Turing, The chemical basis of morphogenesis, Phil. Trans. Royal. Soc. B 237 (1952) 37–72. [17] L.A. Segel, J.L. Jackson, Dissipative structure: an explanation and an ecological example, J. Theor. Biol. 37 (1972) 545–559. [18] M. Banerjee, S. Petrovskii, Self-organised spatial patterns and chaos in a ratio-dependent predator–prey system, Theor. Ecol. 4 (2011) 37–53.

C. Tian, Z. Ling and L. Zhang / Applied Mathematics and Computation 367 (2020) 124778

15

[19] M. Baurmann, T. Gross, U. Feudel, Instabilities in spatially extended predator–prey systems: spatio-temporal patterns in the neighborhood of Turing–Hopf bifurcations, J. Theor. Biol. 245 (2007) 220–229. [20] H.B. Shi, S.G. Ruan, Y. Su, J.F. Zhang, Spatiotemporal dynamics of a diffusive Leslie–Gower predator–prey model with ratio-dependent functional response, Int. J. Bifurc. Chaos 25 (2015) 1530014. [21] E. Barbera, G. Consolo, G. Valenti, Spread of infectious diseases in a hyperbolic reaction–diffusion susceptible-infected-removed model, Phys. Rev. E 88 (2013) 052719. [22] E.P. Zemskov, W. Horsthemke, Diffusive instabilities in hyperbolic reaction–diffusion equations, Phys. Rev. E 93 (2016) 032211. [23] H. Nakao, A.S. Mikhailov, Turing patterns in network-organized activator-inhibitor systems, Nat. Phys. 6 (2010) 544–550. [24] J. Petit, T. Carletti, M. Asllani, D. Fanelli, Delay-induced turing-like waves for one-species reaction–diffusion model on a network, Europhys. Lett. 111 (2015) 58002. [25] A. Berman, R.J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Academic Press, New York, 1979. [26] V. Volterra, Lecons sur la Théorie Mathématique de la Lutte par la Vie, Gauthier-Villars, Paris, 1931. [27] S. Sen, P. Ghosh, S.S. Riaz, D.S. Ray, Time-delay-induced instabilities in reaction–diffusion systems, Phys. Rev. E 80 (2009) 046212. [28] S. Boccaletti, D. Maza, H. Mancini, R. Genesio, F.T. Arecchi, Control of defects and space-like structures in delayed dynamical systems, Phys. Rev. Lett. 79 (1997) 5246. [29] C. Beta, M. Bertram, A.S. Mikhailov, H.H. Rotermund, G. Ertl, Controlling turbulence in a surface chemical reaction by time-delay autosynchronization, Phys. Rev. E 67 (2003) 046224. [30] A.S. Mikhailov, K. Showalter, Control of waves, patterns and turbulence in chemical systems, Phys. Rep. 425 (2006) 79–194. [31] P. Ghosh, Control of the Hopf–Turing transition by time-delayed global feedback in a reaction–diffusion system, Phys. Rev. E 84 (2011) 016222. [32] M. Bertram, A.S. Mikhailov, Pattern formation in a surface chemical reaction with global delayed feedback, Phys. Rev. E 63 (2001) 066102. [33] C. Tian, L. Zhang, Traveling wave governs the stability of spatial pattern in a model of allelopathic competition interactions, Chaos 22 (2012) 043136. [34] S. Grill, V.S. Zykov, S.C. Müller, Feedback-controlled dynamics of meandering spiral waves, Phys. Rev. Lett. 75 (1995) 3368. [35] J. Petit, M. Asllani, D. Fanelli, B. Lauwens, T. Carletti, Pattern formation in a two-component reaction–diffusion system with delayed processes on a network, Physica A 462 (2016) 230–249. [36] B. Hassard, N. Kazarinoff, Y. Wan, Theory and Application of Hopf Bifurcation, Cambridge University Press, Cambridge, 1981. [37] B.J. Kealy, D.J. Wollkind, A nonlinear stability analysis of vegetative turing pattern formation for an interaction-diffusion plant-surface water model system in an arid flat environment, Bull. Math. Biol. 74 (2012) 803–833. [38] K.A. Galvin, R.S. Reid, R.H. Behnke Jr., N.T. Hobbs (Eds.), Fragmentation in Semi-Arid and Arid Landscapes: Consequences for Human and Natural Systems, Springer, Netherlands, 2008. [39] F.B. Samson, F.L. Knopf, W.R. Ostlie, Great plains ecosystems: past, present, and future, Wildl. Soc. B. 32 (2004) 6–15.