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

TL;DR

This paper analyzes the stability and bifurcation structure of steady states in a bacterial colony model with density-dependent motility, revealing conditions for uniform stability and pattern formation.

Contribution

It provides a new analytical framework for understanding how bacterial growth rates influence pattern formation and stability in density-suppressed motility models.

Findings

Large growth rates lead to stable uniform states.

Lower growth rates induce pattern formation.

Analytical results are supported by numerical simulations.

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.

Equations

$$\displaystyle\left\{\begin{array}[]{ccc}\begin{aligned} &u_{t}=\Delta(r(v)u)+\sigma u(1-u),&x\in\Omega,\ t>0,\\ &v_{t}=D\Delta v-v+u,&x\in\Omega,\ t>0\end{aligned}\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{ccc}\begin{aligned} &u_{t}=\nabla\cdot(r(v)\nabla u+ur^{\prime}(v)\nabla v)+\sigma u(1-u),&x\in\Omega,\ t>0,\\ &v_{t}=D\Delta v-v+u,&x\in\Omega,\ t>0,\end{aligned}\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{ccc}\begin{aligned} &\Delta(r(v)u)+\sigma u(1-u)=0,&x\in\Omega,\\ &D\Delta v-v+u=0,&x\in\Omega,\\ &\nabla u\cdot\nu=\nabla v\cdot\nu=0,&x\in\partial\Omega,\end{aligned}\end{array}\right.$$

$$r(v)\in C^2([0,\infty )),r(v)>0\text{and}{r}^{\prime }(v)<0\text{for}v\in [0,\infty ),\underset{v\to +\infty }{\lim }r(v)=0,$$

$$(u(x),v(x))\in \mathbb{B}=\{(u,v):\frac{1}{B}\le u,v\le B\}\text{for}x\in \overline{\Omega }$$

$$0={\lambda }_0<{\lambda }_1<{\lambda }_2<{\lambda }_3<\cdot \cdot \cdot .$$

$$\mathbf{X}=\underset{i=1}{\overset{\infty }{⨁}}{\mathbf{X}}_i,{\mathbf{X}}_i=\underset{j=1}{\overset{\dim S({\lambda }_i)}{⨁}}{\mathbf{X}}_{ij},$$

$$\displaystyle\left\{\begin{array}[]{ccc}\begin{aligned} &\frac{dU}{dt}=r(1)\Delta U+r^{\prime}(1)\Delta V-\sigma U,&x\in\Omega,\ t>0,\\ &\frac{dV}{dt}=D\Delta V-V+U=0,&x\in\Omega,\ t>0,\\ &\nabla U\cdot\nu=\nabla V\cdot\nu=0,&x\in\partial\Omega,t>0.\end{aligned}\end{array}\right.$$

$${\rho }_i^2+\left[(D+r(1)){\lambda }_i+1+\sigma \right]{\rho }_i+[\sigma +r(1){\lambda }_i](1+D{\lambda }_i)+r^{\prime }(1){\lambda }_i=0,i=0,1,2,\cdot \cdot \cdot ,$$

$$\sigma -\left\{-\left[\frac{r^{\prime }(1)}{1+D{\lambda }_i}+r(1)\right]{\lambda }_i\right\}\stackrel{def}{=}\sigma -{\sigma }_i,i=1,2,\cdot \cdot \cdot .$$

$${\sigma }_i>0\text{for}i\in [1,i^c],{\sigma }_{i^c+j}\le 0\text{for}j=1,2,3,\cdot \cdot \cdot \infty .$$

$$r^{\prime }(1)+r(1)<0.$$

$${\sigma }_a=\underset{1\le i\le i^c}{\max }{\sigma }_i=-\left[\frac{r^{\prime }(1)}{1+D{\lambda }_{i_a}}+r(1)\right]{\lambda }_{i_a},i_a\in [1,i^c].$$

$${\lambda }_i=\frac{1}{D}\left(\sqrt{\frac{-r^{\prime }(1)}{r(1)}}-1\right)$$

$${\sigma }_c=\frac{1}{D}{\left(\sqrt{-r^{\prime }(1)}-\sqrt{r(1)}\right)}^2.$$

$$r^{\prime }(1)+r(1)\ge 0$$

$$r^{\prime }(1)+r(1)<0\text{and}\sigma D>-(r^{\prime }(1)+r(1)).$$

$$0<\sigma <{\sigma }_c$$

$$\displaystyle\left\{\begin{array}[]{ccc}\begin{aligned} &\frac{du}{dt}=(r(v)u)^{\prime\prime}+\sigma u(1-u),&x\in(0,l),\ t>0,\\ &\frac{dv}{dt}=Dv^{\prime\prime}-v+u,&x\in(0,l),\ t>0,\\ &u^{\prime}(0)=u^{\prime}(l)=0,\ v^{\prime}(0)=v^{\prime}(l)=0,&t>0,\end{aligned}\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{ccc}\begin{aligned} &(r(v)u)^{\prime\prime}+\sigma u(1-u)=0,&x\in(0,l),\\ &Dv^{\prime\prime}-v+u=0,&x\in(0,l),\\ &u^{\prime}(0)=u^{\prime}(l)=0,\ v^{\prime}(0)=v^{\prime}(l)=0.\end{aligned}\end{array}\right.$$

$$\left\{\begin{array}[]{l}-\varphi^{\prime\prime}(x)=\lambda\varphi(x),\ \quad\quad\quad\quad x\in(0,l),\\[2.84526pt] \varphi^{\prime}(x)=0,\ \quad\quad\quad\quad\quad\quad\quad x=0,l\end{array}\right.$$

$${\lambda }_j=(\pi j/l{)}^2,j=0,1,2,\cdot \cdot \cdot$$

$$\varphi_{j}(x)=\left\{\begin{array}[]{l}1,\ \quad\quad\quad\quad\quad\quad\quad j=0,\\[2.84526pt] \cos(\pi jx/l),\ \quad\quad\quad j>0.\end{array}\right.$$

$$X=\{(u,v):u,v\in C^2([0,l]),u^{\prime }=v^{\prime }=0\mathrm{at}x=0,l\},$$

$$({\omega }_1,{\omega }_2{)}_Y=(u_1,u_2{)}_{L^2(0,l)}+(v_1,v_2{)}_{L^2(0,l)}$$

$$\left\{\begin{array}[]{l}r^{\prime\prime}(v)v^{\prime 2}u+r^{\prime}(v)v^{\prime\prime}u+2r^{\prime}(v)v^{\prime}u^{\prime}+r(v)u^{\prime\prime}+\sigma u(1-u)=0,\ x\in(0,l),\\[2.84526pt] Dv^{\prime\prime}+u-v=0,\ x\in(0,l),\\[2.84526pt] u^{\prime}(0)=u^{\prime}(l)=0,\ v^{\prime}(0)=v^{\prime}(l)=0.\end{array}\right.$$

$$P(\sigma,\omega)=\left(\begin{array}[]{c}r^{\prime\prime}(v)v^{\prime 2}u+r^{\prime}(v)v^{\prime\prime}u+2r^{\prime}(v)v^{\prime}u^{\prime}+r(v)u^{\prime\prime}+\sigma u(1-u)\\ Dv^{\prime\prime}+u-v\end{array}\right),$$

$$P(\sigma ,{\omega }^{\ast })=0\text{for}\sigma >0.$$

$${\lambda }_j<-\frac{r^{\prime }(1)+r(1)}{Dr(1)},\text{and}{\sigma }_j\ne {\sigma }_k\text{for all integers}k\ne j,$$

$$u(\epsilon )=u^{\ast }+\epsilon a_j{\phi }_j+o(\epsilon ),v(\epsilon )=v^{\ast }+\epsilon {\phi }_j+o(\epsilon ),a_j=1+D{\lambda }_j.$$

Full text

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: MXversion28.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

The following nonlinear reaction-diffusion system

[TABLE]

was first introduced in [2] to describe the dynamical behavior of the bacterial species Vibrio Fischeri’s colonies. Here $\Omega$ is a bounded open domain in $\mathbb{R}^{n}$, $n\geq 1$ is a positive integer, and $\Delta=\sum_{i=1}^{n}\frac{\partial^{2}}{\partial x_{i}^{2}}$. The quantities $u(x,t)$ and $v(x,t)$ stand for the density of bacteria and acy-homoserine lactone (AHL) secreted by Vibrio Fischeri, respectively. The positive constants $\sigma$ and $D$ measure, respectively, the logistic growth rate of the bacteria and the diffusion rate of AHL. The diffusion rate of bacteria is state-dependent on $v$ modeled by the positive motility function $r(v)$. As stated in [5], the model (1.2) can be transformed to

[TABLE]

which is a chemotaxis model of Keller-Segel type proposed in [6] if the cells/ bacteria do not sense the concentration between receptors, more details can be found in [6] and [19].

At the present, for the model (1.2), the rigorous mathematical results are very limited. When $\sigma=0$ (namely bacteria have no growth), the existence of global solutions were obtained in [14, 18], the metastability of non-constant steady states was discussed in [16]. For the case where $\sigma>0$, the mechanism of stripe formation of (1.2) was analyzed in [2] when $r(v)$ is a piecewise decreasing function. In [12], authors discussed the pattern solutions and their stability when the diffusion rate of $u$ has a drop at some critical AHL concentration, that is, $r(v)$ is a step function. The *apriori * $L^{\infty}-$ 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 [5] when the system (1.2) is located in a two-dimensional bounded domain with zero Neumann boundary conditions and some conditions are imposed on the motility function $r(v)$. In [10], authors investigated the boundedness, existence and non-existence of non-constant positive classical solutions to the stationary problem of the model (1.2), that is

[TABLE]

where $\nu$ is the outward unit normal vector on $\partial\Omega$ and the domain $\Omega\subset\mathbb{R}^{3}$ is bounded and has smooth boundary. The boundary condition means that there is no flux of either bacteria or AHL across the boundary of the domain. Under the condition that the motility function $r(v)$ satisfies

[TABLE]

in [10] there is the following result:

Lemma 1.1**.**

Let $\Omega$ be a bounded domain in $\mathbb{R}^{n}(1\leq n\leq 3)$ with smooth boundary. Then for any given constant $D_{0}>0$, there exists a positive constant $B>1$, which depends only on $D_{0}$ and $\Omega$, such that any positive solution $(u,v)$ of (1.6) satisfies

[TABLE]

provided that $D\geq D_{0}$. Furthermore, if $\lim\inf_{v\rightarrow\infty}r(v)v\in(r(0),\infty)$, such a constant $B$ is independent of $D_{0}$ and $\Omega$.

Since $r(v)\rightarrow 0$ as $v\rightarrow\infty$, by Lemma 1.1, the case of degeneracy will not happen here. Due to the assumption $r^{\prime}(v)<0$, we call AHL concentration being of the repressive effect on bacterium motility.

In order to further present the preliminaries, we now give some notations. Let $W^{m,p}(\Omega,\mathbb{R}^{N})$ for $m\geq 1,\ 1<p<+\infty$ be the Sobolev space of $\mathbb{R}^{N}$- valued functions with norm $\|\cdot\|_{m,p}$. When $p=2$, $W^{m,2}(\Omega,\mathbb{R}^{N})$ is written as $H^{m}(\Omega)$. Let $L^{p}(\Omega)(1\leq p\leq\infty)$ denote the usual Lebesgue space in a bounded domain $\Omega\subset\mathbb{R}^{n}$ with norm \|f\|_{p}=\Big{(}\int_{\Omega}|f(x)|^{p}dx\Big{)}^{1/p} for $1\leq p<\infty$ and $\|f\|_{\infty}=\mathrm{ess}\sup\limits_{x\in\Omega}|f(x)|$. When $p\in(n,+\infty)$, $W^{1,p}(\Omega,\mathbb{R}^{2})\hookrightarrow C(\Omega,\mathbb{R}^{2})$ which is the space of $\mathbb{R}^{2}$-valued continuous functions.

The following properties of the negative Laplacian operator $-\Delta$ with zero Neumann boundary condition on $\Omega$ will be used later. There is a sequence of eigenvalues ${\lambda}^{\infty}_{i=0}$ satisfying

[TABLE]

Each $\lambda_{i}$ has multiplicity $m_{i}\geq 1$. Let $\varphi_{ij},\ i\geq 0,\ 1\leq j\leq m_{i}$, be the normalized eigenfunctions corresponding to $\lambda_{i}$. Let $S(\lambda_{i})$ be the eigenspace associated with $\lambda_{i}$ in $H^{1}(\Omega;\mathbb{R}^{2})$. Then the set $\{\varphi_{ij},\ i\geq 0,\ j=1,2,\cdots,\mathrm{dim}S(\lambda_{i})\}$ forms a complete orthogonal basis in $L^{2}(\Omega)$. Let $X=[H^{1}(\Omega)]^{2}$ and $X_{ij}=\{c\varphi_{ij}:1\leq j\leq m_{i},c\in\mathbb{R}^{2}\}$. Then

[TABLE]

where $\bigoplus$ denotes the direct sum of subspaces and $\mathrm{dim}S(\lambda_{i})=m_{i}$ .

It is obvious that the system (1.6) has two constant solutions, i.e., $(u(x),v(x))\equiv(0,0)$ and $(u(x),v(x))\equiv(1,1)$ for all $x\in\Omega$. Linearizing (1.2) with Neumann boundary at $(0,0)$ and $(1,1)$ respectively, by a simple computation, we know that $(0,0)$ is always unstable. The linearized system at the point $(1,1)$ reads

[TABLE]

Let $(\varphi,\psi)e^{\rho t}$ be the solution of (1.12). Then the eigenvalue $\rho=\rho(\lambda_{i})\stackrel{{\scriptstyle def}}{{=}}\rho_{i},$ of (1.12) satisfies

[TABLE]

where $\lambda_{i},i=0,1,2,3,\cdot\cdot\cdot$ is defined in (1.9). By the standard stability theory, for the stability/instability of the steady state $(1,1)$ we have a critical discriminant

[TABLE]

It is easy to check that if there is $i$ such that $\sigma_{i}>0$, then there must exist a positive integer $i^{c}$ such that

[TABLE]

Note that (1.15) implies that

[TABLE]

Set

[TABLE]

Moreover, if we regard $\lambda_{i}$ as any real number and use (1.16), then at

[TABLE]

the maximum of $\sigma_{i}$, denoted by $\sigma_{c}$, is attained as

[TABLE]

It is clear that $\sigma_{c}\geq\sigma_{a}$, and that if $i_{a}$ is such that (1.18) is true, then $\sigma_{c}=\sigma_{a}$. We now have the lemma below.

Lemma 1.2**.**

Suppose that (1.7) and (1.15) hold. Then, for (1.2) with the zero Neumann boundary condition we have the following facts:

$(i)$* The steady state $\omega*=(1,1)$ is linearly stable if either*

[TABLE]

or

[TABLE]

$(ii)$* $\omega*$ is unstable if $0<\sigma<\sigma_{a}$; Usually, we call $k_{a}=\sqrt{\lambda_{i_{a}}}$ admissible wave number.*

$(iii)$* $\omega*$ is linearly stable if $\sigma>\sigma_{c}$.*

Naturally, we may expect the existence of non-constant steady states as the constant solutions are unstable and figure out their structure. The purpose of this paper is to establish the existence and structure of positive solutions of (1.6) in one dimensional space $\Omega=(0,l),l>0$ and to derive the criteria for the stability/unstability of each bifurcation branch.

Throughout this paper, by Lemmas 1.1 and 1.2, we assume that both (1.7) and (1.15) are always true, and that $\sigma$ satisfies

[TABLE]

for fixed constants $D$ and $l$.

This paper is organized as follows. In Section 2, we discuss the local and global bifurcation to describe the structure of positive solutions near the bifurcation points and prove that these bifurcation curves can be prolonged as long as the parameter $\sigma$ is less than the critical value $\sigma_{c}$. In Section 3, we use the asymptotic analysis and the adjoint theory to derive the expression of the steady states. Then the stability/unstability criteria of the bifurcating solutions are given. Numerical simulations are carried out to demonstrate all the theoretical results in Section 4.

2 Local and global bifurcation

With $\Omega=(0,l),l>0$ the system (1.2) with Neumann boundary conditions can be rewritten as

[TABLE]

whose stationary system is (1.6) with one dimensional space, i.e.,

[TABLE]

We know that the eigenvalue problem

[TABLE]

has a sequence of simple eigenvalues

[TABLE]

and their corresponding eigenfunctions are

[TABLE]

Obviously, the set of eigenfunctions constitutes an orthogonal basis in $L^{2}(0,l)$. Let

[TABLE]

then $X$ is a Banach space with the usual $C^{2}$ norm, and $Y=L^{2}(0,l)\times L^{2}(0,l)$ is a Hilbert space with the inner product

[TABLE]

for $\omega_{1}=(u_{1},v_{1})\in Y$, $\omega_{2}=(u_{2},v_{2})\in Y$. By expanding the second-order derivative term in the first equation, we have the system (2.4) in the form of

[TABLE]

Define the map $P:\Lambda\longrightarrow Y$ by

[TABLE]

where $\omega=(u,v)$, and $\Lambda=(0,\sigma_{a})\times\mathbb{B}$ is a bounded set in $(0,\infty)\times X$. Hence, looking for the solutions of (2.4) is exactly equivalent to looking for the zero points of this map. Let $\omega^{*}=(u^{*},v^{*})=(1,1)$, then we have

[TABLE]

We recall that, for a number $\alpha>0$, $(\alpha,\omega^{*})$ is a bifurcation point of the equation $P=0$ with respect to the curve $(\sigma,\omega^{*}),\sigma>0$ if every neighborhood of $(\alpha,\omega^{*})$ contains zeros of $P$ in $(0,\infty)\times X$ not lying on this curve. Then the results on local bifurcation of solutions for (2.4) are as follows.

Theorem 2.1**.**

Suppose that (1.20) is true. If $j$ is a positive integer such that

[TABLE]

then $(\sigma_{j},\omega^{*})$ is a bifurcation point of $P=0$ with respect to the curve $(\sigma,\omega^{*}),\sigma>0$, where $\sigma_{j}$ is defined in (1.14). Furthermore, there is a one-parameter family of non-trivial solutions $\Gamma_{j}(\varepsilon)=(\sigma(\varepsilon),u(\varepsilon),v(\varepsilon))$ of the problem (2.4) for $|\varepsilon|$ sufficiently small, where $\sigma(\varepsilon),u(\varepsilon),v(\varepsilon)$ are continuous functions, $\sigma(0)=\sigma_{j}$ and

[TABLE]

The set of zero-points of $P$ consists of two curves $(\sigma,\omega^{*})$ and $\Gamma_{j}(\varepsilon)$ in a neighborhood of the bifurcation point $(\sigma_{j},\omega^{*})$.

Proof.

Fix $j$, according to Theorem 1.7 of [1], we need to verify the following conditions:

$(1)$   the partial derivatives $P_{\sigma},P_{\omega}$, and $P_{\sigma\omega}$ exist and are continuous,

$(2)$   $\ker P_{\omega}(\sigma_{j},\omega^{*})$ and $Y/R(P_{\omega}(\sigma_{j},\omega^{*}))$ are one-dimensional,

$(3)$   let $\ker P_{\omega}(\sigma_{j},\omega^{*})=span\{\varphi\}$, then $P_{\sigma\omega}(\sigma_{j},\omega^{*})\varphi\notin R(P_{\omega}(\sigma_{j},\omega^{*}))$.

Because we have

[TABLE]

and

[TABLE]

it is clear that the linear operators $P_{\sigma},P_{\omega}$, and $P_{\sigma\omega}$ are continuous. Condition $(1)$ is verified.

Let $\Phi=(\varphi,\psi)\in\ker L$ with $\varphi=\sum_{0\leq i\leq\infty,}{a_{i}\varphi_{i}}$ and $\psi=\sum_{0\leq i\leq\infty,}{b_{i}\varphi_{i}}$. Then we have

[TABLE]

which means that, by the definition of $\varphi_{i}$ in (2.5), all the coefficients must vanish, that is,

[TABLE]

This equation has a nonzero solution provided that

[TABLE]

which holds if and only if

[TABLE]

Obviously, if $i=0$, then $\sigma=0$, which is excluded by the assumption of the theorem. In view of (2.9), the equation (2.13) holds only for $i=j$ and

[TABLE]

Then, by solving the equation (2.12) with $i=j$, we have

[TABLE]

where $a_{j}=1+D\lambda_{j}$ and $\varphi_{j}$ is defined as in (2.7). Moreover, the adjoint operator of $L$ is

[TABLE]

which has

[TABLE]

where $a^{*}_{j}=-\frac{1+D\lambda_{j}}{r^{\prime}(1)\lambda_{j}}$. It is well known that $R(L)=(\ker L^{*})^{\bot}$. Then the codimension of $R(L)$ is the same as $\dim\ker L^{*}=1$. Condition $(2)$ is thus satisfied.

We know that

[TABLE]

and

[TABLE]

Hence, $P_{\sigma\omega}(\sigma_{j},\omega^{*})\Phi\notin R(L)$, and so condition $(3)$ is verified. In addition, (2.14) implies that (2.10) is true. The proof is completed. ∎

Remark 2.1**.**

(a)*  Theorem 2.1 and Lemma 1.2 show that if $\sigma_{j}\in(0,\sigma_{c})$, then $(\sigma_{j},\omega^{*})$ is a bifurcation point with respect to the trivial branch $(\sigma,\omega^{*})$. For given $D$ and $l$, the number of such bifurcation points is equal to the cardinality of indices $j$ such that $\sigma_{j}\in(0,\sigma_{c})$, see (1.15). (b)   By Theorem 2.1, let $\Upsilon$ be the closure of the non-trivial solution set of $P=0$, then $\Gamma_{j}$ is the connected component of $\Upsilon\cup\{(\sigma_{j},\omega^{*})\}$ to which $(\sigma_{j},\omega^{*})$ belongs, and in a neighborhood of the bifurcation point the curve $\Gamma_{j}$ is characterized by the eigenfunction $\varphi_{j}$. In the open interval $(0,l)$ the function $\varphi_{j}$ has exactly $j$ zeros, thus we call the non-constant solutions in $\Gamma_{j}$ mode $j$ steady states.*

We next apply the global bifurcation theory of Rabinowitz, in particular, Corollary 1.12 in [13], and the Leray-Schauder degree for compact operators to give the information on the bifurcating curve $\Gamma_{j}$ far from the trivial equilibrium. Following the idea of [9], we first rewrite the system (2.8) as

[TABLE]

where

[TABLE]

Now shifting the constant state $\omega^{*}=(1,1)$ to $\mathcal{O}=(0,0)$ by setting $\tilde{u}=u-1$,$\tilde{v}=v-1$, then (2.15) is transformed into

[TABLE]

where $\tilde{f}$ and $\tilde{g}$ are higher-order terms of $\tilde{u}$ and $\tilde{v}$,

[TABLE]

and

[TABLE]

Note that $g_{1}<0$, and $f_{0}>0$ since $\sigma<\sigma_{c}<\frac{-r^{\prime}(1)}{D}$, see (1.19). Let $F_{\sigma}$ and $F$ be the inverse operator of $f_{0}-\frac{d^{2}}{dx^{2}}$ and $-g_{1}-\frac{d^{2}}{dx^{2}}$, respectively; moreover, set

[TABLE]

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

[TABLE]

Obviously, for any given $\sigma>0$ the linear operator $M(\sigma)$ is compact on $X$. On closed $\sigma$ sub-intervals of $(0,\infty)$ the operator $H(\sigma,\widetilde{\omega})$ is compact on $X$ and $H(\sigma,\widetilde{\omega})=o(\|\widetilde{\omega}\|)$ for $\widetilde{\omega}$ near zero uniformly.

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

Lemma 2.2**.**

Assume that (1.20) and (2.9) are satisfied. Then $1$ is an eigenvalue of $M(\sigma_{j})$ with algebraic multiplicity one.

Proof.

Let $\Phi=(\varphi,\psi)$,  $\varphi=\sum_{i=0}^{\infty}a_{i}\varphi_{i},\ \psi=\sum_{i=0}^{\infty}b_{i}\varphi_{i}$. The proof of Theorem 2.1 implies that

[TABLE]

has one unique solution \Phi=\left(\begin{array}[]{c}1+D\lambda_{j}\\ 1\end{array}\right)\varphi_{j}, which shows that $1$ is an eigenvalue of $M(\sigma_{j})$ with the unique eigenfunction. Thus, we have $\dim\ker(M(\sigma^{j})-I)=1$. Next we will prove that the eigenvalue $1$ is simple. It is well known that the algebraic multiplicity of $1$ is equal to the dimension of the generalized null space $\bigcup_{i=1}^{\infty}\ker(M(\sigma_{j})-I)^{i}$. So we need only to verify $\ker\left(M(\sigma_{j})-I\right)\cap R\left(M(\sigma_{j})-I\right)=\{0\}$. Let $M^{*}(\sigma_{j})$ be the adjoint operator of $M(\sigma_{j})$. If $(\varphi,\psi)\in\ker(M^{*}(\sigma_{j})-I)$, from (2.17) it follows that

[TABLE]

where

[TABLE]

By the definition of $F_{\sigma}$ and $F$, the system (2.18) can be expanded as

[TABLE]

with

[TABLE]

Again set $\varphi=\sum_{i=0}^{\infty}a_{i}\varphi_{i},\ \psi=\sum_{i=0}^{\infty}b_{i}\varphi_{i}$. By (2.19), we get

[TABLE]

By $\sigma\neq 0$ and (2.9), we know that $\det A_{i}^{*}=0$ if and only if $i=j$ and

[TABLE]

Therefore, $\ker(M^{*}(\sigma_{j})-I)$ is generated by the unique element \Phi^{*}=\left(\begin{array}[]{c}f_{0}^{j}+\lambda_{j}\\ f_{1}^{j}\end{array}\right)\phi_{j}, and we know that $(M(\sigma_{j})-I)\Phi=0$ has one unique solution (up to a constant multiple ) \Phi=\left(\begin{array}[]{c}-f_{1}^{j}\\ f_{0}^{j}\lambda_{j}\end{array}\right)\varphi_{j}, and thus $(\Phi,\Phi^{*})=-2f_{1}^{j}\lambda_{j}\neq 0$. Then $\Phi\notin(Ker(M^{*}(\sigma_{j})-I))^{\perp}=R(M(\sigma_{j})-I)$, which implies that $\ker(M(\sigma_{j})-I)\cap R(M(\sigma_{j})-I)=\{0\}$. Thus, the desired result follows. ∎

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

Theorem 2.3**.**

Suppose that (1.20) and (2.9) are true. Then the projection of the bifurcation curve $\Gamma_{j}$ onto the $\sigma-$axis is an interval $(0,\sigma_{j})$. Furthermore, the system (2.4) has at least one non-constant positive solution if $\sigma\in(0,\sigma_{a})$ and $\sigma\neq\sigma_{k}$ for any positive integer $k$.

Proof.

By the proof of Lemma 2.2, we know that the linear operator $I-M(\sigma):X\rightarrow X$ is a bijection when $\sigma\in(0,\sigma_{a})\setminus\sigma_{j}$ and being located in a small neighborhood of $\sigma_{j}$. For this fixed $\sigma$, let $\mathcal{O}$ be an isolated solution of (2.17). The index of this isolated zero of the map $I-\left(T(\sigma,.\right)$ is given by

[TABLE]

where $\mathcal{B}$ is a sufficiently small ball centered at $\mathcal{O}$, and $\beta$ is the sum of the algebraic multiplicities of the eigenvalues of $M(\sigma)$ that are larger than $1$. For our bifurcation analysis, we necessarily verify that this index changes as the bifurcation parameter $\sigma$ crosses $\sigma_{j}$, i.e., for $\varepsilon>0$ sufficiently small,

[TABLE]

Indeed, if $\varrho$ is an eigenvalue of $M(\sigma)$ corresponding to an eigenfunction $(\varphi,\psi)$, then we have

[TABLE]

Once again let $\varphi=\sum_{i=0}^{\infty}a_{i}\varphi_{i}$ and $\psi=\sum_{i=0}^{\infty}b_{i}\varphi_{i}$, then the above system can be expanded as

[TABLE]

Then the set of eigenvalues of $M(\sigma)$ consists of all $\varrho^{\prime}s$ that solve the characteristic equation

[TABLE]

Taking $\sigma=\sigma_{j}$, if $\varrho=1$ is a root of (2.21), then it is concluded that

[TABLE]

and thus by the assumption (1.20) we have $i=j$. Therefore, if we do not count the eigenvalues corresponding to $i=j$ in (2.21), $M(\sigma)$ has the same number of eigenvalues which are larger than $1$ for all $\sigma$ close to $\sigma_{j}$, and have the same multiplicities. So we need only to consider the case of $i=j$ in (2.21). Let $\varrho(\sigma)$ and $\overline{\varrho}(\sigma)$ be the two roots of (2.21), then we have

[TABLE]

Obviously, for $\sigma$ close to $\sigma_{j}$, $\overline{\varrho}(\sigma)<1$ is always true. Because $\varrho(\sigma)$ is an increasing function of $f_{0}^{j}$ and $f_{0}^{j}$ is a decreasing function in $\sigma$, the function $\varrho(\sigma)$ will increase with the decrease of $\sigma$. Thus, we have

[TABLE]

which implies that $M(\sigma_{j}-\varepsilon)$ has exactly one more eigenvalue larger than $1\,,$ than $M(\sigma_{j}+\varepsilon)$ does, and by using the same method as Lemma 2.2 we can prove that the algebraic multiplicity of this eigenvalue is one. Hence, (2.20) is verified.

Now, by (2.20) and Corollary 1.12 in [13], we conclude that $\Gamma_{j}$ either meets $\partial\Lambda$ or meets $(\sigma_{k},0)$ for some $k\neq j$ and $\sigma_{k}>0$. It is easy to check that the system (2.4) is reflective. Thus, we can follow the idea in [15, 17] and use a reflective and periodic extension method, which is also exactly the same as that in [9], to show that the first alternative must occur. Then, by Lemmas 1.1 and 1.2, 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 $(\widehat{u}(x),\widehat{v}(x))$ bifurcating from $\omega^{*}=(u^{*},v^{*})=(1,1)$ by using the asymptotic analysis and perturbation method. We first look for the asymptotic expression of steady states $(\widehat{u}(x),\widehat{v}(x))$. To proceed, we let

[TABLE]

where $\sigma_{0}$ needs to be ascertained afterwards, $0<\epsilon\leq 1$. Then expand $\widehat{u}(x)$ and $\widehat{v}(x)$ as power series in $\epsilon$, that is,

[TABLE]

Substituting (3.1) and (3.2) into (2.8), and expanding $r(v),r^{\prime}(v)$ and $r^{\prime\prime}(v)$ as Taylor expansion at $(1,1)$, we collect the coefficients of $O(\epsilon)$ and $O(\epsilon^{2})$, respectively, and then have the following two systems

[TABLE]

and

[TABLE]

where

[TABLE]

Directly solving the system (3.3) yields a unique non-constant solution (up to a constant multiple ) for some integer $j\in\{1,2,\cdot\cdot\cdot,i^{c}\}$

[TABLE]

as long as $\sigma_{0}$ is equal to

[TABLE]

where $(\lambda_{j},\varphi_{j})$ is given by (2.6) and (2.7). Here the uniqueness of solution indicates that $\sigma_{0}^{j}\neq\sigma_{0}^{k}$ for any integer $k\neq j$. In fact, here $\sigma_{0}^{j}$ is just $\sigma_{j}$ in Section 2. Then we have

[TABLE]

for some positive integer $i_{a}$, and here $\sigma_{0}^{i_{a}}$ is exactly $\sigma_{a}$ in Section 2. It is observed that $i_{a}$ is the wave mode maximizing $\sigma_{0}^{j}$ such that $\sigma_{0}^{i_{a}}$ is the maximum bifurcation value. We shall call $i_{a}$ the * admissible wave mode* corresponding to the admissible wave number in Section 2.

In order to solve (3.4), we consider its adjoint system

[TABLE]

This system has a non-constant solution

[TABLE]

Here $\varphi_{j}$ is the same as in (3.5). By the Fredholm alternative [7], the equation (3.4) admits a solution if and only if

[TABLE]

Solving this equation yields

[TABLE]

Then $G_{1}$ in (3.4) can be simplified to

[TABLE]

By this, we can set a particular solution of (3.4) as

[TABLE]

Substitution of (3.10) into (3.4) leads to

[TABLE]

On account of $\sigma_{1}=0$, we need to find the expression of $\sigma_{2}$. Again substituting (3.1) and (3.2) into (2.8) and equating the coefficients of $O(\epsilon^{3})$, we have

[TABLE]

where

[TABLE]

Then applying the solvability condition $\int_{0}^{l}\overline{u}_{2}G_{2}dx=0$ of (3.12) yields

[TABLE]

We now know that when the parameter $\sigma$ given by (3.1) lies in the neighborhood of $\sigma_{0}^{j}$, the bifurcating solution $(\widehat{u},\widehat{v})$ is described by (3.2) with $(u_{1},v_{1})$ and $(u_{2},v_{2})$ formulated by (3.5) and (3.10), respectively. To bring out the relationship between the solution $(\widehat{u},\widehat{v})$ and its bifurcation location $\sigma_{0}^{j}$, we denote $(\widehat{u},\widehat{v})$ by $(\widehat{u}_{j},\widehat{v}_{j})$, that is,

[TABLE]

Normally, $\omega^{*}=(u^{*},v^{*})$ is called the* base* term of the non-constant steady state $(\widehat{u},\widehat{v})$ whose shape and amplitude primarily depend on the leading term $(u_{1},v_{1})$ when $\epsilon$ is small, that is,

[TABLE]

which shows the maximum change of the bacterial density from the base term. Taking into account (2.5) and (3.5), the leading term has wave mode $j$. Thus, $j$ is called the principal wave mode of the solution $(\widehat{u}_{j},\widehat{v}_{j})$.

We shall analyze the stability of the solution (3.15) located at the $j$th bifurcating branch by discussing the sign of the principal eigenvalue of linearized system of (2.8) around $(\widehat{u}_{j},\widehat{v}_{j})$. Let

[TABLE]

and substitute it into (2.8), and expand $r(v),r^{\prime}(v)$ and $r^{\prime\prime}(v)$ as Taylor expansion at $(\widehat{u}_{j},\widehat{v}_{j})$. Then the linearized system of (2.8) is

[TABLE]

where

[TABLE]

Moreover, set

[TABLE]

and

[TABLE]

Substituting these equations into (3.17) and equating the coefficients of $O(1)$ lead to the following system

[TABLE]

By (1.10), we replace $(\varphi_{0}^{\prime\prime},\psi_{0}^{\prime\prime})$ by $-\lambda_{m}(\varphi_{0},\psi_{0})$. Then the existence of a non-zero solution $(\varphi_{0},\psi_{0})$ yields the following equation

[TABLE]

where

[TABLE]

where $\sigma_{0}^{j}$ is given by (3.6). By (3.7), if $j\neq i_{a}$, when the positive integer $m=i_{a}$ such that $E<0$, then the equation (3.20) has a positive root $\gamma_{0}>0$ which implies that $(\widehat{u}_{j},\widehat{v}_{j})$ is unstable. So we have a conclusion as follows:

Proposition 3.1**.**

The non-constant steady state $(\widehat{u}_{j},\widehat{v}_{j})$ in (3.15) is unstable when $j\neq i_{a}$. In other words, if $(\widehat{u}_{j},\widehat{v}_{j})$ is stable, then it is necessary that $j=i_{a}$.

We shall derive a sufficient condition for the stability of the non-constant steady states with the *admissible * wave mode $i_{a}$. Through a simple calculation, we find that the principal eigenvalue of (3.19) is $\gamma_{0}=0$ with the eigenfunction

[TABLE]

Next, we compute $\gamma_{1}$. Again carrying out the computation of obtaining (3.19) and equating the $O(\varepsilon)$ terms lead to

[TABLE]

where

[TABLE]

Applying the solvability condition of (3.21), we have

[TABLE]

where $(\overline{u}_{2},\overline{v}_{2})$ is given by (3.9) with $j=i_{a}$. Solving (3.22) for $\gamma_{1}$, we have

[TABLE]

A direct calculation yields

[TABLE]

and

[TABLE]

Due to $\gamma_{1}=0$, we need further to compute $\gamma_{2}$. We first simplify $G_{3}$ as

[TABLE]

By this, a particular solution of (3.21) ($\varphi_{1},\psi_{1})$ is of the following form

[TABLE]

with

[TABLE]

where $d_{i}(i_{a})$ is given by (3.11). Again we use the same computation of obtaining (3.19), but now equate the $O(\epsilon^{2})$ terms, and then get the following system including $\gamma_{2}$

[TABLE]

where

[TABLE]

Again using the solvability condition of (3.24), we have

[TABLE]

and

[TABLE]

where

[TABLE]

In view of (3.23), we know that the stability of $(\widehat{u}_{i_{a}},\widehat{v}_{i_{a}})$ completely depends on the sign of $\eta$, and thus we have the result below.

Theorem 3.2**.**

If the positive integer $i_{a}$ is the admissible wave mode, then the small-amplitude steady state $(\widehat{u}_{i_{a}},\widehat{v}_{i_{a}})$ of system $(\ref{model223})$ is stable provided that

[TABLE]

4 Numerical simulation

This section is devoted to presenting some numerical examples to demonstrate the theoretical results obtained in Sections 2 and 3. The model is solved with the MATLAB pde solver based on the finite difference scheme. For the sake of brevity, only the numerical results of the solution component $u$ are presented here. We take the length of the spatial interval $l=20$, the diffusion coefficient of AHL $D=1$ and the motility function

[TABLE]

which obviously satisfies the condition (1.7). The small parameter $\epsilon=0.01$ is always fixed in this section.

In Section 2, Theorem 2.3 shows that each bifurcation $\Gamma_{j}$ emanating from $(\sigma_{0}^{j},\omega^{*}),j=1,2,\cdot\cdot\cdot,i^{c}$ goes backward and meets the vertical axis (i.e., $\sigma=0$), but we do not know whether $\Gamma_{j}$ directly joins with $\sigma=0$ or meets some bifurcation points and then reaches $\sigma=0$; furthermore, for $\sigma=\sigma_{0}^{j}$ the existence of non-constant steady states is not established in our theorem. Therefore, we can only give the local bifurcation diagram of this example by together with Proposition 3.1 and Theorem 3.2. Through a computation, we have $i^{c}=11$, $\sigma_{c}=0.5$, $i_{a}=6$ such that $\sigma_{max}=\sigma_{0}^{6}=0.4967<\sigma_{c}$, and $\eta=10.3042$; Moreover,

[TABLE]

Therefore, all the bifurcations are backward which clarify the results in Theorem 2.3. By Proposition 3.1 and Theorem 3.2, we have that the sixth bifurcating branch is stable and the remaining ten ones are unstable. Furthermore, all the bifurcation values can be put in order as

[TABLE]

Based on the above discussion and the equation $(\ref{amplitude})$, the bifurcation diagram are presented in Figure 1. In order to make the diagram look cleaner, we only depict the indication of six branches, that is, the bifurcation parameter $\sigma$ is close to the following six bifurcation points

[TABLE]

accordingly,

[TABLE]

Next we numerically verify the asymptotic expressions of patterns in (3.15) and the selection mechanism of stable mode of steady states with small amplitudes (i.e., pattern solutions with small amplitudes) established in Proposition 3.1 and Theorem 3.2. By the analytical results, the stable wave mode is $j=i_{a}=6$, and thus the wave number of stationary pattern is $3$, i.e., there are $3$ peaks; moreover, the second-order approximation of the pattern solution in (3.15) is specified by

[TABLE]

which is plotted in $(a)$ of Figure 2.

In Figure 2, we show the comparison between the stationary state (4.3) predicted by the asymptotical analysis and the stationary states reached starting from different initial functions computed by numerically solving the system (2.2). It is seen that all the three simulation solutions $(b)-(d)$ are qualitatively in very good agreement with the analytical solution $(a)$, that is, all of them have three peaks (each of two boundary peaks evidently only counts as a half peak). The quantitative discrepancy is caused by ignoring the higher order terms of (4.3). It is also demonstrated that the bifurcation branch emanating from the maximum bifurcation point (i.e., the sixth one here) is stable since all the three simulation solutions $(b)-(d)$ tends to the steady state with wave mode $j=6$ after about running time $t=56,1200$ and $350$, respectively. This verifies the stability criterion established in Proposition 3.1 and Theorem 3.2.

In Figure 3, we show the evolution of pattern starting from different initial locations. it is also observed that the sixth bifurcation branch is stable which is precisely predicted by our analytical results in Proposition 3.1 and Theorem 3.2. We see that solutions starting from different perturbation of the constant state $(1,1)$ finally reach the stationary pattern (i.e, the sixth bifurcation branch) which possesses 3 wave peaks (i.e., the wave mode $j=6$).

Figure 4 displays the process of pattern formation from the initial state to the stable steady state of Figure 3 (iii). The initial state has $2$ wave peaks (i.e, the wave mode is $j=4$), see the left-most one in the top line. Since the mode $4$ steady state is unstable, as shown in the middle one of the top line, it is followed by the pattern having $4$ wave peaks (i.e, $j=8$) at about $t=60$. Again due to the unstability of pattern with mode $j=8$ predicted by Proposition 3.1, as we see from the left-most one in the bottom line, at about $t=660$ merging phenomenon appears, then this solution transits to a stable state with mode $j=6$ by reducing the $3$ wave peaks in the interior of the domain to $2$ ones. This transition is complete at about $t=750$, and then the solution keeps the state with mode $j=6$ that is stable. This once again confirms the stability criterion stated in Proposition 3.1 and Theorem 3.2.

5 Conclusion

In this work, we present the detailed information of steady states near the trivial equilibrium and the global structure of the steady states of (2.2). The second-order approximation of non-constant steady states is also derived. Then by using the standard linear stability analysis and analytical technique, we further establish the stability criterion of the bifurcation branches. It is shown that the growth rate of bacteria $\sigma$ substantially influences the dynamical behavior of the model (1.2). Particularly, under conditions (1.7) and (1.15) on the motility function $r(v)$, when the growth rate of bacteria is sufficiently large (i.e, for $\sigma>\sigma_{c}$), the system will keep stabilization around the uniform state $(1,1)$; while for $\sigma\in(0,\sigma_{a})$, the pattern formation must occur, and thus the densities of bacteria and AHL always depend on their location. The principal wave mode of stationary pattern coincides with the admissible one which maximizes bifurcation values. The analytical results are corroborated by direct simulations of the underlying system (2.2) through different stages.

There are various interesting questions arising from our present analytical and numerical studies. The existence of non-constant steady states and the propagation of pattern in a large domain have been rigorously discussed and organized as a separate paper. It should be noticed that the question whether the pattern formation occurs in the cases when $\sigma\in[\sigma_{i_{a}},\sigma_{c}]$ or $\sigma=0$ has not been discussed yet. This question may be investigated by applying the approaches similar to that in [4, 11]. The global attractivity of non-constant steady state with the admissible wave mode still remains open. The bifurcating and emerging process in the pattern formation are numerically presented in Figures 3 and 4, but the mathematical behavior of the merging process is not well understood yet. All these questions are very interesting and challenging. We think it is worthwhile to explore theories and methods of solving them in the future.

Acknowledgement. This work was supported by the National Natural Science Foundation of China (No. 11671359) and the joint mobility project of the National Natural Science Foundation of China and Academy of Finland (No. 11811530145). M. Ma gratefully acknowledges the Department of Mathematics and Statistics at University of Turku for its hospitality and the help of its members during her stay there and the financial support of Academy of Finland. M. Vuorinen’s research at the Zhejiang Sci-tech University in April 2018 was supported by the National Natural Science Foundation of China (No. 11811530145) of Prof M. Ma.