An energy-speed-accuracy relation in complex networks for biological discrimination
AAn energy-speed-accuracy relation in complex networks for biological discrimination
Felix Wong, , Ariel Amir, and Jeremy Gunawardena , ∗ School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA Department of Systems Biology, Harvard Medical School, Boston, MA 02115, USA
Discriminating between correct and incorrect substrates is a core process in biology but how is energy ap-portioned between the conflicting demands of accuracy ( µ ), speed ( σ ) and total entropy production rate ( P )?Previous studies have focussed on biochemical networks with simple structure or relied on simplifying kineticassumptions. Here, we use the linear framework for timescale separation to analytically examine steady-stateprobabilities away from thermodynamic equilibrium for networks of arbitrary complexity. We also introduce amethod of scaling parameters that is inspired by Hopfield’s treatment of kinetic proofreading. Scaling allowsasymptotic exploration of high-dimensional parameter spaces. We identify in this way a broad class of com-plex networks and scalings for which the quantity σ ln( µ ) /P remains asymptotically finite whenever accuracyimproves from equilibrium, so that µ eq /µ → . Scalings exist, however, even for Hopfield’s original network,for which σ ln( µ ) /P is asymptotically infinite, illustrating the parametric complexity. Outside the asymptoticregime, numerical calculations suggest that, under more restrictive parametric assumptions, networks satisfy thebound, σ ln( µ/µ eq ) /P < , and we discuss the biological implications for discrimination by ribosomes andDNA polymerase. The methods introduced here may be more broadly useful for analysing complex networksthat implement other forms of cellular information processing. I. INTRODUCTION
In cellular information processing, a biochemical mecha-nism is coupled to an environment of signals and substratesand carries out tasks such as detection [1–5], amplification[6–8], discrimination [9–24], adaptation [25], searching [26]and learning [27–30]. As Hopfield pointed out in his seminalwork on discrimination [9], systems operating at thermody-namic equilibrium have limited information processing capa-bility and energy must be expended to do better [8].We focus here on the widely-studied task of discriminationbetween correct and incorrect substrates, an essential featureof many core biological processes. The accuracy of discrim-ination may have to be traded off against speed while energyremains a limiting resource [25, 31]. How can energy be ap-portioned between such desirable properties as accuracy andspeed and the inevitable dissipation of heat to the environ-ment? Quantitative insights into this question can help us dis-till the principles underlying cellular information processingdespite the pervasive complexity of the underlying molecularmechanisms.Previous studies have usually focussed on particular sys-tems, such as Hopfield’s original proofreading mechanism[9, 19, 20], McKeithan’s T-cell receptor mechanism [14, 18],minimal feedback mechanisms [25] or ladder mechanisms[15, 16, 21]. Substantial parametric complexity has beenfound in these individual systems, with different relationshipsbetween energy, speed and accuracy in different regions oftheir parameter spaces. Murugan et al. analysed general sys-tems using simplifying assumptions about where energy is ex-pended and showed how discriminatory regimes also dependon the topology of the mechanism [15, 16]. One of the chal-lenges in dealing with general systems away from thermody-namic equilibrium is that of steady-state “history dependence” ∗ Corresponding author: [email protected] (see the Discussion), which gives rise to substantial algebraiccomplexity in steady-state calculations [8, 32] and makes itdifficult to identify any universal behaviour.We address these issues here in two ways. First, we use thepreviously developed “linear framework” for timescale sepa-ration, which allows a graph-based treatment of Markov pro-cesses of arbitrary structure away from thermodynamic equi-librium ( § II) in which steady-state probabilities can be analyt-ically calculated ( § III). The linear framework was originallydeveloped to analyse cellular systems like post-translationalmodification and we find it helpful because it provides ageneral foundation for many types of timescale separationcalculations in molecular biology. Second, we introduce away of exploring parameter space by scaling the parameters.This idea is inspired by Hopfield’s original analysis of ki-netic proofreading, which we revisit here to point out certainsubtleties that are not always appreciated ( § IV). The scalingmethod allows us to calculate the asymptotic behaviour ofsteady-state properties of general systems, despite the com-plexities arising from high-dimensional parameter spaces andhistory dependence. In this way, we are able to exhibit a uni-versal asymptotic relationship between energy, speed and ac-curacy for a broad class of general systems, without simpli-fying assumptions as to where energy is expended ( § VI). Wefurther explore whether this asymptotic relationship also hassignificance for finite parameter values and for actual biolog-ical discrimination mechanisms ( § VIII).
II. THE LINEAR FRAMEWORK FOR TIMESCALESEPARATION
In a timescale separation, the mathematical description ofa system is simplified by assuming that it is operating suf-ficiently fast with respect to its environment that it can betaken to have reached steady state. This is then used to elim-inate the “fast” components within the system in terms ofthe rate constants and the “slow” components in the envi- a r X i v : . [ q - b i o . M N ] O c t ronment. The linear framework is a graph-based methodol-ogy for systematically undertaking such eliminations [33]. Itprovides a foundation for the classical timescale separationsin biochemistry and molecular biology and has been usedto analyse protein post-translational modification [34], cova-lent modification switches [35] and eukaryotic gene regulation[8, 32]. Some aspects of the framework have appeared pre-viously in the biophysics literature, in the work of Hill [36]and Schnakenberg [37], as well as in several other literatures,but the scope and implications for biology have only recentlybeen clarified [33]. Since its use distinguishes this paper fromothers, we provide here a brief overview. For details and his-torical connections, see refs. [33] and [38]; for a review, seeref. [39].In the linear framework, a system is represented by a di-rected graph, G , with labelled edges and no self-loops (Fig.1(a)), hereafter a “graph”. The vertices, , · · · , n , can be in-terpreted as the “fast” components and a labelled edge, i a → j ,as an abstract reaction whose effective rate constant is the la-bel a . Labels can be complicated expressions involving rateconstants and “slow” components. In this way, certain nonlin-ear reactions can be rewritten as linear reactions (edges) withcomplicated labels. Provided “fast” and “slow” componentscan be uncoupled in this way, the system dynamics can berewritten as if the edges are reactions under mass action ki-netics. This yields a linear dynamics, du/dt = L ( G ) u , inwhich u ∈ ( R ≥ ) n is the vector of component concentrationsand L ( G ) is the Laplacian matrix of G . For instance, for thesubgraph G C in Fig. 1(a), L ( G C ) = − ( k (cid:48) C + l (cid:48) C ) k C l C + Wk (cid:48) C − ( k C + m (cid:48) ) ml (cid:48) C m (cid:48) − ( l C + W + m ) . Since the total concentration is conserved, there is a conversa-tion law, (cid:80) i u i ( t ) = u tot .In a microscopic interpretation, vertices can be microstatesand edges can be transitions, with the labels as rates. A typicalMarkov process, M , gives rise to a graph, G M , for whichLaplacian dynamics, dp/dt = L ( G M ) p with (cid:80) i p i ( t ) = 1 ,gives the master equation of the Markov process. Here, p i ( t ) is the probability that M is in microstate i at time t .The language of graph theory accommodates both macro-scopic interpretations of molecular populations in biochem-istry [34, 35] and microscopic interpretations of Markov pro-cesses on single molecules [8, 32]. While the linear Lapla-cian dynamics is universal, the treatment of the “slow” com-ponents in the labels, which carry the nonlinearity, dependson the application and the questions being asked. We will as-sume below that the concentrations of “slow” components areconstant.Elimination of the “fast components” arises from two re-sults. First, if the graph is strongly connected, so that any twovertices can be joined by a path of edges in the same direc-tion, then there is a unique steady state up to a scalar multiple.Second, a representative steady state, ρ ( G ) , can be calculatedin terms of the labels by the Matrix-Tree theorem (MTT): if Θ i ( G ) denotes the set of spanning trees rooted at i (Fig. 1(c)),then ρ i ( G ) is the sum of the product of the labels on the edges of each tree, ρ i ( G ) = (cid:88) T ∈ Θ i ( G ) (cid:89) j a → k ∈ T a . (1)Results equivalent to the MTT were proved independently byHill [36], Schnakenberg [37] and many others in different sci-entific contexts; for a historical overview, see [38]. Schnaken-berg’s work has come back into view, as in the work of Mu-rugan et al. [16], but the problem of history-dependent alge-braic complexity (Discussion) may have discouraged attemptsto exploit equation (1) as we do here.If a system reaches a steady state, u ∗ , then u ∗ = λρ ( G ) ,where λ is the only remaining degree of freedom at steadystate. The “fast” components can then be eliminated alongwith λ using the conservation law, u ∗ i = (cid:18) ρ i ( G ) ρ ( G ) + · · · + ρ n ( G ) (cid:19) u tot . (2)If the steady state is one of thermodynamic equilibrium,so that detailed balance is satisfied, then the framework givesthe same result as equilibrium statistical mechanics, with thedenominator in equation (2) being the partition function (upto a constant factor). But equation (2) is also valid awayfrom equilibrium, so the framework offers a form of non-equilibrium statistical mechanics.In contrast to eigenvalues or determinants, the MTT givesthe steady state analytically in terms of the labels (equation(1)). This makes it feasible to undertake a mathematical anal-ysis, without relying on numerical simulation, whose scopeis necessarily more restricted. Substantial algebraic complex-ity can arise in equation (1) through history dependence awayfrom equilibrium (Discussion) but, as we show here, with theappropriate mathematical language, it is possible to draw rig-orous conclusions about structurally complex systems awayfrom thermodynamic equilibrium. III. STEADY STATES OF A BUTTERFLY GRAPH
Discrimination typically requires a mechanism for choos-ing a correct substrate from among a pool of available sub-strates, as in DNA replication, in which DNA polymerasemust choose at each step one correct deoxynucleoside triphos-phate from among the four available (dATP, dGTP, dCTP,dTTP). We follow Hopfield in assuming a single correct sub-strate, C , and a single incorrect substrate, D , and describe thismechanism by a graph G (e.g. Fig. 1(a)) whose vertices rep-resent the microstates of the discriminatory mechanism, suchas DNA polymerase in the case of replication. This graph isnaturally composed of two subgraphs, G X ( X = C, D ), cor-responding to the states in which substrate X is bound. G C and G D share a common vertex, but no edges, so that G has abutterfly shape.We will denote such a butterfly graph G = G C ⊕ v G D ,where v is the shared vertex. For the task of discrimination,the subgraphs G X are structurally symmetric, with symmetric (a) m’m’m mk’ C k’ D k C k D l’ D l + W C l + W D l’ C (b) ΔΔ energy Δ C CDD G C G D (c) m’k’ C C C m’l’ C C C k C l’ C C C C C D D energy FIG. 1. The Hopfield mechanism and the linear framework. (a) La-belled, directed “butterfly” graph for the original Hopfield mecha-nism in ref. [9], consisting of the subgraph G D for the incorrectsubstrate D (within left-hand dashed circle) and the subgraph G C for the correct substrate C (within right-hand dashed circle), whichshare the common vertex . Cyan and magenta denote correct and in-correct substrate binding, respectively. (b) Hypothetical energy land-scape for the Hopfield mechanism showing where energy may beexpended to drive the proofreading step with label m (cid:48) . (c) The span-ning trees of G C rooted at C (circled) are shown. A spanning tree isa subgraph which includes every vertex (spanning) and has no cycleswhen edge directions are ignored (tree); it is rooted at i if i is theonly vertex with no outgoing edges. Any non-root vertex has only asingle outgoing edge. Using equations (1) and (3), the trees shownhere give the left-hand factor in the denominator of equation (1) andthe remaining factors arise in a similar way. vertices, X , · · · , n X , of which C = 1 D = 1 is shared, andsymmetric edges, i C → j C if, and only if, i D → j D . Thelabels on these corresponding edges may, however, be distinct.The vertices i X with i > are the microstates in which X isbound, while vertex is the empty microstate in which no X is bound. All directed edges are assumed to be structurallyreversible, so that, if i X → j X , then j X → i X . The graphs G C , G D and G are therefore all strongly connected.Let G = G C ⊕ v G D be any butterfly graph. Even if G C and G D are not structurally symmetric, as above, it followsreadily from equation (1) that ρ i ( G C ⊕ v G D ) = (cid:26) ρ i ( G C ) ρ v ( G D ) if i ∈ G C ρ i ( G D ) ρ v ( G C ) if i ∈ G D . (3) IV. THE ERROR FRACTION FOR THE HOPFIELDMECHANISM
The original Hopfield kinetic proofreading mechanism isdescribed by the discriminatory butterfly graph G = G C ⊕ G D in Fig. 1(a). The substrates C and D are treated as “slow”components and assumed to have constant concentration overthe timescale of interest. These concentrations are absorbedinto the “on-rates” k (cid:48) C , k (cid:48) D , l (cid:48) C , l (cid:48) D . The discrimination mecha-nism itself is assumed to have the “fast” components and to beat steady state. The rate W for exit from X ( X = C, D ) cor-responds to product generation and release of X , so that theoverall system is open whenever W > , with C and D beingtransformed into correct and incorrect product, respectively. In this mechanism, discrimination occurs twice, throughbinding and unbinding of X to form X and to form X . It isassumed that unbinding, rather than binding, causes discrim-ination, as is often the case in biology [14], so that l (cid:48) C = l (cid:48) D and k (cid:48) C = k (cid:48) D . The correct substrate has a longer residencetime, so that k C < k D , which reflects the free energy differ-ence of ∆ > between C and D (Fig. 1(b)): if energy ismeasured in units of k B T , where k B is Boltzmann’s constantand T is the absolute temperature, then k D = k C e ∆ . Thereis assumed to be no difference in discrimination between X and X , so that k C /k D = l C /l D = e − ∆ < .Hopfield defines the steady-state error fraction, ε as theprobability ratio of the incorrect to the correct exit microstate,which, using equation (2), is given by ε = ρ D ( G ) /ρ C ( G ) ( ε is the inverse of the accuracy µ in the Abstract; we will workwith the former). Using equations (1) and (3), ε = [ l (cid:48) D ( k D + m (cid:48) ) + m (cid:48) k (cid:48) D ][( k C + m (cid:48) )( W + l C ) + mk C ][ l (cid:48) C ( k C + m (cid:48) ) + m (cid:48) k (cid:48) C ][( k D + m (cid:48) )( W + l D ) + mk D ] . (4)If the overall system remains closed, so that W = 0 , while themechanism operates at thermodynamic equilibrium, then ithas the error fraction, ε = k C /k D = l C /l D = e − ∆ (Supple-mentary Material). If the overall system becomes open, so that W > , while the mechanism remains at equilibrium, then ε increases monotonically with increasing W (SupplementaryMaterial). If the mechanism itself operates away from equi-librium, then ε > ε (cid:18) l C + m + Wl D + m + W (cid:19) > ε (5)for all positive values of the parameters (Supplementary Mate-rial). Hopfield shows that ε approaches the minimal error, ε ,as certain parametric quantities become small (SupplementaryMaterial) and suggests how this could be achieved in practiceby expending energy to drive the transition from X to X through the label m (cid:48) . This is kinetic proofreading.There are two aspects of Hopfield’s analysis which havenot always been fully appreciated. First, increasing m (cid:48) is notsufficient of itself for ε to approach ε . Indeed, it followsfrom equation (1) that, when W = 0 , ε → ε as m (cid:48) → ∞ .Too much energy expenditure can increase the error fraction,which behaves non-monotonically with respect to m (cid:48) . (Simi-lar non-monotonicity has been observed for kinetic proofread-ing with the T-cell receptor mechanism in Supplementary Fig.1 [18].) The parameter m (cid:48) must neither be too high nor toolow for the error fraction to approach ε . Second, parame-ters other than m (cid:48) , m and W must also take adequate val-ues for the accuracy to approach this bound: the “on-rate” for → X must be much larger than that for → X , so that l (cid:48) D /k (cid:48) D = l (cid:48) C /k (cid:48) C → (Supplementary Material). The lowerbound of ε is only reached asymptotically as several param-eters take limiting values.For more complex systems, the appropriate parameterregime for the minimal error is not readily found using Hop-field’s approach. We therefore sought an alternative strategy.If we let x = e ∆ = ε − and substitute k D = xk C and l D = xl C into equation (1), we see that, if no other parameterschange, the error fraction ε behaves like x − as x increases.We reasoned that to approach the minimal error of x − , thefold change in other parameter values should be some func-tion of x . By retaining only the highest-order term in x as x → ∞ , the behaviour of ε could be determined while by-passing the parametric complexity. Let R ( x ) ∼ Q ( x ) meanthat R ( x ) /Q ( x ) → c as x → ∞ , where < c < ∞ . Itcan be seen from equation (1) that if either k (cid:48) D = k (cid:48) C ∼ x or l (cid:48) C = l (cid:48) D ∼ x − , while none of the remaining parameters de-pend on x , then ε ∼ x − = ε . This scaling of the “on-rates”corresponds to what was required in the previous paragraphfor Hopfield’s limiting procedure. This suggests a strategy forexploring parameter space that can be extended to more com-plex systems. We exploit this below to examine the relationbetween energy, speed and accuracy. V. DISSOCIATION-BASED MECHANISMS
We introduce here a class of discrimination mechanisms forwhich such a relation can be determined. We consider a dis-criminatory butterfly graph of the form G = G C ⊕ G D con-sisting of structurally symmetric subgraphs G C and G D ofarbitrary complexity. The vertex n X is taken to be the onlyexit microstate in which product is generated, so that there isa return edge n X → . No further structural assumptions aremade but the product generation rate, W , makes an additivecontribution to the label of the return edge n X → , as in Fig.1(a).Multiple internal microstates and transitions are allowed in G X as well as multiple returns to the empty microstate, , al-though only a single one of these, through n X , also generatesproduct. As in Hopfield’s original scheme, we think of themechanism as coupled to sources and sinks of energy, whichmay alter the edge labels. In Hopfield’s scheme, the labels onedges which do not go to the reference microstate were as-sumed to be the same between C and D (Fig. 1(a)). In otherwords, there was no “internal discrimination” between correctand incorrect substrates. Here, we allow internal discrimina-tion between C and D : when j (cid:54) = 1 , the label on i C → j C may be different from that on i D → j D .Graphs of this form been widely employed in the litera-ture. In addition to the original Hopfield mechanism (Fig.1(a)), they include the “delayed” mechanism [10], the multi-step mechanism [12, 13, 40], the T-cell receptor mechanism[14, 18] and generalised proofreading mechanisms [15–17].We follow Hopfield in using the steady-state error fractionand work from now on with probabilities in the microscopicinterpretation. Let p ∗ be the vector of steady state probabil-ities. The discrimination error fraction, ε , is the steady-stateprobability ratio of the incorrect exit microstate, n D , to thecorrect exit microstate, n C , ε = p ∗ n D p ∗ n C . (6)We will analyse the behaviour of G under the assumptionthat some of the labels are functionally dependent on the non-dimensional variable x ∈ R . A function R ( x ) is said to be al- lowable if it is positive, R ( x ) > for x > , and has a degree, deg( R ) , given by R ( x ) ∼ x deg( R ) as x → ∞ . This is welldefined because x a ∼ x b if, and only if, a = b . The degreedetermines the asymptotics of allowable functions: R ∼ Q if,and only if, deg( R ) = deg( Q ) . Note that deg( R ) = 0 if, andonly, R ( x ) → c as x → ∞ , where c > , which is the case if R does not depend on x .The labels in the graph G are assumed to be allowable func-tions of x . (The product generation rate W couples the mech-anism to the environment and is assumed not to depend on x .)If R and Q are allowable functions, then so are R − , RQ and R + Q and (Supplementary Material) deg( R − ) = − deg( R )deg( RQ ) = deg( R ) + deg( Q )deg( R + Q ) = max(deg( R ) , deg( Q )) . (7)Accordingly, any rational function of the labels with only pos-itive terms, such as p ∗ , which acquires this structure throughequation (2) and equation (1), or ε , which acquires it throughequation (6), becomes in turn an allowable function of x .We define a dissociation-based mechanism to be a generaldiscrimination mechanism for which, for the edges betweenthe exit microstates and , deg( (cid:96) → n D ) = deg( (cid:96) → n C )deg( (cid:96) n D → ) = deg( (cid:96) n C → ) + 1 . (8)Here, we use (cid:96) i → j to denote the label on the edge i → j . Eq.8 is analogous to the assumption l (cid:48) C = l (cid:48) D and xl C = l D forthe Hopfield mechanism. Unlike the Hopfield mechanism, wedo not restrict what happens at non-exit microstates.With such general assumptions on the labels, a dissociation-based mechanism may not reach thermodynamic equilibrium.However, if it can, with W > , so that the overall systemremains open, then equation (8) ensures that the equilibriumerror fraction, ε eq , has a simple form. Since detailed balancerequires that each pair of edges is independently at steady state[33], the exit states, n X , satisfy (cid:96) n X → p ∗ n X = (cid:96) → n X p ∗ , sothat ε eq = p ∗ n D p ∗ n C = (cid:18) (cid:96) → n D (cid:96) → n C (cid:19) (cid:18) (cid:96) n C → (cid:96) n D → (cid:19) . (9)Applying equation (8) and using equation (6), we see that, ifequilibrium is reached, the resulting error fraction, ε eq , satis-fies ε eq ∼ x − . (10) VI. THE ASYMPTOTIC RELATION
We now define the measures of speed and energy expen-diture in terms of which our main result will be stated. Areasonable interpretation for the speed of the mechanism, σ ,is the steady-state flux of correct product [41], σ = W p ∗ n C . (11)As for energy expenditure, this is determined at steady stateby the rate of entropy production. Schnakenberg put forwarda definition of this [37] that has been widely used [2, 7]: fora pair of reversible edges, i (cid:29) j , the steady-state entropyproduction rate, P ( i (cid:29) j ) , is the product of the net flux, J ( i (cid:29) j ) = (cid:96) j → i p ∗ j − (cid:96) i → j p ∗ i , and the thermodynamic force, A ( i (cid:29) j ) = ln( (cid:96) j → i p ∗ j /(cid:96) i → j p ∗ i ) : P ( i (cid:29) j ) = J ( i (cid:29) j ) A ( i (cid:29) j ) . (12)Here, we omitted Boltzmann’s constant k B for convenience,so that P has units of (time) − . If T is the absolute tem-perature, then k B T P ( i (cid:29) j ) is the power irreversibly ex-pended through i (cid:29) j . The total entropy production rate ofthe system is then given by P = (cid:80) i (cid:29) j P ( i (cid:29) j ) . Note that P ( i (cid:29) j ) ≥ (and so also P ≥ ) with equality at ther-modynamic equilibrium when detailed balance implies that J ( i (cid:29) j ) = 0 . Positive entropy production, with P > ,signifies energy expenditure away from thermodynamic equi-librium.Both σ and P are functions of x and σ is evidently allow-able. However, J ( i (cid:29) j ) is not a rational function with pos-itive terms and ln( x ) (cid:54)∼ x α for any α , so P ( i (cid:29) j ) and P are not allowable functions. Nevertheless, the asymptotic be-haviour of P can be estimated. Some further notation is help-ful to do this. If R ( x ) and Q ( x ) are functions which are notnecessarily allowable, then R ≺ Q means that R/Q → as x → ∞ . This relation is transitive, so that, if S ≺ R and R ≺ Q , then S ≺ Q . If both functions are allowable, then R ≺ Q , if, and only if, deg( R ) < deg( Q ) . We will say that R (cid:45) Q if R/Q → c , where ≤ c < ∞ , and correspondingremarks about transitivity and allowable degrees hold for thisrelation. Note that ≺ and (cid:45) dominate over ∼ when formingproducts, so, for instance,if T (cid:45) S and R ∼ Q then RT (cid:45) SQ , (13)which we will make use of below.Each summand P ( i (cid:29) j ) has the form ( R − Q ) ln( R/Q ) ,where R and Q are allowable. Let α = deg( R ) , β = deg( Q ) and c = lim R/x α and c = lim Q/x β as x → ∞ . Bydefinition, c , c > . Note that, if S is allowable, then (Sup-plementary Material) ln( S ) ∼ (cid:26) ln( x ) if deg( S ) > x − ) if deg( S ) < . (14)The asymptotic behaviour of P ( i (cid:29) j ) then falls into thefollowing three cases (Supplementary Material), as specifiedon the right: case 1: α (cid:54) = β ∼ x max( α,β ) ln( x ) case 2: α = β , c (cid:54) = c ∼ x α case 3: α = β , c = c ≺ x α . (15)The third case is awkward because the leading-order asymp-totics are lost, which leads to the ≺ relation instead of ∼ .However, c and c are rational expressions in the parame-ters which do not involve x and the equation c = c de-fines a hypersurface in the space of those parameters. The reversible edges which fall into case 3 therefore determine aset of measure zero in the space of parameters. Provided thisset is avoided, the asymptotic behaviour of the summands in P fall into the first two cases and can be controlled. In particu-lar, suppose that the total entropy production rate P is writtenas P = (cid:80) u P u , where P u is a term coming from a pair of re-versible edges i (cid:29) j , as in equation (12). In Appendix A, weshow that, if P k is any summand in case 1 of equation (15),then, outside the measure-zero set defined by case 3, P k (cid:45) P .Let us now assume, for any dissociation-based mechanismas defined previously, that ε ( x ) ≺ x − . (16)This forces the error fraction to be asymptotically better thanif the system were able to reach equilibrium (equation (10))and thereby ensures that energy expenditure is contributing toan improvement in accuracy. Consider any general discrimi-nation mechanism which is dissociation-based, as described inequation (8). If its error fraction obeys equation (16) then, out-side the measure-zero set in parameter space defined by case 3of equation (15), we show in Appendix B that the mechanismsatisfies the asymptotic relation, σ ln( ε − ) (cid:45) P . (17)The exact asymptotics of σ ln( ε − ) /P are difficult to es-timate for a general dissociation-based mechanism with al-lowable labels because each pair of reversible edges must beexamined. However, for the Hopfield mechanism (Fig. 1(a)),under the conditions described above for which ε ∼ x − , wefind (Supplementary Material) lim x →∞ σ ln( ε − ) P = 2 Wl c + W (18)outside the parametric set of measure-zero noted above. VII. A NON-DISSOCIATION BASED MECHANISM
The requirements in equation (8) for being dissociation-based are necessary for the validity of equation (17). In theSupplementary Material, we consider a discrimination mech-anism with a structure identical to that of the Hopfield mecha-nism (Fig. 1(a)) but with labels that do not follow equation (8)(Supplementary Fig. 3). If the mechanism reaches thermody-namic equilibrium, then it follows from equation (9) that itsequilibrium error fraction satisfies ε eq ∼ x − . However, witha particular choice of allowable functions for the labels, forwhich the mechanism is no longer at equilibrium, its errorfraction improves asymptotically, with ε ∼ x − / , while itsspeed remains constant, σ ∼ , and its entropy production iseither constant or vanishes, P (cid:45) , outside a set of measurezero. This evidently does not obey equation (17) and showsthe existence of a different asymptotic interplay between en-ergy, speed and accuracy. VIII. NUMERICAL CALCULATIONS OUTSIDE THEASYMPTOTIC REGIME
To examine further the energy-speed-accuracy relationfound by the asymptotic analysis above, we used more restric-tive assumptions on the allowable labels to facilitate numeri-cal exploration. We considered discrimination-based mecha-nisms in which the x -dependency was similar to Hopfield’soriginal analysis. For any return edge to from a non-exitmicrostate, we assumed that (cid:96) i D → = (cid:96) i C → x ( i (cid:54) = n ) , (19)with an additive contribution of W in the exit microstate ( i = n ), (cid:96) n D → = ax + W , (cid:96) n C → = a + W with a ∈ R > . Asfor the other edges, we assumed no internal discrimination, sothat the labels were the same for C and D , (cid:96) i D → j D = (cid:96) i C → j C ( j (cid:54) = 1) , (20)with no x -dependency. By equation (9), the equilibrium er-ror function when the system is closed ( W = 0 ) satisfies ε = x − . We set x = e , sampled the values ln( (cid:96) i C → j C ) , ln( a ) , and ln( W ) uniformly in [ − , , and determined (cid:96) i D → j D from equations (19) and (20). We plotted P/σ against ln( ε /ε ) , when ε < ε , for the Hopfield mechanism (Fig.2(a)), the T-cell receptor mechanism (Supplementary Fig. 1)and for a mechanism different from both of these (Supplemen-tary Fig. 2). In each case, the resulting region was confinedto the left of a vertical line (Fig. 2(a), black dashed line) andabove the diagonal (Fig. 2(a), red dashed line). For the Hop-field mechanism, the vertical bound comes from equation (5)and similar bounds on ε exist for the other mechanisms (notshown). The diagonal bound, however, is unexpected and im-plies the bound σ ln( ε /ε ) < P ( for ε < ε ) (21)for finite parameter values. It is possible that equation (21)holds for any discrimination-based mechanism whose edge la-bels satisfy equations (19) and (20).The calculations leading to equation (21) assumed no inter-nal discrimination between correct and incorrect substrate, asspecified in equation (20). We were interested to find that ex-perimental data for ribosomes and DNA polymerase, based onthe original Hopfield mechanism, showed substantial internaldiscrimination, extending even to the product generation rate W [19]. To examine the impact of this, we proceeded as fol-lows. For any return edge to from a non-exit microstate, weintroduced an asymmetry between C and D so that (cid:96) i D → = α i (cid:96) i C → x ( i (cid:54) = n ) . (22)For the exit state ( i = n ), the product generation rate makesan additive contribution, W X , which now depends on the sub-strate X , so that (cid:96) n C → = a + W C and (cid:96) n D → = α n ax + W D , where a ∈ R > and W D = α W W C . For the otheredges, we similarly introduced an asymmetry (cid:96) i D → j D = α ij (cid:96) i C → j C ( j (cid:54) = 1) . (23)
0 20 40 P / σ ln( ε /ε)
10 30 10 20 30 40 (a)
0 20 40 P / σ ln( ε /ε)
10 30 10 20 30 40 (b) P ε ε σ P ε ε σ *** FIG. 2. Numerics for the Hopfield mechanism. (a) Plot of
P/σ against ln( ε /ε ) for the Hopfield mechanism (Fig. 1(a)) for approx-imately points. The sampling and the dashed lines are describedin the text. (b) Similar plot to (a) for the Hopfield mechanism withinternal discrimination between correct and incorrect substrates, asdescribed in the text, with the light blue points having a lower asym-metry range ( A = 1 ) and the dark blue points having a higher range( A = 5 ). The coloured overlays represent values from experimentaldata for ribosomes (orange, red and brown regions) and DNA poly-merase (green point), with the former being samples of values esti-mated for a parameter for which no experimental data exists. Onlythose samples for which ε > ε are shown and the asterisks, *, givethe averages of the plotted values. The inset gives the plotted av-erages (for the ribosome variants) and values (for DNAP) of errorfractions, ε and ε , speed, σ and entropy production rate, P (all inunits of s − ). The data from which these values were calculated areshown in Supplementary Table 1. See [42] and the caption of Sup-plementary Table 1 for more details. The multiplicative factors α ij and α W carry the asymmetrybetween C and D in internal discrimination.In view of the asymmetry in product generation rates, it isnatural to redefine the error fraction as ε = W D p ∗ n D W C p ∗ n C . Using equation (9), the equilibrium error fraction whenthe system is closed ( W = 0 ) is given by ε = (cid:96) → n D a/ ( (cid:96) → n C α n a ) .We chose the asymmetry factors by sampling ln( α ij ) and ln( α W ) uniformly in the range [ − A, A ] , for A = 1 and A = 5 , and chose the other parameters as described previ-ously for Fig. 2(a). Fig. 2(b) shows that both the verticalbound and the diagonal bound in Fig. 2(a) are broken, with theextent of the breach increasing with increase in the asymme-try range from A = 1 (Fig. 2(b), light blue points) to A = 5 (Fig. 2(b), dark blue points). Similar results were found forthe other mechanisms that we numerically calculated (Sup-plementary Figs. 1(c) and 2(c)). We see that the absence ofinternal discrimination is essential for the vertical and diago-nal bounds shown in Fig. 2(a) and Supplementary Figs. 1(b)and 2(b).Banerjee et al. have provided parameter values for the Hop-field mechanism based on experimental data for discrimina-tion in mRNA translation by the E. coli ribosome, includ-ing also an error-prone and a hyperaccurate mutant, and inDNA replication by the bacteriophage T7 DNA polymerase(DNAP) [19]. We used these parameter values to calculateentropy production, speed and accuracy as defined here andoverlaid the resulting ln( ε /ε ) , P/σ points on the previousnumerical calculation (Fig. 2(b)) [42].The data show a striking difference between mRNA trans-lation and DNA replication (Fig. 2(b)). All three ribosomevariants (orange, wild type; brown, error-prone; red, hyper-accurate) have much higher
P/σ values than DNAP (green),with the former lying comfortably above the diagonal boundgiven by (21) and the latter lying well below. Nevertheless,all systems exhibit substantial internal discrimination (Sup-plementary Table 1). As the inset in Fig. 2(b) shows, theseparation between translation and replication arises from adecrease of two orders of magnitude in entropy productionrate and an increase of two orders of magnitude in speed. Fur-thermore, DNAP not only shows the smallest error fraction, ε ,by three orders of magnitude, but also the greatest fold changeover the equilibrium error fraction, ε /ε . In contrast, the ri-bosome variants, while showing the expected differences inerror fraction, have lower fold changes over their equilibriumerror fractions. Evolution seems to have tuned the energy dis-sipation, speed and accuracy of the replication machinery to amuch greater degree than the translation machinery. IX. DISCUSSION
The relationship between energy expenditure and desir-able features, like accuracy and speed in discrimination, havebeen the subject of many studies, as cited in our references.One of the challenges here is what we have called “his-tory dependence” [8, 32]. If a linear framework graph is atthermodynamic equilibrium, then the steady-state probabilityof a microstate can be calculated by multiplying the ratios, (cid:96) ( i → j ) /(cid:96) ( j → i ) , along any path to the microstate from ;detailed balance ensures that all paths give the same value.In this sense, the graph is “history independent” at steadystate. Away from equilibrium, not only does the calculationof steady-state probabilities become history dependent, in thesense that different paths yield different values, but, as equa-tion (1) reveals, every path contributes to the final answer. TheMatrix-Tree Theorem does the bookkeeping for this calcula-tion using rooted spanning trees.The resulting combinatorial explosion in enumerating span-ning trees can be super-exponential in the size of the graph[8]. This difficulty may have been apparent to earlier work-ers like Hill [36] and Schnakenberg [37] and may have dis-couraged a more analytical approach. The combinatorialcomplexity has largely been avoided by focussing on sim-ple or highly-structured examples and by astute use of ap-proximation. It is only recently, with the linear framework[8] and the re-engagement with earlier work [16], that non-equilibrium steady-state calculations have been calculated an-alytically without such restrictions.In this paper, we have developed a way to address thiscomplexity that is inspired by Hopfield’s analysis of kineticproofreading. Here, the minimum error fraction can only bereached asymptotically (equation (5)) and only when multiple labels change their values consistently. This has suggesteda method of exploring parameter space by treating the labelsas allowable functions of a scaling variable x . In this way, asystem of arbitrary structure can be analysed away from equi-librium, with relaxed assumptions on how energy is being de-ployed, while rising above the combinatorial explosion fromhistory dependence.Perhaps the most interesting conclusion from this analysisis the emergence of the quantity σ ln( µ ) /P . Our main result,as expressed in equation (17), says that this quantity is asymp-totically finite, for any graph obeying the dissociation-basedcondition on exit edges (equation (8)) and for any scheme ofallowable scaling through which energy increases ( deg > )or decreases ( deg < ) the rates, provided that the accuracyimproves over equilibrium (equation (16)).The advantage of the asymptotic analysis undertaken hereis that it reveals a universal behaviour in σ ln( µ ) /P that tran-scends network structure and parametric complexity. Inter-estingly, our numerical calculations suggest that universalitymay still be found for finite parameter values, in the form ofthe bound in equation (21), as shown in Fig. 2 and Supple-mentary Figs. 1 and 2. However, this bound depends cruciallyon the absence of internal discrimination between correct andincorrect substrates, in contrast to the asymptotic behaviourin equation (17), for which internal discrimination is allowed.Experimental data shows that evolution has discriminated in-ternally to a substantial extent but with very different effectson this bound. All E. coli ribosome variants for which we havedata comfortably obey the bound, while the T7 DNA poly-merase breaks it. This reflects a striking reduction in
P/σ forthe latter, with far less difference between the ribosomes andthe DNA polymerase in the fold change over their equilibriumerror fractions (Fig. 2(b)). It would be interesting to know ifthese same comparative relationships are maintained for otherribosomes and polymerases. While recent work has shownthat local trade-offs between speed and accuracy can differmarkedly between different parametric regions [19], the quan-tities introduced here may be helpful for more global compar-isons between discriminatory mechanisms.A similar quantity to σ ln( µ ) /P has emerged in the workof Tu and colleagues on sensory adaptation, using quite dif-ferent methods [25], suggesting that it may be significant fora broader context of cellular information processing that in-cludes discrimination and adaptation. Indeed, the analogy be-tween discrimination and adaptation has already been drawn[3]. Because of their generality, the methods used here may beparticularly useful for developing such a broader perspective. ACKNOWLEDGMENTS
F.W. was supported by the National Science Foun-dation (NSF) Graduate Research Fellowship under grantDGE1144152. A.A. was supported by the Alfred P. SloanFoundation. J.G. was supported by NSF grant 1462629. Wethank P.-Y. Ho for discussions.
APPENDIX A: PROOF OF P k (cid:45) P Suppose that P k ∼ x α k ln( x ) . If P i is also in case 1 and P i ∼ x α i ln( x ) , then P i /P k → c , where c = 0 , , ∞ , de-pending on the relative values of α i and α k . Similarly, if P j is in case 2 and P j ∼ x α j , then P j /P k → c , where c = 0 , ∞ ,depending on the relative values of α j and α k . By assump-tion, there are no other cases to consider (if P h were in case 3,we could not estimate lim x →∞ P h /P k ). Since P k is one of thesummands in P , it follows that P/P k → c , where ≤ c ≤ ∞ .Equivalently, P k /P → c − , where ≤ c − ≤ . In particu-lar, P k (cid:45) P , as required. APPENDIX B: PROOF OF THE ASYMPTOTIC RELATION
Suppose first that (cid:29) n C falls into case 1 in equation (15).Let α = deg( (cid:96) n c → p ∗ n C ) and β = deg( (cid:96) → n C p ∗ ) , so that α (cid:54) = β . Then, x α ln( x ) (cid:45) x max( α,β ) ln( x ) ∼ P (1 (cid:29) n C ) .Since the product generation rate, W , appears additively, (cid:96) n C → = W + U ( x ) , for some allowable function U . It fol-lows from equation (7) that α = deg( (cid:96) n C → ) + deg( p ∗ n C ) =max(0 , deg( U )) + deg( p ∗ n C ) ≥ deg( W p ∗ n C ) = deg( σ ) .Hence, σ (cid:45) x α . Furthermore, since equation (16) tellsus that deg( ε ) < − , it follows from equation (14) that ln( ε − ) ∼ ln( x ) . Using equation (13), we deduce that σ ln( ε − ) (cid:45) P (1 (cid:29) n C ) . If (cid:29) n C does not fall into case 1 in equation (15),then α = β . Let us then consider (cid:29) n D . Accord-ing to equations (6) and (16), deg( p ∗ n D ) < deg( p ∗ n C ) − .Using equation (7) to combine this with equation (8.2), wesee that deg( (cid:96) n D → p ∗ n D ) < deg( (cid:96) n C → p ∗ n C ) = α = β =deg( (cid:96) → n C p ∗ ) . But now, by equation (8.1) and equation (7), deg( (cid:96) → n C p ∗ ) = deg( (cid:96) → n D p ∗ ) . (24)It follows that deg( (cid:96) n D → p ∗ n D ) < deg( (cid:96) → n D p ∗ ) , (25)so that (cid:29) n D falls into case 1 even though (cid:29) n C doesnot. Therefore, by equation (15), P (1 (cid:29) n D ) ∼ x γ ln( x ) ,in which, because of equation (25), γ = deg( (cid:96) → n D p ∗ ) . Butaccording to equation (24), γ = deg( (cid:96) → n C p ∗ ) = β = α .Hence, by the same argument as above for (cid:29) n C , we de-duce that σ ln( ε − ) (cid:45) P (1 (cid:29) n D ) . We can now appeal to the result in Appendix A to completethe proof. [1] Berg, H. C. and Purcell, E. M.
Physics of chemoreception , Bio-phys. J. , 193 (1977).[2] Mehta, P. and Schwab, D. J. Energetic costs of cellular compu-tation , Proc. Natl. Acad. Sci. USA , 17978 (2012).[3] Hartich, D., Barato, A. C. and Seifert, U.
Nonequilibrium sens-ing and its analogy to kinetic proofreading , New J. Phys. ,055026 (2015).[4] ten Wolde, P. R., Becker, N. B., Ouldridge, T. E. and Mugler,A. Fundamental limits to cellular sensing , J. Stat. Phys. ,1395 (2016).[5] Singh, V. and Nemenman, I.
Simple biochemical networks al-low accurate sensing of multiple ligands with a single receptor ,PLoS Comput. Biol. , e1005490 (2017).[6] Qian, H. and Cooper, J. A. Temporal cooperativity and sensi-tivity amplification in biological signal transduction , Biochem-istry , 2211 (2008).[7] Tu, Y. The nonequilibrium mechanism for ultrasensitivity in abiological switch: sensing by Maxwell’s demons , Proc. Natl.Acad. Sci. USA , 11737 (2008).[8] Estrada, J., Wong, F., DePace, A. and Gunawardena, J.
Infor-mation integration and energy expenditure in gene regulation ,Cell , 234 (2016).[9] Hopfield, J. J.
Kinetic proofreading: a new mechanism for re-ducing errors in biosynthetic processes requiring high speci-ficity , Proc. Natl. Acad. Sci. USA , 4135 (1974).[10] Ninio, J. Kinetic amplification of enzyme discrimination , Bio-chemie , 587 (1975).[11] Bennett, C. H. Dissipation error tradeoff in proofreading , BioSystems , 85 (1979).[12] Blomberg, C. and Ehrenberg, M. Energy considerations forkinetic proofreading in biosynthesis , J. Theor. Biol. , 631(1981).[13] Savageau, M. A. and Lapointe, D. S. Optimization of kineticproofreading: a general method for derivation of the constraintrelations and an exploration of a specific case , J. Theor. Biol. , 157 (1981).[14] McKeithan, T. W. Kinetic proofreading in T-cell receptor signaltransduction , Proc. Natl. Acad. Sci. USA , 5042 (1995).[15] Murugan, A., Huse, D. A. and Leibler, S. Speed, dissipation,and error in kinetic proofreading , Proc. Natl. Acad. Sci. USA , 12034 (2012).[16] Murugan, A., Huse, D. A. and Leibler, S.
Discriminatory proof-reading regimes in nonequilibrium systems , Phys. Rev. X ,021016 (2014).[17] Sartori, P. and Pigolotti, S. Thermodynamics of error correc-tion , Phys. Rev. X , 041039 (2015).[18] Cui, W. and Mehta, P. Optimality in kinetic proofreading andearly T-cell recognition: revisiting the speed, energy, accuracytrade-off (2017). Arxiv.org/abs/1703.03398v2.[19] Banerjee, K., Kolomeisky, A. B. and Igoshin, O. A.
Elucidatinginterplay of speed and accuracy in biological error correction ,Proc. Natl. Acad. Sci. USA , 5183 (2017).[20] Banerjee, K., Kolomeisky, A. B. and Igoshin, O. A.
Accuracyof substrate selection by enzymes is controlled by kinetic dis-crimination , J. Phys. Chem. Lett. , 1552 (2017).[21] Rao, R. and Peliti, L. Thermodynamics of accuracy in kinetic proofreading: dissipation and efficiency trade-offs , J. Stat.Mech. Theor. Exp. , P06001 (2015).[22] Ehrenberg, M. and Blomberg, C.
Thermodynamic constraintson kinetic proofreading in biosynthetic pathways , Biophys. J. , 333 (1980).[23] Freter, R. R. and Savageau, M. A. Proofreading systems of mul-tiple stages for improved accuracy of biological discrimination ,J. Theor. Biol. , 99 (1980).[24] Sartori, P. and Pigolotti, S. Kinetic versus energetic discrimina-tion in biological copying , Phys. Rev. Lett. , 188101 (2013).[25] Lan, G., Sartori, P., Neumann, S., Sourjik, V. and Tu, Y.
Theenergy-speed-accuracy trade-off in sensory adaptation , Nat.Phys. , 422 (2012).[26] Holy, T. E. and Leibler, S. Dynamic instability of microtubulesas an efficient way to search in space , Proc. Natl. Acad. Sci.USA , 5682 (1994).[27] Lang, A. H., Fisher, C. K., Mora, T. and Mehta, P. Thermody-namics of statistical inference by cells , Phys. Rev. Lett. ,148103 (2014).[28] Sartori, P., Granger, L., Lee, C. F. and Horowitz, J. M.
Thermo-dynamic costs of information processing in sensory adaptation ,PLoS Comp. Biol. , e1003974 (2014).[29] Parrondo, J. M. R., Horowitz, J. M., and Sagawa, T. Thermody-namics of information , Nat. Phys. , 131 (2015).[30] Still, S., Sivak, D. A., Bell, A. J., and Crooks, G. E. Thermody-namics of prediction , Phys. Rev. Lett. , 120604 (2012).[31] Das, J.
Limiting energy dissipation induces glassy kinetics insingle-cell high-precision responses , Biophys. J. , 1180(2016).[32] Ahsendorf, T., Wong, F., Eils, R. and Gunawardena, J.
A frame-work for modelling gene regulation which accommodates non-equilibrium mechanisms , BMC Biol. , 102 (2014).[33] Gunawardena, J. A linear framework for time-scale separa-tion in nonlinear biochemical systems , PLoS ONE , e36321(2012). [34] Thomson, M. and Gunawardena, J. Unlimited multistability inmultisite phosphorylation systems , Nature , 274 (2009).[35] Dasgupta, T., Croll, D. H., Owen, J. A., Vander Heiden, M. G.,Locasale, J. W., Alon, U., Cantley, L. C. and Gunawardena, J.
Afundamental trade off in covalent switching and its circumven-tion by enzyme bifunctionality in glucose homeostasis , J. Biol.Chem. , 13010 (2014).[36] Hill, T. L.
Studies in irreversible thermodynamics IV. Diagram-matic representation of steady state fluxes for unimolecular sys-tems , J. Theoret. Biol. , 442 (1966).[37] Schnakenberg, J. Network theory of microscopic and macro-scopic behaviour of master equation systems , Rev. Mod. Phys. , 571 (1976).[38] Mirzaev, I. and Gunawardena, J. Laplacian dynamics on gen-eral graphs , Bull. Math. Biol. , 2118 (2013).[39] Gunawardena, J. Time-scale separation: Michaelis andMenten’s old idea, still bearing fruit , FEBS J. , 473 (2014).[40] Johansson, M., Lovmar, M. and Ehrenberg, M.
Rate and accu-racy of bacterial protein synthesis revisited , Curr. Opin. Micro-biol. , 141 (2008).[41] Hill, T. L. Free Energy Transduction and Biochemical CycleKinetics (Dover Publications, New York, USA, 2004).[42] We note that Banerjee et al. imposed one additional constrainton their numerical values, following [3], by assuming that thecorrect and incorrect substrates used the same external chemi-cal potential to break thermodynamic equilibrium; see equation(2) in [19]. While this is reasonable, we have followed Hopfieldhere and not made that assumption in our analysis and calcula-tions. We also had to estimate two parameter values for the ri-bosome variants in Fig. 2(b), for which experimental data werenot available, and did so by random sampling, as explained inSupplementary Table 1. upplementary Material: An energy-speed-accuracy relation in complex networks for biologicaldiscrimination
Felix Wong, , Ariel Amir, and Jeremy Gunawardena , ∗ School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA Department of Systems Biology, Harvard Medical School, Boston, MA 02115, USA
We provide proofs here of the mathematical assertionsmade in the main text.
I. EQUILIBRIUM ERROR FRACTION FOR THEHOPFIELD MECHANISM
The error fraction, ε , for the Hopfield mechanism is givenin equation (4) of the main text, and is repeated here for con-venience, ε = [ l (cid:48) D ( k D + m (cid:48) ) + m (cid:48) k (cid:48) D ][( k C + m (cid:48) )( W + l C ) + mk C ][ l (cid:48) C ( k C + m (cid:48) ) + m (cid:48) k (cid:48) C ][( k D + m (cid:48) )( W + l D ) + mk D ] . (1)The background assumptions, as mentioned in the main text,are k (cid:48) C = k (cid:48) D , l (cid:48) C = l (cid:48) D and k C /k D < .If the mechanism is at thermodynamic equilibrium, then de-tailed balance must be satisfied. The equivalent cycle condi-tion [1] applied to the two cycles in Fig. 1(a) of the main textyields m (cid:48) m = l (cid:48) C k C l C k (cid:48) C = l (cid:48) D k D l D k (cid:48) D . (2)Note that W does not appear in equation (2) since, althoughthe mechanism itself is at thermodynamic equilibrium, thesystem remains open, with substrate being converted to prod-uct. Denote by ε eq the value of ε under the equilibrium con-straint in equation (2). Using equation (2), define α , β so that α = l C l (cid:48) C = mk C m (cid:48) k (cid:48) C and β = l D l (cid:48) D = mk D m (cid:48) k (cid:48) D . Using the background assumptions, define the quantity ε ,given by ε = αβ = l C l D = k C k D , (3)to which a physical interpretation will be given shortly. Sub-stituting α and β into equation (1) and rewriting, we see that ε eq = W.A + αW.B + β , where A = k C + m (cid:48) ( k C + m (cid:48) ) l (cid:48) C + m (cid:48) k (cid:48) C , B = k D + m (cid:48) ( k D + m (cid:48) ) l (cid:48) D + m (cid:48) k (cid:48) D . ∗ Corresponding author: [email protected] If W = 0 , then by equation (3), ε eq = ε , which showsthat ε , as defined in equation (3), is the error fraction for theclosed system at thermodynamic equilibrium as defined in themain text.We now want to prove that ε eq increases from ε as W in-creases, for which it is sufficient to show that dε eq /dW > .For this, dε eq dW = Aβ − Bα ( W.B + β ) , so that dε eq /dW > if, and only if, A/B > α/β . We have AB = (cid:18) k C + m (cid:48) k D + m (cid:48) (cid:19) (cid:18) ( k D + m (cid:48) ) l (cid:48) D + m (cid:48) k (cid:48) D ( k C + m (cid:48) ) l (cid:48) C + m (cid:48) k (cid:48) C (cid:19) . (4)The following result is straightforward. Lemma.
Consider the rational function r ( x ) = ( a + x ) / ( b + x ) , where a, b > . If a/b > , then r ( x ) decreases strictlymonotonically from r (0) = a/b to . If a/b < , then r ( x ) increases strictly monotonically from r (0) = a/b to . Applying the Lemma repeatedly to the terms in equation(4), and recalling the background assumptions, we see that k C + m (cid:48) k D + m (cid:48) > k C k D and ( k D + m (cid:48) ) l (cid:48) D + m (cid:48) k (cid:48) D ( k C + m (cid:48) ) l (cid:48) C + m (cid:48) k (cid:48) C > . Hence,
A/B > k C /k D = α/β and so dε eq /dW > . Itfollows that ε eq increases strictly monotonically from ε as W increases from . II. DERIVATION OF EQUATION (5) OF THE MAIN TEXT
The non-equilibrium error fraction in equation (1) can berewritten as ε ( m (cid:48) ) = u ( m (cid:48) ) v ( m (cid:48) ) , where u ( m (cid:48) ) = [ l (cid:48) D k D + ( l (cid:48) D + k (cid:48) D ) m (cid:48) ][ l (cid:48) C k C + ( l (cid:48) C + k (cid:48) C ) m (cid:48) ] , and v ( m (cid:48) ) = [ k C ( l C + m + W ) + ( l C + W ) m (cid:48) ][ k D ( l D + m + W ) + ( l D + W ) m (cid:48) ] . Using the background assumptions and the Lemma, we seethat, as m (cid:48) increases, u ( m (cid:48) ) decreases hyperbolically from u (0) = k D /k C = ε − to while v ( m (cid:48) ) increases hyper-bolically between v (0) = ε (cid:18) l C + m + Wl D + m + W (cid:19) and v ( ∞ ) = (cid:18) l C + Wl D + W (cid:19) . Hence, ε (0) = ( l C + m + W ) / ( l D + m + W ) and ε ( m (cid:48) ) → ( l C + W ) / ( l D + W ) as m (cid:48) → ∞ . Furthermore, ε ( m (cid:48) ) > v (0) = ε (cid:18) l C + m + Wl D + m + W (cid:19) > ε , (5)as required for equation (5) in the main text. III. THE LIMITING ARGUMENT IN KINETICPROOFREADING
The argument for kinetic proofreading given by Hopfield[2] is based on the non-dimensional quantities, δ = l (cid:48) D ( m (cid:48) + k D ) m (cid:48) k (cid:48) D , δ = m (cid:48) k C , δ = ml D , δ = Wl D , which are to be taken very small. Accordingly, we considerthe non-equilibrium error-fraction, ε , as defined in equation(1), in the limit as these four quantities → . Since k D > k C ,we have that δ > l (cid:48) C ( m (cid:48) + k C ) m (cid:48) k (cid:48) C > . (6)If we now divide above and below by m (cid:48) k (cid:48) D = m (cid:48) k (cid:48) C in theexpression ε = uv introduced above, we get ε = (cid:18) δ + 1 l (cid:48) C ( m (cid:48) + k C ) /m (cid:48) k (cid:48) C + 1 (cid:19) v . If we take δ → and formally treat v as a constant, wesee from equation (6) that ε → v as δ → . However, theexpression for δ involves m (cid:48) and this is also a parameter in v . Hence, v is not constant during the limiting process, whichhas coupled m (cid:48) to the values of the other parameters. If weignore this coupling, we can divide above and below in v by k C to get v = (1 + δ )( W + l C ) + m ( k D /k C + δ )( W + l D ) + mk D /k C and if then let δ → , we see that v → ε (cid:18) W + l C + mW + l D + m (cid:19) . If we now divide above and below in this expression by l D ,we get W/l D + l C /l D + m/l D W/l D + 1 + m/l D → ε as δ → and δ → . Hence, putting the sequence of limitstogether, we have, formally, lim δ → lim δ → lim δ → lim δ → ε = ε . This seems to be the interpretation that has been given inthe literature to Hopfield’s assertion that kinetic proofreadingachieves the error fraction of ε . The coupling noted above specifically affects m (cid:48) , which hasto satisfy two conditions. On the one hand m (cid:48) has to be large,in order that u should be close to u ( ∞ ) = 1 . That is the roleof δ → . On the other hand m (cid:48) has to be small, in order that v should be close to v (0) . That is the role of δ → . Theremaining limits for δ and δ are only there to make sure that v (0) → ε ; compare equation (5). The consequence of thecoupling between the δ and δ limits, which arises through m (cid:48) , can be seen by rewriting δ , δ = (cid:18) l (cid:48) D k (cid:48) D (cid:19) (cid:18) (cid:18) k D k C (cid:19) δ (cid:19) . Hence, in the limit as δ → and δ → , l (cid:48) D k (cid:48) D = l (cid:48) C k (cid:48) C → . In order to achieve the proofreading limit, it is necessary forrates other than m (cid:48) to change. Specifically, the “on rates” forthe first discrimination, k (cid:48) C = k (cid:48) D , must become large withrespect to those for the second discrimination, l (cid:48) C = l (cid:48) D . IV. DERIVATION OF EQUATION (7) OF THE MAIN TEXT
Suppose that R ( x ) and Q ( x ) are allowable functions, as de-fined in the main text, with R ( x ) /x α → c and Q ( x ) /x β → c , as x → ∞ , where c , c > .Since R − /x − α → ( c ) − > , it follows that R − isallowable and deg( R − ) = − deg( R ) . Since ( RQ ) /x α + β → c c > , it follows that RQ is allowable and deg( RQ ) =deg( R )+deg( Q ) . Finally, suppose, without loss of generality,that max( α, β ) = α , so that α ≥ β . Then, R + Qx α = (cid:18) Rx α (cid:19) + (cid:18) Qx β (cid:19) x β − α . (7)The limit of this, as x → ∞ , is c > , if α > β , or c + c > , if α = β . In either case, the limit is positive. Hence, R + Q is allowable and deg( R + Q ) = max(deg( R ) , deg( Q )) . Thisproves equation (7) of the main text. V. DERIVATION OF EQUATION (14) OF THE MAIN TEXT
Suppose that S ( x ) is an allowable function and that S/x α → c > , where α = deg( S ) . Since ln is a contin-uous function, ln( S ( x )) − α ln( x ) → ln( c ) . Dividing through by ln( x ) , we see that ln( S ( x ))ln( x ) → α . (8)If deg( S ) = α > , then ln( S ) ∼ ln( x ) , while if deg( S ) < then ln( S ) ∼ ln( x − ) , which proves equation (14) of themain text. VI. DERIVATION OF EQUATION (15) OF THE MAINTEXT
Note that if R ( x ) , Q ( x ) , S ( x ) , T ( x ) are functions, not nec-essarily allowable, and if R ∼ Q and S ∼ T , so that R/Q → c > and S/T → c > , then ( RS ) / ( QT ) → c c > ,so that RS ∼ QT . We will use this without reference below.Following the discussion in the main text, consider P ( i (cid:29) j ) = ( R − Q ) ln( R/Q ) where R and Q are allowable func-tions with R/x α → c > and Q/x β → c > . There arethree cases to consider.Suppose that α (cid:54) = β . If α > β , then the same argument asin equation (7) shows that ( R − Q ) ∼ x α . By equation (7)of the main text, deg( R/Q ) > , so that ln( R/Q ) ∼ ln( x ) .Hence, P ( i (cid:29) j ) ∼ x α ln( x ) . If α < β , then ( R − Q ) ∼ − x β and ln( R/Q ) ∼ − ln( x ) , so that P ( i (cid:29) j ) ∼ x β ln( x ) . Thisproves case 1.Suppose that α = β but c (cid:54) = c . Then, ( R − Q ) /x α → c − c (cid:54) = 0 . Also, RQ = (cid:18) Rx α (cid:19) (cid:18) x α Q (cid:19) → c c > . Hence, ln(
R/Q ) → ln( c /c ) . Since ( c − c ) ln( c /c ) > ,it follows that P ( i (cid:29) j ) ∼ x α , which proves case 2.Suppose that α = β and c = c . Then ( R − Q ) /x α → and R/Q → , so that ln( R/Q ) → . Hence, P ( i (cid:29) j ) /x α → as x → ∞ , so that P ( i (cid:29) j ) ≺ x α , which provescase 3. VII. DERIVATION OF EQUATION (20) OF THE MAINTEXT
For the Hopfield mechanism (Fig. 1(a) of the main text),we described in the main text how the asymptotic error rateof ε ∼ x − could be achieved, by assuming that the labelsare allowable functions of x such that: deg( l (cid:48) D ) = deg( l (cid:48) C ) = − , deg( k D ) = deg( l D ) = 1 and deg( k (cid:48) C ) = deg( k C ) =deg( k (cid:48) D ) = deg( m (cid:48) ) = deg( m ) = deg( l C ) = deg( W ) =0 . Using equations (1) and (3) of the main text, we find that deg( p ∗ ) = deg( p ∗ C ) = deg( p ∗ C ) = 0 , deg( p ∗ D ) = − , and deg( p ∗ D ) = − . It is helpful to introduce the notation R ≈ Q , for functions R ( x ) , Q ( x ) which may not be allowable, tosignify that lim x →∞ R/Q = 1 . We can use this to calculatethe asymptotic behaviour of the terms P ( i (cid:29) j ) in the entropyproduction rate (equation (12) of the main text). For instance, P (1 (cid:29) C ) = (( l C + W ) p ∗ C − l (cid:48) C p ∗ ) ln (cid:18) ( l C + W ) p ∗ C l (cid:48) C p ∗ (cid:19) , The only term in this expression which depends on x is l (cid:48) C for which deg( l (cid:48) C ) = − . Since deg(ln( R )) ≈ deg( R ) ln( x ) (equation (8)), it follows that P (1 (cid:29) C ) ≈ ( l C + W ) p ∗ C ln( x ) . Similar calculations yield P (2 C (cid:29) C ) ≈ C , P (1 C (cid:29) C ) ≈ C , P (1 (cid:29) D ) (cid:46) C , P (2 D (cid:29) D ) (cid:46) C ,and P (1 D (cid:29) D ) ≈ C , where C , C , C , C , and C are constants independent of x . Since deg( ε ) = − , ln( x ) ≈ ln( ε − ) / . Hence, P ≈ ( l C + W )2 p ∗ C ln( ε − ) , so that lim x →∞ Pσ ln( ε − ) = l C + W W . (9)This proves equation (20) of the main text.
VIII. ADDITIONAL NUMERICAL CALCULATIONS
In Supplementary Figs. 1 and 2, we consider two discrim-ination mechanisms under the assumptions of equations (21)and (22) of the main text. Supplementary Fig. 1(a) shows agraph for McKeithan’s T-cell receptor mechanism [3], whileSupplementary Fig. 2(a) shows a graph different from boththis and the Hopfield example. We used previously developed,freely-available software [4] to compute the Matrix-Tree for-mula (equation (1) of the main text) for each mechanism, fromwhich we obtained symbolic expressions for P , ε , and σ . Thegraphs in Supplementary Figs. 1(a) and 2(a) have 441 and64 spanning trees rooted at each vertex, respectively, under-scoring the combinatorial complexity which arises away fromequilibrium (main text, Discussion). (If a graph has reversibleedges, so that i → j if, and only if, j → i , which is thecase for all the graphs discussed here, there is a bijection be-tween the sets of spanning trees rooted at any pair of distinctvertices.) Supplementary Figs. 1(b)-(c) and 2(b)-(c) shownumerical plots undertaken in a similar way to those for theHopfield mechanism (main text, Fig. 2), as described in themain text. Similar vertical and diagonal bounds were foundfor the symmetric cases, while similar observations regardingthe asymmetric cases as those made in the main text apply. IX. ASYMPTOTIC RELATION FOR ANON-DISSOCIATION-BASED MECHANISM
We consider a discrimination mechanism having the graphshown in Supplementary Fig. 3(a). Its structure is identicalto that of the Hopfield mechanism (Fig. 1(a) of the main text)but its labels differ to reflect the energy landscape illustrated inSupplementary Fig. 3(b). If the labels are allowable functionswith deg( l (cid:48) D ) = − and deg( l (cid:48) C ) = deg( l C ) = deg( l D ) = 0 ,then, if the mechanism reaches thermodynamic equilibrium, itfollows from equation (9) of the main text that its equilibriumerror fraction satisfies ε eq ∼ x − . (10)If it is further assumed that deg( m (cid:48) D ) = deg( m (cid:48) C ) = − and deg( m D ) = − deg( m C ) = 1 / , while all other labelshave degree , then the mechanism is no longer at equilibrium.Using equations (1) and (3) of the main text, we find that ρ = [( m (cid:48) C + k )( l C + W ) + km C ][( m (cid:48) D + k )( l D + W ) + km D ] ρ C = [ m C ( l (cid:48) C + k (cid:48) ) + k (cid:48) ( l C + W )][( m (cid:48) D + k )( l D + W ) + km D ] ρ C = [ m (cid:48) C ( k (cid:48) + l (cid:48) C ) + kl (cid:48) C ][( k + m (cid:48) D )( l D + W ) + m D k ] ρ D = [( m (cid:48) C + k )( l C + W ) + km C ][ m D ( l (cid:48) D + k (cid:48) ) + k (cid:48) ( l D + W )] ρ D = [ m (cid:48) D ( k (cid:48) + l (cid:48) D ) + kl (cid:48) D ][( k + m (cid:48) C )( l C + W ) + m C k ] . It follows that deg( p ∗ ) = deg( p ∗ C ) = deg( p ∗ C ) = deg( p ∗ D ) = 0 (11)and that deg( p ∗ D ) = − / , (12)Using equations (7) and (14) of the main text along with equa-tions (11) and (12), we can calculate the asymptotics of theterms in the entropy production rate P (equation (12) of themain text), assuming, as in the proof of the Theorem, that weare outside the measure-zero subset of parameter space arising from case 3 of equation (15) of the main text. We find that P (1 (cid:29) C ) ∼ P (2 C (cid:29) C ) ∼ x − / ln( x ) P (1 (cid:29) C ) ∼ P (1 (cid:29) D ) ∼ P (2 D (cid:29) D ) ∼ x − ln( x ) P (1 (cid:29) D ) ∼ x − ln( x ) . It follows that P (cid:45) , so that the entropy production rate isasymptotically constant or vanishes. Furthermore, it can beshown from equations (11) and (12) that deg( ε ) = − / and deg( σ ) = 0 . Hence, the error rate is asymptotically betterthan at equilibrium, for which deg( ε eq ) = − (equation (10)),while the speed remains asymptotically constant. This reflectsa different asymptotic relation to that in equation (19) of themain text. [1] Gunawardena, J. A linear framework for time-scale separa-tion in nonlinear biochemical systems , PLoS ONE , e36321(2012).[2] Hopfield, J. J. Kinetic proofreading: a new mechanism for re-ducing errors in biosynthetic processes requiring high speci-ficity , Proc. Natl. Acad. Sci. USA , 4135 (1974).[3] McKeithan, T. W. Kinetic proofreading in T-cell receptor signaltransduction , Proc. Natl. Acad. Sci. USA , 5042 (1995).[4] Ahsendorf, T., Wong, F., Eils, R. & Gunawardena, J. A frame-work for modelling gene regulation which accommodates non-equilibrium mechanisms , BMC Biol. , 102 (2014).[5] Banerjee, K., Kolomeisky, A. B. and Igoshin, O. A. Elucidatinginterplay of speed and accuracy in biological error correction ,Proc. Natl. Acad. Sci. USA , 5183 (2017).
0 20 40 60 80 100 20 40 60 80 100 (cid:55) (cid:22) (cid:1107) (b) (a) C C C C D D D D
120 120 ln( (cid:1094)(cid:3)(cid:3)(cid:22)(cid:1094)(cid:16)(cid:3)
0 20 40 60 80 100 20 40 60 80 100 (cid:55) (cid:22) (cid:1107) (c)
120 120 ln( (cid:1094)(cid:3)(cid:3)(cid:22)(cid:1094)(cid:16)(cid:3) FIG. 1. Numerics for the T-cell receptor mechanism. (a) Graph for an instance of McKeithan’s T-cell receptor mechanism [3], with labelnames omitted for clarity. (b) Plot of
P/σ against ln( ε /ε ) for approximately points, with the labels satisfying equations (21) and (22) ofthe main text and numerically sampled as described in the main text. The vertical black dashed line corresponds to the bound ε > ε for thismechanism (calculation not shown) that is analogous to equation (5) of the main text for the Hopfield mechanism. The diagonal red dashedline corresponds to equation (23) of the main text, as discussed further there. (c) Similar plot to (b) but with internal discrimination betweencorrect and incorrect substrates, as described in the text, with the light blue points having a lower asymmetry range ( A = 1 ) and the dark bluepoints having a higher range ( A = 5 ). (b) (a) C C D D C D
0 20 40 P / σ
10 30 10 20 30 40 ln( ε /ε) (c)
0 20 40 P / σ
10 30 10 20 30 40 ln( ε /ε) FIG. 2. Numerics for another discrimination mechanism. (a) Graph for a discrimination mechanism that is different from both the Hopfieldand McKeithan mechanisms, with label names omitted for clarity. (b) Points plotted as in Supplementary Fig. 1(b). The vertical black dashedline corresponds to the bound ε > ε for this mechanism (calculation not shown) that is analogous to equation (5) of the main text for theHopfield mechanism. The diagonal red dashed line corresponds to equation (23) of the main text, as discussed further there. (c) Similar plot to(b) but with internal discrimination between correct and incorrect substrates, as described in the text, with the light blue points having a lowerasymmetry range ( A = 1 ) and the dark blue points having a higher range ( A = 5 ). ΔΔ P / σ ln (ε ) energy C C D D (a) (c) energy m’m’m mk’ C k’ D k C k D l’ D l + W C l + W D l’ C C CDD (b)
4 8 12 16 205 1015
Hopfield mechanismnon-dissociation-basedmechanism -1 . .
0 0
FIG. 3. A non-dissociation-based mechanism. (a) Graph with the same structure as that for the Hopfield mechanism (Fig. 1(a) of the main text)but no discrimination between C and D takes place through (cid:29) X , while internal discrimination takes place through X (cid:29) X , as reflectedin the label names. (b) Hypothetical energy landscape for the mechanism shown in (a), illustrating where energy may be expended to drive thesteps with labels m C and m D . (c) Plot of P/σ against ln( ε − ) for a numerical instance of the Hopfield mechanism (Fig. 1(a) of the main text,in blue) and a numerical instance of the non-dissociation-based mechanism in (a) (in red), as x is varied in the range x ∈ [0 , e ] . The numericallabel values have been determined by taking k D = l D = x and l (cid:48) C = l (cid:48) D = x − for the Hopfield mechanism and l (cid:48) D = m (cid:48) D = m (cid:48) C = x − , m D = x / and m C = x − / for the non-dissociation-based mechanism, with all other labels being . label DNAP ribosome (wild type) ribosome (hyperaccurate) ribosome (error-prone) k C
900 0 . .
41 0 . k D
900 47 46 .
002 3 . k (cid:48) C .
001 40 27 37 k (cid:48) D . .
002 36 . m C . .
001 0 .
001 0 . m D . . × − , . − [6 . × − , . − [4 . × − , . − m (cid:48) C
700 25 14 31 m (cid:48) D
700 1 . .
49 3 . l C .
085 0 .
048 0 . l D × − . . . l (cid:48) C
250 0 .
001 0 .
001 0 . l (cid:48) D .
002 [1 . × − , . . . × − , . . . × − , . . W C
250 8 .
415 4 .
752 7 . W D .
012 0 . . . TABLE I. Experimentally measured parameter values, in units of s − , for the Hopfield mechanism in Fig. 1(a) of the main text, shown fordiscrimination during DNA replication by the bacteriophage T7 DNA polymerase (DNAP) and discrimination during mRNA translation bythree E. coli ribosome variants, as annotated. The values were obtained from Tables S1-S4 of [5]. The labels in the first column correspond tothose in Fig. 1(a) of the main text, except that m , m (cid:48) and W now have subscripts C and D , for the correct and incorrect substrates, respectively,to allow for internal discrimination, as explained in the main text. The values of m D and l (cid:48) D were not known for the ribosome variants, sowe chose m D from m C by randomly selecting ln( m D /m C ) from the uniform distribution on [ − , , which is similar to the asymmetryranges of the other parameters, and chose l (cid:48) D to satisfy the external chemical potential constraint used by [5], as explained in footnote [42] ofthe main text. The intervals given for m D and l (cid:48) D indicate the range of sampled values. Some samples have ε < ε0