Ultrasensitivity and fluctuations in the Barkai-Leibler model of chemotaxis receptors in Escherichia coli

A stochastic version of the Barkai-Leibler model of chemotaxis receptors in Escherichia coli is studied here with the goal of elucidating the effects of intrinsic network noise in their conformational dynamics. The model was originally proposed to explain the robust and near-perfect adaptation of E. coli observed across a wide range of spatially uniform attractant/repellent (ligand) concentrations. In the model, a receptor is either active or inactive and can stochastically switch between the two states. The enzyme CheR methylates inactive receptors while CheB demethylates active receptors and the probability for a receptor to be active depends on its level of methylation and ligand occupation. In a simple version of the model with two methylation sites per receptor (M = 2), we show rigorously, under a quasi-steady state approximation, that the mean active fraction of receptors is an ultrasensitive function of [CheR]/[CheB] in the limit of saturating receptor concentration. Hence the model shows zero-order ultrasensitivity (ZOU), similar to the classical two-state model of covalent modification studied by Goldbeter and Koshland (GK). We also find that in the limits of extremely small and extremely large ligand concentrations, the system reduces to two different two-state GK modules. A quantitative measure of the spontaneous fluctuations in activity is provided by the variance σa2 in the active fraction, which is estimated mathematically under linear noise approximation (LNA). It is found that σa2 peaks near the ZOU transition. The variance is a non-monotonic, but weak function of ligand concentration and a decreasing function of receptor concentration. Gillespie simulations are also performed in models with M = 2, 3 and 4. For M = 2, simulations show excellent agreement with analytical results obtained under LNA. Numerical results for M = 3 and M = 4 are qualitatively similar to our mathematical results in M = 2; while all the models show ZOU in mean activity, the variance is found to be smaller for larger M. The magnitude of receptor noise deduced from available experimental data is consistent with our predictions. A simple analysis of the downstream signaling pathway shows that this noise is large enough to affect the motility of the organism, and may have a beneficial effect on it. The response of mean receptor activity to small time-dependent changes in the external ligand concentration is computed within linear response theory, and found to have a bilobe form, in agreement with earlier experimental observations.


Introduction
The pioneering work [1] of Goldbeter and Koshland (GK) brought into light an interesting molecular switch-like transition, now referred to as Zero Order Ultrasensitivity (ZOU) [2], observed in reversible covalent modification (e.g. phosphorylation or methylation) of a protein (substrate), catalyzed by two antagonistic enzymes. This switch-like behavior emerges in the limit where the substrate concentration is exceedingly large compared to the enzyme concentrations as well as their individual Michaelis constants. In this "zero-order" regime, the modified fraction of substrate, in a saturating environment, exhibits a sharp transition at a critical value of the ratio of the total concentrations of the antagonistic enzymes. Various aspects of GK switch has been studied over the years [3][4][5]. The fluctuations associated with the ultrasensitive module has also been studied [3,4]. ZOU was also identified as allosteric cooperativity by Quian (2003) [6] and Ge and Quian (2008) [7].
The bacterium E. coli has thousands of receptors on its cell surface for a precise sensation of the extracellular environment. The two main types of chemotaxis receptors are Tar and Tsr. The receptor protein and the protein kinase CheA are linked by the linker protein CheW. The receptor protein, CheA and CheW function as a single signaling complex, which exists in active or inactive state, and can undergo stochastic switching between the states. The switching is regulated by methylation and ligandation (binding of chemoattractant/repellent) of the receptors. Methylation is carried out by the protein CheR while CheB demethylates receptors. In its active state, the receptor-CheW-CheA complex phosphorylates CheB and the response regulator CheY; while the former demethylates the receptor, the latter induces tumbles by binding to the flagellar motor. In this way, the internal biochemistry (methylation/demethylation) and external stimulus (attractant/repellent) regulates the swimming pattern of the organism (reviewed in [8]).
A remarkable property of chemotaxis in E. coli is perfect adaptation to some chemoattractants. After exposure to a stimulus corresponding to an abrupt rise or fall in ligand concentration in a homogeneous environment, the frequency of tumbles returns to its pre-stimulus level after a time interval *4 seconds [9]. A model of the methylation-demethylation cascade proposed by Barkai and Leibler (BL) [10] successfully explains this property using a few simple assumptions. The model was modified and extended in later years by other authors [11][12][13][14]. The important assumptions in the model are (a) CheR methylates inactive receptors while CheB (or its phosphorylated form, CheBp) demethylates active receptors, (b) the state of activity of a receptor depends on its methylation level, (c) binding of an attractant (ligand) molecule adversely affects the activity of a receptor, (d) ligand binding and dissociation are very fast compared to methylation and demethylation. Mello and Tu [13] showed that perfect adaptation is achieved only if the lowest methylation level is assumed to be always inactive, while the highest level is assumed to be always active. As a consequence of these assumptions, the mean receptor activity in the model in steady state turns out to be independent of the attractant concentration. Further, mathematical modeling has shown that, with suitable extensions, the model also successfully predicts the response of the network to a short-lived spike in attractant concentration [15].
In E. coli, the most abundant receptor, Tar, responding to the attractant methyl aspartate, has four methylation sites per monomer. For modeling purposes, it is convenient to represent the methylation state of a single receptor by an index m, where 0 m M with M being the maximum value (M = 4 for Tar). It is generally assumed that activation/inactivation as well as ligand binding/unbinding of a receptor happens much faster compared to methylation/ demethylation. For M = 1, the BL model is entirely identical with the Goldbeter-Koshland (GK) two-state system; the ligand concentration L plays no role whatsoever in the dynamics.
For M ! 2, the intermediate state(s) which stochastically switches(switch) between active and inactive states modifies the dynamics; however, the model retains the same structure as the GK module (see Fig 1, where we have provided a schematic depiction of the BL model, for M = 2). Therefore, it is pertinent to enquire whether ZOU is present also in the BL model. In the GK system (M = 2 in our model), fluctuations in activity are known to be large near the ZOU transition. In E. coli chemotaxis, intrinsic fluctuations are expected to play an important role. The flagellar motor responds ultrasensitively to [CheYp] (p denoting the phosphorylated form); the sensitive part of the motor response curve is less than 1μM in width [16,17], hence it is likely that spontaneous and random fluctuations in activity are crucial in generating the runand-tumble motion of the organism [12]. For this reason, we believe that a detailed study of ZOU in the BL model, including a systematic treatment of fluctuations will make a useful addition to existing literature on ZOU as well as chemotaxis in E. coli. This is the motivation behind this paper.

Model details
Let us denote by V the volume of a cell having a total of N = A 0 V receptor proteins (substrate), and two enzymes. In E. coli, these correspond to the methyltransferase CheR and methylesterase CheB, the total concentrations of each are denoted by R 0 and B 0 respectively. In the presently accepted model of signal transduction in E. coli, CheR methylates the receptor, while CheB demethylates it. According to the two-state model [18] the receptor exists in a complex form involving two more proteins, CheW and CheA, and the complex may exist in one of two conformational states, which we refer to as active and inactive. In the active state, CheA undergoes autophosphorylation and also phosphorylates the cytoplasmic messenger protein CheY. CheA phosphorylates CheB also, a process crucial in adaptation, but ignored here for the sake of maintaining close structural similarity with the GK system. The probability to be in the active state is an increasing function of the methylation level of the receptor. Chemoattractant molecules, if present, also bind to the receptor, and this process increases the probability to find the receptor in inactive state (with chemorepellents, it is the reverse) [8].
BL proposed a quantitative model [10] of the signal transduction cascade, based on the two-state model of Asakura and Honda [18]. In this model, CheR binds only to inactive receptors, while CheB binds only to active receptors [19]. In general, we denote by M the total number of methylation sites. Let A m symbolically represent a receptor in methylation state m, and A m denotes its enzyme-bound form; the same symbols will also be used to denote the corresponding numbers (out of a total of N). The inactive and active versions are identified by superscripts i and a respectively. The binding rate of CheR to a receptor is k + and k − is the dissociation rate; for CheB, the corresponding quantities are k 0 þ and k 0 À respectively. Attractant binding is relevant only for the intermediate state; we assume that the binding rate is k a , dissociation rate is k d and attractant concentration is L (assumed uniform). ν r and ν b represent the product formation rates from CheR and CheB-bound intermediate states, respectively. The concentrations of the free forms of CheR and CheB are denoted by R f and B f respectively. We also define a set of dissociation constants: If ξ m (t) denote the fraction of receptors in methylation state m at time t, the total active fraction of receptors is given by x a ðtÞ ¼ P M m¼0 a m ðLÞx m ðtÞ, where a m (L) is the probability for a receptor in state m to be active, L being the (uniform) attractant concentration. We assume that the lowest state (m = 0) is always inactive, while m = M is always active; the intermediate state(s) can be either active or inactive, depending on whether it is attractant-bound. A simple mathematical form satisfying these conditions is a m (L) = a m (0)K m /(L + K m ) where K m is the dissociation constant for attractant binding to a receptor in state m [14,15,20]. To satisfy the earlier assumption, we choose K 0 = 0 and K M = 1 (these are ideal values). We further choose a 0 m ¼ m=M [19], a simple mathematical form which satisfies our earlier assumptions.

Fokker-Planck equation for M = 2
We will first carry out a detailed mathematical analysis of the properties of the BL model with M = 2. Fig 1 shows a schematic depiction of this case, showing all the biochemical reactions. In the limit where the modification rates ν r and ν b are very small in comparison with the other rates, it may be visualized as a combination of three weakly-coupled modules (the boxes in Fig 1); the fractional populations of receptors in each module are ξ m with m = 0, 1, 2; hence P 2 m¼0 x m ¼ 1. In general, for a complete stochastic description of the dynamics of this system, we also need to keep track of time evolution of the intra-modular populations viz., x m andx m . For m = 1, there is a further sub-division: as explained in the previous section. For each m, x m ¼ x m þx m . It was shown by us in an earlier paper [2] that these intra-modular fractions can be 'integrated out' using one of the quasi-steady state approximation schemes (QSSA) [21], whereby the replacement x i ! " x i ðξÞ naturally appears in the relevant inter-modular conversion rates. Here, we have introduced the compact vector symbol ξ (ξ 0 , ξ 2 ) (because of the normalization constraint mentioned above, only two of the fractions ξ m are independent, we shall take these to be ξ 0 and ξ 2 purely for reasons of convenience). Therefore, within this approximation, the inter-modular dynamics may be expressed in terms of a probability density P(ξ;t), which satisfies a non-linear Fokker-Planck equation (FPE) in the form of a local equation of continuity [2]: where J {J 0 , J 2 } with J m ¼ v m P À @ @x m ðD m PÞ (for m = 0, 2) denoting the components of the state space probability current "vector". Here, the drift and diffusion coefficients are given by Incorporating the rates ω mn in the above expressions for v m and D m and employing the inter-module dynamics under the standard quasi-steady state approximation (sQSSA) as discussed in Appendix 1, we arrive at the following expressions for the drift and diffusion coefficients: and The free enzyme concentrations R f and B f are also functions of ξ 0 , ξ 2 , and the expressions are given by (see Appendix 1)

Averages as fixed points
From Eq (1), it is seen after the required averaging that the time evolution of the mean receptor fractions is given by the vector equation where vðξÞ ðv 0 ðξÞ; v 2 ðξÞÞ is the drift "vector". In general, this vector vanishes at one or more points in the (ξ 0 , ξ 2 ) space; following van Kampen [22], we shall refer to this point (assuming only one exists) as the fixed point, denoted ξ Ã ¼ ðx Ã 0 ; x Ã 2 Þ. It can be shown (see later) that in the large N-limit, the fixed point becomes identical to the mean value " ξ. Under conditions R f ( K r and B f ( K b , which we shall assume to hold for reasons of simplicity, it follows from Eq (3) that the fixed point is given by where R Ã f R f ðξ Ã Þ and B Ã f B f ðξ Ã Þ are given by Eq (5), with the replacement ξ ! ξ Ã in the right hand side of the equations. Eqs (5) and (7) implicitly give the fixed point ξ Ã .

Evaluation of the covariance matrix from the linear FPE
Next, we expand v m and D m in Taylor Series about the fixed point values. Let us define the deviation ξ 0 = ξ − ξ Ã ; we also define its probability distribution F(ξ 0 ) P(ξ Ã + ξ 0 ), which satisfies the equation where J 0 ðξ 0 Þ Jðξ Ã þ ξ 0 Þ, and may be expanded as follows: n @v m @x n ξ Ã þ ::: Keeping the leading term in each, and noting that vðξ Ã Þ ¼ 0, we arrive at (for m = 0, 2) where b mn ¼ @v m @x n j ξ Ã , D Ã m D m ðξ Ã Þ. Substitution of Eq (10) in Eq (8) leads to the multivariate linear Fokker-Planck equation (LFPE) [22] Let us now define the covariances s mn ¼ hx 0 m x 0 n i for m, n = 0, 2 (for m = n these become the variances), which are evaluated in a convenient way by defining the moment generating function (Fourier transform) where μ = (μ 0 , μ 2 ). The moment generating function admits the Taylor expansion m m m n s mn ðtÞ þ ::::::::: Since F(ξ 0 ) satisfies the LFPE given by Eq (11), the moment generating function satisfies the equation Under stationary conditions (t ! 1), the left hand side of Eq (14) becomes zero; next we substitute Eq (13) in Eq (14), and find that " ξ 0 ¼ 0 in the long-time limit. Therefore, within LNA, the fixed point values are equal to the statistical means of the respective quantities. The steady state covariance matrix satisfies the Lyapunov equation [22] bσ þ σb T where D mn ¼ D Ã m d mn are the elements of the (diagonal) diffusion matrix D. The following expressions follow from Eq (15): Explicit expressions for the coefficients β mn are to be found in Appendix 2.

Mean and fluctuations in activity
The active fraction of receptors in the present model is given by where ℓ = L/K L . Within LNA, the average fractions " x Ã a is given by Eq (17), with the replacements x m ! x Ã m in the right hand side and x Ã m given by Eq (7). The variance of the active fraction is given by a , which is also computed using Eq (17). In terms of the covariances σ mn , this is given by

Ultrasensitivity in the mean activity
We will now explore the large A 0 limit of the fixed point. The fixed point values x Ã 0 and x Ã 2 given by Eq (7) can be expressed alternatively as In order to understand the behavior of the above expressions in the A 0 ! 1 limit, let us conjecture expansions of the form Consider now the expressions for R f and B f as given in Eq (5); after using the expansions given by Eq (20), we find the following asymptotic forms in the large A 0 limit: where The ratio R f /B f , upto Oð1=A 0 Þ, is given by Substituting Eq (23) in Eq (19), we find which give the leading terms in the expansions in Eq (20), and are valid for arbitrary ℓ. In order to take the analysis further, we study the limits ℓ ! 0 and ℓ ! 1 separately.

Small ℓ expansion.
Consider now the limit of small ligand concentrations, ℓ ( 1. We assume that, in this case, the zero th order terms x ð0Þ m in Eq (24) can be expanded as From Eq (22), we have Substituting the small ℓ expansions of x ð0Þ m as given in Eq (25), we obtain Substituting Eq (27) in Eq (24) leads to the equations is the control parameter that characterizes ZOU. The implications are (a) the population in the highest methylation state is vanishingly small in the limits ℓ ! 0 and A 0 ! 1 (b) x 00 0 ¼ 0 or 1 unless α = 1: this points to a jump-like behavior for ξ 0 (and therefore, ξ 1 , considering that x 00 0 þ x 00 1 ¼ 1) in these limits. The 'critical point' α c = 1 is the same as what was derived in the seminal paper of Goldbeter and Koshland [1], and referred to as the GK point in our earlier paper [2].
1.4.2 Large ℓ expansion. We next consider a large ℓ expansion of the zero th order term x ð0Þ m of large A 0 expansion in the following form Substituting these in Eq (26), we obtain Therefore, after substitution of Eq (31) in Eq (24), we find which mirror Eq (28) derived in the opposite, ℓ ! 0 limit. In the present case, it is ξ 2 that displays the jump transition from 1 to 0 as α crosses 1 (and ξ 1 does the reverse), while ξ 0 is vanishingly small at all α. These analytical results, given in Eqs (28) and (32) are supported by numerical simulations, to be discussed in the following section.

Stochastic gillespie simulations
For M = 2, we simulated the reaction scheme given in Fig 1 using Gillespie algorithm [23]. The specific numerical values for various parameters correspond to the methylation-demethylation reactions of chemotaxis receptors in the bacterium E. coli, previously used by various authors [2,15]. Initially, we choose all the receptors to be in inactive, unbound configuration at the lowest methylation/inactive level (m = 0). The system then evolves with the preassigned parameters (refer to Table 1) and eventually reaches its steady state where the mean receptor number in each methylation level becomes time-independent. BL model with M = 2 comprises of 8 conformational states and the total number of possible biochemical reactions is 9 (5 reversible and 4 irreversible reactions). Similarly for M = 3, there are 12 configuration states and 14 biochemical reactions while for M = 4, these numbers become 16 and 19 respectively. Once it is ensured that the system has reached its steady state, we study the average and variance of the active receptor fraction ξ a by varying R 0 , N, L and finally M, while B 0 was kept constant throughout. The fractions of receptors in different methylation levels were also kept track of.

Simulation results for M = 2
For fixed A 0 , the mean active fraction of receptors ξ a undergoes a sharp rise as R 0 is increased, as shown in Fig 2. ξ a is independent of ℓ, as is clear from the figure, and agrees with the analytical prediction (fixed point). The inset shows the variance in activity, which depends on ℓ in a non-monotonic manner (see Fig 3). In all cases, the maximum of the variance occurs close to the point of the steepest rise in the mean. Note also that for A 0 = 5.3μM (Fig 2a), the rise is less steeper compared to A 0 = 13.6μM (Fig 2b); for same ℓ, the variance is more in (a), but has a sharper peak in (b). For fixed ℓ, and changing A 0 , the mean shows the expected ultrasensitive rise as function of R 0 (Fig 4). For the set of parameters used here in Table 1, the critical value for R 0 corresponding to α = 1 Eq (29) is R c 0 ¼ 0:224mM, which agrees with our observations in Fig 4. The variance (Fig 4 inset) develops a sharper peak for larger A 0 . However, note that, unlike our earlier observations [2] in the GK model (M = 1) with only fully inactive and fully active states, the peak value of the variance here does not seem to increase with A 0 .
Why would s 2 a show a non-monotonic variation with ℓ? To find out, we studied the mathematical expressions for the individual variances σ 00 , σ 22 and the covariance σ 02 as functions of R 0 and ℓ, at fixed A 0 . The plots, shown in Fig 5, indicate that the non-monotonicity originates from the covariance σ 02 between the fully inactive and fully active methylation levels. Both σ 00 (Fig 5a) and σ 22 (Fig 5b) have their peaks near R c ; but while the peak value of σ 00 is a decreasing function of ℓ, that of σ 22 is an increasing function of ℓ. σ 02 (Fig 5c and 5d) is negative throughout and has its minimum near R c ; however, in (c) the lowest value changes non-monotonically with ℓ. From Eq (17), it is seen that s 2 a is dominated by σ 00 for small ℓ, while the large ℓ behavior is dominated by σ 22 . In the intermediate ℓ regime, σ 02 also contributes significantly, and therefore the non-monotonic variation of its minimum value with ℓ affects the peak value of s 2 a , which appears to be minimized near ℓ = 2. Next, we try to understand how the receptors distribute themselves among different methylation levels and the roles of ℓ and A 0 in this distribution. For fixed A 0 and very small ℓ, the fraction of receptors in the lowest methylation level (always inactive) state shows an ultrasensitive transition from near-unity to near-zero across R c . With increase in ℓ, the transition becomes less sharper/smoother leading to ξ 0 being almost zero in the entire range of R 0 for  [2] and [15]). The experimental values for the concentration of CheR in E. coli can be found in Table 2 in Appendix 3 but in our simulations, we have varied it in the range 10 −2 − 10μM for the sake of exploring ZOU. For the same reason, A 0 is also varied in the range 5.3 − 27.2μM.  Fig 6c and 6d. With increase in A 0 , the transitions in " x 0 , " x 1 and " x 2 become sharper (Fig 7), and agrees with the predictions   (28) and (32). Therefore, in the ℓ ! 0 as well as ℓ ! 1 regimes, the system effectively separates into two different GK-like two-state modules, with the population mostly shared between m = 0 and m = 1 states in the first case, and between m = 1 and m = 2 in the second case. Next, we investigated whether the reduction to two independent two-state models which was found for M = 2 works well for higher M as well. In Fig 9a-9d, we show, for M = 3, the four mean populations "

Symbol
x 0 ; " x 1 ; " x 2 and " x 3 respectively for a limited range of ℓ and two different values of A 0 . Similar to our observations in M = 2 model, we see that the system effectively reduces to two different 2-state (M = 1) models. For low ℓ (e.g. ℓ = 20), most of the receptors are in m = 0 or m = 3 states; in fact, " x 0 and " x 3 show near-ultrasensitive fall and rise respectively across R c . For high ℓ (e.g. ℓ = 20,000), while both " x 0 and " x 1 become smaller as A 0 is increased, " x 2 and " x 3 become dominant; now, " x 2 falls abruptly as R 0 crosses R c , while " x 3 rises in a similar Ultrasensitivity and fluctuations in Barkai-Leibler model way. Therefore, for M = 3 as well, we observe that the system splits into two different GK modules in the limits ℓ ! 0 and ℓ ! 1. The recent work by Pontius et al. [12] also studied mean and fluctuations in receptor activity in E. coli by varying the enzyme ratio [CheR]/[CheB]. While their model includes more features of chemotaxis receptors like clustering, allosteric interactions and enzyme brachiation, it also suffers from a few drawbacks, in our opinion. (a) In their analytical studies, the total methylation level in a cell is treated as the fundamental stochastic variable, whose dynamics is described by a phenomenological Langevin equation obtained by invoking the LNA; by contrast, our model has the receptor populations in various methylation levels as the basic variables, whose joint probability distribution follows a multivariate LFPE. Attractant concentration was not included as a parameter in their model, while it is explicitly included in ours, and its effect on the ZOU transition has been explored in detail.

Change in mean kinase activity in response to a change in attractant concentration
Bacteria perform chemotaxis by responding to time-dependent changes in attractant/repellent concentration in its immediate environment. In this section, we explore the response characteristics of our 3-state model to a time-dependent change in attractant concentration, within linear response theory, applicable to weak perturbations. Let us consider a situation when a time-dependent change δL(t) in the external ligand concentration is switched on in the extracellular environment of the bacterium at t = 0, such that L(t)!L + δL(t) for t ! 0, while L(t) = L for t < 0. Because ligand binding renders a receptor in the intermediate methylation level inactive, there is a change in the mean net activity at later times, which may be expressed in  where χ a (t) is the response function for mean receptor activity, which can be computed by submitting the system to a (small) step-like change in ligand concentration; δL(t) = δL s Θ(t) where Θ(t) is the Heaviside theta-function and δL s is the size of the step. Let d " x a ðstepÞ ðtÞ be the corresponding response. From Eq (33), it follows that From Eq (17), we find d " where d " x 0 ðtÞ and d " x 2 ðtÞ are the changes in " x 0 and " x 2 , respectively, in response to a change δL (t) in the ligand concentration, and are evaluated using Eq (6). The resulting equations have the form where the coefficients β mn have been defined earlier, following Eq (10), while γ m = @v m /@L| ξ Ã . Explicit expressions for both β mn and γ m are given in Appendix 2.
It is convenient to subject Eqs (35) and (36) to Laplace transforms, after putting δL(t) = δL s Θ(t), and use Eq (34) to compute the response function. For the Laplace transform of χ a (t), we findw a ðsÞ ¼ À wherew m ðsÞ are Laplace transforms of χ m (t) (m = 0, 2) which are defined through the relations d " x m ðtÞ ¼ R t 0 w m ðt À t 0 ÞdLðt 0 Þdt 0 . The explicit expressions turn out to bẽ where d = (s − β 00 )(s − β 22 ) − β 02 β 20 . We have confirmed that the response function curve computed in Eq (37) encloses zero area, i.e.,w a ð0Þ ¼ 0, which is a manifestation of the perfect adaptation (regaining the pre-stimulus state within few seconds) property of the bacterium, reproduced in the BL model (in our notation, this property implies lack of dependence of " x a on ℓ in steady state, recall Fig 2). In Appendix 3, we use the results in Eqs (37) and (38) to derive the chemotactic response function, which determines the fractional change in clockwise bias of the flagellar motor in response to temporal variations in attractant concentration, and consequently the chemotactic drift.

Discussion
The subject of the present study is the BL model of conformational dynamics of chemotaxis receptors in Escherichia coli. According to this model, a receptor undergoes stochastic switching between its active and inactive states, depending on its level of methylation and ligandation. The methylation and demethylation of receptors is carried out by two enzymes, viz. methyltransferase CheR and methylesterase CheB, acting antagonistically; the overall structure of the circuit is similar to the two-state (e.g. active/inactive) covalent modification scheme studied originally by GK [1]. This system is characterized by ZOU, a switch-like transition between all-inactive to all-active phases, in the limit of infinitely large substrate concentration. In the present paper, we investigated ZOU in the BL model within a stochastic formulation where the number N of receptors is treated finite. For analytical tractability, we first study a model with a single intermediate partially active state sandwiched between the fully inactive and fully active states. We show rigorously that the mean activity shows ultrasensitive response to increase in [CheR] (at fixed [CheB]) in the limit N ! 1, independent of the attractant concentration L. At the same time, its variance s 2 a is a non-monotonic function of [CheR] and L, as deduced from linear-noise approximation (LNA) and confirmed in numerical simulations. Interestingly, for very low and very high L, the system effectively reduces to two different "two-state" modules, akin to GK switches. The analysis is extended to models with more intermediate methylation levels (E. coli receptor has 3) via numerical simulations. By and large, qualitatively similar behavior is found in all the cases. For E. coli parameters, we find that σ a ' 10 −2 in the vicinity of the critical point of ZOU, for a large range of L, spanning almost 8 orders of magnitude. This estimate agrees with another recent computational study where explicit receptor clustering was included, and receptor activation-inactivation dynamics was stimulated using an equilibrium MWC model [12]. In Appendix 3, we show that this is also consistent with experimental measurements of [CheYp] fluctuations.
Biochemical noise is likely to be a crucial factor in the motility of E. coli for the following reason. The transition from run to tumble mode of a flagellar motor in E. coli is brought about by a switch in its direction of rotation, from counter-clockwise (CCW) to clockwise (CW). The latter switching is known to occur in an ultrasensitive fashion as a function of [CheYp], the concentration of phosphorylated CheY, which is the response regulator cytoplasmic protein (the phosphorylation being done by active receptors). The ultrasensitive response of flagellar motor for the CW bias with the variation in the concentration of CheYp, was first reported by Cluzel et al. [17] by performing in vivo experiments in single cells in absence of stimulant (attractant/repellent). The authors observed cell to cell variability in the amount of CheYp which was distributed around a mean value (ranging from 0.8 to 6μM) even for a given concentration of the inducer isopropyl β-D-thiogalactoside, in short IPTG (used to express CHEY-GFP in PS2001 E. coli strain). CW bias in various cells, pre-induced with different inducer levels, collapsed onto a single sigmoid curve when plotted against CheYp. The measured CW bias in the range 0.1 to 0.9 could be fitted to a Hill plot with Hill coefficient H * 10.3 ± 1.1, with the sensitive part being less than 3μM in width. The ability of the cell to execute runs and tumbles in succession depends on placing the intracellular [CheYp] somewhat precisely in this window. Cluzel et al. proposed that an additional molecular mechanism such as CheZ acting as the [CheYp] regulator can help in retaining the concentration of the latter in the narrow sensitive window of operation. Later, Yuan et al. [24] found from their experiments with cheRcheB cells that the flagellar motor can partially adapt to changes in [CheYp] over a time scale of several minutes, and thereby laterally shift its operating range by about 0.5μM. Any likelihood of partial adaptation due to dynamic localization of CheZ (suggested by Lipkow [25]) was ruled out by working with cheRcheBcheZ cells and explicitly showing similar results in both the cases. In 2013, Yuan and Berg [16] carried out a new set of experiments with pre-adapted motors, which showed that the response curve in this case is much steeper, compared to adapted motors (as in [17]) with Hill coefficients *16.5±1.1 and 20.7±1.6 for CW biases of 0.8±0.1 and 0.5±0.1 respectively. Because of the higher Hill coefficient, the sensitive window of [CheYp] is reduced to about 1μM in this case, with center around 3μM.
What is the impact of receptor-level noise on the clockwise bias of a single cell? It is clear that spontaneous fluctuations in receptor activity will produce corresponding changes in [CheYp], which would, in turn, affect the CW bias of the motor(s). In Appendix 4 we estimate that the maximum standard deviation in clockwise bias (arising from spontaneous fluctuations in [CheYp]) δP CW lies in the range 0.15-0.36 for Hill coefficient H = 20 [16], with the lower value corresponding to mean CW bias * 0.5 and the upper value for mean CW bias * 0.1. It therefore appears that intrinsic noise could rescue a cell with mean [CheYp] outside the optimal range from being forced to "run forever" without tumbling. The motor-level adaptation modifies these estimates a little, but our detailed analysis presented in Appendix 5 shows that this effect is negligible when [CheYp] lies in the range 2.5μM − 4μM, which more than covers the sensitive window. Outside this range, motor-level adaptation increases the standard deviation in CW bias, (see Appendix 5). We conclude, therefore, that intrinsic noise plays an important role in the run and tumble behavior of E. coli, and by extension, in chemotactic drift as well. A detailed analysis of the latter is presently being carried out and will be reported in the near future.
Being a study focused primarily on exploring the occurrence of ZOU in the BL model, we have not included all the features of signal transduction in E. coli here, even at the level of receptor dynamics. In particular, clustering of receptors and allosteric interactions among them have not been included here, but have been studied recently by other authors [12,26,27]. In [12], the authors show via numerical simulations that receptor clustering and enzyme localization make the mean activity a less sensitive function of [CheR], and thus makes the network robust to variations in protein numbers. At the same time, enzyme localization also leads to larger fluctuations in activity, with almost 10-fold increase in variance. Intrinsic biochemical noise arising from spontaneous fluctuations in receptor activity has been shown to be important in explaining the kinetics of flagellar motor-switching [28,29]. Using a set of parameters extracted from the literature, we show in Appendix 3 that the variance in activity measured in our simulations (as well as predicted using LNA) appear to be sufficiently large in magnitude to generate the [CheYp] fluctuations observed in experiments [28]. It is therefore likely that receptor clustering and enzyme localization may not be necessary to enhance the noise levels, although it could be important in other aspects of signaling, e.g. signal amplification.

Appendix 1 Intra and inter-module dynamics in the 3-state receptor-level biochemical reactions
To evaluate the number of receptors in each of the configurations, we employ sQSSA. The essence of sQSSA lies in the assumption that all the modules are weakly coupled (ν r and ν b assumed to be very small compared to other rates involved in the system) and hence detailed balance exists within each module. Under this condition, the mean numbers of intermediate complexes can be expressed in terms of their unbound counterparts as where R f and B f are the mean concentrations of free/unbound R and B enzymes respectively. After using the module-wise normalization conditions In the last equation, we have used the assumption of infinitely fast ligand binding and dissociation kinetics. Inter-module equilibrium is expressed through the following flux-balance relations which determine the fixed point ξ Ã : The free enzyme concentrations R f and B f are determined by the normalization conditions Using Eq (40) in Eq (42) leads to Eq (5). Perfect adaptation (insensitivity of average activity to L). Although the fixed point ξ Ã depends on L, the mean total active fraction does not; this can be easily shown as follows. The Eq (41) for intermodule dynamics, after summing the left hand side give After substituting the equations for intra-module dynamics Eq (40) (in the limit R f ( K r , B f ( K b ) and the expression for ξ a in terms of ξ m Eq (17) in Eq (43), the following quadratic equation for x Ã a , independent of L, is obtained: The solution to the above equation is

Appendix 2 Evaluation of various parameters
The expressions of β mn can be obtained by differentiating Eq (3). Note that the expressions for v m contain R f and B f which are in turn functions of ξ m . Therefore, To evaluate the linear response function, we have to consider the derivative γ m , the explicit expression of which is The relevant quantities involving derivatives with respect to L are listed below: In the computation of the coefficients β mn and γ m , all the above mentioned quantities are to be evaluated at the fixed point (ξ Ã , R f (ξ Ã ), B f (ξ Ã ) and it has been assumed that R f ( K r and B f ( K b .

Appendix 3 Chemotactic response function
In E. coli, the receptor complex acts as the processing unit for any input signal (e.g. changes in the extracellular environment of the bacterium) while the flagellar motor, which controls the switching kinetics of the flagella, regulate the output (run/tumble motion). While the CCW mode of rotation of flagella corresponds to straight motion of the bacterium, the CW mode corresponds to tumbles. These two key constituents of the signaling network are connected by the response regulator CheY. Attractant binding to the receptors results in suppression of the activity of the latter, hence leads to less phosphorylation of CheY. Phosphorylated CheY (CheYp) binds to flagellar motors and enhances the rate of CCW!CW switching, hence a temporal increase in attractant concentration as sensed by the bacterium would result in elongated run durations when swimming in favorable directions. The dynamics of phosphorylated CheY may be described by the equation where a Y is rate for binding of (non-phosphorylated) CheY to an active receptor, which leads to its phosphorylation, Y 0 is the total concentration of CheY (henceforth assumed much larger than Y(t)) and λ Y is the rate of degradation of CheYp (due to the dephosphatase CheZ). The steady state mean concentration of CheYp can be obtained from Eq (53) by equating left hand side to 0, In response to a change in attractant concentration L ! L + δL(t) for t ! 0, the mean [CheYp] undergoes a shift Y Ã ! Y Ã + δY(t), with δY(t) given by The clockwise bias, i.e., the probability P CW for a flagellar motor to be in clockwise-spinning state (corresponding to tumble) has been found in experiments to be an ultrasensitive function of Y [13,16,17]: where H is the Hill coefficient symbolizing the steepness of the transition. It follows that fractional changes in the CW bias P CW and Y are related as where, in the second part, we have defined the response function χ b (t) for the bias. Using Eqs (55) and (54) in Eq (57), it follows that, in steady state, χ b (t) is related to χ a (t) through From Eq (58), it follows that the Laplace transforms of χ b (t) and χ a (t) are related linearly: w b ðsÞ /w a ðsÞ=ðs þ l Y Þ, hencew b ð0Þ ¼ 0, i.e., the response curve for the bias too encloses zero area [9,15]. After using Eqs (37) and (38) in Eq (58), and performing the inversion, we arrive at the following multi-exponential form for χ b (t): The rate constants A and B, as well as the coefficients X m , Y m and Z m (m = 0, 2) are expressed in terms of β mn and γ m as follows: In Fig 10, we show plots of the mathematical expression for χ b (t) obtained in Eq (59), (a) by varying B 0 and (b) varying ℓ. The curves closely resemble the experimentally measured responses to short-lived stimuli [9,30], for somewhat large values of B 0 , in agreement with earlier results in a mean-field BL model [15]. In (a), with increase in B 0 , there is an overall reduction in time scales and an increase in the depth of the negative lobe. In (b), with increase in ℓ, the depth of the negative lobe decreases with no change in the time scales. For ℓ ) 1, the negative lobe effectively disappears.
(assuming " y %ỹ). So there is a maximum fluctuation of 18% for H = 20 corresponding to a CW bias of 0.5 and [CheYp] *3μM around which the transition becomes ultrasensitive. Therefore, the fluctuations in [CheYp] appear sufficiently large to affect the switching behavior of the motor.
Next, how large the fluctuations in activity need to be, to produce the observed standard deviation in [CheYp]? The following simple picture helps us to obtain an estimate. In E. coli, active receptors undergo autophosphorylation, and then transfer the phosphoryl groups to the proteins CheB and CheY. For simplicity, we ignore phosphorylation of CheB and focus on CheY. The dynamics of Y(t) is expressed in Eq (53). Under steady state conditions, the variance s 2 Y ¼ hY 2 i À hYi 2 is given by the expression where, we have introduced the correlation length λ ξ for the fluctuations in ξ a , defined through the relation hdx a ðtÞdx a ðt 0 Þi ¼ s 2 a exp ðÀ l x jt À t 0 jÞ in steady state, where dx a ðtÞ ¼ x a ðtÞ À " x a is the fluctuation in activity relative to its steady state value. From Eq (65), the inequality follows, the equality being satisfied in the limit λ ξ ( λ Y (which we shall assume henceforth). Various estimates for the parameters in Eq (66) are available in the literature (see Table 2). Using Eq (64), we identify the following ranges for the standard deviation σ a of receptor activity: 4:7 Â 10 À 3 < s a < 10 À 2 ðiÞ 3 Â 10 À 4 < s a < 7:4 Â 10 À 4 ðiiÞ 5:2 Â 10 À 4 < s a < 10 À 3 ðiiiÞ It is also of interest to compare the values of the ZOU parameter α across these sets of parameters, which turn out to be 0.547 (i), 0.076 (ii) and 0.035 (iii) ( Table 2). Note that the estimated variance in (ii) and (iii) are smaller than (i), which is consistent with our expectation of maximum variance near α = 1, over a large range of values of ℓ. The magnitude of receptor noise predicted by Eq (67) also agrees with our predictions (see Fig 3, also Fig 8).

Appendix 5 Effect of flagellar motor level adaptation on fluctuations in clockwise bias
From the experimental data of Yuan and Berg (2013) [16], we propose that the change in clockwise bias in response to a time-dependent change in [CheYp] (denoted δY(t) henceforth) be written in the general linear form where χ m (τ) (defined for τ > 0) is the linear response function for the motor, which we assume to consist of two terms: where a na ¼ @ Y P na CW , and P na CW is the bias of the non-adapted motor, as measured by Yuan and Berg (2013), with Hill coefficient H % 20. The second term represents the slow change in the bias originating from motor-level adaptation, with l À 1 m being the time scale for the same. Now, for a step-like change in Y such that Y ! Y 0 + ΔY at t = 0, let us put δY(t) = ΔYΘ(t); using this in Eq (68) yields the response for the same: As t ! 1, the system will adapt, and hence dP step CW ðtÞ ! a a DY, where a a ¼ @ Y P a CW , and P a CW is the clockwise bias of the adapted motor. From Eq (69), it then follows that Γ = λ(α a − α na ). Let us now compute the effect on the CW bias of random fluctuations in Y, and define the rootmean square change of bias by dP CW where, in steady state, we expect hdYðt 1 ÞdYðt 2 Þi ¼ s 2 Y e À l Y ðt 1 À t 2 Þ for t 1 > t 2 , s 2 Y is the steady state variance of Y as given by Eq (65) and λ Y is the dephosphorylation rate of [CheYp]. Using the expression for χ m (t) from Eq (69) and completing the calculations yields the final result The term outside the square root is precisely the result in Eq (62) in Appendix 3. Explicit evaluation of the term in brackets (see Fig 11) shows, however, that the contribution of motor-level adaptation to the fluctuations in bias is negligible in the sensitive part of the response curve.