Limits of Feedback Control in Bacterial Chemotaxis

Inputs to signaling pathways can have complex statistics that depend on the environment and on the behavioral response to previous stimuli. Such behavioral feedback is particularly important in navigation. Successful navigation relies on proper coupling between sensors, which gather information during motion, and actuators, which control behavior. Because reorientation conditions future inputs, behavioral feedback can place sensors and actuators in an operational regime different from the resting state. How then can organisms maintain proper information transfer through the pathway while navigating diverse environments? In bacterial chemotaxis, robust performance is often attributed to the zero integral feedback control of the sensor, which guarantees that activity returns to resting state when the input remains constant. While this property provides sensitivity over a wide range of signal intensities, it remains unclear how other parameters such as adaptation rate and adapted activity affect chemotactic performance, especially when considering that the swimming behavior of the cell determines the input signal. We examine this issue using analytical models and simulations that incorporate recent experimental evidences about behavioral feedback and flagellar motor adaptation. By focusing on how sensory information carried by the response regulator is best utilized by the motor, we identify an operational regime that maximizes drift velocity along chemical concentration gradients for a wide range of environments and sensor adaptation rates. This optimal regime is outside the dynamic range of the motor response, but maximizes the contrast between run duration up and down gradients. In steep gradients, the feedback from chemotactic drift can push the system through a bifurcation. This creates a non-chemotactic state that traps cells unless the motor is allowed to adapt. Although motor adaptation helps, we find that as the strength of the feedback increases individual phenotypes cannot maintain the optimal operational regime in all environments, suggesting that diversity could be beneficial.


Introduction
Escherichia coli cells navigate their environment by alternating straight runs with direction-changing tumbles to perform a random walk. During a run, the flagellar motors spin counterclockwise (CCW) and propel the cell at constant speed in one direction, which changes slowly due to rotational diffusion. Runs are terminated when one or more motors start rotating clockwise (CW), which causes the cell to tumble [1][2][3]. Cells are able to bias their random walk toward favorable conditions using a twocomponent signal transduction pathway that detects changes in signal intensity during runs and modulates the probability to tumble accordingly, resulting in extended runs in the desired direction and net drift velocity in the direction of the gradient [4].
The sensory module of the chemotaxis pathway ( Figure 1A) consists of large clusters of receptor proteins that bind signal molecules to modulate rapidly (< 0.1s) the activity of an associated histidine kinase, CheA [5][6][7]. The high gain of the receptor cluster is coupled to negative integral feedback control [8][9][10], mediated by slow (~1-30 seconds) methylation and demethylation of the receptors by CheR and CheB, respectively [11][12][13]. This allows the receptors to adapt to a constant background signal while maintaining sensitivity over a wide range of concentrations [14,15]. For example, when cells are stimulated with a step of aspartate, the activity of the receptors returns nearly precisely to its pre-stimulus level after a transient response ( Figure 1B first line). While precise adaptation does not hold when receptors become saturated, adaptation with a precision above 80% has been measured for many relevant signals within the micromolar range [16]. Precise adaptation is an important feature of bacterial chemotaxis because it provides robustness by implementing a "maximin" strategy that guarantees at least minimum chemotactic performance in any environmental condition [17].
The activity of the sensory module is relayed through a diffusible response regulator CheY to the flagellar motors, which act as the actuator, ( Figure 1A). When phosphorylated by CheA, CheY-P binds to the motor subunit FliM and increases the probability of the mo-tor to switch from CCW to CW [18]. Fast dephosphorylation of CheY-P by the phosphatase CheZ ensures rapid transfer of information from the sensor to the actuators. Figure 1. Dynamical coupling between the sensor and the actuator in the bacterial chemotaxis system. A. The bacterial chemotaxis system is composed of a sensor module (receptorkinase complexes) and an actuator module (flagellar motors) coupled through the phosphorylated form of CheY. Both modules are ultra-sensitive and adapt to their respective input signals. Maintaining the output of the sensor within the right range relative to the actuator is critical for chemotaxis performance. B. Diagrams of the CheY-P concentration response to different signals. First line: when cells are immobilized onto a slide, a step stimulus of attractant (e.g. methylaspartate) causes a sudden decrease in CheY-P concentration followed by a slower adaptation. Because of the negative integral feedback architecture of the sensor module, CheY-P adapts back to its pre-stimulus level, the adapted CheY-P concentration, Y 0 . Second line: when immobilized cells are exposed to an exponential ramp in time of the same stimulus, the system, which is log sensing, experiences a constant "force" and adapts towards an operational CheY-P concentration, Y m , lower than the adapted level Y 0 . This deviation of CheY-P activity from Y 0 to Y m changes the coupling between sensor and actuator. Third line: when cells are swimming in a gradient of attractant, their biased random walk causes them to climb the gradient. The average drift velocity of the cell up a chemical gradient affects the average input signal experienced by the cell. This creates a feedback of the behavior onto the input signal, which in turn can significantly affect the operating concentration of CheY-P and thus the coupling between sensor and actuator.
The CW bias of the flagellar motor, which defines the tumbling probability [2], is a sensitive function of the CheY-P concentration (Hill coefficient > 10, Figure 1A) [19,20]. The capability of the system to maintain the CheY-P concentration within the tight dynamic range of the motor CW bias response function ( Figure 1A) is often used to investigate robustness to fluctuations in protein concentrations and receptor activity [21,22]. An important underlying assumption is that performance is maximized when the motor converts small variations in CheY-P into large changes in CW bias.
However, recent experiments and theory suggest that the coupling between sensor and actuator is more complex than previously thought. First, the flagellar motors partially adapt to persistent stimulus [23,24]. Second, the motor CW bias response to CheY-P is steeper than previously reported, further restricting the dynamic range of the motor response to CheY-P fluctuations [20]. Finally, in exponential ramps of chemoattractant, the CheY-P concentration reaches a dynamical equilibrium, Y m , hereafter called operational CheY-P, distinct from the adapted CheY-P concentration, Y 0 , that the cell maintains in constant uniform environments [25,26] (Figure 1B second line). For each of these three findings, the characterization of the internal dynamics of the signaling pathway was performed on immobilized cells using experimentally controlled input signals. However, for cells swimming freely in chemical gradients, the input signal dynamics are determined by the chemotactic response of the cell, creating a feedback of the behavior onto the input signal [27] (Figure 1B third line). Because of this behavioral feedback, it remains unclear how the multiple time scales of the system, from signal detection to motor response, ultimately determine chemotactic drift performance.
Here, we use analytical models and stochastic simulations of individual cells to examine the consequences of these new observations for our understanding of the bacterial chemotaxis strategy. Clonal populations of chemotactic E. coli grown in homogeneous conditions exhibit significant cell-to-cell phenotypic variability, with adaptation times ranging from 1 to 30 seconds [28][29][30][31], and motor clockwise bias ranging from 0.1 to 0.4 [3,31]. Therefore, we consider how different combinations of adaptation times and motor clockwise biases, which define a cell's phenotype, affect individual cell chemotactic drift velocity in different environments. In a phenotypically diverse population, different phenotypes of that population may perform best in different environments.
Focusing on how information transfer from sensor to actuator affects chemotactic performance, we analyze the dynamical relationship between the operational regime of CheY-P, Y m , and the drift velocity, V D , as a function of the phenotype and gradient steepness. We show that there is a unique operational regime of the sensor with respect to the motor that maximizes drift velocity in the direction of the gradient by maximizing the contrast between runs up and down the gradient, and not by maximizing the CW bias response. We characterize the performance trade-off faced by individual cells with different combinations of phenotypic parameters (such as, adapted CheY-P concentration, Y 0 , receptor adaptation time, τ, and cell resting tumble bias, TB 0 ).

Results
Maximizing contrast in run durations rather than CW bias response maximizes chemotactic drift velocity.
Previous studies have examined how E. coli chemotactic drift velocity along a one dimensional gradient depends on the adaptation time [25,32], the shape of the response function of the sensory module [17,33,34] and also behavioral feedback [27]. Instead, we first focus on how the coupling between the sensors and actuators by the response regulator CheY-P ( Figure 1A) affects the chemotactic performance of individual cell phenotypes. What CheY-P concentration maximizes drift velocity of cells navigating exponential gradients of methyl-aspartate? We examine this question using stochastic simulations of individual cells and an analytical model.
Simulations were conducted using a standard model of the chemotaxis pathway in individual cells [2,15] as described in Methods. Receptor-kinase complex activity is modeled as an all-or-none response using quasi-equilibrium dynamics for fast ligand binding, chemoreceptor conformational changes, and phosphorylation cascade [15]. The slower (de)methylation kinetics follow simple negative integral feedback dynamics with adaptation rate τ -1 . The flagellar motor is modeled as an inhomogeneous Poisson process that switches cell behavior between runs and tumbles with rates defined as a function of CheY-P concentration that varies in time, Y(t). The parameters of the motor model are calibrated to recent experimental measurements [19,20,23,24]. While motor adaptation [23] is not included at first, its effects are analyzed later in the paper. During runs, a cell swims with constant speed v = 20 µm -1 in a direction subjected to rotational diffusion (rotation diffusion constant, D r = 0.062 rad 2 s -1 [1]). For simplicity the effects of multiple flagellar motors [2,35] or directional persistence [1] are not included but discussed in the Discussion section. Hence, in this model, motor clockwise bias (CW) and cell tumble bias (TB) are the same. We consider cells containing only Tar receptors and use methyl-aspartate as the ligand. Our results readily extend to more complex receptor cluster configurations.
Three-dimensional trajectories of individual cells were simulated as described in [2] for various cell phenotypes, which are characterized by the receptor adaptation times τ and adapted CheY-P concentrations Y 0 , in gradients of chemoattractant of different steepness, g. Following previous studies [25][26][27], we used exponential gradients (L(x) = L 0 e gx ) so that cells experience an approximately constant "force" from the attractant field, as the chemotaxis system is a fold-change detector (Eq. (2) below). This makes it possible to define a steady-state drift velocity, making the problem analytically tractable. The performance of each cell phenotype, which is defined by a unique adapted CheY-P concentration (Y 0 ) and receptor adaptation time (τ), in each gradient steepness (g) is defined as the drift velocity V D (Y 0 , τ, g) along the gradient direction calculated by averaging the velocity of 10,000 phenotypically identical cells over 4 minutes (Methods). The first simulations were done in a relatively shallow gradient g -1 = 5,000 µm, with adaptation times of τ = 5, 10, and 30 s, and adapted CheY-P concentrations spanning the range Y 0 = 1-4 µM.
Plotting drift velocity as a function of the adapted CheY-P concentration reveals that maximal drift velocity is not achieved for CheY-P concentrations in the linear range of the CW bias response curve, where fluctuations in CheY-P result in large changes in clockwise bias ( Figure 2A). Instead, it occurs when adapted CheY-P is at the lower end of this curve (around 2.4 µM in Figure 2A). To understand the underlying reasons of this result, we derived an analytical relationship between CheY-P concentration and drift velocity along a onedimensional gradient. For simplicity we used a onedimensional analytical representation of bacterial chemotaxis in two or three dimensions [27,33,34,37,38]. In this framework, cells either go up or down the gradient or tumble and the effect of rotational diffusion can be represented as a jump process between runs up and runs down the gradient with transition rate (d-1)D r , where d represents the number of spatial dimension [37].
At quasi-steady state (for time scales longer than single run durations) and with no directional persistence (equal probability to run up or down the gradient), the drift velocity is proportional to the cell swimming speed v times the difference between the expected run durations up The only difference between d = 2 or d = 3 dimensions is a rescaling of the drift velocity and rotational diffusion (factor d in the equations above; see SI text).
The expected run duration up or down the gradient is controlled by the cellular concentration of CheY-P, Y. This quantity is in turn a function of the free energy difference between the inactive and active receptor complexes, F, such that F=ln(α/Y-1), where α is the gain of the phosphorylation cascade. The receptor activity follows simple spring-like dynamics around the adapted free energy difference F 0 with adaptation time τ (Methods, our Results still hold when considering asymmetric methylation/demethylation rates, see SI text and Figure S1): Here s = ± 1 or 0 for cells running up, down the gradient, or tumbling, respectively. As the cell moves along a trajectory x(t) it encounters different concentrations of the ligand L(x(t)).
the magnitude of the change in free energy difference and depends on the local steepness of the gradient at the cell position. Here, v is the speed of the cell when it is swimming, N is the gain of the cooperative receptor system and K i and K a are the dissociation constants between ligand and receptors in the inactive and active conformation. In general, both s and f change as a function of time. For ligand concentra- . If in addition, the gradient of ligand is exponential, L=L 0 e gx , we see that f ≈ vNg becomes constant where g represents the inverse length scale of the gradient. Therefore, in an exponential gradient the free energy difference of the receptors, F, tends to increase at the constant positive rate +vNg when the cell swims up the gradient and to decrease at the constant negative rate -vNg when the cell swims down the gradient. This in turn causes the CheY-P concentration to decrease (increase) when cells swim up (down) the gradient. The exact CheY-P concentration trajectories can be calculated by integrating Eq.
(2) as a function of time for different initial condition F i ( Figure 2B), while a cell is swimming up (s=1 green curve) or down (s=-1 red curve) a gradient.
The expected durations of a run, λ R −1 (Y ) , or a tumble, λ T −1 (Y ) , are plotted as a function of CheY-P concentration in Figure 2B (dashed lines; see definition in Methods). When a cell runs up or down the gradient, the rates of switching from one state to another change as a function of time and the direction of motion because they depend on the CheY-P trajectory. A run up or down the gradient can also be terminated by random reorientation from rotational diffusion with rate (d-1)D r [37]. Altogether, the rate at which a run is terminated by either rotational diffusion or a tumble is thus In a shallow gradient, F deviates little from the adapted value F 0 and the adapted value of τ R provides a good approximation of the expected run duration along a direction (black line in Figure 2B): When swimming up or down the gradient, CheY-P fluctuates ( Figure 2B, green lines for up, red lines for down) and the run lengths are modulated approximately following τ R0 (black line and red and green circles in Figure 2B). According to equation (1), drift velocity is largest where the contrast between run duration up and down the gradient is the largest. Figure 2B reveals that this is the case where the slope of the expected run length as a function of CheY-P concentration is largest, which corresponds to the foot of the motor CW bias curve ( Figure 2A) in agreement with the simulations. In contrast, for higher valued of CheY-P that are within the dynamic range of the CW bias response function, (e.g. Y 0 = 3 µM in Figure 2B) run durations up and down the gradient have a smaller contrast and longer tumble duration (dashed line Figure 2B), resulting in slower drift velocity.
In the limit of shallow gradients, Equation [1] can be linearized around the adapted values F 0 and Y 0 to obtain the drift velocity (Methods and SI Text): Here, TB 0 = λ R0 / (λ T 0 + λ R0 ) represents the tumble bias of the cell as a function of the adapted CheY-P concentration Y 0 and the subscript 0 indicates that the rates λ R and λ T and τ R ′ = dτ R / dF are all evaluated at the adapted state. The integral in Eq. (3) is the time-averaged input over the run durations, which in this approximation are exponentially distributed with characteristic time scale τ R0 . For K i << L <<K a and exponential gradients, the rate of change of the free energy difference s f svNg ≈ is constant during a run. Equation (3) indicates that the drift velocity is proportional to the gradient steepness g and the gain of the receptor cluster N. From Equation (3) we also obtain the chemotaxis coefficient of an individual Plotting the drift velocity as a function of Y 0 on top of the simulation results in Figure 2A shows that Equation (3) provides a good prediction of the drift velocity in shallow gradients and confirms that maximum velocity is reached for CheY-P values at the foot of the CW bias response curve.
In this linear regime, the optimal CheY-P concentration is only weakly dependent on the cell adaptation time and does not depend on the gradient steepness.
Because rotational diffusion imposes an upper bound on the run length along a given direction [27,32], it determines, along with the motor parameters, the optimal range for CheY-P fluctuations. As rotational diffusion becomes smaller, cells are able to maintain their original direction for a longer time. The upper bound on the run length therefore becomes longer (dashed lines in Figure 2C) and the optimal CheY-P concentration becomes smaller (full lines in Figure  2C).
Changes in the switching frequency of the flagellar motor, ω, (see Methods) also affects the optimal CheY-P concentration. This becomes apparent when considering that the rate of switching from run to tumbles scales linearly with the switching rate of the flagellar motor, ω (see Methods). Therefore, the expected run duration of a cell scales like the inverse of ω. The result of this scaling is that for increasing values of ω the inflection point of the expected run length as a function of CheY-P shifts to lower values of CheY-P (dashed lines in Figure 2D). Thus, increasing the switching rate of the flagellar motors tends to decrease the optimal CheY-P concentration (full lines in Figure 2D). It also increases the maximum drift velocity that can be reached.
The analytical model of drift velocity in Eq. (3) is different from previous results [27] in two ways. First, it takes into account both the adaptation time and the tumbling state of the cell. Taking the limit τ → ∞ , 0 0 TB ≈ in our model we recover the previous results. Second, in the previous study, the switching rate of the motor λ R was a steep function of kinase activity centered at the adapted kinase activity level, or equivalently, at the adapted CheY-P concentration Y 0 . This means that changing Y 0 would also change the set point of the motor. However, the adapted CheY-P concentration and the set point of the motor are independent parameters. For this reason here we focused on the relationship between the set point of the motor and CheY-P activity. According to Eq. (3) the flagellar motor has its own sensitivity set point independent of the adapted CheY-P concentration of the sensory system Y 0 (see definition of λ R in Methods; motor adaptation [23,39] is considered below).

A unique operational CheY-P concentration maximizes drift velocity for multiple gradients and adaptation times.
Experiments have shown that when immobilized cells are exposed to an exponential ramp of methylaspartate, CheY-P activity reaches a new steadystate, Y m , which is lower than its adapted activity, Y 0 , because of the relatively slow adaptation rate of the system [9,26] ( Figure 1B second line). When cells are swimming in an exponential gradient ( Figure 1B third line), we expected a similar effect to take place because the average drift of an individual cell up the gradient will cause this cell to experience, on average, an exponential increase in ligand concentration as it makes its way up the gradient. While this effect should be minimal in a shallow gradient, it could be-come important in steep or rapidly changing gradients [27,40], especially for cells with longer adaptation times.
To investigate this issue we simulated cells swimming in a steeper exponential gradient (g -1 = 1,000 µm). After less than one minute of simulation, cell populations (10,000 replicate trajectories for each phenotype) reached a constant steady state drift velocity. We calculated the average ligand concentration that the cells encountered over time ( Figure 3A). This reveals that the swimming cells experience an average exponential increase in ligand concentration over time. This average input is similar to the signal experienced by immobilized cells exposed to temporal exponential ramps ( Figure 1B, second line) [9,26]. However, for the swimming cells the ramp rate is dynamically determined by the average drift velocity in the direction of the gradient ( Figure 3A). Thus, for swimming cells the ramp rate depends on the feedback of the performance onto the input signal ( Figure 1B, third line). Consistent with experimental results obtained with immobilized cells exposed to exponential ramps [26], the average CheY-P concentration in the swimming cells reaches a stable dynamical equilibrium, the operational value Y m , after an initial drop from the adapted CheY-P concentration (Y 0 ) ( Figure 3B).
The fact that the operational CheY-P concentration is not the same as the adapted CheY-P concentration implies that an optimal choice of adapted CheY-P must take into account this behavioral feedback. While a phenotype may for example have an adapted CheY-P concentration equal to the optimal concentration (~2.4 µM), during chemotaxis this level drops to an operational level lower than the optimum, hindering its performance ( Figure 3B, solid black line). This effect is intensified when the adaptation time of the receptor cluster increases ( Figure 3B, grey line). On the other hand, a phenotype with an adapted CheY-P concentration higher than the optimal concentration can approach the optimal operational CheY-P concentration as it reaches its steady-state drift velocity ( Figure 3B, dotted line). The difference between Y 0 and Y m grows larger as drift velocity or the receptor adaptation time increase ( Figure 3C). . Y m is instantaneous CheY-P concentration averaged over the population while drifting between t = 60 and 300 s). Two exponential gradients of methyl-aspartate are considered (g -1 = 1,000 µm (black), 5,000 µm (grey)). Black arrow: cell population in blue in Figure 4C.
To determine whether the adapted or the operational CheY-P concentration is the primary variable that controls the average drift velocity in exponential gradients, we simulated cell populations with different adapted CheY-P concentrations and calculated their respective operational CheY-P concentrations. In a steep gradient (g -1 = 1000 µm), the optimal adapted CheY-P increased to ~2.7 µM compared to ~2.4 µM in a shallow gradient ( Figure 3D). However, the optimal operational CheY-P concentrations for steep and shallow gradients are identical ( Figure   3D). This suggests that a unique operational CheY-P concentration maximizes drift velocity in multiple gradients.  (2) around F m and solve for quasi-steady state: (4) This expression quantifies the deviation between the operational free energy difference F m and the adapted free energy difference F 0 as a function of the drift velocity and is consistent with the results of our simulations ( Figure 3C) and [27]. Equation (4) also makes clear that the strength of the feedback depends on adaptation time, the receptor cluster gain N, and the steepness of the gradient. Behavioral feedback strongly affects performance because it moves Y m away from the optimal operating point relative to the motor. This, in turn, affects the capability of the motors to best use the information carried by CheY-P.
By explicitly taking into account the effect of the behavioral feedback onto the coupling between the operating regime of CheY-P and the motor, Eqs. (3) (with 0→m ) and (4) extend previous studies [17,25,27,[32][33][34]37] and reveal new possible dynamical regimes for the biased random strategy as shown below.
Strong behavioral feedback can push the system through a bifurcation creating two possible chemotactic states for some cell phenotype: a fast drift state and a trapped state.
For a given phenotype (Y 0 , τ) and gradient lengthscale, the steady state drift velocity is determined by the intersection of two curves ( Figure 4A). The first curve (solid line in Figure 4A) describes how the drift velocity depends on the operational CheY-P concentration, Y m . It is defined by Equation (3) (with 0→m ) and its profile can be interpreted as follows. For very low values of CheY-P the cell never tumbles. Thus, the cell diffuses equally in all directions and the net drift along the gradient is zero. For high values of CheY-P, the cell tumbles all the time so drift is zero as well. In between these two extremes, drift velocity is maximized for a specific value of the operational CheY-P concentration. However, Y m is not an independent variable. As we showed above (Eq. (4)), because of the feedback the behavior onto the input, the operational CheY-P concentration is itself a function of the drift velocity (which can also be written as: ) ) ). This equation defines the dashed line in Figure 4A, which intersects the horizontal axis at Y 0 . Because each line in Figure 4A defines a relationship between V D and Y m , the intersection between the two lines fully determines the drift velocity and the operational CheY-P concentration (black circle in Figure 4A) for a given phenotype and gradient.
When the feedback is weak (τNg small, i.e. short adaptation time, small gain, or shallow gradient), the operational CheY-P concentration only exhibits a weak dependency on drift velocity and there is only one steady-state solution (intersection). Therefore, an appropriate adapted CheY-P concentration could be selected to ensure that operational CheY-P concentration is approximately optimal at all times ( Figure  4A).
When the feedback is stronger, drift velocity always acts as a negative feedback onto the operational CheY-P concentration. In contrast, the effect of the operational CheY-P concentration onto drift velocity depends on whether the operational CheY-P concen- tration is below or above the CheY-P concentration that maximizes drift velocity ( Figure 4A). Below this concentration, the system obeys negative feedback dynamics, whereas above it, the system obeys positive feedback dynamics. This positive feedback loops combined with the non-linear decrease of the drift velocity as a function of the operational CheY-P concentration, which arise from the extreme sensitivity of the flagellar motor, can lead to bistability [41]. Indeed, for a stronger feedback (steeper gradient or longer adaptation time) the slope of the feedback curve (dashed line in Figure 4AB), which is proportional to 1/τNg, decreases. Thus, for phenotypes with high enough adapted CheY-P concentration (Y 0 is the intersection of the dashed line with the horizontal ax-is), the two curves can intersect more than once ( Figure 4B). In this case, a single phenotype can now experience three different chemotactic states. Two of these states, Y m1 and Y m2 (filled circles in Figure 4B), are stable and are separated by one unstable state (open circle in Figure 4B). For one of the stable solution, the drift velocity is nearly zero and Y m1 is high and very close to the adapted CheY-P concentration.
For the other stable solution, the drift velocity is large and Y m2 is much smaller than the adapted CheY-P concentration.
This analysis suggests that an individual phenotype might experience two different chemotactic states with dramatically different performance: a fast drifting state and a "trapped" state. To find evidence of these two behaviors, we simulated two cell phenotypes (10,000 replicates for each phenotype) in a steep exponential gradient (g -1 = 1,000 µm). One phenotype was predicted to operate closer to the bifurcation than the other (red and blue dots in Figure  4C, respectively). Although both phenotypes reached the same average operational CheY-P (Y m = 2.3 µM), cells with a phenotype closer to the predicted bifurcation point (Y 0 = 3µM, τ = 30s) exhibited a distribution of behavior (both drift velocity and diffusion) significantly skewed toward the "trapped" state (Figure 4C). Closer examination of the trajectories and CheY-P dynamics of individual cells in this simulation reveals that individual cells transition stochastically back and forth between the "trapped" and fast drifting state ( Figure S2). For cases with higher feedback strength cells spend more and more time within the "trapped" state.
When the feedback is strong and the system becomes multistable, the average includes cells in both the "trapped" and high drift states. This phenomenon explains the decreased average drift velocities observed when adaptation time is increased (above 10 seconds) in a relatively steep gradient (g -1 = 1,000 µm). It also explains the resulting shift of the best operational CheY-P concentration to lower concentrations (from ~2.4 to 2.1µM in Figure 4D), since for phenotypes with lower values of the adapted CheY-P only one stable state exists. Similar results are obtained when asymmetric methylation/demethylation rates of the receptors are taken into account ( Figure  S3).

Motor adaptation partially alleviates the chemotactic "trap".
Recent experiments have shown that the number of FliM monomers in the C-ring of the flagellar motor slowly (~minutes) adapts as a function of the CW bias, affecting both the steepness and the halfmaximum CheY-P concentration of the CW bias motor response curve [24]. To examine the effect of motor adaptation on the relationship between CheY-P concentrations and drift velocity, we added motor adaptation to our stochastic model of an individual chemotactic cell by taking into account recent experimental data [20,23,24] (Methods). The resulting CW bias response curve of the adapted motor agrees well with both recent [20] and earlier [19] experimental measurements. In fact, it matches earlier experiments [19] better than a simple Hill function, suggesting that in these experiments the individual motors measured had adapted to the particular concentration of CheY-P expressed in the corresponding individual cells ( Figure 5A; Methods).
Simulations of cells with motor adaptation in a shallow gradient (g -1 = 5,000 µm) show that motor adaptation changes the shape of the drift velocity curve as a function of operational CheY-P, especially at high CheY-P concentrations (compare Figures 5B and  2A). These results are predicted by the analytical model (Eq. (3) with 0→m ) once modified to include motor adaptation (Methods; lines in Figure 5B). Setting the adapted activity of the motor for a given CheY-P concentration to lower or higher CW biases gives qualitatively equivalent results ( Figure S4 and S5).
How does the motor adaptation affect the bifurcation? In a steep gradient the behavioral feedback (Eq. (4)) must be taken into account ( Figure 5C dashed line). Comparing Figure 5C and 4B we see that motor adaptation enable cells with high adapted CheY-P concentration to avoid the chemotactic trap improving performance (see Figure S6). This should provide a selective advantage because it helps buffer the functional consequences of inevitable cell-to-cell variability in the adapted CheY-P concentration, by increasing the range of CheY expression levels that allows effective chemotaxis.
Motor adaptation also affects the optimal operational CheY-P concentration (compare Figures 5B and 2A), shifting it to lower concentrations. When cells are drifting up a gradient, CheY-P drops to the operational CheY-P, causing the CW bias to drop. With motor adaptation, the lower operating CW bias causes the motor to shift its sensitivity to a lower CheY-P concentration. We see again that maximal performance is reached for Y m at the bottom of the CW bias response curve of the now adapted motor ( Figure  5A). However, the motor can only compensate partially for the shift in operational CheY-P concentration.

An individual phenotype faces a performance trade-off in different gradients.
As long as the system does not undergo bifurcation, maximum drift velocity is achieved by having a long adaptation time while maintaining the operational concentration of CheY-P in the optimal range. Therefore, the optimal adapted CheY-P concentration depends on the gradient length-scale and the adaptation time ( Figure 6).
In shallow gradients, the strength of the feedback is small, as is the difference between operational and adapted CheY-P. Thus, it is possible to select an adapted CheY-P concentration that will perform relatively well for multiple adaptation times ( Figure 6A blue line).In steeper gradients, the feedback is stronger (Eq. (4)) and the difference between Y m and Y 0 grows larger with adaptation time. Maintaining the optimal operational CheY-P concentration requires a higher adapted CheY-P concentration (Figure 5A green and red). The bifurcation of the system imposes an upper bound on the range of Y 0 beyond which a portion of the cells spend a significant amount of time trapped into a non-optimal state even with motor adaptation (Figure 6A dashed lines).  (5)) as a function of the chemoreceptor adaptation time in different exponential gradients (g -1 = 1,000 (red), 2,000 (green), and 5,000 (blue) µm). Dots indicate when the maximal theoretical drift velocities cross the bifurcation point (dotted lines represent the inaccessible optimal state). The optimal operational CheY-P concentration Y m is identical for all gradient length scales (black dashed line). B. Contour plot of drift velocities as a function of adaptation time and the adapted cell tumble bias in different exponential gradients (same colors as A). 75%, 90%, and 95% contours of the maximal theoretical drift velocities for each gradient (colors intensities from light to dark). Black dot: the best cell phenotype that achieves equal relative drift velocities in all three gradients (60% of the maximal V D with τ = 7.5 s and TB 0 = 0.044).
Therefore, the optimal adapted CheY-P concentration is a function of both receptor adaptation time and gradient length-scale, making it difficult for a single phenotype to maximize drift velocity in multiple environments ( Figure 6A).
To characterize the resulting performance trade-off and map it to phenotypic space, we calculated the contours of drift velocity relative to its maximum in each environment, as a function of adaptation time τ and the adapted cell tumble bias ( Figure 6B). In shallow gradients, cells benefit from a relatively long adaptation time and a low adapted CW bias. In steep gradients, cells benefit from a short adaptation time and a higher adapted CW bias. The best generalist phenotype can achieve at most 60% relative performance in all three gradients considered here. Motor adaptation, which was taken into account in generating Figure 6, alleviates only partially the tradeoff faced by single cells.

Discussion
The adaptive response and feedback control of the receptor cluster play a critical role in the robustness of the chemotaxis system [8,10,15,17,33]. However, chemotactic performance also relies on the optimal operation of the flagellar motors, which directly control cell behavior. By focusing on how the CheY-P concentration affects the coupling between sensors and actuators, we revealed the existence of an operational regime for CheY-P concentration, which is distinct from the adapted CheY-P concentration, that maximizes drift velocity in a wide range of gradient length-scales and receptor adaptation times. Fluctuations around the best operational CheY-P concentrations maximize the contrast between run duration up and down the gradients. This occurs outside the most sensitive region of the CW bias response curve of the motor. Thus, chemotactic performance relies on maintaining the operational CheY-P concentration within bounds [21,22,42] around this optimal value.
The best operational CheY-P concentration is also determined by the cell rotational diffusion constant D r , which imposes an upper bound on run durations in any particular direction ( Figure 2C) [27,32]. In a more viscous environment or for longer cells, the lower rotational diffusion will result in a lower optimal operational CheY-P concentration. For an ellipsoid, rotational diffusion is inversely proportional the length of the major axis. Therefore, as cells grow the optimal range will shift to lower CheY-P concentrations. If the cell maintains a constant amount of CheY as it grew, the effective concentration of CheY-P would decrease, resulting in robust chemotactic performance during cell growth.
The switching frequency of the flagellar motors also affects the best operational CheY-P concentration.
Higher switching frequencies tend to increase drift velocity while shifting the maximum to smaller CheY-P concentrations ( Figure 2D). Therefore, the best operational CheY-P concentration is further away from the motor threshold. However, the range of CheY-P concentrations where the drift velocity is high becomes narrower (because the expected run length becomes a steeper function of CheY-P). This tends to increase the performance trade-off between different gradient length-scales. Thus, while selecting a higher switching frequency for the flagellar motors may improve performance of some phenotypes it may be detrimental for the population overall. Another important consideration is that the switching frequency is bounded by the speed at which the motor and associated flagella can switch confirmation [2,43].
Directional persistence (amount by which the swimming direction of a new run is biased towards the swimming direction of previous run) has been shown to affect chemotactic performance in climbing shallow gradients of attractants [1,36,44,45]. However, previous modeling and simulations efforts have been done using cells with non-optimal CheY-P concentrations (usually at 3 µM). In this regime, cells have a high tumbling rate, short run lengths, and low drift velocity. Directional persistence effectively reduces the reorientation rate of cells [45], which is equivalent to reducing the tumbling rate slightly. Therefore, directional persistence will shift the optimal CheY-P concentration to higher concentrations and improve the drift velocity of frequently tumbling cells [44]. On the other hand, when cells operate at or close to the optimal CheY-P concentration, the tumbling rate is low. Therefore, the run length in a given direction is terminated by rotational diffusion and not by tumbles. For optimal phenotypes, the relative effect of directional persistence on chemotactic drift is thus less important.
Previous studies have examined how the adaptation time affects chemotactic performance [12,25,27,32,46]. However, these studies only considered single values for the adapted CheY-P concentration (typically set to a CW bias of 0.5) and concluded that adaptation time should decrease as gradients get steeper to keep the operational CheY-P concentration within the dynamic range of the motor CW bias response. We found that, as long as the cell can maintain the optimal operational CheY-P concentration, longer adaptation time is better because it enhances input signal over the course of a run. However, long adaptation reinforces the feedback from the cell drift velocity on the system and can lead to undesirable bistability. Therefore, the bifurcation boundary imposes an upper limit on adaptation time as a function of the gradient length-scale. Interestingly, the distribution of tumble bias typically observed during exponential growth in single E. coli cells ranges from 0.1 to about 0.4 and not many cells are found that have higher tumble bias [31]. Selection for cells with tumble bias below 0.4 is consistent with our finding that the performance of cells with higher tumble bias will suffer from the existence of the "trapped" chemotaxis state.
Our results also provide a strong justification for the role of the recently-discovered flagellar motor adaptation. Indeed, we found that motor adaptation [23,39] plays a significant role in mitigating the behavioral feedback for cells with high tumble bias. When such feedback was included, cells with high tumble bias could escape the "trap" and gain access to a high drifting state in steep gradient. Our model also resolved an apparent contradiction between two sets of experimental measurements of the CW bias response of the flagellar motor as a function of CheY-P concentration. While one measurements reported a Hill coefficient of n=10 [19], newer experiments reported a Hill coefficient of n=20 [20]. In this paper we used the new value n=20 and showed that the previous measurements are fitted with the same parameter value if one makes the reasonable assumption that the motors had had time to adapt before each individual cell measurement ( Figure 5A).
Because the difference between the operational and adapted CheY-P concentrations depends on the strength of the behavioral feedback, which itself is proportional to gradient steepness, different adapted CheY-P concentrations and adaptation times are required to perform optimally in different gradients. Thus, in conditions where drift velocity is important, cells are faced with a performance trade-off. Even though motor adaptation was included, the best compromising phenotype over the gradient steepness considered in this study achieved at most 60% of the theoretical maximal drift velocity in all gradients. The observed cell-to-cell phenotypic diversity in adaptation time and adapted tumble bias [29,31] in an isogenic population may resolve the performance trade-off faced by single cells to improve the chance of survival of a unique genotype in complex or varying environments. In addition, the negative correlation between tumble bias and adaptation time observed by Park et al. in an isogenic population of E. coli [31], is consistent with our predictions about the most beneficial way to distribute phenotypes ( Figure  6B). At its core, the biased random walk relies on the dynamical control of the probability of reorientation. Overall, our analysis reveals limits to the use of negative integral feedback to control such strategy. Because the biased random walk strategy is used by many organisms, these results will inform our understanding of the constraints faced by other organisms as well.

Model and simulations
We used a standard model of bacterial chemotaxis [15] as described in [2]. For a cell following the trajectory x(t), the output of the sensory module, the CheY-P concentration, is ( ) ( )  [26]. Cells switch between runs, R, and tumble, T, with The motor is modeled as a bistable system with switching frequency ω = 1.3 s -1 (unless otherwise stated) and free energy difference ( ) ε and 3 ε are non-dimensional constants that control the basal rate of switching of the motor when Y=0 and the degree of cooperativity of the motor, respectively. K is the binding constant of CheY-P to FliM at the base of the motor. With ε 2 = ε 3 = 80, and K = 3.06 µM, this coarse-grained motor model fits well recent experimental measurements of CW bias (Hill coefficient 20) and switching frequency [20,23,24]. Motor adaptation is considered below.

Linear expansion
Eq. (2) follows by taking the time derivative of F and using the relations from the previous section. Integration of Eq.
(2) gives: The expected duration of a run along the direction s = ±1 is determined by the integral of the rate τ R −1 of terminating a run along the direction s by tumbling or because of rotational diffusion: Because the average cell drift velocity in the direction of the gradient is determined by the contrast between expected run durations up and down the gradient (Equation (1)), the quantity of interest to calculate from Equation (5) In a shallow gradient, the deviations from the adapted free energy difference (1) and using the solution F(t, s, F i ) we obtain the drift velocity to first order in F Δ (Eq. (3)).

Motor adaptation
The number of FliM molecules in the motor, n, is modeled as a binding and unbinding process with CW bias dependent rates [39]: The constants k on and k off define the rate of adaptation of the motor. n 1 and n 2 are the minimum and maximum FliM ring size that a motor can accommodate. n Δ is an effective half max parameter that guarantees that the effective rates of unbinding and binding to the motor go to zero when n approaches n 1 or n 2 . When n changes it affects the steepness of the motor CW bias response, , which in our case is controlled by ε 3 (see above). We used a simple linear relationship ( ) where ε 3,1 is the slope and n 0 and ε 3,0 are the pre-stimuli level of the number of FliM and motor steepness, respectively. k off = 0.025 s -1 , n 1 = 34, n 1 = 44 from [24]. We choose ε 3,0 = 80, n 0 = 36 to match the Hill coefficient of 20 measured for individual motor response curves [20], and fit ∆n = 2.74, ε 3,1 = 2.31 to reproduce [19] ( Figure 5A). k on = 0.0063 s -1 controls the CW bias that the motor adapts to (0.2 in this case, typical for wild type population of E. coli selected for swimming on agar plates [31]). At steady state, dn/dt = 0 defines CW(n(ε 3 )) (Eq. [S20] in SI text). On the other hand, assuming quasi-equilibrium between the motor and operational CheY-P concentration Y m , we  At steady state, we linearize Eq. (S1) around the mean free energy F m , which is on average lower than the resting free energy F 0 . , The variation in f(t) depends on the direction of motion The above relation shows that as long as which is the same as Eq. (3) in the main text (after the subscript 0 has been replaced by m).

Nonlinear solution
All analytical curves in the paper (lines in Figures 2-5) use the linear approximation around F 0 and F m as described in the main text and Materials and Methods. Here we describe how to solve Eqs. (1)-(3) keeping the nonlinearity of the rates λ R (F) and λ T (F). Eq. (S16) can be integrated numerically to calculate the stopping time (red and green circles) of the red and green trajectories in Figure 2B. All other analytical curves in the paper (lines in Figures 2-6) use the linear approximation around F 0 and F m as described in Materials and Methods.
In the 1D representation, the equation for the free energy difference F can be integrated to get Inverting we also get the time duration as a function of the free energy: At steady state the conditional probability densities of the duration t of runs and tumbles are where s = ±1 and the first two probability densities correspond to runs that terminate into a tumble and into a run of opposite direction, respectively. F iR and F iT are the values of F at the beginning of a run and a tumble, respectively. Noting that The probability density to have free energy F e at the end of a run and tumble cycle is then

Effect of asymmetric methylation/demethylation rates
Experimental data shows the methylation/demethylation rates for receptor adaptation are asymmetric [2]. The rate of change of methylation catalyzed by CheR and CheB is usually described as where V R and V B (a) are the rates of methylation and demethylation; and K R and K B are the constants for each reactions. Experimentally, people found that the asymmetry of methylation/demethylation rates is not significant until a>a B , where 0.78 B a ≈ measured in [2]. V B (a) is a piece-wise linear function: For most of the dynamic range of CheY-P level we are interested in ( ) B a a Θ − will be zero and V B will be approximately constant. Thus variations in the rate of demethylation should not affect much drift velocity and optimal CheY-P level. To verify this, we implemented Eq. (S19) into our stochastic simulations of individual cells, with K R =0.43, K B =0.3 and k B =2.7 [2]. To vary the adapted CheY-P level we varied V R and V B,0 since their ratio determines the adapted CheY-P level. Given these definitions the effective adaptation time scale τ eff obtained by linearizing equation (S19) reads where a 0 is the adapted activity of the receptor (corresponding to adapted CheY-P level Y 0 in this case) . As shown in Fig. S1, the optimal CheY-P level in shallow gradient remains at the same position with respect to the motor response curve as in Fig. 2A. While the cells drifts in the steep gradient with slow methylation and demethylation rates, the behavior feedback will still push the system to a bifurcation (Fig. S2).

Adapted motor response curve
The motor adaptation is considered in this study by assuming that the number of FliM molecules in the motor changes as a function of the CW bias of the motor. At steady state, where dn/dt = 0, the relation between CW bias and number of FliM, n, is given by most visible for Y m between 2.5 and 3.5 µM (Fig. S3) whereas it is between 2 and 3 µM when k on = 0.0063 s -1 (Fig. 4A). We also simulated the case where the CW bias that the motor adapts to is 0.05 (here k on = 0.0013 s -1 ) which results in a flat region of the drift velocity curve as a function of Y m around the optimal operational Supporting Figures Figure S1. Effect of asymmetric methylation/demethylation rates on drift velocity V D in exponential gradient. Simulated drift velocity V D (average velocity of 10,000 cells between t = 60 and 300 s) as a function of operational CheY-P concentration Y m in a shallow gradient (L 0 = 200 µM and g -1 = 5,000 µm) for cells with methylation rates V R = 0.1 s -1 (Red), 0.2 s -1 (Green), and 0.4 s -1 (Blue). Figure S2. A simulated cell can transition in and out of the non-chemotactic state to reach the high drift velocity state when swimming in a steep gradient of methyl-aspartate (g -1 = 1,000 µm) illustrating the bi-stable behavior of this cell phenotype (Y 0 = 3.0 µM and τ = 30 s). A. Single cell drift velocity as a function of its operational CheY-P concentration. When the cell escapes the "trapped" chemotactic state, characterized by a high CheY-P concentration, the behavioral feedback maintains an optimal CheY-P concentration and a high drift velocity. B. Cell position along the gradient as a function of its operational CheY-P concentration. The cell can escape the low drift velocity state and maintain a low CheY-P concentration when running up the gradient. On the other hand the cell can return to the "trapped" state after a long run down the gradient. The CheY-P concentration and drift velocity were calculated over a moving average window of 10 seconds. The time progression along the trajectory is indicated by the color of the stroke from blue to red. Figure S3. Effect of asymmetric methylation/demethylation rates on drift velocity V D in steep exponential gradient. V D from stochastic simulations (methylation rate V R = 0.1s -1 ) as a function of Y 0 (filled circles) and Y m (open circles) in exponential gradient of methyl-aspartate (g -1 = 1,000 µm). Figure S4. Drift velocity as a function of operational CheY-P when the rate of binding between FliM and the motor is k on = 0.025 s -1 . Circles are from simulations. Lines are from analytical solution. Everything is the same as in Figure 5A. The only difference is the value of k on . Figure S5. Drift velocity as a function of operational CheY-P when the rate of binding between FliM and the motor is k on =0.0013 s -1 . In this case, the motor does not adapt fast enough to reach quasi-steady state. The analytical solution (lines) makes the approximation that the system is at steady state. Figure S6. Scatter plot of individual drift velocities (in the direction of the gradient) and root mean square displacements (perpendicular to the gradient) of 10,000 simulated cells with motor adaptation (Black) and without motor adaptation (Red), with adaptation time τ = 30 s adapted CheY-P concentration Y 0 = 3.5 µM, gradient length scale g -1 = 1,000 µm.