Thermodynamics of non-elementary chemical reaction networks

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 often modeled with effective schemes obeying non-mass action kinetics, e.g., the Michaelis and Menten mechanism [1, 2], Hill functions [3] and many others [4]. These models can be thought of as resulting from a coarse graining of elementary reactions by eliminating fast-evolving species and their use is well established [5–11].

However, characterizing the energetics at the level of these effective schemes remains an open challenge especially out-of-equilibrium. Inspired by developments in stochastic thermodynamics [12–14], the nonequilibrium thermodynamics of open chemical networks of elementary reactions is nowadays well established [15–17]. Recently, a thermodynamically consistent coarse graining strategy of the dynamics of biocatalysts was proposed in reference [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 section 2 and of the coarse grained ones in section 3. We then proceed by building the corresponding thermodynamics in section 4. The dynamics of the open networks is considered in section 5 and its thermodynamics in section 6. In section 7 we discuss how our previous findings can be used to characterize the thermodynamics of non-elementary networks when the full elementary description is not known. In section 8 we summarize our results and discuss their implications.

We consider here chemical reaction networks (CRNs) composed of chemical species α = (..., α, ...)^⊺ undergoing elementary reactions [19],

$$\begin{equation}\boldsymbol{\alpha }\cdot {\boldsymbol{\nu }}_{+\rho }\underset{{k}_{-\rho }}{\overset{{k}_{+\rho }}{{\rightleftharpoons}}}\boldsymbol{\alpha }\cdot {\boldsymbol{\nu }}_{-\rho },\end{equation} \tag{ 1 }$$

with ν_±ρ (resp. k_±ρ) the vector of the stoichiometric coefficients (resp. the kinetic constant) of the forward/backward reaction ρ. The state of deterministic CRNs is specified by the concentration vector z(t) = (..., [α](t), ...)^⊺ of all the chemical species α. For closed CRNs, its dynamics follows the rate equation

$$\begin{equation}{\mathrm{d}}_{t}\boldsymbol{z}\left(t\right)=\mathbb{S}\boldsymbol{j}\left(\boldsymbol{z}\left(t\right)\right),\end{equation} \tag{ 2 }$$

where we introduce the stoichiometric matrix $\mathbb{S}$ and the current vector j(z). The stoichiometric matrix $\mathbb{S}$ codifies the topology of the network. Each ρ column ${\mathbb{S}}_{\rho }$ of the stoichiometric matrix specifies the net variation of the number of molecules for each species undergoing the ρ elementary reaction (1), ${\mathbb{S}}_{\rho }={\boldsymbol{\nu }}_{-\rho }-{\boldsymbol{\nu }}_{+\rho }$. The current vector $\boldsymbol{j}\left(\boldsymbol{z}\right)={\left(\dots ,{j}^{\rho }\left(\boldsymbol{z}\right),\dots \right)}^{{\intercal}}$ specifies the net reaction current for every ρ elementary reaction (1) as the difference between the forward j^+ρ(z) and backward reaction current j^−ρ(z)

$$\begin{equation}{j}^{\rho }\left(\boldsymbol{z}\right)={j}^{+\rho }\left(\boldsymbol{z}\right)-{j}^{-\rho }\left(\boldsymbol{z}\right),\end{equation} \tag{ 3 }$$

with j^±ρ(z) satisfying the mass-action kinetics [20–22],

$$\begin{equation}{j}^{{\pm}\rho }\left(\boldsymbol{z}\right)={k}_{{\pm}\rho }{\boldsymbol{z}}^{{\boldsymbol{\nu }}_{{\pm}\rho }}.\end{equation} \tag{ 4 }$$

Note that we use ${\boldsymbol{a}}^{\boldsymbol{b}}={\prod }_{i}{\;a}_{i}^{{b}_{i}}$. In the following, we will refer to equation (2) as the elementary dynamics of closed CRNs.

Example. In all the manuscript, we illustrate our findings using a modified version of the model discussed in reference [10]. It represents the transformation of two identical substrates S into one product P. The process is catalyzed by a membrane enzyme E which interacts with the substrate after it is adsorbed S_m by the membrane:

$$\begin{equation}\begin{aligned}\hfill \text{S}\underset{{k}_{-1}}{\overset{{k}_{+1}}{{\rightleftharpoons}}}{\text{S}}_{\text{m}}\\ \hfill {\text{S}}_{\text{m}}+\;\text{E}\underset{{k}_{-2}}{\overset{{k}_{+2}}{{\rightleftharpoons}}}\text{ES}\\ \hfill {\text{S}}_{\text{m}}+\text{ES}\underset{{k}_{-3}}{\overset{{k}_{+3}}{{\rightleftharpoons}}}\text{ESS}\\ \hfill \text{ESS}\underset{{k}_{-4}}{\overset{{k}_{+4}}{{\rightleftharpoons}}}\text{E}\;+\;\text{P}\end{aligned}\end{equation} \tag{ 5 }$$

For this model, we introduced here the concentration vector and the current vector,

$$\begin{equation}\boldsymbol{z}=\left(\begin{matrix}\hfill \left[\text{E}\right]\hfill \\ \hfill \left[\text{ES}\right]\hfill \\ \hfill \left[\text{ESS}\right]\hfill \\ \hfill \left[\text{P}\right]\hfill \\ \hfill \left[\text{S}\right]\hfill \\ \hfill \left[{\text{S}}_{\text{m}}\right]\hfill \end{matrix}\right),\quad \boldsymbol{j}\left(\boldsymbol{z}\right)=\left(\begin{matrix}\hfill {k}_{+1}\left[\text{S}\right]-{k}_{-1}\left[{\text{S}}_{\text{m}}\right]\hfill \\ \hfill {k}_{+2}\left[{\text{S}}_{\text{m}}\right]\left[\text{E}\right]-{k}_{-2}\left[\text{ES}\right]\hfill \\ \hfill {k}_{+3}\left[{\text{S}}_{\text{m}}\right]\left[\text{ES}\right]-{k}_{-3}\left[\text{ESS}\right]\hfill \\ \hfill {k}_{+4}\left[\text{ESS}\right]-{k}_{-4}\left[\text{E}\right]\left[\text{P}\right]\hfill \end{matrix}\right),\end{equation} \tag{ 6 }$$

as well as the stoichiometric matrix

2.1. Topological properties

As shown in reference [15, 23], the linear independent vectors {ℓ^λ} of the cokernel of the stoichiometric matrix

$$\begin{equation}{\boldsymbol{\ell }}^{\lambda }\cdot \mathbb{S}=0\end{equation} \tag{ 8 }$$

are the conservation laws of equation (2). Indeed, for each vector ℓ^λ, the scalar L^λ(z(t)) ≡ ℓ^λ ⋅ z(t) is a conserved quantity, i.e., ${\mathrm{d}}_{t}{L}^{\lambda }\left(\boldsymbol{z}\left(t\right)\right)={\boldsymbol{\ell }}^{\lambda }\cdot \mathbb{S}\boldsymbol{j}\left(\boldsymbol{z}\left(t\right)\right)=0$. The linear independent vectors $\left\{{\boldsymbol{c}}_{\iota }\right\}$ of the kernel of the stoichiometric matrix

$$\begin{equation}\mathbb{S}{\boldsymbol{c}}_{\iota }=0\end{equation} \tag{ 9 }$$

are the internal cycles of equation (2). They 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., $\mathbb{S}\overline{\boldsymbol{j}}=0$ with $\overline{\boldsymbol{j}}\equiv {\boldsymbol{c}}_{\iota }{\psi }^{\iota }$. 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., a_ibⁱ = ∑_ia_ibⁱ and a_ibⁱ = ∑_ia_ibⁱ.

Example. For the CRN (5), there are two conservation laws,

with a clear physical interpretation. The total concentration of the enzyme is given by L^E = ℓ^E ⋅ z = [E] + [ES] + [ESS], while the total concentration of the substrate is given by L^S = ℓ^S ⋅ z = [ES] + 2[ESS] + 2[P] + [S] + [S_m]. The CRN (5) has no internal cycles.

2.2. Equilibrium

Closed CRNs must be detailed balanced, namely the rate equation (2) admits an equilibrium steady-state z_eq characterized by vanishing reaction currents,

$$\begin{equation}\boldsymbol{j}\left({\boldsymbol{z}}_{\text{eq}}\right)=0.\end{equation} \tag{ 11 }$$

This, together with mass-action kinetics, implies the so-called local detailed balance condition for the kinetic constants k_±ρ of the chemical reactions

$$\begin{equation}\frac{{k}_{+\rho }}{{k}_{-\rho }}={\boldsymbol{z}}_{\text{eq}}^{{\mathbb{S}}_{\rho }}.\end{equation} \tag{ 12 }$$

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].

When the chemical species evolve over two different time scales, they can be divided into two disjoint subsets: the fast-evolving species Q and the slow-evolving species P [7, 8]. We apply the same splitting to the stoichiometric matrix

$$\begin{equation}\mathbb{S}=\left(\begin{matrix}\hfill {\mathbb{S}}^{Q}\hfill \\ \hfill {\mathbb{S}}^{P}\hfill \end{matrix}\right)\end{equation} \tag{ 13 }$$

and to the concentration vector z = (q, p). This also allows us to split the rate equation (2) into

$$\begin{align}\hfill {\mathrm{d}}_{t}\boldsymbol{q}\left(t\right)={\mathbb{S}}^{Q}\boldsymbol{j}\left(\boldsymbol{q}\left(t\right),\boldsymbol{p}\left(t\right)\right),\end{align} \tag{ 14 }$$

$$\begin{align}\hfill {\mathrm{d}}_{t}\boldsymbol{p}\left(t\right)={\mathbb{S}}^{P}\boldsymbol{j}\left(\boldsymbol{q}\left(t\right),\boldsymbol{p}\left(t\right)\right).\end{align} \tag{ 15 }$$

The 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. The first is the existence, for every concentration p of the slow species, of the concentration vector $\hat{\boldsymbol{q}}\left(\boldsymbol{p}\right)$ such that the current

$$\begin{equation}\hat{\boldsymbol{j}}\left(\boldsymbol{p}\right)\equiv \boldsymbol{j}\left(\hat{\boldsymbol{q}}\left(\boldsymbol{p}\right),\boldsymbol{p}\right)\end{equation} \tag{ 16 }$$

is a steady-state current of equation (14) or, equivalently, $\hat{\boldsymbol{j}}\left(\boldsymbol{p}\right)\in \mathrm{ker}\enspace {\mathbb{S}}^{Q}$. This is a topological property of the CRN. The 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., $\boldsymbol{q}\left(t\right)=\hat{\boldsymbol{q}}\left(\boldsymbol{p}\left(t\right)\right)$. This means physically that the concentrations of the fast species relax instantaneously to the steady-state corresponding to the frozen values of p(t). This 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. This 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 $\hat{\boldsymbol{j}}\left(\boldsymbol{p}\right)$ of equation (14). It can now be written as a linear combination of the right null eigenvectors ${\mathbb{S}}^{Q}{\boldsymbol{c}}_{\gamma }=0$,

$$\begin{equation}\hat{\boldsymbol{j}}\left(\boldsymbol{p}\right)={\boldsymbol{c}}_{\gamma }{\psi }^{\gamma }\left(\boldsymbol{p}\right),\end{equation} \tag{ 17 }$$

with p dependent coefficients {ψ^γ}. The specific expression of ψ^γ(p) is not discussed here and is not fundamental in what follows. When the dynamics of the fast species at fixed concentrations of the slow ones is linear, the problem is solved in references [7, 18] using a diagrammatic method developed in references [24, 25]. Note that {c_γ} includes all the internal cycles ${\boldsymbol{c}}_{\iota }$ of equation (9) and, in general, other vectors {c_ɛ} called pseudo-emergent cycles [15, 18, 23]. The former characterize a sequence of reactions which upon completion leaves all species (q, p) unchanged, while the latter only leaves the fast species unchanged. Each coefficient ψ^γ(p) represents a current along the cycle γ. By employing the steady-state current vector (17) in equation (15), we obtain a closed dynamical equation for the slow species,

$$\begin{equation}{\mathrm{d}}_{t}\boldsymbol{p}\left(t\right)={\mathbb{S}}^{P}\hat{\boldsymbol{j}}\left(\boldsymbol{p}\left(t\right)\right)={\mathbb{S}}^{P}{\boldsymbol{c}}_{\varepsilon }{\psi }^{\varepsilon }\left(\boldsymbol{p}\left(t\right)\right),\end{equation} \tag{ 18 }$$

where we used ${\mathbb{S}}^{P}{\boldsymbol{c}}_{\iota }=0$ because of equation (9). This coarse grained dynamics will be denoted as effective for convenience. It can be rewritten in a more compact way,

$$\begin{equation}{\mathrm{d}}_{t}\boldsymbol{p}\left(t\right)=\hat{\mathbb{S}}\boldsymbol{\psi }\left(\boldsymbol{p}\left(t\right)\right),\end{equation} \tag{ 19 }$$

introducing the effective stoichiometric matrix $\hat{\mathbb{S}}$ and the effective current vector ψ(p). The effective stoichiometric matrix $\hat{\mathbb{S}}$ codifies the net stoichiometry of the slow species along the pseudo-emergent cycles. Each ɛ column ${\hat{\mathbb{S}}}_{\varepsilon }={\mathbb{S}}^{P}{\boldsymbol{c}}_{\varepsilon }$ specifies the net variation of the number of molecules for each slow species along the ɛ pseudo-emergent cycle, i.e., the stoichiometry of the effective reactions. The effective current vector ψ(p) collects the pseudo-emergent cycle currents $\boldsymbol{\psi }\left(\boldsymbol{p}\right)={\left(\dots ,{\psi }^{\varepsilon }\left(\boldsymbol{p}\right),\dots \right)}^{{\intercal}}$. Unlike j(z) in equation (2), ψ(p) does not, in general, satisfy mass-action kinetics.

Note that, if there were no pseudo-emergent cycles, $\left(\hat{\boldsymbol{q}}\left(\boldsymbol{p}\right),\boldsymbol{p}\right)$ would be an equilibrium condition of the elementary dynamics (2). The corresponding effective dynamics for the slow species would become trivial: d_tp(t) = 0.

Example. For the CRN (5), we split the chemical species as follows

$$\begin{equation}\underset{Q}{\underbrace{\left\{\text{E},\text{ES},\text{ESS}\right\}}}\cup \underset{P}{\underbrace{\left\{\text{P},\text{S},{\text{S}}_{\text{m}}\right\}}}.\end{equation} \tag{ 20 }$$

Underlying this separation is the assumption that the concentration of the enzyme species changes much more quickly. Since there are no internal cycles of $\mathbb{S}$ of equation (7), all the right null eigenvectors of ${\mathbb{S}}^{Q}$,

are pseudo-emergent cycles. Thus, the effective stoichiometric matrix of equation (19) is given by

Each column of $\hat{\mathbb{S}}$ specifies the net stoichiometry of the following effective reactions for the slow species:

$$\begin{equation}\begin{aligned}\hfill \text{S}\underset{}{\overset{\text{ads}}{{\rightleftharpoons}}}{\text{S}}_{\text{m}}\\ \hfill 2{\text{S}}_{\text{m}}\underset{}{\overset{\text{enz}}{{\rightleftharpoons}}}\text{P}\end{aligned}\end{equation} \tag{ 23 }$$

In figure 1, we compare the elementary dynamics of the CRN (5) for the slow species and the effective one obtained by solving equations (2) and (19), respectively. For this case, the effective current vector ψ(p) is computed according to the procedure introduced in references [24, 25]:

$$\begin{align}\hfill & {\psi }_{\text{ads}}\left(\boldsymbol{p}\right)={k}_{+1}\left[\text{S}\right]-{k}_{-1}\left[{\text{S}}_{\text{m}}\right],\hfill \end{align} \tag{ 24 }$$

$$\begin{align}\hfill & {\psi }_{\text{enz}}\left(\boldsymbol{p}\right)=\frac{{L}^{\text{E}}}{\mathcal{D}}\left\{{k}_{+2}{k}_{+3}{k}_{+4}{\left[{\text{S}}_{\text{m}}\right]}^{2}-{k}_{-2}{k}_{-3}{k}_{-4}\left[\text{P}\right]\right\},\hfill \end{align} \tag{ 25 }$$

with L^E the total concentration of the enzyme and

$$\begin{equation}\begin{aligned}\hfill \mathcal{D}& ={k}_{+2}{k}_{+3}{\left[{\text{S}}_{\text{m}}\right]}^{2}+{k}_{+2}{k}_{+4}\left[{\text{S}}_{\text{m}}\right]+{k}_{+2}{k}_{-3}\left[{\text{S}}_{\text{m}}\right]+{k}_{+3}{k}_{+4}\left[{\text{S}}_{\text{m}}\right]+{k}_{-2}{k}_{-4}\left[\text{P}\right]+{k}_{-3}{k}_{-4}\left[\text{P}\right]\hfill \\ \hfill & \quad +{k}_{+3}{k}_{-4}\left[{\text{S}}_{\text{m}}\right]\left[\text{P}\right]+{k}_{-2}{k}_{-3}+{k}_{-2}{k}_{+4}.\hfill \end{aligned}\end{equation} \tag{ 26 }$$

Zoom In Zoom Out Reset image size

3.1. Topological properties

The conservation laws of the effective dynamics (19) are defined as the linear independent vectors of the cokernel of the effective stoichiometric matrix

$$\begin{equation}{\hat{\boldsymbol{\ell }}}^{\zeta }\cdot \hat{\mathbb{S}}=0,\end{equation} \tag{ 27 }$$

as in equation (8). Indeed, each scalar ${\hat{L}}^{\zeta }\left(\boldsymbol{p}\left(t\right)\right)\equiv {\hat{\boldsymbol{\ell }}}^{\zeta }\cdot \boldsymbol{p}\left(t\right)$ is a conserved quantity of the effective dynamics, i.e., ${\mathrm{d}}_{t}{\hat{L}}^{\zeta }\left(\boldsymbol{p}\left(t\right)\right)={\hat{\boldsymbol{\ell }}}^{\zeta }\cdot \hat{\mathbb{S}}\boldsymbol{\psi }\left(\boldsymbol{p}\left(t\right)\right)=0$.

The conservation laws $\left\{{\hat{\boldsymbol{\ell }}}^{\zeta }\right\}$ of (27) can be specified in terms of the conservation laws {ℓ^λ} of equation (8) as follows. First, we split {ℓ^λ} into two disjoint subsets: the conservation laws with null entries for the slow species {ℓ^ξ} and the other conservation laws {ℓ^ζ}. Second, we consider the projection operator $\mathbb{P}$ onto the space of the slow species: $\mathbb{P}\boldsymbol{z}=\boldsymbol{p}$. Third, we identify each conservation law ${\hat{\boldsymbol{\ell }}}^{\zeta }$ of equation (27) as follows

$$\begin{equation}{\hat{\boldsymbol{\ell }}}^{\zeta }=\mathbb{P}{\boldsymbol{\ell }}^{\zeta }.\end{equation} \tag{ 28 }$$

Indeed $\mathbb{P}{\boldsymbol{\ell }}^{\zeta }\cdot \hat{\mathbb{S}}=0$. This can be verified considering that $\mathbb{P}{\boldsymbol{\ell }}^{\zeta }\cdot {\hat{\mathbb{S}}}_{\varepsilon }={\boldsymbol{\ell }}^{\zeta }\cdot \mathbb{S}{\boldsymbol{c}}_{\varepsilon }=0$, where we only used the definition of pseudo-emergent cycle and conservation law, i.e., ${\mathbb{S}}^{Q}{\boldsymbol{c}}_{\varepsilon }=0$ and ${\boldsymbol{\ell }}^{\zeta }\cdot \mathbb{S}=0$.

Note that all the conservation laws of the effective stoichiometry matrix $\hat{\mathbb{S}}$ can be written as in equation (28). This follows from the rank nullity theorem for $\mathbb{S}$ and $\hat{\mathbb{S}}$, and the absence of cycles for $\hat{\mathbb{S}}$. Indeed, suppose that there is a vector ϕ ≠ 0 such that $\hat{\mathbb{S}}\boldsymbol{\phi }=0$. This means that ${\mathbb{S}}^{P}{\boldsymbol{c}}_{\varepsilon }{\phi }^{\varepsilon }=0$ and, consequently, c_ɛϕ^ɛ is a right null eigenvector of both ${\mathbb{S}}^{P}$ and ${\mathbb{S}}^{Q}$, i.e., an internal cycle. Since c_ɛϕ^ɛ is a linear combination of only pseudo-emergent cycles {c_ɛ}, the effectively stoichiometric matrix cannot admit cycles.

Example. For the CRN (5) and given the splitting between fast and slow species in (20), the only conservation law of $\mathbb{S}$ in (7) with no null entries for the slow species is ℓ^S of equation (10). The corresponding conservation law of the effective dynamics is

Given $\hat{\mathbb{S}}$ in equation (22), one easily verifies that ${\hat{\boldsymbol{\ell }}}^{\text{S}}\cdot \hat{\mathbb{S}}=0$. The corresponding quantity ${\hat{L}}^{\text{S}}={\hat{\boldsymbol{\ell }}}^{\text{S}}\cdot \boldsymbol{p}$ has the clear physical interpretation of the total concentration of substrate: ${\hat{L}}^{\text{S}}=2\left[\text{P}\right]+\left[\text{S}\right]+\left[{\text{S}}_{\text{m}}\right]$.

3.2. Equilibrium

The equilibrium condition p_eq of the effective dynamics (19) is defined by vanishing cycle currents,

$$\begin{equation}\boldsymbol{\psi }\left({\boldsymbol{p}}_{\text{eq}}\right)=0.\end{equation} \tag{ 30 }$$

Its existence is granted by the local detailed balance condition (12). Indeed, the state $\left(\hat{\boldsymbol{q}}\left({\boldsymbol{p}}_{\text{eq}}\right),{\boldsymbol{p}}_{\text{eq}}\right)$ is an equilibrium of the elementary dynamics (2) by definition: it is a steady-state of the rate equation (2) which admits only equilibrium steady-states.

Let us discuss now the equilibrium state to which the elementary and the effective dynamics relax given a unique initial condition z(0) = (q(0), p(0)). The concentrations of the chemical species at equilibrium depend on the value of the conserved quantities. Therefore, to have similar concentrations for the slow species at equilibrium, also the values of L^ζ(z(0)) = l^ζ ⋅ z(0) and ${\hat{L}}^{\zeta }\left(\boldsymbol{p}\left(0\right)\right)={\hat{\boldsymbol{l}}}^{\zeta }\cdot \boldsymbol{p}\left(0\right)$ have to be similar. This means that the concentration of the fast species involved in the ζ conservation laws must be much smaller than the concentration of the slow species. This is consistent with the time scale separation hypothesis requiring much higher abundances for the slow-evolving than for the fast-evolving species.

We now derive the thermodynamics for the effective 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 effective dynamics, without any knowledge of the elementary dynamics. This is the key result of our work.

4.1. Local detailed balance

The 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 specified by their equilibrium form but expressed in terms of nonequilibrium concentrations. The vector of chemical potentials is thus given by

$$\begin{equation}\boldsymbol{\mu }\left(\boldsymbol{z}\right)={\boldsymbol{\mu }}^{{\circ}}+RT\enspace \mathrm{ln}\left(\boldsymbol{z}\right),\end{equation} \tag{ 31 }$$

with μ^◦ the vector of the standard chemical potentials and R the gas constant. The 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 μ^◦):

$$\begin{equation}\frac{{k}_{+\rho }}{{k}_{-\rho }}={\boldsymbol{z}}_{\text{eq}}^{{\mathbb{S}}_{\rho }}=\mathrm{exp}\left(-\frac{{\boldsymbol{\mu }}^{{\circ}}\cdot {\mathbb{S}}_{\rho }}{RT}\right).\end{equation} \tag{ 32 }$$

We now formulate the local detailed balance condition for the effective dynamics (19). This is done by taking the product of the ratio k_+ρ/k_−ρ along each pseudo-emergent cycle ɛ

$$\begin{equation}\prod _{\rho }{\left(\frac{{k}_{+\rho }}{{k}_{-\rho }}\right)}^{{c}_{\varepsilon }^{\rho }}={\boldsymbol{p}}_{\text{eq}}^{{\hat{\mathbb{S}}}_{\varepsilon }}=\mathrm{exp}\left(-\frac{{\hat{\boldsymbol{\mu }}}^{{\circ}}\cdot {\hat{\mathbb{S}}}_{\varepsilon }}{RT}\right),\end{equation} \tag{ 33 }$$

where ${\hat{\boldsymbol{\mu }}}^{{\circ}}=\mathbb{P}{\boldsymbol{\mu }}^{{\circ}}$ is the vector of the standard chemical potentials for the slow species P. It is important to note that first equation (33) establishes the same correspondence between the equilibrium concentrations p_eq of the effective dynamics (19) and the standard chemical potentials ${\hat{\boldsymbol{\mu }}}^{{\circ}}$ as equation (32) does between z_eq and μ^◦ for the elementary reactions. Second, the effective stoichiometric matrix $\hat{\mathbb{S}}$ plays the same role in equation (33) as $\mathbb{S}$ in equation (32). Third, the derivation of equation (33) employs only the topological properties of the CRNs and does not require the time scale separation hypothesis.

4.2. Gibbs free energy of reaction

The Gibbs free energy of the ρ elementary reaction (1), namely the thermodynamic force driving this reaction, is given by

$$\begin{equation}{{\Delta}}_{\rho }G\left(\boldsymbol{z}\right)=\boldsymbol{\mu }\left(\boldsymbol{z}\right)\cdot {\mathbb{S}}_{\rho },\end{equation} \tag{ 34 }$$

which becomes

$$\begin{equation}{{\Delta}}_{\rho }G\left(\boldsymbol{z}\right)=-RT\enspace \mathrm{ln}\enspace \frac{{j}_{+\rho }\left(\boldsymbol{z}\right)}{{j}_{-\rho }\left(\boldsymbol{z}\right)}\end{equation} \tag{ 35 }$$

because of the local detailed balance (32). At equilibrium Δ_ρG_eq ≡ Δ_ρG(z_eq) = 0 since j_+ρ(z_eq) = j_−ρ(z_eq). This means that μ_eq belongs to the cokernel of $\mathbb{S}$.

We define now the Gibbs free energy of each effective reaction Δ_ɛG(p). Collecting the Gibbs free energies of the elementary reactions in the vector ${{\Delta}}_{r}\boldsymbol{G}\equiv {\left(\dots ,{{\Delta}}_{\rho }G,\dots \right)}^{{\intercal}}$, Δ_ɛG is given by

$$\begin{equation}{{\Delta}}_{\varepsilon }G\left(\boldsymbol{p}\right)\equiv {{\Delta}}_{r}\boldsymbol{G}\left(\boldsymbol{z}\right)\cdot {\boldsymbol{c}}_{\varepsilon }=\hat{\boldsymbol{\mu }}\left(\boldsymbol{p}\right)\cdot {\hat{\mathbb{S}}}_{\varepsilon }.\end{equation} \tag{ 36 }$$

We used equation (34), ${\mathbb{S}}^{Q}{\boldsymbol{c}}_{\varepsilon }=0$ and introduced the vector of the chemical potential of the slow species

$$\begin{equation}\hat{\boldsymbol{\mu }}\left(\boldsymbol{p}\right)={\hat{\boldsymbol{\mu }}}^{{\circ}}+RT\enspace \mathrm{ln}\enspace \boldsymbol{p}=\mathbb{P}\boldsymbol{\mu }\left(\boldsymbol{z}\right).\end{equation} \tag{ 37 }$$

The expression of Δ_ɛG in equation (36) is formally the same as Δ_ρG in equation (34). At equilibrium Δ_ɛG_eq ≡ Δ_ɛG(p_eq) = 0, because of the local detailed balance condition (33). This means that ${\hat{\boldsymbol{\mu }}}_{\text{eq}}$ belongs to cokernel of $\hat{\mathbb{S}}$. Unlike Δ_ρG(z), there is no analytical correspondence between Δ_ɛG(p) and ψ^±ɛ(p) and, in general,

$$\begin{equation}{{\Delta}}_{\varepsilon }G\left(\boldsymbol{p}\right)\ne -RT\enspace \mathrm{ln}\enspace \frac{{\psi }^{+\varepsilon }\left(\boldsymbol{p}\right)}{{\psi }^{-\varepsilon }\left(\boldsymbol{p}\right)}.\end{equation} \tag{ 38 }$$

Here, ψ^±ɛ(p) are two positive defined currents such that ψ^ɛ(p) = ψ^+ɛ(p) − ψ^−ɛ(p). This breaks the flux-force relation at the coarse grained level as was already stressed in reference [18].

Finally, we note that, as the local detailed balance, Δ_ɛG is defined using only topological properties of the CRNs.

Example. For the CRN (5) and given the splitting between fast and slow species in (20), the Gibbs free energy of the effective reactions in equation (23) is specified as

$$\begin{equation}\begin{aligned}\hfill {{\Delta}}_{\text{ads}}G& ={\mu }_{{\text{S}}_{\text{m}}}-{\mu }_{\text{S}},\hfill \\ \hfill {{\Delta}}_{\text{enz}}G& ={\mu }_{\text{P}}-2{\mu }_{{\text{S}}_{\text{m}}}.\hfill \end{aligned}\end{equation} \tag{ 39 }$$

4.3. Entropy production rate

The entropy production rate of the elementary dynamics reads

$$\begin{equation}T\dot {{\Sigma}}\left(t\right)=-{{\Delta}}_{r}\boldsymbol{G}\left(\boldsymbol{z}\left(t\right)\right)\cdot \boldsymbol{j}\left(\boldsymbol{z}\left(t\right)\right)=-{{\Delta}}_{\rho }G\left(\boldsymbol{z}\left(t\right)\right){j}^{\rho }\left(\boldsymbol{z}\left(t\right)\right){\geqslant}0.\end{equation} \tag{ 40 }$$

It quantifies the dissipation of the relaxation toward equilibrium. Because of equation (35), the dissipation of each elementary reaction ρ is also non-negative, −Δ_ρG(z(t))j_ρ(z(t)) ⩾ 0 (without summation over ρ).

We define the entropy production rate of the effective dynamics $\hat{\dot {{\Sigma}}}$ using the expression in equation (40), but evaluated along the effective trajectory $\left(\hat{\boldsymbol{q}}\left(\boldsymbol{p}\left(t\right)\right),\boldsymbol{p}\left(t\right)\right)$. Using equations (16), (17) and (36), we obtain that

$$\begin{equation}T\hat{\dot {{\Sigma}}}\left(t\right)=-{{\Delta}}_{\varepsilon }G\left(\boldsymbol{p}\left(t\right)\right){\psi }^{\varepsilon }\left(\boldsymbol{p}\left(t\right)\right)=-{{\Delta}}_{c}\boldsymbol{G}\left(\boldsymbol{p}\left(t\right)\right)\cdot \boldsymbol{\psi }\left(\boldsymbol{p}\left(t\right)\right){\geqslant}0,\end{equation} \tag{ 41 }$$

where we collected the Gibbs free energies of the effective reactions in the vector ${{\Delta}}_{c}\boldsymbol{G}\equiv {\left(\dots ,{{\Delta}}_{\varepsilon }G,\dots \right)}^{{\intercal}}$. The effective entropy production rate $\hat{\dot {{\Sigma}}}\left(t\right)$ is still non-negative by definition: using equation (35) one proves that $\dot {{\Sigma}}\left(t\right){\geqslant}0$ for every j(z) including $\hat{\boldsymbol{j}}\left(\boldsymbol{p}\right)$. However, unlike the elementary dynamics, it is not granted that −Δ_ɛG(p(t))ψ_ɛ(p(t)) ⩾ 0 for every effective reaction ɛ due to the lack of a flux-force relation (38) at the coarse grained level.

If the time scale separation hypothesis holds, namely $\boldsymbol{j}\left(\boldsymbol{z}\right)=\hat{\boldsymbol{j}}\left(\boldsymbol{p}\right)$, the two entropy production rates coincide:

$$\begin{equation}\hat{\dot {{\Sigma}}}\left(t\right)=\dot {{\Sigma}}\left(t\right).\end{equation} \tag{ 42 }$$

When it does not hold, no bound constrains one entropy production rate respect to the other and $\hat{\dot {{\Sigma}}}$ can be lower or greater than $\dot {{\Sigma}}$. This was also observed in reference [26] for driven systems.

Example. For the CRN (5) with the fast and slow species splitting (20), the entropy production rate of the elementary and of the effective dynamics are plotted in figure 2. While $\hat{\dot {{\Sigma}}}\left(t\right){< }\dot {{\Sigma}}\left(t\right)$ for t < 0.5, $\hat{\dot {{\Sigma}}}\left(t\right){ >}\dot {{\Sigma}}\left(t\right)$ for t > 0.5.

Zoom In Zoom Out Reset image size

4.4. Gibbs free energy

The Gibbs free energy of ideal dilute solution is given by

$$\begin{equation}G\left(\boldsymbol{z}\right)=\boldsymbol{\mu }\left(\boldsymbol{z}\right)\cdot \boldsymbol{z}-RT{\Vert}\boldsymbol{z}{\Vert},\end{equation} \tag{ 43 }$$

with ∥a∥   ≡ ∑_ia_i. It is the proper thermodynamic potential of the elementary dynamics of closed CRNs since it has been proven in reference [15] that i) ${\mathrm{d}}_{t}G=-T\dot {{\Sigma}}{\leqslant}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 effective dynamics:

$$\begin{equation}\hat{G}\left(\boldsymbol{p}\right)=\hat{\boldsymbol{\mu }}\left(\boldsymbol{p}\right)\cdot \boldsymbol{p}-RT{\Vert}\boldsymbol{p}{\Vert}.\end{equation} \tag{ 44 }$$

First, we notice that in general Ĝ(p) ≠ G(z). This means that Ĝ does not estimate the exact free energy of CRNs. However, the time derivative of Ĝ according to the effective dynamics (19) reads

$$\begin{equation}\begin{aligned}\hfill {\mathrm{d}}_{t}\hat{G}\left(\boldsymbol{p}\left(t\right)\right)& =\hat{\boldsymbol{\mu }}\left(\boldsymbol{p}\left(t\right)\right)\cdot \hat{\mathbb{S}}\boldsymbol{\psi }\left(\boldsymbol{p}\left(t\right)\right)={{\Delta}}_{c}\boldsymbol{G}\left(\boldsymbol{p}\left(t\right)\right)\cdot \boldsymbol{\psi }\left(\boldsymbol{p}\left(t\right)\right)\hfill \\ \hfill & =-T\hat{\dot {{\Sigma}}}\left(t\right){\leqslant}0\hfill \end{aligned}\end{equation} \tag{ 45 }$$

where we used the definition of the Gibbs free energy of the effective reactions given in equation (36) and the entropy production rate (41). The variation of free energy in a time interval [0, t] thus satisfies

$$\begin{equation}{\Delta}\hat{G}\left(t\right)=-T{\int }_{0}^{t}\mathrm{d}t\enspace \hat{\dot {{\Sigma}}}\left(t\right)=-T{\int }_{0}^{t}\mathrm{d}t\enspace \dot {{\Sigma}}\left(t\right)={\Delta}G\left(t\right),\end{equation} \tag{ 46 }$$

as long as the time scale separation hypothesis holds and, hence, $\hat{\dot {{\Sigma}}}\left(t\right)=\dot {{\Sigma}}\left(t\right)$.

Second, one can prove that Ĝ(p(t)) ⩾ Ĝ_eq ≡ Ĝ(p_eq). Indeed, consider

$$\begin{equation}{\hat{G}}_{\text{eq}}={\hat{\boldsymbol{\mu }}}_{\text{eq}}\cdot {\boldsymbol{p}}_{\text{eq}}-RT{\Vert}{\boldsymbol{p}}_{\text{eq}}{\Vert}.\end{equation} \tag{ 47 }$$

As mentioned in subsection 4.2, ${\hat{\boldsymbol{\mu }}}_{\text{eq}}$ belongs to the cokernel of $\hat{\mathbb{S}}$ and, therefore, it can be expressed as a linear combination of the conservation laws ${\hat{\boldsymbol{l}}}^{\zeta }$:

$$\begin{equation}{\hat{\boldsymbol{\mu }}}_{\text{eq}}={f}_{\zeta }{\hat{\boldsymbol{l}}}^{\zeta }.\end{equation} \tag{ 48 }$$

Thus, the term ${\hat{\boldsymbol{\mu }}}_{\text{eq}}\cdot {\boldsymbol{p}}_{\text{eq}}$ in equation (47) satisfies

$$\begin{equation}{\hat{\boldsymbol{\mu }}}_{\text{eq}}\cdot {\boldsymbol{p}}_{\text{eq}}={f}_{\zeta }\underset{={\hat{L}}^{\zeta }}{\underbrace{{\hat{\boldsymbol{l}}}^{\zeta }\cdot {\boldsymbol{p}}_{\text{eq}}}}={f}_{\zeta }{\hat{\boldsymbol{l}}}^{\zeta }\cdot \boldsymbol{p}\left(t\right)={\hat{\boldsymbol{\mu }}}_{\text{eq}}\cdot \boldsymbol{p}\left(t\right),\end{equation} \tag{ 49 }$$

where we used ${\mathrm{d}}_{t}{\hat{L}}^{\zeta }=0$. This allows us to write the difference between Ĝ(p(t)) and Ĝ_eq as a relative entropy

$$\begin{equation}\hat{G}\left(\boldsymbol{p}\left(t\right)\right)-{\hat{G}}_{\text{eq}}=RT\mathcal{L}\left(\boldsymbol{p}\left(t\right){\Vert}{\boldsymbol{p}}_{\text{eq}}\right),\end{equation} \tag{ 50 }$$

with

$$\begin{equation}\mathcal{L}\left(\boldsymbol{a}{\Vert}\boldsymbol{b}\right)=\sum _{i}{a}_{i}\enspace \mathrm{ln}\left(\frac{{a}_{i}}{{b}_{i}}\right)-\left({a}_{i}-{b}_{i}\right){\geqslant}0,\end{equation} \tag{ 51 }$$

proving that Ĝ(p(t)) ⩾ Ĝ_eq.

In summary, Ĝ in equation (44) is the proper thermodynamic potential of the effective dynamics because of equation (45) and equation (50). It also correctly quantifies the variation of the free energy of the whole CRN because of equation (46).

In open CRNs some species are exchanged with external reservoirs called chemostats. We discuss in parallel the elementary dynamics and the effective one. We consider that the concentrations of chemostatted species are either fixed or slowly driven by the chemostats. These species are therefore treated as slow-evolving species.

We thus split P into two disjoint subsets: the internal species X and the chemostatted species Y. The former (as well as the fast species) evolve in time only because of the chemical reactions and their rate equation is unchanged. The latter evolve in time because of the chemical reactions and of matter flows with the chemostats. These are accounted in the rate equation by introducing the exchange current I(t). For the elementary dynamics of the chemostatted species, the rate equation is given by

$$\begin{equation}{\mathrm{d}}_{t}\boldsymbol{y}\left(t\right)={\mathbb{S}}^{Y}\boldsymbol{j}\left(\boldsymbol{z}\left(t\right)\right)+\boldsymbol{I}\left(t\right),\end{equation} \tag{ 52 }$$

while for the effective dynamics it reads

$$\begin{equation}{\mathrm{d}}_{t}\boldsymbol{y}\left(t\right)={\hat{\mathbb{S}}}^{Y}\boldsymbol{\psi }\left(\boldsymbol{p}\left(t\right)\right)+\boldsymbol{I}\left(t\right).\end{equation} \tag{ 53 }$$

Here we applied the splitting P = (X, Y) to the substoichiometric matrix ${\mathbb{S}}^{P}$ and the effective stoichiometric matrix $\hat{\mathbb{S}}$,

$$\begin{equation}{\mathbb{S}}^{P}=\left(\begin{matrix}\hfill {\mathbb{S}}^{X}\hfill \\ \hfill {\mathbb{S}}^{Y}\hfill \end{matrix}\right),\quad \hat{\mathbb{S}}=\left(\begin{matrix}\hfill {\hat{\mathbb{S}}}^{X}\hfill \\ \hfill {\hat{\mathbb{S}}}^{Y}\hfill \end{matrix}\right),\end{equation} \tag{ 54 }$$

as well as the vector p = (x, y). Note that equations (52) and (53) are merely definitions for the exchange current I(t). At the level of the effective dynamics, a slow variation of y(t) corresponds to a slow variation of I(t). Operationally, one may therefore equally well control I(t) and determine y(t) through equation (53). Both types of external control, on y(t) or I(t), can thus be considered.

Example. For the CRN (5) and given the splitting between fast and slow species in (20), we chemostat the free substrate S and the product P. The adsorbed substrate S_m is now the only internal species. Using the effective stoichiometric matrix in equation (22), we obtain the following two substoichiometric matrix

In figure 3, we compare the elementary dynamics of the CRN (5) for the internal species and the effective one.

Zoom In Zoom Out Reset image size

5.1. 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.

The unbroken conservation laws are those with null entries for the chemostatted species. For the elementary dynamics, the conservation laws {ℓ^ξ} of the fast species (see subsection 3.1) are unbroken by definition. Then, we represent with the vectors $\left\{{\boldsymbol{\ell }}^{{\zeta }_{u}}\right\}\subseteq \left\{{\boldsymbol{\ell }}^{\zeta }\right\}$ the other unbroken conservation laws (if any). Indeed, the quantities ${L}^{{\zeta }_{u}}\left(\boldsymbol{z}\left(t\right)\right)\equiv {\boldsymbol{\ell }}^{{\zeta }_{u}}\cdot \boldsymbol{z}\left(t\right)$, as well as L^ξ(z(t)) ≡ ℓ^ξ ⋅ z(t), are conserved even if the CRN is open:

$$\begin{equation}{\mathrm{d}}_{t}{L}^{{\zeta }_{u}}\left(\boldsymbol{z}\left(t\right)\right)=\underset{=0}{\underbrace{{\boldsymbol{\ell }}^{{\zeta }_{u}}\cdot \mathbb{S}\boldsymbol{j}\left(\boldsymbol{z}\left(t\right)\right)}}+\sum _{\alpha \in Y}\underset{=0}{\underbrace{{\ell }_{\alpha }^{{\zeta }_{u}}}}{I}^{\alpha }\left(t\right)=0.\end{equation} \tag{ 56 }$$

The broken conservation laws are all the other conservation laws, $\left\{{\boldsymbol{\ell }}^{{\zeta }_{b}}\right\}=\left\{{\boldsymbol{\ell }}^{\zeta }\right\}{\backslash}\left\{{\boldsymbol{\ell }}^{{\zeta }_{u}}\right\}$. The corresponding quantities ${L}^{{\zeta }_{b}}\left(\boldsymbol{z}\left(t\right)\right)={\boldsymbol{\ell }}^{{\zeta }_{b}}\cdot \boldsymbol{z}\left(t\right)$ are not conserved:

$$\begin{equation}{\mathrm{d}}_{t}{L}^{{\zeta }_{b}}\left(\boldsymbol{z}\left(t\right)\right)=\underset{=0}{\underbrace{{\boldsymbol{\ell }}^{{\zeta }_{b}}\cdot \mathbb{S}\boldsymbol{j}\left(\boldsymbol{z}\left(t\right)\right)}}+\sum _{\alpha \in Y}\underset{\ne 0}{\underbrace{{\ell }_{\alpha }^{{\zeta }_{b}}}}{I}^{\alpha }\left(t\right)\ne 0.\end{equation} \tag{ 57 }$$

Because of the correspondence in equation (28), the unbroken/broken conservation laws of the effective dynamics are given, respectively, by

$$\begin{equation}{\hat{\boldsymbol{\ell }}}^{{\zeta }_{u}}=\mathbb{P}{\boldsymbol{\ell }}^{{\zeta }_{u}}\quad \text{and}\quad {\hat{\boldsymbol{\ell }}}^{{\zeta }_{b}}=\mathbb{P}{\boldsymbol{\ell }}^{{\zeta }_{b}}.\end{equation} \tag{ 58 }$$

Chemostatting a species does not always break a conservation law [15–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. This allows us to introduce the so-called moieties. They represent parts of (or entire) molecules which are exchanged with the environment through the chemostats. For the elementary dynamics, their concentration vector is specified as

$$\begin{equation}\boldsymbol{m}\left(\boldsymbol{z}\left(t\right)\right)\equiv {\mathbb{M}}^{-1}{\boldsymbol{L}}_{\text{br}}\left(\boldsymbol{z}\left(t\right)\right),\end{equation} \tag{ 59 }$$

while for the effective dynamics it is given by

$$\begin{equation}\hat{\boldsymbol{m}}\left(\boldsymbol{p}\left(t\right)\right)\equiv {\mathbb{M}}^{-1}{\hat{\boldsymbol{L}}}_{\text{br}}\left(\boldsymbol{p}\left(t\right)\right).\end{equation} \tag{ 60 }$$

We introduced the vectors of broken conserved quantities ${\boldsymbol{L}}_{\text{br}}\left(\boldsymbol{z}\left(t\right)\right)={\left(\dots ,{\boldsymbol{\ell }}^{{\zeta }_{b}}\cdot \boldsymbol{z}\left(t\right),\dots \right)}^{{\intercal}}$ and ${\hat{\boldsymbol{L}}}_{\text{br}}\left(\boldsymbol{p}\left(t\right)\right)={\left(\dots ,{\hat{\boldsymbol{\ell }}}^{{\zeta }_{b}}\cdot \boldsymbol{p}\left(t\right),\dots \right)}^{{\intercal}}$, and the matrix $\mathbb{M}$ with entries ${\left\{{\ell }_{\alpha }^{{\zeta }_{b}}\right\}}_{\alpha \in {Y}_{p}}$ (see reference [15–17] for details). The matrix $\mathbb{M}$ is square and nonsingular, and it can be inverted obtaining ${\mathbb{M}}^{-1}$. Comparing equations (59) and (60), one notices that $\boldsymbol{m}\left(\boldsymbol{z}\left(t\right)\right)=\hat{\boldsymbol{m}}\left(\boldsymbol{p}\left(t\right)\right)$ if ${L}^{{\zeta }_{b}}\left(\boldsymbol{z}\left(t\right)\right)={\hat{L}}^{{\zeta }_{b}}\left(\boldsymbol{p}\left(t\right)\right)$ for every broken conservation law, i.e., the slow-evolving species are much more abundant than the fast-evolving species. We stress finally that the number of moieties is equal to the number of Y_p species.

Example. For the CRN (5) and given the splitting between fast and slow species in (20), the conservation law ${\hat{\boldsymbol{\ell }}}^{\text{S}}$ in equation (29) of the effective dynamics is broken when the free product P is chemostatted. If we now chemostat the free substrate S, no conservation law is broken. According to this sequence of chemostatting, the species P belongs to the set Y_p and the matrix $\mathbb{M}$ is simply the scalar

$$\begin{equation}\mathbb{M}=2.\end{equation} \tag{ 61 }$$

The corresponding moiety,

$$\begin{equation}m=\left[\text{P}\right]+\frac{\left[\text{S}\right]}{2}+\frac{\left[{\text{S}}_{\text{m}}\right]}{2},\end{equation} \tag{ 62 }$$

represents the total concentration of product exchanged between the two chemostats.

We consider now the thermodynamic description of open CRNs. The semigrand Gibbs free energy

$$\begin{equation}\mathcal{G}\left(\boldsymbol{z}\right)=G\left(\boldsymbol{z}\right)-{\boldsymbol{\mu }}_{{Y}_{p}}\left(t\right)\cdot \boldsymbol{m}\left(\boldsymbol{z}\right),\end{equation} \tag{ 63 }$$

represents the proper thermodynamic potential for the elementary dynamics [15, 16] since i) ${\mathrm{d}}_{t}\mathcal{G}=-T\dot {{\Sigma}}{\leqslant}0$ for autonomous detailed balanced systems and ii) $\mathcal{G}{\geqslant}{\mathcal{G}}_{\text{eq}}=\mathcal{G}\left({\boldsymbol{z}}_{\text{eq}}\right)$. In analogy to equilibrium thermodynamics, $\mathcal{G}$ is defined from the Gibbs free energy (43) by eliminating the energetic contributions of the matter exchanged with the reservoirs. The latter amounts to concentration of the moieties m(z) of equation (59), times the reference values of their chemical potentials ${\boldsymbol{\mu }}_{{Y}_{p}}\left(t\right)$ which is the vector collecting the values of chemical potentials fixed by the chemostats Y_p. Because of the elementary dynamics (2) and taking into account the exchange current I(t) through equation (52), the evolution of $\mathcal{G}\left(\boldsymbol{z}\right)$ is given by

$$\begin{equation}{\mathrm{d}}_{t}\mathcal{G}\left(\boldsymbol{z}\left(t\right)\right)=-T\dot {{\Sigma}}\left(t\right)+{\dot {w}}_{\text{driv}}\left(t\right)+{\dot {w}}_{\text{nc}}\left(t\right).\end{equation} \tag{ 64 }$$

The entropy production rate $\dot {{\Sigma}}$ is specified in equation (40). The driving work rate _driv accounts for the time dependent manipulation of the chemical potential of the Y_p chemostats,

$$\begin{equation}{\dot {w}}_{\text{driv}}\left(t\right)=-\left({\mathrm{d}}_{t}{\boldsymbol{\mu }}_{{Y}_{p}}\left(t\right)\right)\cdot \boldsymbol{m}\left(\boldsymbol{z}\left(t\right)\right).\end{equation} \tag{ 65 }$$

The nonconservative work rate _nc quantifies the energetic cost of sustaining fluxes of chemical species among the chemostats,

$$\begin{equation}{\dot {w}}_{\text{nc}}\left(t\right)=\boldsymbol{\mathcal{F}}\left(t\right)\cdot \boldsymbol{I}\left(t\right),\end{equation} \tag{ 66 }$$

by means of the force $\boldsymbol{\mathcal{F}}\left(t\right)={\left(\dots ,{\mu }_{\alpha }\left(t\right)-{\left({\boldsymbol{\mu }}_{{Y}_{p}}\left(t\right)\cdot {\mathbb{M}}^{-1}\right)}_{{\zeta }_{b}}{l}_{\alpha }^{{\zeta }_{b}},\dots \right)}_{\alpha \in Y}^{{\intercal}}$. In other words, this is the force keeping the system out of equilibrium (see also appendix A).

For the effective dynamics, we introduce the following semigrand Gibbs free energy:

$$\begin{equation}\hat{\mathcal{G}}\left(\boldsymbol{p}\right)=\hat{G}\left(\boldsymbol{p}\right)-{\boldsymbol{\mu }}_{{Y}_{p}}\left(t\right)\cdot \hat{\boldsymbol{m}}\left(\boldsymbol{p}\right),\end{equation} \tag{ 67 }$$

with $\hat{\boldsymbol{m}}\left(\boldsymbol{p}\right)$ given in equation (60). Note that ${\boldsymbol{\mu }}_{{Y}_{p}}\left(t\right)$ in equation (63) and in equation (67) are exactly the same since they are imposed by the same chemostats. We now show that $\hat{\mathcal{G}}\left(\boldsymbol{p}\right)$ is the proper thermodynamic potential of the effective dynamics of open CRNs.

In general $\hat{\mathcal{G}}\left(\boldsymbol{p}\right)\ne \mathcal{G}\left(\boldsymbol{z}\right)$ since Ĝ(p) ≠ G(z). However, we will now see that their variation in time can be very similar. According to the effective dynamics (19) and taking into account the exchange current I(t) through equation (53), the evolution of $\hat{\mathcal{G}}\left(\boldsymbol{p}\right)$ is given by

$$\begin{equation}{\mathrm{d}}_{t}\hat{\mathcal{G}}\left(\boldsymbol{p}\left(t\right)\right)=-T\hat{\dot {{\Sigma}}}\left(t\right)+{\hat{\dot {w}}}_{\text{driv}}\left(t\right)+{\hat{\dot {w}}}_{\text{nc}}\left(t\right).\end{equation} \tag{ 68 }$$

The entropy production rate $\hat{\dot {{\Sigma}}}$ of equation (41) satisfies $\hat{\dot {{\Sigma}}}=\dot {{\Sigma}}$ as long as the time scale separation hypothesis holds (see subsection 4.3). The driving work rate at the effective level,

$$\begin{equation}{\hat{\dot {w}}}_{\text{driv}}\left(t\right)=-\left({\mathrm{d}}_{t}{\boldsymbol{\mu }}_{{Y}_{p}}\left(t\right)\right)\cdot \hat{\boldsymbol{m}}\left(\boldsymbol{p}\left(t\right)\right),\end{equation} \tag{ 69 }$$

corresponds to _driv of the elementary dynamics when $\hat{\boldsymbol{m}}\left(\boldsymbol{p}\left(t\right)\right)=\boldsymbol{m}\left(\boldsymbol{z}\left(t\right)\right)$. This occurs when ${\hat{L}}^{{\zeta }_{b}}={L}^{{\zeta }_{b}}$ (see subsection 5.1). The nonconservative work rate ${\hat{\dot {w}}}_{\text{nc}}$ has the same expression as the one for the elementary dynamics specified in equation (66). However, the numerical values of ${\hat{\dot {w}}}_{\text{nc}}$ and _nc can be different because different currents I(t) might be necessary to set the same concentrations of the chemostatted species for the elementary and effective 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] satisfies

$$\begin{equation}{\Delta}\hat{\mathcal{G}}={\int }_{0}^{t}\mathrm{d}t\enspace {\mathrm{d}}_{t}\hat{\mathcal{G}}={\int }_{0}^{t}\mathrm{d}t\enspace {\mathrm{d}}_{t}\mathcal{G}={\Delta}\mathcal{G}\end{equation} \tag{ 70 }$$

as long as the time scale separation hypothesis holds (granting also ${\hat{L}}^{{\zeta }_{b}}={L}^{{\zeta }_{b}}$).

If the open system is autonomous (no driving is performed, ${\hat{\dot {w}}}_{\text{driv}}=0$) and detailed balance (all the chemostatted species break a conservation law, ${\hat{\dot {w}}}_{\text{nc}}=0$), then ${\mathrm{d}}_{t}\hat{\mathcal{G}}=-T\hat{\dot {{\Sigma}}}{\leqslant}0$.

Finally, we show that $\hat{\mathcal{G}}\left(\boldsymbol{p}\left(t\right)\right){\geqslant}{\hat{\mathcal{G}}}_{\text{eq}}\equiv \hat{\mathcal{G}}\left({\boldsymbol{p}}_{\text{eq}}\right)$, where the equilibrium state p_eq is set by the Y_p chemostats (see appendix A). To this aim, consider that ${\hat{\boldsymbol{\mu }}}_{\text{eq}}$ belongs to the cokernel of $\hat{\mathbb{S}}$ and can thus be expressed as a linear combination of the conservation laws:

$$\begin{equation}{\hat{\boldsymbol{\mu }}}_{\text{eq}}={f}_{{\zeta }_{u}}{\hat{\boldsymbol{\ell }}}^{{\zeta }_{u}}+{f}_{{\zeta }_{b}}{\hat{\boldsymbol{\ell }}}^{{\zeta }_{b}},\end{equation} \tag{ 71 }$$

where now we distinguish between unbroken and broken conservation laws. Thus, ${\hat{\boldsymbol{\mu }}}_{\text{eq}}\cdot {\boldsymbol{p}}_{\text{eq}}$ satisfies now

$$\begin{equation}\begin{aligned}\hfill {\hat{\boldsymbol{\mu }}}_{\text{eq}}\cdot {\boldsymbol{p}}_{\text{eq}}& ={f}_{{\zeta }_{u}}\underset{={\hat{L}}^{{\zeta }_{u}}}{\underbrace{{\hat{\boldsymbol{\ell }}}^{{\zeta }_{u}}\cdot {\boldsymbol{p}}_{\text{eq}}}}+{f}_{{\zeta }_{b}}{\hat{\boldsymbol{\ell }}}^{{\zeta }_{b}}\cdot {\boldsymbol{p}}_{\text{eq}}\hfill \\ \hfill & ={\hat{\boldsymbol{\mu }}}_{\text{eq}}\cdot \boldsymbol{p}\left(t\right)-{f}_{{\zeta }_{b}}{\hat{\boldsymbol{\ell }}}^{{\zeta }_{b}}\cdot \boldsymbol{p}\left(t\right)+{f}_{{\zeta }_{b}}{\hat{\boldsymbol{\ell }}}^{{\zeta }_{b}}\cdot {\boldsymbol{p}}_{\text{eq}}\hfill \\ \hfill & ={\hat{\boldsymbol{\mu }}}_{\text{eq}}\cdot \boldsymbol{p}\left(t\right)-{\boldsymbol{\mu }}_{{Y}_{p}}\cdot \hat{\boldsymbol{m}}\left(\boldsymbol{p}\left(t\right)\right)+{\boldsymbol{\mu }}_{{Y}_{p}}\cdot {\hat{\boldsymbol{m}}}_{\text{eq}}\hfill \end{aligned}\end{equation} \tag{ 72 }$$

where we used ${\mathrm{d}}_{t}{\hat{L}}^{{\zeta }_{u}}=0$ and the last step is discussed in appendix A. This allows us to write ${\hat{\mathcal{G}}}_{\text{eq}}$ as

$$\begin{equation}{\hat{\mathcal{G}}}_{\text{eq}}={\hat{\boldsymbol{\mu }}}_{\text{eq}}\cdot \boldsymbol{p}\left(t\right)-RT{\Vert}{\boldsymbol{p}}_{\text{eq}}{\Vert}-{\boldsymbol{\mu }}_{{Y}_{p}}\cdot \hat{\boldsymbol{m}}\left(\boldsymbol{p}\left(t\right)\right)\end{equation} \tag{ 73 }$$

and, consequently, the difference between $\hat{\mathcal{G}}\left(\boldsymbol{p}\left(t\right)\right)$ and ${\hat{\mathcal{G}}}_{\text{eq}}$ as a relative entropy

$$\begin{equation}\hat{\mathcal{G}}\left(\boldsymbol{p}\left(t\right)\right)-{\hat{\mathcal{G}}}_{\text{eq}}=RT\mathcal{L}\left(\boldsymbol{p}\left(t\right){\Vert}{\boldsymbol{p}}_{\text{eq}}\right){\geqslant}0.\end{equation} \tag{ 74 }$$

This proves that $\hat{\mathcal{G}}\left(\boldsymbol{p}\left(t\right)\right){\geqslant}{\hat{\mathcal{G}}}_{\text{eq}}$.

Example. For the CRN (5) and given the splitting between fast and slow species in (20), we compare the thermodynamic quantities of the elementary and of the effective dynamics in figure 4.

Zoom In Zoom Out Reset image size

Our thermodynamic framework can be applied directly to effective networks even if the full elementary description is not available. Indeed, all the expressions (e.g., equations (36), (41), (44) and (67)) of the effective thermodynamic quantities require the knowledge of only the effective dynamics (19) and the standard chemical potentials ${\hat{\boldsymbol{\mu }}}^{{\circ}}$ of the slow species. This is the key result of our work.

However, our framework should be applied only to effective 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 effective models must exhibit the following three properties to be thermodynamically consistent.

First, in absence of any chemostatted species (i.e., closed CRNs), the effective model must relax to an equilibrium (30). This 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 effective reactions must be reversible and effective current vector ψ(p) has to admit a nontrivial equilibrium steady-state, i.e., ∃p_eq ≠ 0 such that ψ(p_eq) = 0. If the effective model relaxes to a nonequilibrium steady-state without chemostatted 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 effective network as we proved in subsection 3.1. Third, provided that the time scale separation hypothesis holds, the effective entropy production rate (41) must be greater than or equal to zero. The property $\hat{\dot {{\Sigma}}}\left(t\right){\geqslant}0$ is not granted anymore if the cycle current vector ψ is not specified according to the coarse graining procedure discussed in section 3. In this respect, an effective dynamics giving rise to $\hat{\dot {{\Sigma}}}\left(t\right){< }0$ is necessarily thermodynamically inconsistent.

We now consider two examples, one that can and the other that cannot be characterized using our framework.

7.1. Example 1

Consider the following reactions

$$\begin{equation}\text{S}\underset{}{\overset{{\psi }_{1}}{{\rightleftharpoons}}}I\underset{}{\overset{{\psi }_{2}}{{\rightleftharpoons}}}P,\end{equation} \tag{ 75 }$$

satisfying the dynamics

$$\begin{equation}\begin{aligned}\hfill & {\mathrm{d}}_{t}\left[\text{S}\right]\left(t\right)=-{\psi }_{1}\left(\left[\text{S}\right]\left(t\right),\left[\text{I}\right]\left(t\right)\right)\hfill \\ \hfill & {\mathrm{d}}_{t}\left[\text{I}\right]\left(t\right)={\psi }_{1}\left(\left[\text{S}\right]\left(t\right),\left[\text{I}\right]\left(t\right)\right)-{\psi }_{2}\left(\left[\text{I}\right]\left(t\right),\left[\text{P}\right]\left(t\right)\right)\hfill \\ \hfill & {\mathrm{d}}_{t}\left[\text{P}\right]\left(t\right)={\psi }_{2}\left(\left[\text{I}\right]\left(t\right),\left[\text{P}\right]\left(t\right)\right)\hfill \end{aligned}\end{equation} \tag{ 76 }$$

with

$$\begin{equation}\begin{aligned}\hfill & {\psi }_{1}\left(\left[\text{S}\right],\left[\text{I}\right]\right)=\frac{{a}_{1,s}\left[\text{S}\right]-{a}_{1,i}\left[\text{I}\right]}{{b}_{1,s}\left[\text{S}\right]+{b}_{1,i}\left[\text{I}\right]+{b}_{1,0}},\hfill \\ \hfill & {\psi }_{2}\left(\left[\text{I}\right],\left[\text{P}\right]\right)=\frac{{a}_{2,ii}{\left[\text{I}\right]}^{2}+{a}_{2,i}\left[\text{I}\right]+{a}_{2,ip}\left[\text{I}\right]\left[\text{P}\right]+{a}_{2,p}\left[\text{P}\right]+{a}_{2,pp}{\left[\text{P}\right]}^{2}}{{b}_{2,ii}{\left[\text{I}\right]}^{2}+{b}_{2,i}\left[\text{I}\right]+{b}_{2,ip}\left[\text{I}\right]\left[\text{P}\right]+{b}_{2,p}\left[\text{P}\right]+{b}_{2,pp}{\left[\text{P}\right]}^{2}+{b}_{2,0}}.\hfill \end{aligned}\end{equation} \tag{ 77 }$$

Here {a_•,•} and {b_•,•} are parameters of the effective model.

The first step to study the energetics of the mechanism in equation (75) is to write the dynamical systems (76) as in equation (19). We thus introduce the concentration vector and the stoichiometric matrix

as well as the current vector ψ(p) = (ψ₁([S], [I]), ψ₂([I], [P])). We then verify that the dynamical system (76) does not admit cycles: the stoichiometric matrix in equation (78) has no right null eigenvectors. The left null eigenvector of $\hat{\mathbb{S}}$ is the conservation laws ℓ = (1, 1, 1)^⊺ corresponding to the total concentration, L = ℓ ⋅ p = [S] + [I] + [P]. The equilibrium state of equation (76) can be identified by solving the system of equations ψ(p_eq) = 0 and L = ℓ ⋅ p_eq = ℓ ⋅ p(0).

Assuming that the standard chemical potentials ${\mu }_{\text{S}}^{{\circ}}$, ${\mu }_{\text{I}}^{{\circ}}$ and ${\mu }_{\text{P}}^{{\circ}}$ are known, the second step is to apply the expression of the thermodynamic quantities given in sections 4 and 6. For instance, the dissipation can be quantified with the entropy production rate of equation (41) which reads now

$$\begin{equation}T\hat{\dot {{\Sigma}}}=\left({\mu }_{\text{S}}-{\mu }_{\text{I}}\right){\psi }_{1}+\left({\mu }_{\text{I}}-{\mu }_{\text{P}}\right){\psi }_{2}.\end{equation} \tag{ 79 }$$

By solving the rate equation (76) for a specific set of parameters and initial condition, one can compute the entropy production rate (79) which is shown in figure 5. Notice that $\hat{\dot {{\Sigma}}}\left(t\right){\geqslant}0$ for every time step of the dynamics.With the same strategy, namely, applying the effective expressions given in this work, one can compute also other thermodynamic quantities.

Zoom In Zoom Out Reset image size

Consider now the following elementary mechanism

$$\begin{equation}\begin{aligned}\hfill \text{S}+{\text{E}}_{1}\underset{{k}_{-1}}{\overset{{k}_{+1}}{{\rightleftharpoons}}}{\text{E}}_{1}\text{S}\\ \hfill {\text{E}}_{1}\text{S}\underset{{k}_{-2}}{\overset{{k}_{+2}}{{\rightleftharpoons}}}{\text{E}}_{1}+\text{I}\\ \hfill \text{I}+{\text{E}}_{2}\underset{{k}_{-3}}{\overset{{k}_{+3}}{{\rightleftharpoons}}}{\text{E}}_{2}\text{I}\\ \hfill {\text{E}}_{2}\text{I}\underset{{k}_{-4}}{\overset{{k}_{+4}}{{\rightleftharpoons}}}{\text{E}}_{2}+\text{P}\\ \hfill \text{I}+{\text{E}}_{3}\underset{{k}_{-5}}{\overset{{k}_{+5}}{{\rightleftharpoons}}}{\text{E}}_{3}\text{I}\\ \hfill {\text{E}}_{3}\text{I}\underset{{k}_{-6}}{\overset{{k}_{+6}}{{\rightleftharpoons}}}{\text{E}}_{3}+\text{P}\end{aligned}\end{equation} \tag{ 80 }$$

In the limit of fast-evolving enzymes (E₁, E₂ and E₃) and complexes (E₁S, E₂I and E₃I), the effective dynamics for [S], [I] and [P] provided by the coarse graining procedure discussed in section 3 is consistent with the dynamical system (76). We show in figure 5 the entropy production rate (40) of this elementary mechanism for two different set of initial concentrations of the enzymes and the complexes. In the first 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 effective 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. The initial difference 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.

7.2. Example 2

Consider now the model of gene regulation provided in reference [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. Then, the proteins degrade. The corresponding chemical reaction network is

$$\begin{equation}\begin{aligned}\hfill {G}_{\mathit{\text{A}}}\mathit{\text{A}},\mathrm{e}\mathrm{x}\mathrm{p}{\to }\mathit{\text{A}}\\ \hfill {\mathit{\text{G}}}_{\mathit{\text{B}}}\mathit{\text{B}},\mathrm{e}\mathrm{x}\mathrm{p}{\to }\mathit{\text{B}}\\ \hfill \mathit{\text{A}}\mathit{\text{A}},\mathrm{d}\mathrm{e}\mathrm{g}{\to }\varnothing \\ \hfill \mathit{\text{B}}\mathit{\text{B}},\mathrm{d}\mathrm{e}\mathrm{g}{\to }\varnothing \end{aligned}\end{equation} \tag{ 81 }$$

It evolves according to the following dynamical system

$$\begin{equation}\begin{aligned}\hfill & {\mathrm{d}}_{t}\left[\mathit{\text{A}}\right]\left(t\right)={\psi }_{A,\mathrm{e}\mathrm{x}\mathrm{p}}\left(\left[\mathit{\text{A}}\right]\left(t\right),\left[\mathit{\text{B}}\right]\left(t\right)\right)-{\psi }_{A,\mathrm{d}\mathrm{e}\mathrm{g}}\left(\left[\mathit{\text{A}}\right]\left(t\right)\right)\hfill \\ \hfill & {\mathrm{d}}_{t}\left[\mathit{\text{B}}\right]\left(t\right)={\psi }_{B,\mathrm{e}\mathrm{x}\mathrm{p}}\left(\left[\mathit{\text{A}}\right]\left(t\right),\left[\mathit{\text{B}}\right]\left(t\right)\right)-{\psi }_{B,\mathrm{d}\mathrm{e}\mathrm{g}}\left(\left[\mathit{\text{B}}\right]\left(t\right)\right)\hfill \end{aligned}\end{equation} \tag{ 82 }$$

with the currents

$$\begin{equation}\begin{aligned}\hfill & {\psi }_{A,\mathrm{e}\mathrm{x}\mathrm{p}}\left(\left[\mathit{\text{A}}\right],\left[\mathit{\text{B}}\right]\right)={g}_{\mathit{\text{A}}}+\frac{{a}_{1}{\left[\mathit{\text{A}}\right]}^{4}}{{S}^{4}+{\left[\mathit{\text{A}}\right]}^{4}}+\frac{{b}_{1}{S}^{4}}{{S}^{4}+{\left[\mathit{\text{B}}\right]}^{4}},\hfill \\ \hfill & {\psi }_{B,\mathrm{e}\mathrm{x}\mathrm{p}}\left(\left[\mathit{\text{A}}\right],\left[\mathit{\text{B}}\right]\right)={g}_{\mathit{\text{B}}}+\frac{{a}_{2}{\left[\mathit{\text{B}}\right]}^{4}}{{S}^{4}+{\left[\mathit{\text{B}}\right]}^{4}}+\frac{{b}_{2}{S}^{4}}{{S}^{4}+{\left[\mathit{\text{A}}\right]}^{4}},\hfill \\ \hfill & {\psi }_{A,\mathrm{d}\mathrm{e}\mathrm{g}}\left(\left[\mathit{\text{A}}\right]\right)={k}_{\mathit{\text{B}}}\left[\mathit{\text{A}}\right],\hfill \\ \hfill & {\psi }_{B,\mathrm{d}\mathrm{e}\mathrm{g}}\left(\left[\mathit{\text{B}}\right]\right)={k}_{\mathit{\text{B}}}\left[\mathit{\text{B}}\right].\hfill \end{aligned}\end{equation} \tag{ 83 }$$

Here, g_A, g_B, a₁, a₂, S, b₁, b₂, k_A and k_B are parameters of the model. The concentration of the proteins A and B are the only dynamical variables.

This effective model (82) is designed in such a way that it cannot equilibrate. Indeed, the degradation reactions are irreversible and the currents cannot vanish. The model always relaxes toward 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 defined 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].

In this work, we developed a thermodynamic theory for effective (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 effective models result from underlying elementary (mass-action) networks satisfying the time scale separation hypothesis. This 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 left 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 different species into effective mesostates. These mesostates may then satisfy closed evolution equations. For example, during a catalytic process the complex enzyme–substrate is transformed in many different 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.

This research was funded by the European Research Council project NanoThermo (ERC-2015-CoG Agreement No. 681456).

We discuss here the constraints between the chemical potentials at equilibrium. We consider the case of the effective dynamics, but the same exact reasoning applies to the elementary dynamics (see, for example, appendix A in reference [30]). Because of the local detailed balance (33), the vector of the equilibrium chemical potentials ${\hat{\boldsymbol{\mu }}}_{\text{eq}}$ is a left-null eigenvectors of $\hat{\mathbb{S}}$ and, therefore, it can be expressed as a linear combination of the conservation laws ${\hat{\boldsymbol{l}}}^{\zeta }$:

$$\begin{equation}{\hat{\boldsymbol{\mu }}}_{\text{eq}}^{{\intercal}}=\boldsymbol{f}\cdot \hat{\mathbb{L}},\end{equation} \tag{ A.1 }$$

where we rewrite equation (48) by introducing the vector $\boldsymbol{f}={\left(\dots ,{f}_{\zeta },\dots \right)}^{{\intercal}}$ and the matrix $\hat{\mathbb{L}}$ with entries ${\left\{{\ell }_{\alpha }^{\zeta }\right\}}_{\alpha \in P}$.

We examine the constraints imposed on the chemical potentials by equation (A.1) as if the CRN were open. The chemical species are either internal species X, or chemostatted species Y_p breaking the conservation laws or other chemostatted species Y_f. The set of conservation laws splits into unbroken conservation laws $\left\{{\hat{\boldsymbol{\ell }}}^{{\zeta }_{u}}\right\}$ and broken conservation laws $\left\{{\hat{\boldsymbol{\ell }}}^{{\zeta }_{b}}\right\}$. We show that the equilibrium chemical potentials of the Y_p species set the equilibrium chemical potentials of the Y_f species. We start applying the same splitting to $\hat{\mathbb{L}}$

and f = (f_un, f_br). Here, we introduced the vectors ${\boldsymbol{f}}_{\text{un}}={\left(\dots ,{f}_{{\zeta }_{u}},\dots \right)}^{{\intercal}}$ and ${\boldsymbol{f}}_{\text{br}}={\left(\dots ,{f}_{{\zeta }_{b}},\dots \right)}^{{\intercal}}$, and the matrices ${\mathbb{L}}_{X}^{\text{un}}$ with entries ${\left\{{\ell }_{\alpha }^{{\zeta }_{u}}\right\}}_{\alpha \in X}$, ${\mathbb{L}}_{X}^{\text{br}}$ with entries ${\left\{{\ell }_{\alpha }^{{\zeta }_{b}}\right\}}_{\alpha \in X}$ and ${\mathbb{L}}_{{Y}_{f}}^{\text{br}}$ with entries ${\left\{{\ell }_{\alpha }^{{\zeta }_{b}}\right\}}_{\alpha \in {Y}_{f}}$. The matrix $\mathbb{M}$ was already introduced in subsection 5.1. The zero matrices $\mathbb{0}$ collect the entries of the unbroken conservation laws for chemostatted species (which vanish by definition).

Using equation (A.1) and (A.2), we can recognize that the equilibrium chemical potentials of the Y_p and Y_f species are given by

$$\begin{equation}{\hat{\boldsymbol{\mu }}}_{{Y}_{p},\text{eq}}^{{\intercal}}={\boldsymbol{f}}_{\text{br}}\cdot \mathbb{M},\end{equation} \tag{ A.3 }$$

$$\begin{equation}{\hat{\boldsymbol{\mu }}}_{{Y}_{f},\text{eq}}^{{\intercal}}={\boldsymbol{f}}_{\text{br}}\cdot {\mathbb{L}}_{{Y}_{f}}^{\text{br}},\end{equation} \tag{ A.4 }$$

respectively. The matrix $\mathbb{M}$ is square and nonsingular, and it can be inverted. Hence

$$\begin{equation}{\boldsymbol{f}}_{\text{br}}^{{\intercal}}={\hat{\boldsymbol{\mu }}}_{{Y}_{p},\text{eq}}\cdot {\mathbb{M}}^{-1}.\end{equation} \tag{ A.5 }$$

The equilibrium chemical potentials of the Y_f species thus become

$$\begin{equation}{\hat{\boldsymbol{\mu }}}_{{Y}_{f},\text{eq}}^{{\intercal}}={\hat{\boldsymbol{\mu }}}_{{Y}_{p},\text{eq}}\cdot {\mathbb{M}}^{-1}{\mathbb{L}}_{{Y}_{f}}^{\text{br}},\end{equation} \tag{ A.6 }$$

proving that they depend on the equilibrium chemical potentials of the Y_p species.

In open CRNs, the chemostats fix the chemical potentials of the Y species and, so their concentrations. The Y_p chemostats set a reference equilibrium conditions, meaning

$$\begin{equation}{\hat{\boldsymbol{\mu }}}_{{Y}_{p},\text{eq}}={\boldsymbol{\mu }}_{{Y}_{p}},\end{equation} \tag{ A.7 }$$

$$\begin{equation}{\hat{\boldsymbol{\mu }}}_{{Y}_{f},\text{eq}}^{{\intercal}}={\boldsymbol{\mu }}_{{Y}_{p}}\cdot {\mathbb{M}}^{-1}{\mathbb{L}}_{{Y}_{f}}^{\text{br}}.\end{equation} \tag{ A.8 }$$

Here ${\boldsymbol{\mu }}_{{Y}_{p}}$ is the vector of chemical potentials set by the Y_p chemostats. It is the same whether we consider the elementary dynamics or the effective one. The Y_f chemostats impose the chemical potentials ${\boldsymbol{\mu }}_{{Y}_{f}}$ which, in general, do not satisfy equation (A.8). As a consequence, the nonconservative forces

$$\begin{equation}\boldsymbol{\mathcal{F}}={\left(\dots ,{\mu }_{\alpha }-{\left({\boldsymbol{\mu }}_{{Y}_{p}}\cdot {\mathbb{M}}^{-1}\right)}_{{\zeta }_{b}}{l}_{\alpha }^{{\zeta }_{b}},\dots \enspace \right)}_{\alpha \in Y}^{{\intercal}}\end{equation} \tag{ A.9 }$$

of equation (66) are generated. Indeed, equation (A.9) can be rewritten as

$$\begin{equation}\boldsymbol{\mathcal{F}}={\left(\dots ,{\mu }_{\alpha }-{\mu }_{\alpha ,\text{eq}},\dots \enspace \right)}_{\alpha \in Y}^{{\intercal}},\end{equation} \tag{ A.10 }$$

with μ_α,eq given in equation (A.7) and (A.8). Each entries in $\boldsymbol{\mathcal{F}}$ represents the difference between the chemical potential set by a chemostat and the equilibrium one. Therefore, it represents the force keeping the system out of equilibrium.

Note that, thanks of equation (A.7), we can write ${\boldsymbol{f}}_{\text{br}}^{{\intercal}}={\boldsymbol{\mu }}_{{Y}_{p}}\cdot {\mathbb{M}}^{-1}$. We use this property and the definition of the moieties (60) in the last step of equation (72).