Time-periodic and stable patterns of a two-competing-species Keller-Segel chemotaxis model: effect of cellular growth

This paper investigates the formation of time--periodic and stable patterns of a two--competing--species Keller--Segel chemotaxis model with a focus on the effect of cellular growth. We carry out rigorous Hopf bifurcation analysis to obtain the bifurcation values, spatial profiles and time period associated with these oscillating patterns. Moreover, the stability of the periodic solutions is investigated and it provides a selection mechanism of stable time--periodic mode which suggests that only large domains support the formation of these periodic patterns. Another main result of this paper reveals that cellular growth is responsible for the emergence and stabilization of the oscillating patterns observed in the $3\times3$ system, while the system admits a Lyapunov functional in the absence of cellular growth. Global existence and boundedness of the system in 2D are proved thanks to this Lyapunov functional. Finally, we provide some numerical simulations to illustrate and support our theoretical findings.


Introduction and preliminary results
In this paper, we continue our study in [52] of the following system of (u, v) = (u(x, t), v(x, t)) x ∈ Ω, t > 0, τ w t = ∆w − λw + u + v, x ∈ Ω, t > 0, ∂u ∂n = ∂v ∂n = ∂w ∂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.1) where d i , µ i , a i , i = 1, 2 and λ are positive constants and τ is a nonnegative constant; Ω is a bounded domain in R N , N ≥ 1 with smooth boundary ∂Ω. (1.1) is a chemotaxis model proposed by Tello and Winkler in [49] to study the population dynamics of two competitive biological species attracted by the same nutrition subject to Lotka-Volterra dynamics. Here u(x, t) and v(x, t) represent population densities of the two competing species at space-time location (x, t) ∈ Ω × R + , while w(x, t) denotes concentration of the attracting chemical. It is assumed that both species direct their movement chemotactically along the gradient of chemical concentration over the habitat, hence both χ and ξ are assumed to be positive constants.
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 in which a i , i = 1, 2 measure the levels of inter-specific competition and µ i , i = 1, 2 interpret intrinsic cellular growth. The chemical is produced by both species at the same rate with no saturation and it is consumed by certain enzyme at the rate of λ meanwhile.
Chemotaxis is the oriented movement of cellular organisms towards the high concentration region of a chemical released by the cells. One of the most interesting phenomena in chemotaxis is the cellular aggregation during which initially homogeneously distributed cells aggregate and develop into a fruiting body. For example, during the first phase of its developmental cycle, Dictyostelium discoideum exists as single amoeboid cell, then it differentiates into multiple cells which eventually aggregate into a multicellular organism after the growth phase. In Dictyostelium chemotaxis, it was discovered that the aggregating cells of D. discoideum are attracted by a chemical called cyclic AMP (cAMP), which is synthesized and released by the cells periodically. See the review papers [17,18]. It is of great interest to both biologists and mathematicians to understand the initiation and formation of the self-organized oscillating patterns. Mathematical modeling of chemotaxis dates back to the pioneering works of Patlak [44] and 35,36], where a group of parabolic reaction-diffusion systems have been proposed to describe the spatial-temporal behaviors of cellular distribution and chemical concentration. Diffusion models the random cellular movements and chemical diffusions, the advection term describes the chemotactic cellular movement, and the kinetics interpret the cellular birth-death and chemical degradation-creation. A substantial literature on the modeling and mathematical analysis of chemotaxis has developed since the appearance of works of Keller and Segel. See the survey papers [22,23,24,50] and the references cited therein. We would like to point out that most papers in literature focus on the studies of chemotaxis model with single bacteria and one chemical stimulus.
From the viewpoint of linearized stability analysis, chemotaxis destabilizes spatially homogeneous solutions for reaction-diffusion systems, while diffusion stabilizes homogeneous solutions (unless Turing's instability occurs), therefore one can expect the emergence of spatially inhomogeneous when chemotaxis rate is large. In [52], the authors investigate existence and stability of nontrivial positive steady states to (1.1) over (0, L) through rigorous bifurcation analysis. Moreover, we provide a selection mechanism of stable wavemode to predict the spatial profile of the stable stationary patterns. Numerical simulations in [52] verify the theoretical findings and suggest that system (1.1) over (0, L) also admits time periodic spatial patterns for properly chosen parameters. It is surmised in [52] that stable oscillating patterns emerge since (ū,v,w) loses its stability through Hopf bifurcation. One of the goals of our current work is to rigorously investigate the formation of stable time periodic patterns. In particular we show that chemotaxis and cellular kinetics are responsible for the formation of temporal oscillating patterns. We want to mention that the phenomenon of time periodic oscillations is important not only for reaction-diffusion models in biological and ecological systems ( [15] e.g.), but almost all other dynamical systems of scientific disciplines such as fluid mechanics [27,28,29,46], lasers [21] etc. In particular, for the effect of chemotaxis on bacterial strategies, see the discussions in [48] and the references cited therein.
The rest of this paper is organized as follows. In Section 2, we analyze the linearized stability of (ū,v,w) in terms of chemotaxis rate χ. It is shown in Proposition 2.1 that this homogeneous solution loses its stability as χ surpasses a threshold value χ 0 , which is the minimum of bifurcation values χ S k and χ H k over N + , where χ S k and χ H k are chosen such that the stability matrix (2.3) of (ū,v,w) has zero or purely imaginary eigenvalues respectively. Section 3 is denoted to the rigorous Hopf bifurcation analysis of (1.1) over (0, L), for which time periodic spatial patterns are established-see Theorem 3.1. Our existence results employ Hopf bifurcation theorem for parabolic systems from [2,14] etc. We also investigate the stability of these periodic solutions and establish a selection mechanism for stable oscillating patterns in terms of system parameters. In Section 4, we study the effect of cellular growth on the spatial-temporal dynamics of (1.1). In particular, we find that when µ 1 = µ 2 = 0, (1.1) does not admit time periodic patterns. Our argument is based on the construction of a time-monotone Lyapunov functional to (1.1). Moreover, global existence of this problem over 2D is also established provided that the initial cellular population is not too large-see Theorem 4.2. Section 5 presents various numerical studies that support our theoretical findings. Finally, we discuss our results and propose some open problems for future studies in Section 6.

Linearized stability analysis of homogeneous steady state
In the mathematical analysis of pattern formation in reaction-diffusion systems, the principle of exchange of stability ( [14,47,45] e.g.) is often employed to determine when bifurcation occurs for the family of evolution equations. Generally speaking, the principle states that when a spatially homogeneous solution loses its stability as a parameter crosses a threshold value, there may exist stable spatially inhomogeneous solutions to the system. In particular if the homogeneous solution loses stability through a pair of complex conjugate eigenvalues crossing the imaginary axis, one may expect, under suitable but reasonable technical conditions, that there exist time periodic solutions to the evolution equations. Moreover this principle usually gives a qualitative relationship between the shape of bifurcating curve (such as its turning direction) of solutions and their stabilities.
In this paper we are interested in studying time periodic solutions to (1.1), in contrast to the stable steady states investigated in [52]. For the simplicity of our calculations and without much loss of generality we shall confine our attention to system (1.1) over one-dimensional interval One of our goals is to explore how cellular kinetics effect the formation of spatially inhomogeneous positive solutions to (2.1). For this purpose, we adopt the principle of exchange of stability in the context of Hopf bifurcation, i.e., the bifurcation of a family of time periodic solutions from (ū,v,w). To begin with, we carry out the linearized stability of (ū,v,w) to investigate the spatial-temporal of dynamics to (2.1) around this homogeneous solution. Some of our stability results have been obtained in [52] and more details are needed for our purpose in this paper, therefore we include them here for the completeness and consistency of our arguments.
Linearizing (2.1) by setting U , V and W u =ū + U, v =v + V, w =w + W, with 0 < 1 and substituting these perturbations into (2.1), we obtain Now we look for solutions of (2.2) in the form (U, V, W ) = (C 1 , C 2 , C 3 )e σt+ikx , where k is the wavemode vector and σ is the growth rate of the perturbations respectively; moreover |k| 2 = ( kπ L ) 2 thanks to the L 2 eigen-expansions and C i (i = 1, 2, 3) are constants to be determined. Substituting these solutions into the linearized system above gives us the following problem where the matrices A 0 and D are Or we have the following stability matrix equivalent to (2.2) whose eigen-value is σ According to the Routh-Hurwitz conditions or Corollary 2.2 in [39], the real parts of all eigenvalues to (2.4) are negative (hence (ū,v,w) is locally stable) if and only if for all k ∈ N + , while there exist some eigenvalues with a nonnegative real part if one of the conditions above fails for some k ∈ N + . Moreover since α 2 (k) > 0, we will always have α 1 (χ, k) > 0 whenever α 0 (χ, k) > 0 and α 1 (χ, k)α 2 (k) − α 0 (χ, k) > 0. Therefore, the stability criterion above implies that (ū,v,w) is unstable if there exists k ∈ N + such that either α 0 (χ, k) < 0 or α 1 (χ, k)α 2 (k) − α 0 (χ, k) < 0. The following results are proved in [52].

5)
and We want to point out that Proposition 2.1 holds for multi-dimensional bounded domains R N with ( kπ L ) 2 being replaced by the k-th Neumann eigenvalue of −∆. According to Proposition 2.1, the spatially homogeneous steady state (ū,v,w) loses its stability at χ 0 = min k∈N + {χ S k , χ H k }. It is natural to expect that as χ surpasses this threshold value, this homogeneous solution is driven unstable by spatially inhomogeneous solutions from the view point of principle of exchange of stability. We use the indices S and H on the shoulder of χ k to indicate that the stability is lost through steady state and Hopf bifurcation respectively as χ crosses χ S k and χ H k . One of the goals of this paper is to establish time periodic patterns to (2.1), in contrast to the stable steady states obtained in [52]. Indeed, the authors in [52] carried out rigorous steady state bifurcation analysis on (2.1) which shows that if χ 0 = min k∈N + χ S k , the stability of (ū,v,w) is lost to stable spatially inhomogeneous steady state of (2.1); moreover weakly nonlinear stability analysis near the bifurcating steady states is also performed which provides a wavemode selection mechanism. On the other hand, numerical simulations in [52] suggest that (2.1) admits stable time periodic solutions when χ is around χ 0 = min k∈N + χ H k . Our approach is based on the Hopf bifurcation theorem for which [2,14] and [32,46] are good references. For example, according to Theorem 1 in [2] or Theorem 1.11 in [14], one of the necessary conditions for χ H k to be a bifurcation value of (2.1) is that the stability matrix (2.3) with χ = χ H k has purely imaginary eigenvalues. We first claim that χ H k = χ S k if Hopf bifurcation occurs at χ = χ H k , ∀k ∈ N + . If not and we assume that χ H k = χ S k , then α 0 (χ, k) = α 1 (χ, k) = 0 and (2.3) has three eigen-values σ 1 (k) = −α 2 (k) < 0, σ 2,3 (k) = 0, under which Hopf bifurcation does not occur. To apply the Hopf bifurcation theorem in [2,14], we also need that, if k = j, Therefore we shall assume the following conditions in the rest of our analysis.
From straightforward calculations, we have that The following results are immediate from straightforward calculations.
Lemma 2.1. Letχ 1 k be given by (2.8). Then for each k ∈ N + , we have that either (i) . According to Lemma 2.1 and our discussions above, Hopf bifurcation may occur at (ū,v,w, χ H k ) only when χ H k < χ S k .
Remark 2.1. In general it is very difficult to determine exactly when case (i) or case (ii) occurs in terms of system parameters. However if the interval length L is sufficiently small, we have that 3) has no purely imaginary eigen-values when L is sufficiently small.

Spatially inhomogeneous periodic patterns
In this section we prove the existence of time periodic spatial patterns of (2.1). To be precise, we want to show that, under proper assumptions on system parameters, the constant equilibrium (ū,v,w) loses its stability through Hopf bifurcation as χ surpass χ 0 = min k∈N + {χ H k , χ S k }. According to our analysis in Section 2, the stability matrix (2.3) has a paired purely imaginary eigenvalues if and only if χ = χ H k < χ S k and there does not exist a time periodic solution to (2.1) that bifurcates from (ū,v,w) if χ H k > χ S k . Therefore, we assume that χ H k < χ S k in the sequel in order to perform bifurcation analysis of (2.1) at χ H k .

Hopf bifurcation
For any χ > 0, we denote the eigen-values of (2.3) by σ 1 (χ, k), σ 2 (χ, k) and σ 3 (χ, k). When In order to apply the bifurcation theory from [2] or [14] at point χ H k , we need to verify the eigenvalue crossing condition or the so-called transversality condition.
Before we state our main result, let us introduce the following Sobolev space Our main result on the existence of nontrivial periodic orbits of (2.1) states as follows.
3) has a paired purely imaginary eigen-values where we have skipped the index k in each eige-value without confusing our reader. It is easy to know that that σ 1 , η and ζ are real analytical functions of χ with η(χ H k ) = 0 and ζ(χ H k ) = ζ 0 > 0. Theorem 3.1 follows from Theorem 1 of [2] once we can prove the following transversality condition In particular we will show that η (χ H k ) > 0, where the prime here and in the sequel means the derivative taken with respect to χ. Substituting the eigen-values σ 1 (χ) and η(χ) ± iζ(χ) into (2.4) and equating the real and imaginary parts there, we have that Differentiating the equations above with respect to χ. Since α 2 is independent of χ, we obtain that 2η (χ) + σ 1 (χ) = 0, (3.6) and Since η(χ H k ) = 0 and σ 1 (χ H k ) = −α 2 (k), solving (3.7) with χ = χ H k gives us that (3.6). This verifies all the necessary conditions required in Theorem 1 of [2], from which our Theorem follows. Theorem 3.1 establishes time periodic spatial patterns to (2.1) bifurcating from (ū,v,w). Moreover it determines the exact bifurcation point χ H k and gives the explicit expression of the oscillation patterns, which admit spatial profile of the eigenfunction cos kπx L . The arguments and results in Theorem 3.1 carry over to multi-dimensional bounded domain Ω.
We have to point out that the condition χ H k < χ S k is necessary for the occurrence of Hopf bifurcation at (ū,v,w, χ H k ). Considering the complexity of both terms, it is very hard to determine or evaluate when χ H k < χ S k , however according to Remark 2.1, for each k ∈ N + , if the interval length L is sufficiently small we must have that χ S k < χ H k . This indicates that Hopf bifurcation does not occur in this situation hence (2.1) does not have time periodic solutions bifurcating from (ū,v,w) when the interval length is sufficiently small. Under this condition, it is showed in [52] that the stability of the homogeneous solution is lost through steady state bifurcation at the first bifurcation branch which has stable stationary solutions of (2.1) with wavemode cos πx L . See Theorem 3.2 in [52]. It is worth mentioning that when µ 1 = µ 2 = 0, we will show in Section 4 that χ S k < χ H k for all k ∈ N + . This implies that χ H k is no longer a Hopf bifurcation point hence there does not exist any time periodic solutions to (2.1) that bifurcate from (ū,v,w, χ H k ).

Stability of time periodic bifurcating solutions
We proceed to analyze stability of the time periodic bifurcating solutions on the bifurcation curves ρ k (s), s ∈ (−δ, δ) obtained in Theorem 3.1. By stability here we mean the formal linearized stability of a periodic solution relative to disturbances from ρ k (s). Assume that all the conditions in Theorem 3.1 are satisfied. Suppose that χ H k0 = min k∈N + χ H k < χ S k , ∀k ∈ N + , then our stability results show that ρ k (s), s ∈ (−δ, δ) is asymptotically stable only if k = k 0 , and ρ k (s), s ∈ (−δ, δ) is always unstable for any k = k 0 . Certainly this is a necessary condition for stability. Moreover a rigorous mathematical treatment of stability away from the equilibrium is so far nonexistent Hopf's pioneering work in 1942 established the basic properties of time periodic solutions to ODEs, such as existence and uniqueness, symmetry properties and stability, etc. Since then a considerable amount of work has been done in studying stability of time periodic solutions to Navier-Stokes equations [27,26,45] or abstract evolutions equations [14,41,30,46]. The stability of time periodic solutions refers to the behavior of the Floquet multiplier or Floquet exponent [20,19,26,27] and we refer the reader to [5,20] or [14] for reviews on the Floquet theory.
Denote u k (s, t) = (u k (s, x, t), v k (s, x, t), w k (s, x, t)) and let (u k (s, t), χ k (s)) be the periodic solutions on the branch ρ k (s) obtained in Theorem 3.1. Rewrite (2.1) into the following abstract form and we skip the index k in (u, v, w) without confusing our reader. Differentiating the abstract system against t, writingu = du dt , we have that 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 −κt , where w is a sufficiently small T -periodic function and κ = κ(s) is a continuous function of s, then we have that dw(s, t) dt = G u (u k , χ k (s))w(s, t) + κ(s)w(s, t), (3.8) where G u is the Fréchet derivative with respect to u. Then stability of the bifurcating solutions in the neighborhood of the branch point χ H k can be determined by computing the eigenvalues of this reduced equation. At s = 0 (3.8) is associated with the eigenvalue problem where It is easy to see that the spectrum of G 0 is infinitely dimensional; in particular, we want to point out that G 0 corresponds to the stability matrix of (ū,v,w) given in (2.3).
, then we must have that the real part of one of these eigen-values must be positive, e.g. Re(σ 2 (χ H k , k)) > 0, since (ū,v,w) is unstable if χ > χ 0 according to Proposition 2.1. Therefore for any positive integer k = k 0 , we have that G 0 (k) must have an eigenvalue with positive real part, hence κ(0) < 0 if k = k 0 . By the standard perturbation theory for an eigenvalue of finite multiplicity (see [33] or [20] e.g.), κ(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 . This implies that if a periodic bifurcating solution is stable, it must be on the k 0 -th branch where χ H k achieves its minimum over N + , i.e., it is on the left-most branch, while all the later branches are always unstable.
We proceed to discuss stability of branch ρ k0 (s) around (ū,v,w, χ H k0 ). According to Lemma 2.10 in [14], the eigenvalue κ(s) is a continuous real function of s which is uniquely defined near s = 0. For χ being around χ H k0 , the eigenvalues of A k0 are σ 1 (χ), σ 2,3 (χ) = η(χ) ± iζ(χ). By Theorem 2.13 in [14], κ(s) and sχ k0 (s) have the same zeros in small neighbourhood of s = 0 in which κ(s) and −η (χ H k0 )sχ k0 (s) have the same sign (if they are not zero), and |κ(s) + η (χ H k0 sχ k0 (s)| ≤ |sχ k0 (s)|o (1) Figure 1. In order to evaluate χ k0 (0), one can follow the calculations using the factorization theorem in [31] or the method of integral averaging in [11], or the normal form method and centre manifold theorem from [19] of Hassard et al.. In each method, we need to perform a perturbation analysis in the neighbourhood of the critical bifurcation value χ H k0 , by substituting χ H k (s) and the periodic solution u(s, x, t) as Taylor series of s into (2.1) then we equate the s-terms to find algebraic equations of χ k0 (0) which determines the direction of the Hopf bifurcation if χ k0 (0) = 0. The calculation is routine but extremely complicated, therefore we skip it here.

System without cellular growth
In this section we study the positive solutions to (2.1) with µ 1 = µ 2 = 0, i.e., the following system One of our main results in this section shows that (4.1) and its multi-dimensional counterpart has no positive time periodic solutions, and this indicates that the cellular growth is responsible for the formation of time periodic spatial patterns of system (1.1) and (2.1). According to our results in Section 3 and [52], (ū,v,w) loses its stability to steady state bifurcating solutions when χ 0 = min k∈N + χ S k and to Hopf bifurcating solutions when χ 0 = min k∈N + χ H k . We shall show that Hopf bifurcation does not occur for (2.1) at (ū,v,w) when µ 1 = µ 2 = 0 from the viewpoint of linearized stability analysis. Then we proceed to investigate the effect of cellular kinetics on the dynamics of (4.1). In particular, we shall show that the kinetics are necessary for the formation of periodic patterns of the two competing species chemotaxis system.

Linearized stability of (ū,v,w)
Similar as in Section 2, our starting point is the linearized stability analysis of the constant solution (ū,v,w). To this end we perform some elementary calculations to show that Hopf bifurcation never takes place at (ū,v,w, χ H k ) when µ 1 = µ 2 = 0, under which the stability matrix (2.3) becomes By the same arguments that lead to Proposition 2.1, we have that (ū,v,w) is unstable with respect to (4.1) if χ ≥ χ 0 = min k∈N + {χ S k , χ H k } and it is locally asymptotically stable if χ < χ 0 where χ S k in (2.5) and χ H k in (2.6) become Similar as above, since we consider positive chemotaxis (chemical being chemo-attractive to the cells), we have that χ 0 > 0 hence χ S k > 0 which implies that To see that there does not exist time periodic positive solution to (4.1) that bifurcating from (ū,v,w), we shall prove that χ H k is not a bifurcation value by showing that χ S k < χ H k for all k ∈ N + .
To this end, we denote χ S k − χ H k = J1+J2 J3 , where and Since J 3 > 0, we just need to evaluate J 1 + J 2 in order to determine the sign of F . Simple calculations give rise to Case 2. If d 1 < d 2 , we can rewrite J 1 + J 2 as follows, which gives rise to where we have used the fact ξ < d2(( kπ L ) 2 +λ) v in (4.3), therefore we have that J 1 + J 2 < 0 in both cases hence χ S k < χ H k for each k ∈ N + as claimed. According to Proposition 2.2 and Lemma 2.1, matrix (4.2) does not have purely imaginary eigenvalues for any χ hence (4.1) has no Hopf bifurcating solutions from (ū,v,w).

Remark 4.1.
If Ω is a bounded domain in R N , N ≥ 2, then the constant solution is unstable if χ > χ 0 , with ( kπ L ) 2 being replaced by the k-th Neumann eigenvalue of −∆. Without loss our generality, we assume that (ū,w,w) is the same as given by (1.2). If not, then thanks to the conservation of cellular populations, we must have that then our calculations above are still true with respect to the change of notations.
According to our discussions above, (ū,v,w) loses stability to steady state bifurcating solutions as χ surpasses χ 0 = min k∈N + {χ H k , χ S k }. However, since χ S k < χ H k , we have from Proposition 2.2 and Lemma 2.1 that the stability matrix (4.2) has no purely imaginary eigen-value. Therefore Hopf bifurcation can not occur for system (4.1), which does not admit stable time periodic patterns bifurcating from (ū,v,w). Moreover, according to the results in [52], we know that the stability of (ū,v,w) is lost through steady state bifurcation to spatially inhomogeneous patterns of (4.1), which has a spatial profile cos k0πx L .

Lyapunov functional
The linearized stability analysis of (ū,v,w) suggests that system (4.1) has no time periodic patterns that bifurcate from the constant solution (ū,v,w) and it does not rule out the existence of time periodic patterns of (4.1) since there may exist other oscillating solutions other than those from Hopf bifurcation. However we shall prove that the latter case is indeed impossible by showing the existence of time-monotone Lyapunov functional to (4.1). Our results hold for (4.1) over multi-dimensional domains and we consider the following fully parabolic system where Ω ⊂ R N , N ≥ 1, ∇ is the gradient operator and ∆ is the Laplace operator. The system parameters are the same in (2.1).
We shall show that (4.4) has a time-monotone Lyapunov functional, therefore it admits no time periodic patterns regardless of the space dimension and system parameters as long as there is no cellular growth. In light of this, another main contribution of this paper is the global existence, boundedness and large-time behavior of positive solutions to (4.4). We begin with the verification that (4.4) has a Lyapunov functional in the following form which is non-increasing along the trajectories of (4.4). We have the following Lemma.
Proof. According to the Maximum Principles (e.g. [38]) and positivity of the initial data, both u are v are stricitly positive onΩ × (0, T ). We have from straightforward calculations that To estimate (4.7), we have from the PDEs and the divergence theorem that while the last two terms of (4.7) becomes and In light of (4.8)-(4.11), (4.7) leads us to Therefore F (u, v, w) is always non-increasing in t and (4.6) follows from (4.12).
The time-monotone Lyapunov functional (4.5) indicates that (4.4) can not have timeperiodic solutions, in contrast to system (2.1) which has oscillating solutions as we have from Theorem 3.1. Therefore the formation of oscillating solutions to (1.1) is contributed by the appearance of cellular growth terms.

Global existence and boundedness for N = 2
In [7], the authors investigated parabolic-parabolic-elliptic system of (4.4) with λ = 0 over the whole space R N , N ≥ 2. Their results state that, in loose terms, if the initial data u 0 and v 0 concentrate at some points x i , i ∈ N + , then the solutions to (4.4) blow up in a finite time. In [13], E. Espejo et al studied the parabolic-parabolic-elliptic system of (4.4) with λ = d 2 = 1 over a unit disk in R 2 under homogeneous Dirichlet boundary conditions. They showed that if then there exist global bounded classical positive solutions. In [12], it is proved that if one of the inequalities fails, then the solutions to a similar problem over R 2 blow up. For Ω = R 2 , global existence and large-time behaviors for (4.4) are investigated in [54] provided that (u 0 , v 0 , ∇w 0 ) L 1 (R 2 ) are sufficiently small, following the arguments on invariant sets of (4.4) as in [53].
This section is devoted to study the global existence and boundedness of classical positive solutions to (4.4) as well as their large-time behaviors. Similar as for (1.1), Amann's theories [3,4] guarantee the local existence of (4.4) since it is a normally parabolic triangle system, while L 1 -boundedness of u and v still holds for (4.4) due to the conservation of cellular populations. By the same arguments for Theorem 2.5 in [52] we can prove the local existence and boundedness for (4.4) over (0, T max ) for some T max ∈ (0, ∞]. We are mainly concerned with the global existence and boundedness of (4.4) over Ω ⊂ R 2 . In particular, assuming we show that the positive classical solutions to (4.4) exist globally and are uniformly bounded as follows.
In Figure 2, we plot the numerical simulations to illustrate the evolution of spatiallyinhomogeneous time periodic patterns of (1.1). The numerics there indicate the lack of a stable global attractor to the full system (1.1), at least for the parameter set we choose. Therefore, we are motivated to investigate large-time behavior of positive solutions to (4.4) by establishing the existence time-monotone Lyapunov functional. We show that the classical solutions to (4.4) converge to its stationary states as time goes to infinity. See Theorem 4.6. For example, the first subgraph in Figure 4 plots the spatial-temporal dynamics of (4.4) over (0, 6), where the interior spike is an attractor of the system.
We now pass to present our proof of Theorem 4.2. Our main vehicle is the L 2 -estimates proven in Lemma 4.5 and the application of standard Moser-Alikakos iteration [2]. To derive the L 2 estimates, we shall estimate energy-type functionals u ln u L 1 and v ln v L 1 via a special version of the Moser-Trudinger inequality. It is necessary to remark that the crucial use of the embedding inequalities only applies when N ≤ 2.
and multiplying this inequality by d1 χ Ω u 0 gives rise to Similarly we have for v(x, t) that On the other hand, we apply (4.14) on (4.16) and (4.17) with φ = (1+δ)χi di w, i = 1, 2, respectively, where χ 1 = χ and χ 2 = ξ, to have that where we have applied the boundedness of w L 1 in (4.18). In light of (4.16)-(4.18), we have that where C 3 is a positive constant that depends on u L 1 + v L 1 . Choosing δ > 0 to be sufficiently small, we see from condition (4.13) that which, together with (4.19), implies that Now we can easily see that (4.20) implies that Ω (u + v)w < C 5 .
On the other hand, we have from straightforward calculations therefore both Ω u ln u and Ω v ln v are bounded. This completes the proof of Lemma 4.4.
The following result is an immediate consequence of (4.15) with p = 3 in Lemma 3.5 of [42].

Corollary 1.
Under the same conditions in Theorem 4.2. Let u, v be the classical solutions to (4.4), then for any > 0, there exists a positive constant C( ) such that Next we provide the boundedness of u(·, t) L 2 + v(·, t) L 2 for t ∈ (0, ∞), which suffices to prove the global existence and boundedness of (u, v, w) to (4.4). Proof. In light of the PDEs in (4.4), straightforward calculations involving the integration by parts and Young's inequality lead us to and To estimate (4.24) (similarly (4.25)), we have from the Gagliardo-Ladyzenskaja-Nirenberg interpolation inequality (see [38] e.g.) and Cauchy-Schwartz that, for any u ∈ W 1,4 (Ω), there exists two positive constants C 1 and C 2 such that Moreover we have from Hölder's that in (4.24) and Moreover we have from Corollary 1 in [10] due to Gagliardo-Ladyzenskaja-Nirenberg inequality that for any > 0 there exist C 13 and C 14 such that where C 13 (C 14 ) is a positive constant that depends on , Ω and u 0 L 1 ( v 0 L 1 ). Now we have from (4.29) that and denoting we have from (4.30) that where d = min{d 1 , d 2 }, C 18 = max{2C 11 , 2C 12 } is a positive constant independent of and C 17 is also a positive constant.
To derive the boundedness of y(t), we recall from (4.6) and (4.15) that t 0 w t (·, t) L 2 ds is bounded for all t ∈ (0, ∞). Solving (4.31) gives rise to where y 0 = u 0 2 L 2 + v 0 2 L 2 and C 19 is a positive constant that may depend on . This finishes the proof of Lemma 4.5.
Proof of Theorem 4.2. The proof is exact the same as proof of Theorem 2.5 in [52] where Moser-Alikakos iteration and standard bootstrap arguments are applied, except that the global existence there is established for N = 1, therefore we shall only sketch the main steps here.
First of all, since (4.4) is a triangular system, its local existence follows from the classical results of Amann [3,4] and the regularity of the solutions follows from standard parabolic regularity arguments. By the same estimates for (2.7) in [52], we can find a constant C > 0 such that + v(·, s) L 2 (Ω) + w(·, s) L 2 (Ω) ds , therefore w(·, t) W 1,q (Ω) is uniformly bounded for each q ∈ (1, ∞) since (u, v)(·, t) L 2 (Ω) are bounded. Then we can again apply Gagliardo-Nirenberg interpolation to show that (u, v)(·, t) L 3 are bounded, which implies that w(·, t) W 1,∞ (Ω) is uniformly bounded. Applying the standard Moser-Alikakos iteration we have the uniform boundedness of (u, v) in L ∞ . This completes the proof.

Asymptotic behaviours and nonconstant positive steady states
In this section, we will show that the bounded classical solutions to (4.4) converge to the (probably nontrivial) steady state as t → ∞. Using the monotonicity of the Lyapunov functional F in (4.5), we can prove that the limit of ω-sets of (4.4) is the stationary system. The following theorem can be proved by the same arguments for Lemma 3.1 in [53] thanks to Lemma 4.5. (4.33) According to Remark 4.1, (4.33) has no nonconstant stable steady state if χ is large. It is interesting to study nonconstant positive solutions to (4.33). In particular the steady states with concentrating properties such as boundary or interior spikes can be used to model the aggregation phenomenon for chemotactic cells. For works in this direction, we refer to the results

Numerical simulations
In this section, we perform some numerical studies of stable and time periodic spatially inhomogeneous solutions to system (2.1). To manifest the effect of cellular growth and other parameters on its spatial-temporal dynamics, we fix a 1 = a 2 = 0.5 in all our simulations for which the IBVP which has positive equilibrium (ū,v,w) = ( 2 3 , 2 3 , 4 3λ ), and the initial data are small perturbations from the equilibrium. Then we choose different sets of system parameters to study the initiation and development of spatial patterns to the system.
First of all, we take d 1 = 5, d 2 = 0.1, µ 1 = µ 2 = 1, λ = 5, ξ = 0.1 and consider (2.1) over domain Ω = (0, 6), subject to initial condition (u 0 , v 0 , w 0 ) = (ū,v,w) + 0.001 cos 2πx. According to our stability analysis in Proposition 1.1, (ū,v,w) is unstable if χ > χ 0 = χ H 2 ≈ 63.2 and according to our bifurcation results in Theorem 3.1 and their stability analysis, the homogeneous solution loses its stability to time periodic pattern which has spatial profile cos 2πx L ; moreover its period is approximately given by T = 2π ζ0 ≈ 8 in (3.2). In Figure 2, we choose χ = 80 and plot (u, v, w)(x, t) for t ∈ (0, 100). The initial data has a spatial inhomogeneity of form cos 2πx, but the periodic patterns develop according to the spatial profile cos πx 3 which is the stable wavemode. Moreover the time period of the oscillating patterns matches our theoretical result. Figure 2: Initiation and development of time periodic spatial patterns to (1.1) over (0, 6) with initial data being small perturbations of (ū,v,w). System parameters are chosen to be d 1 = 5, d 2 = 0.1, µ 1 = µ 2 = 1, λ = 5, ξ = 0.1 and χ = 80. Our theoretical results indicate that the homogeneous equilibrium loses its stability at χ 0 = χ H 2 ≈ 63.2 through Hopf bifurcation to a stable time periodic pattern which has spatial profile cos 2πx 6 and period T ≈ 8. Space and time grid sizes are ∆x = 0.02 and ∆t = 0.05. The numerical simulations are in good agreement with our theoretical findings. Figure 4 is devoted to illustrate the effect of cellular growth rates µ 1 and µ 2 on the pattern formations in (2.1). In particular, we choose µ 1 = µ 2 and plot in each subgraph the spatialtemporal behavior of (2.1).
If the interval length L is small, we always have that χ S k < χ H k , ∀k ∈ N + , then (2.1) does not exist time periodic solutions through Hopf bifurcation. Figure 5 includes a set of simulations on the spatial-temporal behaviors of solutions to (2.1) over different intervals, with the same set of system parameters and initial data. We want to point out that when L = 1, our simulation indicates that (ū,v,w) is an attractor to (2.1) when χ is fixed. Indeed, according to Remark 2.1 we know that χ 0 approaches infinity as L approaches zero, therefore χ has to be sufficiently large to support pattern formation when the domain size is small. Our numerics in Figure 6 are devoted to study the effect of χ on the spatial-temporal dynamics of (2.1) when it is far away from χ k0 = 63.2. In each plot we choose χ = 90, 110, 240 and 300 respectively, but the rest parameters and initial data are the same as those in Figure 2. When χ = 90 or 110, we see that the stable time periodic solutions have the same profiles cos πx 3 as described in Theorem 3.1, though a time periodic spatial pattern with mode cos πx 2 is developed for time t up to 60 when χ = 160. We surmise that this oscillating solution is unstable or meta stable, and a nonlinear analysis is required to determine its stability. Moreover (2.1) has periodic patterns when χ = 240 which is far away from χ 0 ≈ 63.2. We surmise that the existence of this periodic solution is not driven by linearized instability of the homogeneous solution but by the nonlinear cellular growth.  Initial data are taken to be small perturbations of (ū,v,w). Space and time grid sizes are ∆x = L/500 = 0.012 and ∆t = 0.05. We observe that the cellular growth rate µ supports the formation of periodic patterns. However the periodic pattern disappear at µ ≈ 2.1, for which we surmise that the oscillating solutions becomes unstable and develop into stable stationary patterns.

Conclusions and Discussions
In this paper, we study the 3 × 3 system (1.1) modeling the spatial-temporal evolution of two competing species and one-attracting chemical. It has been studied in [48,49] that either the unique positive equilibrium or semi-equilibrium is a global attractor to (1.1) if the chemotaxis Figure 5: The effect of domain size on the pattern formation. We choose the system parameters to be the same as those in Figure 3 except that χ is slighter larger than χ k 0 . ∆x = L/500 and ∆t = 0.05 in each graph. Our simulations support our theoretical findings that large domains support periodic patterns with higher modes, however when the domain size is small, therefore does not exist time periodic solutions that bifurcate from the homogeneous solution. coefficients χ and ξ are small compared to the cellular growth rates µ 1 and µ 2 . Nonconstant positive steady states of (1.1) over Ω = (0, L) have been studied in [52] by rigorous analysis. Our results complement the works in [48,49,52] in studying its positive time periodic spatially inhomogeneous solutions observed by the numerical studies in [52]. Periodic patterns have been experimentally observed in chemotaxis of E. coli or Dictyostelium discoideum by many researchers [23,24,22]. Numerical simulations have been performed to investigate the oscillating patterns by various authors [16,43], however very few works are available on the rigorous mathematical analysis of there periodic solutions (the only related work we know is [39]).
The starting point of our mathematical analysis is the linearized stability of the positive constant solution (ū,v,w) to (1.1) over Ω = (0, L), which becomes unstable if χ > χ 0 = min k∈N + {χ S k , χ H k }. It is proved in [52] that if χ 0 = χ S k1 < χ H k , ∀k ∈ N + , then the stability of (ū,v,w) is lost to nonconstant positive stationary solutions of (2.1) at χ S k1 through steady state bifurcation, while all the rest bifurcating solutions around χ S k are unstable if k = k 1 . The first set of results in this paper state that if χ 0 = χ H k0 < χ H k , χ S k , ∀k ∈ N + , then (ū,v,w) loses its stability to time periodic spatial solutions at χ H k0 through Hopf bifurcation, while all the rest Hopf bifurcation branches must be unstable if k = k 0 . These linearized stability results complete the understanding of the local dynamics of (ū,v,w) which show that (ū,v,w) is driven unstable by large chemotaxis rate through either steady state bifurcation if χ 0 = χ S k1 and through Hopf bifurcation if χ 0 = χ H k0 ; moreover out stability results provide a completes wavemode selection mechanism for system (2.1) in the following sense: the only stable bifurcating solution (through Hopf or steady state) must stay on the left-most branch, while all the rest bifurcating branches are academic in that they are all unstable.
We have also investigated the effect of cellular growth on the dynamics of (1.1) over multidimensional bounded domains. In particular we showed that the cellular kinetics are responsible for the formation of time periodic patterns to (1.1), which has no temporal oscillating patterns when µ 1 = µ 2 = 0. Our proof is based on the construction of time-monotone Lyapunov functional. An extra conclusion we have from the Lyapunov functional is that we proved the global existence and boundedness of classical solutions to (1.1) over Ω ⊂ R 2 provided the total cell population is not too large. Considering system (1.1) over (0, L), we have also studied the effect of the domain size on the pattern formation, which shows that small domain supports stationary patterns while large domain supports time periodic patterns. Numerical simulations are implemented to illustrate both stationary and time periodic patterns to (1.1) which support our theoretical findings.
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 destabilize (ū,v,w). Moreover, 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 [49]. This needs an approach totally different from those in this paper.
When (ū,v,w) loses its stability at χ H k0 , a Hopf bifurcation occurs and the stability is lost to a time-periodic solution. Our analysis of the Hopf bifurcation stability states that it depends on the turning direction of the branch around the equilibrium. It seems necessary to evaluate χ (χ H k0 ) for this sake which need some very hard calculations. The global Hopf bifurcation analysis is also a very interesting problem that worths future exploration. It is interesting and important to ask what happens when at a later stage the time periodic solution loses stability? We refer this to the Poincaré map for which D. H. Sattinger [46] is a good reference.
We showed that (1.1) without cellular growth does not have any time periodic solutions by proving the existence of a time monotone Lyapunov functional. The global existence of (1.1) in 2D depends a prior estimate of this functional, where we have assumed the smallness of initial cell populations. From the viewpoint of mathematical analysis, it is an interesting and important question to study the global existence or blow-ups of (1.1) in higher dimensions when µ 1 = µ 2 = 0. In particular, for the global existence on Keller-Segel chemotaxis models, we refer the reader to the very recent survey [6]. It appears that the positivity of χ and ξ is required in the arguments of [49], and in light of our stability analysis, we surmise that the solution of (1.1) 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 [49].