AVOIDANCE BEHAVIOR IN INTRAGUILD PREDATION COMMUNITIES: A CROSS-DIFFUSION MODEL

A cross-diffusion model of an intraguild predation community where the intraguild prey employs a fitness based avoidance strategy is examined. The avoidance strategy employed is to increase motility in response to negative local fitness. Global existence of trajectories and the existence of a compact global attractor is proved. It is shown that if the intraguild prey has positive fitness at any point in the habitat when trying to invade, then it will be uniformly persistent in the system if its avoidance tendency is sufficiently strong. This type of movement strategy can lead to coexistence states where the intraguild prey is marginalized to areas with low resource productivity while the intraguild predator maintains high densities in regions with abundant resources, a pattern observed in many real world intraguild predation systems. Additionally, the effects of fitness based avoidance on eigenvalues in more general systems are discussed.


Introduction
Intraguild predation (IGP) describes a situation where a predator and prey species also compete for a shared resource. It has been observed in a variety of ecological communities including avian, both large and small mammal, reptile, insect, fish and bacteria. In fact, any ecosystem with a complex food web is likely to have examples of intraguild predation within it. In [4] a database of 113 food webs was analyzed for presence of intraguild predation and it was found to be present at high frequencies throughout.
One of the earliest attempts to rigorously model intraguild predation was by Holt and Polis in [11] where basic ODE models were developed for a three species community. One of the conclusions reached in [11] was that their model with strong IGP was particularly prone to species exclusion, even though communities with strong IGP appear to be widespread in nature. Holt and Polis suggested numerous lines of future research on mechanisms to stabilize coexistence states in IGP communities. One of these was to allow for a heterogeneous environment.
Numerous empirical studies have found cases where the intraguild prey (IGPrey) employs nonrandom dispersal in foraging behaviors and habitat selection in an apparent effort to reduce its risk of predation [9], [15], [18], [21], [22], [25]. A common thread in these studies is that the intraguild predator tends to concentrate in the highest quality parts of the habitat (measured by basal resource productivity) while the intraguild prey will concentrate in areas of marginal habitat quality in order to avoid the increased predation pressure present in the higher quality areas. We propose a model where the intraguild prey (IGPrey) employs a fitness based dispersal strategy that accounts for local resource availability and predation risk. We assume that the resource and intraguild predator (IGPredator) will disperse randomly. We analyze what effects this non-random dispersal strategy has on the long-term population dynamics and steady-state distributions of the community.
Amarasekare was the first to model IGP in a heterogeneous environment; first with random movement [5], and then with non-random movement strategies (density, habitat and fitness dependent dispersal were all considered) [6]. Both of these models use an environment consisting of 3 distinct patches, each with a different level of resource productivity. The dynamic equations take the the form of an ODE for each species in each patch. Due to the size of the system (9 equations) all conclusions were based on numerical simulations. This work differs in that we take space to be a continuous variable resulting in a 3-dimensional quasilinear parabolic system of PDEs. We are able to prove the existence of a global attractor for this system and derive conditions for the uniform persistence of the IGPrey.
There has been some past work modeling IGP communities that have incorporated negative penalties from high IGPredator density into the functional response terms of the IGPrey in an effort to model prey vigilance, adaptive foraging, and other anti-predation behavior [12], [17]. These works differ significantly from the model we propose below as they do not explicitly model space and hence do not allow the IGPrey to modify its movement behavior in response to local conditions. We will assume that the habitat, Ω, is a domain in R 2 with smooth boundary. We will use u(x, t) to denote the density of the resource species, v(x, t) for the IGPrey and w(x, t) for the IGPredator. The system we consider is The per-capita growth rates (also referred to as fitness functions), f, g and h are given by The constants d 1 and d 3 are the motility of the resource and IGPredator respectively. The function M (u, w) is the motility of the IGPrey. It is assumed that the IGPrey changes its movement strategy based on the local density of the resource as well as the IGPredator. We assume that M is twice differentiable in u and w and there is a positive constant d such that We will prove global existence of solutions to (1) without specifying a specific form for M , but in Section 3 we will specify M after introducing the concept of fitness based avoidance strategies. The function r(x) is the spatially varying resource productivity which affects both resource growth rate and carrying capacity. We assume that r(x) is C α (Ω) for some α ∈ (0, 1) and that r(x) > 0 on Ω. The parameters µ 1 and µ 2 are the natural mortality rates of the IGPrey and IGPredator and ω 1 , ω 2 and ω 3 are self-limiting/crowding coefficients for each species. The a i parameters are the attack rates, the h i 's are the handling times and the e i 's are the conversion efficiencies of each predation/consumption functional response. We will assume that a i , h i , e i , µ i , ω i and d i are all positive for i = 1, 2, 3. The combined predation term in (5), e 2 a 2 u + e 3 a 3 v 1 + h 2 a 2 u + h 3 a 3 v , arises from assuming that the IGPredator will indiscriminately attack the resource or IGPrey when encountered while searching and does not a priori partition its search time between searching for resources and searching for IGPrey. Note that we have imposed Neumann conditions in (2). This is equivalent to zero flux conditions (replacing ∂v We see from (7) that if ∂v ∂n = 0 in addition to the conditions already imposed on u and w then we will have zero flux for the v-component equation. Conversely, in order to achieve zero flux in the v equation and maintain zero flux in the u and w equations we would need to impose ∂v ∂n = 0.
We will now prove that for nonnegative initial conditions in the Sobolev space W 1,p (Ω) with p > 2, the system (1)-(2) has unique solutions that are global (exist for t ∈ [0, ∞)).

Global Existence and the Compact Global Attractor
For p > n we have W 1,p (Ω) → C(Ω), so in this case we may define [W 1,p + (Ω)] 3 to be the cone of nonnegative triples in [W 1,p (Ω)] 3 .
In [1], [2] and [3], Amann proved existence results for a class of quasilinear parabolic equations that includes (1)- (2). Theorem 1 of [2] implies that for initial conditions in [W 1,p + (Ω)] 3 with p > n there is a unique classical solution to (1)-(2) with a corresponding maximal interval of existence, J (which may be [0, ∞)). Theorem 3 of [2] implies that if the L ∞ -norms of all solution components are bounded for t ∈ J, then the solution exists globally in time, i.e. J = [0, ∞). Standard comparison principles for single parabolic equations with coefficients that depend on time and space can be applied to these classical solutions to conclude that they remain nonnegative on Ω for all t ∈ J (see [14]).
Le proved stronger results for a two component system with additional assumptions bounding the growth of the flux and reaction terms in [13]. Le's system included one component with cross-diffusion and one without. We will use Theorem 2.2 of [13] to prove that for any initial conditions in [W 1,p + (Ω)] 3 , p > 2, the system (1)-(2) has solutions that exist globally in time. Furthermore, when viewed as a semiflow on [W 1,p + (Ω)] 3 the system will have a compact global attractor. In order to state this theorem, we need to define the class of ultimately uniformly bounded functions in [W 1,p + (Ω)] 3 . Definition 1. Let X be a complete metric space and ω : [0, b) × X → R where b ∈ (0, ∞]. We say that ω is ultimately uniformly bounded with respect to X if there exists a function C 0 : R + → R + such that Let P be the set of ultimately uniformly bounded functions with respect to [W 1,p + (Ω)] 3 . For the two component system in [13], Le showed that if the L n -norm of the component with cross-diffusion is in P and the L ∞ -norm of the component without cross-diffusion is in P, then solutions exist globally and there exists an α > 0 such that the C 1+α (Ω) norm of each component is also in P. Although (1)-(2) has three components, Le's main result, Theorem 2.2 of [13], generalizes to this system (the key element being that at least one of the components is a standard reaction-diffusion equation with no density dependence in the flux).
Applying Theorem 2.2 of [13] requires that the reaction terms satisfy conditions analogous to condition (H.2) in [13]. For (1)-(2) this means we must have a continuous function C(u, w) such that Returning to (3)-(5) we see that these conditions are satisfied since v does not appear in any term of (3) or (5) with an exponent greater than one and the term involving v in (4) is negative. Note that these growth conditions would also be satisfied for Lotka-Volterra functional response terms (i.e. h i = 0 for i = 1, 2, 3) in (3)-(5). Solutions to (1)-(2) with initial conditions in [W 1,p + (Ω)] 3 are twice differentiable in space for t ∈J, hence for any Banach space X with C 2 (Ω) → X, we can view the solution norms in X, u(·, t) X , v(·, t) X and w(·, t) X , as functions from J ×[W 1,p + (Ω)] 3 to R, and talk about these norms belonging to P. To keep the notation compact, we will write u X (or u(t) X when the time variable needs explicit mention) instead of u(·, t) X to indicate these functions. In order to apply Theorem 2.2 of [13] to system (1)-(2), we need to prove that u ∞ , v 2 and w ∞ ∈ P.
The u and w component equations involve only standard diffusion operators so comparison principles for parabolic equations can be effectively utilized. Because the predation terms of the u equation in (1) are negative, the u component of a solution to (1)-(2) is a subsolution to the logistic diffusion equation Let r = max x∈Ω r(x). A well known result (see [8]) for (11) is hence u ∞ ∈ P. The same reasoning applies to w ∞ upon noting that It is worth noting that establishing this w ∞ bound is the only place where it is necessary to have a saturating functional response. The analyses of this paper remain valid with h 1 and h 2 = 0; however, it is required that h 3 > 0 to bound w before establishing an upper bound on v.
The main difficulty comes in proving v 2 ∈ P. We can adopt the argument used to prove Theorem 3.1 of [13] to establish this bound. The system Le considered as the subject of Theorem 3.1 was a two component Lotka-Volterra competition system with logistic self-limiting on both terms. The important properties of our system that let the argument carry over are that the IGPrey has a logistic self-limiting term and that we are able to establish the L ∞ bounds on u and w prior to bounding v. We will now give a brief outline of the proof of v 2 ∈ P for our system.
First we establish that v 1 and t+1 t v(s) 2 2 ds are both in P using the same method as Lemma 3.3 of [13] (where the IGPrey self-limiting plays a critical role). We will use the Uniform Gronwall Lemma as found in [24] to pass from t+1 t v(s) 2 2 ds ∈ P to v 2 ∈ P; but first, we must establish an inequality of the form Manipulations analogous to equations (3.6)−(3.10) of [13] will yield such an inequality with η 1 (t) = C( ∇u 4 4 + ∇w 4 4 ) and η 2 (t) = C(1 + ∇u 2 4 + ∇w 2 4 ). To arrive at this bound Gagliardo-Nirenberg inequalities (as found in [16]) are used, which depend on the dimension of the domain Ω and are valid in this case for n ≤ 2.
It is then left to show that t+1 t ∇u(s) 4 4 ds and t+1 t ∇w(s) 4 4 ds ∈ P. This is accomplished by applying a second Gagliardo-Nirenberg inequality (again using the fact that n ≤ 2) to bound ∇u 4 4 in terms of ∆u 2 2 and other L p norms of u already shown to be in P. To this end, we need t+1 t ∆u 2 2 ds ∈ P. To show t+1 t ∆u 2 2 ds ∈ P, we can trace the steps of Lemma 3.4 in [13] and first prove that ∇u 2 ∈ P and then use this to establish t+1 t ∂u ∂t (s) 2 2 ds ∈ P. Then, rearranging the u component equation of (1) and integrating over Ω and from t to t + 1 yields From the bound (8) on the size of |u f (x, u, v, w)| and the fact that u ∞ and w ∞ ∈ P there exists an η(t) ∈ P such that We have already established that v 1 and t+1 t v(s) 2 2 ds ∈ P, therefore t+1 t u(s)f (s) 2 2 ds ∈ P. We then conclude from (14) that t+1 t ∆u 2 2 ds ∈ P and hence t+1 t ∇u 4 4 ds ∈ P. The same argument is then used to show t+1 t ∇w 4 4 ds ∈ P. Finally, the Uniform Gronwall Lemma can be applied to inequality (13) to conclude that v 2 ∈ P. 3 . This solution exists for all t ≥ 0 and, furthermore, there exists an α > 0 such that The system (1)-(2) defines a semiflow on [W 1,p + (Ω)] 3 and this semiflow possesses a compact global attractor.

Motility Functions Modeling Avoidance.
In this section we set forth conditions that define a class of motility functions, M (u, w), that model the intraguild prey employing a fitness based avoidance strategy, i.e. avoiding areas with a bad resource/predation risk trade-off. We then establish a sufficient condition based on this class of motility functions for the IGPrey to be uniformly persistent in the system (asymptotically bounded above and below by a positive constant that does not depend on initial conditions). We will assume that the IGPrey is able to assess the local density of resources and frequency of predator attacks. Both of these assumptions are reasonable for a variety of species [9], [18], [22], [25]. The IGPrey will use the resource availability and frequency of predator attacks as a way of judging the local environmental quality and increase its motility in regions that it judges to be bad while maintaining a lower base level of motility (or perhaps even decreasing it) in regions that it judges to be good.
The per capita growth rate (fitness) function, g(u, v, w), is a good candidate to represent the IGPrey's assessment of local environmental quality. Where g is positive, the abundance of resources outweighs any present predation risk (good areas) and where g is negative, the resources are not sufficient to warrant the risk (bad areas). However, because of the logistic self limiting term, g depends on the local intraspecific density v, and so we cannot choose M = M (g) within the framework we have developed thus far. Instead, we will use the function g * (u, w) = g(u, 0, w) and assume M (u, w) = M (g * (u, w)). The main thrust of our analysis deals with the case when the IGPrey is rare and invading the system, so this should be a good approximation to the actual fitness in these cases. We will refer to g * as the IGPrey's linearized fitness.
We will use an additional parameter, λ, to represent the degree to which the IGPrey is employing this avoidance strategy. More precisely, we will consider a family of motility functions, {M λ (g * )} λ≥0 such that M λ (g * ) ≥ d 2 for all λ and g * < 0 and, Condition (19) states that when the linearized fitness is negative, λ increasing toward infinity will result in ever higher motility for the IGPrey. In this sense, higher values of λ make the IGPrey increase its avoidance response to areas that it judges to be bad. We also require that M λ is twice differentiable with respect to g * for all λ ≥ 0, which implies from (17) and (18) that M λ (0) = d 2 for all λ. (1)-(2), and A the global attractor in Y that was shown to exist in Theorem 1. In order to study the asymptotic dynamics of π it suffices to examine orbits originating in a small neighborhood of A denoted by B(A, ε). WriteY for the interior of Y , i.e. triples of strictly positive functions on Ω, and ∂Y for its boundary, nonnegative triples where at least one function is zero at some point in Ω. Take t 1 , t 2 > 0 and setX = π(B(A, ε), t 1 ), X = π(X, t 2 ), and S = X ∩ ∂Y . Then X, S and X \ S are compact and forward invariant under π (see Theorem 4.1 of [8]). By virtue of the maximum principle for single parabolic equations, S consists of elements of X where at least one component is identically zero on all of Ω, and X \ S consists of strictly positive functions.

Acylclicity Argument for Uniform
We will derive a condition for uniform persistence of the IGPrey by considering its ability to invade the system when absent. Define S uw to be the subset of S where v ≡ 0 and A uw to be A∩S uw . To investigate the structure of A uw we need to examine the u − w subsystem This is a standard system of the type considered in [8], Section 4.5. This system has two boundary equilibria: (u * , 0) and (0, 0) where u * (x) is the unique positive solution of the diffusive logistic (steady-state) equation Any initial condition of the form (u 0 , 0), u 0 not identically zero, converges to (u * , 0); whereas, any initial condition of the form (0, w 0 ) converges to (0, 0). Furthermore, if the principal eigenvalue, σ 1 , of (22) is positive, then this subsystem is permanent and there exists a compact invariant set, A 1 , that is bounded away from the boundaries and attracts all initial data of the form (u 0 , w 0 ), neither u 0 nor w 0 identically zero. However, if the principal eigenvalue σ 1 is negative, then (u * , 0) attracts all such initial data of the subsystem. Therefore, we have A uw = {(0, 0, 0), (u * , 0, 0)} ∪ A 1 (using the convention A 1 = ∅ in the case of σ 1 < 0). We will now establish some lemmas that will be useful in proving our main result.
For a point x ∈ X and a set U ⊆ X define the distance from x to U in the standard way and the unstable set of U by The compactness of U implies that if x ∈ W S (U ) then lim t→∞ d(π(x, t), U ) = 0.
Lemma 1. Let U ⊆ S uw be a compact set. Suppose there exists a continuous function b(x) such that b(x) ≤ g(u,0,w) M (u,w) for all (u, 0, w) ∈ U and the principal eigenvalue σ 2 of for all (u, v, w) ∈ B(U, ξ). Take t 0 such that d ((u(t), v(t), w(t)), U ) ≤ ξ for all t > t 0 . Multiply the v component equation of (1) by v 2 and multiply (25) by M (u, w)v and integrate over Ω to obtain and (28) Apply the divergence theorem twice to the first term in (28) and subtract (28) from (27) to get for all t > t 0 , and we have also assumed that M is such that M (u, w) ≥ d > 0 for all u, w ≥ 0.
Utilizing the fitness based avoidance strategy, M = M λ (g * (u, w)), we will derive conditions that guarantee a b(x) satisfying the requirements of Lemma 1 exists. Given any U that is a compact subset of [C(Ω)] 3 , define and b λ,U (x) = min These functions are both continuous on Ω by the following Lemma: Lemma 2. Suppose U is a compact subset of [C(Ω)] n and φ : R n → R is a continuous function. Then, the function F : Ω → R defined by is continuous.
Proof. F is well-defined (the minimum is obtained by using a function in U at each point in Ω) because U is compact and uses the max-norm. Fix x 0 ∈ Ω and suppose F is not continuous at x 0 . Then there is an ε > 0 and a sequence {x n } ∞ n=1 ⊆ Ω such that x n → x 0 and |F (x n ) − F (x 0 )| > ε for all n. Let {f n } ⊆ U be such that F (x n ) = φ(f n (x n )), for n = 0, 1, . . .
There is a subsequence, {f n k } that converges in [C(Ω)] n to a functionf ∈ U . Choose K sufficiently large so that and which contradicts the minimality of F (x n k ) for f ∈ U . Now, suppose F (x 0 ) > F (x n k ) + ε. Then (33) and (34) imply that which contradicts the minimality of F (x 0 ) for f ∈ U . Therefore, we must have that F (x) is continuous on Ω.
Proof. If U 1 and U 2 are two compact invariant subsets of S, we will say that U 1 is chained to U 2 in S and write U 1 → U 2 if there exists a u ∈ S \(U 1 ∪U 2 ) such that u ∈ W u (U 1 ) ∩ W s (U 2 ), i.e. there is a full orbit, γ(u), passing through u with that connects to U 1 as t → −∞ and connects to U 2 as t → ∞. We say a collection of compact invariant sets, {U 1 , . . . , U m } forms a chain if U 1 → U 2 → . . . → U m , and we say the collection forms a cycle if U m = U 1 .
Dissecting the proof the acyclicity theorem of [10], we can derive the following sufficient condition for v to be uniformly persistent in (1)-(2) for all initial data with neither u 0 nor v 0 identically zero. This will be the case if all of the following are true: (1) No subset of C = {(0, 0, 0), (u * , 0, 0), A 1 } can form a cycle in S.
Addressing item (1): no subset of the three isolated components of A uw can be chained together in S uw to form a cycle. We have (0, 0, 0) → (u * , 0, 0), and if A 1 is not empty we also have (0, 0, 0) → A 1 and (u * , 0, 0) → A 1 ; but, that is all of the existing chains, none of which can be linked to form a cycle.
To address item (2), consider the case when w 0 ≡ 0. Let S uv be the subset of S where w ≡ 0. Applying the acyclicity theorem of [10] to the u − v subsystem, we see that it will be permanent as long as W s ((u * , 0)) does not intersect the interior of S uv (the subset of S uv where u and v are strictly positive). Applying analogues of Lemma 1 and Lemma 3 to guarantees that this is the case for λ sufficiently large. Therefore, there is a compact subset, A 2 , in the interior of S uv that attracts all initial data of the form (u 0 , v 0 , 0), neither u 0 nor v 0 identically zero.
Now, consider what happens for initial data with u 0 ≡ 0. The resulting subsystem is Even though the reaction terms in the first equation of (40) are always strictly negative, the nonlinear diffusion prevents us from using a comparison principle argument to conclude that solutions converge to (0, 0, 0). Instead, multiply the first equation in (40) by e 3 , add it to the second equation and integrate over Ω (applying the divergence theorem to eliminate the Laplacian terms) to obtain We can drop the quadratic terms from the right hand side and set k = min{µ 1 , µ 2 } to get the inequality and then conclude that e 3 v + w 1 → 0 as t → ∞, hence v 1 , w 1 → 0 individually. The semiflow restricted to S, π S (·, t) is compact for t > 0, so for any sequence {t n } with t n → ∞, the sequence of solutions to (40) evaluated at these points in time, {(0, v n , w n )}, has a convergent subsequence in S. However, from (41) we know that this convergent subsequence must be converging to (0, 0, 0) in L 1 (Ω) 3 which implies that it is converging to (0, 0, 0) in S as well. Since every sequence of solution points has a subsequence converging to (0, 0, 0) in S, the solution trajectory must be converging to (0, 0, 0) in S as well. Therefore, the only isolated component of A∩(S \ S uw ) is A 2 (the interior global attractor for the u−v subsystem) which is not chained to any element of C in S. Hence condition (2) is satisfied.

Using Lemma 1 and Lemma 3 with
and λ sufficiently large yields W s (M ) ∩ X \ S = ∅ for M ∈ {(u * , 0, 0), A 1 }. A standard argument that only considers the u-component can be used to show that W s ((0, 0, 0)) ∩ X \ S = ∅ (see [8] for an example). Thus condition (3) is satisfied and we can conclude that the IGPrey is uniformly persistent for initial data with neither the resource nor the IGPrey densities identically zero.
Theorem 2 can be viewed from an ecological perspective as saying that as long as the intraguild prey has one location in its habitat where conditions are favorable for all asymptotically feasible resource-IGPredator configurations, then applying a strong fitness based avoidance strategy will allow it to be uniformly persistent in the ecosystem.
This result is somewhat analogous to the well known result for a randomly diffusing species that as long as there is a point in the domain where linearized fitness is positive, then a species can invade provided its random diffusion rate is sufficiently small. However, the mechanisms are in fact quite different. In the case of small random diffusion, the ability to invade comes from a small flux out of favorable regions; whereas, in the fitness based avoidance case the ability to invade is a result of a high flux out of unfavorable regions. In fact, it would be perfectly reasonable to imagine realistic cases where both of these mechanisms are operating simultaneously.
It is worth noting that the conclusion of Theorem 2 does not directly depend on the IGPrey's random diffusion rate d 2 . The effect of increasing d 2 is that a larger value of λ may be needed (which in turn will result in even higher motility in areas of negative fitness) in order for the IGPrey to be uniformly persistent.
The theoretical insight provided by Theorem 2 is that fast random movement out of "bad" areas can have the same positive effect on long-term survival as slow random movement out of "good" areas. However, in practice, the structure of A 1 and A 2 are generally unknown, so directly verifying the conditions of Theorem 2 for a particular choice of parameters may be impossible. It was shown in [20] that there is a threshold value e * 2 (depending on the other parameters in the system) such that for e 2 > e * 2 the IGPredator can invade the (u * , 0, 0) equilibrium and the u − w subsystem is permanent. Furthermore, for e 2 > e * 2 but sufficiently close to e * 2 , the u−w subsystem has a unique positive equilibrium that is a global attractor for positive initial data, i.e. A 1 is a single point in function space, (û, 0,ŵ). In this case, the conditions of Theorem 2 would simplify to: there exists an x ∈ Ω such that g * (û(x),ŵ(x)) > 0 which is simple to check after computingû andŵ.
A sufficient condition for the uniform persistence of the IGPredator was developed in [20] following an argument along the same lines as the proof of Theorem 2. It requires that the IGPredator is able to invade all points in the A 2 attractor of the u − v subsystem which can guaranteed when a principal eigenvalue of an equation analagous to (25) of Lemma 1 is positive (see [20] for full details).

Numerical Examples
In [11] a three species ODE model for intraguild predation is considered and two necessary conditions for a positive coexistence equilibrium are derived: the IGPredator has sufficient gains from consumption of the IGPrey (e 3 sufficiently large), and the IGPrey is a superior competitor for the basal resource. In this section we will demonstrate by numerical simulation that neither of these conditions is necessary for coexistence in the spatially explicit model (1)- (2). Furthermore, we will see how fitness based avoidance strategies employed by the IGPrey effect the existence of these coexistence states and the population densities at equilibrium.
Our results were proved for a two dimensional domain but are also valid for one dimension. We will perform our numerical experiments using the unit interval for the domain Ω to simplify the visualization of results. Solutions were obtained using the built-in 1-D parabolic PDE solver pdepe of Matlab which implements the discretization scheme described in [23]. For a description of a numerical scheme suitable for a two dimensional domain refer to [20]. For both numerical experiments we will use a resource productivity function, r(x), of the form which is a logistic sigmoid function centered at x = 1/2, saturating at the values r min on the left of the domain (low quality habitat) and r max on the right of the domain (high quality habitat). The parameter α controls the sharpness of the transition of the transition layer around x = 1/2 and will be taken large (α = 40) so that this transition is fairly sharp. Throughout this section we use r min = 0.6 and r max = 1.25. For both examples we take ω 1 = 1 which implies that r(x) is also the local carrying capacity of the resource. For both scenarios we will assume that the resource diffuses slowly, d 1 = 0.0005, so that the effect of the spatially varying habitat quality is not diffused throughout the domain by a fast moving resource. In order to clearly compare the ability of the IGPredator and IGPrey to compete for the shared resource we will assume that attack rates and handling times for the resource and natural mortality parameters are identical between the IGPredator and IGPrey (i.e. a 1 = a 2 = 0.8, h 1 = h 2 = 0.1, µ 1 = µ 2 = 0.2, ω 2 = ω 3 = 0.1) and only the conversion efficiencies are unequal (e 1 = e 2 ).
A simple family of motility functions, {M λ (g * )} λ≥0 , that satisfies conditions (17)- (19) is given by the piecewise linear functions However, this family of motility functions is not twice differentiable at g * = 0 as required. Instead we will use the exponentially smoothed approximation of (43) given by that still satisfy (17)- (19) and are infinitely differentiable at g * = 0. The first scenario we consider is one where the IGPredator does not gain significantly from consumption of the IGPrey (e 3 small or zero). We assume that the IGPredator diffuses slowly (d 3 = 0.001) and is not efficient enough to maintain a high density in the low quality habitat (where r(x) ≈ r min ). The IGPrey is a more efficient competitor for the shared resource (e 1 > e 2 ) and can subsist on the resources in the low quality habitat. Even though the IGPredator does not directly gain significantly from consumption of the IGPrey, we will assume that it is highly aggressive towards the IGPrey (a 3 large). A good representative for this type of scenario might be mammalian carnivore systems with multiple predators of varying size competing for a shared prey [19]. The remaining parameters for this example are: d 2 = 0.05, e 1 = 0.9, e 2 = 0.4, e 3 = 0.1, a 3 = 2, h 3 = 0.5. These parameters result in a u − w subsystem with a positive steadystate that appears to be a unique stable node for the subsystem (all initial conditions tried converged to this steady state). The linearized fitness of the IGPrey at this steady-state is slightly positive in the part of the domain with marginal resources and low IGPredator density but significantly negative in the high resource portion of the domain (see Figure 1). For d 2 = 0.05 the IGPrey is unable to invade this equilibrium using only random diffusion (λ = 0), but can invade and persist using fitness based avoidance (λ = 5). The resulting positive three species steady state after IGPrey invasion is pictured in Figure 1. Note that the IGPrey concentrates in the marginal resource portion of the domain while the IGPredator concentrates in the higher quality portion of the domain; a pattern observed in many actual IGP systems in nature [18], [25]. This type of result is also possible with e 3 = 0, the interspecific killing case (the IGPredator is not consuming the IGPrey, just attacking it).
The second scenario we consider is one where the IGPrey is actually an inferior competitor for the shared resource (e 1 < e 2 ) but is able to invade and persist using fitness based avoidance by exploiting areas where the IG-Predator has under exploited the available resources due to over dispersion.  For this to be possible, we will assume that the IGPredator has a moderate random diffusion rate (d 3 = 0.1) and is only mildly aggressive towards the IGPrey (a 3 small). The parameters used for this scenario are: d 2 = 0.05, e 1 = 0.55, e 2 = 0.7, e 3 = 0.1, a 3 = 0.1, h 3 = 0.5. The results are shown in Figure 2. Due to a higher random diffusion rate, the IGPredator is much more evenly distributed across the entire domain in the u − w equilibrium. The result is a very depressed resource level in the marginal resource area of the habitat. The resulting linearized fitness of the IGPrey is negative in this area but slightly positive in the high quality portion of the habitat. This is the case even though it is an inferior competitor for the resource and it has the extra predation pressure of the IGPredator. As in Scenario 1, the IGPrey cannot invade with d 2 = 0.05 and random diffusion alone (λ = 0); however, with sufficient fitness based avoidance (λ = 5) the IGPrey invades and persists by concentrating in the high quality portion of the habitat. This equilibrium profile might seem odd if observed empirically: a species that is an inferior competitor concentrating in a part of the habitat where its competitor and predator is highest density and coexisting there. In order to understand this type of equilibrium profile it is essential to not only understand the species interactions but also the underlying movement mechanisms of all of the species in the system.

Discussion
Using the example of an intraguild predation system we have demonstrated how a fitness based avoidance strategy, namely increasing motility in areas of negative fitness, can lead to persistence of the IGPrey in the system.
However, regardless of the community structure, this same general principle can always facilitate an invader to increase when rare as long as there is one favorable spot in the habitat by damping negative contributions to the relevant principal eigenvalue from areas in the domain where conditions are unfavorable. As pointed out earlier, this has a similar effect as decreasing a random diffusion rate; but, has very different interpretations. A small diffusion rate allows increase when rare by limiting losses resulting from flux from good regions to bad regions. The mechanism proposed herein allows a species to increase when rare by increasing the flux out of bad regions into good regions.
This might lead one to compare fitness based avoidance to a behavior like advection up fitness gradients, but there are important differences from the perspective of a mechanistic derivation of the movement behavior. The use of fitness based avoidance does not assume that organisms have the ability to track gradients at all, only the ability to assess local conditions for suitability and adjust the rate at which they depart the location in a random direction. The net effect at the aggregate population level can be decomposed into a taxis (gradient tracking) and kinesis (variable rate random motion) component: ∆ [M (g)v] = ∇ · M (g)∇v + vM (g)∇g but, at the organism level all movement is in random directions.
When strong fitness based avoidance is combined with moderate random diffusion, the resulting asymptotic population distribution will tend to have very low density in regions where fitness is negative and concentrate with fairly flat, homogeneous profiles in areas with positive fitness (due to the moderate random diffusion in these areas). For contrast, a species that uses only a very small random movement will also have low density in areas of negative fitness, but the density where fitness is positive will much more closely track the shape of the fitness curve itself. Ecologists doing empirical species distribution studies might be able to exploit the marked differences between these two distributions to more accurately infer underlying dispersal mechanisms.
Given a heterogeneous environment, the mechanism of dispersal for species in an ecosystem can play an important role in the types of asymptotic distributions of population densities observed and whether or not coexistence can occur. We believe that the framework of cross-diffusion systems is a very powerful tool for analyzing the interplay between species interactions with dispersal mechanisms. The system analyzed in this work is fairly simple from the perspective that only one species is dispersing in a non-random fashion. There is ample room to develop this theory further and extend this type of analysis to more complex cases in the future.