Extinction and coexistence of species for a diffusive intraguild predation model with B-D functional response

Extinction and coexistence of species are two fundamental issues in systems with IGP. In this paper, we constructed a mathematical model with IGP by introducing heterogeneous environment and B-D functional response between the predator and prey. First, some sufficient conditions for the extinction and permanence of the time-dependent system were obtained by using comparison principle and upper and lower solution method. Second, we got some necessary and sufficient conditions for the existence of coexistence states by means of the fixed point index theory. In addition, we discussed the uniqueness and stability of coexistence state under some conditions. Finally, we studied the effects of the parameters in system on the spatial distribution of species and obtained some interesting results about the extinction and coexistence of species by using numerical simulations.


1.
Introduction. Intraguild predation, or IGP, is a widespread ecological phenomenon in natural communities, which occurs when two consumers share a common resource and also engage in a predator-prey interaction [35,36,4]. One of the pioneering work to rigorously model intraguild predation was by Holt and Polis in [18], where a basic ODE models were developed for a three species community. It was showed in [18] that IGP significantly influences the distribution, abundance and coexistence of many species.
Following the model in [18], Tanabe and Namba in [41] and Abrams and Fung in [31] mainly obtained some results about the existence of chaos by numerical simulations; Hsu, Ruan and Yang in [19] completely classify the parameter space of model in article [18] to obtain some results about extinction and uniform persistence. On the other hand, a growing number of biological and mathematical models including IGP have been proposed by incorporating some more realistic ecological factors, such as delay [39,10], age or stage structure [45,38,21], adaptive foraging [32,26], defense [20,42,44,23], refuge [29] and additional species [27,17]. The main goal in the studies above is to ascertain mechanisms for extinction and coexistence of species in systems with IGP. We remark that different functional responses between predator and prey have been assumed in models including IGP, such as linear functional response, Holling type II, Holling-Type III and ratio-dependent functional response respectively [1,43,22,15]. We point out that, to our knowledge, few model is constructed with assuming the B-D functional response between the predator and prey to investigate the IGP.
We remark that all the models aforementioned assumed a homogeneous environment with no migration dynamics. As pointed in [18], spatial heterogeneity may have a important effect on the coexistence and extinction of the species in the IGP communities. Amarasekare in [2,3] first model IGP in a heterogeneous environment by using an environment consisting of 3 distinct patches and investigated the effects of different dispersal strategies on the coexistence based on numerical simulations. On the other hand, [37] model IGP in a heterogeneous environment with space to be a continuous variable resulting in a parabolic system of PDEs with Neumann boundary condition. They assumed that the IGprey employs a fitness based avoidance strategy and proved the existence of a global attractor for this system and derived conditions for the uniform persistence of the IGprey.
Motivated by the articles above, we construct a mathematical model with heterogeneous environment by introducing the diffusion and B-D functional response between the predator and prey and obtain the following IGP system: in Ω × (0, ∞), (u(0, x), v(0, x), w(0, x)) = (u 0 (x), v 0 (x), w 0 (x)) ≥ (0, 0, 0) in Ω, (1) where Ω is a bounded domain in R N (N ≥ 1 is an integer) with a smooth boundary ∂Ω and ν is the outward unit rector on ∂Ω. Here we assumed a Robin boundary conditions for the system which means that that the flux of individuals across the boundary point is proportional to the density. One can find more details about the biological meaning of the Robin boundary condition in page 31 in [9]. The three functions u, v and w represent the densities of the resource, intermediate predator (IG prey) and top predator (IG predator) respectively. The terms a1u 1+b1u+e1v , a2u 1+b2u+e2w and a3v 1+b3v+e3w represent the B-D functional response, which turn to Holling-II functional response if e i = 0 and linear functional response if both b i = 0 and e i = 0. Here b i /a i represent the handing time of the predators and e i describes that mutual interference among individual of predators [5,14]. The positive constants a i (i = 1, 2, 3) are the consumption rate and m i (i = 1, 2, 3) represent the conversion rates of the prey to predator and r 1 is the growth rate of resource. We assume that r 3 may be positive or negative. If r 3 > 0, it represents the growth rate of the IGpredator, which implies that the IGpredator is a generalist with a exclusive resources; if r 3 < 0, the death rate of the IGpredator is r 3 − w, which the IGpredator has no exclusive resources. Similarly, we also assume that parameter r 2 may be positive or negative and it represents the death rate of the IGprey if r 2 > 0 and the immigration rate if r 2 < 0. We point out that the model (1) turn to food chain system if a 2 = m 2 = 0 which was studied in [47]; while if a 3 = m 3 = 0, model (1) turn to a system with two predators competing for one prey, which was investigated in [46].
In our work here, one of the main purposes is to study the existence of positive stationary solutions of (1) by using fixed point index theory, which are the positive solutions of in Ω, The positive solutions of (2) are also called coexistence states. The rest of this paper is organized as follows. In Section 2, some sufficient conditions for the extinction and permanence of species in system (1) are obtained by using upper and lower solution method. In Section 3, the existence and nonexistence of coexistence states of model (2) are studied by using some degree theorems developed. In Section 4, we further discuss the uniqueness and stability of coexistence state under some conditions. In Section 5, some interesting results about the extinction and coexistence of species are obtained by using numerical simulations. Finally in Section 6, we give some discussions about the conclusions obtained above.
2. Extinction and permanence. In this section, we would like to investigate the asymptotic behavior of the time-dependent solutions of system (1). We first give some definitions and theorems which will be used in the following.

Extinction.
In this subsection, we study the extinction of species in system (1). We can obtain the following sufficient conditions for the extinction of species in system (1).
be a positive solution of (1).

2.3.
Permanence and global attractor. In this subsection, we would like to investigate the permanence and global attractor of the time dependent system (1) by using upper and lower solution method. We first give the following definition about ordered upper and lower solution of (2). Definition 2.4. A pair of functions (u, v, w) and (u, v, w) in C(Ω) ∩ C 2 (Ω) are called ordered upper and lower solution of (2) if they satisfy the relation u ≥ u, v ≥ v, w ≥ w and the following inequalities: in Ω, Let v be the unique positive solution of the following problem in Ω, and v be the unique positive solution of the following problem Let w be the unique positive solution of the following problem in Ω, and w be the unique positive solution of the following problem in Ω, In order to give the main results, we introduce the following assumptions: Remark 2.5. By Theorem 2.2, it is easy to see that the existence and uniqueness of v , v , w , w follow from the assumptions (19).
The following theorem provides sufficient conditions for permanence of the timedependent system (1).
Theorem 2.6. Assume that the conditions in (19) hold. Then, there exist a pair of functions (u + , v + , w + ) and and satisfy the following relations Θ k1 ( is a positive global attractor of (1).
Proof. It is easy to show that (Θ k1 (r 1 ), v , w ) and (Θ k1 (r 1 − a1 e1 − a2 e2 ), v , w ) are a pair of ordered upper and lower solution of (1) under assumption (19). From the definition 8.1 in page 425 in [33], we know that the reaction functions of the system (2) are quasimonotone. Then the existence of (u + , v + , w + ) and (u − , v − , w − ) can be proved easily by using the iteration scheme (10.3) or Theorem 10.1 in page 438 in [33].
Next, we prove that is a positive global attractor of (1). Note that (u − , v − , w − ) is positive in Ω by maximum principle, it suffice to prove that Let be sufficiently small such that It can be seen that the assumption (19) implies the positivity of the middle terms in (A) and (B). The 2nd inequalities in (A) and (B) follow from the monotonicity of the expression. So we may choose small enough such that the inequalities in (A) and (B) hold. First, we can obtain from the first equation of (1) that and Then by using Theorem 2.2, comparison principle and assumption (19), we can get from (20) and (21) and Then, it follows from (22), (23) and the second equation of (1) that
The next theorem gives sufficient conditions for a global attractor in the case that exactly one species is dying out. This can be proved similarly as in the above theorem. So we omit the proof.
is a global attractor of (1).
, then there exists a pair of quasisolution (u + , w + ) and (u − , w − ) of the following system in Ω, is a global attractor of (1).
3. Coexistence states. In this section, we shall investigate the existence of coexistence states for (2) by using the fixed point index theory. We first give some theorems to calculate the indexes at the trivial and semi-trivial states of (2).
Theorem 3.1. (see [16,28].) Assume h(x) ∈ C α (Ω)(0 < α < 1) and M is a sufficiently large number such that M > h(x) for all x ∈ Ω. Define a positive and compact operator L : From Theorem 3.1, we can see that it is crucial to know the sign of the eigenvalue λ 1,k (h) to determine the spectral radius of L. The following theorem is established to determine the sign of the principle eigenvalue λ 1,k (h).
Theorem 3.2. (see [7,8] Now, we state the fixed point index theory which plays an crucial role in getting the sufficient conditions for the existence of positive solution of model (2).
Let E be a real Banach space and W ⊂ E be the natural positive cone of E. For y ∈ W, define W y = {x ∈ E : y + η ∈ W for some η > 0} and S y = {x ∈ W y : −x ∈ W y }. Then W y is a wedge containing W, y, −y, while S y is a closed subset of E containing y. Let T be a compact linear operator on E which satisfies T (W y ) ⊂ W y . We say that T has property α on W y if there is a t ∈ (0, 1) and an ω ∈ W y \S y such that (I − tT )ω ∈ S y . Let A : W → W is a compact operator with a fixed point y ∈ W and A is a Fréchet differentiable at y. Let L = A (y) be the Fréchet derivative of A at y. Then L maps W y into itself. We denote by where S only contains discrete points. Then the following theorem can be obtained. [47,46,24] where σ is the sum of algebraic multiplicities of the eigenvalues of L which are greater than 1.
Finally, we introduce the following theorem about degree calculations, which was introduced by E. N. Dancer and Y. H. Du in [13] and we state here for convenience.
Assume that E 1 and E 2 are ordered Banach spaces with positive cones W 1 and W 2 , respectively. Let E = E 1 ⊕ E 2 and W = W 1 ⊕ W 2 . Then E is an ordered Banach space with positive cone W. Let D be an open set in W containing 0 and Then the following theorem can be obtained.
(i) If for any u ∈ S, the spectral radius γ( 2 (u, 0)) | W2 > 1 and 1 is not an eigenvalue of 2 (u, 0) | W2 corresponding to a positive eigenvector, then deg Nonexistence of coexistence states. We start from the following lemma which give the priori bounds for the coexistence states of (2). Lemma 3.5. If m 1 r 1 > r 2 (1 + b 1 r 1 ), then any coexistence state (u, v, w) of (2) has an a priori bounds:

EXTINCTION AND COEXISTENCE OF SPECIES 3767
Proof. From the first equation of (2), we can obtain Then by (32), the maximum principle and Theorem 2.2, we know From the second equation of (2), we have Hence by (33) Then from the third equation of (2), we have By using (34), the maximum principle and Theorem 2.2, we have The proof is completed.
Remark 3.6. From the proof of Lemma 3.5, we can see that m 1 r 1 > r 2 (1 + r 1 ) is the necessary condition for problem (2) have coexistence states. So, throughout this subsection, we assume that: (H) m 1 r 1 > r 2 (1 + r 1 ).
Now we can obtain the following necessary conditions for the existence of coexistence states for system (2).
are necessary conditions for the existence of positive solutions to (2).
Proof. (i) Assume (u, v, w) is a positive solution of (2). Then it follows from the first equation and the comparison principle of eigenvalues that which is a contradiction.
(ii) Let (u, v, w) be a positive solution of (2). Then r 1 > λ 1,k1 from (i) and u(x) ≤ Θ k1 (r 1 ) ≤ r 1 , v < K 2 , w < K 3 from Lemma 3.5. We can get from the second equation and the comparison principle of eigenvalues that Similarly, it follows from the third equation that be a positive solution of (2). Then r 1 > λ 1,k1 from (i) and u(x) ≤ Θ k1 (r 1 ) ≤ r 1 . Since r 3 > λ 1,k3 , we know from the third equation of (2) that w ≥ Θ k3 (r 3 ). We can get from the first equation and the comparison principle of eigenvalues that and from the second equation that .
The proof is completed.
3.3. Existence of coexistence states. In the rest of this section, we shall give sufficient conditions for (2) has coexistence states by using the fixed point index theory. Now, we introduce the following notations: From Lemma 3.5, we can see that the coexistence state of (2) must be in D. Take q sufficiently large such that u( 1+b3v+e3w − r 2 v + qv and r 3 w − w 2 + m2uw 1+b2u+e3w + m3uw 1+b3u+e3w + qw are respectively monotone increasing with respect to u, v and w for all (u, v, w) ∈ [0,

EXTINCTION AND COEXISTENCE OF SPECIES 3769
Define a positive and compact operator : Remark 3.8. Observe that (2) is equivalent to (u, v, w) = (u, v, w), and then it is sufficient to prove that has a non-constant positive fixed point in D to show that (2) has a coexistence state.
From the remark above, we can see that it is necessary to calculate the degree of I − in D relative to W and the fixed point index of at (0, 0, 0) relative to W. The following lemma give the corresponding results about deg W (I − , D) and index W ( , (0, 0, 0)).
To study the other semi-trivial solutions of (2), consider the following three possible subsystems: in Ω, and in Ω, It is easy to show that system (43) has only trivial solution (v, w) = (0, 0) and semi-trivial solution (v, w) = (0, Θ k3 (r 3 )). About the positive solutions of systems (41) and (42), we can obtain some results from [6,30] and [16] respectively and simple comparison arguments for elliptic problems. We point out that the corresponding main results are still valid under Robin boundary conditions even if the results in the above references were obtained under Dirichlet boundaries.
Lemma 3.14. Assume that r 1 > λ 1,k1 and −r 2 > λ 1,k2 − Proof. Recall the definitions of E, W and D in section 3, we define Then we are in the framework to use Theorem 3.4.
Similarly, we can obtain the following lemma by using the similar methods as above.
Lemma 3.15. Assume that one of the following conditions hold: and r 3 > λ 1,k3 .
Based on above analysis, we can give the following results about the existence of coexistence states of (2). Theorem 3.16. If one of the following conditions hold, then (2) has at least one coexistence state.
The proof is completed. 4. Stability and uniqueness of coexistence states. In this section, we investigate the stability and uniqueness of coexistence states of (2). We first introduce some notations to give our main results.
Denote v † be the unique positive solution of the following system in Ω, 1+b1Θ k 1 (r1) . Let w † denote the unique positive solution of the following system in Ω, We remark that the existence and uniqueness of positive solution of (53), (54) and (55) follows from Theorem 4. We can also note that the condition (ii) in Theorem 3.16 implies the existence and uniqueness of positive solution of (53), (54) and (55).
The following theorem shows the uniqueness, non-degeneracy and linear stability of the coexistence states for (2) under some assumptions. (i) The coexistence state of (2) converge to (u † , v † , w † ) as a i → 0 for i = 1, 2, 3.
(ii) There exists a positive constantã =ã (a 1 , a 2 , a 3 ) such that for a 1 , a 2 , a 3 <ã, (2) has exactly one coexistence state which is non-degenerate and linearly stable.
Proof. (i) It is easy to see that the compact operator (u, v, w) defined in Section 3 converges to the operator as a i → 0 for i = 1, 2, 3. So the fixed points of (2) converge to the fixed points of (u, v, w). It can be shown that (u † , v † , w † ) is the only fixed point of (u, v, w), the conclusion follows.
If η = 0, we have which is a contradiction. Hence we have ξ = ζ = η = 0. This contradiction indicates that the coexistence state of (2) is non-degenerate and linearly stable. Next, we prove the uniqueness of coexistence state of (2). First the compactness implies that has at most finitely many positive fixed points in the region N defined in Section 4.3 and denote them by (ui, vi, wi) for i = 1, ..., k. Then for sufficiently small a1, a2 and a3, as done in the proof of Lemma 3.9, it can be shown that I − (ui, vi, wi) is invertible and (ui, vi, wi) has no property α on W (u i ,v i ,w i ) . In addition, (ui, vi, wi) does not have a real eigenvalue which is greater than or equal to 1. Thus it follows from Theorem 3.3 that index( , (ui, vi, wi)) = (−1) 0 = 1 for i = 1, ..., k. Then by using the additivity property of the degree, we have The uniqueness is obtained.
(iii) It is easy to see that the compact operator (u, v, w) defined in Section 4.3 converges to the operator as ei → ∞ for i = 1, 2, 3. Hence the fixed points of (2) converge to the fixed points of (u, v, w). It can be shown that (Θ k 1 (r1), 0, Θ k 3 (r3)) is the only fixed point of (u, v, w) as in [24] and the conclusion follows.
(iv) As in the proof of (ii), we can derive a contradiction to prove that the coexistence state is non-degenerate and linearly stable and get the uniqueness by using the additivity property of the degree. So we omit the proof.
5. Numerical simulation. In this section, some numerical results on spatial patterns to the system (1) are provided to verify and complement the analytic results in section 4. For simplicity, the simulations are performed on a one dimensional spatial grid with length of L under Dirichlet boundary condition.
(58) The initial conditions consist of random perturbations about the homogeneous state of u 0 = 1, v 0 = 1, w 0 = 0.5. Here we point out that the characteristics of stationary state remain unchanged even the initial conditions are displace by the other forms such as u 0 = |sinx|, v 0 = |cosx|, w 0 = 0.5. The system (1) will be solved with the approach of discretizing the diffusion operator in space, applying forward Euler integration to the finite-difference equations with time steps ∆t = 0.002. We run the simulations until they reach a stationary state or until they show a behavior that does not seem to change its characteristics anymore. 5.1. Effect of growth rate. Here we first investigate the effects of r 1 , the growth rate of resource, on distribution of species. To proceed, we vary the values of r 1 with the other parameters fixed. Fig.1 illustrates the spatial distribution of resource, IGprey and IGpredator in a one dimension region with different levels of growth rates r 1 for resource. One can see that, the population IGprey is in the state of extinction when the growth rate of resource is low ( Fig.1(a)). With the increase of growth rate (bigger r 1 ), the population IGprey can coexist with resource and IGpredator ( Fig.1(b)). One can also note that larger growth rate of resource leads to larger mass of IGprey and IGpredator ( Fig.1(c)). Fig.2 illustrates the spatial  distribution of resource, IGprey and IGpredator in a one dimension region with different levels of r 2 . Here we recall that the parameter r 2 describes the death rate of IGprey if r 2 > 0 and the immigration rate of IGprey if r 2 < 0. One can see that, the population IGprey is driven to extinction when the death rate r 2 is high (Fig.2(a)). With the decrease of death rate of IGprey (smaller r 2 ), the population IGprey can coexist with resource and IGpredator ( Fig.2(b)). When there is immigration of IGprey, that is r 2 < 0, one can note that larger mass of IG prey can coexist with resource and IGpredator (Fig.2(c)). Fig.3 illustrates the  spatial distribution of resource, IGprey and IGpredator in a one dimension region with different levels of r 3 . Here we recall that the parameter r 3 describes exclusive resource of IGpredator if r 3 > 0 and the death rate of IGpredator if r 3 < 0. When there is no exclusive resource of IGpredator, that is r 3 < 0 and describe death rate, one can see that the population IGpredator is driven to extinction when the death rate is high (Fig.3(a)). With the decrease of death rate of IGprey (r 3 < 0 with smaller absolute value), the population IGpredator can coexist with source and IGpredator (Fig.3(b)). When there is exclusive resources for IGpredator, that is r 3 > 0, one can note that the IGprey is driven to extinction by the growth of IGpredator (Fig.3(c)).
One can conclude from Fig.1 that, under some circumstances, the increase of the growth rate r 1 of resource mediate the coexistence of species in system, while it follows from Fig.2 that the decrease of death rate r 2 of IGprey mediate the coexistence. The effects of the parameter r 3 can be complex and it leads to the extinction of IGpredator if it is small enough, leads to the extinction of IGprey if it is large enough and to coexistence of species when it lies in the mediate space.

Effect of interference of predators.
Here we first investigate the effects of e 2 , which describes that mutual interference among individual of IGpredators [29,30], on distribution of species. To proceed, we vary the values of e 2 with the other parameters fixed. Fig.4 illustrates the spatial distribution of source, IGprey and IGpredator in a one dimension region with different interference levels e 2 among IGpredators. One can see that, the population IGprey is in the state of extinction when interference level of predators is very low (Fig.4(a)). With the increase of interference levels among IGpredators (bigger e 2 ), the population IGprey can coexist with resource and IGpredator ( Fig.4(b)). One can also note that stronger interference among predators leads to larger mass of IGprey and smaller mass of IGpredator (Fig.4(c)). Similarly, Fig.5 illustrates the spatial distribution of resource, IGprey and IGpredator in a one dimension region with different interference levels e 3 among IGpredators. Fig.6 illustrates the spatial distribution of source, IGprey and IGpredator in a one dimension region with different interference levels among IGprey. One can see that, the population IGprey can coexist with resource and IGpredator when interference level among IGprey is very low (Fig.6(a)). With the     increase of interference level among IGprey(bigger e 1 ), the mass of population IGprey decrease (Fig.6(b)) and turn to extinction at last (Fig.6(c)). Fig.7 illustrates the spatial distribution of source, IGprey and IGpredator in a one dimension region with different interference levels e 2 and e 3 among predators. One can see that, the population IGpredator can coexist with resource and IGprey when interference levels among IGpredators is very low (Fig.7(a)). With the increase of interference levels of IGpredators (bigger e 2 and e 3 ), the mass of population IGpredator decreases ( Fig.7(b)) and turn to extinction at last (Fig.7(c)).  One can conclude from Fig.4 and Fig.5 that different interference levels e 2 and e 3 among predators can lead to coexistence of species under some circumstances; while from Fig.6 and Fig.7, we see that they can also lead to extinction of species under some other circumstances.

Effect of intraguild predation.
Here we investigate the effects of a 2 and m 2 , which describes Intraguild predation in a food chain model, on distribution of species. To proceed, we vary the values of a 2 and m 2 with the other parameters fixed. Fig.8 illustrates the spatial distribution of resource, IGprey and IGpreda-  tor in a one dimension region with different intraguild predation strengths a 2 for IGpredator. One note that the parameter r 3 is supposed to be negative, which means the IGpredator has no exclusive resource and is a specialist predator. One can see that, when there is no IGpredation in system, that is a 2 = 0, the population IGpredator is in the state of extinction ( Fig.8(a)). With the introduction of IGpredation (bigger a 2 ), the population IGpredator turn to coexist with resource and IGprey (Fig.8(b)). One can also note that stronger IGpredation leads to larger mass of IGpredator and smaller mass of IGprey (Fig.8(c)). Fig.9 illustrates the spatial distribution of resource, IGprey and IGpredator in a one dimension region with different intraguild predation strengths a 2 for IGpredator. Here we point out the parameters used in Fig.9 are the same with Fig.8 except the interference levels e 1 , e 2 and e 3 . One can see that, when there is no IGpredation, that is a 2 = 0, the population IGprey coexist with resource and IGpredator ( Fig.9(a)). With the introduction of IGpredation (bigger a 2 ), the mass of the population IGprey decrease ( Fig.9(b)) and turn to the state of extinction at last (Fig.9(c)).
From Fig.8 and Fig.9 above, we can see that the introduction of IGpredation can lead to coexistence or extinction of species under some circumstances. 6. Conclusions and discussions. In recent years, the study of the dynamic relationship between predator and prey with IGP has long been one of the most important themes in population dynamics because of its universal existence in nature. Extinction and coexistence of species are two important issues in systems with IGP. Most studies about predator and prey models with IGP were assumed in a homogeneous environment in space and different ODE models were set up. In this paper, we constructed a mathematical model with IGP in the presence of diffusion and B-D functional response between the predator and prey. Here we assumed that the IGpredator may be a generalist with a exclusive resources or just a specialist. We obtained some sufficient conditions for the extinction and permanence of the time-dependent system. We also obtained some necessary and sufficient conditions for the existence of coexistence states and further discussed the uniqueness and stability of coexistence state under some special conditions. The existence of a non-constant time-independent positive solution, also called stationary pattern, always indicates the dynamical richness of the system. We found that the dynamic relationship between predator and prey with IGP may be more complex than that of the models without IGP such as food chain or two predators one prey models.
Generally, it is not enough to just know some conditions on extinction and existence of coexistence states and we hope to know how the parameters in system affect the distribution of species. Obviously, it is difficulty to obtain these results by just using the analytic methods. Thus by using numerical simulations, we studied the effects of the parameters in model on the spatial distribution of species and obtain some interesting results about the extinction and coexistence of species. For example, we found that different interference levels among predators can both lead to coexistence of species under some circumstances and to extinction of species under some other circumstances. We also found that the introduction of IGpredation in a food chain model and the growth rate (or death rate) of IGpredator have similar effects on the extinction and coexistence of species. These results verify and complement the analytic results in section 4.