Universal thermodynamic bounds on symmetry breaking in biochemical systems

Living systems are maintained out-of-equilibrium by external driving forces. At stationarity, they exhibit emergent selection phenomena that break equilibrium symmetries and originate from the expansion of the accessible chemical space due to non-equilibrium conditions. Here, we use the matrix-tree theorem to derive universal thermodynamic bounds on these symmetry-breaking features in biochemical systems. Our bounds are independent of the kinetics and hold for both closed and open reaction networks. We also extend our results to master equations in the chemical space. Using our framework, we recover the thermodynamic constraints in kinetic proofreading. Finally, we show that the contrast of reaction-diffusion patterns can be bounded only by the non-equilibrium driving force. Our results provide a general framework for understanding the role of non-equilibrium conditions in shaping the steady-state properties of biochemical systems.

Introduction.-Livingsystems operate away from thermodynamic equilibrium, continuously harvesting and consuming energy [1].A key feature of out-of-equilibrium living matter is the emergence of various selection phenomena, i.e., biochemical states are populated not only according to their energy (as at equilibrium), but also to their kinetic features [2][3][4][5][6][7].As a consequence, nonequilibrium settings can break equilibrium energetic symmetries, even if the intimate connection between symmetry breaking and dissipation is still under investigation.Only in recent years, a growing amount of works is focusing on highlighting such a link in determined contexts.For example, non-equilibrium conditions have been proven to be crucial to achieve low-error copying mechanisms that preserve DNA and RNA information with an efficiency constrained by the driving force [8][9][10][11].Furthermore, biological systems attain high sensitivity during sensing through non-equilibrium processes [12,13], and the bi-stability of biochemical states, which is crucial for signaling, can be exploited only away from thermodynamic equilibrium [14].Also, the newlydiscovered mechanism of ultra-affinity in chaperone activity is genuinely out-of-equilibrium and bounded by thermodynamic properties [15,16].On a larger scale, spatial symmetry may be broken by the onset of reactiondiffusion patterns, whose emergence is associated with dissipation in the form of energy and information [17,18].
From a broader perspective, the celebrated thermodynamic uncertainty relations set a universal link between dissipation and performance of biological systems [19,20].Indeed, they set a bound for the accuracy of any stochastic current in terms of the average entropy production.However, interesting features of biochemical systems usually are not limited to the statistics of their currents, as for most of the cases mentioned above.
A more general, yet not exhaustive, approach to the quest of understanding the role of dissipation in biochemical systems has been inspired by network theory.Kirch-hoff's law on currents has been applied to study thermodynamic constraints in metabolic networks [21][22][23][24], and the representation of stationary states in terms of spanning trees -presented in [25,26] -has been shown to be a powerful tool to derive thermodynamic and structural constraints on non-equilibrium response [27,28].
In this Letter, we give a unified description of thermodynamic bounds on non-equilibrium symmetry breaking in biochemical systems.By taking advantage of the network representation of the stationary states, we derive universal bounds for the ratio of probabilities of any two states (or set of states), which is a measure of selection and, as such, a quantification of the breaking of equilibrium symmetries.We show that these bounds depend solely on network geometry and equilibrium properties, and notably not on the kinetics, as they can be determined by individual spanning trees (for which detailed balance always holds).The derivation of an upper and lower bound allows us to write down a non-equilibrium phase space, i.e., a region of the entire chemical space in which any feasible stationary solution must lie.These results hold for both closed and open chemical reaction networks (CRN) described in terms of rate equations.In the case of small number of molecules or formation of complexes, our findings can be extended to chemical master equations, providing upper and lower bounds to the ratio of correlation functions of any order.
The universal thermodynamic bounds can be applied to an immense variety of symmetry-breaking mechanisms in biochemical systems.First, we re-derive in a straightforward way the constraints on kinetic proofreading [8,9].Then, we show that the contrast of reaction-diffusion (RD) patterns is limited by the thermodynamic driving force, highlighting a fascinating direct link between observations and non-equilibrium conditions.whose probabilities follow a rate equation: where p i is the probability to find the system in the state i and is equal to the concentration of the species i, c i , divided by the total concentration.This description holds for uni-molecular, catalytic, and auto-catalytic reactions, noticing that the transition rate from j to i might be non-linear, i.e., kij = ω ij (p)k ij , with ω ij (p) = ω ji (p) and p = {p 1 , . . ., p N }.Here, ω ij (p) encodes all the catalytic effects and k ij satisfies generalized detailed balance, the energy difference between the states i and j, and F ij is the non-equilibrium driving force on the edge (ij) of the CRN.At equilibrium, F ij = 0 for all edges, and the stationary solution coincides with the equilibrium Boltzmann distribution, p eq i ∝ e −βEi .Out-of-equilibrium, the steady-state solution, p ss i , depends on the non-linear kinetics and topology of the system, not only on the energies of the states.Notice that the non-equilibrium driving can originate from any source, e.g., chemical potential difference [8,9], thermal gradient [30,31].
An elegant way to write down p ss i is through its diagrammatic representation in terms of spanning trees [25,26].The solution reads as follows: where N is the normalization factor, and A i (T µ ; p ss ) the product of the rates of the µ-th spanning tree, T µ , directed towards the state i, with the sum running over all possible spanning trees.Notice that A i (T µ ; p ss ) generally depends on the solution itself, p ss , as some transition rates might be non-linear.Hence, Eq. ( 2) is closed only for uni-molecular reaction networks but is self-consistent for all other cases governed by Eq. (1).
To gain intuition about the relevance of spanning trees, consider the ratio K eq ij (T µ ) ≡ A i (T µ ; p ss )/A j (T µ ; p ss ).It corresponds to the product of the ratios between forward and backward rates along each edge of the reaction pathway connecting i and j and belonging to spanning tree T µ .Indeed, all the edges outside this reaction pathway cancel out in the ratio.Moreover, for each spanning tree T µ , the reaction pathway connecting any two states is unique.We remark that any ratio of forward and backward transition rates along the same edge is klm / kml = k lm /k ml since the non-linear factor is symmetric.Finally, due to the generalized detailed balance, K eq ij (T µ ) is a pseudo-equilibrium quantity, as it does not depend on the kinetics but only on thermodynamic properties.K eq ij (T µ ) also quantified the ratio of the probabilities of i and j if they were connected only through T µ .In Fig. 1, we derive it in a simple example.The main result of this Letter gives an upper and lower bound on the ratio of probabilities of any two states, i and j, in terms of pseudo-equilibrium quantities [32]: Naively speaking, the maximization (minimization) over all possible spanning trees amounts to maximize (minimize) over all possible reaction pathways connecting i to j.The ratio of steady-state probabilities is thus bounded from both sides by the extremal pseudo-equilibrium values it can take, which are sheer thermodynamic quantities independent from non-linearities of reactions and non-equilibrium kinetics [33].As expected, at equilibrium all reaction pathways give the same contribution due to the fact that detailed balance holds, thus upper and lower bounds are both saturated.This result generalizes similar findings obtained in [34] for a specific case, and in [35] only for linear reaction networks.Below, we show that it can be further extended to sets of states, open CRN, and chemical master equations.In Fig. 2a, we show a pictorial representation of a CRN, where we highlight the spanning tree determining the upper bound on the symmetry breaking between p ss i and p ss j due to non-equilibrium driving, namely T max for simplicity.Fig. 2b presents (in green) the resulting nonequilibrium phase space.At equilibrium, this region collapses into a single line whose slope is determined by the Boltzmann factor e −β∆Eij .Out-of-equilibrium, we have K eq ij (T max ) ≡ K max ij , and analogously we define K min ij .These bounds determine a region in which steady-state probabilities must lie independently of their kinetic features.Noticeably, for a high-dimensional space, the nonequilibrium phase space is a cone whose boundaries are given by bounds on ratios of all pairs of states.Conservation laws, such as mass conservation, shall be introduced as hyperplanes that reduce the accessible region and dimensionality of the phase space.
Eq. ( 3) can be readily generalized to evaluate the selection between two groups of states S = {s 1 , . . ., s N } and Σ = {σ 1 , . . ., σ M }, hence quantifying symmetry breaking between coarse-grained chemical macro-states.For the upper bound, we have [32]: where we (arbitrarily) selected s N and σ M as reference states.The lower bound can be obtained by maximizing the inverse ratio in the l.h.s. of Eq. ( 4).
Open CRN.-In the most general case, chemical networks are connected to external reservoirs that keep the concentrations of some chemicals cm fixed.These couplings lead to an additional term in the rate equation: where c i indicates the concentrations of the species i.Before p i = c i /( j c j ).The spanning tree representation cannot be directly applied, as the chemostatted species are not dynamical variables.However, by merging them into a single moiety ce = m cm , we have that the effective transition rate from e to i through the transition channel m is w im cm /c e [32].Using this identification, the resulting rate equation can be solved with the spanning tree approach, substituting the normalization condition with a constraint on the fixed concentration ce .As a consequence, the bounds shown above can be equally derived for open chemical systems.Kinetic proofreading.-Livingsystems encode information in the form of RNA and DNA.Genome duplica- tion, translation and transcription are processes that use this information to select the correct substrate (e.g., nucleotides, codons) with high fidelity, ensuring the survival of the organism [36].Since the accuracy of this selection is severely limited by equilibrium constraints, back in 1974 Hopfield proposed the existence of energyconsuming intermediate steps that favor kinetic discrimination.The simplest non-equilibrium model for proofreading considers an enzyme, E, that can bind wrong (W ) or right (R) substrates, forming complexes.These can be in a passive (EW and ER) or active (EW * and ER * ) state, e.g., due to ATP hydrolysis.Therefore, the transition from passive to active complexes is driven outof-equilibrium by a chemical potential ∆µ (Fig. 3a).It has been shown that the energetic cost of performing a cycle [8,9] constrains the proofreading performance, i.e., η ≡ p EW * /p ER * ≥ exp (−β( + ∆µ)).In its original formulation, this bound results from an optimization process on the reaction rates.Here, by employing the universal thermodynamic bounds, we have: where η ε = e −βε is the equilibrium error rate due to energetic discrimination (see Fig. 3b).The driving force can either increase or suppress the error rate, depending on the kinetic properties of the network, as shown in Fig. 3c.Our approach recovers the infinite-energy Hopfield rate and can be readily used to find thermodynamic bounds on multi-stage proofreading [32], a process with too many parameters to be bounded via optimization.RD patterns.-Originallyillustrated in the seminal work of Turing [37], the instability of the homogeneous fixed point in reaction-diffusion systems, with the conse-quent emergence of spatial patterns, attracted huge attention in different contexts [38][39][40].RD patterns have also been analyzed from a thermodynamic perspective and their onset classified as a non-equilibrium phase transition [17].Recently, Brauns et al. [41] proposed a phasespace geometric approach to study pattern formation of mass-conserving reaction-diffusion systems.Consider the following dynamics for concentrations u and v: The reaction term f (u, v) contains all the transitions connecting u and v, i.e., f (u, v) , with (i) indicating the reaction channel, while D u and D v are diffusion coefficients.In [41], it is shown that all the values of concentrations explored by a spatially extended pattern must be embedded between the intersection of the reactive nullcline, f (u, v) = 0, and the flux-balance subspace, D u u + D v v = η 0 , where η 0 is a constant determined by the total turnover balance.Mass conservation is a constraint that reduces the degrees of freedom.Here, we study which limits to RD patterns may be set only according to the laws of thermodynamics.
In Fig. 4a, we show the system under investigation.Fig. 4b shows, in red, the maximum and minimum accessible values of u (and v by mass conservation) in the pattern.The universal thermodynamic bounds on the ratios between u and v at stationarity determine the nonequilibrium phase space to which the concentrations must belong.The intersections between these bounds and the flux-balance subspace define the extremal accessible values for u, u * max/min (see Fig. 4b).In [32], we show that: where β∆µ is the non-equilibrium driving force provided by ATP hydrolysis.The rightmost equality is reached when the thermally activated pathway is so slow that the system is dominated by the non-linear rates, while the leftmost equality is attained when At this point, we define the pattern contrast (or visibility) as C x = (x max −x min )/(x max +x min ), with x = u, v From Eq. ( 8), we readily obtain [42]: This inequality immediately tells that pattern formation is a sheer consequence of non-equilibrium conditions, as at equilibrium the non-equilibrium phase-space shrinks into a line and C u goes to zero.Moreover, we only need to know the thermodynamic force driving the system out of equilibrium to bound the contrast of a pattern, independently of the details of the kinetics.In Fig. 4c, we numerically check Eq. (9) in our example.A tighter thermo-kinetic bound for C u can be found by adding information about diffusion coefficients [32].Although we inspected a simple two-state system, the presented approach is valid for any mass-conserving RD systems.Chemical ME.-CRNs with small number of molecules or formation of complexes are impossible to describe using Eq. ( 1).Nevertheless, it is always possible to write down a chemical master equation, where each state is n = {n 1 , . . ., n N }, i.e., the number of molecules of each of the N species.Thus, the solution can be found using the spanning tree method in the entire chemical space.Employing the same methodology leading to Eq.s (3) and (4), we can bound the ratio between correlation functions of any order [32].For l-th order correlations, we have: n a l K eq (j l ,a l ),(j l ,L) K eq (i l ,L),(j l ,L) , (10) where the dependence on T µ has been omitted for clarity.Here, i l = (i 1 , . . ., i l ), a l = (a 1 , . . ., a l ), and K eq (i l ,a l ),(i l ,L) is evaluated from the pathway connecting the states in which i l = a l and i l = L, an (arbitrary) reference state.In [32], we show that Eq. ( 10) can be upper bounded by one single pathway.When l = 1, Eq. ( 10) bounds the ratio of average concentrations, as in Eq. (3).When l = 2, the thermodynamic bound applies to ratio of correlations and might be useful to study how couplings inferred from correlation matrix (e.g., using a maximum entropy approach [43][44][45]) depend on non-equilibrium conditions.Eq. ( 10) is manifestly similar to Eq. ( 4) and the lower bound can be obtained from the upper bound of the inverse ratio.A detailed analysis of applications and criticalities of these universal bounds for Chemical ME is left for future studies, even if they are spiritually identical to the bounds derived above for rate equations.
Discussion.-In this Letter, we find universal thermodynamic bounds on different quantities that might be used to quantify how much equilibrium symmetries are broken by non-equilibrium drivings.We presented two simple systems in which our framework can provide interesting results, in particular highlighting the bound on pattern contrast solely as a function of the thermodynamic force.The general advantage of our findings is that they are independent of non-linear reaction terms and kinetics, usually very important to understand systems operating away from equilibrium.As such, unraveling a reformulation of these bounds in the contexts of information thermodynamics might pave the way to understand the thermodynamic constraints on chemical informationprocessing devices [46,47].
As said above, countless examples of selection phenomena have proven to be relevant in as many different contexts.Surely enough, the generality of our results will be a powerful tool to understand potentialities and limitations of non-equilibrium conditions in shaping the steadystate properties of biochemical systems.Indeed, nonequilibrium drivings expand the accessible chemical space, enabling living systems to adapt to various nonlinear mechanisms and exhibit complex behaviours.For example, the emergence of chiral symmetry breaking in thermodynamically consistent reaction networks is a promising open problem in which the role of a non-equilibrium driving has not been thoroughly explored yet [48].

DERIVATION OF THE UNIVERSAL THERMODYNAMIC BOUNDS
On ratios of steady-state probabilities Any (auto-)catalytic reaction can be absorbed in a non-linear master equation as follows with each pair of transition rates satisfying the generalized detailed balance relation The non-linearity of the rates is due to catalytic reactions and encoded in the prefactor ω ij , i.e., kij = ω ij (p)k ij and kji = ω ij (p)k ji .These non-linear factors have to be the same for forward and backward reaction, ω ij (p) = ω ji (p), so that the general detailed balance relation, Eq. ( 12), is guaranteed.This condition stems from the fact that every chemical reaction must be reversible and ensures thermodynamic consistency.The stationary solution of the non-linear reaction network cannot be explicitly obtained in general.However, the final stationary state have to satisfy, by definition, the following equation: This condition can be rewritten in terms of the spanning trees of the reaction network as follows: In the linear case, i.e., without (auto-)catalytic reactions, the r.h.s. of this equation does not depend on the state of the system and it provides a closed solution for the steady-state probabilities.In the present scenario, the r.h.s.does depend on p ss itself and the equation has to be solved self-consistently, thus being not particularly helpful to find steady-state probabilities.In Fig. 5(a-b), we show the spanning tree decomposition in a particularly simple case.Fig. 5(c) shows the decomposition between thermodynamic and kinetic (non-linear) parts, while Fig. 5(d) highlights that the symmetric features of ω ij makes the ratio of A i and A j , for any i and j, independent from p ss .We can make use of Eq. ( 14) to express any ratio of probabilities: Noting that the following mathematical inequalities always hold (it can be proven by induction, see also Fig. 6): the bounds presented in the main text immediately follow by applying Eq. ( 16) to Eq. ( 15).Indeed, as explained in the main text and shown in Fig. 5(d) above, the ratio between any A i and A j cannot depend on p ss .

On ratios of coarse-grained steady-state probabilities
Sometimes we only have access to the concentrations of coarse-grained chemical states, i.e., a set of states impossible to resolve from each other.The universal thermodynamic bounds on the selection parameter can be obtained even in this case.To this aim, consider two sets of states, S and Σ.We arbitrarily choose a reference state, say s N for S and σ M for Σ.Then, we compute all the ratios within each set with respect to this state:  14).(c) Decomposition into thermodynamic and kinetic part of a specific tree.(d) Proof of the fact that any ratio between Ai and Aj does not depend on the kinetic term.FIG. 6. Visualization of the inequality in Eq. ( 16), i.e., mini(ai/bi) ≤ i ai/ i bi ≤ maxi(ai/bi), which is used to prove the universal thermodynamic bounds presented in the main text.
Noting again that all ratios do not depend on steady probabilities, we can apply the inequality in Eq. ( 16) and readily obtain the upper bound presented in the main text.The lower bound can be analogously derived.Despite a more complex structure with respect to the previous case, these bounds again depend solely on the thermodynamic forces acting on the underlying reaction networks.

BOUNDS FOR OPEN CHEMICAL REACTION NETWORKS
Here, we present the derivation of the bound in the case of open CRNs.The rate equation for these systems is: Here, cm are the chemostatted species whose concentrations are kept constant by external reservoirs.In Fig. 7(a), we present a simple example of this setting.Of course, the spanning tree method cannot be straightforwardly applied, as cm are not dynamical variables, being fixed in time by external constraints.
As graphically shown in Fig. 7(b), we can merge all cm into a fictitious chemical species of concentration ce = m cm .The resulting rate equation in terms of this reduced set of states is: Now, each m identifies a different reaction channel and the spanning tree method can be applied just by changing the normalization condition.Notice that, instead of fixing one of the concentrations so that the sum of them will be c tot , the total concentration, as for closed CRNs, here we have to fix ce to be equal the value externally imposed by chemostats.Thus, the universal thermodynamic bounds can be derived as shown above.In Fig. 7(c-d), we show how to identify trees and reaction pathways, highlighting that each reaction channel has to be considered separately.

KINETIC PROOFREADING
In this section, we discuss how to reach the thermodynamic bounds for proofreading schemes.For convenience, we set β = 1 in the following discussion.In Fig. 8a, we illustrate a basic one-step proofreading network, also presented in the main text.At thermodynamic equilibrium, i.e. ∆µ = ln k1k2k3 k−1k−2k−3 = = 0, there is no proofreading and the two states EW * and ER * are discriminated by their energy difference so that the error rate is η ε = e − .Hopfield showed that better discrimination can be achieved by adding a proofreading step and driving the system out of equilibrium [36].The Hopfield error rate, η h , is the square of the equilibrium error rate, i.e., η h = e −2 .It is easy to show that the Hopfield error rate can be recovered employing a kinetic symmetry constrain k i = k i for all i, and considering the infinitely far-from-equilibrium limit ∆µ → ∞, i.e., an infinite amount of available energy.In the most general case, the thermodynamic bounds are shown in Fig. 8a and reported in the main text.To saturate these bounds, the kinetic symmetry constraint is loosened to allow a dominating reaction pathway to be kinetically favored.This leads to a significant kinetic discrimination effect when the transition from E and EW is negligible.
The dominating pathway argument can be generalized to analyze the thermodynamic cost of discrimination in networks with arbitrary topology.The most efficient discrimination pathway can be identified by comparing the thermodynamic driving forces along different pathways.One simple example is an N-step proofreading network as illustrated in Fig. 8b, which contains n intermediate states before reaching the final releasing states EW * or ER * .Energy is injected to drive the switching between the intermediate states.The pathways that maximize (minimize) the driving force between the two final states determines the upper (lower) universal thermodynamic bound in Fig. 8b.The minimal error is reached by choosing kinetic parameters so that the kinetic rates on the chemical pathways in gray (see Fig. 8b) are negligible.This choice allows us to reach the optimal driving force from EW * to ER * .

Bounds on visibility of reaction-diffusion pattern
Reaction-diffusion pattern originate from chemical reactions with bi-stability.For mass-conserving reaction-diffusion systems, a phase-space geometry approach can be applied to study various properties of reaction-diffusion pattern [41].Our thermodynamic bounds determine the accessible phase-space as a function of the available energy, and thus set limits on such properties.To illustrate this, we consider a simple reaction-diffusion system with two species U and V .The boundaries of the non-equilibrium phase space are determined by two reaction pathways.The autocatalytic reaction scheme shapes the reactive nullclines, enabling switching between these pathways.(c) Stationary state u and v patterns.The fast-diffusion limit of species V results in a uniform v-pattern, while the bi-stable nature of the autocatalytic reaction creates the two plateaus of u-pattern.

Non-equilibrium phase space
An active (high-order) autocatalytic reaction is one of the simplest chemical reaction models that can show multistability.In such a reaction system, we consider two reaction pathways, one describing spontaneous transitions and the other actively fueling the autocatalytic reaction.The reactions are schematized as follows (and shown in Fig. 9a): Here the two reaction pathways lead to different equilibrium states.As illustrated in Fig. 9b, the presence of both branches allows the system to evolve towards a non-equilibrium stationary state.The non-linear term controls the switch between the reaction branches.When the total concentration is small, the catalytic reaction is negligible and the system mainly stays in the equilibrium state determined by the spontaneous branch.In the other limit, with high concentration, the catalytic branch dominates and the system approaches the (effective) equilibrium state compatible with the energy balance of the coupled ATP hydrolysis.Notice that here we are ignoring the presence of the nucleotide-exchange reactions to avoid complications that will not change the main message.
This system can be easily written in the form of a non-linear rate equation with two distinct reaction pathways between U and V .Then, the complete reaction-diffusion equations for the concentrations u and v are with the effective rates k( 1) . These rates satisfy the generalized detailed balance condition in Eq. ( 12).Here, the linear reaction pathway describes spontaneous transitions and is characterized by the following equilibrium constant: where ∆E uv is the energy difference between the two states.In contrast, the active autocatalytic branch is driven by ATP hydrolysis with a driving force ∆µ = ∆µ ), thus its equilibrium constant is: The system is globally out of equilibrium due to the presence of a non-zero thermodynamic affinity along the cycle involving both reaction branches, that is A = ln(K vu /K (2) vu ) = β∆µ.Clearly, A is proportional to the available energy from ATP hydrolysis, encoded in the driving force ∆µ.
Due to the unbalance between the two reaction pathways (Eq.s ( 21) and ( 22)), when the system is driven out of equilibrium, the non-equilibrium phase space expands in the (u, v) plane (see Fig. 9b).Notice that in equilibrium conditions, i.e., ∆µ = 0, this collapses into a line with a slope dictated by the energy difference between the two states.The concentration-dependent nature of the autocatalytic reaction makes reactive nullclines to have a non-trivial s-shape in this plane represented by the blue line in Fig. 9b.This allows for multiple intersections with the fluxbalance subspace (dashed orange line in Fig. 9b).They exactly determine the range of possible concentrations in a stationary pattern.Nevertheless, the equilibrium constants associated with each branch can provide a bound for such concentrations that is reached when the flux-balance subspace is horizontal and the reactive nullcline is attached to the boundaries of the non-equilibrium phase space.In particular, the intersections between the flux-balance subspace and the phase space boundaries give the maximum and minimum u, respectively u * max and u * min .Their ratio is upper bounded by the available energy, i.e., u * max /u * min ≤ e β∆µ , which is saturated only for a horizontal flux-balance subspace.The derivation simply follows from the geometric construction and considering the flux-balance subspace can only have a negative slope as both D u and D v have to be positive.The upper bound of the visibility immediately follows: where B th is the name for this thermodynamic bound.To verify this bound, we did numerical simulations of 1D diffusion-reaction system and plotted the results in the Fig. 4c of the main text.In the simulation, we use the 3-rd order autocatalytic reaction with the parameters

Thermo-kinetic tighter bound
The thermodynamic bound is universal and independent of the kinetics of the system.However, we can get a tighter bound, which we call the thermo-kinetic bound, by including information about the kinetics.The bound on contrast comes from the range of concentration defined by the intersections between the flux-balance subspace, η 0 = D u u + D v v, and the two boundaries of the non-equilibrium phase space, u/v = K max/min vu .Instead of providing an upper limit for the ratio u * max /u * min , we can explicitly solve these intersections as follows: where R D = D u /D v and η = η 0 /D v .The solution give a thermo-kinetic bound on pattern contrast as where γ u = ln . It is evident that γ u is positive and approaches zero in R D → 0 limit.In Fig. 10a, we show several 2D pattern with different values of model parameters and in Fig. 10b, we show that the thermo-kinetic bound provides a better estimate of the actual visibility of reaction-diffusion patterns.Therefore the thermodynamic bound can be recovered when the diffusion of V is infinitely faster than the diffusion of U : This limit agrees with our previous analysis that the maximum possible contrast of a pattern is saturated when the other species is in a uniform distribution due to fast diffusion.Analogously, Eq. ( 23) is saturated and the thermodynamic bound, B th , coincides with B th−kin when the flux-balance subspace is a horizontal line.

BOUNDS FOR CHEMICAL MASTER EQUATION
Every chemical system can be written in the form of a chemical master equation (CME).The state of a CME is defined by the number of molecules of each chemical species.For example, if we have N chemical species, a generic state is n = {n 1 , n 2 , . . ., n N } with constraints coming from the fixed number of molecules per species.For example, if only the total number of molecules is fixed to M tot , then i n i = M tot .Conversely, if A and B can form a complex but cannot convert one into another, two constraints have to be put on the total number of molecules in both species.At any rate, we have a CME defined in terms the probability to be in a given state n, i.e., p( n), plus the normalization condition.A CME cannot always be written as a rate equation with only catalytic reactions, e.g., when complexes are involved and when we cannot take the limit of large number of molecules.However, for non-complex species, it is reasonable to define the average probability to observe a given state, k, as follows: In the absence of complexes, P(k) satisfies a rate equation in the limit of a very high number of molecules, otherwise fluctuations play an important role.For a multi-molecular complex, the definition of an analogous quantity involves the use a normalization factor at the denominator that might be different from M tot , depending on whether or not the interest is to quantity the average amount of complexes with respect to total monomers.

Bounds on the ratio of averages
We can define a symmetry-breaking index as we did before: where the notation n/n i1 indicates that the sum is performed over all indices except for n i1 .Notice that the space of possible states is way larger than before, since we are tracking the number of molecules in each chemical species.Since we are now dealing with monomers in quantifying the symmetry breaking, both n i1 and n j1 can be at least 0 and at most M tot .First, we define i m = (i 1 , . . ., i m ) as an m-dimensional vector whose values range from 1 to N , indicating the species, and a m = (a 1 , . . ., a m ) as an m-dimensional vector whose elements span all possible number of molecules (per species) compatible with the chemical space.Thus, we introduce the following short notation: where i1,a1 n/ni 1 indicates a summation over all n but n i1 , which is fixed to a 1 instead.Moreover, since p( n) is governed by a Master Equation, its steady-state solution can be expressed in terms of spanning trees, as before.Thus, we have: where A a1 i1 (T µ ) is a short notation indicating the product of all the rates belonging to the spanning tree T µ , oriented towards the state n where n i1 = a 1 .Applying the inequality shown above by maximizing only over T µ , we obtain: s ≤ max {Tµ} n a 1 K eq (i1,a1),(i1,L) (T µ ) n a 1 K eq (j1,a1),(j1,L) (T µ ) K eq (i1,L),(j1,L) (T µ ) , by indicating with K eq (i1,a1),(i1,L) (T µ ) = A a1 i1 (T µ )/A L i1 (T µ ) and selecting a 1 = L as a reference state for both i 1 and j 1 for simplicity.The derivation follows the same steps shown above to derive the universal thermodynamic bounds for set of states.Clearly, we can maximize also over a 1 , obtaining: s ≤ max {Tµ},a1 n/ni 1 K eq (i1,a1),(i1,L) (T µ ) n/nj 1 K eq (j1,a1),(j1,L) (T µ ) K eq (i1,L),(j1,L) (T µ ) , (32) Results.-Consider a chemical network of N species arXiv:2212.12074v2[cond-mat.stat-mech]21 May 2023

FIG. 1 .
FIG. 1.(a) A four-state reaction network.(b) All spanning trees, and the directed tree to state 1 constructed from the third spanning tree.(c) Steady-state representation using directed trees.(d) The ratio of two directed trees (to state 1 and 4) obtained from the same spanning tree gives a unique reaction pathway connecting these states.

FIG. 3 .
FIG. 3. (a) Proofreading network driven out of equilibrium by a chemical potential ∆µ associated with a fuel-to-waste conversion.(b) Spanning trees defining the upper and lower bound of the error rate, η ≡ pEW * /pER * .(c) Blue points are the error rates of the proofreading network with random kinetic parameters but the same universal thermodynamic bounds (limiting red and blue solid lines).The dashed-yellow line shows the error rate under equilibrium discrimination.
FIG.4.Thermodynamic bound on the contrast of RD patterns.(a) Model system under investigation, driven out of equilibrium by ATP hydrolysis.(b) Universal thermodynamic bounds (in blue) determine the non-equilibrium phase space and contain the actual range of the observed pattern (in red).(c) Pattern visibility is upper bounded by Eq. (9), two corresponding 1D patterns are plotted in the subpanel.

FIG. 5 .
FIG. 5. (a) Spanning tree decomposition for a simple network.(b) Contribution to p ss 1 , i.e., numerator of the r.h.s. of Eq. (14).(c) Decomposition into thermodynamic and kinetic part of a specific tree.(d) Proof of the fact that any ratio between Ai and Aj does not depend on the kinetic term.

1 FIG. 7 .
FIG. 7. Network decomposition of an open chemical reaction network.(a) An open chemical reaction network with two chemostatted species.(b) The two chemostatted species can be merged into a single moiety with concentration ce = c3 + c4.(c) Spanning tree decomposition of the open chemical reaction network.(d) The reaction pathways between state 1 and state 2 are identified from the spanning trees.
FIG.8.Proofreading networks and the corresponding thermodynamic bounds.

Flux-balance
FIG. 9. (a) Active autocatalytic model.(b)The boundaries of the non-equilibrium phase space are determined by two reaction pathways.The autocatalytic reaction scheme shapes the reactive nullclines, enabling switching between these pathways.(c) Stationary state u and v patterns.The fast-diffusion limit of species V results in a uniform v-pattern, while the bi-stable nature of the autocatalytic reaction creates the two plateaus of u-pattern.
25 and system size L = 50.The pattern in Fig.9cis obtained with the same parameters and β∆µ = 2.6.

FIG. 10 .
FIG. 10.(a) 2D reaction-diffusion patterns under varying driving forces.(b) Thermodynamic and thermo-kinetic bounds on the visibility of the 2D pattern.Adding information about the kinetics clearly improves the estimate.Here, Du = 0.1, Dv = 5, systems size 100 × 100, and all the other kinetic parameters are the same as in Fig. 9.