Boundary spikes of a Keller-Segel chemotaxis system with saturated logarithmic sensitivity

In this paper, we study the nonconstant positive steady states of a Keller-Segel chemotaxis system over a bounded domain $\Omega\subset \mathbb{R}^N$, $N\geq 1$. The sensitivity function is chosen to be $\phi(v)=\ln (v+c)$ where $c$ is a positive constant. For the chemical diffusion rate being small, we construct positive solutions with a boundary spike supported on a platform. Moreover, this spike approaches the most curved part of the boundary of the domain as the chemical diffusion rate shrinks to zero. We also conduct extensive numerical simulations to illustrate the formation of stable boundary and interior spikes of the system. These spiky solutions can be used to model the self--organized cell aggregation phenomenon in chemotaxis.


Introduction and main result.
Chemotaxis is the oriented movement of cells along the gradient of certain chemicals in their environment. One of the most interesting phenomena in chemotaxis is the cell aggregation and it is an important mechanism for the formation of a fruiting body from single cells. See [7,14] and [15] for detailed discussions. This phenomenon has been discovered in bacteria, in particular E. coli, and also in the cell slime molds such as Dictyostelium discoideum. See [2,7] and the references therein. One of the main reasons of the directed movements of the cells is that they can navigate within a complex environment by detecting and processing of a variety of internal and external signals. Certain chemicals, most frequently inorganic salt, amino acids and some proteins called chemokines, can induce chemotaxis in motile cells.
In this paper we study the existence of nonconstant positive solutions (u, v) = (u(x), v(x)) to the following system    ∇ · (∇u − pu∇ ln(v + c)) = 0, x ∈ Ω, 2 ∆v − v + βu/α = 0, x ∈ Ω, ∂u ∂n = ∂v ∂n = 0, where Ω is a bounded domain in R N , N ≥ 1, with smooth boundary ∂Ω. n is the unit outer normal on the boundary. The parameters p, c, , α and β are assumed to be positive constants. Let u be the cellular population density and v be the chemical concentration, then system (1) can be used to model the steady state of a chemotaxis system with saturated logarithmic sensitivity function. We 1232 QI WANG are concerned with boundary spike solutions to (1) and our main theorem goes as follows.
The mathematical modeling of chemotaxis was initiated by Keller and Segel in their pioneering works [16,17,18] during 1970s. Denote by u(x, t) the cell population density and by v(x, t) the chemical concentration at space-time (x, t). Then the general Keller-Segel model reads as follows where d 1 , d 2 > 0 are the diffusion rates of cells and chemicals respectively and χ(u, v) is a function that measures the chemotactic response of cells to the chemical. Moreover, χ is positive if the chemical is chemo-attractive (sugar, nutrition, etc.) to the cells and it is negative if the chemical is chemo-repulsive (poison, toxin for example). φ(v) > 0 is the so-called (signal-dependent) sensitivity function and it reflects the variance of cellular sensitivity with respect to the chemical concentration. k(u, v) is the creation and degradation rate of the chemical. The non-flux boundary conditions mean that the domain Ω is enclosed thus it inhibits cell immigration and chemical diffusion across the boundary. The nonnegative initial data u 0 , v 0 are not identically zero. Keller-Segel system (3) and its variants have been studied extensively over the past few decades. These systems are very rich in mathematical dynamics and they can induce many interesting and striking properties such as global solutions, finite time blow-ups, traveling waves, etc. Furthermore, the intuitively simple system (3), even in its simplest form, successfully demonstrates its ability in presenting solutions with self-organized spatial patterns with concentration properties such as boundary spike or interior spike. These spiky solutions can be used to model the aforementioned cellular aggregation and concentration phenomena. See the survey paper [12] for detailed discussions.
From the viewpoint of mathematical modeling, we say that cells aggregate if there exists a solution u(x, t) that converges to a nonconstant positive function u(x), In this paper, we focus on the interplay between the sensitivity function φ(v) and the terms d 1 , d 2 and χ on the mathematical modeling of cell aggregations phenomenon. To make our goal evident, we assume that d 1 and d 2 are positive constants and k(u, v) = −αv + βu for some positive constants α and β. Under these settings, Keller-Segel model (3) becomes and we want to find nonconstant positive steady states to systems in the form of (4). Time-dependent system (4) can model cell aggregation if L ∞ -norm of the solutions goes to infinity, then the aggregation is simulated by a δ-function or a linear combination of δ-functions that represents the chemical concentration. The formal study of blow-up solutions of Keller-Segel models was initiated in [5,25] and the local profile of the blow-ups was rigourously constructed in [11] over a 2D domain. Though such blow-up solutions are evidently connected to the phenomenon of cellular aggregation, a δ-function is not the optimal choice from the viewpoint of mathematical modeling. It challenges the rationality that chemical concentration can not be infinity at a spot in the domain and it also brings difficulties to numerical simulations. Moreover, rigourous analysis of the states after blow-up is impossible. An alternative approach available in the literature is to show that positive solutions of (4) exist globally in time and converge to bounded steady states as time increases. Then the steady states with concentrating properties such as spikes, transition layers, etc. can be used to describe the cell aggregation phenomenon. We shall take the second approach in our paper and it is the goal of the current work to investigate the existence of nonconstant positive steady states of system (4).
For N = 1, it has been well known that solutions to (4) and a collection of its variants are global and bounded in time. See [13,30] for example. However, the global existence is not always available for higher dimensional domains. It turns out that the sensitivity function φ(v) and system parameters in (4) play essential roles in existence of its global-in-time solutions. Two most commonly used sensitivity functions are φ(v) = v and φ = ln v, which lead (4) to the so-called minimal model and logarithmic model respectively. See the survey paper of Hillen and Painter [12].
For φ(v) = ln v, Nagai, et al. [24] established the global solutions of (4) with d 1 = d 2 = α = β = 1 over Ω ⊂ R 2 provided that either χ < 1 or χ < 5/2, Ω is a disk and (u 0 (x), v 0 (x)) is radial in x. Winkler [40] studied the same problem and established the global classical solutions over Ω ⊂ R N , N ≥ 2 if χ < 2/N and global weak solutions if χ < (N + 2)/(3N − 4). Furthermore, Nagai and Senba [22] proved the global existence and boundedness for the parabolic-elliptic system of (4) when Ω ⊂ R 2 is a disk and (u 0 (x), v 0 (x)) is radial, or χ < 2, Ω ⊂ R 3 is a ball and (u 0 (x), v 0 (x)) is radial. An exploration of literatures suggests that both the smallness of initial data and of the chemo-attraction rate χ tend to prevent solutions of (4) from blowing up in finite or infinite time.
Attempts have been made to show that the steady states of (4) have concentrating properties such as boundary spikes, interior spikes, etc., when the chemical diffusion rate is sufficiently small. This approach was initiated by Lin, Ni and Takagi in a series of remarkable papers [20,28,29] around 1990s. To be precise, their pioneering works focus on the logarithmic case φ(v) = ln v and can be summarized as follows. Through the transformation introduced in [31] and proper parameter scalings, the stationary system of (4) reduces to the following single equation where p = χ d1 . Equation (5) is associated with an energy functional with v + = max{v, 0}. Suppose that p ∈ (1, ∞) if N = 1, 2 and p ∈ (1, N +2 N −2 ) if N ≥ 3, Lin, Ni and Takagi [20] applied variational techniques and min-max arguments to prove that if d 2 is small, (5) has only nonconstant positive least-energy solution, which is a critical point of J d2 in H 1 (Ω). Furthermore, Ni and Takagi [28,29] showed that the least-energy solutions must concentrate at a boundary spike if d 2 is sufficiently small and the spike moves to the place of ∂Ω with the largest mean curvature as d 2 goes to zero. Motivated by their works, Grossi, Pistoia and Wei [9] obtained solutions of (5) with a single interior spike that stays in the most centered region of the domain (in the sense of a critical point of the distance function) for small d 2 . By the Lyapunov-Schmidt method, Gui and Wei [10] construct multi-peak solutions with m interior spikes and n boundary spikes, for any nonnegative integers m, n, as long as d 2 is sufficiently small. The methods and techniques used in these brilliant works began with good guesses on the structures of the spike(s). Moreover, the stability analysis has also been carried out through the infinite dimension Lyapunov-Schmidt reduction and the method of NLEP (nonlocal eigenvalue problem) recently developed by J. Wei, et al. See [32] and the surveys in [26,27] for recent development in this direction.
It is worthwhile to point out that a large chemotactic coefficient χ can also destabilize the constant solution of (4) like a small diffusion rate. Therefore, nontrivial steady states may emerge through large chemotactic coefficient χ. For example, by taking χ as the bifurcation parameter, Wang, et al. [4,36,37] applied global bifurcation methods to investigate nontrivial steady states of (4) over one-dimensional finite interval. Furthermore, these solutions are shown to admit boundary spikes or interior transition layers as χ approaches infinity.
The logarithmic sensitivity function φ(v) = ln v was chosen largely due to the Weber-Fechner's law describing cellular behaviors: the subjective sensation is proportional to the logarithm of the stimulus intensity [12]. For the logarithmic model (4) with φ(v) = ln v, we observe that the dynamics of cellular movements are dominated by the taxis flux χ v ∇v which may be unbounded for v ≈ 0. We want to note that this singularity is an important mechanism in the formation of traveling wave solutions as showed in [16] or discussed in the survey paper [38]. See [21,34] for more works on the logarithmic model. However, from the viewpoint of mathematical modeling, it is not reasonable to assume that low chemical concentration elicits significant chemotactic responses from the cell's motility, hence we modify this problem by taking φ(v) = ln(v + c) for a positive constant c which saturates φ (v) at v = 0. Then the steady state of (4) becomes where d 1 , d 2 , χ, c, α and β are positive constants. First of all, we claim that the total cell population in the time-dependent counterpart of (6) is conserved. Indeed, integrating the first equation over Ω, we have from the Neumann boundary conditions that and this implies that for all t > 0, Throughout this paper, we assume the cell population in (7) is a fixed positive constant denoted by Ω u(x)dx = M . Moreover, since small chemical diffusion rate supports nontrivial positive solutions to (6) as we shall show, we denote for the simplicity of notations that then it is easy to see that (6) leads us to (1). According to our results in Theorem 1.1, v (x) concentrates at the boundary point P if the chemical diffusion rate is sufficiently small. Compared to the spiky solutions in [20,26,27] and [32], v (x) has a boundary spike which is supported on a platform with height t * in (2). The results in [20,26,27] and [32] are based on the application of the classical variational methods due to Ambrosetti and Rabinowitz [1] on the following equation, equipped with the energy functional

QI WANG
It is assumed, among other conditions, that f (t) ≡ 0 for t ≤ 0 and f (t) = o(t) as t → 0 + . As we shall see later, system (1) can be reduced to the single equation (10), however, f (v) there apparently does not satisfy these conditions, therefore it inhibits the direct application of the results in [20,26,27]. On the other hand, one needs some nontrivial modifications of the variational methods in [1] in order to apply them for (10). Furthermore, the nonlocal structure in (10) can not be scaled out as in [20] and this makes the application of previous results on (10) difficult. Therefore the adaption of the ideas in [20], etc. to (1) requires several technical tricks. Compared to the results obtained in [28,29] and [32], where the sensitivity function is φ(v) = ln v, our boundary spike is supported on a platform with height t * > 0. If c = 0, we see from (2) that our results coincide with those obtained in [28,29], except that the assumption on the initial cell population becomes M ≤ 0, which apparently can not be true. However, we can consider an approximation problem with c ≈ 0, then we can choose the parameters properly, for example α being large, such that both the assumption on the initial cell population M and the smallness of can be achieved.
We want to point out that the assumptions on M and p are made out of technical reasons. As we have discussed earlier, either small initial total cell mass or small chemo-attraction coefficient χ = pd 1 tends to prevent solutions of (4) from blowing up in finite or infinite time, therefore it makes our analysis of pattern formations in the steady state realistic. It is also worthwhile to mention that the global existence of a parabolic-elliptic type system (4) with φ(v) = ln(v + c) was studied in [3], and the fully parabolic system has been investigated in [35] for d 1 = d 2 = α provided that χ/d 1 = p < 2/N , N ≥ 2. However, the full understanding about its global existence or finite time blow-ups is away from complete and there are many problems widely left open. See the surveys [12,14,15]. Though our assumption on p = χ/d 1 in Theorem 1.1 does not guarantee the global existence of the timedependent systems, we surmise that the smallness assumption of the cell population might contribute in this matter.
The rest of this paper is organized as follows. In Section 2, we convert system (1) into a single equation with a nonlocal structure. Then we establish the existence of nontrivial positive solutions of the corresponding local problem. Section 3 is denoted to study the asymptotic behaviors of the positive solutions obtained in Section 2 as → 0. The proof of Theorem 1.1 is included at the end of this section. Finally, we present numerical simulations in Section 4 to illustrate and verify our numerical findings. Throughout the rest of this paper, we assume that C and C i , i = 1, 2, ... are generic positive constants that may vary from line to line.

2.
Existence of nonconstant positive solutions. First of all, we observe that system (1) can be reduced to a single equation with an integral constraint. Indeed, we see that u-equation in (1) is equivalent to Testing it by ln u − p ln(v + c) over Ω by parts, we have that therefore, there exists a positive constant C such that Integrating (8) over Ω leads us to where M = Ω u(x)dx is a fixed positive constant representing the total cell population. Denoting m = βM/α, we readily see that v(x) satisfies the following nonlocal equation, which serves as a prototype of (1). Indeed, if v(x) is a positive solution to (10), then it solves (1) together with u(x) obtained through (9). Two main difficulties emerge in proving the existence of positive nonconstant solutions of (10). First of all, the nonlocal structure Ω (v + c) p dx inhibits the direct applications of some standard arguments such as degree-theory etc. If c = 0, the equation (10) can be transformed into the following problem without a nonlocal term 2 ∆w − w + w p = 0, through the scaling v = mw( Ω w p ) −1 as in [28], however, such linear scaling is not available to eliminate the nonlocal structure in (10). Moreover, it is well known that (10) is naturally associated with an energy functional and this motivates one to apply variational methods to obtain nontrivial solutions, e.g, the Mountain Pass Theorem, [1,6], etc. However, the presence of c in (10) requires nontrivial modification of these variational methods to this end. Motivated by the ideas initiated by Ni and Takagi in [28,29] and developed by Sleeman et al. in [32], we first investigate the solutions of (10) with Ω (v + c) p dx replaced by a positive constant δ, while we have to treat the dependence of v(x) on δ, then we proceed to find a solution v(x) such that the integral condition Ω (v + c) p dx = δ is satisfied. To be precise, we first consider the following problem with δ > 0 being an arbitrary but fixed constant. The following results indicate that δ can not be arbitrarily small if (11) has a positive solution.

Proof.
Let v(x) be a positive solution of (11), then we integrate the first equation over Ω to obtain that Since R δ (0) > 0, we must have that R δ (t) is negative at its minimum point over (0, ∞). In particular, we see that and R δ (t) achieves its minimum at the critical point Therefore we have that R δ (t * ) < 0 or equivalently −t * + m(t * +c) p δ < 0. On the other hand, we observe that (t * + c) p−1 = δ mp , then the inequality above implies that After substituting the formula of t * into t * > c p−1 , we readily have the desired lower bound δ 0 in (12). To establish the upper estimate on t 2 , we see that then t 2 < δ m 1 p−1 follows from the inequality above. To show the monotonicity of t 1 (δ) and t 2 (δ) in δ, we note that 0 < t 1 (δ) < t * (δ) < t 2 (δ). Differentiating R δ (t) = 0 with respect to δ gives rise to then this leads us to the desired monotonicity instantaneously thanks to the fact that t 1 (δ) < t * < t 2 (δ). For all δ ≥ δ 0 , we see that t 1 (δ) ≤ t 1 (δ 0 ) and t 1 (δ) is uniformly bounded in δ. Finally, we have The proof of this proposition completes.
We are interested in positive solutions to (11), hence we shall assume δ ≥ δ 0 from now on. In particular, (11) allows only constant positive solutionv ≡ c p−1 if δ = δ 0 . For the sake of simplicity, we skip the index δ in t 1 and t 2 and denote them as the first and the second root of R δ (t) correspondingly, unless otherwise noticed.
To study positive solutions to (11), we introduce the transformation where t 1 is obtained in Proposition 1, then (11) becomes where c δ = 1 − mp δ t 1 + c p−1 ∈ (0, 1), and the nonlinear term reads To obtain nonconstant positive solutions w ,δ (x) of (15), we want to take the variational approach, i.e., to find nontrivial critical points of an energy functional associated with (15) in certain Hilbert space. To this end, we endow Sobolev space and assign (15) the energy functional where It is well known that any critical point of J ,δ is a weak solution of (15) in H 1 (Ω); moreover, since p is subcritical, the weak solution is classical by the standard elliptic regularity theories involving the bootstrap argument.
To obtain nonconstant critical points of J ,δ , we shall apply the well-known Mountain Pass Theorem on the energy functional J ,δ (17) due to Ambrosetti and Rabinowitz in [1] and Ni, et al. in [6,28] which states that, and e(x) ≡ 0 is non-negative with J ,δ (e) = 0. Moreover, it is known that C ,δ is independent of the choice of e. The critical point w ,δ (x) of (17) in H 1 (Ω) is called a least-energy solution of (15) and J ,δ (w ,δ ) = C ,δ is called the least-energy value.
To apply the Mountain Pass Theorem, we first observe that J ,δ ∈ C 1 (H 1 (Ω); R), J ,δ (0) = 0 and J ,δ satisfies the Palais-Smale condition thanks to (19) and the proof of Lemma 3.6 in [1]; moreover, the argument that leads to Lemma 3.1 in [1] together with (19) yields that, there exist some positive constants r and C such that J ,δ (w ,δ ) > 0 if 0 < w ,δ < r and J ,δ (w ,δ ) ≥ C > 0 if w ,δ = r. We now need to show, for > 0 being sufficiently small, that there exist a nonnegative function e(x) ∈ H 1 (Ω) and a positive constant t 0 such that J ,δ (t 0 e) = 0. To this end, we choose a test function following (2.10) in [20].
Without loss of generality, we assume that Ω contains the origin (otherwise, we can shift Ω without changing the shape of this test function) and define

QI WANG
First of all, we have the following lemma which estimates the energy of this test function.
Lemma 2.1. There exist two positive constants C 1 and C 2 with C 1 < C 2 which are independent of and δ such that Proof. We apply the arguments in Lemma 2.4 in [20]. To this end, we shall need to verify that f δ (t) in (16) satisfies (h3 ) in [20]: and there exists some positive constants a 1 and a 2 such that with p ∈ (1, ∞) if N = 1, 2 and p ∈ (1, (N + 2)/(N − 2)) if N ≥ 3. In particular, we need to show that a 1 and a 2 are uniform in δ.
In order to prove (22), we first see that the second inequality holds thanks to (17) and the definition of F δ (t). To show the first inequality, we test (15) by w ,δ and then integrate it over Ω by parts In light of the energy functional defined in (17), we see that therefore we obtain from (23) that and this concludes the proof of Lemma 2.2.
Now we are ready to present the following existence and nonexistence of nonconstant positive solutions to (15).
Proof. Since f δ (t) satisfies (19) and p is subcritical, similar as the proof of Theorem 2 in [20], we can show that the least-energy value C ,δ is achieved at w (x), which is a H 1 solution to (15). Moreover, it follows from the standard elliptic regularity argument that w (x) is a classical solution. On the other hand, similar as the analysis of Section 2 in [20], we can show that C ,δ ≥ C 1 N and w (x) is uniformly bounded in both and δ. Moreover, nonconstant solution w (x) must be strictly positive overΩ according to the strong maximum principle and Hopf's boundary lemma. Furthermore, C ,δ ≤ C 2 N follows from Lemma 2.1 and the fact that C ,δ ≤ sup t≥0 J ,δ (te) .
The proof of (ii) is classical and we will leave it to the reader. To prove (i), we shall show that J −1 ,δ (C ,δ ), which may consist of constant solutions of (15), admits only nonconstant positive solutions if is small independent of δ. To this end, we first observe thatw ≡ (t 2 − t 1 )δ − 1 p−1 is the unique positive solution to (15). To rule out it as a least-energy solution of (15) for small , we have from (22) that and if is small J ,δ (w) > J (w ,δ ) = C N . Thereforew = t 2 − t 1 can not be a least-energy solution and this finishes the proof of Proposition 2.
Corollary 1. Let w (x) be a solution of (15), then we have that where C is a positive constant independent of and δ.
Proof. This follows from (23) and (25). 3. Single boundary spike on a platform. In this section, we construct a positive solution v (x) of (1) that has single boundary spike supported on a platform for small. To this end, we first introduce the equation in the entire space, which we shall use to approximate least-energy solution w ,δ (x) of (15).
Lemma 3.1. For each δ ∈ (δ 0 , ∞), there exists a unique solution w δ to the following problem where c δ and f δ are the same as in (15). Moreover, w δ satisfies the followings: (i) The solution w δ is radially symmetric such that, w δ (x) = w δ (|x|); (ii) w δ < 0 for r > 0, with r = |x|; (iii) w δ (r) ≤ Ce −µr , r > 1, where C and µ are positive constants independent of δ. Furthermore, the system is associated with an energy functional I δ (w), which is positive and uniformly bounded in δ.
Proof. First of all, we see that for all δ ∈ (δ 0 , ∞) then by the celebrated Gidas-Ni-Nirenberg [8] symmetry theorem and the uniqueness result from Kwong and Zhang [19], there exists a unique solution w δ (x) to (27) which is radially symmetric. Choosing = 1 in the test function e(x) defined in (18), we can show that I δ (w δ ) is uniformly bounded in δ by the similar proof as for (25), moreover there exists a positive constant C independent of δ such that Then we can apply the radial lemma of Strauss [33] to obtain that where C is a positive constant which is independent of δ. After applying Proposition 4.1 in [8] to (27), we have that and this concludes the proof of the lemma.

QI WANG
To prove (29), we first see that it is equivalent to find a δ such that Indeed, we integrate the v-equation in (10) over Ω and collect that which shows that (29) is equivalent to Ω v ,δ dx = m as claimed. We now present the proof of our main result.
Proof of Theorem 1.1. Let w ,δ (x) be a least-energy solution obtained in Proposition 3. Following [32], we denote the set of least-energy solutions to (15) as It is easy to see that S δ is nonempty; moreover S δ is compact since where C is uniform in and δ. Taking we have from (14) that where t 1 is given in Proposition 2. By the compactness of S δ , inf w ,δ ∈S δ Ω w ,δ (x)dx is a well-defined and continuous function of δ, hence ρ(δ) is a continuous function of δ. Furthermore we have from Proposition 3 that for any w ,δ ∈ S δ , where w δ is the unique ground state of (27) and o(1) is independent of δ. R N w δ (x) dz is uniformly bounded in δ. According to Proposition 2, we have that where the inequality is due to the assumptions on M and c. On the other hand, since t 1 (δ) → 0 as δ → ∞ uniformly in , we can find δ 1 large such that then it follows from (14) and (31) that Therefore, we can take small enough independent of δ such that Together with (32) and (34), we observe from the Intermediate value theorem that there exists δ ∈ (δ 0 , δ 1 ) such that Now by taking δ = δ , we see that Proposition 3 implies all but part (ii) and (2) of To prove that the maximum of v (x) has a positive lower bound, we have that where the last identity follows from the integral constraint To show (2), we put t * = t 1 (δ ) and conclude that t * = t 1 (δ ) < t 1 (δ 0 ) = c p−1 , where the inequality and the last identity follows from Proposition 2 and its remark. On the other hand, we integrate v (x) over Ω and collect that where the last identity follows from (31) or Proposition 3 with δ = δ . Sending to zero in (35), we can readily conclude that t * approaches to βM α|Ω| as claimed. Thus we have proved (2) and this completes the proof of Theorem 1.1.

4.
Numerical simulations of spiky solutions. In this section, we present numerical results on the the formation and evolution of boundary spike of (4) over Ω = (0, 1) × (0, 1) to illustrate our theoretical results. Putting α = β = 1, and φ(v) = ln(v + c), we perform extensive numerical simulations on (4) by choosing different sets of values for , c and the initial data. We shall see that multi-spike solutions also arise through the system. It is worthwhile to mention that rigorous analysis is needed to fully understand the dynamics such as the large time behaviors and stability of these structures, which is beyond the scope of this paper.
For c = 0.1, numerical simulations of system (4) on the unit square are plotted in Figure 1. The graphes in the first line represent the spatial-temporal behaviors of the cellular population density u(x, y, t) and the graphes in the second line illustrate the behaviors of the chemical concentration v(x, y, t). Our simulations show that u quickly develops an interior spike which gradually moves to the corner. Eventually a stable spike is formed at the boundary point corner (0, 0), where the boundary is most curved. In Figure 2, we choose c = 5 and take the same initial data and system parameters as in Figure 1. For this set of values, solutions quickly evolves into a boundary spike at t = 100 compared to that in Figure 1. It eventually develops into a stable boundary spike concerned at (0, 0). We find through our numerical results that the stable boundary spike stays on a platform with height 0.02, compared with those obtained in Figure 1. According to our steady state analysis, we surmise that this platform is contributed by the large magnitude of c. We need to point out that our numerical results are ambiguous since the height of the platform for Figure 2 is t * = βM α|Ω| = 3, which is greatly larger than 0.02. Figure 1. The formation of boundary spike of system (4) at corner (0, 0) of the unit square in 2-D, where initial data are taken to be u 0 (x, y) = 3 − cos(πx) cos(πy), v 0 (x, y) = 3 + cos(π(x − 1/4)) cos(π(y − 1/4)) + cos(π(x − 1/2)) cos(π(y − 1/2)). The parameters are d 1 = 1, d 2 = 0.01, χ = 3, α = β = 1 and c = 0.1.
However, the numerical simulations support our analysis qualitatively and rigorous stability analysis is needed to fully understand the dynamics of the spiky solutions.
We also perform numerical simulations in Figure 3 which show that the constant solution is a global attractor of (4) if c is sufficiently large, independent of the initial data. This corresponds to the fact that the saturating coefficient c in (1) reduces the chemo-attraction magnitude and stabilizes the positive constant steady state. Figure 2. The initial data (u 0 , v 0 ) and parameters are chosen to be the same as in Figure 1 except that c = 5. We see that a stable single boundary spike develops at (0, 0).
We proceed to simulate the formations of stable interior spike and multi-spikes in Figure 4. Our goal is to illustrate the effect of initial data on the localization of the spikes, i.e., the spot where a spike is formed. To this end, we choose different initial data in Figure 4. The initial data on the left are taken to be u 0 (x, y) = 3 − cos(πx) cos(πy)−cos(π(x−1)) cos(π(y −1)), v 0 (x, y) = 3+cos(π(x−1/2)) cos(π(y − 1/2)), where the initial chemical concentration v 0 attains its maximum at point  Figure 1 except that c = 10. Initial data are u 0 = 3 + cos 2πx cos 2πy and v 0 = 3 − cos 2πx cos 2πy. We see that solutions converge to the global attractor (3,3).
(1/2, 1/2). We see that a stable interior spike is formed at this point. The initial data on the right are u 0 (x, y) = 3−cos(πx) cos(πy)−cos(π(x−1)) cos(π(y−1)), v 0 (x, y) = 3 + cos(πx) cos(πy) + cos(π(x − 1)) cos(π(y − 1)) and v 0 achieves its maximum at two boundary points (0, 0) and (1, 1). Then a stable double boundary spike is developed at both maximum points. Our numerical simulations indicate that the spike of (4) with saturated logarithmic sensitivity can localize at the spot where the initial chemical concentration reaches its largest value, independent of the initial cellular population distribution. This complies with the well-accepted theoretical results that chemotaxis dominates the dynamics of the cellular movements over the domain when χ is not too small or c is not too large. However, a rigorous analysis for this purpose requires totally non-trivial mathematical techniques. Our results indicate that spikes usually develop at the spots where initial chemical concentration maximizes.

QI WANG
Finally, we present Figure 5 to illustrate the emergence of multi-spikes in model (4) with parameters taken to be d 1 = 1, d 2 = 0.001, χ = 3 and α = β = c = 1. The initial data are u 0 (x, y) = 3 + 0.1 cos(πx) cos(πy) and v 0 (x, y) = 3+cos(2πx) cos(2πy). We observe that both u and v develop multi-spikes (t = 100) at the four corners and the center very quickly. The structures of these multi-spikes keep well-preserved for quite a long time period (through t = 100 to t = 300), then four spikes at the corner disappear and both u and v develop into a stable interior spike at the center eventually. Our numerical simulations illustrate that the multispikes can arise from system (4) and they are always unstable or meta-stable. The rigorous stability analysis is beyond the scope of this paper. Figure 5 also suggests that the formation of stable patterns, i.e., solutions with spikes or stripes, etc., depends on the initial data and in particular the initial distribution of chemical concentration. Figure 5. The formation and evolution of multi-spikes. The initial data are u 0 (x, y) = 3 + 0.1 cos(πx) cos(πy), v 0 (x, y) = 3 + cos(2πx) cos(2πy). The parameters are chosen d 1 = 1, d 2 = 0.001, χ = 3 and α = β = c = 1.