Global bifurcation and stability of steady states for a bacterial colony model with density-suppressed motility

# Global bifurcation and stability of steady states for a bacterial colony model with density-suppressed motility

Manjun Ma Department of Mathematics, School of Sciences, Zhejiang Sci-tech University, Hangzhou, Zhejiang, 310018, China; mjunm9@zstu.edu.cn    Peng Xia Department of Mathematics, School of Sciences, Zhejiang Sci-tech University, Hangzhou, Zhejiang, 310018, China; 624598686@qq.com    Qifeng Zhang Department of Mathematics, School of Sciences, Zhejiang Sci-tech University, Hangzhou, Zhejiang, 310018, China; zhangqifeng0504@163.com    Matti Vuorinen Department of Mathematics and Statistics, University of Turku, FIN-20014 Turku, Finland; vuorinen@utu.fi

FILE: MXversion23.tex

Abstract: We investigate the structure and stability of the steady states for a bacterial colony model with density-suppressed motility. We treat the growth rate of bacteria as a bifurcation parameter to explore the local and global structure of the steady states. Relying on asymptotic analysis and the theory of Fredholm solvability, we derive the second-order approximate expression of the steady states. We analytically establish the stability criterion of the bifurcation solutions, and show that sufficiently large growth rate of bacteria leads to a stable uniform steady state. While the growth rate of bacteria is less than some certain value, there is pattern formation with the admissible wave mode. All the analytical results are corroborated by numerical simulations from different stages.

Keywords: Density-suppressed motility, reaction-diffusion model, global bifurcation, stability analysis

MR Subject Classification: 35K55, 35K45, 35K57.

## 1 Introduction and preliminaries

We consider the following nonlinear reaction-diffusion system

 (1.1)

This model was first introduced in [1] to describe the dynamical behavior of the bacterial species Vibrio fischeri’s colonies. Here is a bounded open domain in , is a positive integer, and . The quantities and stand for the density of bacteria and acy-homoserine lactone (AHL) secreted by Vibrio fischeri, respectively. The positive constants and measure, respectively, the logistic growth rate of the bacteria and the diffusion rate of AHL. The diffusion rate of bacteria is state-dependent on modeled by the positive motility function . At the present, for the model (1.1), the rigorous mathematical results are very limited. The mechanism of stripe formation of (1.1) was analyzed in [1] when is a piecewise decreasing function. In [2], authors discussed the pattern solutions and their stability when the diffusion rate of has a drop at some critical AHL concentration, that is, is a step function. The apriori bound, the global existence of classical solutions, the non-existence of pattern solutions and the numerical results of pattern formation and wave propagation were established in [3] when the system (1.1) is located in a two-dimensional bounded domain with zero Neumann boundary conditions and some conditions are imposed on the motility function . When (namely bacteria have no growth), some results on the existence of global solutions were obtained in [5, 6].

As stated in [3], the model (1.1) can be transformed to

 (1.2)

which is a chemotaxis model of Keller-Segel type proposed in [7] if the cells/ bacteria do not sense the concentration between receptors, more details can be found in [7] and [8]. The problem of how diffusion, logistic growth and chemotaxis can interact and how such interaction influences the dynamics of some particular systems has been extensively investigated in the literature, see [9]-[11] and references therein.

In this paper, let (1.1) be subject to the Neumann boundary conditions

 ∇u⋅ν=∇v⋅ν=0, x∈∂Ω,  t>0, (1.3)

where is the outward unit normal vector on and the domain is bounded and has smooth boundary, the initial data

 (u0,v0)∈[W1,∞(Ω)]2 and u0(x), v0(x)≥0, x∈¯¯¯¯Ω, (1.4)

there exists at least one point such that and , furthermore, let the motility function satisfy

 r(v)∈C3([0,∞)), r(v)>0 and r′(v)<0 for v∈[0,∞),limv→+∞r(v)=0. (1.5)

The condition (1.3) means that there is no flux of either bacteria or AHL across the boundary of the domain, and thus (1.1) is a closed system.

Since as , by the apriori bound of established in [3], the case of degeneracy will not happen here. Due to the assumption , we call AHL concentration being of the repressive effect on bacterium motility.

For the readers’ convenience, some notations used in this paper are presented here. Let for be Sobolev space of - valued functions with norm . When , is written as . Let denote the usual Lebesgue space in a bounded domain with norm for and . When , which is the space of -valued continuous functions.

The well-posedness of the parabolic system (1.1) with (1.3)-(1.4) was already proved in [3], see the following lemma.

###### Lemma 1.1.

(Theorem 1.1 in [3]) Assume that possesses smooth boundary and (1.5) holds. Then the problem (1.1) with (1.3)-(1.4) has a unique nonnegative global classical solution satisfying

 ∥u(⋅,t)∥L∞(Ω)+∥v(⋅,t)∥W1,∞)(Ω)≤C for all t>0,

where is a constant independent of .

The following properties of the negative Laplacian operator with zero Neumann boundary condition on will be used later. There is a sequence of eigenvalues satisfying

 0=λ0<λ1<λ2<λ3<⋅⋅⋅. (1.6)

Each has multiplicity . Let , be the normalized eigenfunctions corresponding to . Let be the eigenspace associated with in . Then the set forms a complete orthogonal basis in . Let and . Then

 X=∞⨁i=1Xi, Xi=dimS(λi)⨁j=1Xij, (1.7)

where denotes the direct sum of subspaces and .

Clearly, the non-constant steady states of the system (1.1) with (1.3)-(1.4) are necessarily positive non-constant solutions to the elliptic system

 ⎧⎪⎨⎪⎩Δ(r(v)u)+σu(1−u)=0,x∈Ω,DΔv−v+u=0,x∈Ω,∇u⋅ν=∇v⋅ν=0,x∈∂Ω. (1.8)

The boundedness of positive classical solutions of (1.8) was proved in [12].

###### Lemma 1.2.

Assume that (1.5) is true. Let be a positive classical solution of the system (1.8). For any given constant , there exists a positive constant such that

 (u,v)∈B={(u,v):1B≤u,v≤B}  for  x∈¯¯¯¯Ω (1.9)

provided that , where is independent of and .

It is obvious that the system (1.8) has two constant solutions, i.e., and for all . Linearizing (1.8) at and respectively, by a simple computation, we know that is always unstable. For the linearized system around the point , the eigenvalue satisfies

 ρ2i+[(D+r(1))λi+1+σ]ρi+[σ+r(1)λi](1+Dλi)+r′(1)λi=0, (1.10)

where and is defined in (1.6). By the standard stability theory, for the stability/instability of the steady state we have a critical discriminant

 σ−{−[r′(1)1+Dλi+r(1)]λi}\lx@stackreldef=σ−σi, i=1,2,⋅⋅⋅. (1.11)

It is easy to check that if there is such that , then there must exist a positive integer such that

 σi>0  for  i∈[1,ic], σic+j≤0  for  j=1,2,3,⋅⋅⋅∞. (1.12)

Note that (1.12) implies that

 r′(1)+r(1)<0. (1.13)

Set

 σa=max1≤i≤icσi=−[r′(1)1+Dλia+r(1)]λia,  ia∈[1,ic]. (1.14)

Moreover, if we regard as a real number and use (1.13), then at

 λi=1D(√−r′(1)r(1)−1) (1.15)

the maximum of , denoted by , is attained as

 σc=1D(√−r′(1)−√r(1))2. (1.16)

It is clear that , and that if such that (1.15) is true, then . We now have the lemma below.

###### Lemma 1.3.

Suppose that (1.5) and (1.12) hold. Then we have the follows:

The steady state is linearly stable if either

 r′(1)+r(1)≥0

or

 r′(1)+r(1)<0  and  σD>−(r′(1)+r(1)).

is unstable if ; Usually, we call the admissible wave number.

is linearly stable if .

Naturally, we may expect the existence of non-constant steady states as the constant solutions are unstable and figure out their structure. The aim of this paper is to study the existence and structure of positive solutions of (1.8) in one dimensional space by applying the bifurcation theory and to establish the stability/unstability of each bifurcation branch. By Lemma 1.3, we shall fix constants and , and regard as a bifurcation parameter satisfying

 0<σ<σa. (1.17)

The local bifurcation theory will be used to precisely describe the structure of positive solutions near the bifurcation points. The global bifurcation theory is then employed to prove that these bifurcation curves can be prolonged as long as is less than a certain critical values. The asymptotic analysis and the adjoint theory are applied to derive the expression of the steady states. Furthermore, the linearized stability theory and some analytical techniques are used to give the stability criterion of the bifurcating solutions. Numerical simulations are carried out to demonstrate all the theoretical results.

## 2 Global bifurcation

With the system (1.8) reads as

 ⎧⎪⎨⎪⎩(r(v)u)′′+σu(1−u)=0,x∈(0,l),Dv′′−v+u=0,x∈(0,l),u′(0)=u′(l)=0, v′(0)=v′(l)=0. (2.1)

The eigenvalue problem

 {−φ′′(x)=λφ(x), x∈(0,l),φ′(x)=0, x=0,l, (2.2)

has a sequence of simple eigenvalues

 λj=(πj/l)2, j=0,1,2,⋅⋅⋅, (2.3)

whose corresponding eigenfunctions are

 φj(x)={1, j=0,cos(πjx/l), j>0. (2.4)

Obviously, the set of eigenfunctions constitutes an orthogonal basis in . Let

 X={(u,v):u,v∈C2([0,l]),u′=v′=0at x=0, l},

then is a Banach space with the usual norm, and is a Hilbert space with the inner product

 (ω1,ω2)Y=(u1,u2)L2(0,l)+(v1,v2)L2(0,l)

for , . By expanding the second-order derivative term in the first equation, we have the system (2.1) in the form of

 ⎧⎪⎨⎪⎩r′′(v)v′2u+r′(v)v′′u+2r′(v)v′u′+r(v)u′′+σu(1−u)=0, x∈(0,l),Dv′′+u−v=0, x∈(0,l),u′(0)=u′(l)=0, v′(0)=v′(l)=0. (2.5)

Define the map by

 P(σ,ω)=(r′′(v)v′2u+r′(v)v′′u+2r′(v)v′u′+r(v)u′′+σu(1−u)Dv′′+u−v.),

where and is a bounded set in . Hence, looking for the solutions of (2.1) is exactly equivalent to looking for the zero points of this map. Let , then we have

 P(σ,ω∗)=0 for σ>0.

We recall that, for a number , is a bifurcation point of the equation with respect to the curve if every neighborhood of contains zeros of in not lying on this curve. Then the results on local bifurcation of solutions for (2.1) are as follows.

###### Theorem 2.1.

Suppose that (1.5), (1.12) and (1.17) are true. If is a positive integer such that

 λj<−r′(1)+r(1)Dr(1), and σj≠σk for all integers k≠j, (2.6)

then is a bifurcation point of with respect to the curve , where is defined in (1.11). Furthermore, there is a one-parameter family of non-trivial solutions of the problem (2.1) for sufficiently small, where are continuous functions, and

 u(ε)=u∗+εajφj+o(ε), v(ε)=v∗+εφj+o(ε), aj=1+Dλj. (2.7)

The set of zero-points of consists of two curves and in a neighborhood of the bifurcation point .

###### Proof.

For the fixed , according to Theorem 1.7 of [13], we need to verify the following conditions:

the partial derivatives , and exist and are continuous,

and are one-dimensional,

let , then .

Because we have

 Pσ=(u(1−u)0), Pσω=(1−2u000),

and

 L=Pω(ω∗)=⎛⎜⎝r(1)∂2∂x2−σr′(1)∂2∂x21D∂2∂x2−1⎞⎟⎠,

it is clear that the linear operators , and are continuous. Condition is verified.

Let and . Then we have

 ∞∑i=0(−λir(1)−σ−λir′(1)1−λiD−1)(aibi)φi=0, (2.8)

which means that, by the definition of in (2.4), all coefficients must vanish, that is,

 (−λir(1)−σ−λir′(1)1−Dλi−1)(aibi)=0, i=0,1,2,⋅⋅⋅,∞. (2.9)

This equation has a nonzero solution provided that

 det(−λir(1)−σ−λir′(1)1−Dλi−1)=0, (2.10)

which holds if and only if

 σ=−[r′(1)1+Dλi+r(1)]λi\lx@stackreldef=σi, i=0,1,2,⋅⋅⋅,∞.

Obviously, if , then , which is excluded by the assumption of the theorem. In view of (2.6), the equation (2.10) holds only for and

 σ=σj=−[r′(1)1+Dλj+r(1)]λj.

Then, by solving the equation (2.9) with , we have

 kerL=span{Φ}, Φ=(aj1)φj, (2.11)

where . The adjoint operator of is

 L∗=⎛⎜⎝r(1)∂2∂x2−σ1r′(1)∂2∂x2D∂2∂x2−1⎞⎟⎠,

which has

 kerL∗=span{Φ∗},  Φ∗=(a∗j1)φj,

where . It is well known that . Then the codimension of is the same as . Condition is thus satisfied.

We know that

 Pσω(σj,ω∗)Φ=(−(1+Dλj)φj0)

and

 Y=(1+Dλj)2r′(1)λj<0.

Hence, , and so condition is verified. The proof is completed. ∎

###### Remark 2.1.

(a)  Theorem 2.1 and Lemma 1.3 show that if , then is a bifurcation point with respect to the trivial branch . The number of such bifurcation points is equal to the number of for which , that is for given and , see the expression (1.12). (b)  By Theorem 2.1, let be the closure of the non-trivial solution set of , then is the connected component of to which belongs, and in a neighborhood of the bifurcation point the curve is characterized by the eigenfunction . In the open interval the function has exactly zeros, thus we call the non-constant solutions in mode steady states.

We next apply the global bifurcation theory of Rabinowitz, in particular, Corollary 1.12 in [14], and the Leray-Schauder degree for compact operators to give the information on the bifurcating curve far from the trivial equilibrium. Following the idea of [15], we first rewrite the system (2.5) as

 ⎧⎪⎨⎪⎩−u′′=f(u,v),  x∈(0,l),−v′′=g(u,v),  x∈(0,l),u′(0)=u′(l)=0, v′(0)=v′(l)=0, (2.12)

where

 f(u,v)=1r(v)[r′′(v)v′2u+r′(v)v′′u+2r′(v)v′u′+σu(1−u)],  g(u,v)=1D(u−v).

Now shifting the constant state to by setting ,, then (2.12) is transformed into

 ⎧⎪ ⎪⎨⎪ ⎪⎩−~u′′=f0~u+f1~v+~f(~u,~v), x∈(0,l),−~v′′=g0~u+g1~v+~g(~u,~v), x∈(0,l),~u′(0)=~u′(l)=0, ~v(0)=~v′(l)=0, (2.13)

where and are higher-order terms of and ,

 f0=fu(1,1)=−r′(1)Dr(1)−σr(1),  f1=fv(1,1)=r′(1)Dr(1),

and

 g0=gu(1,1)=1D,  g1=gv(1,1)=−1D.

Note that and . Let and be the inverse operator of and , respectively; moreover, set

Then the boundary value problem (2.13) can be rewritten in the matrix form

 ˜ω=M(σ)˜ω+H(σ,˜ω)\lx@stackreldef=T(σ,˜ω), M(σ)=(2f0Fσf1Fσg0F0),  ˜ω∈X. (2.14)

Obviously, for any given the linear operator is compact on . On closed sub-intervals of the operator is compact on and for near zero uniformly.

The result below will play a critical role in the proof of the global bifurcation of the solutions to (2.1).

###### Lemma 2.2.

Assume that (1.5), (1.12), (1.17) and (2.6) are satisfied. Then is an eigenvalue of with algebraic multiplicity one.

###### Proof.

Let ,  . The proof of Theorem 2.1 implies that

 (M(σj)−I)Φ=⎛⎜⎝r(1)d2dx2−σjr′(1)d2dx21Dd2dx2−1⎞⎟⎠Φ=0,

has one unique solution , which shows that is an eigenvalue of with the unique eigenfunction. Thus, we have . Next we will prove that the eigenvalue is simple. It is well known that the algebraic multiplicity of is equal to the dimension of the generalized null space . So we need only to verify . Let be the adjoint operator of . If , from (2.14) it follows that

 {2fj0Fσ(φ)+gj0F(ψ)=φ,fj1Fσ(φ)=ψ, (2.15)

where

 fj0=−r′(1)Dr(1)−σjr(1),  fj1=r′(1)Dr(1),  gj0=1D.

By the definition of and , the system (2.15) can be expanded as

 {−fj1φ′′=fφφ+fψψ,−φ′′=fj1φ−fj0ψ (2.16)

with

 fφ=fj1gj1+2fj0fj1,  fψ=fj1gj0−2fj0gj1−2(fj0)2.

Again set . By (2.16), we get

 ∞∑i=0A∗i(aibi)ϕi=0, A∗i=(fϕ−fj1λifψfj1−λi−fj0).

By and (2.6), we know that if and only if and

 A∗j=(00fj1−λj−fj0).

Therefore, is generated by the unique element , and we know that has one unique solution (up to a constant multiple ) , and thus . Then , which implies that . Thus, the desired result follows. ∎

We now state the results on the global bifurcation of the boundary value problem (2.1).

###### Theorem 2.3.

Suppose that (1.5), (1.12), (1.17) and (2.6) are true. Then the projection of the bifurcation curve onto the axis is an interval . Furthermore, the system (2.1) has at least one non-constant positive solution if and for any positive integer .

###### Proof.

By the proof of Lemma 1.1, we know that the linear operator is a bijection when , and lies in a small neighborhood of as well as is an isolated solution of (2.14) for this fixed . The index of this isolated zero of the map is given by

 index(I−T(σ,.),(σ,O))=deg(I−M(σ),B,O)=(−1)β,

where is a sufficiently small ball centered at , and is the sum of the algebraic multiplicities of the eigenvalues of that are larger than . For our bifurcation analysis, we necessarily verify that this index changes as the bifurcation parameter crosses , i.e., for sufficiently small,

 index(I−T(σj−ε,⋅),(σj−ε,O))≠index(I−T(σj+ε),(σj+ε,O)). (2.17)

Indeed, if is an eigenvalue of corresponding to an eigenfunction , then we have

 {−ϱφ′′=(2−ϱ)f0φ+f1ψ,−ϱψ′′=g0φ+ϱg1ψ.

Once again let and , then the above system can be expanded as

 ∞∑i=0((2−ϱ)f0−ϱλif1g0(g1−λi)ϱ)(aibi)φi=0.

Then the set of eigenvalues of consists of all that solve the characteristic equation

 (f0+λi)(λi−g1)ϱ2−2f0(λi−g1)ϱ−g0f1=0,i=0,1,2,⋅⋅⋅. (2.18)

Taking , if is a root of (2.18), then it is concluded that

 σj=−[r′(1)(1+Dλi)+r(1)]λi=σi,

and thus by the assumption we have . Therefore, if we do not count the eigenvalues corresponding to in (2.18), has the same number of eigenvalues which are larger than for all close to , and have the same multiplicities. So we need only to consider the case of in (2.18). Let and be the two roots of (2.18), then we have

 ϱ(σj)=1, ¯¯¯ϱ(σj)=fj0−λjfj0+λj<1.

Obviously, for close to , is always true. Because is an increasing function of and is a decreasing function in , the function will increase with the decrease of . Thus, we have

 ϱ(σj−ε)>1,  ϱ(σj+ε)<1,

which implies that has exactly one more eigenvalue larger than than does, and by using the same method as Lemma 1.1 we can prove that the algebraic multiplicity of this eigenvalue is one. Hence, (2.17) is verified.

Now by the argument in the proof of Theorem 1.3 and Corollary 1.12 in [14], we conclude that either meets or meets for some and . It is easy to check that the system (2.1) is reflective. Thus, we can follow the idea in [16]-[17] and use a reflective and periodic extension method exactly same as that in [15] to show that the first alternative must occur. Then, by Lemmas 1.2 and 1.3, we know that the desired results are true. The proof is completed. ∎

## 3 Stability of bifurcating branches

In this section we shall study the stability of steady states bifurcating from by using the asymptotic analysis and perturbation method. We first look for the asymptotic expression of steady states . To proceed, we let

 σ=σ0+∞∑k=1εkσk, (3.1)

where . Then expand and as power series in , that is,

 {ˆu=u∗+∑∞k=1εkuk,ˆv=v∗+∑∞k=1εkvk. (3.2)

Substituting (3.1) and (3.2) into (2.5), we collect the coefficients of and , respectively, and then have the following two systems