Statistical parameter inference of bacterial swimming strategies

A typical trajectory of bacteria such as E.coli is described by a run-and-tumble random walk. During the run phase the bacterium moves forward along a nearly straight line, only rotational thermal noise affects its persistence. During the tumble phase the bacterium's speed is reduced and it reorients strongly into a new direction. The angle between the orientations before and after the tumble event is the tumble angle β. We express the velocity of the bacterium in two dimensions as the product of speed v(t) and unit vector ${\bf{e}}(t)=(\cos {\rm{\Theta }},\sin {\rm{\Theta }})$,

$$\begin{eqnarray}&&\dot{{\bf{r}}}(t)=v(t){\bf{e}}(t),\end{eqnarray} \tag{ 1 }$$

where the orientation angle Θ is measured with respect to the x axis. Figure 1(a) illustrates the different quantities.

We set up two overdamped Langevin equations for speed and orientation angle, which fully describe the bacterial motion,

$$\begin{eqnarray}&&\dot{v}(t)=r[{v}_{0}-v(t)]+{\xi }_{{\rm{sp}}}(t)+q(t),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\dot{{\rm{\Theta }}}(t)=\sqrt{2{D}_{\mathrm{rot}}(t)}\,{\xi }_{{\rm{an}}}(t).\end{eqnarray} \tag{ 3 }$$

We introduce both Langevin equations in more detail.

(1) The equation for speed v(t) contains three terms, which are associated with drift, diffusion, and jumps. We start with the last term

$$\begin{eqnarray}&&q(t)=-\displaystyle \sum _{i=1}^{{N}^{\lambda }}\eta v(t)\delta (t-{t}_{i}).\end{eqnarray} \tag{ 4 }$$

It initiates each tumble event at time t_i by a shot-noise process, while the occurrence of times t_i follows a Poisson process with tumble rate λ. At the beginning of each tumble, the bacterial speed is reduced by the relative jump height η to $(1-\eta ){v}_{t}$ and N^λ is the actual number of tumble events. The first and second term represent a conventional Ornstein–Uhlenbeck process. After a tumble event the speed relaxes with relaxation rate r towards the swimming speed v₀ of the run phase. Thus r⁻¹ is the mean duration of a tumble event, which we call tumble time in the following. The Gaussian white noise term is fully determined by $\langle {\xi }_{{\rm{sp}}}\rangle =0$ and $\langle {\xi }_{{\rm{sp}}}(s){\xi }_{{\rm{sp}}}(t)\rangle ={\sigma }^{2}\delta (t-s)$, where we introduce the white noise strength σ. It describes the ubiquitous noise due to internal noise of the swimming mechanism and variations between individual bacteria.

Note that the actual tumble time of a bacterium is exponentially distributed. In our model the white noise term also induces stochastic fluctuations in the duration of the tumble events as visible in figure 2(b). Altogether, the stochastic speed process is determined by five parameters: $\{\lambda ,r,{v}_{0},\sigma ,\eta \}$.

Zoom In Zoom Out Reset image size

(2) The stochastic equation for the orientation angle Θ is fully described by rotational diffusion, where the white noise process is defined by $\langle {\xi }_{{\rm{an}}}\rangle =0$ and $\langle {\xi }_{{\rm{an}}}(s){\xi }_{{\rm{an}}}(t)\rangle =\delta (t-s)$. Following [31], we model tumbles as a random walk on a unit sphere with enhanced rotational diffusion. Thus, the rotational diffusion coefficient ${D}_{{\rm{rot}}}(t)$ is no longer a constant but alternates between two values: the thermal rotational diffusion coefficient D₀ during run phases and an enhanced value D_T during tumble phases. We describe each transition between the two states by a Poisson process and thus obtain a telegraph process. The transition rate from the run to the tumble phase is the tumble rate λ, whereas the transition rate in the opposite direction is the speed relaxation rate r or the inverse tumble time. A full definition and basic properties of the telegraph process are given in the appendix B.2 or can be found in [25].

To link the telegraph process to the shot-noise process for the speed value in equation (4), the diffusion coefficient switches at the same times t_i from the thermal (D₀) to the enhanced (D_T) value. Note, while the speed process allows a second tumble although the first one is not finished yet, this is not possible in the telegraph process for rotational diffusion. However, for bacteria like E.coli the time between tumble events is typically one order of magnitude larger than the tumble time r⁻¹. This makes these double events very rare and tumble events in both speed and angular processes coincide. All in all, we have four parameters governing the stochastic process for the orientation angle: $\{\lambda ,r,{D}_{0},{D}_{T}\}$.

Figure 2 shows a typical simulated trajectory (a) and the corresponding time series for speed and angular displacement ΔΘ during time step Δt = 0.1 s (b). It has to be compared to the experimental time series of both quantities in (c). Note that $\tfrac{{\rm{\Delta }}{\rm{\Theta }}}{{\rm{\Delta }}t}$ represents the turning rate of the bacterium. In the following we will always work with the angular displacement ΔΘ.

In this section we state moments, stationary distributions, and time autocorrelation functions for the stochastic processes of speed and orientation angle in equations (2) and (3). They depend on the swimming parameters introduced above. Matching the theoretical expressions of these quantities to the values determined by averaging over all individual tracks of the experiments, we are able to infer the mean swimming parameters of an E.coli population. We refer to appendix B for details of the derivations and only state the final expressions in the following. Finally we note, while the calculated moments and distribution functions are valid for any point process with arbitrary run-time distribution, the formulas for the autocorrelation functions are only valid for a Poisson process, which gives an exponential distribution for the run times. In appendix G we demonstrate this for a Γ distribution.

2.2.1. Speed

The moments ${m}_{n}^{V}=\langle v{(t)}^{n}\rangle$ of equation (2), where the average is taken over all times t and all tracks in the long-time limit, can be calculated as a function of the reduced parameter set $(\lambda /r,\eta ,{v}_{0},{\sigma }^{2}/r)$. For the first moment, the mean speed, we obtain

$$\begin{eqnarray}&&{m}_{1}^{V}\left(\displaystyle \frac{\lambda }{r},\eta ,{v}_{0},\displaystyle \frac{{\sigma }^{2}}{r}\right)=\displaystyle \frac{{v}_{0}}{1+\eta \lambda /r}\,.\end{eqnarray} \tag{ 5 }$$

The mean speed is smaller than the swimming speed v₀ since during the tumble phase speed is reduced by a factor η. More generally, a recursive formula for the nth moment is given by

$$\begin{eqnarray}&&{m}_{n}^{V}\left(\displaystyle \frac{\lambda }{r},\eta ,{v}_{0},\displaystyle \frac{{\sigma }^{2}}{r}\right)=\displaystyle \frac{{v}_{0}\,{m}_{n-1}^{V}\,+\,\tfrac{1}{2}(n-1)\tfrac{{\sigma }^{2}}{r}\,{m}_{n-2}^{V}}{1+\tfrac{\lambda }{{nr}}-\tfrac{\lambda }{{nr}}{(1-\eta )}^{n}}\,,\end{eqnarray} \tag{ 6 }$$

where the zeroth moment is m₀ = 1 due to normalization. We now have access to all the speed moments. As an example, figure 3(a) shows a histogram for the distribution of speed values recorded in an experiment, from which the speed moments can be calculated. The orange line represents the distribution obtained from numerically solving the speed equation (2) using the actual parameters inferred from this experiment. The two distributions nicely agree, which is an a posteriori verification of our Langevin equation.

Zoom In Zoom Out Reset image size

From the moments we can only infer the ratios λ / r and σ²/r. In order to determine the full set of parameters of equation (2), we also use the speed autocorrelation function for our model. It has an exponential form with relaxation rate r + ηλ,

$$\begin{eqnarray}&&{g}_{V}(\tau )=\langle [v(t+\tau )-{m}_{1}^{V}][v(t)-{m}_{1}^{V}]\rangle ={{\rm{\Delta }}}^{2}{{v}{\rm{e}}}^{-(r+\eta \lambda )\tau }\,,\end{eqnarray} \tag{ 7 }$$

where we have introduced the variance ${{\rm{\Delta }}}^{2}v=\langle {(v-{m}_{1}^{V})}^{2}\rangle$. Figure 3(b) shows the autocorrelation function for the experimental data of E.coli. Indeed, the curve is well-fitted by an exponential over two decades up to $\tau \simeq 1s$, which is around half the mean track length. This agreement supports the validity of our stochastic description of the speed process in equation (2).

2.2.2. Angle

Here, we work directly with the steady-state probability distribution $p(| {\rm{\Delta }}{\rm{\Theta }}| )$ for the absolute angular displacement $| {\rm{\Delta }}{\rm{\Theta }}|$ during a finite time step Δt. We determine $p(| {\rm{\Delta }}{\rm{\Theta }}| )$ from equation (3) for the orientation angle as a function of the reduced parameter set $(\lambda /r,{D}_{0},{D}_{T})$. In the long-time limit the probability distribution $p(| {\rm{\Delta }}{\rm{\Theta }}| )$ becomes stationary and is given by

$$\begin{eqnarray}&&p(| {\rm{\Delta }}{\rm{\Theta }}| )=\displaystyle \frac{r}{\lambda +r}{ \mathcal N }(0,\sqrt{2{D}_{0}{\rm{\Delta }}t})+\displaystyle \frac{\lambda }{\lambda +r}{ \mathcal N }(0,\sqrt{2{D}_{T}{\rm{\Delta }}t}),\end{eqnarray} \tag{ 8 }$$

where ${ \mathcal N }(0,\sigma )$ denotes the normal distribution with zero mean and standard deviation σ. For our parameter inference we use the same time step ${\rm{\Delta }}t=0.1{{\rm{s}}}^{-1}$ as in [31].

Figure 4(a) presents a histogram for all angular displacements in time step Δt recorded in the experiment. It shows a deviation from the theoretical distribution of equation (8) in the tail at angles larger than π/2, which is visible only in the semi-logarithmic plot. Note that the region $| \beta | \gt \pi /2$ only represents roughly 3% of all angular displacements. There are two possible reasons for this deviation: First, we record angular displacements Θ = π +  as a displacement −(π − ) since we cannot distinguish between tumbles to the right and left during one time step. Second, it is also possible that the diffusion model for tumbling does not apply for such large angles.

Zoom In Zoom Out Reset image size

For completeness we also give the nth moment of the absolute angular displacement, ${m}_{n}^{{\rm{\Delta }}{\rm{\Theta }}}=\langle | {\rm{\Delta }}{\rm{\Theta }}{| }^{n}\rangle$. It follows directly from the probability distribution of equation (8):

$$\begin{eqnarray}\begin{array}{rcl}{m}_{n}^{{\rm{\Delta }}{\rm{\Theta }}}(\lambda ,r,{D}_{0},{D}_{T}) & = & \left(\displaystyle \frac{{(2{D}_{0}{\rm{\Delta }}t)}^{\tfrac{n}{2}}}{1+\lambda /r}+\displaystyle \frac{{(2{D}_{T}{\rm{\Delta }}t)}^{\tfrac{n}{2}}}{1+r/\lambda }\right)(n-1)!!\\ & & \cdot \left\{\begin{array}{ll}\sqrt{\displaystyle \frac{2}{\pi }}\,\,\, & {\rm{if}}\,{\rm{n}}\,{\rm{is}}\,{\rm{odd}}\\ 1\,\,\, & {\rm{if}}\,{\rm{n}}\,{\rm{is}}\,{\rm{even}}\end{array},\right.\end{array}\end{eqnarray} \tag{ 9 }$$

where $n!!$ denotes the double factorial.

Similar to the speed process, we can only infer the ratio λ / r from fits to the probability distribution $p(| {\rm{\Delta }}{\rm{\Theta }}| )$ of equation (8). In order to determine the full set of parameters of equation (3), we use again the autocorrelation function of our model, now for the swimming direction e(t). Numerical investigations of our model (see appendix B.3) suggest that it has a simple exponential form with relaxation rate α_Θ for parameters relevant to the experiments:

$$\begin{eqnarray}&&{g}_{{\rm{\Theta }}}(\tau )=\langle {\bf{e}}(t+\tau )\cdot {\bf{e}}(t)\rangle \propto {{\rm{e}}}^{-{\alpha }_{{\rm{\Theta }}}\tau }.\end{eqnarray} \tag{ 10 }$$

Analytically, we are not able to calculate this exponential form. However, in the time interval ${(\lambda +r)}^{-1}\lt \tau \lt \langle {D}_{\mathrm{rot}}{\rangle }^{-1}$ relevant to the experiments, we can derive the linear approximation

$$\begin{eqnarray}&&{g}_{{\rm{\Theta }}}(\tau )\approx 1-{\alpha }_{{\rm{\Theta }}}\tau \approx 1-\left(\langle {D}_{\mathrm{rot}}\rangle -\displaystyle \frac{{{\rm{\Delta }}}^{2}{D}_{\mathrm{rot}}}{\lambda +r}\right)\tau \end{eqnarray} \tag{ 11 }$$

and thereby obtain an expression for the relaxation rate α_Θ. Here we have introduced the respective mean $\langle {D}_{\mathrm{rot}}\rangle$ and variance ${{\rm{\Delta }}}^{2}{D}_{\mathrm{rot}}$ of the telegraph process ${D}_{\mathrm{rot}}(t)$,

$$\begin{eqnarray}\begin{array}{rcl}\langle {D}_{\mathrm{rot}}\rangle & = & \displaystyle \frac{{D}_{0}}{1+\lambda /r}+\displaystyle \frac{{D}_{T}}{1+r/\lambda }\\ {{\rm{\Delta }}}^{2}{D}_{\mathrm{rot}} & = & \langle {({D}_{\mathrm{rot}}-\langle {D}_{\mathrm{rot}}\rangle )}^{2}\rangle =\displaystyle \frac{{({D}_{0}-{D}_{T})}^{2}\lambda /r}{{(1+\lambda /r)}^{2}}.\end{array}\end{eqnarray} \tag{ 12 }$$

Figure 4(b) shows the directional autocorrelation function for the experimental data of E.coli moving in a uniform buffer medium. Indeed, the curve is well-fitted by an exponential up to τ ≃ 5 s excluding the first point. This agreement supports the validity of our stochastic description of the angle process in equation (3). The deviation in the experimental data for the first point is caused by the offset for angular displacements larger than π/2, where the experimental distribution function in figure 4(a) deviates from theory. For two and more time steps the influence of this offset becomes smaller and smaller.

E.coli AW405 strain was cultured overnight in liquid Tryptone Broth (TB) (10 g l⁻¹ Difco Bacto^TM-Tryptone and 5 g l⁻¹ NaCl) at 37 °C on a rotary shaker at 300 rpm. The cell suspension was diluted 1:100 into fresh TB, and grown to mid-exponential phase (OD₆₀₀ = 0.5). Then the bacterial suspension was washed two times in order to remove any remaining medium from the growth medium following [3, 32] and resuspended in motility buffer (11.2 g l⁻¹ K₂HPO₄, 4.8 g l⁻¹ KH₂PO₄, 3.93 g l⁻¹ NaCl, 0.029 g l⁻¹ EDTA and 0.5 g l⁻¹ glucose; pH 7.0). Afterward, the cell suspension was divided into two fractions. One was centrifuged and resuspended in the same motility buffer, and the other was centrifuged and resuspended in motility buffer supplemented with the chemoattractant α-methyl-aspartate (Sigma-Aldrich, USA) in a final concentration of 0.5 mM. In both cases, the final OD₆₀₀ of the cell suspensions was 0.07 before filling them into chemotaxis chambers.

In this study, a µ-Slide Chemotaxis 3D (ibidi, Martinsried, Germany) was used in order to maintain a stable linear gradient of the chemoattractant α-methyl-aspartate. This chemotaxis chamber consists of two large reservoirs connected to a central observation area (see figure 1(b)). For the chemotaxis assay, the cell suspension with chemoattractant was filled into the reservoir on the right-hand side and the chemoattractant-free cell suspension into the reservoir on the left-hand side. The central observation area was filled with motility buffer (see appendix A). A stable linear chemoattractant gradient is generated by diffusion in the observation area and maintained for several hours [33]. For the control assay, both reservoirs were filled with chemoattractant-free cell suspension. In this case, a homogeneous environment without any gradient was established in the observation area.

An IX71 inverted microscope with a 20x UPLFN-PH objective (both Olympus, Germany) in phase contrast mode was used for imaging cell trajectories. Five image sequences were taken with 10 min intervals between them using a Orca Flash 4.0 CMOS camera (Hamamatsu Photonics, Japan). For each sequence, the images were acquired at 20 frames per second for 30 s. The field of the view was placed in the center of the gradient region at 30 μm above the bottom of the chamber (total height in the observation area was 70 μm).

A custom Matlab program based on the Image Processing Toolbox (version R2015a, The MathWorks, USA) was used to process the image sequences automatically. For each image sequence, a background image was calculated by pixel-wise time average projection. It was subtracted from each frame to eliminate non-motile objects and shading effects. The built-in Matlab function imerode was then applied for morphological erosion (with a disk of radius 0.6 μm) to reduce the background noise. The putative bacterial cells are distinguished from background using the maximum entropy thresholding algorithm by Kapur et al [34]. The threshold was calculated for each image in the sequence separately. The median of all threshold values was used to segment the whole sequence. The binary images were further processed with the morphological operations, imopen and imclose (with a disk of radius 0.3 μm) to eliminate any noise caused by segmentation. The builtin function bwconncomp was used to find all connected objects in the binary images. Size and centroid of the objects were determined using the regionprops function. Afterwards, particles with an area between 1 and 15.6 μm² were considered as single bacterial cell. Finally, trajectories were obtained employing the tracking algorithms by Crocker and Grier [35].

To avoid tracking artifacts caused by tumble events when cells enter and leave the focal plane, the first and last 0.5 s of each track were removed. We use the curvatures of the trajectories to filter out bacteria that circle close to the cover slip. The curvature at each data point is calculated and after that the median curvatures for each single trajectory. Finally, the trajectories with the highest 5% of median values are removed. We checked that our results do not change if we go up to 20%. Also, tracks with a total displacement <10 μm were eliminated, as they most likely result from damaged flagella. The minimal track length is 0.5 s and the maximal length is 19.35 s. The control data set consists of 769 tracks with a total length of 1629 s. The gradient data set consists of 3498 tracks with a total length of 7206 s.

The trajectories were smoothed using a second-order Savitztky–Golay filter with a window size of 5 data points corresponding to 250 ms [36]. Instantaneous speed $v=\tfrac{{\rm{\Delta }}s}{{\rm{\Delta }}t}$, direction of propagation θ, and turning rate $\omega =\tfrac{{\rm{\Delta }}\theta }{{\rm{\Delta }}t}$ were evaluated on the smoothed tracks. The tumble events were detected as described previously [9, 22, 37]. Briefly, in the time series of speed and turning rate, local minima and local maxima were detected, respectively, to identify tumble events. Four parameters, two for the speed and two for the turning rate, were adjusted such that the recognition of tumble events was correct as checked by visual examination (threshold parameters α = 3 and β = 6.5 and tumble duration parameters 0.55 × Δv and 0.65 × Δω, see the supporting information S5 in [22]).

Figures (3) and (4) show distributions and autocorrelation functions for speed and angular displacements recorded for E.coli when swimming in a homogeneous buffer without any chemical gradient. Note that speed and angle inference are performed separately from each other but they are linked by the tumble rate λ and the inverse tumble time r.

4.1.1. Speed inference

4.1.2. Angle inference

4.1.3. Inferred parameters

Table 1 gives an overview of the inferred swimming parameters for the two stochastic processes for speed and angle. The two inferred tumble rates λ are very close together and the inverse tumble times r agree within the error bars. Our results are in good agreement with tumble rate λ = 0.84 s⁻¹ and swimming velocity ${v}_{0}=20.7\,\mu {\mathrm{ms}}^{-1}$ determined with a heuristic tumble recognizer (see section 3.4 and [22]). This validates our inference method. Moreover, our findings are in good agreement with previously measured tumble rates [3, 9, 22] and swimming speeds [38]. The inferred value for the thermal rotational diffusivity D₀ agrees with previously reported values in the literature, which range from 0.06 s⁻¹ [20, 22] to 0.18 s⁻¹ [39].

Table 1.  Inferred parameters for the stochastic processes of speed and angle for E.coli moving in a buffer medium without a chemical gradient (control experiment).

Speed			Angle
λ	0.83 ± 0.04 s−1		λ	0.84 ± 0.02 s−1
r	4.41 ± 0.3 s−1		r	3.81 ± 0.3 s−1
v0	20.8 ± 0.2 μm s−1		D0	0.09 ± 0.002 s−1
$\sqrt{\tfrac{{\sigma }^{2}}{r}}\,$	5.11 ± 0.07 μm s−1		DT	2.31 ± 0.12 s−1
η	0.85 ± 0.01

We use the enhanced rotational diffusion coefficient ${D}_{T}=2.31\,{{\rm{s}}}^{-1}$ and the inverse tumble time r = 3.81 s⁻¹ of the angle stochastic process to determine the distribution function of absolute tumble angles, $P(| \beta | )$, by recording the angular displacement for exponentially distributed tumble times with mean r⁻¹. The corresponding three-dimensional distribution function is obtained by multiplying the two-dimensional quantity with $\sin \beta$ from the solid angle element. The resulting distribution is shown in orange in figure 5 for $| \beta | \lt \pi$. It has a maximum at β_max = 0. 78 = 45° and the mean tumble angle is $\langle | \beta | \rangle =1.06=61^\circ$, which are remarkably close to the values β_max = 45° and $\langle | \beta | \rangle =62^\circ$ from [3]. The shape of the distribution function is similar to the one obtained with the heuristic tumble recognizer (blue bars). Also, the maximum values are very close. While the main characteristics of the two curves agree well, the heuristic tumble recognizer determines more tumbles for angles close to π. As a result, it finds a larger mean tumble angle $\langle | \beta | \rangle =1.43=82^\circ$. This might be explained as follows. Some tumbles occur only in one time interval, where one cannot distinguish between a leftward tumble angle $\tilde{\beta }$ and a rightward tumble $| \tilde{\beta }-2\pi |$. Thus, the heuristic tumble recognizer chooses always the smaller angle and, therefore, the distribution of tumble angles close to π is enhanced. In contrast, our inference for the angle process only uses angular displacements up to π/2 and thereby avoids differences of P(ΔΘ) between theory and experiment for values larger than π/2 , as discussed in connection with figure 4(a). Therefore, we think that our inferred tumble angle distribution gives a better account of $P(| \beta | )$.

Zoom In Zoom Out Reset image size

Compared to literature we define the tumble time differently by setting ${\tau }_{t}={r}^{-1}$. Usually, one employs a tumble recognizer and identifies the tumble state when the angular displacement (per time step) exceeds a threshold value [3, 7, 8, 16]. The duration of this period is then the tumble time (see also figure 6(a)), for which values of τ_t = 0.12 s and 0.14 s were measured using different thresholds [3, 8]. However, this procedure underestimates the duration of a tumble event, which starts when a flagellum leaves the bundle and ends when it returns to the bundle. At the beginning and end of this period the angular displacement (per time step) can of course be below the given threshold value. Indeed, [16] showed that the duration of a tumble event obtained from visualizing the flagellar dynamics during tumbling is significantly larger than the time determined by tumble recognizers.

Zoom In Zoom Out Reset image size

In contrast to tumble recognizers, our method defines the tumble time as the inverse relaxation rate ${\tau }_{t}={r}^{-1}$. This is a more rational quantification of the tumble time without the need of an a priori threshold value. Tumbles are initiated when the speed jumps below the swimming speed and they end when the speed has relaxed back to the swimming speed. We argue that the higher value τ_t = 0.23 s obtained by our method describes the tumble process more precisely.

Next, we apply our method to experimental data of E.coli recorded in a constant gradient of a chemoattractant concentration. Conditioning the analysis on the swimming direction, we are able to determine how the swimming parameters depend on the orientation or swimming angle θ. Thus, we divide the experimental data into eight subsets or sectors each spanning a range of orientation angles centered at ${\theta }_{n}=2\pi n/8$ for n = {0, 1, ..., 7 }. Here, θ = 0, 2π means swimming up the gradient and θ = π against the gradient. In practice, instead of dividing the data for the orientation angle into 8 disjunct sectors, we use smooth weighting based on Gaussian kernels as in [22] (for further details see appendix D).

Figure 7 shows the results from applying our inference method to the moments of speed and to the distribution of angular displacements. Graph (a) plots the tumble bias λ/r, the ratio of tumble time to run time, versus orientation angle. It is lowered when swimming up the gradient (θ = 0, 2π) and increased when swimming down the gradient (θ = π). This confirms the classical chemotaxis strategy. The curves from angle inference (orange) and speed inference (blue) show good agreement. Again, we recognize that both inference strategies give coherent results, even though they are performed independently from each other. In figure 7(b) the rotational diffusion coefficient D_T during tumbling also depends on the swimming direction. It is lowered when swimming up the gradient and increased when swimming down the gradient. This suggests angular persistence or a reduced mean tumble angle, when swimming in a favorable direction, as a chemotaxis strategy. It was already reported in [8, 22]. We will comment more on this strategy in the following.

Zoom In Zoom Out Reset image size

Adding the speed autocorrelation function to the parameter inference, we investigate whether tumble rate λ and tumble time r⁻¹ are separately modulated during chemotaxis. Figure 8 shows the results for the speed parameters λ, r, v₀, σ. Indeed, we recover the classical chemotaxis strategy in plot (a) with a strong reduction of the tumble rate when swimming up the chemical gradient. The tumble rate for θ = 0 is less than half of the tumble rate for θ = π. The same trend occurs for the tumble time r⁻¹, which increases when swimming down the gradient. This bias in tumble time together with the same trend for the diffusion coefficient D_T found above confirms a bias in the mean tumble angle $\langle \beta \rangle$. It is enhanced when swimming in an unfavorable direction, which confirms the alternative chemotaxis strategy identified in [8, 22]. No significant modulations are visible for the swimming speed v₀ plotted in (c). So there is no chemokinesis. The same applies to the jump height η, which is shown in figure E1 of appendix E. In the last plot (d) we report a bias in speed fluctuations. The swimming speed is significantly more volatile when swimming up a chemical gradient compared to swimming against it. We are not aware that this has been reported so far and it is also not clear to us what purpose it serves.

Zoom In Zoom Out Reset image size

In order to perform the derivations, we rewrite Langevin equation (2) as a stochastic differential equation (SDE) using mathematical notation:

$$\begin{eqnarray}&&{{\rm{d}}{v}}_{t}=r({v}_{0}-{v}_{t}){\rm{d}}{t}+\sigma {{\rm{d}}{W}}_{t}-\eta {v}_{t}{{\rm{d}}{N}}_{t}^{\lambda }\,.\end{eqnarray} \tag{ B1 }$$

Here, we define the Poisson process where ${{\rm{d}}{N}}_{t}^{\lambda }=1$ occurs with probability λ dt for each time step indicating the start of a tumble and ${{\rm{d}}{N}}_{t}^{\lambda }=0$ otherwise. Moreover, we introduce the Wiener process dW_t. Integrating equation (B1) and splitting the Poisson process into a deterministic part and a fluctuating part ${{\rm{d}}{N}}_{t}^{\lambda }=\lambda {\rm{d}}{t}+{\rm{d}}{\tilde{N}}_{t}^{\lambda }$ [43] yields

$$\begin{eqnarray}&&{v}_{t}={\int }_{0}^{t}r({v}_{0}-{v}_{s}){\rm{d}}{s}+{\int }_{0}^{t}\sigma {{\rm{d}}{W}}_{s}-{\int }_{0}^{t}\eta {v}_{s}\lambda {\rm{d}}{s}-{\int }_{0}^{t}\eta {v}_{s}{\rm{d}}{\tilde{N}}_{t}^{\lambda }.\end{eqnarray} \tag{ B2 }$$

Note that the second and fourth term on the rhs are martingales [43]. Thus, their expectation values vanish. We will use this property when calculating the moments and autocorrelation function of the speed variable. Taking the expectation value $\langle \ldots \rangle$ on both sides, we obtain the first moment:

$$\begin{eqnarray}&&{m}_{1}=\langle {v}_{t}\rangle ={{rv}}_{0}t-{\int }_{0}^{t}(r+\eta \lambda )\langle {v}_{s}\rangle {\rm{d}}{s}.\end{eqnarray} \tag{ B3 }$$

To ease the notation, we dropped the superscript V from the main text. Taking the time derivative on both sides, we obtain a non-homogeneous ordinary differential equation (ODE):

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{m}}_{1}}{{\rm{d}}{t}}={{rv}}_{0}-(r+\eta \lambda ){m}_{1}.\end{eqnarray} \tag{ B4 }$$

Its full solution with initial value C at time t₀ reads

$$\begin{eqnarray}&&{m}_{1}(t)=\displaystyle \frac{{v}_{0}}{1+\eta \tfrac{\lambda }{r}}+{{\rm{e}}}^{-(r+\eta \lambda )(t-{t}_{0})}\left(C-\displaystyle \frac{{v}_{0}}{1+\eta \tfrac{\lambda }{r}}\right)\,.\end{eqnarray} \tag{ B5 }$$

Taking the long-time limit $t\to \infty$, we recover the equation (5) from the main text.

Next, we calculate the nth moment ${m}_{n}=\langle {v}^{n}\rangle$. Using Ito's lemma [44], we first formulate a SDE for an arbitrary function f(v_t) of the speed variable:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{d}}{f}({v}_{t}) & = & \left(f^{\prime} ({v}_{t})r({v}_{0}-{v}_{t})+\displaystyle \frac{1}{2}f^{\prime\prime} ({v}_{t}){\sigma }^{2}\right){\rm{d}}{t}\\ & & +f^{\prime} ({v}_{t})\sigma {{\rm{d}}{W}}_{t}+[f({v}_{t-}-\eta {v}_{t-})-f({v}_{t-})]{{\rm{d}}{N}}_{t}^{\lambda }.\end{array}\end{eqnarray} \tag{ B6 }$$

Here, ${v}_{t-}$ denotes the value right before a jump. Setting $f({v}_{t})={v}_{t}^{n}$, integrating equation (B6), and taking the expectation value on both sides yields:

$$\begin{eqnarray}\begin{array}{rcl}\langle {v}_{t}^{n}\rangle & = & {\displaystyle \int }_{0}^{t}({{nrv}}_{0}\langle {v}_{s}^{n-1}\rangle +(\lambda [{(1-\eta )}^{n}-1)]-{nr})\langle {v}_{s}^{n}\rangle \\ & & +\left.\displaystyle \frac{n(n-1)}{2}{\sigma }^{2}\langle {v}_{s}^{n-2}\rangle \right){\rm{d}}{s},\end{array}\end{eqnarray} \tag{ B7 }$$

where we again extracted the deterministic part of the Poisson process and all martingales dropped out. Taking the time derivative on both sides, we obtain an ODE, which also contains the lower-order moments ${m}_{n-1}$ and ${m}_{n-2}$:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{\rm{d}}{m}}_{n}}{{\rm{d}}{t}} & = & \,{{nrv}}_{0}\,{m}_{n-1}+(\lambda [{(1-\eta )}^{n}-1]-{nr})\,{m}_{n}\\ & & +\displaystyle \frac{n(n-1)}{2}{\sigma }^{2}\,{m}_{n-2}.\end{array}\end{eqnarray} \tag{ B8 }$$

The solution of this ODE in the long-time limit $t\to \infty$, where ${{\rm{d}}{m}}_{n}/{\rm{d}}{t}=0$, yields equation (6) in the main text

$$\begin{eqnarray}&&{m}_{n}=\displaystyle \frac{{v}_{0}\,{m}_{n-1}\,+\,\tfrac{1}{2}(n-1)\tfrac{{\sigma }^{2}}{r}\,{m}_{n-2}}{1+\tfrac{\lambda }{{nr}}-\tfrac{\lambda }{{nr}}{(1-\eta )}^{n}}\,.\end{eqnarray} \tag{ B9 }$$

Finally, we calculate the speed autocorrelation function $g(s,t)=\langle ({v}_{s}-{m}_{1})({v}_{t}-{m}_{1})\rangle$ of equation (B1). We define the probability distributions for the speed process $P(v^{\prime} )$ and the conditional probability $P(v,t| v^{\prime} ,s)$ of having v at time t given that we have $v^{\prime}$ at time s and obtain

$$\begin{eqnarray}\begin{array}{rcl}g(s,t) & = & \displaystyle \int \int [(v-{m}_{1})(v^{\prime} -{m}_{1})P(v,t| v^{\prime} ,s)P(v^{\prime} )]{\rm{d}}{v}{\rm{d}}{v}^{\prime} \\ & = & \displaystyle \int [\langle v(t)-{m}_{1}| [v^{\prime} ,s]\rangle (v^{\prime} -{m}_{1})P(v^{\prime} )]{\rm{d}}{v}^{\prime} \\ & = & \displaystyle \int [((v^{\prime} -{m}_{1}){{\rm{e}}}^{-(r+\eta \lambda )| t-s| })(v^{\prime} -{m}_{1})P(v^{\prime} )]{\rm{d}}{v}^{\prime} \\ & = & {\rm{\Delta }}{v}^{2}\,{{\rm{e}}}^{-(r+\eta \lambda )| t-s| }\,,\end{array}\end{eqnarray} \tag{ B10 }$$

where we have have used equation (B5) with $C=v^{\prime}$ in the second last step. We recover equation (7) after setting s = t + τ. Identifying the relaxation rate α_V, we can write the following formulas for λ and r:

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{{\alpha }_{V}}{\eta +{(\lambda /r)}^{-1}}\,,\end{eqnarray} \tag{ B11 }$$

$$\begin{eqnarray}&&r=\displaystyle \frac{{\alpha }_{V}}{1+\eta \lambda /r}\,.\end{eqnarray} \tag{ B12 }$$

We rewrite the Langevin equation (3) from the main text as a SDE using mathematical notation:

$$\begin{eqnarray}&&{\rm{d}}{{\rm{\Theta }}}_{t}=\sqrt{2{D}_{t}}{{\rm{d}}{W}}_{t}.\end{eqnarray} \tag{ B13 }$$

The SDE contains two stochastic processes: the telegraph process D_t, where we drop here the subscript $\mathrm{rot}$ used in the main text, and the white noise process dW_t. These two processes are stochastically independent of each other. Thus, the moments for the angular displacement during time step Δt factorize into contributions from each process

$$\begin{eqnarray}&&\langle | {\rm{\Delta }}{\rm{\Theta }}{| }^{n}\rangle =\langle {[2{D}_{t}]}^{\tfrac{n}{2}}\rangle \langle | {\rm{\Delta }}{W}_{t}{| }^{n}\rangle \,.\end{eqnarray} \tag{ B14 }$$

The pdf $p({\rm{\Delta }}{W}_{t})$ and the absolute moments of the white noise increments ΔW_t during time step Δt are given by

$$\begin{eqnarray}&&p({\rm{\Delta }}{W}_{t})={ \mathcal N }(0,\sqrt{{\rm{\Delta }}t})\,,\end{eqnarray} \tag{ B15 }$$

$$\begin{eqnarray}\langle | {\rm{\Delta }}W(t){| }^{n}\rangle ={({\rm{\Delta }}t)}^{\tfrac{n}{2}}(n-1)!!\left\{\begin{array}{ll}\sqrt{\displaystyle \frac{2}{\pi }} & \mathrm{if}\ \ n\ \ \mathrm{is}\ \ \mathrm{odd}\\ 1 & \mathrm{if}\ \ n\ \ \mathrm{is}\ \ \mathrm{even}\end{array}\right.,\end{eqnarray} \tag{ B16 }$$

where ${ \mathcal N }(0,\sigma )$ denotes the normal distribution with zero mean and standard deviation σ and $n!!$ denotes the double factorial.

For the telegraph process D_t with states D₀ and D_T, the two probabilities for being in one of the states at time t obey the following master equations:

$$\begin{eqnarray*}\begin{array}{rcl}{\partial }_{t}P({D}_{0},t| C,{t}_{0}) & = & -\lambda P({D}_{0},t| C,{t}_{0})+{rP}({D}_{T},t| C,{t}_{0}),\\ {\partial }_{t}P({D}_{T},t| C,{t}_{0}) & = & \,\,\,\lambda P({D}_{0},t| C,{t}_{0})-{rP}({D}_{T},t| C,{t}_{0}).\end{array}\end{eqnarray*}$$

Here, λ is the transition rate from D₀ to D_T and r the transition rate for the reverse process. The variable C indicates the initial condition at time t₀. We first state the pdf p(D) in the long-time limit $t\to \infty$ as well as the auto-correlation function $\langle {D}_{t}{D}_{s}\rangle$ from literature [25]:

$$\begin{eqnarray}&&p({D}_{0})=\displaystyle \frac{r}{\lambda +r}\,,\end{eqnarray} \tag{ B17 }$$

$$\begin{eqnarray}&&p({D}_{T})=\displaystyle \frac{\lambda }{\lambda +r}\,,\end{eqnarray} \tag{ B18 }$$

$$\begin{eqnarray}&&\langle {D}_{t}{D}_{s}\rangle =\langle D{\rangle }^{2}+{{\rm{\Delta }}}^{2}D\,{{\rm{e}}}^{-(\lambda +r)| t-s| }\,.\end{eqnarray} \tag{ B19 }$$

In the last equation we have introduced the mean $\langle D\rangle$ and the variance Δ²D in the long time limit. They are given by

$$\begin{eqnarray}&&\langle D\rangle =\displaystyle \frac{{D}_{0}r+{D}_{T}\lambda }{\lambda +r}\,,\end{eqnarray} \tag{ B20 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{2}D=\displaystyle \frac{{({D}_{0}-{D}_{T})}^{2}\lambda r}{{(\lambda +r)}^{2}}\,.\end{eqnarray} \tag{ B21 }$$

The mean value of D_t for any time t with initial condition C at time t₀ is given by

$$\begin{eqnarray}&&\langle {D}_{t}\rangle =\langle D\rangle +{{\rm{e}}}^{-(\lambda +r)(t-{t}_{0})}(C-\langle D\rangle )\,.\end{eqnarray} \tag{ B22 }$$

We can use the pdf p(D) to calculate the first factor on the RHS of equation (B14) in the long time limit

$$\begin{eqnarray}&&\langle {[2{D}_{t}]}^{\tfrac{n}{2}}\rangle =\displaystyle \frac{{(2{D}_{0})}^{\tfrac{n}{2}}}{1+\lambda /r}+\displaystyle \frac{{(2{D}_{T})}^{\tfrac{n}{2}}}{1+r/\lambda }\,.\end{eqnarray} \tag{ B23 }$$

Inserting this expression and equations (B16) in (B14) leads to equation (9) stated in the main text.

The pdf of the absolute angular displacement $p(| {\rm{\Delta }}{\rm{\Theta }}| )$ can be calculated straightforwardly. Using the independence of the two stochastic processes and combining equations (B15), (B17), and (B18), we obtain

$$\begin{eqnarray}&&p(| {\rm{\Delta }}{\rm{\Theta }}| )=\displaystyle \frac{r}{\lambda +r}{ \mathcal N }(0,\sqrt{2{D}_{0}{\rm{\Delta }}t})+\displaystyle \frac{\lambda }{\lambda +r}{ \mathcal N }(0,\sqrt{2{D}_{T}{\rm{\Delta }}t})\,.\end{eqnarray} \tag{ B24 }$$

This agrees with equation (8) from the main text.

Finally, we calculate the directional autocorrelation function $g(\tau )=\langle {\bf{e}}(\tau )\cdot {\bf{e}}(0)\rangle =\langle \cos ({\rm{\Theta }}(\tau )-{\rm{\Theta }}(0))\rangle$. Integrating equation (3) and using the real part ${\mathfrak{R}}$ of the Euler identity ${{\rm{e}}}^{{\rm{i}}{x}}=\cos (x)+{\rm{i}}\sin (x)$ yields:

$$\begin{eqnarray}&&g(\tau )={\mathfrak{R}}\,\left\langle {{\rm{e}}}^{{\rm{i}}{\int }_{0}^{\tau }\sqrt{2{D}_{s}}{{\rm{d}}{W}}_{s}}\right\rangle .\end{eqnarray} \tag{ B25 }$$

The term in the real part operator can be interpreted as the characteristic function of the random variable $X(\tau )={\int }_{0}^{\tau }\sqrt{2{D}_{s}}{{\rm{d}}{W}}_{s}$ for wavenumber k = 1. Using the moment representation of the characteristic function, we obtain

$$\begin{eqnarray}&&g(\tau )={\mathfrak{R}}\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{{i}^{n}}{n!}{m}_{n}(\tau ),\end{eqnarray} \tag{ B26 }$$

where we have defined the moments ${m}_{n}=\langle {X}^{n}\rangle$. For symmetry reasons, the odd moments vanish

$$\begin{eqnarray}&&{m}_{2n+1}=0,\end{eqnarray} \tag{ B27 }$$

and the real part operator can be skipped. First, we calculate m₂, where we use again the independence of the two stochastic processes dW_t and D_t in the second line

$$\begin{eqnarray}\begin{array}{rcl}{m}_{2}(\tau ) & = & \left\langle {\displaystyle \int }_{0}^{\tau }\sqrt{2{D}_{{s}_{1}}}{{\rm{d}}{W}}_{{s}_{1}}{\displaystyle \int }_{0}^{\tau }\sqrt{2{D}_{{s}_{2}}}{{\rm{d}}{W}}_{{s}_{2}}\right\rangle \\ & = & {\displaystyle \int }_{0}^{\tau }{\displaystyle \int }_{0}^{\tau }\langle \sqrt{2{D}_{{s}_{1}}}\sqrt{2{D}_{{s}_{2}}}\rangle \langle {{\rm{d}}{W}}_{{s}_{1}}{{\rm{d}}{W}}_{{s}_{2}}\rangle \\ & = & {\displaystyle \int }_{0}^{\tau }{\displaystyle \int }_{0}^{\tau }\langle \sqrt{2{D}_{{s}_{1}}}\sqrt{2{D}_{{s}_{2}}}\rangle \delta ({s}_{1}-{s}_{2}){{\rm{d}}{s}}_{1}{{\rm{d}}{s}}_{2}\\ & = & 2{\displaystyle \int }_{0}^{\tau }\langle D\rangle {{\rm{d}}{s}}_{1}\\ & = & 2\langle D\rangle \tau .\end{array}\end{eqnarray} \tag{ B28 }$$

Next, we calculate the fourth moment m₄, where we use the correlation function of equation (B19) in the fourth line:

$$\begin{eqnarray}\begin{array}{rcl}{m}_{4}(\tau ) & = & {\displaystyle \int }_{0}^{\tau }{\displaystyle \int }_{0}^{\tau }{\displaystyle \int }_{0}^{\tau }{\displaystyle \int }_{0}^{\tau }\langle \sqrt{2{D}_{{s}_{1}}}\sqrt{2{D}_{{s}_{2}}}\sqrt{2{D}_{{s}_{3}}}\sqrt{2{D}_{{s}_{4}}}\rangle \langle {{\rm{d}}{W}}_{{s}_{1}}{{\rm{d}}{W}}_{{s}_{2}}{{\rm{d}}{W}}_{{s}_{3}}{{\rm{d}}{W}}_{{s}_{4}}\rangle \\ & = & {\displaystyle \int }_{0}^{\tau }{\displaystyle \int }_{0}^{\tau }{\displaystyle \int }_{0}^{\tau }{\displaystyle \int }_{0}^{\tau }\langle \sqrt{2{D}_{{s}_{1}}}\sqrt{2{D}_{{s}_{2}}}\sqrt{2{D}_{{s}_{3}}}\sqrt{2{D}_{{s}_{4}}}\rangle 3\delta ({s}_{1}-{s}_{2})\delta ({s}_{3}-{s}_{4}){{\rm{d}}{s}}_{1}{{\rm{d}}{s}}_{2}{{\rm{d}}{s}}_{3}{{\rm{d}}{s}}_{4}\\ & = & 12{\displaystyle \int }_{0}^{\tau }{\displaystyle \int }_{0}^{\tau }\langle {D}_{{s}_{1}}{D}_{{s}_{2}}\rangle {{\rm{d}}{s}}_{1}{{\rm{d}}{s}}_{2}\\ & = & 12{\displaystyle \int }_{0}^{t}{\displaystyle \int }_{0}^{t}\langle D{\rangle }^{2}+\displaystyle \frac{{({D}_{0}-{D}_{T})}^{2}r\lambda }{{(\lambda +r)}^{2}}{{\rm{e}}}^{-(\lambda +r)| {s}_{1}-{s}_{2}| }{{\rm{d}}{s}}_{1}{{\rm{d}}{s}}_{2}\\ & = & 12(\langle D{\rangle }^{2}{\tau }^{2}+I)\\ & = & 12\left(\langle D{\rangle }^{2}{\tau }^{2}+2\displaystyle \frac{{{\rm{\Delta }}}^{2}D}{\lambda +r}\left[\tau +\displaystyle \frac{{{\rm{e}}}^{-(\lambda +r)(\tau -0)}}{\lambda +r}-\displaystyle \frac{1}{\lambda +r}\right]\right).\end{array}\end{eqnarray} \tag{ B29 }$$

Replacing in the double integral ${\int }_{0}^{\tau }\ldots {{\rm{d}}{s}}_{1}$ by $2{\int }_{0}^{{s}_{2}}\ldots {{\rm{d}}{s}}_{1}$, the integral I is calculated as follows:

$$\begin{eqnarray}\begin{array}{rcl}I & = & 2{\displaystyle \int }_{0}^{\tau }{\displaystyle \int }_{0}^{{s}_{2}}{{\rm{\Delta }}}^{2}{{D}{\rm{e}}}^{-(\lambda +r)({s}_{2}-{s}_{1})}{{\rm{d}}{s}}_{1}{{\rm{d}}{s}}_{2}\\ & = & 2{\displaystyle \int }_{0}^{\tau }{\left[\displaystyle \frac{{{\rm{\Delta }}}^{2}D}{\lambda +r}{{\rm{e}}}^{-(\lambda +r)({s}_{2}-{s}_{1})}\right|}_{{s}_{1}=0}^{{s}_{1}={s}_{2}}{{\rm{d}}{s}}_{2}\\ & = & 2\displaystyle \frac{{{\rm{\Delta }}}^{2}D}{\lambda +r}{\displaystyle \int }_{0}^{\tau }1-{{\rm{e}}}^{-(\lambda +r){s}_{2}}{{\rm{d}}{s}}_{2}\\ & = & 2\displaystyle \frac{{{\rm{\Delta }}}^{2}D}{\lambda +r}{\left[{s}_{2}+\displaystyle \frac{{{\rm{e}}}^{-(\lambda +r){s}_{2}}}{\lambda +r}\right|}_{{s}_{2}=0}^{{s}_{2}=\tau }\\ & = & 2\displaystyle \frac{{{\rm{\Delta }}}^{2}D}{\lambda +r}\left[\tau +\displaystyle \frac{{{\rm{e}}}^{-(\lambda +r)\tau }}{\lambda +r}-\displaystyle \frac{1}{\lambda +r}\right].\end{array}\end{eqnarray} \tag{ B30 }$$

Truncating the sum of equation (B26) for n > 4 , we finally obtain:

$$\begin{eqnarray}\begin{array}{rcl}{g}_{{\rm{\Theta }}}(\tau ) & = & 1-\langle D\rangle \tau +\langle D{\rangle }^{2}{\tau }^{2}/2\\ & & +\displaystyle \frac{{{\rm{\Delta }}}^{2}D}{\lambda +r}\left(\tau +\displaystyle \frac{{{\rm{e}}}^{-(\lambda +r)\tau }}{\lambda +r}-\displaystyle \frac{1}{\lambda +r}\right).\end{array}\end{eqnarray} \tag{ B31 }$$

This form suggests a slope $-\langle D\rangle$ of the correlation function for times $\tau \lt {(\lambda +r)}^{-1}$, which in our case means τ < 0.2 s and is just valid for the very initial time range of the correlation function. From equation (B31) we can extract another linear approximation by concentrating on the time range ${(\lambda +r)}^{-1}\lt \tau \lt \langle D{\rangle }^{-1}$. It gives equation (11) from the main text

$$\begin{eqnarray}&&{g}_{{\rm{\Theta }}}(\tau )=1-\left(\langle D\rangle -\displaystyle \frac{{{\rm{\Delta }}}^{2}D}{\lambda +r}\right)\tau \,,\end{eqnarray} \tag{ B32 }$$

from which we obtain an expression for the relaxation rate α_Θ measured in experiments. It is determined by $\langle D\rangle$ and the second term in the brackets is a correction. But it is sufficient to determine separate values for r and λ, when r / λ is known from the analysis of the pdf $p(| {\rm{\Delta }}{\rm{\Theta }}| )$. Solving the equation for α_Θ for either λ or r, we obtain the formulas

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{{{\rm{\Delta }}}^{2}D}{(1+{(\lambda /r)}^{-1})(\langle D\rangle -{\alpha }_{{\rm{\Theta }}})}\,,\end{eqnarray} \tag{ B33 }$$

$$\begin{eqnarray}&&r=\displaystyle \frac{{{\rm{\Delta }}}^{2}D}{(1+\lambda /r)(\langle D\rangle -{\alpha }_{{\rm{\Theta }}})}\,.\end{eqnarray} \tag{ B34 }$$

The directional autocorrelation function ${g}_{{\rm{\Theta }}}(\tau )=\langle {\bf{e}}(\tau )\cdot {\bf{e}}(0)\rangle$ has an exponential form in experiments up to ca. 2s (see figure 4). Here, we validate this dependence by numerically solving equation (3) with the inferred parameters of table 1. The semi-logarithmic plot in figure B1(a) shows the resulting autocorrelation function (blue data points). It is in good agreement with the exponential decay of equation (10) using the relacation rate α_Θ from equation (11), which we derived in the previous section in equation (B32). This validates our proposition for the relaxation rate.

Zoom In Zoom Out Reset image size

Moreover, we can further validate the exponential fit to the experimental directional autocorrelation function using the theoretical value for the relaxation rate. After having inferred the reduced parameter set $(\lambda /r,{D}_{0},{D}_{T})$ as described in the main text using the pdf $p(| {\rm{\Theta }}| )$, we determine the directional autocorrelation function by simulating the angle process with the reduced parameter set for different values of the parameter λ. Figure B1(b) shows the mean squared error Σ of the simulated autocorrelation function compared to the experimental function plotted versus the tumble rate λ. The best match is for a λ very close to the value shown in table 1, which was determined using the theoretical prediction of equation (11) for the relaxation rate α_Θ.

We apply our method to experimental data where E.coli swims in a uniform concentration c₀ of the chemoattractant L-serine. Cells were harvested and washed two times in the motility buffer enriched with the chemoattractant L-serine (Sigma-Aldrich, USA) at final concentrations c₀ = 1, 10, 100, and 1000 μM. Figure H1(a) shows the normalized tumble rate λ / λ* and figure H1(b) shows the normalized swimming speed ${v}_{0}/{v}_{0}^{* }$ both as a function of the serine concentration. Here, λ* and ${v}_{0}^{* }$ denote the reference values from the control experiment without any chemoattractant. First of all, the parameters determined with the help of the tumble recognizer and our inference method agree well within the error bars. Our results for the tumble rate are in qualitative agreement with previous measurements, see figure H1(a). Except for the data point at c₀ = 100 μM our data shows a decreasing tumble rate in higher concentrations of serine as in [3, 48]. Currently, we have no explanation for the outlier at c₀ = 100 μM, which needs further investigation. Within error bars, our results for the swimming speed shown in figure H1(b) are similar to those obtained for the E.coli RP437 strain [48, 49], in particular we do not see the 40% increase in swimming speed of E.coli AW405 in serine reported by [3], neither. We emphasize that both analysis methods, heuristic tumble recognizer and our inference method, do not show this increase. Thus, we conclude that our method is also applicable for experiments in serine. The lack of the increase in swimming speed is an experimental issue and goes beyond the scope of this paper. A potential reason could be the different motility buffer used in our experiments compared to the one in [3].

Zoom In Zoom Out Reset image size

We apply our method to experimental data from our previous work [22]. A similar experimental set-up to the gradient experiment from the main text was used except for the growth medium, which was Lysogeny broth (LB) instead of TB. Our method infers a tumble rate of λ = 0.37 ± 0.02 s⁻¹ and a swimming speed of ${v}_{0}=18.1\pm 0.2\,\mu {\rm{m}}\,{{\rm{s}}}^{-1}$. This is in agreement with the tumble rate determined in [22] (λ = 0.39 ± 0.03 s⁻¹) and the swimming speed, which we obtained using the heuristic tumble recognizer (${v}_{0}=18.0\pm 0.1\,\mu {{\rm{m}}{\rm{s}}}^{-1}$). Thus, our method also reproduces essential parameters in a different growth medium, where bacteria tumble less often.