## Figures

## Abstract

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 in the active fraction, which is estimated mathematically under linear noise approximation (LNA). It is found that 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.

**Citation: **Roy U, Gopalakrishnan M (2017) Ultrasensitivity and fluctuations in the Barkai-Leibler model of chemotaxis receptors in *Escherichia coli*. PLoS ONE 12(4):
e0175309.
https://doi.org/10.1371/journal.pone.0175309

**Editor: **Ramon Grima,
University of Edinburgh, UNITED KINGDOM

**Received: **October 18, 2016; **Accepted: **March 23, 2017; **Published: ** April 13, 2017

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

**Data Availability: **All relevant data are within the Manuscript.

**Funding: **The author(s) received no specific funding for this work.

**Competing interests: ** The authors have declared that no competing interests exist.

## 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–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–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 run-and-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.

The red boxes denote the three states; *ξ*_{0}, *ξ*_{1} and *ξ*_{2} being the fractional concentrations of receptors in each. Blue dotted arrows with rates *ν*_{r} show the methylation process resulting in the increment in the state index while the green ones with rates *ν*_{b} depict demethylation with the decrement in the same.

## Materials and methods

### 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 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 and 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 , 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} = ∞ (these are ideal values). We further choose [19], a simple mathematical form which satisfies our earlier assumptions.

## 1 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 . 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} and . For *m* = 1, there is a further sub-division: and , as explained in the previous section. For each *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 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]: (1) where

**J**≡ {

*J*

_{0},

*J*

_{2}} with (for

*m*= 0, 2) denoting the components of the state space probability current “vector”. Here, the drift and diffusion coefficients are given by (2) with

*δ*=

*N*

^{−1}, , , and .

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:
(3)
and
(4)

The free enzyme concentrations *R*_{f} and *B*_{f} are also functions of *ξ*_{0}, *ξ*_{2}, and the expressions are given by (see Appendix 1)
(5)

## Results

### 1.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
(6)
where 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 . 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
(7)
where and 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*

**ξ****.*

**ξ**### 1.2 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 * ξ*′ =

*−*

**ξ****; we also define its probability distribution Φ(*

**ξ***′)≡*

**ξ***P*(

** +*

**ξ***′), which satisfies the equation (8) where , and may be expanded as follows: (9) Keeping the leading term in each, and noting that , we arrive at (for*

**ξ***m*= 0, 2) (10) where , . Substitution of Eq (10) in Eq (8) leads to the multivariate linear Fokker-Planck equation (LFPE) [22] (11)

Let us now define the covariances 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)
(12)
where ** μ** = (

*μ*

_{0},

*μ*

_{2}). The moment generating function admits the Taylor expansion (13)

Since Φ(* ξ*′) satisfies the LFPE given by Eq (11), the moment generating function satisfies the equation
(14)

Under stationary conditions (*t* → ∞), the left hand side of Eq (14) becomes zero; next we substitute Eq (13) in Eq (14), and find that 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]
(15)
where are the elements of the (diagonal) diffusion matrix * D*. The following expressions follow from Eq (15):
(16)
where . Explicit expressions for the coefficients

*β*

_{mn}are to be found in Appendix 2.

### 1.3 Mean and fluctuations in activity

The active fraction of receptors in the present model is given by
(17)
where *ℓ* = *L*/*K*_{L}. Within LNA, the average fractions , hence , where is given by Eq (17), with the replacements in the right hand side and given by Eq (7).

The variance of the active fraction is given by , which is also computed using Eq (17). In terms of the covariances *σ*_{mn}, this is given by
(18)

### 1.4 Ultrasensitivity in the mean activity

We will now explore the large *A*_{0} limit of the fixed point. The fixed point values and given by Eq (7) can be expressed alternatively as
(19)

In order to understand the behavior of the above expressions in the *A*_{0} → ∞ limit, let us conjecture expansions of the form
(20)

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:
(21)
where
(22)

The ratio *R*_{f}/*B*_{f}, upto , is given by
(23)

Substituting Eq (23) in Eq (19), we find
(24)
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 *ℓ* → ∞ separately.

#### 1.4.1 Small *ℓ* expansion.

Consider now the limit of small ligand concentrations, *ℓ* ≪ 1. We assume that, in this case, the *zero*^{th} order terms in Eq (24) can be expanded as
(25)

From Eq (22), we have (26)

Substituting the small *ℓ* expansions of as given in Eq (25), we obtain
(27)

Substituting Eq (27) in Eq (24) leads to the equations
(28)
where
(29)
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} → ∞ (b) or 1 unless *α* = 1: this points to a jump-like behavior for *ξ*_{0} (and therefore, *ξ*_{1}, considering that ) 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 of large *A*_{0} expansion in the following form
(30)

Substituting these in Eq (26), we obtain (31)

Therefore, after substitution of Eq (31) in Eq (24), we find
(32)
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.

## 2 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).

(previously used in [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*.

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.

### 2.1 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 , 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}.

The figures show the results for two values of *A*_{0}, 5.3*μ*M (a) and 13.6*μ*M (b). In both the figures, the insets show the corresponding variances. Note that the maximum of the variance occurs at the same *R*_{0} for various *ℓ*.

Interestingly, the maximum is a non-monotonic function of *ℓ*, and appears to have a minimum around *ℓ* = 2.

Why would 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 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 , which appears to be minimized near *ℓ* = 2.

For fixed *A*_{0}, near the critical point, *σ*_{00} is a decreasing function of *ℓ*, while *σ*_{22} is an increasing function of the same. By contrast, the covariance *σ*_{02} (c) displays non-monotonic change with *ℓ* near the critical point, with a minimum reached near *ℓ* = 2 (compare with Fig 3). But *σ*_{02} (d) varies monotonically with *A*_{0} for fixed *ℓ*.

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 very large *ℓ*, as given in Fig 6a. On the other hand, the fraction in the highest methylation level (always active) state is almost zero for very small *ℓ* and then increases in the upper-critical regime with increase in *ℓ*, finally leading to an ultrasensitive switch from near-zero to near-unity for very large *ℓ* Fig 6b. But the fraction in the intermediate methylation level exhibits switch-like behaviors of opposite nature at very low and high *ℓ* values; while for intermediate *ℓ*, it shows non-monotonic change with *R*_{0} with a peak near *R*_{c} Fig 6c and 6d. With increase in *A*_{0}, the transitions in , and become sharper (Fig 7), and agrees with the predictions made in Eqs (28) and (32). Therefore, in the *ℓ* → 0 as well as *ℓ* → ∞ 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.

### 2.2 Simulation results for *M* > 2

Numerical simulations for models with more intermediate methylation levels, namely *M* = 3 and *M* = 4 were also carried out for a limited range of values of *ℓ* and *A*_{0}. Fig 8 shows a comparison of the mean and variance in active fraction between models with *M* = 2, 3 and 4, for *ℓ* = 20 and *A*_{0} = 13.6*μ*M. For larger *M*, the mean shows sharper rise near *R*_{0} = *R*_{c}, but the peak value of the variance for *M* = 4 is only about a third of the corresponding number for *M* = 2. In general, models with larger number of intermediate states has smaller fluctuations in activity, presumably due to the dominant role of inter-state covariances.

The attractant concentration is *ℓ* = 20. The inset displays the corresponding variances. The numerical results establish that the qualitative behavior of the BL model, in particular, the ZOU transition, remains unchanged in presence of additional intermediate methylation levels.

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 and 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, and show near-ultrasensitive fall and rise respectively across *R*_{c}. For high *ℓ* (*e.g.* *ℓ* = 20,000), while both and become smaller as *A*_{0} is increased, and become dominant; now, falls abruptly as *R*_{0} crosses *R*_{c}, while rises in a similar way. Therefore, for *M* = 3 as well, we observe that the system splits into two different GK modules in the limits *ℓ* → 0 and *ℓ* → ∞.

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. (b) No direct quantitative comparisons between model predictions and stochastic simulations were done in their paper, while our model predictions are shown to agree well with simulations. (c) 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.

### 2.3 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 the general form
(33)
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 be the corresponding response. From Eq (33), it follows that
(34)

From Eq (17), we find
(35)
where and are the changes in and , respectively, in response to a change *δL*(*t*) in the ligand concentration, and are evaluated using Eq (6). The resulting equations have the form
(36)
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 find
(37)
where are Laplace transforms of *χ*_{m}(*t*) (*m* = 0, 2) which are defined through the relations . The explicit expressions turn out to be
(38)
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., , 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 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* → ∞, independent of the attractant concentration *L*. At the same time, its variance 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
(39)
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 , , and , we find
(40)

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 * ξ**:
(41)

The free enzyme concentrations *R*_{f} and *B*_{f} are determined by the normalization conditions
(42)

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 (43)

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 , independent of *L*, is obtained:
(44)

### 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,
(46)
where
(47) (48) (49)
To evaluate the linear response function, we have to consider the derivative *γ*_{m}, the explicit expression of which is
(50)
The relevant quantities involving derivatives with respect to *L* are listed below:
(51) (52)
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
(53)
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,
(54)
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
(55)

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]:
(56)
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
(57)
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
(58)
From Eq (58), it follows that the Laplace transforms of *χ*_{b}(*t*) and *χ*_{a}(*t*) are related linearly: , hence , 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*):
(59)
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:
(60) (61)
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.

### Appendix 4 On the role of biochemical noise in the run and tumble motion of *E. coli*

From Eq (56), it follows that a fluctuation in [CheYp] equal to its standard deviation in magnitude would cause the following change in *P*_{CW}:
(62)
where *y* = *Y*/*K*_{Y} and is its mean value in steady state. The function has a single maximum (and no minimum) at , with maximum value 1/4 for *H* ≫ 1. Using this property in Eq (62) leads to the following upper bound on *δP*_{CW}:
(63)

The standard deviation *σ*_{Y} has been measured in experiments, and the range of variation has been estimated [28] as
(64)
with the actual value depending on the mean clockwise bias, . The dissociation constant *K*_{Y} ≃ 3*μM* and the Hill coefficient *H* ≃ 20 (for pre-adapted motors) have also been measured in single cell experiments [16]. Using this estimate in Eq (63), we find *δP*_{CW} ≤ 0.15 − 0.367 (assuming ). 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 is given by the expression
(65)
where, we have introduced the correlation length *λ*_{ξ} for the fluctuations in *ξ*_{a}, defined through the relation in steady state, where is the fluctuation in activity relative to its steady state value. From Eq (65), the inequality
(66)
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:
(67)

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
(68)
where *χ*_{m}(*τ*) (defined for *τ* > 0) is the linear response function for the motor, which we assume to consist of two terms:
(69)
where , and 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 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:
(70)
As *t* → ∞, the system will adapt, and hence , where , and 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 root-mean square change of bias by , with the average to be carried out over different realizations of the fluctuations in [CheYp]. From Eq (68), it turns out that
(71)
where, in steady state, we expect for *t*_{1} > *t*_{2}, 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
(72)
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.

Note the flat region, nearly 1.5*μM* in width, where the correction arising from motor-level adaptation is negligible. *α*_{a} and *α*_{na} have been computed from Eq (56) with *H*_{a} = 10 [17] and *H*_{na} = 20 [16]. We have taken *λ*_{m} = 0.0167*s*^{−1} considering that motor adaptation takes place over a minute [24] and *λ*_{Y} = 30*s*^{−1} [14, 31].

## Acknowledgments

We thank P.G. Senapathy Centre for Computing Resources, IIT Madras, for providing us with computing facilities in the Virgo cluster.

## Author Contributions

**Conceptualization:**MG UR.**Data curation:**UR.**Formal analysis:**UR MG.**Investigation:**UR MG.**Methodology:**UR MG.**Project administration:**UR MG.**Resources:**UR MG.**Software:**UR.**Supervision:**MG.**Validation:**UR.**Visualization:**UR MG.**Writing – original draft:**UR MG.**Writing – review & editing:**UR MG.

## References

- 1. Goldbeter A, Koshland DE Jr. An amplified sensitivity arising from covalent modification in biological systems. Proc. Natl. Acad. Sci. 1981 73(11):6840–6844.
- 2. Jithinraj PK, Roy U, Gopalakrishnan M. Zero-order ultrasensitivity: A study of criticality and fluctuations under the total quasi-steady state approximation in the linear noise regime. J. Theor. Biol. 2014 344:1–11. pmid:24309434
- 3. Berg OG, Paulsson J, Ehrenberg M. Fluctuations and quality of control in biological cells: zero-order ultrasensitivity reinvestigated. Biophys. J. 2000 79:1228–1236. pmid:10968987
- 4. Elf J, Ehrenberg M. Fast evaluation of fluctuations in biochemical networks with the linear noise approximation. Genome Res. 2003 13:2475–2484. pmid:14597656
- 5. Xu Y, Gunawardena J. Realistic enzymology for post-translational modification: zero-order ultrasensitivity revisited. J. Theor. Biol. 2012 311:139–152. pmid:22828569
- 6. Qian H. Thermodynamic and kinetic analysis of sensitivity amplification in biological signal transduction. Biophys. Chem. 2003 105:585–593. pmid:14499920
- 7. Ge H, Qian M. Sensitivity amplification in the phosphorylation–dephosphorylation cycle: nonequilibrium steady states, chemical master equation, and temporal cooperativity. J. Chem. Phys. 2008 129:015104. pmid:18624503
- 8. Tindall MJ, Porter SL, Maini PK, Gaglia G, Armitage JP. Overview of mathematical approaches used to model bacterial chemotaxis I: the single cell. Bull. Math. Biol. 2008 70: 1525–1569. pmid:18642048
- 9. Segall JE, Block SM, Berg HC. Temporal comparisons in bacterial chemotaxis. Proc. Natl. Acad. Sci. USA 1986 83:8987–8991. pmid:3024160
- 10. Barkai N, Leibler S. Robustness in simple biochemical networks. Nature. 1997 387: 913–917. pmid:9202124
- 11.
Rao CV, Kirby JR, Arkin AP. Design and Diversity in Bacterial Chemotaxis: A Comparative Study in
*Escherichia coli*and*Bacillus subtilis*. PLoS Biol. 2004 2:E49. pmid:14966542 - 12. Pontius W, Sneddon MW, Emonet T. Adaptation dynamics in densely clustered chemoreceptors. PLoS Comp. Biol. 2013 9(9): e1003230
- 13. Mello D, Tu Y. Perfect and near-perfect adaptation in a model of bacterial chemotaxis. Biophys. J. 2003 84:2943–2956. pmid:12719226
- 14. Kollmann M, Lovdok L, Bartholome K, Timmer J, Sourjik V. (2005) Design principles of a bacterial signalling network. Nature 2005 438: 504–507
- 15. Reneaux M, Gopalakrishnan M. Theoretical results for chemotactic response and drift of E. coli in a weak attractant gradient. J. Theor. Biol. 2010 266:99–106. pmid:20558183
- 16. Yuan J, Berg HC. Ultrasensitivity of an adaptive bacterial motor. J. Mol. Biol. 2013 425: 1760–1764 pmid:23454041
- 17. Cluzel P, Surette M, Leibler S. An ultrasensitive bacterial motor revealed by monitoring signaling proteins in single cells. Science. 2000. 287: 1652–1655. pmid:10698740
- 18. Asakura S, Honda M. Two-state model for bacterial chemoreceptor proteins: the role of multiple methylation. J. Math. Biol. 1984 176:349–367. pmid:6748079
- 19. Emonet T, Cluzel P. Relationship between cellular response and behavioural variability in bacterial chemotaxis. Proc. Nat. Acad. Sci. USA 2008 105:3304–3309. pmid:18299569
- 20. Sourjik V, Berg HC. Receptor sensitivity in bacterial chemotaxis. Proc. Nat. Acad. Sci. USA 2002 99: 123–127. pmid:11742065
- 21. Schnoerr D, Sanguinetti G, Grima R. Approximation and inference methods for stochastic biochemical kinetics—a tutorial review. J. Phys. A: Math. Theor. 2016 10.1088: 1751–8121
- 22.
van Kampen NG. Stochastic Processes in Physics and Chemistry. Elsevier. 2007
- 23. Gillespie DT. Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 1977 81(25):2340–2351.
- 24. Yuan J, Branch RW, Hosu BG, Berg HC. Adaptation at the output of the chemotaxis signalling pathway. Nature 2012 484: 233–236 pmid:22498629
- 25. Lipkow K. Changing cellular localization of CheZ predicted by molecular simulations. PLoS Comput. Biol. 2006 2, e39 pmid:16683020
- 26. Endres RG, Wingreen NS Precise adaptation in bacterial chemotaxis through “assistance neighbourhoods”. Proc. Natl. Acad. Sci. USA 2006 103: 13040–13044 pmid:16924119
- 27. Hansen CH, Endres RG, Wingreen NS Chemotaxis in Escherichia coli: a molecular model for robust precise adaptation. PLos Comp. Biol. 2008 4: e1 pmid:18179279
- 28. Park H, Pontius W, Guet CC, Marko JF, Emonet T, Cluzel P. Interdependence of behavioural variability and response to small stimuli in bacteria. Nature 2010 468:819–823 pmid:21076396
- 29. Park H, Oikonomou P, Guet CC, Cluzel P. Noise underlies switching behaviour of the bacterial flagellum. Biophys. J. 2011 101(10): 2336–2340 pmid:22098731
- 30. Masson JB, Voisinnea G, Wong-Nga J, Celani A and Vergassola M Noninvasive inference of the molecular chemotactic response using bacterial trajectories. Proc. Natl. Acad. Sci. USA 2012 109: 1802–1807 pmid:22307649
- 31. Rao CV, Arkin AP. Stochastic chemical kinetics and the quasi-steady-state assumption: Application to the Gillespie algorithm. J. Chem. Phys. 2003 118: 4999–5010
- 32. Morton-Firth CJ, Shimizu TS, Bray D. A free-energy-based stochastic simulation of the Tar receptor complex. J. Mol. Biol. 1999 286: 1059–1074 pmid:10047482