PATTERN FORMATION OF THE ATTRACTION-REPULSION KELLER-SEGEL SYSTEM

In this paper, the pattern formation of the attraction-repulsion Keller-Segel (ARKS) system is studied analytically and numerically. By the Hopf bifurcation theorem as well as the local and global bifurcation theorem, we rigorously establish the existence of time-periodic patterns and steady state patterns for the ARKS model in the full parameter regimes, which are identified by a linear stability analysis. We also show that when the chemotactic attraction is strong, a spiky steady state pattern can develop. Explicit time-periodic rippling wave patterns and spiky steady state patterns are obtained numerically by carefully selecting parameter values based on our theoretical results. The study in the paper asserts that chemotactic competitive interaction between attraction and repulsion can produce periodic patterns which are impossible for the chemotaxis model with a single chemical (either chemo-attractant or chemo-repellent).


1.
Introduction. Chemotaxis is characterized by the sustained migration of cells in the direction of an increasing concentration of chemo-attractant or decreasing concentration of chemo-repellent, where the former is referred to as attractive chemotaxis and the latter to repulsive chemotaxis. Chemotactic pattern formation observed both experimentally and numerically exhibits rich interesting coherent structures such as aggregate, fruiting bodies, clusters, spirals, spots, rings, labyrinthine pattern, stripes and prorogating waves [1,4,11,14,16,31,55]. The patterns generated by chemotaxis models have been advocated to explain various biological processes/pattern formations, such as fish pigmentation patterning [36], angiogenesis in tumor progression [5], primitive streak formation [37], blood vessel formation [12], wound healing [40], and so on. The prototypical chemotaxis model was well-known as the Keller-Segel system [23] as follows u t = ∆u − ∇ · (λu∇v), x ∈ Ω, t > 0, v t = ∆v + αu − βv, x ∈ Ω, t > 0, (1.1) where u = u(x, t) denotes the cell density, v = v(x, t) represents the chemical concentration, and Ω is a bounded domain in R n . Here the attractive (resp. repulsive) chemotaxis corresponds to λ > 0 (resp. λ < 0), and |λ| ∈ R\{0} measures the strength of chemotactic response; parameters α, β > 0 denote the production and degradation rates of the chemical. The Keller-Segel system (1.1) describes the chemotactic interaction between cells and one chemical signal (either attractive or repulsive), and it has triggered a tremendous amount of studies in the past four decades, cf. [16,17,38,47] and references therein. However, in many biological processes, cells interact with combined attractive and repulsive chemical signals to achieve certain biological functions, for instance the leukocyte trafficking and neuronal pathfinding are governed by gradients of both chemical attractants and repellents [19]; see more examples in [6,13,28,35].
In this paper, we shall consider the spatiotemporal pattern formation of the following dimensionless attraction-repulsion Keller-Segel (ARKS) system x ∈ Ω, t > 0 v t = ∆v + αu − βv, x ∈ Ω, t > 0, w t = ∆w + γu − δw, x ∈ Ω, t > 0, ∂u ∂ν = ∂v ∂ν = ∂w ∂ν = 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) where u(x, t) denotes the cell density, v(x, t) represents the concentration of an attractive signal (i.e. chemo-attractant), and w(x, t) is the concentration of a repulsive signal (i.e. chemo-repellent); Ω ⊂ R n (n ≥ 1) is a bounded connected domain with a smooth boundary ∂Ω, and each of u, v and w satisfies a no-flux boundary condition; The parameters satisfy χ ≥ 0, ξ ≥ 0, α, β, γ, δ > 0, where χ and ξ measure the strength of the attraction and repulsion, respectively, and α, β, γ, δ represent the production and degradation rates of the two chemicals, respectively. The first cross-diffusive term in the first equation implies that the cell movement is directed up the chemo-attractant concentration, whereas the second cross-diffusive term indicates that cells move down the chemo-repellent concentration. The second and third equations in (1.2) state that chemo-attractant and chemo-repellent are released by cells and undergo natural decay. The attraction-repulsion chemotaxis model (1.2) was proposed in [28] to describe the aggregation of microglia observed in Alzheimer's disease and in [35] to address the quorum effect in the chemotactic process. The system (1.2) can also be regarded as a particularized version of multi-species models proposed in a recent work [18], but the analysis there does not apply to (1.2) due to the lack of a Lyapunov functional. The global existence, asymptotic behavior and steady states of solutions to the model (1.2) were established in [25] for a one dimensional spatial domain and in [48] for the multi-dimensional case in some parameter regimes, where the existence of steady states was obtained only for the case β = δ. It should be mentioned that when β = δ, the system (1.2) can be reduced to the classical Keller-Segel model (1.1) and hence conventional methods can be employed, see details in [48]. Traveling wave solutions of an attraction-repulsion chemotaxis model with volume-filling effect were recently studied in [39].
The purpose of this paper is to employ local and global bifurcation theory to study the time-periodic solutions and steady state solutions of (1.2) in the whole parameter regimes. We find that if β = δ, the interaction between attraction and repulsion may lead to time-periodic rippling patterns which were not possible for the classical Keller-Segel model (1.1) due to the existence of a time-monotone Lyapunov functional. We note that the parameter regimes for the pattern formation of the ARKS model (1.2) was also identified in [28] by a linear stability analysis without rigorous proof, where time-periodic patterns were found under the argument of short ranged attraction and long ranged repulsion (i.e. the spatial range (length) of the attractant is longer than that of the repellent). In this paper, we study the pattern formation of the ARKS model (1.2) where the attractant and repellent have the same spatial range and identify the full parameter regimes in which the existence of time-periodic and steady state patterns are rigorously proved. Our results indicate that time-periodic patterns still can develop if the spatial range of attractant and repellent are the same. Moreover we explicitly find the numerical time-periodic rippling pattern which was not shown in [28]. To the best of our knowledge, this is the first time that the existence of time-periodic patterns is rigorously established for reaction-diffusion systems with chemotactic effect. For reaction-diffusion systems without chemotactic effect, time-periodic patterns have been found in many well-known systems via Hopf bifurcations [26,56,58]. In [30], Hopf bifurcation and time-periodic patterns were obtained for a three-dimensional diffusive competition system with cross-diffusion. Note that time-periodic pattern is one of the spatiotemporal patterns predicted by Turing in his seminal work [49] sixty years ago, and Turing also pointed out that such pattern can only occur with three or more morphogens, which resonates the fact that Hopf bifurcation cannot occur for the Keller-Segel system with only one chemical. In this paper, we assume that all three species u, v, and w have the same diffusion rate equal to one for simplicity and brevity; however our analysis and results directly extend to the case where the diffusion rates of u, v, and w are different. We shall discuss a general attraction-repulsion Keller-Segel model with different diffusion rates in a forthcoming paper.
In numerical simulation of (1.2) with parameters in the time-periodic pattern regime, propagating rippling waves which pass through one another without interference have been found. It is worth underlining that such patterns have been discovered experimentally in the early stage of starvation-induced fruiting body development of the starved myxobacteria ( [20,21,57]) where the contact interaction between cells was postulated to be the underlying mechanism. However our results indicate that the chemotactic interaction of attraction and repulsion may offer an alternative mechanism, as indicated in [45].
For the steady state solutions, we recall that the existence of steady states to the Keller-Segel system was first analyzed in [42] by using local and global bifurcation techniques, and bifurcation of spatially heterogeneous steady state solutions were also obtained in [15,29,32]. Recently a global bifurcation analysis for steady states was carried out in [53] for a class of regularized models of the classical Keller-Segel models surveyed in [16], see also related work in [44,51,52]. On the other hand, the existence of spiky steady state solutions for chemotactic systems in the small diffusion coefficient case has been shown in [22,24] and many other work, see related surveys [33,34]. The steady state solutions of the attraction-repulsion Keller-Segel model has not been considered in afore-mentioned works. Here we prove the existence of non-constant steady state solutions with bifurcation theory, and numerical simulation in the paper also shows spiky pattern although the diffusion coefficients are not small.
By integrating the first equation with the boundary conditions, one immediately observes that the cell mass is preserved for all t > 0. Hence the value ofū = 1 |Ω| Ω u(x, t)dx, will be prescribed in our study. Then the steady state solutions of (1.2) satisfy It is clear that (ū,v,w) is a constant steady state (i.e. equilibrium) solution of (1.3), wherev = αū/β,w = γū/δ. In the paper we shall show, givenū > 0, the system (1.3) has positive solutions for a suitable range of χ. The solutions will be the ones bifurcating from the constant steady states (ū,v,w). The chemotactic coefficient χ and ξ will be used as the bifurcation parameter, and the bifurcating solutions have arbitrarily fixed total cell mass. The remaining part of the paper is organized as follows. In Section 2, we present the basic linear stability analysis for the constant steady state (ū,v,w); In Sections 3 and 4, we consider Hopf bifurcations and steady state bifurcations respectively; Discussion on pattern formation and numerical simulations are given in Section 5; and Section 6 is the appendix which collects the general Hopf bifurcation theorem, the local and global steady state bifurcation theorem, and the Routh-Hurwitz criterion for cubic equations, which are the analytic tools used in the paper.
2. Linearized stability analysis. Before proceeding, we introduce some important notations used in what follows. By N, we denote the set of all the positive integers, and N 0 = N ∪ {0}. We use the Sobolev spaces X = u ∈ W 2,p (Ω) : ∂u ∂ν = 0 , Y = L p (Ω), Y 0 = u ∈ Y : Ω u(x)dx = 0 , and X 0 = X Y 0 . We always assume that p > n.
Linearizing the system (1.3) about the constant equilibrium solution (ū,v,w), we obtain an eigenvalue problem First we show that eigenvalue problem (2.1) can be reduced to a sequence of matrix eigenvalue problems.
Lemma 2.1. Let {λ n } be the sequence of eigenvalues of −∆ with Neumann boundary condition, such that 0 = λ 0 < λ 1 ≤ λ 2 ≤ · · · , and let y n (x) be its corresponding eigenfunction, n ∈ N. Let (ū,v,w) be a positive constant equilibrium. Define Then 1. If µ is an eigenvalue of (2.1), then there exists n ∈ N, such that µ is an eigenvalue of (2.2); if (a n , b n , c n ) is an eigenvector associated with µ for (2.2), then (a n , b n , c n )y n is an eigenvector associated with µ for (2.1). if there exists an n ∈ N, such that A n has at least one eigenvalue with nonnegative real part.
In view of Corollary 2.2, we define the quantity (2.7) Applying Corollary 2.2 and (2.4)-(2.7), we have the following rough results regarding the stability/instablity of the constant equilibrium (ū,v,w).
Remark 2.4. The stability condition (2.8) in Proposition 2.3 generalizes a result in [48], which asserts that in the absence of attraction (χ = 0) or when the repulsion dominates over attraction (in the sense that ξγ ≥ αχ) under an additional assumption that δ = β (same death rate), the constant equilibrium is stable and hence there is no pattern formation. This is consistent with biological intuition since the repulsion acts as diffusion which helps to stabilize the system. However if the attraction dominates over the repulsion, spatiotemporal pattern formation may be expected.
For that purpose we define the set to be the steady state bifurcation curve, and to be the Hopf bifurcation curve (see [50]). The curves S and H are where the constant equilibrium is neutrally stable, i.e. the linearized equation has an eigenvalue withy zero real part. Since the expression of a 0 (χ, p) and T (χ, p) are linear in χ, then S and H are graphs of functions defined as follows: and Then the functions χ S (p) and χ H (p) have the following properties.

PING LIU, JUNPING SHI AND ZHI-AN WANG
Next we differentiate χ H (p) and obtain (2.14) Then Hence any critical point of χ H (p) must be a local minimum, and furthermore it must be unique since lim p→0 + χ H (p) = ∞ and lim p→∞ χ H (p) = ∞. Therefore this critical point must be a global minimum point.
Now we can state a more precise stability result for the constant equilibrium (ū,v,w).
If χ * = min n∈N χ S (λ n ), then the linearized equation (2.1) has a zero eigenvalue at χ = χ * , and a bifurcation of nontrivial steady state solutions of (1.2) occurs if some transversality conditions are met; similarly if χ * = min n∈N χ H (λ n ), then the linearized equation (2.1) has a pair of purely imaginary eigenvalues at χ = χ * , and a Hopf bifurcation generating a family of periodic orbits of (1.2) occurs if some transversality conditions are met. We shall give more detailed analysis for both types of bifurcations in the next two sections. However one can see that whether the constant equilibrium loses stability to a spatial stationary pattern through a steady state bifurcation or to a spatial time-periodic pattern through a Hopf bifurcation is a delicate problem depending on the shape of the functions χ = χ S (p) and χ = χ H (p).
3. Hopf bifurcation. In this section we prove the existence of periodic orbits of (1.2) for a certain parameter range using the Hopf bifurcation theorem. For the reader's convenience, a review of the Hopf bifurcation theorem for quasilinear parabolic systems is given in the appendix (subsection 6.2).
First we identify the parameter range of α, β, δ, γ, ξ so that (2.1) could have purely imaginary eigenvalues. In the following lemma, we first show that if β > δ and ξγū is large, then χ H (p) < χ S (p) for some p values, which indicates the possibility of a Hopf bifurcation.
Lemma 3.1. Let (ū,v,w) be a positive constant equilibrium and define wherep is the unique positive root of the equation Proof. By (2.12) and (2.13), we have Then it is apparent that when δ ≥ β, χ H (p) > χ S (p) for all p > 0. For the case β > δ, with A = ξγū, χ H − χ S can be expressed as Hence χ H (p) − χ S (p) = (<)0 for p > 0 if and only if Q 1 (p) = (<)0. Since Q 2 (p) is a convex and increasing cubic function for p > 0, and Q 3 (p) is linear, then Q 1 (p) = 0 has either 0, 1 or 2 roots when A <, = or > A * , for a critical value A * . At the critical value A = A * , we have a unique rootp such that Q 2 (p) = Q 3 (p) and Q 2 (p) = Q 3 (p). Hence at p =p, Simplifying (3.6) we obtain Since (3.7) has a unique positive rootp, then the critical value A * is Hence when A > A * , Q 1 (p) = 0 has two roots, denoted by p 1 and p 2 . Noticing that Q 1 (p) is convex, then the stated results follow. Now we can determine the potential Hopf or steady state bifurcation points by using the Hopf bifurcation curve and steady state bifurcation curve.
Proposition 3.2. Recall that {λ n } are the eigenvalues of −∆ with Neumann boundary condition, such that 0 = λ 0 < λ 1 ≤ λ 2 ≤ · · · , and y n (x) is the corresponding eigenfunction for λ n . Let χ S (p) and χ H (p) be defined as in (2.12) and (2.13), respectively. For n ∈ N, define Then for fixed α, β, δ, γ, ξ, 1. the eigenvalue problem (2.1) has an eigenvalue µ = 0 if and only if χ = χ S n for some n ∈ N, and the corresponding eigenfunction is V n y n , where V n satisfies A n V n = 0 with A n defined as in (2.2); 2. the eigenvalue problem (2.1) has a pair of purely imaginary eigenvalues µ = ±iν if and only if (3.5) is satisfied, χ = χ H n for some n ∈ N and p 1 < λ n < p 2 where p 1 and p 2 are defined in Lemma 3.1, and the corresponding eigenfunction is V ± n y n , where V ± n satisfies A n V ± n = ±iνV ± n . Proof. The first part follows directly from part 1 of Lemma 2.1 and the fact that a zero eigenvalue of A n is equivalent to a 0 (χ, λ n ) = 0. For the second part, (2.1) has a pair of purely imaginary eigenvalues if and only if for some n ∈ N, A n has a pair of purely imaginary conjugate eigenvalues. Suppose that A n has eigenvalues iν, −iν and θ where ν, θ ∈ R and ν = 0. Then for a 0 , a 1 and a 2 defined in (2.4)-(2.6), we have from the Routh-Hurwitz criterion (see Appendix 6.3) Since from (2.4), a 2 (χ, λ n ) = 3λ n + β + δ > 0, hence θ < 0 and a 0 (χ, λ n ) > 0. On the other hand, it is clear that T (χ, λ n ) = a 2 (χ, λ n )a 1 (χ, λ n ) − a 0 (χ, λ n ) = 0. Hence χ must satisfy χ H (λ n ) = 0 and χ S (λ n ) > 0. From Lemma 3.1, this can only occur when (3.4) is satisfied, χ = χ H n for some n ∈ N and p 1 < λ n < p 2 . Finally at a potential Hopf bifurcation point χ H n , we show that a nondegeneracy condition (see (H3) in Theorem 6.1) holds. Lemma 3.3. Suppose that (3.5) is satisfied and p 1 < λ n < p 2 where p 1 and p 2 are defined as in Lemma 3.1, and χ H n is defined as in Proposition 3.2. Then for χ near χ H n , A n has a unique eigenvalue σ(χ) + iν(χ) such that σ(χ H n ) = 0 andν(χ H n ) = ν 0 > 0, and σ (χ H n ) > 0. Proof. The existence of a unique eigenvalue σ(χ)+iν(χ) follows from the continuous dependence of eigenvalues on the parameter χ. For χ near χ H n , let the three roots (3.12) Now we state our main result in this section on the existence of nontrivial periodic orbits of (1.2).
Theorem 3.4. Assume that the parameters α, β, δ, γ, ξ,ū satisfy (3.4), and let χ H n and χ S n be defined as in (3.8). Suppose that (A1) for some j ∈ N, λ j is a simple eigenvalue of −∆ in Ω with Neumann boundary condition, and the corresponding eigenfunction is y j (x); (A2) The eigenvalue λ j satisfies
Proof. We follow the general framework given in Section 6.1. Denote U = (u, v, w). Define (3.17) Then the chemotaxis system (1.2) is in a form of (6.6) with a and a 0 defined as in (3.17), a j ≡ 0, b ≡ 0, f ≡ 0 and g ≡ 0. To apply Theorem 6.1, we make a change of variables (u, v, w) = (ū +ũ,v +ṽ,w +w). DenoteŨ = (ũ,ṽ,w) andŪ = (ū,v,w). Then after the transformation, (1.2) becomes x ∈ ∂Ω, t > 0, (3.18) and for χ ≥ 0, (ũ,ṽ,w) = (0, 0, 0) is an equilibrium. The linearized equation (6.7) atŪ is exactly (2.1). Hence to apply Theorem 6.1, we only need to verify conditions (H1)-(H3) are satisfied. From Proposition 3.2, at χ = χ H j , the matrix A j has a pair of purely imaginary eigenvalues ±iν 0 , with ν 0 given by (3.16) following (3.9) and (2.5). Since λ n is a simple eigenvalue of −∆, and χ H j = χ H n for n = j, then µ = ±iν 0 is a pair of simple eigenvalues of (2.1), and hence (H1) and (H2) (for k ∈ Z\{±1, 0}) are satisfied. Since χ H j = χ S n for any n ∈ N, then 0 cannot be an eigenvalue for (2.1) when χ = χ H j , thus (H2) is also satisfied for k = 0 . Finally Lemma 3.3 implies that (H3) is satisfied. Therefore we can apply Theorem 6.1 to obtain the results stated. The form (3.15) can be obtained by solving the linearized parabolic system. Theorem 6.1 shows that the occurrence of a Hopf bifurcation in (1.2) depends on several factors: (i) parameters need to satisfy (3.4); (ii) the spatial domain Ω (or indeed the eigenvalue λ j ) need to satisfy conditions (A1)-(A3). The parameter choice (3.4) can be achieved when all parameters are properly chosen, and one example is shown in Fig. 1. The conditions (A1) and (A3) hold for generic spatial domains (i.e. all domains except a zero measure set), but (A2) means that the size and shape of the domain have to be chosen properly. For a parameter set satisfying (3.4), if the domain Ω is so small that λ 1 > p 2 , then (A2) is not met.
When the Hopf bifurcations do occur in (1.2), results in this section give precise algorithm to determine the exact bifurcation points and the nature of the Hopf bifurcation. The parameter selection is achieved by using (3.4), and to satisfy (A2), one can use (3.1) and (3.2) and properly choose λ j . The expression (3.15) gives the oscillation frequency, and the eigenfunction y j (x) gives the spatial profile of the oscillation. We also emphasize that Theorem 6.1 holds for domain Ω in any dimension n ≥ 1. When n = 1 and Ω = (0, l), the eigenfunction is y j (x) = cos(jπx/l) which shows the time-periodic orbit is also periodic in space. 4. Steady state bifurcation. In this section we prove the existence of nontrivial steady state solutions of (1.2), or equivalently solutions of (1.3). We use a bifurcation approach similar to the one in [53]. We first will consider the case of a general bounded domain Ω ⊂ R n (n ≥ 1), and then we give some more specific results for the case Ω = (0, l) and n = 1. We first prove some nonexistence results of (1.3) for small χ ≥ 0. Proof. The first part is indeed proved in [48, Proposition 2.3], so we only prove the second part. When χ = 0, the system (1.3) is reduced to (4.1) Note that the first equation of (4.1) can be written as ∇[u∇(ln u + ξw)] = 0, then it follows that ln u(x) + ξw(x) ≡ C 1 or u(x) ≡ Ce −ξw(x) for a constant C ∈ R. By using the mass conservation equality and the equation of w, we obtain that w(x) satisfies It is easy to verify that for a fixed A > 0, g(w)(w −w) < 0 for any w > 0 and w =w = γū/δ. From the maximum principle, it follows that 0 ≤ w(x) ≤w for x ∈ Ω. Integrating the equation in (4.2) yields Ω g(w(x))dx = 0. Hence w(x) ≡w and consequently (u, v, w) ≡ (ū,v,w).
The main result for global bifurcation of steady state solutions in the attractionrepulsion chemotaxis system is as follows.
Proof. Define a mapping F : We apply Theorem 6.2 to the equation F (u, v, w, χ) = 0 at (ū,v,w, χ S j ). Clearly, F (ū,v,w, χ S j ) = 0, and F is continuously differentiable. We verify the conditions in Theorem 6.2 in the following steps.
Step 1. F U (ū,v,w, χ S j ) is a Fredholm operator with index zero, and the kernel space N (F U (ū,v,w, χ S j )) is one-dimensional, where U = (u, v, w). By using the same argument as in Lemma 2.3 in [53], one can show that the linear operator F U (ū,v,w, χ S j ) : . Then from Lemma 2.1, there exists j ∈ N such that zero is an eigenvalue of A j (defined in (2.2)), and the corresponding eigenvector is (ā j ,b j ,c j )y j = λ j + β, α, γ(λ j + β) λ j + δ y j . (4.7) From the condition (A1), the eigenvector is unique up to a constant multiple. Thus one has N (F U (ū,v,w, χ S j )) = span{(ā j ,b j ,c j )y j }, which is one-dimensional.
Step 2. F χU (ū,v,w, χ S j )[(ā j ,b j ,c j )y j ] ∈ R(F U (ū,v,w, χ S j )). We claim that the range space R(F U (ū,v,w, χ S j )) can be characterized as follows: where (a * j , b * j , c * j ) is a non-zero eigenvector for the eigenvalue µ = 0 of A T n (the transpose of A n defined in (2.2)): and its adjoint operator Then we have where ·, · is the inner product in [L 2 (Ω)] 3 . This proves that if Since (4.12) defines a codimension-1 set in Y 0 × Y 2 × R, and we know that then R(F U (ū,v,w, χ S j )) must be in form of (4.8). Now noting that ). Finally we can use the arguments in Lemma 2.3 of [53] to show that for any (u, v, w, χ) ∈ X 3 × R, the operator F U (u, v, w, χ) is a Fredholm operator with index zero. Then all conditions in Theorem 6.2 are verified, and hence we obtain that solutions bifurcated from (ū,v,w, χ S j ) are on a connected component C j of the set of nontrivial solutions of (1.3). We shall show that all solutions on C j are positive. This is apparently true for solutions near the bifurcation point (ū,v,w, χ S j ) sincē u > 0,v > 0 andw > 0. From the equation of v and w in (1.3), one can see that if u is nonnegative, then v and w must be positive in Ω. Hence if not all solutions on C j are positive, then there exists (u, v, w, χ) ∈ C j , such that v(x) > 0, w(x) > 0 in Ω, and u(x) = 0 for some x ∈ Ω. But from the equation of u in (1.3) (which is linear in u), and strong maximum principle, we have a contradiction. Hence all solutions on C j are positive.
For the last part, since the only nonnegative solution of (1.3) is (ū,v,w) when χ = 0 from Lemma 4.1, and χ = 0 cannot be a bifurcation point from Lemma 2.1, then the continuum C j cannot be extended to χ = 0. Hence proj χ C j ⊂ (0, ∞). The result for the case of β = δ is from the other result in Lemma 4.1. This completes the proof.
It is difficult to obtain more detailed information of the global branch C j in the general setting as in Theorem 4.2. When the spatial domain is one-dimensional, a more precise description of C j can be obtained. In the remaining part of this section, we consider (1.3) in the case n = 1 and Ω = (0, l), that is To obtain a more precise global bifurcation diagram, we shall prove the following a priori estimates of positive solutions of (4.13) (the proof is similar to the one in [53]): Lemma 4.3. Suppose that α, β, γ, δ > 0 and χ, ξ ≥ 0, and (u(x), v(x), w(x)) is a non-negative solution of (4.13). Then either (u( (4.14) Proof. Clearly the constant (ū,v,w) is a solution to (4.13). Next we assume that (u(x), v(x), w(x)) is a non-constant nonnegative solution of (4.13). Assume that v(x 0 ) = min Integrating the v-equation in (4.13) and using the boundary condition, we have

Proof.
We can verify that all conditions in Theorem 6.3 are also satisfied, then Theorem 6.3 implies that C 1 has a subset C + 1 , which contains the s > 0 part of the curve Γ 1 in Theorem 4.2 but excludes the s < 0 part. Note that it is easy to verify that χ S 1 can be calculated as in (4.15). From Lemma 4.1 and β = δ, we know that proj χ C + 1 ⊂ (ξγ/α, ∞). On the other hand, from Lemma 4.3, all positive solutions (u, v, w) are bounded for χ in a compact subset of (0, ∞). Hence when C + 1 is not compact, we must have proj U C + 1 ⊃ (χ S 1 , ∞). Furthermore when β = δ, one has that χ S 1 > ξγ/α from Lemma 2.5 and (4.13) is then reduced to with w = (γ/α)v. Then we can apply Theorem 3.1 in [53] to obtain other assertions in the theorem.
We conjecture that the results in Theorem 4.4 hold at least for the case δ ≥ β, but a rigorous proof is difficult to obtain. Note that at the same bifurcation point χ = χ S 1 , another family of steady state solutions (u, v, w) satisfying u , v , w > 0 also bifurcates from the line of trivial solutions {(ū,v,w, χ)} by taking the portion of Γ 1 with s < 0, hence a pitchfork bifurcation for (4.13) with β = δ occurs at The stability of the bifurcating steady state solutions can be determined by calculation using formulas in [27,43], but the calculation is lengthy hence not included here. From the numerical simulations in Section 5, the monotone steady state near χ = χ S 1 appears to be locally asymptotically stable when δ ≥ β. Finally we remark that in the results of Sections 2-4, we have used the attraction strength parameter χ as the bifurcation parameter. Indeed χ and ξ are both important parameters which measure the strength of the attraction and repulsion respectively. In this section we show that the bifurcation analysis in Sections 2-4 is still valid when we use ξ as the bifurcation parameter instead of χ. Since the proofs are all similar we only state the results for the stability analysis without proof. Similar to the functions χ S (p) and χ H (p) in Section 2, from (2.10) and (2.11) we can define and Then the properties of ξ S (p) and ξ H (p) are as follows: Lemma 4.5. Suppose that α, β, γ, δ, χ,ū > 0 are fixed, and let ξ S (p) and ξ H (p) be defined as in (4.17) and (4.18).
We notice that ξ S (p * ) can be explicitly calculated: where Both ξ S (p * ) and ξ H (p * ) can be computed if the parameter values are given. In Fig.  2, graphs of ξ = ξ S (p) and ξ = ξ H (p) for several different sets of parameters are shown. In particular if the graphic scenario in the right panel of Fig. 2 occurs, Hopf bifurcations are possible as we have shown in Section 3. Indeed similar to Lemma 3.1, we have the following result characterizing the parameter regimes for timeperiodic patterns.
Lemma 4.7. Let (ū,v,w) be a positive constant equilibrium and define wherep is the unique positive root of the equation 1. If either δ ≥ β or β > δ and χαū < A * , then ξ S (p) > ξ H (p) for all p > 0.
5. Numerical simulations. The classical Keller-Segel model with one chemical does not admit time-periodic pattern solutions due to the existence of a timemonotone Lyapunov functional. In this paper, we prove that the attraction-repulsion Keller-Segel (ARKS) model (1.2), in which two competing chemicals (i.e. chemoattractant and chemo-repellent) interact with one species, may produce time-periodic patterns. Our study shows that the death rates of two competing chemicals, denoted by β and δ, are critical in creating periodic patterns which arise only if the chemo-attractant has larger death rate than the chemo-repellent (i.e. δ < β) and repulsion is large in the sense that ξγū > A * , where A * is given in Lemma 3.  Fig.  1 (c). The initial conditions are set as a small random perturbation of the homogeneous steady state (20, 5/4, 20).
In this section, we numerically explore the spatiotemporal patterns generated by the ARKS models (1.2) in different parameter regimes identified in the paper. The PDEPE computing package in Matlab is implemented to solve the system (1.2) in one dimension. Lemma 3.1 presents complete parameter regimes in which various bifurcations occur, which are illustrated in Fig. 1. We first choose parameter values as in Fig.  1 (c) in an interval Ω = (0, l), then we obtain that Fig. 1 (c), we see that if we select χ slightly larger than the smallest Hopf bifurcation critical value χ H j given in (2.13) with p = (jπ/l) 2 , then only one Hopf bifurcation occurs and the constant steady state loses the stability through the Hopf bifurcation, and a time-periodic patterns is expected in this case (see Theorem 3.4). Indeed, in the simulation shown in Fig. 3, we choose n = 1 and Ω = (0, 3), and then it can be easily calculated that the smallest Hopf bifurcation value is χ H 1 = χ H (π 2 /9) ≈ 2.97, while the next Hopf bifurcation value is χ H 2 = χ H (4π 2 /9) ≈ 3.27.  Here we choose χ = 3.07 and set the initial condition as a small random perturbation of the homogeneous steady state (u s , v s , w s ) = (20, 5/4, 20). Then the spatiotemporal pattern of the cell density u is plotted in Fig. 3 (a), in which we observe highly organized rippling periodic waves arise from an initially almost homogeneous colony and periodically reverse their directions. In a three dimensional visualization shown in Fig. 3 (b), we see that two wave crests pass through one another without interference when colliding, different from the typical chemotactic waves which coalesce upon collision (cf. [16,54]). In the numerical simulation plotted in Fig. 3-5, time-periodic patterns are found by a careful selection of parameter values according to Fig. 1 (c) and Lemma 3.1 such that only Hopf bifurcations occur. It is also of interest to consider the pattern formation where only steady state bifurcations occur or both bifurcations occur. In the second simulation, we consider the parameter values given in Fig. 1 (a) so that χ S (p) = p + 2, χ H (p) = 4p + 7 + 2 p .
In this case, no Hopf bifurcations can occur according to Lemma 3.1 and Theorem 3.4, but a steady state bifurcation occurs at every χ S j = χ S (λ j ). Again we choose Ω = (0, 3). Then χ S 1 = χ S (π 2 /9) = π 2 /9 + 2 ≈ 3.10. In the numerical simulation shown in Fig. 6 (a), we choose χ = 3.20 and set the initial condition as a small random perturbation of the homogeneous steady state (1, 1, 1). We observe that a spatially monotone decreasing steady state pattern arises from the initial uniform colony. But the pattern solution u shown in Fig. 7 (a) is not spiky since the value of χ is not large (just slightly bigger than the bifurcation value). If we increase the value of χ, we find from the numerical simulation shown in Fig. 6 (b) that the spiky steady state pattern emerge, and the solution profile is monotone decreasing in space as shown Fig. 7 (b). These numerical simulations are entirely consistent with our theoretical results in Theorem 4.4.
6.1. Hopf bifurcation for quasilinear parabolic systems. We recall the Hopf bifurcation theorem for a quasilinear parabolic system from Amann [3] (see also Crandall and Rabinowitz [8], Da Prato and Lunardi [9]). For simplicity we only consider flux type boundary conditions not Dirichlet ones.
The theorem here is essentially the same as the theorem in [3] with some rephrasing. As mentioned in [3], the C ∞ smoothness assumption can be considerably weakened, which we do not need here as our system (1.2) does satisfy the C ∞ smoothness assumption.
6.2. Bifurcation of steady state solutions. We first recall some bifurcation theorems for an abstract equation where F : R × X → Y is a nonlinear differentiable mapping, and X, Y are Banach spaces. In the following N (L) and R(L) are the null space and the range of a linear operator L; F u is the Fréchet partial derivatives of F with respect to argument u, and F λu is the mixed Fréchet partial derivatives of F with respect to arguments u and λ.
This result follows directly from [44,Theorem 3.3], since the Augmon's condition is satisfied from Case 1 of [44,Remark 2.5.5]. Hence the key condition on the Fredholm property of the linear operator in Theorems 6.2 and 6.3 is satisfied in this setting. We also remark that the differentiability of the norm function is also satisfied for L p space (see the remark after [44,Theorem 4.4]).
6.3. Cubic roots. Next we consider the roots of the cubic equation with constants A, B and C λ 3 + Aλ 2 + Bλ + C = 0, A > 0.
Then the Routh-Hurwitz criterion (e.g. see [31]) says that all the roots are negative or all the roots have negative real parts if and only if A > 0, C > 0, AB − C 2 > 0.
Since A > 0, there is at least one root with negative real part. There are two cases to consider. Case 1. C = −λ 1 λ 2 λ 3 < 0.
(1) If three roots are all real with λ 1 < 0, then λ 2 λ 3 < 0. This indicates that either λ 2 or λ 3 is a positive number. Hence the equation has a positive root.