Anomalous diffusion in run-and-tumble motion

TL;DR

This paper models bacterial run-and-tumble motion using a combined Brownian and Lévy walk scheme, analyzing how different dwelling time distributions lead to various diffusion behaviors.

Contribution

It introduces a continuous-time random walk model capturing diverse diffusion regimes in bacterial motion based on dwelling time distributions.

Findings

Normal diffusion, superdiffusion, and ballistic spreading can result from different dwelling time distributions.

The model explains how microscopic phase durations influence macroscopic diffusion behavior.

Long-time and short-time behaviors of mean squared displacement are characterized.

Abstract

A random walk scheme, consisting of alternating phases of regular Brownian motion and L\'evy walks, is proposed as a model for run-and-tumble bacterial motion. Within the continuous-time random walk approach we obtain the long-time and short-time behavior of the mean squared displacement of the walker as depending on the properties of dwelling time distribution in each phase. Depending on these distributions, normal diffusion, superdiffusion and ballistic spreading may arise.

Equations

$${\lambda }_1(\mathbf{x}|t)=\frac{1}{{\left(4\pi Dt\right)}^{\frac{d}{2}}}\exp \left(-\frac{{\mathbf{x}}^2}{4Dt}\right).$$

$${\lambda }_2(\mathbf{x}|t)=\frac{1}{S_{d-1}(ct{)}^{d-1}}\delta \left(|\mathbf{x}|-ct\right).$$

$$⟨{\mathbf{x}}^2(t)⟩=\int {\mathbf{x}}^{2}G(\mathbf{x},t|0,0)d\mathbf{x}.$$

$${\Lambda }_i(\mathbf{x},t)={\lambda }_i(\mathbf{x}|t){\int }_t^{\infty }{\psi }_i(t^{\prime })d{t}^{\prime }.$$

$$G_1$$

$$=\sum _{m=0}^{\infty }[{\sigma }_1(\mathbf{k},s){\sigma }_2(\mathbf{k},s){]}^m[{\Lambda }_1(\mathbf{k},s)+{\sigma }_1(\mathbf{k},s){\Lambda }_2(\mathbf{k},s)],$$

$$G_1(\mathbf{k},s)=\frac{{\Lambda }_1(\mathbf{k},s)+{\sigma }_1(\mathbf{k},s){\Lambda }_2(\mathbf{k},s)}{1-{\sigma }_1(\mathbf{k},s){\sigma }_2(\mathbf{k},s)}.$$

$$G(\mathbf{k},s)=\frac{{\Lambda }_1+{\Lambda }_2+{\sigma }_1{\Lambda }_2+{\sigma }_2{\Lambda }_1}{2(1-{\sigma }_1{\sigma }_2)},$$

$$⟨{\mathbf{x}}^2(s)⟩=-{\Delta }_{\mathbf{k}}{G(\mathbf{k},s)|}_{\mathbf{k}=0}$$

$${\sigma }_1(\mathbf{k},s)$$

$$={\int }_0^{\infty }dt{e}^{-st-D{\mathbf{k}}^2}{\psi }_1(t)={\psi }_1(s+D{\mathbf{k}}^2)$$

$${\Lambda }_1(\mathbf{k},s)=\frac{1-{\psi }_1(s+D{\mathbf{k}}^2)}{s+D{\mathbf{k}}^2}.$$

$${\sigma }_2(\mathbf{k},s)=$$

$$={\int }_{{\Omega }_{d-1}}d\mathbf{e}{\int }_0^{\infty }dr{r}^{d-1}{\int }_0^{\infty }dt\frac{\delta (r-ct){\psi }_2(t)}{S_{d-1}(ct{)}^{d-1}}{e}^{-st+ir\mathbf{ke}}.$$

$${\sigma }_2(\mathbf{k},s)=\frac{1}{S_{d-1}}{\int }_{{\Omega }_{d-1}}d\mathbf{e}{\psi }_2(s-ic\mathbf{ke})$$

$${\Lambda }_2(\mathbf{k},s)=\frac{1}{S_{d-1}}{\int }_{{\Omega }_{d-1}}d\mathbf{e}\frac{1-{\psi }_2(s-ic\mathbf{ke})}{s-ic\mathbf{ke}}.$$

$$⟨{\mathbf{x}}^2(s)⟩$$

$$+\frac{c^2(1+{\psi }_1)}{s\left(1-{\psi }_1{\psi }_2\right)}\left(\frac{{\psi }_2^{\prime }}{s}+\frac{1-{\psi }_2}{s^2}\right),$$

$$⟨{\mathbf{x}}^2(s)⟩=$$

$$=\frac{dD{\tau }_1(s{\tau }_2+2)}{s^2(s{\tau }_1{\tau }_2+{\tau }_1+{\tau }_2)}+\frac{c^2{\tau }_2^2(s{\tau }_1+2)}{s^2(s{\tau }_2+1)(s{\tau }_1{\tau }_2+{\tau }_1+{\tau }_2)},$$

$$⟨{\mathbf{x}}^2(t)⟩\simeq dDt$$

$$⟨{\mathbf{x}}^2(t)⟩\simeq 2\frac{dD{\tau }_1+c^2{\tau }_2^2}{{\tau }_1+{\tau }_2}t$$

$${\psi }_1(s)$$

$${\psi }_2(s)$$

$$⟨{\mathbf{x}}^2(s)⟩\approx \frac{2A_2{c}^2\left(1-\beta \right)s^{\beta -2}}{A_1{s}^{\alpha +1}+A_2{s}^{\beta +1}}.$$

$$\langle\mathbf{x}^{2}(t)\rangle\simeq c^{2}(1-\beta)\left\{\begin{array}[]{l l}t^{2}&\mbox{for}\;\beta<\alpha\\ \frac{A_{2}}{A_{1}+A_{2}}t^{2}&\mbox{for}\;\beta=\alpha\\ \frac{2}{\Gamma(3+\alpha-\beta)}\frac{A_{2}}{A_{1}}t^{2-(\beta-\alpha)}&\mbox{for}\;\beta>\alpha\end{array}\right.$$

$${\psi }_2=1-A_{2}s+B_2{s}^{\beta }+\cdots \text{.}$$

$$⟨{\mathbf{x}}^2(s)⟩\approx 2\frac{dD{A}_1{s}^{\alpha -1}+B_2{c}^2\beta s^{\beta -2}}{A_1{s}^{\alpha +1}+A_2{s}^2}$$

$$\langle\mathbf{x}^{2}(t)\rangle\simeq\left\{\begin{array}[]{l l}\frac{2(\beta-1)c^{2}}{\Gamma(3+\alpha-\beta)}\frac{B_{2}}{A_{1}}t^{2-(\beta-\alpha)}&\mbox{for}\;\beta-1<\alpha<1\\ 2\left(dD+c^{2}(\beta-1)\frac{B_{2}}{A_{1}}\right)t&\mbox{for}\;\beta-1=\alpha<1\\ 2dDt&\mbox{for}\;\alpha<\beta-1<1\end{array}\right.$$

Full text

We propose a qualitative model for the movement of motile bacteria like E. Coli. These organisms have several flagella, by which they are propelled. The motion of E. Coli can be subdivided in two different phases: tumble- and run-mode. In the tumbling-mode the flagella rotate in different directions, causing an erratic movement and a fast change of orientation. In the run-mode all flagella rotate in the same direction and entangle, resulting in a straight forward movement with a constant velocity. The bacteria use the tumbling phase for chemotaxis, swimming in the direction of nutrients or away from poisonous areas. The whole process is described for instance in Berg2004 and its numerical and theoretical description is a subject of recent discussion, see e.g. Othmer1988 ; Visser2003 ; NicolauJr2009 .

Random walks are the model of choice in description of motion of living organisms Berg2004 ; Othmer1988 ; NicolauJr2009 ; Meysman2008 . Run and tumble motion is often described either by a Lévy-walk or a continuous time random walk (CTRW), also called “velocity” and “jumping process” Othmer1988 . CTRWs describe a random process consisting in instantaneous jumps separated by periods of rest (waiting times) following a given probability distribution as first discussed in Montroll1965 . A Lévy-walk is a motion at a constant speed where the direction of velocity is changed after random epochs Shlesinger1987 . It is essentially a CTRW in which the waiting period corresponds to a motion with constant velocity and jumps lead to a change of the direction of motion. Although Lévy-walks are commonly used to describe the (random) motion of bacteria, the origin of this kind of motion is a topic of recent studies (see e.g. Matthaeus2009 ).

In this paper we combine both processes and model the run-and-tumble motion as an alternating random walk consisting of phases of simple diffusion (tumble phase) and of Lévy walk (run phase). A sketch of the overall process is shown in Figure 1. The speed of motion in the run phase will be denoted by $c$ and the diffusion coefficient in the tumbling phase by $D$. The dwelling time $t$ in each phase is given by waiting time probability density functions (PDFs) $\psi_{1}(t)$ and $\psi_{2}(t)$ for the tumbling and run phases, respectively. All other properties in the run and tumbling phases will be denoted by subscripts $1$ and $2$ in what follows. The whole motion takes place in the space of $d$ dimensions, with $d=2$ or $3$. For the sake of simplicity, we assume that the observation starts at a beginning of a run or of a tumble phase (in order to avoid the discussion of aging effects). Although the alternating setup is not a new process in the theory of random walks, and the general discussion follows the standard lines (see Klafter2011 ), the model discussed below is new and far from being trivial.

The displacements $\mathbf{x}$ in a tumbling phase are given by the transition probability density of Brownian motion

[TABLE]

It depends on the time $t$ the walker has spent in this phase. Respectively, if $t$ is the elapsed time in the run phase, the transition probability density for a walk over distance $x$ reads:

[TABLE]

Here $S_{d-1}$ is the (hyper)surface of a $d$-dimensional unit ball which is given by: $S_{d-1}=(2\pi^{d/2})/\Gamma(d/2)$, with $\Gamma(x)$ being the Gamma function. The corresponding probability densities are isotropic as long as no chemotactic effects are included.

We discuss the properties of the propagator $G(\mathbf{x}_{1},t_{1}|\mathbf{x}_{0},t_{0})$ of the complete process starting at $t=0$ and stopped (observed) at time $t$, and concentrate on the mean square displacement (MSD)

[TABLE]

Typically, the MSD follows a power law pattern, $\langle\mathbf{x}^{2}(t)\rangle\propto t^{\gamma}$; the value $\gamma=1$ corresponds to normal diffusion and cases with $\gamma>1$ to superdiffusion. The special limit $\gamma=2$ is referred to as ballistic spreading. Subdiffusion, with $0<\gamma<1$, is impossible in our setting. Since our situation is homogeneous both in space and in time, it is easier to use the Fourier representation for the coordinate and the Laplace representation for time in which $G$ is given by an algebraic expression.

Typically, the process is stopped not at a time it changes the phase. Hence, we need two additional transition probability densities, which are

[TABLE]

They describe the displacement in a prematurely finished phase ($i=1,2$). The integral in this expression gives the probability that the dwelling time in phase $i$ exceeds $t$. Its Laplace transform is $(1-\psi_{i}(s))/s$. Thus, the $\Lambda_{i}$ are the complete analogons to the regular CTRWs waiting time PDFs after the last jump. Eq. (4) takes into account the fact that the walker continuously moves during the waiting period.

Let us assume that walks start with the beginning of a tumble phase. To find the propagator, we note that the corresponding process may have completed some number $m\geq 0$ of full tumble and run cycles before the observation at time $t$ took place, and that the time is stopped either during the tumble or during the following run phase of the last, incomplete, cycle. A propagator corresponding to a full cycle is simply a convolution of $\sigma_{1}(\mathbf{x},t)=\lambda_{1}(\mathbf{x}|t)\psi_{1}(t)$ and $\sigma_{2}(\mathbf{x},t)=\lambda_{2}(\mathbf{x}|t)\psi_{2}(t)$. In a Fourier-Laplace domain $\sigma_{i}(\mathbf{x},t)$ transforms to $\sigma_{i}(\mathbf{k},s)$. The Fourier-Laplace representation of the propagator for $m$ completed cycles is given by $[\sigma_{1}(\mathbf{k},s)\sigma_{2}(\mathbf{k},s)]^{m}$. The remaining last cycle can be stopped either in the tumble phase or in the run phase. In the first case, the joint PDF of displacement and time in this last phase will be given by $\Lambda_{1}(\mathbf{k},s)$, or, in the Fourier-Laplace representation, by $\Lambda_{1}(\mathbf{k},s)$. If the incomplete cycle finishes in the run phase, this PDF is given by the convolution of $\sigma_{1}(\mathbf{x},t)$ and $\Lambda_{2}(\mathbf{x},t)$, corresponding to $\sigma_{1}(\mathbf{k},s)\Lambda_{2}(\mathbf{k},s)$ in the Fourier-Laplace representation. Summing up all the possibilities, we get

[TABLE]

yielding for walkers starting with a tumbling phase

[TABLE]

Since the walker could have started in both phases, we average over the two possibilities. The averaged propagator reads

[TABLE]

where we have omitted the arguments of the corresponding functions. Due to its symmetry, this expression is easier to handle and will be used for the further analysis.

The MSD in the Laplace domain is given by

[TABLE]

where $\Delta_{\mathbf{k}}$ denotes the Laplacian in the Fourier space. To relate this MSD to the PDFs of the waiting times and to the transition probabilities, we perform explicitly the Fourier-Laplace transforms of $\sigma_{1,2}$ and $\Lambda_{1,2}$. It gives in case of tumbling:

[TABLE]

where $\psi_{1}(...)$ is the Laplace representation of $\psi_{1}(t)$. The same procedure applied to $\Lambda_{1}(\mathbf{x},t)$ gives

[TABLE]

Let us now turn to the Lévy walk phase. It reads

[TABLE]

Performing integration in $r$ we obtain

[TABLE]

Here $\Omega_{d-1}$ denotes the hypersurface of a $d$-dimensional unit ball, and $\mathbf{e}$ is an unit vector defining the point on this surface. Integration in $\mathbf{e}$ corresponds in $d=2,3$ to the angle-integration or over the solid angle. Analogously,

[TABLE]

Evaluating the Laplacian of the propagator one obtains a lengthy expression involving $\sigma_{1,2}$, $\Lambda_{1,2}$ and their first and second partial derivatives taken at $\mathbf{k}=0$. All items are expressed through $\psi_{1,2}$ and their first and second derivatives with respect to their argument. For the running phase, they have additionally to be integrated over the surface of the unit ball. All first derivatives of $\sigma_{1,2}$ and $\Lambda_{1,2}$ vanish due to spacial symmetry. Terms containing a second derivative of $\psi_{1,2}$ enter expressions for $\sigma_{1,2}$ and for $\Lambda_{1,2}$ with opposite signs and cancel. The compact final result consists of two parts corresponding to the running and the tumbling phase. It reads:

[TABLE]

where the derivative $\psi^{\prime}$ is taken with respect to the Laplace-variable $s$.

Exponential Waiting Time PDFs. We now proceed by making specific assumptions about the waiting time PDFs $\psi_{i}$. We first assume that all waiting time PDFs take exponential forms so that $\psi_{i}(s)=1/(s\tau_{i}+1)$ with $\tau_{i}$ being mean dwelling times in the corresponding phases. By plugging these expressions into (12), we can get an expression for the second moment in the spectral domain

[TABLE]

and from that, the asymptotic scaling in the time domain

[TABLE]

In order to obtain this result, we used Tauberian theorems, which state that the small $s$ limit corresponds to a large $t$ limit in original domain and vice versa. In both limits $t\rightarrow 0$ and $t\rightarrow\infty$, one observes normal diffusion, albeit with different diffusion coefficients.

Power Law PDFs. We will now consider waiting time distributions asymptotically following power laws for $t\rightarrow\infty$: $\psi_{1}(t)\propto t^{-\left(1+\alpha\right)}$ and $\psi_{2}(t)\propto t^{-\left(1+\beta\right)}$. Interesting new features arise if $\alpha\leq 1$ and $\beta\leq 2$. If the exponents are equal or larger than unity, the first moments of the waiting times are finite but the second ones diverge. If the exponents are less than unity, i.e. $\alpha,\beta<1$, even the means do not exist. As we proceed to show, the diffusive behavior is governed by the first non-analytic term of $\psi_{1}(s)$, resp. $\psi_{2}(s)$.

Let us first discuss the case when the first moments of both dwelling time distributions diverge. The corresponding expansions of the Laplace-transforms in the limit $s\rightarrow 0$ read

[TABLE]

with $\alpha,\beta<1$. In this case Eq.(12) gives in the leading order

[TABLE]

Remarkably, this asymptotic expression does not explicitly depend on the diffusion coefficient $D$ of the tumbling phase. Transformed back to the time domain, the diffusive scaling depends on the relation between $\alpha$ and $\beta$:

[TABLE]

Hence, we obtain superdiffusion in the long-time limit, which is always ballistic for $\beta\leq\alpha$. It is independent on the properties of the tumbling phase for $\beta<\alpha$, and depends on a prefactors $A_{1}$ and $A_{2}$ if $\alpha=\beta$. The first two regimes are dominated by Lévy walks. The third superdiffusive one with $\beta>\alpha$ is close to a sequence of Lévy walks interrupted by rests. The tumbling periods hardly contribute to the displacement, and, therefore, a dependence on $D$ is still absent.

As one can see, the asymptotic expressions (19) fail for $\beta\geq 1$. In this case, subleading orders in Eqs.(17) must be taken into account. With $\beta>1$, the Laplace transform of the waiting time density in the run phase reads

[TABLE]

$A_{2}$ stands for the finite mean waiting time in this phase. The case $\beta=1$ leads to logarithmic corrections, which do not change the qualitative behaviour. We now consider the case $1<\beta\leq 2$ when the first moment of the dwelling time in the run phase does exist, but the second one diverges except for the limiting value $\beta=2$. Again we let $0<\alpha\leq 1$ for the tumble phase. Repeating the procedure, we get

[TABLE]

wherein an explicit dependence on $D$ appears. The behavior in the time domain depends on relation between $\alpha$ and $\beta$:

[TABLE]

If $\beta-1<\alpha$, the behavior of the walker still corresponds to superdiffusion, and the $D$-dependence is suppressed. This case again is similar to Lévy walks interrupted by rests where the tumbling does not contribute to the displacement. In the two following situations with $\beta-1\geq\alpha$, the walker exhibits a normal diffusive behavior. In the last one with $\beta>1+\alpha$, the MSD is dominated by the random motion in the tumbling phase. The running phase has no effect on the value of the diffusion coefficient. For short times $t\to 0$ the motion is always dominated by normal diffusion.

In the limit $\alpha\to 1$ the PDF of waiting times in the tumbling phase gets a finite first moment. The corresponding scalings are obtained from Eqs. (19) and (22). If $\beta<1$, i.e. if the waiting times in the run phase do not possess neither a first nor a second moment, the motion is ballistic. With growing $\beta>1$ the motion tends to a superdiffusive one. The special limit is when $\alpha=1$ and $\beta=2$ corresponds to the purely diffusive case. Then $\psi_{1}$ does possess the first moment $\tau_{1}$ corresponding to $A_{1}$ in the expansion (17) and $\psi_{2}$ does possess the second moment with value $2\tau_{2}^{2}$ expressed by $2B_{2}$ in the Laplace transform (17). The MSD can be obtained from the expression (21) by taking now into account the contribution of the second term in the denominator. This result coincides in the long time asymptotics with the one for the exponentially distributed waiting times Eq.(15).

We have seen that by choosing proper waiting time distributions of an alternating CTRW process as described above one can describe different propagation regimes in the long-time limit ranging from normal diffusion up to ballistic spreading. In the case of power-law waiting time PDFs normal diffusive behavior can only be achieved if the distribution for the Brownian part has a heavier tail than the one for runs ($\alpha\leq\beta-1\leq 1$), which has to possess a first moment, i.e. the PDF of the Brownian tumbling phase has to have significantly more mass in the tail than its Lévy-counterpart. Such power-law-like behavior of a waiting time distribution is no rarity and occurs in several other (particularly biological) contexts Mejias2010 .

We therefore expect a transition from normal to superdiffusion in cases with $\alpha>\beta-1>0$ or with $\beta<1$. Fig. 2 shows what diffusive regime may be expected for different values of $\alpha$ and $\beta$.

We note that if the observation does not start at time of a phase change, the model with the considered power-law distributed waiting times shows aging effects Klafter2011 . It may result in changes of prefactors in the corresponding asymptotic expressions but won’t lead to different asymptotics, so that the classification given in Fig. 2 still will hold true.

Let us summarize our findings. We proposed a phenomenological model for bacteria performing a run-and-tumble motion in the absence of chemotaxis. We have determined the possible diffusive regimes of such motion. Depending on the distribution of times in the two phases of motion, the model shows transitions from normal diffusive to superdiffusive and to ballistic behavior. Our model may prove to be useful to model situations that exhibit such transitions, such as the one described in Metzner2007 . It also may be suitable to describe non-biological situations like the motion of nanorods, which may also exhibit superdiffusion Campos2009 ; Golestanian2009 .

The authors acknowledge financial support by DFG within IRTG 1740 research and training group project.