Resolving the binding-kinase discrepancy in bacterial chemotaxis: A nonequilibrium allosteric model and the role of energy dissipation

To start, we focus on the in vitro measurements by Amin and Hazelbauer on Tar receptors embedded in native membrane vesicles [27]. The receptors were fixed at different levels of methylation by substituting glutamates (E) with glutamine (Q) at the receptor methyl-accepting sites. By comparing dose-response curves of receptors with zero or three modifications, it was found that methylation shifts the kinase response curve significantly but only induces a much smaller shift in ligand binding (see Fig. 1B).

Before introducing the nonequilibrium allosteric model, we first examine whether the behaviors (especially ligand binding and kinase activity) of the chemoreceptors can be consistently described by an equilibrium allosteric model. To this end, we present a parametric test that systematically detects inconsistency between the measurements and the MWC model without fitting the data to any specific function. Other equilibrium models, such as Ising models, can be ruled out using a similar approach (see SI Appendix).

In the MWC model, the receptor complex has $N$ identical binding sites, whose occupancy is denoted by $\sigma_{i}$ $(i=1,2\dots N)$; $\sigma_{i}=0,1$ represents a vacant or occupied receptor, respectively. The receptor methylation level is denoted by $m$. The MWC model assumes all-or-none behavior for the whole complex. Namely, the receptor activity $s$ is either on ($s=1$) or off ($s=0$). The energy (Hamiltonian) of the receptor complex is given by

with $\mu=\log\quantity([L]/K_{i})$ being the chemical potential for binding, which depends on ligand concentration $[L]$ and a (methylation-dependent) dissociation constant $K_{i}(m)$. $E_{s}$ is the energy difference between the active and inactive receptor states in the absence of a ligand. Each occupied binding site increases this energy difference by $E_{0}>0$, thereby suppressing activity at high occupancy.

Given the MWC Hamiltonian, the average receptor activity is

and the average binding is

In the limit of large $E_{s}$, we have $\expectationvalue{s}\ll 1$, so the binding curve becomes a Hill function with Hill coefficient $n_{H}=1$ (no cooperative binding), consistent with experiments using vesicle-bound chemoreceptors [27]. The maximum activity $s_{\max}=\quantity(1+\exp(E_{s}))^{-1}$ is reached at $[L]=0$, and the minimum activity $s_{\min}=\quantity(1+\exp(E_{s}+NE_{0}))^{-1}$ occurs at $[L]=\infty$.

The MWC model has seen great success in capturing the methylation-dependence of the downstream CheY-P response [7, 25], even in the presence of time-varying stimuli [31]. These studies suggest the following basic mechanism: methylation affects the energy for the active state, $E_{s}(m)$, and thereby shifts the kinase response curves across orders of magnitude in ligand concentration. However, this mechanism does not take into account measurements of ligand binding.

Here, we examine whether the binding and activity curves measured at the same methylation level are both compatible with the equilibrium MWC model. This can be done by eliminating the concentration variable $[L]$ from Eqs. (2) and (3) to obtain a parametric relation between $\expectationvalue{s}$ and $\expectationvalue{\sigma}$. For example, solving Eq. (2) for $[L]$ and substituting the solution into Eq. (3) gives the mean occupancy as a function of the receptor activity,

where $s_{\max}$ and $s_{\min}$ are the maximum and minimum activity mentioned above, which can be determined from experimental measurements of kinase activity. For any fixed $N$, $s_{\max}$ and $s_{\min}$ uniquely determine the value of the energy parameters $E_{0}$ and $E_{s}$, which gives $\tilde{\sigma}$ as a function of $(\expectationvalue{s},s_{\max},s_{\min})$. The expression of $\tilde{\sigma}$ for finite $N$ can be complicated, but its behavior can be illustrated in small and large $N$ limits. For $N=1$, the equilibrium model simplifies to $\tilde{\sigma}=(s_{\max}-\langle s\rangle)/(s_{\max}-s_{\min})$, which equates ligand binding with the normalized kinase activity. On the other hand, for large $N$, the inferred receptor occupancy converges to

The curves for intermediate $N$ lie between these two extreme limits.

To test the validity of the equilibrium MWC model, we transform the measured kinase response curves for two different methylation levels [27] according to Eq. (4) and plot the inferred occupancy ($\tilde{\sigma}$) against the corresponding measurements of the receptor occupancy ($\sigma_{\mathrm{exp}}$) ¹¹1Following previous studies, we assume the measured CheY-P concentration is proportional to activity in the model: $[\text{CheY-P}]=A_{0}\langle s\rangle$, where $A_{0}$ is the saturating concentration of CheY-P. We use $A_{0}=60$pM chosen to be slightly above the maximum CheY-P concentration measured in this set of experiments (the results are robust to increasing the value of $A_{0}$).. If the equilibrium model is valid, all the points should collapse onto the diagonal $\sigma=\tilde{\sigma}$ for some choice of $N$. However, the data for both methylation levels ($m=0,3$) lie well off the diagonal for any choice of $N$ (Fig. 1A), indicating that the system (even for a single methylation level) cannot be described by the equilibrium MWC model with any $N$.

The intuition behind this inconsistency is that the MWC model cannot capture the relative shift between binding and response curves. To demonstrate this, we fit the MWC model only to the measured kinase activity (CheY-P level) and compare the resulting binding and activity curves with experiments (Fig. 1B). The MWC model successfully captures kinase activity (upper panel) but does not produce the correct binding curves (lower panel). Instead, it predicts that the increase in binding and the decrease in kinase activity should occur at around the same ligand concentration. This prediction is inconsistent with experiments, which found the sharpest change in kinase activity occurring at a concentration that is either much lower ($m=0$) or higher ($m=3$) than that of binding. Moreover, as the methylation level changes from $m=0$ to $m=3$, the shift in activity curves is much greater than the shift in binding curves. These results suggest that the MWC model is unable to capture the data either at a single methylation level or the change between different methylation levels. Similarly, fitting the MWC model to the measured binding curves results in discrepancies with the measured kinase activity (see SI Appendix). Note that these fittings are only shown to build intuition. They are all encompassed by the parametric test (Fig. 1A), which offers stronger evidence by revealing inconsistency with the MWC model without relying on fitting with any assumptions or any particular formulation of the loss function.

Similar results hold for more complex equilibrium models: for example, if we treat the measured kinase activity $a$ as a separate degree of freedom from the receptor activity $s$. Assuming the two activities are coupled by equilibrium mechanisms, with the active receptor promoting kinase phosphorylation, we add the following terms to the MWC Hamiltonian,

where $F_{0}$ is the energy of the active kinase when the receptor is inactive and $F_{1}=F_{0}+\Delta F<F_{0}$ is the energy when the receptor is active. The average activities of the joint Hamiltonian $H+H_{a}$ are linearly related,

Therefore, our analysis above also applies to equilibrium models with additional binary degrees of freedom. In the SI Appendix, we show that this linear relation holds for arbitrary chains of binary variables coupled via equilibrium interactions. Hence, these types of models are also inconsistent with the experiments.

Beyond MWC models, the binding and kinase response cannot be simultaneously captured by equilibrium Ising-type models (for various spatial structures), where the receptor activities $s_{i}$ are variable across the receptor complex (see SI Appendix).

The fact that equilibrium models cannot capture the binding and kinase response measurements indicates that the E. coli chemotaxis signaling pathway operates out of equilibrium. We now present a nonequilibrium allosteric model that extends the MWC model and can simultaneously capture the measured binding and response curves as well as their asymmetric shifts due to methylation.

In this model (pictured in Fig. 2A), the receptor conformation is again represented as a binary variable $s$ with inactive (0) and active (1) states. It is coupled to ligand occupancy $\sigma_{i}$ by an equilibrium MWC model. The phosphorylation state of the response regulator CheY is represented by a new binary variable $a$, with dephosphorylated (0) and phosphorylated (1) states. The measured CheY-P concentration is proportional to the average phosphorylation variable $\expectationvalue{a}$. The transitions between states $a=0,1$ involve phosphorylation-dephosphorylation (PdP) cycles (blue boxes), whose reaction rates depend on both the conformation $s$ and the methylation $m$. The system operates out of equilibrium sustained by continuous ATP hydrolysis in the PdP cycle, and the average receptor behavior is determined by the probability distribution $P(a,s)$ of the nonequilibrium steady state.

When the receptor is active ($s=1$), it catalyzes the phosphorylation of CheY through the autophosphorylation of the histidine kinase CheA, which subsequently transfers the phosphate group to CheY. In our model, we do not consider all the detailed steps of phosphotransfer from ATP to CheY and assign an overall rate $k_{1}$ to describe the phosphorylation of CheY. For simplicity, we also assume that the phosphorylation rate is negligible when the receptor is inactive ($s=0$). The dephosphorylation of CheY-P is catalyzed by CheZ, with a dephosphorylation rate $k_{z}$ independent of the receptor state. To ensure thermodynamic consistency, all reactions are reversible with reverse rates $k_{1}^{\prime}=k_{1}e^{-G_{1}}$ and $k_{z}^{\prime}=k_{z}e^{-G_{z}}$. The PdP cycle is driven by ATP hydrolysis with free energy $G_{\mathrm{ATP}}$ (measured in units of $k_{B}T$), which is dissipated per PdP cycle per CheY molecule. The rates obey the thermodynamic constraint:

where $G_{1}$ and $G_{z}$ are the free energy released by phosphorylation and dephosphorylation respectively. Note that $\frac{k_{1}k_{z}}{k_{1}^{\prime}k_{z}^{\prime}}=e^{G_{\mathrm{ATP}}}\gg 1$, which clearly indicates a violation of detailed balance and the far-from-equilibrium nature of our model.

Motivated by structural details of the bacterial chemoreceptors (see Discussion), we assume that receptor methylation can affect kinase activity by changing the energy barrier in the phosphorylation reaction, which scales $k_{1}$ and $k_{1}^{\prime}$ by the same factor. On the other hand, the dephosphorylation rate $k_{z}$ remains independent of $m$. The ratio of the two rates can be written as: $k_{1}/k_{z}=\exp\quantity[E_{p}-\Delta E_{B}(m)]$, where $E_{p}$ is a constant and $\Delta E_{B}(m)$ captures the allosteric effect of receptor methylation on the energy barrier for phosphorylation. Here, we adopt the simplest model where the energy barrier has a linear dependence on the methylation level, i.e. $\Delta E_{B}(m)=-\epsilon_{m}Nm$. The energy difference $G_{1}=\log(k_{1}/k_{1}^{\prime})$ is not affected by methylation.

The receptor conformation $s$ can switch between ON and OFF states at rates $\omega$ and $\omega^{\prime}$, whose ratio is defined as $\alpha=\omega/\omega^{\prime}$. Throughout our analysis, switching is assumed to be much faster than the PdP reactions ($\omega,\omega^{\prime}\gg k_{1},k_{z}$), which means $s$ is in fast equilibrium with $\sigma$. Since the coupling between $s$ and $\sigma$ does not involve any nonequilibrium driving, we can describe it with the MWC free energy (Eq. 1). Matching the switching dynamics to the MWC model gives the rate ratio:

The measured CheY-P concentration is proportional to the steady-state average $\expectationvalue{a}=(1+\mathcal{P})^{-1}$, where $\mathcal{P}$ is the probability ratio of being not phosphorylated versus phosphorylated:

In the large dissipation (strong nonequilibrium driving) limit, the reverse reaction rates are negligible, i.e. $k_{i}^{\prime}/k_{i}\to 0$. With an additional assumption $e^{E_{s}}\gg 1$, the activity reduces to,

This expression has the same form as the MWC activity Eq. (2) with a new effective energy for the active state $E_{s}^{\text{eff}}=E_{s}-E_{p}-\epsilon Nm$, which depends on the methylation level $m$. Given that the MWC model [Eq. (2)] has been successfully employed to fit the CheY-P curves (Fig. 1B), we can expect a similar, if not better, fit from the nonequilibrium model. Indeed, Fig. 2B shows that the model successfully captures the CheY-P level of vesicle-bound Tar receptors at all five methylation levels [27].

Since there is no feedback from the CheY phosphorylation state to the receptor, the average receptor conformation $\expectationvalue{s}$ and ligand occupancy $\expectationvalue{\sigma}$ are given by the MWC model, Eqs. (2) and (3) respectively. For large $E_{s}$, we again have $\langle s\rangle\ll 1$ so that the binding curve closely resembles a Hill function with Hill coefficient $n_{H}=1$. Fig. 2B shows that the same set of parameters used to fit kinase response also produces the correct binding curves, which is a significant improvement from the MWC model.

Why is the nonequilibrium model able to capture the asymmetric shift of binding and kinase response curves due to methylation while equilibrium models cannot? In equilibrium models, methylation can only affect receptor behavior through thermodynamic control, i.e. by changing the energy difference between different states through $E_{s}$ and $E_{0}$. Since equilibrium interactions are always reciprocal (symmetric), the fact that ligand binding controls receptor activity means that there has to be an equally strong feedback from the receptor state to the ligand occupancy. Therefore, this type of control shifts the binding and kinase response curves by similar amounts, inconsistent with experimental observations. In the nonequilibrium model, however, the interaction between ligand binding and kinase activity can be non-reciprocal (asymmetric): the receptor complex acts as an enzyme that exerts kinetic control by changing the phosphorylation energy barrier $\Delta E_{B}(m)$. In the absence of strong feedback from the substrate (CheY) to the receptor complex, these changes in the energy barrier enable amplified shifts in the kinase response while maintaining modest shifts in binding response. For this mechanism to function, however, the system must be driven out of equilibrium (in this case by continuous ATP hydrolysis). Energy dissipation is required to enable kinetic control, which has no impact on the steady-state distribution if the system is in equilibrium, and is necessary to maintain the asymmetric coupling between the phosphorylation state $a$ and the receptor state $s$. As shown later, the response amplitude vanishes in the absence of energy dissipation.

In addition to explaining the asymmetric shifts in binding ($\sigma$) and kinase response ($a$) curves, our model also predicts different dose-response curves for the receptor conformation ($s$) and kinase response ($a$). Vaknin and Berg measured Tar receptor conformation in vivo by fusing fluorescent proteins to the C-termini of chemoreceptors and measuring the intensity of fluorescence resonance energy transfer (FRET) between receptor homodimers [29]. They measured the receptor conformation for mutants with fixed methylation states EEEE ($m=0$) and QQQQ ($m=4$) and the corresponding kinase response for QQQQ ($m=4$) and $cheR^{+}cheB^{+}$ mutants. The $cheR^{+}cheB^{+}$ cells do not have a fixed methylation state, but their kinase response is known to be similar to QEEE mutants ($m=1$) [31]. The shift of the kinase response curve from $m=1$ to $m=4$ was found to be much more significant than the shift of the receptor conformation curve from $m=0$ to $m=4$ (Fig. 3A, circles versus triangles). Note that the conformation measurements were carried out in the absence of kinase CheA and linker protein CheW, which leads to noncooperative responses ($N=1$). Similar to the in vitro case discussed above, the discrepancy between the dose-response curves suggests the necessity for a nonequilibrium allosteric model. Indeed, as shown in Fig. 3A, the model simultaneously fits both the receptor conformation and kinase response curves for all mutants (solid and dashed lines).

Another abundant chemoreceptor in E. coli is the serine receptor Tsr, which regulates CheY phosphorylation using the same microscopic mechanisms as the Tar receptor. Levit and Stock [28] measured the binding and kinase response for Tsr receptor complexes with three different receptor methylation levels (Fig. 3B). This was achieved by expressing the WT receptor without CheR or CheB to fix its methylation level and subsequently increasing the methylation level by adding CheR and S-adenosylmethionine (SAM) or decreasing the methylation level by adding CheB_c, the catalytic domain of CheB. As the methylation level varied between these three conditions, the serine concentration required to inhibit kinase activity varied by more than two orders of magnitude, while the binding affinity only changed by two-fold (from $K_{d}=10\mu$M to $K_{d}=20\mu$M). The asymmetric shift in binding and kinase response curves of Tsr receptors is similar to that found for Tar receptors, which is inconsistent with equilibrium models. Once again, the discrepancy between these shifts can be fully captured by the nonequilibrium allosteric model (Fig. 3B, solid and dashed lines).

Our results above show that strong nonequilibrium driving, which is provided by the free energy released during ATP hydrolysis ($G_{\mathrm{ATP}}$), is necessary to produce the experimentally observed large shifts in kinase response curves with changes in receptor methylation level (or other forms of modification to the receptor) while leaving ligand binding relatively unchanged. This raises the question of how much energy dissipation is required for the system to exhibit behaviors that match the experimentally observed ligand binding and kinase response.

To address this question, we evaluate the difference between the average phosphorylation level $\expectationvalue{a}$ (proportional to the CheY-P concentration) and its value in the infinite dissipation limit $\expectationvalue{a}_{\infty}$ [given by Eq. (11)] by expanding the full solution to the nonequilibrium model [Eq. (10)] in terms of the reverse reaction rates $k_{1}^{\prime}$ and $k_{z}^{\prime}$. The leading order correction $\delta\langle a\rangle\equiv\expectationvalue{a}-\expectationvalue{a}_{\infty}$ is

where $P_{\text{ON}}=\alpha/(1+\alpha)$ is the probability of the receptor being in the ON state ($s=1$). The two terms in the parenthesis have clear physical meanings: they quantify the irreversibility of the PdP cycle by comparing the reaction rates along or against the nonequilibrium driving (yellow arrow in Fig. 2A). The first term $k_{z}^{\prime}/k_{1}$ is the ratio between the reverse dephosphorylation rate $k_{z}^{\prime}$ (counterclockwise transition) and the phosphorylation rate (clockwise transition) originating from the CheY state ($a=0$). Similarly, the second term $k_{1}^{\prime}/k_{z}$ is the ratio between reverse phosphorylation (counterclockwise) and dephosphorylation (clockwise) reactions from the CheY-P state ($a=1$). Since phosphorylation only occurs in the ON state, while dephosphorylation occurs independent of the receptor state, each term is weighted by an appropriate factor of $P_{\text{ON}}$. To make each of these terms small requires large $G_{z}$ and $G_{1}$ (and hence large $G_{\mathrm{ATP}}$) to drive the PdP cycle sufficiently irreversibly so that the system behaves close to the infinite-dissipation limit.

To quantify the minimum required dissipation $G_{\mathrm{ATP}}$, we start with Eq. (12) and focus on the extreme cases $m=0$ and $m=4$, whose response curves envelop those of intermediate $m$. For $m=0$, the CheY-P response is suppressed: the phosphorylation rate $k_{1}=k_{z}\exp(E_{p})$ is much less than the dephosphorylation rate $k_{z}$ ($E_{p}\approx-4$ for our best fit in Fig. 2), so the first term in Eq. (12) dominates. We want the error to be small compared to the maximal infinite-dissipation activity, i.e., $|\delta\expectationvalue{a}/\max(\expectationvalue{a}_{\infty})|<d_{th}$ for some loss threshold $d_{th}$, where $\max(\expectationvalue{a}_{\infty})=(1+e^{Es}k_{z}/k_{1})^{-1}$. Since $E_{s}$ is large and therefore $\alpha$ is small, this leads to $k_{z}^{\prime}/(k_{1}\alpha)<d_{th}$, or

Conversely, for $m=4$, the CheY-P response is amplified by the kinetic barrier shift due to methylation, and we have $k_{1}=k_{z}\exp(E_{p}+4\epsilon_{m}N)\gg k_{z}$. Therefore, the second term in Eq. (12) dominates. Again requiring this term is small compared to the maximal infinite-dissipation activity leads to $k_{1}^{\prime}/k_{1}<d_{th}$, or

Adding these two inequalities gives an operation threshold for $G_{\mathrm{ATP}}$,

From our fits, $E_{p}-E_{s}\approx-4$, so using a threshold of $d_{th}=10^{-2}$ leads to $G_{\min}=13.2$, which is the minimum dissipation energy needed to drive the system to have the experimentally observed behaviors shown in Fig. 2B.

As given in Eq. (8), the total dissipation energy ($G_{\mathrm{ATP}}$) can be decomposed into two parts: $G_{\mathrm{ATP}}=G_{1}+G_{z}$ where $G_{1}$ and $G_{z}$ are used to suppress the reverse reactions for the phosphorylation and the dephosphorylation processes, respectively. In the infinite dissipation limit with $G_{1},G_{z}\rightarrow\infty$, all the reverse reactions can be neglected, and we obtain the exact expression for the kinase activity given in Eq. (11), which agrees quantitatively with the experiments. For finite values of $G_{1}$ and $G_{z}$, the behavior of the model can be different. The difference between the phosphorylation levels at finite and infinite dissipation can be quantified by

where the weight $w_{m}=1/\expectationvalue{a}_{\infty}(m,0)$ normalizes curves by the infinite-dissipation maximal activity at the corresponding methylation state, and the average is calculated by sampling uniformly over the activity range of $\expectationvalue{a}_{\infty}$. Fig. 4A shows $d(G_{1},G_{z})$ for different values of $G_{1}$ and $G_{z}$. As expected, the difference is small in the upper right region when both $G_{1}$ and $G_{z}$ are sufficiently large, which defines the operational regime of the chemoreceptor sensory system. The estimated thresholds Eqs. (13) and (14), shown by white dashed lines, accurately predict when the models have quantitatively similar kinase activity. Moreover, they are indeed limited by methylation levels $m=0$ and $m=4$ as shown by the contours of the output discrepancy for individual methylation levels (white solid lines). The blue dashed line shows the operational threshold $G_{\min}=13.2$, where the difference is only small for a particular pair of $(G_{1},G_{z})$. In contrast, for a typical physiological value $G_{\mathrm{ATP}}=20$ [33] (blue solid line), $(G_{1},G_{z})$ is allowed to vary within a much larger range without deviating from the experimentally observed kinase response.

Next, we investigate the role of energy dissipation in the signaling pathway. In particular, we find that it enhances the amplitude and the sensitivity range of the response. To illustrate this, we tune the nonequilibrium driving $G_{\mathrm{ATP}}$ and track how it affects the CheY-P response amplitude and the half-maximum aspartate concentration ($k_{1/2}$) for each methylation level. To demonstrate their typical behaviors, we consider two paths shown by the purple arrows in Fig. 4A, with the corresponding responses shown in Fig. 4B and C. Both paths start from $G_{\mathrm{ATP}}=20$ ($G_{1}=G_{z}=10$) and move towards the equilibrium limit ($G_{\mathrm{ATP}}=0$) by decreasing $G_{1}$ with fixed $G_{z}$ and vice versa. In each case, the amplitude decreases to zero as the dissipation approaches zero. Indeed, the receptor state affects CheY phosphorylation through kinetic control, which has no effect in the equilibrium limit. Furthermore, non-zero dissipation enables the large spread in CheY-P response curves. Because our model has no feedback from the phosphorylation state to the receptor, and therefore no response in equilibrium, the response curves snap together suddenly when $G_{\text{ATP}}=0$. For the response to be nonzero in the equilibrium limit, there has to be feedback from CheY to the receptor. Thus, the relation Eq. (7) holds and the spread between kinase response curves must be comparable to that for binding curves. Interestingly, the $k_{1/2}$ scaling is qualitatively different depending on the path taken to equilibrium: varying $G_{1}$ leads to both non-monotonic behavior, with $k_{1/2}$ initially increasing for each $m$ before decreasing dramatically near equilibrium. Conversely, varying $G_{z}$ produces monotonic scaling with an initial gradual decrease in the sensitivity range before a similar collapse at equilibrium.

The dissipation energies $G_{1}$ and $G_{z}$ are connected to biochemical parameters via the following relationships,

and

where $G^{0}_{\mathrm{ATP}}$ and $G^{0}_{\mathrm{CheY-P}}$ are the phosphate bond energies for ATP and CheY-P, respectively and $[\mathrm{ATP}]_{0}$, $[\mathrm{ADP}]_{0}$ and $[\mathrm{P_{i}}]_{0}$ are the reference concentrations.

A possible scheme to test our model could be to measure kinase dose-response curves while tuning the [ATP]/[ADP] ratio to change $G_{1}$. For a fixed value of $G_{z}$ within the operational regime, the response amplitude and sensitivity range will decrease as the [ATP]/[ADP] ratio decreases. For example, if $G_{z}=10$, as is the case in Fig. 4B, the response amplitude for $m=4$ will decrease by about 5% when the [ATP]/[ADP] ratio decreases by three orders of magnitude ²²2Changes in ATP concentration of nearly this magnitude (400-fold) have been used in studies of the processive motion of molecular motors [63]. Note that this estimate is nearly the worst-case scenario, since our choice of initial $(G_{1},G_{z})$ has $G_{z}$ close to its threshold value. If the true $G_{z}$ is larger than 10, a much larger reduction in amplitude will occur for the same decrease in the [ATP]/[ADP] ratio (e.g., 51% for $G_{z}=13$). Furthermore, since our model is coarse-grained (for example, it ignores intermediate steps in the phosphotransfer between CheA and CheY), the predicted dissipation rate is smaller than that of the real system [35, 36]. Therefore, the real kinase response may be more susceptible to changes in ATP concentration than the model predicts. Another possible scheme to test our model is to vary $G_{z}$ by controlling the inorganic phosphate concentration $[P_{i}]$ as shown in Fig. 4C. Quantitative measurements of the amplitude and sensitivity range using one of the proposed schemes above may help more precisely pinpoint the operational regime of chemoreceptor signaling complexes.

The proposed microscopic mechanism underlying the nonequilibrium allosteric model is summarized in Fig. 5. The average phosphorylation state of the response regulator CheY is controlled by the methyl-accepting chemotaxis proteins (MCP, blue box) through CheA and CheW (green box). The MCP can undergo multiple distinct conformational changes (denoted by $s$ and $\tilde{s}$), each of which affects the phosphorylation rate of CheY in different ways.

The conformation $s$ is predominantly controlled by ligand binding ($\sigma$) via an equilibrium mechanism captured by the MWC energy function $H_{\mathrm{MWC}}$ (methylation may have a minor effect on $s$ through $K_{i}$, $E_{s}$, or both). When $s=0$, the receptor is in the OFF state with almost no kinase activity, which is equivalent to having a very large energy barrier ($E_{B}\gg 1$). When $s=1$, the receptor is in the ON state with a finite (but not unique) kinase activity, whose intensity depends on the methylation $m$. Given that the receptor methylation sites are away from the kinase CheA, we hypothesize that a different conformation state $\tilde{s}$ mediates the allosteric control of methylation on the kinase activity. In particular, increasing $m$ induces a conformational change $\tilde{s}(m)$ that lowers the energy barrier $E_{B}$ and thereby increases the phosphorylation level of the ON state. The change in the barrier height is exactly $\Delta E_{B}(m)$ as defined in our nonequilibrium allosteric model. It is important to note that while $s$ acts like a binary switch that can be described by a two-state model, $\tilde{s}$ may represent multiple or even a continuous spectrum of conformations, which lead to different kinetic rates or equivalently different barrier heights in the CheY phosphorylation reaction.

The idea that ligand binding and receptor methylation induce different conformational changes in the chemoreceptor has been suggested in the experimental literature [44, 45] and is supported by recent in vivo crosslinking measurements of Tsr receptors [46], which found differing receptor structure changes in response to ligand binding and methylation. In particular, the change in conformation (quantified by the fraction of crosslinking products) induced by ligand binding is about twice as large as the corresponding change between methylation levels $m=0$ and $m=4$.

Another line of evidence comes from the dynamical properties of the cytoplasmic helical domains of nanodisc-inserted receptors, which have been measured using electron paramagnetic resonance (EPR) spectroscopy [47, 48]. It was found that increasing methylation reduces the mobility of the helical structure but preserves the mobility difference between companion helices and among different functional regions. In contrast, the conformation change induced by ligand binding has no detectable effect on helical mobility. Since helical mobility is primarily affected by methylation, similar to the conformation $\tilde{s}$ in the model, we conjecture that helical dynamics may play a role in tuning the phosphorylation energy barrier in the receptor ON state.

These findings directly support the idea of distinct conformational changes ($s$ and $\tilde{s}$) used in our model, which leads to the following emerging picture for chemoreceptors operation: ligand binding induces a large conformational change that switches the receptor from active to inactive; methylation leads to more subtle conformational and dynamical changes that modulate the kinetic properties of the active state. More structural insight into the conformation states is needed to more concretely disentangle the differing effects of ligand binding and methylation. For example, repeating recent FRET measurements of nanodisc-inserted Tar receptors [49] at different methylation levels would help to clarify this picture.

Beyond chemotaxis, our model should be broadly applicable to other two-component signaling systems, which have been studied extensively in E. coli and other bacteria [50]. Here, we illustrate these applications with the E. coli oxygen sensor protein DosP [51, 52, 30]. Increasing oxygen concentration promotes DosP’s ability to hydrolyze cyclic di-GMP, an important bacterial second messenger that triggers downstream responses [53, 54]. Previous work measuring both oxygen binding and phosphodiesterase activity of DosP as functions of oxygen concentration reported that the two measurements cannot be consistently explained by an equilibrium model [30]. Indeed, inconsistency with equilibrium MWC models can be more systematically demonstrated using the parametric test introduced in this work (Fig. 6A). A slightly generalized version of the nonequilibrium model (see SI Appendix) is able to capture both the binding and the phosphodiesterase activity curves. As shown in Fig. 6B, it captures both the sharpness and sensitivity range difference between the binding and activity curves. In contrast, as established by the parametric test, the MWC model can only capture either but not both response curves (see SI Appendix).

Since many of these two-component systems involve the hydrolysis of energy-rich molecules such as NTP, we expect simultaneous measurements of ligand binding and downstream activity dose-response curves to reveal potential inconsistencies with equilibrium models that are widely adopted. Indeed, inconsistency with equilibrium models has also been reported for the FixL/FixJ system [55]. Large discrepancies between the half-maximal concentrations for binding and activity have also been observed for the PhoP/PhoQ system [56], suggesting that it operates out of equilibrium. Overall, the theoretical framework presented here enables a deeper understanding of the mechanism of two-component systems by fitting simultaneously to measurements of binding and enzymatic activity.

More broadly, we expect that the nonequilibrium allosteric model developed here may be useful for understanding other biological signaling systems such as the G-protein coupled receptor (GPCR) signaling pathways if the receptor binding and the output activity can be measured and modeled simultaneously.

Besides application to other two-component systems, another extension of this work is to better incorporate structural information. In the introduction, we noted that the MWC model has two major simplifications: it is equilibrium and it neglects the rich spatial structure of large protein complexes. We have addressed the first of these shortcomings by developing a minimal nonequilibrium model of allostery that is sufficient to explain most of the existing measurements of chemotactic signaling. Looking forward, incorporating spatial structure into the nonequilibrium model will allow us to probe one of the most important questions in biology: how does structure determine function? As the structural information for these complexes becomes available thanks to the recent development of high-resolution techniques such as cryo-electron tomography, the challenge is how this structural information can be used to understand the functions of these large complexes. E. coli chemoreceptor clusters form a beautiful hexagonal lattice structure by tiling core units made up by two trimers of receptor dimers bound with a CheA dimer and two CheW [58, 12]. Developing a nonequilibirum model containing these spatial details will be essential for understanding a variety of recent experimental discoveries including asymmetric activity switching times [59, 60], emergent asymmetric coupling in mixed receptor complexes [61], and slow logarithmic decay in activity on long timescales in saturating ligand [62].