Nonequilibrium Thermodynamics of Chemical Reaction Networks: Wisdom from Stochastic Thermodynamics

We build a rigorous nonequilibrium thermodynamic description for open chemical reaction networks of elementary reactions. Their dynamics is described by deterministic rate equations satisfying mass action law. Our most general framework considers open networks driven by time-dependent chemostats. The energy and entropy balances are established and a nonequilibrium Gibbs free energy is introduced. The difference between this latter and its equilibrium form represents the minimal work done by the chemostats to bring the network in its nonequilibrium state. It is minimized in nondriven detailed-balanced networks (i.e. networks which relax to equilibrium states) and has an interesting information-theoretic interpretation. We further show that the entropy production of complex balanced networks (i.e. networks which relax to special kinds of nonequilibrium steady states) splits into two non-negative contributions. One charaterizing the dissipation of the nonequilibrium steady state and the other the transients due to relaxation and driving. Our theory lays the path to study time-dependent energy and information transduction in biochemical networks.


I. INTRODUCTION
Thermodynamics of chemical reactions has a long history. The second half of the 19th century witnessed the dawn of the modern studies on thermodynamics of chemical mixtures. It is indeed at that time that J. W. Gibbs introduced the concept of chemical potential and used it to define the thermodynamic potentials of non-interacting mixtures [1]. Several decades later, this enabled T. de Donder to approach the study of chemical reacting mixtures from a thermodynamic standpoint. He proposed the concept of affinity to characterize the chemical force irreversibly driving chemical reactions and related it to the thermodynamic properties of mixtures established by Gibbs [2]. I. Prigogine, who perpetuated the Brussels School founded by de Donder, introduced the assumption of local equilibrium to describe irreversible processes in terms of equilibrium quantities [3,4]. In doing so, he pioneered the connections between thermodynamics and kinetics of chemical reacting mixtures [5].
During the second half of the 20th century, part of the attention moved to systems with small particle numbers which are ill-described by "deterministic" rate equations. The Brussels School, as well as other groups, produced various studies on the nonequilibrium thermodynamics of chemical systems [6][7][8][9][10][11] using a stochastic description based on the (Chemical) Master Equation [12,13]. These studies played an important role during the first decade of the 21st century for the development of Stochastic Thermodynamics, a theory that systematically establishes a nonequilibrium thermodynamic description for systems obeying stochastic dynamics [14][15][16][17], including chemical reaction networks (CRNs) [18][19][20][21][22].
Another significant part of the attention moved to the thermodynamic description of biochemical reactions in terms of deterministic rate equations [23,24]. This is not so surprising since living systems are the paramount example of nonequilibrium processes and they are powered by chemical reactions. The fact that metabolic processes involve thousands of coupled reactions also emphasized the importance of a network description [25][26][27]. While complex dynamical behaviors such as oscillations were analyzed in small CRNs [28,29], most studies on large biochemical networks focused on the steady-state dynamics. Very few studies considered the thermodynamic properties of CRNs [30][31][32][33]. One of the first nonequilibrium thermodynamic description of open biochemical networks was proposed in Ref. [34]. However, it did not take advantage of Chemical Reaction Network Theory which connects the network topology to its dynamical behavior and which was extensively studied by mathematicians during the seventies [35][36][37] (this theory was also later extended to stochastic dynamics [38][39][40][41]). As far as we know, the first and single study that related the nonequilibrium thermodynamics of CRNs to their topology is Ref. [22], still restricting itself to steady states.
In this paper, we consider the most general setting for the study of CRNs, namely open networks driven by chemostatted concentrations which may change over time. To the best of our knowledge, this was never considered before. In this way, steady-state properties as well as transient ones are captured. Hence, in the same way that Stochastic Thermodynamics is built on top of stochastic dynamics, we systematically build a nonequilibrium thermodynamic description of CRNs on top of deterministic chemical rate equations. In doing so, we establish the energy and entropy balance and introduce the nonequilibrium entropy of the CRN as well as its nonequilibrium Gibbs free energy. We show this latter to bear an information-theoretical interpretation similar to that of Stochastic Thermodynamics [42][43][44][45] and to be related to the dynamical potentials derived by mathemati-cians. We also show the relation between the minimal chemical work necessary to manipulate the CRNs far from equilibrium and the nonequilibrium Gibbs free energy. Our theory embeds both the Prigoginian approach to thermodynamics of irreversible processes [5] and the thermodynamics of biochemical reactions [23]. Making full use of the mathematical Chemical Reaction Network Theory, we further analyze the thermodynamic behavior of two important classes of CRNs: detailed-balanced networks and complex-balanced networks. In absence of time-dependent driving, the former converges to thermodynamic equilibrium by minimizing their nonequilibrium Gibbs free energy. In contrast, the latter converges to a specific class of nonequilibrium steady states and always allow for an adiabatic-nonadiabatic separation of their entropy production, which is analogous to that found in Stochastic Thermodynamics [46][47][48][49][50]. We note that while finalizing this paper, a result similar to the latter was independently found in Ref. [51].

Outline and Notation
The paper is organized as follows. After introducing the necessary concepts in chemical kinetics and Chemical Reaction Network Theory, sec. II, the nonequilibrium thermodynamic description is established, sec. III. As in Stochastic Thermodynamics, we build it on top of the dynamics and formulate the entropy and energy balance, § III D and § III E. Chemical work and nonequilibrium Gibbs free energy are also defined and the informationtheoretic content of the latter is discussed. The special properties of detailed-balance and of complex-balanced networks are considered in sec. V and IV, respectively. Conclusions and perspectives are drawn in sec. VI, while some technical derivations are detailed in the appendices. We now proceed by fixing the notation. We consider a system composed of reacting chemical species X σ , each of which is identified by an index σ ∈ S, where S is the set of all indices/species. The species populations change due to elementary reactions, i.e. all reacting species and reactions must be resolved (none can be hidden), and all reactions must be reversible, i.e. each forward reaction +ρ has a corresponding backward reaction −ρ. Each pair of forward-backward reactions is a reaction pathway denoted by ρ ∈ R. The orientation of the set of reaction pathways R is arbitrary. Hence, a generic CRN is represented as The constants k +ρ (k −ρ ) are the rate constants of the forward (backward) reactions. The stoichiometric coefficients −∇ σ +ρ and ∇ σ −ρ identify the number of molecules of X σ involved in each forward reaction +ρ (the stoichiometric coefficients of the backward reactions have opposite signs). Once stacked into two non-negative ma- FIG. 1: Representation of a closed CRN. The chemical species are {Xa, · · · , Xe}. The two reaction pathways are labeled by 1 and 2. The nonzero stoichiometric coefficients are −∇ a +1 = −1, ∇ b −1 = 2 and ∇ c −1 = 1 for the first forward reaction and −∇ c +2 = −1, −∇ d +2 = −1 and ∇ e −2 = 1 for the second one. Since the network is closed, no chemical species is exchanged with the environment.

integer-valued stoichiometric matrix
The reason for the choice of the symbol "∇" will become clear later.

Example 1.
The stoichiometric matrix of the CRN depicted in fig. 1 is Physical quantities associated to species and reactions are represented in upper-lower indices vectorial notation. Upper and lower indexed quantities have the same physical values, e.g. Z i = Z i , ∀ i. We use the Einstein summation notation: repeated upper-lower indices implies the summation over all the allowed values for those indices-e.g. σ ∈ S for species and ρ ∈ R for reactions. Given two arbitrary vectorial quantities a = {a i } and b = {b i }, the following notation will be used Finally, given the matrix C, whose elements are {C i j }, the elements of the transposed matrix C T are {C j i }. The time derivative of a physical quantity A is denoted by d t A, its steady state value by an overbar,Ā, and its equilibrium value by A eq or A eq . We reserve the overdot, A, to denote the rate of change of quantities which are not exact time derivatives.

II. DYNAMICS OF CRNS
In this section, we formulate the mathematical description of CRNs [52,53] in a suitable way for a thermodynamic analysis. We introduce closed and open CRNs and show how to drive these latter in a time-dependent way. We then define conservation laws and cycles and review the dynamical properties of two important classes of CRNs: detailed-balanced networks and complex-balanced networks.
We consider a chemical system in which the reacting species {X σ } are part of a homogeneous and ideal dilute solution: the reactions proceed slowly compared to diffusion and the solvent is much more abundant than the reacting species. Temperature T and pressure p are kept constant. Since the volume of the solution V is overwhelmingly dominated by the solvent, it is assumed constant. The species abundances are large enough so that the molecules discreteness can be neglected. Thus, at any time t, the system state is well-described by the molar concentration distribution The reaction kinetics is controlled by the reaction rate functions J ±ρ {Z σ } , which measure the rate of occurrence of reactions and satisfy the mass action kinetics [52,54,55] The net concentration current along a reaction pathway ρ is thus given by Example 2. For the CRN in fig. 1 the currents are

A. Closed CRNs
A closed CRN does not exchange any chemical species with the environment. Hence, the species concentrations vary solely due to chemical reactions and satisfy the rate equations Since rate equations are nonlinear, complex dynamical behaviors may emerge [29]. The fact that the rate equations (7) can be thought of as a continuity equation for the concentration, where the stoichiometric matrix ∇ (2) acts as a discrete differential operator, explains the choice of the symbol "∇" for the stoichiometric matrix [56].

B. Driven CRNs
In open CRNs, matter is exchanged with the environment via reservoirs which control the concentrations of some specific species, fig. 2. These externally controlled species are said to be chemostatted, while the reservoirs controlling them are called chemostats. The chemostatting procedure may mimic various types of controls by the environment. For instance, a direct control could be implemented via external reactions (not belonging to the CRN) or via abundant species whose concentrations are negligibly affected by the CRN reactions within relevant time scales. An indirect control may be achieved via semipermeable membranes or by controlled injection of chemicals in continuous stirred-tank reactors.
Among the chemical species, the chemostatted ones are denoted by the indices σ y ∈ S y , and the internal ones by σ x ∈ S x (S ≡ S x ∪S y ). Also, the part of the stoichiometric matrix related to the internal (resp. chemostatted) species is denoted by ∇ X = ∇ σx ρ (resp. ∇ Y = ∇ In nondriven open CRNs the chemostatted species have constant concentrations, i.e. {d t Z σy = 0}. In driven open CRNs the chemostatted concentrations change over time according to some time-dependent protocol π(t): {Z σy ≡ Z σy (π(t))}. The changes of the internal species are solely due to reactions and satisfy the rate equations Instead, the changes of chemostatted species {d t Z σy } are not only given by the species formation rates {∇ σy ρ J ρ } but must in addition contain the external currents {I σy }, which quantify the rate at which chemostatted species enter into the CRN (negative if chemostatted species leave the CRN), This latter equation is not a differential equation since the chemostatted concentrations {Z σy } are not dynamical variables. It shows that the external control of the chemostatted concentration is not necessarily direct, via the chemostatted concentrations, but can also be indirectly controlled via the external currents. We note that (10) is the dynamical expression of the decomposition of changes of species populations in internal-external introduced by de Donder (see [ for given chemostatted concentrations {Z σy }.

C. Conservation Laws
In a closed CRN, a conservation law = { σ } is a left null eigenvector of the stoichiometric matrix ∇ [23,25] Conservation laws identify conserved quantities L ≡ σ Z σ , called components [23,25], which satisfy We denote a set of independent conservation laws of the closed network by λ and the corresponding components by L λ ≡ λ σ Z σ . The choice of this set is not unique, and different choices have different physical meanings. This set is never empty since the total mass is always conserved. Physically, conservation laws are often related to parts of molecules, called moieties [58], which are exchanged between different species and/or subject to isomerization (see example 4).
In an open CRN, since only {Z σx } are dynamical variables, the conservation laws become the left null eigenvectors of the stoichiometric matrix of the internal species ∇ X . Stated differently, when starting from the closed CRN, the chemostatting procedure may break a subset of the conservation laws of the closed network { λ } [56]. E.g. when the first chemostat is introduced the total mass conservation law is always broken. Within the set λ , we label the broken ones by λ b and the unbroken ones by λ u . The broken conservation laws are characterized by system environment k +1 k -1 where the first term is nonvanishing for at least one ρ ∈ R.
The broken components {L λ b ≡ λ b σ Z σ } are no longer constant over time. On the other hand, the unbroken conservation laws are characterized by where the first term vanishes for all ρ ∈ R. Therefore, the unbroken components {L λu ≡ λu σ Z σ } remain constant over time. Without loss of generality, we choose the set { λ } such that the entries related to the chemostatted species vanish, λu σy = 0, ∀λ u , σ y .
When chemostatting as in fig. 2, the first two conservation laws break while the last one remains unbroken. We also note that this set is chosen so that the unbroken conservation law satisfies 3 a = 3 e = 0. When considering the specific implementation in fig. 3 of the CRN in fig. 2, we see that the first two conservation laws in (16) represent the conservation of the concentrations of the moiety H and C, respectively. Instead, the third conservation law in (16) does not have a straightforward interpretation. It is related to the fact that when the species H or C are produced also O must be produced and vice versa.

D. Detailed-Balanced Networks
A steady state (11) is said to be an equilibrium state, {Z σ eq }, if it satisfies the detailed balance property [57, § 9.4], i.e. all concentration currents (5) vanish, For open networks, this means that the external currents, eq. (11b), must also vanish {I σy eq = 0}. By virtue of mass action kinetics, eq. (4), the detailed balance property (17) can be rewritten as A CRN is said to be detailed-balanced if, for given kinetics {k ±ρ } and chemostatting {Z σy }, its dynamics exhibits an equilibrium steady state (17). For each set of unbroken components {L λu }-which are given by the initial condition and constrain the space where the dynamics dwells-the equilibrium distribution is globally stable [59]. Equivalently, detailed-balanced networks always relax to an equilibrium state, which for a given kinetics and chemostatting is unique and depends on the unbroken components only, see also sec. V.
Closed CRNs must be detailed-balanced. This statement can be seen as the zeroth law for CRNs. Consequently, rather than considering eq. (18) as a property of the equilibrium distribution, we impose it as a property that the rate constants must satisfy and call it local detailed-balance property. It is a universal property of elementary reactions which holds regardless of the network state. Indeed, while the equilibrium distribution depends on the components, the r.h.s. of eq. (18) does not. This point will become explicit after introducing the thermodynamic structure, eq. (88) in sec.V. The local detailed-balance property will be rewritten in a thermodynamic form in § III B, eq. (50).
In open nondriven CRNs, the chemostatting procedure may prevent the system from reaching an equilibrium state. To express this scenario algebraically we now introduce the concepts of emergent cycle and cycle affinity.
A cyclec = {c ρ } is a right null eigenvector of the stoichiometric matrix [56], namely Since ∇ is integer-valued,c can always be rescaled to only contain integer coefficients. In this representation, its entries denote the number of times each reaction occurs (negative signs identify reactions occurring in backward direction) along a transformation which overall leaves the concentration distributions {Z σ } unchanged, see example 5. We denote by {c α } a set of linearly independent cycles. An emergent cycle c = {c ρ } is defined algebraically as [56] ∇ σx In its integer-valued representation, the entries of c denote the number of times each reaction occurs along a transformation which overall leaves the concentrations of the internal species {Z σx } unchanged while changing the concentrations of the chemostatted species by an amount ∇ σy ρ c ρ . These latter are however immediately restored to their prior values due to the injection of −∇ σy ρ c ρ molecules of X σy performed by the chemostats. Emergent cycles are thus pathways transferring chemicals across chemostats while leaving the internal state of the CRN unchanged. We denote by {c ε } a set of linearly independent emergent cycles.
When chemostatting an initially closed CRN, for each species which is chemostatted, either a conservation law breaks-as mentioned in § II C-or an independent emergent cycle arises [56]. This follows from the rank nullity theorem for the stoichiometric matrices ∇ and ∇ X , which ensures that the number of chemostatted species |S y | equals the number of broken conservation laws |λ b | plus the number of independent emergent cycles |ε|: |S y | = |λ b | + |ε|. Importantly, the rise of emergent cycles is a topological feature: it depends on the species which are chemostatted, but not on the chemostatted concentrations. We also note that emergent cycles are modeled as "flux modes" in the context of metabolic networks [60][61][62].

Example 5.
To illustrate the concepts of cycles and emergent cycles, we use the following CRN [56] whose Y 1 and Y 2 species are chemostatted. The stoichiometric matrix decomposes as The set of linearly independent cycles (19) consists of only one cycle, which can be written as As the CRN is chemostatted one linearly independent emergent cycle (20) arises We now see that if each reaction occurs a number of times given by the entry of the cycle (23), the CRN goes back to the initial state, no matter which one it is. On the other hand, when the emergent cycle (24) is performed, the state of the internal species does not change, while two molecules of Y 1 are annihilated and two of Y 2 are created. However, since the chemostats restore their initial values, the overall result of c is to transfer two Y 1 , transformed in Y 2 , from the first to the second chemostat.
The closed version of this CRN has two independent conservation laws, the first of which, 1 , is broken following the chemostatting of any of the two species Y 1 or Y 2 . The other chemostatted species, instead, gives rise to the emergent cycle (24), so that the relationship |S y | = |λ b | + |ε| is satisfied.
Any cyclec α and emergent cycle c ε bears a cycle affinity From the definition of cycle (19) and current (5), and the local detailed balance (18), it follows that the cycle affinities along the cycles (19) vanish, {Ã α = 0}, and that the cycle affinities along the emergent cycles only depend on the chemostatted concentrations Since emergent cycles are pathways connecting different chemostats, the emergent affinities quantify the chemical forces acting along the cycles. This point will become clearer later, when the thermodynamic expressions of the emergent cycle affinities {A ε } will be given, eq. (49). A CRN is detailed-balanced if and only if all the emergent cycle affinities {A ε } vanish. This condition is equivalent to the Wegscheider's condition [59]. This happens when the chemostatted concentrations fit an equilibrium distribution. As a special case, unconditionally detailedbalanced networks are open CRNs with no emergent cycle. Therefore, they are detailed-balanced for any choice of the chemostatted concentrations. Consequently, even when a time-dependent driving acts on such a CRN and prevents it from reaching an equilibrium state, a well-defined equilibrium state exists at any time: the equilibrium state to which the CRN would relax if the time-dependent driving were stopped. Example 6. Any CRN with one chemostatted species only (|S y | = 1) is unconditionally detailed-balanced. Indeed, as mentioned in § II C, the first chemostatted species always breaks the mass conservation law, |λ b | = 1 and thus no emergent cycle arises, The open CRN in fig. 2 is an example of unconditionally detailed-balanced network with two chemostatted species, since the chemostatting breaks two conservation laws, see example 4. Indeed, a nonequilibrium steady state would require a continuous injection of Y a and ejection of Y e (or vice versa). But this would necessary result in a continuous production of X b and consumption of X d which is in contradiction with the steady state assumption.
Finally, a tacit assumption in the above discussion is that the network involves a finite number of species and reactions, i.e. the CRN is finite-dimensional. Infinitedimensional CRNs can exhibits long-time behaviors different from equilibrium even in absence of emergent cycles [63].

E. Complex-Balanced Networks
To discuss complex-balanced networks and complexbalanced distributions, we first introduce the notion of complex in open CRNs.
A complex is a group of species which combines in a reaction as products or as reactants. Each side of eq. (1) defines a complex but different reactions might involve the same complex. We label complexes by γ ∈ C, where C is the set of complexes.
The set of complexes is and the complex 2X b is involved in both the second and third reaction.
The notion of complex allows us to decompose the stoichiometric matrix ∇ as We call Γ = {Γ σ γ } the composition matrix [35,37]. Its entries Γ σ γ are the stoichiometric number of species X σ in the complex γ. The composition matrix encodes the structure of each complex in terms of species, see example 8. The matrix ∂ = {∂ γ ρ } denotes the incidence matrix of the CRN, whose entries are given by The incidence matrix encodes the structure of the network at the level of complexes, i.e. how complexes are connected by reactions. If we think of complexes as network nodes, the incidence matrix associates an edge to each reaction pathway and the resulting topological structure is a reaction graph, e.g. fig. 1 and eqs. (21) and (29). The stoichiometric matrix instead encodes the structure of the network at the level of species. If we think of species as the network nodes, the stoichiometric matrix does not define a graph since reaction connects more than a pair of species, in general. The structure originating is rather a hyper-graph [56,65], or equivalently a Petri net [66,67].
Example 8. The composition matrix and the incidence matrix of the CRN in (29) are where the complexes are ordered as in example 7. The corresponding reaction hyper-graph is where only the forward reactions are depicted.
In an open CRN, we regroup all complexes γ ∈ C of the closed CRN which have the same stoichiometry for the internal species (i.e. all complexes with the same internal part of the composition matrix Γ X γ regardless of the chemostatted part Γ Y γ ) in sets denoted by C j , for j = 1, 2, . . . . Complexes of the closed network made solely of chemostatted species in the open CRN are all regrouped in the same complex C 0 . This allows to decompose the internal species stoichiometric matrix as where {Γ σx j ≡ Γ σx γ , for γ ∈ C j } are the entries of the composition matrix corresponding to the internal species, and {∂ j ρ ≡ γ∈Cj ∂ γ ρ } are the entries of the incidence matrix describing the network of regrouped complexes. This regrouping corresponds to the-equivalent-CRN only made of internal species with effective rate constant {k ±ρ Z σy ∇ ±ρ σy } ruling each reaction.

Example 9.
Let us consider the CRN (29) where the species X a and X c are chemostatted. The five complexes of the closed network, see example 7, are regrouped as In terms of these groups of complexes, the composition matrix and incidence matrix are which corresponds to the effective representation A steady-state distribution {Z σx } (11) is said to be complex-balanced if the net current flowing in each group of complexes C j vanishes, i.e. if the currents {J ρ } satisfy Complex-balance steady states are therefore a subclass of steady states (11a) which include equilibrium ones (17) as a special case While for generic steady states only the internal species formation rates vanish, for complex-balanced ones the complex formation rates also vanish. For a fixed kinetics ({k ±ρ }) and chemostatting (S y and {Z σy }), a CRN is complex-balanced if its dynamics exhibits a complex-balanced steady state (37) [35,36]. The complex-balanced distribution (37) depends on the unbroken components {L λu }, which can be inferred from the initial conditions, and is always globally stable [68]. Hence, complex-balanced networks always relax to a-complexbalanced-steady state. Detailed-balanced networks are a subclass of complex-balanced networks.
Whether or not a CRN is complex-balanced depends on the network topology (∇), the kinetics ({k ±ρ }) and the chemostatting (S y and {Z σy }). For any given network topology and set of chemostatted species S y , one can always find a set of effective rate constants {k ±ρ Z σy ∇ ±ρ σy } which makes that CRN complex-balanced [37]. However, for some CRNs, this set coincides with the one which makes the CRN detailed-balanced [69]. A characterization of the set of effective rate constants which make a CRN complex balanced is reported in Refs. [37,69].
Deficiency-zero CRNs are a class of CRN which is complex-balanced irrespective of the effective kinetics {k ±ρ Z σy ∇ ±ρ σy } [35][36][37]. The network deficiency is a topological property of the CRN which we briefly discuss in app. D, see Refs. [22,52,53] for more details. Consequently, regardless of the way in which a deficiency-zero CRN is driven in time, it will always remain complexbalanced. Throughout the paper, we will refer to these CRNs as unconditionally complex-balanced, as in the seminal work [35].

Example 10.
The open CRN (36) has a single steady stateZ b for any given set of rate constants and chemostatted concentrations Z a and Z c [64], defined by (11a) If the stronger condition (37) holds which is equivalent to the steady state is complex-balanced. Yet, if the steadystate currents are all independently vanishing, i.e., eq. (41) is equal to zero, then the steady state is detailed-balanced. When, for simplicity, all rate constants are taken as 1, the complex-balanced set of quadratic equations (41) The former case corresponds to a genuine complex-balanced state,Z b = 1 with currents J 1 =J 2 =J 3 = 1 − Z c , while the second to a detailedbalance stateZ b = √ Z c with vanishing currents. When, for example, Z a = 1 and Z c = 4, neither of the two previous conditions holds: the nonequilibrium steady state

Example 11. Let us now consider the following open
where the species Y a and Y e are chemostatted. Out of the four complexes of the closed network, {Y a , X b , X c + X d , Y e }, two are grouped into C 0 = {Y a , Y e } and the other two remain This network is deficiency-zero and hence unconditionally complex-balanced [22]. Therefore, given any set of rate constants k ±1 , k ±2 , and k ±3 , and the chemostatted concentrations Z a and Z e , the steady state of this CRN is complex-balanced, i.e. the steady state always satisfies a set of condition like those in eq. (40). Indeed, contrary to example 10, steady state currents {J 1 ,J 2 ,J 3 } different from each other cannot exist since they would induce a growth or decrease of some concentrations.

III. THERMODYNAMICS OF CHEMICAL NETWORKS
Using local equilibrium, we here build the connection between the dynamics and the nonequilibrium thermodynamics for arbitrary CRNs. In the spirit of Stochastic Thermodynamics, we derive an energy and entropy balance, and express the dissipation of the CRN as the difference between the chemical work done by the reservoirs on the CRN and its change in nonequilibrium free energy. We finally discuss the information-theoretical content of the nonequilibrium free energy and its relation to the dynamical potentials used in Chemical Reaction Network Theory.

A. Local Equilibrium
Since we consider homogeneous reaction mixtures in ideal dilute solutions, the assumption of local equilibrium [57, § 15.1] [70] means that the equilibration following any reaction event is much faster than any reaction time scale. Thus, what is assumed is that the nonequilibrium nature of the thermodynamic description is solely due to the reaction mechanisms. If all reactions could be instantaneously shut down, the state of the whole CRN would immediately become an equilibrated ideal mixture of species. As a result, all the intensive thermodynamic variables are well-defined and equal everywhere in the system. The temperature T is set by the solvent, which acts as a thermal bath, while the pressure p is set by the environment the solution is exposed to. As a result, each chemical species is characterized by a chemical potential [23, § 3.1] where R denotes the gas constant and {µ • σ ≡ µ • σ (T )} are the standard-state chemical potentials, which depend on the temperature and on the nature of the solvent. The total concentration of the solution is denoted by Z tot = σ Z σ + Z 0 , where Z 0 is the concentration of the solvent. We assume for simplicity that the solvent does not react with the solutes. In case it does, our results still hold provided one treats the solvent as a nondriven chemostatted species, as discussed in app. A. Since the solvent is much more abundant than the solutes, the total concentration is almost equal to that of the solvent which is a constant, Z tot Z 0 . Without loss of generality, the constant term −RT ln Z tot −RT ln Z 0 in eq. (45) is absorbed in the standard-state chemical potentials. Consequently, many equations appear with non-matching dimensions. We also emphasize that standard state quantities, denoted with " • ", are defined as those measured in ideal conditions, at standard pressure (p • = 100 kPa) and molar concentration (Z • σ = 1 mol/dm 3 ), but not at a standard temperature [71, p. 61].
Due to the assumption of local equilibrium and homogeneous reaction mixture, the densities of all extensive thermodynamic quantities are well-defined and equal everywhere in the system. With a slight abuse of notation, we use the same symbol and name for densities as for their corresponding extensive quantity. E.g. S is the molar entropy divided by the volume of the solution, but we denote it as entropy. We apply the same logic to rates of change. E.g. we call entropy production rate the molar entropy production density rate.

B. Affinities, Emergent Affinities and Local Detailed Balance
The thermodynamic forces driving reactions are given by differences of chemical potential (45) The local detailed-balance (18) allows us to express these thermodynamic forces in terms of reaction affinities, gives the external thermodynamic forces the network is coupled to, as we shall see in eq. (61), and thus provide a thermodynamic meaning to the cycle affinities (28).
Combining the detailed-balance property (18) and the equilibrium condition on the affinities A eq ρ = 0 (46), we can relate the Gibbs free energies of reaction to the rate constants This relation is the thermodynamic counterpart of the local detailed balance (18). It plays the same role as in Stochastic Thermodynamics, namely connecting the thermodynamic description to the stochastic dynamics. We emphasize that the local detailed balance property as well as the local equilibrium assumption by no mean imply that the CRN operates close to equilibrium. Their importance is to assign well defined equilibrium potentials to the states of the CRN, which are then connected by the nonequilibrium mechanisms, i.e. reactions.

C. Enthalpies and Entropies of Reaction
To identify the heat produced by the CRN, we need to distinguish the enthalpic change produced by each reaction from the entropic one. We consider the decomposition of the standard state chemical potentials [23, § 3.2]: The standard enthalpies of formation {h • σ } take into account the enthalpic contributions carried by each species [ The entropies of formation {s σ ≡ s • σ − R ln Z σ } account for the entropic contribution of each species in the CRN [23, § 3.2]. Entropy changes along reactions are given by called entropies of reaction [23, § 3.2].

Entropy Production Rate
The entropy production rate is a non-negative measure of the break of detailed balance in each chemical reaction. Its typical form is given by [57, § 9.5
Furthermore it can be rewritten in a thermodynamically appealing way using (48): It can be further expressed as the sum of two distinct contributions [56]: The first term is due to changes in the internal species and thus vanishes at steady state. The second term is due to the chemostats. It takes into account both the exchange of chemostatted species and the time-dependent driving of their concentration. If the system reaches a nonequilibrium steady state, the external currents {Ī σy } do not vanish and the entropy production reads which clearly emphasizes the crucial role of emergent cycles in steady state dissipation.

Entropy Flow Rate
The entropy flow rate measures the reversible entropy changes in the environment due to exchange processes with the system [57]. Using the expressions for the enthalpy of reaction (52) and entropy of formation (53), we express the entropy flow rate as The first contribution is the heat flow rate (positive if heat is absorbed by the system). When divided by temperature, it measures minus the entropy changes in the thermal bath. The second contribution accounts for minus the entropy change in the chemostats.

System Entropy
The entropy of the ideal dilute solution constituting the CRN is given by (see app. A) The total concentration term and the constant S 0 together represent the entropic contribution of the solvent. S 0 may also account for the entropy of chemical species not involved in the reactions. We also prove in app. B that the entropy (63) can be obtained as a large particle limit of the stochastic entropy of CRNs. S would be an equilibrium entropy if the reactions could be all shut down. But in presence of reactions, it becomes the nonequilibrium entropy of the CRN. Indeed, using eqs. (53), (58) and (62), we find that its change can be expressed as This relation is the nonequilibrium formulation of the Second Law of Thermodynamics for CRNs. It demonstrates that the non-negative entropy production (55) measures the entropy changes in the system plus those in the reservoirs (thermal and chemostats) [57].

First Law of Thermodynamics
Since the CRN is kept at constant pressure p, its enthalpy is equal to the CRN internal energy, up to a constant. Indeed, the enthalpy H is a density which, when written in terms of the internal energy (density) U , reads H = U + p.
Using the rate equations (9) and (10), the enthalpy rate of change can be expressed as the sum of the heat flow rate, defined in eq. (62), and the enthalpy of formation exchange rate The last term on the r.h.s. of eq. (68) is the free energy exchanged with the chemostats. It represents the chemical work rate performed by the chemostats on the CRN [21,23]Ẇ c ≡ I σy µ σy .
Either eq. (67) or (68) may be considered as the nonequilibrium formulation of the First Law of Thermodynamics for CRNs. The former has the advantage to solely focus on energy exchanges. The latter contains entropic contributions but is appealing because it involves the chemical work (69).

Nonequilibrium Gibbs Free Energy
We are now in the position to introduce the thermodynamic potential regulating CRNs. The Gibbs free energy of ideal dilute solutions reads As for entropy, the total concentration term −RT Z S and the constant G 0 represent the contribution of the solvent (see app. A). Furthermore, in presence of reactions, G becomes the nonequilibrium Gibbs free energy of CRNs.
We will now show that the nonequilibrium Gibbs free energy of a closed CRN is always greater or equal than its corresponding equilibrium form. A generic nonequilibrium concentration distribution {Z σ } is characterized by the set of components {L λ = λ σ Z σ }. Let {Z σ eq } be the corresponding equilibrium distribution defined by the detailed-balance property (18) and characterized by the same set of components {L λ } (a formal expression for the equilibrium distribution will be given in eq. (88)). At equilibrium, the Gibbs free energy (70) reads As discussed in § III B, the equilibrium chemical potentials must satisfy ∇ σ ρ µ eq σ = 0. We deduce that µ eq σ must be a linear combination of the closed system conservation laws (12) where {f λ } are real coefficients. Thus, we can write the equilibrium Gibbs free energy as In this form, the first term of the Gibbs free energy appears as a bilinear form of components {L λ } and conjugated generalized forces {f λ } [23, § 3.3], which can be thought of as chemical potentials of the components. From eq. (72) and the properties of components (13), the equality Z σ eq µ eq σ = Z σ µ eq σ follows. Hence, using the definition of chemical potential (45), the nonequilibrium Gibbs free energy G of the generic distribution {Z σ } defined above is related to G eq (73) by where we introduced the relative entropy for nonnormalized concentration distributions, also called Shear Lyapunov Function, or pseudo-Helmholtz function [35,73,74] This quantity is a natural generalization of the relative entropy, or Kullback-Leibler divergence, used to compare two normalized probability distributions [75]. For simplicity, we still refer to it as relative entropy. It quantifies the distance between two distributions: it is always positive and vanishes only if the two distributions are identical, {Z σ } = {Z σ }. Hence, eq. (74) proves that the nonequilibrium Gibbs free energy of a closed CRN is always greater or equal than its corresponding equilibrium form, G ≥ G eq . We now proceed to show that the nonequilibrium Gibbs free energy is minimized by the dynamics in closed CRNs, viz. G-or equivalently L {Z σ }|{Z σ eq } [59, 76]-acts as a Lyapunov function in closed CRNs. Indeed, the time derivative of G (70) always reads When using the rate equation for closed CRNs (7) we find that d t G = −J ρ ∇ σ ρ µ σ . Using eq. (74) together with eqs. (46) and (58), we get which proves the aforementioned result.

Chemical Work
In arbitrary CRNs, the rate of change of nonequilibrium Gibbs free energy (76) can be related to the entropy production rate (59) using the rate equations of open CRN (9) and (10) and the chemical work rate (69), This important results shows that the positivity of the entropy production sets an intrinsic limit on the chemical work that the chemostats must perform on the CRN to change its concentration distribution. The equality sign is achieved for quasi-static transformations (Ṡ i 0). If we now integrate eq. (78) along a transformation generated by an arbitrary time dependent protocol π(t), which drives the CRN from an initial concentration distribution {Z σ i } to a final one {Z σ f }, we find where ∆G = G f − G i is the difference of nonequilibrium Gibbs free energies between the final and the initial state. Let us also consider the equilibrium state {Z σ eq i } (resp. {Z σ eq f }) obtained from {Z σ i } (resp. {Z σ f }) if one closes the network (i.e. interrupt the chemostatting procedure) and let it relax to equilibrium, as illustrated in fig. 4. The Gibbs free energy difference between these two equilibrium distributions, ∆G eq = G eq f − G eq i , is related to ∆G via the difference of relative entropies, eq. (74), where Thus, the chemical work (79) can be rewritten as which is a key result of our paper. ∆G eq represents the reversible work needed to reversibly transform the CRN from {Z σ eq i } to {Z σ eq f }. Implementing such a reversible transformation may be difficult to achieve in practice. However, it allows us to interpret the difference W irr c ≡ W c − ∆G eq in eq. (82) as the chemical work dissipated during the nonequilibrium transformation, i.e. the irreversible chemical work. The positivity of the entropy production implies that This relation sets limits on the irreversible chemical work involved in arbitrary far-from-equilibrium transformations. For transformations connecting two equilibrium distributions, we get the expected inequality W irr c ≥ 0. More interestingly, eq. (83) tells us how much chemical work the chemostat need to provide to create a nonequilibrium distribution from an equilibrium one. It also tells us how much chemical work can be extracted from a CRN relaxing to equilibrium.
The conceptual analogue of (82) in Stochastic Thermodynamics (where probability distributions replace nonnormalized concentration distributions) is called the nonequilibrium Landauer's principle [42,43] (see also [77][78][79]). It has been shown to play a crucial role to analyze the thermodynamic cost of information processing (e.g. for Maxwell demons, feedback control or proofreading). The inequality (83) is therefore a nonequilibrium Landauer's principle for CRN.

IV. THERMODYNAMICS OF COMPLEX-BALANCED NETWORKS
In this section, we focus on unconditionally complexbalanced networks. We shall see that the thermodynamics of these networks bears remarkable similarities with Stochastic Thermodynamics. Let us first observe that whenever a CRN displays a well-defined steady-state distribution {Z σx }, the entropy production rate (55) can be formally decomposed as the sum of an adiabatic and a nonadiabatic contribution in analogy to what was done in Stochastic Thermodynamics [46][47][48][49][50]. As discussed in § II E, unconditionally complex-balanced networks have a unique steady-state distribution Z σx ≡Z σx (π(t)) , eq. (37), for any value of the chemostatted concentrations {Z σy ≡ Z σy (π(t))} and of the fixed unbroken components {L λu }. The decomposition (84) is thus well-defined at any time, for any protocol π(t). As a central result, we prove in app. C that the adiabatic and nonadiabatic contribution are non-negative for unconditionally complex-balanced networks as well as for complex-balanced networks without time-dependent driving.
The adiabatic entropy production rate encodes the dissipation of the steady state Z σx . It can be rewritten in terms of the steady state Gibbs free energy of reaction {∆ rḠρ } (48) as This inequality highlights the fact that the transient dynamics-generating the currents {J ρ }-is constrained by the thermodynamics of the complex balanced steady state, i.e. by {∆ rḠρ }. The nonadiabatic entropy production rate characterizes the dissipation of the transient dynamics. It can be decomposed as whereZ Sx = σx∈SxZ σx (see Refs. [46,48] for the analogous decomposition in the stochastic context). The first term is proportional to the time derivative of the relative entropy (75) between the nonequilibrium concentration distribution at time t and the corresponding complexbalanced steady-state distribution. Hence, it describes the dissipation of the relaxation towards the steady state. The second term, TṠ d , is related to the time-dependent driving performed via the chemostatted species and thus denoted driving entropy production rate [46]. It vanishes in nondriven networks where we obtaiṅ This result shows the role of the relative entropy L {Z σx }|{Z σx } as a Lyapunov function in nondriven complex-balanced networks with mass action kinetics. It was known in the mathematical literature [35,80], but we provide a clear thermodynamic interpretation to this result by demonstrating that it derives from the nonadiabatic entropy production rate. We mention that an alternative derivation of the adiabatic-nonadiabatic decomposition for nondriven complex-balanced networks with mass action kinetics was found in Ref. [51], while we were finalizing our paper.

V. THERMODYNAMICS OF OPEN DETAILED-BALANCED NETWORKS
We finish our study by considering detailed-balanced networks. We discuss the equilibrium distribution, introduce a new class of nonequilibrium potentials and derive a new work inequality.
Let us also emphasize that open detailed-balanced CRNs are a special class of open complex-balanced CRNs for which the adiabatic entropy production rate vanishes (since the steady state is detailed balanced) and thus the nonadiabatic entropy production characterizes the entire dissipation.

A. Equilibrium Distribution
As discussed in § II D, for given kinetics {k ±ρ }, chemostatting {Z σy } and unbroken components {L λu }, detailed-balanced networks always relax to a unique equilibrium distribution. Since the equilibrium chemical potentials can be expressed as a linear combination of conservation laws, eq. (72), we can express the equilibrium distribution as inverting the expression for the chemical potentials (45).
Since the independent set of unbroken conservation laws, { λu }, are such that λu σy = 0, ∀λ u , σ y , see § II C, we have that We thus conclude that the |λ b | broken generalized forces {f λ b } only depend on the chemostatted concentrations {Z σy }. Instead, the remaining |λ u | unbroken generalized forces f λu can be determined by inverting the nonlinear set of equations L λu = λu σx Z σx eq . They therefore depend on both {Z σy } and {L λu }.

B. Open nondriven networks
As a consequence of the break of conservation laws, the nonequilibrium Gibbs free energy G (70) is no longer minimized at equilibrium in open detailed-balanced networks. In analogy to equilibrium thermodynamics [23], the proper thermodynamic potential is obtained from G by subtracting the energetic contribution of the broken conservation laws. This transformed nonequilibrium Gibbs free energy reads We proceed to show that G is minimized by the dynamics in nondriven open detailed-balanced networks. Let {Z σx } be a generic concentration distribution in a detailed-balanced network characterized by {L λu } and {Z σy }, and let {Z σx eq } be its corresponding equilibrium. Using the relation between equilibrium chemical potentials and conservation laws (72), the transformed Gibbs free energy (90) at equilibrium reads Yet, combinig eq. (72) and the properties of unbroken components, one can readily show that Z σ eq µ eq The relation between the nonequilibrium G and the corresponding equilibrium value thus follows (we show in app. A the derivation of the latter in presence of reacting solvent). The non-negativity of the relative entropy for concentration distributions L {Z σx }|{Z σx eq } ensures that the nonequilibrium transformed Gibbs free energy is always greater or equal to its equilibrium value, G ≥ G eq . Since entropy production and nonadiabatic entropy production coincide, using eqs. (87) and (92), we obtain which demonstrates the role of G as a Lyapunov function. The relative entropy L ({Z σx }|{Z σx }) was known to be a Lyapunov function for detailed-balanced networks [76,81], but we provided its clear connection to the transformed nonequilibrium Gibbs free energy. To summarize, instead of minimizing the nonequilibrium Gibbs free energy G (70) as in closed CRNs, the dynamics minimizes the transformed nonequilibrium Gibbs free energy G in open nondriven detailed-balanced CRNs.

C. Open driven networks
We now consider unconditionally detailed-balanced CRNs. As discussed in § II D, they are characterized by a unique equilibrium distribution Z σx eq ≡ Z σx eq (π(t)) , defined by eq. (18), for any value of the chemostatted concentrations {Z σy = Z σy (π(t))}.
We start by showing that the external fluxes {I σy } can be expressed as influx rate of moieties. Since the CRN is open and unconditionally detailed balanced, each chemostatted species broke a conservation law (no emergent cycle is created, § II D). Therefore, the matrix whose entries are { λ b σy } in eq. (89) is square and also nonsingular [82]. We can thus invert eq. (89) to get where {ˆ σy λ b } denote the entries of the inverse matrix of that with entries { λ b σy }. Hence, using the definition of broken component, From the rate equations for the chemostatted concentrations (10), we find that We can thus interpret M σy as the concentration of a moiety which is exchanged with the environment only through the chemostatted species X σy . Eq. (95) shows that the energetic contribution of the broken components can be expressed as the Gibbs free energy carried by these specific moieties.  fig. 2, whose conservation laws given in example 4, the concentrations of the exchanged moieties are For the specific implementation of that CRN, fig. 3, the first term (resp. second term) is the total number of moiety 2H (resp. C) in the system, which can be exchanged with the environment only via the chemostatted species H 2 O (resp. CO).
We now turn to the new work relation. From the general work relation (78), using (90) and (95), we find where the driving work due to the time-dependent driving of the chemostatted species is obtained using the chemical work rate (69) together with eqs. (95) and (96) Equivalently, the driving work rate (99) can be defined as the rate of change of the transformed Gibbs free energy (90) due to the time dependent driving only, i.e.
To relate this alternative definition to eq. (99), all {Z σy } must be expressed in terms of {µ eq σy } using the definition of chemical potential (45).
The driving work rateẆ d vanishes in nondriven CRNs, where (98) reduces to (93). After demonstrating that the entropy production rate is always proportional the difference between the chemical work rate and the change of nonequilibrium Gibbs free energy in eq. (79), we showed that, for unconditionally detailed-balanced CRNs, it is also proportional the difference between the driving work rate and the change in transformed nonequilibrium Gibbs free energy, eq. (98).
We end by formulating a nonequilibrium Landauer's principle for the driving work instead the chemical work done in § III E 3. We consider a time-dependent transformation driving the unconditionally detailed-balanced {Z σ eq f }) denotes the equilibrium distribution obtained from {Z σ i } (resp. {Z σ f }) by stopping the time-dependent driving and letting the system relax towards the equilibrium, fig. 4. Note that this reference equilibrium state is different from the one obtained by closing the network in § III E 3. Integrating (98) over time and using (92), we get where ∆G eq represents the reversible driving work, and the irreversible driving work satisfies the inequality This central relation sets limits on the irreversible work spent to manipulate nonequilibrium distributions. It is a nonequilibrium Landauer's principle for the driving work by the same reasons why inequality (83) is a nonequilibrium Landauer's principle for the chemical work. The key difference is that the choice of the reference equilibrium state is different in the two cases. The above discussed inequality (103) only holds for unconditionally detailed-balanced CRNs while eq. (83) is valid for any CRNs.

VI. CONCLUSIONS AND PERSPECTIVES
Following a strategy reminiscent of Stochastic Thermodynamics, we systematically build a nonequilibrium thermodynamic description for open driven CRNs made of elementary reactions in homogeneous ideal dilute solutions. The dynamics is described by deterministic rate equations whose kinetics satisfies mass action law. Our framework in not restricted to steady states and allows to treat transients as well as time-dependent drivings performed by externally controlled chemostatted concentrations. Our theory embeds the nonequilibrium thermodynamic framework of irreversible processes established by the Brussels School of Thermodynamics.
We now summarize our results. Starting from the expression for the entropy production rate, we established a nonequilibrium formulation of the first and second law of thermodynamics for CRNs. The resulting expression for the system entropy is that of an ideal dilute solution. The clear separation between chemostatted and internal species allowed us to identify the chemical work done by the chemostats on the CRN and to relate it to the nonequilibrium Gibbs potential. We were also able to express the minimal chemical work necessary to change the nonequilibrium distribution of species in the CRN as a difference of relative entropies for non-normalized distributions. These latter measure the distance of the initial and final concentration distributions from their corresponding equilibrium ones, obtained by closing the network. This result is reminiscent of the nonequilibrium Landauer's principle derived in Stochastic Thermodynamics [43] and which proved very useful to study the energetic cost of information processing [45]. We also highlighted the deep relationship between the topology of CRNs, their dynamics and their thermodynamics. Closed CRNs (resp. nondriven open detailed-balanced networks) always relax to a unique equilibrium by minimizing their nonequilibrium Gibbs free energy (resp. transformed nonequilibrium Gibbs free energy). This latter is given, up to a constant, by the relative entropy between the nonequilibrium and equilibrium concentration distribution. Non-driven complex-balanced networks relax to complex-balanced nonequilibrium steady states by minimizing the relative entropy between the nonequilibrium and steady state concentration distribution. In all these cases, even in presence of driving, we showed how the rate of change of the relative entropy relates to the CRN dissipation. For complex-balanced networks, we also demonstrated that the entropy production rate can be decomposed, as in Stochastic Thermodynamics, in its adiabatic and nonadiabatic contributions quantifying respectively the dissipation of the steady state and of the transient dynamics.
Our framework could be used to shed new light on a broad range of problems. We mention only a few.
Our theory could also be used to study metabolic networks. However, these require some care, since complex enzymatic reaction mechanisms are involved [106]. Nevertheless our framework provides a basis to build effective coarse-graining procedures for enzymatic reactions [107]. For instance, proofreading mechanisms operating in metabolic processes could be considered [108]. We foresee an increasing use of thermodynamics to improve the modeling of metabolic networks, as recently shown in Refs. [30,32,33].
Since our framework accounts for time-dependent drivings and transient dynamics, it could be used to represent the transmission of signals through CRNs or their response to external modulations in the environment. These features become crucial when considering problems such as signal transduction and biochemical switches [24,109,110], biochemical oscillations [28,111], growth and self-organization in evolving bio-systems [112,113], or sensory mechanisms [85,87,[114][115][116][117]. Also, since transient signals in CRNs can be used for computation [118] and have been shown to display Turing-universality [119][120][121][122], one should be able to study the thermodynamic cost of chemical computing [123].
As closing words, we believe that our results constitute an important contribution to the theoretical study of CRNs. It does for nonlinear chemical kinetics what Stochastic Thermodynamics has done for stochastic dynamics, namely build a systematic nonequilibrium thermodynamics on top of the dynamics. It also opens many new perspectives and builds bridges between approaches from different communities studying CRNs: mathematicians who study CRNs as dynamical systems, physicists who study them as nonequilibrium complex systems, and biochemists as well as bioengineers who aim for accurate models of metabolic networks.
We show that the nonequilibrium Gibbs free energy (70) is the Gibbs free energy of an ideal dilute solution [131,Ch. 7] (see also [51]). We also show that in open detailed-balanced networks in which the solvent reacts with the solutes, the expression of the transformed Gibbs free energy (92) is recovered by treating the solvent as a special chemostatted species.
The Gibbs free energy (density) of an ideal dilute mixture of chemical compounds kept at constant temperature and pressure reads where the labels σ ∈ S refer to the solutes and 0 to the solvent. The chemical potentials of each species (45) read Since the solution is dilute, Z tot = σ∈S Z σ + Z 0 Z 0 and the standard state chemical potentials {µ • σ } depend on the nature of the solvent. Hence, the chemical potentials of the solutes read while that of the solvent where Z S ≡ σ∈S Z σ . Therefore, the Gibbs free energy (A1) reads which is eq. (70) in the main text, where G 0 is equal to Z 0 µ • 0 plus possibly the Gibbs free energy of solutes which do not react.
We now consider the case where the solvent reacts with the solutes. We assume that both the solutes and the solvent react according to the stoichiometric matrix where the first row refers to the solvent, the second block of rows to the internal species and the last one to the chemostatted species. The solvent is treated as a chemostatted species such that d t Z 0 = 0. In order to recover the expression for the transformed Gibbs free energy (92) in unconditionally detailed balanced networks, we observe that, at equilibrium ∇ σ ρ µ eq σ + ∇ 0 ρ µ eq 0 = 0 .
Therefore, the equilibrium chemical potentials are a linear combination of the conservation laws of ∇ (A6) µ eq σ = f λ λ σ µ eq 0 = f λ λ 0 . (A8) As mentioned in the main text, § II C, the chemostatting procedure breaks some conservation laws, which are labeled by λ b . The unbroken ones are labeled by λ u . The transformed Gibbs free energy is defined as in eq. (90), reported here for convenience where G reads as in eq. (A1), {L λ b } are the broken components and {f λ b } are here interpreted as the conjugated generalized forces. Adding and subtracting the term Z σ µ eq σ + Z 0 µ eq 0 from the last equation and using eq. (A8) we obtain G = G eq + Z σ (µ σ − µ eq σ ) + Z 0 (µ 0 − µ eq 0 ) , (A10) where G eq = f λu L λu .
From eqs. (A3) and (A4) and the fact that Z σy = Z σy eq and Z 0 = Z eq 0 we obtain in agreement with the expression derived in the main text, eq. (92).

Appendix B: Entropy of CRNs
We show how the nonequilibrium entropy (63) can be obtained as a large particle limit of the stochastic entropy. We point out that while finalizing the paper similar derivations for other thermodynamic quantities have obtained in Refs. [51,132].
In the stochastic description of CRNs, the state is characterized by the population vector n = {n σ }. The probability to find the network is in state n at time t is denoted p t (n). The stochastic entropy of that state reads [21,107], up to constants, S(n) = −k B ln p t (n) + s(n) . (B1) The first term is a Shannon-like contribution while the second term is the configurational entropy s • σ is the standard entropy of one single X σ molecule, and n 0 is the very large number of solvent molecules. We now assume that the probability becomes very narrow in the large particle limit n σ 1 and behaves as a discrete delta function p t (n) δ(n −n(t)). The vector n(t) ≡ {n σ } denotes the most probable and macroscopic amount of chemical species, such that Z σ =n σ /(V N A ). Hence, the average entropy becomes  Dividing by V and using the relation R = N A k B we finally get the macroscopic entropy density (63) where the (molar) standard entropies of formation s • σ reads Mindful of the information-theoretical interpretation of the entropy [133], we note that the uncertainty due to the stochasticity of the state disappears (the first term on the r.h.s. of eq. (B1)). However, the uncertainty due to the indistinguishability of the molecules of the same species-quantified by the configurational entropy (B2)-remains and contributes to the whole deterministic entropy function (63).
where ψ eq γ ≡ Z eq σ Γ σ γ . In order to prove the non-negativity of the adiabatic term (84), we rewrite it aṡ The detailed balance property is used in the first equality, and the decomposition of the stoichiometric matrix (30) in the second one. Also, the constant RT is taken equal to one. Using (C3), eq. (C5) can be rewritten aṡ From the log inequality − ln x ≥ 1 − x and the detailed balance property (C4), we obtaiṅ The last equality follows from the assumption of complexbalanced steady state (C2), the properties of the groups of complexes {C j } ( § II E), and the fact that {Z σy =Z σy }. Indeed Concerning the nonadiabatic term (84), using the rate equations (9) and the fact that {Z σy =Z σy }, we can rewrite it aṡ Because of (C3), we further get thaṫ From the log inequality − ln x ≥ 1 − x and from (C4) The last equality again follows from the assumption of complex-balance steady state (C2) as in eq. (C8).