Exponential Signaling Gain at the Receptor Level Enhances Signal-to-Noise Ratio in Bacterial Chemotaxis

Cellular signaling systems show astonishing precision in their response to external stimuli despite strong fluctuations in the molecular components that determine pathway activity. To control the effects of noise on signaling most efficiently, living cells employ compensatory mechanisms that reach from simple negative feedback loops to robustly designed signaling architectures. Here, we report on a novel control mechanism that allows living cells to keep precision in their signaling characteristics – stationary pathway output, response amplitude, and relaxation time – in the presence of strong intracellular perturbations. The concept relies on the surprising fact that for systems showing perfect adaptation an exponential signal amplification at the receptor level suffices to eliminate slowly varying multiplicative noise. To show this mechanism at work in living systems, we quantified the response dynamics of the E. coli chemotaxis network after genetically perturbing the information flux between upstream and downstream signaling components. We give strong evidence that this signaling system results in dynamic invariance of the activated response regulator against multiplicative intracellular noise. We further demonstrate that for environmental conditions, for which precision in chemosensing is crucial, the invariant response behavior results in highest chemotactic efficiency. Our results resolve several puzzling features of the chemotaxis pathway that are widely conserved across prokaryotes but so far could not be attributed any functional role.


Quantitative chemotaxis pathway model
The signal signaling gain of the E. coli chemotaxis pathway is modelled in this work by g = 18 receptor homodimers that form independent allosteric units and change activity in unison [1]. Here, the probability to find an active receptor complex takes the form P (t) = [1 + exp[−g∆F]] −1 , with ∆F = − + m − f (L(t)) ( = 1.0, = 0.44obtained from a linear fit to the data shown in Ref. [1]). ∆F is a function of the average methylation level per receptor dimer, m , and the free energy contribution of attractant binding, f (L(t)) = i g i /g[ln(1 + L(t)/K of f i ) − ln(1 + L(t)/K on i )], with L(t) a specific ligand concentration at time t and K on i and K of f i the corresponding dissociation constants of ligand binding to the g i receptors of type i within an allosteric unit in the on and off state, respectively [1]. The global methylation dynamics follows the rate equation˙ with the activity dependent functions F i , i ∈ {1, 2, 3}. Measurements of the methylation rate in the regime of low receptor activity are consistent with the approximation F 1 [P (t)] = const [2]. For demethylation, the situation is more complicated, as phosphorylation of CheB shows more than one order of magnitude higher demethylation activity but is non-essential for adaptation. Thus, F 3 is expected to show a monotone increasing dependency on P (t), even in the low activity regime. However, the contribution F 3 can be neglected to good approximation for the native expression level of CheY [3] as the ∼ 70-fold more active phosphorylated form of CheB dominates the demethylation kinetics [4]. This is confirmed by a non-phosphorylatable mutant of CheB, expressed at the estimated native level, showing an order of magnitude higher adapted kinase activity than what is expected for wild type conditions ( Fig. 4F and Fig. S3). Note that the high induction of the CheY-YFP/CheZ-CFP FRET pair (c.f. Fig. S5) leads to a significantly higher kinase activity, for which the F 3 contribution cannot be neglected any more. As the phosphoflux from CheA to CheB, v B , determines the amount of phosphorylated CheB, P (t) ∼ v B ∼ B p , a linear dependency of the methylation rate on the mean receptor activity is the consequence,˙ m ∼ const − P (t), if the F 2 term is approximately constant. This linear relation has been experimentally confirmed in the low activity regime [2].

Pathway response dynamics
Bacterial chemotaxis results from a positive (negative) correlation between the average duration of a straight swimming run and the temporal change of attractant (repellent) concentration. The alternating periods of swimming (counter clockwise flagella rotation) and tumbling (clockwise flagella rotation) of a motile E. coli cell can be described to good approximation by a straight swimming cell that is subject to rotational diffusion. This is possible as the measured directional changes, ϕ, between swim runs are moderate, with standard deviation ϕ 2 ≈ 58 • [5]. We can therefore write the total rotational diffusion coefficient of the cell, D R = D B + D ϕ , as a superposition of Brownian motion, D B ≈ 0.068 rad 2 s −1 (30 • directional change per second in agar, with viscosity 2.7 cp at 32 • C [6]), and stochastic tumbling events D ϕ = ϕ 2 /(4τ (t)) that are separated on average by the time span τ (t) ( Here, k B is the Boltzmann constant, T , the temperature, η the solvent viscosity, and a ≈ 1µm the effective hydrodynamics radius of an E. coli cell. For moderate chemoeffector gradients the average duration of a swim run, τ (t), shows only small changes for cells that have moved up or down the gradient [6]. This allows us to linearize around the average duration of a mean swim run, τ , in presence of a stationary chemoeffector gradient to give τ (t) = τ + δτ , with δτ = 0. The averaging, indicated by the brackets . , runs over a sufficiently large amount of tumbling events and defines the quasi-stationary state as indicated by the overline. With the definitions D ϕ = ϕ 2 /(4τ ) and D R = D B + D ϕ we get for Eq. (2) the expression The temporal change in the probability to find a receptor in the active state, P (t), depends on the encountered changes in chemoligand. In case of a stationary chemoeffector landscape this probability depends exclusively on the history of the swimming trajectory weighted by the loss of memory due to adaptation with rate γ (see Methods section in the main text). From Eq. (9) of the main text we find after linearization around the stationary states y = z(t)H(u) and f (L) the expression which results for the gain function P (t) : with , δy = y(t) − y, δf (L) = f (L) − f (L), and γ := N r. Here, ] is the per receptor free energy contribution of ligand binding that together with the free energy contribution from chemoreceptor methylation, F, determines the activity state of the receptor cluster [1]. In Eq. (6) we have defined the stationary receptor activity by P := P a + δP (t) , where P a denotes the adapted receptor activity. We emphasize that the definition of a stationary receptor activity, P , makes only sense in case of a sufficiently slowly varying extracellular chemoeffector concentration. The fraction of time the cell spends in clockwise rotation (CW) can be described to good approximation by the expression [7] with n = 10 the Hill coefficient of the motor response curve and Y p the concentration of phosphorylated CheY. With the definition of the clockwise bias CW = τ R /(τ + τ R )with τ R the average tumbling time -we arrive with a linear dependency between receptor activity, P (t), and phosphorylated CheY, Y p (t) ≈ y(t) = z(t)P (t), at the relation Here, we have introduced with τ d the response time needed for an extracellular change in chemoligand concentration to affect the swimming behavior of the cells. Note that z(t − τ d ) ≈ z(t) by definition. For small fluctuation δP (t) around the stationary activity state of the receptor clusters, P , we find for the relative change in the mean swim duration the relation Here, we included the reasonable assumption that the signaling time delay τ d is to good approximation the consequence of an exponentially decay of the phosphorylation and binding reactions. To calculate the drift velocity in a linear gradient we employ the approximation where the term L (t) arises as δτ is measured relative to the stationary swim run duration τ . Here, L(t * ) is the chemoattractant concentration at time t * ∈ [t 0 , t] for a cell located in vicinity of the position x(t * ) in direction of the gradient. The expansion coefficient, h, is given by .
For stationary attractant gradients, L(t * ) =L(x(t * )), we can use the expression with the attractant gradient at position x given by ∇L(x) and V D the drift velocity of the cell.

Drift in linear gradients
In the following we denote the current angle of the swim direction to the attractant gradient by θ(t). This allows us to rewrite the actual position of the cell, x(t), in direction of an attractant gradient as with V the effective swim velocity along the trajectory. Thus, the expression for changes in the average swimming time can be rewritten with the help of Eq. (11) as Further note that the quantity we are after -the net drift in gradient direction V D -is given by which requires the calculation of probability distribution, P (θ, t), to find the cell in orientation θ at time t. This probability distribution satisfies the rotational diffusion equation with L the classical analog to the quantum mechanical angular momentum operator on the unit sphere L = −ir × ∇ = i θ (sin θ) −1 ∂ φ −φ∂ θ . A formal solution to Eq. (19) is given by with T the time ordering operator. In the following we expand the solution of P (θ, t) to first order in the smallness parameter D ϕ δτ (t)/τ by employing Eq. (3) and using Dyson series This expansion allows us to give an accurate, explicit expression for the chemotaxis drift velocity, V D , for moderate gradients Due to interruption by tumbling events, a reduced effective swim velocity, V , is the consequence, which is given to leading order by V = τ τ +τ R V 0 , with V 0 the swimming velocity of the cell. We thus arrive at the final expression for the drift velocity, with τ R = 0.14 s the average tumbling time. Our analysis shows that the time a cells spends in the tumbling mode should be minimized to arrive at highest chemotactic drift. However, minimizing tumbling times means that eventually rotational Brownian motion will dominate directional changes, overriding the chemotaxis strategy to prolong swimming when moving towards favorable environments and to shorten swim run duration when moving away from favorable environments. Therefore, swim run duration should dominate the tumbling duration but should remain far below the average time needed for a cell to reverse direction due to Brownian motion, τ B . As the adaptation time should be long enough to allow for reliable comparison of chemoeffector concentrations at different times but short enough to cut off any information that is not correlated to the actual swimming we arrive at the relation τ R < τ , γ −1 < τ B for highly chemotactic cells.

Comparison with computer simulations
To show the regimes of validity of the analytical approximation for the net drift velocity in direction of an attractant gradient we employed computer simulations to generate swimming trajectories of individual cells. To this end we integrated Eq. (4) numerically using equidistant time steps of 0.005 seconds. After every time step the probability to enter or leave a tumbling and swimming mode is calculated. The mean duration of swim run in the adapted state is fixed to 0.8072s and the adaptation rate is set to 0.4022s −1 . The later corresponds to an effective methyltransferase activity of k R R = 0.05674s −1 . In the computer simulations, swim and tumbling duration are taken from exponential distributions, to generate similar statistics as experimentally observed [8,9]. The tumbling angle between swim runs is chosen to be Gaussian distributed with a standard deviation of 58 • and rotational diffusion changes swimming direction by 30 • per second. To allow an easy comparison we set the signaling time to zero, τ d = 0.

Molecular description of the chemotaxis pathway
The central argument why mass action equations are sufficient to describe cellular signal pathways (see below) is a consequence of intracellular noise reduction. As the effect of noise is always linked to both molecule copy number and the reaction time scale, we expect that the fast response times of the pathway, τ ∼ 10 −1 s, requires sufficiently high abundance of the proteins per cell (strain RP437, [3]) of CheY (∼ 8200 copies), CheZ (∼ 3200 copies), and CheA (long isoform) (∼ 4500 copies). In contrast the more than one order of magnitude slower adaptation kinetics can buffer the small concentration for CheB (∼ 240 copies) and CheR (∼ 140 copies) [3]. We emphasize that noise due to low copy number of the protein CheR is measurable but influences only to small extend the chemotactic efficiency of a population [10].

Methylation dynamics
The molecular mechanisms of how the average cellular receptor activity, P (t), drives the methylation rates are currently not fully understood. However the leading order functional dependencies of the adaptation dynamics can be determined from response measurements of this work. In the following we make use of the fact that the phosphorylation dynamics is much faster than the methylation dynamics and thus can be treated as equilibrated on the adaptation time scale. The mass action equations for the methylation dynamics of receptors of type α ∈ {1, 2, 3, 4}, (e.g. Tar, Tap, Tsr, Trg) in the methylation state m, are given by with H 1 and G 1 monotone increasing functions of the total concentrations of the methyltransferase CheR, R T , and methylesterase CheB, B T , respectively, and H 2 and G 2 the unknown dependencies on the probability to find a receptors in the active state, P (t), at time t. Averaging Eq. (31) over all receptor types and methylation states results in with m = α,m m T m , assuming that the highest methylation state is not significantly populated within the physiological regime. Note that there must not exist any dependency on the methylation state in Eq. (33) as a consequence of perfect adaptation [11]. The functional form of H i and G i , i ∈ {1, 2} can be inferred from experiments, shown in Fig. 2 in the main text. Here, an increase in CheR abundance leads to both an increase in the activity [12] and the adaptation rate, whereas an increase in CheB reduces kinase activity but leaves the adaptation rate unchanged. As the probability to find a receptor in the active state, P (t), equilibrates rapidly and is determined by the chemoattractant concentration, L, and the average methylation level, m , the receptor activity has on the mean field level a functional time dependence given by P (t) = P ( m (t), L(t)). By taking the time derivative, ∂ t P (t) = ∂ m P ∂ t m (t) + ∂ L P ∂ t L(t), we can rewrite Eq. (33) as From Fig. 2A in the main text follows that multiplication of Eq. (34) by the total concentration of CheB, B T , shows identical adaptation dynamics over one oder of magnitude in CheB abundance. Therefore a transformation of Eq. (34) to a scaled activity, P B (t) := P (t)B T , must not show any explicit dependency on B T to leading order, Independence of the solution P B (t) on the CheB concentration requires the functional forms H 2 = 1 and with a yet unknown exponent α, which results in where we have introduced the proportionality constant b. For the adapted state (∂ t P B (t) = 0 and ∂ t L(t) = 0) follows that . As the adapted state is unchanged for concerted overexpression of CheR and CheB (data not shown) we obtain the equation where we introduced the proportionality constant r. The yet unspecified derivatives of P B (t) with respect to m and L must be functions that depend uniquely on P B (t), , because of the dynamic scaling property Fig. 2B in the main text. We therefore arrive at the equation As the adaptation rate increases linear with R T [12] (see also Fig. 2E in the main text) and P a B ∝ R T (Fig. S2) to good approximation in the low activity regime, we observe to leading order invariance of the functionP B (τ ) := P B (τ /R T )/R T , with τ := t R T , for different concentrations of CheR and stepwise changes in attractant concentration. Multiplication of Eq. (38) by (R T ) −2 gives Note that a step change of attractant with magnitude ∆L at time t = 0 results in the expression ∂ t L(t) = ∆L δ(t) = ∆L δ(τ /R T ) = R T ∆L δ(τ ), with δ(t) the delta distribution. Elimination of all dependencies on CheR from Eq. (39) requires the remaining unknown functionals to take the forms f m (x) ∝ x 2−α and f L (x) = l x. A least squares fit (see Fig. 2 in the main text) reveals that the choice α = 1 reproduces accurately the response of E. coli to chemoeffectors where we have to introduce a negative proportionality constant, l < 0, if L(t) reflects the ambient attractant concentration. Resubstitution gives the final result This logistic differential equation reproduces accurately the response dynamic as shown in Fig. 2 in the main text by the solid lines. The solution of Eq. (41) reads This result shows that the time needed to recover the pre-stimulus value depends on the methyltransferase concentration, R T , but not on the methylesterase concentration, B T . So far we have given very strong evidence that the integral feedback that leads to adaptation is mediated by the activity dependent term k B b B T P (t). For this term there exists experimental evidence that CheB phosphorylation is the dominant but not the only contribution. This is because E. coli cells show adaptation to chemoeffectors for a CheB mutant that cannot be phosphorylated. We conjecture that there is second contribution resulting an integral feedback where CheB bind with weak affinity to active receptors. Thus there is evidence that the integral feedback derived from the scaling arguments above is biochemically realized by two mechanisms, resulting in

Phosphorylation dynamics
The phosphorylation dynamics of CheY has be extensively studied in E. coli. Phosphorylation takes place at the kinase CheA, which consists of a P 2 binding domain and a P 1 domain that transfers phospho-groups to CheY and CheB. Both CheY and CheB are most effectively phosphorylated by binding to the P 2 domain but can also receive phospho-groups by a direct phosphotransfer from the P 1 domain. Transfer of phospho-groups from the ATP binding pocket to the P 1 domain is believed to be the rate limiting step in the phosphorylation cascade and dephosphorylation by CheZ is the by far strongest contribution to CheYp hydrolysis. This gives rise to the equation with Z T the total concentration of CheZ phosphatases, k Z the dephosphorylation rate, K Z the corresponding Michaelis-Menten constant, A the concentration of non occupied kinases associated to fully functional receptor complexes and P (t) denotes the average receptor activity. The whole phosphorylation network is presented in the following. We assume that autophosphorylation of the kinase CheA depends on the average receptor activity P (t), and thus gives rise to the reaction equation Phosphorylation of the response regulators CheY and CheB predominantly takes place by binding to the P 2 domain of CheA and subsequent phosphotransfer from the phospho-rylated P 1 domain to CheY and CheB, respectively. This is summarized in the reaction equations Finally, dephosphorylation of CheYp is mediated by the phosphatase CheZ whereas phosphorylated CheB autodephosphorylates Applying the law of mass action, we obtain the set of rate equations where complexes are denoted by brackets. These equations are completed by the conservation laws Here, A c denotes the total concentration of kinases associated to fully functional receptor complexes. Due to the small amount of CheB, i.e. B T A c and B T Y T , we can neglect all interactions of CheB with CheA in the differential equation governing the dynamics of the phosphorylation state of CheA. Hence, instead of Eqs. (49) and (55) we consider and we get the set of differential algebraic equations where we used the set of conservation laws (56)-(58) and (60) and introduced the Michaelis-Menten constants Phosphorylation dynamics is much faster than methylation of the receptor. Thus we can assume that phosphorylation is equilibrated on the timescale the average receptor activity P (t) changes due to methylation. To obtain this quasi steady state, we have to solve the set of algebraic equations The phosphotransfer rate from CheA to CheB to leading order in P (t) is given by In the adapted state, i.e.Ṗ (t) = 0, it follows from Eq. (6) in the main text that holds and thus, using Eq. (77), we get with P a the adapted receptor activity. Solving for P a yields From Eqs. (70) and (71) follows. Making the same assumption as above, i.e.
If we assume that k A k Y and since by definition P (t) ≤ 1 holds, we get The complex [Y pZ] can be written as a function of Y p and Z T , Using Eq. (87) and (88) we get for (85) We obtain for the adapted state by using the approximation for P a derived in (84) independent of k A and A c .

Absence of CheZ dependence on kinase activity
To show that CheZ does not alter kinase activity we employed the fact that phosphorylated CheB strongly localizes at the receptor complexes [13]. By substituting native CheB by a catalytic inactive mutant that is fluorescently tagged, CheBS164C-YFP, we were able to measure CheBS164C-YFP localization at the receptors by employing the aspartate receptor fusion construct Tar-CFP. Localization is proportional to the FRET signal Tar-CFP/CheBS164C-YFP and can be compared, after stimulation with methylaspartate, with the localization of unphosphorylated CheBS164C-YFP. We employed a CheY deleted strains to avoid any interference between the two phospho-receivers CheB and CheY. The FRET amplitude arising from the difference between the FRET signal before and after the stimulus shows no significant difference if CheZ is present at native level or absent: 0.0154 ± 0.0027 for the ∆CheY strain and 0.0162 ± 0.0032 for the ∆CheYZ strain. Kinase activities for three independent experiments for each strain are listed in the table below.
Exp. No. kinase activity (AU), ∆CheY kinase activity (AU), ∆CheYZ 1 0.0163 0.0151 2 0.0175 0.0136 3 0.0124 0.0198 Table S1 5 Cell-to-cell variability of receptor clusters  Figure S4: Histogram of the cell-to-cell variability of the total amount of receptors that are organized in clusters. The variability is significantly higher than expected from gene expression noise alone, as additional noise contributions arise from stochasticity in receptor cluster assembly and uneven distribution of cluster among daughter cells.

Response behavior under exponential attractant ramps
To test the quality of the mathematical model shown in the Box of the main text we reproduced the response behavior of the flagellar rotation bias upon temporal changes in attractant (α-methyl-dl-aspartate). The flagella rotation bias is a highly significant measure of the kinase activity and has the advantage not to rely on fluorescent fusion constructs. Fig. 1 shows the change in counter clockwise (CCW) bias in exponentially increasing (upper red area) or exponentially decreasing (lower red area) attractant concentrations over time. For the calculation we employed the linearized version of the free energy contribution of ligand bindingṠ(t) ≈ h(L)L(t). The red areas arise from inserting the whole range of possible expansion parametersL ∈ [43.2µM, 448µM ]. Here, 43.2µM is the initial attractant concentration and 448µM is the final attractant concentrations. The attractant concentration at time t is given by L(t) ∝ exp[a t], with a the ramp rate [8,9]. We employed an adaption rate γ = 0.4 s −1 for 36 • C that has been inferred from FRET measurements of the kinase activity. Note that the adaptation rate is highly temperature dependent and decreases approximately eight fold at 20 • C,    Table S2. The errors are given by the 68% confidence intervals.
In order to calculate absolute concentrations of CheY-YFP from FACS, two samples, one from a population with 10 µM and the other with 25 µM IPTG induction were measured with a fluorometer and FACS. The fluorometer was calibrated by measuring a series of solutions of known YFP concentrations. Furthermore, the number of cells per volume was determined by two different methods, i.e. plating and Neubauer chamber. We calculated for each IPTG-induction a weighted mean for the cell density and estimated the error based on the observed range of measured values. In this way we obtained for the molecule numbers per cell

Calibration of Imaging Scales
Another method to measure expression levels is given by analyzing microscopy images. This methods yields single cell data and is important for quantifying expression levels in swarm assays, where bacteria are placed on an agar medium in a petri-dish and are allowed to perform chemotaxis for several hours. Depending on the chemotaxis efficiencies of the bacteria strains, a characteristic pattern emerges: bacteria with low chemotaxis efficiency reside in the middle of the petri dish, whereas the better performing phenotypes are found in an outer ring. In order to calibrate the fluorescence scales for YFP and CFP, we have to assume, that the bacteria in the inner region are equivalent to those grown in liquid medium. The latter ones are used for FACS, Fluorometer and FRET-experiments. At first we check whether the linear scaling for CFP-and YFP also holds in this case. Calculating the mean of the YFP-and CFP-fluorescence for different IPTG-levels shows, that linear scaling applies as shown in Fig. 6. We further assume, that the expression levels are a similar function of the IPTG-concentration as given by Eq. (91). However, fluorescent intensities measured by imaging differ by a scaling factor α I , which is determined by fitting this expression to the data obtained from swarm assays, using only α I as a free parameter. The remaining parameter values are the same as given in Eq. (91). The resulting curve for YFP is shown in Fig. 7. The obtained scaling factors are given by α YFP Fluorescent intensity I is a measure for molecule numbers per cell, whether measured by FACS, fluorometer or microscope. However, biologically more important are molecular concentrations. This was accounted for by determining also the length L of bacteria, which allowed us to calculate a concentration measure, I/L. Until now we assumed, that it is valid to use I/L = I / L . This however does not have to be true, since L and I could be correlated. First of all, L does not seem to be a function of the used IPTG concentrations, as depicted in Fig. 8. The mean length is given by L = 50 ± 13, which was calculated by using all data. Plotting I/L and I / L as a function of the IPTG concentrations give similar results, see Fig. 9. Hence, it seems to be valid to assume I/L = I / L . This allows us to determine the "length concentrations" I/L as a function of IPTG, yielding This can now be used to calibrate the imaging scales in units of µM, but the estimated error seems to be too small, since it does not take into account the uncertainty of determining the functional dependence of FACS on IPTG.