Universal Thermodynamic Bounds on Nonequilibrium Response with Biochemical Applications

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 common biochemical motifs.


I. INTRODUCTION
One of the basic characteristics of any physical system is its response to small perturbations [1]. For instance, response is used to quantify everything from material properties-such as conductivity [2] and viscoelasticity [3]-to the sensing capability of cells [4,5] and the accuracy of biomolecular processes [6][7][8]. Near thermodynamic equilibrium, response is completely determined by the nature of spontaneous fluctuations, according to the fluctuation-dissipation theorem (FDT) [2]. 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 [1]. 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 [1] 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 [9][10][11][12][13]. While these predictions offer fundamental theoretical insight, the necessary correlations are often prohibitively difficult to measure except in simple single-particle systems [14][15][16]. In certain special cases, such as under stalling conditions, the FDT holds unmodified [17]. 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 [18,19], a useful prediction that is utilized to measure dissipation and efficiency in molecular motors [20]. Alternatively, violations of the FDT can be used to fit model parameters, as suggested for models of biomolecular processes [21,22]. More often FDT violations are framed in terms of system-specific "effective temperatures" [23][24][25], 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 [26,27], we show here that generic nonequilibrium steadystate 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; …; 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 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 ¼ − P j≠i W ji are fixed by the conservation of probability. Time reversibility of the underlying microscopic dynamics implies that W ij ≠ 0 only if W ji ≠ 0 [28]. We additionally suppose that, for any two states, there is some sequence of allowed transitions (W ij ≠ 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 π that satisfies P N j¼1 W ij π j ¼ 0. This distribution π, 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 Eq. (1) playing out on a transition graph-a weighted directed graph G, as in Fig. 1(a), where the vertices fig represent the states and directed edges fe mn g represent possible transitions, weighted by the rates W mn . Note that, by assumption, every edge in G has a reverse, so we 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: 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 [29,30]: These cycle forces are linear combinations of thermodynamic affinities multiplied by their conjugate distancesfor 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 correspondence means that the cycle forces depend on macroscopically tunable parameters-such as the 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 [31] characteristic of thermodynamic equilibrium.

III. STATIC RESPONSE TO PERTURBATIONS
Now, suppose the transition rates W ij ðλÞ depend on a control parameter λ, 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 hQi π ¼ P j Q j π j respond to small changes in λ: At thermal equilibrium, the steady state π eq i ∝ e −βϵ i ðλÞ depends only on the underlying (free) energy landscape ϵ i ðλÞ, irrespective of the precise form of the transition rates, where β ¼ 1=k B T with k B Boltzmann's constant and T the temperature. This simplifying fact immediately implies the static FDT, which equates the static response to an equilibrium correlation function: FIG. 1. Transition graphs and cycles. (a) Representative transition graph for a four-state system, which acts as a recurring illustrative example. (b) Cycles around which the cycle forces F C drive the system out of equilibrium. ∂ λ hQi eq ¼ βC eq ðQ; VÞ; ð4Þ where C eq ðQ; VÞ ¼ hQVi eq − hQi eq hVi eq and the subscript "eq" emphasizes that averages are taken with respect to the equilibrium distribution [2]. Here, V ¼ −∂ λ ϵ is known as the coordinate conjugate to λ and represents the displacement induced by λ-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 [9,10], but the relevant "conjugate coordinate" ∂ λ ln π requires a knowledge of the parameter dependence of the full nonequilibrium steady-state distribution, which can be challenging to calculate or measure. Still, this response formula is given thermodynamic meaning by relating the nonequilibrium conjugate coordinate to the stochastic entropy production rate [13] as well as to the time-reversal symmetry properties of the path action [11].

IV. PARAMETERIZING PERTURBATIONS
Here, we turn our attention to the variations of the steady-state distribution with the transition rates, ∂π i =∂W kl , constraining them in terms of π 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 are able to provide clear and measurable thermodynamic constraints on the static response.
To classify perturbations, it proves fruitful to parameterize the rate matrix as 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 nonuniquely. To see this fact, consider the following program for identifying such a parameterization: Choose the vertex parameters fE 1 ; …; E N g arbitrarily and then set The nonuniqueness 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 , are equivalent to scaling all of the rates out of state j: Similarly, derivatives with respect to B ij and F ij multiplicatively scale transition rates associated with a single edge: These facts are illustrated in Fig. 2. We note that the righthand sides of Eqs. (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 Eq. (6) does 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 ¼ P e ij ∈C F ij . Furthermore, if all the F ij ¼ 0, the steady-state distribution has the Boltzmann form π i ∝ 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 FIG. 2. Parameterizing perturbations. Any perturbation of rates can be decomposed into the variation of some combination of (a) vertex parameters, (b) edge parameters, and (c) asymmetric edge parameters. Affected rates are highlighted.
The freedom in the rate parameterization makes this decomposition nonunique, and the tightness of our inequalities depends on the specific decomposition. For example, one choice could force ∂ λ 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 the Appendix A, which gives an exact algebraic expression for the steady-state probabilities π i in terms of the rates W ij [29,32]. 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 Appendixes.

V. VERTEX PERTURBATIONS
Our first main result is the exact expression for the response to a vertex perturbation (Appendix B) 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 ≠ j, 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 a farfrom-equilibrium response has an equilibriumlike 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 ϵ i ðλÞ and that the rates depend on the energy as W ij ¼ ω ij e βϵ j , with arbitrary energy-independent ω ij . A comparison with Eq. (6) shows that variations in the energy ϵ i in this case can be parameterized as vertex parameters E i . Then, Eq. (13) implies that arbitrarily far from equilibrium the response maintains the equilibriumlike form of the FDT, with the response proportional to the nonequilibrium steady-state correlation with the coordinate conjugate to the energy V ¼ −∂ λ ϵ [cf. Eqs. (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 modifies 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, a response to edge perturbations is constrained by thermodynamic affinities through the cycle forces. For proofs of the results in this section, see Appendix C.

A. 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: where 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 equilibriumthen F 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 π i ð1 − π i Þ as the variance of the occupation fluctuations of state i; thus, we see a manifestation of how thermodynamics shapes the interplay between the response and fluctuations. These inequalities, applying to all discrete stochastic dynamics, significantly generalize a bound for two-state systems derived by Hartich, Barato, and Seifert in a model of nonequilibrium sensing [33]. The conditions for equality in Eqs. (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 Eq. (17): a single cycle with strong timescale separation so that the system effectively has only 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.

B. Multiple edges
The response to a perturbation of multiple edge parameters can be bounded by combining Eq. (17) with the triangle inequality. For example, for any set S of jSj edges, 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 π 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 Eq. (19). Suppose we vary the edge parameters associated to the edges S of a subgraph H ≠ 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, where F i↔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 a self-intersection).
Whenever there is only one path through state space between i and j, and in all cases at thermodynamic equilibrium, F i↔j ¼ 0.
We finally note that our results (16), (17), and (20) admit generalization to the response of a ratio of positive observables hAi=hBi. In this general case, the bounds remain unchanged, except that F i↔j is replaced by its maximum over all pairs of vertices i, j.

VII. ASYMMETRIC EDGE PERTURBATIONS
Last, the MTT allows us to bound the response to asymmetric edge perturbations as (Appendix F) which is related to, but distinct from, inequalities established in Ref. [34]. This result is a consequence of an identical inequality that holds for general rate perturbations ∂π i =∂ 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 introduce, 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.
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 FIG. 3. Thermodynamics and topology bound response. (a) Transition graph for our representative four-state system with a single perturbed edge connecting states 1 and 2, e 12 , highlighted in blue. (b) Only cycle forces F C around cycles that contain the perturbed edge e 12 , highlighted in blue, constrain the static response. (c) The maximum response max j j∂π j =∂B 12 j to the perturbation of the edge parameter B 12 as a function of the maximum cycle force around cycles containing e 12 for 15 000 randomly sampled rate matrices (gray dots). All samples fall below the predicted bound ð1=4Þ tanhðF max =4Þ (red line). 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 Δμ ¼ μ ATP − μ ADP − μ 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 [35].
In the other limit-that of a single substrate moleculeour results (17) limit the sensitivity of the ratio π S Ã =π S for a particular substrate molecule to changes in the enzyme concentration (Appendix G): where F max ¼ Δμ=k 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 [37,38]. 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 multiedge bound (20) then applies with F S↔S Ã ¼ Δμ=k 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 (Δμ ¼ 0), it is clear that 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 [36,40,41]. Our results provide a general quantification of this requirement.

B. Biochemical sensing
The covalent modification cycle, discussed in the previous section, is an integral component of numerous biochemical models for cellular sensing [5,39,[42][43][44]. 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 number 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 μ s ¼ Nπ S Ã and variance σ 2 s ¼ Nπ S Ã ð1 − π 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 ares molecules of S Ã . Supposing ½E 2 and all rate constants assume known, fixed values, one can produce an estimate c ½E 1 of ½E 1 by choosing c ½E 1 to be the value of ½E 1 that gives μ s ð c ½E 1 Þ ¼s. The variance of the estimate c ½E 1 so constructed is often well approximated, when the noise is small, by [5,39] 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 π S Ã ð1 − π S Ã Þ ≤ 1=4 then leads to the following bound on the relative error: This result interpolates between bounds on error established by Govern and ten Wolde [5,39] 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.

C. Kinetic proofreading
As a third application, we turn to the effectiveness of kinetic proofreading [45,46]. 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ð−ΔÞ, where k B TΔ 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 the 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, Huse, and Leibler [47]: At equilibrium, ν ¼ 1. The simplest nonequilibrium scheme to increase discrimination is the single-cycle network illustrated in Fig. 5(a). Note that we suppose the binding energy Δ appears exclusively in the unbinding rates. Hopfield observes that, in a certain nonequilibrium limit of the rates, ν → 2 [45].
Our results lead to a constraint on ν that interpolates between the equilibrium case and this limit.
In the single-cycle network, the variation of the binding energy Δ is equivalent to the variation of two vertex parameters (that of ES and ES Ã ) and a single-edge parameter (ES ↔ ES Ã ), leading to the inequality (Appendix G) where F max ¼ Δμ=k B T is the chemical driving around the cycle. This bound, which can be saturated, reduces correctly to ν ¼ 1 at equilibrium and is consistent with ν → 2 in the limit of strong driving Δμ → ∞.
We can also bound ν in the case of a more general kinetic proofreading scheme [7,47] in which there are m complexes that can dissociate. Each of the dissociation transitions can be thought of as crossing a "discriminatory fence" [47], its rate depending on the binding energy Δ, as in Fig. 6. We suppose that the dissociation transitions are the only ones that depend on Δ. We make no assumptions about the structure of the transition graph on either side of the fence. In such a network, perturbing Δ 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) where F E↔ES is the maximum entropy produced to go from E to ES and back again. Notably, we recover the result jν − 1j ≤ m − 1 of [47].

IX. CONCLUSION
In this work, we have developed a series of universal bounds on a 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 the 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 or 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 the response, because the response to the scaling of any one rate is bounded by 1=4. At the same time, our multiedge 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 [47] and biochemical copy processes [48].
Our results point to numerous other extensions, including bounds on the response of currents with implications for the Green-Kubo and Einstein relations [49][50][51]. 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 strikes us as a promising direction for future work.

ACKNOWLEDGMENTS
We thank Alexandre Solon and Matteo Polettini for very useful discussions. J. A. O. thanks 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 ¼ fi → j; k → l; …g, we define the weight wðSÞ to be the product of the weights of the edges: 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 toward" r, which we then call the "root." The resulting directed graph, which we write as T r , is a rooted spanning tree of G. The steady-state distribution π is given explicitly by the MTT [29,[52][53][54][55][56] in terms of weights of rooted spanning trees of G.
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 where N ¼ P N k¼1 P 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 [52] and has been rediscovered repeatedly in different literature; see, e.g., Refs. [29,[53][54][55][56] for further discussion.
The MTT offers a graphical representation of the steadystate distribution that provides a convenient method for organizing the structure of the solution. We illustrate this result in Fig. 8.

APPENDIX B: VERTEX PERTURBATIONS
Proof.-The matrix-tree theorem implies that π i can be expressed as the ratio of sums of weights of rooted spanning trees. So, to evaluate ∂π i =∂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 steadystate distribution in a convenient manner as illustrated in Fig. 9.
Thus, for i ¼ k, the matrix-tree theorem implies that where a ¼ 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 ≠ k, the MTT yields by a similar argument with The theorem now follows by differentiating these expressions. When i ≠ k, and similarly for i ¼ k.
Proof.-First, note that Now we apply Theorem 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)   where F max is the magnitude of the cycle force that is the largest in magnitude, among all those associated to cycles containing the distinguished edge m ↔ n (in either direction). Our proof relies on the following technical lemma, which we prove in Appendix D.
Lemma 1: "Tree surgery."-Let E mn be the set of spanning trees of G containing the distinguished (undirected) edge m ↔ n. Then, for any two distinct vertices i, j of G, 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 E mn .
Then, by the matrix-tree theorem, we can write where 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.

Differentiating yields
Now note that by the inequality of arithmetic and geometric means (AM-GM) the denominator is bounded as To complete the proof, we need to bound the ratio b 0 a 1 =a 0 b 1 by expðF max Þ. To establish this bound, we match up terms above and below, writing the fraction as The desired result is now a consequence of the inequality P n i¼1 x i P n i¼1 y i followed by Lemma 1. ▪ From Theorem 2, we readily obtain our bounds on steady-state response. For Eq. (17), Proof.-Choose the observables in Theorem 2 to be A l ¼ δ il and B l ¼ δ kl , where δ ij is the Kronecker delta. ▪ We also have Corollary 3.-Let π X ¼ P k∈X π k be the total probability of a set of states X. Then, Proof.-Choose the observables in Theorem 2 to be A i ¼ δ i ðXÞ and B i ¼ 1 − δ i ðXÞ, where the indicator δ i ðXÞ ¼ 1 if i ∈ X and δ i ðXÞ ¼ 0 otherwise. Note that we then have hAi ¼ π X and hBi ¼ 1 − π X .
▪ If X ¼ fig consists of only a single state, we recover the bound (16).

Multiple edges
In this section, we derive our inequality for the response to perturbations by multiple edge parameters Eq. (20). The proof proceeds in two steps. We first prove a bound on an arbitrary set of edges S from which Eq. (20) and other results are ready corollaries.
Here, the magnitude of response is bounded by a different function F i↔j of cycle forces. The quantity F i↔j is defined for any graph G and vertices i and j to be where P i→j is a (non-self-intersecting) path from i to j, P j→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 max to be the size of the largest intersection S has with any spanning tree of G. Similarly, define c min to be the size of the smallest such intersection. Then, The appearance of F i↔j in this result stems from Lemma 2, 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, We also rely on the following lemma, which generalizes the first part of the proof of Theorem 2.
Lemma 3.-For any symbols fa n g; fb n g, ð P j n¼i na n Þ P j n¼i b n − ð P j n¼i nb n Þ P j n¼i a n P j n¼i a n P j n¼i b n ≤ X j m¼iþ1 tanh 1 4 ln P j n¼m a n P m−1 n¼i b n P j n¼m b n P m−1 n¼i a n : Proof.-First, note that we can rearrange the sum as X j n¼i na n ¼ X j n¼i X n m¼1 a n ¼ i X j n¼i a n þ X j m¼iþ1 X j n¼m a n ; ðC18Þ which is illustrated in Fig. 11. As a result, we have ð P j n¼i na n Þ P j n¼i b n − ð P j n¼i nb n Þ P j n¼i a n P j n¼i a n P j n¼i b n ¼ X j m¼iþ1 ð P j n¼m a n Þ P j n¼i b n − ð P j n¼m b n Þ P j n¼i a n P j n¼i a n P j n¼i b n ¼ X j m¼iþ1 ð P j n¼m a n Þ P m−1 n¼i b n − ð P j n¼m b n Þ P m−1 n¼i a n ð P m−1 n¼i a n þ P j n¼m a n Þð so that we have, for all c, By the matrix-tree theorem, the derivative we wish to bound can be written in terms of these quantities as Expanding the derivative and applying Lemma 3 yields P c max n¼m a n P m−1 n¼c min b n P c max n¼m b n P m−1 n¼c min a n : ðC23Þ To prove the theorem, all that remains is to demonstrate that P c max n¼m a n P m−1 n¼c min b n P c max n¼m b n P m−1 n¼c min a n holds for all m. This result follows by an application of Lemma 2. So we have 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, Proof.-Let m be the number of edges in G. The largest possible intersection of a spanning tree and S cannot exceed jSj in size, so we have c max ≤ jSj. 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 min ≥ jSj − r. Therefore, c max − c min ≤ 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, 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 jWj component trees. To see this fact, suppose it had jWjþ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 jWj, there can be at most jWj, as components cannot be connected to each other). But, in that case, T (as a subgraph of G) is disconnected. This is a contradiction, because as a spanning tree, T must be connected.
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 fact 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 max or n − jWj ¼ c 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 come across ratios of products of spanning tree weights, such as P T∈E mn which we bound using Lemmas 1 and 2, yielding our theorems. Here, we present proofs of these key lemmas.
The arguments 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 Þ → ðT 0 j ; S 0 i Þ. We construct these mappings explicitly, but first we set out some relevant notation and definitions.
First, we 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, is denoted Hnh, and the graph formed by adding an edge h to H is denoted H ∪ 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 toward the root i (since a tree has no cycles, this direction is defined unambiguously). Generalizing from this fact, 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 result 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 toward i if between k and i and toward 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. Lemma 1: "Tree surgery."-Let E mn be the set of spanning trees of G containing the distinguished (undirected) edge m ↔ n. Then, for any two distinct vertices i, j of G, where F max is the magnitude of the cycle force that is the largest in magnitude, among all those associated to cycles containing the distinguished edge m ↔ 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 → n or n → m). In that case, apply cycle flip in that step, so that the distinguished edge is not exchanged (Fig. 13).
It is guaranteed that this iterative process eventually terminates, 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 0 j ; S 0 i Þ, whose associated weight wðT 0 j ÞwðS 0 i Þ appears in the denominator of Eq. (D2). This transformation 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 argue above. Therefore, as long as it is possible to uniquely determine which is applied at each step, the whole sequence of operations is invertible. But this identification 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 is, in fact, applied in that step. Namely, cycle flip is applied if f is the distinguished edge, and edge swap is 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 0 j ÞwðS 0 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 apply cycle flip only if the path in T b from tðeÞ to its root contains m ↔ n, C is always a cycle containing m ↔ 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 → n or n → 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
for some cycle C that contains the edge m ↔ n, as desired.
To prove the inequality (D2), 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: ▪ Lemma 2: "Cycle flip only."-For any spanning trees T, S and vertices i, j of G, where F i↔j is the largest possible value of ln½wðP i→j ∪ P j→i Þ=wðP Ã i→j ∪ P Ã j→i Þ, where P i→j is a (non-self-intersecting) path from i to j and P j→i is a (nonself-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→w from any vertex v to any other vertex w. "Rerooting" a tree changes its weight as follows: which implies The fraction appearing here is of the form wðP i→j ∪ P j→i Þ= wðP Ã i→j ∪ P Ã j→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 a 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 max . While we are unable to prove this example 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 Eq. (17), as this bound 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 where we use the notation The first inequality in Eq. (E1) 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 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 for all T ∈ E mn and S ∉ E mn . Thus, every term in the numerator is proportional to a term in the denominator with the same proportionality constant: Thus, the tree surgery inequality is exactly satisfied in this case. Equality in the AM-GM inequality is reached when While there are numerous choices for the rates that cause this equality to be satisfied, we 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 Eq. (E7). 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 : Solving Eq. (E7) for the symmetric edge parameters then leads to the relation Thus, it is possible to saturate our inequality for the response of the ratio lnðπ m =π n Þ to perturbations of B mn . This case may seem rather special, 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 "singlecycle" behavior. To see this possibility, 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ðπ i =π j Þ, again the dynamics can effectively reproduce the situation discussed above, where we focus on the vertices adjacent to e mn . This result is possible 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 a 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.
Proof.-By the matrix-tree theorem, we can write where a, b, c, and d are non-negative quantities formed from sums of weights of rooted spanning trees that do not depend on W ij . By normalization of probability π k ≤ 1, so we have c ≥ a, d ≥ b. Differentiating these expressions yields which after rearranging gives but the value of each of the two fractions on the right-hand side is not smaller than zero or greater than one, which means their difference is not greater than one in magnitude, implying the result.
Proof.-The asymmetric edge parameter F ij appears in two rates, W ij and W ji . This fact implies, by the chain rule, which implies, by the triangle inequality, 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 wellknown 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, Δμ ¼ μ ATP − μ ADP − μ Pi for ATP hydrolysis).
There are several ways that studying a biochemical network might lead us to consider a linear time evolution equation like Eq. (1), with P 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 ; …Þ of chemical species A; B; …, is of this form. 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 [57,58], this generic situation can arise from strong timescale 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 [56,59]. Whichever interpretation we take, the mathematics that arises is the same, and our results can be put to work.

Covalent modification cycle
Goldbeter and Koshland [35] study 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: The total substrate concentration S tot ¼ ½S þ ½S Ã þ ½E 1 S þ ½E 2 S is conserved in these reactions, as are the enzyme totals E 1;tot ¼½E 1 þ½E 1 S and E 2;tot ¼ ½E 2 þ ½E 2 S. In the limit of saturating substrate S tot ≫ E 1;tot ; E 2;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;tot and E 2;tot [35]. Sensitivity of the steady state to changes in enzyme concentrations is possible only out of equilibrium [40]. In Eq. (G2), 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 tot ≪ E 1;tot ; E 2;tot , we have that ½E 1 ≈ E 1;tot and ½E 2 ≈ E 2;tot , and the nonlinear mass-action dynamics implied by Eq. (G2) 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 π S Ã =π 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 and ADP in the case of phosphorylation or dephosphorylation) are held at fixed values. These choices yield a system of the form we study in the preceding sections, with linear dynamics of the form (1), held out of equilibrium by the cycle force F C ¼lnð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 ¼ Δμ=k B T. Our results proved above then imply a bound, in terms of Δμ, on the sensitivity s of the steadystate ratio π S Ã =π S to a change in ½E 1 .
Perturbing ½E 1 is equivalent to perturbing two edge parameters and a vertex parameter: Now we can apply Corollaries 1 and 4 to bound the sensitivity Remarkably, the form of the bound (G4) remains unchanged even if the assumption that catalysis proceeds via a single intermediate complex is completely relaxed. In particular, following Gunawardena et al. [37,38], 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: 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 fðE 1 SÞ i g containing E 1 , together with S and S Ã , and whose edges E are all the edges between the vertices V. Scaling ½E 1 is equivalent to decreasing all the edge parameters associated to edges in E and the vertex parameters associated to vertices in the set V I ¼ VnfS; S Ã g. This decomposition yields the result  14. (a) The transition graph G arising from the Goldbeter-Koshland model in the low-substrate limit, with product rebinding. Two transition rates (red) depend on the (assumed constant) free enzyme concentration ½E 1 that we vary. Scaling ½E 1 is equivalent to a perturbation of two edge parameters and one vertex parameter (blue). (b) State numbers and rate labels we use in this subsection. Key equivalences are 1 ¼ S, 3 ¼ S Ã , where the last line follows from Corollary 5 with W ¼ fS; S Ã g, jWj ¼ 2.

Biochemical sensing and the Govern-ten Wolde trade-off
Now we show how our sensing bound (24) arises and how it reduces to the results of Govern and ten Wolde [39] in the appropriate limits.
To arrive at Eq. (24), we first employ the approximation (23): Now we recognize the derivative in the denominator as being an instance of the kinds we have considered already. In particular, the application of Eq. (20), in the form of a ratio of general observables, together with our vertex perturbation results, yields ½E 1 ∂π S Ã ∂½E 1 þ π S Ã π Y ≤ π S Ã ð1 − π S Ã Þ tanhðΔμ=4k B TÞ; ðG8Þ where π Y is the fraction of the total substrate bound up in complexes involving E 1 . In Ref. [39], these enzymatic intermediates are neglected (e.g., Eq. [S22] of Ref. [39]), and to proceed we do the same here, supposing that π Y is very small. We then get σ e ½E To make contact between our low potential limit Δμ ≪ k B T (G11) and the results of Govern and ten Wolde, we now review the context of their results. In that paper, the authors study the error δc=c in an estimate of the concentration c of an extracellular ligand L, which binds to receptor R, forming a complex: R þ L ↔ 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 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 Eq. (G12) relates the error of these estimates: δc The authors of Ref. [39] compare the sensing error ðδc=cÞ 2 not to the chemical potential difference Δμ, 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 the steady state, R. So, in our notation, where J is the net current around the cycle, per substrate molecule. The arguments of Govern and ten Wolde show that, in the limit Δμ ≪ k B T, δc c This result is slightly tighter than the inequality (S26) they write in Ref. [39] but also follows from their argument. As a consequence of Eqs. (G13) and (G15), we then get To show that this bound coincides with our result (G11), we need to show that for a fixed small force the smallest value achieved by 4=w is 64ðk B TÞ 2 =NðΔμÞ 2 . If this fact 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 Eqs. (G11) and (G16) do coincide for a fixed small force, we use the inequality which holds (for a two-state system), by the same algebra (C6)-(C8) that leads to our symmetric edge perturbation result. See also Malaguti and ten Wolde [60] [Eqs. (S112)-(S114)], who give an explicit expression for jJj=R. Plugging Eq. (G17) into the definition (G14), we get Additionally, this inequality can be saturated in the limit Δμ ≪ k B T, because there is a near-equilibrium regime in which jJj=R ≈ Δμ=16k B T. This result establishes the desired equivalence between our results.

Kinetic proofreading
In our presentation and analysis here, we follow closely the papers of Murugan, Huse, and Leibler [7,47]. Our results generalize bounds on the discriminatory index ν 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 [45,46]. A perturbation of the binding energy Δ can be decomposed as a linear combination of vertex and symmetric edge parameter perturbations. In terms of the notation we introduce in Fig. 15 23 : Now we can apply Corollaries 1 and 2 to derive the bound In the case of the more general kinetic proofreading scheme [7,47], where m complexes can dissociate, described in Fig. 6, perturbing Δ is equivalent to perturbing the edge and vertex parameters associated to the edges E and vertices V on one side of the fence. We then have  15. (a) The transition graph G for the single-loop kinetic proofreading mechanism of Hopfield. The dissociation rates of the complexes ES Ã and ES are the only rates that depend on the binding energy Δ. Perturbing Δ is equivalent to vertex and edge perturbations (blue). (b) State numbers and rate labels we use in this subsection. Key equivalences are 1 ¼ E, 3 ¼ ES, W 12 ¼ k 1 e Δ , and W 13 ¼ k 2 e Δ .