Existence and metastability of non-constant steady states in a Keller-Segel model with density-suppressed motility

We are concerned with stationary solutions of a Keller-SegelModel with density-suppressed motility and without cell proliferation. we establish the existence and the analytical approximation of non-constant stationary solutions by applying the phase plane analysis and bifurcation analysis. We show that the one-step solutions is stable and two or more-step solutions are always unstable. Then we further show that two or more-step solutions possess metastability. Our analytical results are corroborated by direct simulations of the underlying system.


Introduction
Stripe pattern formation was observed in the experiment of the engineered E. coli strains with the behavior of density suppressing motility in an isolated apparatus (see [9]), which showed that spatio-temporal patterns could be driven by a "self-trapping" mechanism besides diffusion-driven and chemotaxis-driven instabilities [10]. In order to describe the essential features of the stripe pattern formation driven by the density-suppressed motility, in [1] authors proposed the following model x ∈ Ω, t > 0, v t = D∆v + ηu − βv, x ∈ Ω, t > 0, ∂u ∂ν = ∂v ∂ν = 0, x ∈ ∂Ω, t > 0, u(x, 0) = u 0 (x), v(x, 0) = v 0 (x), x ∈ Ω, (1.1) where the domain Ω ⊂ R N , N ≥ 1 is bounded and has a smooth boundary ∂Ω, ν is the outward unit normal vector on ∂Ω; u(x, t) and v(x, t) denote the densities of E. coli cells and chemical substance acyl-homoserine lactone (AHL), respectively; The chemical substance AHL is produced by E. coli cells with a rate η > 0, degraded with a rate β > 0, and diffused with a rate D. E. coli cells have a logistic growth with an intrinsic rate σ ≥ 0 saturated at the normalized density 1 and a non-random diffusion with the diffusion rate r(v). The motility function r(v) decreases as the density of AHL increases, i.e., dr dv < 0, which implies the suppressing effect of AHL concentration on cell's motility.
When σ > 0, as far as we know, there have been the following study to (1.1). The existence of global classical solutions and the stability of constant steady state for Ω ⊂ R 2 were investigated in [2].
The global existence of classical solutions was recently extended to higher dimensions (N ≥ 3) under appropriate conditions in [14]. In [12] the authors studied the dynamics of interface of discontinuity of solutions when r(v) is a piecewise constant function.
For the case where σ = 0, global classical solutions in two dimensions and global weak solutions in three dimensions were established in [13] by supposing that r(v) has a positive lower and upper bounds. In [15], the authors studied the global existence of classical solutions, the stability of constant steady states and the existence of non-constant solutions in any dimensions for the motility function r(v) given by r(v) = c 0 /v p , p > 0, c 0 > 0 and small enough. (1.2) Specifically, the reference [15] dealt with the model as follows: x ∈ Ω, t > 0 (1. 3) with the following boundary condition and initial value ∂u where the initial functions u 0 (x) ≥ 0 and v 0 (x) > 0 are smooth. The first equation of (1.3) is also in the form of Obviously, the system (1.3) is a special form of original Keller-Segel model [3] It is seen that the case of L = 0 corresponds to the model (1.3) when both η(v) and β(v) are positive constants. For the biological interpretation of L = 0 and L ∈ (0, 1) as well as other more details about (1.3), we refer interested readers to [15]. The Keller-Segel models [3] have been extensively investigated for a cell aggregation phenomenon and the global existence and boundedness of solutions in various forms, for example, see [4,5,6,16,17]. The aim of this paper is to establish the existence of non-constant steady states of (1.3)-(1.4) by using a different method from [15] and to derive conditions for their stability and metastability in one-dimensional space. For some different models the metastability was discussed, for example, see [7,8,11]. Throughout the paper, we assume (1.2) is true.
Our presentation is structured as follows. In Section 2 we give some properties of the negative Laplace operator used later, the interval and the number of unstable modes and the sufficient conditions for the instability of constant steady state, and the expression of the most unstable mode. In Section 3 we establish the existence of non-constant steady states by applying the phase plane analysis and the third-order approximate expression of the local bifurcations. In Section 4, we analyze the stability and metastability of the stationary solutions with small amplitudes. Full numerical solutions of the original system are also carried out to corroborate the results of our analytical analysis. Finally, in Section 5 we conclude our work and bring some forward problems for further study.

Preliminaries
In this section we give some results which will play an important role in our later discussion. We first present one known property of the negative Laplace operator −∆ in the interval [0, l], where l is a positive real number. The eigenvalue problem has a sequence of simple eigenvalues λ j = (jπ/l) 2 , j = 0, 1, 2, · · ·, (2.8) whose corresponding eigenfunctions are The steady states of (1.3)-(1.4) with N = 1 satisfy the elliptical boundary-value problem (2.10) Obviously, system(1.3) possesses one conserved quantity where M is an implicit positive parameter. Then (2.10) has the constant solution (M, η β M ). Naturally, we first study its stability to establish the unstable mode band and its relationship with the system parameters. (ii) the number of the unstable Fourier modes j max is equal to the greatest j satisfying (2.12), i.e., Here [·] denotes the integer part.
Then substitute this into (1.3) to obtain the linearized system of U and It is well know that the jth mode cos( πjx l ) grows at the exponential function exp(µt), where µ is the larger eigenvalue of the matrix This gives the characteristic equation of µ Thus, the discriminant of (2.14) is (2.15) By this, we have .
(2. 16) and µ > 0 if and only if detA < 0, i.e., 0 < λ < λ * , which implies that on a bounded domain [0, l] there is a finite set of unstable modes j satisfies Then (i) and (ii) are proved. Next, we prove (iii). We regard µ as a continuous function of λ and assume that the maximum of µ(λ) attains at λ u . To find the value of λ u we have to solve the equation dµ dλ = 0. By (2.16), we introduce a new function h(µ) = 2µ − trA, hence dµ dλ = 0 is equivalent to Then taking the square of both sides and substituting the expression of h 2 and dh dλ leads to To simplify the expression we let .
Then (2.18) can be rewritten as which has a unique positive root .
It is easy to check that Then (iii) follows.

Existence and analytical approximation
3.1. Existence of non-constant steady states. This subsection is devoted to the discussion of the existence of nonconstant solutions to the stationary system (2.10) by using the method of phase plane analysis. To this end, we apply the first equation of (2.10) with the given boundary conditions to get Applying the conserved quantity M , we know that the integration constant v 0 is determined by Integrating the second equation of (2.10) from 0 to l and using the boundary conditions yield an additional information for v(x) Substituting (3.19) into the the second equation of (2.10) leads to the single second-order equation In [18] the authors established the conditions for the existence and the attractivity of positive nonconstant solutions to (3.23).
Let v x = ω and regard v 0 (v) in (3.19) as a given constant. Then (3.23) is equivalent to the first-order system    v x = w, which is a Hamiltonian system if letting the variable x act as the time. Its Hamiltonian function is Thus system (3.25) has only two types of fixed points, i.e., saddles and centers. It is easy to observe that (3.25) has two or three fixed points denoted by (v k , 0) satisfying and the eigenvalues ρ of the fixed point (v k , 0) solve the equation The type and the number of the fixed points explicitly depend on the values of the parameter p. By a straightforward computation, we have the proposition below.  To obtain spatially positive and inhomogeneous solutions (i.e., non-constant positive solutions) of (3.23)-(3.24) , it is well known to require that the trajectory of (3.25) on the phase plane meets two conditions: (i) it begins and ends at the line ω = 0 to satisfy the boundary value conditions; (ii) the transition time x between this two points is equal to l. Then, based on the results of Bendixson and Poincaré, only some of the periodic trajectories circling the center can be candidates. Therefore, by Proposition 3.1, (3.23)-(3.24) (equivalent to (2.10)) has no non-constant positive solutions when 0 < p ≤ 1.
Next we discuss the case where p > 1. Notice of (B) in Proposition 3.1, there are two types of phase portraits for four different values of p, which are shown in Fig.1. Let v * ∈ [v 0 , v 1 ] and (v * , 0) is the point where the trajectory touches the V −axis. We use l(v * ) to represent the length of a half circle which ends at (v * , 0). If v * approaches v 0 , then the corresponding orbit approaches a homoclinic, then If v * approaches v 1 , then by the linearized system of (3.25) around the center v 1 , the length of the half circle is lim Thus, when the interval length l > l * , there is at least one non-constant steady state in (3.23)-(3.24). Through the above discuss, we conclude the following result on solutions to (2.10).
Theorem 3.2. If 0 < p ≤ 1, then (2.10) has only constant positive solutions. Let the positive parameters D, β be fixed, and assume that l > l * . Then (2.10) has at least one non-constant solution provided that p > 1. . This result implies that the occurrence of non-constant solutions of (2.10) requires β D to be large for some fixed p. This can be done by giving a sufficiently small diffusion rate D and sufficiently large degradation rate β of the chemical substance AHL.

3.2.
Analytical expression of local bifurcation. Next under the condition that p > 1, we shall establish the third-order approximations of non-constant solutions with small amplitudes of (2.10) (equivalent to (3.23)-(3.24)). The above phase plane analysis shows that non-constant solutions bifurcate from the spatially homogenous solution (M, η β M ) having at least one unstable mode. If we treat β as a bifurcation parameter, then, by Theorem 3.2, the bifurcation points are where, for two given positive integers j 1 and j 2 , Notice of (2.12), from (3.28) it follows that β 1 0 is the smallest bifurcation point β min . Let (u * , v * ) be a non-constant solutions of (2.10). Then v * solves (3.23)-(3.24)). we now make an asymptotic analysis for v * with a small amplitude by assuming with β 0 = β j 0 , β 1 = β j 1 · · · , j = 1, 2, · · · and u n = n j=1 a nj cos πjx l , n = 1, 2, · · ·.
Here the positive parameter 0 < ε 1 keeps the parameter β is in the small neighborhood of the bifurcation location β 0 so that the corresponding small amplitude solution (u * , v * ) can be bifurcated at this location from the constant steady state (M, η β M ). Substituting (3.30) into (3.23)-(3.24), under the condition (3.29) we have β 1 = 0 and as well as and where d 1 (j) and d 2 (j) are defined in (3.35). Fig.2 compares the long-time numerical steady states (see left panels) with the prediction (3.36) from our asymptotic analysis (see right panels) for the bifurcations with the principal wave modes j = 3 and j = 4, respectively. For the sake of brevity, only the numerical results of the solution component v are presented here. As it can be noticed from the figure, there is a qualitative agreement between the full numerics and the analytical prediction. The variation in amplitude originates from omitting higher order terms in the analysis. Therefore, we shall use the expression in (3.36) to analyze the stability of (u * , v * ) by estimating the sign of the principal eigenvalue.

Stability and metastability
This section is devoted to the analysis of stability and metastability for non-constant steady states with small amplitudes.

Stability analysis.
Following the method used in [19], we first indicate the relationship between the solution (u * , v * ) and its bifurcation location β j 0 by relabelling (u * , v * ) as (u * j , v * j ). For system (4.37) Substitution of (4.37) into (1.3)-(1.4) yields a linear system Then substituting the asymptotic expansions of u * j , v * j and β in (3.30) and into (4.38) and equating the O(1) terms lead to The sign of γ 0 determines the stability of the stationary solution (u * j , v * j ). To solve the eigenvalue problem (4.40) for γ 0 , in view of (2.7)-(2.9) we replace (ϕ 0 , ψ 0 ) with −λ m (ϕ 0 , ψ 0 ) for some integer m ≥ 0. Thus, the characteristic equation is It is easy to see that β j 0 = β min when j = 1, and then there exists a integer m = 1 such that τ < 0. Hence (4.41) has a positive root γ 0 . We now have the following result.
In other words, the stable non-constant steady state with small-amplitude is always located on the first bifurcation.

4.2.
Metastability of multi-step solutions. Solutions located on the jth bifurcation possess the principal wave mode j. Hence we call them j-step solutions. As obtained in the previous sections, under the condition (4.50) the one-step solution is stable and multi-step ones are unstable. By the equation (4.41), we know that the principal eigenvalues of multi-step solutions correspond to m = 1, that is, where σ(j) = λ 1 r( ηM Through a simple analysis of (4.51), we know that γ 0 (j, l) is an increasing function of j for some fixed l, and if the wave mode j is fixed, then γ 0 (j, l) is decreasing with the increase of length l; When the length l is sufficiently large, it is obvious that, for all j the principal eigenvalue γ 0 is close to zero. These dependency are also explicitly shown in Fig.4. Thus, two or more-step solutions have metastability; Moreover, the less steps the stationary solution has, the more stable it is.  We now give a numerical example to demonstrate our theoretical results. System parameters are taken as c 0 = 0.1, p = 2, D = 0.5, β = 0.4 and η = 0.4. The interval length is chosen as l = 20. The mean density is set as M = 0.5. By a computation, we have that (1) the upper bound of unstable modes is λ * = 0.8; (2) the number of unstable Fourier modes is j max = 6; (3) the most unstable mode is j = 3 corresponding to λ u = 0.2492. Fig.5 demonstrates all of above analytical results obtained by applying Lemma2.1, Theorem 3.2, Theorem 4.2 and our metastability analysis. The initial data is a random perturbation of the spatially homogeneous background (M, ηM β ) = (0.5, 0.5). As observed in Fig.5, the most unstable three-step solution appears first at about t = 60. Its left step disappears at about t = 1.1 × 10 3 , and then a two-step solution develops. The two-step solution finally becomes a stable one-step solution at about t = 1.1 × 10 5 . Hence, the three-step solution persists for about 10 3 time units and the two-step solution exits for about 10 6 time units. When we vary system parameters, especially increase the length of interval, we do observe that every transient period of duration will become longer, which is not shown here. Since such non-constant steady states that stays almost unchanged for a rather long period may not be distinguishable from true stable steady sates, we call them metastable steady states.

Conclusion
In this work, in one-dimensional space we obtain the conditions for the existence of non-constant steady states of (1.3) by using the phase plane analysis. Then relying on the bifurcation analysis by treating the decay rate of chemical substance as a bifurcation parameter, we derive the thirdorder approximate expression of non-constant steady states with small amplitudes. Furthermore, the stability of one-step solutions and the metastability of two or more-step solutions are established. The analytical results are corroborated by the numerical computation as well as by numerical simulations of the underlying Keller-Segel system. The metastable states are the phenomenon that can be seen in experiments. Thus the establishment of Metastability of stationary solutions will provide theoretical . The initial data is (u 0 (x), v 0 (x)) = (M + rand(1), ηM β + rand (1)). The blue and solid curve is u(x). The red and dashed curve is v(x).
basis for experimental study related to the model (1.3). The analytical results show that E. coli cells and chemical substance AHL, under the condition that other environmental factors remain unchanged, will be inhomogenously distributed when the diffusion rate of E. coli cells is sufficiently small and the degradation rate of AHL is sufficiently large (see Theorem 3.2 ), and their state always appears to be stable when the spatial domain is large enough (see (4.51)). This makes good biological sense.
However, how to explain the formation of metastability and how to understand mergings and dissolvings of steps of non-constant steady states are not explored here. This is an interesting problem. Another challenging possibility is to consider the metastability of (1.3) in two or higher dimensional spaces.