On a Chemotaxis Model with Saturated Chemotactic Flux

We propose a PDE chemotaxis model, which can be viewed as a regularization of the Patlak-Keller-Segel (PKS) system. Our modification is based on a fundamental physical property of the chemotactic flux—its boundedness. This means that the cell velocity is proportional to the magnitude of the chemoattractant gradient only when the latter is small, while when the chemoattractant gradient tends to infinity the cell velocity saturates. Unlike the original PKS system, the solutions of the modified model do not blow up in either finite or infinite time in any number of spatial dimensions, thus making it possible to use bounded spiky steady states to model cell aggregation. After obtaining local and global existence results, we use the local and global bifurcation theories to show the existence of one-dimensional spiky steady states; we also study the stability of bifurcating steady states. Finally, we numerically verify these analytical results, and then demonstrate that solutions of the two-dimensional model with nonlinear saturated chemotactic flux typically develop very complicated spiky structures.


Introduction
Chemotaxis is a phenomenon of collective movement of microorganisms in the direction of increasing chemical concentration. The simplest and classical PDE model of chemotaxis was introduced in [28] and [20,21]. In this model, which we will refer to as the Patlak-Keller-Segel (PKS) model, the cell density ρ(x, t) and the chemoattractant concentration c(x, t) are governed by the following system of convection-diffusion-reaction equations: ρ t + χ∇ · (ρ ∇c) = ν∆ρ, c t = ∆c − γ c c + γ ρ ρ. (1.1) The most important phenomenon in chemotaxis is the aggregation of cells, namely, the concentration of ρ as t increases. The biological phenomenon behind this is that even when the cells are initially distributed almost evenly over the habitat Ω, later on they, being chemotactic to a chemical released by themselves, start to aggregate in a number of "centers" ( [1, 4-6, 9, 32, 37]). In the literature, two ways of mathematically modeling cell aggregation have been proposed: (i) solutions of (1.1) blow up in finite time and at the blowup time, ρ is a linear combination of several δ-functions, plus a regular part-see [8,14,25]; (ii) time-dependent solutions converge to bounded but spiky steady states-see [24,26,35,36] for such results on several modifications (regularizations) of (1.1). See also the survey papers [18,19].
While the blowup and the formation of the δ-function are not an unreasonable modeling of the cell aggregation phenomenon, they create enormous, and also unnecessary, challenges to numerics and analysis (be it formal or rigorous). Thus we prefer chemotaxis models that only have bounded, global-in-time solutions that approach spiky steady states as time increases. Such models may be obtained by regularizing the PKS system. A variety of regularizations has been proposed over the past decades, see the review papers [17][18][19], the monograph [31] and references therein.
In this paper, we consider a regularization of the PKS model, which is based on a fundamental physical property of the chemotactic flux-its boundedness (this feature is almost always lost in weakly nonlinear, small gradients expansions, underlying the derivation of most continuum models). To derive the modified system we replace the linear chemotactic flux ρ∇c by a nonlinear saturated one, ρQ(∇c), which is proportional to the magnitude of the chemoattractant gradient only when the latter is small and is bounded when the chemoattractant gradient tends to infinity. The regularized model then reads: ρ t + χ∇ · (ρ Q (∇c)) = ν∆ρ, c t = ∆c − γ c c + γ ρ ρ, (1.2) where a smooth saturated chemotactic flux Q(u 1 , . . . , u d ) = Q(u) = (Q 1 (u), . . . , Q d (u)) satisfies the following properties: where C i are constants. Without loss of generality, one may assume that max The synthesized form of the saturated flux is a Pade approximate which connects universal features present at both very small and very large gradients. There is a certain arbitrariness in the choice of the chemotactic flux function Q. A typical example of a saturated chemotactic flux, which is used in all of our numerical experiments, is if |∇c| ≤ s * , |∇c| − s * 1 + (|∇c| − s * ) 2 + s * ∇c |∇c| , otherwise, (1.4) where s * is a switching parameter, which defines a small gradient values, for which the system (1.2) reduces to the original PKS system (1.1) so that the effect of saturated chemotactic flux is felt at large gradient regimes only. Note that when s * = 0, the flux (1.4) becomes a mean curvature type flux: (1.5) The chemotaxis system (1.2) studied here is similar to the regularization (M7) from [17], where a specific form of the bounded function Q was considered (see also [2,27,33]). In [16], a more general type of the cell density equation was considered and an a-priori L ∞ -bound on its solutions was established. The result from [16] also applies to the ρ-equation in (1.2). In §3, we give an alternative proof of the L ∞ bounds on both ρ and c, which also applies to the time independent version of the system (1.2). We then proceed in §4 with the proof of the local existence result for a more general chemotaxis system (4.1)-to the best of our knowledge, this general local existence result is new. This result directly applies to the system (1.2) and together with the a-priori bounds obtained in §3 leads to the global existence of the solution of (1.2). This result is supported by our numerical experiments, in which we compare the blowing up solutions of the PKS system (1.1) with the spiky, but bounded solutions of (1.2), (1.4).
In §5, we study one-dimensional (1-D) steady-state solutions and use both the local [11] and global [30,34] bifurcation theories to investigate existence of nontrivial steady states. We then show that when the chemotaxis coefficient χ is large, the ρ component of the steady state is spiky, while the c component is close to a bounded function given explicitly; we also use [10] to study the stability of bifurcating steady states. Our extensive numerical experiments support the analytical findings. Two-dimensional (2-D) steady states are considered in §6, where formation and stability of the spiky solutions are illustrated numerically.
All numerical results reported in this paper have been obtained using a second-order positivity preserving central-upwind scheme developed in [7]. A brief description of the 2-D version of the scheme is presented in Appendix A.

Several Basic Definitions and Results
In this section, we provide the reader with a briefly review of several basic definitions and results, which will be used afterward.
We begin with Young's inequality, which states that for all positive real numbers a, b, p and q, such that 1/p + 1/q = 1, This inequality also gives rise to the so-called Young's inequality with ε (valid for any ε > 0), We will use the following Sobolev and Hölder spaces:
Let Ω be a bounded domain in R d with a smooth boundary. Then for any u ∈ H 1 (Ω) , and M is a constant depending only on d and Ω.

Corollary 2.1
Let Ω be a bounded domain in R d with a smooth boundary. Then for any u ∈ H 1 (Ω) and any ε > 0

4)
where K ≥ 1 depends only on d and Ω.

A Priori Estimates
Consider the following initial-boundary value problem (IBVP): where Ω is a bounded domain in R d with a smooth boundary ∂Ω and n is the outer normal vector field on ∂Ω.
We first prove that positive solutions of the above IBVP remain bounded for all times.
Theorem 3.1 Let (ρ(x, t), c(x, t)) be a positive classical solution of the IBVP (3.1) with bounded nonnegative initial data. Then, for all x ∈Ω and t ≥ 0,

3)
where C = C(d, Ω) is a constant, which depends on d and Ω only.
The last term in (3.4) is estimated using the inequality (2.4) with u = ρ s 2 and ε such that This results in (3.5) We then fix T ∈ (0, ∞), multiply both sides of the last inequality by the integrating factor exp χ 2 s(s − 1) 2ν t and integrate over the time interval [0, t] for t ∈ [0, T ] to obtain the following estimate: Let us now define the function

7)
Remark 3.2 By using the L ∞ -bounds established in Theorem 3.1, parabolic boundary L pestimates and Schauder estimates (see, e.g., [22]), one can obtain that ρ t , c t and all spatial partial derivatives of ρ and c up to order two are bounded onΩ × [0, ∞).

Existence Results
In this section, we consider the IBVP for a more general chemotaxis system, for which we shall establish a local existence result. This result together with the a-priori estimates proved in §3 will lead to a global existence result for the IBVP (3.1).

Notations, Definitions, Introductory Results
We will use the semigroup theory, for which [13] and [29] are good references. Let X 1 = C(Ω) and X 2 = L p (Ω). We define the linear operator A 1 = −ν∆ + I in X 1 with the domain D(A 1 ) = u ∈ W 2,q (Ω), ∀ q > d ∆u ∈ C(Ω), ∂u ∂n ∂Ω = 0 and the linear operator A 2 = −∆ + I in X 2 with the domain D(A 2 ) = u ∈ W 2,p (Ω) ∂u ∂n ∂Ω = 0 . Then, if the boundary ∂Ω is C 2+ω smooth for some ω > 0, −A i is a generator of an analytic semigroup in X i , denoted by e −A i t , i = 1, 2. For any u 0 ∈ X i , e −A i t u 0 is the solution of the linear parabolic equation u t + A i u = 0 on Ω × (0, ∞) with the homogeneous Neuman boundary condition and the initial condition u(x, 0) = u 0 (x).
Since the smallest eigenvalue of both A 1 and A 2 is 1 and the operators are the generators of the analytic semigroups, we can define their fractional powers A a i for a ≥ 0, with domain D(A a i ) and target space X i . The basic facts are: We shall repeatedly use the following facts: and that for i = 1, 2 where C a , C * a and σ ∈ (0, 1) are constants. We now formally convert the IBVP (4.1) into the following system of integral equations: ∇ · e −A 1 (t−s) (P (ρ(·, s), c(·, s), ∇c(·, s))) ds , c(·, s))) ds. (4.5) Note that for technical reasons, in the first integral of the ρ-equation, we write the product of the operator e −A 1 (t−s) and the divergence operator in the reversed order. These two operators commute only if everything else in that integral is smooth enough; thus at this moment, (4.1) and (4.5) are only formally equivalent. Nonetheless, in this paper we call a solution of (4.5) a weak solution of (4.1).

Local Existence
We shall first prove that (4.5) has a unique local solution (Lemma 4.1), then show that this local solution is smooth and thus it is a classical solution of IBVP (4.1) (Lemmas 4.2 and 4.3).
Proof: For any τ 0 ∈ (0, 1), we define the following space: and for each fixed (ρ 0 , c 0 ) ∈ C(Ω) × D(A a 2 ), we define its subset S by We also define the operator T 1 by with the right-hand side (RHS) equal to the RHS of (4.5).
We shall first verify that there exists a small τ 0 > 0 such that T 1 : S → S . To this end, we establish several bounds.
For any (ρ(·, t), c(·, t)) ∈ S , (4.2) and (4.6) imply that there exists a constant K 0 depending only on ρ 0 C(Ω) and c 0 D(A a 2 ) , such that for any t ∈ [0, Thus, for any t ∈ [0, τ 0 ] and a fixed b ∈ ( 1 2 , 1), we use the bounds (4.3) and (4.7) to obtain and A a 2 e −A 2 (t−s) (c(·, s) + S(ρ(·, s), c(·, s))) L p (Ω) ds where C 1 and C 2 are constants independent of t and τ 0 . Using the fact that lim t→0 + e −A 1 t ρ 0 = ρ 0 in C(Ω) and lim t→0 + e −A 2 t A a 2 c 0 = A a 2 c 0 in L p (Ω), it follows from (4.8) and (4.9) that there exists a small τ 0 (depending only on ρ 0 and c 0 ) such that (ϕ, ψ) ∈ T and In the following we shall prove that for a sufficiently small τ 0 , the operator T 1 is a contracting mapping from S to S . To this end, we first note that for any (ρ 1 , c 1 ), (ρ 2 , c 2 ) ∈ S , we have where K * is a constant depending only on ρ 0 C(Ω) and c 0 D(A a 2 ) . Similarly to (4.8) and (4.9), we obtain that for any t ∈ [0, τ 0 ] with a sufficiently small τ 0 > 0, the following estimate holds: We then use the contracting mapping theorem to conclude that there exists a unique local solution of the system (4.5) in C([0, τ 0 ], C(Ω)) × C([0, τ 0 ], D(A a 2 )).
Next, we treat the ρ-component in a way similar to the treatment of the c-component in the beginning of this proof. Namely, we use [13, Theorem 3.2.2] to infer that ρ is a strong solution (ρ ∈ C 1 ((0, τ 0 ], C(Ω)) ∩ C((0, τ 0 ], D(A 1 ))) of the ρ-equation in (4.1) satisfying the corresponding boundary and initial conditions. We now apply the parabolic Schauder regularity theory to conclude that ρ(x, t) We are now ready to state and prove the following main results on the existence of weak and classical maximal solution of IBVP (4.1). (i) For any given ρ 0 ∈ C(Ω) and c 0 ∈ W 2,p (Ω) with p > d and ∂c 0 ∂n ∂Ω = 0 (or more generally, (ii) Suppose that, in addition to the conditions in (i), ρ 0 ∈ C 1+α (Ω) for some α > 0 and is also a classical solution of the IBVP (4.1) satisfying for some small δ > 0.
Proof: The existence and uniqueness of the maximal solution of (4.5) satisfying (4.18) follows directly from Lemmas 4.1 and 4.2 and from the standard extension argument (see [13,Theorem 3.3.4]); the same extension argument also implies that if T max < ∞, then lim sup In the following, we shall prove (by contradiction) that if (4.21) is true then (4.19) holds. To this end, we assume that (4.21) is satisfied while (4.19) is not, that is, there exists a constant C 1 such that Using (4.22) and the estimates similar to those in (4.9), one can easily obtain that there exists a constant C 2 depending only on C 1 , T max and c 0 D(A a 2 ) such that which, together with (4.22), contradicts (4.21). This completes the proof of (i).
It follows from Lemma 4.3 and its proof that if one can prove the non-existence of a τ 1 ∈ (τ 0 , T max ) such that lim sup Using (4.17), (4.24) and (4.25), we obtain Finally, we apply the general Gronwall inequality ( [13, Lemma 7.1.1]) to (4.26) and conclude that there exists a constant C depending only on K 1 , K 2 , τ 0 and τ 1 such that , which, together with (4.24), contradicts (4.23). This completes the proof of (ii) and of Theorem 4.1.

Global Existence
In this section, we combine the local existence results from §4.2 and the a-priori estimates established in §3 to obtain the existence of the global classical solution of the chemotaxis system with a saturated chemotactic flux. Proof: First, Theorem 4.1 implies that there exists a unique classical solution of the IBVP (3.1) defined on some maximal existence interval [0, T max ). At the same time, Theorem 3.1 ensures that the solution is uniformly bounded there, hence T max = ∞, which means that the classical solution is global. The positivity of the solution follows from the strong maximum principle and Hopf boundary point lemma.

Numerical Examples
We now illustrate the existence of global solutions of (1.2), (1.4) and their stability. In the numerical experiments below, we consider the system (1.2), (1.4) subject to the initial condition    phenomenon inherent in the original PKS model (1.1). Note that the time snapshots of the solution of the PKS system are given at much smaller times than for the system with the saturated flux, since by the last shown time (t = 10 −6 ) the solution has already blown up. To numerically verify this we have performed a mesh refinement study, which clearly shows that, unlike the previous (saturated) case, the height of the spike increases by a factor of 4 when the number of grid cells is doubled. This is consistent with the fact that the magnitude of the "numerical" δ-function is proportional to 1/(∆x∆y).

One-Dimensional Steady-State Solutions
In this section, we consider the 1-D version of the system (1.2) subject to the homogeneous Neumann boundary conditions and seek its steady-state solutions satisfying the following boundary value problem (BVP): where Q is a C 2 -smooth bounded increasing function satisfying (compare with (1.3)) Since in the time-dependent case the total mass of ρ is preserved, the value of is prescribed and, as we will show below, given a m ρ > 0, the system (5.1) has positive solutions for a suitable range of χ.
We shall establish the existence of a solution of (5.1) by using the local bifurcation theory from [11] and the global bifurcation theory for nonlinear Fredholm mappings from [30,34]. To this end, we write (5.1), (5.4) in the abstract form and state several results that are basic for the bifurcation theory to apply.
Lemma 5.1 The operator F , defined in (5.6), satisfies the following properties: (iii) For any fixed (ρ 0 , c 0 ) ∈ X × X , the Frechet derivative is given by Proof: Properties (i)-(iii) can be verified by straightforward calculations, which are left to the reader. Property (iv) can be proved as follows. We rewrite (5.7) as Obviously, T 4 : X × X → Y 0 × Y × R is linear and compact. On the other hand, by Theorem 4.4 from [3] or Remark 2.5 (case 3) from [34], the differential system in T 3 is elliptic and satisfies "Agmon's condition". Then, by Theorem 3.3 and Remark 3.4 from [34], Equipped with Lemma 5.1, we proceed as follows. By property (i), (m ρ , m c , χ) is a constant (trivial) solution of (5.5) for all χ ∈ R. We seek now nontrivial solutions of (5.5) bifurcating from these trivial solutions. The necessary condition for bifurcation to occur at (m ρ , m c , χ) is This BVP will have a nontrivial solution if and only if γ ρ χm ρ Q ′ (0) ν − γ c is one of the Neumann eigenvalues for the interval (0, L), that is, γ ρ χm ρ Q ′ (0) ν − γ c = k 2 π 2 L 2 for some nonnegative integer k. Then, c is the corresponding eigenfunctionc k (x) = cos kπx L . Notice that the case k = 0 can be excluded because L 0 c(x) dx = 0. We have thus obtained that (5.8) is satisfied only for Theorem 5.1 Assume that the function Q satisfies (5.2). Then, for each k ∈ N, there exists an interval (−δ, δ) and continuous functions: and (ρ k (s), c k (s), χ k (s)) is a solution of (5.5) (and (5.1), (5.4)). Moreover, all nontrivial solutions of (5.5) near the bifurcation point (m ρ , m c ,χ k ) are on the curve (ρ k (s), c k (s), χ k (s)).
Proof: Our proof is based on Theorem 1.7 from [11]. In the previous discussion, we have checked all but the following "transvesality condition" required by the theorem: If this condition fails, then the following BVP has a solution (ρ, c). Similarly, as in the discussion for (5.9), we have which we then substitute into the second equation in (5.13) and obtain   Sincec k is a solution of the corresponding homogeneous BVP (5.10), we reach a contradiction (to Fredholm Alternative). This completes the proof of (5.12) and hence that of the theorem.
Next, we wish to extend globally the local bifurcation curves in Theorem 5.1. We shall do so for the first curve (k = 1) by using the global bifurcation theory; for the existence of nonmonotone solutions, we prefer to use the reflection (in x) argument and the monotone solutions (that stay on the first bifurcation branch). The following result guarantees the existence of positive monotone solutions of (5.1) with (5.4).
Finally, we proceed with the proof of (d). By Theorem 4.4 of [34], each of P ± satisfies one of the following three alternatives: (i) It is not compact in X × X × R; (ii) It contains a point (m ρ , m c , χ * ) with χ * =χ 1 ; (iii) It contains a point (m ρ +ρ, m c +c, χ), where 0 = (ρ,c) ∈ Z with Z being a closed complement of N D (ρ,c) F (m ρ , m c ,χ 1 ) = span{(ρ 1 ,c 1 )} in X × X . We take (see (5.11)) If alternative (ii) occurs, then χ * is a bifurcation value, which is impossible (this can be proved as in the proof of (c) above). If alternative (iii) occurs, then we integrate by parts: and reach a contradiction. Thus, according to alternative (i), P ± are not compact in X ×X ×R. This means that they are unbounded in X ×X ×R by the elliptic regularity theory. By Remark 3.1, if χ is bounded, then (ρ, c) is also bounded. Thus, the projection of each of P ± on the χ-axis is unbounded. Since P ± are connected, the projection must be an interval of the form [a, ∞) with a ≤χ 1 . This completes the proof of (d) and thus of Theorem 5.2.

Theorem 5.3
As χ → ∞, ρ concentrates at x = 0, that is, in the sense of distribution, and uniformly on [0, L].
On the other hand, after integrating the second equation in (5.1), we have Therefore, by the Arzela-Ascoli theorem, after passing to a subsequence of χ → ∞, c(x) converges to some function, which we denote by c ∞ (x), uniformly on [0, L]. If we show that c ∞ is given by (5.22), then the convergence in (5.22) holds without passing to a subsequence.
Next, we integrate the first equation in ( If there exists x 0 ∈ (0, L] such that ρ ∞ (x 0 ) = 0, then γρχρQ(c ′ ) Then, by the Sturm oscillation theorem applied to (5.24), c ′ changes sign on (0, x 0 ) for large χ, which contradicts the monotonicity assumption. We have thus shown that ρ ∞ ≡ 0 on (0, 1]. Finally, we integrate (5.23) to obtain Taking the limit in this equation as χ → ∞ and using the Lebesque dominated convergence theorem, we have Hence, This and the fact that

Numerical Examples
The purpose of this section is to illustrate the statement of Theorem 5.3. To this end, we consider the 1-D version of the system (1.2), We shall now show that the bifurcating solutions of (5.1), (5.4) (mentioned in Theorem 5.1) are asymptotically stable provided the coefficient γ c is not too large. Before doing so, we need to determine the direction in which the bifurcation curve turns. By Theorem 1.7 from [11], where Z is defined in (5.20). Furthermore, if Q is C 5 -smooth, then the operator F defined in (5.6) is C 4 -smooth and hence, by Theorem 1.18 of [11], (ρ 1 , c 1 , χ 1 ) is a C 3 -smooth function of s. Thus, we can write the following expansions: where (ψ 1 , ϕ 1 ) ∈ Z , (ψ 2 , ϕ 2 ) ∈ Z , and the o(s 4 ) in the expansions for ρ 1 and c 1 are with respect to the H 2 -norm, and K 2 and K 3 are constants. Since We also make additional assumptions on the chemotactic flux Q (satisfied, for example, by (5.3)): Q ′′ (0) = 0, Q ′′′ (0) < 0, (5.28) so that its Taylor expansion reads and substitute (5.27) into the ρ-equation in (5.1) to obtain Collecting the O(s) terms in (5.29), we obtain On the other hand, recall that (ψ 1 , ϕ 1 ) ∈ Z , which together with the first formula in (5.11) leads to which, in turn, together with (5.33) implies Now multiplying (5.30) by cos πx L , integrating with respect to x and using (5.34), we obtain from which it follows that K 2 = 0. Next, gathering the O(s 2 ) terms in (5.29), we have Note that (5.34) also holds with ψ 1 and ϕ 1 replaced by ψ 2 and ϕ 2 , respectively. This enables us to eliminate ψ 2 and ϕ 2 from (5.35) after the multiplication by cos πx L and integration by parts, so that we obtain Observe that integration by parts yields Multiplying (5.30) by cos 2πx L , integrating with respect to x and taking into account that K 2 = 0, we have Integrating by parts the above integrals and evaluating the first integral on the RHS, yields Multiplying (5.32) by cos 2πx L and integrating by parts, we have This together with (5.39) leads to Finally, combining (5.40) with (5.36)-(5.38), evaluating the second and the last integrals in (5.36), and using the relation (5.11), we obtain Our assumptions on Q, (5.2) and (5.28), guarantee that the last term on the RHS of (5.41) is positive. Therefore, the cubic polynomial in γ c inside the bracket in (5.41) has only one root, which we denote by γ * c = γ * c (L, γ ρ , m ρ , Q ′ (0), Q ′′′ (0)) > 0, and we conclude that sgn(K 3 ) = sgn(γ * c − γ c ) Thus, we have proven the following theorem: Theorem 5.4 Assume that the function Q is C 5 -smooth and satisfies (5.2) and (5.28). Then, the bifurcation curve (studied in Theorem 5.2) at (m ρ , m c ,χ 1 ) turns to the right if γ c ∈ (0, γ * c ) and to the left if γ c ∈ (γ * c , ∞). Also, χ ′ 1 (0) = 0 and sgn(χ ′′ 1 (0)) = sgn(γ * c − γ c ). Two typical bifurcation curves, illustrating the results of Theorem 5.4, are schematically presented in Figure 5.3.
We now study the stability of the bifurcation steady states (ρ 1 (s, x), c 1 (s, x)) for s ∈ (−δ, δ), δ > 0. In order to prove that they are asymptotically stable, we need to show that the real part of any eigenvalue λ of the following eigenvalue problem is positive: x ∈ (0, L), Note that when s = 0, 0 is an eigenvalue and the corresponding eigenspace is N D (ρ,c) F (m ρ , m c ,χ 1 ) = span {(ρ 1 ,c 1 )}. To study small eigenvalues of (5.42) when s = 0, we need to first use the eigenvalue perturbation result of [10,Corollary 1.13]. To fit into the abstract framework of [10], we define K : so that K is linear and bounded, and consider the following eigenvalue problem: Obviously, (5.42) and (5.43) are equivalent (see (5.7)). We first prove that λ = 0 is a "simple eigenvalue" of the pair D (ρ,c) F (m ρ , m c ,χ 1 ), K , which, according to [10], means that is Fredholm with zero index and has onedimensional null space (this has been proved in Lemma 5.1); Proof of (ii): Assume that (ii) is not true. Then, there exists (ρ, c) ∈ X × X such that Integrating the first equation in the above system and using the formula forρ 1 in (5.11), we obtain which contradicts Fredholm Alternative.
To show the stability in the case of γ c ∈ (0, γ * c ), we just need to show that given a small neighborhood of the origin of the complex plane, (5.42) does not have an eigenvalue which has a negative real part and which is outside of this neighborhood. This would follow from the standard eigenvalue perturbation theory if we can show that the limit of (5.42) as s → 0, that is, the eigenvalue problem has no nonzero eigenvalues with nonpositive real parts. To this end, we expand ρ and c as follows: where a 0 = 0 because L 0 ρ(x) dx = 0. Then and λ is an eigenvalue of (5.46) if and only if there exist k ≥ 0 and (a k , b k ) = (0, 0) such that (5.47) holds, which is equivalent to Thus, any nonzero eigenvalue λ must have a positive real part. We surmise that if γ c ∈ (0, γ * c ), then all nontrivial solutions of (5.1), (5.4) on the bifurcation curve P are asymptotically stable, while if γ c ∈ (γ * c , ∞), then the nontrivial solutions of (5.1), (5.4) on P, which are not near the bifurcation point (m ρ , m c ,χ 1 ) are also asymptotically stable, see

Nonmonotone steady states on (0, L)
In this section, we provide a few examples of how to construct multi-spike solutions of (5.1), (5.4) using the monotone solutions, discussed above, and the reflection argument. We also numerically demonstrate that these solutions may emerge as steady-state solutions of a time-dependent IBVP for the system (5.25) and conduct their experimental stability study.
Obviously, we can produce steady states with arbitrary many boundary and interior spikes. The stability of each of these non-monotone steady steady states is the same as the stability of the monotone ones (which we use as building blocks), within the class of functions that have the same symmetry as the non-monotone steady state. However, if the perturbation does not have the same symmetry as the corresponding steady-state solution, the limiting (as t → ∞) solution may have a completely different structure. To illustrate these results numerically, we consider a particular example, in which the system (5.25), (5.3) with different values of the chemosensitivity constant χ and subject to the homogeneous Neumann boundary conditions is numerically studied on the interval [0, 1]. Here, Q is given by (5.3), and ν = γ c = γ ρ = 1.
We first consider the initial data ρ(x, 0) = 1 + sχ 2 cos(2πx), c(x, 0) = 1 + s cos(2πx), withχ 2 = 1 + 4π 2 (see (5.11)) and the small parameter s = ±0.01. We take either the small χ = 1.2χ 2 or the intermediate χ = 12χ 2 . Taking s = 0.01 leads to the convergence towards the numerical steady-state solution, whose ρ-component contains two boundary spikes, see Figure  5.4. Such solution corresponds to an analytical steady state described in Example 1. Switching the sign of s to s = −0.01 leads to the convergence to a different steady state with only one interior spike in the ρ-component of the solution, see Figure 5.5. This corresponds to the single interior spike steady-state solution described in Example 2. We note that in both cases, the obtained non-monotone steady-state solutions are stable under small numerical perturbations present in every numerical computation.  We then consider the following initial data, ρ(x, 0) = 1 + 0.01χ 3 cos(3πx), c(x, 0) = 1 + 0.01 cos(3πx), where ρ(x, 0) has two local maxima and two local minima on the interval [0, 1]. We take a small chemotaxis sensitivity constant χ = 1. the symmetry of the initial datum. The other solution (the one computed with ∆x = 1/200) is presented in Figure 5.7. In this case, the solution preserves its initial symmetry for some time (in fact, if one stops the computation at time t = 1, one may conclude that the numerical solution had already reached its symmetric steady state by that time). However, at later times, the symmetry gets destroyed and the resulting steady-state solution has only one boundary spike at x = 1. The reason why the two numerical solutions are so different is in the fact that 201 is divisible by 3 while 200 is not. Notice that in the latter case, the grid does not correspond to the data symmetry, which leads to the lack of symmetry in numerical perturbations, which in turn causes the break of the solution symmetry for large t.

Multidimensional Steady-State Solutions
In the d-dimensional case, d ≥ 2, we can also obtain a local bifurcation result, which is an analogue of Theorem 5.1. This can be done by slightly modifying the preliminaries and the proof of Theorem 5.1. First, we take X = {u ∈ W 2,p (Ω) | ∂u ∂n | ∂Ω = 0} and Y = L p (Ω) with a fixed p > d. Then, N D (ρ,c) F (m ρ , m c , χ) = {0} if and only if γρχmρQ ′ (0) ν − γ c is one of the positive Neumann eigenvalues, which we denote by µ k , k = 1, 2, . . .. Next, letc k (x) be a Neumann eigenfunction associated with µ k , Then, Theorem 5.1 with cos kπx L replaced byc k (x) holds for any k ≥ 1 such that the Neumann eigenvalue corresponding to µ k is the 1-D one.
In the remaining part of this section, we perform an extensive numerical study of 2-D steadystate solutions. We consider the IBVP (3.1) with ν = γ c = γ ρ = 1, Q given by (1.5), and choose different values of χ and different sets of initial data.
First, we numerically study the case corresponding to the Neumann eigenvalue µ 3 = 5π 2 4 (notice that for smaller eigenvalues the corresponding steady states are quasi 1-D, see (6.1)). We take χ = 10µ 3 and the initial data ρ(x, y, 0) = 1 − 0.01µ 3 cos πx 2 cos(πy), c(x, y, 0) = 1 − 0.01 cos πx 2 cos(πy), (6.2) prescribed in the domain Ω = [0, 2]×[0, 1] (see Figure 6.1). The solution of this IBVP is expected to develop two spikes at two opposite corners of Ω, (2, 0) and (0, 1). Indeed, the numerical steady state computed on the uniform grid with ∆x = ∆y = 1/100, shown in Figure 6.1, meets the expectations.  In the next experiment, we take χ = 100 and consider the IBVP (3.1), (1.5) in the same domain Ω = [0, 2] × [0, 1], but subject to different initial conditions, ρ(x, y, 0) ≡ c(x, y, 0) = 100 + cos πx 2 sin(πy), (6.4) shown in Figure 6.3. In this case, one can also expect the solution to develop spikes at the locations of initial local extrema. However, the numerical solution, computed on the uniform grid with ∆x = ∆y = 1/100, behaves in an unpredictable way: it develops one large interior spike and two smaller spikes at the opposite sides of Ω, and an additional smaller spike at the middle of the third side of Ω, see Figure 6.3. It is instructive to compare the above solution with another numerical solution, computed on a slightly finer uniform grid with ∆x = 2/201 and ∆y = 1/101. The large spike seems to be located at the same interior point (its height is slightly different though), but we now have four additional smaller spikes: two of them are located at the opposite sides of Ω, and two (even smaller) ones emerge at the corners (2, 0) and (2,1). This five-spike steady-state solution is presented in Figure 6.4. The obtained results demonstrate that small numerical perturbations of initial data may lead to quite different solution structures.
The last example suggests that steady states may have quite complicated spiky structures. To further investigate this, we numerically solve the IBVP (3.1), (1.5) with χ = 50 in a larger square domain Ω = [0, 10] × [0, 10] with the initial data ρ(x, y, 0) = 1 + σ, c(x, y, 0) = 0, (6.5) where σ is a random variable uniformly distributed in [−0.1, 0.1]. The numerical solution develops a stable multi-spike structure, see the obtained numerical steady state plotted in Figure 6.5. This result suggests that the studied chemotaxis system with a saturated chemotactic flux can be used to model solutions with multiple spikes (both interior and boundary ones) appearing in real biological systems.
Remark 6.1 We would like to point out that our numerical experiments, in which the multispike solutions arise, provide only preliminary results. A further rigorous analysis of the 2-D chemotaxis system with a saturated chemotactic flux is required to fully understand the development, evolution and stability of such complicated solutions. Figure 6.5: Multi-spike structure emerging out of the random initial data (6.5).

Remark A.3
The semi-discrete scheme (A.1) is a system of time-dependent ODEs, which has to be integrated numerically using a stable and accurate ODE solver. In this paper, we have used a third-order strong stability preserving (SSP) Runge-Kutta method from [12]. The efficiency of the fully discrete method can be improved by applying an SSP implicit-explicit Runge-Kutta method (see, e.g., [15] and references therein), as discussed in [7].
Remark A.4 The 1-D version of the numerical method used in the paper can be easily reduced from the 2-D scheme described above.