Universal thermodynamic bounds on nonequilibrium response with biochemical applications

TL;DR

This paper derives universal thermodynamic bounds that constrain the response of nonequilibrium systems, extending fluctuation-dissipation concepts far from equilibrium, with applications to biochemical systems.

Contribution

It introduces a set of equalities and inequalities that apply to nonequilibrium steady states, providing new tools to analyze their response characteristics.

Findings

Bounds on nonequilibrium response derived

Application to biochemical motifs demonstrated

Framework extends fluctuation-dissipation relations

Abstract

Diverse physical systems are characterized by their response to small perturbations. Near thermodynamic equilibrium, the fluctuation-dissipation theorem provides a powerful theoretical and experimental tool to determine the nature of response by observing spontaneous equilibrium fluctuations. In this spirit, we derive here a collection of equalities and inequalities valid arbitrarily far from equilibrium that constrain the response of nonequilibrium steady states in terms of the strength of nonequilibrium driving. Our work opens new avenues for characterizing nonequilibrium response. As illustrations, we show how our results rationalize the energetic requirements of two common biochemical motifs.

Equations

$${\dot{p}}_i(t)=\sum _{j=1}^N{W}_{ij}{p}_j(t),$$

$$F_C=\ln \left(\frac{W_{i_0{i}_m}\cdots W_{i_2{i}_1}{W}_{i_1{i}_0}}{W_{i_m{i}_0}\cdots W_{i_1{i}_2}{W}_{i_0{i}_1}}\right).$$

$${\partial }_{\lambda }⟨Q{⟩}_{\pi }=\sum _i{Q}_i{\partial }_{\lambda }{\pi }_i=\sum _i{Q}_i\sum _{kl}\frac{\partial W_{kl}}{\partial \lambda }\frac{\partial {\pi }_i}{\partial W_{kl}}.$$

$${\partial }_{\lambda }⟨Q{⟩}_{\mathrm{eq}}=\beta C_{\mathrm{eq}}(Q,V),$$

$${\partial }_{\lambda }⟨Q{⟩}_{\pi }=C_{\pi }(Q,{\partial }_{\lambda }\ln \pi ),$$

$$W_{ij}=\exp \left[-(B_{ij}-E_j-F_{ij}/2)\right],$$

$$B_{ij}=\frac{1}{2}\left[(E_j-\ln W_{ij})+(E_i-\ln W_{ji})\right]$$

$$F_{ij}=-\left[(E_j-\ln W_{ij})-(E_i-\ln W_{ji})\right].$$

$$\frac{\partial {\pi }_k}{\partial E_j}=\sum _i{W}_{ij}\frac{\partial {\pi }_k}{\partial W_{ij}}.$$

$$\frac{\partial {\pi }_k}{\partial B_{ij}}=-W_{ij}\frac{\partial {\pi }_k}{\partial W_{ij}}-W_{ji}\frac{\partial {\pi }_k}{\partial W_{ji}}$$

$$\frac{\partial {\pi }_k}{\partial F_{ij}}=\frac{1}{2}\left(W_{ij}\frac{\partial {\pi }_k}{\partial W_{ij}}-W_{ji}\frac{\partial {\pi }_k}{\partial W_{ji}}\right).$$

$$\frac{\partial {\pi }_k}{\partial \lambda }=\sum _j\frac{\partial E_j}{\partial \lambda }\frac{\partial {\pi }_k}{\partial E_j}+\sum _{i>j}\frac{\partial B_{ij}}{\partial \lambda }\frac{\partial {\pi }_k}{\partial B_{ij}}+\sum _{i>j}\frac{\partial F_{ij}}{\partial \lambda }\frac{\partial {\pi }_k}{\partial F_{ij}}$$

$$\frac{\partial {\pi }_i}{\partial E_k}=\{\begin{array}{ll}-{\pi }_i(1-{\pi }_i)& \text{mbox}ifi=k\\ {\pi }_i{\pi }_k& \text{mbox}ifi\ne k\end{array}.$$

$$\frac{\partial \ln \left({\pi }_i/{\pi }_j\right)}{\partial E_k}=\{\begin{array}{ll}-1& \text{mbox}ifi=k\\ 1& \text{mbox}ifj=k\\ 0& \text{otherwise}\end{array},$$

$${\partial }_{\lambda }⟨Q{⟩}_{\pi }=\beta C_{\pi }(Q,V),$$

$$\left|\frac{\partial {\pi }_i}{\partial B_{mn}}\right|\le {\pi }_i(1-{\pi }_i)\tanh \left(F_{\max }/4\right)$$

$$\left|\frac{\partial ({\pi }_i/{\pi }_j)}{\partial B_{mn}}\right|\le \left(\frac{{\pi }_i}{{\pi }_j}\right)\tanh \left(F_{\max }/4\right),$$

$$F_{\max }=\underset{C\ni e_{mn}}{\max }|F_C|$$

$$\begin{array}{rl}\left|\sum _{e_{mn}\in S}\frac{\partial \ln ({\pi }_i/{\pi }_j)}{\partial B_{mn}}\right|& \le \sum _{e_{mn}\in S}\left|\frac{\partial \ln ({\pi }_i/{\pi }_j)}{\partial B_{mn}}\right|\\ & \le |S|\tanh (F_{\max }/4).\end{array}$$

$$\left|\sum _{e_{mn}\in S}\frac{\partial }{\partial B_{mn}}\ln \left(\frac{{\pi }_i}{{\pi }_j}\right)\right|\le \left(|W|-1\right)\tanh \left(\frac{F_{i\leftrightarrow j}}{4}\right).$$

$$\left|\frac{\partial {\pi }_i}{\partial F_{mn}}\right|\le {\pi }_i(1-{\pi }_i)\le \frac{1}{4},$$

$$s=\left|\frac{\partial \ln \left({\pi }_{S^{\ast }}/{\pi }_S\right)}{\partial \ln [E_1]}\right|\le \tanh (\Delta \mu /4k_{\mathrm{B}}T),$$

$${\sigma }_e^2\approx \frac{{\sigma }_s^2}{(\partial {\mu }_s/\partial [E_1]{)}^2},$$

$${\left(\frac{{\sigma }_e}{[E_1]}\right)}^2\gtrsim \frac{4}{N\tanh (\Delta \mu /4k_{\mathrm{B}}T{)}^2}.$$

$$\nu =-\frac{\partial \ln ({\pi }_E/{\pi }_{ES})}{\partial \Delta }.$$

$$|\nu -1|\le \tanh (\Delta \mu /4k_{\mathrm{B}}T),$$

$$\left|\nu -1\right|\le (m-1)\tanh (F_{E\leftrightarrow ES}/4),$$

$$w(S)=W_{ji}{W}_{lk}\cdots .$$

$${\pi }_k=\frac{1}{\mathcal{N}}\sum _{\begin{array}{c}\text{spanning trees}\\ \text{T of G}\end{array}}w(T_k),$$

$$\frac{\partial {\pi }_i}{\partial E_k}=\{\begin{array}{ll}-{\pi }_k(1-{\pi }_k)& \text{mbox}ifi=k\\ {\pi }_k{\pi }_i& \text{mbox}ifi\ne k\end{array}.$$

$${\pi }_k=\frac{a}{a+b{e}^{E_k}},$$

Full text

I Introduction

One of the basic characteristics of any physical system is its response to small perturbations Chaiken and Lubensky (1995). For instance, response is used to quantify everything from material properties—such as conductivity Kubo et al. (1985) and viscoelasticity Mason and Weitz (1995)—to the sensing capability of cells Qian (2012); Govern and ten Wolde (2014a) and the accuracy of biomolecular processes Bintu et al. (2005); Murugan et al. (2011); Estrada et al. (2016). Near thermodynamic equilibrium, response is completely determined by the nature of spontaneous fluctuations, according to the fluctuation-dissipation theorem (FDT) Kubo et al. (1985). This deep connection between response and fluctuations is not only of theoretical interest, but also finds practical application. The FDT forms the basis of powerful experimental techniques, such as microrheology, spectroscopy, and dynamic light scattering Chaiken and Lubensky (1995). It additionally has implications for the design of mesoscopic devices: highly-responsive equilibrium devices are always plagued by noise. To combine low fluctuations with high sensitivity, a device must be driven away from equilibrium.

The great utility of the FDT near equilibrium Chaiken and Lubensky (1995) has led to significant interest in expanding its validity and developing generalizations for nonequilibrium situations. Generically, response can be related to some formal nonequilibrium correlation functions Agarwal (1972); Risken (1984); Baiesi et al. (2009); Prost et al. (2009); Seifert and Speck (2010). While these predictions offer fundamental theoretical insight, the necessary correlations are often prohibitively difficult to measure except in simple single-particle systems Gomez-Solano et al. (2009); Mehl et al. (2010); Bohec et al. (2013). In certain special cases, such as under stalling conditions, the FDT holds unmodified Altaner et al. (2016). More commonly, however, the study of nonequilibrium response has focused on how the FDT is violated. For example, the violation of the velocity FDT for Brownian particles can be related to the steady-state heat dissipation through the Harada-Sasa equality Harada and Sasa (2005, 2006); a useful prediction that has been utilized to measure dissipation and efficiency in molecular motors Toyabe et al. (2010). Alternatively, violations of the FDT can be used to fit model parameters, as has been suggested for models of biomolecular processes Yan and Hsu (2013); Sato et al. (2003). More often FDT violations are framed in terms of system-specific “effective temperatures” Cugliandolo (2011); Ben-Isaac et al. (2011); Dieterich et al. (2015), whose time-dependence under some circumstances can reveal information about effective equilibrium descriptions.

Inspired by the recent demonstration of thermodynamic bounds on far-from-equilibrium dynamical fluctuations Barato and Seifert (2015); Gingrich et al. (2016), we show here that generic nonequilibrium steady-state response can be constrained in terms of experimentally-accessible thermodynamic quantities. In particular, we present equalities and inequalities—akin to the FDT but valid arbitrarily far from equilibrium—that link static response to the strength of nonequilibrium driving. Our results open new possibilities to experimentally characterize away-from-equilibrium response and suggest design principles for high-sensitivity, low-noise devices. As illustrations, we show how our results rationalize the energetic requirements of biochemical switches, biochemical sensing, and kinetic proofreading.

II Modeling nonequilibrium steady states

Nonequilibrium steady states are characterized by the constant and irreversible exchange of energy and matter between a system and its environment. These flows are driven by thermodynamic affinities—quantities like temperature gradients, chemical potential differences and nonconservative mechanical forces. The underlying dynamics leading to the establishment of such steady states are often well-modeled as a continuous-time Markov jump process on a finite set of states $i=1,\dots,N$, which represent (coarse-grained) physical configurations. The probability $p_{i}(t)$ to find the system in state $i$ at time $t$ then evolves according to the master equation

[TABLE]

where the off-diagonal entries of the transition rate matrix $W_{ij}$ specify the probability per unit time to jump from $j$ to $i$, and diagonal entries $W_{ii}=-\sum_{j\neq i}W_{ji}$ are fixed by the conservation of probability. Time-reversibility of the underlying microscopic dynamics implies that $W_{ij}\neq 0$ only if $W_{ji}\neq 0$ Seifert (2012). We will additionally suppose that for any two states, there is some sequence of allowed transitions ($W_{ij}\neq 0$) connecting them, a property known as irreducibility. Under this assumption, no matter the initial condition, the solution of the master equation (1) converges at long times to the unique steady state distribution $\pi$ that satisfies $\sum_{j=1}^{N}W_{ij}\pi_{j}=0$. This distribution $\pi$, and in particular its dependence on physical quantities through the transition rates $W_{ij}$, serves as a general model of a nonequilibrium steady state.

Many key properties of this nonequilibrium steady state, including its thermodynamics, come to light by picturing the stochastic dynamics described by (1) playing out on a transition graph—a weighted directed graph $G$, as in Fig. 1$(a)$, where the vertices $\{i\}$ represent the states and directed edges $\{e_{mn}\}$ represent possible transitions, weighted by the rates $W_{mn}$.

Note that, by assumption, every edge in $G$ has a reverse, so we will represent and discuss the transition graph as if it were an undirected graph, with the understanding that every undirected edge represents two opposing directed edges.

The cycles in the graph, like those in Fig. 1$(b)$, play a central role in the thermodynamics of the steady state. A cycle is a sequence of directed edges and vertices connecting the initial vertex to itself without self-intersecting, $C=\{i_{0}\xrightarrow{e_{i_{1}i_{0}}}i_{1}\to\cdots\to i_{m}\xrightarrow{e_{i_{0}i_{m}}}i_{0}\}$. The asymmetry of the rates around these cycles then encodes the thermodynamic affinities driving the system out of equilibrium through the cycle forces—the log of the product of rates around the cycle divided by the product of rates in the reverse orientation Schnakenberg (1976); Andrieux and Gaspard (2006):

[TABLE]

These cycle forces are linear combinations of thermodynamic affinities multiplied by their conjugate distances—for example a chemical potential gradient times a change in particle number. As such, the cycle forces equal the dissipation (entropy production) in the environment accrued every time the system flows around the cycle $C$. This means that the cycle forces depend on macroscopically tunable parameters—such as environmental temperature or chemical potential—that characterize how strongly the system is driven away from equilibrium. If all the cycle forces vanish, the system satisfies detailed balance, a statistical time-reversal symmetry Zia and Schmittmann (2007) characteristic of thermodynamic equilibrium.

III Static response to perturbations

Now, suppose the transition rates $W_{ij}(\lambda)$ depend on a control parameter $\lambda$, which could represent, say, the strength of an applied electric field, a temperature, or even a microscopic kinetic parameter such as a reaction barrier. In this work, we focus on the response to static perturbations, that is how steady-state averages $\langle Q\rangle_{\pi}=\sum_{j}Q_{j}\pi_{j}$ respond to small changes in $\lambda$:

[TABLE]

At thermal equilibrium, the steady state $\pi^{\rm eq}_{i}\propto e^{-\beta\epsilon_{i}(\lambda)}$ depends only on the underlying (free) energy landscape $\epsilon_{i}(\lambda)$, irrespective of the precise form of the transition rates, where $\beta=1/k_{\rm B}T$ with $k_{\rm B}$ Boltzmann’s constant and $T$ temperature. This simplifying fact immediately implies the static FDT, which equates the static response to an equilibrium correlation function,

[TABLE]

where $C_{\rm eq}(Q,V)=\langle QV\rangle_{\rm eq}-\langle Q\rangle_{\rm eq}\langle V\rangle_{\rm eq}$ and the subscript “eq” emphasizes that averages are taken with respect to the equilibrium distribution Kubo et al. (1985). Here, $V=-\partial_{\lambda}\epsilon$ is known as the coordinate conjugate to $\lambda$ and represents the displacement induced by $\lambda$—for example, volume is conjugate to pressure and particle number is conjugate to chemical potential. The FDT’s utility in part stems from the fact that we often know the conjugate coordinate from basic physical reasoning and it is easily measured.

Away from equilibrium, the steady-state distribution generally has a complicated dependence on the rates. Nevertheless, response can always be related to a nonequilibrium correlation function Agarwal (1972); Risken (1984),

[TABLE]

but the relevant “conjugate coordinate” $\partial_{\lambda}\ln\pi$ requires knowledge of the parameter-dependence of the full nonequilibrium steady-state distribution, which can be challenging to calculate or measure. Still, this response formula has been given thermodynamic meaning by relating the nonequilibrium conjugate coordinate to the stochastic entropy production rate Seifert and Speck (2010) as well as to the time-reversal symmetry properties of the path action Baiesi et al. (2009).

IV Parameterizing Perturbations

Here, we turn our attention to the variations of the steady-state distribution with the transition rates, $\partial\pi_{i}/\partial W_{kl}$, constraining them in terms of $\pi$ and the cycle forces $F_{C}$. We accomplish this goal by breaking any general perturbation into a linear combination of three special types of perturbations that change rates in a coordinated way. By focusing on these subclasses of perturbations, we will be able to provide clear and measurable thermodynamic constraints on static response.

To classify perturbations, it will prove fruitful to parameterize the rate matrix as

[TABLE]

introducing the vertex parameters $E_{j}$, (symmetric) edge parameters $B_{ij}=B_{ji}$, and asymmetric edge parameters $F_{ij}=-F_{ji}$, which can all be varied independently. Any rate matrix can be cast in this form, albeit non-uniquely. To see this, consider the following program for identifying such a parameterization: choose the vertex parameters $\{E_{1},\dots,E_{N}\}$ arbitrarily, then set

[TABLE]

The non-uniqueness of the parameterization is manifest in this construction because of the freedom to choose the $E_{i}$. For example, we could choose $E_{i}=0$ for all $i$. We emphasize, however, that no matter the choice, the derivatives with respect to $E_{j}$, $B_{ij}$ or $F_{ij}$ are independent of the values of the parameterization. Variations of a vertex parameter, say $E_{j}$, is equivalent to scaling all of the rates out of state $j$

[TABLE]

Similarly, derivatives with respect to $B_{ij}$ and $F_{ij}$ multiplicatively scale transition rates associated with a single edge

[TABLE]

[TABLE]

These facts are illustrated in Fig. 2. We note that the right hand sides of equations (9)–(11) do not depend on the choice of parameterization (6), even though the parameterization is not unique.

Our parameterization is reminiscent of the Arrhenius expression for transition rates for a system evolving in an energy landscape with wells of depth $E_{i}$ and barriers of height $B_{ij}$ driven by forces $F_{ij}$. While we stress that (6) will not in general support such an interpretation, the analogy is suggestive in several ways. For example, the asymmetric edge parameters $F_{ij}$ are the sole contributors to the cycle forces (affinities) $F_{C}=\sum_{e_{ij}\in C}F_{ij}$. Furthermore, if all the $F_{ij}=0$, the steady-state distribution has the Boltzmann form $\pi_{i}\propto\exp(-E_{i})$, with the $E_{i}$ acting as a dimensionless energy.

Our main results are a series of simple thermodynamic equalities and inequalities for how the steady state responds to perturbations of the $E_{j}$, $B_{ij}$ and $F_{ij}$. By combining these results, we can constrain the response to any arbitrary perturbation of the rates through a decomposition of the form

[TABLE]

The freedom in the rate parameterization makes this decomposition non-unique, and the tightness of our inequalities will depend on the specific decomposition. For example, one choice could force $\partial_{\lambda}E_{i}=0$ for all $i$. We are not presently aware of a good strategy to identify the decomposition which yields optimally tight inequalities. In this work, we show that simple decompositions can nevertheless yield interesting bounds.

In deriving our response results, the basic mathematical tool we rely on is the matrix-tree theorem (MTT), presented in Appendix A, which gives an exact algebraic expression for the steady-state probabilities $\pi_{i}$ in terms of the rates $W_{ij}$ Schnakenberg (1976); Hill (1977). All our results—presented in the following sections—are obtained by reasoning about the result of differentiating the expression given by the MTT. Proofs are given in the appendices.

V Vertex perturbations

Our first main result is the exact expression for the response to a vertex perturbation (Appendix B)

[TABLE]

We stress that the $B_{ij}$ and $F_{ij}$ are unrestricted, so this equality holds even for nonequilibrium steady states. As an immediate consequence, we also find that for $i\neq j$,

[TABLE]

which implies that the relative probability between two states is insensitive to the vertex parameters elsewhere in the graph.

Remarkably, these equalities are exactly equivalent to the response of a Boltzmann distribution to energy perturbations, which leads to the surprising conclusion that far-from-equilibrium response has an equilibrium-like structure if the perturbation leaves the $B_{ij}$ and $F_{ij}$ fixed. To leverage this observation, let us assume that we vary the rates only through the system’s energy function $\epsilon_{i}(\lambda)$ and that the rates depend on the energy as $W_{ij}=\omega_{ij}e^{\beta\epsilon_{j}}$, with arbitrary energy-independent $\omega_{ij}$. Comparison with (6), shows that variations in the energy $\epsilon_{i}$ in this case can be parameterized as vertex parameters $E_{i}$. Then (13) implies that arbitrarily far from equilibrium the response maintains the equilibrium-like form of the FDT,

[TABLE]

with the response proportional to the nonequilibrium steady-state correlation with the coordinate conjugate to the energy $V=-\partial_{\lambda}\epsilon$ [cf. (4) and (5)]. This prediction implies that experimental verification of the static FDT is not sufficient to conclude that a system is in equilibrium.

VI Symmetric edge perturbations

More generally, a perturbation will modify not only the vertex parameters $E_{i}$, but also the edge parameters $B_{ij}$. While at equilibrium the steady state is independent of edge parameters $B_{ij}$, this is generically not the case out of equilibrium. In this section, we demonstrate that in fact response to edge perturbations is constrained by thermodynamic affinities through the cycle forces. For proofs of the results in this section, see Appendix C.

VI.1 Single edge

Our second main result is that the response to the perturbation of a symmetric edge parameter $B_{mn}$, associated to a single edge $e_{mn}$, is constrained by the cycle forces:

[TABLE]

where

[TABLE]

is the maximum cycle force over all cycles that contain the (undirected) perturbed edge $e_{mn}$ (illustrated in Fig. 3). If the cycle forces all equal zero—as they must at equilibrium—then $F_{\rm max}=0$, and the response is zero, as expected. In addition, only perturbations of an edge contained in a cycle can induce a response: perturbations of edges whose removal would disconnect $G$ therefore cannot alter the steady state. Equation (16), furthermore, has the character of the FDT, once we recognize $\pi_{i}(1-\pi_{i})$ as the variance of the occupation fluctuations of state $i$; thus, we see a manifestation of how thermodynamics shapes the interplay between response and fluctuations. These inequalities, applying to all discrete stochastic dynamics, significantly generalize a bound for two-state systems derived by Hartich et al. in a model of nonequilibrium sensing Hartich et al. (2015).

The conditions for equality in (16) and (17) would suggest methods for designing optimized or highly responsive devices. As detailed in Appendix E, we can exhibit at least one scenario that does saturate (17): a single cycle with strong time-scale separation so that the system effectively only has two states. This limiting scenario suggests that small single cycle systems are ideal for optimizing response. Systematically deducing the system parameters that saturate our inequalities in general remains for future work.

VI.2 Multiple edges

The response to a perturbation of multiple edge parameters can be bounded by combining (17) with the triangle inequality. For example, for any set $S$ of $|S|$ edges,

[TABLE]

It is clear, however, that this inequality is not always the best we can do. Consider for example the case where $S$ consists of every edge in $G$. In this case, increasing all the edge parameters by the same amount (which is what the sum above amounts to) is like rescaling time, which cannot affect $\pi$ and therefore has zero response.

In this section, we provide a different bound on the response to a perturbation of multiple edge parameters that in many cases improves on (19). Suppose we vary the edge parameters associated to the edges $S$ of a subgraph $H\neq G$. Let $W$ be the set of vertices of $H$ that connect it to the rest of the graph (i.e. the set of vertices of $H$ incident to an edge not in $H$). Then,

[TABLE]

where $F_{i\leftrightarrow j}$, defined precisely in Appendix C, can be physically identified as the largest possible entropy produced in the environment when the system goes from $i$ to $j$ and back again (along paths without self-intersection). Whenever there is only one path through state space between $i$ and $j$, and in all cases at thermodynamic equilibrium, $F_{i\leftrightarrow j}=0$.

We finally note that our results (16), (17), and (20) admit generalization to the response of a ratio of positive observables $\langle A\rangle/\langle B\rangle$. In this general case, the bounds remain unchanged, except that $F_{i\leftrightarrow j}$ is replaced by its maximum over all pairs of vertices $i,j$.

VII Asymmetric edge perturbations

Lastly, the MTT allows us to bound the response to asymmetric edge perturbations as (Appendix F)

[TABLE]

which is related to, but distinct from, inequalities established in Thiede et al. (2015). This result is a consequence of an identical inequality that holds for general rate perturbations $\partial\pi_{i}/\partial\ln W_{jk}$.

Any perturbation of rates can be decomposed into a linear combination of perturbations of the vertex and edge parameters $E_{i}$, $B_{ij}$, and $F_{ij}$ we have introduced, and so the response can be bounded using our inequalities (via the triangle inequality). What our results in this section show is that even for a general perturbations, there is a universal bound—the response to the variation of a single rate is bounded by a constant independent of the structure of $G$ or rates of other transitions, meaning that high sensitivity always requires many different transitions to be perturbed, their cumulative effect generating a response that can greatly exceed $1/4$.

VIII Biochemical applications

In this section, we illustrate the use of our main results by detailing applications to well-studied motifs appearing in biochemical networks.

VIII.1 Covalent modification cycle

First, we consider a well-studied model Goldbeter and Koshland (1981); Qian (2007); Xu and Gunawardena (2012); Dasgupta et al. (2014) of a biological switch—the modification/demodification cycle depicted in Fig. 4, also known as the Goldbeter–Koshland loop Goldbeter and Koshland (1981); Xu and Gunawardena (2012), or “push–pull” network Govern and ten Wolde (2014b).

The network consists of a substrate with two forms, an “unmodified” $S$ and “modified” $S^{*}$, along with enzymes, $E_{1}$ and $E_{2}$, that actively catalyze its modification and demodification, respectively. For example, if $E_{1}$ is a kinase, $E_{2}$ a phosphatase, and $S^{*}$ a singly-phosphorylated form of $S$, then the system is driven by the chemical potential gradient $\Delta\mu=\mu_{\rm ATP}-\mu_{\rm ADP}-\mu_{{\rm P}_{i}}$ for ATP hydrolysis. In the limit in which the substrate is very abundant compared to its modifying enzymes, it is well known that such a system can exhibit unlimited sensitivity to changes in the ratio of the concentrations of the modifying and demodifying enzymes Goldbeter and Koshland (1981).

In the other limit—that of a single substrate molecule—our results (17) limit the sensitivity of the ratio $\pi_{S^{*}}/\pi_{S}$ for a particular substrate molecule to changes in the enzyme concentration (Appendix G):

[TABLE]

where $F_{\rm max}=\Delta\mu/k_{\rm B}T$ is the single chemical driving force. For the simple cycle in Fig. 4(a) where each enzyme has a single intermediate and we assume mass-action kinetics, this result arises from unraveling a change in $[E_{1}]$ as a change in the vertex parameter associated to $E_{1}S$, together with changes in the parameters of the edges connecting $E_{1}S$ to $S$ and $E_{1}S$ to $S^{*}$.

In fact, inequality (22) turns out to hold under assumptions much more general than those of Fig. 4(a). It remains true even if catalysis by $E_{1}$ and $E_{2}$ proceeds via any number of intermediate complexes with arbitrary rates as in Fig. 4(b), as long as there is no irreversible formation of a dead-end complex and the chemical driving is the same around every cycle in which $E_{1}$ makes the modification of $S$ and $E_{2}$ removes it Xu and Gunawardena (2012); Dasgupta et al. (2014). In this general case, the many perturbed vertex and edge parameters (Fig. 4(b)) form a subgraph that acts effectively as a single edge perturbation. Our multi-edge bound (20) then applies with $F_{S\leftrightarrow S^{*}}=\Delta\mu/k_{\rm B}T$ being the maximum entropy produced to go from the unmodified $S$ to the modified $S^{*}$ form and back again.

In the absence of nonequilibrium drive ($\Delta\mu=0$), it is clear this switch cannot work, because it operates by varying the kinetics via an enzyme concentration, and at equilibrium the steady state is independent of kinetics. It has long been known that switches require energy Goldbeter and Koshland (1987); Qian (2007); Tu (2008). Our results provide a general quantification of this requirement.

VIII.2 Biochemical sensing

The covalent modification cycle, discussed in the previous section, is an integral component of numerous biochemical models for cellular sensing  Berg and Purcell (1977); Bialek and Setayeshgar (2005); Kaizu et al. (2014); Govern and ten Wolde (2014a, b). So far, we have described a single substrate molecule stochastically switching between its two forms due to the action of abundant enzymes $E_{1}$ and $E_{2}$. Here, we imagine there are $N$ substrate molecules, which act as $N$ independent copies of the system studied above, as long as the numbers of both enzymes greatly exceeds $N$. Then the number $s$ of modified substrate molecules $S^{*}$ can be interpreted as a noisy readout of the enzyme concentration $[E_{1}]$. The random variable $s$ is binomially distributed, with mean $\mu_{s}=N\pi_{S^{*}}$ and variance $\sigma_{s}^{2}=N\pi_{S^{*}}(1-\pi_{S^{*}})$, which implicitly depend on $[E_{1}]$. Thus, this scenario provides a mechanism to measure a chemical concentration, by exploiting the relation between $s$ and $[E_{1}]$.

Now suppose one makes the observation at some time that there are $\tilde{s}$ molecules of $S^{*}$. Supposing $[E_{2}]$ and all rate constants assume known, fixed values, one can produce an estimate $\hat{[E_{1}]}$ of $[E_{1}]$ by choosing $\hat{[E_{1}]}$ to be the value of $[E_{1}]$ that gives $\mu_{s}(\hat{[E_{1}]})=\tilde{s}$. The variance of the estimate $\hat{[E_{1}]}$ so constructed is often well-approximated, when the noise is small, by Govern and ten Wolde (2014a, b)

[TABLE]

where the quantities on the right hand side should be evaluated at the true concentration $[E_{1}]$. Our result (16) combined with the probabilistic inequality $\pi_{S^{*}}(1-\pi_{S^{*}})\leq 1/4$ then leads to the following bound on the relative error

[TABLE]

This result interpolates between bounds on error established by Govern and ten Wolde Govern and ten Wolde (2014a, b) in two limits of resource limitation in cellular sensing systems. That work studies a model of sensing in which cell surface receptors bind to an extracellular ligand whose concentration the cell needs to determine. The ligand-bound receptors then participate in a modification/demodification cycle like the one we study here, playing the role of $E_{1}$. See Appendix G.2 for details.

VIII.3 Kinetic proofreading

As a third application, we turn to the effectiveness of kinetic proofreading Hopfield (1974); Ninio (1975). A common challenge faced by biomolecular processes is that of discriminating between two very similar chemical species. At equilibrium, the probability of an enzyme $E$ being bound to a substrate $S$, divided by the probability of that enzyme being free is $\exp(-\Delta)$, where $k_{B}T\Delta$ is the binding (free) energy of the enzyme-substrate complex $ES$. Two substrates with very similar binding energies are constrained to be bound by the enzyme a similar fraction of time.

Kinetic proofreading is a scheme to use nonequilibrium driving to improve discrimination based on binding energy. One way to quantify the discriminatory ability of a kinetic network is using the discriminatory index introduced by Murugan et al. Murugan et al. (2014),

[TABLE]

At equilibrium, $\nu=1$. The simplest nonequilibrium scheme to improve on this is the single-cycle network illustrated in Fig. 5(a). Note that we have supposed the binding energy $\Delta$ appears exclusively in the unbinding rates. Hopfield observed that in a certain nonequilibrium limit of the rates, $\nu\to 2$ Hopfield (1974). Our results lead to a constraint on $\nu$ that interpolates between the equilibrium case and this limit.

In the single-cycle network, the variation of the binding energy $\Delta$ is equivalent to the variation of two vertex parameters (that of $ES$ and $ES^{*}$) and a single edge parameter ($ES\leftrightarrow ES^{*}$), leading to the inequality (Appendix G),

[TABLE]

where $F_{\mathrm{max}}=\Delta\mu/k_{\rm B}T$ is the chemical driving around the cycle. This bound, which can be saturated, reduces correctly to $\nu=1$ at equilibrium and is consistent with $\nu\to 2$ in the limit of strong driving $\Delta\mu\to\infty$.

We can also bound $\nu$ in the case of a more general kinetic proofreading scheme Murugan et al. (2012, 2014) in which there are $m$ complexes that can dissociate. Each of the dissociation transitions can be thought of as crossing a “discriminatory fence” Murugan et al. (2014), its rate depending on the binding energy $\Delta$, as in Fig. 6.

We suppose that the dissociation transitions are the only ones that depend on $\Delta$. We make no assumptions about the structure of the transition graph on either side of the fence. In such a network, perturbing $\Delta$ is equivalent to perturbing the edge and vertex parameters on one side of the “fence” forming the subgraph highlighted in blue in Fig. 6. We then have (Appendix G)

[TABLE]

where $F_{E\leftrightarrow ES}$ is the maximum entropy produced to go from $E$ to $ES$ and back again. Notably, we recover the result $\left|\nu-1\right|\leq m-1$ of Murugan et al. (2014).

IX Conclusion

In this work, we have developed a series of universal bounds on nonequilibrium response in terms of the strength of the nonequilibrium driving. We show that for a large class of static perturbations, a result equivalent to the FDT continues to hold out of equilibrium. For many other perturbations, we bound the response in terms of the dimensionless thermodynamic forces, which quantify departure from equilibrium.

The illustrations detailed in the previous section demonstrate the potential of our results to unify long-standing observations about the importance of energy “expenditure” in many different models. The tasks of making a sharp molecular switch, a good sensor, or discriminating between two similar ligands, all have in common the need for a large response to a small perturbation. We find new bounds interpolating between known limits in these systems, and show how they all descend from our results on vertex and edge perturbations.

A more detailed analysis of the conditions under which our bounds are saturated would lead to design principles for optimal response. Our preliminary investigation identified single cycles as ideal when a single edge parameter is varied. We expect that for more complex perturbations, the most highly responsive systems may have a more complicated structure.

An important theme highlighted by our work is that sensitivity is limited not only by nonequilibrium driving, but also, very strictly, by network size and structure. The total number of transitions in a biochemical network limits response, because the response to the scaling of any one rate is bounded by $1/4$. At the same time, our multi-edge results show how many enlargements or complications of networks (e.g. departure from Michaelis-Menten assumptions in the covalent modification cycle), do not confer any advantage. In this sense, our results build on the work of others who studied similar questions in the context of kinetic proofreading Murugan et al. (2014) and biochemical copy processes Ouldridge et al. (2017).

Our results point to numerous other extensions, including bounds on the response of currents with implications for the Green-Kubo and Einstein relations Seifert (2010); Baiesi et al. (2011); Dechant and Sasa . We have also focused on results that hold in general, not taking into account possible characteristic structures in the graph of states and transitions, which are present for example in many natural examples, such as chemical reaction networks. The study of such extensions and special cases strike us as promising directions for future work.

Acknowledgments

We would like to thank Alexandre Solon and Matteo Polettini for very useful discussions. JAO would like to thank Jeremy England for advice and support.

Appendix A Matrix-tree theorem

The key tool that we apply in our analysis of nonequilibrium response is the matrix-tree theorem (MTT). To state the theorem, we must introduce some additional notation and concepts.

For any set of directed edges $S=\{i\to j,k\to l,\dots\}$, we define the weight $w(S)$ to be the product of the weights of the edges,

[TABLE]

The weight $w(H)$ of any subgraph $H$ we define to be the weight of its edge set.

We also need to introduce spanning trees, which are connected subgraphs of a graph $G$ that contain every vertex, but have no cycles, see Fig. 7.

Every graph that is connected (as is, by assumption, the transition graph of our system) has at least one spanning tree. For any spanning tree $T$ and vertex $r$ of $G$, there is a unique way to direct the edges of $T$ so that they all “point towards” $r$, which we then call the “root”. The resulting directed graph, which we write $T_{r}$, is a rooted spanning tree of $G$. The steady-state distribution $\pi$ is given explicitly by the matrix-tree theorem (MTT) Tutte (1948); Hill (1966); Shubert (1975); Schnakenberg (1976); Leighton and Rivest (1986); Mirzaev and Gunawardena (2013) in terms of weights of rooted spanning trees of $G$.

Theorem** (matrix-tree theorem).**

Let $W$ be the transition rate matrix of an irreducible continuous-time Markov chain with $N$ states. Then the unique steady-state distribution is

[TABLE]

where ${\mathcal{N}}=\sum_{k=1}^{N}\sum_{T}\ w(T_{k})$ is the normalization constant.

This theorem, also known as the Markov chain tree theorem, is a consequence of a result of Tutte Tutte (1948), and has been rediscovered repeatedly in different literatures, see e.g. Hill (1966); Shubert (1975); Schnakenberg (1976); Leighton and Rivest (1986) and Mirzaev and Gunawardena (2013) for further discussion.

The MTT offers a graphical representation of the steady-state distribution that provides a convenient method for organizing the structure of the solution. We illustrate this result in Fig. 8.

Appendix B Vertex perturbations

Theorem 1**.**

[TABLE]

Proof.

The matrix-tree theorem implies that $\pi_{i}$ can be expressed as the ratio of sums of weights of rooted spanning trees. So to evaluate $\partial\pi_{i}/\partial E_{k}$, we need to understand in which spanning trees, and in what form, $E_{k}$ appears. The only rates that depend on $E_{k}$ are rates of transitions out of $k$, $W_{*k}=\exp(E_{k}-B_{*k}+F_{*k}/2)$, see Fig. 2. Furthermore, any rooted spanning tree has exactly one edge directed out of $k$, unless the tree is rooted at $k$, in which case it has none. These observations allow us to group spanning trees in the MTT expression for the steady-state distribution in a convenient manner as illustrated in Fig 9.

Thus, for $i=k$, the matrix-tree theorem implies that

[TABLE]

where

[TABLE]

Here, $a$ is the sum of weights of all spanning trees rooted at $k$—these do not depend on $E_{k}$ since they have no edge directed out of $k$—and $be^{E_{k}}$ is the sum of weights of all spanning trees not rooted at $k$—each of these has exactly one factor of $E_{k}$, making $b$ independent of $E_{k}$.

If $i\neq k$, the MTT yields by a similar argument

[TABLE]

with

[TABLE]

The theorem now follows by differentiating these expressions. When $i\neq k$,

[TABLE]

and similarly for $i=k$. ∎

Corollary 1**.**

If $i\neq j$,

[TABLE]

Proof.

First, note that

[TABLE]

Now we apply Theorem 1. If $i=k$, then $j\neq k$, and $\frac{\partial\ln\left(\pi_{i}/\pi_{j}\right)}{\partial E_{k}}=-(1-\pi_{k})-\pi_{k}=-1$. If $j=k$, then $i\neq k$, and $\frac{\partial\ln\left(\pi_{i}/\pi_{j}\right)}{\partial E_{k}}=\pi_{k}+(1-\pi_{k})=1$. And if neither $i$ nor $j$ equal $k$, then $\frac{\partial\ln\left(\pi_{i}/\pi_{j}\right)}{\partial E_{k}}=\pi_{k}-\pi_{k}=0$. ∎

Appendix C Symmetric edge perturbations

C.1 Single edge

In this appendix, we bound the response to the perturbation of a single symmetric edge parameter in terms of the cycle forces driving the system out of equilibrium.

First, we prove a general bound on the response of a ratio of observables. Equations (16) and (17) will then follow as corollaries by choosing suitable observables.

Theorem 2**.**

Consider any two observables $A,B\in\mathbb{R}_{\geq 0}^{N}$ with at least one positive entry. Then,

[TABLE]

where $F_{\mathrm{max}}$ is the magnitude of the cycle force that is largest in magnitude, among all those associated to cycles containing the distinguished edge $m\leftrightarrow n$ (in either direction).

Our proof relies on the following technical lemma, which we prove in Appendix D.

Lemma 1** (“Tree surgery”).**

Let ${\mathcal{E}}_{mn}$ be the set of spanning trees of $G$ containing the distinguished (undirected) edge $m\leftrightarrow n$. Then for any two distinct vertices $i,j$ of $G$,

[TABLE]

Proof of Theorem 2.

The matrix-tree theorem offers a graphical representation of the steady-state distribution in terms of rooted spanning trees. This observation suggests that we can segregate those contributions to steady-state averages that contain $B_{mn}$ by selecting those (undirected) spanning trees in $G$ that contain the edge $e_{mn}$. Let us call this set ${\mathcal{E}}_{mn}$.

Then by the matrix-tree theorem, we can write

[TABLE]

where

[TABLE]

where $a_{1}$ and $b_{1}$ are linear in $\exp(-B_{mn})$, since they contain edge $e_{mn}$, whereas $a_{0}$ and $b_{0}$ are independent of $B_{mn}$. An illustrative example is presented in Fig. 10.

This implies

[TABLE]

Now note that by the AM-GM inequality the denominator is bounded as

[TABLE]

Since the numerator $b_{0}a_{1}-a_{0}b_{1}=\left(\sqrt{b_{0}a_{1}}-\sqrt{a_{0}b_{1}}\right)\left(\sqrt{b_{0}a_{1}}+\sqrt{a_{0}b_{1}}\right)$, the bound (44) implies,

[TABLE]

To complete the proof, we need to bound the ratio $b_{0}a_{1}/a_{0}b_{1}$ by $\exp\left(F_{\rm max}\right)$. To do this, we match up terms above and below, writing the fraction as

[TABLE]

The desired result is now a consequence of the inequality

[TABLE]

to give

[TABLE]

followed by Lemma 1.

∎

From Theorem 2 we readily obtain our bounds on steady-state response. For (17):

Corollary 2**.**

[TABLE]

Proof.

Choose the observables in Theorem 2 to be $A_{l}=\delta_{il}$ and $B_{l}=\delta_{kl}$, where $\delta_{ij}$ is the Kronecker delta. ∎

We also have:

Corollary 3**.**

Let be $\pi_{X}=\sum_{k\in X}\pi_{k}$ be the total probability of a set of states $X$. Then,

[TABLE]

Proof.

Choose the observables in Theorem 2 to be $A_{i}=\delta_{i}(X)$ and $B_{i}=1-\delta_{i}(X)$, where the indicator $\delta_{i}(X)=1$ if $i\in X$ and $\delta_{i}(X)=0$ otherwise. Note that we then have $\langle A\rangle=\pi_{X}$ and $\langle B\rangle=1-\pi_{X}$. ∎

If $X=\{i\}$ consists of only a single state we recover the bound (16).

C.2 Multiple edges

In this section, we derive our inequality for the response to perturbations by multiple edge parameters (20). The proof proceeds in two steps. We first prove a bound on an arbitrary set of edges $S$ from which (20) and other results are ready corollaries.

Here, the magnitude of response is bounded by a different function $F_{i\leftrightarrow j}$ of cycle forces. The quantity $F_{i\leftrightarrow j}$ is defined for any graph $G$ and vertices $i$ and $j$ to be

[TABLE]

where $P_{i\to j}$ is a (non-self-intersecting) path from $i$ to $j$, $P_{j\to i}$ is a (non-self-intersecting) path from $j$ to $i$, and the superscript ‘$*$’ denotes the reverse path.

Theorem 3**.**

Let $S$ be a set of edges, and define $c_{\mathrm{max}}$ to be the size of the largest intersection $S$ has with any spanning tree of $G$. Similarly, define $c_{\mathrm{min}}$ to be the size of the smallest such intersection. Then,

[TABLE]

The appearance of $F_{i\leftrightarrow j}$ in this result stems from this lemma, that we rely on here and prove in Appendix D.

Lemma 2** (“Cycle flip only”).**

For any spanning trees $T,S$ and vertices $i,j$ of $G$,

[TABLE]

We will also rely on the following lemma, which generalizes the first part of the proof of Theorem 2.

Lemma 3**.**

For any symbols $\{a_{n}\},\{b_{n}\}$,

[TABLE]

Proof.

First note we can rearrange the sum as

[TABLE]

which is illustrated in Fig. 11.

As a result, we have

[TABLE]

Comparison with equations (43) through (45) in the proof of Theorem 2, together with the triangle inequality, establishes the desired result. ∎

Now we are ready to proceed with the proof of the theorem.

Proof of Theorem 3.

Define

[TABLE]

so that we have, for all $c$

[TABLE]

By the matrix-tree theorem, the derivative we wish to bound can be written in terms of these quantities as

[TABLE]

Expanding the derivative and applying Lemma 3 yields

[TABLE]

To prove the theorem, all that remains is to demonstrate that

[TABLE]

holds for all $m$. This follows by an application of Lemma 2. So we have

[TABLE]

as desired. ∎

Theorem 3 has a number of simple corollaries.

Corollary 4**.**

If $G$ has $r$ independent cycles, then for any set $S$ of edges,

[TABLE]

Proof.

Let $m$ be the number of edges in $G$. The largest possible intersection of a spanning tree and $S$ cannot exceed $|S|$ in size, so we have $c_{\mathrm{max}}\leq|S|$. Furthermore, each spanning tree of $G$ has exactly $m-r$ edges. So the smallest possible intersection is realized if all $r$ edges a spanning tree excludes are edges in the set $S$, which means $c_{\mathrm{min}}\geq|S|-r$. Therefore, $c_{\mathrm{max}}-c_{\mathrm{min}}\leq r$, and the corollary follows from Theorem 3. ∎

Corollary 5**.**

Let $H$ be a subgraph of $G$, and write $W$ for the set of vertices of $H$ incident to an edge not in $H$. Let $S$ be the edge set of $H$. Then,

[TABLE]

Proof.

Consider a spanning tree $T$ of $G$. Viewed as a subgraph of $H$, $T$ is still at least a spanning forest (i.e. it may no longer be connected, but still has no cycles and includes every vertex of $H$), with no more than $|W|$ component trees. To see this, suppose it had $|W|+1$ component trees. In this case, one component would have to be disconnected from all the vertices in $W$ (if every component is connected to a vertex in $|W|$, there can be at most $|W|$, as components cannot be connected to each other). But in that case, $T$ (as a subgraph of $G$) was disconnected—it was never a spanning tree at all.

Let $n$ be the number of vertices in $H$. The number of edges in a spanning forest is always the number of vertices in the forest minus the number of components (trees in the forest). This means that for our graph $G$, the size of the intersection of $S$ and the edge set of $T$ is restricted to lie between $n-1=c_{\mathrm{max}}$ or $n-|W|=c_{\mathrm{min}}$. By Theorem 3, this implies the result. ∎

Appendix D Proofs of the root-swapping lemmas

In the course of proving our results above we came across ratios of products of spanning tree weights, such as

[TABLE]

which we bounded using Lemmas 1 and 2, yielding our theorems. Here, we present proofs of these key lemmas. The arguments will depend on the existence of invertible mappings between the pairs of spanning trees in the numerator to pairs of spanning trees in the denominator, which have their roots “swapped”: $(T_{i},S_{j})\to(T^{\prime}_{j},S^{\prime}_{i})$. We will construct these mappings explicitly, but first we set out some relevant notation and definitions.

First, we will find it helpful in this section to use the standard notation $s(e)$ (the source) for the vertex at the tail of a directed edge $e$ and $t(e)$ (the target) for the vertex at the head of $e$, where the arrow points. In addition, the graph formed by the removal of the edge $h$ from a graph $H$, i.e. by the deletion of $h$, will be denoted $H\setminus h$, and the graph formed by adding an edge $h$ to $H$ will be denoted $H\cup h$.

Second, we need to define a new kind of spanning tree. We have already introduced spanning trees, as well as the notion of a spanning tree $T_{i}$ rooted at vertex $i$. Recall that in a rooted spanning tree, every edge is directed towards the root $i$ (since a tree has no cycles, this direction is defined unambiguously). Generalizing from this, we define a doubly-rooted spanning tree, schematically depicted in Fig. 12(a). We start with a spanning tree $S$ and two vertices $i$ and $j$. We first note that all the edges in the rooted trees $S_{i}$ and $S_{j}$ are oriented in the same direction except for those edges along the unique path between $i$ and $j$. This inspires us to pick a vertex $k$ on this path and define a doubly-rooted spanning tree $S_{ij,k}$ with branch point $k$ to be the spanning tree $S$ with every edge directed as it is in $S_{i}$ and $S_{j}$—when those directions are the same—and otherwise directed towards $i$ if between $k$ and $i$, and towards $j$ if between $k$ and $j$. One can view a (singly) rooted spanning tree $S_{j}$ as a sort of “degenerate” doubly-rooted tree $S_{ij,i}$ with branch point $i$.

Our mappings are then built from repeated applications of the following operations on pairs of the form $(T_{b},S_{mn,b})$, where $T_{b}$ is a spanning tree rooted at some vertex $b$, and $S_{mn,b}$ is a doubly-rooted spanning tree with roots $m$ and $n$ and branch point $b$:

• •

Cycle flip. Consider the unique edge $e$ pointing out of $b$ towards $n$ in $S_{mn,b}$. Reroot the tree $T_{b}$ to the target $t(e)$ to form $T_{t(e)}$ and flip the edge $e\to e^{*}$ to form $S_{mn,t(e)}=\left(S_{mn,b}\setminus e\right)\cup e^{*}$. Output $(T_{t(e)},S_{mn,t(e)})$.

• •

Edge swap. Consider the unique edge $e$ pointing out of $b$ towards $n$ in $S_{mn,b}$. Let $f$ be first edge, along the unique directed path in $T_{b}$ from $t(e)$ to $b$, that reconnects $S_{mn,b}\setminus e$. Swap these edges to form $T^{\prime}_{s(f)}=(T_{b}\setminus f)\cup e$ and $S^{\prime}_{mn,s(f)}=(S_{mn,b}\setminus e)\cup f$. Output $(T^{\prime}_{s(f)},S^{\prime}_{mn,s(f)})$.

The output of each of these operations is another pair $(T^{\prime}_{b^{\prime}},S^{\prime}_{mn,b^{\prime}})$ consisting of a tree rooted at $b^{\prime}$ and a doubly-rooted tree with branch point $b^{\prime}$ (see Figure 12 for an illustration of this in the case of edge swap). Furthermore, no edges are reoriented in the edge swap, although edges are exchanged between $T$ and $S$. In a general cycle flip, no edges are exchanged, and the edges that are reoriented form the single cycle obtained from the union of $e$ with the unique path in $T$ from $t(e)$ to $s(e)$. In the degenerate case where the path in $T$ from $t(e)$ to $s(e)$ consists of the single edge $e^{*}$, cycle flip and edge flip are equivalent.

Notably, both of these operations are invertible, in the sense that given the output of either, and knowledge of which was applied, we can uniquely recover the original pair $(T_{b},S_{mn,b})$ from $(T^{\prime}_{b^{\prime}},S^{\prime}_{mn,b^{\prime}})$.

• •

To invert the cycle flip all we need is to identify the original edge $e$—it is the reverse of the unique edge pointing out of $b^{\prime}$ towards $m$. Note that $s(e)=b$, the original branch point of $S$ and root of $T$.

• •

To invert the edge swap, we need to identify the original edges $e$ and $f$. The unique edge pointing out of $b^{\prime}$ towards $m$ is $f$. The original $e$ is the first edge, going back along the path in $T^{\prime}_{b^{\prime}}$ from $s(f)=b^{\prime}$ to $t(f)$, that reconnects $S^{\prime}_{mn,b^{\prime}}\setminus f$.

Lemma 1** (“Tree surgery”).**

Let ${\mathcal{E}}_{mn}$ be the set of spanning trees of $G$ containing the distinguished (undirected) edge $m\leftrightarrow n$. Then for any two distinct vertices $i,j$ of $G$,

[TABLE]

where $F_{\mathrm{max}}$ is the magnitude of the cycle force that is largest in magnitude, among all those associated to cycles containing the distinguished edge $m\leftrightarrow n$ (in either direction).

Proof.

To prove this result, it is sufficient to find a bijection between terms in the numerator and those in the denominator, such that each term and its partner are equal or differ by a factor of $\exp(F_{C})$, where $F_{C}$ is the cycle force associated to a cycle $C$ that contains the distinguished edge.

Consider any term $w(T_{i})w(S_{j})$ in the numerator. We map it to a term in the denominator as follows. Starting with the pair $(T_{i},S_{j})$, viewing $S_{j}$ as a doubly-rooted tree $S_{ij,i}$, we repeatedly apply edge swap until the root of the rooted tree (equivalently the branch point of the doubly-rooted tree) equals $j$, unless the edge $f$ that would be removed from $T_{b}$ in the process is the distinguished edge ($m\to n$ or $n\to m$). In that case, apply cycle flip in that step, so that the distinguished edge is not exchanged.

It is guaranteed that this iterative process will eventually terminate, because at every step, the branch point of the doubly-rooted tree $S_{ij,b}$ moves closer to $j$, and the part of $S_{ij,b}$ rooted at $i$ grows. Eventually, the branch point hits $j$, and the edge swap and cycle flip operations cannot be applied.

At the end of this iterative process the initial pair $(T_{i},S_{j})$ has been transformed into a pair $(T^{\prime}_{j},S^{\prime}_{i})$, whose associated weight $w(T^{\prime}_{j})w(S^{\prime}_{i})$ appears in the denominator of (66). This defines a bijection between terms in the numerator and terms in the denominator. To see that the map is invertible, we note that every step along the way (an application of edge swap or cycle flip) is invertible, as we argued above. Therefore, as long as it is possible to uniquely determine which was applied at each step, the whole sequence of operations is invertible. But this is possible, because when inverting a step, we can find the edge $f$ that would have been removed from $T$ by edge swap in that step, and that determines whether or not edge swap or cycle flip was in fact applied in that step. Namely, cycle flip was applied if $f$ was the distinguished edge, and edge swap was applied otherwise.

Having found this bijection between terms, it remains only for us to ask what is the ratio of the terms $w(T_{i})w(S_{j})$ and $w(T^{\prime}_{j})w(S^{\prime}_{i})$? The operation edge swap has no effect on this product of weights, since it merely moves edges between $T$ and $S$. However, when cycle flip is applied, edges change the way they are directed, and the weight $w(T_{i})w(S_{j})$ changes by a factor of $\exp(F_{C})$ where $C$ is the (directed) cycle that gets flipped. Since we only apply cycle flip if the path in $T_{b}$ from $t(e)$ to its root contains $m\leftrightarrow n$, $C$ is always a cycle containing $m\leftrightarrow n$. Furthermore, in the iteration described above, cycle flip is applied at most once. To see this, note that the original tree $T_{i}$ contains either $m\to n$ or $n\to m$, never both. Furthermore, the edge $f$ that comes up in edge swap always points from the part of $S$ rooted at $j$ to the part rooted $i$. Thus, if cycle flip flips the distinguished edge to point the other way, it will never come up as $f$ in edge swap again, because the part of $S$ rooted at $i$ only ever grows during this algorithm.

So we have

[TABLE]

for some cycle $C$ that contains the edge $m\leftrightarrow n$, as desired.

To prove the inequality (66), we now match up terms above and below using this bijective “tree surgery”, putting them in an order such that each term above has the same position as its partner (e.g. image) below. The Lemma then follows from the inequality:

[TABLE]

∎

Lemma 2** (“Cycle flip only”).**

For any spanning trees $T,S$ and vertices $i,j$ of $G$

[TABLE]

where $F_{i\leftrightarrow j}$ is largest possible value of $\ln[w(P_{i\to j}\cup P_{j\to i})/w(P_{i\to j}^{*}\cup P_{j\to i}^{*})]$ where $P_{i\to j}$ is a (non-self-intersecting) path from $i$ to $j$ and $P_{j\to i}$ is a (non-self-intersecting) path from $j$ to $i$, and the superscript ‘$\,*$’ denotes the reverse orientation.

Proof.

As above, we consider the pair $(T_{i},S_{j})$ but this time just apply cycle flip to it repeatedly until it can no longer be applied (because the branch point of $S$ has become $j$). The effect of these steps is to “swap the roots” of the two trees $T_{i}$ and $S_{j}$, changing the directions of edges without changing the underlying (undirected) spanning trees. Along any undirected spanning tree $T$, there is a unique directed path $T_{v\to w}$ from any vertex $v$ to any other vertex $w$. “Re-rooting” a tree changes its weight as follows

[TABLE]

which implies

[TABLE]

The fraction appearing here is of the form $w(P_{i\to j}\cup P_{j\to i})/w(P_{i\to j}^{*}\cup P_{j\to i}^{*})$, as required in the statement, establishing the result. ∎

It is important to note that neither Lemma 1 nor Lemma 2 implies the other, although their proofs can be viewed as depending on a common technique.

Appendix E Saturating the inequalities

We have established a number of thermodynamic bounds on steady-state response to edge perturbations. It remains an open question whether we can saturate these inequalities. In this section, we exhibit one example where we can saturate our bounds—the case of a system whose transition graph $G$ consists of a single cycle $C$ with cycle force $F_{C}=F_{\rm max}$. While we are unable to prove this is the only way to saturate our inequalities, we do argue for its general relevance.

To keep the discussion as straightforward and precise as possible, we focus on the ratio bound in (17), as this turns out to be the simplest to investigate. We first specialize to the case where we vary the edge parameter $B_{mn}$ associated to the edge $e_{mn}$, and ask for the response of the ratio of steady-state probabilities of the adjacent states $m$ and $n$. In this case, the series of inequalities that lead to our bound can be summarized as

[TABLE]

where we use the notation

[TABLE]

The first inequality in Eq. (72) is an application of the AM-GM inequality and the second comes about from our “tree surgery” argument of Lemma 1. We address each in turn.

Let us begin with the tree surgery inequality, which comes about from analyzing the ratio

[TABLE]

The tree surgery provides an invertible mapping between the terms in the numerator and those in the denominator. For the case of a single cycle with vertices $m,n$ adjacent to the distinguished edge $e_{mn}$, we have

[TABLE]

for all $T\in{\mathcal{E}}_{mn}$ and $S\notin{\mathcal{E}}_{mn}$. Thus, every term in the numerator is proportional to a term in the denominator with the same proportionality constant:

[TABLE]

Thus, the “tree surgery” inequality is exactly satisfied in this case.

Equality in the AM-GM inequality is reached when

[TABLE]

While there are numerous choices for the rates that cause this equality to be satisfied, we will just exhibit a particular one to show that it is possible. To do so, we first make a simplifying observation: each term on both sides of the equality is a product of the weight of a spanning tree rooted at $m$ and one that is rooted at $n$. Therefore, each term has exactly the same dependence on the vertex parameters $E_{j}$, so we can cancel all the $E_{j}$ on both sides of (78). Thus, all we need to do is fix the symmetric and asymmetric edge parameters. We first fix the asymmetric edge parameters by choosing all the weight of the cycle force to be on the perturbed edge $e_{mn}$,

[TABLE]

Solving Eq. (78) for the symmetric edge parameters then leads to the relation

[TABLE]

Thus, it is possible to saturate our inequality for the response of the ratio $\ln(\pi_{m}/\pi_{n})$ to perturbations of $B_{mn}$.

This may seem like a rather special case, but we believe the situation is more general than it first appears, since it is possible for the dynamics on more complicated graphs $G$ (e.g. with multiple cycles) to effectively have this “single-cycle” behavior. To see this, note that if the rates of transitions in $G$ are very small, apart from those around a single cycle containing the perturbed edge $e_{mn}$, then the graph is effectively composed of a single cycle, for the purposes of understanding response of the states on the cycle. In addition, if we look at the response of ratios of arbitrary states on the cycle, such as $\ln(\pi_{i}/\pi_{j})$, again the dynamics can effectively reproduce the situation discussed above, where we focused on the vertices adjacent to $e_{mn}$. This is because if the rates along the unique paths from $i$ to $m$ and $j$ to $n$ on the cycle are extremely fast, the states along these paths rapidly reach a local steady-state distribution. The two paths then act as two “effective states” adjacent to the perturbed edge $e_{mn}$.

These arguments suggest that, for a general graph $G$, there are limits of the rates that give rise to response approaching arbitrarily closely the bound set by Corollary 2.

Appendix F Asymmetric edge perturbation inequality

Our asymmetric edge perturbation bound follows from a more general inequality for arbitrary perturbations of a single rate:

Proposition 1**.**

[TABLE]

Proof.

By the matrix-tree theorem, we can write

[TABLE]

where $a,b$, $c$ and $d$ are nonnegative quantities formed from sums of weights of rooted spanning trees that do not depend on $W_{ij}$. By normalization of probability $\pi_{k}\leq 1$, so we have $c\geq a$, $d\geq b$. Differentiating these expressions yields

[TABLE]

which after re-arranging gives

[TABLE]

but the value of each of the two fractions on the right hand side is not smaller than zero or greater than one. This means their difference is not greater than one in magnitude, implying the result. ∎

Corollary 6**.**

[TABLE]

Proof.

The asymmetric edge parameter $F_{ij}$ appears in two rates, $W_{ij}$ and $W_{ji}$. This implies, by the chain rule

[TABLE]

which implies, by the triangle inequality,

[TABLE]

Now applying Proposition 1 establishes the desired result. ∎

Appendix G Biochemical applications

So far, we have stated and proved equalities and inequalities about the response to perturbations of physical systems whose dynamics are well-modeled as continuous in time and Markovian over a finite state space. In this section, we describe specializations of these general results to two well-known motifs found in biochemical networks. In each case, we find an inequality relating some figure of merit to a chemical potential difference driving the network out of equilibrium (for example, $\Delta\mu=\mu_{\mathrm{ATP}}-\mu_{\mathrm{ADP}}-\mu_{\mathrm{Pi}}$ for ATP hydrolysis).

There are several ways that studying a biochemical network might lead us to consider a linear time evolution equation like (1),

[TABLE]

with $\sum_{i}W_{ij}=0$ for all $j$. First, the chemical master equation, which governs the evolution of the distribution over counts $(n_{A},n_{B},\dots)$ of chemical species $A,B,\dots$, is of this form. However, for chemical systems with many particles the number of states $N$ in such a description is enormous.

However, for some chemical reaction networks, the linear equation (1) arises as the rate equation governing the deterministic evolution of the concentrations of chemical species. As emphasized by Gunawardena Gunawardena (2012, 2014), this is a generic situation that can arise from strong time-scale separation. When the rate equation of a reaction network is of the form (1), we can equivalently view it as the master equation of a continuous-time Markov chain describing the stochastic transitions of a single molecule subject to a set of effectively monomolecular reactions Mirzaev and Gunawardena (2013); Wong et al. (2017). Whichever interpretation we take, the mathematics that arises is the same, and our results can be put to work.

G.1 Covalent modification cycle

Goldbeter and Koshland Goldbeter and Koshland (1981) studied a model of the covalent modification and demodification of a substrate by two enzymes, assuming the action of both enzymes obeys mass-action kinetics with a single intermediate complex and no product rebinding:

[TABLE]

The total substrate concentration $S_{\mathrm{tot}}=[S]+[S^{*}]+[E_{1}S]+[E_{2}S]$ is conserved in these reactions, as are the enzyme totals $E_{1,\mathrm{tot}}=[E_{1}]+[E_{1}S]$, $E_{2,\mathrm{tot}}=[E_{2}]+[E_{2}S]$. In the limit of saturating substrate $S_{\mathrm{tot}}\gg E_{1,\mathrm{tot}},E_{2,\mathrm{tot}}$, the kinetics are effectively Michaelis-Menten in form, and the steady-state ratio $[S^{*}]/[S]$ can exhibit unlimited sensitivity to changes in $E_{1,\mathrm{tot}}$ and $E_{2,\mathrm{tot}}$ Goldbeter and Koshland (1981).

Sensitivity of the steady state to changes in enzyme concentrations is only possible out of equilibrium Goldbeter and Koshland (1987). In (89), the nonequilibrium nature of the system is reflected in the combination of the irreversible product release reactions with the overall reversibility of the modification of $S$.

In the regime of low substrate $S_{\mathrm{tot}}\ll E_{1,\mathrm{tot}},E_{2,\mathrm{tot}}$, we have that $[E_{1}]\approx E_{1,\mathrm{tot}}$ and $[E_{2}]\approx E_{2,\mathrm{tot}}$, and the nonlinear mass-action dynamics implied by (89) reduce to linear kinetics, with the enzyme concentrations “absorbed” into the rate constants (see Fig. 14).

In this work, we consider the low-substrate limit, and study the relative probability $\pi_{S^{*}}/\pi_{S}$ for a particular substrate molecule to be modified. For thermodynamic consistency, all reactions must be reversible, so we must include product rebinding. We further suppose that concentrations of other participants in these reactions (e.g. ATP, ADP in the case of phosphorylation/dephosphorylation) are held at fixed values. These choices yield a system of the form we have studied in the preceding sections, with linear dynamics of the form (1), held out of equilibrium by the cycle force $F_{C}=\ln(\frac{W_{21}W_{32}W_{43}W_{14}}{W_{12}W_{23}W_{34}W_{41}})$. For a system such as this one, driven chemically, we can identify $F_{C}=\Delta\mu/k_{\mathrm{B}}T$. Our results proved above then imply a bound, in terms of $\Delta\mu$, on the sensitivity $s$ of the steady-state ratio $\pi_{S^{*}}/\pi_{S}$ to a change in $[E_{1}]$.

Perturbing $[E_{1}]$ is equivalent to perturbing two edge parameters and a vertex parameter:

[TABLE]

Now we can apply Corollaries 1 and 4 to bound the sensitivity

[TABLE]

Remarkably, the form of the bound (91) remains unchanged even if the assumption that catalysis proceeds via a single intermediate complex is completely relaxed. In particular, following Gunawardena et al. Xu and Gunawardena (2012); Dasgupta et al. (2014), we consider an arbitrary reaction network built out of a collection of any number of reactions of the following form, which include an arbitrary number of intermediates and reactions between them:

[TABLE]

A general network of this form is schematically represented in Fig. 4(b). In any such network, consider the subgraph whose vertices $V$ are all the intermediates $\{(E_{1}S)_{i}\}$ containing $E_{1}$, together with $S$ and $S^{*}$, and whose edges $\mathcal{E}$ are all the edges between the vertices $V$. Scaling $[E_{1}]$ is equivalent to decreasing all the edge parameters associated to edges in $\mathcal{E}$, and the vertex parameters associated to vertices in the set $V_{I}=V\setminus\{S,S^{*}\}$. This decomposition yields the result

[TABLE]

where the last line follows from Corollary 5 with $W=\{S,S^{*}\}$, $|W|=2$.

G.2 Biochemical sensing and the Govern–ten Wolde trade-off

Now we will show how our sensing bound (24) arises, and how it reduces to the results of Govern and ten Wolde Govern and ten Wolde (2014b) in the appropriate limits.

To arrive at (24), we first employ the approximation (23)

[TABLE]

Now we recognize the derivative in the denominator as being an instance of the kinds we have considered already. In particular, application of (20), in the form of a ratio of general observables, together with our vertex perturbation results, yields

[TABLE]

where $\pi_{Y}$ is the fraction of the total substrate bound up in complexes involving $E_{1}$. In Govern and ten Wolde (2014b), these enzymatic intermediates are neglected (e.g. equation [S22] of Govern and ten Wolde (2014b)), and to proceed we will do the same here, supposing that $\pi_{Y}$ is very small. We then get

[TABLE]

as desired. In the limit $\Delta\mu\gg k_{\mathrm{B}}T$,

[TABLE]

whereas in the limit $\Delta\mu\ll k_{\mathrm{B}}T$, this yields

[TABLE]

To make contact between our low potential limit $\Delta\mu\ll k_{\mathrm{B}}T$ (98) and the results of Govern and ten Wolde, we now review the context of their results. In that paper, the authors study the error $\delta c/c$ in an estimate of the concentration $c$ of an extracellular ligand $L$, which binds to receptor $R$, forming a complex: $R+L\leftrightarrow RL$. The complex $RL$ then plays the role of $E_{1}$ in our discussion above, participating in a covalent modification cycle. The concentration of the ligand-bound receptor (in our notation $E_{1}$, in their notation $RL$) is given by

[TABLE]

where $K$ is the dissociation constant and $R_{T}$ is the total concentration of receptors.

An estimate of $c$ can be constructed, just as we describe in the main text producing an estimate of $[E_{1}]$. Using the approximation (23) together with the equation (G.2) relates the error of these estimates,

[TABLE]

The authors of Govern and ten Wolde (2014b) compare the sensing error $(\delta c/c)^{2}$ not to the chemical potential difference $\Delta\mu$, but to a quantity $w$, which is the dissipation rate of the whole system, normalized by the sum of all the rates (forward and reverse) around the modification cycle (in the case studied by the authors, consisting of only two states, in our notation, $S$ and $S^{*}$, and no intermediates). We shall call this quantity, which is the rate of relaxation to steady state, $\mathcal{R}$. So in our notation,

[TABLE]

where $J$ is the net current around the cycle, per substrate molecule. The arguments of Govern and ten Wolde show that in the limit $\Delta\mu\ll k_{\mathrm{B}}T$,

[TABLE]

This is slightly tighter than the inequality (S26) they write in Govern and ten Wolde (2014b), but also follows from their argument. As a consequence of (100) and (102), we then get

[TABLE]

To show that this bound coincides with our result (98), we need to show that for fixed small force the smallest value achieved by $4/w$ is $64(k_{\mathrm{B}}T)^{2}/N(\Delta\mu)^{2}$. If this were not so, it would imply that either our result was weaker than that of Govern and ten Wolde, or vice versa.

To show that (98) and (103) do coincide for fixed small force, we use the inequality,

[TABLE]

which holds (for a two-state system), by the same algebra (43)–(45) that led to our symmetric edge perturbation result. See also Malaguti and ten Wolde Malaguti and ten Wolde (2019) (equations S112 to S114), who give an explicit expression for $|J|/R$.

Plugging (104) into the definition (101), we get

[TABLE]

Additionally, this inequality can be saturated in the limit $\Delta\mu\ll k_{\mathrm{B}}T$, because there is a near equilibrium regime in which $|J|/R\approx\Delta\mu/16k_{\mathrm{B}}T$. This establishes the desired equivalence between our results.

G.3 Kinetic proofreading

In our presentation and analysis here, we follow closely the papers of Murugan et al. Murugan et al. (2012, 2014). Our results generalize bounds on the discriminatory index $\nu$ found in those works.

First, we consider the single-loop, three-state network (see Fig. 15) equivalent to the system studied by Hopfield and Ninio Hopfield (1974); Ninio (1975).

A perturbation of the binding energy $\Delta$ can be decomposed as a linear combination of vertex and symmetric edge parameter perturbations. In terms of the notation we introduce in Fig. 15(b), we have

[TABLE]

Now we can apply Corollaries 1 and 2 to derive the bound

[TABLE]

In the case of the more general kinetic proofreading scheme Murugan et al. (2012, 2014) where $m$ complexes can dissociate, described in Fig. 6, perturbing $\Delta$ is equivalent to perturbing the edge and vertex parameters associated to the edges $\mathcal{E}$ and vertices $V$ on one side of the “fence”. We then have

[TABLE]

where in the last line we have applied Corollary 5.