Traveling concentration pulses of bacteria in a generalized Keller-Segel model

TL;DR

This paper develops a generalized Keller-Segel model for bacteria that accounts for run-and-tumble dynamics, deriving new equations that predict traveling bacterial pulses and their maximum capacity, aligning well with experimental data.

Contribution

It introduces a polarization extended model incorporating angle bias and provides a bounded chemotactic drift, improving upon classical Keller-Segel equations.

Findings

Derived a Markovian response theory for run-and-tumble particles.

Identified a maximum bacterial capacity for traveling pulses.

Achieved good agreement with experimental observations.

Abstract

We formulate the Smoluchowski equation for a run-and-tumble particle. It includes the mean tumble rate in a chemical field, for which we derive a Markovian response theory. Using a multipole expansion and a reaction-diffusion equation for the chemoattractant field, we derive a polarization extended model, which also includes the recently discovered angle bias. In the adiabatic limit we recover generalized Keller-Segel equations with diffusion and chemotactic coefficients that depend on the microscopic swimming parameters. By requiring the tumble rate to be positive, our model possesses an upper bound of the chemotactic drift velocity, which is no longer singular as in the original Keller-Segel equations. Using the Keller-Segel model, we present an extensive study of traveling bacterial concentration pulses demonstrating how speed, width, and height of the pulse depend on the microscopic parameters. Most importantly, we discover a maximum number of bacteria that the pulse can sustain - the maximum carrying capacity. Finally, we obtain a remarkably good match to experimental results on the traveling bacterial pulse. It does not require a second, signaling chemical field nor a singular chemotactic drift velocity.

Tables

Table 1. Table 1: List of parameters used for the reference pulse in Fig. 2 and for the match to the experimental system shown in Fig.  6 .

Parameter	Value Fig. 2	Value Fig. 6	Ref.
$D$	$0.2\mu m^2{s}^{-1}$	same	Berg (2008)
$D_{\mathrm{rot}}$	$0.06s^{-1}$	same	Pohl et al. (2017)
$r$	$0$	$1.67\text{×}{10}^{-4}{s}^{-1}$	Saragosti et al. (2011)
$k$	$3.35\text{×}{10}^6{s}^{-1}$	same	Adler (1969)
${\lambda }_{\mathrm{equ}}$	$3s^{-1}$	same	Saragosti et al. (2011)
$v_0$	$25\mu {\mathrm{ms}}^{-1}$	same	Saragosti et al. (2011)
$⟨\cos \beta ⟩$	$0.392$	same	Saragosti et al. (2011)
$D_c$	$8\text{×}{10}^2\mu m^2{s}^{-1}$	same	Saragosti et al. (2011)
${\chi }_0$	$0.64{\lambda }_{\mathrm{equ}}\delta v_0^{}{}_{}{}^{-1}$	same	Saragosti et al. (2011)
$\delta$	$600\mu m$	same	Saragosti et al. (2011)
$c_0$	$1.26\text{×}{10}^6\mu m^{-3}$	$2.61\text{×}{10}^6\mu m^{-3}$	—
$c_t$	${10}^{-12}{c}_0$	${10}^{-1}{c}_0$	—
$A$	$5\text{×}{10}^4\mu m{}^2$	same	Saragosti et al. (2011)
$x_0$	$50\mu m$	same	Saragosti et al. (2011)
$N_0$	$1.5\text{×}{10}^5$	$1.5\text{×}{10}^5$	Saragosti et al. (2011)

Equations

$$\lambda (t)={\lambda }_{\mathrm{equ}}-\underset{-\infty }{\overset{t}{\int }}R(t-t^{\prime })c(\mathbf{r}(t^{\prime }))d{t}^{\prime }\text{leavevmode},$$

$$\lambda (\mathbf{r},\mathbf{e})={\lambda }_{\mathrm{equ}}-{\chi }_0{v}_0\mathbf{e}\cdot \frac{\text{gradient}c(\mathbf{r})}{c(\mathbf{r})}.$$

$$\lambda (\mathbf{r},\theta )={\lambda }_{\mathrm{equ}}+{\chi }_0{v}_0\frac{|\text{gradient}c(\mathbf{r})|}{c(\mathbf{r})}\cos \theta .$$

$$\lambda (\mathbf{r},\mathbf{e})={\lambda }_{\mathrm{equ}}-{\chi }_0\frac{v_0}{\delta }\tanh \left(\frac{c}{c_t}\right)\tanh \left(\delta \frac{|\text{gradient}c|}{c}\right)\mathbf{e}\cdot \hat{\mathbf{s}}.$$

$$\chi \left(\frac{|\text{gradient}c|}{c}\right):={\chi }_0\frac{v_0}{\delta }\tanh \left(\frac{c}{c_t}\right)\tanh \left(\delta \frac{|\text{gradient}c|}{c}\right).$$

$$\begin{array}{rl}\frac{\partial \psi }{\partial t}=& -\text{divergence}{\mathbf{J}}_{\mathrm{trans}}+-\mathcal{R}\cdot {\mathbf{J}}_{\mathrm{rot}}\\ & +\mathcal{F}\{\psi \}+\frac{r}{S_d}\int \psi (\mathbf{r},{\mathbf{e}}^{\prime },t)d{\mathbf{e}}^{\prime }.\end{array}$$

$$\begin{array}{rl}\mathcal{F}\{\psi \}=& -\lambda (\mathbf{r},\mathbf{e})\psi \\ & +\int P(\mathbf{e}-{\mathbf{e}}^{\prime },{\mathbf{e}}^{\prime })\lambda (\mathbf{r},{\mathbf{e}}^{\prime })\psi (\mathbf{r},{\mathbf{e}}^{\prime },t)d{\mathbf{e}}^{\prime }.\end{array}$$

$$\begin{array}{rl}\frac{\partial \psi }{\partial t}=& -\text{divergence}(v_0\mathbf{e}\psi )+D{\nabla }^2\psi +D_{\mathrm{rot}}{\mathcal{R}}^2\psi \\ & -\lambda (\mathbf{r},\mathbf{e})\psi +\int P(\mathbf{e}-{\mathbf{e}}^{\prime },{\mathbf{e}}^{\prime })\lambda (\mathbf{r},{\mathbf{e}}^{\prime })\psi (\mathbf{r},{\mathbf{e}}^{\prime },t)d{\mathbf{e}}^{\prime }\\ & +\frac{r}{S_d}\int \psi (\mathbf{r},{\mathbf{e}}^{\prime },t)d{\mathbf{e}}^{\prime },\end{array}$$

$$\begin{array}{r}\frac{\partial c}{\partial t}=D_c{\nabla }^{2}c-k\int \psi (\mathbf{r},{\mathbf{e}}^{\prime },t)d{\mathbf{e}}^{\prime }.\end{array}$$

$$\begin{array}{r}\int P(\mathbf{e}-{\mathbf{e}}^{\prime })d\mathbf{e}=1,\end{array}$$

$$\begin{array}{r}\int \mathbf{e}P(\mathbf{e}-{\mathbf{e}}^{\prime })d\mathbf{e}=⟨\cos \beta ⟩{\mathbf{e}}^{\prime },\end{array}$$

$$\frac{\partial \rho }{\partial t}=-\text{divergence}(v_0\mathbf{P})+D{\nabla }^2\rho +r\rho .$$

$$\begin{array}{rl}\frac{\partial \mathbf{P}}{\partial t}=& -v_0\nabla \mathbf{Q}-\frac{v_0}{d}\text{gradient}\rho +D{\nabla }^2\mathbf{P}\\ & -[D_{\mathrm{rot}}(d-1)+(1-⟨\cos \beta ⟩){\lambda }_{\mathrm{equ}}]\mathbf{P}\\ & +(1-⟨\cos \beta ⟩)\chi \left(\frac{|\text{gradient}c|}{c}\right)\mathbf{Q}\hat{\mathbf{s}}\\ & +\frac{1-⟨\cos \beta ⟩}{d}\chi \left(\frac{|\text{gradient}c|}{c}\right)\rho \hat{\mathbf{s}},\end{array}$$

$$\omega =D_{\mathrm{rot}}(d-1)+(1-⟨\cos \beta ⟩){\lambda }_{\mathrm{equ}},$$

$$\begin{array}{rl}\frac{\partial \mathbf{P}}{\partial t}=& -\omega \mathbf{P}+D{\nabla }^2\mathbf{P}-\frac{v_0}{d}\text{gradient}\rho \\ & +\frac{1-⟨\cos \beta ⟩}{d}\chi \left(\frac{|\text{gradient}c|}{c}\right)\rho \hat{\mathbf{s}}.\end{array}$$

$$\frac{\partial c}{\partial t}=D_c{\nabla }^{2}c-k\rho .$$

$$\mathbf{P}=-\frac{v_0}{\omega d}\text{gradient}\rho +\frac{1-⟨\cos \beta ⟩}{\omega d}\chi \left(\frac{|\text{gradient}c|}{c}\right)\rho \hat{\mathbf{s}}.$$

$$\begin{array}{rl}\frac{\partial \rho }{\partial t}=& D_{\mathrm{eff}}{\nabla }^2\rho +r\rho \\ & -\frac{v_0(1-⟨\cos \beta ⟩)}{\omega d}\text{divergence}\left[\chi \left(\frac{|\text{gradient}c|}{c}\right)\rho \hat{\mathbf{s}}\right],\end{array}$$

$$\begin{array}{r}\frac{\partial c}{\partial t}=D_c{\nabla }^{2}c-k\rho ,\end{array}$$

$${\mathbf{v}}_{\text{ch}}=\frac{v_0(1-⟨\cos \beta ⟩)}{\omega d}\chi \left(\frac{|\text{gradient}c|}{c}\right)\hat{\mathbf{s}}.$$

$$v_{\text{ch}}\le \frac{v_0}{d}.$$

$$\begin{array}{r}\int P(\mathbf{e}-{\mathbf{e}}^{\prime },{\mathbf{e}}^{\prime })d\mathbf{e}=1,\end{array}$$

$$\begin{array}{r}\int \mathbf{e}P(\mathbf{e}-{\mathbf{e}}^{\prime },{\mathbf{e}}^{\prime })d\mathbf{e}=⟨\cos \beta ⟩({\mathbf{e}}^{\prime }){\mathbf{e}}^{\prime }.\end{array}$$

$$\begin{array}{r}⟨\cos \beta ⟩({\mathbf{e}}^{\prime })=⟨\cos \beta {⟩}_0+\sigma \left(\frac{|\text{gradient}c|}{c}\right){\mathbf{e}}^{\prime }\cdot \hat{\mathbf{s}},\end{array}$$

$$\displaystyle\begin{aligned} \frac{\partial\mathbf{P}}{\partial t}=&-\bigg{\{}\omega+\left[\frac{1}{d+2}\chi\left(\frac{|\gradient{c}|}{c}\right)\sigma\left(\frac{|\gradient{c}|}{c}\right)\right]\\ &\qquad\times(\mathbf{1}+2\hat{\mathbf{s}}\otimes\hat{\mathbf{s}})\bigg{\}}\mathbf{P}+D\nabla^{2}\mathbf{P}-\frac{v_{0}}{d}\gradient{\rho}\\ &+\bigg{[}\frac{1-\langle\cos{\beta}\rangle_{0}}{d}\chi\left(\frac{|\gradient{c}|}{c}\right)+\frac{\lambda_{\mathrm{equ}}}{d}\sigma\left(\frac{|\gradient{c}|}{c}\right)\bigg{]}\rho\ \hat{\mathbf{s}}\,.\\ \end{aligned}$$

$$\begin{array}{rl}\frac{\partial \tilde{\rho }}{\partial \tilde{t}}=& {\tilde{\nabla }}^2\tilde{\rho }+\tilde{r}\tilde{\rho }-{\chi }_0(1-⟨\cos \beta ⟩)/\tilde{\delta }\\ & \times \widetilde{\text{divergence}}\left[\tanh \left(\frac{\tilde{c}}{{\tilde{c}}_t}\right)\tanh \left(\tilde{\delta }\frac{|\widetilde{\text{gradient}}\tilde{c}|}{\tilde{c}}\right)\tilde{\rho }\hat{\mathbf{s}}\right]\end{array}$$

$$\begin{array}{r}\frac{\partial \tilde{c}}{\partial \tilde{t}}=\frac{D_c}{D_{\text{eff}}}{\tilde{\nabla }}^2\tilde{c}-\frac{{\rho }_0}{c_0}\tilde{\rho }.\end{array}$$

$$v_p^{th}=\frac{N_{p}k}{A{c}_0}.$$

$$\lambda (t)={\lambda }_{\mathrm{equ}}-\underset{-\infty }{\overset{t}{\int }}R(t-t^{\prime })c(\mathbf{r}(t^{\prime }))d{t}^{\prime }\text{leavevmode}.$$

$$c({\mathbf{r}}^{\prime })=c({\mathbf{r}}_t)+\text{gradient}c\cdot ({\mathbf{r}}^{\prime }-{\mathbf{r}}_t)\text{leavevmode},$$

$$\begin{array}{rl}\lambda (t)=& {\lambda }_{\mathrm{equ}}-\text{gradient}c/c_a\cdot \underset{t-{\tau }_m}{\overset{t}{\int }}\tilde{R}(t-t^{\prime })\mathbf{r}(t^{\prime })d{t}^{\prime }\text{leavevmode},\end{array}$$

$$\mathbf{r}(t^{\prime })=\mathbf{r}(t)+v_0{\int }_t^{t^{\prime }}\mathbf{e}(t^{\prime \prime })d{t}^{\prime \prime }\text{leavevmode},$$

Full text

I Introduction

Collective motion of biological and artificial microswimmers is a most appealing research field as demonstrated in several review articles Ramaswamy (2010); Romanczuk et al. (2012); Marchetti et al. (2013); Zöttl and Stark (2016); Bechinger et al. (2016). The formation of various patterns and clustering have been investigated both experimentally and theoretically in systems of bacteria Adler (1966a); Budrene and Berg (1991); Mittal et al. (2003); Sokolov et al. (2007); Tailleur and Cates (2008); Zhang et al. (2010); Dombrowski et al. (2004); Reinken et al. (2018), of eukaryotic cells such as Dictyostelium discoideum or human sperm Rothschild (1948); Bonner and Savage (1947); Palsson and Othmer (2000); Ben-Jacob et al. (2000); Strassmann et al. (2000); Kaiser (2003); Riedel et al. (2005); Friedl and Gilmour (2009); Schoeller and Keaveny (2018), as well as in suspensions of active colloids Saintillan and Shelley (2008); Theurkauff et al. (2012); Buttinoni et al. (2013); Zöttl and Stark (2014); Pohl and Stark (2014); Speck et al. (2014); Liebchen et al. (2015); Pohl and Stark (2015); Zöttl and Stark (2016); Blaschke et al. (2016); Kuhr et al. (2017). In this article we study the collective behavior of a bacterial population, which in the concentration field of a chemoattractant forms a traveling solitary pulse.

The motility mechanism of the run-and-tumble bacterium E.coli has been extensively studied Berg and Brown (1972); Block et al. (1982); Berg (2008); Saragosti et al. (2011); Sourjik and Wingreen (2012); Masson et al. (2012); Pohl et al. (2017); Seyrich et al. (2018). Bacteria perform chemotaxis, the ability to sense and respond to chemical gradients in order to find better living conditions. They realize the chemotactic drift motion along a chemical gradient by elongated run phases if the environment becomes more favorable while runs are shortened in the opposite case. The internal chemotaxis machinery of the bacterium senses and compares the nutrient concentration in time, which is rationalized in a linear response theory for the tumble rate Macnab and Koshland (1972); Block et al. (1982); De Gennes (2004); Eisenbach (2004). More recently, a second chemotaxis strategy, called angle bias, has been reported Saragosti et al. (2011); Pohl et al. (2017); Seyrich et al. (2018). The mean reorientation angle during tumbling is reduced if the bacterium swims along a chemical gradient and increased in the opposite case. This also generates a net drift motion in the favorable direction. Finally, using logarithmic sensing, E.coli is able to perform chemotaxis in concentration fields varying by many orders of magnitude Dahlquist et al. (1972); Kalinin et al. (2009); Lazova et al. (2011). Such an ability is commonly described by Weber’s law in different physical areas Weber (1834); Fechner (1860).

A very interesting collective phenomenon in a bacterial population is a concentration pulse that travels along a capillary tube with almost no dispersion nearly like a soliton Adler (1966a); Saragosti et al. (2011), most recently also observed in a population with non-genetic variations Fu et al. (2018). The pulse is initiated in an initially uniform environment of a chemoattractant. A bacterial population concentrated in space eats the nutrient and thereby creates a chemical gradient along which it drifts towards untouched regions. Moreover, Adler in his experiments also observed that not all bacteria travel with the pulse but are left behind at the initial location Adler (1966a), which indicates a finite carrying capacity of the traveling pulse. Further chemoattractants present in Adler’s experiments then initiated further pulses emerging from the bacteria left behind.

A very prominent theoretical approach to describe the traveling bacterial pulse is the celebrated Keller-Segel model Keller and Segel (1971), originally introduced for the aggregation of slime molds Keller and Segel (1970). It couples a diffusion-drift equation for the bacterial density to a reaction equation for the nutrient. However, the Keller-Segel model has two drawbacks. First, a soliton solution (classified as unstable Nagai and Ikeda (1991)) only occurs if the chemotactic drift velocity diverges for vanishing nutrient concentration. Second, nutrient diffusion was neglected. Later, based on analytic arguments, Ref. Nagai and Ikeda (1991) demonstrated that traveling pulses also exist in the presence of nutrient diffusion. More importantly, Brenner et al. showed that the singularity in the chemotactic drift velocity is not necessary if one introduces a second chemoattractant, which the bacteria excrete themselves Brenner et al. (1998). Reference Saragosti et al. (2011) followed this approach to formulate a kinetic model (inspired by Ref. Othmer et al. (1988)), which describes traveling pulses in their experiments. Finally, a modification of this kinetic model has recently been used to investigate pulse propagation in the presence of two E.coli populations Emako et al. (2016). The Keller-Segel equations find wide applications in modeling bacterial chemotaxis as reviewed in Ref. Tindall et al. (2008). They have also been derived for active Brownian particles, which propel by self-diffusiophoresis, and for quorum-sensing run-and-tumble particles Rein et al. (2016).

Multipole expansions have frequently been applied to microswimmers in order to approximate the Smoluchowski equation for the full distribution function in the microswimmer’s position and orientation Golestanian (2012); Pohl and Stark (2014); Stark (2016); Rein et al. (2016); Reinken et al. (2018). Besides for density such expansions also provide an additional dynamic equation for the polarization, which unraveled interesting collective behavior of Janus particles Liebchen et al. (2015) and which also allowed to investigate steady-state distributions of run-and tumble particles Schnitzer (1993). Our derivation is inspired by the approach of the latter reference but extends it by introducing the concentration field of a chemoattractant.

In this article we formulate a Markovian response theory for the tumble rate. It includes logarithmic sensing for which we introduce an upper threshold. We use the tumble rate in the Smoluchowski equation and derive a polarization extended model (PE) to treat chemotaxis of non-interacting E.coli bacteria. The PE model contains equations for the bacterial density, the bacterial polarization, and the chemical concentration field. In a second step, we also include the recently discovered angle bias. In the adiabatic limit the PE model simplifies to a generalized Keller-Segel model (KS) where the coefficients for diffusion and chemotactic drift velocity depend on the microscopic swimming parameters of the bacterium. In particular, the chemotactic coefficient is not singular in the chemical concentration. We numerically solve both models for an initially uniform chemoattractant and a bacterial population concentrated in space using parameters that are realistic for the E.coli bacterium. The traveling bacterial pulse generated by both the PE and KS model are identical thus the KS model is a valid approximation of the full kinetic formalism. We present a detailed parameter study of the traveling pulse and identify a maximum carrying capacity as a consequence of the bounded chemotactic drift velocity, which has not been mentioned so far. It means that the pulse can only sustain a finite number of bacteria. Finally, we tune our parameters to match the experimental realization of the bacterial pulse in Ref. Saragosti et al. (2011). Hence, our generalized Keller-Segel model is able to describe traveling bacterial pulses without the need neither for a singular chemotactic drift velocity nor for a second chemoattractant.

The remainder of the article is organized as follows. We present the Markovian response theory for the tumble rate in Sect. II.1. We use it to derive the polarization extended model (PE) and the generalized Keller-Segel model (KS) in Sects. II.2 and II.3. We also incorporate the angle bias and formulate a non-dimensional version of the KS model in Sects. II.4 and II.5. Details of the numerical solution scheme are given in Sect. III and Sect. IV presents our detailed numerical study. We close with conclusions and an outlook in Sect. V.

II Model

II.1 Markovian response theory for tumble rate

Bacteria tumble less when moving up a chemical gradient. Based on the established linear-response theory, we formulate an equation for the tumble rate $\lambda(\mathbf{r},\mathbf{e})$ as a function of the swimming direction $\mathbf{e}$. Below, we will relate it to the angle $\theta$ relative to the local gradient $\gradient c$ of a chemoattractant with density $c$.

We start with the linear-response theory. It gives the tumble rate $\lambda(t)$ as a function of time and depends on the bacterium’s past trajectory $\mathbf{r}(t^{\prime})$,

[TABLE]

where we have introduced the response kernel $R(t)$ and $\lambda_{\mathrm{equ}}$ is the tumble rate without any chemical gradient. Note that Eq. (1) describes a non-Markovian process. We convert it to a Markovian process with $\lambda$ depending on location $\mathbf{r}$ and swimming direction $\mathbf{e}$ by averaging over all possible histories, which gives the mean history for a given swimming direction $\mathbf{e}$. To do so, we split the integral on the R.H.S. of Eq. (1) into contributions from individual runs, during which the according swimming directions $\mathbf{e}_{i}$ are assumed to be constant. Averaging over the history of all possible paths, we can show that each of these contributions gives a term proportional to the scalar product $\mathbf{e}\cdot\frac{\gradient{c}}{c}$. We also used here that the response kernel is proportional to the inverse background conversation, $R\propto 1/c$, which was indeed measured in experiments Masson et al. (2012). For the detailed derivation we refer to appendix A and present the final result,

[TABLE]

Here, $v_{0}$ is the swimming velocity of the bacterium and $\chi_{0}$ is a unitless measure of the chemotactic strength. It depends on integrals over the response function $R$ and moments of the tumble angle distribution $P(\beta)$. Note that we obtain here $\lambda\propto\gradient(\ln c)$ commonly known as logarithmic sensing and Weber’s law. It was measured, for example, in Ref. Lazova et al. (2011).

When we define the orientation angle relative to the negative chemical gradient, $\cos\theta=-\mathbf{e}\cdot\frac{\gradient{c}}{|\gradient{c}|}$, the tumble rate becomes

[TABLE]

In Fig. 1 adapted from Ref. Pohl et al. (2017), the red points show experimental data for the mean tumble rate $\lambda(\mathbf{r},\theta)$. It was obtained by averaging over a population of around $1000$ individual bacteria in a linear gradient. The appropriate cosine fit (red line) confirms our theoretically derived result of Eq. (2).

In general, Eq. (2) can produce a negative tumble rate for a sufficiently large gradient of $\log c$, which is even singular at $c=0$. The reason is the relation $R\propto 1/c$ for the response function mentioned earlier. While this dependence was measured for a wide range of background concentrations Masson et al. (2012), clearly the deviation from $\lambda_{\mathrm{equ}}$ in Eq. (2) has to saturate to a value smaller than $\lambda_{\mathrm{equ}}$. Furthermore, it is known that the bacterium needs a small threshold concentration $c_{t}$ to perform chemotaxis Adler (1969). To implement both constraints we use the hyperbolic tangent function in Eq. (2) and write,

[TABLE]

Here, we have introduced $\hat{\mathbf{s}}=\frac{\gradient{c}}{|\gradient{c}|}$ and do not explicitly state the space dependence of the concentration $c$. This expression recovers Eq. (2) for $\delta|\gradient{c}|/c<1$ and $c>c_{t}$, while it smoothly approaches the minimum tumble rate $\lambda_{\mathrm{equ}}-\chi_{0}\tfrac{v_{0}}{\delta}$ at a large relative chemical gradient or it becomes $\lambda_{\mathrm{equ}}$ for $c<c_{t}$. The chemotactic length $\delta$ quantifies the strength of the logarithmic derivative of $c(\mathbf{r})$. In the following, we use the notation

[TABLE]

II.2 Polarization extended model (PE)

Smoluchowski equation. – We first construct dynamic equations for the evolution of the one-particle distribution function $\psi(\mathbf{r},\mathbf{e},t)$ of position $\mathbf{r}$ and orientation $\mathbf{e}$ at time $t$ and the concentration of chemoattractant, $c(\mathbf{r},t)$. We begin with a generalized Smoluchowski equation for $\psi$, which contains the usual contributions from translational and rotational currents, $\mathbf{J}_{\mathrm{trans}}$ and $\mathbf{J}_{\mathrm{rot}}$, but also contributions from tumble events represented by $\mathcal{F}\{\psi\}$ and from cell division and death, $\frac{r}{S_{d}}\int\psi(\mathbf{r},\mathbf{e}^{\prime},t)d\mathbf{e}^{\prime}$:

[TABLE]

Here, $\mathcal{R}=\mathbf{e}\times\partial_{\mathbf{e}}$ and $r=\ln{(2)}/\tau$ is the net creation rate with $\tau$ being the mean doubling time of bacterial cells Saragosti et al. (2011). Here, we assume that the net creation of cells does not depend of their directon $\mathbf{e}$ and $S_{d}$ is the surface area of a $d$ dimensional unit sphere (full solid angle). For the translational current we include active motion and translational diffusion, $\mathbf{J}_{\mathrm{trans}}=v_{0}\mathbf{e}\psi-D\gradient{\psi}$, where $D$ is the translational diffusion coefficient and $v_{0}$ is the bacterial swimming speed. The rotational current is purely diffusive, $\mathbf{J}_{\mathrm{rot}}=-D_{\mathrm{rot}}\mathcal{R}\psi$, where $D_{\mathrm{rot}}$ is the rotational diffusion coefficient. According to Ref. Berg and Brown (1972) we take a Poisson distribution for the run times and so write the term for the tumble events as

[TABLE]

We introduced the tumble rate $\lambda(\mathbf{r},\mathbf{e})$ and $P(\mathbf{e}-\mathbf{e}^{\prime},\mathbf{e}^{\prime})$ is the probability of a bacterium to reorient from orientation $\mathbf{e}^{\prime}$ to $\mathbf{e}$. In Eq. (7) the first term on the R.H.S. represents events, which cause bacteria with orientation $\mathbf{e}$ to tumble into any orientation, and the second term represents all events, which cause bacteria with other orientations to tumble into orientation $\mathbf{e}$. The complete Smoluchowski equation for the evolution of $\psi$ now reads

[TABLE]

For completeness, we also write a reaction-diffusion equation for the chemoattractant concentration $c$, which is also consumed by bacteria with constant rate $k$,

[TABLE]

Multipole expansion. – In order to make progress with Eqs. (8) and (9), we assume that the probability distribution for a specific tumble event does not depend on the initial orientation of the bacterium, $P(\mathbf{e}-\mathbf{e}^{\prime},\mathbf{e}^{\prime})=P(\mathbf{e}-\mathbf{e}^{\prime})$. Therefore, we can write for the zeroth and first moment,

[TABLE]

where $\langle\cos{\beta}\rangle$ is the mean of the cosine of the reorientation angle $\beta$ Schnitzer (1993). We take the tumble rate $\lambda(\mathbf{r},\mathbf{e})$ to vary as in Eq. (4). Now, we integrate Eq. (8) over all orientations $\mathbf{e}$ and define the bacterial density $\rho(\mathbf{r},t)=\int\psi(\mathbf{r},\mathbf{e},t)d\mathbf{e}$ and polarization $\mathbf{P}(\mathbf{r},t)=\int\mathbf{e}\psi(\mathbf{r},\mathbf{e},t)d\mathbf{e}$, which correspond to the zeroth and first moments of $\psi$ respectively, to obtain

[TABLE]

We have also used $\int\mathcal{R}^{2}\psi d\mathbf{e}=0$ and the normalization condition Eq. (10) to show that tumbling does not contribute to Eq. (12), as it should be.

In order to derive a dynamic equation for the polarization, we compute $\int\mathbf{e}\,$Eq. (8)$\,d\mathbf{e}$ and introduce the quadrupole moment ${\mathit{\mathbf{Q}}}=\int(\mathbf{e}\otimes\mathbf{e}-\frac{\mathbf{1}}{d})\psi(\mathbf{r},\mathbf{e},t)d\mathbf{e}$, with $d$ being the number of spatial dimensions. This gives

[TABLE]

where we used $\mathcal{R}^{2}\mathbf{e}=-(d-1)\mathbf{e}$ and Eq. (11). To truncate the multipole expansion, we neglect the quadrupole moment ${\mathit{\mathbf{Q}}}$ and define the relaxation rate

[TABLE]

on which polar order relaxes or decorrelates in time. Thus, we ultimately obtain

[TABLE]

Finally, with our definition of bacterial density $\rho$, we can write Eq. (9) in a simpler form,

[TABLE]

II.3 The Keller-Segel model as adiabatic limit

In the case of high Peclet numbers ($\mathrm{Pe}=av_{0}/D\gg 1$), where we can neglect translational diffusion, and on large times $t\gg\frac{1}{\omega}$, where the adiabatic limit $\tfrac{\partial\mathbf{P}}{\partial t}\approx 0$ applies, the polarization from Eq. (15) becomes

[TABLE]

We remind that $\hat{\mathbf{s}}=\gradient{c}/|\gradient{c}|$. Substituting Eq. (17) into Eq. (12), we obtain the generalized Keller-Segel model

[TABLE]

where $D_{\mathrm{eff}}=v_{0}^{2}/\omega d+D$ is the typical translational diffusion coefficient of an active particle, the orientation of which decorrelates on the characteristic time $\omega^{-1}$.

From the third term on the R.H.S. of Eq. (18) we read off the chemotactic velocity along the chemical gradient,

[TABLE]

Taking $\chi=\chi_{0}v_{0}|\gradient{c}|/c$, we recover the model suggested by Keller and Segel with $\mathbf{v}_{\text{ch}}\propto\gradient{c}/c$ Keller and Segel (1971), and the chemotactic drift velocity is determined by a combination of microscopic parameters, $\tfrac{\chi_{0}v_{0}^{2}(1-\langle\cos\beta\rangle)}{\omega d}$, called the chemotactic constant. However, as stated earlier, according to Eq. (2) the form for $\chi$ implies negative tumble rates for small $c$ and sufficiently steep chemoattractant gradient. The maximum value that $\chi$ can physically assume is $\lambda_{\mathrm{equ}}$, where the tumble rate becomes negative. As a result, the chemotactic speed $v_{\text{ch}}=|\mathbf{v_{\textrm{ch}}}|$ is also bounded. Taking $\chi=\lambda_{\mathrm{equ}}$ and approximating $\omega\approx(1-\langle\cos{\beta}\rangle)\lambda_{\mathrm{equ}}$ since $D_{\mathrm{rot}}$ in Eq. (14) is usually much smaller than $\lambda_{\mathrm{equ}}$, we find

[TABLE]

This shows that an appropriately bounded tumble rate is closely linked to a physically bounded chemotactic drift.

II.4 Bias of tumble angles

Up to this point we have ignored the effect of a bias in the tumble angle towards smaller mean values when swimming up the chemical gradient. This has recently been observed in experiments Saragosti et al. (2011); Pohl et al. (2017). We now introduce it by allowing the probability distribution for a specific tumble event to explicitly depend on the initial orientation of the bacterium, $\mathbf{e^{\prime}}$. Hence Eqs. (10) and (11) become

[TABLE]

Equation (22) states, whatever the initial orientation of the bacterium the respective distribution is normalized. In Eq. (23) the value of the mean cosine of the tumble angle $\langle\cos{\beta}\rangle(\mathbf{e}^{\prime})$ now explicitly depends on the initial orientation $\mathbf{e}^{\prime}$ before the tumble event.

The effect of an angle bias is to lower the mean tumble angle when the bacterium aligns with the chemoattractant gradient, hence the value of $\langle\cos{\beta}\rangle(\mathbf{e}^{\prime})$ will increase for stronger alignment. Expanding $\langle\cos{\beta}\rangle(\mathbf{e}^{\prime})$ up to the first Legendre polynomial, thus taking into account the leading polar correction, yields

[TABLE]

where $\sigma$ is a positive and monotonically increasing function. It is bounded such that its maximum value $\sigma_{\mathrm{max}}\leq 1-\langle\cos{\beta}\rangle_{0}$, with $\langle\cos{\beta}\rangle_{0}$ being the mean cosine of the tumble angle, when the angle bias is not taken into account. Using Eqs. (2), (22), (23), and (24), we can retrace the steps of the multipole expansion (see appendix B for details) to obtain an extended form for Eq. (15) with Eqs. (12) and (16) remaining unchanged,

[TABLE]

One immediately recognizes that the angle bias renormalizes the relaxation rate of the polarization and makes it anisotropic. Thus polarizations along and perpendicular to the chemical gradient relax with different rates. In the adiabatic limit $\frac{\partial\mathbf{P}}{\partial t}\approx 0$ and for large $\mathrm{Pe}$ we can again solve for the polarization by inverting the matrix in front of $\mathbf{P}$ in Eq. (25). Substituting the resulting equation into Eq. (12), we again obtain a generalized Keller-Segel equation and a chemotactic velocity $\mathbf{v}_{\text{ch}}$ along the chemical gradient. It now also depends on the angle bias quantified by $\sigma$. We refrain from giving the lengthy expression.

II.5 Rescaling the Keller-Segel equations

In order to identify essential parameters especially in the generalized Keller-Segel equations (18) and (19), we introduce unitless quantities. First, we rescale time with the chemical consumption rate, $\tilde{t}=kt$, lengths by the distance $l=(\tfrac{v_{0}^{2}}{\omega dk})^{1/2}$, by which a bacterium diffuses in time $k^{-1}$, $\tilde{\mathbf{r}}=\mathbf{r}/l$, and the net creation by $k$, $\tilde{r}=r/k$. Second, we refer the bacterial and chemical densities to their initial values, $\tilde{\rho}=\rho/\rho_{0}$ and $\tilde{c}=c/c_{0}$, respectively. Finally, we introduce the rescaled chemotactic length $\tilde{\delta}=\delta/l$ and the rescaled threshold density $\tilde{c}_{t}=c_{t}/c_{0}$. This allows us to write the generalized Keller-Segel equations (18) and (19), where chemotactic response is bounded by the hyperbolic tangents, in rescaled form:

[TABLE]

To arrive at this rescaling, we used $D_{\mathrm{eff}}\approx v_{0}^{2}/\omega d$, where we neglected the thermal contribution to $D_{\mathrm{eff}}$. The rescaling shows that the generalized Keller-Segel equations are described by a set of six relevant parameters: $\{\tfrac{D_{c}}{D_{\mathrm{eff}}},\tilde{r},\chi_{0}(1-\langle\cos\beta\rangle),\tfrac{\rho_{0}}{c_{0}},\tilde{\delta},\tilde{c}_{t}\}$.

III Details of Numerical Solution Scheme

In the following we study in detail traveling bacterial pulses in an initally uniform density field of a chemoattractant by numerically solving both the polarization extended model (PE) of Eqs. (12), (15) and (16) and the generalized Keller-Segel model (KS) of Eqs. (18) and (19). The experiments in Ref. Saragosti et al. (2011) are performed in microchannels with cross section $A$. Neglecting any effects at the channel walls, we take all densities to just depend on the $x$ coordinate along the channel. To solve the respective system of equations, we apply a predictor-corrector method at any given time step to efficiently propagate the field variables in time Press et al. (2007). As initial field values we choose an exponentially distributed bacterial density, $\rho(x,t=0)=\rho_{0}\exp(-x/x_{0})$, a uniform density of the chemoattractant, $c(x,t=0)=c_{0}$, and zero polarization $P_{x}(x,t=0)=0$. During time integration no-flux boundary conditions are employed at $x=0$ and a sufficiently large $x_{\infty}$, $\partial_{x}\rho\rvert_{0,x_{\infty}}=0$ and $\partial_{x}c\rvert_{0,x_{\infty}}=0$, while polarization stays zero, $P_{x}\rvert_{0,x_{\infty}}=0$. Finally, when integrating Eq. (16) the sink term can produce negative concentrations of the chemoattractant Rosen (1974); Rosen and Baloga (1975). To avoid this, we set the concentration $c$ to zero whenever it would become negative. This allows the bacteria to fully degrade the chemoattractant without producing negative values for $c$.

When we solve our equations with real parameters, we rely on Ref. Saragosti et al. (2011) and take the channel length $x_{\infty}=10^{5}\mathrm{\mu m}$ and the initial decay length of the bacterial density as $x_{0}=$50\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}. This ensures that at $t=0$ $99\%$ of the bacteria can be found within $200\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ at the channel end at $x_{0}=0$. We also assume a channel cross section $A=$500\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$\times$100\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m} to calculate the initial number of bacteria $N_{0}=\rho_{0}x_{0}A$, which we use as a parameter instead of $\rho_{0}$ in the following. We divide the channel length into $5\text{\times}{10}^{4}\text{\,}$ grid points so that the grid length is $2\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ and use the time step $0.01\text{\,}\mathrm{s}$ for integrating our equations in time. All the relevant parameters are given in Table 1. Finally, we will also numerically solve the rescaled Keller-Segel equations (26) and (27) in order to explore the dependence on some of the relevant dimensionless parameters.

IV Traveling concentration pulses of bacteria

We first introduce the traveling bacterial pulse for a reference system using two values for the initial number of bacteria, then perform a systematic parameter study, and finally demonstrate a perfect match with the experimentally observed bacterial pulse reported in Ref. Saragosti et al. (2011).

IV.1 Reference system

Figure 2(a) shows a series of snapshots of the bacterial density profile $\rho(x,t)$, the concentration field $c(x,t)$, and the polarization $P_{x}(x,t)$ at equally spaced times for realistic parameters of the E.coli bacterium listed in Table 1. Solid lines represent numerical solutions of the polarization-extended model (PE) and dashed lines the adiabatic approximation leading to the generalized Keller-Segel model (KS). Video S1 of the supplemental material shows an animation of the propagating profiles.

Clearly, while eating up the chemoattractant completely, a traveling pulse in the bacterial density forms that propagates with constant pulse speed $v_{p}=$4.68\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}\mathrm{s}^{-1}$$. It has a comparable width to the traveling step in the chemoattractant profile. In contrast to the bacterial solitons derived from the original KS model in Ref. Keller and Segel (1971), our bacterial pulse shows a small dispersion visible from the slight decrease of the pulse height. It is caused by bacteria that cannot follow the pulse at small chemoattractant concentrations since in our model the chemotactic drift velocity $\mathbf{v}_{\text{ch}}$ of eq. (20) has an upper bound. In contrast, in the original KS model $\mathbf{v}_{\text{ch}}$ diverges at small chemoattractant concentrations Keller and Segel (1971), which allows all bacteria to stay in the traveling pulse. Thus, we demonstrate when one allows dispersion a singular chemotactic drift velocity is no longer necessary for traveling pulses to occur. Note, already Ref. Scribner et al. (1974) used the KS model with a bounded chemotactic drift and observed traveling pulses which was ignored in more recent works Brenner et al. (1998); Saragosti et al. (2011). Instead, a second chemoattractant was proposed as explained in the introduction.

One realizes that the profiles generated from the KS and PE model are identical except in the beginning. We start with $P_{x}=0$ in the PE model whereas in the KS model a non-zero polarization is calculated from eq. (17). It is due to the initial gradient in the bacterial density. Thus, we conclude that the adiabatic limit $\frac{\partial P}{\partial t}\approx 0$ as a condition for deriving the generalized KS model is fulfilled. Indeed, the decorrelation or decay time $\omega^{-1}=$0.49\text{,}\mathrm{s}$$ is much smaller than the characteristic time for the pulse propagation. For the latter, we approximately find $200\text{\,}\mathrm{s}$ for the pulse to travel its own width and the polarization always assumes its stationary value. Our finding also means that the kinetic models of Refs. Saragosti et al. (2011, 2012); Emako et al. (2016), which work with the full one-particle distribution function $\psi(\mathbf{r},\mathbf{e},t)$, are not necessary to describe pulse propagation. They can be reduced to the Keller-Segel equations.

In Fig. 2(a) not all bacteria travel with the pulse but some remain at the initial location. This also occurs in the experiments of Ref. Adler (1966a). However, the remaining bacteria perform chemotaxis in oxygen as a second chemoattractant and thereby initiate a secondary pulse. Since we do not incorporate another chemoattractant, the bacterial distribution at the initial location only broadens by diffusion. Finally, we mention previous numerical work on the KS model that also showed the bacteria left behind Scribner et al. (1974).

In their original work Keller and Segel derived an analytic expression for the speed of the bacterial soliton as a function of the number of bacteria in the soliton $N_{p}$, the consumption rate $k$, the cross section $A$, and the initial chemoattractant concentration $c_{0}$ Keller and Segel (1971),

[TABLE]

By integrating the bacterial density along the $x$ direction at time $t=$2000\text{,}\mathrm{s}, we obtain $N_{p}=A\int\rho dx=$0.88\text{\times}{10}^{5}\text{\,} bacteria in the pulse and for the number of bacteria left behind close to $x=0$, $N_{c}=$0.57\text{\times}{10}^{5}\text{,}. Note, $N_{p}$ and $N_{c}$ do not add up to $N_{0}$ since there are also bacteria in the trail between the initial location and the pulse. Using $N_{p}$ and the parameter values of the reference system from Table [1](#S3.T1) in expression ([28](#S4.E28)), we obtain $v_{p}^{th}=$4.69\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}\mathrm{s}^{-1}, which is in very good agreement with our numerical value of $v_{p}=$4.68\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}\mathrm{s}^{-1}$$.

In Fig. 2(b) we lower the number of bacteria $N_{0}$ by a factor of three. Now, all bacteria travel in the pulse and none are left behind. This suggests that the traveling pulse can only carry a certain number of bacteria and thus has a maximum carrying capacity. We will investigate it in more detail in the following parameter studies, where we use the two traveling pulses from Fig. 2 as a reference. For the traveling pulse in Fig. 2(b) we determine a smaller pulse speed of $v_{p}=$2.66\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}\mathrm{s}^{-1}. It matches very well with the theoretical prediction from Eq. ([28](#S4.E28)) using $N_{p}=N_{0}=$0.5\text{\times}{10}^{5}\text{\,}. Video S2 of the supplemental material shows an animation of the traveling profiles.

IV.2 Parameter studies

IV.2.1 Influence of bounded chemotactic drift

To keep the tumble rate positive, we introduced the chemotactic length $\delta$ in eq. (4), which prevents the chemotactic drift velocity in eq. (20) to become arbitrarily large. Furthermore, for the chemotactic response the lower chemical threshold $c_{t}$ was introduced. In Fig. 3 we explore the influence of both parameters on the traveling bacterial pulse. The chemotactic length $\delta$ increases from left to right and the threshold concentration $c_{t}$ from top to bottom. The reference system of Fig. 2(a) is in the center.

A smaller chemotactic length $\delta$ means that the bacterium can sense larger chemical gradients and that the drift velocity $\mathbf{v}_{\text{ch}}$ saturates at a larger value proportional to $\delta^{-1}$. However, it cannot become arbitrarily small since then the tumble rate in eq. (4) becomes negative. In this case our numerical solution scheme is unstable and the bacterial density becomes negative.

The length $\delta=$384\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ is close to the minimal value. In particular for small $c_{t}$ nearly all bacteria travel in the pulse, none are left at the origin. As a result, the pulse travels the fastest. Increasing $\delta$ to $600\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, bacteria left behind are clearly visible. The pulse contains less bacteria and, therefore, is slower. This is also in agreement with the smaller chemotactic drift velocity. Interestingly, for the smallest $c_{t}$ we observe a second propagating pulse strongly decreasing in height. Finally, if we increase $\delta$ by a factor of $10$ to $6000\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, the majority of bacteria stay close to the initial location while only a smaller number travels in the pulse (note the 10 times smaller range of the vertical axis). Thus, the pulse speed is small and the pulse has not yet separated from the non-propagating bacteria. In conclusion, increasing $\delta$ decreases the maximum carrying capacity significantly and makes the pulse slower.

Increasing the threshold concentration $c_{t}$ from nearly zero to $c_{t}=0.1c_{0}$ has three effects. First, the dispersion of the pulse increases (especially strongly from the second to the third row) which slows down the pulse. Second, the shape of the pulse becomes more asymmetric as bacteria at the rear flank cannot follow the pulse. Third, the number of bacteria left behind at the initial location increases slightly. In the upper middle plot we recognize that the threshold $c_{t}$ is so low that the remaining bacteria can still travel by chemotaxis, although with a stronger dispersion as the first pulse. Finally, even for the vanishing threshold of the upper left plot a slight dispersion is visible. This again suggests that a true propagating soliton, for which the pulse shape does not vary in time, is not possible as long as the chemotactic drift velocity is bounded.

IV.2.2 Quantitative study of the rescaled Keller-Segel equations

We now consider the rescaled Keller-Segel equations and study the propagating bacterial pulse in detail by plotting the rescaled pulse speed $v_{p}$, pulse full width at half maximum $\Delta x_{h}$, and pulse amplitude $\rho_{\textrm{max}}/\rho_{0}$ as a function of the remaining parameters $\rho_{0}/c_{0}$, $D_{c}/D_{\textrm{eff}}$ and $\chi_{0}(1-\langle\cos\beta\rangle)$. Again, we neglect bacterial growth by setting $\tilde{r}=0$. Figure 4 shows all results and the relevant parameters are given in the figure caption.

In Fig. 4(a) we see that the pulse speed depends linearly on $\rho_{0}/c_{0}$ in agreement with the Keller-Segel prediction of Eq. (28) but then saturates at a constant value. The reference pulse from Fig. 2(b) (purple disc), where no bacteria are left behind at the origin, is located in the linear regime, while the reference pulse from Fig. 2(a) (orange star), where some bacteria remain close to the origin, propagates in the saturated regime. Thus, in the first case adding more bacteria to the system increases the number of bacteria in the traveling pulse and speeds it up. In contrast, in the second case additional bacteria remain close to the initial location. Thus, the traveling concentration pulse has a maximum carrying capacity $N^{\ast}$ with respect to the amount of bacteria it can hold and all other bacteria are left behind. The transition between both regimes occurs at the critical ratio $(\rho_{0}/c_{0})^{*}$. This allows to discuss the following scenario, for example. Lowering $c_{0}$ at constant $\rho_{0}$ speeds up the pulse in the linear regime since bacteria degregate the chemoattractant faster. However, once $(\rho_{0}/c_{0})^{*}$ is reached the pulse looses bacteria to keep the pulse velocity constant. Thus the carrying capacity of the pulse decreases for lower $c_{0}$.

The ratio $\rho_{0}/c_{0}$ also influences the pulse shape. In the linear regime with increasing $\rho_{0}/c_{0}$ the pulse becomes narrower while its absolute height $\rho_{\mathrm{max}}$ roughly increases with $\rho_{0}^{2}$. When reaching the saturated regime, the pulse width stays constant as should $\rho_{\mathrm{max}}$. Thus for the relative height we find $\rho_{\mathrm{max}}/\rho_{0}\propto(\rho_{0}/c_{0})^{-1}$.

In Fig. 4(b) we show the pulse speed does not significantly depend on the ratio of diffusion constants, $D_{c}/D_{\textrm{eff}}$, for both study cases a (blue stars) and b (green circles). This is in contrast to Ref. Holz and Chen (1979) where the authors proposed a correction term to Eq. (28), which predicts a decrease of the pulse speed with increasing $D_{c}$. However, when examining the bacterial pulse profile, we observe that for larger $D_{c}/D_{\textrm{eff}}$ the pulse needs longer to form. It needs longer to eat up all the chemoattractant at the origin due to the larger diffusive flux of chemoattractant into the depleted areas. But once the bacteria have fully degraded the chemoattractant, the pulse propagates with the same speed $v_{p}$ independent of $D_{c}$. For increasing $D_{c}/D_{\textrm{eff}}$ the width of the pulse also increases while the amplitude decreases. However, both trends are very small since the quantities do not even change by a factor of two while varying $D_{c}/D_{\textrm{eff}}$ over four orders of magnitude. Finally, noting the relevant length scale $l=\sqrt{D_{\textrm{eff}}/k}$ to depend on the effective bacterial diffusion constant $D_{\textrm{eff}}$, we find $v_{p}\propto\sqrt{D_{\textrm{eff}}}$. Moreover, the pulse width increases significantly with $D_{\textrm{eff}}$ while the absolute height increases only slightly.

Figure 4(c) shows the results for varying the chemotactic parameter $\chi_{0}(1-\langle\cos\beta\rangle)$. For values larger than $50$ the numerical scheme becomes unstable akin to the instability in the chemotactic length $\delta$ already discussed in Sec. IV.2.1. The pulse speed increases linearly in the chemotactic parameter for the study case a (blue stars) and also for the study case b (green circles) in the range $\chi_{0}(1-\langle\cos\beta\rangle)<20$. To understand this finding, we looked in detail at the bacterial profiles. In study case a we find that with increasing $\chi_{0}(1-\langle\cos\beta\rangle)$ more and more bacteria from the vicinity of the initial location enter the pulse, which according to Eq.  (28) then speeds up. Thus we conclude for the maximum carrying capacity of the pulse, $N^{*}\propto\chi_{0}(1-\langle\cos\beta\rangle$. A similar observation in connection with the scenario of Fig. 4(b) gives $N^{*}\propto\sqrt{D_{\textrm{eff}}}$. Note that our results are in contrast to Ref. Rosen (1975) which found $v_{p}\propto\sqrt{\chi_{0}}$. In study case b (green circles) we start with a smaller number of bacteria. Thus, at $\chi_{0}(1-\langle\cos\beta\rangle)\approx 20$ all bacteria have entered the pulse, which then travels with constant speed for further increase in $\chi_{0}(1-\langle\cos\beta\rangle)$. With increasing chemotactic parameter the width of the traveling pulse decreases in the study case a (blue stars). The curve of study case b (green circles) follows this trend until the pulse speed becomes constant and then steadily increases.

For the pulse amplitude of the two study cases the behavior is inverted. The curves are well fitted by $y(x)=Ax^{B}$, where the constants $A$ differ by approximately a factor of three, which is the factor by which the density ratios $\rho_{0}/c_{0}$ of both cases differ. It becomes visible since we plot the reduced amplitude $\rho_{\mathrm{max}}/\rho_{0}$. The exponents are nearly the same. The amplitude of the study case b (green circles) decreases for $\chi_{0}(1-\langle\cos\beta\rangle>20$ and thereby compensates the increasing pulse width as the number of bacteria in the pulse is constant.

IV.2.3 Influence of growth rate $r$

Finally, we investigate the influence of the growth term in the Keller-Segel equation (18). Figure 5 shows propagating pulses for three different growth rates $r$, while the other parameters are chosen as in the reference system but with the reduced initial number of bacteria $N_{0}=$0.5\text{\times}{10}^{5}\text{,}$$. It is below the maximum carrying capacity of the pulse and was used in Fig. 2(b).

Consequently, in the upper panel the pulse grows due to the non-zero growth rate and speeds up in time until the maximum carrying capacity is reached at around $2600\text{\,}\mathrm{s}$. Then, the pulse propagates with constant shape like a perfect soliton. However, in our case the pulse leaves a trail of bacteria behind, which originate from the continuous bacterial growth.

In the middle and lower panel, the maximum carrying capacity is reached after $1000\text{\,}\mathrm{s}$ and $50\text{\,}\mathrm{s}$, respectively. Interestingly, the pulse does no longer separate from the broad distribution of bacteria originating from the initial location but rather sits on top of its right flank. In the lower panel the pulse is fastest and its amplitude is highest. This comes from the fact that the broad distribution of the remaining bacteria is much larger compared to the middle panel and more bacteria actively take part in the degregation process of the chemoattractant. As a consequence, the pulse propagates faster.

Last, we observe that with increasing growth rate the pulse becomes more peaked. This is reminiscent to Fig. 4(a) and (c), where a faster pulse has a smaller pulse width.

IV.3 Matching the experimental pulse

Figure 6 (upper panel) shows the traveling bacterial concentration pulse recorded in the experiments of Ref. Saragosti et al. (2011) and compares it to the numerical solution of our polarization extended model. Both propagating pulses agree very well in shape and in speed $v_{p}=$3.8\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}\mathrm{s}^{-1}$$. We extracted the experimental data from Fig. 1(b) in Ref. Saragosti et al. (2011) and in our model mainly used parameters from the same publication including a non-zero growth rate but also added missing values from Refs. Pohl et al. (2017) and Adler (1969).

Moreover, a realistic value for the sensing threshold $c_{t}=$2.61\text{\times}{10}^{5}\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{-3}$=0.1c_{0}$ was chosen Adler (1969). This was necessary to match the asymmetry and dispersion of the pulse. The full parameter set is given in Table 1.

V Discussion and Conclusions

In the first part of this article we derived the celebrated Keller-Segel equation from first principles using a generalized Smoluchowski equation for the full distribution function in position and orientation and a multipole expansion. An important ingredient is the bacterial tumble rate in a chemical field, for which we derived a Markovian response theory from the classical chemotaxis strategy based on temporal sensing. Our expression for the tumble rate includes logarithmic sensing, a lower chemical threshold, and bounds the chemotactic response to keep the tumble rate positive. The multipole expansion provides a polarization extended model (PE) from which we derived the Keller-Segel equation in the adiabatic limit, where the bacterial polarization instantly follows variations in the density. We thereby obtain microscopic expressions for the diffusion coefficient and the chemotactic drift velocity. Due to the bounded chemotactic response the inherent and unrealistic singularity in the drift velocity is removed.

Our detailed study of the traveling bacterial pulse shows that its characteristic time is much larger than the relaxation time of the bacterial polarization. Thus, PE and KS model provide identical results except for the initial fields and we conclude that the full Smoluchowski equation as used for example in Refs. Saragosti et al. (2011); Emako et al. (2016) is not necessary. This drastically reduces the computational effort and allowed us to perform extensive numerical studies of the bacterial pulse propagation.

We find that due to the upper bound of the chemotactic velocity the traveling pulse can only carry a limited number of bacteria. To the best of our knowledge such a maximum carrying capacity has not yet been reported in the context of traveling bacterial pulses. In particular, it is not predicted by the analytic soliton solution of the original Keller-Segel model Keller and Segel (1971). Another consequence of the upper bound of the chemotactic velocity is an effective dispersion of the pulse. While propagating, the pulse leaves a trail of bacteria behind and hence the pulse height decreases and the width expands. This is consistent with results from Ref. Scribner et al. (1974). The loss of bacteria can be compensated by a non-zero growth rate and we have seen that soliton-like pulses, which propagate with constant shape, are possible.

Exploiting a rescaled version of our KS model, we quantify how pulse speed, pulse width, and pulse amplitude depend on the different unitless parameters. We mention some key results. First, throughout our parameter study we find that the analytic soliton solution of the original KS model still provides a correct estimate for the pulse speed as a function of the number of bacteria in the pulse. Second, we find the maximum carrying capacity to be proportional to the chemotactic strength $\chi_{0}(1-\langle\cos\beta\rangle)/\tilde{\delta}$ and $\sqrt{D_{\textrm{eff}}}$. As a consequence these parameters affect the pulse speed as long as there are sufficient bacteria in the system so that the maximum carrying capacity is reached. Third, the diffusion coefficient of the chemoattractant does not influence the pulse speed as predicted in a theoretic model in Ref. Holz and Chen (1979). The pulse only takes longer to eat up all the chemoattractant at the origin due to the larger diffusive flux of chemoattractant into depleted areas.

Finally, we show that our simulated pulse propagation is able to match quantitatively the traveling bacterial pulse in the exeriments of Ref. Saragosti et al. (2011) in speed and shape. In contrast to the models used in Refs. Brenner et al. (1998) and Saragosti et al. (2011), we do not need a second chemoattractant to generate a traveling concentration pulse as a solution of our generalized KS model.

We mention four directions into which our approach can be extended. First, so far we did not explore the full PE model. In future works it would be interesting to explore the possibility of having the bacterial polarization as an independent field and its potential to induce complex dynamics. For example, the alignment or polarization of magnetotatic bacteria can be controlled by an external magnetic field Blakemore (1975); Faivre and Schuler (2008), which offers the possibility to address polarization as an independent field variable, e.g., by a time-varying external stimulus. The dependence of cell characteristics on polarization could also evoke a feedback loop in highly non-linear equations. For example, it was shown that the nutrient uptake of bacteria depends on cell shape Young (2007), meaning that the consumption rate may depend on the polarization, which then influences chemotaxis Desai et al. (2018). For Janus colloids with effective phoretic repulsion this can generate interesting collective dynamics on times much smaller than the characteristic time scale of the bacterial pulse Liebchen et al. (2015). Finally, in complex geometries with characteristic lengths similar to the persistence length of the bacterium, we expect the polarization equation also to become important.

Second, to describe the multiple pulses that have been observed in experiments with several nutrients Adler (1966b, a), one can extend our model by coupling the bacterial density to several nutrient fields. Bacteria that are left behind by the first pulse can then perform chemotaxis in a second nutrient field and thereby create a second pulse.

Third, in our generalized Smoluchowski equation the swimming speed is a constant as it is commonly done for the E.coli bacterium also during chemotaxis Berg and Brown (1972); Berg (2008); Seyrich et al. (2018). However, some bacteria are known to couple their swimming speed to the concentration of a chemical field Barbara and Mitchell (2003); Garren et al. (2014); Son et al. (2016), a strategy which is called chemokinesis. Reference Rein et al. (2016) derived coupled equations for bacterial polarization and density from a Smoluchowski equation where the swimming velocity depends on the chemical concentration. It is certainly interesting to extend our theory in order to investigate the combined effect of chemokinesis and chemotaxis. A cell that performs both strategies is the sperm cell Eisenbach (2004).

Finally, it would be interesting to extend our approach to chemoattractants to which E.coli is not perfectly adapted such as serine Berg and Brown (1972); Wong-Ng et al. (2016); Seyrich et al. (2018). For this chemoattractant the mean tumble rate drops as the concentration of serine increases. Thus, when swimming up a chemical gradient, the chemotactic velocity increases Wong-Ng et al. (2016). In our generalized KS model, the effective diffusion coefficient of bacteria $D_{\textrm{eff}}$, which directly depends on the tumble rate $\lambda_{\mathrm{equ}}$, now is enhanced in front of the pulse as runs are longer, while it is smaller in the back of the pulse where runs are shorter. This should affect the pulse propagation and indeed, experiments with serine showed that below a certain strength of the chemical gradient traveling pulses do not form Wang and Chen (1986).

ACKNOWLEDGMENTS

Fruitful discussions with J.-T. Kuhr are acknowledged. This work was supported by the German Research Foundation (DFG) within the research training group GRK 1558 and the RISE Germany Program of DAAD.

Appendix A Markovian response theory for tumble rate

Equation (1) describes the non-Markovian linear response of the tumble rate on the concentration history of the bacterial trajectory $c(\mathbf{r}(t^{\prime}))$. Now, by averaging over all possible trajectories that arrive at location $\mathbf{r}$ with orientation $\mathbf{e}$, we are able to derive a Markovian response theory for the mean tumble rate. We make use of the fact that rightward and leftward tumbling is equally probable, thus the tumble angle distribution is an even function, $P\left(\mathbf{e}-\mathbf{e^{\prime}}\right)=P\left(\mathbf{e^{\prime}}-\mathbf{e}\right)$. This has indeed been measured for E.coli in experiments Berg and Brown (1972); Pohl et al. (2017); Seyrich et al. (2018). In the following, we derive Eq. (2) from the main text.

We start with repeating the expression for the tumble rate from the linear response theory:

[TABLE]

In the following, we will use three key properties of the response function $R(\tau)$ that were measured in experiments Block et al. (1982); Segall et al. (1986); Masson et al. (2012). First, starting from $\tau=0$ it is non-zero over a time interval $\tau_{m}\lessapprox$15\text{,}\mathrm{s}$$, which we call the memory time. Second, it fulfills $\int_{-\infty}^{0}R(\tau)d\tau=0$, which means the tumble rate does not depend on the absolute chemical concentration (perfect adaptation). Third, it is inversely proportional to the adaptation concentration, $R(t-t^{\prime})=\widetilde{R}(t-t^{\prime})/c_{a}$. Adaption occurs during the memory time, thus we can set $c_{a}=c(\mathbf{r}(t))$. Taking these properties into account, we will use the approximation $\int_{-\infty}^{0}R(\tau)f(\tau)dt\approx 1/c_{a}\int_{-\tau_{m}}^{0}\tilde{R}(\tau)f(\tau)d\tau$. In particular, due to perfect adaption any additive constant in $f(\tau)$ will not contribute to the integral.

To evaluate Eq. (29), we need an expression for $c(\mathbf{r}(t^{\prime}))$. Therefore, we perform a Taylor expansion around the current position $\mathbf{r}(t)=\mathbf{r}_{t}$ and locally approximate the chemical field by

[TABLE]

where we used $\mathbf{r}^{\prime}=\mathbf{r}(t^{\prime})$. In the following derivation, all locations $\mathbf{r}^{\prime}$ that contribute to the integral in Eq. (29) should be close to the current location $\mathbf{r}_{t}$ so that the linear approximation is valid. Moreover, we assume that temporal variations of the chemical field are negligible within the memory time $\tau_{m}$ so that $\gradient{c}$ is constant. Both requirements are justified for the bacterial pulse. On the one hand, the mean run length of bacteria is much smaller than the width of the step in the chemoattractant concentration, and on the other hand on times comparable to $\tau_{m}$ the step hardly moves.

Using Eq. (30) in Eq. (30), we obtain

[TABLE]

where we applied the property of perfect adaption to set $\int_{-\tau_{m}}^{t}R(t-t^{\prime})\,\left[c(\mathbf{r}_{t})-\gradient{c}(t)\cdot\mathbf{r}_{t}\right]dt^{\prime}=0$ and that $\gradient{c}$ is constant within the memory time $\tau_{m}$.

To proceed, we write the trajectory of a bacterium that swims with constant velocity $v_{0}$ along the direction given by unit vector $\mathbf{e}(t)$ as

[TABLE]

where the bacterium has been at location $\mathbf{r}(t^{\prime})$ before reaching $\mathbf{r}(t)$, thus $t^{\prime}\leq t$. Changes in the swimming direction due to rotational diffusion are much smaller than due to tumbling. So the trajectory of the bacterium is a sequence of straight runs and instantaneous tumble events. Denoting tumble events by index $i$, the tumble time $t_{i}$, and $\mathbf{e}_{i}$ the direction prior to tumble event $i$, we can write the orientation vector in Eq. (32) at any time $t^{\prime\prime}<t$ with the help of a telescope sum:

[TABLE]

Here we set $\mathbf{e}_{0}=\mathbf{e}(t)$ for the current direction after the last tumble event $i=1$ and $\mathbf{e}_{-1}=\mathbf{0}$ is used. The number of tumble events in the time interval $t-t^{\prime\prime}$ is $n$ and we number the tumble events backwards in time.

Now, we determine the mean tumble rate $\langle\lambda(t)\rangle$ by averaging the right-hand side of Eq. (31) over an ensemble of bacterial trajectories $\mathbf{r}(t^{\prime})$ that all reach the position $\mathbf{r}(t)$ with swimming direction $\mathbf{e}(t)$. For this, we first have to evaluate $\langle\mathbf{e}(t^{\prime\prime})\rangle$ in Eq. (32) by averaging over $n$ independent tumble events and considering that $n$ is a random variable. It is determined by the probability distribution $P(n,t-t^{\prime\prime})$ of having $n$ tumble events in the time interval $t-t^{\prime\prime}$. We can thus write

[TABLE]

Note, for a constant tumble rate $P(n,t-t^{\prime\prime})$ becomes a Poisson distribution. To calculate the mean tumble direction $\langle\mathbf{e}_{i}-\mathbf{e}_{i-1}\rangle$ we use the probability distribution $P(\mathbf{e}_{i-1}-\mathbf{e}_{i})$ from the main text and calculate the first moment as in Eq. (11) but now with respect to the incoming direction $\mathbf{e}_{i}$ of the tumble event. This gives $\langle\mathbf{e}_{i}\rangle=\int\mathbf{e}_{i}P(\mathbf{e}_{i-1}-\mathbf{e}_{i})d\mathbf{e}_{i}=\langle\cos\beta\rangle\mathbf{e}_{i-1}$, where the tumble angle is determined by $\cos\beta=\mathbf{e}_{i-1}\cdot\mathbf{e}_{i}$ and we used that $P(\mathbf{e}_{i-1}-\mathbf{e}_{i})$ is an even function meaning that left- and rightward tumbles are equally probable. Repeating the formula for $\langle\mathbf{e}_{i}\rangle$ for the whole sequence of tumble events, we finally have $\langle\mathbf{e}_{i}\rangle=\langle\cos\beta\rangle^{i}\mathbf{e}(t)$ and the telescope sum in Eq. (34) becomes

[TABLE]

Combining the last two equations yields

[TABLE]

where the only remaining orientation vector is the current one, $\mathbf{e}(t)$.

Now, with the last formula and Eq. (32) we can formulate the average location

[TABLE]

Using it in the tumble rate (31), we finally obtain

[TABLE]

with

[TABLE]

Setting $c_{a}=c(\mathbf{r}$) we then recover Eq. (2) from the main text.

Appendix B Multipole expansion with bias in the tumble angle

Following the same steps as in the multipole expansion without angle bias, we average Eq. (8) over all orientations $\mathbf{e}$ using Eq. (22) and obtain

[TABLE]

Similarly, we compute the polarization $\int\mathbf{e}\,$Eq. (8)$\,d\mathbf{e}$ using Eqs. (23) and (24) and with $\hat{\mathbf{s}}=\frac{\gradient{c}}{|\gradient{c}|}$ we obtain

[TABLE]

The last term in Eq. (41) is part of the octupole moment which is defined as

[TABLE]

and represents the interplay between tumble rate variation and tumble angle bias. By neglecting all moments above the first and again defining a relaxation rate $\omega=D_{\mathrm{rot}}(d-1)+(1-\bar{\Theta})\lambda_{\mathrm{equ}}-r$ we arrive at

[TABLE]