Stationary and time periodic patterns of two-predator and one-prey systems with prey-taxis

This paper concerns pattern formation in a class of reaction-advection-diffusion systems modeling the population dynamics of two predators and one prey. We consider the biological situation that both predators forage along the population density gradient of the preys which can defend themselves as a group. We prove the global existence and uniform boundedness of positive classical solutions for the fully parabolic system over a bounded domain with space dimension $N=1,2$ and for the parabolic- -parabolic-elliptic system over higher space dimensions. Linearized stability analysis shows that prey-taxis stabilizes the positive constant equilibrium if there is no group defense while it destabilizes the equilibrium otherwise. Then we obtain stationary and time-periodic nontrivial solutions of the system that bifurcate from the positive constant equilibrium. Moreover, the stability of these solutions is also analyzed in detail which provides a wave mode selection mechanism of nontrivial patterns for this strongly coupled system. Finally, we perform numerical simulations to illustrate and support our theoretical results.


Introduction
One of the central problems in the study of ecological systems is to understand the spatialtemporal behaviors of population distributions of interacting species. Over the past few decades, many mathematical models have been proposed and developed to investigate the collective influence of individual's dispersals on the spatial-temporal population distribution at a group level. For example, reaction-diffusion equations have been applied to model the densities and spatial distributions of single or multiple interacting species in [13,39,40,44] etc. In particular, the formation of nontrivial patterns in these systems such as diffusion-driven heterogeneity or Turing's pattern, especially those with concentration properties, can be used to describe the aggregation or segregation phenomena.
In this paper, we study the following 3 × 3 reaction-advection-diffusion system x ∈ Ω, where Ω is a bounded domain in R N , N ≥ 1 with smooth boundary ∂Ω and unit outer normal n. ∇ is the gradient operator and ∆ is the Laplace operator. u, v and w are functions of space-time variable (x, t). d i , α i , i = 1, 2, 3, and β i , β 3i , i = 1, 2, are positive constants. χ and ξ are also assumed to be positive constants and φ(w) is a smooth function.
(1.1) models the population dynamics of three interacting species subject to Lotka-Volterra kinetics, where u and v are population densities of two predators at location x and time t, and w is the population density of the prey. It is assumed that both predators consume the same prey species and disperse over the habitat by a combination of random diffusion and directed movement (prey-taxis) along the gradient of prey population density, while the preys move in the habitat randomly. Diffusion rates d i , i = 1, 2, 3, measure the intensity of random dispersals of the species. Here prey-taxis is the phenomena that predators with the ability to perceive the heterogeneity of prey distribution approach the patches with high preys density. The positive constants χ and ξ measure the intensity of the directed movement of each predator in response to the prey-taxis. The prey-dependent function φ reflects the strength of prey-tactic movement of the predators with respect to the variation of prey population density. Various specific forms of φ can be chosen, depending on the biological situation that one tries to model. In particular, χuφ(w) > 0 represents the biological situation that predator u forage the preys while χuφ(w) < 0 means that predator u retreat the habitat of the preys which can defend themselves as a group when the population density is high. See [15,38,64] and the references therein for more examples and further discussions about antipredation of preys through group defense. The population kinetics are of classical Lotka-Volterra type, where α i are the intrinsic growth rates of the species which reproduce logistically, β i are the growth rates of the predators and β 3i are the death rates of the preys due to predation.
In search of preys, prey-taxis is the process that predators move preferentially towards patches with high density of prey. All predators forage preys in accordance with spatial distributions of prey density. Though individual predator tends to forage the vicinity of recent captures, swarms of predators exhibit prey-taxis behaviors. This uneven searching effort may result in both larger predation rate and aggregation of predators where preys are abundant. For example, insects can find preys which are concentrated in some areas and sparse in others, though they lack long-distance sensory perception. See [26] and the references therein for detailed discussions and field experiments. This mechanism is the same as bacterial chemotaxis in which cellular organisms sense and move along the concentration gradient of stimulating chemicals in the environment [18,19,20,28].
Various reaction-diffusion systems have been proposed to describe spatial predatorprey distributions under directed movements. In [26], Kareiva and Odell proposed a mechanistic approach, formulated as partial differential equations with spatially varying dispersals and advection, to demonstrate and explain that area-restricted search does create predator aggregation. Since then, various reaction-diffusion systems have been proposed to model predator-prey dynamics with prey-taxis. Sapoukhina et al. [49] assumed that such directed movement is determined by the velocity variation (i.e. the acceleration). They investigated the effect of prey-taxis on the predator's ability to maintain pest population below a certain economic threshold value. It is assumed in [12,14,56] etc. that the directed movement of predator is due to the advective velocity. We refer to [18,28,46,57] etc. for the derivation or justification of (1.1).
Reaction-diffusion models with prey-taxis subject to different population dynamics have been studied by various authors extensively. Ainseba et al. [2] established the existence and uniqueness of weak solutions to prey-taxis models with volume-filling effect in prey-tactic sensitivity function. Global existence and boundedness of classical solutions are obtained in [16,53], while nonconstant positive steady states are investigated in [58, 59,61] and [34] using Crandall-Rabinowitz bifurcation theory and fixed point theory, respectively. Lee et al. [33] studied pattern formation in prey-taxis system under a variety of nonlinear functional responses. Travelling wave solutions for prey-taxis models are studied by the same trio in [32]. Effects of prey-taxis on predator-prey models are investigated numerically in [7] which suggest that both response functions and initial data play important roles in pattern formation; moreover, predator-prey models admit chaotic patterns when the prey-taxis coefficient is sufficiently large. We also want to mention that for predator-prey systems without prey-taxis, Okubo and Levin [44] pointed out that an Allee effect in the functional response and a density-dependent death rate of the predator are necessary for the pattern formation. See the book of Murray [40] for more works on prey-taxis models. We also want to mention that there are also many works on the predator-prey models with density-dependent diffusion or the so-called cross-diffusions, such as [29,30,43,48,65] etc.
All the aforementioned works are devoted to studying prey-taxis models with one predator and one prey. There are some works on two-predator and one-prey model without prey-taxis [22,37,45,55]. Moreover, we refer the reader to [23,51,57] etc. for the studies on the dynamics of competitive reaction-diffusion systems with cross-diffusion or advection. In [35], Lin et al. considered (1.1) with χ = ξ = 0 over Ω = (0, L) with or without diffusion. They investigated the global dynamics of all equilibria of the ODE system and obtained traveling wave solutions of the PDE system. In this paper, we consider system (1.1) with two predators and one prey, both predators foraging the preys prey-tactically. We are motivated to study the effect of prey-taxis and group defense on the spatial-temporal dynamics of (1.1), in particular, its positive solutions that exhibit stationary or time-oscillating spatial structures.
There are several scientific goals of our paper and the rest part of this paper is organized as follows. In Section 2, global existence of positive classical solutions are obtained for the full parabolic system in Theorem 2.6 over 1D and 2D bounded domains and for the parabolic-parabolic-elliptic system over arbitrary space dimensions in Theorem 2.7 respectively. We also prove that these classical solutions are uniformly bounded in time.
In Section 3, we study the linearized stability of the positive equilibrium to (1.1). It is shown that prey-taxis destabilizes the positive equilibrium if there is group defense in preys and it stabilizes the equilibrium otherwise. Section 4 and Section 5 are devoted to the steady state and Hopf bifurcation analysis of (1.1) over Ω = (0, L). Existence and stability of stationary and time-periodic spatial patterns are established in Theorem 4.2, Theorem 4.3 and Theorem 5.2. Though there have been many works devoted to the study of nonconstant steady states of the 2×2 prey-taxis systems, regular time-periodic spatial patterns are quite new to the literature. Finally, we perform extensive numerical studies in Section 6 to illustrate and support our theoretical findings. We would like to note that, C and C i are assumed to be positive constants that may vary from line to line in the sequel.

Existence and boundedness of global solutions
In this section, we study the global existence and boundedness of classical positive solutions to (1.1). There are two well-established methods in proving the global boundedness for reaction-diffusion systems in the literature. One is to construct its time-monotone Lyapunov-functional and the other is to use Gagliardo-Nirenberg type estimate through Moser-L p iteration. As we shall see in our mathematical analysis and numerical simulations, (1.1) admits time-periodic spatially inhomogeneous solutions for properly chosen parameters and hence lacks time-monotone Lyapunov-functional and maximum principle. Our proof of global existence and boundedness of classical positive solutions to (1.1) is based on the local theory of Amann [4,6] and the Moser-Alikakos L p iteration technique [3].

Local existence and preliminary results
We first obtain the local existence and uniqueness of positive classical solutions to (1.1) and their extensibility criterion based on Amann's theory [6]. To this end, we convert it into the following triangular form subject to nonnegative initial data u 0 , v 0 , w 0 and homogeneous Neumann boundary conditions. Since all the eigenvalues of D 0 are positive, system (2.1) is normally parabolic, and we have from the standard parabolic maximum principles that u, v, w ≥ 0 in Ω × R + . Moreover, the following results are evident from Theorem 7.3 and Theorem 9.3 of [5], Theorem 5.2 in [6] and the standard parabolic regularity arguments.
Let Ω be a bounded domain in R N , N ≥ 1 with smooth boundary ∂Ω. Assume that d i , α i for i = 1, 2, 3, and β i , β 3i for i = 1, 2 are positive constants. Suppose that the initial data (u 0 , v 0 , w 0 ) ∈ C 1 (Ω) × C 1 (Ω) × W 2,p for some p > N , and According to Theorem 2.1, the local solutions (u, v, w) are global if their L ∞ -norms are bounded in time. To establish the global existence, we state some basic properties of the local solutions.
over Ω × (0, T max ). Then there exist positive constants C 1 and C 2 such that Proof. We first prove (2.3). Letw(t) be the solution of then it is easy to see thatw(t) is uniformly bounded and it is a super-solution to the third equation of (1.1). Hence we have that w(x, t) ≤w(t) from the maximum principle and this gives rise to (2.3). Moreover, if w 0 (x) ∈ (0, 1),ŵ(t) ≡ 1 is a super-solution and this implies that w(x, t) ∈ (0, 1). We next prove the boundedness of u(·, t) L 1 and the same argument applies for v(·, t) L 1 . Integrating the first equation in (1.1) over Ω, we have from the boundedness of w(·, t) L ∞ that here and in the rest of this section we skipped dx. In light of ( Ω u(x, t)) 2 ≤ |Ω| Ω u 2 (x, t), we have Solving this ordinary differential inequality gives In order to estimate L p -norms of u and v, we need to provide an a priori estimate on ∇w L ∞ . The following lemma is due to the well-known smoothing properties of operator −d 3 ∆ + 1 and embeddings between the analytic semigroups generated by {e t∆ } t>0 . We refer the reader to [21,63] for references. Lemma 2.3. Suppose that Ω is a bounded domain in R N , N ≥ 1. Assume that all the conditions in Theorem 2.1 hold and let (u, v, w) be the classical solution of (1.1) over (0, T max ). Then there exists a constant C dependent on w 0 W 1,q (Ω) and |Ω| such that for We first write w-equation into the following variation-of-constants formula (2.7) where ν is the first Neumann eigenvalue of −d 3 ∆ + 1. On the other hand, we know from the gamma function that therefore (2.5) follows from (2.7).
In order to apply the Moser-iteration, it is sufficient to prove the boundedness of ∇w(·, t) L 2(N +1) thanks to the quadratic decay kinetics in (1.1).

Proposition 2.1.
Let Ω be a bounded domain in R N , N ≥ 1, and (u, v, w) be the classical solutions of (1.1) over (0, T max ). If there exists a constant C 1 such that Therefore (u, v, w) is a global solution to (1.1) over Ω × (0, ∞).
Proof. For each p ≥ 1, we have from straightforward calculations To estimate the second integral in (2.8), we have from Young's inequality that, for any constant > 0, there exists C( ) > 0 such that On the other hand, for each p > 0 there exists C 0 > 0 such that For the simplicity of notations and without loss of our generality, we assume that φ(w) L ∞ ≤ 1 in light of (2.3). Choosing p = N in (2.8), thanks to the boundedness of ∇w(·, t) L 2(N +1) we have where is chosen to be sufficiently small. Then we can have from (2.9) that u(·, t) L N ≤ C, ∀t ∈ (0, T max ) and this, in light of Lemma 2.3, implies that ∇w(·, t) L q ≤ C for any q ∈ [1, ∞).
In particular, we note that ∇w(·, t) L 2(N +2) is bounded, therefore we can choose p = N + 1 in (2.8) and perform the same calculations as in (2.9) to show that u(·, t) L N +1 is uniformly bounded. Once again by applying Lemma 2.3, we can show that ∇w(·, t) L ∞ is bounded.
Without losing the generality of our analysis, we assume that φ(w)∇w L ∞ ≤ 1, and then we have from (2.8) by taking smaller than d 1 pχ . Letting p go to ∞ and applying the standard Moser-Alikakos L p -iteration [3] on (2.10), we can show that there exists a constant C > 0 such that u(·, t) L ∞ ≤ C, ∀t ∈ (0, T max ). Similarly, we can prove the uniform boundedness of v(·, t) L ∞ . The proof of this proposition completes.
2.2 A priori estimates of fully parabolic system for N ≤ 2 According to Lemma 2.3 and Proposition 2.1, in order to obtain the global existence of (1.1), it is sufficient to establish the boundedness of u and v in their L p -norms for some p ≥ N . Lemma 2.2 already gives L 1 boundedness from which global existence of (1.1) in 1D follows. In this section, we restrict our attention to study (1.1) over 2D and we want to prove the boundedness of (u, v)(·, t) L 2 . First of all, we introduce several entropy-type inequalities to estimate weighted functions u ln u L 1 , v ln v L 1 and ∇w L 2 , etc. The a priori estimates rely on Gagliardo-Nirenberg interpolation inequality for which [41,54] are good references. We now give the following lemma. Proof. Using u-equation of (1.1) and the boundedness of φ(w) L ∞ , we have from integration by parts and Young's inequality where C 1 , C 2 and C 3 are positive constants. Similarly we obtain from the v-equation On the other hand, we have from straightforward calculations where C i 's are positive constants. In light of (2.3), we have from Gagliardo-Nirenberg interpolation that there exists C * > 0 such that Multiplying (2.15) by µ 1 and then adding it to (2.12) and (2.13), we have where C i are positive constants. Choose µ 1 large such that d 3 µ 1 2 ≥ C * (C 2 + C 4 + 1). We apply the fact that s ln s ≤ s 2 ln s + C , ∀s > 0, and Young's inequality to obtain Proof. We have from the integration by parts, Hölder's inequality and (2.3) that where to derive the inequality we have assumed that φ(w) L ∞ , φ (w) L ∞ ≤ 1 without loss of our generality.
To estimate (2.18), we apply Lemma 3.5 of [41] with p = 3 to have that for any . This fact and (2.11) imply that Thanks to the Gagliardo-Nirenberg inequality [54] and (2.11), we have and and that therefore (2.18) gives us By the same arguments, we can show that On the other hand, we operate ∇ to w-equation in (1.1). Then it follows from Young's inequality and integration by parts that where we used the inequalities Multiplying (2.24) by µ 2 sufficiently small with = d 3 µ 2 4(2+µ 2 ) , we add it with (2.22) and (2.23) to have d dt (2.25) Therefore (2.17) follows from (2.25) thanks to Grönwall's inequality.

Existence and boundedness of global solutions
We now present our main results on the global existence and uniform boundedness of positive classical solutions to (1.1).
Suppose that all the parameters in (1.1) are positive and φ ∈ C 3 (R; R). Then for any nonnegative initial data Proof. According to Theorem 2.1 and Lemma 2.3, we only need to show that (u(·, t), v(·, t)) L ∞ is uniformly bounded for all t ∈ (0, T max ), then T max = ∞ and the existence part of Theorem 2.6 follows. Moreover, one can apply the standard parabolic boundary L p -estimates and Schauder estimates in [31] to verify that u t , v t , w t and all spatial partial derivatives of u, v and w up to second order are bounded onΩ × [0, ∞), therefore (u, v, w) have the regularities stated in Theorem 2.6 Thanks to Lemma 2.3 and Lemma 2.5, we have that ∇w(·, t) L p is bounded for any p ∈ (1, ∞), then the uniform boundedness of (u, v)(·, t) L ∞ follows from Proposition 2.1. Together with (2.3), we conclude the proof of Theorem 2.6.
Our proof of the boundedness of (u, v) L ∞ hence the global existence of (1.1) is restricted to the domain with dimension N = 1, 2 due to technical reasons. It is well known that, in order to prove the boundedness of ∇w L ∞ , it is sufficient to prove the L p -boundedness of u and v for some p > 2(N + 1) for a wide class of reaction-diffusion system with biological relevance, for instance the Keller-Segel chemotaxis models with no cellular growth. The logistic decay in (1.1) contributes so one only needs to prove L p -boundedness of u and v for some p > N . Whether or not the logistic dynamics are sufficient to prevent finite or infinite time blowups for (1.1) over higher dimensional domains is unclear in the literature and it is beyond the scope of this paper. However, we can show the global existence and boundedness for the following parabolic-parabolicelliptic system of (1.1) in the following system regardless of space dimensions To be precise, we have the following results.
Let Ω be a bounded domain in R N , N ≥ 3 with smooth boundary ∂Ω. Suppose that φ (w) ≥ 0, ∀w ≥ 0. Then (2.26) has a unique classical solution (u, v, w) which exists globally in time and satisfies the following estimate with a positive constant Proof. By the same analysis as in the proof of Proposition 2.1, we shall only need to show the uniform boundedness of u(·, t) L ∞ and v(·, t) L ∞ . Similar to (2.8), we can obtain from the integration by parts where C 1 and C 2 are positive constants that depend on the system parameters and w L ∞ . Applying the Moser-iteration gives rise to the boundedness of u(·, t) L ∞ . By the same way we can prove the boundedness of v(·, t) L ∞ . The proof of this Theorem completes.
The assumption φ (w) ≥ 0 corresponds to the biological situation that the intensity of the directed dispersals of predator species increases as prey density increases, therefore there is no group defense in prey species. We have to leave it for future works on the global existence of (2.26) and its fully parabolic counter-part over higher-dimensions.

Linearized stability of homogeneous equilibrium
From the viewpoint of mathematical modeling, it is interesting and meaningful to investigate (1.1) for its nontrivial solutions describing the spatial distributions of the interacting species over the habitat. We shall see from the mathematical analysis that the formation of such nontrivial patterns is due to the joint effect of prey-taxis and sensitivity function.
To illustrate the effect of prey-taxis on the pattern formation and for the simplicity of calculations, we restrict our attention to (1.1) over one dimension and consider the following system where denotes the derivative taken with respect to x and all the parameters in (3.1) are the same as those in (1.1). We can find that (3.1) has six equilibria and one of them is which is positive if and only if α 3 > β 31 + β 32 . To explore the existence and stability of stationary and oscillatory nonconstant positive solutions to (3.1), our starting point is the linear stability of (ū,v,w) and we shall assume that it is positive from now on. Linearizing (3.1) around the constant equilibrium (ū,v,w) and letting (u, v, w) = (ū,v,w) + (U, V, W ), U, V, W being small perturbations from (ū,v,w), we obtain the following system of (U, V, W ) According to the standard linearized stability principle ([52] e.g.), the stability of (ū,v,w) is determined by eigenvalues to the following matrix We have the following result.
Proposition 3.1. Assume that all the parameters in (3.1) are positive and the condition α 3 > β 31 + β 32 is satisfied. Suppose that φ(w) < 0, then the positive equilibrium (ū,v,w) is locally asymptotically stable if χ < χ 0 and it is unstable if χ > χ 0 , where and where Proof. The characteristic equation for an eigenvalue σ of the stability matrix (3.2) is is unstable if one of the conditions above fails for some k ∈ N + . Note that we always have that η 2 (χ, k) > 0 for each k ∈ N + ; moreover, On the other hand, since φ(w) < 0, it follows from straightforward calculations that , if χ is larger than the minimum of χ S k and χ H k over k ∈ N. Similarly we can show that (ū,v,w) is locally asymptotically stable if χ < χ 0 . This finishes the proof.
According to Proposition 3.1, (ū,v,w) becomes unstable as χ surpasses the threshold value χ 0 if φ(w) < 0 which we shall assume from now on. We would like to point out that biologically φ(w) < 0 describes situation that a huge amount of preys can aggregate to form group defense and keep predators away from the habitat when the prey population density surpassesw. This amounts to a switch from prey-attraction to prey-repulsion in the predation. See [1,15,38,64] and the references therein for the detailed description of prey group defense.
We will see in our coming mathematical analysis and numerical simulations that (ū,v,w) loses its stability to time-periodic patterns through Hopf bifurcation if χ = χ H k and to stationary patterns through steady state bifurcation if χ = χ S k . Therefore we have used the indices H and S to denote Hopf and steady state bifurcation respectively. Here and in the rest part of this paper, we study the effect of prey-taxis on the formation of nontrivial patterns to (3.1). Without losing the generality of our analysis, we treat χ as the variable parameter and fix all the rest parameters, while similarly, we have that Proposition 3.1 also holds for ξ > ξ 0 = min k∈N + {ξ H k , ξ S k }, where ξ H k and ξ S k are functions of χ.
Corollary 1. Suppose that φ(w) > 0 and all the rest conditions in Proposition 3.1 are satisfied, the homogeneous equilibrium (ū,v,w) is locally asymptotically stable if χ >χ 0 and it is unstable if χ <χ 0 , whereχ When φ(w) > 0, both χ S k and χ H k are negative for any k ∈ N + henceχ 0 < 0, therefore (ū,v,w) is always stable according to Corollary 1 since χ > 0. This result corresponds to the widely held belief ( [33] e.g.) that prey-taxis stabilizes homogeneous equilibrium and inhibits the formation of spatial patterns for one-predator and one-prey system. It states that the same holds true for two-predator and one-prey model as long as φ(w) > 0. However, if φ(w) < 0, the prey-taxis destabilizes homogeneous equilibrium which becomes unstable as χ surpasses χ 0 . Therefore, in order to investigate the formation of nontrivial patterns in (3.1), we shall assume that φ(w) < 0 in the coming analysis. Remark 3.1. As we shall see in our coming analysis, the occurrence of Hopf or steady state bifurcation at (ū,v,w) depends on whether χ 0 in (3.3) is achieved at min k∈N + χ S k or min k∈N + χ H k . We divide our discussions into the following cases: (1). If χ 0 = χ S k 0 < min k∈N + χ H k , then η 0 (χ 0 , k 0 ) = 0 and the eigenvalues of (3.2) at χ 0 are Since (ū,v,w) is unstable for all χ > χ 0 , we know that (3.2) has at least one eigenvalue with positive real part when χ = χ S k , k = k 0 . This fact will be applied in the proof of the stability of steady state bifurcating solutions to (3.1) around (ū,v,w, χ S k ). In particular, it implies that the only stable bifurcating solutions must be on the branch around (ū,v,w, χ S k 0 ) that turns to the right while all the rest branches are always unstable. Moreover, Hopf bifurcation does not occur at (ū,v,w) in this case.
In this case, (3.2) has three eigenvalues σ 1 = −η 2 (χ, k) < 0, σ 2 = σ 3 = 0, and linear stability of the (ū,v,w) is lost since there are two zero eigenvalues. This inhibits the application of our steady state and Hopf bifurcation analysis which requires the null space of (3.2) to be one-dimensional. Therefore we assume that χ S k = χ H k , ∀k ∈ N + in our bifurcation analysis.
In general, it is not obvious to determine when case (1) or case (2) in Remark 3.1 occurs. However, if the interval length L is sufficiently small, we have since φ(w) < 0. This implies that for small intervals, χ 0 = min k∈N + χ S k = χ S 1 and (ū,v,w) loses its stability to the steady state bifurcating solution as χ surpasses χ 0 . Since k 0 = 1, the bifurcating solution has a stable wave mode cos πx L which is monotone in x. Moreover, we shall observe from our numerics that k 0 is increasing or nondecreasing in L. Therefore, small domain only supports monotone stable solutions, while large interval supports non-monotone solutions, at least when χ is around

Nonconstant positive steady states
This section is devoted to studying nonconstant positive steady states to system (3.1), i.e., nonconstant solutions to the following system where u, v and w are functions of x and all the parameters are the same as those in (3.1). We assume that α 3 > β 31 +β 32 so that (ū,v,w) is the unique positive equilibrium to (4.1).
In order to look for nonconstant positive solutions to (4.1), we shall perform steady state bifurcation analysis at (ū,v,w). When φ(w) < 0, we already know that prey-taxis χ destabilizes (ū,v,w) which becomes unstable when χ surpasses χ 0 , therefore we are concerned with the conditions under which the spatially inhomogeneous solutions emerge through bifurcation as χ increases. We refer these as prey-taxis induced patterns in analogy to Turing's instability.

Steady state bifurcation
To apply the bifurcation theory of Crandall-Rabinowitz [9,10] with χ being the bifurcation parameter, we introduce the spaces and convert (4.1) into the following abstract equation It is easy to see that F(ū,v,w, χ) = 0 for any χ ∈ R and F : X ×X ×X ×R → Y ×Y ×Y is analytic. Moreover, for any fixed (û,v,ŵ) ∈ X × X × X , the Fréchet derivative of F is given by where We collect the following facts about F.
To seek non-trivial solutions of (4.1) that bifurcate from equilibrium (ū,v,w), we first check the following necessary condition 4) where N denotes the null space. Taking (û,v,ŵ) = (ū,v,w) in (4.3), we see that the null space in (4.4) consists of solutions to the following problem In order to verify (4.4), we substitute the following eigen-expansions into (4.5) t k , s k and r k constants and collect k = 0 can be easily ruled out since α 3 > β 31 + β 32 . For each k ∈ N + , (4.5) has nonzero solutions (t k , s k , r k ) if and only if the coefficient matrix of (4.6) is singular or equivalently where H 1 , H 2 and H 3 are given in Proposition 3.1. Note that χ S k in (4.7) is the same as (3.4). χ S k > 0 if and only if Condition (4.4) is satisfied if χ = χ S k and N (DF(ū,v,w, χ S k )) = span{(ū k ,v k ,w k )} which is of one dimension where and (4.10) Having the potential bifurcation value χ S k in (4.7), we now prove in the following theorem that the steady state bifurcation occurs at (ū,v,w, χ S k ) for each k ∈ N + , which establishes existence of nonconstant positive solutions to (4.1).
Theorem 4.2. Assume that α 3 > β 31 + β 32 and φ(w) < 0. Suppose that for positive integers k, j ∈ N + , where χ S k and χ H k are given by (3.4) and (3.5) respectively. Then for each k ∈ N + , there exist a positive constant δ and a unique one-parameter curve Γ k (s) = {(u k (s, x), v k (s, x), w k (s, x), χ k (s)) : s ∈ (−δ, δ)} of spatially inhomogeneous solutions (u, v, w, χ) ∈ X × X × X × R to (4.1) that bifurcate from (ū,v,w) at χ = χ S k . Moreover, the solutions are smooth functions of s such that and where (ū k ,v k ,w k ) is given by (4.8) and O(s 2 ) ∈ Z is in the closed complement of N (DF(ū,v,w, χ)) defined by (4.14) Proof. All the necessary conditions except the following have been verified in order to apply the Crandall-Rabinowitz local theory in [9] d dχ (DF(ū,v,w, χ))(ū k ,v k ,w k )| χ=χ S k ∈ R(DF(ū,v,w, χ)), (4.15) where R is the range of the operator. We argue by contradiction and suppose that condition (4.15) fails, then there exists a nontrivial solution (u, v, w) that satisfies Multiplying equations in (4.16) by cos kπx L and integrating them over (0, L) by parts, we obtain that The coefficient matrix is singular because of (4.7), then we reach a contradiction and this completes the proof of condition (4.15). Finally the statements in Theorem 4.2 follow from Theorem 1.7 of [9].

Stability of bifurcating solutions near (ū,v,w, χ S k )
Now we proceed to study the stability of the spatially inhomogeneous solution (u k (s, x), v k (s, x), w k (s, x)) established in Theorem 4.2. Here the stability or instability is that of the bifurcation solution regarded as an equilibrium of system (3.1). To this end, we want to determine the turning direction of the bifurcation branch Γ k (s) around each bifurcation point χ S k . It is easy to see that the operator F is C 4 -smooth if φ is C 5 -smooth, therefore according to Theorem 1.8 in [10], we can write the following expansions where (ϕ i , ψ i , γ i ) ∈ Z in (4.14) and K i are constants for i = 1, 2. o(s 3 ) are taken with respect to the X -topology and o(s 2 ) is a constant. Moreover we have the Taylor expansion First of all, we claim that the bifurcation branch Γ k is of pitch-fork type by showing K 1 = 0. Substituting (4.18) into (4.1) and collecting s 2 terms, we obtain the following system where Multiplying the equations in (4.20) by cos kπx L and integrating them over (0, L) by parts give usū Since (ϕ 1 , ψ 1 , γ 1 ) ∈ Z, we infer from (4.14) that where P k , Q k are given by (4.9) and (4.10). Combining (4.22)-(4.24) gives The determinant of the coefficient matrix M to the system above is and this implies that K 1 = 0 in (4.21). Thus the bifurcation branch Γ k (s) is of pitch-fork, i.e., being one sided. Now we present another main result of this paper which states that the stability of the bifurcating solutions depends on the sign of K 2 .  (s, k), v k (s, k), w k (s, k), χ k (s))} be the bifurcation branch given by (4.12)-(4.13). Denote χ 0 = min k∈N + {χ S k , χ H k } as in (3.3). Then it holds that: (i) If χ 0 = χ S k 0 < min k∈N + χ H k , then Γ k 0 (s) around (ū,v,w, χ S k 0 ) is asymptotically stable when K 2 > 0 and it is unstable when K 2 < 0, while Γ k (s) around (ū,v,w, χ S k ) is always unstable for each k = k 0 ; (ii) If The bifurcation curves Γ k (s) in case (i) are illustrated in Figure 1 schematically. Our results suggest that (ū,v,w, χ S k ) loses its stability to stable steady state bifurcating solution with wave mode number k 0 for which χ S k achieves its minimum over N + . When case (ii) occurs, we surmise that the stability of the homogeneous solution is lost to stable Hopf bifurcating solutions. This is rigorously verified in Section 5. We would like to mention that, K 2 can be evaluated in terms of system parameters and we give the detailed calculations in Appendix. Proof of Theorem 4.3. Our proof follows the approaches in [57,60] based on slight modifications in the arguments for Corollary 1.13 of [10], or Theorem 3.2 of [57], Theorem 5.5, Theorem 5.6 of [8]. We shall only prove case (ii) and case (i) can be treated similarly.
For each k ∈ N + , we linearize (4.1) around (u k (s, x), v k (s, x), w k (s, x), χ k (s)) and obtain the following eigenvalue problem (s, x), v k (s, x), w k (s, x), χ k (s)) is asymptotically stable if and only if the real part of eigenvalue σ(s) is negative.
Sending s → 0, we know from the proof of (4.15) thatσ = 0 is a simple eigenvalue of DF(ū,v,w, χ S k ) = σ(u, v, w) or equivalently which has one-dimensional eigen-space N DF(ū,v,w, χ S k ) = {(P k , Q k , 1) cos kπx L }. Multiplying the system above by cos kπx L and integrating them over (0, L) by parts, we have that σ = 0 is an eigenvalue of (3.2) with χ = χ S k which reads If χ 0 = min k∈N + χ H k < χ S k for all k ∈ N + , or χ 0 = min k∈N + χ S k < χ H k for k = k 0 , we have from the proof of Proposition 3.1 that this matrix always has an eigenvalue σ with positive real part. From the standard eigenvalue perturbation theory in [27], for s being small, there exists an eigenvalue σ(s) to the linearized problem above that has a positive real part and therefore (u k (s, x), v k (s, x), w k (s, x), χ k (s)) is unstable for s ∈ (−δ, δ).

According to Theorem 4.3, the only stable bifurcation branch must be
loses its stability only to nonconstant steady state with wave mode cos k 0 πx L . This gives a wave mode selection mechanism for system (1.1) when χ is around the bifurcation value. In general it is very difficult to determine whether χ 0 is achieved at χ S k or χ H k . According to the discussions after Remark 3.1, if the interval length L is sufficiently small, χ 0 = χ S 1 < min k∈N + χ H k and the only stable bifurcating solution has wave mode cos πx L which is spatially monotone. The wave mode section mechanism given in Theorem 4.3 is verified and illustrated in our numerical studies of (1.1) in Section 6.

Time-periodic positive solutions
In this section, we study the periodic orbits of (3.1) that bifurcate from (ū,v,w) at χ = χ H k . We want to show that under proper assumptions on system parameters, the constant equilibrium (ū,v,w) loses its stability through Hopf bifurcation as χ comes across χ 0 = min k∈N + {χ S k , χ H k }. To apply the bifurcation theory for (3.1) at point χ = χ H k , we need to verify that the real part of eigenvalue crosses the imaginary axis at χ H k . According to the discussions in Section 3, Hopf bifurcation occurs for (3.1) at (ū,v,w) only if χ = χ H k and η 1 (χ, k) > 0, when the stability matrix (3.2) has purely imaginary eigenvalues given by To determine when η 1 (χ H k , k) > 0, we let χ M k be the unique root of η 1 (χ, k) = 0 which is given explicitly in the following form We first give the following fact which will be used in our coming analysis.
Lemma 5.1. Let χ M k be given as above, then for each k ∈ N + ,it holds that either According to Lemma 5.1 and discussions in Section 3, the stability matrix (3.2) has a pair of purely imaginary eigenvalues if and only if χ = χ H k < χ S k , therefore Hopf bifurcation may occur at (ū,v,w, χ H k ) only when χ H k < χ S k . We shall always assume this condition in the coming Hopf bifurcation analysis.

Hopf bifurcation
In this subsection, we prove the existence of Hopf bifurcation of (3.1) assuming that χ H k < χ S k . We recall the notation of Sobolev space X = {u ∈ H 2 (0, L)|u (0) = u (L) = 0} from Section 4. According to the proof of Theorem 2.1, we know that (3.1) is normally parabolic, therefore we can apply the Hopf bifurcation theory from [5] (or Theorem 1.11 from [11], Theorem 6.1 from [36]). Our main result on the existence of nontrivial periodic orbits of (3.1) states as follows.
Proof. We follow the approach in the proof of Theorem 5.2 in [59] or Theorem 3.4 in [36]. According to Proposition 3.1 and Remark 3.1, the stability matrix (3.2) with χ = χ H k has a pair of purely imaginary eigenvalues σ H 2,3 (k) = ± η 1 (χ H k )i; moreover since χ H k = χ H j for ∀j = k, matrix (3.2) has no eigenvalue of the form mτ 0 i for m ∈ N + \{±1}.
Theorem 5.2 implies that system (3.1) admits time-periodic spatial patterns that bifurcate from (ū,v,w, χ H k ) if and only if χ H k < χ S k . Furthermore, it gives the explicit expression of the time-periodic spatial patterns as ϑ k mentioned above with the spatial profile of eigen-function cos kπx L . As we have discussed in Section 3, it is very difficult to determine the necessary condition χ H k < χ S k in terms of system parameters, however if the interval L is sufficiently small, we always have χ H k > χ S k for each k ∈ N + , and therefore this indicates that Hopf bifurcation dose not occur for (3.1) when the interval length is sufficiently small. Indeed, in this case, we already know from the discussions after the proof of Theorem 4.3 that the stability of the homogeneous solution (ū,v,w) is lost through the steady state bifurcation at the first bifurcation branch (ū,v,w, χ S 1 ), which contains stable stationary solutions of (3.1) with eigenfunction cos( πx L ).

Stability of time-periodic bifurcating solutions
We continue to explore the stability of the time-periodic bifurcating solutions on the bifurcation curves ϑ k (s) obtained in Theorem 5.2. The stability here we mean is the formal linearized stability of a periodic solution relative to perturbations from ϑ k (s). Suppose that χ H k 0 = min k∈N + χ H k < χ S k , ∀k ∈ N + , and assume that all the conditions in Theorem 5.2 are satisfied here, then our stability results show that ϑ k (s), s ∈ (−δ, δ) is asymptotically stable only if χ = χ H k 0 . Denote u k (s, x, t) = (u k (s, x, t), v k (s, x, t), w k (s, x, t)) and let (u k (s, x, t), T k (s), χ k (s)) be the periodic solutions on the branch ϑ k (s) obtained in Theorem 5.2. Then we can rewrite (3.1) into the following form Differentiating the system against t, writingu k = du dt , we have du k dt = G u (u k , χ k (s))u k , then we observe that 0 is a Floquet exponent and 1 is a Floquet multiplier for u k . Linearize the periodic solution around the bifurcation branch ϑ k (s) by substituting the perturbed solution u k + we −lt , where w is a sufficiently small T -periodic function and l = l(s) is a continuous function of s, then we have that where G u is the Fréchet derivative with respect to u given by The stability of the bifurcating solutions around χ H k can be determined by computing the eigenvalues of this reduced equation. When s = 0, (5.7) is associated with the eigenvalue problem G 0 (k)w = l(0)w, the spectrum of which is infinitely dimensional. Moreover G 0 corresponds to the stability matrix (3.2).
Suppose that min k∈N + {χ H k , χ S k } = χ H k 0 for some k 0 ∈ N + . We first show that ϑ k (s) around χ H k is unstable for any k = k 0 . Denote the eigenvalues of A k (χ H k ) by σ H 1 (χ H k , k), σ H 2 (χ H k , k) and σ H 3 (χ H k , k). According to the Proposition 3.1, there exists at least one eigenvalue with positive real part if χ > χ 0 . Therefore for any positive integer k = k 0 , we have that G 0 (k) must have an eigenvalue with positive real part hence l(0) < 0 if k = k 0 . By the standard perturbation theory for an eigenvalue of finite multiplicity [17,27], l(s) < 0 for s being small if k = k 0 , therefore all the bifurcation branches ϑ k (s) around (ū,v,w) are unstable if k = k 0 . The result indicates that if a periodic bifurcation solution is stable, it must be on the ϑ k 0 (s) branch where χ H k 0 < min k∈N + χ S k , i.e., it is on the left-most branch, while the branches on its right hand side are always unstable. Now we proceed to discuss the stability of branch ϑ k 0 (s) around (ū,v,w, χ H k 0 ). According to Lemma 2.10 in [11] (or [24,25]), the eigenvalue l(k) is a continuous real function of s near the origin. For χ being around χ H k 0 , the eigenvalue of (5.9) are σ 1 (χ, k) and σ 2,3 (χ, k) = λ(χ, k) ± iτ (χ, k). According to Theorem 2.13 in [11], l(s) and sχ k 0 (s) have the same zeros in small neighbourhood of s = 0 where l(k) and −λ(χ H k 0 )sχ k 0 s(s) have the same sign −λ(χ H k 0 )sχ k 0 s(s) = 0, l(s) = 0, and |l(s) + λ((χ H k 0 )sχ k 0 s(s)| ≤ |sχ k 0 (s)|o(1), as s → 0.
According to Theorem 8.2.3 in [17], if l(s) > 0, the periodic bifurcation solutions are orbitally asymptotically stable and, if l(s) < 0, the periodic bifurcation solutions are orbitally unstable. We have proved that λ(χ H k 0 ) < 0 and l(s) and sχ k 0 (s) has the same sign. Therefore, assuming that χ k 0 (0) = 0, if the branching solutions appear supercritical, they are stable and if they appear subcritical, they are unstable. Therefore, one need to compute χ k 0 and/or χ k 0 similarly as in Section 4. The calculations are straightforward but complicated and we skip them here for simplicity.

Numerical simulations
This section is devoted to the numerical studies of system (3.1). We are motivated to investigate the effects of prey-taxis on the formation of nontrivial patterns to this system. In particular, we show that the one-dimensional system admits both stationary and timeperiodic solutions emerging from bifurcations. Moreover, we shall see that, when the prey-taxis rate χ is taken to be greatly larger than the critical bifurcation value χ 0 , (3.1) can develop various interesting patterns with striking structures such as spikes, propagation, coarsening, etc.
(ū,v,w)+(0.01, 0.01, 0.01) cos πx, which are small perturbations from the homogeneous equilibrium. We see that the initial data have spatial profiles in the form of cos πx, but the spatial-temporal patterns develop according to the stable wave mode cos 6πx 7 . Numerical simulations in Figure 3 are devoted to verifying that (ū,v,w) loses its stability to steady state bifurcations when χ 0 is achieved at χ S k 0 for this set of system parameters. These simulations support the results on the stability of the bifurcating solutions in Theorem 4.3. . System parameters in all the graphes are taken to be the same as in Table 1 except that χ = 8, which is larger than χ S k 0 ≈ 6.05 given in Table 1. Initial data are (u 0 , v 0 , w 0 ) = (ū,v,w) + (0.01, 0.01, 0.01) cos πx, while the stable pattern has wave mode cos 6πx 7 . These graphes support our stability analysis of the bifurcating solutions. 6.04 6.03 6.04 6.04 6.04 6.04 6.04 6.04 Table 2: Stable wave mode numbers and the corresponding bifurcation values χ 0 for different interval lengthes. System parameters are chosen to be the same as those in Table 1. We see that the threshold value χ 0 is always achieved at the steady state bifurcation point χ S k 0 . This table also indicates that larger intervals support higher wave modes.

Time periodic patterns
Our next set of numerical results are provided to demonstrate that (3.1) admits timeperiodic patterns through Hopf bifurcations. To this end, we set system parameters to be d 1 = d 3 = 1, d 2 = 0.1, α 1 = 0.02, α 2 = 0.04, α 3 = 8 and β 1 = 0.05, β 2 = β 31 = β 32 = 0.5, while the sensitivity function is chosen φ(w) = 0.05w(0.2 − w). We shall show that the equilibrium (ū,v,w) = (2.13, 6.65, 0.45) loses its stability to time-periodic orbits. Table 3 lists the values of χ S k and χ H k when the interval length is L = 7. We see that the threshold value χ 0 is achieved at χ H 3 ≈ 92.57. In Figure 4, we plot the numerical solutions of (3.1) subject to initial data (u 0 , v 0 , w 0 ) = (ū,v,w) + (0.01, 0.01, 0.01) cos πx, small perturbations from the homogeneous equilibrium. The initial data have spatial profiles in the form of cos πx, but the spatial-temporal patterns develop according to the stable Figure 3: Formation of stationary patterns of (3.1) over intervals with lengthes L = 9, 11, 13 and 15. System parameters here are taken to be the same as those in Table 1 except that χ = 8, which is slightly larger than χ S k 0 ≈ 6.04 given in Table 2. Initial data are (u 0 , v 0 , w 0 ) = (ū,v,w) + (0.01, 0.01, 0.01) cos πx. These graphes support our stability analysis of the bifurcating solutions and indicate that large intervals support more aggregates than small intervals. . We see that min k∈N + {χ S k , χ H k } = χ H 3 . Therefore (ū,v,w) loses its stability to the time-periodic solutions with wave mode cos 3πx 7 . This is numerically verified in Figure 4.  Table 4: Stable wave mode numbers and the corresponding bifurcation values χ 0 for different interval lengthes, where system parameters are chosen to be the same as those in Table 3. We see that the threshold value χ 0 is always achieved at the Hopf bifurcation point χ H k 0 .
Numerical simulations in Figure 5 are devoted to verifying that the (ū,v,w) loses its stability to Hopf bifurcations when χ 0 is achieved at χ H k 0 , where system parameters are . System parameters in all the graphes are taken to be the same as in Table 3 except that χ = 120, which is slightly larger than χ H k 0 ≈ 92.57 given in Table 3. Initial data are (u 0 , v 0 , w 0 ) = (ū,v,w) + (0.01, 0.01, 0.01) cos πx, however the stable oscillating patterns have spatial profile cos 3πx 7 , which emerge periodically. These plots support our stability analysis in Section 5.
taken to be the same as those in Table 3. In particular, we select the interval lengths to be L = 9, 11, 13 and 15 respectively. These simulations support our results on the stability of the Hopf bifurcating solutions obtained in Theorem 5.2. Figure 5: Formation of time-periodic spatial patterns of (3.1) over intervals with lengthes L = 9, 11, 13 and 15 respectively. System parameters in all the graphes are taken to be the same as in Table 3 except that χ = 120, which is slightly larger than χ H k 0 given in Table 4. Initial data are (u 0 , v 0 , w 0 ) = (ū,v,w) + (0.01, 0.01, 0.01) cos πx. These graphes support our stability analysis of the bifurcating solutions.

Other interesting patterns
In Figure 6, we plot the formation of stable boundary spikes of (3.1) through traveling wave over Ω = (0, 7). System parameters are taken to be d 1 = 5, d 2 = 0.5, d 3 = 1, α 1 = α 2 = 0.05, α 3 = 4, β 1 = 0.3, β 2 = 0.5, β 31 = 1, β 32 = 1 3 , ξ = 0.1 and φ(w) = w(0.1 − w). Prey-taxis rate χ = 3000 is greatly larger than the critical bifurcation value χ 0 = 956.79. We observe that the boundary spike is developed through traveling wave solution. However, rigorous analysis of qualitative properties of the propagating solutions is out of the scope of our paper. Finally, we present numerical simulations in Figure 7 to show that when the prey-taxis is much larger than χ 0 , (3.1) admits some other interesting and striking dynamics such as merging and emerging of spikes, irregular spatial-temporal oscillations etc. For example, Subplot (i) of Figure 7 shows that there occurs a coarsening process in (3.1) in which interior spikes of u(x, t) shift to the boundary or the center to merge into another stable spike. We also observe the spontaneous emergence of stable interior spikes at time t ≈ 100. All the parameters and the initial data in (3.1) are chosen to be the same as in Figure 2, except that χ = 124, which is much far away from χ 0 ≈ 6.05. In subplot (ii), when the system parameters and the initial data are chosen to be the same as in Figure 4, except that χ = 160, (3.1) initially develops time-periodic spatial patterns which are metastable. Then the oscillating patterns develop into stable stationary spikes. Time-periodic patterns and spontaneous initiation of interior multiple spikes are observed in subplots (iii) and (iv).

Conclusions and discussions
Our paper investigates population dynamics of a two-predator and one-prey model with prey-taxis, given by a 3×3 reaction-advection-diffusion system. It is proved that the system admits positive classical solution which is global and uniformly bounded in time over 1D or 2D bounded domains. The same results are obtained for its parabolic-parabolicelliptic counterpart for domains of arbitrary space dimension.
Stability of the unique positive equilibrium is studied when the domain is a finite interval. It shows that both prey-taxis χ and sensitivity function φ determines the linearized stability of this equilibrium. It is known (see [33] e.g.) that, in contrast to chemotaxis [18,19] or advection for competition system [57], prey-taxis stabilizes constant equilibrium for one-predator and one-prey system. However, our result reveals that this is true only when there is no group defense in the preys, i.e., a huge amount of preys can aggregate and keep their predators away from the habitat. If the predators retreat from the habitat, which can be modeled by choosing φ(w) < 0, prey-taxis destabilizes the constant equilibrium, which becomes unstable as χ surpasses χ 0 = min k∈N {χ S k , χ H k } given by (3.3). Therefore group defense is an important mechanism in the formation of nontrivial patterns in (1.1).
We have obtained both stationary and time-periodic spatial patterns to the system over 1D bounded interval through steady state bifurcation at χ = χ S k and Hopf bifurcation at χ = χ H k respectively. Stabilities of these bifurcating solutions are also investigated rigorously. It is proved that steady state bifurcation occurs at (ū,v,w, χ S k ) for each k ∈ N + , however, only the k 0 -branch that turns to the right is stable if χ 0 = χ S k 0 < min k∈N χ H k . In other words, if the steady state bifurcation curve is stable, it must be on the left most branch on the χ-axis. On the other hand, Hopf bifurcation (ū,v,w, χ H k ) only if χ H k < χ S k , while only the left most branch can be stable. Moreover, our analysis indicates that small intervals only supports steady state bifurcations while large intervals may lead to Hopf bifurcation when the system parameters are chosen properly. Extensive numerical simulations are performed to illustrate and support our theoretical findings. Apparently, the formation of these nontrivial patterns is due to the effect of large prey-taxis and prey group defense effect.
Global existence and bounded are obtained for (1.1) over 2D and it is interesting to ask the same question for the system over higher dimensions. Logistic decays in the kinetics help to prevent finite or infinite time blow-ups, however, whether or not they are sufficient over higher dimensions, in particular when the prey-taxis rate is large, is unknown in the literature.
Our bifurcation analysis is based on the local versions in [5,9] etc. From the viewpoint of mathematical analysis, it is interesting to investigate the behavior or shape of these local branches, in particular in the study of positive steady states when large prey-taxis may lead to striking structures such as spikes and layers, etc. For example, according to the global theory of Rabinowitz [47] and its developed version in [50], global continuum of Γ k (s) either intersects with the χ-axis at another bifurcating point, or extends to infinity, or intersects with a singular point. Populations growth terms in (3.1) inhibits the application of topology argument developed in [8,62] etc.
When the prey-taxis rate is around bifurcation values, our results provide almost a complete understanding of the spatial-temporal dynamics of (1.1) over 1D. Further research is needed on its pattern formations when χ is away from χ 0 and in particular when it is sufficiently large. For example, rigorous analysis of the profile of the spikes obtained in numerical simulations can be an interesting problem to probe in the future. There are also some interesting problems such as the investigation of chaotic dynamics in (1.1) or bifurcation analysis of (1.1) over higher dimensions. It is also meaningful to ask about the biologically realistic traveling wave solutions to (3.1), compared to those for the system without prey-taxis obtained in [35].