Ultrasensitivity and Fluctuations in the Barkai-Leibler Model of Chemotaxis Receptors in {\it Escherichia coli}

TL;DR

This paper investigates how intrinsic noise affects the ultrasensitive response of E. coli chemotaxis receptors, revealing that fluctuations peak near the ultrasensitivity transition and can influence bacterial motility.

Contribution

It introduces a stochastic model of E. coli chemotaxis receptors demonstrating ultrasensitivity and quantifies receptor noise, linking it to bacterial behavior.

Findings

Receptor activity exhibits zero-order ultrasensitivity similar to Goldbeter-Koshland model.

Fluctuation variance peaks near the ultrasensitivity transition.

Receptor noise can positively influence bacterial motility.

Abstract

A stochastic version of the Barkai-Leibler model of chemotaxis receptors in {\it E. coli} is studied here to elucidate the effects of intrinsic network noise in their conformational dynamics. It was originally proposed to explain the robust and near-perfect adaptation of {\it E. coli} observed across a wide range of spatially uniform attractant/repellent (ligand) concentrations. A receptor is either active or inactive and can stochastically switch between the two states. Enzyme CheR methylates inactive receptors while CheB demethylates active ones and the probability for it to be active depends on its level of methylation and ligandation. A simple version of the model with two methylation sites per receptor (M=2) shows zero-order ultrasensitivity (ZOU) akin to the classical 2-state model of covalent modification studied by Goldbeter and Koshland (GK). For extremely small and large ligand concentrations, the system reduces to two 2-state GK modules. A quantitative measure of the spontaneous fluctuations in activity (variance) estimated mathematically under linear noise approximation (LNA) is found to peak near the ZOU transition. The variance is a weak, non-monotonic and decreasing functions of ligand and receptor concentrations. Gillespie simulations for M=2 show excellent agreement with analytical results obtained under LNA. Numerical results for M=2, 3 and 4 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 have a beneficial effect on the motility of the organism. The response of mean receptor activity to small time-dependent changes in the external ligand concentration, computed within linear response theory, is found to have a bilobe form.

Figures

Tables

Table 1. Table 1: A list of values of the various parameters used in numerical simulations (previously used in [ 2 ] and [ 15 ] ). The experimental values for the concentration of CheR in E.coli can be found in Table 2 in 6.3 but in our simulations, we have varied it in the range 10 − 2 − 10 ​ μ superscript 10 2 10 𝜇 10^{-2}-10\mu M for the sake of exploring ZOU. For the same reason, A 0 subscript 𝐴 0 A_{0} is also varied in the range 5.3-27.2 μ 𝜇 \mu M.

Symbol	Quantity	Estimated value
$V$	Cell Volume	${10}^{-15}$ L (602.3 $\mu$M-1)
$B_0$	CheB concentration	0.28 $\mu$M
$K_r$	CheR dissociation constant	0.39 $\mu$M
$K_b$	CheB dissociation constant	0.54 $\mu$M
$K_L$	Ligand binding constant	0.1 $\mu$M
${\nu }_r$	Methylation rate	0.75 $s^{-1}$
${\nu }_b$	Demethylation rate	0.6 $s^{-1}$

Table 2. Table 2: A list of relevant parameter values collected from the literature.

		Numerical values
Symbol	Quantity	Morton-Firth et al. [32]	Rao et al.[31]	Kollmann et al.[14]
$a_Y$	CheY phosphorylation rate	3 $\mu M^{-1}{s}^{-1}$	100 $\mu M^{-1}{s}^{-1}$	100 $\mu M^{-1}{s}^{-1}$
${\lambda }_Y$	CheYp dephosphorylation rate	14.15 $s^{-1}$	30.1 $s^{-1}$	30.1 $s^{-1}$
$Y_0$	CheY concentration	18 $\mu$M	17.9 $\mu$M	9.7 $\mu$M
$A_0$	CheA concentration	5 $\mu$M	5 $\mu$M	5.3 $\mu$M
$R_0$	CheR concentration	0.235 $\mu$M	0.3 $\mu$M	0.16 $\mu$M
$B_0$	CheB concentration	2.27 $\mu$M	2 $\mu$M	0.28 $\mu$M
${\nu }_r$	Methylation rate	0.819 $s^{-1}$	0.255 $s^{-1}$	0.39 $s^{-1}$
${\nu }_b$	Demethylation rate	0.155 $s^{-1}$	0.5 $s^{-1}$	6.3 $s^{-1}$
$\alpha \left(\text{<a><span>29</span></a>}\right)$	ZOU parameter	0.547	0.076	0.035

Equations

$$\frac{\partial P(\mathbf{\xi },t)}{\partial t}=-\frac{\partial }{\partial \mathbf{\xi }}\cdot \mathbf{J}$$

$$v_m=\delta ({\omega }_{1m}-{\omega }_{m1});D_m=\frac{{\delta }^2}{2}({\omega }_{1m}+{\omega }_{m1});m=0,2$$

$$v_0(\mathbf{\xi })$$

$$v_2(\mathbf{\xi })$$

$$D_0(\mathbf{\xi })$$

$$D_2(\mathbf{\xi })$$

$$R_f(\mathbf{\xi })$$

$$B_f(\mathbf{\xi })$$

$$\frac{d\overline{\xi }}{dt}=\overline{\mathbf{v}(\mathbf{\xi })},$$

$${\xi }_0^{\ast }=\frac{B_f^{\ast 2}{K}_r^2{\nu }_b^2{K}_L}{(B_f^{\ast }{K}_r{\nu }_b+R_f^{\ast }{K}_b{\nu }_r)(B_f^{\ast }{K}_r{\nu }_b{K}_L+R_f^{\ast }{K}_b{\nu }_{r}L)}$$

$${\xi }_2^{\ast }=\frac{R_f^{\ast 2}{K}_b^2{\nu }_r^{2}L}{(B_f^{\ast }{K}_r{\nu }_b+R_f^{\ast }{K}_b{\nu }_r)(B_f^{\ast }{K}_r{\nu }_b{K}_L+R_f^{\ast }{K}_b{\nu }_{r}L)},$$

$$\frac{\partial \Phi ({\mathbf{\xi }}^{\prime },t)}{\partial t}=-\frac{\partial }{\partial {\mathbf{\xi }}^{\prime }}\cdot {\mathbf{J}}^{\prime }$$

$$J^{\prime}_{m}(\boldsymbol{\xi}^{\prime})=\Phi(\boldsymbol{\xi}^{\prime})\Big{[}v_{m}(\boldsymbol{\xi}^{*})+\sum_{n}\xi^{\prime}_{n}\frac{\partial v_{m}}{\partial\xi_{n}}\Big{\rvert}_{\boldsymbol{\xi}^{*}}+...\Big{]}-\frac{\partial}{\partial\xi_{m}^{\prime}}\Big{[}\Phi(\boldsymbol{\xi}^{\prime})\Big{(}D_{m}(\boldsymbol{\xi}^{*})+\sum_{l}\xi^{\prime}_{l}\frac{\partial D_{m}}{\partial\xi_{l}}\Big{\rvert}_{\boldsymbol{\xi}^{*}}+...\Big{)}\Big{]}$$

$$J_m^{\prime }\simeq \Phi ({\mathbf{\xi }}^{\prime })\sum _{n=0,2}{\beta }_{mn}{\xi }_n^{\prime }-D_m^{\ast }\frac{\partial \Phi ({\mathbf{\xi }}^{\prime })}{\partial {\xi }_m^{\prime }}$$

$$\frac{\partial\Phi(\boldsymbol{\xi}^{\prime},t)}{\partial t}=-\sum_{m=0,2}\frac{\partial}{\partial\xi_{m}^{\prime}}\Big{[}\Phi\sum_{n=0,2}\beta_{mn}\xi^{\prime}_{n}\Big{]}+\sum_{m=0,2}D_{m}^{*}\frac{\partial^{2}\Phi}{\partial{\xi_{m}^{\prime}}^{2}},$$

$$G(\mathbf{\mu },t)={\int }_{-\infty }^{\infty }d{\xi }_0^{\prime }d{\xi }_2^{\prime }{e}^{-i\mathbf{\mu }\cdot {\mathbf{\xi }}^{\prime }}\Phi ({\mathbf{\xi }}^{\prime },t)$$

$$G(\mathbf{\mu },t)=1-i\mathbf{\mu }\cdot \overline{{\mathbf{\xi }}^{\prime }(t)}-\frac{1}{2}\sum _{m,n=0,2}{\mu }_m{\mu }_n{\sigma }_{mn}(t)+.........$$

$$\frac{\partial G(\mathbf{\mu },t)}{\partial t}=\sum _{m,n=0,2}{\beta }_{mn}{\mu }_m\frac{\partial G}{\partial {\mu }_n}-G\sum _{m=0,2}{D}_m^{\ast }{\mu }_m^2$$

$$\mathbf{\beta }\mathbf{\sigma }+\mathbf{\sigma }{\mathbf{\beta }}^T+2\mathbf{D}=0$$

$${\sigma }_{00}$$

$${\sigma }_{22}$$

$${\sigma }_{02}$$

$${\xi }_a={\xi }_1^a+{\xi }_2=\frac{(1-{\xi }_0)}{1+\ell }+\frac{\ell }{1+\ell }{\xi }_2.$$

$${\sigma }_a^2=\frac{{\sigma }_{00}+{\ell }^2{\sigma }_{22}-2\ell {\sigma }_{02}}{(1+\ell {)}^2}.$$

$$\displaystyle\xi_{0}^{*}=\left[1+\Big{(}\frac{R_{f}^{*}}{B_{f}^{*}}\Big{)}^{2}\frac{K_{b}^{2}\nu_{r}^{2}\ell}{K_{r}^{2}\nu_{b}^{2}}+\Big{(}\frac{R_{f}^{*}}{B_{f}^{*}}\Big{)}\frac{K_{b}\nu_{r}(\ell+1)}{K_{r}\nu_{b}}\right]^{-1}$$

$$\displaystyle\xi_{2}^{*}=\left[1+\Big{(}\frac{B_{f}^{*}}{R_{f}^{*}}\Big{)}^{2}\frac{K_{r}^{2}\nu_{b}^{2}}{K_{b}^{2}\nu_{r}^{2}\ell}+\Big{(}\frac{B_{f}^{*}}{R_{f}^{*}}\Big{)}\frac{K_{r}\nu_{b}(\ell+1)}{K_{b}\nu_{r}\ell}\right]^{-1}$$

$${\xi }_m^{\ast }={\xi }_m^{(0)}+\frac{1}{A_0}{\xi }_m^{(1)}+\frac{1}{A_0^2}{\xi }_m^{(2)}+...$$

$$R_f$$

$$B_f$$

$$R_f^{(1)}$$

$$R_f^{(2)}$$

$$B_f^{(1)}$$

$$B_f^{(2)}$$

$$\frac{R_f}{B_f}$$

$${\xi }_0^{(0)}$$

$${\xi }_2^{(0)}$$

Full text

1 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 $\sim$ 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\leq m\leq 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\geq 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$\mu$M in width [17, 16], 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.

2 Materials and Methods

2.1 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 $\tilde{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_{+}^{\prime}$ and $k_{-}^{\prime}$ 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). $\nu_{r}$ and $\nu_{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: $K_{r}=k_{-}/k_{+},K_{b}=k_{-}^{\prime}/k_{+}^{\prime},K_{L}=k_{d}/k_{a}$.

If $\xi_{m}(t)$ denote the fraction of receptors in methylation state $m$ at time $t$, the total active fraction of receptors is given by $\xi_{a}(t)=\sum_{m=0}^{M}a_{m}(L)\xi_{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$ [20, 14, 15]. To satisfy the earlier assumption, we choose $K_{0}=0$ and $K_{M}=\infty$ (these are ideal values). We further choose $a_{m}^{0}=m/M$[19], a simple mathematical form which satisfies our earlier assumptions.

2.2 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 $\nu_{r}$ and $\nu_{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 $\xi_{m}$ with $m=0,1,2$; hence $\sum_{m=0}^{2}\xi_{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}$ and ${\tilde{x}}_{m}$. For $m=1$, there is a further sub-division: $x_{1}=x_{1}^{a}+x_{1}^{i}$ and ${\tilde{x}}_{1}={\tilde{x}}_{1}^{a}+{\tilde{x}}_{1}^{i}$, as explained in the previous section. For each $m$, $\xi_{m}=x_{m}+{\tilde{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}\to\overline{x_{i}}(\boldsymbol{\xi})$ naturally appears in the relevant inter-modular conversion rates. Here, we have introduced the compact vector symbol $\boldsymbol{\xi}\equiv(\xi_{0},\xi_{2})$ (because of the normalization constraint mentioned above, only two of the fractions $\xi_{m}$ are independent, we shall take these to be $\xi_{0}$ and $\xi_{2}$ purely for reasons of convenience).

Therefore, within this approximation, the inter-modular dynamics may be expressed in terms of a probability density $P(\boldsymbol{\xi};t)$, which satisfies a non-linear Fokker-Planck equation (FPE) in the form of a local equation of continuity[2]:

[TABLE]

where $\textbf{J}\equiv\{J_{0},J_{2}\}$ with $J_{m}=v_{m}P-\frac{\partial}{\partial\xi_{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

[TABLE]

with $\delta=N^{-1}$, $\omega_{10}=N\nu_{b}\tilde{x}_{1}^{a}$, $\omega_{01}=N\nu_{r}\tilde{x}_{0}$, $\omega_{12}=N\nu_{r}\tilde{x}_{1}^{i}$ and $\omega_{21}=N\nu_{b}\tilde{x}_{2}$.

Incorporating the rates $\omega_{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 6.1, we arrive at the following expressions for the drift and diffusion coefficients:

[TABLE]

and

[TABLE]

The free enzyme concentrations $R_{f}$ and $B_{f}$ are also functions of $\xi_{0},\xi_{2}$, and the expressions are given by (see 6.1)

[TABLE]

3 Results

3.1 Averages as fixed points

From (1), it is seen after the required averaging that the time evolution of the mean receptor fractions is given by the vector equation

[TABLE]

where ${\bf v}(\boldsymbol{\xi})\equiv(v_{0}(\boldsymbol{\xi}),v_{2}(\boldsymbol{\xi}))$ is the drift “vector”. In general, this vector vanishes at one or more points in the $(\xi_{0},\xi_{2})$ space; following van Kampen[22], we shall refer to this point (assuming only one exists) as the fixed point, denoted $\boldsymbol{\xi}^{*}=(\xi_{0}^{*},\xi_{2}^{*})$. It can be shown (see later) that in the large $N$-limit, the fixed point becomes identical to the mean value $\overline{{\bf\xi}}$.

Under conditions $R_{f}\ll K_{r}$ and $B_{f}\ll K_{b}$, which we shall assume to hold for reasons of simplicity, it follows from (2.2) that the fixed point is given by

[TABLE]

where $R_{f}^{*}\equiv R_{f}(\boldsymbol{\xi}^{*})$ and $B_{f}^{*}\equiv B_{f}(\boldsymbol{\xi}^{*})$ are given by (2.2), with the replacement $\boldsymbol{\xi}\to\boldsymbol{\xi}^{*}$ in the right hand side of the equations. (2.2) and (3.1) implicitly give the fixed point $\boldsymbol{\xi}^{*}$.

3.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 $\boldsymbol{\xi}^{\prime}=\boldsymbol{\xi}-\boldsymbol{\xi}^{*}$; we also define its probability distribution $\Phi(\boldsymbol{\xi}^{\prime})\equiv P(\boldsymbol{\xi}^{*}+\boldsymbol{\xi}^{\prime})$, which satisfies the equation

[TABLE]

where $\bf{J}^{\prime}(\boldsymbol{\xi}^{\prime})\equiv\bf{J}(\boldsymbol{\xi}^{*}+\boldsymbol{\xi}^{\prime})$, and may be expanded as follows:

[TABLE]

Keeping the leading term in each, and noting that ${\bf v}(\boldsymbol{\xi}^{*})=0$, we arrive at (for $m$=0,2)

[TABLE]

where \beta_{mn}=\frac{\partial v_{m}}{\partial\xi_{n}}\Big{|}_{\boldsymbol{\xi}^{*}}, $D_{m}^{*}\equiv D_{m}(\boldsymbol{\xi}^{*})$. Substitution of (10) in (8) leads to the multivariate linear Fokker-Planck equation (LFPE) [22]

[TABLE]

Let us now define the covariances $\sigma_{mn}=\langle\xi^{\prime}_{m}\xi^{\prime}_{n}\rangle$ 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)

[TABLE]

where $\boldsymbol{\mu}=(\mu_{0},\mu_{2})$.The moment generating function admits the Taylor expansion

[TABLE]

Since $\Phi(\boldsymbol{\xi}^{\prime})$ satisfies the LFPE given by (11), the moment generating function satisfies the equation

[TABLE]

Under stationary conditions ($t\to\infty$), the left hand side of (14) becomes zero; next we substitute (13) in (14), and find that $\overline{\boldsymbol{\xi}^{\prime}}=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]

[TABLE]

where $D_{mn}=D^{*}_{m}\delta_{mn}$ are the elements of the (diagonal) diffusion matrix $\boldsymbol{D}$. The following expressions follow from (15):

[TABLE]

where ${\mathcal{D}}=(\beta_{00}+\beta_{22})(\beta_{00}\beta_{22}-\beta_{02}\beta_{20})$. Explicit expressions for the coefficients $\beta_{mn}$ are to be found in 6.2.

3.3 Mean and fluctuations in Activity

The active fraction of receptors in the present model is given by

[TABLE]

where $\ell=L/K_{L}$. Within LNA, the average fractions $\overline{\xi}_{m}=\xi_{m}^{*}$, hence $\overline{\xi}_{a}=\xi_{a}^{*}$, where $\xi_{a}^{*}$ is given by (17), with the replacements $\xi_{m}\to\xi_{m}^{*}$ in the right hand side and $\xi_{m}^{*}$ given by (3.1).

The variance of the active fraction is given by $\sigma_{a}^{2}=\overline{\xi_{a}^{2}}-\overline{\xi}_{a}^{2}$, which is also computed using (17). In terms of the covariances $\sigma_{mn}$, this is given by

[TABLE]

3.4 Ultrasensitivity in the mean activity

We will now explore the large $A_{0}$ limit of the fixed point. The fixed point values $\xi_{0}^{*}$ and $\xi_{2}^{*}$ given by (3.1) can be expressed alternatively as

[TABLE]

In order to understand the behavior of the above expressions in the $A_{0}\to\infty$ limit, let us conjecture expansions of the form

[TABLE]

Consider now the expressions for $R_{f}$ and $B_{f}$ as given in (2.2); after using the expansions given by (20), we find the following asymptotic forms in the large $A_{0}$ limit:

[TABLE]

where

[TABLE]

The ratio $R_{f}/B_{f}$, upto ${\mathcal{O}}(1/A_{0})$, is given by

[TABLE]

Substituting (23) in (3.4), we find

[TABLE]

which give the leading terms in the expansions in (20), and are valid for arbitrary $\ell$. In order to take the analysis further, we study the limits $\ell\to 0$ and $\ell\to\infty$ separately.

3.4.1 Small $\ell$ expansion

Consider now the limit of small ligand concentrations, $\ell\ll 1$. We assume that, in this case, the $zero^{th}$ order terms $\xi_{m}^{(0)}$ in (3.4) can be expanded as

[TABLE]

From (3.4), we have

[TABLE]

Substituting the small $\ell$ expansions of $\xi_{m}^{(0)}$ as given in (25), we obtain

[TABLE]

Substituting (27) in (3.4) leads to the equations

[TABLE]

where

[TABLE]

is the control parameter that characterizes ZOU. The implications are (a) the population in the highest methylation state is vanishingly small in the limits $\ell\to 0$ and $A_{0}\to\infty$ (b) $\xi_{0}^{00}=0$ or 1 unless $\alpha=1$: this points to a jump-like behavior for $\xi_{0}$ (and therefore, $\xi_{1}$, considering that $\xi_{0}^{00}+\xi_{1}^{00}=1$) in these limits. The ’critical point’ $\alpha_{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].

3.4.2 Large $\ell$ expansion

We next consider a large $\ell$ expansion of the $zero^{th}$ order term $\xi_{m}^{(0)}$ of large $A_{0}$ expansion in the following form

[TABLE]

Substituting these in (26), we obtain

[TABLE]

Therefore, after substitution of (31) in (3.4), we find

[TABLE]

which mirror (28) derived in the opposite, $\ell\to 0$ limit. In the present case, it is $\xi_{2}$ that displays the jump transition from 1 to 0 as $\alpha$ crosses 1 (and $\xi_{1}$ does the reverse), while $\xi_{0}$ is vanishingly small at all $\alpha$. These analytical results, given in equations (28) and (32) are supported by numerical simulations, to be discussed in the following section.

4 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 $\xi_{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.

4.1 Simulation results for $M=2$

For fixed $A_{0}$, the mean active fraction of receptors $\xi_{a}$ undergoes a sharp rise as $R_{0}$ is increased, as shown in Fig (2). $\xi_{a}$ is independent of $\ell$, as is clear from the figure, and agrees with the analytical prediction (fixed point). The inset shows the variance in activity, which depends on $\ell$ 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\mu$M (Fig 2a), the rise is less steeper compared to $A_{0}=13.6\mu$M (Fig 2b); for same $\ell$, the variance is more in (a), but has a sharper peak in (b). For fixed $\ell$, and changing $A_{0}$, the mean shows the expected ultrasensitive rise as function of $R_{0}$ (Fig4). For the set of parameters used here in Table 1, the critical value for $R_{0}$ corresponding to $\alpha=1$ (29) is $R_{0}^{c}=0.224\mu$M, which agrees with our observations in Fig4. 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 $\sigma_{a}^{2}$ show a non-monotonic variation with $\ell$? To find out, we studied the mathematical expressions for the individual variances $\sigma_{00},\sigma_{22}$ and the covariance $\sigma_{02}$ as functions of $R_{0}$ and $\ell$, at fixed $A_{0}$. The plots, shown in Fig 5, indicate that the non-monotonicity originates from the covariance $\sigma_{02}$ between the fully inactive and fully active methylation levels. Both $\sigma_{00}$ (Fig 5a) and $\sigma_{22}$ (Fig 5b) have their peaks near $R_{c}$; but while the peak value of $\sigma_{00}$ is a decreasing function of $\ell$, that of $\sigma_{22}$ is an increasing function of $\ell$. $\sigma_{02}$ (Fig 5c and d) is negative throughout and has its minimum near $R_{c}$; however, in (c) the lowest value changes non-monotonically with $\ell$. From (17), it is seen that $\sigma_{a}^{2}$ is dominated by $\sigma_{00}$ for small $\ell$, while the large $\ell$ behavior is dominated by $\sigma_{22}$. In the intermediate $\ell$ regime, $\sigma_{02}$ also contributes significantly, and therefore the non-monotonic variation of its minimum value with $\ell$ affects the peak value of $\sigma_{a}^{2}$, which appears to be minimized near $\ell=2$.

Next, we try to understand how the receptors distribute themselves among different methylation levels and the roles of $\ell$ and $A_{0}$ in this distribution. For fixed $A_{0}$ and very small $\ell$, 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 $\ell$, the transition becomes less sharper/smoother leading to $\xi_{0}$ being almost zero in the entire range of $R_{0}$ for very large $\ell$, 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 $\ell$ and then increases in the upper-critical regime with increase in $\ell$, finally leading to an ultrasensitive switch from near-zero to near-unity for very large $\ell$ Fig (6b). But the fraction in the intermediate methylation level exhibits switch-like behaviors of opposite nature at very low and high $\ell$ values; while for intermediate $\ell$, it shows non-monotonic change with $R_{0}$ with a peak near $R_{c}$ Fig (6c and d). With increase in $A_{0}$, the transitions in $\overline{\xi}_{0}$, $\overline{\xi}_{1}$ and $\overline{\xi}_{2}$ become sharper (Fig 7), and agrees with the predictions made in (28) and (32). Therefore, in the $\ell\to 0$ as well as $\ell\to\infty$ 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.

4.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 $\ell$ 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 $\ell=20$ and $A_{0}=13.6\mu$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.

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-d, we show, for $M=3$, the four mean populations $\overline{\xi}_{0},\overline{\xi}_{1},\overline{\xi}_{2}$ and $\overline{\xi}_{3}$ respectively for a limited range of $\ell$ 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 $\ell$ (e.g. $\ell=20$), most of the receptors are in $m=0$ or $m=3$ states; in fact, $\overline{\xi}_{0}$ and $\overline{\xi}_{3}$ show near-ultrasensitive fall and rise respectively across $R_{c}$. For high $\ell$ (e.g. $\ell=20,000$), while both $\overline{\xi}_{0}$ and $\overline{\xi}_{1}$ become smaller as $A_{0}$ is increased, $\overline{\xi}_{2}$ and $\overline{\xi}_{3}$ become dominant; now, $\overline{\xi}_{2}$ falls abruptly as $R_{0}$ crosses $R_{c}$, while $\overline{\xi}_{3}$ 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 $\ell\to 0$ and $\ell\to\infty$.

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.

4.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 $\delta L(t)$ in the external ligand concentration is switched on in the extracellular environment of the bacterium at $t=0$, such that $L(t)\to L+\delta L(t)$ for $t\geq 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

[TABLE]

where $\chi_{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; $\delta L(t)=\delta L_{s}\Theta(t)$ where $\Theta(t)$ is the Heaviside theta-function and $\delta L_{s}$ is the size of the step. Let $\delta\overline{\xi_{a}}^{\rm(step)}(t)$ be the corresponding response. From (33), it follows that

[TABLE]

From (17), we find

[TABLE]

where $\delta\overline{\xi}_{0}(t)$ and $\delta\overline{\xi}_{2}(t)$ are the changes in $\overline{\xi}_{0}$ and $\overline{\xi}_{2}$, respectively, in response to a change $\delta L(t)$ in the ligand concentration, and are evaluated using (6). The resulting equations have the form

[TABLE]

where the coefficients $\beta_{mn}$ have been defined earlier, following (10), while $\gamma_{m}=\partial v_{m}/\partial L|_{\boldsymbol{\xi}^{*}}$. Explicit expressions for both $\beta_{mn}$ and $\gamma_{m}$ are given in 6.2.

It is convenient to subject (35) and (36) to Laplace transforms, after putting $\delta L(t)=\delta L_{s}\Theta(t)$, and use (34) to compute the response function. For the Laplace transform of $\chi_{a}(t)$, we find

[TABLE]

where $\tilde{\chi}_{m}(s)$ are Laplace transforms of $\chi_{m}(t)$ ($m=0,2$) which are defined through the relations $\delta\overline{\xi}_{m}(t)=\int_{0}^{t}\chi_{m}(t-t^{\prime})\delta L(t^{\prime})dt^{\prime}$. The explicit expressions turn out to be

[TABLE]

where $d=(s-\beta_{00})(s-\beta_{22})-\beta_{02}\beta_{20}$. We have confirmed that the response function curve computed in (37) encloses zero area, i.e., $\tilde{\chi}_{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 $\overline{\xi}_{a}$ on $\ell$ in steady state, recall Fig (2)). In 6.3, we use the results in (37) and (4.3) 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.

5 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\to\infty$, independent of the attractant concentration $L$. At the same time, its variance $\sigma_{a}^{2}$ 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 $\sigma_{a}\simeq 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 6.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\mu M$) even for a given concentration of the inducer isopropyl $\beta$-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\sim 10.3\pm 1.1$, with the sensitive part being less than $3\mu$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$\mu$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 $\sim 16.5\pm 1.1$ and $20.7\pm 1.6$ for CW biases of $0.8\pm 0.1$ and $0.5\pm 0.1$ respectively. Because of the higher Hill coefficient, the sensitive window of [CheYp] is reduced to about 1$\mu$M in this case, with center around 3$\mu$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 6.4 we estimate that the maximum standard deviation in clockwise bias (arising from spontaneous fluctuations in [CheYp]) $\delta 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 $\sim$ 0.5 and the upper value for mean CW bias $\sim$ 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 6.5 shows that this effect is negligible when [CheYp] lies in the range $2.5\mu M-4\mu M$, which more than covers the sensitive window. Outside this range, motor-level adaptation increases the standard deviation in CW bias, (see figure in 6.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 6.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.

6 Appendices

6.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 ($\nu_{r}$ and $\nu_{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

[TABLE]

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 $x_{0}+\tilde{x}_{0}=\xi_{0}$, $x_{1}^{i}+\tilde{x}_{1}^{i}=\xi_{1}^{i}$, $x_{1}^{a}+\tilde{x}_{1}^{a}=\xi_{1}^{a}$ and $x_{2}+\tilde{x}_{2}=\xi_{2}$, we find

[TABLE]

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 $\boldsymbol{\xi}^{*}$:

[TABLE]

The free enzyme concentrations $R_{f}$ and $B_{f}$ are determined by the normalization conditions

[TABLE]

Using (6.1) in (6.1) leads to (2.2).

Perfect adaptation (insensitivity of average activity to L):

Although the fixed point $\boldsymbol{\xi}^{*}$ depends on $L$, the mean total active fraction does not; this can be easily shown as follows. The equations (41) for intermodule dynamics, after summing the left hand side give

[TABLE]

After substituting the equations for intra-module dynamics (6.1) (in the limit $R_{f}\ll K_{r}$, $B_{f}\ll K_{b}$) and the expression for $\xi_{a}$ in terms of $\xi_{m}$ (17) in (43), the following quadratic equation for $\xi_{a}^{*}$, independent of $L$, is obtained:

[TABLE]

The solution to the above equation is

[TABLE]

6.2 Evaluation of various parameters

The expressions of $\beta_{mn}$ can be obtained by differentiating (2.2). Note that the expressions for $v_{m}$ contain $R_{f}$ and $B_{f}$ which are in turn functions of $\xi_{m}$. Therefore,

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

To evaluate the linear response function, we have to consider the derivative $\gamma_{m}$, the explicit expression of which is

[TABLE]

The relevant quantities involving derivatives with respect to $L$ are listed below:

[TABLE]

[TABLE]

In the computation of the coefficients $\beta_{mn}$ and $\gamma_{m}$, all the above mentioned quantities are to be evaluated at the fixed point $(\boldsymbol{\xi}^{*},R_{f}(\boldsymbol{\xi}^{*}),B_{f}(\boldsymbol{\xi}^{*})$ and it has been assumed that $R_{f}\ll K_{r}$ and $B_{f}\ll K_{b}$ .

6.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$\to$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

[TABLE]

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 $\lambda_{Y}$ is the rate of degradation of CheYp (due to the dephosphatase CheZ). The steady state mean concentration of CheYp can be obtained from (53) by equating left hand side to [math],

[TABLE]

In response to a change in attractant concentration $L\to L+\delta L(t)$ for $t\geq 0$, the mean [CheYp] undergoes a shift $Y^{*}\to Y^{*}+\delta Y(t)$, with $\delta Y(t)$ given by

[TABLE]

The clockwise bias, i.e., the probability $P_{\rm{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$[16, 13, 17]:

[TABLE]

where $H$ is the Hill coefficient symbolizing the steepness of the transition. It follows that fractional changes in the CW bias $P_{\rm{CW}}$ and $Y$ are related as

[TABLE]

where, in the second part, we have defined the response function $\chi_{b}(t)$ for the bias. Using (55) and (54) in (57), it follows that, in steady state, $\chi_{b}(t)$ is related to $\chi_{a}(t)$ through

[TABLE]

From (58), it follows that the Laplace transforms of $\chi_{b}(t)$ and $\chi_{a}(t)$ are related linearly: $\tilde{\chi_{b}}(s)\propto\tilde{\chi_{a}}(s)/(s+\lambda_{Y})$, hence $\tilde{\chi_{b}}(0)=0$, i.e., the response curve for the bias too encloses zero area[9, 15]. After using (37) and (4.3) in (58), and performing the inversion, we arrive at the following multi-exponential form for $\chi_{b}(t)$:

[TABLE]

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 $\beta_{mn}$ and $\gamma_{m}$ as follows:

[TABLE]

[TABLE]

In Fig.10, we show plots of the mathematical expression for $\chi_{b}(t)$ obtained in (59), (a) by varying $B_{0}$ and (b) varying $\ell$. 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 $\ell$, the depth of the negative lobe decreases wih no change in the time scales. For $\ell\gg 1$, the negative lobe effectively disappears.

6.4 On the role of biochemical noise in the run and tumble motion of E.coli

From (56), it follows that a fluctuation in [CheYp] equal to its standard deviation in magnitude would cause the following change in $P_{\rm{CW}}$:

[TABLE]

where $y=Y/K_{Y}$ and $\overline{y}$ is its mean value in steady state. The function $\overline{y}^{H-1}/(1+\overline{y}^{H})^{2}$ has a single maximum (and no minimum) at $\overline{y}=[(H-1)/(H+1)]^{1/H}\equiv\tilde{y}$, with maximum value 1/4 for $H\gg 1$. Using this property in (62) leads to the following upper bound on $\delta P_{\rm{CW}}$:

[TABLE]

The standard deviation $\sigma_{Y}$ has been measured in experiments, and the range of variation has been estimated[28] as

[TABLE]

with the actual value depending on the mean clockwise bias, ${\overline{P}}_{\rm CW}$. The dissociation constant $K_{Y}\simeq 3\mu M$ and the Hill coefficient $H\simeq 20$ (for pre-adapted motors) have also been measured in single cell experiments [16]. Using this estimate in (63), we find $\delta P_{\rm CW}\leq 0.15-0.367$ (assuming $\overline{y}\approx\tilde{y}$). So there is a maximum fluctuation of $18\%$ for $H=20$ corresponding to a CW bias of $0.5$ and [CheYp] $\sim 3\mu 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 (53). Under steady state conditions, the variance $\sigma_{Y}^{2}=\langle Y^{2}\rangle-\langle Y\rangle^{2}$ is given by the expression

[TABLE]

where, we have introduced the correlation length $\lambda_{\xi}$ for the fluctuations in $\xi_{a}$, defined through the relation $\langle\delta\xi_{a}(t)\delta\xi_{a}(t^{\prime})\rangle=\sigma_{a}^{2}\exp(-\lambda_{\xi}|t-t^{\prime}|)$ in steady state, where $\delta\xi_{a}(t)=\xi_{a}(t)-\overline{\xi_{a}}$ is the fluctuation in activity relative to its steady state value. From (65), the inequality

[TABLE]

follows, the equality being satisfied in the limit $\lambda_{\xi}\ll\lambda_{Y}$ (which we shall assume henceforth). Various estimates for the parameters in (66) are available in the literature (see Table2). Using (64), we identify the following ranges for the standard deviation $\sigma_{a}$ of receptor activity:

[TABLE]

It is also of interest to compare the values of the ZOU parameter $\alpha$ across these sets of parameters, which turn out to be 0.547 (i), 0.076 (ii) and 0.035 (iii) (Table2). Note that the estimated variance in (ii) and (iii) are smaller than (i), which is consistent with our expectation of maximum variance near $\alpha=1$, over a large range of values of $\ell$. The magnitude of receptor noise predicted by (6.4) also agrees with our predictions (see Fig (3), also Fig (8)).

6.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 $\delta Y(t)$ henceforth) be written in the general linear form

[TABLE]

where $\chi_{m}(\tau)$ (defined for $\tau>0$) is the linear response function for the motor, which we assume to consist of two terms:

[TABLE]

where $\alpha_{\rm na}=\partial_{Y}P_{\rm CW}^{\rm na}$, and $P_{\rm CW}^{\rm na}$ is the bias of the non-adapted motor, as measured by Yuan and Berg (2013), with Hill coefficient $H\approx 20$. The second term represents the slow change in the bias originating from motor-level adaptation, with $\lambda_{m}^{-1}$ being the time scale for the same. Now, for a step-like change in $Y$ such that $Y\to Y_{0}+\Delta Y$ at $t=0$, let us put $\delta Y(t)=\Delta Y\Theta(t)$; using this in (68) yields the response for the same:

[TABLE]

As $t\to\infty$, the system will adapt, and hence $\delta P_{\rm CW}^{\rm step}(t)\to\alpha_{\rm a}\Delta Y$, where $\alpha_{\rm a}=\partial_{Y}P_{\rm CW}^{\rm a}$, and $P_{\rm CW}^{\rm a}$ is the clockwise bias of the adapted motor. From (69), it then follows that $\Gamma=\lambda(\alpha_{\rm a}-\alpha_{\rm 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 $\delta P_{\rm CW}\equiv\sqrt{\langle\delta P_{\rm CW}(t)\rangle^{2}}$, with the average to be carried out over different realizations of the fluctuations in [CheYp]. From (68), it turns out that

[TABLE]

where, in steady state, we expect $\langle\delta Y(t_{1})\delta Y(t_{2})\rangle=\sigma_{Y}^{2}e^{-\lambda_{Y}(t_{1}-t_{2})}$ for $t_{1}>t_{2}$, $\sigma^{2}_{Y}$ is the steady state variance of $Y$ as given by (65) and $\lambda_{Y}$ is the dephosphorylation rate of [CheYp]. Using the expression for $\chi_{m}(t)$ from (69) and completing the calculations yields the final result

[TABLE]

The term outside the square root is precisely the result in (62) in 6.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.

Acknowledgments

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