Thermodynamics of Non-Elementary Chemical Reaction Networks

We develop a thermodynamic framework for closed and open chemical networks applicable to non-elementary reactions that do not need to obey mass action kinetics. It only requires the knowledge of the kinetics and of the standard chemical potentials, and makes use of the topological properties of the network (conservation laws and cycles). Our approach is proven to be exact if the network results from a bigger network of elementary reactions where the fast-evolving species have been coarse grained. Our work should be particularly relevant for energetic considerations in biosystems where the characterization of the elementary dynamics is seldomly achieved.


I. INTRODUCTION
Many processes in biology result from the combined effect of numerous elementary chemical reactions obeying mass-action kinetics. Well-known examples are enzymatic reactions, signaling cascades, and gene regulations. Since their detailed kinetic characterization at the elementary level is hard to achieve, they are o en modeled with e ective schemes obeying non-mass action kinetics, e.g., the Michaelis and Menten mechanism [1,2], Hill functions [3] and many others [4]. ese models can be thought of as resulting from a coarse graining of elementary reactions by eliminating fastevolving species and their use is well established [5][6][7][8][9][10][11].
However, characterizing the energetics at the level of these e ective schemes remains an open challenge especially outof-equilibrium. Inspired by developments in stochastic thermodynamics [12][13][14], the nonequilibrium thermodynamics of open chemical networks of elementary reactions is nowadays well established [15][16][17]. Recently, a thermodynamically consistent coarse graining strategy of the dynamics of biocatalysts was proposed in Ref. [18]. At steady-state this approach reproduces the correct dissipation at the coarse grained level.
Building on this approach, in the present work we develop a thermodynamics directly applicable to (closed and open) chemical networks of non-elementary reactions described by deterministic non-mass action rate equations. Crucially, this approach can be validated using thermodynamics of elementary reactions. Indeed, if the non-elementary network can be constructed from a network of elementary reactions by coarse graining the fast evolving species, one can show that the thermodynamic quantities evaluated in both are equivalent.
Our work is structured as follows. Starting from a network of elementary reactions and assuming time scale separation, we derive a thermodynamics for the network obtained after coarse graining the fast evolving species. We start by considering the dynamics of the closed network of the elementary reactions in Sec. II and of the coarse grained ones in * francesco.avanzini@uni.lu † gianmaria.falasco@uni.lu ‡ massimiliano.esposito@uni.lu Sec. III. We then proceed by building the corresponding thermodynamics in Sec. IV. e dynamics of the open networks is considered in Sec. V and its thermodynamics in Sec. VI. In Sec. VII we discuss how our previous ndings can be used to characterize the thermodynamics of non-elementary networks when the full elementary description is not known. In Sec. VIII we summarize our results and discuss their implications.

II. ELEMENTARY DYNAMICS OF CLOSED CRNS
We consider here chemical reaction networks (CRNs) composed of chemical species α = (. . . , α, . . . ) undergoing elementary reactions [19], with ν ±ρ (resp. k ±ρ ) the vector of the stoichiometric coecients (resp. the kinetic constant) of the forward/backward reaction ρ. e state of deterministic CRNs is speci ed by the concentration vector z(t) = (. . . , [α](t), . . . ) of all the chemical species α . For closed CRNs, its dynamics follows the rate equation where we introduce the stoichiometric matrix S and the current vector j(z). e stoichiometric matrix S codi es the topology of the network. Each ρ column S ρ of the stoichiometric matrix speci es the net variation of the number of molecules for each species undergoing the ρ elementary reaction (1), S ρ = ν −ρ − ν +ρ . e current vector j(z) = (. . . , j ρ (z), . . . ) speci es the net reaction current for every ρ elementary reaction (1) as the di erence between the forward j +ρ (z) and backward reaction current j −ρ (z) with j ±ρ (z) satisfying the mass-action kinetics [20][21][22], Note that we use a b = i a b i i . In the following, we will refer to Eq. (2) as the elementary dynamics of closed CRNs. arXiv:2008.10495v1 [cond-mat.stat-mech] 24 Aug 2020 Example. In all the manuscript, we illustrate our ndings using a modi ed version of the model discussed in Ref. [10]. It represents the transformation of two identical substrates S into one product P. e process is catalyzed by a membrane enzyme E which interacts with the substrate a er it is adsorbed S m by the membrane: For this model, we introduced here the concentration vector and the current vector, as well as the stoichiometric matrix

A. Topological Properties
As shown in Ref. [15,23], the linear independent vectors { λ } of the cokernel of the stoichiometric matrix are the conservation laws of Eq. (2). Indeed, for each vector λ , the scalar L λ (z(t)) ≡ λ · z(t) is a conserved quantity, i.e., d t L λ (z(t)) = λ · Sj(z(t)) = 0. e linear independent vectors {c ι } of the kernel of the stoichiometric matrix are the internal cycles of Eq. (2). ey are sequence of reactions that leave the state of CRNs invariant. Any linear combination of internal cycles gives a steady-state current vector, i.e., Sj = 0 with j ≡ c ι ψ ι . Note that in all the work we use the Einstein notation: repeated upper-lower indices implies the summation over all the allowed values for the indices, i.e., Example. For the CRN (5), there are two conservation laws, with a clear physical interpretation. e total concentration of the enzyme is given by , while the total concentration of the substrate is given by (5) has no internal cycles.

B. Equilibrium
Closed CRNs must be detailed balanced, namely the rate equation (2) admits an equilibrium steady-state z eq characterized by vanishing reaction currents, is, together with mass-action kinetics, implies the so-called local detailed balance condition for the kinetic constants k ±ρ of the chemical reactions It ensures that i) z eq is the only steady-state of the rate equation (2) ii) z(t) relaxes to z eq and hence the thermodynamic consistency [23].

III. COARSE GRAINED DYNAMICS OF CLOSED CRNS
When the chemical species evolve over two di erent time scales, they can be divided into two disjoint subsets: the fastevolving species Q and the slow-evolving species P [7,8]. We apply the same spli ing to the stoichiometric matrix and to the concentration vector z = (q, p). is also allows us to split the rate equation (2) into d t p(t) = S P j(q(t), p(t)) .
e coarse graining procedure provides a closed dynamical equation for the slow species only. It is based on two assumptions collectively called time scale separation hypothesis or quasi-steady-state assumption.
e rst is the existence, for every concentration p of the slow species, of the concentration vectorq(p) such that the current is a steady-state current of Eq. (14) or, equivalently,ĵ(p) ∈ ker S Q . is is a topological property of the CRN. e second is the equivalence, during the evolution of z(t) = (q(t), p(t)) according to the rates equations (14) and (15), between the actual concentration vector of the fast species q(t) and the steady-state one, i.e., q(t) =q(p(t)). is means physically that the concentrations of the fast species relax instantaneously to the steady-state corresponding to the frozen values of p(t). is occurs when i) the reactions involving only the P species are much slower than the reactions involving only Q species and ii) the abundances of the P species are much higher than the abundances of the Q species. is latter condition ensures that when the P and Q species are coupled by a reaction, on the same time scale the concentrations of the Q species can dramatically change, the concentrations of the P species remain almost constant.
Hence, j(q(t), p(t)) in (15) becomes the steady-state currentĵ(p) of Eq. (14). It can now be wri en as a linear combination of the right null eigenvectors S Q c γ = 0, with p dependent coe cients {ψ γ }. e speci c expression of ψ γ (p) is not discussed here and is not fundamental in what follows. When the dynamics of the fast species at xed concentrations of the slow ones is linear, the problem is solved in Refs. [7,18] using a diagrammatic method developed in Ref. [24,25]. Note that {c γ } includes all the internal cycles {c ι } of Eq. (9) and, in general, other vectors {c ε } called pseudo-emergent cycles [15,18,23]. e former characterize a sequence of reactions which upon completion leaves all species (q, p) unchanged, while the la er only leaves the fast species unchanged. Each coe cient ψ γ (p) represents a current along the cycle γ . By employing the steady-state current vector (17) in Eq. (15), we obtain a closed dynamical equation for the slow species, where we used S P c ι = 0 because of Eq. (9). is coarse grained dynamics will be denoted as e ective for convenience. It can be rewri en in a more compact way, introducing the e ective stoichiometric matrixŜ and the effective current vector ψ(p). e e ective stoichiometric ma-trixŜ codi es the net stoichiometry of the slow species along the pseudo-emergent cycles. Each ε columnŜ ε = S P c ε speci es the net variation of the number of molecules for each slow species along the ε pseudo-emergent cycle, i.e., the stoichiometry of the e ective reactions. e e ective current vector ψ(p) collects the pseudo-emergent cycle currents ψ(p) = (. . . ,ψ ε (p), . . . ) . Unlike j(z) in Eq. (2), ψ(p) does not, in general, satisfy mass-action kinetics.
Note that, if there were no pseudo-emergent cycles, (q(p), p) would be an equilibrium condition of the elementary dynamics (2). e corresponding e ective dynamics for the slow species would become trivial: d t p(t) = 0.
Example. For the CRN (5), we split the chemical species as follows Underlying this separation is the assumption that the concentration of the enzyme species changes much more quickly.
Since there are no internal cycles of S of Eq. (7), all the right null eigenvectors of S Q , c ads = are pseudo-emergent cycles. us, the e ective stoichiometric matrix of Eq. (19) is given bŷ Each column ofŜ speci es the net stoichiometry of the following e ective reactions for the slow species: In Fig. 1, we compare the elementary dynamics of the CRN (5) for the slow species and the e ective one obtained by solving Eq. (2) and Eq. (19), respectively. For this case, the e ective current vector ψ(p) is computed according to the procedure introduced in Refs. [24,25]: with L E the total concentration of the enzyme and A. Topological Properties e conservation laws of the e ective dynamics (19) are de ned as the linear independent vectors of the cokernel of the e ective stoichiometric matrix [S](t ) el.
Note that all the conservation laws of the e ective stoichiometry matrixŜ can be wri en as in Eq. (28). is follows from the rank nullity theorem for S andŜ, and the absence of cycles forŜ. Indeed, suppose that there is a vector ϕ 0 such thatŜϕ = 0. is means that S P c ε ϕ ε = 0 and, consequently, c ε ϕ ε is a right null eigenvector of both S P and S Q , i.e., an internal cycle. Since c ε ϕ ε is a linear combination of only pseudo-emergent cycles {c ε }, the e ectively stoichiometric matrix cannot admit cycles.
Example. For the CRN (5) and given the spli ing between fast and slow species in (20), the only conservation law of S in (7) with no null entries for the slow species is S of Eq. (10). e corresponding conservation law of the e ective dynam- GivenŜ in Eq. (22), one easily veri es thatˆ S ·Ŝ = 0. e corresponding quantityL S =ˆ S · p has the clear physical interpretation of the total concentration of substrate: B. Equilibrium e equilibrium condition p eq of the e ective dynamics (19) is de ned by vanishing cycle currents, Its existence is granted by the local detailed balance condition (12). Indeed, the state (q(p eq ), p eq ) is an equilibrium of the elementary dynamics (2) by de nition: it is a steadystate of the rate equation (2) which admits only equilibrium steady-states. Let us discuss now the equilibrium state to which the elementary and the e ective dynamics relax given a unique initial condition z(0) = (q(0), p(0)).
e concentrations of the chemical species at equilibrium depend on the value of the conserved quantities. erefore, to have similar concentrations for the slow species at equilibrium, also the values of L ζ (z(0)) = l ζ · z(0) andL ζ (p(0)) =l ζ · p(0) have to be similar. is means that the concentration of the fast species involved in the ζ conservation laws must be much smaller than the concentration of the slow species. is is consistent with the time scale separation hypothesis requiring much higher abundances for the slow-evolving than for the fast-evolving species.

IV. THERMODYNAMICS OF CLOSED CRNS
We now derive the thermodynamics for the e ective dynamics starting from the full elementary network. We emphasize that all the resulting thermodynamic quantities (e.g., entropy production and free energy) can be evaluated at the level of the e ective dynamics, without any knowledge of the elementary dynamics. is is the key result of our work.
A. Local Detailed Balance e thermodynamic theory of CRNs presumes that all degrees of freedom other than concentrations are equilibrated at temperature T and pressure of the solvent. In this way, thermodynamic state functions can be speci ed by their equilibrium form but expressed in terms of nonequilibrium concentrations.
e vector of chemical potentials is thus given by with µ • the vector of the standard chemical potentials and R the gas constant. e local detailed balance condition (12) can be then restated to establish a correspondence between the elementary dynamics (i.e., kinetic constants k ±ρ ) and the thermodynamics (i.e., standard chemical potentials µ • ): We now formulate the local detailed balance condition for the e ective dynamics (19). is is done by taking the product of the ratio k +ρ /k −ρ along each pseudo-emergent cycle ε whereμ • = Pµ • is the vector of the standard chemical potentials for the slow species P. It is important to note that rst Eq. (33) establishes the same correspondence between the equilibrium concentrations p eq of the e ective dynamics (19) and the standard chemical potentialsμ • as Eq. (32) does between z eq and µ • for the elementary reactions. Second, the e ective stoichiometric matrixŜ plays the same role in Eq. (33) as S in Eq. (32).
ird, the derivation of Eq. (33) employs only the topological properties of the CRNs and does not require the time scale separation hypothesis.

B. Gibbs Free Energy of Reaction
e Gibbs free energy of the ρ elementary reaction (1), namely the thermodynamic force driving this reaction, is given by which becomes because of the local detailed balance (32). At equilibrium . is means that µ eq belongs to the cokernel of S. We de ne now the Gibbs free energy of each e ective reaction ∆ ε G(p). Collecting the Gibbs free energies of the elementary reactions in the vector ∆ r G ≡ (. . . , ∆ ρ G, . . . ) , ∆ ε G is given by We used Eq. (34), S Q c ε = 0 and introduced the vector of the chemical potential of the slow specieŝ e expression of ∆ ε G in Eq. (36) is formally the same as ∆ ρ G in Eq. (34). At equilibrium ∆ ε G eq ≡ ∆ ε G(p eq ) = 0, because of the local detailed balance condition (33). is means that µ eq belongs to cokernel ofŜ. Unlike ∆ ρ G(z), there is no analytical correspondence between ∆ ε G(p) and ψ ±ε (p) and, in general, Here, ψ ±ε (p) are two positive de ned currents such that ψ ε (p) = ψ +ε (p)−ψ −ε (p). is breaks the ux-force relation at the coarse grained level as was already stressed in Ref. [18]. Finally, we note that, as the local detailed balance, ∆ ε G is de ned using only topological properties of the CRNs.
Example. For the CRN (5) and given the spli ing between fast and slow species in (20), the Gibbs free energy of the e ective reactions in Eq. (23) is speci ed as C. Entropy Production Rate e entropy production rate of the elementary dynamics reads It quanti es the dissipation of the relaxation towards equilibrium. Because of Eq. (35), the dissipation of each elementary reaction ρ is also non-negative, −∆ ρ G(z(t))j ρ (z(t)) ≥ 0 (without summation over ρ).
We de ne the entropy production rate of the e ective dy-namicsˆ Σ using the expression in Eq. (40), but evaluated along the e ective trajectory (q(p(t)), p(t)). Using Eqs. (16), (17) and (36), we obtain that where we collected the Gibbs free energies of the e ective reactions in the vector ∆ c G ≡ (. . . , ∆ ε G, . . . ) . e e ective entropy production rateˆ Σ(t) is still non-negative by de nition: using Eq. (35) one proves that Σ(t) ≥ 0 for every j(z) includingĵ(p). However, unlike the elementary dynamics, it is not granted that −∆ ε G(p(t))ψ ε (p(t)) ≥ 0 for every e ective reaction ε due to the lack of a ux-force relation (38) at the coarse grained level.
If the time scale separation hypothesis holds, namely j(z) = j(p), the two entropy production rates coincide: When it does not hold, no bound constrains one entropy production rate respect to the other andˆ Σ can be lower or greater than Σ. is was also observed in Ref. [26] for driven systems. Example. For the CRN (5) with the fast and slow species spli ing (20), the entropy production rate of the elementary and of the e ective dynamics are plo ed in Fig. 2. Whilê

D. Gibbs Free Energy
e Gibbs free energy of ideal dilute solution is given by with a ≡ i a i . It is the proper thermodynamic potential of the elementary dynamics of closed CRNs since it has been proven in Ref. [15] that i) d t G = −T Σ ≤ 0 and ii) G ≥ G eq ≡ G(z eq ). We now explain why the following Gibbs free energy is the proper thermodynamic potential of the e ective dynamics: First, we notice that in generalĜ(p) G(z). is means thatĜ does not estimate the exact free energy of CRNs. However, the time derivative ofĜ according to the e ective dynamics (19) reads where we used the de nition of the Gibbs free energy of the e ective reactions given in Eq. (36) and the entropy production rate (41). e variation of free energy in a time interval [0, t] thus satis es as long as the time scale separation hypothesis holds and, hence,ˆ Σ(t) = Σ(t). Second, one can prove thatĜ(p(t)) ≥Ĝ eq ≡Ĝ(p eq ). Indeed, considerĜ eq =μ eq · p eq − RT p eq . (47) As mentioned in Subs. IV B,μ eq belongs to the cokernel ofŜ and, therefore, it can be expressed as a linear combination of the conservation lawsl ζ :μ eq = f ζl ζ .
(48) us, the termμ eq · p eq in Eq. (47) satis eŝ where we used d tL ζ = 0. is allows us to write the di erence betweenĜ(p(t)) andĜ eq as a relative entropŷ proving thatĜ(p(t)) ≥Ĝ eq . In summary,Ĝ in Eq. (44) is the proper thermodynamic potential of the e ective dynamics because of Eq. (45) and Eq. (50). It also correctly quanti es the variation of the free energy of the whole CRN because of Eq. (46).

V. DYNAMICS OF OPEN CRNS
In open CRNs some species are exchanged with external reservoirs called chemostats. We discuss in parallel the elementary dynamics and the e ective one. We consider that the concentrations of chemosta ed species are either xed or slowly driven by the chemostats. ese species are therefore treated as slow-evolving species.
We thus split P into two disjoint subsets: the internal species X and the chemosta ed species Y . e former (as well as the fast species) evolve in time only because of the chemical reactions and their rate equation is unchanged. e la er evolve in time because of the chemical reactions and of ma er ows with the chemostats. ese are accounted in the rate equation by introducing the exchange current I (t).
For the elementary dynamics of the chemosta ed species, the rate equation is given by while for the e ective dynamics it reads [S] chem.
[P] chem. Here we applied the spli ing P = (X , Y ) to the substoichiometric matrix S P and the e ective stoichiometric matrixŜ, as well as the vector p = (x, ). Note that Eq. (52) and (53) are merely de nitions for the exchange current I (t). At the level of the e ective dynamics, a slow variation of (t) corresponds to a slow variation of I (t). Operationally, one may therefore equally well control I (t) and determine (t) through Eq. (53). Both types of external control, on (t) or I (t), can thus be considered.
Example. For the CRN (5) and given the spli ing between fast and slow species in (20), we chemostat the free substrate S and the product P. e adsorbed substrate S m is now the only internal species. Using the e ective stoichiometric matrix in Eq. (22), we obtain the following two substoichiometric matrix In Fig. 3, we compare the elementary dynamics of the CRN (5) for the internal species and the e ective one.

A. Topological Properties
When CRNs are open, some conservation laws do not correspond anymore to conserved quantities. Hence, we split the set of conservation laws into two disjoint subsets: the unbroken and the broken conservation laws. e unbroken conservation laws are those with null entries for the chemosta ed species. For the elementary dynamics, the conservation laws { ξ } of the fast species (see Subs. III A) are unbroken by de nition. en, we represent with the vectors { ζ u } ⊆ { ζ } the other unbroken conservation laws (if any). Indeed, the quantities L ζ u (z(t)) ≡ ζ u · z(t), as well as L ξ (z(t)) ≡ ξ · z(t), are conserved even if the CRN is open: (56) e broken conservation laws are all the other conservation laws, e corresponding quantities Because of the correspondence in Eq. (28), the unbroken/broken conservation laws of the e ective dynamics are given, respectively, bŷ Chemosta ing a species does not always break a conservation law [15][16][17]27]. We thus distinguish the set of controlled species Y p ⊆ Y breaking the conservation laws from the others Y f = Y \ Y p . Note that the number of Y p species is equal to the number of broken conservation laws.
is allows us to introduce the so-called moieties. ey represent parts of (or entire) molecules which are exchanged with the environment through the chemostats. For the elementary dynamics, their concentration vector is speci ed as while for the e ective dynamics it is given bŷ We introduced the vectors of broken conserved quantities L br (z(t)) = (. . . , ζ b · z(t), . . . ) andL br (p(t)) = (. . . ,ˆ ζ b · p(t), . . . ) , and the matrix M with entries { ζ b α } α ∈Y p (see Ref. [15][16][17] for details). e matrix M is square and nonsingular, and it can be inverted obtaining M −1 . Comparing Eq. (59) and Eq. (60), one notices that m(z(t)) =m(p(t)) if for every broken conservation law, i.e., the slow-evolving species are much more abundant than the fast-evolving species. We stress nally that the number of moieties is equal to the number of Y p species.
Example. For the CRN (5) and given the spli ing between fast and slow species in (20), the conservation lawˆ S in Eq. (29) of the e ective dynamics is broken when the free product P is chemosta ed. If we now chemostat the free substrate S, no conservation law is broken. According to this sequence of chemosta ing, the species P belongs to the set Y p and the matrix M is simply the scalar e corresponding moiety, represents the total concentration of product exchanged between the two chemostats.

VI. THERMODYNAMICS OF OPEN CRNS
We consider now the thermodynamic description of open CRNs. e semigrand Gibbs free energy represents the proper thermodynamic potential for the elementary dynamics [15,16] since i) d t G = −T Σ ≤ 0 for autonomous detailed balanced systems and ii) G ≥ G eq = G(z eq ). In analogy to equilibrium thermodynamics, G is dened from the Gibbs free energy (43) e entropy production rate Σ is speci ed in Eq. (40). e the driving work rate w driv accounts for the time dependent manipulation of the chemical potential of the Y p chemostats, e nonconservative work rate w nc quanti es the energetic cost of sustaining uxes of chemical species among the chemostats, by means of the force In other words, this is the force keeping the system out of equilibrium (see also App. A).
For the e ective dynamics, we introduce the following semigrand Gibbs free energy: withm(p) given in Eq. (60). Note that µ Y p (t) in Eq. (63) and in Eq. (67) are exactly the same since they are imposed by the same chemostats. We now show thatĜ(p) is the proper thermodynamic potential of the e ective dynamics of open CRNs.
In generalĜ(p) G(z) sinceĜ(p) G(z). However, we will now see that their variation in time can be very similar. According to the e ective dynamics (19) and taking into account the exchange current I (t) through Eq. (53), the evolution ofĜ(p) is given by d tĜ (p(t)) = −Tˆ Σ(t) +ˆ w driv (t) +ˆ w nc (t) .
(68) e entropy production rateˆ Σ of Eq. (41) satis esˆ Σ = Σ as long as the time scale separation hypothesis holds (see Subs. IV C). e driving work rate at the e ective level, corresponds to w driv of the elementary dynamics when m(p(t)) = m(z(t)). is occurs whenL ζ b = L ζ b (see Subs. V A). e nonconservative work rateˆ w nc has the same expression as the one for the elementary dynamics speci ed in Eq. (66). However, the numerical values ofˆ w nc and w nc can be di erent because di erent currents I (t) might be necessary to set the same concentrations of the chemosta ed species for the elementary and e ective dynamics if the time scale separation hypothesis does not hold perfectly. We can thus conclude that the variation of the semigrand Gibbs free energy in a time interval [0, t] satis es as long as the time scale separation hypothesis holds (granting . If the open system is autonomous (no driving is performed, w driv = 0) and detailed balance (all the chemosta ed species break a conservation law,ˆ w nc = 0), then d tĜ = −Tˆ Σ ≤ 0. Finally, we show thatĜ(p(t)) ≥Ĝ eq ≡Ĝ(p eq ), where the equilibrium state p eq is set by the Y p chemostats (see App. (A)). To this aim, consider thatμ eq belongs to the cokernel ofŜ and can thus be expressed as a linear combination of the conservation laws: where now we distinguish between unbroken and broken conservation laws. us,μ eq · p eq satis es noŵ where we used d tL ζ u = 0 and the last step is discussed in App. (A). is allows us to writeĜ eq aŝ and, consequently, the di erence betweenĜ(p(t)) andĜ eq as a relative entropŷ G(p(t)) −Ĝ eq = RT L(p(t) p eq ) ≥ 0 .
FIG. 4. ermodynamic quantities of the elementary dynamics and e ective dynamics in Fig. 3.
We use RT (k +1 ) 2 /k +2 as units of measure for the thermodynamic quantities.
Example. For the CRN (5) and given the spli ing between fast and slow species in (20), we compare the thermodynamic quantities of the elementary and of the e ective dynamics in Fig. 4.

VII. EFFECTIVE NETWORKS WITHOUT ELEMENTARY COUNTERPART
Our thermodynamic framework can be applied directly to e ective networks even if the full elementary description is not available. Indeed, all the expressions (e.g., Eqs. (36), (41), (44) and (67)) of the e ective thermodynamic quantities require the knowledge of only the e ective dynamics (19) and the standard chemical potentialsμ • of the slow species. is is the key result of our work.
However, our framework should be applied only to e ective models which are compatible with an underlying elementary network satisfying the time scale separation hypothesis. Indeed, we proved the thermodynamic consistency only in this case. On the one hand, one should have some physical evidence supporting that the models result from elementary reactions satisfying the time scale separation hypothesis. On the other hand, the e ective models must exhibit the following three properties to be thermodynamically consistent.
First, in absence of any chemosta ed species (i.e., closed CRNs), the e ective model must relax to an equilibrium (30).
is is fundamental for energetic considerations. Physically, it means that every system has to reach an equilibrium when there are no energy sources (the chemostats) balancing the dissipation. Mathematically, it means that the e ective reactions must be reversible and e ective current vector ψ(p) has to admit a nontrivial equilibrium steady-state, i.e., ∃p eq 0 such that ψ(p eq ) = 0. If the e ective model relaxes to a nonequilibrium steady-state without chemosta ed species, then it means that it does not account for hidden sources of energy and any energetic consideration becomes meaningless. Second, provided that the time scale separation hypothesis holds, there must be no cycles in the e ective network as we proved in Subs. III A. ird, provided that the time scale separation hypothesis holds, the e ective entropy production rate (41) must be greater than or equal to zero. e propertŷ Σ(t) ≥ 0 is not granted anymore if the cycle current vector ψ is not speci ed according to the coarse graining procedure discussed in Sec. III. In this respect, an e ective dynamics giving rise toˆ Σ(t) < 0 is necessarily thermodynamically inconsistent.
We now consider two examples, one that can and the other that cannot be characterized using our framework.

A. Example 1
Consider the following reactions satisfying the dynamics We then verify that the dynamical system (76) does not admit cycles: the stoichiometric matrix in Eq. (78) has no right null eigenvectors. e le null eigenvector ofŜ is the conservation laws = (1, 1, 1) corresponding to the total concentration, L = · p = [S] + [I] + [P]. e equilibrium state of Eq. (76) can be identi ed by solving the system of equations ψ(p eq ) = 0 and L = · p eq = · p(0).
Assuming that the standard chemical potentials µ • S , µ • I and µ • P are known, the second step is to apply the expression of the thermodynamic quantities given in Sec. IV and VI. For For simplicity, we use quantities scaled by some arbitrary reference unit.
instance, the dissipation can be quanti ed with the entropy production rate of Eq. (41) which reads now By solving the rate equation (76) for a speci c set of parameters and initial condition, one can compute the entropy production rate (79) which is shown in Fig. 5. Notice that Σ(t) ≥ 0 for every time step of the dynamics. With the same strategy, namely, applying the e ective expressions given in this work, one can compute also other thermodynamic quantities. Consider now the following elementary mechanism In the limit of fast-evolving enzymes (E 1 , E 2 and E 3 ) and complexes (E 1 S, E 2 I and E 3 I), the e ective dynamics for [S], [I] and [P] provided by the coarse graining procedure discussed in Sec. III is consistent with the dynamical system (76). We show in Fig. 5 the entropy production rate (40) of this elementary mechanism for two di erent set of initial concentrations of the enzymes and the complexes. In the rst case, the concentrations are not small enough and the time scale separation hypothesis does not hold. Hence, the entropy production rate at the elementary level does not correspond to that of the e ective model (75). In the second case, the time scale separation hypothesis holds and the entropy production rate at the elementary level is well approximated by the effective one. e initial di erence between the two is due to the relaxation of the initial state of the elementary dynamics to the corresponding quasi-steady-state for the fast-evolving species.

B. Example 2
Consider now the model of gene regulation provided in Ref. [28]. It represents the synthesis of two proteins A and B via the expression of the two genes G A and G B . Each protein promotes its synthesis and represses the synthesis of the other. en, the proteins degrade. e corresponding chemical reaction network is Here, A , B , a 1 , a 2 , S, b 1 , b 2 , k A and k B are parameters of the model. e concentration of the proteins A and B are the only dynamical variables. is e ective model (82) is designed in such a way that it cannot equilibrate. Indeed, the degradation reactions are irreversible and the currents cannot vanish. e model always relaxes towards a nonequilibrium steady state, but the major issue is that energy sources preventing the equilibration cannot be accounted for. While our thermodynamic quantities can be formally de ned for this model (they just require the dynamical system and the standard chemical potential), they are meaningless.
We note that the use of irreversible reactions and/or nonvanishing currents is a very common feature of kinetic models for biology and does not preclude per se a consistent thermodynamic analysis as recently illustrated for the irreversible Michaelis-Menten enzymatic scheme [29].

VIII. CONCLUSIONS
In this work, we developed a thermodynamic theory for e ective (non-mass-action) models of both closed and open CRNs. We focused here only on deterministic models. Our approach provides the exact thermodynamic quantities when the e ective models result from underlying elementary (mass-action) networks satisfying the time scale separation hypothesis.
is was proven by exploiting the topological properties of the CRNs. Our time scale separation hypothesis can be considered as a zero-order expansion in the ratio between the fast and the slow time scale. Exploring higher order corrections and whether they can be used to bound deviations in entropy production rates is le for future work.
Similar approaches might be employed in other frameworks. First, the topological properties could be used to derive a thermodynamically consistent coarse-graining of stochastic CRNs. Second, one can exploit weakly broken conservation laws to collect di erent species into e ective mesostates. ese mesostates may then satisfy closed evolution equations. For example, during a catalytic process the complex enzyme-substrate is transformed in many di erent species that can be considered as a unique mesostate using the conservation of the total concentration of the enzyme. We leave also these points to future investigations. the chemical potentials µ Y f which, in general, do not satisfy Eq. (A8). As a consequence, the nonconservative forces with µ α,eq given in Eq. (A7) and (A8). Each entries in F represents the di erence between the chemical potential set by a chemostat and the equilibrium one. erefore, it represents the force keeping the system out of equilibrium.
Note that, thanks of Eq. (A7), we can write f br = µ Y p ·M −1 . We use this property and the de nition of the moieties (60) in the last step of Eq. (72).