∎ 11institutetext: Graduate School of Information Science, University of
Hyogo, Kobe 650-0047, Japan.
11email<EMAIL_ADDRESS>
#  Numerical study of the volcano effect in chemotactic aggregation based on a
kinetic transport equation with non-instantaneous tumbling

Shugo Yasuda

###### Abstract

Aggregation of chemotactic bacteria under a unimodal distribution of chemical
cues was investigated by Monte Carlo (MC) simulation based on a kinetic
transport equation, which considers an internal adaptation dynamics as well as
a finite tumbling duration.

It was found that there exist two different regimes of the adaptation time,
between which the effect of the adaptation time on the aggregation behavior is
reversed; that is, when the adaptation time is as small as the running
duration, the aggregation becomes increasingly steeper as the adaptation time
increases, while, when the adaptation time is as large as the diffusion time
of the population density, the aggregation becomes more diffusive as the
adaptation time increases. Moreover, the aggregation profile becomes bimodal
(volcano) at the large adaptation-time regime when the tumbling duration is
sufficiently large while it is always unimodal at the small adaptation-time
regime.

A remarkable result of this study is the identification of the parameter
regime and scaling for the volcano effect. That is, by comparing the results
of MC simulations to the continuum-limit models obtained at each of the small
and large adaptation-time scalings, it is clarified that the volcano effect
arises due to the coupling of diffusion, adaptation, and finite tumbling
duration, which occurs at the large adaptation-time scaling.

###### Keywords:

chemotaxis, aggregation, volcano effect, kinetic transport equation, non-
instantaneous interaction, Monte Carlo simulation

##  1 Introduction

Chemotactic bacteria, such as Escherichia coli, migrate by alternate running
and tumbling, where the bacteria run straightly for a certain duration
(typically approximately 1 second) and subsequently tumble very quickly
(typically, the tumbling duration is approximately 0.1 seconds) BB1972 . Since
the length of the run is modulated according to the temporal variation of
extracellular chemical cues sensed by the bacteria along their moving pathway,
collective dynamics such as traveling waves and aggregations occurs at the
population level when the spatial gradients of external chemical cues are
sufficiently large.

Kinetic transport equations have been proposed to model the run-and-tumble
motion of chemotactic bacteria ODA1998 ; HO2000 ; EO2004 ; CMPS2004 ; DS2005
and successfully utilized to elucidate the mathematics and physics behind the
complicated collective dynamics SCBPBS2011 ; EGBAV2016 ; XXT2018 ; C2020 .
Since the duration of tumbling is much shorter than the running duration, the
tumbling duration is usually ignored, and instead, the instantaneous velocity
jump process is considered in the kinetic transport equations. However, it is
fundamental and important to know how that small but finite duration of
tumbling affects the collective dynamics.

Recently, studies on non-instantaneous interactions in biological phenomena
have become increasingly prevalent. For example,  in Ref. XXT2018 , concentric
stripe patterns created by engineering E. coli, which have longer tumbling
duration than the wild-type E. coli, was investigated based on a kinetic
transport model with non-instantaneous tumbling process.  Additionally, in
Ref. LTLTH2022 , the sub-diffusion of hydration water molecules that stay
attached to proteins between jumps was investigated. Moreover, in Ref. KST2022
, a novel kinetic transport model with non-instantaneous collision operator
was proposed to consider the contact inhibition of movement upon cell-cell
collisions or collisions under cell-cell adhesion. The present study also aims
to contribute to a fundamental understanding of the effects of the non-
instantaneous interaction (i.e., the non-instantaneous tumbling in this study)
on a simple chemotactic aggregation of bacteria.

In this paper, we regard the kinetic transport equation, which considers the
finite tumbling duration as well as the internal adaptation dynamics EO2004 ,
and clarify their effects on chemotactic aggregation via numerical
simulations.  Especially, we are concerned with the volcano effect (i.e.,
bimodal aggregation of chemotactic bacteria, which was first observed in an
micro-scale experiment MBBO2003 ). The volcano effect was previously
investigated numerically based on a similar kinetic transport model SM2011
and an individual-based simulation JJP2018 . However, the biological and
physical conditions under which the volcano effect takes place and the
multiscale mechanisms behind the volcano effect are still largely unknown. We
tackle with these problems by means of a Monte Carlo (MC) simulation of the
kinetic transport model. The results of MC simulations are also compared with
the continuum-limit equations obtained at different scalings of the adaptation
time, which were previously derived in literatures, e.g., Refs. CMPS2004 ;
EO2004 ; HO2000 ; PTV2016 ; PSTY2020 ; XXT2018 . An important result of this
study is the identification of the parameter regime and scaling for the
volcano effect to arise. That is, it will be clarified that the volcano effect
arises due to the coupling of diffusion, adaptation, and finite tumbling
duration, which occurs at a certain large adaptation-time scaling. These
results contribute to advance the elucidation of the mathematics behind the
volcano effect.

Incidentally, the effect of the adaptation time on the chemotactic aggregation
was investigated in our previous study Y2021  based on the kinetic transport
equation without the finite tumbling duration. In the previous study, the
volcano effect was never observed, but intead, the trapezoidal aggregation,
where the aggregation profile in the central region is rather flat, was
observed at the large adaptation-time regime. The present study is an
extension of the previous study to include the finite tumbling duration.
Interestingly, as it will be seen in Sec. 3, although the modification of the
continuum-limit equation due to the finite tumbling duration is very small, it
is crucial to produce the volcano effect.

The rest of the paper is organized as follows: In Sec. 2, we present the basic
kinetic transport model, which considers both the internal adaptation dynamics
and the tumbling process. The nondimensionalization of the kinetic transport
model is also given. In Sec. 3,  asymptotic equations of the kinetic transport
model, which are utilized to compare with the MC results obtained in Sec. 4,
are summarized.  In Sec. 4, numerical simulations of chemotactic aggregation
under an exponential distribution of chemical cues are implemented for one-
and two-dimensional spaces by using the MC code of the kinetic transport
model. From the numerical results, the emergence of bimodal aggregation and
its microscopic and macroscopic mechanism are discussed in detail. Finally, we
provide concluding remarks in Sec. 5.

##  2 Forumulation

We consider the run-and-tumble bacteria, where the bacteria modulate their run
length according to the memory of the external chemical cue sensed along their
moving pathway, while in the tumbling phase, they alter their moving direction
in a uniformly random way in a short period of time.

The memory of the external chemical cue is described by the internal state
determined via intracellular chemical signal transduction, which has been
extensively studied by many authors (for instance, see Refs BL1997 ; KLBTS2005
; TSB2008 ; JOT2010 ; H2012 ). The mathematical model of intracellular signal
transduction is very complicated in general, but the fundamental property
necessary for the chemotactic response is described by the excitation and the
adaptation process SPO1997 . In the present paper, we consider the following
simple adaptation dynamics of the internal state $m\in R$:

$\dot{m}=\frac{M(S)-m}{\tau_{a}},$

where $\tau_{a}>0$ is the adaptation time and $M(S)$ denotes the local
equilibrium of the internal state according to the local concentration of the
external chemical cue $S$.

The bacteria modulate their tumbling frequency according to the deviation of
the internal state $m$ from the local equilibrium state $M(S)$. More
specifically, the bacteria temporarily decrease (or increase) the tumbling
frequency when the local equilibrium state $M(S)$ is higher (or lower) than
the current internal state $m$. The modulated tumbling frequency returns to
the basal state when the internal state approaches the local equilibrium state
via the adaptation dynamics. Thus, the adaptation process allows the bacteria
to further respond to a subsequent change in the external chemical cue.

Since the modulation of the tumbling frequency is stiff and bounded BSB1983 ,
we write it as follows:

$\Lambda_{\delta}(M(S)-m)=1-F\left(\frac{M(S)-m}{\delta}\right),$ (1)

where $\delta$ represents the stiffness of the chemotactic response and the
response function $F(X)$ has the following property:

$F(0)=0,F^{\prime}(X)>0,F(x\rightarrow\pm\infty)\rightarrow\pm\chi,$

where $\chi$ ($0<\chi<1$) represents the modulation amplitude.

When we write the population density of running cells with a velocity
$v\in[v,v+dv]$ and an internal state $m\in[m,m+dm]$ at time $t>0$ and at space
$x\in R^{d}$ (where $d$ is the dimension of space) as
$d\rho_{f}=f(t,x,v,m)dvdm$ and that of the tumbling cells with an internal
state $m\in[m,m+dm]$ as $d\rho_{g}=g(t,x,m)dm$, the time evolution of the
densities $f(t,x,v,m)$ and $g(t,x,m)$ is described as follows:

$\partial_{t}f+v\cdot\nabla_{x}f+\partial_{m}\left\\{\left(\frac{M(S)-m}{\tau_{a}}\right)f\right\\}=\mu\frac{g}{||V||}-\lambda\Lambda_{\delta}(M(S)-m)f$
(2a)
$\partial_{t}g+\partial_{m}\left\\{\left(\frac{M(S)-m}{\tau_{a}}\right)g\right\\}=\lambda\Lambda_{\delta}(M(S)-m)\int_{V}f(t,x,v,m)dv-\mu
g,$ (2b)

where $\lambda>0$ is the mean tumbling frequency at the basal state, $\mu>0$
is the mean frequency when the tumbling cells change to the running cells, and
$||V||$ is the volume of the velocity space, i.e., $||V||=\int_{V}dv$. Here,
the velocity space is the surface of the ball (i.e., $V=|v|\mathbb{S}^{d}$).
We also introduce the notation $\nu=\mu^{-1}$, which denotes the mean tumbling
duration.

We also remark that since the duration of tumbling is governed by a Poisson
process with a constant rate $\mu$ in Eq. (2), very short tumbling events also
occur. Then it is quite an idealization that the post-tumbling direction is
uniformly randomized. The extension to the non-uniform and non-instantaneous
tumbling kernel would be an important future study.

The total population density $\rho$, the population density of running cells
$\rho_{f}$, and the population density of tumbling cells $\rho_{g}$ are given
as follows:

$\rho(t,x)=\rho_{f}(t,x)+\rho_{g}(t,x),$ (3a)
$\rho_{f}(t,x)=\int_{\mathbb{R}}\int_{V}f(t,x,v,m)dvdm,$ (3b)
$\rho_{g}(t,x)=\int_{\mathbb{R}}g(t,x,m)dm.$ (3c)

###  2.1 Nondimensionalization

We introduce the nondimensional quantities as follows:

$\displaystyle\hat{f}=f/(\rho_{\mathrm{c}}/||V||),\quad\hat{g}=g/\rho_{\mathrm{c}},\quad\hat{t}=t/t_{\mathrm{c}},\quad\hat{x}=x/L_{\mathrm{c}},\quad\hat{v}=v/\mathrm{v_{c}}.$
(4)

Here, the subscript “c” denotes the characteristic quantities, and we set
$\mathrm{v_{c}}=|v|$ in the following. Then, we can rewrite (2) as follows:

$\sigma\partial_{\hat{t}}\hat{f}+\hat{v}\cdot\nabla_{\hat{x}}\hat{f}+\partial_{m}\left\\{\left(\frac{M(S)-m}{\hat{\tau}}\right)\hat{f}\right\\}=\frac{1}{\varepsilon}\left[\hat{\mu}\hat{g}-\Lambda_{\delta}(M(S)-m)\hat{f}\right],$
(5a)
$\sigma\partial_{\hat{t}}\hat{g}+\partial_{m}\left\\{\left(\frac{M(S)-m}{\hat{\tau}}\right)\hat{g}\right\\}=\frac{1}{\varepsilon}\left[\Lambda_{\delta}(M(S)-m)<\hat{f}>-\hat{\mu}\hat{g}\right],$
(5b)

where $<\hat{f}>$ is the average of $\hat{f}$ over the velocity space
$\hat{V}$, which is defined as follows:

$<f>=\frac{1}{||\hat{V}||}\int_{\hat{V}}f(\hat{t},\hat{x},\hat{v},m)d\hat{v},$
(6)

with $||\hat{V}||=\int_{\hat{V}}d\hat{v}$. Here, we also introduce the
following nondimensional parameters:

$\displaystyle\sigma=L_{\mathrm{c}}/(\mathrm{v}_{\mathrm{c}}t_{\mathrm{c}}),\quad\varepsilon=\mathrm{v}_{\mathrm{c}}/(\lambda
L_{\mathrm{c}}),\quad\hat{\tau}=\tau_{a}/(L_{\mathrm{c}}/\mathrm{v}_{\mathrm{c}}),\quad\hat{\mu}=\mu/\lambda.$
(7)

Here, $\sigma$ is the time parameter and $\varepsilon$ represents the mean run
length at the reference state.

The population densities defined in Eq. (3) are written as follows:

$\hat{\rho}(\hat{t},\hat{x})=\hat{\rho}_{f}(\hat{t},\hat{x})+\hat{\rho}_{g}(\hat{t},\hat{x}),$
(8a)
$\hat{\rho}_{f}(\hat{t},\hat{x})=\int_{\mathbb{R}}<\hat{f}>(\hat{t},\hat{x},m)dm,$
(8b)
$\hat{\rho}_{g}(\hat{t},\hat{x})=\int_{\mathbb{R}}\hat{g}(\hat{t},\hat{x},m)dm.$
(8c)

We note that when the tumbling duration is negligibly small compared to the
running duration (i.e., $\hat{\nu}=\hat{\mu}^{-1}\ll 1$), the density of
tumbling cells $\hat{g}$ becomes negligibly small $\hat{g}\ll 1$, and thus,
Eq. (5) is reduced as follows:

$\sigma\partial_{t}\hat{f}+v\cdot\nabla_{x}\hat{f}+\partial_{m}\left\\{\left(\frac{M(S)-m}{\tau}\right)\hat{f}\right\\}=\frac{1}{\varepsilon}\Lambda_{\delta}(M(S)-m)\left(<\hat{f}>-\hat{f}\right).$
(9)

This equation was used in a previous study Y2021  where instantaneous tumbling
events were considered. Thus, this study is an extension of the previous study
to consider the non-instantaneous tumbling events in the chemotactic
aggregations.

##  3 Continuum-limit equations

It has been proved that different types of continuum-limit (i.e.,
$\varepsilon\rightarrow 0$) equations are obtained by the asymptotic analysis
of the kinetic transport equation at different scalings of the adaptation time
EO2004 ; PTV2016 ; XXT2018 ; PSTY2020 . In this study, we utilize the
continuum-limit equations to confirm the asymptotic behaviors of the MC
simulations of the kinetic transport equation (5), which will be given in the
next section. In this section, we summarize the continuum-limit equations
obtained at two different scalings of the adaptation time. An asymptotic
relation between the two different continuum-limit equations is also briefly
explained.

We consider small and large adaptation-time scalings, i.e.,
$\hat{\tau}=O(\varepsilon)$ and $\hat{\tau}=O(1/\varepsilon)$, respectively.
We also consider the diffusive time scale $\sigma=\varepsilon$.  These
settings of the time scale parameters are physically interpreted as follows:
The time scale parameter $\sigma=\varepsilon$ reads that the characteristic
time $t_{\mathrm{c}}$ corresponds to the diffusion time of the population
density, i.e., $t_{\mathrm{c}}=t_{d}$, where the diffusion time $t_{d}$ is
defined as

$t_{d}=L_{\mathrm{c}}^{2}/D_{\rho}$ (10)

with the diffusion constant defined as $D_{\rho}=\mathrm{v_{c}}^{2}/\lambda$.
The small adaptation time scaling $\hat{\tau}=O(\varepsilon)$ indicates that
the adaptation time is comparable to the running duration (i.e.,
$\tau_{a}\sim\lambda^{-1}$), while the large adaptation-time scaling
$\hat{\tau}=O(1/\varepsilon)$ indicates that the adaptation time is comparable
to the diffusion time (i.e., $\tau_{a}\sim t_{d}$).

In the following, we only consider the case where the stiffness of the
modulation function is moderate, i.e., $\delta=O(1)$ in Eq. (1) (although the
stiff chemotactic response, such as those considered in the previous studies
PTV2016 ; PY2018 ; PSTY2020  are more realistic in general). The asymptotic
analysis for the stiff chemotactic response with the finite tumbling duration
should be an important future work.

Hereafter, we write the nondimensional quantities without ”$\hat{\quad}$” for
simplicity unless otherwise stated.

###  3.1 Small adaptation-time scaling

We consider the small adaptation-time scaling in Eq. (5) as follows:

$\tau=\alpha\varepsilon,\quad\sigma=\varepsilon,$ (11)

where the parameter $\alpha=O(1)$ denotes the ratio of the adaptation time to
the mean running duration (i.e., $\alpha=\tau_{a}/\lambda^{-1}$).

Then, Eq. (5) is written as follows:

$\varepsilon^{2}\partial_{t}f_{\varepsilon}+\varepsilon
v\cdot\nabla_{x}f_{\varepsilon}+\partial_{m}\left\\{\left(\frac{M_{\varepsilon}-m}{\alpha}\right)f_{\varepsilon}\right\\}=\mu
g_{\varepsilon}-\Lambda(M_{\varepsilon}-m)f_{\varepsilon},$ (12a)
$\varepsilon^{2}\partial_{t}g_{\varepsilon}+{\color[rgb]{0,0,0}\partial_{m}}\left\\{\left(\frac{M_{\varepsilon}-m}{\alpha}\right)g_{\varepsilon}\right\\}=\Lambda(M_{\varepsilon}-m)<f_{\varepsilon}>-\mu
g_{\varepsilon},$ (12b)

where the subscript $\varepsilon$ represents the expansion of the quantity
with respect to $\varepsilon$, e.g., $f_{\varepsilon}=f_{0}+\varepsilon
f_{1}+\varepsilon^{2}f_{2}\cdots$. Here, we also write $M(S_{\varepsilon})$ as
$M_{\varepsilon}=M(S_{\varepsilon})=M_{0}+\varepsilon M_{1}+\cdots$.

The asymptotic analysis of the above equation gives the following standard KS
equation for the total population density $\rho$ at the continuum limit
$\varepsilon\rightarrow 0$ as follows:

$\sigma_{\nu}\partial_{t}\rho_{0}-\nabla_{x}\cdot
c_{d}\left[\nabla_{x}\rho_{0}+\frac{\Lambda^{\prime}(0)\alpha\rho_{0}}{1+\alpha}\nabla_{x}M_{0}\right]=0,$
(13)

where $\sigma_{\nu}$ is the time-scale parameter defined as
$\sigma_{\nu}=1+\nu$ and $c_{d}$ is the diffusion constant calculated as
$c_{d}=1/d$ for $d=1$, 2, and 3. Here, we note again that the parameter
$\nu=\mu^{-1}$ denotes the relative mean tumbling duration to the mean running
duration. The population densities of the running and tumbling cells are
obtained as follows:

$\displaystyle\rho_{f}=\frac{1}{1+\nu}\rho_{0},\quad\rho_{g}=\frac{\nu}{1+\nu}\rho_{0}.$
(14)

The formal derivation of the KS equation (13) is concisely described in
Appendix A.1.

Equation (13) shows that the tumbling duration $\nu$ only affects the time
scale, but the spatial distribution of the population density $\rho_{0}$ in
the steady state are not affected by the tumbling duration at the continuum
limit $\varepsilon\rightarrow 0$ in the small adaptation-time scaling (11).

###  3.2 Large adaptation-time scaling

We consider the large adaptation time scaling at Eq. (5) as follows:

$\tau=\beta/\varepsilon,\quad\sigma=\varepsilon,$ (15)

where the parameter $\beta=O(1)$ denotes the ratio of the adaptation time to
the diffusion time (i.e., $\beta=\tau_{a}/t_{d}$).

Then, Eq. (5) is written as follows:

$\varepsilon^{2}\partial_{t}f_{\varepsilon}+\varepsilon
v\cdot\nabla_{x}f_{\varepsilon}+\varepsilon^{2}\partial_{m}\left\\{\left(\frac{M_{\varepsilon}-m}{\beta}\right)f_{\varepsilon}\right\\}=\mu
g_{\varepsilon}-\Lambda(M_{\varepsilon}-m)f_{\varepsilon},$ (16a)
$\varepsilon^{2}\partial_{t}g_{\varepsilon}+\varepsilon^{2}{\color[rgb]{0,0,0}\partial_{m}}\left\\{\left(\frac{M_{\varepsilon}-m}{\beta}\right)g_{\varepsilon}\right\\}=\Lambda(M_{\varepsilon}-m)<f_{\varepsilon}>-\mu
g_{\varepsilon}.$ (16b)

The asymptotic analysis of Eq. (16) gives the continuum-limit equation at
$\varepsilon\rightarrow 0$ as follows:

$\partial_{t}h_{0}-\nabla_{x}\cdot\left[\frac{c_{d}}{\Lambda(M_{0}-m)}\nabla_{x}\left(\frac{h_{0}}{1+\nu\Lambda(M_{0}-m)}\right)\right]+\partial_{m}\left[\left(\frac{M_{0}-m}{\beta}\right)h_{0}\right]=0,$
(17)

where $h_{0}$ denotes the density of the cells with internal state $m$ and is
defined as follows:

$h_{0}(t,x,m)=<f_{0}>(t,x,m)+g_{0}(t,x,m).$

Although the above continuum-limit equation (17) was previously derived in
Ref. XXT2018 , we concisely describe the formal derivation of Eq. (17) in
Appendix A.2 for the completeness of the present paper. We also remark that
Eq. (17) with $\nu=0$ was also derived in Ref. Y2021  from the kinetic
transport equation without finite tumbling duration (9). However,
interestingly, the volcano effect was not observed in the previous study
while, as it will be seen in Sec. 4, the volcano effect arises both in MC
simulations and numerical results of Eq. (17) with the finite tumbling
duration $\nu\neq 0$. This indicates that the small modification introduced in
Eq. (17) with the parameter $\nu$ enables to produce the volcano effect.

The total population density of cells $\rho_{0}(t,x)$ is given by the
integration of $h_{0}$ with respect to the internal state $m$, i.e.,

$\rho_{0}(t,x)=\int_{\mathrm{R}}h_{0}(t,x,m)dm.$

In the following text, we call Eq. (17) the extended Keller-Segel (ExKS) model
because the consistency with the standard Keller-Segel model is confirmed at
$\beta\rightarrow 0$, as shown in the next subsection.

###  3.3 Consistency between the KS and ExKS models

It is easily seen that by taking the limit as $\alpha\rightarrow\infty$ in Eq.
(13), the KS equation is written as follows:

$\sigma_{\mu}\partial_{t}\rho_{0}-\nabla_{x}\cdot
c_{d}\left[\nabla_{x}\rho_{0}+\Lambda^{\prime}(0)\rho_{0}\nabla_{x}M_{0}\right]=0.$

The same equation is also obtained by taking the limit $\beta\rightarrow 0$ in
Eq. (17). This can be seen as follows:

When taking the limit $\beta\rightarrow 0$ at Eq. (17), we have the equation
as follows:

$\partial_{m}\left[(M_{0}-m)h_{0}\right]=0.$

Thus, the solution $h_{0}$ at $\beta\rightarrow 0$ is written as follows:

$h_{0}(t,x,m)=\rho_{0}(t,x)\delta(M_{0}-m),$ (18)

where $\delta(m)$ is the Dirac delta function.

On the other hand, by integrating Eq. (17) with respect to $m$, we obtain the
equation as follows:

$\partial_{t}\rho_{0}-\nabla_{x}\cdot j_{1}=0,$ (19)

where

$j_{1}=c_{d}\int_{\mathbb{R}}B_{1}(t,x,m)dm,$ (20)

and $B_{1}$ is defined as Eq. (A.16b).

By substituting Eq. (18) into Eq. (A.16b), we obtain the equation as follows:

$\displaystyle B_{1}$
$\displaystyle=-\frac{\mu}{\Lambda(M_{0}-m)}\nabla_{x}\left(\frac{\rho_{0}\delta(M_{0}-m)}{\mu+\Lambda(M_{0}-m)}\right)$
$\displaystyle=-\frac{\mu\delta(M_{0}-m)\nabla_{x}\rho_{0}}{\Lambda(M_{0}-m)(\mu+\Lambda(M_{0}-m))}-\frac{\mu\rho_{0}\nabla_{x}M_{0}}{\Lambda(M_{0}-m)}\left[\frac{\delta^{\prime}(M_{0}-m)}{\mu+\Lambda(M_{0}-m)}-\frac{\Lambda^{\prime}(M_{0}-m)\delta(M_{0}-m)}{(\mu+\Lambda(M_{0}-m))^{2}}\right].$
(21)

Thus, the flux $j_{1}$ (Eq. (20)) at $\beta\rightarrow 0$ is written as
follows:

$\displaystyle j_{1}$ $\displaystyle=-\frac{\mu
c_{d}\nabla_{x}\rho_{0}}{\Lambda(0)(\mu+\Lambda(0))}-\mu
c_{d}\rho_{0}\nabla_{x}M_{{\color[rgb]{0,0,0}0}}\left[\frac{\Lambda^{\prime}(0)(\mu+\Lambda(0))+\Lambda(0)\Lambda^{\prime}(0)}{\Lambda^{2}(0)(\mu+\Lambda(0))^{2}}-\frac{\Lambda^{\prime}(0)}{\Lambda(0)(\mu+\Lambda(0))^{2}}\right]$
$\displaystyle=-\frac{\mu
c_{d}}{\Lambda(0)(\mu+\Lambda(0))}\left[\nabla_{x}\rho_{0}+\frac{\Lambda^{\prime}(0)}{\Lambda(0)}\rho_{0}\nabla_{x}M_{{\color[rgb]{0,0,0}0}}\right].$
(22)

By using $\Lambda(0)=1$ and $\nu=\mu^{-1}$, it is seen that Eq. (19) with Eq.
(3.3) provides the KS model at $\alpha\rightarrow\infty$. Due to this
consistency, we can say that the ExKS model (17) is an extension of the
standard KS model to involve the effects of the internal adaptation dynamics
at the large adaptation-time scaling.

##  4 Numerical analysis

###  4.1 Problem and method

We consider the aggregation of chemotactic bacteria in the one-dimensional
interval $-L/2\leq x\leq L/2$ and two-dimensional square
$(x_{1},x_{2})\in[-L/2,L/2]\times[-L/2,L/2]$ with the periodic boundary
condition under the following distribution of the external chemical cues
$S(x)$:

$S(x)={\color[rgb]{0,0,0}\left\\{\begin{array}[]{cc}\exp(-|x|),&(\mathrm{for}\quad
d=1),\\\ \exp(-r),&(\mathrm{for}\quad d=2),\end{array}\right.}$ (23)

where $r=\sqrt{x_{1}^{2}+x_{2}^{2}}$ is the radial distance in the two-
dimensional space $d=2$. We also consider the logarithmic sensing of the
external chemical cue KJTW2009 , where the local equilibrium of the internal
state is modeled as $M(S)=\log S$. Thus, for Eq. (23), $M(S)$ is described as
follows:

$M(S)={\color[rgb]{0,0,0}\left\\{\begin{array}[]{cc}-|x|,&(d=1),\\\
-r,&(d=2).\end{array}\right.}$ (24)

The modulation of tumbling frequency in the chemotactic response is determined
by the deviation of the internal state $m$ from the local equilibrium state as
described in Eq. (1). In the following numerical simulations, we consider the
response function $F(X)$ in Eq. (1), which is described as follows:

$F\left(X\right)=\frac{\chi X}{\sqrt{1+X^{2}}}.$ (25)

The kinetic transport equation (5) is numerically solved for one- and two-
dimensional spaces by the Monte Carlo method described in Appendix B, while
the ExKS model (17) is calculated only for one-dimensional space by a standard
finite-volume (FV) scheme, which is similar to that described in Appendix B.1
in Ref. Y2021 . In the MC simulations, we set the number of mesh intervals as
$I=100$, the total number of MC particles as $N=720000$, and the time step
size as $\Delta t=2\times 10^{-4}$ for one-dimensional problem. For two-
dimensional problem, we set the number of lattice cells as $50\times 50$
(i.e., $I=50$ in Appendix B), the total number of MC particles as $N=1.8\times
10^{7}$, and the time step size as $\Delta t=2\times 10^{-4}$.  Here, it
should be noted that in the MC method, the characteristic time
$t_{\mathrm{c}}$ is fixed as $t_{\mathrm{c}}=L_{\mathrm{c}}/\mathrm{v_{c}}$
independent of $\varepsilon$, which reads $\sigma=1$ at Eq. (5), but the time
period of the MC simulation is varied according to $\varepsilon$ as $0\leq
t\leq\frac{2L^{2}}{\varepsilon}$ in order to compare the MC results with the
continuum-limit models obtained in the diffusion time scaling, i.e.,
$\sigma=\varepsilon$. The MC results are also time-averaged over the time
period $\delta t=\frac{0.1L^{2}}{\varepsilon}$. In the FV method of the ExKS
model, we set the number of mesh intervals as $I=100$ in space and $K=800$ in
the internal state, where the internal state variable $m$ is discretized as
$m_{k}=-Y+k\Delta y$ ($k=0,\cdots,K$) with $\Delta y=2Y/K$ and $Y=5$. The
time-step size is set as $\Delta t=1\times 10^{-4}$, and the simulations are
performed over the time period $0<t\leq 25$.

Numerical simulations are performed for various values of the adaptation time
$\tau$, run length $\varepsilon$, stiffness $\delta$, and modulation amplitude
$\chi$, while the spatial extent $L=10$ is fixed.  In this paper, most of the
numerical results are given for the spatially one-dimensional problem. Only a
few results are presented for the two-dimensional problem in order to
demonstrate that the volcano effect is certainly reproduced by the kinetic
transport equation (5) in two-dimensional space.

###  4.2 Occurrence of bimodal aggregation

Figure 1: Spatial distributions of the total population density $\rho$ in the
steady state at different adaptation times $\tau$ for the one-dimensional
problem. Figure (a) shows the results at the small adaptation-time scaling,
i.e., $\tau=\alpha\varepsilon$, and Figure (b) shows the results at the large
adaptation-time scaling, i.e., $\tau=\beta/\varepsilon$. The parameters
$\varepsilon=0.1$, $\nu=0.3$, $\delta=1.25$, and $\chi=0.7$ are fixed. The
colored solid lines show the results of the MC simulations in both figures,
while the black dashed lines, which almost overlap with the colored solid
lines, show the results of the KS model (13) in Fig. (a) and those of the ExKS
model (17) in Fig. (b).  Figure 2: Spatial distributions of the population
densities of running and tumbling cells in the steady states, i.e., $\rho_{f}$
[in (a)] and $\rho_{g}$ [in (b)] at the large adaptation time scaling
$\tau=\beta/\varepsilon$ with $\beta$=0.2, 0.5, 1.0, and 2.0 for the one-
dimensional problem. The parameter values of $\varepsilon$, $\nu$, $\delta$,
and $\chi$ are the same as those in Fig. 1. For the line types, see the
caption in Fig. 1.

Figure 1 shows the comparison of the aggregation profiles of the total
population density $\rho$ at different values of the adaptation time $\tau$.
It is clear that there is a nonmonotonic dependence of the adaptation time on
the aggregation profile between the small and the large adaptation-time
scalings; at the small adaptation-time scaling $\tau=O(\varepsilon)$ [in Fig
1(a)], the peak of the aggregation profile increases as $\tau$, while at the
large adaptation-time scaling $\tau=O(\varepsilon^{-1})$ [in Fig. 1(b)], it
decreases as $\tau$.

This nonmonotonic dependence was also observed in a previous study Y2021 ,
where the self-organized aggregation of bacteria without the tumbling phase
(i.e., $\nu=0$) was considered. However, a remarkable difference between the
cases with and without the tumbling phase can be observed at the large-
adaptation time scaling. In Fig. 1(b), the volcano effect (i.e., bimodal
aggregation) is observed at $\beta=$0.5, 1, and 2.

To clarify the contribution of the tumbling cells to bimodal aggregation, we
show the spatial distributions of the running and tumbling cells in Fig. 2(a)
and Fig. 2(b), respectively. It is evident that at the large adaptation-time
regime, the distribution of running cells $\rho_{f}$ is unimodal, while that
of the tumbling cells $\rho_{g}$ is bimodal, where the hollow at the central
region becomes increasingly larger as $\beta$ increases. At the small
adaptation-time regime, say $\tau\lesssim 0.1$, the distribution of tumbling
cells is unimodal, which can be seen from the red-colored triangles in Fig.
3(a). Thus, the bimodal distribution of the total population density $\rho$
observed in Fig. 1(b) is due to the local decrease of the tumbling cells at
the central region.

Figure 3: The diagrams of bimodal aggregation with respect to the adaptation
time $\tau$ [in (a)] and the tumbling duration $\nu$ [in (b)] for the one-
dimensional problem. The second-order derivatives of the total population
density $\rho$ and the population density of the tumbling cells $\rho_{g}$ at
$x=0$, $\rho^{\prime\prime}$ and $\rho^{\prime\prime}_{g}$ are shown. Thus,
$\rho^{\prime\prime}>0$ and $\rho^{\prime\prime}_{g}>0$ represent the
occurrence of bimodal aggregation. In both figures, the parameters
$\varepsilon=0.1$, $\delta=1.25$, and $\chi=0.7$ are commonly fixed, while
$\nu=0.3$ is set in (a) and $\tau=10.0$ is set in (b). The results obtained by
the MC simulation and ExKS model are shown in each figure.

Figure 3 is the diagram of the occurrence of bimodal aggregation with respect
to the adaptation time $\tau$ and the tumbling duration $\nu$. Here, the
second derivatives of the total population density at $x=0$ are calculated
from the numerical results as follows:

$\rho^{\prime\prime}=\frac{\rho_{\frac{I}{2}+1}-\rho_{\frac{I}{2}}-\rho_{\frac{I}{2}-1}+\rho_{\frac{I}{2}-2}}{(\Delta
x)^{2}},$

where $\rho_{i}$ is the averaged population density in the interval
$x\in[x_{i},x_{i+1}]$ with $x_{i}=i\Delta x$ and $x_{I}=L$.

It can be seen that bimodal aggregation occurs when both the adaptation time
$\tau$ and the tumbling duration $\nu$ are sufficiently large, i.e.,
$\tau\gtrsim 5$ and $\nu\gtrsim 0.1$. The effect of the tumbling duration
$\nu$ on the bimodal aggregation is also seen in SI 3 in the supplemental
information (see Online Resource, SI.pdf). It can be seen that the hollow at
the central region becomes increasingly deeper as $\nu$ increases.

Notably, it is also seen that the ExKS model (17) can reproduce the transient
behavior between unimodal and bimodal aggregations described by the kinetic
transport equation (although the deviation of the ExKS model from the MC
results increases as the tumbling duration increases). This indicates that the
ExKS model (17) inherits some essential mathematical structure that should be
necessary to describe the volcano effect.

Figure 4 shows the MC results for the two-dimensional problem. It is clearly
seen that the volcano effect arises in the two-dimensional space at the large
adaptation-time regime. As in the one-dimensional problem, the population
density of running cells is unimodal while that of tumbling cells is bimodal.
The spatial distribution of the local mean run length $\bar{\xi}$, which is
defined as Eq. (26), sharply increases as the radial distance $r$ approaches
to the center of the aggregation. This feature is also similar to that
observed in the local mean run length in the one-dimensional problem, which
will be discussed in detail in the next section.

Figure 4:  Volcano aggregation in two-dimensional space. Spatial distributions
of the total population density $\rho$ (in (a)), population density of running
cells $\rho_{f}$ (in (b)), population density of tumbling cells $\rho_{g}$ (in
(c)), and local mean run length $\bar{\xi}$ in the two-dimensional square with
$L=10$ are shown. Each inset show the y-distributions at the intersection
$x=5$ (which are shown in sold black lines on the surface of the two-
dimensional distributions). The parameters are set as $\varepsilon=0.1$,
$\tau=10$, $\nu=0.3$, $\delta=0.1$, and $\chi=0.9$.

We note that MC results for the two-dimensional problem at different values of
$\varepsilon$ and $\tau$ are also given in SI 4-6 in the supplemental
information (see Online Resource, SI.pdf).

###  4.3 Why volcano occurs

In this subsection, we further investigate the microscopic mechanism of the
volcano effect with using the numerical results for the one-dimensional
problem. Figure 5 shows the spatial distributions of the local mean run length
of the bacteria defined by the equation as follows:

$\xi^{\pm}(t,x)=\int_{R}\frac{\varepsilon}{\Lambda(M(S)-m)}<f>^{\pm}(t,x,m)dm,$
(26)

where

$<f>^{\pm}=\frac{2}{||V||}\int_{\pm v\cdot\nabla_{x}S>0}f(v)dv.$

Thus, $\xi^{\pm}$ represent the local mean run lengths of the bacteria,
climbing and descending the gradient of chemical cues, respectively, and
$\bar{\xi}=\frac{\xi^{+}+\xi^{-}}{2}$ is the local mean run length for all of
the moving bacteria.

Figure 5: Spatial distributions of the local mean run length defined by Eq.
(26) at $\tau=1$ (a), $\tau=5$ (b) and $\tau=10$ (c). The parameters
$\varepsilon=0.1$, $\nu=0.3$, $\delta=0.25$, and $\chi=0.7$ are fixed. The
downward arrows on the horizontal axis in each figure show the position at
$x=\sqrt{\varepsilon\tau}$, which represents the thickness of the diffusion
layer within the adaptation time. The colored lines show the results obtained
by the MC simulations, while the black dashed lines show the results of
$\bar{\xi}$ obtained by the ExKS model (17).

Since the local equilibrium of internal state $M(S)$ is described as Eq. (24),
the spatial gradient of $M(S)$, which is sensed by moving bacteria along their
pathway, is constant except the boundaries at $x=0$ and $x=\pm 5$, but it
changes stepwise at the boundaries. The mean run length substantially changes
only within the layers near the boundaries with thickness proportional to the
diffusion length of the population density during the adaptation time, i.e.,
$|x|\lesssim{\color[rgb]{0,0,0}\sqrt{\varepsilon\tau}}$ and
$|x-5|\lesssim{\color[rgb]{0,0,0}\sqrt{\varepsilon\tau}}$. We call this layer
the diffusion layer since it is created due to the diffusion coupled with the
internal adaptation dynamics as explained below. Here, we note that the non-
dimensional form of the diffusion constant defined below Eq. (10) is written
as $\hat{D}_{\rho}=\varepsilon$. Thus, $\sqrt{\varepsilon\tau}$ represents the
thickness of the diffusion layer within the adaptation time.

In the middle region except the diffusion layers near the boundaries, the
bacteria uniformly create a biased random motion according to the local
gradient of $M(S)$, where the bacteria climbing the gradient are more likely
to have longer run lengths (see $\xi^{+}$ in Fig. 5), while those descending
the gradient are more likely to have shorter run lengths (see $\xi^{-}$ in
Fig. 5). However, in the close vicinity of the boundary at $x=0$, the bacteria
moving across the boundary at $x=0$ with positive velocity are more likely to
have longer run lengths than the bacteria moving with positive velocity in the
middle region, since they have been climbing the gradient of $M(S)$ in $x<0$.
Thus, the bacteria coming from the left side of the boundary contribute to
increasing the local mean run length in the close vicinity of the boundary at
$x=0$.

Furthermore, since the bacteria keep the memory of the internal state during
the adaptation time $\tau$, which is comparable to the diffusion time in the
characteristic length at the large adaptation-time scaling (see the second
paragraph in Sec. 3), the population of the bacteria having a longer run
length in the close vicinity at $x=0$ diffuses via the random motions of
individual bacteria during the adaptation time. Thus, the diffusion layer with
thickness $\sqrt{\varepsilon\tau}$ is formed near the boundary at $x=0$.
Notably, the increase in the local mean run length in the diffusion layer
$|x|<{\color[rgb]{0,0,0}\sqrt{\varepsilon\tau}}$ reduces the local population
of the tumbling cells $\rho_{g}$, as shown in Fig. 2(b).

The same argument also holds in the right-side diffusion layer; the bacteria
moving across the boundary at $x$=5 with a negative velocity are more likely
to have shorter run lengths than those of the bacteria moving with negative
velocity in the middle region since they have been descending the gradient of
the chemical cue in $x>5$ to ensure that they contribute to decreasing the
local mean run length. Thus, as shown in Fig. 2, the population of tumbling
cells slightly increases in the diffusion layer at $|x|=5$.

The dependence of the adaptation time on the spatial profile is more clearly
observed in Fig. 6, where the spatial distributions of $\rho$, $\rho_{g}$, and
$\bar{\xi}$ are shown in the scaled coordinate $x/\sqrt{\beta}$. Notably, the
consistency of this scaling property with the ExKS model (17) is confirmed;
that is, when the equilibrium of the internal state is the linear function
described as Eq. (24), the solution of Eq. (17) has the scaling property as
follows:

$h_{\delta,\beta}(t,x,m)\propto h_{a\delta,a^{2}\beta}(a^{2}t,ax,am),$ (27)

where $a>0$ is an arbitrary constant and $h_{\delta,\beta}(t,x,m)$ is the
solution of Eq. (17) with the stiffness $\delta$ and the parameter of the
adaptation time $\beta$. Thus, setting $a=1/\sqrt{\beta}$, we obtain the
equation as follows:

$h_{\delta,\beta}(t,x,m)\propto
h_{\frac{\delta}{\sqrt{\beta}},1}\left(\frac{t}{\beta},\frac{x}{\sqrt{\beta}},\frac{m}{\sqrt{\beta}}\right).$
(28)

By integrating the above relation with respect to $m$, we obtain the equation
as follows:

$\rho_{\delta,\beta}(x)=\rho_{\frac{\delta}{\sqrt{\beta}},1}\left(\frac{x}{\sqrt{\beta}}\right).$
(29)

This scaling property indicates that the spatial profile is mainly determined
by the parameter $\beta$ when the stiffness $\delta$ does not substantially
affect the spatial profile. Indeed, as seen in SI 1 and SI 2 in the
supplemental information (see Online Resource, SI.pdf), the variation in the
stiffness $\delta$ does not considerably affect the spatial scale of the
distribution of population density. Thus, the scaling property (29) and the
numerical results in SI 1 and SI 2 are consistent with the observation about
the dependence of $\beta$ on the spatial profile in Fig. 6.

One may think that the mean run length of the bacteria $\varepsilon$ should
directly affect the volcano profile. Figure 7 shows the effect of the mean run
length $\varepsilon$ on the volcano at different values of $\beta$. The
aggregation profile becomes more diffusive as the mean run length
$\varepsilon$ increases. However, the peak position of the aggregation is not
affected by the mean run length $\varepsilon$, but it is clearly affected by
the parameter $\beta$, as is already shown in Fig. 6 and is explained by the
scaling property (29). Thus, due to the above numerical results, we conclude
that the volcano observed at the large adaptation-time scaling is generated
due to the coupling of the internal adaptation and the diffusive motion by the
runs and tumbles of bacteria.

Figure 6: Spatial distributions of the total population density $\rho$ (in
(a)), population density of tumbling cells $\rho_{g}$ (in (b)), and local mean
run length $\bar{\xi}$ in the re-scaled spatial axis $x/\sqrt{\beta}$ at the
large adaptation time scaling $\tau=\beta/\varepsilon$. The parameters
$\nu=0.3$, $\chi=0.7$ and $\delta=0.25$ are fixed. The solid lines show the
results of the ExKS model, while the dashed line shows the results obtained by
MC simulations at $\varepsilon=0.1$.  Figure 7: The effect of the run length
$\varepsilon$ on the bimodal aggregation and the asymptotic behavior to the
continuum limit $\varepsilon\rightarrow 0$. The solid lines show the results
obtained by the MC simulations, where the adaptation time $\tau$ is set as
$\tau=\beta/\varepsilon$ with $\beta=0.2$ in (a), $\beta=0.5$ in (b), and
$\beta=1$ in (c), while the black dashed lines show the results of the ExKS
model (17). The parameters $\nu=0.3$, $\chi=0.7$ and $\delta=0.25$ are fixed.

It can also be confirmed from the numerical results in Fig. 7 that the ExKS
model is a good approximation of the kinetic transport model (16) at the large
adaptation-time scaling $\tau=\beta/\varepsilon$ when $\varepsilon$ is
moderately small, say $\varepsilon\lesssim 0.2$.

##  5 Concluding remarks

The volcano effect of run-and-tumble chemotactic bacteria was investigated via
MC simulations based on the kinetic transport equation (5), which considers a
non-instantaneous tumbling as well as an internal adaptation dynamics.  The MC
simulations were performed for one- and two-dimensional problems for a wide
range of parameters. Especially, in order to identify the parameter regime and
scaling for the volcano effect to arise, the results of one-dimensional MC
simulations were compared with the numerical results of continuum-limit
equations obtained at different scalings of the adaptation time.

MC simulations uncovered that the distribution of running cells is always
unimodal irrespective to the adaptation time $\tau$, while those of tumbling
cells and total population density become bimodal when both the adaptation
time $\tau$ and the tumbling duration $\nu$ are sufficiently large, i.e.,
$\tau=O(1/\varepsilon)$ and $\nu>0.1$.

In order to clarify the microscopic behavior of bacteria forming volcano, the
distribution of local mean run length was investigated. It is clarified that
the decrease of the tumbling cells at the central region is caused by the
increase of the local mean run length in the central region, where the
diffusion layer is created due to the coupling of diffusion and internal
adaptation of bacteria. More concretely, the local mean run length of bacteria
increases at the close vicinity of the boundary at $x=0$ due to the mirror
symmetry of chemical cues. The bacteria who have longer run lengths at the
central region diffuse due to their individual run-and-tumble motions with
keeping the memory of their internal state during the adaptation time. Thus,
the layer in which the local run length is larger than that in the middle
region extends to the thickness of the diffusion length in the adaptation time
$\sqrt{\varepsilon\tau}$. The increase of the mean run length indicates the
decrease of the tumbling cells, so that the tumbling cells considerably
decrease in the diffusion layer and form volcano when the adaptation time is
comparable to the diffusion time. Thus, the coupling of diffusion, adaptation,
and tumbling duration is crucial for the formation of volcano.

The fact that the ExKS model (17), which is obtained at the large adaptation-
time scaling $\tau=O(1/\varepsilon)$, well approximates the volcano
aggregations obtained by MC simulations also convinces that the diffusion of
the memory is crucial for the formation of volcano in run-and-tumble bacteria.
The scaling property found in the ExKS model (29) also clearly shows that the
spatial profile of the population density is mainly determined by the
diffusion length coupled with the internal adaptation dynamics.

One may think that such a large adaptation time $\tau=O(1/\varepsilon)$ is
biologically unrealistic. However, in the experiment of volcano MBBO2003 , the
system size is so small as $L_{c}\sim 100\,\mu\mathrm{m}$, so that the non-
dimensional run length $\varepsilon$ is estimated as $\varepsilon=0.2$ by
using the typical run length of E. coli (i.e., $v_{c}/\lambda$=20 [$\mu$m]).
Furthermore, by using the typical quantities of run duration $\lambda^{-1}$=1
[s], tumbling duration $\mu^{-1}$=0.1 [s], and adaptation time $\tau_{a}$=10
[s] for E. coli, the other non-dimensional parameters defined in Eq. (7) are
calculated as $\sigma$=0.2, $\tau$=2, and $\nu$=0.1, where $t_{c}=t_{d}$ is
used in the calculation of $\sigma$. Thus, the large adaptation-time scaling
$\tau=O(1/\varepsilon)$ is not unrealistic but is rather relevant for volcano.

In conclusion we remark that diffusion, adaptation, and non-instantaneous
tumbling are key ingredients to describe the non-monotonic aggregation of
chemotactic bacteria observed at scales of some tens to hundreds of micro-
meters such as volcano, and the kinetic transport model, which considers the
non-instantaneous tumbling as well as the internal adaptation dynamics would
be useful to elucidate the mathematics behind the complicated aggregation
behaviors.

##  Appendix A Derivation of the continuum-limit model

The continuum-limit models, which are utilized to compared with the MC results
in the main text, were previously derived in literatures, e.g., Refs. CMPS2004
; EO2004 ; HO2000 ; PTV2016 ; PSTY2020 ; XXT2018 . In this appendix, we
briefly describe the derivation of the models for the completeness of this
paper.

###  A.1 Small adaptation-time scaling

By following the procedure in Ref. PSTY2020 , we change the variables of Eq.
(12) as follows:

$f_{\varepsilon}(t,x,v,m)=p_{\varepsilon}\left(t,x,v,y=\frac{M(S_{\varepsilon})-m}{\varepsilon}\right),\quad
g_{\varepsilon}(t,x,m)=q_{\varepsilon}\left(t,x,y=\frac{M(S_{\varepsilon})-m}{\varepsilon}\right).$

Then, we have

$\varepsilon^{2}\partial_{t}p_{\varepsilon}+\varepsilon
v\cdot\nabla_{x}p_{\varepsilon}+\partial_{y}\left\\{\left(v\cdot
G_{\varepsilon}-\frac{y}{\alpha}\right)p_{\varepsilon}\right\\}=\mu
q_{\varepsilon}-\Lambda(\varepsilon y)p_{\varepsilon},$ (A.1a)
$\varepsilon^{2}\partial_{t}q_{\varepsilon}+\partial_{y}\left\\{-\frac{y}{\alpha}q_{\varepsilon}\right\\}=\Lambda(\varepsilon
y)<p_{\varepsilon}>-\mu q_{\varepsilon},$ (A.1b)

where we write $G_{\varepsilon}=\nabla_{x}M(S_{\varepsilon})$. By integrating
the sum of the above equations with respect to $v$ and $y$, we obtain the
following conservation law:

$\partial_{t}\rho_{\varepsilon}+\nabla_{x}\cdot\left(\frac{j_{\varepsilon}}{\varepsilon}\right)=0,$
(A.2)

where $\rho_{\varepsilon}$ is the total population density, i.e.,
$\rho_{\varepsilon}=\rho_{p_{\varepsilon}}+\rho_{q_{\varepsilon}}$ with
$\rho_{p_{\varepsilon}}=\int<p_{\varepsilon}>dy$ and
$\rho_{q_{\varepsilon}}=\int q_{\varepsilon}{\color[rgb]{0,0,0}d}y$, and the
flux $j_{\varepsilon}$ is defined as

$j_{\varepsilon}=\int<vp_{\varepsilon}>dy.$ (A.3)

As can be seen below, the KS equation is obtained from Eq. (A.2) at the
continuum limit ($\varepsilon\rightarrow 0$).

We assume $p_{\varepsilon}$ and $q_{\varepsilon}$ are compactly supported with
respect to $y$ and $\Lambda(\varepsilon y)$ in Eq. (A.1) can be expanded as
$\Lambda(\varepsilon y)=1+\varepsilon\Lambda^{\prime}(0)y+O(\varepsilon^{2})$.
Then, from the leading-order terms of Eq. (A.1), we obtain the leading-order
equation as follows:

$\partial_{y}\left\\{\left(v\cdot
G_{0}-\frac{y}{\alpha}\right)p_{0}\right\\}=\mu q_{0}-p_{0},$ (A.4a)
$\partial_{y}\left\\{-\frac{y}{\alpha}q_{0}\right\\}=<p_{0}>-\mu q_{0}.$
(A.4b)

By integrating each of Eqs. (A.4) w.r.t $y$, we obtain the following relation
between the leading-order population densities:

$\int p_{0}dy=\mu\rho_{q_{0}}=\rho_{p_{0}}=\frac{\mu}{1+\mu}\rho_{0}.$ (A.5)

Furthermore, by taking the moment of the above equation multiplied by $v$, we
obtain the flux $j_{0}$ as follows:

$j_{0}=\int<vp_{0}>dy=0.$ (A.6)

From the first-order terms of Eq. (A.1), we have the equation as follows:

$\nabla_{x}\cdot(vp_{0})+\partial_{y}\left\\{\left(v\cdot
G_{0}-\frac{y}{\alpha}\right)p_{1}\right\\}+\partial_{y}(v\cdot
G_{1}p_{0})=\mu q_{1}-p_{1}-\Lambda^{\prime}(0)yp_{0},$ (A.7a)
$\partial_{y}\left\\{\left(-\frac{y}{\alpha}\right)q_{1}\right\\}=<p_{1}>-\mu
q_{1}+\Lambda^{\prime}(0)y<p_{0}>.$ (A.7b)

By integrating the sum of above equations multiplied by $v$ with respect to
$v$ and $y$, we obtain

$\begin{split}\nabla_{x}\cdot\int<v\otimes
vp_{0}>dy=-j_{1}-\Lambda^{\prime}(0)\int y<vp_{0}>dy,\\\
\frac{\mu}{1+\mu}c_{d}\nabla_{x}\rho_{0}=-j_{1}-\Lambda^{\prime}(0)\int
y<vp_{0}>dy,\end{split}$

where we use Eq. (A.5) and $<v\otimes v>=c_{d}I$ with $c_{d}=1/d$ for
$d=1,2,3$. Here, $I$ is the identity matrix.  The last term of the above
equation is obtained by integrating Eq. (A.4a) multiplied by $vy$ with respect
to $v$ and $y$:

$\begin{split}-\int y\partial\left\\{<v\otimes vp_{0}>\cdot
G-\frac{y}{\alpha}<vp_{0}>\right\\}dy=\int y<vp_{0}>dy,\\\ \int<v\otimes
vp_{0}>dy\cdot G-\frac{1}{\alpha}\int y<vp_{0}>dy=\int y<vp_{0}>dy,\\\ \int
y<vp_{0}>dy=\frac{\mu}{1+\mu}\frac{\alpha}{1+\alpha}c_{d}\rho_{0}G.\end{split}$

Hence, the flux $j_{1}$ is written as follows:

$j_{1}=-\frac{\mu}{1+\mu}c_{d}\left[\nabla_{x}\rho_{0}+\frac{\alpha}{1+\alpha}\Lambda^{\prime}(0)G\rho_{0}\right].$
(A.8)

Thus, by taking the limit $\varepsilon\rightarrow 0$ at Eq. (A.2), we obtain
the KS equation (13).

###  A.2 The Extended Keller-Segel model

The derivation of Eq. (17) is as follows. In the following, we write the
average of $f_{\varepsilon}$ over the velocity space as
$A_{\varepsilon}=<f_{\varepsilon}>$.

From the leading-order terms of Eq. (16), we can write the leading-order
solution as follows:

$f_{0}=A_{0}(t,x,m),\quad g_{0}=\frac{\Lambda(M_{0}-m)}{\mu}A_{0}(t,x,m).$
(A.9)

Here, $A_{0}(t,x,m)$ is an unknown function independent of the velocity $v$.

From the $\varepsilon^{1}$ terms of Eq. (16), we obtain

$v\cdot\nabla_{x}A_{0}=\mu
g_{1}-\Lambda(M_{0}-m)f_{1}-\Lambda^{\prime}(M_{0}-m)M_{1}A_{0},$ (A.10a)
$0=\Lambda(M_{0}-m)A_{1}+\Lambda^{\prime}(M_{0}-m)M_{1}A_{0}-\mu g_{1}.$
(A.10b)

By taking the sum of the above equations, we obtain the following equation,

$v\cdot\nabla_{x}A_{0}=\Lambda(M_{0}-m)(A_{1}-f_{1}).$

Hence, $f_{1}$ can be written in the form

$f_{1}=A_{1}(t,x,m)+v\cdot B_{1}(t,x,m),$ (A.11a) with
$B_{1}=-\frac{\nabla_{x}A_{0}}{\Lambda(M_{0}-m)}.$ (A.11b)

From Eq. (A.10b), we can also write $g_{1}$ as follows:

$g_{1}=\frac{1}{\mu}\left(\Lambda(M_{0}-m)A_{1}+\Lambda^{\prime}(M_{0}-m)M_{1}A_{0}\right).$
(A.12)

Subsequently, from the $\varepsilon^{2}$ terms of Eq. (16), we obtain the
equation as follows:

$\begin{split}\partial_{t}A_{0}+v\cdot\nabla_{x}(A_{1}+v\cdot
B_{1})+\partial_{m}\left\\{\left(\frac{M_{0}-m}{\beta}\right)A_{0}\right\\}=\mu
g_{2}-\Lambda(M_{0}-m)f_{2}\\\ -\Lambda^{\prime}(M_{0}-m)M_{1}(A_{1}+v\cdot
B_{1})-\frac{1}{2}\Lambda^{\prime\prime}(M_{0}-m)M_{2}A_{0},\end{split}$
(A.13a)
$\begin{split}\frac{\Lambda(M_{0}-m)}{\mu}\partial_{t}A_{0}+\partial_{m}\left\\{\left(\frac{M_{0}-m}{\beta}\right)\frac{\Lambda(M_{0}-m)}{\mu}A_{0}\right\\}=\Lambda(M_{0}-m)A_{2}-\mu
g_{2}\\\
+\Lambda^{\prime}(M_{0}-m)M_{1}A_{1}{\color[rgb]{0,0,0}+}\frac{1}{2}\Lambda^{\prime\prime}(M_{0}-m)M_{2}A_{0}.\end{split}$
(A.13b)

By integrating the sum of the above equations with respect to $v$, we obtain
the equation to determine the leading-order solution $A_{0}$ as follows:

$\left(1+\frac{\Lambda(M_{0}-m)}{\mu}\right)\partial_{t}A_{0}+\nabla_{x}\cdot(c_{d}B_{1})+\partial_{m}\left\\{\left(\frac{M_{0}-m}{{\color[rgb]{0,0,0}\beta}}\right)\left(1+\frac{\Lambda(M_{0}-m)}{\mu}\right)A_{0}\right\\}=0.$
(A.14)

When we write the total density of cells with an internal state $m$ as
$h_{0}=<f_{0}>+g_{0}$, i.e., from Eq. (A.9),

$h_{0}(t,x,m)=\left(1+\frac{\Lambda(M_{0}-m)}{\mu}\right)A_{0}(t,x,m),$ (A.15)

we can rewrite Eq. (A.14) as follows:

$\partial_{t}h_{0}+\nabla_{x}\cdot(c_{d}B_{1})+\partial_{m}\left\\{\left(\frac{M_{0}-m}{{\color[rgb]{0,0,0}\beta}}\right)h_{0}\right\\}=0,$
(A.16a) with
$B_{1}=-\frac{1}{\Lambda(M_{0}-m)}\nabla_{x}\left(\frac{h_{0}}{1+\frac{\Lambda(M_{0}-m)}{\mu}}\right).$
(A.16b)

Thus, we obtain the ExKS model (17).

##  Appendix B Monte Carlo method

We extend the Monte Carlo (MC) method developed in Refs. Y2017 ; PY2018 ;
Y2021  to include the tumbling duration. In the MC method, we use the time
scale $t_{0}=L_{0}/V_{0}$ (which reads $\sigma=1$ in Eq. (5)). We also use the
variable $y=M(S)-m$ instead of the internal state $m$ itself, since this
change of variable rewrites the internal adaptation dynamics into the
formulation involving the material derivative of the chemical cue along the
moving pathway of bacteria, which can be calculated in a straightforward way
in the MC code.

By changing the variable as $f(t,x,v,m)=p(t,x,v,y=M(S)-m)$ and
$g(t,x,m)=q(t,x,y=M(S)-m)$, Eq. (5) reads as follows:

$\partial_{t}p+v\cdot\nabla_{x}p+\partial_{y}\left\\{\left(D_{t}M(S)-\frac{y}{\tau}\right)p\right\\}=\frac{1}{\varepsilon}\left[\mu
q-\Lambda(y)p\right],$ (B.1a)
$\partial_{t}q+\partial_{y}\left\\{\left(\partial_{t}M(S)-\frac{y}{\tau}\right)q\right\\}=\frac{1}{\varepsilon}\left[\Lambda(y)<p>-\mu
q\right],$ (B.1b)

where $D_{t}$ denotes the material derivative along the moving pathway, i.e.,
$D_{t}=\partial_{t}+v\cdot\nabla_{x}$. The $y$-derivative term in Eq. (B.1)
denotes the internal dynamics of individual cells, where the internal state
$y$ changes according to the temporal variation of the external chemical cue
sensed by individual cells along their moving pathways, i.e., for the
logarithmic sensing $M(S)=\log S$, the equation is as follows:

$\dot{y}=\frac{D_{t}S}{S}-\frac{y}{\tau}.$ (B.2)

For the one-dimensional problem ($d=1$), the velocity space is composed of
three discrete velocities, i.e., $v=\\{-1,0,1\\}$ and the one-dimensional
space $-L/2\leq x\leq L/2$ is divided into the uniform mesh cells
$C_{i}=[i\Delta x,(i+1)\Delta x]$ ($i=0,\cdots,I-1$), where $\Delta x=L/I$ is
the width of the mesh interval.  For the two-dimensional problem ($d=2$), the
velocity space is composed of $v=0$ and $v=(\cos\theta,\sin\theta)$, where the
angle $\theta$ is uniformly distributed in $[0,2\pi]$, and the two-dimensional
space $[-L/2,L/2]\times[-L/2,L/2]$ is divided into the uniform lattice mesh
cells $C_{i}=[i_{1}\Delta x,(i_{1}+1)\Delta x]\times[i_{2}\Delta
x,(i_{2}+1)\Delta x]$ ($i_{1},i_{2}=0,\cdots,I-1$) with $i=i_{1}+i_{2}\times
I$.  Initially, MC particles are uniformly distributed in each mesh cell with
the equilibrium internal state at $y=0$. The initial velocities of each MC
particle are randomly determined from $v=\\{-1,+1\\}$ for $d=1$ while, for
$d=2$, the direction $\theta$ of the individual initial velocities
$v=(\cos\theta,\sin\theta)$ are uniformly randomly determined.  Thus, the
initial condition for Eq. (B.1) is described as $p=\delta(y)$ and $q=0$ at
$t=0$.

Then, the position $r^{k}_{l}$, velocity $v_{l}^{k}$, and internal state
$y^{k}_{l}$ of the $l$th MC particle at time $t=k\Delta t$ ($k\geq 1$) are
determined as follows:

  1. 1.

Each MC particle moves as follows:

$r^{k}_{l}=r^{k-1}_{l}+v_{l}^{k-1}\Delta t.$ (B.3)

  2. 2.

The population density in the $i$th mesh cell $C_{i}$, $\rho^{k}_{i}$ are
calculated as follows:

$\rho_{i}^{k}=\frac{1}{\bar{N}}\sum_{l=0}^{N}\int_{C_{i}}\delta(x-r_{l}^{k})dx,$
(B.4)

where $\bar{N}$ is the average number of the MC particles in the mesh cell,
i.e., $\bar{N}=N/I$ for $d=1$ and $\bar{N}=N/I^{2}$ for $d=2$.

  3. 3.

The internal state of the $l$th MC particle, $y^{k}_{l}$, is updated according
to Eq. (B.2) as follows:

$\frac{y^{k}_{l}-y^{k-1}_{l}}{\Delta t}=\frac{S^{k}_{l}-S^{k-1}_{l}}{\Delta
tS^{k-1}_{l}}-\frac{y^{k}_{l}}{\tau},$ (B.5)

where $S^{k}_{l}$ denotes the local concentration of the chemical cue at the
position of the $l$th MC particle, i.e., $S^{k}_{l}=-\exp|r_{l}^{k}|$. We note
that in Eq. (B.5), the pathway derivative $D_{t}S$ in Eq. (B.2) is given by
the rate of change of $S$ sensed by each MC particle, i.e.,
$(S_{l}^{k}-S_{l}^{k-1})/\Delta t$.

  4. 4.

The running or tumbling state of the $l$th MC particle is stochastically
determined according to the right-hand side of Eq. (B.1). When
$v^{k-1}_{l}\neq 0$, the velocity is changed to $v^{k}_{l}=0$ by the
probability $\frac{\Delta t\Lambda(y_{l}^{k})}{\varepsilon}$, while when
$v^{k-1}_{l}=0$, the velocity is changed to $|v^{k}_{l}|=1$ by the probability
$\frac{\Delta t\mu}{\varepsilon}$, where the moving direction, i.e.,
$v=\\{-1,+1\\}$ for $d=1$ and the angle $\theta\in[0,2\pi]$ of
$v=(\cos\theta,\sin\theta)$ for $d=2$, is uniformly randomly determined. The
particles that are not selected to change their velocities retain their
states, i.e., $v_{l}^{k}=v_{l}^{k-1}$.

  5. 5.

Return to 1.

###### Acknowledgements.

The author would like to acknowledge partial funding from the Japan-France
Integrated action Program PHC SAKURA, Grant No. JPJSBP120193219.

## References

  * (1) Barkai, N., Leibler, S.: Robustness in simple biochemical networks.  Nature 387(6636), 913–917 (1997)
  * (2) Berg, H., Brown, D.: Chemotaxis in escherichia coli analysed by three-dimensional tracking.  Nature 239(5374), 500–504 (1972)
  * (3) Block SM Segall JE, B.H.: Adaptation kinetics in bacterial chemotaxis.  journal of bacteriology 154(1), 312–323 (1983)
  * (4) Calvez, V.: Chemotactic waves of bacteria at mesoscale.  Journal of European Mathematical Society 22(2), 593–668 (2020)
  * (5) Calvez, V., Perthame, B., Yasuda, S.: Traveling wave and aggregation in a flux-limited keller-segel model.  Kinetic and Related Models 11(4), 891–909 (2018)
  * (6) Chalub, F., Markowich, P., Perthame, B., C., S.: Kinetic models for chemotaxis and their drift-diffusion limits.  Monatsh. Math. 142, 123–141 (2004)
  * (7) Dolak, Y., Schmeiser, C.: Kinetic models for chemotaxis: Hydrodynamic limits and spatio-temporal mechanisms.  J. Math. Biol. 51, 595–615 (2005)
  * (8) Emako, C., Gayrard, C., Buguin, A., Almeida, N., Vauchelet, N.: Traveling pulses for a two-species chemotaxis model.  PLoS Computational Biology 12(4) (2016)
  * (9) Erban, R., Othmer, H.: From individual to collective behaviour in bacterial chemotaxis.  SIAM J. Appl. Math. 65(2), 361–391 (2004)
  * (10) G., H.: Bacterial chemotaxis: the early years of molecular studies.  Annu. Rev. Microbiol 66, 285–303 (2012)
  * (11) Hillen, T., Othmer, H.: The diffusion limit of transport equations derived from velocity jump processes.  SIAM J. Appl. Math. 61(3), 751–775 (2000)
  * (12) Javen, G., Jashnsaz, H., Pressé, S.: From the volcano effect to banding: a minimal model for bacterial behavioral transitions near chemoattractant sources.  Physical Biology 15(4), 046002 (2018)
  * (13) Jiang, L., Ouyang, Q., Tu, Y.: Quantitative modeling of escherichia coli chemotactic motion in environments varying in space and time.  PLoS Computational Biology 6(4) (2010)
  * (14) Kalinin, Y.V., Jiang, L., Tu, Y., Wu, M.: Logarithmic sensing in escherichia coli bacterial chemotaxis.  Biophys J. 96(6), 2439–2448 (2009)
  * (15) Kanzler, L., Schmeiser, C., Tora, V.: Two kinetic models for non-instantaneous binary alignment collisions (2022).  DOI 10.48550/ARXIV.2203.15711.  URL https://arxiv.org/abs/2203.15711
  * (16) Keller, E., Segel, L.: Initiation of slime mold aggregation viewed as an instability.  J. Theor. Biol. 26, 399–415 (1970)
  * (17) Keller, E., Segel, L.: Model for chemotaxis.  J. Theor. Biol. 30, 225–234 (1971)
  * (18) Kollmann, M., Løvdok, L., Bartholomé, K., Timmer, J., Sourjik, V.: Design principles of a bacterial signalling network.  Nature 438(7067), 504–507 (2005)
  * (19) Li, S., Tan, P., Li, J., Tang, M., Hong, L.: Many-body effect renders universal subdiffusion to water on different proteins ((in print))
  * (20) Mittal, N., Budrene, E., Brenner, M., Van Oudenaarden, A.: Motility of escherichia coli cells in clusters formed by chemotactic aggregation.  Proceedings of the National Academy of Sciences of the United States of America 100(23), 13259–13263 (2003)
  * (21) Othmer, H., Dunbar, S., Alt, W.: Models of dispersal in biological systems.  Journal of Mathematical Biology 26(3), 263–298 (1988)
  * (22) Perthame, B. and Tang, M. and Vauchelet, N.: Derivation of the bacterial run-and-tumble kinetic equation from a model with biological pathway  Journal of Mathmatical Biology 73, 1161–1178 (2016)
  * (23) Perthame, B., Sun, W., Tang, M., Yasuda, S.: Multiple asymptotics of kinetic equations with internal states.  Mathematical Models and Methods in Applied Sciences 30(06), 1041–1073 (2020)
  * (24) Perthame, B., Yasuda, S.: Stiff-response-induced instability for chemotactic bacteria and flux-limited Keller-Segel equation.  Nonlinearity 31(9), 4065–4089 (2018)
  * (25) Saragosti, J., Calvez, V., Bournaveas, N., Perthame, B., Buguin, A., Silberzan, P.: Directional persistence of chemotactic bacteria in a traveling concentration wave.  Proceedings of the National Academy of Sciences 108(39), 16235–16240 (2011)
  * (26) Simons, J., Milewski, P.: The volcano effect in bacterial chemotaxis.  Mathematical and Computer Modelling 53(7-8), 1374–1388 (2011)
  * (27) Spiro, P., Parkinson, J., Othmer, H.: A model of excitation and adaptation in bacterial chemotaxis.  Proceedings of the National Academy of Sciences 94(14), 7263–7268 (1997)
  * (28) Tu, Y., Shimizu, T., Berg, H.: Modeling the chemotactic response of escherichia coli to time-varying stimuli.  Proceedings of the National Academy of Sciences 105(39), 14855–14860 (2008)
  * (29) Xue, C., Othmer, H.: Multiscale models of taxis-driven patterning in bacterial populations.  SIAM Journal on Applied Mathematics 70(1), 133–167 (2009)
  * (30) Xue, X., Xue, C., Tang, M.: The role of intracellular signaling in the stripe formation in engineered escherichia coli populations.  PLoS Computational Biology 14(6) (2018)
  * (31) Yasuda, S.: Monte carlo simulation for kinetic chemotaxis model: An application to the traveling population wave.  Journal of Computational Physics 330, 1022 – 1042 (2017)
  * (32) Yasuda, S.: Effects of internal dynamics on chemotactic aggregation of bacteria.  Physical Biology 18(6), 066001 (2021)