Home | About | Journals | Submit | Contact Us | Français |

**|**PLoS Comput Biol**|**v.13(3); 2017 March**|**PMC5358899

Formats

Article sections

Authors

Related links

PLoS Comput Biol. 2017 March; 13(3): e1005429.

Published online 2017 March 6. doi: 10.1371/journal.pcbi.1005429

PMCID: PMC5358899

Oleg A Igoshin, Editor^{}

Rice University, UNITED STATES

The authors have declared that no competing interests exist.

**Conceptualization:**TE SWZ.**Data curation:**JL.**Formal analysis:**JL TE.**Funding acquisition:**JL TE SWZ.**Investigation:**JL TE.**Methodology:**JL TE.**Project administration:**TE.**Software:**JL TE.**Supervision:**TE SWZ.**Validation:**JL.**Visualization:**JL TE SWZ.**Writing – original draft:**JL TE SWZ.**Writing – review & editing:**JL TE SWZ.

* E-mail: ude.elay@tenome.yrreiht

Received 2016 December 13; Accepted 2017 February 28.

Copyright © 2017 Long et al

This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Many organisms navigate gradients by alternating straight motions (runs) with random reorientations (tumbles), transiently suppressing tumbles whenever attractant signal increases. This induces a functional coupling between movement and sensation, since tumbling probability is controlled by the internal state of the organism which, in turn, depends on previous signal levels. Although a negative feedback tends to maintain this internal state close to adapted levels, positive feedback can arise when motion up the gradient reduces tumbling probability, further boosting drift up the gradient. Importantly, such positive feedback can drive large fluctuations in the internal state, complicating analytical approaches. Previous studies focused on what happens when the negative feedback dominates the dynamics. By contrast, we show here that there is a large portion of physiologically-relevant parameter space where the positive feedback can dominate, even when gradients are relatively shallow. We demonstrate how large transients emerge because of non-normal dynamics (non-orthogonal eigenvectors near a stable fixed point) inherent in the positive feedback, and further identify a fundamental nonlinearity that strongly amplifies their effect. Most importantly, this amplification is asymmetric, elongating runs in favorable directions and abbreviating others. The result is a “ratchet-like” gradient climbing behavior with drift speeds that can approach half the maximum run speed of the organism. Our results thus show that the classical drawback of run-and-tumble navigation—wasteful runs in the wrong direction—can be mitigated by exploiting the non-normal dynamics implicit in the run-and-tumble strategy.

Countless bacteria, larvae and even larger organisms (and robots) navigate gradients by alternating periods of straight motion (runs) with random reorientation events (tumbles). Control of the tumble probability is based on previously-encountered signals. A drawback of this run-and-tumble strategy is that occasional runs in the wrong direction are wasteful. Here we show that there is an operating regime within the organism’s internal parameter space where run-and-tumble navigation can be extremely efficient. We characterize how the positive feedback between behavior and sensed signal results in a type of non-equilibrium dynamics, with the organism rapidly tumbling after moving in the wrong direction and extending motion in the right ones. For a distant source, then, the organism can find it fast.

Navigation up a gradient succeeds by finding those directions in which signals of interest increase. This can be difficult when the size of the navigator is small compared to the length scale of the gradient because local directional information becomes unreliable. In this case, cells [1, 2], worms [3], larvae [4], and even robots [5] often adopt a run-and-tumble strategy to navigate. During runs the organism moves approximately straight, collecting differential sensor information in one direction. Tumbles, or reorientations at zero speed, enable the organism to explore other directions. Signal levels are transduced rapidly to the motility apparatus through an internal state variable, so that increases in attractant transiently raise the probability to run longer (and tumble less) before a negative integral feedback adapts it back [6]. Classically, averaging over many runs and tumbles results in a net drift up the gradient, although this is usually rather modest because of the occasional runs in the wrong direction. We here focus on the positive feedback inherent to this strategy wherein motion up the gradient lowers the probability to tumble, which further boosts drift up the gradient. Our analysis reveals an unstudied regime in which rapid progress can be achieved. Small fluctuations in the speed of the organism along the gradient grow into large transients in the correct direction but small ones otherwise. We show that this asymmetric amplification arises from the positive feedback, which causes the eigenvectors near the adapted state of the dynamical system to become non-orthogonal, therefore leading to non-normal dynamics. The resulting large transient are further boosted by a nonlinearity that is intrinsic to the positive feedback. Such non-normal dynamics were first discovered in fluid mechanics where they were shown to play an important role in the onset of turbulence in the absence of unstable modes [7, 8].

Past theoretical studies of run-and-tumble navigation have mostly focused on what happens when adaptation dominates the dynamics (e.g. [9–13]). In this regime, the internal state of the organism exhibits small fluctuations around its mean, and mean field theory (MFT) can be applied to make predictions. This approach has been used to describe the motile behavior or populations of *E. coli* bacteria in exponential ramps [13–15] and oscillating gradients [16]. Beyond the well-understood negative feedback-dominated regime there is a large portion of physiologically relevant parameter space where the positive feedback between movement and sensation dominates the run-and-tumble dynamics. Agent-based simulations have shown that, in this case, large transient fluctuations can emerge in the internal state of an individual organism climbing a gradient, precluding the use of mean field approaches [13]. While systems of partial differential equations (PDEs) can be integrated numerically to reproduce these dynamics [17], a precise understanding of the role of the positive feedback in generating such large fluctuations and the impact of those on the performance of a biased random walk are fundamental questions that remain largely unanswered because of difficulties in obtaining analytical results.

Here we develop an analytical model of run-and-tumble gradient ascent that preserves the rich nonlinearity of the problem and incorporates the internal state, 3D-direction of motion, and position of the organism as stochastic variables. We find that large fluctuations in the internal state originate from two key mechanisms: (i) the non-normal dynamical structure of the positive feedback that enables small fluctuations to grow, and (ii) a quadratic nonlinearity in the speed along the gradient that further amplifies such transients asymmetrically. Utilizing phase space analysis and stochastic simulations, we show how these two effects combine to generate a highly effective “ratchet-like” gradient-climbing mode that strongly mitigates the classic drawback of biased random walks: wasteful runs in wrong directions. In this new regime an organism should be able to achieve drift speeds on the order of the maximum swim speed. Our results are general in that they apply to a large class of biased random walk strategies, where run speed and sampling of new directions may be modulated based on previously encountered signals.

Consider a random walker with an internal state variable *F* that follows linear relaxation towards the adapted state *F*_{0} over the timescale *t*_{M}, which represents the memory duration of the random walker. We assume that the perceived signal, *ϕ*(**X**, *t*) = *ϕ*(*C*(**X**, *t*)), at position **X** and time *t* (here *C* represents the signal), is rapidly transduced to determine the value of an internal state variable *F* via a receptor with gain *N*:

$$\dot{F}=-\frac{F-{F}_{0}}{{t}_{M}}+N\left(\frac{\partial}{\partial t}+\dot{\mathbf{X}}\xb7\nabla \right)\varphi (\mathbf{X},t).$$

(1)

Stochastic switching between runs and tumbles depends on *F* and follows inhomogeneous Poisson statistics with probability to run *r*(*F*) = λ_{T}/(λ_{R} + λ_{T}) = 1/(1 + exp(−*HF*)), where *H* is the gain of the motor, and λ_{R}(*F*) and λ_{T}(*F*) are the transition rates from run to tumble and vice versa [15, 18].

During runs the speed is constant $\parallel \dot{\mathbf{X}}\parallel ={v}_{0}$ and the direction of motion is subject to rotational Brownian motion with diffusion coefficient *D*_{R}. During tumbles the speed is nil and reorientation follows rotational diffusion *D*_{T} > *D*_{R} to account for persistence effects [19]. Taken together, these two processes cause the random walker to lose its original direction at the expected rate ${t}_{D}^{-1}=(n-1)(r{D}_{R}+(1-r{D}_{T})$ where *n* = 2, 3 for two- and three-dimensional motion respectively. Note that, in this minimal model we ignore possible internal signaling noise [20, 21], and all randomness comes from the rotational diffusions *D*_{R} and *D*_{T} as well as the stochastic switchings with rates λ_{R}(*f*) and λ_{T}(*f*). The effect of signaling noise is considered below using agent-based simulations. Since *ϕ*(*C*) can be nonlinear, Eq (1) includes possible effects of saturation of the sensory system.

Consider a static one-dimensional gradient and define the length scale of the perceived gradient as *L*(**X**) = 1/||∇*ϕ*(**X**)|| and the direction of motion as $s=\widehat{\mathbf{u}}\xb7\widehat{\mathbf{X}}$. Then from Eq (1) the internal dynamics satisfies the following equations during runs and tumbles, respectively:

$$\begin{array}{cccc}& \dot{F}{|}_{run}\hfill & \hfill =& -\frac{F-{F}_{0}}{{t}_{M}}+\frac{N{v}_{0}}{L\left(\mathbf{X}\right)}s\hfill \\ & \dot{F}{|}_{tumble}\hfill & \hfill =& -\frac{F-{F}_{0}}{{t}_{M}}.\hfill \end{array}$$

(2)

We are interested in the displacement of the random walker along the gradient over timescales longer than individual runs and tumbles. In the limit where the switching timescale *t*_{S} = 1/(λ_{R} + λ_{T}) is much shorter than the other timescales we derive from a two-state stochastic model and Eq (2) (Methods Eqs (13) and (14)):

$$\dot{r}=r\left(1-r\right)\left(\underset{\text{negative}\phantom{\rule{4.pt}{0ex}}\text{feedback}}{\underbrace{-\frac{f\left(r\right)-{f}_{0}}{{t}_{M}}}}+\underset{\text{positive}\phantom{\rule{4.pt}{0ex}}\text{feedback}}{\underbrace{\frac{r\phantom{\rule{0.166667em}{0ex}}s}{L\left(\mathbf{X}\right)/\left(NH{v}_{0}\right)}}}\right),$$

(3)

where *f* = *HF*. The first term is the negative feedback towards the adapted run probability *r*_{0} = *r*(*f*_{0}). The second term shows how motion up the gradient (*s* > 0) causes the probability to run *r* to feed back on itself—when the organism is oriented up the gradient (*s* > 0), *F* increases only during runs (Eq (2)), and this increase in turn raises *r*(*F*) so that the probability that the dynamics of *F* follows $\dot{F}{|}_{run}$ rather than $\dot{F}{|}_{tumble}$ is increased, and so on. A positive feedback is thereby created with characteristic timescale *t*_{E} = *L*/(*NHv*_{0}). Steeper gradient (smaller *L*), stronger receptor gain *N* or motor gain *H*, or faster speed *v*_{0}, all lead to stronger positive feedback (shorter *t*_{E}). This important timescale, *t*_{E}, together with *t*_{M} (memory duration) and *t*_{D} (direction decorrelation time), effectively determines the dynamics.

Expressing time in units of *t*_{M}, we introduce the following two non-dimensional timescales:

$$\begin{array}{cc}\hfill {\tau}_{E}(\mathbf{X})& =\frac{{t}_{E}(\mathbf{X})}{{t}_{M}}=\frac{L(\mathbf{X})}{{t}_{M}NH{v}_{0}}\hfill \\ \hfill {\tau}_{D}\left(f\right)& =\frac{{t}_{D}}{{t}_{M}}=\frac{1/{t}_{M}}{(n-1)\left(r\left(f\right){D}_{R}+\left(1-r\left(f\right)\right){D}_{T}\right)}.\hfill \end{array}$$

(4)

Here *τ*_{E} quantifies the ratio between the negative and positive feedbacks. (See Table 1 for a summary of the symbols used.) From above, we expect that the dynamics will depend on how *τ*_{E} and *τ*_{D} compare with one.

To explore how run-and-tumble dynamics depend on *τ*_{E} and *τ*_{D}, we used a previously published stochastic agent-based simulator of the bacteria *E. coli* that reproduces well available experimental data on the wild-type laboratory strain RP437 ([15, 22] and S1 Appendix). In this case the internal state *F* represents the free energy of the chemoreceptors. Since *E. coli* approximately detects log-concentrations (S1 Appendix Eq (S11)), we simulated an exponential gradient so that *τ*_{E} is a constant. In this case the cells reach steady state with a constant drift speed *V*_{D}. Calculating *V*_{D} from 10^{4} simulated trajectories for a range of *τ*_{E} and *τ*_{D0} = *τ*_{D}(*r*_{0}) values reveals that cells climb the gradient much faster when the positive feedback dominates (*τ*_{E} < 1) (Fig 1A). The trajectories of individual cells resembled that of a ratchet that moves almost only in one direction (Fig 1B green). In contrast, when the negative feedback dominates (*τ*_{E} > 1) the trajectories exhibit both up and down runs of similar although slightly biased lengths (Fig 1B red). *V*_{D} also depends on *τ*_{D} and peaks when the direction decorrelation time is on the same order as the memory duration (*τ*_{D} 1), consistent with previous studies [11, 23, 24]. In these simulations the adapted probability to run *r*_{0} = 0.8 and the ratio *D*_{T}/*D*_{R} = 37 were kept constant. Changing these values did not change the main results (S1A–S1C Fig).

In a wild type population, individual isogenic cells will have different values of *τ*_{E} and *τ*_{D0} due to cell-to-cell variabilities in swimming speed and in the abundance of chemotaxis proteins [22, 25, 26]. In a recent experimental study, the phenotype and performance of individual wild type cells (RP437 strain) was quantified by tracking cells swimming up a static quasi-linear gradient of methyl-aspartate (varying from 0 to 1 *mM* over 10 *mm*). This experiment revealed large differences among the performances of individual cells within the isogenic population [22], which could be reproduced by complementing the model of bacterial chemotaxis just described with a simple model of noisy gene expression (Fig 2 in [22]). To examine in which region of the (*τ*_{E}, *τ*_{D0}) space these cells might have been operating, we used this same model (complemented with diversity in rotational diffusion coefficients *D*_{R} and *D*_{T} due to variations in cell length; see S1 Appendix) with the same parameter values to run simulations of 16,000 cells climbing the experimentally measured (Fig 1A). We find that even in this relatively shallow gradient some cells might have been operating in the positive-feedback-dominated regime, especially near the bottom of the gradient (black dots). As the cells climb the gradient, *τ*_{E} becomes larger (white dots) because, as concentration increases, the log-sensing cells in the quasi-linear gradient face a shallower gradient, and thus weaker positive feedback.

To better understand the origin of the fast drift speed and its associated “ratchet-like” behavior, we examine the relationship between the drift speed *V*_{D} and the statistics of the internal state *f*. Using *t*_{M} as the unit of time and *v*_{0}*t*_{M} as the unit of length we derive a Fokker-Planck equation for the probability *P*(*x*, *f*, *s*, *τ*) that at time *τ* = *t*/*t*_{M} the cell is at position *x* = *X*/(*v*_{0}*t*_{M}) with internal state *f* and orientation *s* (Methods
Eq (14)):

$${\partial}_{\tau}P=-{\partial}_{f}\left(\left(-\left(f-{f}_{0}\right)+\frac{r\left(f\right)s}{{\tau}_{E}\left(x\right)}\right)P\right)+\frac{{\widehat{L}}_{s}P}{(n-1){\tau}_{D}\left(f\right)}-{\partial}_{x}\left(r\left(f\right)sP\right).$$

(5)

Here ${\widehat{L}}_{s}={(1-{s}^{2})}^{\frac{3-n}{2}}{\partial}_{s}\left({(1-{s}^{2})}^{\frac{n-1}{2}}{\partial}_{s}\right)$ is the rotational diffusion operator on the (*n* − 1)-sphere. All symbols used are summarized in Table 1.

For simplicity we consider a log-sensing organism swimming in a static exponential gradient. In this case, *τ*_{E}(*x*) = *τ*_{E} is constant (more complex gradient profiles and the effect of receptor saturation are considered later in the paper). Therefore the positive feedback becomes independent of position and the system can reach a steady state drift speed. Separating the variable *x* and integrating over *x* we obtain (Methods
Eq (17))

$$\begin{array}{c}\hfill {V}_{D}={\tau}_{E}\overline{\u27e8f-{f}_{0}\u27e9},\end{array}$$

(6)

where represents averaging over *f* and *s* and the bar indicates steady state. Eq (6) indicates that the drift speed is determined by the steady state marginal distribution $\overline{p}\left(f\right)$. To find an analytical expression for $\overline{p}\left(f\right)$, we expand the steady state joint distribution $\overline{P}(f,s)$ in orthonormal eigenfunctions of the angular operator ${\widehat{L}}_{s}$—the first two coefficients are the marginal distribution $\overline{p}\left(f\right)$ and the first angular moment ${\overline{p}}_{1}\left(f\right)/\sqrt{n}=\int \overline{P}s\mathrm{d}s$—and discard higher orders to obtain a closed system of equations. The analytical solution for the steady state marginal distribution $\overline{p}\left(f\right)$ reads

$$\overline{p}\left(f\right)=\frac{1}{W}\frac{\frac{r\left(f\right)}{{\tau}_{E}}}{\frac{1}{n}{\left(\frac{r\left(f\right)}{{\tau}_{E}}\right)}^{2}-{\left(f-{f}_{0}\right)}^{2}}exp\left(-{\int}^{f}\frac{\frac{{f}_{1}-{f}_{0}}{{\tau}_{D}\left({f}_{1}\right)}}{\frac{1}{n}{\left(\frac{r\left({f}_{1}\right)}{{\tau}_{E}}\right)}^{2}-{\left({f}_{1}-{f}_{0}\right)}^{2}}\mathrm{d}{f}_{1}\right),$$

(7)

where *W* is a normalization constant. The full derivation is provided in Methods Eqs (25)–(27), together with an interpretation of the distribution as a potential solution $\overline{p}\left(f\right)\propto exp(-V\left(f\right))$ where *V*(*f*) is the “potential”. We also examine how the shape of the potential depends on *τ*_{E} and *τ*_{D}.

The solution $\overline{p}\left(f\right)$ is plotted in Fig 1C. When the negative feedback dominates (*τ _{E}* 1) the distribution is sharply peaked and nearly Gaussian with variance ${\sigma}^{2}={\tau}_{D0}{r}_{0}^{2}/n{\tau}_{E}^{2}$ (Methods
Eq (32)) and its mean barely deviates from the adapted state

In the fast gradient climbing regime (*τ*_{E} 1) trajectories resemble that of a ratchet (Fig 1B). To gain mechanistic insight into this striking efficiency we examined the Langevin system equivalent to the Fokker-Planck Eq (5). Defining *v* = *rs* as the normalized run speed projected along the gradient, we change variables from (*f*, *s*) to (*r*, *v*) and obtain (Methods Eqs (47)–(52)):

$$\begin{array}{cc}\hfill \frac{\mathrm{d}r}{\mathrm{d}\tau}=& \phantom{\rule{3.33333pt}{0ex}}r\left(1-r\right)\left(-\left(f\left(r\right)-{f}_{0}\right)+\frac{v}{{\tau}_{E}}\right)\hfill \\ \hfill \frac{\mathrm{d}v}{\mathrm{d}\tau}=& \phantom{\rule{3.33333pt}{0ex}}v\left(1-r\right)\left(-\left(f\left(r\right)-{f}_{0}\right)+\frac{v}{{\tau}_{E}}\right)-\frac{v}{{\tau}_{D}\left(r\right)}+\sqrt{\frac{\left({r}^{2}-{v}^{2}\right)}{{\tau}_{D}\left(r\right)}}\eta \left(\tau \right),\hfill \end{array}$$

(8)

where *v* = *dx*/*dτ* and *η*(*τ*) denotes delta-correlated Gaussian white noise. The nullclines of the system (Fig 2A and 2C) intersect at the only stable fixed point (*r*, *v*) = (*r*_{0}, 0) of Eq (8) where the eigenvalues of the relaxation matrix

$$\left[\begin{array}{cc}-1& \left(1-{r}_{0}\right){r}_{0}/{\tau}_{E}\\ 0& -1/{\tau}_{D0}\end{array}\right]$$

(9)

are both negative (Methods
Eq (53)). Stochastic fluctuations due to rotational diffusions *D*_{R} and *D*_{T} (heat maps in Fig 2A and 2C) continuously push the system away from the fixed point. The magnitude of these fluctuations is large near the fixed point, causing the system to quickly move away. Fluctuations are smaller near *r* = 1 and *v* = 1, enabling the organism to climb the gradient at high speed for a longer time. Net drift results from spending more time in the region where *v* > 0.

Stochastic excursions in the (*r*, *v*)-plane away from the fixed point exhibit distinctive trajectories depending on the value of *τ*_{E}. When the positive feedback dominates (*τ*_{E} 1; Fig 2A) the eigenvectors of the relaxation matrix, (1, 0)^{T} and ${(\frac{(1-{r}_{0}){r}_{0}}{{\tau}_{E}}\frac{{\tau}_{D0}}{{\tau}_{D0}-1},1)}^{T}$, are highly non-orthogonal. This defines a non-normal dynamics that enables linear deviations to grow transiently [7, 8] to feed the nonlinear positive feedback (*v*^{2} term second line in Eq (8)) leading to large deviations. Importantly, this only happens for runs that start in the correct direction. If the run is in the wrong direction the linear deviation does not grow (Fig 2B; see also S1 Movie). Asymmetry arises because the *v*^{2} term is always positive. Similar selective amplification properties are observed in neuronal networks, where non-normal dynamics enables the network to respond to certain signals while ignoring others (including noise) [27, 28]. Thus, a random walker running in the correct direction is aided by the positive feedback, which pushes its internal dynamics towards the upper right corner of the phase plane where *r* = 1 and *v* = 1. If, instead, the run is in the wrong direction (*v* < 0), the nonlinearity pushes the system back into the high noise region near the fixed point where it will rapidly pick a new direction (Fig 2B).

In contrast, when the negative feedback dominates (*τ*_{E} 1; Fig 2C), the eigenvectors become nearly orthogonal. Linear deviations from the fixed point simply relax to the fixed point regardless of the initial direction of the run. Thus runs up and down the gradient are only marginally different in length, resulting in a small net drift (Fig 2D). This key difference between the positive and negative feedback regimes is reflected in the flow field (white curves in Fig 2A and 2C).

For simplicity in our analytical derivations we assumed the environment was a constant exponential gradient with concentrations in the (log-sensing) sensitivity range of the organism. Here we explore what happens when the organism encounters concentrations beyond its sensitivity range. For wild type *E. coli* the change in the free energy of the chemorecetor cluster due to ligand binding is proportional to ln((1 + *C*/*K*_{i})/(1 + *C*/*K*_{a})) (S1 Appendix Eq (S8)). Therefore the receptor is log-sensing to methyl-aspartate only for concentrations between *K*_{i} *C* *K*_{a}, where *K*_{i} = 0.0182 *mM* and *K*_{a} = 3 *mM* are the dissociation constants of the inactive and active states of the receptor. When *C* < *K*_{i} the receptor senses linear concentration [29], whereas when *C* > *K*_{a} the receptors saturate [30–33]: as a cell approaches a high concentration source its sensitivity decreases (S1 Appendix Eq (S10)). This in turn increases the value of *τ*_{E}. Simulations in an exponential gradient show that this effect results in an eventual slow-down as the cell approaches the source (Fig 3A–3C).

Realistic gradients are typically limited in spatial extent and often are not exponential, in which case *L* and therefore *τ*_{E} are different in different regions. *L* is long near the source in a linear gradient, for example, and decreases linearly with distance from the source. Simulations show that the cell initially climbs the gradient fast but later slows down as the gradient length scale *L* increases and *τ*_{E} increases (Fig 3D–3F). On the contrary, for a static localized source in three dimensions, *L* is short near the source but increases linearly with distance from it (Methods). Thus, *τ*_{E} decreases and the cell accelerates as it approaches the source (Fig 3G–3I). Comparing cells in various dynamical regimes (different values of *τ*_{E}) across these different gradients suggests that a lower value of *τ*_{E} results in faster gradient ascent.

When entering a food gradient, it is natural to try to climb as fast as possible. However this strategy could create a problem: the longer runs implied by the positive feedback mechanism could propel the accelerating *E. coli* beyond the nutrient source. This is the case in Fig 3E, where cells with the lowest *τ*_{E} (green) reach the source first but overshoot slightly; they settle, on average, at a further distance than those with intermediate *τ*_{E} (blue). Thus there is a trade-off between transient gradient climbing and long-term aggregating, as previously observed [13, 15, 23]. In nature, as chemotactic bacteria live in swarms, chasing and eating nutrient patches driven by flows and diffusions while new plumes of nutrients are constantly created by other organisms [2, 34], the actual environments experienced by bacteria are far more complex. The trade-off we found here hints that in these random small fluctuating gradients [11, 16, 35–37] the bacteria should not aim for maximal drift speed but need to deal with this trade-off to avoid overshoot. In general, natural environments will be complex, with a variety of different sources and gradients, implying different parameter domains will be optimal for *E. coli* at different times. Such phenotypic diversity may well confer an advantage [15, 37–41].

Our results illustrate the surprisingly new capabilities that can emerge when living systems exploit the full nonlinearity inherent within an otherwise simple and widely used strategy. For the particular case of bacterial chemotaxis we showed that cells that swim fast, have long memory (adaptation time), or large signal amplification, are likely to exhibit “ratchet-like” climbing behavior in a positive-feedback-dominated regime, even in shallow gradients. As we showed from simulations using a model that fits experimental data, this regime should be accessible to wild-type bacterial populations. Actually identifying these “ratcheting” cells from experimental trajectories would require observing them for a sufficient time (*T* *t*_{M}, *t*_{D}, *t*_{E}) and in a sufficiently steep gradient over the distance traveled (Δ*X* = *V*_{D}*T* ~ 0.5*v*_{0}*T*). Using parameter estimates from [13, 20], for *t*_{E} < *t*_{D} *t*_{M} ~ 10 *s* we take *T* ~ 200 *s*, and for *v*_{0} = 20 *μm*/*s* we get Δ*X* = 2 *mm*. To see how this compares with existing experimental setup with a quasi-linear gradient varying from 0 to 1 *mM* over 10 *mm* [22], we note that the black dots in Fig 1A show that some cells located 1.5 *mm* away from 0 concentration can operate in the positive-feedback-dominated regime. Thus, using the same setup as in [22] these requirements would be satisfied near the bottom of the gradient if the source concentration was increased to 3 *mM*.

It is common to make simplifying assumptions to facilitate analysis, but we do not believe that ours are limiting. We showed with simulations that our results hold (S1 Fig for details) when we take into account: (i) different values of *r*_{0} and *D*_{T}/*D*_{R}; (ii) the limited range of the receptor sensitivity [15, 18] (S1 Appendix Eq (S10)); (iii) possible nonlinearities (S1 Appendix Eq (S4)) and asymmetries of adaptation rates [14, 42]. A hallmark of *E. coli* chemotaxis is that, in the absence of a gradient, run-and-tumble behavior adapts back to prestimulus statistics [6, 43]. These robust properties of integral feedback control [6] remain in place in our study because the transients originate from non-normal dynamics around the stable fixed point. The boost from positive feedback described here is independent from other mechanisms that can enhance drift up a gradient such as imperfect adaptation in the response to some amino acids [44] and stochastic fluctuations in the adaptation mechanism [20, 21]. The latter has been shown to enhance chemotactic performance in shallow gradients by transiently pushing the system into a regime of slower direction changing provided it is running up the gradient. There are some similarities between the effect of signaling noise and the positive feedback mechanism presented here: both can affect drift speed by causing long-lasting asymmetries in the internal state when running up the gradient. In S3 Fig we show using simulations that signaling noise in the adaptation mechanism does not change our conclusion that the drift speed is maximal in the positive-feedback-dominated regime. Depending on the region of the (*τ*_{E}, *τ*_{D}) parameter space, the signaling noise can either enhance the drift speed by less than 10% or reduce it by up to 30%.

The fact that non-normal dynamics might be exploited to boost runs in the correct direction parallels recent findings in neuroscience [27] that suggest neuronal networks use similar strategies to selectively amplify neural activity patterns of interest. Thus, non-normal dynamics could be a feature that is selected for in living dynamical systems. Although we used bacterial chemotaxis as an example, our results do not depend on the specific form of the functions *r*(*f*) and *t*_{D}(*f*), provided they are increasing. Therefore our findings should be applicable to a large class of biased random walk strategies exhibited by organisms when local directional information is unreliable. In essence, any stochastic navigation strategy requires a memory, *t*_{M}, to make temporal comparisons, a reorientation mechanism, *t*_{D}, to sample new directions, and external information, *t*_{E}, relayed to decision-making circuitry through motion and signal amplification. Our theoretical contribution showed the (surprisingly) diverse behavioral repertoire that is possible by having these work in concert. In retrospect, perhaps this should not be surprising given the diverse environments in which running-and-tumbling organisms can thrive.

A detailed description of the chemotaxis model and agent-based simulations is provided in S1 Appendix (parameters in S1 Table). Briefly, the agent-based simulations were performed using Euler’s method to integrate a standard model of *E. coli* chemotaxis [12–14, 18, 20] in which the cell relays information from the external environment to the flagellar motor through a signaling cascade triggered by ligand-binding receptors. The receptors are described by a two-state model where the activity *a* is determined by the free energy difference *F* between the active and inactive states, which is determined by both the ligand concentration *C* and the receptors’ methylation level *m*. At each time step, the cell moves forward or stays in place according to its motility state (run or tumble), which also determines whether its direction changes with rotational diffusion coefficients *D*_{R} or *D*_{T}. At the new position changes in ligand concentration *C* cause changes in free energy *F* and thus activity *a*, and the methylation state adapts to compensate that change to maintain a constant activity. The updated value of the free energy *F* then determines the switching rates between the clockwise and counter-clockwise rotation of the flagellar motor state, which in turn determines the motility state of the cell according to rules and parameters in [20], completing one time step.

In Fig 1A we considered a wild type population in the scatter plot. To generate a population with realistic parameters, we used a recent model [22] of phenotypic diversity in *E. coli* chemotaxis that reproduces available experimental data on the laboratory strain RP437 climbing a linear gradient of methyl-aspartate. In this model individual cells have different abundances of the chemotaxis proteins (CheRBYZAW) and receptors (Tar, Tsr). These molecular abundances then determine the memory time *t*_{M} and the adapted probability to run *r*_{0} [15]. The run speed was different among cells and sampled from a Gaussian distribution to match experimental observations [22]. Rotational diffusion coefficients were also distributed to reflect differences in cell length.

We define ${P}_{R}(\mathbf{X},\widehat{\mathbf{u}},F,t)$ and ${P}_{T}(\mathbf{X},\widehat{\mathbf{u}},F,t)$ as the probability distributions at time *t* to be running or tumbling at position **X** in direction $\widehat{\mathbf{u}}$ with internal variable *F*. As described, there is Poisson switching between runs and tumbles with rates λ_{R}(*F*) and λ_{T}(*F*), runs and tumbles follow rotational diffusion with *D*_{R} and *D*_{T}, and motion is constant in runs and 0 in tumbles. Thus we construct a two-state stochastic master equation model [45]

$$\begin{array}{cc}\hfill {\partial}_{t}{P}_{R}=& -{\partial}_{F}\left(\dot{F}{|}_{run}{P}_{R}\right)-\nabla \xb7\left({v}_{0}\widehat{\mathbf{u}}{P}_{R}\right)+{\nabla}_{\widehat{\mathbf{u}}}^{2}\left({D}_{R}{P}_{R}\right)-{\lambda}_{R}{P}_{R}+{\lambda}_{T}{P}_{T}\hfill \\ \hfill {\partial}_{t}{P}_{T}=& -{\partial}_{F}\left(\dot{F}{|}_{tumble}{P}_{T}\right)+{\nabla}_{\widehat{\mathbf{u}}}^{2}\left({D}_{T}{P}_{T}\right)+{\lambda}_{R}{P}_{R}-{\lambda}_{T}{P}_{T},\hfill \end{array}$$

(10)

where $\dot{F}{|}_{run,tumble}$ are defined in Eq (2).

Since the gradient varies in one direction only we focus on motion in the gradient direction and integrate the probability over all other directions. Thus $\nabla \xb7\widehat{\mathbf{u}}=s{\partial}_{X}$ and ${\nabla}_{\widehat{\mathbf{u}}}^{2}={\widehat{L}}_{s}$, the polar angle part of the rotational diffusion operator on the (*n* − 1)-sphere. To derive the analytical form of ${\widehat{L}}_{s}$ we note in *n*-dimensional space we can iteratively write down the Laplace-Beltrami operator [46] as

$$\begin{array}{c}\hfill {\nabla}_{{S}^{n-1}}^{2}={\left(sin\theta \right)}^{2-n}{\partial}_{\theta}\left({\left(sin\theta \right)}^{n-2}{\partial}_{\theta}\right)+{\left(sin\theta \right)}^{-2}{\nabla}_{{S}^{n-2}}^{2},\end{array}$$

(11)

where 0 < *θ* < *π* is the polar angle. In a one-dimensional gradient we define the gradient direction as the polar axis, thus $s=\widehat{\mathbf{u}}\xb7\widehat{\mathbf{X}}=cos\theta $. We can write $sin\theta =\sqrt{1-{s}^{2}}$ and ${\partial}_{\theta}=-\sqrt{1-{s}^{2}}{\partial}_{s}$. Then the polar angle part is

$$\begin{array}{c}\hfill {L}_{s}={\left(1-{s}^{2}\right)}^{\frac{3-n}{2}}\widehat{{\partial}_{s}}\left({\left(1-{s}^{2}\right)}^{\frac{n-1}{2}}{\partial}_{s}\right).\end{array}$$

(12)

Using the definitions of the normalized internal state *f* = *HF*, of the timescale of switching between runs and tumbles *t*_{S} = 1/(λ_{R} + λ_{T}) [45], and of the probability to run *r* = λ_{T}/(λ_{R} + λ_{T}), we obtain

$$\begin{array}{cc}\hfill {\partial}_{t}{P}_{R}=& -{\partial}_{f}\left(\left(-\frac{f-{f}_{0}}{{t}_{M}}+\frac{NH{v}_{0}}{L}s\right){P}_{R}\right)+{D}_{R}{\widehat{L}}_{s}{P}_{R}\hfill \\ & -\frac{1-r}{{t}_{S}}{P}_{R}+\frac{r}{{t}_{S}}{P}_{T}-s{\partial}_{X}\left({v}_{0}{P}_{R}\right)\hfill \\ \hfill {\partial}_{t}{P}_{T}=& -{\partial}_{f}\left(-\frac{f-{f}_{0}}{{t}_{M}}{P}_{T}\right)+{D}_{T}{\widehat{L}}_{s}{P}_{T}\hfill \\ & +\frac{1-r}{{t}_{S}}{P}_{R}-\frac{r}{{t}_{S}}{P}_{T}.\hfill \end{array}$$

(13)

If we assume the switching terms with ${t}_{S}^{-1}$ in Eq (13) dominate, the probabilities to be running and tumbling equilibrate on a much faster timescale than the other ones. Therefore we can let *P* = *P*_{R} + *P*_{T} and can approximate the actual probability to run as *P*_{R}/*P* ≈ *r*. Adding the two equations above yields the Fokker-Planck equation:

$$\begin{array}{cc}\hfill {\partial}_{t}P\approx & -{\partial}_{f}\left(\left(-\frac{f-{f}_{0}}{{t}_{M}}+\frac{rs}{L/NH{v}_{0}}\right)P\right)\hfill \\ & +\left(r{D}_{R}+\left(1-r\right){D}_{T}\right){\widehat{L}}_{s}P-rs{\partial}_{X}\left({v}_{0}P\right).\hfill \end{array}$$

(14)

This is equivalent to a system of Langevin equations. Considering *dr*/*df* = *r*(1 − *r*) the internal variable dynamics (the first term on the right) gives Eq (3) which defines *t*_{E}. The angular dynamics (the second term on the right) defines *t*_{D}.

Using the time scale definitions in Eq (4) and non-dimensionalizing time *τ* = *t*/*t*_{M} and position *x* = *X*/(*v*_{0}*t*_{M}), we obtain the Fokker-Planck Eq (5).

From the Fokker-Planck Eq (5) we consider the steady state so that _{τ} = 0. For a log-sensing organism moving in an exponential gradient *τ*_{E} does not depend on *x*. We can therefore integrate over *x* to get an equation for the marginal steady state distribution $\overline{P}(f,s)$—this removes the _{x} term. Integrating over *s* gives

$$0=-{\partial}_{f}\left(-\left(f-{f}_{0}\right)\int \overline{P}w\left(s\right)\mathrm{d}s+\frac{r\left(f\right)}{{\tau}_{E}}\int s\overline{P}w\left(s\right)\mathrm{d}s\right),$$

(15)

where the bar indicates steady state. By the boundary conditions that *P* → 0 at ±∞, we must have

$$\begin{array}{c}\hfill r\left(f\right)\int s\overline{P}w\left(s\right)\mathrm{d}s={\tau}_{E}\left(f-{f}_{0}\right)\int \overline{P}w\left(s\right)\mathrm{d}s.\end{array}$$

(16)

From the −_{x}(*rsP*) term of the Fokker-Planck Eq (5), the spatial flux is *r*(*f*)*s* and the drift speed is its average over the distribution. Thus we get the drift speed as Eq (6)

$${V}_{D}=\overline{\u27e8rs\u27e9}=\int \int \phantom{\rule{0.277778em}{0ex}}r\left(f\right)s\overline{P}w\left(s\right)\mathrm{d}s\mathrm{d}f={\tau}_{E}\phantom{\rule{0.277778em}{0ex}}\int \int \phantom{\rule{0.277778em}{0ex}}\left(f-{f}_{0}\right)\overline{P}w\left(s\right)\mathrm{d}s\mathrm{d}f={\tau}_{E}\overline{\u27e8f-{f}_{0}\u27e9}.$$

(17)

Here we use separation of variables and expand the solution to the Fokker-Planck Eq (5) as a sum of eigenfunctions of the operator ${\widehat{L}}_{s}$ on *s*. We then ignore high order terms assuming *τ*_{D0} 1 and derive an approximate analytical solution.

The eigenvalue problem of the angular operator ${\widehat{L}}_{s}$, defined in Eq (12), is

(1 - *s*^{2})*y*^{″} - (*n* - 1)*s**y*^{′} = λ*y*.

(18)

We identify this as the Gegenbauer differential equation [47], with eigenfunctions the Gegenbauer polynomials ${C}_{k}^{(n/2-1)}\left(s\right)$ and the corresponding eigenvalues ${\lambda}_{k}^{(n/2-1)}=-k(k+n-2)$. When *n* = 3 they are Legendre polynomials with eigenvalues ${\lambda}_{k}^{(1/2)}=-k(k+1)$. The first few Gegenbauer polynomials are

$$\begin{array}{cc}\hfill {C}_{0}^{(n/2-1)}\left(s\right)& =1\hfill \\ \hfill {C}_{1}^{(n/2-1)}\left(s\right)& =\left(n-2\right)s\hfill \\ \hfill {C}_{2}^{(n/2-1)}\left(s\right)& =\frac{n-2}{2}\left(n{s}^{2}-1\right).\hfill \end{array}$$

(19)

They are orthogonal in the sense that

$${\int}_{-1}^{1}{C}_{k}^{(n/2-1)}}(s){C}_{l}^{(n/2-1)}(s){\left(1-{s}^{2}\right)}^{\frac{n-3}{2}}\text{d}s={N}_{k}^{(n/2-1)},$$

(20)

where the normalization constants are ${N}_{k}^{(n/2-1)}=\frac{\pi {2}^{4-n}(k+n-3)!}{k!(2k+n-2){(\Gamma (n/2-1\left)\right)}^{2}}$. When *n* = 3 they are ${N}_{k}^{(1/2)}=\frac{2}{2k+1}$, those of Legendre polynomials.

The weight in the integration above is consistent with the geometry on an (*n* − 1)-sphere *S*^{n−1}, whose the volume element are iteratively defined [46] as

d_{Sn-1}*ω* = (sin*θ*)^{n-2}d*θ*d_{Sn-2}*ω*.

(21)

After a change of variable *s* = cos *θ* and integrating over all remaining dimensions, we see that any integration of *s* should carry a weight

$$\begin{array}{c}\hfill w\left(s\right)\text{ds}={\left(1-{\mathrm{s}}^{2}\right)}^{\frac{\mathrm{n}-3}{2}}\mathrm{d}\mathrm{s}.\end{array}$$

(22)

From orthogonality and completeness, we write any function of *s*, in particular the probability distribution *P*, as a series of Gegenbauer polynomials. When *n* = 3 this is the Fourier-Legendre Series.

$$\begin{array}{cc}\hfill P(x,f,s,\tau )\phantom{\rule{3.33333pt}{0ex}}& =\sum _{k=0}^{\infty}{p}_{k}(x,f,\tau )\frac{{C}_{k}^{(n/2-1)}\left(s\right)}{\sqrt{{N}_{0}^{(n/2-1)}{N}_{k}^{(n/2-1)}}}\hfill \\ & =\frac{1}{{N}_{0}^{(n/2-1)}}\left({p}_{0}+{p}_{1}\sqrt{n}s+{p}_{2}\sqrt{\frac{n+2}{n-1}}\frac{n{s}^{2}-1}{2}+\cdots \right),\hfill \\ \hfill {p}_{k}(x,f,\tau )\phantom{\rule{3.33333pt}{0ex}}& ={\int}_{-1}^{1}\sqrt{\frac{{N}_{0}^{(n/2-1)}}{{N}_{k}^{(n/2-1)}}}{C}_{k}^{(n/2-1)}\left(s\right)P(x,f,s,\tau ){\left(1-{s}^{2}\right)}^{\frac{n-3}{2}}\mathrm{d}s,\hfill \end{array}$$

(23)

where we normalize the definitions to ensure ${p}_{0}={\int}_{-1}^{1}P{(1-{s}^{2})}^{\frac{n-3}{2}}\mathrm{d}s$ is the same as the marginal distribution. When *n* = 3, the above is

$$\begin{array}{cc}\hfill P(x,f,s,\tau )\phantom{\rule{3.33333pt}{0ex}}& =\sum _{k=0}^{\infty}{p}_{k}(x,f,\tau )\frac{\sqrt{2k+1}}{2}{C}_{k}^{(1/2)}\left(s\right)\hfill \\ & =\frac{1}{2}\left({p}_{0}+{p}_{1}\sqrt{3}s+{p}_{2}\sqrt{\frac{5}{2}}\frac{3{s}^{2}-1}{2}+\cdots \right),\hfill \\ \hfill {p}_{k}(x,f,\tau )\phantom{\rule{3.33333pt}{0ex}}& ={\int}_{-1}^{1}\sqrt{2k+1}{C}_{k}^{(1/2)}\left(s\right)P(x,f,s,\tau )\mathrm{d}s.\hfill \end{array}$$

(24)

From now on we denote the marginal distribution *p*(*f*) = *p*_{0}(*f*). Also, from this definition ${p}_{1}=\sqrt{n}{\int}_{-1}^{1}sP{(1-{s}^{2})}^{\frac{n-3}{2}}\mathrm{d}s$.

Substitute the expansion Eq (23) into the Fokker-Planck Eq (5) and use the orthogonality Eq (20), we obtain

$${\partial}_{\tau}{p}_{k}=-{\partial}_{f}\left(-\left(f-{f}_{0}\right){p}_{k}+\frac{r\left(f\right)}{{\tau}_{E}}{\widehat{s}}_{kl}{p}_{l}\right)+\frac{{\lambda}_{k}^{(n/2-1)}}{(n-1){\tau}_{D}\left(f\right)}{p}_{k}-{\partial}_{x}\left({\widehat{s}}_{kl}{p}_{l}\right),$$

(25)

where ${\widehat{s}}_{kl}=\sqrt{\frac{k(k+n-3)}{(2k+n-4)(2k+n-2)}}{\delta}_{k-1,l}+\sqrt{\frac{(k+1)(k+n-2)}{(2k+n-2)(2k+n)}}{\delta}_{k+1,l}$ (summation over *l* implied) is an operator relating neighboring orders. It comes from the positive feedback term. When *n* = 3 it is ${\widehat{s}}_{kl}=\frac{k}{\sqrt{4{k}^{2}-1}}{\delta}_{k-1,l}+\frac{k+1}{\sqrt{4{(k+1)}^{2}-1}}{\delta}_{k+1,l}$. The first few equations are

$$\begin{array}{cc}\hfill {\partial}_{\tau}p=& -{\partial}_{f}\left(-\left(f-{f}_{0}\right)p+\frac{r\left(f\right)}{{\tau}_{E}}\frac{1}{\sqrt{n}}{p}_{1}\right)-r\left(f\right){\partial}_{x}\frac{1}{\sqrt{n}}{p}_{1}\hfill \\ \hfill {\partial}_{\tau}{p}_{1}=& -{\partial}_{f}\left(-\left(f-{f}_{0}\right){p}_{1}+\frac{r\left(f\right)}{{\tau}_{E}}\left(\frac{1}{\sqrt{n}}p+\sqrt{\frac{2(n-1)}{n\left(n+2\right)}}{p}_{2}\right)\right)\hfill \\ & -\frac{1}{{\tau}_{D}\left(f\right)}{p}_{1}-r\left(f\right){\partial}_{x}\left(\frac{1}{\sqrt{n}}p+\sqrt{\frac{2(n-1)}{n\left(n+2\right)}}{p}_{2}\right)\hfill \\ \hfill {\partial}_{\tau}{p}_{2}=& -{\partial}_{f}\left(-\left(f-{f}_{0}\right){p}_{2}+\frac{r\left(f\right)}{{\tau}_{E}}\left(\sqrt{\frac{2(n-1)}{n\left(n+2\right)}}{p}_{1}+\sqrt{\frac{3n}{\left(n+2\right)\left(n+4\right)}}{p}_{3}\right)\right)\hfill \\ & -\frac{2n}{(n-1){\tau}_{D}\left(f\right)}{p}_{2}-r\left(f\right){\partial}_{x}\left(\sqrt{\frac{2(n-1)}{n\left(n+2\right)}}{p}_{1}+\sqrt{\frac{3n}{\left(n+2\right)\left(n+4\right)}}{p}_{3}\right).\hfill \end{array}$$

(26)

In the definition of ${\widehat{s}}_{kl}$, when *k* 1 the non-zero entries approach a constant 1/2. This means for large *k* the coefficients *p*_{k} in Eq (25) evolve similarly except that higher orders decay with faster rates *k*(*k* + *n* − 2)/(*n* − 1)*τ*_{D}. Therefore when *τ*_{D0} 1 we can neglect the 2nd and higher orders, which closes the infinite series of moment equations and leaves two equations concerning the zeroth and first marginal moments in *s*, *p*(*x*, *f*, *τ*) and *p*_{1}(*x*, *f*, *τ*) respectively. At steady state the approximation gives the analytical solution

$$\overline{p}\left(f\right)=\frac{1}{W}\frac{\frac{r\left(f\right)}{{\tau}_{E}}}{\frac{1}{n}{\left(\frac{r\left(f\right)}{{\tau}_{E}}\right)}^{2}-{\left(f-{f}_{0}\right)}^{2}}exp\left(-{\int}^{f}\frac{\frac{{f}_{1}-{f}_{0}}{{\tau}_{D}\left({f}_{1}\right)}}{\frac{1}{n}{\left(\frac{r\left({f}_{1}\right)}{{\tau}_{E}}\right)}^{2}-{\left({f}_{1}-{f}_{0}\right)}^{2}}\mathrm{d}{f}_{1}\right),$$

(27)

where *W* is a normalization constant. Eq (27) is the same as Eq (7) in the main text.

We can interpret the steady state distribution as a potential solution $\overline{p}\left(f\right)\propto exp(-V\left(f\right))$ where *V*(*f*) is the “potential”. In this case the equivalent “force” in internal state is

$$\begin{array}{cc}\hfill F\left(f\right)=& -{V}^{\prime}\left(f\right)=\frac{\mathrm{d}ln\overline{p}\left(f\right)}{\mathrm{d}f}\hfill \\ \hfill =& \frac{\mathrm{d}}{\mathrm{d}f}ln\frac{\frac{r\left(f\right)}{{\tau}_{E}}}{\frac{1}{n}{\left(\frac{r\left(f\right)}{{\tau}_{E}}\right)}^{2}-{\left(f-{f}_{0}\right)}^{2}}-\frac{\frac{f-{f}_{0}}{{\tau}_{D}\left(f\right)}}{\frac{1}{n}{\left(\frac{r\left(f\right)}{{\tau}_{E}}\right)}^{2}-{\left(f-{f}_{0}\right)}^{2}}.\hfill \end{array}$$

(28)

Since *τ*_{D0} 1 the second term dominates, making the “force” a spring-like system, with spring constant

$$k\left(f\right)=\frac{\frac{1}{{\tau}_{D}\left(f\right)}}{\frac{1}{n}{\left(\frac{r\left(f\right)}{{\tau}_{E}}\right)}^{2}-{\left(f-{f}_{0}\right)}^{2}}.$$

(29)

Three observations can be made from this spring constant in intuitively understanding the steady state distribution $\overline{p}\left(f\right)$. (i) *k*(*f*)→∞, i.e. the “spring” becomes infinitely “stiff”, when the denominator approaches 0. Therefore, the bounds of the distribution $\overline{p}\left(f\right)$ are proportional to 1/*τ*_{E}, the ratio between the positive and negative feedbacks (Eq (4)). Intuitively, a stronger positive feedback (smaller *τ*_{E}) drives the internal state *f* further away from *f*_{0}, so the spring constant *k*(*f*) is smaller and the distribution $\overline{p}\left(f\right)$ is wider. (ii) A slower change in direction (smaller *τ*_{D}) leads to a larger spring constant *k*(*f*) 1/*τ*_{D}(*f*), and thus the distribution $\overline{p}\left(f\right)$ is more concentrated near the “origin” *f*_{0}. Intuitively, a shorter direction correlation time *τ*_{D} inhibits coherent motion in a single direction, which is required by the positive feedback to consistently drive the internal state *f* away. Thus the distribution $\overline{p}\left(f\right)$ is more concentrated. (iii) Asymmetries are created by the functional dependencies of *r*(*f*) and *τ*_{D}(*f*), both increasing in *f*—a “weaker spring” for higher values of *f* shifts the distribution $\overline{p}\left(f\right)$ there. Intuitively, more positive feedback *r*(*f*) and more coherent motion *τ*_{D}(*f*) in the positive direction asymmetrically drives the internal state towards higher values. These 3 observations can all be found in Fig 1C.

We expand the steady state solution Eq (27) in orders of $\frac{1}{{\tau}_{E}}\u2aa11$ and *τ*_{D0} 1 and obtain a near-Gaussian approximation, from which we integrate using Eq (6) to obtain MFT results.

First, we write the steady state distribution Eq (27) as

$$\begin{array}{c}\hfill \overline{p\left(f\right)}=\frac{1}{W}B\left(f\right)exp\left(-{\int}^{f}A\left({f}_{1}\right)\mathrm{d}{f}_{1}\right).\end{array}$$

(30)

From the Taylor expansion of the integrand in the exponent

$$\begin{array}{cc}\hfill A\left(f\right)=& \frac{\frac{f-{f}_{0}}{{\tau}_{D}\left(f\right)}}{\frac{1}{n}{\left(\frac{r\left(f\right)}{{\tau}_{E}}\right)}^{2}-{\left(f-{f}_{0}\right)}^{2}}\hfill \\ \hfill =& \frac{n{\tau}_{E}^{2}}{{r}_{0}^{2}{\tau}_{D0}}\left(f-{f}_{0}\right)+\left(1+O\left(\frac{1}{{\tau}_{E}^{2}}\right)\right){\Sigma}_{m=1}^{\infty}\frac{{n}^{m+1}{\tau}_{E}^{2m+2}}{{r}_{0}^{2m+2}{\tau}_{D0}}{\left(f-{f}_{0}\right)}^{2m+1}\hfill \\ & -\left(1+O\left(\frac{1}{{\tau}_{E}^{2}}\right)\right){\Sigma}_{m=1}^{\infty}\frac{{n}^{m}{\tau}_{E}^{2m}{r}_{0}^{\prime}}{{r}_{0}^{2m+1}{\tau}_{D0}}\left(2m+\frac{{r}_{0}{\tau}_{D0}^{\prime}}{{r}_{0}^{\prime}{\tau}_{D0}}\right){\left(f-{f}_{0}\right)}^{2m},\hfill \end{array}$$

(31)

where ′ = *d*/*df*, we see that if we define

$${\sigma}^{2}=\frac{{r}_{0}^{2}{\tau}_{D0}}{n{\tau}_{E}^{2}},$$

(32)

the first term in *A*(*f*) will give $-\int A\left(f\right)\mathrm{d}f=-\frac{{(f-{f}_{0})}^{2}}{2{\sigma}^{2}}+...$. If we can show that the rest of the terms are small when $\frac{1}{{\tau}_{E}}<1$ and *τ*_{D0} < 1, we can write $\overline{p}\left(f\right)$ as a small deviation from a Gaussian.

Indeed, if we consider the integration range $|f-{f}_{0}|\sim \sigma \sim O(\sqrt{{\tau}_{D0}}/{\tau}_{E})$, we can write

$$\begin{array}{cc}\hfill -{\int}_{{f}_{0}}^{f}A\left({f}_{1}\right)\mathrm{d}{f}_{1}=& -\frac{{\left(f-{f}_{0}\right)}^{2}}{2{\sigma}^{2}}-{\Sigma}_{m=1}^{\infty}\frac{{\tau}_{D0}^{m}}{2m+2}\frac{{\left(f-{f}_{0}\right)}^{2m+2}}{{\sigma}^{2m+2}}\hfill \\ & +{\Sigma}_{m=1}^{\infty}\frac{{r}_{0}^{\prime}}{{r}_{0}}\frac{{\tau}_{D0}^{m-1}}{2m+1}\left(2m+\frac{{r}_{0}{\tau}_{D0}^{\prime}}{{r}_{0}^{\prime}{\tau}_{D0}}\right)\frac{{\left(f-{f}_{0}\right)}^{2m+1}}{{\sigma}^{2m}}+O\left(\frac{1}{{\tau}_{E}^{2}}\right)\hfill \\ \hfill =& -\frac{{\left(f-{f}_{0}\right)}^{2}}{2{\sigma}^{2}}+\frac{{r}_{0}^{\prime}}{{r}_{0}}\frac{1}{3}\left(2+\frac{{r}_{0}{\tau}_{D0}^{\prime}}{{r}_{0}^{\prime}{\tau}_{D0}}\right)\frac{{\left(f-{f}_{0}\right)}^{3}}{{\sigma}^{2}}-\frac{{\tau}_{D0}}{4}\frac{{\left(f-{f}_{0}\right)}^{4}}{{\sigma}^{4}}\hfill \\ & +\frac{{r}_{0}^{\prime}}{{r}_{0}}\frac{{\tau}_{D0}}{5}\left(4+\frac{{r}_{0}{\tau}_{D0}^{\prime}}{{r}_{0}^{\prime}{\tau}_{D0}}\right)\frac{{\left(f-{f}_{0}\right)}^{5}}{{\sigma}^{4}}+O\left(\frac{1}{{\tau}_{E}^{2}}\right)+O\left({\tau}_{D0}^{2}\right),\hfill \end{array}$$

(33)

Similarly, the prefactor is

$$B\left(f\right)=\frac{n{\tau}_{E}}{{r}_{0}}(1-\frac{{r}_{0}^{\prime}}{{r}_{0}}\left(f-{f}_{0}\right)+{\tau}_{D0}\frac{{\left(f-{f}_{0}\right)}^{2}}{{\sigma}^{2}}-3{\tau}_{D0}\frac{{r}_{0}^{\prime}}{{r}_{0}}\frac{{\left(f-{f}_{0}\right)}^{3}}{{\sigma}^{2}}+O\left(\frac{1}{{\tau}_{E}^{2}}\right)+O\left({\tau}_{D0}^{2}\right)).$$

(34)

Substitute Eqs (33) and (34) back into Eq (30) and taking care of the orders of all cross terms, we obtain

$$\begin{array}{cc}\hfill \overline{p}\left(f\right)& =\frac{1}{Z}\frac{{\mathrm{e}}^{-\frac{{\left(f-{f}_{0}\right)}^{2}}{2{\sigma}^{2}}}}{\sqrt{2\pi {\sigma}^{2}}}(1-\frac{{r}_{0}^{\prime}}{{r}_{0}}\left(f-{f}_{0}\right)+{\tau}_{D0}\frac{{\left(f-{f}_{0}\right)}^{2}}{{\sigma}^{2}}\hfill \\ & +\frac{{r}_{0}^{\prime}}{3{r}_{0}}\left(2-9{\tau}_{D0}+\frac{{r}_{0}{\tau}_{D0}^{\prime}}{{r}_{0}^{\prime}{\tau}_{D0}}\right)\frac{{\left(f-{f}_{0}\right)}^{3}}{{\sigma}^{2}}-\frac{{\tau}_{D0}}{4}\frac{{\left(f-{f}_{0}\right)}^{4}}{{\sigma}^{4}}\hfill \\ & +\frac{{r}_{0}^{\prime}}{60{r}_{0}}\left(103+32\frac{{r}_{0}{\tau}_{D0}^{\prime}}{{r}_{0}^{\prime}{\tau}_{D0}}\right){\tau}_{D0}\frac{{\left(f-{f}_{0}\right)}^{5}}{{\sigma}^{4}}\hfill \\ & -\frac{{r}_{0}^{\prime}}{12{r}_{0}}\left(2+\frac{{r}_{0}{\tau}_{D0}^{\prime}}{{r}_{0}^{\prime}{\tau}_{D0}}\right){\tau}_{D0}\frac{{\left(f-{f}_{0}\right)}^{7}}{{\sigma}^{6}}+O\left(\frac{1}{{\tau}_{E}^{3}}\right)+O\left({\tau}_{D0}^{\frac{5}{2}}\right)).\hfill \end{array}$$

(35)

with normalization constant *Z*.

We notice from Eq (30) that the range of distribution is bounded by *f*_{L} and *f*_{U}, defined by

$$\begin{array}{cc}\hfill {f}_{L}-{f}_{0}=& -\frac{1}{\sqrt{n}}\frac{r\left({f}_{L}\right)}{{\tau}_{E}},\hfill \\ \hfill {f}_{U}-{f}_{0}=& \frac{1}{\sqrt{n}}\frac{r\left({f}_{U}\right)}{{\tau}_{E}}.\hfill \end{array}$$

(36)

Since $\sigma =\frac{{r}_{0}\sqrt{{\tau}_{D0}}}{\sqrt{n}{\tau}_{E}}\u2aa1\frac{{r}_{0}}{\sqrt{n}{\tau}_{E}}$, we see that the integration range is much larger than the standard deviation of the Gaussian factor, and thus can be considered from −∞ to ∞. Therefore we get the normalization constant

$$Z=1+\frac{{\tau}_{D0}}{4}+O\left(\frac{1}{{\tau}_{E}^{2}}\right)+O\left({\tau}_{D0}^{2}\right).$$

(37)

Substitute Eqs (35) and (37) into Eq (6) and carry out the integrals

$${V}_{D}=\frac{{r}_{0}{\tau}_{D0}}{n{\tau}_{E}}\frac{\left(1-\frac{3}{4}{\tau}_{D0}\right)\left({r}_{0}^{\prime}+{r}_{0}\frac{{\tau}_{D0}^{\prime}}{{\tau}_{D0}}\right)+O\left(\frac{1}{{\tau}_{E}}\right)+O\left({\tau}_{D0}^{\frac{3}{2}}\right)}{1+\frac{{\tau}_{D0}}{4}+O\left(\frac{1}{{\tau}_{E}^{2}}\right)+O\left({\tau}_{D0}^{2}\right)}.$$

(38)

Finally, noticing that by the definition of *τ*_{D} in Eq (4)

$${\tau}_{D}\left(f\right)={\tau}_{D0}\frac{{r}_{0}{D}_{R}+\left(1-{r}_{0}\right){D}_{T}}{r\left(f\right){D}_{R}+\left(1-r\left(f\right)\right){D}_{T}},$$

(39)

we can get

$${\tau}_{D0}^{\prime}={\tau}_{D0}\frac{{D}_{T}-{D}_{R}}{{r}_{0}{D}_{R}+\left(1-{r}_{0}\right){D}_{T}}{r}_{0}^{\prime}.$$

(40)

Therefore

$$\begin{array}{cc}\hfill {r}_{0}^{\prime}+{r}_{0}\frac{{\tau}_{D0}^{\prime}}{{\tau}_{D0}}& =\frac{{\tau}_{D0}^{\prime}}{{\tau}_{D0}}\left(\frac{{r}_{0}{D}_{R}+\left(1-{r}_{0}\right){D}_{T}}{{D}_{T}-{D}_{R}}+{r}_{0}\right)\hfill \\ & =\frac{{\tau}_{D0}^{\prime}}{{\tau}_{D0}}\frac{{D}_{T}}{{D}_{T}-{D}_{R}}.\hfill \end{array}$$

(41)

Taking *D*_{T} *D*_{R}, we put this back into Eq (38) and get

$$\begin{array}{cc}\hfill {V}_{D}& =\frac{{r}_{0}{\tau}_{D0}^{\prime}}{n{\tau}_{E}}\frac{1-\frac{3}{4}{\tau}_{D0}+O\left(\frac{1}{{\tau}_{E}}\right)+O\left({\tau}_{D0}^{\frac{3}{2}}\right)}{1+\frac{{\tau}_{D0}}{4}+O\left(\frac{1}{{\tau}_{E}^{2}}\right)+O\left({\tau}_{D0}^{2}\right)}\hfill \\ & =\frac{{r}_{0}{\tau}_{D0}^{\prime}}{n{\tau}_{E}\left(1+{\tau}_{D0}\right)}\left(1+O\left(\frac{1}{{\tau}_{E}}\right)+O\left({\tau}_{D0}^{\frac{3}{2}}\right)\right).\hfill \end{array}$$

(42)

When converted back to real units (*t* instead of *τ* = *t*/*t*_{M}), the highest-order term is identical, except for notations, to Eq (3) in Dufour *et al.* [13] obtained from a different approach. It can also be reduced to Eq (12) in Si *et al.* [12] by assuming a high running probability and a long memory. It agrees with Eq (6.24) in Erban & Othmer [48] and Eq (16) in Franz *et al.* [49] with appropriate inclusion of rotational diffusion.

In Eq (38) we expanded the distribution as a near-Gaussian around *f*_{0}. From Eq (6) we see the mean internal state ${f}_{m}=\overline{\langle f\rangle}$ has a slight shift, so it’s more accurate to expand around *f*_{m}. From Eqs (38) and (6) in the main text, we see $\overline{\langle f-{f}_{0}\rangle}\sim O(1/{\tau}_{E}^{2})$. Thus considering the shift in *f*_{m} the resulting *V*_{D} has the same form compared to Eq (42):

$${V}_{D}=\frac{{r}_{m}{\tau}_{Dm}^{\prime}}{n{\tau}_{E}\left(1+{\tau}_{Dm}\right)}\left(1+O\left(\frac{1}{{\tau}_{E}}\right)+O\left({\tau}_{Dm}^{\frac{3}{2}}\right)\right).$$

(43)

The first term in Eq (5) says the flux in *f*-space is non-negative provided −(*f* − *f*_{0}) + *rs*/*τ*_{E} > 0, or, noting *s* ≤ 1

(44)

Thus the upper bound *f*_{U} of the distribution *p*(*f*) is achieved at equality. Similarly, the lower bound *f*_{L} is achieved when we take equal signs of

(45)

noting *s* ≥ 1.

When *τ*_{E} becomes small we note *f*_{U,L} deviates far away from *f*_{0} as 1/*τ*_{E} → ∞. Using the definition *r* = 1/(1 + exp(−*f*)), we write

$$\begin{array}{c}\hfill {f}_{L,U}=\mp \frac{1}{{\tau}_{E}\left(1+exp\left(-{f}_{L,U}\right)\right)}.\end{array}$$

(46)

The plus sign gives exp(−*f*_{U}) 1 and *f*_{U} ≈ 1/*τ*_{E}. The minus sign gives exp(−*f*_{L}) 1 and *f*_{L} = −exp(*f*_{L})/*τ*_{E}. Taking logarithm, the latter gives *f*_{L} = ln (|*f*_{L}|*τ*_{E}) ≈ ln *τ*_{E}.

To derive Langevin equations from the Fokker-Planck equation we need to consider the geometric weight factor *w*(*s*) in Eq (22) for anglular integration. In deriving *s*-dynamics, we start with the angular part of the Fokker-Planck Eq (5)

$${\partial}_{\tau}P=\frac{{\widehat{L}}_{s}P}{(n-1){\tau}_{D}\left(f\right)}+\dots .$$

(47)

Multiplying an arbitrary function *A*(*s*) and integrating over all dimensions, we obtain

$$\begin{array}{cc}& \int \mathrm{d}x\int \mathrm{d}f{\int}_{-1}^{1}w\left(s\right)\mathrm{d}sA\left(s\right){\partial}_{\tau}P\hfill \\ \hfill =& \int \mathrm{d}x\int \mathrm{d}f{\int}_{-1}^{1}w\left(s\right)\mathrm{d}sA\left(s\right)\frac{{\left(1-{s}^{2}\right)}^{\frac{3-n}{2}}{\partial}_{s}\left({\left(1-{s}^{2}\right)}^{\frac{n-1}{2}}{\partial}_{s}P\right)}{(n-1){\tau}_{D}\left(f\right)}\hfill \\ \hfill =& -\int \mathrm{d}x\int \mathrm{d}f{\int}_{-1}^{1}w\left(s\right)\mathrm{d}s\frac{sP}{{\tau}_{D}\left(f\right)}{\partial}_{s}A\left(s\right)\hfill \\ & +\int \mathrm{d}x\int \mathrm{d}f{\int}_{-1}^{1}w\left(s\right)\mathrm{d}s\frac{\left(1-{s}^{2}\right)P}{(n-1){\tau}_{D}\left(f\right)}{\partial}_{s}^{2}A\left(s\right).\hfill \end{array}$$

(48)

To apply the standard result of equivalence between Fokker-Planck equations and Langevin equations, we need to change the measure in *s*-space to unity. This prompts the definition *Q*(*s*, *t*) = *w*(*s*)*P*(*y*, *f*, *s*, *t*)d*x*d*f* so that the above becomes

$$\begin{array}{cc}\hfill {\int}_{-1}^{1}\mathrm{d}sA\left(s\right){\partial}_{\tau}Q=& -{\int}_{-1}^{1}\mathrm{d}s\frac{sQ}{{\tau}_{D}\left(f\right)}{\partial}_{s}A\left(s\right)+{\int}_{-1}^{1}\mathrm{d}s\frac{\left(1-{s}^{2}\right)Q}{(n-1){\tau}_{D}\left(f\right)}{\partial}_{s}^{2}A\left(s\right)\hfill \\ \hfill =& {\int}_{-1}^{1}\mathrm{d}sA\left(s\right){\partial}_{s}\frac{sQ}{{\tau}_{D}\left(f\right)}+{\int}_{-1}^{1}\mathrm{d}sA\left(s\right){\partial}_{s}^{2}\frac{\left(1-{s}^{2}\right)Q}{(n-1){\tau}_{D}\left(f\right)},\hfill \end{array}$$

(49)

where we integrated by parts and discarded boundary terms. Since *A*(*s*) is an arbitrary function, we can write down the Fokker-Planck equation

$${\partial}_{\tau}Q={\partial}_{s}\frac{sQ}{{\tau}_{D}\left(f\right)}+{\partial}_{s}^{2}\frac{\left(1-{s}^{2}\right)Q}{(n-1){\tau}_{D}\left(f\right)},$$

(50)

which is equivalent [45] to the Langevin equation

$$\frac{\mathrm{d}s}{\mathrm{d}\tau}=-\frac{s}{{\tau}_{D}\left(r\right)}+\sqrt{\frac{2\left(1-{s}^{2}\right)}{(n-1){\tau}_{D}\left(r\right)}}\eta \left(\tau \right).$$

(51)

where *η*(*τ*) denotes the Gaussian white noise with *η*(*τ*_{1})*η*(*τ*_{2}) = *δ*(*τ*_{1} − *τ*_{2}).

The other two variables follow standard results [45] from the Fokker-Planck Eq (5) in the main text

$$\begin{array}{cc}\hfill \frac{\mathrm{d}f}{\mathrm{d}\tau}& =-\left(f-{f}_{0}\right)+\frac{r\left(f\right)s}{{\tau}_{E}},\hfill \\ \hfill \frac{\mathrm{d}x}{\mathrm{d}\tau}& =r\left(f\right)s.\hfill \end{array}$$

(52)

Now we change variables according to the definitions *r*(*f*) = 1/(1 + exp(−*f*)) and *v* = *rs*, and derive from the above dynamics in Eqs (51) and (52) to get the Langevin Eq (8).

Near the fixed point (*r*_{0}, 0), the eigenvectors and eigenvalues of the linearized Langevin Eq (8) are:

$$\begin{array}{cc}\hfill \left[\begin{array}{c}1\\ 0\end{array}\right]& \phantom{\rule{4.pt}{0ex}}\phantom{\rule{4.pt}{0ex}}\phantom{\rule{4.pt}{0ex}}\text{for}\phantom{\rule{4.pt}{0ex}}\text{eigenvalue}\phantom{\rule{4.pt}{0ex}}\phantom{\rule{4.pt}{0ex}}\phantom{\rule{4.pt}{0ex}}-1;\hfill \\ \hfill \left[\begin{array}{c}\frac{\left(1-{r}_{0}\right){r}_{0}}{{\tau}_{E}}\frac{{\tau}_{D0}}{{\tau}_{D0}-1}\\ 1\end{array}\right]& \phantom{\rule{4.pt}{0ex}}\phantom{\rule{4.pt}{0ex}}\phantom{\rule{4.pt}{0ex}}\text{for}\phantom{\rule{4.pt}{0ex}}\text{eigenvalue}\phantom{\rule{4.pt}{0ex}}\phantom{\rule{4.pt}{0ex}}\phantom{\rule{4.pt}{0ex}}-\frac{1}{{\tau}_{D0}}.\hfill \end{array}$$

(53)

When *τ*_{E} is large, $\frac{(1-{r}_{0}){r}_{0}}{{\tau}_{E}}\frac{{\tau}_{D0}}{{\tau}_{D0}-1}\u2aa11$ and the eigenvectors are almost orthogonal. When *τ*_{E} is small, $\frac{(1-{r}_{0}){r}_{0}}{{\tau}_{E}}\frac{{\tau}_{D0}}{{\tau}_{D0}-1}\gg 1$ and the eigenvectors are not orthogonal.

In Fig 1A heat map the drift speed *V*_{D} was calculated by fitting the linear part of the mean trajectory. In Fig 1B the first 50 *s* were removed to avoid the start up transient. In Fig 1C, the steady state $\overline{p}\left(f\right)$ from agent-based simulations was calculated from the histogram of all the internal values of the 10^{4} simulated cells between *τ* = 10 and *τ* = 20, sampled at regular steps of *τ* = 0.01. Numerical solutions of the Fokker-Planck Eq (5) were obtained by expanding the distribution in angles, as in Eq (25), and keeping the first 10 orders. The steady state $\overline{p}\left(f\right)$ was found by solving an initial value problem using the NDSolve function in Mathematica, with 10^{4} spatial points and integration time up to *τ* = 10. Further orders, finer grid, and longer integration times were checked to ensure solution accuracy. In Fig 1D, *V*_{D} from agent-based and Fokker-Planck were calculated by plugging into Eq (6)
$\overline{p}\left(f\right)$ obtained from those methods in C. MFT was calculated by combining Eq (42) with Eq (6) to find both ${f}_{m}=\overline{\langle f\rangle}$ and *V*_{D} [12, 13]. In the inset, the black curves show the approximate distribution in Eq (35).

In Fig 2B and 2D the Langevin trajectories were generated using Euler’s method to integrate Eq (8).

In Fig 3C, 3F and 3I the *τ*_{E} calculation considered receptor saturation as well as the varying gradient length scales, with *C* and *L* evaluated at mean positions. Note this is not the average *τ*_{E} over the population.

(PDF)

Click here for additional data file.^{(160K, pdf)}

(A) Same as Fig 1A (where *r*_{0} = 0.8, Table S1 Table) except with *r*_{0} = 0.5. (B) *r*_{0} = 0.7. (C) Same as Fig 1A (where *D*_{T}/*D*_{R} ≈ 37, Table S1 Table) except with *D*_{T}/*D*_{R} = 5. (D) Same as Fig 1A except without assuming receptors in log-sensing range, i.e. Eq (S7) was used rather than Eq (S11). (E) Same as D, but additionally implements adaptation asymmetry, where the adaptation rate ${t}_{M}^{-1}$ in equation Eq (S4) depends on *m* [42]. Here the adaptation is 3 times faster when *m* > *m*(*C*) than when *m* < *m*(*C*), and *t*_{M} is defined as the time scale when *m* < *m*(*C*). (F) Same as D but with nonlinear adaptation rate Eq (S4). Values for the parameters are: *a*_{B} = 0.74, *r*_{B} = 4.0, *K*_{R} = 0.32, *K*_{B} = 0.30, and *V*_{R} and *V*_{B}(0) chosen to ensure $\frac{\mathrm{d}m}{\mathrm{d}t}=0$ when *a* = *a*_{0} and the adaptation time is *t*_{M} when linearized.

(TIF)

Click here for additional data file.^{(11M, tif)}

5 sample trajectories (thin solid curves) and the mean over 10^{4} sample trajectories (thick lines) of non-dimensionalized position *x* = *X*/(*v*_{0}*t*_{M}) as a function of time *τ* = *t*/*t*_{M}. Colors correspond to cells with matching (green *τ*_{D0} = 1) and non-matching (orange *τ*_{D0} = 0.1) reorientation as in Fig 1 (same methods).

(TIF)

Click here for additional data file.^{(906K, tif)}

(A) Same as Fig 1A but adding a signaling noise term ${\sigma}_{m}\sqrt{2/{t}_{M}}\Gamma \left(t\right)$ in Eq (S6) where Γ(*t*) is the standard Wiener process (see Eq [4] in [20]). From Eqs (S1)-(S2) and *Y* = *αa* we obtain *dY*/*dm* = *Y*(1 − *a*)*ϵ*_{1}. Then plugging in *σ*_{Y}/*Y* = 0.1 [20, 21], *a* ≈ 0.5 and *ϵ*_{1} = −1, we obtain *σ*_{m} = 0.2. (B) The absolute difference in drift speed between A (with signaling noise) and Fig 1A (without signaling noise), Δ*V*_{D} = *V*_{D}|_{noise} − *V*_{D}|_{nonoise}, shows how signaling noise can either enhance or reduce the drift speed depending on (*τ*_{E}, *τ*_{D}). Note the colorbar is different. Black contour lines show level sets of Δ*V*_{D} at the colorbar ticks (-0.05, -0.025, 0, 0.025, and 0.05). (C) The relative difference in drift speed (B divided by the drift speed without signaling noise). Again, the color scale is different and black contour lines show level sets at the colorbar ticks.

(TIF)

Click here for additional data file.^{(5.4M, tif)}

(PDF)

Click here for additional data file.^{(116K, pdf)}

(AVI)

Click here for additional data file.^{(9.0M, avi)}

We thank D Clark, Y Dufour, N Frankel, X Fu, S Kato, N Olsman, DC Vural, and A Waite for discussions. This work was supported by the HPC facilities operated by, and the staff of, the Yale Center for Research Computing.

JL received support from the Natural Sciences and Engineering Research Council of Canada (NSERC) Postgraduate Scholarships-Doctoral Program (http://www.nserc-crsng.gc.ca/Students-Etudiants/PG-CS/BellandPostgrad-BelletSuperieures_eng.asp) PGSD2-471587-2015. TE received support from the National Institute of General Medical Sciences (www.nigms.nih.gov) grant 4R01GM106189-04. TE, JL and SWZ were supported by the Allen Distinguished Investigator Program (grant 11562) through the Paul G. Allen Frontiers Group (www.pgafamilyfoundation.org/programs/investigators-fellows). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability

All relevant data are within the paper and its Supporting Information files.

1.
Berg HC, Brown DA. Chemotaxis in *Escherichia coli* analysed by three-dimensional tracking. Nature. 1972.
October;239:500–504. doi: 10.1038/239500a0
[PubMed]

2.
Stocker R. Marine microbes see a sea of gradients. Science. 2012.
November;338:628–633. doi: 10.1126/science.1208929
[PubMed]

3.
Albrecht DR, Bargmann CI. High-content behavioral analysis of *Caenorhabditis elegans* in precise spatiotemporal chemical environments. Nat Methods. 2011.
July;8(7):599–606. doi: 10.1038/nmeth.1630
[PMC free article] [PubMed]

4.
Gomez-Marin A, Stephens GJ, Louis M. Active sampling and decision making in *Drosophila* chemotaxis. Nat Commun. 2011.
March;2:441
doi: 10.1038/ncomms1455
[PMC free article] [PubMed]

5.
Taylor-King JP, Franz B, Yates CA, Erban R. Mathematical modelling of turning delays in swarming robots. IMA J Appl Math. 2015.
March;80:1454–1474.

6.
Yi TM, Huang Y, Simon MI, Doyle J. Robust perfect adaptation in bacterial chemotaxis through integral feedback control. Proc Natl Acad Sci U S A. 2000.
April;97(9):4649–4653. doi: 10.1073/pnas.97.9.4649
[PubMed]

7.
Trefethen LN, Trefethen AE, Reddy SC, Driscoll TA. Hydrodynamic Stability Without Eigenvalues. Science. 1993.
July;261:578–584. doi: 10.1126/science.261.5121.578
[PubMed]

8.
Schmid PJ. Nonmodal Stability Theory. Annu Rev Fluid Mech. 2007.
January;39:129–162. doi: 10.1146/annurev.fluid.38.050304.092139

9.
Keller EF, Segel LA. Traveling bands of chemotactic bacteria: a theoretical analysis. J Theor Biol. 1971;30:235–248. doi: 10.1016/0022-5193(71)90051-8
[PubMed]

10.
Schnitzer MJ. Theory of continuum random walks and application to chemotaxis. Phys Rev E. 1993.
October;48(4):2553–2568. doi: 10.1103/PhysRevE.48.2553
[PubMed]

11.
Celani A, Vergassola M. Bacterial strategies for chemotaxis response. Proc Natl Acad Sci U S A. 2010.
January;107(4):1391–1396. doi: 10.1073/pnas.0909673107
[PubMed]

12.
Si G, Wu T, Ouyang Q, Tu Y. Pathway-based mean-field model for *Escherichia coli* chemotaxis. Phys Rev Lett. 2012.
July;109(4):048101
doi: 10.1103/PhysRevLett.109.048101
[PMC free article] [PubMed]

13.
Dufour YS, Fu X, Hernandez-Nunez L, Emonet T. Limits of feedback control in bacterial chemotaxis. PLoS Comput Biol. 2014.
June;10:e1003694
doi: 10.1371/journal.pcbi.1003694
[PMC free article] [PubMed]

14.
Shimizu TS, Tu Y, Berg HW. A modular gradient-sensing network for chemotaxis in *Escherichia coli* revealed by responses to time-varying stimuli. Mol Syst Biol. 2010.
June;6:382–395. doi: 10.1038/msb.2010.37
[PMC free article] [PubMed]

15.
Frankel NW, Pontius W, Dufour YS, Long J, Hernandez-Nunez L, et al.
Adaptability of non-genetic diversity in bacterial chemotaxis. eLife. 2014.
October; doi: 10.7554/eLife.03526
[PMC free article] [PubMed]

16.
Zhu X, Si G, Deng N, Ouyang Q, Wu T, et al.
Frequency-dependent *Escherichia coli* chemotactic behavior. Phys Rev Lett. 2012.
March;108(12):128101
doi: 10.1103/PhysRevLett.108.128101
[PMC free article] [PubMed]

17.
Xue C, Yang X. Moment-flux models for bacterial chemotaxis in large signal gradients. J Math Biol. 2016.
February;:1–24. doi: 10.1007/s00285-016-0981-9
[PubMed]

18.
Tu Y. Quantitative modeling of bacterial chemotaxis signal amplification and accurate adaptation. Annu Rev Biophys. 2013.
February;42:337–359. doi: 10.1146/annurev-biophys-083012-130358
[PMC free article] [PubMed]

19.
Saragosti J, Silberzan P, Buguin A. Modeling *E. coli* tumbles by rotational diffusion. Implications for chemotaxis. PLoS ONE. 2012.
April;7(4):e35412
doi: 10.1371/journal.pone.0035412
[PMC free article] [PubMed]

20.
Sneddon MW, Pontius W, Emonet T. Stochastic coordination of multisple actuators reduce latency and improves chemotactic response in bacteria. Proc Natl Acad Sci U S A. 2012.
January;109(2):805–810. doi: 10.1073/pnas.1113706109
[PubMed]

21.
Flores M, Shimizu TS, ten Wolde PR, Tostevin F. Signaling Noise Enhances Chemotactic Drift of *E. coli*. Phys Rev Lett. 2012.
October;109:148101
doi: 10.1103/PhysRevLett.109.148101
[PubMed]

22.
Waite AJ, Frankel NW, Dufour YS, Johnston JF, Long J, Emonet T. Non-genetic diversity modulates population performance. Mol Syst Biol. 2016. accepted. doi: 10.15252/msb.20167044
[PMC free article] [PubMed]

23.
Clark DA, Grant LC. The bacterial chemotactic response reflects a compromise between transient and steady-state behavior. Proc Natl Acad Sci U S A. 2005.
June;102(26):9150–9155. doi: 10.1073/pnas.0407659102
[PubMed]

24.
Kafri Y, da Silveira RA. Steady-state chemotaxis in *Escherichia coli*. Phys Rev Lett. 2008.
June;100(23):238101
doi: 10.1103/PhysRevLett.100.238101
[PubMed]

25.
Dufour YS, Gillet S, Frankel NW, Weibel DB, Emonet T. Direct correlation between motile behaviors and protein abundance in single cells. PLoS Comput Biol. 2016.
September;12(9):e1005041
doi: 10.1371/journal.pcbi.1005041
[PMC free article] [PubMed]

26.
Spudich JL, Koshland DE Jr. Non-genetic individuality: chance in the single cell. Nature. 1976.
August;262:467–471. doi: 10.1038/262467a0
[PubMed]

27.
Murphy BK, Miller KD. Balanced amplification: a new mechanism of selective amplification of neural activity patterns. Neuron. 2009.
February;61:635–648. doi: 10.1016/j.neuron.2009.02.005
[PMC free article] [PubMed]

28.
Hennequin G, Vogels TP, Gerstner W. Non-normal amplification in random balanced neural networks. Phys Rev E. 2012.
July;86(1):011909
doi: 10.1103/PhysRevE.86.011909
[PubMed]

29.
Colin R, Zhang R, Wilson LG. Fast, hight-throughput measurement of collective behavior in a bacterial population. J R Soc Interface. 2014.
July;11:20140486
doi: 10.1098/rsif.2014.0486
[PMC free article] [PubMed]

30.
Mello BA, Tu Y. Quantitative modeling of sensitivity in bacterial chemotaxis: The role of coupling among different chemoreceptor species. Proc Natl Acad Sci U S A. 2003.
July; 100(14):8223–8228. doi: 10.1073/pnas.1330839100
[PubMed]

31.
Sourjik V, Berg HC. Functional interactions between receptors in bacterial chemotaxis. Nature. 2004.
March;428:437–441. doi: 10.1038/nature02406
[PubMed]

32.
Keymer JE, Endres RG, Skoge M, Meir Y, Wingreen NS. Chemosensing in *Escherichia coli*: Two regimes of two-state receptors. Proc Natl Acad Sci U S A. 2006.
February;103(6):1786–1791. doi: 10.1073/pnas.0507438103
[PubMed]

33.
Hansen CH, Endres RG, Wingreen NS. Chemotaxis in *Escherichia coli*: A Molecular Model for Robust Precise Adaptation. PLoS Comput Biol. 2008.
January;4(1):e1
doi: 10.1371/journal.pcbi.0040001
[PubMed]

34.
Blackburn N, Fenchel T, Mitchell J. Microscale Nutrient Patches in Planktonic Habitats Shown by Chemotactic Bacteria. Science. 1989.
December;282:2254–2256. doi: 10.1126/science.282.5397.2254
[PubMed]

35.
Clausznitzer D, Micali G, Neumann S, Sourjik V, Endres RG. Predicting Chemical Environments of Bacteria from Receptor Signaling. PLoS Comput Biol. 2014.
October;10(10):e1003870
doi: 10.1371/journal.pcbi.1003870
[PMC free article] [PubMed]

36.
Hein AM, Brumley DR, Carrara F, Stocker R, Levin SA. Physical limits on bacterial navigation in dynamic environments. J R Soc Interface. 2016.
January;13:20150844
doi: 10.1098/rsif.2015.0844
[PMC free article] [PubMed]

37.
Edgington MP, Tindall MJ. Understanding the link between single cell and population scale responses of *Escherichia coli* in differing ligand gradients. Comput Struct Biotechnol J. 2015.
October;13:528–538. doi: 10.1016/j.csbj.2015.09.003
[PMC free article] [PubMed]

38.
Elowitz MB, Levine AJ, Siggia ED, Swain PS. Stochastic Gene Expression in a Single Cell. Science. 2002.
August;297:1183–1186. doi: 10.1126/science.1070919
[PubMed]

39.
Kussell E, Leibler S. Phenotypic Diversity, Population Growth, and Information in Fluctuating Environments. Science. 2005.
September;309:2075–2078. doi: 10.1126/science.1114383
[PubMed]

40.
Acar M, Mettetal JT, van Oudenaarden A. Stochastic switching as a survival strategy in fluctuating environments. Nat Genet. 2008.
April;40(4):471–475. doi: 10.1038/ng.110
[PubMed]

41.
Ackermann M. A functional perspective on phenotypic heterogeneity in microorganisms. Nat Rev Microbiol. 2015.
August;13:497–508. doi: 10.1038/nrmicro3491
[PubMed]

42.
Clausznitzer D, Oleksiuk O, Løvdok L, Sourjik V, Endres RG. Chemotactic Response and Adaptation Dynamics in *Escherichia coli*. PLoS Comput Biol. 2010.
May;6(5):e1000784
doi: 10.1371/journal.pcbi.1000784
[PMC free article] [PubMed]

43.
Block SM, Segall JE, Berg HC. Impulse Responses in Bacterial Chemotaxis. Cell. 1982.
November;31:215–226. doi: 10.1016/0092-8674(82)90421-4
[PubMed]

44.
Wong-Ng J, Melbinger A, Celani A, Vergassola. The Role of Adaptation in Bacterial Speed Races. PLoS Comput Biol. 2016.
June;12(6):e1004974
doi: 10.1371/journal.pcbi.1004974
[PMC free article] [PubMed]

45.
Gardiner CW. Handbook of stochastic methods for physics, chemistry and the natural sciences. Berlin: Springer Berlin Heidelberg; 2004.
doi: 10.1007/978-3-662-05389-8

46.
Jost J. Riemannian Geometry and Geometric Analysis. Berlin: Springer Berlin Heidelberg; 2011.
doi: 10.1007/978-3-642-21298-7

47.
Abramowitz M, Stegun I. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Mineola: Dover Publications; 1964.

48.
Erban R, Othmer HG. From individual to collective behavior in bacteria chemotaxis. SIAM J Appl Math. 2004.
December;65(2):361–391.

49.
Franz B, Xue C, Painter KJ, Erban R. Travelling Waves in Hybrid Chemotaxis Models
Bull Math Biol. 2014.
February;76(2):377–400. doi: 10.1007/s11538-013-9924-4
[PubMed]

Articles from PLoS Computational Biology are provided here courtesy of **Public Library of Science**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |