Nonlocal interaction driven pattern formation in a prey–predator model

Applied Mathematics and Computation 308 (2017) 73–83

Contents lists available at ScienceDirect

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

Nonlocal interaction driven pattern formation in a prey–predator model Canrong Tian a, Zhi Ling b,∗, Lai Zhang c a

Department of Basic Sciences, Yancheng Institute of Technology, Yancheng 224003, China School of Mathematical Science, Yangzhou University, Yangzhou 225002, China c Department of Mathematics and Mathematical Statistics, Umeå University, Umeå SE-90187, Sweden b

a r t i c l e

i n f o

MSC: 35B36 35B40 92D25 Keywords: Nonlocal interactions Pattern formation Prey–predator model

a b s t r a c t A widely observed scenario in ecological systems is that populations interact not only with those living in the same spatial location but also with those in spatially adjacent locations, a phenomenon called nonlocal interaction. In this paper, we explore the role of nonlocal interaction in the emergence of spatial patterns in a prey–predator model under the reaction–diffusion framework, which is described by two coupled integro-differential equations. We first prove the existence and uniqueness of the global solution by means of the contraction mapping theory and then conduct stability analysis of the positive equilibrium. We find that nonlocal interaction can induce Turing bifurcation and drive the formation of stationary spatial patterns. Finally we carry out numerical simulations to demonstrate our analytical findings. © 2017 Elsevier Inc. All rights reserved.

1. Introduction The theory of pattern formation dates back to the pioneering work by Turing [1] in 1952, who studied a reaction– diffusion model consisting of two chemical species, an activator and an inhibitor. The activator stimulates the production of the other species, which, in turn, inhibits the generation of the activator. Turing found that as time advances the concentrations of the two species can evolve from an initial near-homogeneity into a spatially inhomogeneous pattern when the diffusion of the inhibitor is greater than that of the activator. This finding implies that the equilibrium of the nonlinear system, which is stable in the absence of diffusion, can be destabilised by diffusion, a counter-intuitive observation, as diffusion is usually considered to exhibit a stabilising effect. This scenario is termed as Turing instability and the emergent spatial patterns are called Turing patterns. Since Turing’s work, pattern formation has attracted great attention from both theorists and experimentalists during the last few decades. Levin and Segel [2] for the first time extended Turing idea from chemistry to ecology. They added diffusion to a planktonic system and theoretically investigated how diffusion leads to emergence of spatial patterns. Their method turns out to be a routine framework for studying pattern formation. Much research effort has been devoted to the mechanisms behind the emergence of spatial patterns and a variety of mechanisms have been found. These mechanisms can basically be classified into two categories. The first category is that spatial patterns emerge from locally stable kinetic dynamics due to the presence of self-diffusion [1–4], cross-diffusion [5–7], and advection [8–10]. The other category is that spatial patterns emerge from locally unstable kinetic dynamics (such as limit ∗

Corresponding author. E-mail address: [email protected] (Z. Ling).

http://dx.doi.org/10.1016/j.amc.2017.03.017 0 096-30 03/© 2017 Elsevier Inc. All rights reserved.

74

C. Tian et al. / Applied Mathematics and Computation 308 (2017) 73–83

cycle, chaos) through spatially inhomogeneous perturbations. Thus, factors that can destabilise the local kinetic dynamics can potentially be mechanisms for pattern formation, for instance, time delay [11–13] and single generation cycles [14]. It is worth pointing out that there is generally a requirement on the diffusion rates of involving species to trigger pattern formation in the former category but not in the latter [15]. Recently, population dynamics with free boundary (i.e., size of spatial domain can change as time advances) in a heterogeneous environment have been considered and it was shown that the free boundary can also affect the spatial distributions of population [16,17]. A common feature of above mentioned works is that a given individual interacts only with those living in the same spatial location (i.e., local interaction). In reality interactions can possibly occur within individuals living in adjacent locations (i.e., nonlocal interaction), which has been receiving growing attentions [18–20]. Moreover, some biological experiments demonstrated that bacteria cultures in Petri dishes release toxic substances or nutrients, which cause nonlocal interactions [21,22]. In contrast to local interactions where spatial and temporal scales of movements and interactions are comparable, nonlocal interactions assume a large spatial scale of movement but a small temporal scale of movement compared to other processes. Nevertheless, the difference between spatial and temporal scales is not too big to be indistinguishable [23]. The biological mechanism of the nonlocal interactions can be attributed to the effect of mobility [18]. When predators hunt sparse preys, their motion direction should depend not only on the prey population at the same spatial location but also on some weighted average of the population from other locations [18]. Recently researchers have considered spatio-temporal models of interacting populations (such as preypredator and competition models) with nonlocal interactions and investigated the existence of traveling wave, periodic traveling wave and modulated traveling wave [19,20,24–26]. In [27], Tanzy et al. considered in a competition model the effect of asymmetric kernel function of nonlocal interaction on the resulting patterns. In [28,29] the authors studied in a predator–prey model the effect of nonlocal interaction on the oscillatory and chaotic patterns as well as the existence of multiple stationary branches, in particular, the transition from one type of pattern to another. However, all of these studies focused on one-dimensional spatial domain, which raises a question to which degree the achieved results apply to higher dimensional space. More importantly, we wonder whether the emergent spatial patterns whenever feasible are solely due to the nonlocal interactions or combined effects of nonlocal interactions with other factors such as alternative spatially inhomogeneous steady state and/or Hopf bifurcation. Motivated by these questions, we consider the following prey–predator model

⎧ ∂ u1 /∂ t − d1 u1 = u1 (a1 − b11 φδ ∗ u1 − b12 φδ ∗ u2 ), ⎪ ⎪ ⎨∂ u /∂ t − d u = u (−a + b φ ∗ u − b φ ∗ u ), 2 2 2 2 2 21 δ 1 22 δ 2 ⎪∂ u1 /∂ν = ∂ u2 /∂ν = 0, ⎪ ⎩ u1 (0, x ) = u10 (x ), u2 (0, x ) = u20 (x ),

t > 0,

x ∈ ,

t > 0,

x ∈ ,

t > 0,

x ∈ ∂ ,

(1.1)

x ∈ .

Here u1 and u2 are population densities of the prey and predator at time t and space x, respectively.  is a bounded domain in Rn with a smooth boundary ∂ , and ν is the outward unit normal on ∂ , with a zero-flux on the boundaries. di (i = 1, 2 ) are the diffusive rates of prey and predator, respectively. a1 is the intrinsic growth rate of the prey while a2 the background death rate of predator. bii (i = 1, 2 ) are parameters related to nonlocal intraspecific competition, while b12 and b21 are parameters related to nonlocal predation. The nonlocal interactions are given by the presentation of convolution between a spatial kernel function φ δ and population density. The spatial kernel function measures the extent to which a given individual can interact with others at different spatial locations. For instance, the effect of predation on per capital growth rate of prey is given by the term −b12 u1 φδ ∗ u2 , while the contribution of prey on the per capital growth rate of predator is given by the term b21 u2 φ δ ∗u1 . All above appearing parameters are positive, and we consider positive initial conditions, that is, ui0 (x) ≥ 0 for i, j = 1, 2. Moreover, in order to ensure the existence of the classical solution, we assume that ui0 (x) ∈ Cα (). The convolution terms in (1.1) are mathematically represented as

φδ ∗ u i =



Rn

φδ (x − y )ui (t, y )dy,

(1.2)

where φ δ is a spatial kernel function satisfying



Rn

φδ (y )dy = 1.

(1.3)

We also assume that φ δ is an even, nonnegative function with support in the interval [−δ, δ ]n . These conditions on φ δ imply that φ δ approaches a δ -function as δ → 0. Thus, the interactions turn out to be local in this limit. Several different kernel functions have been studied in the literature. The simplest one is the step function [30–32], which is symmetric, unweighted average over a specified interval, φδ (x ) = 21δ , x ∈ (−δ, δ ), whilst φδ (x ) = 0, x ∈ (−δ, δ ). A function with a continuous, rather than abrupt, transition from unity to zero is also possible, but if the transition is sufficiently rapid, the transition region can be neglected, such as the kernel function parabolic segments [18] φδ (x ) = 43δ 2 (δ 2 − x2 ), x ∈ (−δ, δ ), whilst φδ (x ) = 0, x ∈

(−δ, δ ), and elevated triangle [32] φδ (x ) =

αδ −|x| , x ∈ (−δ, δ ), whilst (2α −1 )δ 2

φδ (x ) = 0, x ∈ (−δ, δ ).

The paper is structured as follows. We prove the existence and uniqueness of the solution to the problem (1.1) in the next section. The difficulty is to raise the regularity of the solutions to (1.1) including the spatial convolution. To overcome

C. Tian et al. / Applied Mathematics and Computation 308 (2017) 73–83

75

this difficulty, we construct an iteration sequence that converges to the solution. The ladder argument of the integral representation ensures the regularity of the solutions [33]. In Section 3, we explore the role of the spatially nonlocal interactions in the emergence of spatial patterns, and further analytically derive the conditions for spatial patterns to arise. We analytically show that spatial patterns do not arise in system (1.1) with local interactions but emerge in system (1.1) with nonlocal interactions when the width of the kernel φ δ exceeds a threshold value. Numerical simulations are carried out in Section 4, the results of which well confirm our analytical predictions. The paper is closed with a brief conclusion. 2. Existence and uniqueness of the solution Let  be the closure of , and for any finite T > 0, we set

D T = ( 0 , T ] × ,

DT = [0, T ] × ,

S T = ( 0 , T ] × ∂ .

Denote by Cα (DT ) the set of Hölder continuous functions in DT with exponent α ∈ (0, 1), and by C(DT ) the set of continuous functions in DT . Denote also by C1, 2 (DT ) the set of functions that are once continuously differentiable in t and twice continuously differentiable in x. The above spaces for vector-valued functions are denoted by C α (DT ), C (DT ) and C 1,2 (DT ), respectively. For simplicity, throughout this paper, we denote

f1 (u1 , u2 ) ≡ u1 (a1 − b11 φδ ∗ u1 − b12 φδ ∗ u2 ),

f2 (u1 , u2 ) ≡ u2 (−a2 + b21 φδ ∗ u1 − b22 φδ ∗ u2 ). ˜ = (u˜1 , u˜2 ) and u ˆ = (uˆ1 , uˆ2 ) ∈ C (DT ) ∩ C 1,2 (DT ) is called coupled upper and Definition 2.1. A pair of nonnegative functions u ˆ and if ˜ ≥u lower solutions of (1.1) if u

∂ u˜1 /∂ t − d1 u˜1 ≥ u˜1 (a1 − b11 φδ ∗ u˜1 − b12 φδ ∗ uˆ2 ) ∂ u˜2 /∂ t − d2 u˜2 ≥ u˜2 (−a2 + b21 φδ ∗ u˜1 − b22 φδ ∗ u˜2 ) ∂ uˆ1 /∂ t − d1 uˆ1 ≤ uˆ1 (a1 − b11 φδ ∗ uˆ1 − b12 φδ ∗ u˜2 ) ∂ uˆ2 /∂ t − d2 uˆ2 ≤ uˆ2 (−a2 + b21 φδ ∗ uˆ1 − b22 φδ ∗ uˆ2 ) ∂ uˆi /∂ν ≤ 0 ≤ ∂ u˜i /∂ν for i = 1, 2, on ST , u˜i (0, x ) ≥ ui0 (x ), uˆi (0, x ) ≤ ui0 (x ) for i = 1, 2,

in DT , in DT , in DT ,

(2.1)

in DT , in .

ˆ means for i = 1, 2, u˜i ≥ uˆi for (t, x) ∈ DT . ˜ ≥u Here u ˜ and u ˆ , we set For a given pair of coupled upper and lower solutions u

uˆi , u˜i  = {ui ∈ C (DT ) : uˆi ≤ u ≤ u˜i }, uˆ , u˜  = {u ∈ C (DT ) : uˆ ≤ u ≤ u˜ }. ˆ, u ˜  if there exists a constant ki > 0, for any u, v ∈ u ˆ, u ˜ , Definition 2.2. fi is called Lipschitz continuous with respect to u such that

| f i ( u 1 , u 2 , φδ ∗ u 1 , φδ ∗ u 2 ) − f i ( v 1 , v 2 , φδ ∗ v 1 , φδ ∗ v 2 ) | ≤ k i ( | u 1 − v 1 | + | u 2 − v 2 | + φδ ∗ | u 1 − v 1 | + φδ ∗ | u 2 − v 2 | ) .

(2.2)

ˆ, u ˜ , we call f = ( f1 , f2 ) is Lipschitz continuous with Furthermore, if f1 and f2 are Lipschitz continuous with respect to u ˆ, u ˜ . respect to u ˜ and u ˆ are bounded, in terms of (1.3), direct calculations show that there exist constants k1 and k2 , where If u

k1 = max{a1 + b11 (||u1 ||0 + ||v1 ||0 ) + b12 ||u2 ||0 , b12 ||u1 ||0 }, k2 = max{a2 + b22 (||u2 ||0 + ||v2 ||0 ) + b21 ||u1 ||0 , b21 ||u2 ||0 }, ˆ, u ˜ . and ||ui ||0 = supD |ui |, such that f is Lipschitz continuous in u T In addition, for the sake of convenience, we define the following operators

∂ u1 /∂ t − d1 u1 + k1 u1 , L2 u2 = ∂ u2 /∂ t − d2 u2 + k2 u2 , F1 (u1 , u2 ) = k1 u1 + f1 (u1 , u2 ), F2 (u1 , u2 ) = k2 u2 + f2 (u1 , u2 ). L1 u1 =

Then the system (1.1) can be reformulated as follows

Li ui = Fi (u1 , u2 ) for i = 1, 2, ∂ u1 /∂ν = ∂ u2 /∂ν = 0, u1 (x, 0 ) = u10 (x ), u2 (x, 0 ) = u20 (x ),

in DT , on ST , in .

(2.3)

Now we are in position to show that the system (2.3) has a unique global solution. To this aim we construct a sequence

{u(m) } ≡ {u1(m) , u2(m) } according to the following iteration process Li ui(m ) = Fi (u1(m−1) , u2(m−1) ) for i = 1, 2, ∂ u1(m) /∂ν = ∂ u2(m) /∂ν = 0, u1(m ) (x, 0 ) = u10 (x ), u2(m ) (x, 0 ) = u20 (x )

in DT , on ST , in .

76

C. Tian et al. / Applied Mathematics and Computation 308 (2017) 73–83

with u0 ∈ C α (DT ) ∩ C (DT ). To show the convergence of the sequence {u(m) }, we perform the transformation wi = e−γ t ui (i = 1, 2 ) for some positive constant γ > 0 to ensure contraction property. The system (2.3) is equivalent to the following system

Li wi + γ wi = Hi (w1 , w2 ) for i = 1, 2, ∂ w1 /∂ν = ∂ w2 /∂ν = 0, w1 (x, 0 ) = w10 (x ), w2 (x, 0 ) = w20 (x ),

in DT , on ST , in ,

(2.4)

where

w i 0 ( x ) = e −γ t u i 0 ( x )

for i = 1, 2,

H1 (w1 , w2 ) = k1 w1 + w1 (a1 − b11 eγ t φδ ∗ w1 − b12 eγ t φδ ∗ w2 ),

H2 (w1 , w2 ) = k2 w2 + w2 (−a2 + b21 eγ t φδ ∗ w1 − b22 eγ t φδ ∗ w2 ). According to (2.4), we can construct sequences (m )

(m )

(m−1 )

{w(m) }

(m−1 )

Li wi + γ wi = Hi (w1 , w2 ) ∂ w1(m) /∂ν = ∂ w2(m) /∂ν = 0, w1(m ) (x, 0 ) = w10 (x ), w2(m ) (x, 0 ) = w20 (x ),

(2.5)

via the following iteration process

for i = 1, 2, in DT , on ST , in .

In term of the integral representation theory for linear parabolic boundary-value problems, the sequence {w(m) } can be expressed as

wi(m ) (t, x ) =



t 0

dτ 

+

0

t

 

dτ

i (t, x; τ , ξ )(Hi (w(m−1) ))(τ , ξ )dξ  ∂

i (t, x; τ , ξ )ψi(m−1) (τ , ξ )dξ +

 

i (t, x; 0, ξ )wi0 (ξ )dξ ,

(2.6)

where  i is the fundamental solution of the linear parabolic operator Li + γi , and ψ i is the single-layer potential [33]. We now show that the sequence {w(m) } converges in C (DT ) to a unique solution of the associated integral Eq. (2.6). In Lemma 2.1, we give some properties of the function Hi (w), defined in (2.5), in the space X ≡ X1 × X2 , where

Xi ≡ {wi ∈ C α (DT ) ∩ C (DT ) : wi (0, x ) = wi0 (x ) in }, for i = 1, 2. Lemma 2.1. If w, w ∈ X , then Hi (w ) ∈ C α (DT ) ∩ C (DT ) and

|Hi (w ) − Hi (w )| ≤ 3ki (|w1 − w 1 | + |w2 − w 2 | ) for i = 1, 2.

(2.7)

Proof. We first show Hi (w ) ∈ C α (DT ) ∩ C (DT ). Without loss of generality we only need to treat the case of spatial convolution because we will see later that the plus and the multiplication of Eq. (2.5) do not change the Hölder continuous property of the functions. By (1.2),

φδ ∗ wi (t, x ) =



Rn

φδ (x − y )wi (t, y )dy =



Rn

φδ (y )wi (t, x − y )dy.

Let (t, x), (t , x ) ∈ DT . By the property of φ δ and the Hölder continuous property of wi , we have

|φδ ∗ wi (t, x ) − φδ ∗ wi (t , x )| = ≤



Rn



Rn

φδ (y )|wi (t, x − y ) − wi (t , x − y )|dy ci (|t − t |α + |x − x |α )φδ (y )dy

≤ ci (|t − t |α + |x − x |α ) for some constant ci . This shows that φδ ∗ wi (t, x ) is Hölder continuous in DT for i = 1, 2. It follows from the Hölder continuity of Hi (w1 , w2 , φδ ∗ w1 , φδ ∗ w2 ) in (t, x) ∈ DT and the Lipschitz condition (2.2) that Hi (w) is Hölder continuous in DT . To show the relation (2.7) we observe from (2.2) that

|Hi (w ) − Hi (w )| = |(ki wi + e−γ t fi (eγ t w1 , eγ t w2 , eγ t φδ ∗ w1 , eγ t φδ ∗ w2 )) −(ki w i + e−γ t fi (eγ t w 1 , eγ t w 2 , eγ t φδ ∗ w 1 , eγ t φδ ∗ w 2 ))| ≤ ki (|wi − w i | ) + e−γ t ki (|w1 − w 1 | + |w2 − w 2 | ) + ki (|w1 − w 1 | + |w2 − w 2 | ) ≤ 3ki (|w1 − w 1 | + |w2 − w 2 | ). This proves the lemma.



Based on the global Lipschitz property (2.7), we have the following existence result for the system (2.3).

C. Tian et al. / Applied Mathematics and Computation 308 (2017) 73–83

77

˜, u ˆ ) be a pair of coupled upper and lower solutions of (1.1). Then the system (1.1) has a unique solution Theorem 2.1. Let (u ˆ, u ˜ . Moreover, for any u(0 ) ∈ C α (DT ) ∩ C (DT ) with u(0 ) (0, x ) = (u10 (x ), u20 (x )) in , the sequence {u(m) } u∗ (t, x) and u∗ ∈ u obtained from (2.4) converges to u∗ as m → ∞. Proof. We prove the theorem using the contraction mapping theorem in C (DT ) for the transformed Eqs. (2.4). For each i = 1, 2, we define operators Ai : D(Ai ) → R(Ai ) and Hi : X → C α (DT ) ∩ C (DT ) by

Ai (wi ) ≡ Li wi + γ wi (wi ∈ D(Ai )), Hi (w ) ≡ Hi (w1 , w2 )

( w ∈ X ),

(2.8)

where D(Ai ) is the domain of Ai given by

D(Ai ) = {wi ∈ C 1,2 (DT ) ∩ C (DT ) :

∂ wi /∂ ν = 0 on ST , wi (0, x ) = wi0 (x ) in }.

R(Ai ) is the range of Ai , and Hi (w1 , w2 ) is given by (2.5). It follows from Lemma 2.1 that Hi maps X into C α (DT ) ∩ C (DT ). In terms of the operators Ai and Hi the iteration process (2.6) can be written as

Ai wi(m ) = Hi (w(m ) )

(wi(m−1) ∈ D(Ai )) for i = 1, 2,

and in vector form it becomes

Aw(m ) = H (w(m−1) )

(w(m) ∈ D(A )).

(2.9)

From the standard parabolic theorem the inverse operator

||A w − A −1

−1

w ||0 ≤

( γ + k3 )

−1

||w − w ||

A−1

exists and possesses the property

α 0 for w, w ∈ C (DT ) ∩ C (DT )

(2.10)

where k3 = min{k1 , k2 }. This implies that Eq. (2.9) is equivalent to

w(m ) = A−1 H (w(m−1) )

(w(m) ∈ D(A )),

(2.11)

which can be considered as a compact form for the integral representation (2.6) in the space C α (DT ) ∩ C (DT ). In terms of (2.7), there exists a constant k, independent of γ , such that

||H(w ) − H(w )||0 ≤ k||w − w ||0 for w, w ∈ X .

(2.12)

Combining (2.10) and (2.12), we have

||A−1 H(w ) − A−1 H(w )||0 ≤ k(γ + k3 )−1 ||w − w ||0 for w, w ∈ X . ||A−1 H(w )

− A−1 H (w )||

(2.13)

−1 0 ≤ ||w − w ||0 for w, w ∈ X . Thus, the operator A H possesses a

By choosing γ > k, we have contraction property in X . This ensures that the sequence {w(m) } converges in C (DT ) to some w∗ ∈ C (DT ). By the equivalence between (2.11) and (2.6) the sequence {wi(m ) } given by (2.6) converges in C (DT ) to w∗i for i = 1, 2. To show that w∗ is the solution of (2.4), we need to raise the regularity of w∗ . Letting m → ∞ in (2.6), the analogous argument used in [33] shows that the single-layer potential ψim converges to some continuous function ψi∗ , and w∗i satisfies the integral equation

w∗i (t, x ) =



t 0

+

dτ



t 0

 

dτ

i (t, x; τ , ξ )(Hi (w∗ ))(τ , ξ )dξ

 ∂

i (t, x; τ , ξ )ψi∗ (τ , ξ )dξ +

 

i (t, x; 0, ξ )wi0 (ξ )dξ .

Since w∗i is continuous in DT , a ladder argument as that in [33] shows that w∗i ∈ C 1,2 (DT ) ∩ C (DT ). This implies that w∗ is the unique solution of (2.4). Since u(m ) = eγ t w(m ) , the sequence u(m) governed by (2.4) converges to a unique solution u∗ = eγ t w∗ to the Eq. (2.3). By the equivalence between (2.3) and (1.1), u∗ is the unique solution of the system (1.1). This proves the theorem.  Theorem 2.1 demonstrates that to show the existence and uniqueness of the solution to the system (1.1), we only need ˜ and u ˆ , which satisfies (2.2). If we choose u ˜ and u ˆ being constant to seek a pair of coupled upper and lower solutions u vectors c˜ and cˆ , these constants need to satisfy

0 ≥ c˜1 (a1 − b11 c˜1 − b12 cˆ2 ), 0 ≥ c˜2 (−a2 + b21 c˜1 − b22 c˜2 ), 0 ≤ cˆ1 (a1 − b11 cˆ1 − b12 c˜2 ), 0 ≤ cˆ2 (−a2 + b21 cˆ1 − b22 cˆ2 ), c˜i ≥ sup ui (0, x ), 

cˆi ≤ inf ui (0, x ) for i = 1, 2. 

Thus, we have the existence result for the system (1.1) and the bounded estimate.

78

C. Tian et al. / Applied Mathematics and Computation 308 (2017) 73–83

Theorem 2.2. If the initial function ui0 (x ) ∈ C α () ∩ C () and ui0 (x) ≥ 0 for i = 1, 2, then the system (1.1) has a unique global solution u∗ (t, x) for (t, x) ∈ (0, ∞) × , and for (t, x) ∈ (0, ∞) × ,

0 ≤ u∗1 (t, x ) ≤ M1 , 0 ≤ u∗2 (t, x ) ≤ M2 ,

(2.14)

where

M1 = max{sup u10 (x ), a1 /b11 }, 

M2 = max{sup u20 (x ), (b21 M1 − a2 )/b22 }. 

Proof. From the above argument, we only need to choose cˆ1 = cˆ2 = 0, c˜1 = max{sup u10 (x ), a1 /b11 }, c˜2 = max{sup u20 (x ), (b21 c˜1 − a2 )/b22 }. Thus, the inequalities in (2.14) hold, which implies that c˜ and cˆ are a pair of coupled upper and lower solutions of (1.1). Applying the Theorem 2.1, the proof completes.  Remark 2.1. The result obtained in this section applies to the bounded domain in Rn , but to derive the specific conditions driving the emergence of spatial patterns by nonlocal interactions, we confine the spatial domain in R2 in the coming section. Furthermore, from a biological point of view, two dimensional space coincides with the real-world ecological community. 3. Nonlocal interaction driven pattern formation In this section, we aim to derive the conditions on the spatial convolution such that the positive spatially uniform equilibrium, which is actually unique when feasible, is locally stable in the case of δ = 0 but unstable for some δ > 0. This phenomenon is called nonlocal interaction driven instability.  Note that the spatial averaging kernel satisfies  φδ (y )dy = 1, which means that the nonlocal coupling has no effect on spatially uniform solutions. Thus, the uniform equilibria of the system (1.1) satisfy the algebraic system

u1 (a1 − b11 u1 − b12 u2 ) = 0, u2 (−a2 + b21 u1 − b22 u2 ) = 0.

(3.1)

It is clear that system (3.1) admits a trivial equilibrium (0, 0), corresponding to extinction of both species. Since the net reproductive rates a1 and a2 are both positive, the extinction equilibrium (0, 0) is unstable. Furthermore, there exists in system (3.1) a semi-trivial equilibrium, i.e., (a1 /b11 , 0), corresponding to extinction of the predator. In this situation, the system (1.1) is reduced to a scalar logistic population model with nonlocal coupling. For the scalar logistic population model, it was shown that the nonlocal interactions destabilize the uniform equilibrium, which is stable for the system with local iterations [30,31]. We are interested in the positive equilibrium, corresponding to the coexistence of prey and predator. In fact, there is a unique positive equilibrium in system (3.1)



a1 b22 − a2 b12 a2 b11 + a1 b21 , b11 b22 + b12 b21 b11 b22 + b12 b21

(u , v ) = ∗

∗



(3.2)

under the following parametric condition

b11 a1 < . b21 a2

(3.3)

Furthermore, straightforward computation shows that condition (3.3) is necessary and sufficient for the coexistence equilibrium (u∗ , v∗ ) being stable when interactions are local, i.e., when δ → 0. Next we consider the case where (u∗ , v∗ ) is destabilized by the nonlocal interactions (1.2). Taking into account the fact that the

spatial kernel functions are supported in the interval [−δ, δ ]n , we obtain that φˆ is only dependent on β = kδ, here k ≡ k2 + k2 is the wave number, k = (kx , ky ) δ

x

y

is the wave-vector. In this paper, we take φˆ δ as a function only of β . Theorem 3.1. If (3.3) holds and φˆ δ (β ) > 0 for all β , then the positive equilibrium (u∗ , v∗ ) is locally asymptotically stable. Proof. We carry out the linear stability analysis of (1.1) around (u∗ , v∗ ). To this end, we set the perturbations u = u1 − u∗ , v = u2 − v∗ , and plug these perturbations into (1.1). Retaining only the linear terms in u and v, we have the linear integrodifferential system

∂ u/∂ t − d1 u = −b11 u∗ φδ ∗ u − b12 u∗ φδ ∗ v, ∂v/∂ t − d2 v = b21 v∗ φδ ∗ u − b22 v∗ φδ ∗ v. Applying the Fourier transform in x yields the following matrix form

d dt

 uˆ vˆ



=

−k2 d1 − b11 u∗ φˆ δ (k ) b21 v∗ φˆ δ (k )



−b12 u∗ φˆ δ (k ) 2 −k d2 − b22 v∗ φˆ δ (k )

uˆ vˆ

≡M



uˆ , vˆ

(3.4)

C. Tian et al. / Applied Mathematics and Computation 308 (2017) 73–83

79

where k is the wavenumber, and uˆ, vˆ and φˆ δ denote the Fourier transform of u, v and φ δ , respectively. In view of the spatial kernel function φ δ being even, φˆ δ (k ) is real for all k. Moreover, (1.3) implies that φˆ δ (0 ) = 1. In term of β = kδ, the system (3.4) shows that (u∗ , v∗ ) is locally asymptotically stable if and only if the eigenvalues of matrix M have negative real parts, which is equivalent to trace(M) < 0 and det(M ) > 0, i.e.,

− ( d1 + d2 ) d1 d2

β2 − (b11 u∗ + b22 v∗ )φˆ δ (β ) < 0, δ2

(3.5)

β4 β2 + (b11 d2 u∗ + b22 d1 v∗ ) 2 φˆ δ (β ) + (b11 b22 + b12 b21 )u∗ v∗ φˆ δ2 (β ) > 0. 4 δ δ

(3.6)

Since φˆ δ (β ) > 0 for all β according to our assumption, the inequalities (3.5) and (3.6) are satisfied, which completes the proof.  Remark 3.1. When δ → 0, φˆ δ (0 ) = 1, which implies that φˆ δ (β ) > 0 for all β is naturally satisfied. Then from Theorem 3.1 the positive equilibrium is locally asymptotically stable. Theorem 3.1 covers the stability result of classical prey–predator model with local interactions. Theorem 3.2. Assume that (3.3) and

(b11 d2 u∗ − b22 d1 v∗ )2 > 4b12 b21 d1 d2 u∗ v∗

(3.7)

hold. We take δ as a bifurcation parameter. If there exists a critical point δ 0 such that

1

δ02

≤ min S(β ),

(3.8)

β >0

where

S (β ) =

−φˆ δ (β ) (b11 d2 u∗ + b22 d1 v∗ − (b11 d2 u∗ − b22 d1 v∗ )2 − 4b12 b21 d1 d2 u∗ v∗ ), 2 2d1 d2 β

(3.9)

then the positive equilibrium (u∗ , v∗ ) is unstable. Furthermore, if d1 = d2 , the system (1.1) undergoes Turing bifurcation at δ 0 . Proof. Notice that the violations of relations (3.5) and (3.6), respectively, correspond to the onset of temporal Hopf bifurcation, and spatial Turing bifurcation. To show spatial Turing bifurcation, we only need to verify that (3.6) fails while 2 (3.5) holds. Considering (3.6) as a binomial expression of β and φˆ (β ), we have δ2

δ

β β + (b11 d2 u∗ + b22 d1 v∗ ) 2 φˆ δ (β ) + (b11 b22 + b12 b21 )u∗ v∗ φˆ δ2 (β ) δ4 δ  2 1 β 2 b11 d2 u∗ + b22 d1 v∗ ˆ = d1 d2 2 + φδ ( β ) d1 d2 2 δ 4

2

d1 d2

−

φˆ δ2 (β )  4d1 d2

(b11 d2 u∗ + b22 d1 v∗ )2 − 4d1 d2 (b11 b22 + b12 b21 )u∗ v∗ .

(3.10)

It is easy to see the condition (3.7) is necessary for (3.10) being smaller than 0. If additionally

1

δ2

≤

−φˆ δ (β )  ∗ ∗ b d u + b d v − (b11 d2 u∗ − b22 d1 v∗ )2 − 4b12 b21 d1 d2 u∗ v∗ 11 2 22 1 2d1 d2 β 2

(3.11)

is satisfied, which is the condition (3.8), then (u∗ , v∗ ) is unstable, which completes the proof. To show that the system (1.1) undergoes spatial Turing bifurcation at δ 0 , we need to verify spatial Turing bifurcation to occur prior to the temporal Hopf bifurcation as δ increases to δ 0 . From the above argument, we only need to show that if (3.5) fails then (3.6) must have already failed as δ increases. When (3.5) fails and additionally d1 = d2 , we have

1

δ2

=−

b11 u∗ + b22 v∗ . 2d1 β 2

(3.12)

Plugging (3.12) into (3.10), we see that the square term in (3.10) vanishes, and

β4 β2 + (b11 d2 u∗ + b22 d1 v∗ ) 2 φˆ δ (β ) + (b11 b22 + b12 b21 )u∗ v∗ φˆ δ2 (β ) 4 δ δ φˆ δ2 (β )  = − (b11 d2 u∗ + b22 d1 v∗ )2 − 4d1 d2 (b11 b22 + b12 b21 )u∗ v∗ ≤ 0.

d1 d2

4d1 d2

Thus (3.6) does not hold, which completes the proof.

(3.13)



Remark 3.2. From Theorem 3.2, φˆ δ < 0 is necessary for spatial Turing instability. This case can happen if φ δ has a discontinuity, e.g., the spatial kernel is a step function (see next section).

80

C. Tian et al. / Applied Mathematics and Computation 308 (2017) 73–83

0.02 0.01 0 −0.01 S(β) 2

1/(δ ) with δ=10

−0.02

2

1/(δ ) with δ=15 −0.03 −0.04 −0.05 2

4

6

8

10

β

12

14

16

18

20

Fig. 1. S(β ) and 1/δ 2 for different δ . Parameters are given in (4.1). The positive equilibrium (4.2) is stable for δ = 10 but unstable for δ = 15 at β = 4 corresponding to the norm of wavenumber k = 3.75.

4. Numerical results In this section, we carry out numerical simulations to demonstrate our analytical findings. The following set of parameters is employed

a1 = 1, a2 = 0.2, b11 = 0.1, b12 = 0.1, b21 = 0.05, b22 = 0.1, d1 = 0.1, d2 = 0.1, Lx = 100, Ly = 100,

(4.1)

which gives rise to a unique positive equilibrium

( u∗ , v∗ ) = ( 8, 2 ).

(4.2)

It is easy to check that conditions (3.3) and (3.7) are satisfied under this set of parameters, which means that the positive equilibrium is locally stable with spatially local interactions but can be unstable with spatially nonlocal interactions. For nonlocal interaction driven pattern formation, Theorem 3.2 shows that the spatial kernel function must be φˆ δ < 0 for some δ . To fulfil this condition, we consider the following step function

φδ (x, y ) =

1 , 4δ 2

0,

(x, y ) ∈ (−δ, δ ) × (−δ, δ ), (x, y ) ∈ (−δ, δ ) × (−δ, δ ),

whose Fourier transform is

1

φˆ δ (k ) = √

2π

sin(kδ ) , kδ

where k is the spatial wavenumber. From conditions (3.8) and (3.9), we obtain a spatial bifurcation curve, defined by S(β ), which is graphically illustrated in Fig. 1. The figure shows that the positive equilibrium (4.2) is unstable for those δ such that the curve S(β ) intersects with the line 1/δ 2 . For given parameter values in (4.1), we see that the minimum δ required for the emergence of Turing instability is approximately δ = 15. Below we present details for numerical simulation of spatial patterns. We consider an initial data as a random perturbation (white noise) around the positive equilibrium (u∗ , v∗ ) in , with a variance lower than the amplitude of the positive equilibrium. More precisely,

u10 (x ) = u∗ + η1 (x ),

u20 (x ) = v∗ + η2 (x ),

where η1 ∈ [−1.5, 1.5] and η2 ∈ [−2, 2]. Moreover, we consider a rectangular domain  = [−L/2, L/2] × [−L/2, L/2] ⊂ R2 , and resolve system (1.1) on a grid with L × L sites with a space step of x = y = 1. The wave-vector for this two-dimensional domain is thereby

k = 2π (m/L, n/L ), k = 2π



L 2

L 2

L 2

(m/L )2 + (n/L )2 , m, n = − , − + 1, . . . , .

C. Tian et al. / Applied Mathematics and Computation 308 (2017) 73–83

81

25

20

u

1

15

10

5

0

4

6

8

10

12

δ

14

16

18

20

22

Fig. 2. Bifurcation diagram of u1 with respect to δ . The amount of vertical axis is the spatially averaged prey density.

For the evolution in time, we apply the first order backward Euler time advancing scheme with a time step t = 0.005. Discretizing the Laplacian in the grid with lattice sites denoted by (i, j) yields the following the nine-point formula

u|(i, j ) =

1 [4al (i, j )u(i − 1, j ) + 4ar (i, j )u(i + 1, j ) + 4ad (i, j )u(i, j − 1 ) 6x2 + 4au (i, j )u(i, j + 1 ) + al (i, j )u(i − 1, j + 1 ) + au (i, j )u(i + 1, j + 1 ) + ad (i, j )u(i − 1, j − 1 ) + ar (i, j )u(i + 1, j − 1 ) − 20u(i, j )],

where the matrix elements of al , ar , ad , au are unity except at the boundary. When (i, j) is at the left boundary, that is i = 0. We define al (i, j )u(i − 1, j ) ≡ u(i + 1, j ), which guarantees zero-flux of reactants in the left boundary. Similarly we define ar (i, j), ad (i, j), au (i, j) such that the boundary is zero-flux. The nine-point formula for the Laplacian has an error of O(x4 ). Fig. 2 shows the spatially averaged density of the prey species before and after the onset of Turing bifurcation. Results are qualitatively similar for u2 , and hence omitted. For δ less than 15, i.e., the Turing instability does not occur, we see that the density of u1 is spatially homogeneous. For δ greater than 15, i.e., spatial Turing bifurcation happens, we see that the density of u1 is spatially inhomogeneous. In our numerical simulation, we set δ = 15 and L = 100, such that the width of the kernel does not encounter the boundary of the domain. Sampled emergent spatial patterns are presented in Fig. 3. We observe only regular spot patterns with a spot width increasing with δ . We do not observe other spatial patterns such as strip pattern, although different initial values are tried. In fact, the selection of strip pattern or spot pattern depends upon the non-linearities of the reaction kinetics. Specifically, it has been shown that the presence of quadratic nonlinearities in the reaction kinetics leads to spot pattern, but the absence of quadratic terms leads to strip pattern [34]. The Taylor series expansion of the system (1.1) is

∂ u1 /∂ t − d1 u1 = (−b11 u∗ u1 − b12 u∗ u2 ) − (b11 u1 φδ ∗ u1 + b12 u1 φδ ∗ u2 ), ∂ u2 /∂ t − d2 u2 = (b21 v∗ u1 − b22 v∗ u2 ) + (b21 u2 φδ ∗ u1 − b22 u2 φδ ∗ u2 ).

(4.3)

Noticing that the reaction kinetics of (4.3) has only quadratic nonlinearities. Thus, in view of the theory of pattern selection [34], we can only observe spot patterns in our model. Moreover, we compare our numerical experiment with the pattern selection for competing population model proposed by Segal et al. [20]. They considered one dimensional space and found by numerical simulations that the emergent spatial pattern is a cosine-like curve and for most observed patterns, population density at the valleys of the curve is near extinction. Our numerical experiments exhibited a qualitatively similar outcome in two dimensional space that is, spot pattern, but we also observed a quantitative difference in the population density which is far from extinction even at the valleys (Fig. 3).

82

C. Tian et al. / Applied Mathematics and Computation 308 (2017) 73–83

δ = 16

δ = 18 10

10

9

9

8

8

7

7

δ = 22

δ = 20 10

10

9

9

8

8

7

7

Fig. 3. Stationary spatial patterns for varied δ under the default parameters in (4.1).

5. Conclusion In this paper, we studied pattern formation in a prey–predator model with nonlocal intraspecific competition and nonlocal predation. We proved that our model (1.1) admits a unique global bounded solution, which is locally stable. This means that in case of local interactions there are no spatial patterns including Turing patterns. However introducing nonlocal interactions breaks up the homogeneous distribution of population densities, resulting in emergence of spot spatial patterns (Fig. 3), as long as the width of nonlocal interaction kernel goes beyond a threshold value (Fig. 1). Our results unravel a new mechanism for the emergence of spatial patterns, that is, nonlocal interaction driven spatial patterns. Generally in spatially homogeneous environments, spatial patterns arise either due to stable kinetic dynamic with Turing instability or due to unstable kinetic dynamics (via temporal Hopf bifurcation) with spatially inhomogeneous perturbations. However, we found that, on the one hand, when interactions are local (i.e., δ = 0) the unique equilibrium in our model (1.1) is locally stable and there is no Turing instability, and on the other hand, there exists pattern formation with d1 = d2 = 1 when nonlocal interactions are introduced. Thus, the emergence of spatial patterns cannot be attributed to Turing instability since d1 = d2 because Turing instability requires sufficiently large difference in terms of the diffusion rate between predator and prey [15]. Nor is it due to the temporal Hopf bifurcation because the equilibrium is locally stable. Taken together, it is because of nonlocal interactions induced spatial Turing bifurcation that triggers the formation of spatial patterns. Since we only observe spot patterns in this paper, an interesting extension is to explore how nonlocal interactions affect spatial temporal dynamics in a system with diverse spatial patterns such as shown in [3,4].

C. Tian et al. / Applied Mathematics and Computation 308 (2017) 73–83

83

Acknowledgments Canrong Tian is grateful for the financial support from the Oversea Study Programme of Talent Young Teachers of Jiangsu Province. Zhi Ling and Lai Zhang are supported by the PRC Grant NSFC 11571301. Zhi Ling is supported by the NSF of Jiangsu Province BK20151305 also. Lai Zhang further acknowledges the financial support from the Swedish Strategic Research Programme eSSENCE. 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]

A. Turing, The chemical basis of morphogenesis, Philos. Trans. R. Soc. B 237 (1952) 37–72. S.A. Levin, L.A. Segel, An hypothesis to explain the origin of planktonic patchiness, Nature 259 (1976) 659. M. Banerjee, S. Petrovskii, Self-organised spatial patterns and chaos in a ratio-dependent predator–prey system, Theor. Ecol. 4 (2011) 37–53. 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. Y. Lou, W.M. Ni, Diffusion, self-diffusion and cross-diffusion, J. Differ. Equ. 131 (1996) 79–131. C. Tian, Turing patterns created by cross-diffusion for a Holling II and Leslie–Gower type three species food chain model, J. Math. Chem. 49 (2011) 1128–1150. Z. Ling, L. Zhang, Z.G. Lin, Turing pattern formation in a predator–prey system with cross diffusion, Appl. Math. Model. 38 (2014) 5022–5032. A.B. Rovinsky, M. Menzinger, Self-organization induced by the differential-flow of activator and inhibitor, Phys. Rev. Lett. 70 (1993) 778–781. R.A. Satnoianu, J.H. Merkin, S.K. Scott, Pattern formation in a differential-flow reactor model, Chem. Eng. Sci. 55 (20 0 0) 461–469. R.A. Satnoianu, P.K. Maini, M. Menzinger, Parameter space analysis, pattern sensitivity and model comparison for turing and stationary flow-distributed waves (FDS), Phys. D 160 (2001) 79–102. C. Tian, L. Zhang, Delay-driven irregular spatiaotemporal patterns in a plankton system, Phys. Rev. E 88 (2013) 012713. S. Sen, P. Ghosh, S.S. Riaz, D.S. Ray, Time-delay-induced instabilities in reaction–diffusion systems, Phys. Rev. E 80 (2009) 046212. M. Banerjee, L. Zhang, Influence of discrete delay on pattern formation in a ratio-dependent prey–predator model, Chaos, Solitons Fractals 67 (2014) 73–81. L. Zhang, U.H. Thygensen, M. Banerjee, Size-dependent diffusion promotes the emergence of spatiotemporal patterns, Phys. Rev. E. 90 (2014) 012904. S.V. Petrovskii, H. Malchow, A minimal model of pattern formation in a prey–predator system, Math. Comput. Model. 29 (1999) 49–63. J. Ge, K. Kim, Z. Lin, H. Zhu, A SIS reaction–diffusion-advection model in a low-risk and high-risk domain, J. Differ. Equ. 259 (2015) 5486–5509. C. Lei, Z. Lin, Q. Zhang, The spreading front of invasive species in favorable habitat or unfavorable habitat, J. Differ. Equ. 257 (2014) 145–166. N.F. Britton, Aggregation and the competitive exclusion principle, J. Theor. Biol. 136 (1989) 57–66. N.F. Britton, Spatial structures and periodic travelling waves in an integrodifferential reaction–diffusion population model, SIAM J. Appl. Math. 50 (1990) 1663–1688. B.L. Segal, V.A. Volpert, A. Bayliss, Pattern formation in a model of competing populations with nonlocal interactions, Phys. D 253 (2013) 12–22. E. Ben-Jacob, I. Cohen, H. Levine, Cooperative self-organization of microorganisms, Adv. Phys. 49 (20 0 0) 395–554. M.G. Clerc, D. Escaff, V.M. Kenkre, Analytical studies of fronts, colonies, and patterns: combination of the Allee effect and nonlocal competition interactions, Phys. Rev. E 82 (2010) 036210. C.T. Lee, M.F. Hoopes, J. Diehl, W. Gilliland, G. Huxel, E.V. Leaver, K. McCann, J. Umbanhowar, A. Mogilner, Nonlocal concepts and models in biology, J. Theor. Biol. 201 (2001) 201–219. N. Apreutesei, N. Bessonov, V. Volpert, V. Vougalter, Spatial structures and generalized travelling waves for an integro-differential equation, Discr. Contin. Dyn. Syst. B 13 (2010) 537–557. M. Banerjee, V. Volpert, Prey–predator model with a nonlocal consumption of prey, Chaos 26 (2016) 083120. J.A. Sherratt, Periodic traveling waves in integro-differential equations for nonlocal dispersal, SIAM J. Appl. Dyn. Syst. 13 (2014) 1517–1541. M.C. Tanzy, V.A. Volpert, A. Bayliss, M.E. Nehrkorn, Stability and pattern formation for competing populations with asymmetric nonlocal coupling, Math. Biosci. 246 (2013) 14–26. M. Banerjee, V. Volpert, Spatio-temporal pattern formation in Rosenzweig–Macarthur model: effect of nonlocal interactions, Ecol. Complex. (2017), doi:10.1016/j.ecocom.2016.12.002. (in press) M. Banerjee, L. Zhang, Stabilizing role of nonlocal interaction on spatio-temporal pattern formation, Math. Model. Nat. Phenom. 11 (2016) 103–118. S. Genieys, V. Volpert, P. Auger, Pattern and waves for a model in population dynamics with nonlocal consumption of resources, Math. Model. Nat. Phenom. 1 (2006) 63–80. M.A. Fuentes, M.N. Kuperman, V.M. Kenkre, Nonlocal interaction effects on pattern formation in population dynamics, Phys. Rev. Lett. 91 (2003) 158104. M.A. Fuentes, M.N. Kuperman, V.M. Kenkre, Analytical considerations in the study of spatial patterns arising from nonlocal interaction effects, J. Phys. Chem. B 108 (2004) 10505–10508. C.V. Pao, Nonlinear Parabolic and Elliptic Equation, Plenum, New York, 1992. B. Ermentrout, Stripes or spots? Non-linear effects in bifurcation of reaction–diffusion equations on the square, Proc. R. Soc. Lond. A 434 (1991) 413–417.