PATTERN FORMATION IN A PREDATOR-MEDIATED COEXISTENCE MODEL WITH PREY-TAXIS

. Can foraging by predators or a repulsive prey defense mechanism upset predator-mediated coexistence? This paper investigates one scenario in- volving a prey-taxis by a prey species. We study a system of three populations involving two competing species with a common predator. All three popu- lations are mobile via random dispersal within a bounded spatial domain Ω, but the predator’s movement is inﬂuenced by one prey’s gradient representing a repulsive eﬀect on the predator. We prove existence of positive solutions, and investigate pattern formation through bifurcation analysis and numerical simulation.


1.
Introduction. It has been known for decades that predation can enhance species diversity, especially within aquatic environments. That is, the presence of a common predator can increase species diversity in competitive situations in what is referred to as predator-mediated coexistence (see for example, [13,20,30,32,33]). Questions of the relationship between species competition and defense remain of current interest (see, for example [36]). The impact of the predator on coexistence need not be the same for all prey species and depends on mechanisms involved in predation (see [10] and references there in). One question we address is how prey-taxis, either attractive or repulsive, may impact the emergence of predator mediated coexistence. Literature for predator-prey systems and two species in competition suggest that attractive prey-taxis does not destabilize a coexistence steady state [37,16]. As such we focus our motivation on exploration of how the behavior response of avoidance of self harm (repulsive prey-taxis) may impact predator mediated coexistence. We will also ask if foraging behavior by the predator (attractive prey-taxis) that does not destabilize coexistence in a predator-prey system may destabilize a predator mediated coexistence state.
In general there is a tendency for prey to resist damage from a predator with many prey species having adapted anti-predator mechanisms to decrease vulnerability to a predator. Some anti-predator behaviors such as the release of deterrent chemicals or noises, or mimicry may result in avoidance behavior by the predator. However, it remains unclear what is the impact on predator-mediated coexistence if one of the prey species were to have such a defense mechanism that leads to avoidance behavior by the predator. If the predator were to be repelled, would the dynamics of competition reign and lead to potential competitive exclusion; even if the weaker competitor is the prey equipped with a defense mechanism? In this article we analyze a model for species in competition with a common predator when one of the prey species is equipped with a defense mechanism to repel the predator. We propose a model of a chemotaxis type prey defense mechanism where the predator is detered by changes in the abundance of a chemical released by one prey. We study a model of a prey-taxis type prey defense mechanism where the predator is detered by higher prey densities due to the prey defenses. This latter model is considered as an approximation of the former in the case where there is rapid spread of the chemical.
Assume we have two competing species that are prey to a single predator, with all three species being mobile within a bounded spatial domain via simple diffusion. Assume one of the prey species produces a chemical that acts as a repulsive prey defense mechanism to the predator species. Then one question would be how does prey producing the chemorepellent-like mechanism fare relative to the other prey in the presence of the predator species. If, in a connected, bounded spatial region Ω, u and v represent two competing population densities, and w represents a common predator to u and v, then a reasonable model would be 3 vw w t = ∇ · (d 3 ∇w + χ 0 w∇z) + ew − κw 2 + ε 1 α 3 uw + ε 2 β 3 vw z t = d 4 z − µ z z + δu Here K u , K v are carrying capacities for these species that are critical densities beyond which there is insufficient environmental resources to support further population growth. The α 1 , β 1 are intrinsic growth rate parameters, and α 2 , α 3 , β 2 , β 3 are conversion rate parameters, that is, the rate per interaction of mortality on the prey or competitor, or growth benefit to the predator due to the prey-predator interaction. The efficiency of the predator in receiving a growth benefit from predation is given by ε 1 and ε 2 . Also, e is an intrinsic growth (death) rate for the predator, κ is a self-limiting or crowding coefficient for the predator, and χ 0 is a positive chemotactic sensitivity parameter, whose term represents a repulsive mechanism toward the predator due to production of chemical signal z by species u (at rate δ per individual). The chemical can become degraded, or lost by other means, hence the µ z term. While all parameters in the model are positive constants, e and χ 0 are the only parameters in this model that can reasonably be treated as a negative constant. A negative value for χ 0 would correspond to an attractive mechanism. We will provide a global existence proof for e as either positive or negative and provide examples where e is positive.
If we instead assume the chemical acts quickly so a quasi-steady state assumption can be employed, and that the chemical does not diffuse, then an alternative model formulation would be This can be considered a direct repulsive reaction of the predator species with respect to the presence of the u species, while the first model would represent the reaction being indirect based on a chemosensory mechanism as a defense scent put out by the u prey species. In this paper we will concentrate on analysis of system (1) only. First, we scale the system (1) to place it in a non-dimensional form. Letũ = u/K u ,ṽ = v/K v ,w = w/w char . We do not need to specify w char , but typically it could be w char = e/κ. With a characteristic length scale L = d 3 /α 1 , let the characteristic time scale be T = L 2 /d 3 = 1/α 1 , so thatt = t/T andx = x/L,ỹ = y/L. Here the length scale is associated with the predator migration rate per unit time, which makes the characteristic time scale tied to the intrinsic growth rate of the u prey population.
Then substituting everything into system (1) and dropping the tilde notation, we have We assume this system is defined on the connected, bounded domain Ω ⊂ R n , which has smooth boundary, and that u, v, w satisfy no-flux boundary conditions and initial conditions where u 0 , v 0 , w 0 ≥ 0 on Ω. Here ∂/∂ν = ν · ∇, where ν is the unit outward normal vector on ∂Ω.
In the next section we mention some relevant previous literature, but not trying to be comprehensive since the population literature is so extensive. In section 3 we establish local and global solvability of system (2). In section 4 we discuss the stability of the constant, positive steady state solutions, and conditions for convergence of solutions to (2)-(4) to the steady state solution. We also treat the "prey-tactic" sensitivity parameter χ as a bifurcation parameter and establish conditions for pattern formation. We then fix χ and give a brief discussion of destabilizing the system by considering the prey diffusivity d as a bifurcation parameter. In section 5 we present some numerical results for our model system. In section 5.1 we outline our numerical method and in section 5.2 numerical results concentrating on two scenarios, one where we would expect prey u to "win", with prey v going to extinction, and one where we would expect v to "win." We summarize and briefly discuss unresolved issues in section 6.
As far as the authors know, this investigation, while mainly using standard analysis techniques, is the first report incorporating a chemotaxis-type term into the model ecosystem consisting of two competing species with a common predator, and exploring pattern formation via prey-tactic and diffusive mechanisms.
Although we have restricted ourselves to Lotka-Volterra type dynamics, which is a common, but narrow, modeling viewpoint of species interactions, results below show sufficiently rich pattern behavior to justify investigation of the model.
2. Some background. Understanding spatial and temporal behaviors of interacting species in an ecological system is an important problem in population ecology. Our model falls into the realm of a reaction-diffusion-advection system. One survey of such systems, with emphasis on ecological considerations, is [12]. A couple of reviews of taxis mechanism model systems include [18], and [17]. Older papers involving one predator and two prey species, with diffusion, but without chemotaxis, include [15], and [9]. [15] shows existence of positive steady state solutions, their stability, and some limiting behavior. [9] add a crowding effects term and gives conditions for persistence of the populations. For steady state solutions in this context, there is [25], and [23]. [25] uses degree theory in cones, positive operator techniques, and sub-and super-solution techniques.
[1] models a predator-prey dynamic with prey-taxis, that is direct movement of predators in response to a prey gradient, and prove existence of a weak solution by fixed-point methods, and uniqueness by a duality technique. They also study linear stability around the equilibrium. Preytaxis models were derived by [22] and extended by [26] to discuss various functional forms in the predator-prey system and explore necessary conditions for pattern formation, and include some numerical simulations. See also [42], and [7] and their references. [16] with their competing population model with chemotaxis establish a global existence result for classical solutions in one space dimension, and global existence for higher dimensions when one species is at a quasi-steady state. They also use bifurcation theory to establish non-constant steady state solutions, and explore their stability. Tello and Wrzosek [34,35] also examine a competing species model with taxis, namely a chemorepellent. Some interesting variants using similar analytical techniques include one species producing both a chemoattractant and a chemorepellent (see [21], and [28]), and two (non-competing) species producing the same chemoattractant (see [24]). [39] though it involves a single species producing a chemoattractant, they introduce a population crowding (volume filling) mechanism (see also [41]) and explore pattern formation via stability and bifurcation analysis, and numerical simulation. Although we were originally motivated by the work of Tello and Wrzosek [34,35], perhaps the closest work to that presented here is that of [39] and [27]. [27] have two (non-competing) prey, one common predator, with prey-taxis in the sense of the movement of the predator depending on the gradient of both prey species.
In the recent work of [37] a single predator-prey system with a repulsive prey-taxis mechanism similar to the one we employ is considered. The prey dynamics have a logistic growth rate like our model, but the interaction term is of Holling type II and the diffusivities of the two populations are equal. Their analysis concentrates on a global bifurcation result, and stability of the resulting nontrivial stationary states.
3. Global existence in one spatial dimension. We first establish local existence of classical solutions to model system (2)-(4) in arbitrary dimension n ≥ 1. Since later numerical simulations are done in one space dimension, we complete a global existence result for n = 1. Recall that Ω ⊂ R n is connected, bounded, with smooth boundary, and outer unit normal defined on the boundary ∂Ω given by ν. Let p > n ≥ 1; then W 1,p (Ω; R n ) is continuously embedded in the continuous function space C(Ω; R n ). Let , depending on the initial conditions (u 0 , v 0 , w 0 ), such that problem (2) -(4) has a unique maximal classical solution Proof. Let z = (u, v, w) ∈ R 3 . Then we can write (2), (3), (4) in the form where Since the eigenvalues of A are positive, (5) is normally parabolic [5,4]. Thus, local existence in (i) follows from Theorem 7.3 and Corollary 9.3 in [5]. Because (5) is (lower-) triangular, (iii) follows from Theorem 5.2 in [4]. For (ii) (u, v, w) ≡ (0, 0, 0) is a strict sub-solution of the problem, then by the maximum principle and Hopf's boundary point lemma, the positivity of the solution follows.
To establish global in time solutions we need to show that u, v, and w are uniformly bounded in L ∞ (Ω). This is done in Theorem 3.5. A number of preliminary lemmas are needed for this theorem and are shown below.
In the rest of this section we assume b 1 ≥ c 1 , b 2 ≥ c 2 . This is equivalent to ε 1 (K u /w char ), ε 2 (K v /w char ) ≤ 1. So, if w char is greater than the carrying capacities of u and v, we are assuming the efficiency in converting consumed prey to new predators is sufficiently low.
From the comparison principle we have that Remark: By Jensen's inequality and Theorem 1.

Solving for
Hence, taken together with theorem 3.1 (ii) we have Corollary 1. Under the conditions of the above Lemma we have ||u(·, t)|| L 1 (Ω) and ||v(·, t)|| L 1 (Ω) are bounded for t ≥ 0, with the bound depending only on |Ω|, Because our numerical work below is all spatially one-dimensional, we assume in the rest of the paper that Ω = (0, L) is a bounded interval.
Here we modify and apply an argument given in [16] for species in competition.

Lemma 3.3.
Let Ω = (0, L) be a bounded interval in R. Assume that the nonnegative initial data (u 0 , v 0 , w 0 ) ∈ W 1,2 (Ω; R 3 ). Then there exists a positive constant C dependent on the parameters of (2) such that, for any p ∈ (1, ∞) Proof. For any dimension we have the integral form appealing to [40]. We have a similar expression for v(·, t). We also can write down comparable integral forms for u x (·, t) and v x (·, t), so any argument we make for u x (·, t) will have an analogous argument for v x (·, t); hence, we only discuss the u case. Taking the x derivative of the u equation, writing out the equation for u x (·, t) taking the L p norm, and appealing to Lemma 1.3, [40], we have for the one dimensional case, where λ is the first nonzero eigenvalue of − on our domain, with Neumann boundary conditions. Since Suppose the assumptions of Lemma 3.3 hold and let a positive classical solution of (2) be given by (u(x, t), v(x, t), w(x, t)). Then for each p ∈ (2, ∞), there exists a positive constant C(p) such that Proof. For p > 2 multiply the w equation of (2) by w p−1 and integrate over Ω = (0, L): (let (·) mean Ω (·)dx, and denote norm || · || L p (0,L) as || · || p , p ≤ ∞) where C 0 = c 1 ||u|| ∞ + c 2 ||v|| ∞ + |µ|. By Hölder's inequality and Lemma 3.4, and by Young's inequality If we let ε = 2/(pχC(2(p + 1))), then With an application of Young's inequality, Then the differential inequality becomes Since ||w|| 1 is uniformly bounded, let ||w|| Hence, another application of Hölder's inequality gives = w p , our differential inequality becomes Since α > 1, the solution y(t) is bounded for all t ∈ (0, T max ). Now the following corollary follows: where C is a positive constant dependent on the parameters of (2).

Theorem 3.5.
Let Ω = (0, L), then for any positive initial data Proof. To prove this theorem, we need to show that ||w(·, t) L ∞ is uniformly bounded for all t ∈ (0, T max ), and show that T max = ∞ ; Theorem 3.5 then follows from the Theorem 3.1(iii). Through the remainder of the proof we use the fact that solutions are nonnegative (i.e. u ≥ 0, v ≥ 0, w ≥ 0). For p > 2 we multiply the w equation of (2) by w p−1 and integrate over Ω; that is, proceeding as before, we have Without loss of generality we can set The following inequality follows from the Gagliardo-Ladyženskaja-Nirenberg inequality and Young's inequality that for any w ∈ W 1,2 (Ω) and ε > 0 where K only depends on |Ω| = L [11,16]. Choose ε = 2(p−1) and note p−1 and solve the differential inequality for each t ∈ (0, T max ).
Now, we employ Moser-Alikakos iteration [2] to establish an L ∞ (Ω) estimate for w. Define Noting that w 0 > 0 from (6) we have: Where C is a constant that only depends on |Ω| = L, and M (2 i0 + 2 −i+i0 ) is bounded for i > i 0 and for all t ∈ (0, ∞) in light of Lemma 3.2 and Lemma 3.4 that give L 1 and L p , p > 2 bounds respectively for t ∈ (0, T max ) for w.
T max ) and the main part of our operator is lower-triangular; thus, it follows from [4] that the uniform L ∞ bound ensures the extendability criterion for the maximal solution to get T max = +∞. 4. Constant steady state solutions. For the rest of the paper assume c 1 = b 1 , c 2 = b 2 , and η = 1. This is done for notational convenience only, reducing the number of parameters. Numerical simulations indicate no significant pattern formation change by employing this assumption. For system (2), noting this assumption, besides the trivial steady states (u, v, w) = (0, 0, 0), (0, 1, 0), (1, 0, 0), (0, 0, µ ), we have the "degenerate" steady state solutions, namely The requirement of positivity of u s , v s , w s places constraints on the system parameters (see next subsection).

Comments on the ODE systems. For the ODE probleṁ
we have the equilibrium states (7), (8). ODE theory gives global existence, positivity, and uniform boundedness of solutions for this system. For a good discussion on predator-prey and competitive species models, see [31]. Let us first consider the case of two competing species without the presence of a predator. In this case we have the semi-trivial steady states (0, 0), (1, 0), (0, 1) and Here are principle cases: • If a 1 < 1, a 2 < r, then (ū s ,v s ) exists, and is stable (it is a stable node); • If a 1 > 1, a 2 > r, then (ū s ,v s ) exists, but is unstable (it is a saddle point); • If a 1 < 1, a 2 > r, then (u, v) → (1, 0) as t → ∞; that is, "u wins"; • If a 1 > 1, a 2 < r, then (u, v) → (0, 1) as t → ∞; that is, "v wins" (u is the weaker competitor).
The last two cases are examples of Gause's principle or principle of competitive exclusion [31]. Including the predator in the dynamics, for linear stability about the equilibrium state (8) of the ODE system, consider The characteristic polynomial is λ 3 + B 2 λ 2 + B 1 λ + B 0 = 0. For (u s , v s , w s ) to be linearly stable, it is necessary and sufficient to have , so a necessary condition for stability is for B 1 > 0, and this is guaranteed when r − a 1 a 2 > 0. For the sufficient and necessary condition, we need the parameters to be such that Note that r − a 1 a 2 is only assured to be positive in the stable coexistence case; for all other cases, (10) requires that r − a 1 a 2 is not too negative.
Assume that for consistency with (10) that r − a 1 a 2 is not too negative. Further, we assume that Given (11), in order to have u s > 0, v s > 0, and w s > 0, respectively, assume Expression (10) is the condition for B 1 B 2 − B 0 > 0. Expression (11) and (12) gives existence of the steady state solution (u s , v s , w s ) in the positive octant. From (11), θ > 0, so B 0 > 0, and since B 2 > 0, (10) guarantees (u s , v s , w s ) is stable. Note that the case where (ū s ,v s ) is unstable places conditions on a 1 and a 2 such that r − a 1 a 2 is negative. However, r −a 1 a 2 can be negative for stability of the coexistence steady state (u s , v s , w s ).

Now let
One can check that A is nonempty with (0.2, 0.8, 0.5, 0.1, 0.96, 0.5, 0.1) as an element of A. With these parameters, the coexistence steady state is (u s , v s , w s ) = (0.5, 0.5, 0.8). Hence, with parameters satisfying the above assumptions in the set A express the predator-mediated coexistence situation we consider.  If (a 1 , a 2 , b 1 , b 2 , r, µ, ) ∈ A, then (u s , v s , w s ) is the globally stable positive coexistence solution to (9) for positive initial conditions.
Remark. We assume (10), (11), and (12) hold and restrict ourselves to analysis for χ > 0 in the rest of the paper. Biological interpretation of these assumptions requires returning to the dimensional form of system (1). For example the nondimensional assumption of r − a 1 a 2 is not too negative in dimensional form is: . This states that the relative intrinsic growth rate of species v is sufficiently large relative to the strength of competition and carrying capacities. Similarly, assumption (11) is a statement that the relative growth rate of species v and the predation conversion rates for u and v are sufficiently high for the predator to be able to both control the size of the prey species population and the prey species are able to grow sufficiently fast to not go extinct [19].

Pattern formation.
Pattern formation refers to the process that by changing a (bifurcation) parameter, the spatially homogeneous steady state loses stability to spatially inhomogeneous perturbations, and stable inhomogeneous solutions arise.
Here we want to investigate pattern formation for our system, particularly with regard to the prey-tactic sensitivity parameter, χ. We also investigate possible patterns when random dispersal (diffusion) parameter, d, is changed. , t), and w = w s + W (x, t), and substitute into the u, v, w system and linearize; then we have We follow a stability analysis like that given in [31], namely consider (U, V, W ) T = e λt+i k· x Ξ, k = | k|, then we obtain the Jacobian matrix The dispersion relation associated with our linearized system is det(J ) = 0, that is, where k is the wave number, and λ is the eigenvalue, which determines the temporal growth and depends on the wave number. Here The dispersion relation (13) gives λ(k) as a function of wave number k at the steady state (u s , v s , w s ). Now, Re λ(k) < 0 if, and only if a(k 2 ) > 0, c(k 2 , χ) > 0, and a(k 2 )b(k 2 , χ) − c(k 2 , χ) > 0. For the steady state to be unstable we need Re λ(k) to be positive for some sufficiently large range of k = 0. With a(k 2 ) > 0, one way for this to happen is if c(k 2 , χ) < 0. Define By (11) p b (r) > 0; however in the v "wins" case for 0 < d < 1, p a (r) can be negative. We assume for this exposition that p a (r) > 0. Thus, c(k 2 , χ) < 0 if χ > χ c , where By (10) χ c > 0 for k 2 large enough, that is ). The potential for k 2 min to be positive implies that pattern formation may potentially be possible in the χ < 0 (attracting case). From now on we only consider k 2 > k 2 min with χ > 0 and will comment later on the corresponding case for χ < 0.

4.3.
Bifurcations with prey-tactic sensitivity χ. For χ = 0, c(k 2 , χ) remains positive because of (11); that is, given parameters in A. Then, as χ is increased, c(k 2 , χ) will become negative (χ > χ c ). Figure 1 shows an example of the graph of χ c (red curve) in the k 2 , χ plane. For values of (k 2 , χ) above (respectively, below) the graph of χ c , c(k 2 , χ) < 0 (respectively, c(k 2 , χ) > 0). For k 2 min = vs(a1b2−rb1) b1d > 0 (so r < a 1 b 2 /b 1 ), then and for k 2 min = 0, either r = a 1 b 2 /b 1 , or r > a 1 b 2 /b 1 . In the first case and in the second case In either case for k 2 1, χ c ∼ d b1usws k 2 as seen in figure 1. Pattern formation can also potentially occur when a(k 2 )b(k 2 , χ)− c(k 2 , χ) < 0. The blue dashed line curve in figure 1 gives values of (k 2 , χ) where this expression is zero and will be negative (positive) for values lying above (below) this curve. Note that since a(k 2 ) > 0, whenever c(k 2 , χ) < 0 there will be a real positive root of the dispersion relation even if a(k 2 )b(k 2 , χ) − c(k 2 , χ) < 0. However when a(k 2 )b(k 2 , χ) − c(k 2 , χ) < 0 and c(k 2 , χ) > 0, there may not be any roots of the dispersion relation with positive real part leading to the potential that no pattern formation is possible in this case.
Defining χ 0 = min k>kmin χ c (k 2 ) then χ > χ 0 does not guarantee pattern formation since the domain is bounded; the allowable wave numbers are discrete. But, from the asymptotics, for a fixed χ > χ 0 , the interval where c(k 2 , χ) < 0 has ap- min . Let (k 2 1 , k 2 2 ) be such an interval; that is, (k 2 1 (χ), k 2 2 (χ)) = {k 2 |c(k 2 , χ) < 0}. Now consider, for a convenient, though generalizable case, the one-dimensional domain Ω = (0, L). Then the wave numbers correspond to the eigenvalues of the operator − with homogeneous Neumann boundary conditions, namely, k = nπ/L, where n = 0, ±1, ±2, . . .. Identify k 1 with n 1 π/L, and k 2 with n 2 π/L, for some real numbers 0 < n 1 < n 2 . Then we need χ large enough for the interval to contain a discrete unstable mode; that is, for n 1 ≤ n ≤ n 2 . Note that n 1 is monotone decreasing, and n 2 is monotone increasing with respect to χ, as long as χ > χ 0 . χ can be increased to the point where either n 1 = n 1 (χ) or n 2 = n 2 (χ), or both, are integers. Let χ m be the smallest such χ bigger than χ c such that condition holds. Then this can be our bifurcation value.
A similar analysis follows for the attractive (χ < 0) case when k 2 min > 0 where 0 > χ c > −k 2 min . We summarize these observations in the following (preytaxisinduced instability) statement.
Theorem 4.2. Assume (a 1 , a 2 , b 1 , b 2 , r, , µ) ∈ A. In the absence of our taxis mechanism, that is, χ = 0, the steady state (u s , v s , w s ) is locally stable. Let χ m > χ c be the first value of χ such that either k 1 L/π or k 2 L/π is an integer. Then χ m is a bifurcation number and χ > χ m is a necessary and sufficient condition for pattern formation of system(2), (3).
Note that the repulsive case requires a sufficiently strong prey defense (χ > χ c > 0) for pattern formation to develop. In the attractive case, however, the attraction needs to be sufficiently weak (0 > χ > −χ c ) for pattern formation to develop.

4.4.
Bifurcations with diffusivity parameter d. Here we focus on the repulsive, χ > 0, case and provide an alternative way to consider c(k 2 , χ) is (11). That is, Note that for all k 2 , χ c1 (k 2 , r) > 0, so that A(χ, k 2 ) < 0 if, and only if χ > χ c1 (k 2 , r). The sign of χ c2 (k 2 , r) depends on the sign of rb 1 − a 1 b 2 and the value of r − a 1 a 2 , given (11), along with the range of k 2 . If, for example, rb 1 − a 1 b 2 = 0, then B is independent of χ and depends on the sign and size of r − a 1 a 2 ; that is, Another case is when rb 1 − a 1 b 2 > 0 and r − a 1 a 2 > 0, where figure 2 shows a case where χ c1 and χ c2 cross at somê k 2 , providing four regions of (k 2 , χ) space. Region I (respectively, region II) is for χ c1 < χ < χ c2 , for k 2 <k 2 (resp., for χ > χ c1 , χ c2 ). In both regions (u s , v s , w s ) is unstable for 0 < d < d + , where d + is the largest root of the quadratic (with respect to d), c(k 2 , χ) = 0. In region III, where k 2 >k 2 , and χ c2 < χ < χ c1 , (u s , v s , w s ) is unstable for 0 < d − < d < d + (d ± are the two real, positive roots of the quadratic). In region IV, where χ < χ c1 , χ c2 , (u s , v s , w s ) is always stable. However, there are a number of cases to consider to determine whether a range of d exists giving c(k 2 , χ) < 0. This depends on the size of χ and conditions on k 2 , and we leave details to the interested reader.
One approach is to consider a roots argument. Principle cases would be (i) B > A 2 /4k 6 ; (ii) 0 < B < A 2 /4k 6 ; (iii) B < 0. For case (i) the roots to c(k 2 , χ) = 0, namely d ± , are complex and so c(k 2 , χ) > 0 for all arguments. So consider the other two cases in more detail.
• Case (ii) (0 < B < A 2 /4k 6 ): Both roots are real but their sign depends on the sign of A. If χ < χ c1 , then A > 0, so d ± < 0 and this implies c(k 2 , χ) > 0 for all d > 0. Respectively, if χ > χ c1 , then A < 0, so d ± > 0 and there exists an interval, 0 < d − < d < d + , of length where c(k 2 , χ) < 0. This interval must be long enough to contain at least one unstable eigenvalue. Write the discriminant in the form (14) implies χ must be sufficiently large. • Case (iii) (B < 0): Both roots are real but are of opposite signs. Then c(k 2 , χ) < 0 for 0 < d < d + = (−A + √ A 2 − 4k 6 B)/2k 6 , and d + must be large enough to contain at least one unstable eigenvalue. This again requires χ to be sufficiently large.

4.5.
Global stability of the co-existence steady state. Since, by above analysis, we have L ∞ norms on u, v, w that are bounded, then we have a constant upper solution to the system (2). Note that (u, v, w) = (0, 0, 0) provides a lower solution for the system. Hence, we have is a positive invariant set for system (2) for some finite constant M > 0 whose existence is guaranteed by Theorem 3.5. Let V is a Lyapunov functional (see, e.g. [26,38]). If (u, v, w) is a solution to (2), then So, for notational convenience, let f (u, v, w) .
Note that with similar expressions for the corresponding v and w terms. So, Note that by adding and subtracting the appropriate quantities, we have where we have used the fact that f (u s , v s , w s ) = g(u s , v s , w s ) = h(u s , v s , w s ) = 0. Hence, r then we can write Since M is real symmetric, I 2 will be negative if M is positive definite; this will be so if both eigenvalues are positive. Now tr(M ) = 1+r > 0, and det(M ) = r−( a1+a2 2 ) 2 . The eigenvalues of M are guaranteed positive if det(M ) > 0, that is, r > ( a1+a2 2 ) 2 . Now consider I 1 . Let us first motivate our calculation by considering the one dimensional case. Let The twodimensional case will carry over to arbitrary dimensions. Let Y 2 be the vector with components u y , w y , then in the two-dimensional case we can write As in the case with matrix M , we can guarantee I 1 is negative if N has only positive eigenvalues. Note that the tr(N ) = d us u 2 + ws w 2 > 0 and det(N ) = d usws (uw) 2 − χ 2 4 ( ws w ) 2 . Therefore, the det(N ) is positive if χ 2 < 4d u 2 us ws . If the uniform bound on u is given by u ≤ 1, then a sufficient condition on the positivity of eigenvalues would be χ 2 < 4dus ws . Given the conditions for positive determinents of matrices M, N , the only possibility for d dt V (u, v, w) = 0 is when (u, v, w) = (u s , v s , w s ). By the LaSalle invariance principle , we obtain the global stability of (u s , v s , w s ). Theorem 4.3. Assume (u 0 , v 0 , w 0 ) ∈ H. If r > (a 1 + a 2 ) 2 /4, and χ 2 < 4du s /w s hold, then (u s , v s , w s ) is globally asymptotically stable.
5. Numerical simulation . In this section we describe a numerical scheme for simulation of model system (2) -(4) and provide numerical simulations to illustrate some of the pattern formation possible from small perturbations of the homogeneous predator-mediated coexistence steady state using the prey-tactic sensitivity parameter, χ, as a bifurcation parameter. We concentrate on the cases where in the pure competition system in the absence of a predator either u "wins" or v "wins" and the predator mediated coexistence is expected in the presence of the predator in the absence of a repulsive prey defense mechanism.
We implement a discontinuous Galerkin method described in section 5.1 for a one dimensional spatial domain using the FeNiCs software package [3,29]. Although we have performed simulations in two dimensional spatial domain, we have not learned anything particularly new over the one dimensional case, so simulations shown in section 5.2 are all one dimensional spatially. Also, in all simulations the initial data take the form u 0 (x) = u s +ε sin(nπx/L+φ), where for the given parameter set, u s is the u component of the constant steady state, as given in (8), and the perturbation amplitude ε is 10% of the u s value, and similarly for v 0 (x) and w 0 (x). We have generally kept n = 1 or n = 2.
5.1. Discontinuous Galerkin method . We implement a discontinuous Galerkin method that is derived from schemes that are well suited for advection-diffusionreaction equations [6,14,8]. Methods such as these that are well suited for strongly convective systems should also be reasonable methods for quasi-linear parabolic systems such as system (2)- (4).
For ease of presentation in stating our discontinuous Galerkin method we begin by looking at just the equation for the predator, w, in steady state where h(u, v) = µ − εw + b 1 u + b 2 v. We will state the bilinear forms that are used in seeking a weak solution to (15) then apply these forms to the full system (2)-(4) in steady state and ultimately place in a time stepping scheme for numerical solution of the original system. The corresponding weak formulation to (15) is The discrete formulation is to find in a broken polynomial space defined over the For a continuous finite element method, one could use the bilinear forms a diff (w h , ξ h ) and a adv (w h , ξ h ; ∇ h u h ) as defined by (16). For the discontinuous Galerkin method illustrated here, the diffusion term is implemented with a symmetric interior penalty method for diffusion a diff (w, ξ) = a SIP (w, ξ) and the advection term is implemented with an upwind scheme a adv (w, ξ; u) = a upw (w, ξ; u) given by and where E o h is the set of all interior edges for the mesh, |e| is the mesh element diameter, α > 0 is a penalty parameter, and n + is the unit normal to e such that (−χ∇ h u h ) · n + > 0. [[·]] and {·} are the jump and average operators respectively.
Returning to the steady state version of the full system (2) the discrete formulation becomes to find functions To obtain a numerical solution to the original system (2) -(4) we implement a backward Euler method. Rewriting (17) , for a given temporal step size ∆t.

5.2.
Case studies in one dimension . Pattern formation may occur when the coefficient c(k 2 , χ) < 0 in the dispersion relation (13). Theorem 4.2 gives conditions for pattern formation to occur when c(k 2 , χ) < 0. In figures 3 and 4 we exemplify stationary spatial pattern formation for cases where u "wins" and v "wins" respectively. In figures 5 and 6 we exemplify spatio-temporal pattern formation for cases where u "wins" and v "wins" respectively. In all cases χ is chosen so that c(k 2 , χ) < 0 when k = π/L or 2π/L where L = 5 is the length of the domain. Figure 3 shows a case where u "wins" in the pure competition setting. In the absence of an effective repulsive prey defense mechanism, the expected homogeneous steady state is (u s , v s , w s ) = (0.15, 0.48, 0.97). Note that in making an initial condition that is a small periodic perturbation of this coexistence state, u starts at an extreme disadvantage as u 0 (x) is much smaller than both v 0 (x) and w 0 (x). The predator is locally repelled by u, which allows for the competition dynamics to dominate locally. Hence, coexistence is disrupted and v is driven to extinction : v "wins": In this example the parameters are chosen so that in the pure competition system v "wins". Note that while parameters are such that in the absence of a repulsive prey defense mechanism by u, a predator mediated coexistence is expected; yet, once this state is disrupted, the prey coexist with peak densities out of phase with the predator in the presence of this taxis mechanism. Parameter set for this figure is: (a 1 , a 2 , b 1 , b 2 , r, ε, µ, χ, d) = (1.1, 0.1, 1.4, 0.3, 0.4, 0.4, 0.2, 50, 0.25) by u. Then the resultant pattern formation is what would be expected for the predator-prey (u, w) system. In this simulation χ is just enough larger than χ c so that one unstable mode is captured. Figure 4 shows a case where v "wins" in the pure competition setting. In the absence of an effective repulsive prey defense mechanism, the expected homogeneous steady state is (u s , v s , w s ) = (0.02, 0.10, 0.63). Note that in making an initial condition that is a small periodic perturbation of this coexistence state, u starts at an extreme disadvantage as u 0 (x) is much smaller than both v 0 (x) and w 0 (x). However, the repulsion of the predator by u does not result in v driving u to extinction as it tries to "win" the competition. Rather, in this case as v out competes u, the ability of u to repel w is diminished and the two prey species coexist with peaks in their respective densities that are out of phase with the predator density peak.
In figures 5 and 6 we revisit the v "wins" and u "wins" cases respectively for parameters where we have observed spatio-temporal patterns. In each figure the first three rows show surface plots of u(x, t), v(x, t), and w(x, t) respectively. The first column shows the initial stages where a non-homogenous spatial pattern appears to be emerging. The second column shows the prey and predator densities at a later time once they have settled into a steady spatio-temporal pattern. The last row of the figures shows the densities after the pattern has emerged at a fixed time in the first column and for a fixed position and short time period in the second column. In both cases the competing species coexist and are out of phase with the predator density as in figures 3 and 4. Note that in the u "wins" case of figure 6 that v is not Figure 5. Spatio-temporal patterns: v "wins" case. The first column shows the initial dynamics from the perturbation of the homogeneous steady state and the second column shows the spatiotemporal pattern that is achieved. In the last row u(x, t), v(x, t), and w(x, t) are given by the blue, green, and red lines respectively. The parameter set for this figure is: (a 1 , a 2 , b 1 , b 2 , r, ε, µ, χ, d) = (1.25, 1.0, 0.25, 0.5, 1.25, 0.5, 0.25, 25, 0.1) driven to extinction. Here the growth rate, r, for v is relatively high and hence v is able to survive in the regions where the predator density is large. Then v diffuses into the regions where thte predator density is low and maintains the competition with u.
We have commented that in our model there is potential for pattern formation in the case of attracting (χ < 0) prey-taxis. Figure 7 shows an example of pattern formation in the case where v "wins". However there is a notable difference in the steady state from the repulsive case (χ > 0) shown in figure 4. In this example, Figure 6. Spatio-temporal patterns: u "wins" case. The first column shows the initial dynamics from the perturbation of the homogeneous steady state and the second column shows the spatiotemporal pattern that is achieved. In the last row u(x, t), v(x, t), and w(x, t) are given by the blue, green, and red lines respectively. The parameter set for this figure is: (a 1 , a 2 , b 1 , b 2 , r, ε, µ, χ, d) = (0.9, 1.1, 0.2, 0.2, 1.09, 0.5, 1.0, 16, 0.25) the prey species u attracts the predator speices w resulting in u and w peaking at the same location. However, in area where u is low and hence so is w, the other prey species v is able to thrive and reaches its peak out of phase from u and w. Thus a predator that forages for a prey that is the weaker competitor, the stronger competitor is locally benefited. 6. Conclusions. We have considered the interaction of predator mediated effects of predator mediated coexistence and prey-taxis with an emphasis on prey defense Figure 7. Attracting (χ < 0) case: In this example the parameters are chosen so that in the pure competition system v "wins". Note that while parameters are such that in the absence of a attracing prey-taxis mechanism by u, a predator mediated coexistence is expected; yet, once this state is disrupted, the prey coexist with peak densities out of phase with the predator in the presence of this taxis mechanism. Parameter set for this figure is: (a 1 , a 2 , b 1 , b 2 , r, ε, µ, χ, d) = (1.25, 1.00, 0.25, 0.50, 1.25, 0.5, 0.25, −0. 1,10) in a three species model consisting of two prey species in competition sharing a common predator species. While other studies have investigated some "tactic" mechanism such as chemotaxis or prey-taxis in predator-prey or two species competition systems, the work here is the first the authors know of incorporating a prey-taxis mechanism for one prey species in a competing species environment with a common predator. The taxis "breaks" the symmetry of behavior that would normally exist in an ecosystem of prey with a common predator, leading to qualitative differences between v "wins" versus u "wins" scenarios. While we treat the mechanism as a generic repulsive form of "prey-taxis," this mechanism may be representative of any number of prey behaviors including release of chemical odors, loud noises, or mimicry of a larger predator to give a few examples. We have shown that the model has non-negative, bounded, global in time solutions and demonstrated the global asymptotic stability of the homogeneous coexistence steady state. We have primarily focused in this report on pattern formation, treating the prey-tactic sensitivity parameter as a bifurcation parameter, and have illustrated this analysis with numerical simulations.
The ability of a predator to establish a coexistence between competing prey is long known from both theoretical mathematical modeling studies and experimental field studies. However the impact of behavioral responses of predator and prey on predator mediated coexistence is little studied. Here we have looked explicitly at cases where in the presence of a predator, one of the prey engages a defense mechanism that elicits an avoidance behavior in the predator. We observe in the presence of a generalist predator and prey defense mechanism a difference in coexistence behavior between the competitors that depends on the expected winner of the competition. When the weaker competitor is provided the defense mechanism, we see a predator mediated coexistence arise between the competitor prey. Moreover, the weaker competitor provides defense for the stronger competitor in the sense that the prey species coexist out of phase from the predator population. However when the stronger competitor is also able to defend against the predator, we see a lack of coexistence between the competitors; that is, the stronger competitor continues to win the competition. We also observe that the surviving prey and predator densities are out of phase of one another. While not shown, we see similar results with a specialist predator that preys only on the species with defense mechanism.
Previous literature indicates that in predator-prey and two species competition systems that pattern formation requires a repulsive mechanism. As such we did not expect to find pattern formation for χ < 0. Surprisingly, we found that in cases where v "wins" in the competition setting and u is weakly attracting to the predator, then pattern formation is possible. Through numerical experimentation we did not find cases for u "wins" where pattern formation occurs and conjecture that for χ < 0, a 1 < 1, and a 2 > r pattern formation is not possible.
A predator defense mechanism such as releasing a chemical odor that is toxic or otherwise alarming, or attracting to a predator would follow a chemotaxis type model where the chemical diffuses through the environment after release by the prey species. The prey-taxis type model presented here can be considered as an approximation to the chemotaxis type model where the chemical is assumed to act quickly so a quasi-steady state assumption where the chemical does not diffuse can be employed. It is a separate consideration beyond the scope of this paper to qualitatively compare the differences between this prey-taxis type model and its chemotaxis type model variant where the chemical acts quickly leading to very little diffusion.
The observed pattern formation in the numerical simulations gives rise to a number of interesting questions. One question is the determination of the requirements on the parameters for stationary non-homogeneous patterns and spatio-temporal patterns. However, such a detailed bifurcation analysis is beyond the scope of this initial exploration of the interaction of predator mediated effects and is reserved for future work. Additionally, we are intrigued by the role of the growth rate, r, of the prey without the defense mechanism on the dynamics of pattern formation. We conjecture that values of r separate the formation of stationary patterns from spatio-temporal patterns where a relatively large (small) value of r is required for the formation of spatio-temporal (stationary) patterns as the relatively rapid growth of v will aid it in competing with u which may diminish the repelling gradient of u allowing w to invade and decrease v and allow the repelling gradient of u to grow.