Delay induced spatiotemporal patterns in a diffusive intraguild predation model with Beddington-DeAngelis functional response.

A diffusive intraguild predation model with delay and Beddington-DeAngelis functional response is considered. Dynamics including stability and Hopf bifurcation near the spatially homogeneous steady states are investigated in detail. Further, it is numerically demonstrated that delay can trigger the emergence of irregular spatial patterns including chaos. The impacts of diffusion and functional response on the model's dynamics are also numerically explored.


1.
Introduction. Competition and predation are two fundamental ecological relationships among species and have been widely studied [1]. Recently, it is been recognized that intraguild predation (IGP), which is a combination of competition and predation, has significant impacts on the distribution, abundance, persistence and evolution of the species involved [2]. As a result, growing attention has been paid to IGP models [3,4,5,6,7,8,9].
The general framework of IGP described below was established by Holt and Polis [5]      where R(t), N (t), P (t) represent the densities of basal resource, IG prey and IG predator, respectively. The quantities ρ 2 (R, N, P )R and ρ 3 (R, N, P )N are functional responses of the IG predator to the resource and IG prey, respectively; ρ 1 (R, N, P )R is the functional response of the IG prey to the basal resource; and m 1 and m 2 are density-independent morality rates. The parameters e 1 and e 2 are the conversion rates of resource consumption into reproduction for the IG prey and IG predator, respectively; the parameter e 3 denotes the conversion rate of the IG predator from its consumption of IG prey; Rϕ(R) is recruitment of the basal resource. Functional response describes how the consumption rate of individual consumers varies with respect to resource density and is often used to model predator-prey interactions. For IGP models, several functional response functions have been studied. For instance, Velazquez et al. [10] and Hsu et al. [11] investigated the case with a linear functional response. Abrams and Fung [12] considered Holling type-II functional response. Verdy and Amarasekare [13] and Freeze et al. [14] investigated Holling type-II and ratio-dependent functional responses, respectively. Kang and Wedekin [15] considered Holling-III functional response.
Note that the reproduction of predator following the consumption of prey is not instantaneous, but rather is mediated by some reaction-time lag required for gestation. Time delay plays an important role in ecology and it can induce very complex dynamical behaviors [16,17,18,19,20,21,22]. For IGP models, it has been shown that a time delay greatly impacts their dynamics [23,24]. In [24] where r is the growth rate of R in the absence of N and P , K is the carrying capacity of resource. c 1 is the predation rate of IG prey for resource, c 2 is the predation rate of IG predator for resource, c 3 is the consumption rate of IG predator to IG prey and all other parameters have the same meanings as those in (1).
Note that for each species, individuals tend to migrate towards regions with lower population densities. Hence the species are distributed over space and interact with each other within their spatial domains. To take spatial effects into consideration, reaction diffusion equations become a natural choice [25,26,27,28,29,30,31,32,33,34,35,36]. In this work, we consider a reaction diffusion IGP model with delay and Beddington-DeAngelis functional response.
Suppose Ω ⊂ R n is a bounded domain with smooth boundary ∂Ω. Let R(t, x), N (t, x), P (t, x) represent the densities of basal resource, IG prey and IG predator at time t and location x, respectively. The basal resource is assumed to grow logistically. We assume the basal resource is consumed by the IG prey at a rate c 1 R(t, x)N (t, x), and the IG prey is consumed by the IG predator is c 3 N (t, x)P (t, x) at time t and location x. In this paper, we will assume the functional response takes the Beddington-DeAngelis (B-D) form, i.e., the consumption of the resource by the IG predator is characterized by c2P (t,x)R(t,x) 1+a1R(t,x)+a2P (t,x) . The reproduction of IG prey from consuming the basal resource is e 1 c 1 R(t − τ, x)N (t − τ, x), where the time-lag parameter is introduced in a manner analogous to the treatment in [24]. We further assume the populations cannot cross the boundary of Ω. Our model then reads as where d 1 , d 2 , d 3 denote the diffusion coefficients of the three species, respectively; ∆ is the Laplacian operator in the n dimensional space, ν is the outward unit normal vector on ∂Ω, and the homogeneous Neumann boundary conditions reflect the situation where the population cannot across the boundary of Ω. The meanings and units of the parameters of model (3) are summarized in Table 1.
For rescalling, we let

RENJI HAN, BINXIANG DAI AND LIN WANG
Throughout the paper, we denoteΩ = Ω ∂Ω, The space of continuous functions inD T is denoted by C(D T ). For vector-value functions we use the product spaces with the usual inner product ·, · . The rest of the paper is organized as follows. In Section 2, we study the existence and uniqueness of solution of (4) and estimate the solution's priori bounds. In Section 3, we discuss the existence of nonnegative spatially homogeneous steady states. In Section 4, we carry out stability analysis and Hopf bifurcation analysis about the unique positive spatially homogeneous steady state of System (4). Numerical simulations are presented in Section 5 to illustrate the impacts of delay, diffusion and the functional response on the dynamics of our IGP model. We conclude this paper with a brief summary and discussion in Section 6.

2.
Existence of solution of System (4) and priori bound estimation. Theorem 2.1. Consider System (4), we have the following conclusions.
To estimate the priori bounds of u 2 (t, x) and u 3 (t, x), we set

RENJI HAN, BINXIANG DAI AND LIN WANG
Then From the Neumann boundary conditions, we further obtain Similarly, there exists Thus, for t ≥ T 2 + τ, we have and hence lim sup For the case with d 1 = d 2 = d 3 = d and τ = 0, we can similarly show that for any ε > 0, there exists Thus for t > T 3 , we have Consider the system   The comparison argument implies that This completes the proof.
3. Spatially homogeneous steady states of System (4). Same as in [24], we denote by R i = βi γi = eiciK mi (i = 1, 2) the reproduction numbers for the IG prey and IG predator, respectively. Consider (4), we easily have the following existence results on trivial and semi-trivial spatially homogeneous steady states. (ii) There is a weakly semi-trivial steady state in the absence of IG Prey and IG Predator E 1 = (1, 0, 0). (iii) The IG Prey-only strong semi-trivial steady state E 10 :

System (4) admits a positive steady state E
if E * is a positive solution to the following three equations: It follows from (10) that Combining (9), (11) and (12), we obtain where .
For the distribution of roots of Eq. (13), we have the following results.
According to Lemma 3.1, the following proposition is valid.

Dynamics of System
we have the following lemma from [37].

4.1.
Stability of E 0 and E 1 . In the following, we consider the stability of E 0 and E 1 . (ii) The semi-trivial steady state E 1 is locally asymptotically stable if R 1 < 1 and Proof. Linearizing System (4) at a constant steady state where and From Lemma 4.1, we know that the eigenvalues of the System (4) is confined on the subspace X i , and λ is an eigenvalue of (14) on X i if and only if it is an eigenvalue Then the characteristic equation of (14) is where I 3 stands for the 3 × 3 identity matrix.
(ii) If u ⋄ = E 1 , then we obtain the following characteristic equation Thus it follows from [24, Lemma 6] that the eigenvalues λ 2i have negative real parts for all µ i . It follows from [40, Corollary 1.11] that the equilibrium E 1 is locally asymptotically stable for R 1 < 1 and R 2 < 1 + b. If R 1 > 1 then there exists µ 1 = 0 such that d 2 µ 1 + γ 1 < β 1 , so it follows from the [24, Lemma 6] that at least one of the eigenvalues λ 2i has a positive real part. If R 2 > 1 + b then there exists µ 1 = 0 such that λ 31 > 0. From [40,Corollary 1.11], the steady state E 1 is unstable in either case.

4.2.
Stability of E 10 and E 01 . Next, we consider the stability of the two strong semi-trivial steady states: E 10 and E 01 . For the IG prey-only strong semi-trivial steady state E 10 , we have the following result.
Proof. For E 10 , the characteristic equation is with R1+b + β(1 − 1 R1 ) − γ 2 > 0 holds. Therefore, m 1 (λ) has at least one zero with a positive real part and the characteristic Eq. (20) has at least one positive root with a positive real part. The proof of (i) is complete. DenoteĒ Then we have SinceĒ 1 +L 1 > 0 andĒ 1L1 +J 1 > 0, it is easy to see that Eq. (21) has no positive zeros for all µ i .
Next, we discuss the distribution of positive zeros of G(u). Clearly F (0) = Thus, G(u) has no positive zeros for all µ i . It follows from [24,Lemma 11] that E 10 is locally asymptotically stable for all τ ≥ 0. This completes the proof of (ii).
For the IG predator-only strong semi-trivial steady state E 01 , we have the following result.
, then E 01 is locally asymptotically stable for all τ ≥ 0. Hereû 1 is defined in Proposition 3.1.

4.3.
Stability of the unique positive spatially homogenous steady state E * .

Stability and Hopf bifurcation.
In this subsection, by taking τ as the bifurcation parameter, we investigate the stability and the Hopf bifurcation near the unique positive spatially homogeneous steady state E * . For this purpose, we always assume (H 1 ) q < 0 and R 2 < The above assumption guarantees the uniqueness of the positive spatially homogeneous steady state E * . For the spatially homogeneous positive steady state E * = (u * 1 , u * 2 , u * 3 ), the characteristic equation is given below: where
Next, we discuss the effect of the delay τ = 0 on the stability of the positive steady state E * . We first determine critical values of τ at which a pair of simple purely imaginary eigenvalues appears.
Let λ = iω (ω > 0) be a root of Eq. (24). Substituting λ = iω into Eq. (24) yields It follows from (26) and (27) that From (24), we get In the following, we need to seek conditions required for Eq. (29) to have at least one positive root. For this purpose, we further make the following hypotheses: (H 5 ) According to the above analysis, we have the following lemma. Proof. It follows from (H 2 ) and (H 4 ) that b 0i + c 0i > 0 for all i ∈ N. Thus r i < 0 if and only if b 0i − c 0i < 0. It follows from (H 2 ), (H 5 ) and (H 6 ) that there exists µ i (i = 1, 2, · · · , N 2 ) such that r i < 0. Since lim s→∞ h(s) = ∞ for fixed µ i , then (29) has at least one positive root for each i ∈ {1, 2, · · · , N 2 }. This completes proof of part (i). Similarly, we can prove part (ii).
Differentiating both sides of (33) with respect to ω yields Substituting iω * into (32) yields If iω * is not simple, then iω * must satisfy (34) and (35), we have This is a contradiction. Thus ±iω * is a pair of simple purely imaginary roots of Eq. (24).
Since ±iω * are simple purely imaginary roots and we may consider λ = λ(τ ) to be a differentiable function. Differentiating (22) with respect to τ yields .
Using (32) again, we obtain Thus, from (36), we have This completes the proof.
Let us define C * = C([0, 1], R 3 * ), where R 3 * is the 3−dimensional vector space of row vectors, A * with domain dense in C * and range in C * . Let P and P * be the center subspace, that is, the generalized eigenspace of A and A * associated with In view of the definition of the two infinitesimal generators A and A * , we have the following conclusions.
It follows from Lemma 4.12 that Thus In what follows, in order to determine the bifurcation direction and stability, we compute the coordinates to describe the center manifold C 0 . As the formulas to be developed for the bifurcation direction and stability are all relative to µ = 0 only, we set µ = 0 in Eq. (4.24) and obtain a center manifold with the range in Q. The flow of Eq. (4.24) on the center manifold can be written as Moreover, from (39), z satisfieṡ where g(z,z) = g 20 z 2 2 + g 11 zz + g 02z By the Taylor expansion where Noting that (40), we get Substituting (38) into (42) and combining (41) yield  check that the hypotheses (H 1 ) − (H 6 ) hold and h ′ ((ω * ) 2 ) = 0.0938 > 0. By Theorem 4.9, local Hopf bifurcation occurs at τ * ≈ 0.7859. We use the forward Euler method to find numerical solutions to System (4) with τ = 0.7 < τ * and τ = 1.2 > τ * , respectively. As illustrated in Figures 1 and 2, when τ = 0.7 < τ * ,

5.2.
The spatiotemporal dynamics in two-dimensional space. Consider System (4) with u 1 = u 1 (t, x, y), u 2 = u 2 (t, x, y), u 3 = u 3 (t, x, y) and ∆ = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 . For this purpose, the domain of System (4) is confined to a fixed spatial domain Ω = [0, L] × [0, L] ⊂ R 2 with L = 400. we solve System (4) on a grid with 400 × 400 sites by a simple Euler method with a time step size of δt = 0.01 and a space step size of δx = δy = 1. The discretization of the Laplacian term takes the form where (i, j) denote the lattice sites and s = 1 is the lattice constant. The matrix elements of A l , A r , A d , A u are unity except at the boundary. When (i, j) lies on the left boundary, that is i = 0, we define A l (i, j)u(i − 1, j) ≡ u(i + 1, j), which guarantees zero-flux of individuals in the left boundary. Similarly we define A r (i, j), A d (i, j), A u (i, j) such that the zero-flux boundary condition is satisfied. With the given Neumann boundary conditions, the eigenvalues of −∆ on Ω are µ i = π 2 L 2 (n 2 + m 2 ), n, m ∈ Z, where Z represents the integer set. In order to discuss the impacts of delay and diffusion on the dynamics of System (4), we will compare the temporal model (that is, System (4) without diffusion) with System (4). We take parameters as (P 2 ) : α = 0.7, β = 0.9, β 1 = 1.95, β 2 = 1.85, γ 1 = 0.2, γ 2 = 0.8, b = 2.5, c = 5, d 1 = 1, d 2 = 2, d 3 = 4.
We consider two types of different initial conditions: (IC 2 ) : φ 1 (t, x, y) = u * 1 + 0.001 sin x sin y, φ 2 (t, x, y) = u * 2 + 0.001 sin x sin y, φ 3 (t, x, y) = u * 3 + 0.001 sin x sin y, and for (t, x, y) ∈ [−τ, 0] × Ω. Here ε 1 = 2 × 10 −7 , ε 2 = 3 × 10 −5 , ε 3 = 1.2 × 10 −3 , ε 4 = 6 × 10 −5 , ε 5 = 3 × 10 −5 . Under (P 2 ) and (IC 2 ), it is easy to check that all conditions of Theorem 4.9 are satisfied and there is a unique positive steady state E * ≈ (0.1754, 0.7419, 0.2029). The corresponding Hopf bifurcation value is computed as τ 1 1,0 ≈ 0.8028. Figure 3 depicts the population dynamics of the temporal model and the spatiotemporal model at τ = 1. Both the temporal model and the spatiotemporal If we increase τ to τ = 1.5 and use initial condition (IC ′ 2 ), as shown in Figure  4, the temporal model still exhibits regular oscillations, while the spatiotemporal model exhibits irregular oscillations and the calculated Lyapunov exponent is 0.0011 > 0 (By the method proposed in [42]), which indicates the occurrence of chaos. Figure 5 depicts the snapshots of the contour maps of specie u 1 for the temporal model and the spatiotemporal model at time t = 5000. The temporal model exhibits the spiral wave pattern. However, the spatiotemporal model presents the chaotic wave pattern. To have a better understanding on the evolution process of the  spatiotemporal pattern, in Figure 6, we present the snapshots of contour maps of the basal resource u 1 at t=200,500, 1000, 1200, 1500, 2500, respectively. As pointed out in [43], the spirals are usually observed under suitable parametric conditions near Turing-Hopf bifurcation threshold. In addition, in the spatially extended system the existence of a stable limit cycle normally results in the formation of chaotic spatiotemporal patterns [44]. As can be seen from Figure 6, as time t increases, an chaotic wave spatial pattern is gradually formed starting from a regular spiral wave pattern.

The effect of delay.
To explore the impact of delay, in Figure 7, we take the snapshots of the contour maps of specie u 1 at time t = 1500 for several different values of τ . As can be seen in Figure 7, the time delay can lead to the formation of an irregular spatial pattern from a regular spiral pattern in the whole domain as the time delay increases and surpasses some critical value.

5.2.2.
The effect of diffusion. As seen from Figure 6, System (4) has a regular spiral wave pattern when τ = 1.5, d 1 = 1, d 2 = 2 and d 3 = 4 at t = 1500. To numerically examine how the diffusion affects the pattern, we take the snapshots of contour maps of u 1 at t = 1500 with several different choices of the diffusion coefficients. As shown in Figure 8, spiral wave pattern emerges firstly when d 1 = d 2 = d 3 = 0, then as the three diffusion coefficients change to d 1 = d 2 = 0.01, d 3 = 0.04, the spiral wave structure disappears around the center of the spirals wave, with the increase in these diffusion coefficients, it grows steadily, and eventually the chaotic wave pattern dominates the whole domain. Differing from the instability mechanism in Figure 7, the spirals wave loses its stability due to the Doppler effect [45].

5.2.3.
The impact of the prey saturation constant b. Figure 9 demonstrates how the prey saturation constant constant b affects the pattern formation of the basal resource u 1 at time t = 1500 with τ = 1.5 : when b is small, we observe a pattern with stripes firstly; as the constant b increases, the spiral wave pattern emerges, then it grows steadily, as b goes beyond a certain value, chaotic wave pattern appears.

5.2.4.
The impact of the predator interference constant c. To see how the predator interference constant c influences the spatiotemporal pattern, we numerically simulation (4) with different values of c and plot the snapshots of the contour maps of the basal resource u 1 at t = 1500 in Figure 10. It is illustrated in Figure 10 that the predator interference constatn c can also lead to the formation of chaotic wave spatial pattern, which can be preceded from the evolution of a regular spiral patterns as the predator interference constant c decreases. 6. A summary and discussion. In this work, we have investigated the spatiotemporal dynamics of a diffusive IGP model with delay and the Beddington-DeAngelis functional response. we have established locally asymptotically stability results of the trivial, semi-trivial and strong semi-trivial steady states. In the case that there is a unique positive spatially homogeneous steady state E * , we have carried out the Hopf bifurcation analysis. Unlike competition models with monotone response functions ( [17]) where delay does not induce sustained oscillations, in our IGP models, delay promotes complex dynamics including bistability, and the emergence of spiral wave pattern and chaotic wave pattern.
Compared with the temporal model in [24], we also observe bistability is possible in System (4). In addition, the diffusion also has impacts on the formation of spatiotemporal patterns as it can change the distribution of characteristic roots of the corresponding characteristic equations, and hence has an important effect on the dynamics for the constant steady state of System (4). This has been illustrated via numerical simulations as well (See Figure 8). Moreover, we have observed that the functional response can also influence the formation of complex patterns. As demonstrated in Figures 9 and 10, the functional responses can also trigger the emergence of spiral wave pattern and chaotic wave spatial pattern.