Stiff-response-induced instability for chemotactic bacteria and flux-limited Keller-Segel equation

TL;DR

This paper introduces a new instability mechanism in chemotactic bacteria due to stiff responses, leading to pattern formation, supported by theoretical proofs and numerical simulations aligning with biological observations.

Contribution

It presents a novel instability mechanism based on stiff chemotactic responses, with proofs for both microscopic and macroscopic models, and demonstrates pattern formation through simulations.

Findings

Instability occurs in both microscopic and macroscopic models.

Unstable frequencies are bounded, similar to Turing instability.

Numerical simulations show formation of periodic patterns.

Abstract

Collective motion of chemotactic bacteria as E. Coli relies, at the individual level, on a continuous reorientation by runs and tumbles. It has been established that the length of run is decided by a stiff response to a temporal sensingof chemical cues along the pathway.We describe a novel mechanism for pattern formation stemming from the stiffness of chemotactic response relying on a kinetic chemotaxis model which includes a recently discovered formalism for the bacterial chemotaxis. We prove instability both for amicroscopic description in the space-velocity space and for the macroscopic equation, a flux-limited Keller-Segel equation, which has attracted much attention recently.A remarkable property is that the unstable frequencies remain bounded, as it is the case in Turing instability. Numerical illustrations based on a powerful Monte Carlo method show that the stationary homogeneous state of population density isdestabilized and periodic patterns are generated in realistic ranges of parameters. These theoretical developments are in accordance with several biological observations.

Tables

Table 1. Table 1: The parameter sets with which the Monte Carlo simulations are performed. The parameter sets A, B, C, and D corresponds to the marks A, B, C, and D in Fig.  1 , respectively. The values of the round brackets show the values of d / k 𝑑 𝑘 d/k , χ / k 𝜒 𝑘 \chi/\sqrt{k} , and k ​ δ 𝑘 𝛿 \sqrt{k}\delta , i.e., ( d / k 𝑑 𝑘 d/k , χ / k 𝜒 𝑘 \chi/\sqrt{k} , k ​ δ 𝑘 𝛿 \sqrt{k}\delta ). The black squares satisfy the kinetic instability condition Eq. ( 7 ) while white squares do not satisfy it.

$k$	A(1, 0.5, 0.05)	B(1, 0.5, 0.0625)	C(0.7, 0.5, 0.0625)	D(1, 0.5, 0.1)
0.1	–	$■$	–	–
1.0	$■$	$\square$	$■$	$\square$
2.0	$\square$	–	$\square$	–

Table 2. Table 2: Experimental values for mean tumbling frequency ψ 0 subscript 𝜓 0 \psi_{0} , doubling time t 0 subscript 𝑡 0 t_{0} , running speed V 0 subscript 𝑉 0 V_{0} , diffusion coefficient of chemoattractant d ~ ~ 𝑑 \tilde{d} , and stiffness δ ~ − 1 superscript ~ 𝛿 1 \tilde{\delta}^{-1} and modulation χ 𝜒 \chi in chemotactic response obtained in Ref.  [ 34 ] .

${\psi }_0$	3.0 [s]
$t_0$	4500 [s]
$V_0$	25 [$\mu$m/s]
$\tilde{d}$	8$\times {10}^{-6}$ [cm2/s]
${\tilde{\delta }}^{-1}$	20 [s]
$\chi$	0.2

Equations

$$\begin{array}{rl}& {\partial }_{t}f(t,\mathbf{x},\mathbf{v})+\mathbf{v}\cdot \nabla f=\\ & \frac{1}{k}\left\{\frac{1}{4\pi }{\int }_{V}K[D_t\log S{|}_{{\mathbf{v}}^{\prime }}]f({\mathbf{v}}^{\prime })d\Omega ({\mathbf{v}}^{\prime })-K[D_t\log S{|}_{\mathbf{v}}]f(\mathbf{v})\right\}+P[\rho ]f(\mathbf{v}).\end{array}$$

$$-d\Delta S(t,\mathbf{x})+S(t,\mathbf{x})=\rho (t,\mathbf{x}),$$

$$\rho (t,\mathbf{x})=\frac{1}{4\pi }{\int }_{V}f(t,\mathbf{x},\mathbf{v})d\Omega (\mathbf{v}),$$

$$P[\rho ]>0,\text{mbox}for0<\rho <1,$$

$$P[\rho ]<0,\text{mbox}for\rho >1,$$

$$P[\rho ]\simeq 1-\rho ,\text{mbox}for\rho \simeq 1.$$

$$K[X]=1-F[X],$$

$$F[0]=0,$$

$$\frac{dF[X]}{dX}>0,$$

$$F[X]\to \pm \chi \text{mbox}asX\to \pm \infty .$$

$$\frac{F^{\prime }[0]}{k}>\underset{\lambda }{\inf }\left[1+\frac{k}{\frac{k\lambda }{\arctan (k\lambda )}-1}\right](1+d{\lambda }^2).$$

$$f(t,\mathbf{x},\mathbf{v})$$

$$S(t,\mathbf{x})$$

$$\rho (t,\mathbf{x})$$

$$P[1+{\rho }_g{e}^{\mu t}]=-{\rho }_g{e}^{\mu t},$$

$$K[D_t\log (1+S_g{e}^{\mu t}){|}_{\mathbf{v}}]=1-F^{\prime }[0](\mu S_g+\mathbf{v}\cdot \nabla S_g)e^{\mu t}.$$

$$k(\mu g(\mathbf{x},\mathbf{v})$$

$$=\frac{1}{4\pi }{\int }_V(1-F^{\prime }[0](\mu S_g(\mathbf{x})+{\mathbf{v}}^{\prime }\cdot \nabla S_g(\mathbf{x}))e^{\mu t}+g(\mathbf{x},{\mathbf{v}}^{\prime })e^{\mu t})d\Omega ({\mathbf{v}}^{\prime })$$

$$-(1-F^{\prime }[0](\mu S_g(\mathbf{x})+\mathbf{v}\cdot \nabla S_g(\mathbf{x}))e^{\mu t}+g(\mathbf{x},\mathbf{v})e^{\mu t})-k{\rho }_g(\mathbf{x}),$$

$$k\mu g(\mathbf{x},\mathbf{v})$$

$$=\frac{1}{4\pi }{\int }_{V}g(\mathbf{x},{\mathbf{v}}^{\prime })d\Omega ({\mathbf{v}}^{\prime })+F^{\prime }[0]\mathbf{v}\cdot \nabla S_g(\mathbf{x})-g(\mathbf{x},\mathbf{v})-k{\rho }_g(\mathbf{x}).$$

$$\hat{g}(\mathbf{\lambda },\mathbf{v})=\frac{1-k+\mathrm{i}\frac{F^{\prime }[0]}{1+d|\mathbf{\lambda }{|}^2}\mathbf{\lambda }\cdot \mathbf{v}}{1+k\mu +\mathrm{i}k\mathbf{\lambda }\cdot \mathbf{v}}{\hat{\rho }}_g(\mathbf{\lambda }).$$

$${\hat{\rho }}_g(\mathbf{\lambda })=\frac{1}{2}{\int }_{-1}^1\frac{1-k+\mathrm{i}\frac{F^{\prime }[0]\lambda v}{1+d{\lambda }^2}}{1+k\mu +\mathrm{i}k\lambda v}dv{\hat{\rho }}_g(\mathbf{\lambda }),$$

$${\hat{\rho }}_g(\mathbf{\lambda })=\frac{1}{2}{\int }_{-1}^1\frac{\left(1-k+\mathrm{i}\frac{F^{\prime }[0]\lambda v}{1+d{\lambda }^2}\right)\left(1+k{\mu }_1-\mathrm{i}k\lambda ({\mu }_2+v)\right)}{(1+k{\mu }_1{)}^2+k^2{\lambda }^2({\mu }_2+v{)}^2}dv{\rho }_g(\mathbf{\lambda }),$$

$${\int }_{-1}^1\frac{(1-k)(1+k{\mu }_1)+\frac{F^{\prime }[0]k{\lambda }^{2}v}{1+d{\lambda }^2}({\mu }_2+v)}{(1+k{\mu }_1{)}^2+k^2{\lambda }^2({\mu }_2+v{)}^2}dv=2,$$

$${\int }_{-1}^1\frac{(1-k)k\lambda ({\mu }_2+v)-\frac{F^{\prime }[0]\lambda v}{1+d{\lambda }^2}(1+k{\mu }_1)}{(1+k{\mu }_1{)}^2+k^2{\lambda }^2({\mu }_2+v{)}^2}dv=0.$$

$$\begin{array}{rl}\left(\alpha -\frac{\beta }{\xi }\right)\left[\arctan \left(\xi ({\mu }_2+1)\right)-\arctan \left(\xi ({\mu }_2-1)\right)\right]& \\ +{\mu }_2\beta \log \left(\frac{{\xi }^{-2}+({\mu }_2-1{)}^2}{{\xi }^{-2}+({\mu }_2+1{)}^2}\right)& =2-2\beta ,\end{array}$$

$$\begin{array}{rl}{\mu }_2\beta \left[\arctan \left(\xi ({\mu }_2+1)\right)-\arctan \left(\xi ({\mu }_2-1)\right)\right]& \\ +\frac{1}{2}\left(\alpha -\frac{\beta }{\xi }\right)\log \left(1+\frac{4{\mu }_2}{{\xi }^{-2}+({\mu }_2-1{)}^2}\right)& =0,\end{array}$$

$$\alpha =\frac{1-k}{k\lambda },\beta =\frac{F^{\prime }[0]}{k(1+d{\lambda }^2)},\xi =\frac{k\lambda }{1+k{\mu }_1}.$$

$$\begin{array}{r}\left|\arctan (\xi ({\mu }_2+1))-\arctan (\xi ({\mu }_2-1))\right|=\left|\arctan \left(\frac{2\xi }{1+{\xi }^2({\mu }_2^2-1)}\right)\right|\\ <\left|\arctan \left(\frac{2\xi }{1-{\xi }^2}\right)\right|\to |2\xi +\mathcal{O}({\xi }^2)|,\end{array}$$

$$\left(\alpha \xi -\beta \right)\frac{\arctan (\xi )}{\xi }=1-\beta .$$

$${\mu }_1=\frac{\lambda }{\xi }-\frac{1}{k},$$

$$\varphi (\xi )=\frac{\arctan (\xi )}{\xi },$$

$$\psi (\xi )=\frac{1-\beta }{\alpha \xi -\beta },$$

$$\varphi (k\lambda )<\psi (k\lambda ).$$

$$0<\frac{k}{1-\beta }<1.$$

Full text

We describe a novel mechanism for pattern formation stemming from the stiffness of chemotactic response relying on a kinetic chemotaxis model which includes a recently discovered formalism for the bacterial chemotaxis. We prove instability both for a microscopic description in the space-velocity space and for the macroscopic equation, a flux-limited Keller-Segel equation, which has attracted much attention recently.

A remarkable property is that the unstable frequencies remain bounded, as it is the case in Turing instability. Numerical illustrations based on a powerful Monte Carlo method show that the stationary homogeneous state of population density is destabilized and periodic patterns are generated in realistic ranges of parameters. These theoretical developments are in accordance with several biological observations.

1 Introduction

Collective motion of chemotactic bacteria as E. Coli relies, at the individual level, on a continuous reorientation by runs and tumbles sensing extracellular chemoattractants produced by themselves [1, 2, 3, 4]. It has been established that the length of run is decided by a stiff response to temporal sensing of chemical cues along pathway, i.e., bacteria reduce their tumbling frequency and extend the run length as they sense an increase in concentrations of chemoattractants along the pathway. Thus, the modulation of tumbling frequency in the chemotactic response is an essential mechanism for bacterial communities self-organization.

This paper is concerned with the pattern formation of the population density of run-and-tumble chemotactic bacteria as E. Coli. We describe a novel self-organized pattern formation mechanism stemming from a modulation of tumbling frequency with stiffness in chemotactic response. Our analysis relies on a solid mathematical analysis and simulations using a unique Monte Carlo code.

In order to investigate the multiscale nature in this new self-organized pattern formation mechanism, we rely on a mesoscopic description, i.e., a kinetic reaction-transport equation for the chemotactic bacteria coupled with a reaction-diffusion equation for the chemoattractants. The microscopic dynamic properties such as tumbling rate, modulation in stiff response, and proliferation (division/death) rate are included at the individual level. We consider the following three main ingredients in the pattern formation: (i) the random run-and-tumble motion of bacteria, where the bacteria run linearly with a constant speed when rotating their flagella in counter-clockwise direction, but occasionally change the running directions (tumbling) when rotating their flagella in clockwise direction; (ii) the stiff and bounded signal response to the logarithmic sensing of chemoattractants along the pathway of bacterium, which generates the biased random motion searching for the higher-concentration region of chemoattractant; and (iii) the division/death of bacteria, where the population-growth rate depends on the local population density of bacteria. The kinetic transport equation considered in this paper describes all these ingredients at the microscopic (individual) level.

The pattern formations in the chemotaxis with population growth have been investigated at the macroscopic level by the Keller-Segel type equations [5, 6, 7]; for example, in [9, 10, 11, 12], the pattern formations induced by the properties of chemotaxis, i.e., so-called chemotaxis-induced instability, are demonstrated both theoretically and numerically. Our paper is also concerned with the chemotaxis-induced instability, but which is based on the kinetic transport equation, which up to our knowledge, has not been carried out so far.

The kinetic approach has a distinctive advantage in studying the multiscale mechanism and mathematical hierarchy between the individual dynamics and macroscopic phenomena. It has a long history and was first proposed in [13, 14] and then further developed toward involving the spatiotemporal variation of the chemoattractant along the pathway of bacterium [16], the internal dynamics of the cellular states [15], and the multi-cellular interactions [17, 18]. The mathematical foundations for the kinetic chemotaxis model have been strengthened involving the mathematical hierarchy between kinetic and continuum equations and the existence of solution for the kinetic chemotaxis equation [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33]. The numerical methods to solve the kinetic chemotaxis equations have been also developed in [25, 37, 38, 39, 40, 41, 42]. The use of the kinetic chemotaxis model is also advanced due to the development of experimental technologies, which allow experimentalists to measure the individual velocities and turning angles in the collective motions of bacteria and give access to time scale measurements. For example, in [34, 35, 36], the advantage of the kinetic modeling is demonstrated by the comparison of the numerical and experimental results.

Our analysis also applies to the flux-limited Keller-Segel system, which is a very active research subject nowadays[26, 43, 44, 45, 46]. The flux-limited Keller-Segel system is derived as the asymptotic limit of the kinetic chemotaxis model mentioned earlier in the so-called “diffusion limit” [47]. It incorporates a saturation of the chemotactic flux which avoids the blow-up of solutions. When diffusion is ignored, it has the property of finite speed of propagation. Thus the flux-limited Keller-Segel system can describe collective behaviours observed in various biological systems more realistically. Our instability result is also a new observation for the flux-limited Keller-Segel system.

In this paper, we propose a new mechanism leading to the linear instability of a kinetic chemotaxis equation coupled with a diffusion-reaction equation for the chemoattractant. The kinetic chemotaxis equation involves a population-growth term, which depends on the local population density of bacteria, as well as a chemotactic response function in the integral kernel, which depends on the spatiotemporal variation of the chemoattractant along the pathway of each bacterium. We obtain a linear instability condition based on the stiffness of the response. The stationary homogeneous state of the population density of bacteria becomes linearly unstable and periodic patterns are generated. We also numerically demonstrate the pattern formations by Monte Carlo method in which we vary the parameters involved in the linear instability condition. A theoretical foundation for the Monte Carlo method is also presented.

2 Main result

Since the observation of the run-and-tumble movement of bacteria [48, 49], the kinetic chemotaxis equation has been proposed as an accurate description [14, 15, 16]. In this study, we include a recently advocated formalism for bacterial chemotaxis, i.e., a logarithmic sensing [50] and a stiff and bounded signal response [51]. That is,

[TABLE]

Here $f(t,\bm{x},\bm{v})$ is the microscopic population density of bacteria with a velocity $\bm{v}\in V$ at position $\bm{x}\in\mathbb{R}$ and time $t\geq 0$. The velocity space of $\bm{v}$, $V$ is the surface of the unit ball, i.e., $|\bm{v}|=1$ and $\Omega(\bm{v})$ is the unit measure on $V$. The tumbling kernel $K[D_{t}\log S|_{\bm{v}}]$ represents the stiff and bounded signal response (which is explained in Eq. (5)) to the logarithmic sensing of chemical attractant $S(t,{\bm{x}})$ along the pathway with velocity ${\bm{v}}$. Here $D_{t}X|_{\bm{v}}$ is the material derivative, i.e., $D_{t}X|_{\bm{v}}=\partial_{t}X+\bm{v}\cdot\nabla X$. The concentration of the chemical attractant $S(t,\bm{x})$ and macroscopic population density of bacteria $\rho(t,\bm{x})$ are calculated as, respectively,

[TABLE]

[TABLE]

where $d$ is the molecular diffusion constant. In Eq. (1), $P[\rho]$ is the population-growth rate of bacteria which depends on the local population density $\rho$ as

[TABLE]

Thus, the bacteria may divide when the local population density is lower than unity and the new born bacteria have the same velocities as the parents, but they may die when the local population density is larger than the unity.

We remark that the existence of the traveling wave in the kinetic transport equations with population growth (but without chemotactic responses, i.e., $K[X]=$const.) are proved in [55, 56, 57], where new born particles may choose their velocities according to a prescribed equilibrium velocity distribution. In this paper, we use the simplest population-growth model among those for which the existence of the traveling wave is proved mathematically and the logistic population-growth term is recovered in the continuum limit.

The tumbling kernel $K[X]$ in Eq. (1) is a decreasing function of $X$ and we choose it as

[TABLE]

where $F[X]$ is a smooth and bounded function which satisfies the following properties,

[TABLE]

Here, $F[X]$ represents the chemotactic response of the bacteria, say the response function, and $\chi$ represents the amplitude of modulation in the chemotactic response and takes a constant value between $0\leq\chi<1$.

In Eqs. (1)–(4), all quantities are nondimensionalized by the following characteristic quantities; the characteristic time $t_{0}$ is defined as $t_{0}=|(\frac{d\tilde{P}}{d\rho})_{\rho=1}|^{-1}$, where $\tilde{P}$ represents the population-growth rate in the dimensional form, the characteristic length $L_{0}$ is defined as $L_{0}=t_{0}V_{0}$, where $V_{0}$ is a constant speed of the bacteria. The nondimensional parameter $k$ is defined as $k=1/(t_{0}\psi_{0}$), where $\psi_{0}$ is a mean tumbling frequency of the bacteria. The population density is scaled by that in the uniform stationary state $\rho_{0}$ and the concentration of the chemoattractant is scaled by $(a_{0}/b_{0})\rho_{0}$, where $a_{0}$ is the production rate of chemoattractant by the bacteria and $b_{0}$ is the degradation rate of chemoattractant.

It is easily seen that Eqs. (1)–(4) have a constant uniform solution with $f(t,\bm{x},\bm{v})=1$, $S(t,\bm{x})=1$, and $\rho(t,\bm{x})=1$. In our main result, this uniform solution gives to Turing-like instability [58, 59]. That is, the uniform solution is linearly unstable if the stiffness of the response function $F^{\prime}[0]$ is sufficiently large as

[TABLE]

In addition, the unstable eigenmodes are bounded, i.e., no high frequency oscillations occur and thus patterns are formed.

Furthermore, Eq. (7) includes the linear stability condition of a flux-limited Keller-Segel equation obtained by the asymptotic analysis of the kinetic chemotaxis equation Eq. (1) under a diffusion scaling introduced in Eq. (28). This proves that our instability condition is sharp in the continuum limit ($k\rightarrow 0$).

Figure 1 shows the linear instability diagram. We numerically calculate the minimum values of the right-hand side of Eq. (7) with variation in the Fourier mode $\lambda$ for given values of $k$ and $d$. One can observe that as increasing the stiffness of the response function $F^{\prime}[0]$, the stationary homogeneous state with $f=S=1$ is destabilized (the chemotaxis-induced instability). The critical lines increase monotonically as the diffusion coefficient $d$ increases, so that the stationary homogeneous state is more destabilized when the diffusion coefficient $d$ becomes smaller. It is also seen that as decreasing $k$, the critical line decreases but converges to that for the flux-limited Keller-Segel equation. The stationary homogeneous state is more likely destabilized as decreasing $k$.

3 Linear instability analysis

3.1 Linearization

It is easily seen that Eq. (1) has the uniform solution, $f(t,\bm{x},\bm{v})=S(t,\bm{x})=1$. We carry out the linear instability analysis about this uniform solution. We consider a small perturbation around the uniform solution as

[TABLE]

where $\mu$ is a constant which distinguishes stability ($\mathrm{Re}(\mu)<0$) and instability ($\mathrm{Re}(\mu)>0$). From Eqs. (4)–(6), we can linearize the population-growth rate $P[\rho]$ and tumbling kernel $K[D_{t}\log S|_{\bm{v}}]$ in Eq. (1) as, respectively,

[TABLE]

and

[TABLE]

Thus, Eq. (1) is linearized as

[TABLE]

By taking the Fourier transform of Eqs. (2) and (11), we obtain

[TABLE]

Hereafter, $\mathrm{i}$ represents the imaginary unit, $\bm{\lambda}\in\mathbb{R}$ the wave vector, and $\hat{a}(\bm{\lambda})$ the Fourier transform of the function $a(\bm{x})$ as $\hat{a}(\bm{\lambda})=\int_{\mathbb{R}}a(\bm{x})e^{-\mathrm{i}\bm{\lambda}\cdot\bm{x}}d\bm{x}$. By integrating the above equation as to $\bm{v}$ on $V$, we obtain an equation for $\hat{\rho}_{g}(\bm{\lambda})$ as

[TABLE]

where $\mu_{1}$ is the real part of $\mu$, i.e., $\mu_{1}=\mathrm{Re}(\mu)$, and $\mu_{2}$ is defined as $\mu_{2}=\mathrm{Im}(\mu)/\lambda$. Here we write $\lambda=|\bm{\lambda}|$.

Thus, in order to obtain a non-trivial solution of $\hat{\rho}_{g}(\bm{\lambda})$, the following equations must be simultaneously satisfied, i.e.,

[TABLE]

and

[TABLE]

Further, Eqs. (14) and (15) are analytically calculated as, respectively,

[TABLE]

and

[TABLE]

where

[TABLE]

Note that Eqs. (16) and (17) are symmetric as to the sign of $\mu_{2}$ and $\mu_{2}$=0 always satisfies Eq. (17) irrespective of the values of $\alpha$, $\beta$, and $\xi$. The eigenvalue $\mu\,(=\mu_{1}+\mathrm{i}\mu_{2}\lambda)$ is obtained by solving Eqs. (16) and (17) simultaneously, and the sign of the growing rate $\mu_{1}$ determines the instability of the uniform solution of the kinetic equation, Eq. (1).

Before we consider the instability condition, we first verify that Eq. (16) is never satisfied as $\lambda\rightarrow\infty$. This is explained as follows.

We consider the first term of the left hand side of Eq. (16). In Eq. (18), $\alpha$ and $\beta$ vanishes as $\lambda\rightarrow\infty$ while the limiting values of $\xi$ and $\mu_{2}$ are unknown. In the case that $\xi$ converges to a finite value or diverges as $\lambda\rightarrow\infty$, the first term of the L.H.S of Eq. (16) vanishes because the first factor vanishes while the second factor is bounded. On the other hand, if $\xi$ converges to zero as $\lambda\rightarrow\infty$, the second factor of the first term of the L.H.S of Eq. (16) is estimated as

[TABLE]

so that the first term of the L.H.S of Eq. (16) vanishes as $\lambda\rightarrow\infty$. Thus, the first term of the L.H.S of Eq. (16) vanishes as $\lambda\rightarrow\infty$ regardless the limiting values of $\xi$ and $\mu_{2}$. However, the second term of the L.H.S of Eq. (16) is always non-positive while the R.H.S of Eq. (16) converges to 2 as $\lambda\rightarrow\infty$. Thus, Eq. (16) cannot be satisfied in the limit of large $\lambda$ whatever the limiting values of $\xi$ and $\mu_{2}$ take, so that eigenmodes cannot exhibit large oscillations.

3.2 Instability condition

We now assume the imaginary part of eigenvalue is zero, i.e., $\mu_{2}=0$, and consider a linear instability condition, i.e., $\mu_{1}>0$. Because $\mu_{2}=0$ satisfies Eq. (17), we only consider Eq. (16) with $\mu_{2}$=0, i.e.,

[TABLE]

We can write $\mu_{1}$ as a function of $\xi$, i.e.,

[TABLE]

so that $\mu_{1}$ takes a positive value if and only if $0<\xi<k\lambda$. Thus, the instability condition, i.e., $\mu_{1}>0$, is equivalent to the condition that Eq. (20) has a solution with $0<\xi<k\lambda$.

First, we consider Eq. (20) with $\beta=1$, which is obtained at $\lambda=\sqrt{\frac{1}{d}(F^{\prime}[0]/k-1)}$ when $F^{\prime}[0]/k\geq 1$. It is immediate that Eq. (20) gives $\mu_{1}=-1$, so that $\beta=1$ does not hold the instability condition. (We also note that $\mu_{1}=-1$ is also obtained as $\lambda$=0.)

For $\beta\neq 1$, we consider the condition that the intersection of the following functions,

[TABLE]

and

[TABLE]

exists in $0<\xi<k\lambda$. Note that $\phi(\xi)$ is monotonically decreasing for $\xi\geq 0$ and $\phi(0)=1$ and $\phi(+\infty)=0$, while the behavior of $\psi(\xi)$ depends on the coefficients $\alpha$ and $\beta$.

Case I: $\alpha>0$ and $\beta>1$

The coefficient of $(\xi-\beta/\alpha)^{-1}$ in Eq. (23) is negative, so that $\psi(\xi)$ takes positive values only when $\xi$ is smaller than the asymptotic line, $\xi=\beta/\alpha$. See also Fig. 2(a). Thus, we only consider the intersection of $\phi(\xi)$ and $\psi(\xi)$ in $0<\xi<\beta/\alpha$. It is found that $\psi(\xi)$ monotonically increases in $0<\xi<k\lambda$ because the asymptotic line, $\xi=\beta/\alpha$, is always larger than $\xi=k\lambda$, because $\frac{\beta}{\alpha}>\frac{1}{\alpha}=\frac{k\lambda}{1-k}>k\lambda$. It is also found that $\psi(\xi)$ is always smaller than $\phi(\xi)$ at $\xi=0$, because $\psi(0)=1-\frac{1}{\beta}<1$. Thus, the condition that $\phi(\xi)$ and $\psi(\xi)$ intersect in $0<\xi<k\lambda$ is given by

[TABLE]

Case II: $\alpha>0$ and $0<\beta<1$

The coefficient of $(\xi-\beta/\alpha)^{-1}$ in Eq. (23) is positive, so that $\psi(\xi)$ takes positive values only when $\xi>\beta/\alpha$. It is obvious that if the asymptotic line $\xi=\beta/\alpha$ is larger than or equal to $\xi=k\lambda$, $\phi(\xi)$ and $\psi(\xi)$ do not intersect in $0<\xi<k\lambda$. See also Fig. 2(b). Thus, we only consider the case $\beta/\alpha<k\lambda$, in which $\psi(\xi)$ is positive and monotonically decreases in $\beta/\alpha<\xi<k\lambda$. However, $\psi(\xi)$ is always larger than the unity in this range. This is explained as follows.

From the condition $\beta/\alpha<k\lambda$ and Eq. (18), we obtain

[TABLE]

Since $\psi(\xi)$ monotonically decreases in $\beta/\alpha<\xi<k\lambda$, we get

[TABLE]

Thus, we can conclude that $\phi(\xi)$ and $\psi(\xi)$ do not intersect in $0<\xi<k\lambda$ in this case.

Case III: $\alpha=0$, i.e., $k=1$.

In this case, $\psi(\xi)=1-\frac{1}{\beta}$ is constant, so that the condition that $\phi(\xi)$ and $\psi(\xi)$ intersect in $0<\xi<k\lambda$ is given by Eq. (24).

Case IV: $\alpha<0$ and $\beta>1$.

The coefficient of $(\xi-\beta/\alpha)^{-1}$ in Eq. (23) is positive and the asymptotic line $\xi=\beta/\alpha$ is in the negative region, so that $\psi(\xi)$ takes finite positive values and monotonically decreases in $0<\xi<k\lambda$. See also Fig. 2(b). In addition, $\psi(\xi)$ is smaller than $\phi(\xi)$ at $\xi=0$, because $\psi(0)=1-1/\beta<1$. Thus, the condition that $\phi(\xi)$ and $\psi(\xi)$ intersect in $0<\xi<k\lambda$ is given by Eq. (24).

Case V: $\alpha<0$ and $0<\beta<1$.

The coefficient of $(\xi-\beta/\alpha)^{-1}$ in Eq. (23) is negative and the asymptotic line $\xi=\beta/\alpha$ is in the negative region, so that $\psi(\xi)$ is always negative in $0<\xi<k\lambda$. Thus, $\phi(\xi)$ and $\psi(\xi)$ do not intersect in $0<\xi<k\lambda$.

In summary, the condition that $\phi(\xi)$ and $\psi(\xi)$ intersect in $0<\xi<k\lambda$ is solely given by Eq. (24). Thus, the instability condition of the perturbation with mode $\lambda$ for Eqs. (1)–(4) is written as

[TABLE]

We also remark that the auxiliary condition $\beta>1$, i.e., $0<\lambda<\sqrt{(F^{\prime}[0]/k-1)/d}$, is automatically satisfied under the above condition. Thus, when $F^{\prime}[0]/k$ exceeds the minimum of the right-hand side of Eq. (27) with variation in $0<\lambda<\sqrt{(F^{\prime}[0]/k-1)/d}$, the stationary homogeneous state in the population density becomes linearly unstable. Furthermore, stationary periodic patterns are generated because no unstable eigenmodes exist in the limit of large $\lambda$.

3.3 Continuum limit

We introduce a small parameter $\varepsilon$ and scale $k$, $x$, and $d$ in Eqs. (1) and (2) as follows,

[TABLE]

We also suppose the gradient of the response function is scaled as

[TABLE]

Then, by the asymptotic analysis of small $\varepsilon$ for Eqs. (1) and (2), we obtain the following flux-limited Keller-Segel equation in the limit $\varepsilon\rightarrow 0$,

[TABLE]

[TABLE]

where $U[\log S]$ is defined as

[TABLE]

We remark that the flux $U[\log S]$ is bounded as $|U[\log S]|\leq\frac{\chi}{2}$ because of the bounded signal response Eq. (6c).

When we carry out the linear instability analysis for Eqs. (30)–(32) as in the previous subsection, we obtain the growth rate

[TABLE]

where $\hat{\lambda}$ is the Fourier variable as to $\hat{x}$, i.e., $\hat{\lambda}=\varepsilon\lambda$. Eq. (33) achieves its maximum at $\hat{\lambda}^{2}=(\sqrt{\hat{F}^{\prime}[0]}-1)/\hat{d}$, and the sign of the maximum value of the growth rate $\mu_{1}$ determines the linear stability. Thus, the instability condition for Eqs. (30)–(32) is written as

[TABLE]

Remarkably the above equation is consistent with the linear stability condition of the Keller-Segel system obtained earlier by Nadin, et. al., Ref. [11]. We also note that Eq. (34) can be obtained by using Eq. (7) with the same scaling as Eq. (28). This means Eq. (7) is not only sufficient but also necessary condition for linear instability in the continuum limit.

4 Numerical analysis

The numerical simulations are performed for Eqs. (1)–(3) with the uniform initial density, i.e., $f(0,x,v)=S(0,x)=1$, and periodic boundary condition in the one-dimensional interval $x=[0,L]$, i.e., $f(t,0,v)=f(t,L,v)$ and $S(t,0)=S(t,L)$. The kinetic equation Eq. (1) is solved by the Monte Carlo (MC) method, which has been recently developed in Ref. [42], coupled with the finite volume (FV) scheme for Eq. (2). The details of the MC method and FV scheme is presented in Sec. 5. In the MC simulations, the one-dimensional interval $L$ is set as $L$=100 and divided into 2000 cubic lattice boxes with a side length $\Delta x$=0.05. The 1$\times 10^{7}$ simulation particles are used initially as a total in the whole lattice boxes, and they are distributed randomly into each lattice box. See also Fig. 6. The time-step size $\Delta t$ is set as $\Delta t$=5$\times 10^{-3}$. In the previous study, it has been verified that these simulation parameters produce accurate numerical solutions. Eq. (2) is discretized using the FV scheme on the uniform lattice mesh system with a mesh interval $\Delta x$ and solved implicitly at each time step.

For the population-growth rate $P[\rho]$ and response function $F[X]$, which satisfies Eq. (4) and Eq. (6), respectively, we consider

[TABLE]

and

[TABLE]

Here, $\chi$ and $\delta^{-1}$ represent the amplitude of modulation and stiffness of the response function, respectively. Note that for Eq. (36), $F^{\prime}[0]$ is given as $F^{\prime}[0]=\chi/\delta$. The values of parameters $\chi$, $\delta$, and $d$ are set so as to correspond to the marks A, B, C in Fig. 1 with four values of $k$, i.e., $k$=0.1, 1, 2, and 10. More specifically, we set the parameters as $(\chi/\sqrt{k},\sqrt{k}\delta,d/k)$=(0.5, 0.05, 1) for A, (0.5, 0.0625, 1) for B, (0.5, 0.0625, 0.7) for C, and (0.5, 0.1, 1) for D. We perform the MC simulations for the parameters listed in Table 1 and investigate the compatibility and sharpness of the kinetic instability condition Eq. (7) numerically. Table 1 also shows the prediction of the linear instability by Eq. (7). For example, the parameter set C with $k=1.0$ is slightly above the critical line in Fig. 1, so that the kinetic instability condition Eq. (7) affirms the occurrence of periodic pattern formation. However, for the parameter set C with $k=2.0$, which is slightly below (but very close to) the critical line in Fig. 1, we cannot confirm neither the linear instability nor homogeneous state theoretically because Eq. (7) is a sufficient condition of the linear instability. However, the MC simulation can numerically demonstrate how sharply the kinetic instability condition Eq. (7) can predict the pattern formations.

Figure 3 shows the time progress of the population density $\rho$. It is obviously seen that the results for the black squares in Table 1 exhibit the stationary periodic patterns after some transient period. On the other hand, those for the white squares in Table 1 do not show any distinct patterns. In order to quantify the patterns, we also calculate the power spectra of the Fourier transforms of population density profiles between $t$=[400,500]. Figure 4 shows the results of the power spectra. Here, we take the averages of the snap shots of power spectra obtained at $t$=[400,500] with a time interval 4. The power spectrum for the continuum equations Eqs. (30)–(32) is calculated from the snap shot of the population density at $t=200$. For the parameter sets shown as the black squares $\blacksquare$ in Table 1 (See Fig. 4(a)), steep peaks are observed and the first peaks appear in $0<\hat{\lambda}<\sqrt{(\hat{F}^{\prime}[0]-1)/\hat{d}}$ for each parameter set. The second and third peaks also appear as the non-linear effects although they are much smaller than the first peaks. The power spectra decreases as the wave number $\hat{\lambda}$ decreases from the first-peak position, while they neither grow nor damp at the large wave number, so that a plateau regime appears at the large wave numbers. The peaks of the power spectrum for the continuum equations with the parameter set B coincide with those for the parameter set B with $k=0.1$. However, no plateau regime appears for the continuum equations. This result confirms that there are no eigenmodes of the linearized kinetic equation.

For the parameter sets shown as the white squares $\square$ in Table 1 (See Fig. 4(b)), the behaviors of power spectra are similar to those for the black squares in Table 1, i.e., Fig. 4(a), except the peak behaviors. In Fig. 4(b), we cannot see the steep first peaks as seen in Fig. 4(a). The small peaks of A, B, and C in Fig. 4(b) are even smaller than the second peaks appearing in the figure (a). For the parameter set C, we cannot see any peaks. Thus, we cannot observe the stationary periodic patterns evidently in Fig. 3 for the parameter sets shown as white squares in Table 1.

The parameter set A and B are very close to the critical lines for $k$=2.0 and $k$=0.1 in Fig. 1, respectively. (The parameter set A is only slightly lower and B is only slightly upper than each critical line.) However, the result of the parameter set B shows the stationary periodic pattern evidently but the result of the parameter set A does not show it. These numerical results demonstrate that the critical lines in the instability diagram Fig. 1 can predict sharply the occurrence of the linear instability and pattern formations.

Finally, we show the effect of modulation amplitude $\chi$ on the instability pattern profile. Figure 5 shows the pattern profile for the parameter set B with $k=0.1$ in Table 1 and that obtained when the modulation amplitude $\chi$ is about four times larger than the parameter set B for $k=0.1$. The pattern profile obtained for the former parameter is periodic oscillation around the initial uniform state $\rho=1$. The formation of periodic oscillatory patter applies to all of the instability patterns obtained in the parameter sets in Table 1, where $\chi$ is fixed as $\chi/\sqrt{k}=0.5$. However, in the case of a large modulation amplitude, i.e., $\chi/\sqrt{k}=2.06$, the population of bacteria becomes localized due to a strong chemotactic response so that it forms periodical bounded spikes. The boundedness in instability pattern formation stems from the flux-limited property in the non-linear stiff response function. Incidentally, the boundedness property is also inherent in the flux-limited Keller-Segel equation which is obtained by the asymptotic analysis of the kinetic chemotaxis equation. The variety of solution types is observed in the flux-limited Keller-Segel equation, which will be addressed in our forthcoming paper.

5 Monte Carlo method

The motions of the chemotactic bacteria are simulated by using the Monte Carlo particles which follows the process described by the kinetic chemotaxis equation Eq. (1) coupled with the reaction-diffusion equation for the chemoattractant Eq. (2). The one-dimensional space interval $x\in[0,L]$ is divided into $I$ cubic lattices with a uniform side length $\Delta x$, i.e., $L=I\Delta x$. See also Fig. 6. The reaction-diffusion equation Eq. (2) is implicitly solved on the uniform lattice system by using a finite volume scheme as

[TABLE]

where $S_{i}^{n}$ and $\rho_{i}^{n}$ are the concentration of chemical attractant and population density of bacteria in the $i$th lattice site $[i\Delta x,(i+1)\Delta x]\times\Delta x^{2}$ at a time $t=n\Delta t$, respectively. Here we also set $S^{n}_{-1}=S^{n}_{I-1}$ and $S^{n}_{I}=S^{n}_{0}$ according to the periodic condition. Hereafter the superscript represents the time-step number, the subscript without parenthesis represents the lattice-site number, and the subscript in parenthesis represents the index of each MC particle. The population density of bacteria $\rho^{n}_{i}$ is calculated from the number of the MC particles involved in the $i$th lattice site, $M_{i}^{n}$, as

[TABLE]

where $M$ is the number of MC particles involved in one lattice site in the reference state, i.e, $M=N/I$, where $N$ is the total number of MC particles in the initial state.

The MC simulation is conducted using the following steps. Hereafter, the position and velocity of the $l$th particle are expressed as ${\bm{r}}_{(l)}^{n}$ and ${\bm{v}}_{(l)}^{n}$, respectively.

At $\hat{t}=0$, MC particles are distributed according to the initial density. In each lattice site, MC particles are distributed uniformly at random positions and their velocities ${\bm{v}}$ are determined by the probability density $f_{i}^{0}({\bm{v}})/(4\pi\rho_{i}^{0})$. 2. 1.

Particles move with their velocities for a duration $\Delta\hat{t}$:

[TABLE]

where $N^{n}$ is the total number of simulation particles at a time step $n$, i.e., $N^{n}=\sum_{i=0}^{I-1}M_{i}^{n}$. The particles that move beyond the boundaries are inserted at the opposite boundaries according to the periodic boundary conditions. 3. 2.

At each lattice site, the macroscopic population densities $\rho^{n+1}_{i}$ and concentrations of chemical cues $S^{n+1}_{i}$ ($i=0,\cdots,I-1$) are calculated by Eqs. (38) and (37), respectively. 4. 3.

The tumbling of each particle is calculated using the scattering kernel in Eq. (1). The tumbling of the $l$th particle may occur with a probability

[TABLE]

where $D_{t}\log S_{(l)}$ represents the temporal variation of the chemical cue experienced by the $l$th MC particle along the pathway, and is defined by the following forward difference,

[TABLE]

The local concentration of chemical cue at the position of the $l$th MC particle is calculated by using linear interpolation between the neighboring lattice sites, i.e.,

[TABLE]

This generates a chemoattractant gradient and MC particles that stay at the same lattice site after a single time step can sense the chemoattractant gradient along their pathways.

For the particle that is judged to tumble, say the $l_{t}$th particle, a new velocity after the tumbling, ${\bm{v}}^{n+1}_{(l_{t})}$, is determined randomly as,

[TABLE]

Here $U_{1}$ and $U_{2}$ are the uniform random variables between 0 and 1. 5. 4.

The divisions/deaths are judged for all MC particles. The division (or death) occurs with a probability $|P[\rho^{n}_{i}]|\Delta t$, if $|P[\rho_{i}^{n}]|$ is positive (or negative), where $\rho_{i}^{n}$ is the local population density at the lattice site where each MC particle is involved. For a particle that is judged to undergo division, e.g., the $l$th particle, a new particle with the same velocity ${\bm{v}}_{(l)}$ is created at a random position within the same lattice site. The numbers of MC particles involved in each lattice site are counted, $M_{i}^{n+1}$ ($i=0,\cdots,I-1$), and the total number of simulation particles is updated as $N^{n+1}$. 6. 5.

Return to step 1 with the obtained ${{\bm{r}}}_{(l)}$, ${{\bm{e}}}_{(l)}$, $S_{(l)}$ ($l$=1,$\cdots$,$M$) at the new time step.

Weak formulation The overall procedure corresponds to the first-order time difference equation of the kinetic chemotaxis equation in the weak formulation. We consider the following moment equation,

[TABLE]

where $\Phi(\bm{x},\bm{v})$ is an arbitrary smooth function which vanishes outside the computational domain on $\bm{x}$.

Here $<\quad,\quad>$ defines the integration of the arbitrary functions $a(\bm{x},\bm{v})$ and $b(\bm{x},\bm{v})$ as

[TABLE]

We consider the functions $f^{A}$, $f^{B}$, and $f^{C}$ which are determined, respectively, as

[TABLE]

Then, $f(t+\Delta t,\bm{x},\bm{v})$ is obtained by the sum of three functions:

[TABLE]

In the MC method, the microscopic population density $f$ at the $n$th time step is approximated as

[TABLE]

By substituting Eq. (48) into Eq. (46a) we obtain the following equation:

[TABLE]

It is easily seen that the function $f^{A}$ corresponds to the distribution that is obtained by the moving process, say ${\cal S}_{1}$ in the MC method, i.e.,

[TABLE]

where ${\cal S}_{1}\{\quad\}$ represents the operator of the process 1 in the MC method.

The tumbling process and division/death process in the MC method, say ${\cal S}_{3}$ and ${\cal S}_{4}$ respectively, are performed independently in each lattice site. Thus, we consider the processes in a fixed lattice site, say the $i$th lattice site, with using a test function written as $\Phi(\bm{x},\bm{v})=\phi_{i}(\bm{x})\psi(\bm{v})$, where $\phi_{i}(\bm{x})$ is an arbitrary smooth function which vanishes outside the $i$th lattice site. Then, by substituting Eq. (48) into Eqs. (46b) and (46c) we obtain the following equations:

[TABLE]

[TABLE]

where the subscript $m$ ($m=1,\cdots,M_{i}^{n}$) counts the MC particles involved in the $i$th lattice site at the $n$th time step. Here $D_{t}\log S_{(m)}^{n}$ is defined in Eq. (41), $\bar{\psi}=\frac{1}{4\pi}\int\psi(\bm{v})d\Omega(\bm{v})$, and $\bar{\phi}_{i}=\phi(\bm{x}_{i+\frac{1}{2}})$, where $\bm{x}_{i+\frac{1}{2}}$ represents the center of the $i$th lattice site. In the derivation of Eq. (51), we approximate the tumbling kernel as

[TABLE]

We note that the second equality in Eq. (5) is obtained under the assumption of the uniform distribution of large number of particles in each lattice site.

In the tumbling process ${\cal S}_{3}$, the particle which creates the tumbling, say the $m_{t}$th particle, changes its velocity randomly as $\bm{v}^{n}_{m_{t}}\rightarrow\bm{v}^{n+1}_{m_{t}}=\bm{u}_{t}$, where the random velocity $\bm{u}_{t}$ is given by Eq. (43). We now suppose that $M_{t}$ particles make the tumbling in the $i$th lattice site, then the microscopic population density $f^{n}$ in the $i$th lattice site changes as

[TABLE]

and the moment is written as

[TABLE]

We introduce the stochastic variables $Z_{m}$ and $Z$ which are defined as

[TABLE]

and

[TABLE]

If the tumbling occurs for $M_{t}$ particles in the $i$th lattice site, the realized value of $Z$ is written as

[TABLE]

Here we use Eq. (55). On the other hand, the expected value of $Z_{m}$ is written as

[TABLE]

and the expected value of $Z$ is written by the moment Eq. (51) as

[TABLE]

From Eqs. (58) and (60), it is seen that the microscopic population density obtained by the tumbling process on $f^{n}$, ${\cal S}_{3}\{f^{n}\}$ is approximated by the sum of $f^{n}$ and $f^{B}$ under the assumption of the law of large numbers, i.e., $\hat{Z}\simeq E_{p}(Z)$,

[TABLE]

In the division/death process ${\cal S}_{4}$, the particles in the $i$th lattice site may divide (or die) if the local macroscopic population density $\rho_{i}$ is smaller (or larger) than unity. In the divisions, where $P[\rho_{i}]>0$, the particle creates a new particle with the same velocity at a random position within the same lattice site. Thus, when the divisions occur for $M_{c}$ particles for the distribution Eq. (48), the microscopic population density in the $i$th lattice site changes as

[TABLE]

where $\bm{w}$ is a random vector whose components are uniform random numbers in [$-\frac{1}{2}$,$\frac{1}{2}$]. The moment is written as

[TABLE]

In deaths, where $P[\rho_{i}]<0$, the particle is just removed from the lattice site. Thus, when the deaths occurs for $M_{c}$ particles in the $i$th lattice site, the microscopic population density obtained after the death process and its moment are written as changing the sign of the second term and replacing $\bm{x}_{i+\frac{1}{2}}+\Delta x\bm{w_{c}}$ with $\bm{r}_{(m_{c})}^{n}$ in Eqs. (62) and (63), respectively. Hence, under the assumption of the uniform distribution of large number of particles in each lattice site, the moment for the population density after the division/death process ${\cal S}_{4}$ can be written as

[TABLE]

We remark that $M_{c}/M$ is estimated as $M_{c}/M\sim{\cal O}(\Delta t)$.

We introduce the stochastic variable $\Xi_{m}$ and $\Xi$ which are defined as

[TABLE]

and

[TABLE]

If the divisions(or deaths) occur for $M_{c}$ particles, the realized value of $\Xi$ is written as

[TABLE]

On the other hand, the expected value of $\Xi$ is written by Eq. (5) as

[TABLE]

Thus, the microscopic population density obtained by the division/death process on $f^{n}$, ${\cal S}_{4}\{f^{n}\}$ is approximated by the sum of $f^{n}$ and $f^{C}$ under the assumption of the law of large numbers, i.e., $\hat{\Xi}\simeq E_{p}(\Xi)$,

[TABLE]

We remark that the error of the above equation is estimated at most ${\cal O}(\Delta t^{2})$ when $\Delta t$ is at most the second order of $\Delta x$, i.e., $\Delta t\propto\Delta x^{\alpha}$ s.t. $\alpha\leq 2$.

In the MC simulation, the processes ${\cal S}_{1}$, ${\cal S}_{3}$, and ${\cal S}_{4}$ are successively conducted. For example, the process ${\cal S}_{2}$ is performed on the distribution obtained by the process ${\cal S}_{1}$ on $f^{n}$;

[TABLE]

where $f^{B^{\prime}}$ is obtained by replacing $f^{n}$ with $f^{A}$ in Eq. (46b). However, the difference of $f^{A}$ and $f^{n}$ is ${\cal O}(\Delta t)$, so that $f^{B^{\prime}}$ can be replaced with $f^{B}$ within the difference of ${\cal O}(\Delta t^{2})$. Similarly, it is seen that the microscopic population density obtained by the successive three processes of ${\cal S}_{1}$, ${\cal S}_{3}$, and ${\cal S}_{4}$ approximates the microscopic population density at the next time step which satisfies the weak formulation Eq. (5) within the difference of ${\cal O}(\Delta t^{2})$, i.e.,

[TABLE]

6 Concluding remarks

We studied the self-organized pattern formation of chemotactic bacteria based on a kinetic chemotaxis model which includes a recently advocated formalism for bacterial chemotaxis, i.e., the logarithmic sensing of chemical cues along the pathway of bacterium and stiff and bounded signal response. We have discovered a novel linear instability condition Eq. (7) stemming from the stiffness of chemotactic response. Apart from the macroscopic description, we have been able to uncover the instability mechanism at the microscopic level. The stationary homogeneous state of the macroscopic population density becomes linearly unstable and stationary periodic patterns are generated under the linear instability condition. A remarkable property is that no eigenmodes exist in the large-oscillation limit in the linearized kinetic equation, which explains that pattern formations occur as observed in experiments. Our new dispersion relation for instability also turns out to be sharp in the macroscopic limit, i.e., the flux-limited Keller-Segel equation.

MC simulations rigorously based on the kinetic chemotaxis model are performed with changing the parameters involved in the linear instability condition. The numerical results demonstrate that the obtained linear instability condition is compatible and even sharply predicts the occurrence of the periodic pattern formations. See Fig. 3. The power spectra of the macroscopic population density show the plateau regime at the large wave numbers, where the perturbations neither grow nor damp irrespective of the linear instability condition. See Fig. 4 This observation is compatible with the fact that no eigenmodes exist in the large oscillation limit in the linearized kinetic chemotaxis equation. Unexpectedly, the instability pattern undergoes transitions from the periodic oscillation around the uniform state $\rho=1$ to the periodic localized spikes over the zero-density state $\rho=0$ as increasing the modulation amplitude in chemotactic response. See Fig. 5.

The obtained instability condition Eq. (7) includes three control parameters, i.e., $k$, $d/k$, and $F^{\prime}[0]/k$, which are written in the dimensional form as

[TABLE]

Here $\psi_{0}$ is the mean tumbling frequency, $V_{0}$ is the running speed of bacteria, $b_{0}$ is the degradation rate of chemoattractant, and $t_{0}$ is the characteristic time which corresponds to the doubling time in cell division for Eq. (35). See also the paragraph above Eq. (7). Here $\tilde{d}$ and $\tilde{\delta}^{-1}$ are the diffusion coefficient of chemoattractant and stiffness of response function, respectively, in the dimensional form. The values of control parameters are estimated from experimental data in Ref. [34] (See Table 2) as $k=7.5\times 10^{-5}$, $d/k=3.8\times(t_{0}b_{0})^{-1}$, and $F^{\prime}[0]/k=12$. It is difficult to measure the degradation rate of chemoattractant $b_{0}$ in experiments. In some references [36, 62, 63], the value of $b_{0}$ is estimated as $b_{0}=4\times 10^{-3}\sim 5\times 10^{-2}$ from the comparison between experimental and numerical results, so that we may estimate $d/k$ as $d/k\lesssim 1$. Thus, from our analysis, we can expect the stationary homogeneous state becomes destabilized and pattern formation occurs for chemotactic bacteria. Furthermore, for example, in Ref. [2], it is argued that the pattern formation is suppressed by reduction of chemotactic sensitivity. This argument is also consistent with our instability condition. Although it remains to be assessed how quantitatively our instability condition explains the experimental results in terms of pattern formation, the present study convinces us that the self-organized pattern formation occurs due to the modulation of stiff response in chemotaxis in a realistic range of parameters.

Finally, our powerful MC method derived rigorously here has a possible advantage that can be extended to include internal states stemming from an intra-cellular chemical pathway. The kinetic chemotaxis model used in this study is based on a simplified model where the intra-cellular adaptation dynamics in chemotactic response is ignored and replaced by a instantaneous material derivative of chemical cue along the pathway of bacterium. The time scale of adaptation, say $\tau_{M}$, is larger than the inverse of tumbling rate $\psi_{0}^{-1}$ but is much smaller than the doubling time $t_{0}$, i.e., $\psi_{0}^{-1}/t_{0}<\tau_{M}/t_{0}\ll 1$ [60, 61]. Thus, the adaptation dynamics may be significant for the pattern formation. In order to consider the adaptation dynamics internal states have to be taken into account in the tumbling kernel[15, 29, 31]. In a forthcoming paper, we plan to extend the MC method toward this direction and thus be able to challenge problems with another time scale related to the internal-state dynamics.

Acknowledgements

This study was financially supported by JSPS KAKENHI Grant Number 15KT0110 and 16K17554 and Institut Henri Poincaré RIP program.