Global existence and steady states of a two competing species Keller-Segel chemotaxis model

We study an one{dimensional quasilinear system proposed by J. Tello and M. Winkler [19] which models the population dynamics of two competing species attracted by the same chemical. The kinetics terms of the interacting species are chosen to be the Lotka{Volterra type. We prove the existence of global bounded and classical solutions for all chemoattraction rates. Under homogeneous Neumann boundary conditions, we establish the existence of nonconstant steady states by local bifurcation theory. The stability of the bifurcating solutions is also obtained when the diffusivity of both species is large. Finally, we perform extensive numerical studies to demonstrate the formation of stable positive steady states with various interesting spatial structures.


Introduction
This paper is devoted to study the global solutions and nonconstant positive steady states of the following reaction-advection-diffusion system, x ∈ (0, L), t > 0 w t = w xx − λw + u + v, x ∈ (0, L), t > 0 u x (x) = v x (x) = w x (x) = 0, x = 0, L, t > 0, u(x, 0) = u 0 (x), v(x, 0) = v 0 (x), w(x, 0) = w 0 (x), x ∈ (0, L), (1.1) where u, v and w are functions of space x and time t. d i , µ i , a i and λ are positive constants for i = 1, 2; χ and ξ are constants and τ is nonnegative. (1.1) is the one-dimensional version of the following model x ∈ Ω, t > 0, τ w t = ∆w − λw + u + v, x ∈ Ω, t > 0, ∂u ∂n = ∂u ∂n = ∂u ∂n = 0, x ∈ ∂Ω, t > 0, u(x, 0) = u 0 (x), v(x, 0) = v 0 (x), w(x, 0) = w 0 (x), x ∈ Ω, (1.2) which was considered by J. Tello and M. Winkler [19] to study the spatial-temporal behaviors of two competing species attracted by the same chemical. It is assumed in [19] that Ω ⊂ R N , N ≥ 1, is a bounded domain with smooth boundary ∂Ω. u(x, t) and v(x, t) represent the population densities of two competing species at space-time location (x, t) ∈ Ω × R + , while w(x, t) denotes the concentration of the attracting chemical. It is assumed that both species u and v direct their movements chemotactically along the gradient of the chemical concentration over the habitat. This is modeled by taking both χ > 0 and ξ > 0. Biologically, χ and ξ measure the strength of chemical attraction to species u and v respectively. The kinetics of the species are assumed to be of the classical Lotka-Volterra type and a 1 , a 2 interpret the levels of inter-specific competitions. Moreover, the chemical is produced by both species at the same rate with no saturation effect and is consumed by certain enzyme at the rate of λ meanwhile.
It is easy to see that system (1.2) has a positive constant steady state provided with 0 ≤ a 1 < 1, 0 ≤ a 2 < 1.
(1.4) J. Tello and M. Winkler investigated the global existence of parabolic-parabolic-elliptic system (1.2) and the global asymptotic stability of the constant equilibrium (ū,v,w). Their results can be summarized as follows.
Theorem 1.1 suggests that, for any given µ 1 , µ 2 > 0, small chemoattraction rates χ, ξ and competition rates a 1 , a 2 stabilize the positive equilibrium (ū,v,w), which is globally asymptotically stable. Therefore, the stationary system of (1.2) has no positive solutions other than (ū,v,w) provided with condition (1.5). From the view point of mathematical modeling, it is interesting to obtain the existence of nonconstant steady states of the two species chemotaxis model (1.2). On the other hand, it is also important to study the evolutions of species population distributions through inter-specific competition and chemotaxis processes. It is the motivation of this paper to study the range of the parameters χ and ξ in the one-dimensional model (1.1) admit nonconstant positive steady states. In particular, we are concerned with the existence of nonconstant positive steady states of the stationary system of (1.1).
Some other two-species chemotaxis systems closely related to (1.2) are proposed and studied by various authors. For µ 1 = µ 2 = τ = 0 and λw being replaced by a positive constant, Espejo et al. [10] investigated simultaneous finite-time blow-ups of (1.2) when Ω is a circle in R 2 . For µ 1 = µ 2 = λ = τ = 0, P. Biler et al. [4] obtained the blow-up properties of (1.2) with Ω = R N , N ≥ 2. Similar blow-up mechanisms, in particular related to the initial data size, have been studied in [5] and [6] for Ω = R 2 . X. Wang and Y. Wu [20] studied the qualitative behaviors of two competing species chemotaxis models with cellular growth kinetics different from those in (1.1) In particular, they investigated the large time behaviors of global solutions and the existence of nonconstant positive steady states. Moreover, the effect of cell motility and chemotaxis on population growth have also been studied. Horstmann [12] proposed some multi-species chemotaxis models with attraction and repulsion between interacting species.
A linearized stability analysis of the constant steady state is performed and the existence of Lyapunov functional is also obtained for some of the proposed new models. P. Liu et al. [15] study the pattern formation of an attraction-repulsion Keller-Segel system with two chemicals and one species.
In this paper, we study the reaction-advection-diffusion system (1.1) and its stationary system. We are concerned with the global existence of system (1.1) and the formation of its spatially inhomogeneous steady states when the smallness condition (1.5) on χ and ξ is relaxed. Our results state that, (1.1) admits a unique pair of classical solution (u, v, w) which is uniformly bounded in L ∞ for all t ∈ (0, ∞) for all χ, ξ ∈ R-see Theorem 2.5. In section 3, we establish the existence and stability of nonconstant positive steady states by rigorous bifurcation analysissee Theorem 3.1 and Theorem 3.2. In loose terms, large chemotaxis rate χ (the same for ξ) tends to destabilize homogeneous positive steady state and nonconstant positive steady states emerge to form stable patterns. In section 4, we present some numerical simulations to illustrate the formation of self-organized spatial patterns of (1.1) that have complicated and interesting structures. Finally, we include some concluding remarks and propose interesting future questions in section 5.

Existence of global-in-time solutions
In this section, we study the global existence of uniformly bounded and classical solutions (u, v, w) of (1.1). We shall first apply the well-known results of Amann [2,3] to obtain the local existence, then we obtain the global existence once we can find the L ∞ -bounds of (u, v, w). We shall study (1.1) with τ > 0 throughout the rest of this section. Without loss our generality, we choose τ = 1 from now on.

Local existence
We present the existence and uniqueness of the positive local classical solutions of (1.1) as follows.
Theorem 2.1. Assume that the constants τ > 0 and a i , d i , µ i > 0 for i = 1, 2. Suppose that the initial data u 0 , v 0 and w 0 ∈ H 1 (0, L) satisfy u 0 , v 0 and w 0 ≥ 0, ≡ 0 on [0, L]. Then for any χ, ξ ∈ R, the following statements hold true: (i) There exists a unique solution (u(x, t), v(x, t), w(x, t)) of (1.1) which is nonnegative on System (2.1) is normally parabolic since all the eigenvalues of A 0 are positive. Then (i) follows from Theorem 7.3 and Theorem 9.3 of [2]. Moreover, (ii) follows from Theorem 5.2 in [3] since (2.1) is a triangular system. The nonnegativity of (u, v, w) follows from parabolic Maximum Principles-see [14].

A-prior estimates and global existence
We now proceed to establish the L ∞ -bounds of (u, v, w) under the same conditions in Theorem 2.1. This is a consequence of several lemmas.
Lemma 2.2. Let µ 1 , µ 2 > 0 and (u, v, w) be the unique positive solution of (1.1), then there exists a constant C > 0 dependent on (u 0 , v 0 , w 0 ) L 1 and L such that Proof. We integrate the u-equation of (1.1) over (0, L) and have that then it follows from the Gronwall's lemma that Moreover, integrating the w-equation over (0, L), we easily see that w(·, t) L 1 (0,L) is uniformly bounded for all t ∈ (0, ∞). This completes the proof of lemma2.2.
Lemma 2.2 gives us the L 1 -bounds of u, v and w. To obtain their L ∞ -bounds, we shall see that it is necessary to obtain the boundedness of w x (·, t) L ∞ . For this purpose, we first convert the w-equation into the following abstract form where ∆ = d 2 dx 2 . To estimate w(x, t) in (2.3), we want to apply the well-known smoothing properties of the operator −∆ + 1 and the estimates between the linear analytic semigroups generated by e t∆ ,t ≥ 0, for which [11], [13] and [21] are good references. For example, we have from Lemma 1.3 with N = 1 in [21], for all 1 ≤ p ≤ q ≤ ∞, that there exists a positive constant C dependent on d 3 , µ 1 , µ 2 , and w 0 W 1,q (0,L) such that for all T ∈ (0, ∞) and any t ∈ (0, T ), where ν is the first Neumann eigenvalue of −∆. We have the following a-prior estimate on w(·, t) W 1,p (0,L) . Proof. Choosing p = 1 in (2.4), we have that + v(·, s) L 1 (0,L) + w(·, s) L 1 (0,L) ds , where C depends on (u 0 , v 0 , w 0 ) W 1,q (0,L) . On the other hand, we can easily have that there exists a constant C 0 > 0 such that for any q ∈ (1, ∞), therefore, we conclude from (2.6) that Finally, (2.5) is an immediate consequence of (2.2) and (2.7).
Proof. We only show that u(·, t) L p ≤ C(p) and v(·, t) L p ≤ C(p) can be proved by the same arguments. For p > 2, we multiply the first equation of (1.1) by u p−1 and integrate it over (0, L) by parts, then it follows from simple calculations that where C 1 is a positive constant that depends on p. On the other hand, it follows that the Hölder's inequality and the Young's inequality that where we have applied the facts that u and w x L 2(p+1) ≤ C due to (2.5).
In light of (2.10), we obtain from 2.9 that Denoting y p (t) = L 0 u p (x, t)dx, one readily sees that (2.11) implies p , y p (0) = u 0 p L p . Solving this differential inequality, we conclude that y p (t) ≤ C(p) for all t ∈ (0, ∞). Similarly we can show that L 0 v p (x, t)dx ≤ C(p). This finishes the proof of Lemma 2.4. By taking p sufficiently large but fixed in (2.4), we can easily obtain the uniform boundedness of w x (·, t) L ∞ thanks to Lemma 2.4. Corollary 1. Under the conditions in Lemma 2.2, we have that w(·, t) W 1,∞ ≤ C, ∀t ∈ (0, T max ), (2.12) We are now ready to present the following results on the global existence of uniformly bounded and classical positive solutions to (1.1).
Theorem 2.5. Let a i , µ i , i = 1, 2, and λ be positive constants. Then given positive initial data (u 0 , v 0 , w 0 ) ∈ H 1 (0, L) × H 1 (0, L) × H 1 (0, L) and any constants χ, ξ ∈ R, (1.1) has a unique bounded positive solution (u( Proof. According to (ii) of Theorem 2.1 and Corollary 1, we only need to show that (u(·, t), v(·, t)) L ∞ is uniformly bounded for all t ∈ (0, T max ), then we must have that T max = ∞ and the existence part of Theorem 2.5 follows. Moreover, one can apply parabolic boundary L p estimates and Schauder estimates to show that u t , v t , w t and all spatial partial derivatives of u, v and w up to order two are bounded on [0, L] × [0, ∞), therefore (u, v, w) have the regularities stated in Theorem 2.5. Without loss of our generality, we assume in (2.12) that w x (., t) L ∞ ≤ 1. Through the same calculations that lead to the proof of Lemma 2.4 and using the fact u ≥ 0, v ≥ 0, we obtain where we have used Young's inequality in third line of (2.13). We shall need the following estimate from P. 63 in [14] and Corollary 1 in [7] with d = 1 due to Gagliardo-Ladyzhenskaya-Nirenberg inequality for any u ∈ H 1 (0, L) and any > 0, there exists some C 0 > 0 which only depends on L such that (2.14) Thanks to (2.15) and (2.13) gives rise to On the other hand, we can choose p 0 large such that for all p ≥ p 0 , = (2.16) . For each T ∈ (0, ∞), we solve the differential inequality (2.16) to have that for all p ≥ p 0 (2.17) We now use the Moser-Alikakos iteration [1] to establish L ∞ -estimate of u through (2.17). To this end, we let Taking p = 2 i , i = 1, 2, ... and choosing p 0 = 2 i0 for i 0 large but fixed, we have where C is a constant that only depends on L and M (2 i0 ) is bounded for all t ∈ (0, ∞) in light of (2.8). Sending i → ∞ in (2.18), we finally conclude from Lemma 2.4 and (2.18) that By the same calculations we can show that v(·, t) L ∞ is uniformly bounded for all t ∈ (0, ∞). This completes the proof of Theorem 2.5.

Existence and stability of nonconstant positive steady states
In this section, we consider the stationary system of (1.1) in the following form where denotes the derivative taken with respect to x. From the view point of mathematical modeling, it is interesting to study whether or not the chemotaxis can be responsible for cell aggregation through interspecific competitions. Cellar aggregation can be modeled by nonconstant positive solution of 3.1. For this purpose, we study the emergence of nonconstant steady states in (3.1). In particular, we are concerned with the effect of the chemotactic sensitivities on the pattern formations. Unlike random movements, a chemotactic process is anti-diffusion and it has the effect of destabilizing the spatially homogeneous solutions. Then spatially inhomogeneous solutions may arise through bifurcation as the homogeneous solution becomes unstable.

Linearized stability analysis of the homogeneous solution (ū,v,w)
To study the mechanism through which spatially inhomogeneous solutions of (3.1) emerge, we carry out implement the standard linearized stability of (ū,v,w) viewed as an equilibrium of (1.1). To this end, we take (u, v, w) = (ū,v,w) + (U, V, W ), where U , V and W are small perturbations from (ū,v,w), then it is easy to see that (U, V, W ) satisfies the following system According to the standard linearized stability analysis-see [18] for example, we have that the stability of (ū,v,w) can be determined by the following matrix, where Λ = ( kπ L ) 2 > 0, k = 1, 2, ..., are the eigenvalues of − d 2 dx 2 on (0, L) under the Neumann boundary conditions. We have the following result on the linearized instability of (ū,v,w).

3)
and locally asymptotically stable if χ < χ 0 = min k∈N + {χ k ,χ k }, wherẽ Proof. By the principle of the linearized stability-see Theorem 5.2 in [18], (ū,v,w) is asymptotically stable with respect to (1.1) if and only if the real parts of all eigenvalues of the matrix (3.2) is negative. To this end, we first see that its characteristic polynomial reads It is easy to see that α 2 > 0 for all positive χ and ξ. According to the Routh-Hurwitz conditions, or Corollary 2.2 in [15], we have that the real parts of all σ are negative if and only if while the real parts of some σ are positive if one of the conditions above fails. Moreover, we can show by simple calculations that α 1 > 0 is a necessary condition whenever α 0 > 0 and α 1 α 2 − α 0 > 0. Therefore, in order to obtain the unstable condition in terms of χ, we just need one of the following conditions to be satisfied, If α 0 < 0, we have from simple calculations that Therefore (ū,v,w) is unstable if χ > χ 0 according to (S1) of Corollary 2.2 in [15]. Similarly, we can show that (ū,v,w) is locally asymptotically stable if χ < χ 0 .
We see that the homogeneous steady state (ū,v,w) loses its stability at χ = χ 0 . The critical value χ 0 = χ 0 (ξ) decreases as ξ increases; moreover χ 0 < 0 if ξ is sufficiently large. Therefore only one of χ and ξ is needed to be large to destabilize (ū,v,w). If ξ < 0, i.e., species v is repulsive to the chemical gradient, then χ needs to be large to stabilize (ū,v,w). Above all, the local stability analysis suggests that chemo-attraction destabilizes constant steady states and the chemo-repulsion stabilizes constant steady states. The constant solution is always stable both when χ < 0 and ξ < 0. Therefore, we surmise that (ū,v,w) is also a global attractor of (1.1) if χ < 0 and ξ < 0, though the positiveness of χ and ξ is required in [19]. We want to remark that Proposition 1 also holds for multi-dimensional domain Ω ⊆ R N with N ≥ 2 with ( kπ L ) 2 being replaced by the k-eigenvalue of −∆ under the Neumann boundary condition.

Bifurcation of steady states
We now proceed to perform rigourous steady state bifurcation analysis to seek non-constant positive solutions to the stationary system (3.1). As we have shown above, the chemoattraction rate χ and ξ has the effect of destabilizing (ū,v,w), which becomes unstable if χ crosses χ 0 . The linearized instability of (ū,v,w) in (1.1) is insufficient to guarantee the existence of spatially inhomogeneous steady states. Then we are concerned if a stable spatially inhomogeneous steady state of (3.1) may emerge through bifurcations as χ increases. We shall see that the constant solution (ū,v,w) loses its stability to spatially inhomogeneous solutions through steady state bifurcation. Clearly, the emergence of spatially inhomogeneous solutions is due to the effect of large chemoattraction rate χ. We refer this as advection-induced patterns in the sense of Turing's instability. We shall apply the bifurcation theory of Crandall and Rabinowitz [8] to seek non-constant positive solutions of (3.1). To this end, we introduce the Hilbert space and convert (3.1) into the following abstract form It is easy to see that F(ū,v,w, χ) = 0 for any χ ∈ R and F : X × X × X × R → Y × Y × Y is analytic for Y = L 2 (0, L); moreover, we can get from straight calculations, for any fixed (û,v,ŵ) ∈ X × X × X , that the Fréchet derivative of F is given by To see this, we denote u = (u, v, w) T and write (3.6) into It is obvious that operator (3.6) is strictly elliptic since all eigenvalues of A 0 are positive; moreover it satisfies the Agmon's condition according to Remark 2.5 of case 2 in Shi and Wang [17] with N = 1 . Therefore, D (u,v,w) F(û,v,ŵ, χ) is a Fredholm operator with zero index thanks to Theorem 3.3 and Remark 3.4 of [17].
For local bifurcations to occur at (ū,v,w, χ), we need the Implicit Function Theorem to fail on F, or equivalently the following non-triviality condition is satisfied where N denotes the null space. Taking (û,v,ŵ) = (ū,v,w) in (3.6), we have that and its null space consists of solutions to the following system Write (u, v, w) into their eigen-expansions and substitute them into the system above. We have that (t k , s k , r k ) satisfies First of all, k = 0 can be easily ruled out in (3.8) if (t k , s k , r k ) is nonzero. For k ∈ N + , the coefficient matrix of (3.8) is singular and it follows from simple calculations that It is easy to see that χ k > 0 if and only if Moreover, the null space D (u,v,w) F(ū,v,w, χ) has dimension 1 and it is spanned by We want to point out that the minimum bifurcation value χ k in (3.9) is greater or equal to χ 0 in Proposition 1, where (ū,v,w) loses its stability. If χ 0 = min k∈N + χ k , then we shall see that (ū,v,w) loses its stability to spatially inhomogeneous solutions emerge through steady state bifurcations. If χ 0 = min k∈N +χ k < min k∈N + χ k , then a Hopf bifurcation may occur at (ū,v,w) to create spatial time-periodic patterns, which is demonstrated numerically in section 4. However, the rigourous analysis is out of the scope of this paper and we refer the reader to the discussions in Section 2 of [15].
Having the potential bifurcation value χ k in (3.9), we now verify that local bifurcation does occur at (ū,v,w, χ k ) for each k ∈ N + in the following theorem, which establishes the existence of nonconstant positive solutions to (3.1).
Proof. To apply the local theory of Crandall and Rabinowitz in [8], we now only need to verify the following transversality condition where R(·) denotes the range of the operator. If (3.17) fails, there must exist a nontrivial solution (u, v, w) to the following problem Multiplying the both equations in (3.18) by cos kπx L and integrating them over (0, L) by parts, we obtain We readily see that the coefficient matrix of this system is singular because of (3.9), therefore we reach a contradiction in (3.19) and this verifies (3.17). On the other hand, since χ k = χ j for all positive integers k = j, (3.13) follows from straightforward but lengthy calculations, which is left to the reader. This finishes the proof of Theorem 3.1.

Stability of bifurcating solutions
Now we proceed to study the stabilities of the bifurcating solutions established in Theorem 3.1.
Here the stability or instability is that of u k (s, x), v k (s, x), w k (s, x) viewed as equilibrium of system (1.1). To this end, we first determine the turning direction of each bifurcation branch Γ k (s) around (ū,v,w, χ k ), k ∈ N + . It is easy to see that F is C 4 -smooth, therefore according to Theorem 1.18 from [8], the bifurcating solution (u k (s, k), v k (s, k), w k (s, k), χ k (s)) is C 3 -smooth of s and we can write the following expansions where for i = 1, 2, (ϕ i , ψ i , γ i ) ∈ Z defined in (3.16) and K i is a constant to be determined. o(s 3 ) -terms in (3.20) are taken with respect to the X -topology.
Before we perform the stability analysis, we claim that K 1 = 0 and the bifurcation is of pitch-fork type. Moreover the sign of K 2 determines the stability of the bifurcating solution if K 2 = 0. To be precise, we have the following proposition.
The bifurcation curves Γ k (s) are schematically illustrated in Figure 1.
Γ k (s) turns to the right and the bifurcating solutions is stable around (ū,v,w, χ k ) if K 2 > 0; Γ k (s) turns to the right and the bifurcating solutions is stable around (ū,v,w, χ k ) if K 2 > 0. Proof. We only show that K 1 = 0 for the completeness and our coming analysis, while the rest part of this proposition can be proved by slightly modifying the arguments in Corollary 1.13 of [9], or Theorem 5.5, Theorem 5.6 of [7].
Substituting (3.20) into (3.1) and equating the s 2 -terms, we collect the following system Multiplying the first three equations of (3.21) by cos kπx L and then integrating them over (0, L) by parts, we have from simple calculations that On the other hand, since (ϕ 1 , ψ 1 , γ 1 ) ∈ Z, we have from (3.16) that where P k , Q k are defined in (3.11) and (3.12). From (3.23)-(3.25), we arrive at the following system We claim that the coefficient matrix M is nonsingular. Indeed, we have from straightforward calculations that its determinant equals where the last identity follows by direct substitutions of (3.11) and (3.12). Hence we can quickly have that We now proceed to evaluate K 2 which determines the stability of (u k (s, x), v k (s, x), w k (s, x)) for s ∈ (−δ, δ). Equating the s 3 -terms in (3.1), we collect the following system where we have substituted χ k in (3.9) above and used the fact that P k + Q k = ( kπ L ) 2 + λ. On the other hand, we can test the second and the third equation of (3.28) by cos kπx L over (0, L) to obtain Moreover, it follows from (ϕ 2 , ψ 2 , γ 2 ) ∈ Z defined in (3.16) that, Putting (3.30)-(3.32) together, we arrive at the following system, where we have used the following notation Solving the algebraic system (3.33) gives us where we have denoted

38)
and Integrate (3.21) over (0, L) by parts, then thanks to the non-flux boundary conditions and K 1 = 0, we obtain that

44)
and Similarly, we test (3.21) by cos 2kπx L and obtain that Solving system (3.47), we have that where we have applied the following notations in (3.48) and (3.49), and We recall from (3.29) that K 2 consists of integrals (3.35)-(3.36), (3.42) and (3.48)-(3.49). As we can see from the formulas above, the calculations to evaluate K 2 are extremely complicated and lengthy, therefore, for the simplicity of calculations, we suppose that both d 1 and d 2 are large through out the rest part of our calculations. Biologically, this interprets the scenario that both species u and v diffuse much faster than the chemical. Moreover, we shall see in our numerical simulations that, besides the small amplitudes solutions established through bifurcations, large diffusion rates d 1 and d 2 also support formation of nontrivial patterns that have interesting structures, when other system parameters are properly chosen.
First of all, we have the following asymptotic expansions from (3.9) and (3.11)-(3.12) as d 1 → ∞ and d 2 = O(d 1 ) These asymptotic formulas be extensively used through the rest of this section. Substituting (3.54)-(3.55) into (3.50) and (3.51), we have that and Together with (3.56)-(3.57), we have from straightforward calculations that We are ready to present the stability of the nonconstant steady states (u k (s, x), v k (s, x), w k (s, x)) constructed in Theorem 3.1.
Theorem 3.2. For each k ∈ N + , the bifurcation branch Γ k (s) turns to the right and the nonconstant positive steady state (u k (s, x), v k (s, x), w k (s, x)) obtained in Theorem 3.1 is asymptotically stable if K 2 > 0; Γ k (s) turns to the left and (u k (s, x), v k (s, x), w k (s, x)) is unstable if K 2 < 0; moreover, for each k ∈ N + , there exists a largeD such that for all d 1 and d 2 larger thanD, Proof. Following the same calculations that lead to (3.58), we obtain from (3.50), (3.52) and (3.53) that and Moreover, it follows from (3.37)-(3.40) that and Finally, we conclude from (3.43)-(3.45) that and Putting (3.58)-(3.65) all together, we finally obtain that 48(1−a1a2) 2 d 1 > 0 in light of (1.4). Therefore, for d 1 and d 2 being sufficiently large, we have that Theorem 3.2 follows immediately from straightforward calculations in (3.67) thanks to Proposition 2.
Theorem 3.2 asserts that, when both species diffuse fast, the small amplitude bifurcating solution (u k (s, x), v k (s, x), w k (s, x)) is unstable if the chemical decay rate is strong, or the domain size is large. In this scenario, it is natural to expect the formation of stable steady states of (1.1) with large amplitude, such as boundary layer, interior spikes, etc., or the formation of time periodic spatial patterns. However, rigourous analysis for these solutions requires nontrivial mathematical tools which is out of the scope of this paper.

Numerical simulations
In this section, we perform extensive numerical studies to show the spatial-temporal behaviors of the IBVP (1.1). We choose different sets of parameters to study the regime under which spatially inhomogeneous steady states with striking patterns can develop. Our results also demonstrate that small initial data of being perturbations from the homogeneous equilibrium may lead to quite different structures of the solutions.
In Figure 2 and Figure 3, we plot numerical simulations of the evolutions of cell densities and chemical concentration to the stable steady state that has a single boundary spike. The boundary spike corresponds to the phenomenon of cellular aggregation. Figure 2 shows that the competing species u and v aggregates at the same location when chemical w is chemoattractive to both species. Figures 3 shows that u converges to a boundary spike at x = L and v to an inverted spike. This models the segregation phenomenon between species u and v through interspecific competitions, when u is chemo-attractive to w and v is chemo-repulsive. We observe in Figure 2 that the chemical gradient is a dominating mechanism for the coexistence of the competing species, while Figure 3 suggests that chemo-repulsion can be a leading mechanism that initiates the segregation. Figure 4 plots the emergence of multi-spiky solutions. We also observe the spatial segregation of the competing species. From the results in Figure 2 and Figure 3, one may expect the solution to develop spike at the locations where initial chemical concentration attains its maximum.
Motivated by the theoretical results in [16], we demonstrate the variations of pattern evolutions as domain size L changes in the set of simulations in Figure 5. Similar as in [16], we find that the large domain size L tends to increase the number of spikes or cell aggregates. It is also observed that multi-spiky solutions undergoes a coarsening process in which two interiors spikes merge into a single spike and new spike develops and eventually merges with another spike. The instability of the interior spikes may be responsible for spike merging, however, in contrast to the classical Keller-Segel model where the interior spike is known to be unstable, some interior spikes are stabilized for some proper length L-see the stable multi-spike steady states for L = 15. Moreover, rigourous analysis is required if one wants to find the critical threshold responsible for spike merging and stabilization. See the analysis in section 5 of [16] for the calculations and selection of wave modes for the single species chemotaxis model.  Finally, we present in Figure 6 the dynamical behaviors of some time-periodic spatially inhomogeneous patterns. According to our numerical studies, we surmise that the combination of chemo-attraction (χ > 0) and chemo-repulsion (ξ < 0) is necessary for the production of such time-periodic patterns. It is also interesting and important to find the period and rigour analysis and calculations are required for the justification.

Conclusions and Discussions
In this paper, we study the 3 × 3 chemotaxis system (1.1) of two competing species with Lotka-Volterra type kinetics. It is assumed that both species u and v move chemotactically to the same chemical w, which is also produced by the species. Tello and Winkler [19] showed that the homogeneous equilibrium (ū,v,w) is the global attractor of (1.1) if the chemoattraction rates χ Formations of stable single boundary spike on the left with L = 0.5, double spikes in the middle with L = 5 and shifting of interior spike to the boundary on the right with L = 7.
Interior spikes merge and shifts to the boundary on the left with L = 10, coarsening processes occur in the middle with L = 12 and merging of interiors spikes to form stable steady states on the right with L = 15. Figure 5: Spatial-temporal behaviors of u(x, t) of (1.1) for domains with different lengthes. The parameters are d 1 = 0.2, d 2 = 0.3, χ = 20, ξ = 50; initial data are (u 0 , v 0 , w 0 ) = (ū,v,w) + (0.01, 0.001, 0.02) cos(2.4πx). We observe that large interval length supports more interior spikes and merging or emerging of spikes occurs for large L. and ξ is small compared to their intrinsic growth rates µ 1 and µ 2 as in (1.5).
We obtain the global existence of classical solutions (u(x, t), v(x, t), w(x, t)) of (1.1) which are uniformly bounded for all t ∈ (0, ∞). For the stationary system (3.1), we show that the homogeneous equilibrium (ū,v,w) becomes unstable for χ > min k∈N + {χ k ,χ k } in the sense of Turing's instability driven by chemotaxis/ or advection. Then we apply the Crandall-Rabinowitz bifurcation theories to establish its nonconstant positive steady states. The stability or instability of these bifurcation solutions has also been obtained when the diffusion rates d 1 and d 2 are sufficiently large. Our main results asserts that, when both species diffuse fast, the small perturbation bifurcating solutions are unstable if the chemical decay rate is small or the domain size if large. Numerical simulations have been carried to illustrate the dynamical behaviors of solutions to (1.1) that exhibit complex spatio-temporal patterns. Our numerical results also include time periodic spatial patterns that has interesting structures. We also observe the merging and emerging of spikes in the spatial-temporal patterns.
It is an interesting and also important question to ask about the global existence of (1.2) for Ω ⊂ R N , N ≥ 2, when the assumption (1.5) is relaxed. The literature suggests that blow-ups can be inhibited by the degradation in the cellular kinetics, however, whether or not this is sufficient when chemo-attraction rates χ and ξ are large remains open so far, in particular for high space dimensions. It appears that the positivity of χ and ξ is required in the arguments of [19], and in light of our stability analysis, we surmise that the solution of (1.2) is always global if χ < 0 and ξ < 0, since chemo-repulsion has smoothing effect like diffusions. To prove this, one needs an approach totally different from that in [19].
There are also some other unsolved questions concerning the steady states of (3.1). First of all, the emergence of time periodic spatial pattern suggests that the homogeneous equilibrium loses its stability through Hopf bifurcation. Detailed calculations can be performed to study the structures of the bifurcating solutions. If χ 0 = min k∈N +χ k < min k∈N +χ k , (ū,v,w) loses its stability at χ 0 , but steady state bifurcation does not occur until χ passes min k∈N +χ k . It is a natural question ask about the existence of nonconstant steady states of (3.1) and also its driving mechanism when χ is between min k∈N +χ k and min k∈N +χ k .
The existence of non-existence of nontrivial solutions in multi-dimensional domain deserves exploring. In particular, rigourous construction of boundary spike and interior spike is another challenging problem that worths future attentions, even in one-dimensional domain. The stability of these spiky solutions is an important but very challenging problem that one can pursue in the future. We also want to mention that d 1 and d 2 are required large in the stability analysis, however, some extremely complicated and difficult calculations need to be done to remove this constraint.