Mathematical Biosciences 259 (2014) 1–11

Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Population dynamics of intraguild predation in a lattice gas system Yuanshi Wang, Hong Wu∗ School of Mathematics and Computational Science, Sun Yat-sen University, Guangzhou 510275, PR China

a r t i c l e

i n f o

Article history: Received 17 June 2014 Revised 1 November 2014 Accepted 4 November 2014 Available online 11 November 2014 Keywords: Mean-field theory Parasitism Competition Persistence Bifurcation

a b s t r a c t In the system of intraguild predation (IGP) we are concerned with, species that are in a predator–prey relationship, also compete for shared resources (space or food). While several models have been established to characterize IGP, mechanisms by which IG prey and IG predator can coexist in IGP systems with spatial competition, have not been shown. This paper considers an IGP model, which is derived from reactions on lattice and has a form similar to that of Lotka–Volterra equations. Dynamics of the model demonstrate properties of IGP and mechanisms by which the IGP leads to coexistence of species and occurrence of alternative states. Intermediate predation is shown to lead to persistence of the predator, while extremely big predation can lead to extinction of one/both species and extremely small predation can lead to extinction of the predator. Numerical computations confirm and extend our results. While empirical observations typically exhibit coexistence of IG predator and IG prey, theoretical analysis in this work demonstrates exact conditions under which this coexistence can occur. © 2014 Elsevier Inc. All rights reserved.

1. Introduction Intraguild predation (IGP) occurs between species in the same community which utilize similar resources (space or food), and thus there is competition between them. IGP is classified as asymmetrical or symmetrical. In asymmetrical IGP, one species consistently preys upon the other. But in symmetrical IGP, either species can prey upon the other. In this work, we are concerned with asymmetrical IGP. For convenience, we call it as IGP in the following discussions. IGP has been widely observed in both terrestrial and aquatic communities [27]. For examples, there exists IGP between large mammalian carnivores. Large canines and felines are the mammal groups often involved in IGP, with larger species such as lions and gray wolves preying upon smaller species such as foxes and lynx [14]. Coyotes function as predators on gray foxes and bobcats in North America [3]. Since empirical observations typically exhibit coexistence of IG predator and IG prey, an interesting question has been focused on that under which conditions the coexistence can occur [16,17]. Several models have been established to characterize IGP. Holt and Polis [5] formed a three-specie model of IGP with Holling Type I functional response, in which two species that have a predator–prey relationship, also compete for a shared resource. They also formed a two-species system, which is an extension of the exploitative competition model introduced by Schoener [19,20]. Theoretical analysis demonstrated a general criterion for coexistence in IGP that IG prey

∗

Corresponding author. Tel.: +86 2084034848. E-mail address: [email protected] (H. Wu).

http://dx.doi.org/10.1016/j.mbs.2014.11.001 0025-5564/© 2014 Elsevier Inc. All rights reserved.

should be superior at exploitive competition for the shared resource, whereas the IG predator should gain significantly from the IG prey. Local stability analysis and numerical computations also showed that at intermediate levels of environmental productivity, there exist alternative states (either IG prey dominance or IG predator dominance, or either the IG predator dominance or coexistence). Ruggieri and Schreiber [18] considered the Schoener–Polis–Holt model, in which IG prey and IG predator compete for resources (food). A global analysis exhibits six dynamics of the model. Okuyama [28] studied IGP in a spatial setting by establishing a lattice IGP model. For homogeneous resources, pair approximation was used to study the effect of spatially structured species interactions. The qualitative results of the pair approximation model predicted coexistence of the species over a wider range of parameters than the non-spatial model. Takimoto et al. [22,23] analyzed models of IGP with three and four dimensions, and demonstrated complex but systematic sequences of alternative states along a productivity gradient, where sufficient conditions that determine which sequences to occur, are clarified. Takimoto et al. [24] presented a Levins-type patch occupancy model of IGP. Theoretical analysis exhibited conditions for feasibility of each equilibrium, while numerical computations displayed both equilibrium stability and food-chain length with different strength of local IGP. The model showed that ecosystem size can promote coexistence and increase food-chain length even when local IGP is strong, which is consistent with empirical patterns. Kang and Wedekin [7] formed a model of IGP with Holling Type III functional response. Sufficient conditions are provided for persistence and extinction of species in all possible situations and multiple attractors and periodic oscillations are exhibited. For other relevant work, we refer to

2

Y. Wang, H. Wu / Mathematical Biosciences 259 (2014) 1–11

Refs. [1–3,10,13]. As far as we know, mechanisms by which IG prey and IG predator can coexist in IGP models with spatial competition, have not been shown. Therefore, forming an appropriate model of IGP with spatial competition and demonstrating its properties are necessary. In this work, we use a lattice gas model to describe IGP. The lattice gas model is individual-based and is different from the classical “lattice model”. Since the lattice model characterizes local interactions occurring between adjacent sites on a lattice, the lattice gas model describes interactions between any pair of lattice sites and the interactions occur randomly and independently [21]. While dynamics of lattice models cannot be described by mathematical equations, dynamics of lattice gas models are usually depicted by differential equations when the lattice is sufficiently large, which are called the mean-field theory of lattice model. These equations have been widely applied in characterizing competition, predation and mutualism in biology for years. For example, a lattice version of predator–prey model was presented by Nakagiri et al. [11] and dynamics of the model exhibit conditions under which the species can coexist and under which the predator would be excluded by competition of the prey. A typical lattice version of competition model has also been studied by several authors [8,9,12]. In a recent study, a lattice version of mutualismcompetition model was introduced by Iwata et al. [6]. In the twospecies system, spatial competition is considered and the space where the two species live, is assumed to be divided into many sites and is regarded as a lattice. Since each site can be occupied by one individual of the populations, the species are competitive for sites on the lattice. The species are also mutualistic since each of them produces resources for the other. The benefit of the mutualism is represented by the increased reproduction rate of both species populations. Dynamics of the mixed mutualism-competition model demonstrate that when the mutualistic effects vary, interaction outcomes between the species can change among mutualism, parasitism and competition in a smooth fashion. Inspired by the work of Nakagiri et al. [11] and Iwata et al. [6], we apply a lattice version of predation-competition model to describe IGP. The model focuses on spatial competition. The space occupied by the population is regarded as a lattice and each site can be occupied by one individual of the species. Thus, the species are competitive for sites on the lattice. The species are also involved in a predator–prey relationship since predation can occur as a predator meets the prey. Similar to Ref. [6], the benefit of the predation is represented by the increased reproduction rate of the predator. Dynamics of the model demonstrate properties of IGP and mechanisms by which the IGP leads to persistence/extinction of species: (i) In IGP where IG predator cannot survive in the absence of IG prey, when the efficiency of predator in converting its consumption into fitness is high, the predator can persist. When the efficiency is intermediate, the predator can persist if it has a high initial density; if the initial density is low, the predator goes to extinction. When the efficiency is low, the predator will go to extinction. (ii) In IGP where IG predator can survive in the absence of IG prey, when the efficiency of predator is high, the IGP can enhance population density of predator. When the efficiency is low, the predator can persist if its population density is large. Otherwise, the predator will be driven into extinction by its prey through competition. (iii) Intermediate predation can lead to persistence of the predator, while extremely big predation can lead to extinction of one/both species and extremely small predation can lead to extinction of the predator. Saddle-node bifurcation and pitchfork bifurcation in the model are demonstrated, while numerical computations confirm and extend our results. The paper is organized as follows. The model is described in Section 2. Sections 3 and 4 show dynamics of the model. Discussion is in Section 5.

2. Model In this section, we form the lattice gas model of IGP and exhibit boundedness of solutions and nonexistence of periodic orbit of the model. Let X and Y represent the prey and predator, respectively. A site on a lattice is labeled by X (or Y) if it is occupied by an individual of species X (or Y). When a site is empty, it is labeled by O. On the lattice, any pair of sites can interact randomly and independently. When there is only one species on the lattice, the interactions can be described by contact process [4]. For example, in the system of species X, if a site is occupied by X, then it will become O in a mortality rate mX . If sites X and O interact, then the site O will become X in a birth rate BX . A similar discussion can be given for the one-species system of Y. Moreover, when the prey X and predator Y emerge on the same lattice, predation can occur: X would be killed and consumed if it meets Y, which promotes the growth of Y. Therefore, reactions on the lattice of species X and Y can be depicted as follows:

X→O

with mortality rate mX ,

Y→O

with mortality rate mY ,

X + O → 2X

with birth rate BX ,

Y + O → 2Y

with birth rate BY ,

X+Y →O+Y

with killing rate k,

(2.1)

where the first and second reactions respectively describe the mortality processes of species X and Y, while the third and fourth reactions respectively characterize their birth processes. The fifth reaction represents the predation of species Y on X, in which species X is killed by Y and the site occupied by X becomes empty. When the lattice size is large, reactions of (2.1) are usually described by differential equations, which are called lattice models of mean-field theory [6,26]:

dx = −mX x + BX x(1 − x − y) − kxy, dt dy = −mY y + BY y(1 − x − y), dt

(2.2)

where x and y represent fractions of sites occupied by species X and Y, respectively. The factor (1 − x − y) is the fraction of empty sites. Since one site can be occupied by one individual, x and y represent the sizes of population densities of species X and Y. For convenience, x and y are called densities of the two species. The first and second terms in the righthand side of each equation come from the mortality and birth processes, while the third term in the first equation comes from the predation process. Thus, the lattice gas model (2.2) has the same form as Lotka–Volterra equations. Parameters mX and mY respectively represent mortality rates of species X and Y, and the birth rates are defined by

BX = rX ,

¯ BY = rY + ekx,

where rX (rY ) denotes the birth rate of species X (Y) in the absence of the other. Parameter k represents the killing rate of species Y on X, while e¯ denotes the efficiency of Y in converting the consumption into fitness. In the following discussion, we denote

¯ e = ek. In system (2.2), we assume rY ≥ 0 and all other parameters are positive. We consider solutions (x(t), y(t)) of (2.2) with initial values x(0) > 0, y(0) > 0. Thus we have x(t) > 0, y(t) > 0 as t > 0. When rX ≤ mX , we have dx/dt < 0 by the first equation of (2.2). Then species X goes to extinction and system (2.2) becomes a one-species system

Y. Wang, H. Wu / Mathematical Biosciences 259 (2014) 1–11

3

(a) The trivial equilibrium O(0, 0) always exists and has eigenvalues r1 (1 − m1 ), r2 (1 − m2 ). (b) The semi-trivial equilibrium P1 (1 − m1 , 0) always exists since 1 − m1 > 0. P2 (0, 1 − m2 ) exists if 1 − m2 > 0. The eigenvalues of Pi are −ri (1 − mi ), −rj Cj with

without intraguild predation. Therefore, we assume in this work

rX > mX . The following result demonstrates boundedness of solutions and nonexistence of periodic orbit of (2.2).

C1 = m2 − [1 + α(1 − m1 )]m1 ,

Theorem 2.1. Solutions of (2.2) are bounded and system (2.2) admits no periodic orbit.

C2 = m1 − 1 + (1 + β)(1 − m2 ),

Proof. When x + y ≥ 1, we have dx/dt < 0 and dy/dt < 0, which implies that all solutions of (2.2) satisfy x(t) + y(t) ≤ 1 as t is sufficiently large. Thus solutions of (2.2) are bounded. Let φ(x, y) = 1/(xy), which is called the Dulac function. Let f (x, y) and g(x, y) be the right-hand sides of equations in (2.2), respectively. Then we have

i = j, i, j = 1, 2.

(c) There are at most two interior equilibria P − (x− , y− ) with

x± =

√ −B ±  , 2A

y± =

(3.4)

P + (x+ , y+ )

1 − m1 − x± , 1+β

and

(3.5)

where x± are roots of the equation

rX rY + ex ∂(φ f ) ∂(φ g) + = −k − − 0, y > 0. It follows from Bendixson–Dulac theorem [15] that system (2.2) admits no periodic orbit.

A = αβ , B = β − αβ − α m1 ,

(3.6)

with

C = m2 (1 + β) − β − m1 ,

¯  = B − 4AC = A¯β 2 + B¯ β + C, 2

It follows from Theorem 2.1 that all solutions of (2.2) converge to equilibria. In the following analysis, we consider cases of rY > 0 and rY = 0, which imply that species Y can/cannot reproduce in the absence of X, respectively. Remark 2.2. Since X and Y compete for empty space, Y may replace X after killing it. In a lattice gas system, it is assumed that one individual of each species can only occupy one site. Thus, if the site occupied by X is replaced by Y, then the previous site of Y will become empty and the reaction becomes X + Y → Y + O, which does not affect the model (2.2). 3. The case of rY > 0 In this section, we consider the case of rY > 0, which implies that species Y can reproduce in the absence of X. Let

mX , rX

m1 =

m2 =

mY , rY

β=

k , rX

α=

e , rY

r1 = rX , r 2 = rY , (3.1)

then model (2.2) becomes

dx = r1 x[1 − m1 − x − (1 + β)y], dt dy = r2 y[−m2 + (1 + α x)(1 − x − y)]. dt

with

A¯ = (1 + α)2 − 4α m2 , C¯ = α 2 m21 .

B¯ = 2α(m1 + α m1 − 2m2 ),

The following result shows stability of equilibria P ± . Lemma 3.1. If P + and P − are interior equilibria of (3.2), then P+ is a saddle point and P − is locally asymptotically stable. If P + and P − coincide, then the Jacobian matrix of (3.2) at P+ has a simple zero eigenvalue. Proof. A direct computation shows that the Jacobian matrix A of (3.2) at P ± satisfies

trA|P± = −r1 x± − r2 y± (1 + α x± ) < 0,

√ detA|P± = ∓ r1 r2 β x± y± .

Thus det A|P− > 0 and det A|P+ < 0, which implies that when P − is an interior equilibrium of (3.2), it is locally asymptotically stable; when P + is an interior equilibrium of (3.2), it is unstable and is a saddle point. When  = 0, we have x+ = x− and det A|P+ = 0, which implies the two equilibria P + and P − coincide. Since trA|P+ < 0, the Jacobian matrix A at P + has a simple zero eigenvalue as  = 0. Dynamics of system (3.2) can be determined by the relative position of isoclines l1 and l2 of (3.2). Denote

l1 : y = f (β , x) = (3.2)

For simplicity, in the following discussion of this work, parameter β is regarded as the killing rate of species Y on X, and α as the efficiency of species Y in converting the consumption into fitness. From rX > mX , we have

l2 :

1 − m1 − x , 1+β

(1 + α x)(1 − x − y) = m2 with y = g(α , x) = 1 − x −

m2 . 1 + αx

The x-isocline l1 and the x-axis intersect at P1 (1 − m1 , 0). l1 is a line with

lim f (β , x) = 1 − m1 − x,

lim f (β , x) = 0.

(3.7)

1 − m1 > 0.

β →0

In order to analyze stability of equilibria, we compute the Jacobian matrix A of (3.2) by

The y-isocline l2 and the y-axis intersect at P2 (0, 1 − m2 ). l2 is a hyperbola with asymptotes 1 − x − y = 0 and x = −1/α , as shown in Fig. 1. Since

A=

 a11

a12

a21

a22



,

(3.3)

β →+∞

∂ 2 g/∂ x2 = −2m2 α 2 (1 + α y)−3 < 0,

where

¯ y¯ ) l2 is convex upward with a vertex (x,

a11 = r1 [1 − m1 − x − (1 + β)y] − r1 x,

¯ y¯ ) = − (x,

a12 = −r1 (1 + β)x,



a21 = r2 y[α(1 − x − y) − 1 − α x],

a22 = r2 [−m2 + (1 + α x)(1 − x − y)] − r2 y(1 + α x). The equilibria of (3.2) are considered as follows, while their local stability is determined by eigenvalues of Jacobian matrix A at the corresponding equilibria.

1

α



+

m2

α

,1 +

1

α



−2

m2

α

 .

Since ∂ g/∂α = m2 x/(1 + α x)2 > 0, we have g(α¯ , x) > g(α , x) when α¯ > α and x > 0, as shown in Fig. 1. Since we focus on nonnegative solutions of (3.2), l2 ∩ intR2+ = ∅ is necessary. It follows from the convexity of l2 that l2 ∩ intR2+ = ∅ if and only if l2 and the positive x-axis have at least one intersection point

4

Y. Wang, H. Wu / Mathematical Biosciences 259 (2014) 1–11

Hyperbola l2 with m2=0.9 α =1.8 α = 4.4 α = 20

y

1−x−y=0

x Fig. 1. The y-isocline l2 . We fix m2 = 0.9 and let α vary. When α increases from 1.8, 4.4 to 20, l2 increases monotonically. If α tends to infinity, then l2 converges to its asymptote 1 − x − y = 0, which is the dotted line in the figure.

(except the case of being tangent). That is, the following equation has at least one positive root x:

α x2 − (α − 1)x + m2 − 1 = 0. Thus, if 1 − m2 > 0, then l2 ∩ intR2+ = ∅. If 1 − m2 = 0, then l2 ∩ intR2+ = ∅ as α > 1. If 1 − m2 < 0, then l2 ∩ intR2+ = ∅ as α > α0 with



α0 = 2m2 − 1 + 2 m2 (m2 − 1). We consider dynamics of (3.2) in cases of 1 − m2 ≤ 0 and 1 − m2 > 0. When 1 − m2 ≤ 0, we have the following Theorem 3.2, while the proof is in Appendix A. Denote

α1 β2

m2 − m1 = , m1 (1 − m1 )  ¯ −B¯ +  , = ¯ 2A

β1

m2 (1 − m1 ) = 2 − 1, m1 + (1 − 2m1 )m2

¯ = B¯ 2 − 4A¯C. ¯ Then when α > α1 , P1 is below l2 and C1 < 0. where  β1 is defined to compare slopes of l1 and l2 at P1 . When β = β2 , l1 and l2 are tangent in the first quadrant as α > α0 . Theorem 3.2. Assume 1 − m2 ≤ 0. (i) If α > α1 , then equilibrium P − (x− , y− ) is globally asymptotically stable, as shown in Fig. 2a. (ii) Let α < α1 . Then system (3.2) experiences a saddle-node bifurcation at β = β2 , as shown in Fig. 2b–d. If α > α0 and β > β2 , then system (3.2) has two interior equilibria P − and P + . P − is asymptotically stable while P+ is a saddle point, as shown in Fig. 2b. If α > α0 and β = β2 , then equilibria P− and P+ coincide and form a saddle-node point, as shown in Fig. 2c. The separatrices of the saddle point divide the first quadrant into two regions, one is the basin of attraction of P1 while the other is that of P − . In other

cases, all positive solutions of (3.2) converge to equilibrium P1 , as shown in Fig. 2d. (iii) Let α = α1 . Then system (3.2) experiences a pitchfork bifurcation at α = α1 and experiences a degenerate transcritical bifurcation at β = β1 : If 2m1 − 1 ≤ 0 and β > β1 , or, if 2m1 − 1 > 0, m2 < m21 /(2m1 − 1) and β > β1 , then equilibrium P − (x− , y− ) is globally asymptotically stable. In other cases, all positive solutions of (3.2) converge to equilibrium P1 (1 − m1 , 0). When 1 − m2 > 0, we have the following Theorem 3.3, while the proof is in Appendix B. Denote

α2 =

m2 − m1 , m2 (1 − m1 )

β3 =

m2 − m1 . 1 − m2

Thus, when β > β3 , P2 is above l1 and C2 > 0. When α = α2 and β = β3 , l1 and l2 are tangent at P2 . Theorem 3.3. Assume 1 − m2 > 0. (i) Let β > β3 . If α ≥ α1 , then equilibrium P2 is globally asymptotically stable, as shown in Fig. 4a. If α < α1 , then equilibrium P + is a saddle point. The separatrices of P + divide the first quadrant into two regions, one is the basin of attraction of P1 while the other is that of P2 , as shown in Fig. 4b. (ii) Let β = β3 . (iia) If α ≥ α1 , then P2 is globally asymptotically stable. (iib) If α2 < α < α1 , then P + is an interior equilibrium of (3.2) and is a saddle point. The separatrices of P + divide the first quadrant into two regions, one is the basin of attraction of P1 while the other is that of P2 , as shown in Fig. 4c. (iic) If α ≤ α2 , then equilibrium P1 is globally asymptotically stable, as shown in Fig. 4d. (iii) Let β < β3 . (iiia) If α ≥ α1 , then equilibrium P − is globally asymptotically stable, as shown in Fig. 5a.

Y. Wang, H. Wu / Mathematical Biosciences 259 (2014) 1–11

5

(b) α = 10, β = 0.6

(c) α = 10, β = 0.3

(d) α = 10, β = 0.1

x

x

y

y

(a) α = 25, β = 0.01

Fig. 2. Phase-plane diagrams of (3.2). Red and blue curves are the isoclines for X and Y, respectively. The black lines are separatrices of saddle points. The separatrices subdivide the plane into two regions, which are the basins of attraction of the corresponding stable equilibria. The stable and unstable equilibria are denoted by filled and open circles, respectively. Gray arrows display the direction and strength of the vector fields in the phase-plane space. We fix m1 = 0.1, m2 = 1.2 and let α , β vary. (a) When the efficiency of Y is high (α = 25, β = 0.01), Y can persist through predating X. (b) and (c) When the efficiency of Y is low but predation of Y on X is big (α = 10, β = 0.6 or 0.3), Y can persist if its population density is large and that of X is small. Otherwise, Y goes to extinction. (d) When the efficiency of Y becomes low (α = 10, β = 0.1), Y cannot survive in the system. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

(iiib) Let α < α1 . Then system (3.2) experiences a saddle-node bifurcation at β = β2 . If β > β2 , then system (3.2) has two interior equilibria P − and P + . P − is asymptotically stable while P + is a saddle point, as shown in Fig. 5b. If β = β2 , then equilibria P − and P + coincide and form a saddle-node point, as shown in Fig. 5c. The separatrices of the saddle point divide the first quadrant into two regions, one is the basin of attraction of P1 while the other is that of P − . If β < β2 , then all positive solutions of (3.2) converge to P1 , as shown in Fig. 5d.

on theoretical analysis on the model, while the relevant ecological meaning is similar to that of Theorem 3.2 and is omitted. Let

Remark 3.4. The condition of m1 < m2 is unnecessary in this work, which is explained as follows. Theorem 3.2 considers the case of 1 − m2 ≤ 0. Since we assume 1 − m1 > 0 in the whole work (see Section 2), we have m1 < m2 . Theorem 3.3 considers the case of 1 − m2 > 0. If m1 ≥ m2 , then α1 < 0, β3 < 0. Thus, the condition in Theorem 3.3(i) is naturally satisfied and P2 (0, 1 − d2 ) is globally asymptotically stable. The ecological meaning is that species X is driven into extinction by Y. Therefore, the condition of m1 < m2 is necessary for persistence of IG-prey. Recall that m1 = mX /rX , m2 = mY /rY . The meaning of m1 < m2 is that the IG-prey is better competitor than IG-predator when there is no IGP interaction, which is consistent with the criterion for coexistence in IGP that IG-prey should be superior at exploitive competition for the shared resource, as shown by Holt and Polis [5].

(4.1)

4. The case of rY = 0 In this section, we consider model (2.2) with rY = 0, which implies that predator Y cannot reproduce in the absence of X. We focus

m1 =

mX , rX

β=

m2 = mY ,

k , rX

α = e, r1 = rX ,

then model (2.2) becomes

dx = r1 x[1 − m1 − x − (1 + β)y], dt dy = y[−m2 + α x(1 − x − y)]. dt

where parameters are positive and have the same meanings as those in model (2.2). Dynamics of system (3.2) can be determined by the relative position of isoclines l1 and l¯2 of (3.2). Denote

l1 : y =

1 − m1 − x , 1+β

l¯2 :

α x(1 − x − y) = m2 .

Since we focus on nonnegative solutions of (4.1), l¯2 ∩ intR2+ = ∅ is necessary. It follows from the convexity of l¯2 that l¯2 ∩ intR2+ = ∅ if and only if l¯2 and the positive x-axis have at least one intersection point (except the case of being tangent). That is, the following equation has at least one positive root x:

α x2 − α x + m2 = 0. Thus, l¯2 ∩ intR2+ = ∅ if α > α¯ 0 with

α¯ 0 = 4m2 .

6

Y. Wang, H. Wu / Mathematical Biosciences 259 (2014) 1–11

Equilibria of (4.1) are considered as follows, while their local stability is determined by eigenvalues of Jacobian matrix of (4.1) at the corresponding equilibria. (a) The trivial equilibrium O(0, 0) always exists and has eigenvalues r1 (1 − m1 ), −m2 . Thus, O is a saddle point. (b) The semi-trivial equilibrium P1 (1 − m1 , 0) exists and has eigenvalues −r1 (1 − m1 ), −m2 + α m1 (1 − m1 ). Denote

α¯ 1 =

m2 . m1 (1 − m1 )

Then if α > α¯ 1 , P1 is below l¯2 and is a saddle point. (c) There are at most two interior equilibria P + (x+ , y+ ) and P − (x− , y− ) with

 −B1 ± 1 x = , 2A1 ±

y± =

1 − m1 − x± . 1+β

(4.2)

x± are roots of the equation

A1 x2 + B1 x + C1 = 0, with

A1 = αβ ,

1 =

B21

B1 = −α(β + m1 ),

C1 = m2 (1 + β),

− 4A1 C1 = A¯1 β 2 + B¯ 1 β + C¯ 1 ,

with

A¯ 1 = α − 4m2 ,

B¯ 1 = 2α m1 − 4m2 ,

C¯ 1 = α m21 .

When α > α¯ 1 , P1 is below l¯2 . It follows from the convexity of l¯2 that there is a unique equilibrium P − of (4.1). When α ≤ α¯ 1 , P1 is above l¯2 . Since B1 < 0 and C1 > 0, it follows from the convexity of l¯2 that there exist interior equilibria of (4.1) if and only if α > α¯ 0 and β ≥ β¯ 2 with

β¯ 1 =

m1 , 1 − 2m1

 ¯1 −B¯ 1 +  ¯ β2 = , ¯ 2 A1

¯ 1 = B¯ 2 − 4A¯ 1 C¯ 1 = 16m1 m2 (1 − m1 )(α¯ 1 − α) ≥ 0.  1 By a proof similar to that of Theorem 3.2, we conclude the following result. Theorem 4.1. (i) If α > α¯ 1 , then equilibrium P − (x− , y− ) is globally asymptotically stable. (ii) Let α = α¯ 1 . If 2m1 − 1 < 0 and β > β¯ 1 , then equilibrium P − (x− , y− ) is globally asymptotically stable. In other cases, all positive solutions of (4.1) converge to equilibrium P1 (1 − m1 , 0). (iii) Let α < α¯ 1 . Then system (4.1) experiences a saddle-node bifurcation at P + . If α > α¯ 0 and β > β¯ 2 , then system (4.1) has two interior equilibria P − and P + . P − is asymptotically stable while P + is a saddle point. If α > α¯ 0 and β = β¯ 2 , then equilibria P − and P + coincide and form a saddle-node point. The separatrices of the saddle point divide the first quadrant into two regions, one is the basin of attraction of P1 while the other is that of P − . In other cases, all positive solutions of (4.1) converge to equilibrium P1 . 5. Discussion In this paper, we consider a lattice gas model of IGP. Dynamics of the model demonstrate properties of IGP. First, intraguild predation can promote growth of the predator. In the situation considered by Theorem 3.2, species Y cannot survive in the absence of X while X can persist alone. As shown in Theorem 3.2(i), predator Y can survive when its efficiency in converting the consumption into fitness is high (α > α1 ). The reason is that the predation of species Y on X provides more food for Y, which leads to survival of

Y who has high efficiency. Since predator Y cannot persist alone, it is the intraguild predation that leads to its survival. In the situation considered by Theorem 3.3, each species can persist in the absence of the other. When there is no intraguild predation, species Y will approach its carrying capacity 1 − d2 . However, when there exists such predation and the efficiency of Y is high (α > α1 ), Theorem 3.3(iiia) exhibits that species Y can approach a density y− , which is larger than its carrying capacity, as shown in Fig. 5a. The reason is that the predation provides more food for Y, which promotes its growth. Thus, intraguild predation can enhance population density of predators. Second, initial densities of species can play an important role in survival of species. In Theorem 3.2(ii), the efficiency of species Y is intermediate (α0 < α < α1 ). When the initial density of species Y is large, Y can persist. Otherwise, Y goes to extinction. Therefore, it is the initial density that determines survival of predator in Theorem 3.2(ii). In Theorem 3.3(i), when the efficiency of predator Y is low (α < α1 ), Y can persist if its population density is large. The reason is that species Y with the large density can occupy more empty sites, which leads to its survival in spatial competition with X. If the density of Y is small, Y goes to extinction, which displays a typical phenomenon of bi-stability in Lotka–Volterra competitive model. Since Y can persist in the absence X, a small density would lead to its extinction in the presence of X. Thus, population density is important in the survival of species. A similar discussion can be given for Theorem 3.3(ii) and (iii). Third, intermediate predation can lead to persistence of the predator. As shown in Theorem 3.2(ii) and Fig. 2b, it follows from the convexity of l2 that there exists an intermediate predation β¯ (β2 < β¯ < +∞) such that when β = β¯ , species Y can approach a maximal density y¯ − . Over-predation (β > β¯ ) and under-predation (β < β¯ ) cannot lead to the maximum, which we can see in Fig. 2b by varying β . The reason is that over-predation leads to a decreased density of species X, which provides decreased food for Y. Under-predation leads to an increased density of species X, which forms an increased competition with Y on lattice. A similar discussion can be given for Theorem 3.2(i) and (ii). Therefore, intermediate predation is the best for predators when in coexistence. A similar discussion can demonstrate that intermediate efficiency (α ) is beneficial to predators when in coexistence. In Theorem 3.3, the slope of l2 at P2 is k¯ 2 = α m2 − 1. When k¯ 2 > 0 ¯ y¯ ) (α > 1/m2 ), it follows from the convexity of l2 that its vertex P¯ (x, is in the first quadrant as shown in Theorem 3.3(iiia) and Fig. 5a. Thus, it follows from the convexity of l2 that there exists an intermediate predation β¯ (0 < β¯ < β3 ) such that when β = β¯ , species Y can approach a maximal density y¯ − . When y¯ = y¯ − , that is, P¯ = P − , over-predation (β > β¯ ) and under-predation (β < β¯ ) are not the best for predators. The reason is that over-predation leads to a decreased density of species X, which provides decreased food for Y. Underpredation leads to an increased density of species X, which forms an increased competition with Y on lattice. A similar discussion can be given for Theorem 3.3(iiib). Therefore, intermediate predation is the best for predators when in coexistence. A similar discussion can demonstrate that intermediate efficiency (α ) is beneficial to predators when in coexistence. Fourth, IG predator can be driven into extinction by IG prey. In Theorem 3.3(iiib), when the efficiency of predator Y is low (α < α1 ) and predation is small (β < β2 ), Y goes to extinction. The reason is that under these conditions, predator Y fails in spatial competition with X. Since Y can persist in the absence of X, it is the low efficiency and small predation that leads to its extinction. A similar discussion can be given for Theorem 3.3(i) and (ii). Fifth, intermediate levels of environmental productivity (rX and rY ) can lead to coexistence of IG prey (X) and IG predator (Y), while the high levels would lead to extinction of species. (a) In Theorem 3.2, species Y cannot survive in the absence of X while X can persist alone. That is, rY < mY and rX > mX by (3.1). Theorem 3.2(i) exhibits that

Y. Wang, H. Wu / Mathematical Biosciences 259 (2014) 1–11

when α > α1 , species X and Y coexist. Since α > α1 can be written as mX < rX < α mX by (3.1), an intermediate level of environmental productivity (rX ) can lead to the coexistence. However, when the level is high (rX > α mX ), we have α < α1 , which implies extinction of species Y as shown in Theorem 3.2(ii) and Fig. 2d. (b) In Theorem 3.3, each species can survive in the absence of the other. That is, rX > mX and rY > mY by (3.1). Theorem 3.3(iiia) exhibits that when α ≥ α1 and β < β3 , species X and Y coexist. Since β < β3 can be written as

rX >

k(rY − mY ) + rY mX mY

an

intermediate level of environmental productivity (i.e., < rX ≤ α mX ) can lead to the coexistence. However, when the level is high (rX > α mX ), we have α < α1 , which implies extinction of species Y as shown in Theorem 3.3(iiib) and Fig. 5d. A similar discussion can be given for the extinction of species X when rY is high by Theorem 3.3(i) and Fig. 4a. Finally, extremely big predation can lead to extinction of one/both species and extremely small predation can lead to extinction of the predator. In Theorem 3.2, when the predation is extremely big (β → +∞), isocline l1 tends to the x-axis by (3.7) and we have limβ →+∞ y− = 0 in equilibrium P − (x− , y− ), which implies that species Y goes to extinction. The reason is that under the extremely big predation, species X cannot approach a density that could lead to survival of Y. In particular, extremely big predation can lead to extinction of both X and Y if 1 − m2 = 0. Indeed, when 1 − m2 = 0 in Theorem 3.2, we have limβ →+∞ x− = y− = 0 in equilibrium P − (x− , y− ), as shown in Fig. 3. This is because isocline l2 passes through O(0, 0) as 1 − m2 = 0, and the intersection point of isoclines l1 and l2 tends to O(0, 0) as k(rY −mY )+rY mX mY

7

β → +∞. Moreover, when efficiency of species Y is extremely high (α → +∞), P − (x− , y− ) tends to the y-axis and we have limβ →+∞ x− = 0 in equilibrium P − (x− , y− ), which implies extinction of species X. Thus species Y will go to extinction since it cannot survive alone. Therefore, extremely high efficiency of species Y (α → +∞) will lead to extinction of both species. The ecological reason is that under the extremely high efficiency, species Y grows explosively and drives X into extinction through both predation and competition, which eventually leads to extinction of itself. On the other hand, when the predation β is extremely small (β < β2 ), Y cannot survive by Theorem 3.2(iiib). Hence, extremely small predation will lead to extinction of Y. A similar discussion can show that extremely low efficiency of species Y (α < α0 ) will lead to extinction of Y. The reason is that under these conditions, species Y cannot approach a reproduction success rate. A similar discussion can be given for Theorem 4.1. In Theorem 3.3(i), when the predation is extremely big (β → +∞), species X is driven into extinction by Y: (a) When α ≥ α1 , P2 is globally asymptotically stable as β → +∞, which implies extinction of X. (b) When α < α1 , P + tends to x-axis as β → +∞. Thus, separatrices of P + tends to x-axis, which implies extinction of X, as shown in Fig. 4b. The reason is that under these conditions, species X cannot approach a reproduction success rate because of the extremely big predation. Hence, extremely big predation can lead to extinction of species X. A similar discussion can show that extremely high efficiency of species Y can lead to extinction of X. On the other hand, in the situation considered by Theorem 3.3(iiib), the efficiency of species Y is low. When the predation is extremely small (β < β2 ), species Y is driven into extinction by X. The reason is that under these conditions, species Y cannot approach a reproduction success rate because of the competition from X. Since species Y can persist alone, it is the extremely small

m1=0.1, m2=1.0, α = 25 1 β = 0.2 β = 1.0 β = 3.0 β = 20

0.9

0.8

0.7

y

0.6

0.5

0.4

0.3

0.2

0.1

0

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

Fig. 3. Extinction of both IG prey X and IG predator Y as predation rate β tends to infinity in Theorem 3.2. We fix m1 = 0.1, m2 = 1.0, α = 25, and let β vary. When β increases from 0.2, 1, 3 to 20, the intersection point (stable equilibrium) decreases from (0.143,0.618), (0.0425,0.431), (0.0219,0.223), to (0.0012,0.041). Therefore, the point tends to O(0, 0) as β → ∞, which implies the extinction of X and Y.

8

Y. Wang, H. Wu / Mathematical Biosciences 259 (2014) 1–11

(b) α = 1.6, β = 1.6

(c) α = 1.6, β = 1.35

(d) α = 0.6, β = 1.35

x

x

y

y

(a) α = 5, β = 1.6

Fig. 4. Phase-plane diagrams of (3.2). Red and blue curves are the isoclines for X and Y, respectively. The black lines are separatrices of saddle points. The separatrices subdivide the plane into two regions, which are the basins of attraction of the corresponding stable equilibria. The stable and unstable equilibria are denoted by filled and open circles, respectively. Gray arrows display the direction and strength of the vector fields in the phase-plane space. We fix m1 = 0.3, m2 = 0.7 and let α , β vary. (a) When the efficiency of Y is high and predation of Y on X is big (α = 5, β = 1.6), X is driven into extinction by the predation of Y. (b) and (c) When the efficiency of Y is intermediate (α = 1.6, β = 1.6 or 1.35), the species with large density persists, while the other is driven into extinction. (d) When the efficiency of Y becomes low (α = 0.6, β = 1.35), Y is driven into extinction by the competition of X. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

predation (i.e., the competition from X) that leads to its extinction. A similar discussion can show that extremely low efficiency can lead to extinction of Y. Numerical computations confirm and extend our results. In Fig. 2, we fix m1 = 0.1, m2 = 1.2, and let α , β vary. Thus, species X can survive in the absence of Y, while Y cannot persist alone. When the efficiency of predator Y in converting the consumption into fitness is high (α = 25, β = 0.01), Y can persist at y = 0.475 as shown in Fig. 2a. When the efficiency of Y is intermediate, it can persist at y = 0.405, 0.215 if it has a large density and the predation on X is big (α = 10, β = 0.6 or 0.3), as shown in Fig. 2b–c. Otherwise, if the predation on X is small (α = 10, β = 0.1), Y goes to extinction as shown in Fig. 2d. Moreover, Fig. 2d displays that under-predation results in extinction of Y, and Fig. 2b exhibits that intermediate predation leads to the maximal density of Y, which extends the result of Theorem 3.2. In Fig. 3, we fix m1 = 0.1, m2 = 1.0, α = 25, and let β vary. When β increases from 0.2, 1, 3 to 20, the intersection point (stable equilibrium) decreases from (0.143,0.618), (0.0425,0.431), (0.0219,0.223), to (0.0012,0.041). Therefore, the point tends to O(0, 0) as β → ∞, which implies the extinction of X and Y. In Fig. 4, we fix m1 = 0.3, m2 = 0.7, and let α , β vary. Thus, either species persists in the absence of the other. When the efficiency of Y is high (α = 5, β = 1.6), X is driven into extinction by Y as shown in Fig. 4a. When the efficiency of Y is intermediate (α = 1.6, β = 1.6 or 1.35), the species with large density persists while the other is driven into extinction, as shown in Fig. 4b–c. When the efficiency of Y is low (α = 0.6, β = 1.35), Y is driven into extinction by X as shown in Fig. 4d.

In Fig. 5, we fix m1 = 0.3, m2 = 0.7, and let α , β vary. Thus, either species persists in the absence of the other. When the efficiency of Y is high and the predation on X is big (α = 4, β = 1), the species coexist at (0.025,0.338) and Y can approach a density larger than its carrying capacity, as shown in Fig. 5a. When the efficiency of Y becomes intermediate (α = 1.6, β = 1 or 0.72), it can persist at (0.11, 0.295) if it has a large density. Otherwise, Y is driven into extinction by X as shown in Fig. 5b–c. When the predation on X becomes small (α = 1.6, β = 0.6), Y goes to extinction as shown in Fig. 5d. The present lattice gas model is different from Levins-type patch occupancy model, which is shown as follows. In a lattice gas model, it is assumed that one individual can only occupy one site while one site can be occupied by only one individual. Thus, when an individual is killed, there will exist an empty site, which can be occupied by either species in a random way. However, in a Levins-type patch occupancy model (e.g., Ref. [25]), when a patch of species 2 is excluded by species 1, the patch will belong to species 1 even when there is no individual of species 1 to occupy it. On the other hand, one patch can be occupied by two species in the Levins-type patch occupancy model. For IGP systems with spatial competition, the result in this work is different from those in previous papers. Previous works showed properties of the IGP by local stability analysis and numerical computations, while this work theoretically exhibits mechanisms by which IG prey and IG predator can coexist. That is, varying one parameter or population density can lead to transitions of interaction outcomes in a smooth fashion. For example, in the situation where IG predator cannot survive alone in Theorem 3.2, (a) When the efficiency of IG predator is large (α > α1 ), it can coexist with IG prey. (b) When

Y. Wang, H. Wu / Mathematical Biosciences 259 (2014) 1–11

9

(b) α = 1.6, β = 1

(c) α = 1.6, β = 0.72

(d) α = 1.6, β = 0.6

x

x

y

y

(a) α = 4, β = 1

Fig. 5. Phase-plane diagrams of (3.2). Red and blue curves are the isoclines for X and Y, respectively. The black lines are separatrices of saddle points. The separatrices subdivide the plane into two regions, which are the basins of attraction of the corresponding stable equilibria. The stable and unstable equilibria are denoted by filled and open circles, respectively. Gray arrows display the direction and strength of the vector fields in the phase-plane space. We fix m1 = 0.3, m2 = 0.7 and let α , β vary. (a) When the efficiency of Y is high and predation of Y on X is big (α = 4, β = 1), Y can approach a density larger than its carrying capacity. (b) and (c) When the efficiency of Y is intermediate (α = 1.6, β = 1 or 0.72), the species with large density persists, while the other is driven into extinction. (d) When predation of Y on X becomes small (α = 1.6, β = 0.6), Y is driven into extinction by the competition of X. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the efficiency is intermediate (α0 < α < α1 ), IG predator can persist if its predation is big (β > β2 ) and its initial density is high (above the stable manifold of P + in Fig. 2b). Otherwise, IG predator goes to extinction. (c) When the efficiency becomes small (α < α0 ), IG predator cannot survive. (d) Intermediate predation (β = β¯ ) can lead to the maximal density of IG predator, while extremely big predation can lead to extinction of both species. A similar discussion can be given for other situations in Theorems 3.3 and 4.1. In the present model, we use linear functional response to describe interaction between IG prey and IG predator. When other functional responses (e.g., Holling Type functional responses) are used, more complex dynamics such as periodic oscillation may occur. Thus we will pursue it in a future work. Acknowledgments We would like to thank D.L. DeAngelis for providing useful references and comments, and the three anonymous reviewers for their helpful comments on the manuscript. This work was supported by NSF of Guangdong S2012010010320, 1414050000636 and STF of Guangzhou 1563000413. Appendix A. Proof of Theorem 3.2 (i) Since α > α1 , P1 is below l2 and is a saddle point. It follows from the convexity of l2 that there is a unique interior equilibrium P − of (3.2), which is globally asymptotically stable in intR2+ by Theorem 2.1.

(ii) Since α < α1 , P1 is above l2 and C1 > 0. Thus P1 is asymptotically stable. From α > α0 we have α > 1. By (3.6),  = 0 can be rewritten as

A¯ β 2 + B¯ β + C¯ = 0. It follows from convexity of l2 that when β > β2 , we have  > 0, and l1 and l2 have two intersection points in the first quadrant, which implies that system (3.2) has two interior equilibria P − and P + in (3.6). By Lemma 3.1, P− is asymptotically stable while P + is a saddle point. By Theorem 2.1, the separatrices of the saddle point divide the first quadrant into two regions, one is the basin of attraction of P1 while the other is that of P − . When β = β2 , we have  = 0 and equilibria P − and P + coincide. Since the Jacobian matrix A of (3.2) at P + satisfies trA < 0 and det A = 0 as  = 0, there is a simple eigenvalue λ = 0 of A as shown in Lemma 3.1. By (3.3), the matrix A has a right eigenvector v and left eigenvector w with

v = (a12 , −a11 )T ,

w = (a21 , −a11 )T ,

where

a11 = −r1 x|P+ ,

a12 = −r1 (1 + β)x|P+ ,

a21 = r2 y[α(1 − x − y) − 1 − α x]|P+ , a22 = −r2 y(1 + α x)|P+ . A direct computation (1 + β)|P+ < 0.

shows

that

a21 = −r2 y(1 + α x)/

10

Y. Wang, H. Wu / Mathematical Biosciences 259 (2014) 1–11

Let F = (f1 , f2 )T be the righthand side of (3.2). From the definition of D2 F (x0 ) [15, page 69], we have

∂ 2 F (P+ ) ∂ 2 F (P+ ) − 2a11 a12 2 ∂ x∂ y ∂x 2 + ∂ F ( P ) . + a211 ∂ y2

D2 F (P + )(v, v) = a212

equilibrium P1 of (3.2) is globally asymptotically stable in the first and fourth quadrants. However, when there is a small variation of β near β1 , degenerate transcritical bifurcation occurs: If β < β1 , there are two equilibria: P1 becomes unstable, P − is in the first quadrant and is asymptotically stable. If β > β1 , there are two equilibria: P1 becomes unstable, P + is in the fourth quadrant and is asymptotically stable.

A direct computation shows that

  −r1 xy , 0 P+   0 [D2 F (P + )(v, v)] = , −2r1 r2 a11 α(1 + β)xy P+ Fβ (P + ) =

and

wT Fβ = −a21 r1 xy|P+ = 0, wT [D2 F (P + )(v, v)] = 2r1 r2 a211 α(1 + β)xy|P+ = 0. Thus, at β = β2 , Sotomayor’s theorem implies the existence of a saddle-node bifurcation for (3.2) [15, p. 338]. In other cases, there is no interior equilibrium of (3.2). It follows from Theorem 2.1 that all positive solutions of (3.2) converge to equilibrium P1 . (iii) Since α = α1 , P1 is on l2 . The slopes of l1 and l2 at P1 are

k1 = −

1 , 1+β

k2 =

(2m1 − 1)m2 − m21 . m2 (1 − m1 )

Thus when k2 < 0 and |k1 | < |k2 |, there is a unique interior equilibrium P − of (3.2), which is globally asymptotically stable. Here, k2 < 0 can be written as 2m1 − 1 ≤ 0, or 2m1 − 1 > 0, m2 < m21 /(2m1 − 1), while |k1 | < |k2 | can be written as β > β1 . In other cases, there is no interior equilibrium of (3.2). Thus all solutions of (3.2) converge to P1 by Theorem 2.1. System (3.2) experiences a pitchfork bifurcation at equilibrium P1 as α = α1 . Indeed, when β = β1 , equilibria P − and P + coincide at P1 . Since the Jacobian matrix A of (3.2) at P1 satisfies trA = −r1 (1 − m1 ) < 0 and det A = 0, there is a simple eigenvalue λ = 0 of A. By (3.3), the matrix A has a right eigenvector v and left eigenvector w with

v = (a12 , −a11 )T ,

w = (a21 , −a11 )T ,

where

a11 = −r1 (1 − m1 ), a21 = 0,

a12 = −r1 (1 + β)(1 − m1 ),

a22 = 0.

Let F = (f1 , f2 )T be the righthand side of (3.2). From the definition of DFα (x0 ), D2 F (x0 ), D3 F (x0 ) [15, page 69], a direct computation shows that



Fα =



0

−r2 xy(1 − x − y)

wT Fα (P1 ) = 0,

, P1

wT [DFα (P1 )v] = r2 m1 (1 − m1 )a211 = 0,

wT [D2 F (P1 )(v, v)] = 0, and

wT [D3 F (P1 )(v, v, v)] = −2r1 r2 α1 (2 + β1 )(1 − m1 )a211 a12 = 0. Thus, at α = α1 , Sotomayor’s theorem implies the existence of a pitchfork bifurcation for (3.2) [15, p. 338)]. That is, when α ≤ α1 , there is a unique equilibrium P1 of (3.2), which is globally asymptotically stable. when α > α1 , there are three equilibria: P1 becomes unstable, P − and P + are asymptotically stable while P + is in the fourth quadrant. System (3.2) also experiences a degenerate transcritical bifurcation at equilibrium P1 as β = β1 . Indeed, when β = β1 ,

Appendix B. Proof of Theorem 3.3 (i) Since β > β3 , P2 is above l1 and C2 > 0. Thus P2 is asymptotically stable. If α ≥ α1 , then P1 is below l2 and is a saddle point. It follows from the convexity of l2 that there is no interior equilibrium of (3.2), which implies that P2 is globally asymptotically stable in intR2+ by Theorem 2.1. If α < α1 , then P1 is above l2 and is asymptotically stable. It follows from the convexity of l2 that there is a unique interior equilibrium P + of (3.2). By Lemma 3.1, P + is a saddle point. By Theorem 2.1, the separatrices of P + divide the first quadrant into two regions, one is the basin of attraction of P1 while the other is that of P2 . (ii) Since β = β3 , P2 is on l1 . When α ≥ α1 , there is no interior equilibrium of (3.2). Thus all solutions of (3.2) converge to P2 by Theorem 2.1. When α < α1 , P1 is above l2 . The slopes of l1 and l2 at P2 are

k¯ 1 = −

1 , 1+β

k¯ 2 = α m2 − 1.

If α > α2 , then k¯ 2 > k¯ 1 and there is a unique interior equilibrium P + of (3.2), which is a saddle point. By Theorem 2.1, the separatrices of P + divide the first quadrant into two regions, one is the basin of attraction of P1 while the other is that of P2 . If α ≤ α2 , then k¯ 2 ≤ k¯ 1 and there is no interior equilibrium of (3.2). Thus all solutions of (3.2) converge to P1 by Theorem 2.1. (iii) Since β < β3 , P2 is below l1 and C2 < 0. Thus P2 is a saddle point. When α ≥ α1 , P1 is below l2 and is a saddle point. It follows from the convexity of l2 that there is a unique interior equilibrium P − , which is globally asymptotically stable by Theorem 2.1. When α < α1 , P1 is above l2 and is asymptotically stable. By a proof similar to that of Theorem 3.2(ii) that system (3.2) experiences a saddle-node bifurcation at β = β2 . Thus the result of Theorem 3.3(iii) is proved. References [1] P. Amarasekare, Spatial dynamics of communities with intraguild predation: the role of dispersal strategies, Am. Nat. 170 (2007) 819–831. [2] P. Amarasekare, Coexistence of intraguild predators and prey in resource-rich environments, Ecology, 89 (2008) 2786–2797. [3] J.M. Fedriani, T.K. Fuller, R.M. Sauvajot, E.C. York, Competition and intraguild predation among three sympatric carnivores, Oecologia 125 (2000) 258–270 [4] T.E. Harris, Contact interaction on a lattice, Ann. Probab. 2 (1974) 969–988. [5] R.D. Holt G.A. Polis, A theoretical framework for intraguild predation, Am. Nat. 149 (1997) 745–764. [6] S. Iwata, K. Kobayashi, S. Higa, J. Yoshimura, K. Tainaka, A simple population theory for mutualism by the use of lattice gas model, Ecol. Model. 222 (2011) 2042–2048. [7] Y. Kang, L. Wedekin, Dynamics of a intraguild predation model with generalist or specialist predator, J. Math. Biol. 67 (2013) 1227–1259. [8] T. Kawai, Y. Tadokoro, K. Tainaka, T. Hayashi, J. Yoshimura, A lattice model of fashion propagation with correlation analysis, Int. J. Syst. Sci. 39 (2008) 947–957. [9] H. Matsuda, N. Ogita, A. Sasaki, K. Sato, Statistical mechanics of population: the lattice Lotka-Volterra model, Prog. Theor. Phys. 88 (1992) 1035–1049. [10] S.D. Mylius, K. Klumpers, A.M. de Roos, L. Persson, Impact of intraguild predation and stage structure on simple communities along a productivity gradient, Am. Nat. 158 (2001) 259–276. [11] N. Nakagiri, K. Tainaka, T. Tao, Indirect relation between extinction and habitat destruction, Ecol. Model. 137 (2001) 109–118. [12] C. Neuhauser, Ergodic theorems for the multitype contact process, Probab. Theory Relat. Fields 91 (1992) 467–506.

Y. Wang, H. Wu / Mathematical Biosciences 259 (2014) 1–11 [13] S. Nishijima, G. Takimoto, T. Miyashita, Roles of alternative prey for mesopredators on trophic cascades in intraguild predation systems: a theoretical perspective, Am. Nat. 183 (2014) 625–637. [14] F. Palomares, T.M. Claro, Interspecific killing among mammalian carnivores, Am. Nat. 153 (1999) 492–508. [15] L. Perko, Differential Equations and Dynamical Systems, Springer-Verlag, New York, 2001. [16] G.A. Polis, C.A. Myers, R.D. Holt, The ecology and evolution of intraguild predation: potential competitors that eat each other, Ann. Rev. Ecol. Syst. 20 (1989) 297–330. [17] G.A. Polis, R.D. Holt, Intraguild predation: the dynamics of complex trophic interactions, Trends Ecol. Evol. 7 (1992) 151–154. [18] E. Ruggieri, S.J. Schreiber, The dynamics of the Schoener-Polis-Holt model of intraguild predation, Math. Biosci. Eng. 2 (2005) 279–288. [19] T.W. Schoener, Competition and the form of habitat shift, Theor. Popul. Biol. 6 (1974) 265–307. [20] T.W. Schoener, Alternatives to Lotka-Volterra competition: models of intermediate complexity, Theor. Popul. Biol. 10(1976) 309–333.

11

[21] K. Tainaka, Stationary pattern of vortices or strings in biological systems: lattice version of the Lotka-Volterra model, Phys. Rev. Lett. 63 (1989) 2688–2691. [22] G. Takimoto, T. Miki, M. Kagami, Intraguild predation promotes complex alternative states along a productivity gradient, Theor. Popul. Biol. 72 (2007) 264–273 [23] G. Takimoto, T. Miki, M. Kagami, Corrigendum to “Intraguild predation promotes complex alternative states along a productivity gradient”, Theor. Popul. Biol. 73 (2008) 460. [24] G. Takimoto, D.M. Post, D.A. Spiller, R.D. Holt, Effects of productivity, disturbance, and ecosystem size on food-chain length: insights from a metacommunity model of intraguild predation, Ecol. Res. 27 (2012) 481–493. [25] D.E. Taneyhill, Metapopulation dynamics of multiple species: the geometry of competition in a fragmented habitat, Ecol. Monogr. 70 (2000) 495–516. [26] Y. Wang, Dynamics of a lattice gas system, J. Math. Anal. Appl. 419 (2014) 1303–1321. [27] R. Hall, Intraguild predation in the presence of a shared natural enemy, Ecology 92 (2011) 352–361. [28] T. Okuyama, Intraguild predation with spatially structured interactions, Basic and Applied Ecology 9 (2008) 135–144.