Mathematics in Applied Sciences and Engineering https://doi.org/10.5206/mase/8120 Volume 1, Number 1, March 2020, pp.1-15 https://ojs.lib.uwo.ca/mase EXISTENCE AND METASTABILITY OF NON-CONSTANT STEADY STATES IN A KELLER-SEGEL MODEL WITH DENSITY-SUPPRESSED MOTILITY PENG XIA, YAZHOU HAN, JICHENG TAO, AND MANJUN MA Abstract. We are concerned with stationary solutions of a Keller-Segel Model 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 numerical simulations of the underlying system. 1. 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  ut = ∆(r(v)u) + σu(1 −u), x ∈ Ω, t > 0, vt = D∆v + ηu−βv, x ∈ Ω, t > 0, ∂u ∂ν = ∂v ∂ν = 0, x ∈ ∂Ω, t > 0, u(x, 0) = u0(x),v(x, 0) = v0(x), x ∈ Ω, (1.1) where the domain Ω ⊂ RN,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 Ω ⊂ R2 were investigated in [2]. Received by the editors 1 July 2019; revised 27 September 2019; accepted 27 September 2019; published online 30 September 2019. 2000 Mathematics Subject Classification. 35K55, 35K45, 35K57, 35K50, 92C15, 92C17. Key words and phrases. Keller-Segel Model, density-suppressed motility, metastability, nonconstant steady states. Manjun Ma and Peng Xia were supported by the National Natural Science Foundation of China (No. 11671359). Yazhou Han was supported by the provincial Natural Science Foundation of Zhejiang (No. LY18A010013)). Jicheng Tao was supported by the provincial Natural Science Foundation of Zhejiang (No. LY16A010009)). 1 2 PENG XIA, YAZHOU HAN, JICHENG TAO, AND MANJUN MA 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) = c0/v p, p > 0, c0 > 0 and small enough. (1.2) Specifically, the reference [15] dealt with the model as follows:{ ut = ∇· ( r(v) ( ∇u− p v u∇v )) , x ∈ Ω, t > 0, vt = D∆v + ηu−βv, x ∈ Ω, t > 0 (1.3) with the following boundary condition and initial value{ ∂u ∂ν = ∂v ∂ν = 0, x ∈ ∂Ω, t > 0, u(x, 0) = u0(x),v(x, 0) = v0(x), x ∈ Ω, (1.4) where the initial functions u0(x) ≥ 0 and v0(x) > 0 are smooth. The first equation of (1.3) is also in the form of ut = ∇· (r(v)∇u + r′(v)u∇v). Obviously, the system (1.3) is a special form of original Keller-Segel model [3]{ ut = ∇· (γ(v)∇u−χ(v)u∇v) , x ∈ Ω, t > 0, vt = D∆v + η(v)u−β(v)v, x ∈ Ω, t > 0 (1.5) with χ(v) = (L− 1)γ′(v), γ′(v) < 0, L ∈ [0, 1). (1.6) 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. EXISTENCE AND METASTABILITY OF NON-CONSTANT STEADY STATES 3 2. 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  −φ′′ = λφ, x ∈ (0, l), φ′ = 0 at x = 0, l, (2.7) has a sequence of simple eigenvalues λj = (jπ/l) 2, j = 0, 1, 2, · · ·, (2.8) whose corresponding eigenfunctions are φj(x) = { 1, j = 0, cos(jπx/l), j > 0. (2.9) The steady states of (1.3)-(1.4) with N = 1 satisfy the elliptical boundary-value problem  ( r(v) ( u′ − p v uv′ ))′ = 0, x ∈ (0, l), Dv′′ + ηu−βv = 0, x ∈ (0, l), u′(0) = u′(l) = 0, v′(0) = v′(l) = 0. (2.10) Obviously, system(1.3) possesses one conserved quantity M = 1 l ∫ l 0 u(x,t)dx = 1 l ∫ l 0 u(x, 0)dx, (2.11) 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. Lemma 2.1. If there exists j ∈ N satisfies 0 < ( πj l )2 < λ∗, λ∗ = β D (p− 1), (2.12) then (i) the constant steady state solution of (1.3) is linearly unstable. (ii) the number of the unstable Fourier modes jmax is equal to the greatest j satisfying (2.12), i.e., jmax = l √ λ∗ π − 1 if l √ λ∗ π is a positive integer; otherwise, jmax = [ l √ λ∗ π ] . (iii) the most unstable mode ju = [ l √ λu π ] or ju = [ l √ λu π ] + 1, where λu > √ r(ηM β )D(√ r(ηM β ) + √ D )2 λ∗; Furthermore, the wave number √ λu is monotone increasing in β if 1 < p ≤ 2; When p > 2, √ λu is monotone decreasing if β ∈ [β∗, +∞) or is monotone increasing if β ∈ (0,β∗), where β∗ = ηMc −1 p 0 ( (p + 2) √ D p− 2 )2 p . Here [·] denotes the integer part. 4 PENG XIA, YAZHOU HAN, JICHENG TAO, AND MANJUN MA Proof. Let u = M + U(x,t), v = ηM β + V (x,t). Then substitute this into (1.3) to obtain the linearized system of U and V   Ut = r( ηM β )U′′ + r′( ηM β )MV ′′, x ∈ Ω, t > 0, Vt = DV ′′ −βV + ηU, x ∈ Ω, t > 0, Ux = Vx = 0,x = 0, l. (2.13) 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 A(λ) = ( −λr(ηM β ) −r′(ηM β )Mλ η −Dλ−β ) , λ = ( πj l )2. This gives the characteristic equation of µ µ2 − trAµ + detA = 0, (2.14) where trA = − ( r( ηM β ) + D ) λ−β < 0, det A = −Dr( ηM β )λ(λ∗ −λ). Thus, the discriminant of (2.14) is D =(trA)2 − 4detA = [( r( ηM β ) + D ) λ + β ]2 + 4Dr( ηM β )λ(λ∗ −λ) = ( r( ηM β ) −D )2 λ2 + 2 ( r( ηM β ) + D ) βλ + β2 + 4Dr( ηM β )λλ∗ ≥ ( r( ηM β ) −D )2 λ2 − 2 ( r( ηM β ) −D ) βλ + β2 = [( r( ηM β ) −D ) λ−β ]2 ≥ 0. (2.15) By this, we have µ = trA + √ D 2 . (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 0 < j < l √ λ∗ π . 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 dh dλ = − trA dλ = r( ηM β ) + D. (2.17) Applying (2.15) and (2.16), we have h2 = (trA)2 −4detA. Differentiation of this expression by λ yields 2hhλ = 2trA(trA)λ − 4(detA)λ. EXISTENCE AND METASTABILITY OF NON-CONSTANT STEADY STATES 5 Then taking the square of both sides and substituting the expression of h2 and dh dλ leads to 4 [( (r( ηM β ) + D ) λ + β)2 + 4Dr( ηM β )λ(λ∗ −µ) ]( r( ηM β ) + D )2 = [ 2 ( (r( ηM β ) + D ) λ + β) ( r( ηM β ) + D ) + 4Dr( ηM β )(λ∗ − 2λ) ]2 . (2.18) To simplify the expression we let τ = r(ηM β ) + D Dr(ηM β ) . Then (2.18) can be rewritten as( τ2Dr( ηM β ) − 4 ) λ2 + 2λ(τβ + 2λ∗) −λ∗(λ∗ + τβ) = 0, which has a unique positive root λu = −2(τβ + 2λ∗) + √ (4(τβ + 2λ∗)2 + 4(τ2Dr(ηM β ) − 4)λ∗(λ∗ + τβ))) 2(Dr(ηM β )τ2 − 4) . It is easy to check that d( λuλ∗ ) dλ∗ < 0 so that λu λ∗ > lim λ∗→+∞ λu λ∗ = 1 2 + τ √ Dr(ηM β ) = √ r(ηM β )D(√ r(ηM β ) + √ D )2 . Then (iii) follows. � 3. 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 u = v0 r(v) . (3.19) Applying the conserved quantity M, we know that the integration constant v0 is determined by v0 = lMc0∫ l 0 vpdx > 0. (3.20) Integrating the second equation of (2.10) from 0 to l and using the boundary conditions yield an additional information for v(x) ∫ l 0 vdx = lMη β . (3.21) Combining (3.20) with (3.21), we conclude that value range of v0 is 0 < v0 ≤ r( ηM β )M when p > 1. (3.22) The maximum v0 = r( ηM β )M corresponds to the positive constant steady state (M, η β M). In the cases where the steady states are not constant, v0 are gotten by solving (3.20). Substituting (3.19) into the the second equation of (2.10) leads to the single second-order equation Dvxx + ηv0 c0 vp −βv = 0, (3.23) 6 PENG XIA, YAZHOU HAN, JICHENG TAO, AND MANJUN MA with the boundary conditions vx(0) = vx(l) = 0. (3.24) In [18] the authors established the conditions for the existence and the attractivity of positive non- constant solutions to (3.23). Let vx = ω and regard v0(v) in (3.19) as a given constant. Then (3.23) is equivalent to the first-order system  vx = w, Dwx = β(v −g(v)), g(v) = ηv0 βr(v) = ηv0 βc0 vp, (3.25) which is a Hamiltonian system if letting the variable x act as the time. Its Hamiltonian function is H(w,v) = 1 2 w2 + ηv0 D(p + 1)c0 vp+1 − β 2D v2. 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 (vk, 0) satisfying v 1−p k = ηv0 βc0 , (3.26) and the eigenvalues ρ of the fixed point (vk, 0) solve the equation Dρ2 = β(1 −g′(v̄k) = β ( 1 − ηv0 βc0 pv p−1 k ) . (3.27) 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. Proposition 3.1. The following statements are true: (A) For all p > 0, system (3.25) has always a fixed point v0 = (0, 0). it is a saddle if p > 1; while it is a non-differentiable point if 0 < p < 1. (B) Suppose that p > 1 is a integer or p−1 = s q , where s and q are coprime positive integers. Then we have (b1) if p is odd or s is even, then system (3.25) has two nontrivial fixed points: v1 = (( βc0 ηv0 ) 1 p−1 , 0 ) and v2 = ( − ( βc0 ηv0 ) 1 p−1 , 0 ) , and both of them are centers. (b2) if p is even or s is odd, then system (3.25) has one nontrivial fixed point v1 = (( βc0 ηv0 ) 1 p−1 , 0 ) , and it is a center. (C) Assume that 0 < p < 1 and 1 −p = s q , where s,q are coprime integers. Then we have (c1) if s is even, then system (3.25) has two nontrivial fixed points: v3 = (( ηv0 βc0 ) 1 1−p , 0 ) and v4 = ( − ( ηv0 βc0 ) 1 1−p , 0 ) ; moreover, they are saddles. (c2) if s is odd, then system (3.25) has one nontrivial fixed point v3 = (( ηv0 βc0 ) 1 1−p , 0 ) , and it is a saddle. (D) if p = 1, then there is one nontrivial fixed point (M, η β M) which is a node. For the boundary value problem (3.23)-(3.24), fixed points are its spatially homogeneous solutions (i.e., constant solutions). Even if there are more than one fixed points, only one of them satisfies (3.20) and thus it is a solution of (3.23)-(3.24). Obviously, corresponding to this constant solution, we have vk = η β M for all p > 0 and v0 = Mr( ηM β ). EXISTENCE AND METASTABILITY OF NON-CONSTANT STEADY STATES 7 (a) (b) (c) (d) Figure 1. Typical phase portraits for (3.25) with D = 1,β = 1,η = 1 and M = 2. (a) p = 2; (b) p = 3; (c) p = 7/3; (d) p = 8/3. 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∗ ∈ [v0,v1] 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 v0, then the corresponding orbit approaches a homoclinic, then lim v∗→v0 l(v∗) = +∞. If v∗ approaches v1, then by the linearized system of (3.25) around the center v1, the length of the half circle is lim v∗→v1 l(v∗) = π√ β D (g′(v1) − 1) = π√ β D (p− 1) def = l∗. 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 parame- ters D,β be fixed, and assume that l > l∗. Then (2.10) has at least one non-constant solution provided that p > 1. 8 PENG XIA, YAZHOU HAN, JICHENG TAO, AND MANJUN MA Remark 3.1. Theorem 3.2 is just the result of Proposition 3.1 and Theorem 3.3 with N = 1 in [15]. 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 bifur- cate 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 β0 = Dλj p− 1 = D(πj l )2 p− 1 def = β j 0,j = 1, 2, 3, · · ·, (3.28) where, for two given positive integers j1 and j2, D(πj1 l )2 p− 1 6= D(πj2 l )2 p− 1 if j1 6= j2. (3.29) Notice of (2.12), from (3.28) it follows that β10 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  v∗ = ηM β + εv1 + ε 2v2 + ε 3v3 + · · ·, v0 r(v) = u∗ = M + εu1 + ε 2u2 + ε 3u3 + · · ·, β = β0 + εβ1 + ε 2β2 + ε 3β3 + · · · (3.30) with β0 = β j 0,β1 = β j 1 · · · ,j = 1, 2, · · · and un = n∑ j=1 anj 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{ u1 = c(j) cos( πjx l ), v1 = cos( πjx l ), (3.31) where c(j) = λjD + β0 η = pDλj (p− 1)η > 0 (3.32) as well as β2 = β30 (p 2 + 3p) 12(ηM)2 > 0, (3.33) and { u2 = d1(j) cos( 2πjx l ), v2 = d2(j) cos( 2πjx l ), (3.34) where   d1(j) = 4Dλj+β0 η d2(j) = (4p−3)Dλj (p−1)η d2(j), d2(j) = − 2r′( ηM β0 )c(j)+r′′( ηM β0 )M 4r( ηM β0 ) ( (4p−3)Dλj (p−1)η ) +4r′( ηM β0 )M . (3.35) EXISTENCE AND METASTABILITY OF NON-CONSTANT STEADY STATES 9 0 5 10 15 20 Space x 13.7 13.75 13.8 13.85 v (x , t) t = 100000 0 5 10 15 20 Space x 12.49 12.495 12.5 12.505 12.51 12.515 A p p ro x im a ti o n o f v 0 5 10 15 20 Space x 13.65 13.7 13.75 13.8 v (x , t) t = 20000 0 5 10 15 20 Space x 12.49 12.495 12.5 12.505 12.51 12.515 A p p ro x im a ti o n o f v Figure 2. Non-constant steady states to system (1.3) in the interval (0, 20) for M = 5 and different wave modes. System parameters are taken as D = 1,c0 = 0.01,p = 2,η = 1 and β = 0.4. Left: the numerical solutions obtained by integrating the full system (1.3) with the initial condition (u0(x),v0(x)) = ( M + rand(1), ηM β + rand(1) ) . Right: the analytical solutions (3.36) obtained by our asymptotic analysis with ε = 0.01 . The top line is the three-step pattern at t = 105, i.e., j = 3. The bottom line is the four-step pattern at t = 2 × 104, i.e., j = 4. By (3.30), (3.31) and (3.34), the third-order approximation of non-constant solutions with small amplitudes to system (2.10) reads{ u∗ = M + ε pDλj (p−1)η cos( πjx l ) + ε2d1(j) cos( 2πjx l ) + O(ε3), v∗ = ηM β + ε cos(πjx l ) + ε2d2(j) cos( 2πjx l ) + O(ε3), (3.36) where d1(j) and d2(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. 4. Stability and metastability This section is devoted to the analysis of stability and metastability for non-constant steady states with small amplitudes. 4.1. 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 10 PENG XIA, YAZHOU HAN, JICHENG TAO, AND MANJUN MA (1.3)-(1.4) with N = 1, we set { u = u∗j + ϕ(x)e γt, v = v∗j + ψ(x)e γt. (4.37) Substitution of (4.37) into (1.3)-(1.4) yields a linear system  r(v∗j )ϕ ′′ + r′(v∗j )u ∗ jψ ′′ + 2r′(v∗j )v ∗ j ′ϕ′ + Q1ψ ′ + Q2ψ + Q3ϕ = γϕ, x ∈ (0, l), Dψ′′ + ηϕ−βψ = γψ, x ∈ (0, l), ϕ′(0) = ϕ′(l) = 0, ψ′(0) = ψ′(l) = 0, (4.38) where Q1 = 2r ′′(v∗j )v ∗ j ′ u∗j + 2r ′(v∗j )v ∗ j ′ , (4.39) Q2 = r ′′′(v∗j )(v ∗ j ′ )2u∗j + r ′′(v∗j )(v ∗ j ) ′′u∗j + 2r ′′(v∗j )v ∗ j ′ u∗j ′ + r′(v∗j )u ∗ j ′′ , Q3 = r ′′(v∗j )(v ∗ j ′ )2 + r′(v∗j )v ∗ j ′′ . Then substituting the asymptotic expansions of u∗j , v ∗ j and β in (3.30) and  γ = γ j 0 + �γ j 1 + � 2γ j 2 + · · · , ϕ = ϕ0 + �ϕ1 + � 2ϕ2 + · · · , ψ = ψ0 + �ψ1 + � 2ψ2 + · · · , into (4.38) and equating the O(1) terms lead to  r(ηM β0 )ϕ′′0 + r ′(ηM β0 )Mψ′′0 = γ0ϕ0, x ∈ (0, l), Dψ′′0 + ηϕ0 −β0ψ0 = γ0ψ0, x ∈ (0, l), ϕ′0(0) = ϕ ′ 0(l) = 0, ψ′0(0) = ψ ′ 0(l) = 0. (4.40) 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 γ20 + σγ0 + τ = 0 (4.41) with σ = ( λmr( ηM β j 0 ) + Dλm + β j 0 ) , τ = λm(p− 1) c0(β j 0) p (ηM)p ( βm0 −β j 0 ) . It is easy to see that β j 0 6= βmin when j 6= 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. Proposition 4.1. For system (1.3)-(1.4) with N = 1, non-constant steady state (u∗j,v ∗ j ) is unstable if j ≥ 2. In other words, the stable non-constant steady state with small-amplitude is always located on the first bifurcation. Next, we look for the sufficient condition of stability of the first bifurcation. Obviously, when j = 1, the principal eigenvalue of (4.41) is γ0 = 0 corresponding to m = 1 with the eigenvector (ϕ0,ψ0) = ( pDπ2 (p− 1)l2 cos( πx l ), cos( πx l ) ) . (4.42) EXISTENCE AND METASTABILITY OF NON-CONSTANT STEADY STATES 11 Thus, we have to find the value of γ1. Now equating the O(ε) terms in (4.38) gives  r(ηM β0 )ϕ′′1 − r′( ηM β0 )Mη D ϕ1 + r′( ηM β0 )Mβ0 D ψ1 = γ1 ( ϕ0 − r′( ηM β0 )M D ψ0 ) −G1, x ∈ (0, l), Dψ′′1 + ηϕ1 −β0ψ1 = γ1ψ0, x ∈ (0, l), ϕ′1(0) = ϕ ′ 1(l) = 0, ψ′1(0) = ψ ′ 1(l) = 0, (4.43) where G1 =r ′( ηM β0 ) (v1ϕ ′′ 0 + u1ψ ′′ 0 + +2v ′ 1ϕ ′ 0 + 2u ′ 1ψ ′ 0 + u ′′ 1ψ0 + v ′′ 1ϕ0) + r′′( ηM β0 )M (v1ψ ′′ 0 + 2v ′ 1ψ ′ 0 + v ′′ 1ψ0) . Simplifying G1 by applying (3.20) and (4.42), we have G1 = −2 ( 2r′( ηM β0 )c(1)λ1 + r ′′( ηM β0 )Mλ1 ) cos( 2πx l ). (4.44) Since the solution to the adjoint system of the homogeneous system corresponding to (4.43) is{ ϕ = c(1) cos(πx l ), ψ = cos(πx l ), where c(1) = pD r′(ηM β0 )M < 0, (4.45) by solvability condition of (4.43), we have γ1 = ∫ l 0 G1ϕdx∫ l 0 [ ϕ0ϕ− r′( ηM β0 )M D ψ0ϕ + ψ0ψ ] dx = 0. Thus, we have to compute the value of γ2. In view of (4.44), we set the particular solution of (4.43) as{ ϕ1 = d3(1) cos( 2πx l ), ψ1 = d4(1) cos( 2πx l ). (4.46) Substitution of this into (4.43) yields  d3(1) = (4p−3)Dλ1 (p−1)η d4(1), d4(1) = − 2r′( ηM β0 )c(1)+r′′( ηM β0 )M 2r( ηM β0 ) ( (4p−3)Dλ1 (p−1)η ) +2r′( ηM β0 )M . (4.47) By equating the O(ε2) terms in (4.38), we have  r(ηM β0 )ϕ′′2 − r′( ηM β0 )Mη D ϕ2 + r′( ηM β0 )Mβ0 D ψ2 = γ2(ϕ0 − r′( ηM β0 )M D ψ0) − r′( ηM β0 )M D ψ0β2 −G2, x ∈ (0, l), Dψ′′2 + ηϕ2 −β0ψ2 = γ2ψ0 + ψ0β2, x ∈ (0, l), ϕ′2(0) = ϕ ′ 2(l) = 0, ψ′2(0) = ψ ′ 2(l) = 0, (4.48) where the explicit expression of G2 is omitted here since it is too cumbersome. Using the the solvability condition once more, we obtain γ2 = (p− 1)2β2 1 2 (c(1)c(1) − (p− 1)) , (4.49) 12 PENG XIA, YAZHOU HAN, JICHENG TAO, AND MANJUN MA where c(1) and c(1) are defined in (3.32) and (4.45), respectively. Thus, by this, (3.33) and Proposition 4.1, we have the following theorem. Theorem 4.2. Let the positive parameters c0,D,η and l be fixed. Assume that p > 1, β > D p− 1 (π l )2 . (4.50) Then, for any positive constant M the first bifurcation (i.e., (u∗1,v ∗ 1 )) of (2.10) is supercritical and linearly stable. All other bifurcations (i.e., (u∗j,v ∗ j ),j ≥ 2) are also supercritical but unstable. Taking into account Proposition 4.1 and Theorem 4.2, we give the bifurcation diagram in Fig.3. (a) (b) 0 0.05 0.1 0.15 0.2 0.25 0.3 β 0 50 100 150 200 250 300 v m a x ,m in 0.02 0.022 0.024 0.026 0.028 0.03 β 100 150 200 250 300 v m a x ,m in Figure 3. Bifurcation diagram of vmax,min as a function of β. The mean density is set as M = 5. Other parameters are taken as D = 1,η = 1,p = 2 and l = 20. The blue and thick curve shows the constant solution v = ηM β . The black and thin curves correspond to bifurcations. Solid curves correspond to stable solutions, dashed to unstable solutions. The constant solution loses its stability at the first bifurcation point where the non-constant solution with the principal wave mode 1 appears. (b) is an enlargement of the first bifurcation point from (a). 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, γ0(j, l) = −σ(j) + √ σ(j)2 − 4τ(j) 2 > 0, j = 2, 3, · · ·, (4.51) where σ(j) = ( λ1r( ηM β j 0 ) + Dλ1 + β j 0 ) = ( λ1r( ηM β j 0 ) + Dλ1 + Dλj p− 1 ) , τ(j) = λ1(p− 1)r( ηM β j 0 ) ( β10 −β j 0 ) = λ1Dr( ηM β j 0 ) (λ1 −λj) , λj = ( πj l )2 . 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 EXISTENCE AND METASTABILITY OF NON-CONSTANT STEADY STATES 13 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. 10 20 30 40 50 l 0 0.2 0.4 0.6 0.8 1 γ 0 ×10 -8 j=2 j=3 j=4 j=5 Figure 4. Examples of dependence of principal eigenvalues γ0 on l and j for 2-5-step solutions. The mean density is M = 5 and the interval length is l = 50. System parameters are set as D = 0.1,c0 = 0.01,η = 1 and p = 2. We now give a numerical example to demonstrate our theoretical results. System parameters are taken as c0 = 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 jmax = 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, The- orem 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 × 103, and then a two-step solution develops. The two-step solution finally becomes a stable one-step solution at about t = 1.1×105. Hence, the three-step solution persists for about 103 time units and the two-step solution exits for about 106 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. 5. 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 third- order 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 14 PENG XIA, YAZHOU HAN, JICHENG TAO, AND MANJUN MA 0 5 10 15 20 Space x 0.54 0.545 0.55 0.555 u , v t = 1 0 5 10 15 20 Space x 0.54 0.545 0.55 0.555 0.56 u , v t = 50 0 5 10 15 20 Space x 0.52 0.54 0.56 0.58 u , v t = 60 0 5 10 15 20 Space x 0 0.5 1 1.5 2 2.5 u , v t = 300 0 5 10 15 20 Space x 0 1 2 3 u , v t = 1000 0 5 10 15 20 Space x 0 1 2 3 4 u , v t = 1100 0 5 10 15 20 Space x 0 1 2 3 4 u , v t = 30000 0 5 10 15 20 Space x 0 1 2 3 4 u , v t = 110000 0 5 10 15 20 Space x 0 2 4 6 8 u , v t = 115000 Figure 5. Example of metastable steady states in (1.3). The initial data is (u0(x),v0(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 dis- solvings 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. Acknowledgement. We would like to thank the anonymous reviewers for their valuable comments and suggestions which greatly improved the exposure of this manuscript. EXISTENCE AND METASTABILITY OF NON-CONSTANT STEADY STATES 15 References [1] X. Fu, L.-H. Tang, C. Liu, J. -D. Huang, T. Hwa, P. Lenz, Stripe formation in bacterial systems with density -suppressed motility, Physical Review Letters 108 (2012): 198102. [2] H. Y. Jin, Y. J. Kim, and Z. A. Wang, Boundedness, stabilization, and pattern formation driven by density- sup- pressed motility , SIAM J. Appl. Math. 78 (2018), 1632-1657. [3] E. F. Keller and L. A. Segel, Initiation of slime mold aggregation viewed as an instability, J. Theor. Biol. 26 (1970), 399-415. [4] T. Hillen and K.J. Painter, A users guide to PDE models for chemotaxis, J. Math. Biol. 58 (2009), 183-217. [5] M. A. Herrero and J. J. Velázquez, A blow-up mechanism for a chemotaxis model, Ann. Sc. Norm. Super. Pisa, Cl. Sci. Ser. IV 24 (1997), 633-683. [6] D. Horstmann, From 1970 until present: the Keller-Segel model in chemotaxis and its consequences , I Jahresberichte DMV 105 (2003), 103-165. [7] T. Kolokolnikov, T. Erneu and J. Wei, Mesa-type patterns in the one-dimensional Brusselator and their stability, Physica D: Nonlinear Phenomena 214 (2006), 63-77. [8] K. Kang, T. Kolokolnikov and M. J. Ward, The stability and dynamics of a spike in the one-dimensional Keller-Segel model, IMA Journal of Applied Mathematics 72 (2007), 140-162. [9] C. Liu, X. Fu, L. Liu, X. Ren, C. K. Chau, S. Li, L. Xiang, H. Zeng, G.Chen, L.-H. Tang, et al., Sequential establishment of stripe patterns in an expanding cell population, Science 334 (2011), 238-241. [10] J. Murray, Mathematical Biology II: Spatial Models and Biomedical Applications. Springer, 2003. [11] A. Potapov and T. Hillen, Metastability in chemotaxis models, J. Dyna. Diff. Eqns. 17 (2005), 293-330. [12] J. Smith-Roberge, D. Iron and T. Kolokolnikov, Pattern formation in bacterial colonies with density-dependent diffusion, Europ. J. Appl. Math. 30 (2019), 196-218. [13] Y. Tao and M. Winkler, Effects of signal-dependent motilities in a Keller-Segel-type reaction-diffusion system, Math. Models Meth. Appl. Sci. 27 (2017), 1645-1683. [14] J. Wang and M. Wang, Boundedness in the higher-dimensional Keller-Segel model with signal-dependent motility and logistic growth, J. Math. Phys. 60 (2019): 011507. [15] C. Yoon and Y.-J. Kim, Global existence and aggregation in a Keller-Segel model with Fokker-Planck diffusion, Acta Application Mathematics 149 (2017), 101-123. [16] C. Yoon and Y.-J. Kim, Bacterial chemotaxis without gradient sensing, J. Math. Biol. 70 (2015), 1359-1380. [17] T. Senba and T. Suzuki, Parabolic system of chemotaxis: blowup in a finite and the infinite time, Meth. Appl. Anal. 8 (2001), 349-367. [18] C. S. Lin, W. M. Ni, and I. Takagi, Large amplitude stationary solutions to a chemotaxis system, J. Diff. Eqns. 72 (1988), 1-27. [19] M. J. Ma, C.H. Ou, and Z.A. Wang. Stationary solutions of a volume filling chemotaxis model with logistic growth and their stability, SIAM J. Appl. Math. 72 (2012), 740-766. Department of Mathematics, School of Sciences, Zhejiang Sci-tech University, Hangzhou, Zhejiang, 310018, China E-mail address: 624598686@qq.com Department of Mathematics, College of Science, China Jiliang University, Hangzhou, Zhejiang 310018, China E-mail address: taoer558@163.com Department of Mathematics, College of Science, China Jiliang University, Hangzhou, Zhejiang 310018, China E-mail address: 02a0802031@cjlu.edu.cn Corresponding author. Department of Mathematics, School of Sciences, Zhejiang Sci-tech University, Hangzhou, Zhejiang, 310018, China E-mail address: mjunm9@zstu.edu.cn