How the Motility Pattern of Bacteria Affects Their Dispersal and Chemotaxis

Most bacteria at certain stages of their life cycle are able to move actively; they can swim in a liquid or crawl on various surfaces. A typical path of the moving cell often resembles the trajectory of a random walk. However, bacteria are capable of modifying their apparently random motion in response to changing environmental conditions. As a result, bacteria can migrate towards the source of nutrients or away from harmful chemicals. Surprisingly, many bacterial species that were studied have several distinct motility patterns, which can be theoretically modeled by a unifying random walk approach. We use this approach to quantify the process of cell dispersal in a homogeneous environment and show how the bacterial drift velocity towards the source of attracting chemicals is affected by the motility pattern of the bacteria. Our results open up the possibility of accessing additional information about the intrinsic response of the cells using macroscopic observations of bacteria moving in inhomogeneous environments.

For the case of the power law distributed run times from Eq. (17), we proceed with the asymptotic analysis, assuming t 1 , t 2 τ 0 , meaning s, p → 0. The leading terms in this expansion after the inverse Laplace transform yield Eq. (18).
In the following, we derive our main result for the velocity autocorrelation function, MSD, and diffusion coefficient of the run-reverse-flick pattern. As defined in Eq. (19), we denote the probability density function of turning angles at reversals and flicks by g r (ϕ) and g f (ϕ), respectively. The probability densities of run times after reversals and flicks are denoted as f r (τ ) and f f (τ ), respectively. The corresponding survival probabilities are F r,f (τ ) = 1 − τ 0 dt f r,f (t). At this stage, we consider a CTRW in the angle ϕ(t) without rotational diffusion during the run events; due to Eq. (6), the latter can be included in the velocity autocorrelation function afterwards. To formulate the CTRW, we introduce the auxiliary probability densities ν r,f (ϕ, t), with ν r,f (ϕ, t)dϕ dt denoting the probability that a reversal (r) or flick (f ) occurs during the time interval [t, t + dt], while the cell moves along the direction [ϕ, ϕ + dϕ]. We also refer to ν r,f (ϕ, t) as the frequency of reversals and flicks. Similar to the standard CTRW model, these densities satisfy the integral equations (see Ref. [33] for example) where N 0 r,f (ϕ, τ ) is the initial condition for the probability density that a cell has direction ϕ at time t 1 after the reversal (flick) occured time τ ago. The probability densities P r,f (ϕ, t) for a cell to be in the run mode after a reversal or flick with direction ϕ at time t is then expressed as and P (ϕ, t) = P r (ϕ, t) + P f (ϕ, t) gives the total probability density that a cell has direction ϕ at time t.
To complete the set of CTRW equations, we couple the densities P r,f (ϕ, t) to the frequencies ν r,f (ϕ, t) via [33] The initial condition N 0 r,f (ϕ, τ ) is of most general form. For our purpose of calculating P (ϕ 1 , t 1 ; ∆ϕ, t 2 ), it will be given as a result of the previous evolution during the time t 1 . In that case, also t 2 = t 1 + t and ϕ = ∆ϕ. We introduce the notation and after a Fourier-Laplace transform and some algebraic manipulations, we obtain For the densities P r,f (k, s), we find To proceed, we have to take into account the correct initial conditions (previous evolution during t 1 ). We restrict to the simpler and experimentally relevant case of exponential run time distributions, Now, we have to find the time evolution during t 1 . Although we can implement an arbitrary initial condition at t = −t 1 , we will proceed in a different way. It is easy to show that is the steady state density and for t ∈ [−t 1 , 0], P f,r (k 1 = 0, t) = P 0, eq f,r = const. Therefore, we use the following, now intermediate conditions,ν 0 r (s) = P 0, eq ff f (s), ν 0 f (s) = P 0, eq rfr (s), (S11b) and insert them into Eq. (S8). This gives . (S12) Finally, we consider explictly the turning angle distributions,ḡ r,f (k 2 ) = cos (k 2 ϕ r,f ), which givesḡ r (k 1+sτ r,f , and we finally arrive at . (S13) If we combine Eq. (S13) with Eqs. (8,9), we recover the velocity autocorrelation function from Eq. (20).
The MSD displays linear scaling for large times with the diffusion coefficient After including rotational diffusion as well, the diffusion coefficient becomes the expression from Eq. (24).

II. Estimate of the strength of the chemotactic response
Here, we present an estimate for the strength of the chemotactic response W from Eq. (26). We refer to recent measurements by Vuppula et al.

III. Alternative derivation of the velocity autocorrelation function
We now present a different derivation for the velocity autocorrelation function for run-reverse-flick motion, which is valid for two and three spatial dimensions. In this context, we consider a more general random walk with two alternating tumbling events, as depicted in Fig. S1. The angular changes during tumbling are specified by the preferential turning angles ∆ϕ 1 and ∆ϕ 2 and expressed in terms of the persistence parameters α = cos ∆ϕ 1 and β = cos ∆ϕ 2 . It turns out that only the mean cosines of the turning angles, α and β, enter the directional correlation function and the corresponding diffusion coefficient. Note that different turning angle distributions can yield the same persistence parameter. However, in the exact limiting cases of α, β = ±1, only delta-peaked distributions at angles 0 • and 180 • are possible. Figure S1. Pattern of the random walk with two alternating tumbling events. The persistence parameters α = cos ∆ϕ 1 and β = cos ∆ϕ 2 correspond to the two preferential turning angles ∆ϕ 1 and ∆ϕ 2 , respectively. The walker moves with constant speed v along the direction e(t), and the run times are exponentially distributed with mean τ run .
Our goal is to determine the directional part of the velocity autocorrelation function, C rw (t, t ) = v 2 e(t) · e(t ) rw . We first neglect rotational diffusion during the run events and assume an exponential run time distribution f (τ ) = λ exp (−λτ ), where λ is the tumbling rate and the mean run time is τ run = λ −1 . The survival probability that a run that starts at t = 0 is not interrupted before time t is given by . For a random walk whose first turning event is specified by α, the directional correlation function is given by The interpretation of the terms on the right-hand side of Eq. (S16) is as follows. A run may start at t = 0 and is not interrupted till t, such that the direction of motion is constant; a run may also take place from time 0 . . . t 1 , then tumbling of type α occurs, and a second run without tumbling follows during t 1 . . . t; the next term considers the possibility of a run (0 . . . t 1 ), tumbling of type α at t 1 , another run (t 1 . . . t 2 ), tumbling of type β at t 2 , and one more run up to time t, etc. The emergence of products α Nα β N β in the higher-order terms of Eq. (S16) stems from the occurrence of N α and N β tumbling events of type α and β, respectively. It is important to note that Eq. (S16) is valid in two and three dimensions.
To calculate the series in Eq. (S16), we observe that for the exponential distribution Eq. (S16) can be written as where the time-dependent terms e −λt (λt) k k! coincide with the Poisson distribution. Using series expansions for hyperbolic functions results in the expression e(t) · e(0) rw = e −λt α β sinh αβ λt + cosh αβ λt .
Eq. (S19) is not symmetric with respect to interchanging α ↔ β. This reflects the fact that, so far, each random walker interrupts the first run with a tumbling event of type α. However, we are interested in modeling an ensemble of walkers where the first tumbling event is either of type α or β with the same probability. To generalize Eq. (S19) to an ensemble average for arbitrary times, we thus take into account that, when choosing a walker of the ensemble, the probability is 1/2 for the next tumbling event to be α or β. Formally, we symmetrize Eq. (S19) with respect to α and β, and obtain Next, we generalize the result for e(t) · e(t ) if we also take into account rotational diffusion during the runs. Rotational diffusion in three dimensions with rotational diffusion constant D r leads to e(t) · e(t ) rot = exp (−2D r |t − t |). With rotational diffusion, Eq. (S16) is modified according to Compared to Eq. (S16), Eq. (S21) contains an additional factor of exp(−2D r t). Note that this statement also holds true if the run times are not exponentially distributed; it corresponds to the factorization according to Eq. (7). The MSD is given by from which we deduce the diffusion coefficient Setting α = −1 and β = 0 yields the diffusion coefficient for run-reverse-flick motion from Eq. (23); setting α = β provides Eq. (14). Recall that this derivation, based on Eq. (S21), is independent of the spatial dimension. In addition, it does not depend on the details of the distribution of turning angles ∆ϕ 1 and ∆ϕ 2 , but only on the persistence parameters α and β.

IV. Calculation of the chemotactic drift speed for run-tumbleflick motion
Here, we calculate the chemotactic drift speed v d for the motility pattern run-tumble-flick. This is the random walk where tumbling events, specified by persistence parameter α, alternate with flicks, where a new direction of motion is randomly chosen. In the random walk scheme of Fig. S1, this corresponds to setting β = 0. For V. alginolyticus, the reversal between forward and backward runs results in α = −1.
The run times are exponentially distributed, and for simplicity, we assume a single mean run time τ run = λ −1 for forward and backward motion in a homogeneous environment. We consider a small chemical gradient |∇c| in z direction, and the concentration c(t) at the cell's position z(t) becomes c(t) = |∇c| z(t). (S24) Next, we generalize de Gennes' approach from Eq. (25) for the time-dependent rates of forward and backward runs, where the index r denotes the forward run after the flick, and b the backward run after the reversal. Note that we assume the same chemotactic response function after both kinds of turning event; a more general calculation can be found in Ref.
[52]. The adaptive response, ∞ 0 dt R(t) = 0, has allowed us to set a possible constant concentration in Eq. (S24) to zero, as it does not affect the tumbling rates in Eq. (S25). We will use the experimentally measured response function R(t) of E. coli [Eq. (26)], as we are primarily interested in understanding the influence of the random walk pattern on the chemotactic drift speed.
To simplify forthcoming calculations, we make use of the structure of Eq. (S25) and first consider the special case of a delta-response in time R(t) = A δ(t − T ) with delay time T and strength A, as it was done in Refs. [47,48]. We denote the chemotactic drift speed for the delta-response as v δ . Having found this intermediate result for v δ , it is generalized to the chemotactic drift speed v d by integrating v δ with the response function R(t), according to We calculate v δ for a random walk with a flick at t = 0. The flick totally randomizes the direction of the subsequent forward run, which is interrupted at time t r . At t r , the direction of motion is changed with persistence parameter α and the particle performs a backward run of duration t b . As every flick event destroys the directional persistence and memory of the previous random walk steps, v δ is calculated by considering the mean displacement along the z axis during one cycle of the run-tumble-flick process. To obtain v δ in a first-order expansion of |∇c|, we determine the average displacement of a forward run z r and a subsequent backward run z b . As these expressions are first order in |∇c|, the mean duration of the cycle is given by 2τ run , and the chemotactic drift speed for the delta-response becomes The random walker starts with a flick at t = 0 from the initial position z = 0. We need to calculate the displacement along the direction of the gradient, the z axis, during the forward run, t ∈ [0, t r ], and the subsequent backward run, t ∈ [t r , t r + t b ]. The probability density p r (0 → t r ) that a forward run starts at t = 0 and stops at t r is given by and p b (t r → t r + t b ) is the probability density that a backward run starts at t r and stops at t r + t b . The time-dependent rates λ i (t) (i = r, b) due to chemotaxis are given in Eq. (S25). Taking a delta-response R(t) = A i δ(t − T i ) with strength A i and delay time T i , and using Eq. (S24), we can write the rates as and they thus depend on the previous positions z(t − T i ).
The mean displacement z of the forward run z r and subsequent backward run z b reads The averages · · · in Eq. (S30) are taken over all possible paths due to random swimming directions after the flick event. In addition, they represent ensemble averages to account for the rotational diffusion during the runs. Without rotational diffusion, we could write the positions in Eq. (S30b) as z( To proceed, we shift the averages · · · into the integrals in Eq. (S30b). We start to calculate z r : In the first step, a partial integration led us to Eq. (S31b). Then we introduced v z (t r ) = dz(tr) dtr as the speed in z direction, and substituted Eq. (S29) into (S31b); after expanding the exponential, we only kept the first-order terms in |∇c| to arrive at Eq. (S31c). The averages in this equation will be calculated without chemotaxis to remain first-order in |∇c|. We immediately find v z (t r ) = 0, as the direction of motion after the flick at t = 0 is random. Further, we use z(t) = t 0 ds v z (s) and get Without chemotaxis, the velocity distribution is isotropic, which implies and the directional correlation function reads e(s) · e(t r ) = e −2Dr(tr−s) 0 ≤ s < t r , 0 , else.
After inserting Eqs. (S33) and (S34) into Eq. (S32), we perform the resulting integrals and obtain the mean displacement of a forward run, which is proportional to the gradient |∇c| and the response strength A r . To check for consistency, we note that Eq. (S35) agrees with Locsei's intermediate result for α = 0 from Ref. [47]. Now, we turn to the mean displacement of the backward run z b from Eq. (S30). We proceed in the same way as before and obtain up to first order in |∇c|: Again, the expectation values are calculated without chemotaxis, such that v z (t r ) = 0 and Eq. (S33) holds. Be aware that v z (t b ) is the speed in z direction at the absolute time t r + t b . The directional correlation function becomes e(s) · e(t b ) =      e −2Drt b α e −2Dr(tr−s) , 0 ≤ s < t r , e −2Dr(tr+t b −s) , t r ≤ s < t r + t b , 0 , else. (S37) To obtain the first line of Eq. (S37), we made the decomposition e(s) · e(t b ) = e(t b ) · e(t + r ) e(t + r ) · e(t − r ) e(t − r ) · e(s) , see also Ref. [47]. We obtain Note that the dependence of z b on A r reflects a chemotactic coupling between forward and backward runs. Finally, we consider equal strength of the chemotactic response A r = A b = A and T r = T b = T , and combine all preceding results according to Eq. (S27). In first-order in |∇c|, we thus arrive at the chemotactic drift speed for the delta-response, v δ = A|∇c| v 2 λe −(λ+2Dr)T λ[2 + α(1 + λT )] + 2D r [2 + α(λT − 1)] 6(λ + 2D r ) 3 .
The final drift speed v d from Eq. (28) is obtained by integrating v δ with the response function from Eq. (26) according to Eq. (S26).