Open Access

A modular gradient‐sensing network for chemotaxis in Escherichia coli revealed by responses to time‐varying stimuli

Thomas S Shimizu, Yuhai Tu, Howard C Berg

Author Affiliations

  1. Thomas S Shimizu1,,
  2. Yuhai Tu2 and
  3. Howard C Berg*,1
  1. 1 Department of Molecular and Cellular Biology, Harvard University, Cambridge, MA, USA
  2. 2 T. J. Watson Research Center, IBM, Yorktown Heights, NY, USA
  1. *Corresponding author. Department of Molecular and Cellular Biology, Harvard University, 16 Divinity Avenue, Cambridge, MA 02138, USA. Tel.: +1 617 495 0924; Fax: +1 617 496 1114; E-mail: hberg{at}
  • Present address: FOM Institute for Atomic and Molecular Physics (AMOLF), Science Park 104, 1098 XG Amsterdam, The Netherlands


The Escherichia coli chemotaxis‐signaling pathway computes time derivatives of chemoeffector concentrations. This network features modules for signal reception/amplification and robust adaptation, with sensing of chemoeffector gradients determined by the way in which these modules are coupled in vivo. We characterized these modules and their coupling by using fluorescence resonance energy transfer to measure intracellular responses to time‐varying stimuli. Receptor sensitivity was characterized by step stimuli, the gradient sensitivity by exponential ramp stimuli, and the frequency response by exponential sine‐wave stimuli. Analysis of these data revealed the structure of the feedback transfer function linking the amplification and adaptation modules. Feedback near steady state was found to be weak, consistent with strong fluctuations and slow recovery from small perturbations. Gradient sensitivity and frequency response both depended strongly on temperature. We found that time derivatives can be computed by the chemotaxis system for input frequencies below 0.006 Hz at 22°C and below 0.018 Hz at 32°C. Our results show how dynamic input–output measurements, time honored in physiology, can serve as powerful tools in deciphering cell‐signaling mechanisms.

Visual Overview


In searching for better environments, bacteria sample their surroundings by random motility, and make temporal comparisons of experienced sensory cues to bias their movement toward favorable directions (Berg and Brown, 1972). Thus, the problem of sensing spatial gradients is reduced to time‐derivative computations, carried out by a signaling pathway that is well characterized at the molecular level in Escherichia coli. Here, we study the physiology of this signal processing system in vivo by fluorescence resonance energy transfer (FRET) experiments in which live cells are stimulated by time‐varying chemoeffector signals. By measuring FRET between the active response regulator of the pathway CheY‐P and its phosphatase CheZ, each labeled with GFP variants, we obtain a readout that is directly proportional to pathway activity (Sourjik et al, 2007). We analyze the measured response functions in terms of mechanistic models of signaling, and discuss functional consequences of the observed quantitative characteristics.

Experiments are guided by a coarse‐grained modular model (Tu et al, 2008) of the sensory network (Figure 1), in which we identify two important ‘transfer functions’: one corresponding to the receptor–kinase complex, which responds to changes in input ligand concentration on a fast time scale, and another corresponding to the adaptation system, which provides negative feedback, opposing the effect of ligand on a slower time scale. For the receptor module, we calibrate an allosteric MWC‐type model of the receptor–kinase complex by FRET measurements of the ‘open‐loop’ transfer function G([L],m) using step stimuli. This calibration provides a basis for using the same FRET readout (between CheY‐P and CheZ) to further study properties of the adaptation module.

It is well known that adaptation in E. coli's chemotaxis system uses integral feedback, which guarantees exact restoration of the baseline activity after transient responses to step stimuli (Barkai and Leibler, 1997; Yi et al, 2000). However, the output of time‐derivative computations during smoothly varying stimuli depends not only on the presence of integral feedback, but also on what is being integrated. As this integrand can in general be any function of the output, we represent it by a black‐box function F(a) in our model, and set out to determine its shape by experiments with time‐varying stimuli.

We first apply exponential ramp stimuli—waveforms in which the logarithm of the stimulus level varies linearly with time, at a fixed rate r. It was shown many years ago that during such a stimulus, the kinase output of the pathway changes to a new constant value, ac that is dependent on the applied ramp rate, r (Block et al, 1983). A plot of ac versus r (Figure 5A) can thus be considered as an output of time‐derivative computations by the network, and could also be used to study the ‘gradient sensitivity’ of bacteria traveling at constant speeds.

To obtain the feedback transfer function, F(a), we apply a simple coordinate transformation, identified using our model, to the same ramp‐response data (Figure 5B). This function reveals how the temporal rate of change of the feedback signal m depends on the current output signal a. The shape of this function is analyzed using a biochemical reaction scheme, from which in vivo kinetic parameters of the feedback enzymes, CheR and CheB, are extracted. The fitted Michaelis constants for these enzymatic reactions are small compared with the steady‐state abundance of their substrates, thus indicating that these enzymes operate close to saturation in vivo. The slope of the function near steady state can be used to assess the strength of feedback, and to compute the relaxation time of the system, τm. Relaxation is found to be slow (i.e. large τm), consistent with large fluctuations about the steady‐state activity caused by the near‐saturation kinetics of the feedback enzymes (Emonet and Cluzel, 2008).

Finally, exponential sine‐wave stimuli are used to map out the system's frequency response (Figure 5C). The measured data points for both the amplitude and phase of the response are found to be in excellent agreement with model predictions based on parameters from the independently measured step and ramp responses. No curve fitting was required to obtain this agreement. Although the amplitude response as a function of frequency resembles a first‐order high‐pass filter with a well‐defined cutoff frequency, νm, we point out that the chemotaxis pathway is actually a low‐pass filter if the time derivative of the input is viewed as the input signal. In this latter perspective, νm defines an upper bound for the frequency band over which time‐derivative computations can be carried out.

The two types of measurements yield complementary information regarding time‐derivative computations by E. coli. The ramp‐responses characterize the asymptotically constant output when a temporal gradient is held fixed over extended periods. Interestingly, the ramp responses do not depend on receptor cooperativity, but only on properties of the adaptation system, and thus can be used to reveal the in vivo adaptation kinetics, even outside the linear regime of the kinase response. The frequency response is highly relevant in considering spatial searches in the real world, in which experienced gradients are not held fixed in time. The characteristic cutoff frequency νm is found by working within the linear regime of the kinase response, and depends on parameters from both modules (it increases with both cooperativity in the receptor module, and the strength of feedback in the adaptation module).

Both ramp responses and sine‐wave responses were measured at two different temperatures (22 and 32°C), and found to differ significantly. Both the slope of F(a) near steady state, from ramp experiments, and the characteristic cutoff frequency, from sine‐wave experiments, were higher by a factor of ∼3 at 32°C. Fits of the enzymatic model to F(a) suggest that temperature affects the maximal velocity (Vmax) more strongly than the Michaelis constants (Km) for CheR and CheB.

Successful application of inter‐molecular FRET in live cells using GFP variants always requires some degree of serendipity. Genetic fusions to these bulky fluorophores can impair the function of the original proteins, and even when fusions are functional, efficient FRET still requires the fused fluorophores to come within the small (<10 nm) Förster radius on interactions between the labeled proteins. Thus, when a successful FRET pair is identified, it is desirable to make the most of it. We have shown here that combined with careful temporal control of input stimuli, and appropriately calibrated models, a single FRET pair can be used to study the structure of multiple transfer functions within a signaling network.

  • Combining in vivo FRET with time‐varying stimuli, such as steps, ramps, and sinusoids allowed deduction of the molecular mechanisms underlying cellular signal processing.

  • The bacterial chemotaxis pathway can be described as a two‐module feedback circuit, the transfer functions of which we have characterized quantitatively by experiment. Model‐driven experimental design allowed the use of a single FRET pair for measurements of both transfer functions of the pathway.

  • The adaptation module's transfer function revealed that feedback near steady state is weak, consistent with high sensitivity to shallow gradients, but also strong steady‐state fluctuations in pathway output.

  • The measured response to oscillatory stimuli defines the frequency band over which the chemotaxis system can compute time derivatives.


Signaling networks are complex, but typically possess modular architectures (Hartwell et al, 1999). This observation provides hope for understanding the function of the whole in terms of simpler representations of its constituent parts. Much recent work on signaling networks has focused on assigning discrete functions to these parts; however, an essential next step is to determine the nature of the coupling between them.

Here, we have combined theoretical modeling with in vivo measurements of signaling to characterize the transfer functions of the two modules involved in computations of time derivatives by the Escherichia coli chemotaxis‐signaling pathway. The first module, which we call the receptor module, detects changes in environmental conditions to generate an intracellular signal. This module has been shown to amplify signals over a wide dynamic range (Bray, 1998; Sourjik and Berg, 2002b). The second module, which we call the adaptation module, is known to maintain the intracellular signal at a steady state that is indifferent to ambient concentrations of ligand (Berg and Tedesco, 1975; Spudich and Koshland, 1975; Barkai and Leibler, 1997; Alon et al, 1999). Both modules have been characterized extensively by experiment, and theoretical efforts to explain their behavior have enjoyed remarkable success. However, understanding how these modules combine in living cells to produce the biologically important function of gradient sensing requires additional knowledge about the way in which they are coupled in vivo.

To define this link quantitatively, we used a simple theoretical model of the chemotaxis pathway that abstracts many of the known molecular details, but preserves the essential characteristics of the receptor response and adaptation dynamics (Tu et al, 2008). We exploited an important result of this analysis that provides us with a method for inferring the dynamics of one module through input–output measurements of the other. Specifically, we used fluorescence resonance energy transfer (FRET) to monitor the real‐time activity of the sensory receptor–kinase complex, while subjecting cells to time‐varying stimuli that render the output of this receptor module constant in time. Such a stimulus was identified by Block et al (1983), who found that the output, as inferred by the switching statistics of the flagellar motor, reached a steady value during exponential ramps, that is, ligand profiles that increase or decrease exponentially in time. Under these conditions, the feedback signal must exactly cancel the change in input signal, so we can infer the dynamics in the adaptation module through measurements of the receptor module, given a well‐calibrated model of the latter and good control over the temporal profile of ligand input. Fortunately, the receptor module has been characterized in detail by a large body of in vitro experiments (Dunten and Koshland, 1991; Borkovich et al, 1992; Li and Weis, 2000; Bornhorst and Falke, 2001; Antommattei et al, 2004; Lai et al, 2005), in vivo FRET experiments (Sourjik and Berg, 2002b, 2004; Sourjik et al, 2007; Endres et al, 2008), and theoretical modeling (Shi and Duke, 1998; Duke and Bray, 1999; Mello and Tu, 2003, 2005, 2007; Shimizu et al, 2003; Mello et al, 2004; Sourjik and Berg, 2004; Keymer et al, 2006; Skoge et al, 2006). Our model uses a variant of the Monod–Wyman–Changeux (MWC) model for the receptor module, which we calibrated by measuring FRET responses to step stimuli in a series of adaptation‐deficient mutants with systematically varying receptor‐modification levels. This ‘open‐loop’ transfer function of the receptor module defines quantitatively the way in which ligand concentration, the input signal, is balanced by receptor covalent modification, the feedback signal.

This model‐driven approach provided two important advantages. First, we were able to measure the downstream adaptation dynamics through the same FRET method earlier developed to characterize the upstream receptor–kinase response (Sourjik et al, 2007). Second, because we measured the receptor module's output as a proxy for inferring the dynamics of the adaptation module, the measurement automatically yielded a second input–output relation, namely the response in adaptation kinetics to input from the receptor module. In short, we could elucidate the essential dynamical character of the downstream adaptation module without directly monitoring any of its molecular constituents.

In addition to ramp responses, we also measured responses to oscillatory input signals, which revealed the frequency response of the system. Together, these measurements allowed us to quantify two biologically important performance characteristics of the E. coli chemotaxis system, namely the output of time‐derivative computations, and the frequency band over which such computations can be carried out. In contrast to earlier results obtained by monitoring motor switching statistics, no discontinuities in kinase activity were observed. Data collected at two different temperatures suggest a strong dependence of gradient‐sensing performance on the ambient temperature.


We present in Figure 1A a molecular view of the chemotaxis network, and in Figure 1B, the simplified view adopted in our model. At the molecular level, the intracellular signal pathway comprises multiple species of transmembrane receptors (also known as methyl‐accepting chemotaxis proteins), a scaffolding protein CheW (not shown), the histidine kinase CheA, the response regulator CheY, its phosphatase CheZ, and the receptor‐modification enzymes CheR and CheB. The phosphorylated form of CheY, CheY‐P, is the final output of the intracellular pathway that directly interacts with the flagellar motor to modulate swimming behavior. The steady‐state level of this signal is determined by the balance between production of CheY‐P, catalyzed by CheA, and its destruction, catalyzed by CheZ. We denote the activity of the receptor–kinase complex by a, which is feedback regulated by the methyltransferase CheR and the methylesterase/deamidase CheB through the receptor methylation level, m. The latter is defined as the average number of methylated glutamyl residues per receptor monomer—each receptor has a fixed number of these sites subject to reversible modification by CheR and CheB; for the aspartate receptor Tar, there are four, so the maximum value per monomer M≡max(m)=4. Whereas the receptor–kinase‐activity a is known to increase with m, that is ∂a/∂m>0, increased a in turn leads to up‐regulation of CheB (Borczuk et al, 1986), which removes methyl groups, and down‐regulation of CheR (Boldog et al, 2006), which adds methyl groups. Thus, the sense of the feedback is negative, tending to restore a toward its steady‐state value on changes in the input ligand concentration, [L]. CheZ, the phosphatase for CheY‐P, decreases the signal lifetime, thus accelerating the response of the pathway. In our FRET assay, the energy transfer pair is CheY and CheZ, fluorescently labeled by genetic fusion to yellow fluorescent protein (YFP) and cyan fluorescent protein (CFP), respectively. The FRET signal is thus proportional to [CheY‐P‐CheZ], the concentration of the intermediate species in the enzymatic hydrolysis of CheY‐P. At steady state, the production rate of CheY‐P, catalyzed by CheA, is exactly balanced by its destruction rate, which is proportional to [CheY‐P‐CheZ]. Hence, the FRET signal serves as a readout of CheA kinase activity, the output of the receptor module.

Figure 1.

A modular gradient‐sensing network. (A) Molecular view of the chemotaxis network. The linear path from input to output begins with the input ligand concentration, [L], being sensed by the membrane‐associated receptor–kinase complex, A, to regulate its autophosphorylation‐activity, a. A then transfers phosphate to the response regulator, CheY (Y), the phosphorylated form of which (Y‐P) interacts with the flagellar motor (M), to control swimming behavior. The feedback loop is closed by the methyltransferase CheR (R) and the methylesterase/deamidase CheB (B), by regulation of the receptor methylation level, m. CheZ (Z), the phosphatase for CheY‐P, decreases the signal lifetime, thus accelerating the response of the pathway. (B) Modular view of the network. Focusing on the functional modules, rather than the variables, of the network yields this block diagram, in which the variables are viewed as inputs or outputs (represented along wires) of two discrete signal processing modules (represented as boxes). The input–output relation of the receptor module is described by the function G, which takes [L] and m as inputs to produce an output a, which connects to the downstream linear pathway toward motor output. The adaptation module constitutes the feedback loop of the network, in which the output a is converted through F(a) to dm/dt and integrated over time.

In addition to the molecular constituents, we highlight in Figure 1A three important dynamical variables from this network, namely [L], a, and m, which, respectively, correspond to the input, output, and feedback signals. We capture the dynamics of these variables in two simple equations:

Embedded Image

Embedded Image

As the receptor‐modification reactions are much slower than all other reactions in the system, we model the rate of change of the feedback signal m explicity by a differential equation (equation (1)), and denote its dependence on the current signal output, a, by the function F. The receptor–kinase output, a, is known to relax much more rapidly (Sourjik and Berg, 2002a), and so its value at every instant is given by the algebraic equation (equation (2)), in which the function G denotes the dependence on the current ligand input, [L], and feedback signal, m. Thus, equations (1) and (2) describe the adaptation and receptor modules, respectively. This two‐module architecture is emphasized in the block diagram of Figure 1B, which clarifies the function of the two functions appearing in equations (1) and (2) within this signaling circuit: F is a transfer function contributing to the feedback gain of the network, whereas G combines the ligand input [L] and the feedback signal m to produce the network output a. We use for G an allosteric model of the MWC type (Monod et al, 1965), in which the effects on kinase activity of both the ligand input [L] and the feedback signal m contribute additively to the free‐energy balance governing the two‐state output of receptor–kinase complexes with N ligand‐binding subunits.

The properties of the functions F and G would depend on details of the many molecular components, but, importantly, each has only one or two dynamical inputs and both have a single output. The integral sign appearing before F in Figure 1B signifies integral feedback (Yi et al, 2000), which is a well‐known engineering design for robust adaptation, and is a direct consequence of equation (1): as the function F defines the feedback signal m's rate of change in time, m itself must be the time integral of F, that is Embedded Image. Thus, the function F is an important design feature of this signaling circuit, as it not only determines the dynamics of the feedback regulation (equation (1)), but also defines a crucial link between the two modules of the network: from the amplified output of the receptor module, to integral feedback in the adaptation module. Here, we experimentally probe this crucial link and quantitatively determine the feedback strength by measuring the response to time‐varying stimuli in vivo.

Exponential ramps shift the steady‐state kinase activity

We first confirmed that the kinase output during exponential ramps of the form [L](t)=[L]0ert as measured by FRET reaches a constant steady‐state level. To dynamically control the input [L](t), we fabricated a mixing apparatus (described in Materials and methods) based on the design of Block et al (1983). Cells were first adapted to a baseline concentration, [L]0, of the attractant α‐methyl‐dl‐aspartate (MeAsp). The FRET signal was stabilized at FRET0, although some fluctuations were unavoidable because of noise originating in the syringe‐pump‐driven injection of attractant to the mixing chamber. On applying both up ramps (r>0; Figure 2A, C, E, and G), and down ramps (r<0; Figure 2B, D, F, and H), the FRET signal was found to relax to a new, constant value, FRETc, as expected. As the magnitude of this change, ΔFRET (expressed in arbitrary units), is proportional to the change in kinase‐activity, Δa, and the latter, by definition lies on the unit interval, we can convert from FRET units to kinase‐activity (a) units by measuring the full range of ΔFRET. We obtain this by measuring the saturating values for ΔFRET in both directions, negative (ΔFRETsat) and positive (ΔFRETsat+), by addition and removal of a large concentration of attractant, respectively. Thus, the pre‐stimulus kinase activity a0=−ΔFRETsat/(ΔFRETsat+−ΔFRETsat), the change in kinase activity Δa=ΔFRET/(ΔFRETsat+−ΔFRETsat), and the constant kinase activity during exponential ramps ac=a0ac.

Figure 2.

Temporal profiles of pathway activity during exponential ramps. FRET responses (gray points) were observed during stimulation by exponential ramps in concentration of the form [L](t)=[L]0ert of the attractant MeAsp (blue curves). Responses were quantified by fitted functions to the FRET response (red curves). For both ramps up (r>0; A, C, E, G) or down (r<0; B, D, F, H), the FRET readout reached a steady level after an initial transient, which was well fit by a single‐exponential decay. The change in the level of FRET could be converted to units of kinase activity, Δa (see text), which was negative for up ramps, and positive for down ramps, as expected. The traces in the panels here were responses to ramp rates r=±0.001 (A, B), r=±0.005 (C, D), r=±0.01 (E, F), and r=±0.02 (G, H). The gray points in each panel are aligned averages of two to three separately measured responses to identical stimuli. Source data is available for this figure at

Ramp responses determine the feedback function and the gradient sensitivity

The constant response in kinase‐activity ac that is reached during exponential ramps (Figure 2A–H) can be viewed as the output of time‐derivative computations by the chemotaxis network. As the receptor module responds to the logarithmic change of the input signal (Tu et al, 2008; Kalinin et al, 2009), the relevant time derivative is that of the logarithm of input, that is dln[L]/dt, which corresponds in these experiments to the exponential ramp rate r. We therefore conducted exponential ramp‐response measurements of the type depicted in Figure 2A and B over a range of ramp rates r. The asymptotic kinase response, ac, obtained through such measurements, is plotted in Figure 3A as a function of r. The steady‐state activity in the absence of stimuli (i.e. at r=0) was found to be a0≈1/3. This plot reveals the sensitivity of E. coli to temporal gradients of MeAsp, and the overall shape is sigmoidal, with a steep slope (Δacr≈−30 s) near r=0. This implies that the system is tuned to respond sensitively to very shallow gradients, but it has a relatively narrow dynamic range: at greater absolute ramp rates, it becomes largely insensitive to changes in the gradient. If we define the slope Δacr as the gradient sensitivity, its value is large and nearly constant in the small interval near r=0, but decays rapidly outside of it. Importantly, we observed no response thresholds at small ramp rates (Figure 3A, inset), in contrast to Block et al (1983), in which it was found that the ramp‐response magnitude reached zero at low ramp rates (r≈0.005 for up ramps, r≈0.01 for down ramps).

Figure 3.

Gradient sensitivity and the feedback transfer function F(a). (A) By measuring exponential ramp responses in the manner of Figure 2A–H over a range of ramp rates r, we constructed a gradient‐sensitivity curve, relating the kinase‐activity a, to the steepness of the temporal gradient experienced by cells. The results for two FRET strains considered wild type for chemotaxis (VS104, an RP437 derivative, cyan circles; TSS178, an AW405 derivative, dark blue squares) were essentially identical, and collapsed on to a sigmoidal curve with a steep region near r=0 (slope of fitted line, Δar≈−30 s). The steady‐state activity in the absence of stimuli (i.e. at r=0) was found to be a0≈1/3. The inset in (A) is an expanded view about the point (r=0, a=a0), showing the absence of thresholds, at least down to r=±0.001 s−1. (B) Using our model, the data of (A) can be used to map the feedback transfer function F(a). The steady‐state relation rF(ac) implies that we can rescale the r axis by the constant factor α, obtained from our calibration of the receptor‐module transfer function G, and invert the axes about (r=0, a=a0) to obtain F(a). The shallow slope near this origin, F ′(a0)≈−0.01, implies weak negative feedback. The blue curve is a fit of a Michaelis–Menten reaction scheme Embedded Image; see text for interpretation of parameters. Source data is available for this figure at

The data of Figure 3A can be used further to infer the shape of the feedback transfer function F(a). The recorded values of ac are the outputs of the receptor module at which the effects of the time‐varying ligand input, [L](t), and that of the adaptation feedback, m(t), are canceling one another exactly. Therefore, given a model of how these signals are processed within the receptor module, that is if we know the form of the function G([L], m), we can infer the rate of change of the feedback signal, that is dm(t)/dt=F(a), from the temporal input profile [L](t) that is being cancelled. For an MWC‐type receptor module experiencing exponential ramp stimuli, one can show (see Materials and methods) that cancellation of the ligand‐ and methylation‐dependent free energy yields dfL/dt+dfm/dt=r−αF(ac)=0, where α≡−dfm/dm is the free‐energy change per added methyl group (in units of kT). Thus, we can rescale the abscissa of Figure 3A to render it a reading of F(ac)=r/α, and invert the axes about the point (r=0, ac=a0) to obtain a plot of F(a) (Figure 3B). The slope near the fixed point, F ′(a0)≈−0.01, is negative and shallow, implying weak negative feedback.

Frequency–response measurements reveal the time‐derivative computation and its bandwidth

Although the asymptotic gradient sensitivity (Figure 3A) can be viewed as the output of a time‐derivative computation, it does not address what happens if such gradients change over time. For a swimming bacterium in the real world, gradients will rarely be constant, so how rapidly the cell can respond to changing gradients will also be important. We therefore measured the frequency response of the network by measuring responses to oscillatory stimuli Embedded Image over a range of input frequencies, ν (Figure 4A). For small modulation amplitudes AL, the response in FRET, converted again to receptor–kinase‐activity units, was always found to be well fit by a sinusoid of the form a(t)=a0+∣A∣cos (2πνt−φD), where the frequency ν always matched that of the input modulation, whereas the response amplitude ∣A∣ and phase delay φD were found to depend on ν (see below). This is the expected linear response to sine‐wave modulation, consistent with the notion that the receptor module takes the logarithm of ligand concentration (Tu et al, 2008).

Figure 4.

Frequency response and the bandwidth for derivative computations. FRET responses (gray dots) during stimulation by exponential sinusoids of the form Embedded Image (A) with a frequency ν equal to that of the applied stimulus (green line), and the fit by a sinusoid (red line) of the form a(t)=a0+∣A∣cos(2πνt−φD), where ∣A∣ is the amplitude of the response and φD is the phase delay. When the amplitude (B) and phase (C) of responses are plotted against the driving frequency ν, the data are in excellent agreement with the analytical solutions of our model (equations (13) and (14)), plotted here without any free‐fitting parameters—the solution uses parameters obtained separately from the ramp‐response experiments of Figure 2 and 3 and MWC‐model parameters obtained separately in dose–response experiments using step stimuli. These analytical solutions define the characteristic frequency of the response νm≈0.006 Hz, below which the network is able to compute time derivatives. The dashed green curve in (B) is the derivative‐filtering function ∣H∣ obtained by factoring the solution of our linearized model (equation (15)). The characteristic frequency, νm determines the upper boundary of the frequency band over which time derivatives can be computed (shaded region). Presumably, the lower boundary would be determined by a noise floor, which in this figure is taken arbitrarily to be where the response amplitude (black curve in B) falls below 5% of maximum. Source data is available for this figure at

The measured frequency dependence of the response amplitude, ∣A∣, and phase delay, φD, are plotted as points in Figure 4B and C. Given the measured value for F ′(a0) obtained in our ramp experiments, our model can be used to determine the dependence of ∣A∣ and φD on frequency. These predictions, shown as solid curves in Figures 4B and C, are in excellent agreement with the data. This solution yields a characteristic frequency νm defining the upper limit of the frequency band over which the system is able to take time derivatives. Using the measured parameters a0≈1/3 and F ′(a0)≈−0.01 from our ramp‐response measurements above, and N=6, α=2 from our calibration of the receptor module, G (see Materials and methods), we obtain νm≈0.006 Hz. By appropriate factoring of our model solution (see Materials and methods), we also obtain a function ∣H∣, the chemotaxis system's filtering properties for the time derivative of the input signal (Figure 5B, green dashed curve). The shape of ∣H∣ clarifies the low‐pass property of this network for the derivative signal: the derivative signal is passed most efficiently below νm, in which the phase delay φD approaches π/2.

Figure 5.

Effects of temperature on sensitivity to gradients and frequency response. (A) Sensitivity to gradients was markedly decreased at 32°C (red triangles; strain TSS178), in which the steady‐state activity a0≈1/2. For comparison, the data at 22°C (same as Figure 3A) also are plotted (cyan circles; strains VS104 and TSS178). The slope of the linear fit to the 32°C points near a0 (red curve) was Δar≈−11 s. (B) The map of F(a) obtained by conversion of the data in (A) has a similar shape as that at 22°C, but the slope at the zero crossing, F ′(a0)≈−0.03, is approximately threefold steeper, implying stronger negative feedback. The red curve is a fit to the same Michaelis–Menten model as in Figure 3B (see text for parameter values and interpretation). (C) The frequency response is also shifted at 32°C (red triangles). The characteristic cutoff frequency νm≈0.018, obtained from the model fit (black curve), is approximately threefold higher than that at 22°C. For comparison, the 22°C data and the corresponding model prediction from Figure 4A also are reproduced here (cyan circles and blue curve). Source data is available for this figure at

Temperature affects both gradient sensitivity and frequency response

All of the above measurements were made at room temperature (22°C) for consistency with earlier in vivo FRET studies (Sourjik and Berg, 2002b, 2004), but the earlier results of Block et al (1983) were obtained at the higher temperature of 32°C. For comparison, therefore, we also collected ramp‐ and frequency–response data at 32°C. Interestingly, the gradient sensitivity (Figure 5A) was attenuated at this higher temperature, and the slope near a0 (≈1/2 at 32°C) was reduced to Δar≈−11 s. To compare this with the gradient sensitivity found by Block et al (1983), we must convert our response, Δa, measured in kinase‐activity units, into Δ〈CCW〉, the response in counter‐clockwise‐rotational bias of the flagellar motor. Using a Hill function Embedded Image for the motor response with nH=10 (Cluzel et al, 2000), we find Embedded Image, which is within ∼threefold of the Block et al result (Δ〈CCW〉/Δr≈20 s). However, no response thresholds were observed, even at this higher temperature (Figure 5A, inset). By the same procedure used to obtain Figure 3B from 3A, we can convert the gradient‐sensitivity data of Figure 5A to obtain the shape of the feedback transfer function F(a) at 32°C (Figure 5B). The frequency response of the system was also tested at 32°C using exponential sine‐wave stimuli, and was again found to match closely the model prediction with νm≈0.018, based on the slope F ′(a0)≈−0.03 measured at this temperature.

The shape of the gradient response (Figure 5A) obtained from our FRET experiments contrasts strongly from that obtained by Block et al (1983) through motor‐rotation measurements, in which a flat profile near r=0 had suggested the existence of response thresholds at small values of r. These differences remain unexplained, but we speculate that they are attributable to either or both of the following: (a) the difference in the output that has been measured; whereas we have measured the kinase output by FRET between CheY and CheZ, Block et al (1983) measured the motor response, or (b) the amount of data collected in the two studies; whereas our FRET technique has allowed us to efficiently collect a relatively large number of data points, each averaged over hundreds of cells (59 data points, each representing an average over hundreds of cells), the gradient response of Block et al (1983) was obtained using only a few cells (26 data points from 3 cells). The dependence of motor switching on [CheY‐P] was found by Cluzel et al (2000) to be monotonic, steep, and highly reproducible between separately tested individual motors; therefore, there does not seem to be any way in which the continuous variation in [CheY‐P] measured by FRET can give rise to a discontinuous motor response. Thus, (b) seems to be the most plausible explanation.


We have measured the dynamical response of E. coli's chemotactic gradient‐sensing circuit to time‐varying chemical stimuli. Our results reveal this signaling network's sensitivity to temporal gradients, as well as its frequency‐filtering properties. Using an important relationship identified by theoretical analysis (Tu et al, 2008), we derived from our exponential ramp‐response measurements the feedback transfer function F(a), which couples the two modules of this signaling pathway and reveals the dynamics of the adaptation module. Responses to oscillatory stimuli, the power of which has been exemplified in a number of recent studies (Bennett et al, 2008; Hersen et al, 2008; Mettetal et al, 2008), were used to characterize the frequency response of the pathway, which determines the frequency band over which the chemotaxis system can compute time derivatives (Andrews et al, 2006). We found the characteristic frequency of the adaptation system to depend on temperature, but it is in general compatible with weak feedback and large noise amplitudes at steady state. The use of time‐varying stimuli allowed us to probe these dynamical characteristics of the adaptation module without direct measurement of the underlying biochemical reactions of receptor methylation/demethylation. This general approach of combining theory with experiment to probe the dynamics of pathway components not directly accessible to experiment should be applicable to many systems, and aid in the effort to uncover the design principles of biological circuits (Alon, 2007). We discuss here the implications of our findings in the context of earlier studies, and highlight possible directions for future work.

Effects of cell‐to‐cell variability: measurements in single cells and populations

A crucial difference between the measurements reported here and those of Block et al (1983) is that our FRET signals are collected simultaneously from hundreds of cells, whereas Block et al recorded the motor response of individual tethered bacteria. Although our data have much higher signal‐to‐noise ratios than the stochastic binary time series collected in the motor‐response experiments, it is pertinent to ask whether the population‐averaged nature of our FRET measurements might account for some of the differences in our observations. In particular, a striking difference was found in the exponential ramp responses to very shallow gradients (Figure 3A, inset); whereas Block et al could not detect responses to ramp rates in the interval r∈(−0.01, +0.005) s−1, no such thresholds were found in the exponential ramp responses we measured by FRET. Could this difference arise from ensemble averaging the output of individual cells, which individually possess such response thresholds?

We argue here that this is not possible if the typical single cell does indeed possess such response thresholds at the level of the receptor–kinase response. In the context of our preceding analyses, the qualitative question of threshold existence can be recast as a quantitative one of gradient sensitivity: does there exist a region of negligible gradient sensitivity at very low ramp rates (Δacr≈0 near r=0)? To make the distinction between individuals and populations explicit, we denote by aci(r) the steady‐state activity of each individual cell in a population during an exponential ramp with rate r, and by 〈ac〉(r) the population‐averaged activity we measure by FRET. As the latter is just a linear combination of single‐cell activities, Embedded Image, where Ncells is the number of cells in the population, the gradient sensitivity at the population level is just the average of the single‐cell gradient sensitivities Embedded Image. Thus, if the majority of cells possess response thresholds, this should have been evident in our population measurements.

Transient time of ramp responses: inference of cooperativity from dynamics

The gradient sensitivity, which we defined using the constant kinase‐activity ac reached during exponential ramps, can be viewed as the output of time‐derivative computations by the chemotaxis system. Remarkably, the solution for this output in our model, ac=F−1(r/α), is independent of the parameter N, which represents the degree of receptor cooperativity in the MWC model. Thus, although receptor interactions have an amplifying function in the open‐loop response to step and impulsive stimuli (Sourjik and Berg, 2004; Keymer et al, 2006; Mello and Tu, 2007; Tu et al, 2008)—the property traditionally referred to as the ‘gain’ of the chemotaxis system (cf. e.g. Segall et al, 1986)—the amplitude of responses to sustained temporal gradients is dictated only by the adaptation system (through the function F(a)). This does not mean, however, that receptor cooperativity is unimportant within the ‘closed‐loop’ control structure of the chemotactic gradient‐sensing circuit. The function of receptor coupling becomes clear when one considers the time required for the chemotaxis network to take this time derivative, which is an equally important measure of performance for this gradient‐sensing system. Our model provides a basis for investigating the molecular parameters that determine the speed of this computation, and yields the following simple, approximate expression for the time, t1, required for the kinase activity to reach ac during exponential ramps:

Embedded Image

where Δf(ac) is the free energy difference between the steady states with and without the ramp: Δf(ac)=ln(1/ac−1)−ln(1/a0−1). The quantity Δf(ac)/N can be intuitively understood as a characteristic amount of chemical work that the adaptation system must perform before it can catch up with the rate r at which the ligand concentration is being ramped.

In the limit of shallow ramps (small r), where Δac is small, we can linearize equations (1) and (2) to solve for a limiting, minimum response time, t1 → τm=(−αNa0(1−a0) F ′(a0))−1=(2πνm)−1, where νm is the characteristic frequency of the linear response to oscillatory stimuli (see equation (14) of Materials and methods). This solution is independent of the applied (small) ramp rate r. At greater ramp rates, the linear approximation is no longer valid, but we can still estimate t1 by the simple relation presented in equation (3). The justification for this approximation is illustrated in Figure 6A, which is a plot of the time course of the two free‐energy contributions to kinase activity, stemming from ligand binding (fL) and methylation (fm). When an exponential ramp with rate r is started at t=0, fL begins to rise immediately as ∼rt. The methylation‐dependent contribution, −fm, ‘lags’ behind fL, but this lag asymptotically stabilizes as the methylation rate F(a) reaches r/α, from which point on, the free‐energy lag Δf(a), and hence the kinase activity become constant (i.e. Δf (a) → Δf (ac), as a(t) → ac).

Figure 6.

Transient responses to ramps: the time required for derivative computations. (A) Schematic illustration of the changes in activity and free energies during exponential ramps. As the receptor–kinase activity approaches the constant value ac during exponential ramps, the rate of change of ligand‐ and methylation‐dependent free energies (fL and fm, respectively) balance one another, leading to the asymptotic relation rF(ac) between the exponential ramp rate r and the net rate of change in intracellular methylation F(a) (see also Materials and methods). The time required for fm(t) to change by an amount equal to the total free‐energy change Δf(ac) can be estimated as t1≈Δf(ac)/Nr, and provides a robust estimate for the time required for the activity a(t) to reach ac (see text). (B) A plot of rt1 versus Δf(ac) for exponential ramp responses at 22°C. The slope of the fitted curve yields an estimate for the extent of receptor cooperativity, N≈4.9, which is in good agreement with the value of N=6 obtained independently from dose–response curves from experiments using step stimuli (see text).

In the linear regime in which the transient time is independent of r, a(t) will decay exponentially in time as Embedded Image toward ac, and equation (3) is exact. The constant solution t1m in this case can be obtained by letting rF′(a0)(aca0) in equation (3) and taking the limit aca0. In the more general case, t1 will not be independent of r, but equation (3) still provides a useful approximation. In particular, it retains the correct scaling with the extent of receptor cooperativity, N. To compute the points in Figure 6B, t1 was estimated by fitting an exponential decay function of the form Embedded Image to the FRET data recorded during exponential ramp stimuli (such as those of Figure 2) at different values of r, and Δf(ac) was computed from the a0 and ac values from the same measurement.

As Δf (ac) is independent of the degree of receptor cooperativity, N, we can use equation (3) to extract from our experimental data an estimate for the parameter N. Although it recently has been reported (Endres et al, 2008) that the degree of receptor cooperativity seems, under certain conditions, to depend on the modification level of receptors (i.e. that N could depend on m, which is changing during our ramp measurements), the fact that we observe a constant activity during exponential ramps suggests that such changes in the organization of the receptor complex occur on time scales considerably slower than that of our time‐varying stimuli. We, therefore, apply equation (3) to infer a single value of N from the set of t1 values extracted from our ramp‐response data. Although accurate determination of t1 is difficult, because its estimation is more sensitive to experimental noise, the approximate linear scaling of the product rt1N−1 Δf (ac) follows immediately from equation (3), and is more robust to noise than t1 itself (as the factor r conveniently diminishes the error more strongly in which Δac is small, and errors in estimating t1 tend to be large). In Figure 6B, we plot rt1 against Δf (ac), in which the slope of the fitted line gives an estimate N∼4.9, which is in satisfactory agreement with the more reliable value of N=6, estimated through fits to dose–response curves to step stimuli (Mello and Tu, 2007; see also Materials and methods and Figure 7).

Figure 7.

Calibration of the receptor‐module transfer function, G([L],m). (A) FRET responses in kinase activity (points) to step stimuli of the attractant MeAsp were measured in a series of mutants in which the modification state of Tar receptors were fixed by deletions of the cheR and cheB genes. To reveal the response of the Tar receptor, these measurements were conducted in a genetic background in which Tsr and Tap receptors are not expressed. Fits to the allosteric MWC model (equations (6), (7) and (8)), with parameters N=6, KI/KA=0.0062 are shown as solid lines. The gray level of the points and solid curves indicate the modification state of the mutant strain, and is tabulated in the legend (lighter shades of gray for higher modification levels). The blue points and curve, denoted wt in this figure, is the strain VS178, which has the same receptor complement as the other mutants (Tsr−Tap−), but retains the wild‐type genes for CheR and CheB. Kinase activity is shown normalized to the pre‐stimulus value for the strain VS104. (B) Values of fm, obtained for each modification state represented in (A), plotted against the number of modified sites. The data for modification states containing only glutamate (E) and glutamine (Q) residues (black points) fell on a straight line (dotted) when plotted as a function of the number of E → Q transitions (nE → Q; black points). The data for modification states containing only Q's and methylated glutamates (Em) fell on a different straight line (dashed), when plotted as a function of the number of Q → Em transitions (nQ → Em; magenta points). The values for the two extreme methylation levels, EEEE (corresponding to m=0) and EmEmEmEm (corresponding to m=4), were obtained by extrapolation of the two straight lines, as they could not be obtained directly from fits to the FRET data (cells with the Tar population fixed in this state did not respond to any concentration of MeAsp, as is seen in (A)). The solid line connecting these two extreme states reveals the dependence of fm on the E → Em transitions (nE → Em; gray points), and the CheR+CheB+ strain VS178 (blue point, denoted wt) falls on this line, as expected. When fitted by equation (9), this line yields the parameters m0≈0.5 and α≈2 kT used throughout this study. Source data is available for this figure at

Biochemical interpretation of F(a): enzyme saturation and CheB activation

The steady‐state identity rF(ac) derived from our model allowed us to map the feedback transfer function F(a) from the gradient‐response data measured at each temperature. We found that F(a) is a nonlinear, monotonically decreasing function with a shallow slope near the steady‐state fixed point, F(a0)=0, and a much steeper slope near the high‐activity extreme, F(1). As F(a) represents the net rate of change of methylation catalyzed by CheR and CheB, a biochemical interpretation of this curve can be obtained by an enzymatic reaction model,

Embedded Image

where KR and KB are the Michaelis constants for each reaction, and VR and VB(a) represent the velocities of the methylation and demethylation reactions, respectively, when those enzymes are saturated with substrate (all concentrations are normalized, in units of the CheA kinase concentration; for example KR,B are dimensionless, VR,B have units s−1). This expression assumes that CheR can only bind the inactive fraction and CheB can only bind the active fraction of receptor–kinase complexes. Our measurements place strict limits on these end point values, F(0) and F(1). In our ramp‐response data, F(0) can be inferred from the ramp rate rc at which the response amplitude to positive ramp‐stimuli saturates (i.e. in which ac reaches zero) by the relation rcF(0). Reading off values for rc from Figure 5B, and using again α=2 from our receptor‐module calibration (described in Materials and methods), we expect F(0)≈0.01 at 22°C and F(0)≈0.025 at 32°C. The data are noisier at the opposite extreme of activity, but an approximate lower bound for the absolute value of F(1) can readily be inferred by examining the maximal absolute response observed, which yields ∣F(1)∣>∼4F(0) at both temperatures. These constraints, combined with the high curvature near F(1) and the location of a0, precluded a fitting with a constant value for VB. Evidently, what is required is a sharp increase in VB(a) as a approaches unity. The solid curves shown in Figures 3B and 5B use a simple piece‐wise linear form, Embedded Image, where θ(x) is the unit step function (θ(x)=1 for x>0, θ(x)=0 otherwise), aB represents the value of a above which a VB(a) increases, and rB determines the maximum value of VB(a). For the data of Figure 5B, we found that aB=0.74 yields good fits for data at both temperatures. Values for other parameters are as follows: at 22°C, {VR,VB(0)}={0.010, 0.013} s−1, {KR, KB}={0.32,0.30}, and rB=4.0; at 32°C, {VR,VB(0)}={0.030, 0.030} s−1 and {KR, KB}={0.43,0.30}, and rB=2.7.

The comparison between the two sets of experiments offers some insight regarding the temperature dependence of the parameters. The dominant difference is in the saturating velocities for methylation and demethylation, VR and VB(a), which can be interpreted as the product of the active enzyme concentration and the catalytic rate constant kcat. As the growth conditions before the experiments were identical, and measurements were conducted in a medium that does not support protein synthesis, the expression levels of enzymes were the same in the experiments at 22 and 32°C. Thus, we conclude that kcat's for methylation and demethylation are the parameters most sensitive to temperature in the adaptation system.

The mechanism responsible for the strongly nonlinear behavior of VB(a) remains to be determined. A likely candidate is the phosphorylation of CheB, but we note that straightforward mechanistic assumptions, such as a linear, hyperbolic, or quadratic dependence of [CheB‐P] on a would not suffice to produce the observed sharp increase of VB(a) as a approaches unity. Future experiments with phosphorylation‐deficient mutants could shed light on this issue.

Noise and feedback dynamics near steady state

Single‐cell measurements of motor output (Korobkova et al, 2004) have established that the steady‐state activity of the E. coli chemotaxis pathway exhibits large fluctuations over a slow timescale. Through theoretical modeling and simulations, Emonet and Cluzel (2008) have argued that this could be attributable to the adaptation enzymes CheR and CheB operating near saturation, a phenomenon related to zero‐order ultrasensitivity of reversible covalent‐modification systems (Goldbeter and Koshland, 1981). Our data provide an experimental test for this hypothesis, as the biochemical parameters inferred from the above analysis of F(a) can be used to assess the degree of saturation of the in vivo adaptation kinetics, and, in turn, the steady‐state noise.

The values we obtained above of the Michaelis constants, KR and KB, and the pre‐stimulus kinase‐activity a0, imply that >∼70% of both CheR and CheB are bound to their receptor substrates at steady state. Here, we apply the stochastic analysis of Emonet and Cluzel (2008) to our biochemical model for the adaptation system (equation (4) of the main text), to estimate the steady‐state noise expected for kinase activity in vivo. The feedback transfer function F(a) determines the dynamics of the adaptation system through equation (1). In the absence of chemical stimuli, the pathway activity, in turn, is determined completely by the adaptation system; by taking the time derivative of equation (2) and substituting into this equation (1), we find Embedded Image. Applying this result to the stochastic time‐evolution equation obtained by Emonet and Cluzel (2008) by the ‘linear noise approximation’, the displacement from the steady‐state activity, Δa(t)=a(t)−a0, evolves in time as

Embedded Image

where τa=(−αNa0(1−a0)F ′(a0))−1 is the relaxation time of the receptor–kinase‐activity a, η(t) represents temporally uncorrelated white noise, and Embedded Image is the strength of noise originating from the discrete methylation events occurring in Ntot receptor complexes. Computing the value of the relaxation time constant τa using the parameters a0, F ′(a0), N and α from our measurements, we obtain τa≈29 s at 22°C and τa≈11 s at 32°C. The solution for the amplitude of the steady‐state fluctuations, σa, in the Emonet–Cluzel model (equation (5)) is Embedded Image. With parameters for the CheR/CheB reactions inferred from our measured profile of F(a), and assuming Ntot=104/N (i.e. partitioning the ∼104 receptor dimers in each cell into MWC clusters of size N), we obtain σa/a0≈0.087 at 22°C and σa/a0≈0.077 at 32°C.

Thus, our results are consistent with slow steady‐state fluctuations in kinase activity, with a large amplitude of around ∼10% of the mean and relaxation time constants in the range ∼10–30 s. In an analysis of the noise‐spectral data obtained by Korobkova et al (2004), Tu and Grinstein (2005) found that intracellular [CheY‐P] fluctuations must be both slow (the correlation time, τa≫1 s) and of large amplitude (>20% of the mean) to explain the power‐law dwell‐time distribution in motor switching, as well as the 1/f noise observed in the switching‐time power spectrum. According to the above analysis, the requirement for long correlation times (i.e. large τa) is clearly satisfied. Although the calculated noise amplitudes seem to fall somewhat short, we note that we have estimated these values from population measurements. Given the natural variation in protein levels between individual cells, we would expect the reported 1/f behavior to be measurable in some fraction of the population.

Concluding remarks

In this study, we have viewed the chemotaxis‐signaling pathway as a signal processor, and characterized its performance by two phenomenological observables: the ‘gradient sensitivity,’ a quantity we defined using the constant activity reached by the pathway during exponential ramp stimuli (Figures 3A and 5A), and the frequency response, described here by the amplitude and phase of the response to exponential sine waves (Figures 4 and 5C). Calibration of the allosteric MWC model for the receptor module using step stimuli (Figure 7) revealed a linear dependence of receptor‐free energy on the level of covalent modification. This finding admitted a particularly simple solution to our pathway model (equations (1) and (2)), in which the gradient sensitivity depends on the way in which the adaptation and receptor modules are coupled in vivo (through the function F(a) of the adaptation module and the parameter α of the receptor module), but not explicitly on receptor cooperativity (the parameter N of the receptor module). The frequency response, however, depends also on the degree of receptor cooperativity (both the cutoff frequency for time‐derivative computations and the amplitude are proportional to N).

Although the ramp‐response data allowed us to deduce the underlying chemical kinetics of the adaptation enzymes (Figures 3B and 5B), we note that the frequency response, with its dependence on N, is relevant in considering the chemotactic performance of bacteria executing spatial searches in the real world. Swimming bacteria execute random walks by a run/tumble mechanism (Berg and Brown, 1972) to convert spatial gradients of chemoeffectors into temporal ones. In such searches, time‐derivative computations must be made while swimming direction is being randomized. Thus, the sensitivity to sustained temporal ramps cannot be directly interpreted as the sensitivity of swimming cells to spatial gradients. The frequency response is more relevant in considering this problem, as we can view the input fluctuations because of randomized swimming as additional ‘noise’ with its characteristic frequency spectrum (determined by the statistics of tumbling and rotational Brownian motion). The question of how such noise arising from random motility is filtered during spatial searches has been addressed implicitly in a number of recent studies, in which numerical simulations were used to predict the migratory behavior of bacterial populations with different adaptation time scales (Andrews et al, 2006; Bray et al, 2007; Emonet and Cluzel, 2008; Jiang et al, 2010). Our measurements and analyses of the frequency response reported here provide a firm grounding for further analytical treatments of this problem (Sartori and Tu, submitted). Viewed as a signal filter, a crucial characteristic of the frequency response is the cutoff frequency, νm, below which time‐derivative sensing can occur. As νmNF′(a0), the receptor cooperativity N can have a buffering function in keeping νm reasonably high when adaptation feedback is weak (i.e. when the value of F′(a0) is small).

Another important source of noise in chemotactic signaling, namely that due to fluctuations in the intracellular concentrations of signaling species, could be addressed directly by our current findings. The relatively slow relaxation times we obtained for the pathway activity near steady state provide support for the view that the large‐amplitude fluctuations observed by Korobkova et al (2004) could originate in the attenuated feedback near steady state (i.e. small F′(a0)) resulting from enzyme saturation. This result also highlights the importance of using dynamically modulated stimuli to probe the adaptation kinetics near steady state—the often used strategy of linearly interpolating between the limiting adaptation kinetics measured during large step responses (F(0) and F(1) in our model) would have provided a very different estimate for F′(a0), and hence the relaxation time. Moreover, our data and analysis make clear that the strength of feedback affects not only the steady‐state fluctuations, but also the gradient sensitivity and frequency response of the pathway, as evidenced by measurements at two temperatures (Figure 5). It is then reasonable to expect that the pathway's relaxation time scale, which could be tuned by the expression level or kinetic parameters of CheR and CheB, is under selective pressure for its effects on gradient‐sensing parameters such as the cutoff frequency for time‐derivative computations, in addition to its effects on the steady‐state noise.

Materials and methods

Strains and plasmids

All bacterial strains used in this study were derivatives of E. coli K12 strain AW405 (Armstrong et al, 1967). VS104 is a cheYcheZ mutant (Sourjik and Berg, 2002b) of a widely studied AW405 derivative, RP437 (Parkinson and Houts, 1982). TSS178 is an analogous cheYcheZ mutant of the parent strain AW405, constructed for this study. For the step‐response measurements to calibrate the receptor‐module transfer function, G, we used Tsr− Tap− CheR− CheB− derivatives of VS104 in which the tar receptor gene is mutated at its covalent‐modification sites (gifts of V Sourjik). The identity of these strains were, in order of increasing amidation levels, VS144 (EEEE), VS141 (QEEE), VS148 (QEQE), VS150 (QEQQ), SB1 (QQQQ). The strain VS178 is Tsr− Tap−, but retains the wild‐type cheR and cheB genes. Plasmid pVS113 encodes the cheR gene under the control of an arabinose‐inducible promoter, and plasmid pVS88 encodes the fluorescent fusion proteins CheY‐YFP and CheZ‐CFP under control of an isopropyl β‐d‐thiogalactopyranoside‐inducible promoter (Sourjik and Berg, 2004).

In vivo FRET measurements and data analysis

FRET microscopy of bacterial populations were carried out essentially as described earlier (Sourjik et al, 2007). The FRET donor–acceptor pair was CheZ‐CFP and CheY‐YFP, expressed from a plasmid in strain backgrounds that lack the native copies of the cheY and cheZ genes. Cells attached to poly‐l‐lysine‐treated microscope coverslips were seated at the top face of a flow cell (Berg and Block, 1984). The microscope was a Nikon TE300 equipped with a PlanFluor × 40 0.50 n.a. objective. The sample was illuminated by a 75 W super quiet xenon lamp (Hamamatsu, Bridgewater, NJ) through an excitation bandpass filter (Semrock FF01‐438/24‐25) and a dichroic mirror (Semrock FF458‐Di01‐25 × 36), and epifluorescent emission was further split into donor (cyan, C) and acceptor (yellow, Y) channels by a second dichroic mirror (Chroma 515DCXR) and collected through emission bandpass filters (Semrock FF01‐483/32‐25 and FF01‐542/27) by photon‐counting photomultipliers (Hamamatsu H7421‐40).

Signal intensities of the acceptor and donor channels were recorded by a computer running LabView (National Instruments), and the ratio between the two channels (R=Y/C) provided an indicator of FRET activity that is robust to fluctuations in excitation intensity. The change in FRET efficiency, ΔFRET, can be computed at every time point from the recorded data using the ratio change ΔRS≡(R−Rpre) as ΔFRET=(RpreRS−R0)/(RpreRS+∣ΔYC∣)−(RpreR0)/(Rpre+∣ΔYC∣), where R0 is the acceptor‐to‐donor ratio in the absence of FRET, Rpre is pre‐stimulus ratio of acceptor‐ to donor‐channel intensities (measured either in the absence of attractant, or after adapting to a constant value of [L]), and ∣ΔYC∣ is the (constant) absolute ratio between the changes in the acceptor‐ and donor‐channel signals per FRET pair (Sourjik et al, 2007). Under the conditions of the measurements performed here, however, backgrounds in the Y and C channels were negligible, Rpre and ∣ΔYC∣ were practically constant, and we found consistently that Rpre+∣ΔYC∣≫ΔRS. Thus, ΔFRET≈ΔRS/(Rpre+∣ΔYC∣) was essentially proportional to ΔRS (worst case nonlinearity ∼3%) and comparable between measurements. We, therefore, expressed ΔFRET, for simplicity, in arbitrary units of ΔRS throughout this study.

MWC model for the receptor module G([L], m)

In this study, we use a calibrated allosteric model of the receptor module, G. The shape of G has been mapped in detail in recent in vivo FRET experiments of kinase activity through dose–response measurements using step stimuli. This function characterizes the response of the receptor–kinase complex in Figure 1A, taking as input [L] and m to yield the output a. Much recent work has shown this function to be well approximated by a two‐state model,

Embedded Image

wherein the total free‐energy difference between the two output states, ft, is simply proportional to the sum of ligand‐dependent and modification‐dependent parts:

Embedded Image

When the parameter N is an integer greater than unity, it can be interpreted as the number of ligand‐binding subunits in an oligomeric complex with tightly coupled output (i.e. their activity transitions are concerted so that the collective output has only two states), and is equivalent to the MWC model of allostery (Monod et al, 1965). An important assumption of the model is that the ligand dissociation constant for binding sites depends on the activity state of the oligomer. This leads to a simple analytical form for the ligand‐dependent free energy per ligand‐binding subunit, fL:

Embedded Image

where KA and KI are the ligand dissociation constants for the active and inactive states of the oligomer, respectively. The model does not specify the functional form of the modification‐dependent free energy, fm(m), but this can be determined experimentally by fitting the MWC model (equations (6), (7) and (8)) to dose–response curves of kinase activity. In vivo FRET experiments, which measured the response of the aspartate receptor Tar when stimulated by the attractant MeAsp yielded a linear form for fm:

Embedded Image

where α is the free‐energy change per added methyl group (Figure 7); we note that these measurements were made at room temperature, that is ∼22°C). With N=6 and KI/KA=0.0062, α was found to have a value ∼2 kT, and the crossover methylation level, m0 (defined as the methylation level at which fm(m)=0) was found to be ∼0.5. This calibration of the function G provides the basis for our strategy to probe the adaptation module of the pathway through the receptor module, as it establishes the quantitative relationship between effects of methylation and ligand concentration on kinase output, which we monitor by FRET.

Input–output measurements by FRET

We used three classes of time‐varying stimuli to probe the transfer functions of the chemotaxis system. Step stimuli were used to characterize the function G([L], m) for the Tar receptor in an array of Tsr− Tap− CheR− CheB− mutants with fixed receptor‐modification states. As these cells lack the receptor‐modifying enzymes CheR and CheB, the input–output relation of the receptor module can be measured in an open‐loop configuration, without adaptation feedback. The lack of the Tsr and Tap receptors renders the Tar receptor the dominant majority in the remaining receptor population. In this genetic background, step responses to rapid changes in [MeAsp] were measured in ‘modification‐standard strains’ bearing various point mutations at the four modification sites of the Tar receptor, corresponding to the EEEE, QEEE, QEQE, QEQQ, and QQQQ modification states (in which the four‐letter code denotes the state of the amino‐acid residue at the four modification sites; ‘E’ signifies a glutamate, and ‘Q’ a glutamine). As the glutamine (Q) residues differ in their chemical structure from a glutamate (E) only by the addition of an amide group, E → Q substitution mutations have long been used as an experimental proxy for studying the effects of the more biologically relevant methylated glutamate residues (Em), which cannot be encoded genetically. Here, we also studied the effect of Em modifications by overexpressing CheR from a plasmid in the same ‘modification‐standard strains,’ thus yielding, in addition to the above, the modification states QEmQQ, QEmQEm, QEmEmEm, and EmEmEmEm. Dose–response data were collected for strains representing these nine states, and for Tsr− Tap− CheR+ CheB+ cells (Figure 7A). Fitting of these data to the MWC model (see Materials and methods) yielded estimates, in units of free energy, of how strongly each added Q‐ or Em‐modification affects the activity of the Tar receptor (Figure 7B). We found that Em‐modifications have an ∼2.5‐fold stronger effect than Q‐modifications.

We then used exponential ramp stimuli to probe the shape of the feedback transfer function F(a) (equation (1)). The main idea of the measurement hinges on the observation by Block et al (1983) that the pathway output reaches a constant value during exponential ramps of the form [L](t)=[L]0ert (Figure 2A and B). Within our MWC model of the receptor module, the dependence of kinase activity on the total free‐energy ft is monotonic (equation (6)). Therefore, as the activity a(t) approaches a constant, so does the free energy, that is

Embedded Image

In our calibrated MWC model (the function G([L],m) in equation (2)), the ligand‐dependent free energy, fL([L]) (equation (8)), possesses a broad (∼2.5 orders of magnitude) domain in its input ligand concentration [L] over which fL is essentially linear in ln[L]. This implies that an exponential ramp that operates within this range (KI<<[L]<<KA) at rate r will cause a constant change in fL per unit time, at precisely the same ramp rate r (i.e. fL(t)=rt). On the other hand, our receptor‐model calibration also obtained a linear relationship between the modification‐dependent portion of free energy, fm(m), and the methylation level, m (equation (9)), so combined with our model for the rate of change of methylation (equation (1)), the methylation‐dependent free energy follows fm(t)=−αF(a)t. Combining these with the constancy condition for the total free energy (equation (10)), we obtain a simple relationship between the applied ramp rate r and the in vivo feedback transfer function F(a) during exponential ramps,

Embedded Image

Thus, by controlling the input, r while measuring the constant output, ac the functional form of F(a) can be mapped by plotting r/α against the measured values of ac.

The third class of input–output measurements we performed involved measuring the response in kinase activity to exponential sine‐wave stimuli of the form Embedded Image

(Figure 4A). The response was always found to be sinusoidal, a(t)=a0+∣A∣cos (2πνt−φD) where the frequency ν (in Hz) always matched that of the input, as expected for a linear system. We can solve for both the response amplitude ∣A∣ and phase delay φD in our model by linearizing equations (1) and (2). This gives

Embedded Image

Embedded Image

for the response amplitude, and

Embedded Image

for the phase delay, where in equations (12), (13) and (14), νm=−αNa0 (1−a0) F ′(a0)/2π defines a characteristic frequency determined by the adaptation kinetics and receptor parameters. Considering the two limits of equation (12), ϕD → π for ν≫νm and ϕD → π/2 for ν≪νm, makes it clear that νm delineates the range of possible input frequencies ν into two regimes: below νm, the system takes the time derivative of the input (i.e. the output is 90° out of phase), whereas above νm, the system simply follows the input in phase (although the output seems 180° out of phase, because of the negative response of kinase activity to attractant). Finally, by taking the Fourier transform of the signal's time derivative Embedded Image and factoring this out from the general solution (equation (12)), we obtain the system's time‐derivative‐filtering function

Embedded Image

The absolute value of this function, ∣H∣, plotted in Figure 4B, clarifies the low‐pass‐filtering properties of the chemotaxis pathway for the input signal's time derivative.

Dynamic modulation of chemoeffector stimuli

Control of the temporal profile of attractant stimulus was achieved by a mixing apparatus (Figure 8A), based on the design of Block et al (1983). The mixing chamber was a cylinder fabricated of Delrin, with an internal volume Vmix≈100 μl, with two input channels (β and γ) and two output channels (δ and ε). The rate of attractant delivery, β, at concentration [L]β, by a syringe pump (Harvard Apparatus PHD 22/2000) was controlled through the same LabView program collecting the data to facilitate input/output comparisons. A peristaltic pump (Rainin Rabbit) pushed buffer into the chamber through the other input channel at a constant rate, γ. Negative pressure was applied to both output channels of the mixer, the first by a peristaltic‐pump pulling fluid at a constant rate, δ, through the flow cell, and the second by a gravity‐induced hydrostatic force (applied by lowering the height at which the outlet of the connected tubing rested), which pulled fluid to drain at a rate ε. When β is constant, the output concentration [L](t) approaches [L]=[L]ββ/(β+γ) as Embedded Image, where [L]0 is the initial concentration and τmixVmix/(β+γ). In all of our experiments, we set γ≫β so that [L]≈[L]ββ/γ and τmixVmix/γ. Thus, when β is dynamic, [L](t) is proportional to β as long as the latter is varied smoothly over time scales much longer than τmix (Figure 8B).

Figure 8.

Chemical waveform generator: design and performance. (A) The chemical waveform generator, based on the design of Block et al (1983). A mixing chamber of volume Vmix was fed two inputs: attractant at a high concentration at computer‐controlled rate β and buffer at rate γ. The output concentration [L], seen by cells in the flow cell, perfused at rate δ, was proportional to the rate β, provided that γ≫β. Under these conditions, the time for mixing, τmix, was fixed at Vmix/γ. Thus, by programming the pump controlling β, arbitrary waveforms could be generated as a function of time. In our experiments, Vmix≈120 μl and γ>1200 μl/min, so τmix was <6 s. The rate δ does not affect τmix, but affects the exchange time τex of fluid in the flow cell and the polyethylene tubing connecting it to the mixing chamber, so this was also kept relatively high (∼1000 μl/min). (B) We confirmed, by measuring the optical density of bromothymol blue (black points), that temporal gradients of arbitrary function could be programmed (blue curve) provided that the frequency of the programmed signal did not exceed the mixing frequency νmix=(2πτmix)−1. This condition was violated for the blue segment at the left end of the plot.

Conflict of Interest

The authors declare that they have no conflict of interest.

Supplementary Information

Source data for figure 5A [msb201037-sup-0001-SourceData-S1.csv]

Source data for figure 5B [msb201037-sup-0002-SourceData-S2.csv]

Source data for figure 5C [msb201037-sup-0003-SourceData-S3.csv]

Source data for figure 2A [msb201037-sup-0004-SourceData-S4.csv]

Source data for figure 3A [msb201037-sup-0005-SourceData-S5.csv]

Source data for figure 3B [msb201037-sup-0006-SourceData-S6.csv]

Source data for figure 4A [msb201037-sup-0007-SourceData-S7.csv]

Source data for figure 4B [msb201037-sup-0008-SourceData-S8.csv]

Source data for figure 4C [msb201037-sup-0009-SourceData-S9.csv]

Source data for figure 7A [msb201037-sup-0010-SourceData-S10.csv]


We thank Junhua Yuan, Ady Vaknin, Victor Sourjik, Thierry Emonet, and Philippe Cluzel for helpful discussions. This work was supported by National Institutes of Health Grants AI016478 (to HCB) and GM081747 (to YT) and a National Institutes of Health Postdoctoral Fellowship AI063747 (to TSS).


This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.