Out-of-equilibrium thermodynamics of autocatalytic networks

Armand Despons
Gulliver Laboratory, UMR CNRS 7083, PSL Research University, ESPCI, Paris F-75231, France

The thermodynamic implications of the topological features of autocatalytic networks remain to be addressed in general. The purpose of this work is to clarify how the stoichiometric trait of autocatalytic networks, namely their absence of a conservation law, shapes their non-equilibrium behavior. To do so, we consider an autocatalytic network coupled with external species acting as food/waste materials, necessary to fulfill mass conservation. Then, we show that the production of autocatalytic species is associated with a conservative influx of external species. From this, we derive the thermodynamic potential of an autocatalytic sub-network. The latter can be obtained from the usual semigrand free-energy of an open system by taking into account the work associated with the conservative influx of external species fueling the production of autocatalytic species. In the end, we identify the cost dedicated to sustaining the production of autocatalytic species and its efficiency. It reveals that sustaining steady production of species in an autocatalytic network is possible only if they are coupled with the environment.

§ INTRODUCTION

At the heart of cellular metabolism is its capability to replicate essential biomolecules from pre-existing ones, employing elementary compounds as elemental building materials.
This trait emphasizes the vital role of autocatalysis, namely, systems that produce more of themselves, within living organisms.
It results in highly intricate behavior <cit.>, enabling growth <cit.> and self-reproduction <cit.>.
These attributes are believed to have played a central role for abiogenesis <cit.>.
Since the pioneering work of Kauffman on autocatalytic sets <cit.>, the development of the theory of chemical reaction networks (CRN) <cit.> has yielded a formal theory of autocatalytic network, extending Ostwald's original definition to a network scale.
Recent studies, in particular, have elucidated the topological features of autocatalytic networks, allowing for their classification and detection in large network <cit.>.

Concurrently, thermodynamics of chemical systems has also undergone intense investigation since the foundational works of Gibbs, who introduced the concept of chemical potential <cit.>, and de Donder, who proposed the idea of affinity <cit.>, quantifying the irreversible forces driving chemical systems.
Specifically, in the 1970s, particular attention has been drawn to the stochastic dynamics of chemical systems based on the chemical master equation <cit.>.
The latter has allowed for the derivation of the foundational concepts regarding non-equilibrium thermodynamics of chemical systems <cit.>.
Since then, the emergence and growing popularity of stochastic thermodynamics, combined with chemical reaction networks theory, have led to the development of a consistent theory assessing the non-equilibrium thermodynamics of chemical network <cit.>.
Living systems representing arguably the most important example of out-of-equilibrium systems, this theory was quickly applied to study biochemical networks <cit.>.
However, establishing stronger connections between topology and non-equilibrium behavior is still the subject of active researches <cit.>.
To the author's knowledge, linking the particular topology of autocatalytic networks to their non-equilibrium behavior, has not yet been addressed in the general case.

The present paper aims to provide a thermodynamic framework
to study the cost associated to the production of species of an autocatalytic network.
To do so, we rely on the chemical reaction network formalism applied to autocatalytic networks in the deterministic regime, in particular, our approach is based on the recent framework developed by Avanzini et al. in Ref. <cit.>.
Specifically, we consider an arbitrary chemical network in which an autocatalytic sub-network is embedded.
This approach emphasizes the different role of the species: either they belong to the autocatalytic sub-network or they are external.
These are crucial to fulfill mass-conservation, a necessary condition for any (bio)chemical network.

The paper is organized as follows: we first introduce the general framework used throughout the work.
Most of it relies on the block decomposition of the stoichiometric matrix highlighting the presence of an autocatalytic sub-network, consisting of autocatalytic species and autocatalytic reactions, coupled with external species (possibly involved in some additional reactions).
Then, after recalling the usual properties of CRNs, we show how the absence of conservation law in this sub-network influences can be used to decompose the elementary fluxes on a physically relevant basis that tell apart the production of autocatalytic species and the cycles of the sub-network.
This decomposition demonstrates that the steady production of an autocatalytic species is possible only if this species is extracted from the system by the mean a flux coupling the sub-network to its surroundings.
From this we extract the condition for unconditional detailed balance in the sub-network.
The latter is expressed in a topological condition relating the cycle of the sub-network to those of the full network.
Then, we derive the proper thermodynamic potential for the autocatalytic sub-network and its associated decomposition of the entropy production rate.
These constitute the central results of the present work.
From this we extract the cost to sustain production of autocatalytic species.
Finally, we study the efficiencies associated to the production in the autocatalytic network.
For the sake of simplicity, we address the case of ideal solutions for the majority of the text.
In Appendix A, we relax this assumption by specifying the modifications that need to be taken into account in the case of non-ideal solutions.

§ AUTOCATALYTIC CRN

§.§ Setup

We consider a Chemical Reaction Network (CRN) <cit.> where  species are either autocatalytic species, denoted $z \in \Z$, or external species, $e \in \E$ <cit.>.
These species are typically the food and waste species as well as the external catalysts/enzymes.
The two types of species are involved in autocatalytic reactions $\rho \in \R$ and, in addition, we assume that the external species may be reactant/product of some additional reactions $r \in \nR$.
From this, the stoichiometric matrix $\Nabla$ of such CRN can be written as follows:

\begin{equation}
\label{eq:Nabla}
\boldsymbol{\nabla} = \quad \begin{pNiceArray}{m{1em}m{1em}|m{1em}m{1em}}[last-row, first-col]
\Vdots[line-style={solid, <->}, shorten=2pt]_{\small \rotatebox{90}{$\E$} }  & \Block{2-2}{\NablaER} &  & \Block{2-2}{\NablaEnR} & \\
& & & & \\
\Hline
\Vdotsfor{3}[line-style={solid, <->}, shorten=2pt]_{\small \rotatebox{90}{$\Z$} } \hspace*{1em} & \Block{3-2}{\S} & & \Block{3-2}{ \boldsymbol{0} } & \\[-0.5em]
&  &  &  &   \\
& & & & \\
& \Hdotsfor{2}[line-style={solid, <->}, shorten=1pt]_{ \R } &  \Hdotsfor{2}[line-style={solid, <->}, shorten=1pt]_{ \nR }  \\
\end{pNiceArray},
\vspace*{0.5em}
\end{equation}

where the rectangular matrix $\S \equiv \Nabla^{\Z}_{\R}$ is the restriction of $\Nabla$ on the subsets of autocatalytic species and autocatalytic reactions and $\NablaER$ (resp. $\NablaEnR$) is the restriction of $\Nabla$ on the external species and autocatalytic reactions (resp. external reactions).
On top of that, we assume $\S$ that has no conservation law or, equivalently, $\ker \left[ \S^\top \right] = \varnothing$.
Thanks to this condition, the triplet $(\Z, \, \R, \, \S)$ defines an autocatalytic CRN.
For this reason, in what follows, we will refer to this sub-network as the autocatalytic sub-network.

We call $\j$ (resp. $\v$) the elementary fluxes associated to the autocatalytic reactions $$ (resp. additional reactions $$).
The dynamic of the open CRN follows the standard kinetic rate equations:
\begin{align}
\label{eq:dyn_external}
\d_t [\e] & ~=~ \NablaER \cdot \j ~+~ \NablaEnR \cdot \v  ~+~ \I^{\E}, \\[0.5em]
\label{eq:dyn_autocat}
\d_t [\z] & ~=~ \S \cdot \j ~+~ \I^{\Z},
\end{align}
where $[s]$, for $s ∈S$ (with $S = $ or $$), is the vector of chemical concentrations and $^S$ the external fluxes coupling the CRN with the environment.
Furthermore, we assume that, for all chemical reactions, its backward counterpart also exists.
Then, the stoichiometric matrix and the elementary fluxes are the resultant of a backward and a forward contribution:
\begin{gather}
\label{eq:Nabla_pm}
\Nabla = \Nabla_- - \Nabla_+, \\[0.2em]
\label{eq:j_pm}
\j = \j^+ - \j^- , \\[0.2em]
\label{eq:v_pm}
\v = \v^+ - \v^- ,
\end{gather}
where $_-$ (resp. $_+$) describes the products (resp. reactants) stoichiometry and $j^+$/$v^+$ (resp. $j^-$/$v^-$) are the forward (resp. backward) unidirectional fluxes.
Note that the Eq.~\eqref{eq:Nabla_pm} can be written for all the non-vanishing blocks constituting $$ in Eq.~\eqref{eq:Nabla}.
We also assume that a chemical master equation describes the stochastic trajectories in the concentration space.
In the infinite volume limit, the latter results in the deterministic dynamics Eqs.~\eqref{eq:dyn_external}-\eqref{eq:dyn_autocat} with mass-action law for the unidirectional fluxes \cite{Ge2016}:
\begin{gather}
\label{eq:MAL_def}
j^{\pm \rho} = k^{\pm \rho} \, \prod_{s \in \E \sqcup \Z} [s]^{(\nabla_\pm)^s_\rho}, \\[0.5em]
%	& \quad &
v^{\pm r} = \nu^{\pm r} \, \prod_{e \in \E} [e]^{(\nabla_\pm^\E)^s_r},
\end{gather}
where, $ρ∈$ and $r ∈$, and with $k^±ρ$/$ν^±r$ the kinetic rate constants.
We will assume ideal solution for the majority of the text.
In this case, the chemical potential of species $s ∈⊔$ is
\begin{equation}
\label{eq:mu_ideal}
\mu_s = \mu^\circ_s + RT \ln \, [s].
\end{equation}
In Appendix A, we discuss the non-ideal case, which is necessary to assess the case of interacting species.
This is encountered in biological systems, charged species are ubiquitous and interactions between species can be important.
Unless otherwise stated, the main results obtained below remain valid in the non-ideal case.
The local detailed balance condition relates the standard chemical potentials to the kinetic rate constants, ensuring thermodynamic consistency:
\begin{align}
\label{eq:LDB}
\ln \dfrac{
= - \dfrac{
\boldsymbol{\mu}^\circ \cdot \Nabla_\rho
& \quad &
\ln \dfrac{
\nu^{+r}
\nu^{-r}
= - \dfrac{
\boldsymbol{\mu}^\circ \cdot \Nabla_r
\end{align}
where $_ρ$ (resp. $_r$)  corresponds to the column associated to the autocatalytic reaction $ρ∈$ (resp. additional reaction $r ∈$) in the full stoichiometric matrix.
We now turn to application of this setup on the following network that we will keep all along to illustrate the main results.

\NiceMatrixOptions{
code-for-last-row = \color{black!70}\scriptstyle,
code-for-first-col = \color{black!70}\scriptstyle,
\paragraph{\textbf{Example 1.}} We consider the following CRN:
\begin{equation}
\label{eq:CRN_example}
\begin{gathered}
\ce{\sf{F}_0 + \sf{Z}_0 <=>[\sf{1}][] $\sf{Z}_0^{\, \ast}$ <=>[\sf{2}][]  \sf{Z}_0 + \sf{Z}_1 + \sf{W}_0 } \\
\ce{\sf{F}_1 + \sf{Z}_1 <=>[\sf{3}][] $\sf{Z}_1^{\, \ast}$ <=>[\sf{4}][] \sf{Z}_1 + \sf{Z}_0 + \sf{W}_1 } \\[0.5em]
\ce{$\sf{F}_{\mathfrak{a}}$ <=>[\sf{5}][] \sf{Z}_0 + $\sf{W}_{\mathfrak{a}}$ } \\
\ce{$\sf{F}_{1-\mathfrak{a}}$ <=>[\sf{6}][] \sf{Z}_1 + $\sf{W}_{1-\mathfrak{a}}$ } \\[0.5em]
\ce{\sf{W}_0 <=>[\sf{7}][] \sf{F}_0} \\
\ce{\sf{W}_1 <=>[\sf{8}][] \sf{F}_1}
\end{gathered}
\end{equation}
consisting for its four first reactions of the Hinshelwood cycle with two species ($Z_i$) and with intermediates ($Z_i^  ∗$), coupled to fuel ($F_i$) and waste ($W_i$) species.
Moreover, aside from the production of $Z_i$ catalyzed by $Z_1-i$, we also consider its non-catalytic production in reactions $5$, $6$.
The fuel/waste species in these two reactions depend on the value of the parameter $𝔞$ ($𝔞 = 0,  1$), the reason for that will become clear later.
Finally, a waste species can be converted back to its associated fuel species.
This makes network Eq.~\eqref{eq:CRN_example} one the simplest food-generated autocatalytic sets (RAF) \cite{Hordijk2004, Xavier2020}.
Indeed, the autocatalysts $Z_i$ cross-catalyze each other and the two can be produced directly from the food set thank to reaction $5$, $6$.
Note that, often, catalyzed reactions are written without the intermediate species.
In particular, the Hinshelwood cycle with two species is usually written as \ce{\sf{Z}_i  <=>  \sf{Z}_i + \sf{Z}_{1-i}} for $i = 0, 1$.
This last reaction is not elementary, hence, neither mass action kinetics nor the local detailed balance relation is expected to hold, \emph{a priori}.
Recent studies have been interested in overcoming this limitation by building a consistent thermodynamic framework for coarse-grained chemical networks \cite{Wachtel2018, Avanzini2020}.
These go beyond the scope of the present paper and we will limit our analysis to network of elementary reactions. \\
In the full CRN Eq.~\eqref{eq:CRN_example} the subsets of species reads:
\begin{align}
\E = \left\lbrace \sf{F}_0, \, \sf{F}_1, \, \sf{W}_0, \, \sf{W}_1 \right\rbrace,
& &
\Z = \left\lbrace \sf{Z}_0, \, \sf{Z}_0^{\, \ast}, \, \sf{Z}_1, \, \sf{Z}_1^{\, \ast} \right\rbrace,
\end{align}
while the subsets of reactions are
\begin{align}
\R = \left\lbrace \sf{1}, \, \sf{2}, \, \sf{3}, \, \sf{4}, \, \sf{5}, \, \sf{6} \right\rbrace,
& \quad &
\nR = \left\lbrace \sf{7}, \, \sf{8} \right\rbrace.
\end{align}
Hence, the two stoichiometric sub-matrices associated to the external species are:
\begin{gather}
\label{eq:Nabla_E_example}
\NablaER = \hspace{1em}
\begin{pNiceArray}{cccccc}[first-col, last-row]
\sf{F}_0 & -1 & 0 & 0 & 0 & -\delta_{\mathfrak{a},\,0} & -\delta_{\mathfrak{a},\,1} \\[0.3em]
\sf{F}_1 & 0 & 0 & -1 & 0 & -\delta_{\mathfrak{a},\,1} & -\delta_{\mathfrak{a},\,0} \\[0.3em]
\sf{W}_0 & 0 & 1 & 0 & 0 & \delta_{\mathfrak{a},\,0} & \delta_{\mathfrak{a},\,1} \\[0.3em]
\sf{W}_1 & 0 & 0 & 0 & 1 & \delta_{\mathfrak{a},\,1} & \delta_{\mathfrak{a},\,0} \\[0.3em]
& \sf{1} & \sf{2} & \sf{3} & \sf{4} & \sf{5} & \sf{6}
\end{pNiceArray}, \\[0.8em]
%	& \quad &
\NablaEnR = \hspace{1em}
\begin{pNiceArray}{cc}[first-col, last-row]
\sf{F}_0 & 1 & 0 \\[0.3em]
\sf{F}_1 & 0 & 1  \\[0.3em]
\sf{W}_0 & -1 & 0  \\[0.3em]
\sf{W}_1 & 0 & -1  \\[0.3em]
& \sf{7} & \sf{8}
\end{pNiceArray},
\end{gather}
where $δ_𝔞,   i$ stands for the Kronecker symbol.
Furthermore, the stoichiometric matrix of the autocatalytic sub-network is,
\begin{equation}
\label{eq:S_example}
\S = \hspace{1em}
\begin{pNiceArray}{cccccc}[first-col, last-row]
\sf{Z}_0 		& -1 & 1 & 0 & 1 & 1 & 0 \\[0.3em]
\sf{Z}_0^{\, \ast} 	& 1 & -1 & 0 & 0 & 0 & 0 \\[0.3em]
\sf{Z}_1 		& 0 & 1 & -1 & 1 & 0 & 1 \\[0.3em]
\sf{Z}_1^{\, \ast} 	& 0 & 0 & 1 & -1 & 0 & 0 \\[0.3em]
& \sf{1} & \sf{2} & \sf{3} & \sf{4} & \sf{5} & \sf{6}
\end{pNiceArray}.
\end{equation}
The latter has an empty left nullspace thus $(,   ,   )$ has no conservation law and is the autocatalytic sub-network.

\subsection{Conservation laws}

Conservation laws of the full CRN \eqref{eq:Nabla} are the row vectors $ℓ$ that belong to the left nullspace of the full stoichiometric matrix, $ℓ ·= 0$ \cite{Rao2018, Rao2016, Avanzini2021}.
Considering a linear basis ${ℓ^i }$ of the conservation laws, we construct $Ł= {ℓ^i }$, the matrix whose $i$-th row is the conservation law $ℓ^i$.
Taking into account the block decomposition of full stoichiometric matrix Eq.~\eqref{eq:Nabla} we have:
\begin{align}
\L_\E \cdot \Nabla^\E_\R + \L_\Z \cdot \S = 0, \label{eq:L_R} & \quad &
\L_\E \cdot \Nabla^\E_{\nR} = 0.
\end{align}

Furthermore, because the autocatalytic sub-network has no conservation law, its stoichiometric matrix $$ has an empty left nullspace.
As a consequence, the rank of $$ is equal to the number of autocatalytic species $||$ and, then, it admits a (non-unique) right inverse $𝔾$:
\begin{equation}
\S \cdot \invS = \mathbbm{1}.
\end{equation}
In what follows, the column of $𝔾$ associated to species $z ∈$ is denoted $g_z$.
It defines an \emph{elementary mode of production} of species $z$, which is a pathway of autocatalytic reactions producing one unit of species $z$, and letting the other autocatalytic species unchanged.
By a linear combination of these modes, we can construct any state $∈ℝ^||_+$, this implies that the stoichiometric manifold of the autocatalytic sub-network is spanning the whole space.
It makes the dynamics of the autocatalytic species unconstrained when considered independently from the rest of the network.
An important consequence of this is found considering the sum of all the elementary modes, $𝐠 = ∑_z g_z$.
The latter verifies:
\begin{equation}
\S \cdot \mathbf{g} = \boldsymbol{1},
\end{equation}
where $1$ is a vector whose components are equal to 1.
This pathway defines an overall reaction where all the autocatalytic species are increased by one unit, its existence is a sufficient condition for autocatalysis.
Crucially, along an elementary mode of production, mass conservation is violated when one considers only the autocatalytic species.
This underlines the need of the external species bearing the conservation laws in the full network.
Indeed, multiplying both sides of $(<ref>)$ by $$ connects $Ł_$ to $Ł_$:
\begin{equation}
\L_{\Z} = - \L_{\E} \cdot \Nabla^\E_\R \cdot \invS.
\end{equation}
Hence, imposing a non-vanishing external flux on all the external species breaks all the conservation law of full CRN.
More precisely, there exists (a non unique) set of external species, called the \emph{potential species}, such that each of the potential species is associated to a conservation law of the full CRN.
When coupled with the environment, a potential species breaks its associated conservation law.
Formally, calling $⊂$ the subset of potential species, the sub-matrix $Ł_$ is square and, most importantly, non-singular.
Note that linear algebra ensures the existence of such a subset as, by construction, $Ł$ is full rank.
The remaining external species are the \emph{external force species}, which belong to the subset $⊂$, hence, one has $= ⊔$.
As we will see, the reason why potential species are named as such is that the non-conservative work associated with them necessarily cancels out.
Hence, the potential species are associated to a change in the equilibrium state of the full network and, if no external force species are present, no non-conservative force is there to maintain the full network away from this equilibrium state \cite{Avanzini2021, Avanzini2022}.
Splitting Eq.~\eqref{eq:L_R} on these new subsets yields:
\begin{gather}
\L_{\Ep} \cdot \Nabla^{\Ep}_{\R} +  \L_{\Ef} \cdot \Nabla^{\Ef}_{\R} + \L_{\Z } \cdot \S = 0,
%	\L_{\Ep} \cdot \Nabla^{\Ep}_{\nR} +  \L_{\Ef} \cdot \Nabla^{\Ef}_{\nR} = 0.
\end{gather}
and similarly for the additional reactions $$.
From the non-singularity of $Ł_$, we can relate the stoichiometric matrices associated to the potential species to the other stoichiometric sub-matrices,
\begin{gather}
\Nabla^{\Ep}_{\R} = - \M_{\Ef} \cdot \Nabla^{\Ef}_{\R} - \M_{\Z} \cdot \S, \label{eq:Nabla_Ep_R}
%	\Nabla^{\Ep}_{\nR} = - \M_{\Ef} \cdot \Nabla^{\Ef}_{\nR} \label{eq:Nabla_Ep_nR}.
\end{gather}
where we introduced the \emph{moiety matrix} \cite{Avanzini2021}:
\begin{equation}
\label{eq:M_def}
\M = \left(\L_{\Ep}\right)^{-1} \cdot \L = \left( \mathbbm{1}_{|\Ep|}  ~,~ \M_{\Ef} ~,~ \M_{\Z} \right).
\end{equation}
We will make extensive use of the  decomposition Eq.~\eqref{eq:Nabla_Ep_R} in what comes next.
\vspace{2em}

\paragraph{\textbf{Example 2.}}
The network Eq.~\eqref{eq:CRN_example} has two linearly independent conservation laws, ${ℓ^1, ℓ^2 }$, which are independent of the value taken by the parameter $𝔞$:
\begin{equation}
\label{eq:L_example}
\L = \hspace{1em}
\begin{pNiceArray}{cc|cc|cccc}[first-col, last-row]
\boldsymbol{\ell}^1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 1 \\[0.5em]
\boldsymbol{\ell}^2 & 1 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \\[0.5em]

& \sf{F}_0 & \sf{F}_1 & \sf{W}_0 & \sf{W}_1 & \sf{Z}_0 & \sf{Z}_0^{\, \ast} & \sf{Z}_1& \sf{Z}_1^{\, \ast}
\end{pNiceArray}.
\end{equation}
From this we see that the fuel species $F_1$  can be associated with the conservation law $ℓ^1$ such that, when $F_1$ is externally controlled by the mean of a non vanishing external flux $I^F_0$, $ℓ^1$ is broken.
Similarly, the fuel species $F_0$ can be paired with $ℓ^2$.
\begin{align}
\Ep = \left\lbrace
\sf{F}_0, \, \sf{F}_1
\right\rbrace,
& \quad &
\Ef = \left\lbrace
\sf{W}_0, \, \sf{W}_1
\right\rbrace.
\end{align}
Note that this choice is not unique and, apart from ${F_0, W_0 }$ and ${F_1, W_1 }$, any other choice would have been a valid choice for $$.
According to this choice, the \emph{moiety matrix} reads:
\begin{equation}
\label{eq:M_example}
\M = \hspace{1em}
\begin{pNiceArray}{cc|cc|cccc}[first-col, last-row]
\sf{F}_0 & 1 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \\[0.5em]
\sf{F}_1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 1 \\[0.5em]

& \sf{F}_0 & \sf{F}_1 & \sf{W}_0 & \sf{W}_1 & \sf{Z}_0 & \sf{Z}_0^{\, \ast} & \sf{Z}_1& \sf{Z}_1^{\, \ast}
\end{pNiceArray}.
\end{equation}
This matrix captures the intuition that a fragment, or \emph{moiety}, of $F_0$ (resp. $F_1$) should be present in $Z_0^  ∗$ and $W_0$ (resp. $Z_1^  ∗$ and $W_1$).

\begin{figure}[t!]
\centering
\includegraphics{Schema_system_2.eps}
\caption{\label{fig:schema_system}
Schematic representation of the system Eq.~\eqref{eq:Nabla}.
To fulfill mass-conservation, the autocatalytic sub-network must be coupled with external species.
Moreover, the autocatalytic reactions ($\R$) can be further split into two disjoint subsets.
Reactions in $\R_1$ account for the production of autocatalytic species ($\Z$), it requires a conservative influx of external species, $\Delta \e (\z)$, given by Eq.~\eqref{eq:Delta_e}.
On the other hand, reactions belonging to $\R_2$ are associated to the cycles of $\S$, these are sustained thanks to a non-conservative influx of external species.
Finally, all species can be coupled with the environment by the mean of external fluxes.
\end{figure}

\subsection{Cycles}

\subsubsection{Cycles of the sub-network}

Cycles of the full CRN are the column vectors in the right nullspace of the full stoichiometric matrix: $·γ = 0$.
According to the block decomposition Eq.~\eqref{eq:Nabla} we have:
\begin{gather}
\NablaER \cdot \boldsymbol{\gamma}^\R + \NablaEnR \cdot \boldsymbol{\gamma}^{\nR} = \boldsymbol{0}, \label{eq:c_E} \\[0.5em]
\S \cdot \boldsymbol{\gamma}^\R = \boldsymbol{0}, \label{eq:c_Z}
\end{gather}
where $γ^R$ stands for the restriction of the cycle on the subset of reactions $R$.
Along any cycle $γ$, the Wegscheider condition prescribes \cite{Schuster1989}:
\begin{equation}
\label{eq:Wegscheider}
\prod_{\rho \in \R} \left(  \dfrac{k^{+ \rho}}{k^{-\rho}} \right)^{\gamma^\rho}
\prod_{r \in \nR} \left(  \dfrac{\nu^{+ r}}{\nu^{- r}} \right)^{\gamma^r} = 1.
\end{equation}
Further, from Eq.~\eqref{eq:c_Z} we see that if $γ$ is a cycle then, either $γ^= 0$, or $γ^$ is itself a cycle of the autocatalytic sub-network.
Once \emph{all} the external species are externally controlled, and if no external control is exerted on the autocatalytic species, the cycles of $$ become the \emph{emergent cycles} \cite{Polettini2014}.
From the rank-nullity theorem, there is $|| - ||$ independent cycles in the autocatalytic sub-network.
In what follows, we let ${c_ε}_1 ≤ε≤|| - || $ be a linear basis of these cycles.

We can split the set of autocatalytic reactions into two disjoint subsets, $= _1 ⊔_2$,
such that,
\begin{equation}
\S = \left( \S_{\R_1} ~,~ \S_{\R_2} \right),
\end{equation}
with $__1$ being square and non-singular.
Note that, here, autocatalysis (\emph{i.e.} the absence of conservation law for $$) is essential.
Indeed, from this we know that $$ is full rank and, as a consequence, linear algebra ensures that one only needs to remove reactions (columns of $$) to end up with a square non-singular matrix.
Removing one reaction in $_2$ from the network causes exactly one cycle of $$ to disappear, hence, the network $__1$ has no cycle \cite{Cengio2022}.
Consequently, the number of reactions in $_2$ matches the number of cycles associated to $$, \emph{i.e.}, $|_2| = || - ||$.
Importantly, for cycle $c_ε$, its restriction to $_1$, $c_ε^_1$, is totally determined by $c_ε^_2$ \emph{via}
\begin{equation}
\label{eq:c_R1}
\boldsymbol{c}_\varepsilon^{\R_1} = - \left( \S_{\R_1} \right)^{-1} \cdot \S_{\R_2} \cdot \boldsymbol{c}_\varepsilon^{\R_2}.
\end{equation}
This last relation implies that one can the fluxes along the cycles of $$ are completely characterized by the of of reactions $_2$.

\subsubsection{Topological condition for the cycles of the sub-network}

We now examine the condition upon which the cycles of $$ are also cycles for $$ which turns out to play a crucial role for what follows.
For any cycle of $$ we let
\begin{equation}
\widetilde{\boldsymbol{c}_\varepsilon} = \begin{pmatrix} \boldsymbol{c}_\varepsilon \\[0.5em] \boldsymbol{0} \end{pmatrix},
\end{equation}
be its extension for the full CRN, where the null part spans reactions in $/$.
Plugging this expression into Eq.~\eqref{eq:c_E} yields:
\begin{equation}
\left( \Nabla^\E_{\R_2} - \Nabla^\E_{\R_1} \cdot \left( \S_{\R_1} \right)^{-1} \cdot \S_{\R_2} \right) \cdot \boldsymbol{c}_{\varepsilon}^{\R_2} = \boldsymbol{0},
\end{equation}
in which we used Eq.~\eqref{eq:c_R1} to eliminate $c_ε^_1$.
As a consequence, denoting
\begin{equation}
\label{eq:effective_stoichiometric}
\widetilde{\Nabla}^\E_{\R_2} = \Nabla^\E_{\R_2} - \Nabla^\E_{\R_1} \cdot \left( \S_{\R_1} \right)^{-1} \cdot \S_{\R_2},
\end{equation}
if $^__2 = 0$, all the cycles of $$ are also cycles of the full network.
Conversely, any cycle of $$ that is also a cycle of $$ belongs to the right nullspace of $^__2$. \\
This \emph{effective} stoichiometric matrix plays a central role in the behavior of the autocatalytic sub-network.
Indeed, the presence of cycles in the autocatalytic sub-network prevent it from reaching equilibrium when all the external species are externally controlled, except when $^__2 = 0$.
In that case, the autocatalytic CRN can reach an equilibrium state for arbitrary values of the chemical potential of the external species.

\subsubsection{A unique right-inverse}

To tell apart the contribution of the  cycles (reactions in $_2$) from the the modes of production (columns of $𝔾$), we search for a right-inverse such that $𝔾^_2 = 0$.
This additional constraint removes the non-uniqueness of $𝔾$ as the only right-inverse fulfilling this property is,
\begin{equation}
\label{eq:G_matrix}
\mathbb{G} = \quad
\begin{pNiceArray}{m{1.5em}m{1.5em}}[last-row, first-col]
\Vdots[line-style={solid, <->}, shorten=2pt]_{\small \rotatebox{90}{$\R_1$} }  & \Block{2-2}{(\S_{\R_1})^{-1}} & \\
& & \\
\hdottedline
\Vdotsfor{2}[line-style={solid, <->}, shorten=2pt]_{\small \rotatebox{90}{$\R_2$} } \hspace*{1em} & \Block{2-2}{\boldsymbol{0}}  & \\
&  & \\
& \Hdotsfor{2}[line-style={solid, <->}, shorten=1pt]_{ \Z }  \\
\end{pNiceArray}.
\vspace*{0.5em}
\end{equation}
From this, we see that the elementary modes drawn from the above right inverse have vanishing entries in reactions in the subset $_2$, which accounts for the cycles of $$.
Importantly, the elementary modes of production drawn from Eq.~\eqref{eq:G_matrix} are linearly independent.
Furthermore,  gathering these modes with the basis of the cycles of $$, ${c_ε}$, defines a linear basis.
\vspace{2em}

\paragraph{\textbf{Example 3}.}
The CRN Eq.~\eqref{eq:CRN_example} admits two linearly independent cycles $γ_1$ and $γ_2$.
Both depend on the value of $𝔞$, if $𝔞 = 0$,
\begin{equation}
\label{eq:cycle_example_0}
\begin{gathered}
\boldsymbol{\gamma}_1  = \hspace{0.5em}
\begin{pNiceArray}{w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}}[last-row]
0 & 0 & -1 & -1 & 1 & 0 & 1 & -1 \\
\sf{1} & \sf{2} & \sf{3} & \sf{4} & \sf{5} & \sf{6} & \sf{7} & \sf{8}
\end{pNiceArray}^\top, \\[0.8em]
\boldsymbol{\gamma}_2  = \hspace{0.5em}
\begin{pNiceArray}{w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}}[last-row]
1 & 1 & 1 & 1 & -1 & -1 & 0 & 0 \\
\sf{1} & \sf{2} & \sf{3} & \sf{4} & \sf{5} & \sf{6} & \sf{7} & \sf{8}
\end{pNiceArray}^\top ;
\end{gathered}
\end{equation}
while, if $𝔞 = 1$, one has:
\begin{equation}
\label{eq:cycle_example_1}
\begin{gathered}
\boldsymbol{\gamma}_1  = \hspace{0.5em}
\begin{pNiceArray}{w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}}[last-row]
-1 & -1 & 0 & 0 & 0 & 1 & 0 & 0 \\
\sf{1} & \sf{2} & \sf{3} & \sf{4} & \sf{5} & \sf{6} & \sf{7} & \sf{8}
\end{pNiceArray}^\top, \\[0.8em]
\boldsymbol{\gamma}_2  = \hspace{0.5em}
\begin{pNiceArray}{w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}}[last-row]
0 & 0 & -1 & -1 & 1 & 0 & 0 & 0 \\
\sf{1} & \sf{2} & \sf{3} & \sf{4} & \sf{5} & \sf{6} & \sf{7} & \sf{8}
\end{pNiceArray}^\top.
\end{gathered}
\end{equation}
On the other hand, the autocatalytic sub-network has two cycles:
\begin{equation}
\label{eq:emergent_cycles_example}
\begin{gathered}
\boldsymbol{c}_1 = \hspace{0.5em}
\begin{pNiceArray}{w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}}[last-row]
-1 & -1 & 0 & 0 & 0 & 1   \\
\sf{1} & \sf{2} & \sf{3} & \sf{4} & \sf{5} & \sf{6}
\end{pNiceArray}^\top, \\[0.8em]
\boldsymbol{c}_2  = \hspace{0.5em}
\begin{pNiceArray}{w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}w{c}{1.2em}}[last-row]
0 & 0 & -1 & -1 & 1 & 0  \\
\sf{1} & \sf{2} & \sf{3} & \sf{4} & \sf{5} & \sf{6}
\end{pNiceArray}^\top.
\end{gathered}
\end{equation}
Upon the removal of reaction $5$ (resp. $6$), $c_2$ (resp. $c_1$) ceases to be a cycle of $$.
Hence, we can choose these two reactions to constitute the subset $_2$:
\begin{align}
\label{eq:S_R1_example}
\S_{\R_1} = \hspace{.3em}
\begin{pNiceArray}{cccc}[first-col, last-row]
\sf{Z}_0 		& -1 & 1 & 0 & 1   \\[0.3em]
\sf{Z}_0^{\, \ast} 	& 1 & -1 & 0 & 0   \\[0.3em]
\sf{Z}_1 		& 0 & 1 & -1 & 1  \\[0.3em]
\sf{Z}_1^{\, \ast} 	& 0 & 0 & 1 & -1  \\[0.3em]
& \sf{1} & \sf{2} & \sf{3} & \sf{4}
\end{pNiceArray}, \\[.5em]
%	& &
(\S_{\R_1})^{-1} = \hspace{.3em}
\begin{pNiceArray}{cccc}[first-col, last-row]
\sf{1}	& 0 & 1 & 1 & 1   \\[0.3em]
\sf{2}	& 0 & 0 & 1 & 1   \\[0.3em]
\sf{3}	& 1 & 1 & 0 & 1  \\[0.3em]
\sf{4}  & 1 & 1 & 0 & 0  \\[0.3em]
& \sf{Z}_0  & \sf{Z}_0^{\, \ast} & \sf{Z}_1 & 	\sf{Z}_1^{\, \ast}
\end{pNiceArray}.
\end{align}
From Eq.~\eqref{eq:G_matrix}, the inverse $(__1)^-1$ fully characterizes the right-inverse $𝔾$ and its columns are the elementary mode of production of the autocatalytic species.
For example, the first column of $(__1)^-1$ is the elementary modes of production of species $Z_0$, it corresponds to the following overall equation in terms of the autocatalytic species only:
\begin{equation}
\boldsymbol{g}_{\sf{Z}_0}: \hspace{1em} \ce{
\sf{Z}_1 + \sf{Z}_1^{\, \ast} ->[]  \sf{Z}_1 + \sf{Z}_1^{\, \ast} + \sf{Z}_0},
\end{equation}
which produces one unit of  $Z_0$ letting the other autocatalytic species unaffected. \\
Furthermore one notices that, when $𝔞 = 1$, the cycles of the full CRN Eq.~\eqref{eq:cycle_example_1} coincide with those of the autocatalytic sub-network Eq.~\eqref{eq:emergent_cycles_example}.
As a consequence, when $𝔞 = 1$, one has $^__2 = 0$.
On the contrary, when $𝔞 = 0$, the cycles of $$ differ from the cycles of the full CRN Eq.~\eqref{eq:cycle_example_0} hence, in that case, the effective stoichiometric matrix Eq.~\eqref{eq:effective_stoichiometric} has non-vanishing entries:
\begin{equation}
\widetilde{\Nabla}^\E_{\R_2} = \hspace{1em}
\begin{pNiceArray}{cc}[first-col, last-row]
\sf{F}_0 & -1 & 1 \\[0.3em]
\sf{F}_1 & 1 & -1  \\[0.3em]
\sf{W}_0 & 1 & -1 \\[0.3em]
\sf{W}_1 & -1 & 1 \\[0.3em]
& \sf{5} & \sf{6}
\end{pNiceArray}.
\end{equation}
Importantly, the case $𝔞 = 1$ in the network Eq.~\eqref{eq:CRN_example} corresponds to the case where a non-catalyzed reaction can be kinetically favored by the mean of a catalyst.
In this situation, the catalytic pathway, combined with its non-catalytic counterpart define, a cycle for the autocatalytic sub-network and the full network.
In Appendix B, we extend this result in the general case by showing that, if two pathways of reaction in a network are related in a similar way (\emph{i.e.} one being the catalyzed version of the other, with an arbitrary number of linear intermediates), then, a cycle for both the autocatalytic sub-network and the full network can be constructed. \\
When $𝔞 = 0$, however, the external reactant and product in the catalyzed pathways and their non-catalyzed counterparts are different, hence, they are \emph{a priori} two different processes.

\subsection{Dynamics and steady-state}

The dynamics of the autocatalytic species Eq.~\eqref{eq:dyn_autocat} can be written,
\begin{equation}
\label{eq:dyn_free_Z}
\d_t [\z] ~=~ \S_{\R_1} \cdot \j^{\R_1} + \S_{\R_2} \cdot \j^{\R_2} + \I^\Z.
\end{equation}
Hence, using the invertibility of $__1$ we can express the the elementary fluxes of reactions $_1$, $ȷ^_1$, as a function of the other parameters.
%	\j^{\R_1}
%	= \left( \S_{\R_1} \right)^{-1} \cdot \left(
%	\d_t [\z]
%	- \S_{\R_2} \cdot \j^{\R_2}
%	- \I^\Z
%	\right).
Plugging this expression in the dynamics of the external species Eq.~\eqref{eq:dyn_external} reads:
\begin{multline}
\label{eq:dyn_external_2}
\d_t [\e] ~=~
\d_t \Delta \e (\z) + \widetilde{\Nabla}^\E_{\R_2} \cdot \j^{\R_2}
+ \NablaEnR \cdot \v  \\
+ \I^{\E}
- \Nabla^\E_{\R_1} \cdot (\S_{\R_1})^{-1} \cdot \I^\Z,
\end{multline}
\begin{equation}
\label{eq:Delta_e}
\Delta \e ( \z ) = \Nabla^\E_{\R_1} \cdot (\S_{\R_1})^{-1} \cdot [\z].
\end{equation}
Equation \eqref{eq:dyn_external_2} reveals the presence of hidden conservative dynamics for the external species due to production of the autocatalytic species along the elementary modes of Eq.~\eqref{eq:G_matrix}.
The flux along the cycles of $$, on the other hand, requires a non-conservative influx of external species, $ ^__2 ·ȷ^_2$, which vanishes when the cycles of $$ are also cycles for the full network.
These two complementary phenomena will play an important role for the thermodynamics associated to the autocatalytic sub-network, they are represented in Fig.~\ref{fig:schema_system} as the two sources of external species in the autocatalytic sub-network.

\subsubsection{Condition for equilibrium in the sub-network}

Assuming the autocatalytic network is undriven, $^= 0$, and the external species are kept at an arbitrary chemical potential $Μ_$ (\emph{i.e.} chemostatted), we now look for the condition to reach detailed balance in the sub-network.
When this arises the sub-network is unconditionally detailed balanced which means that the chemical potentials of the autocatalytic species, $Μ_^eq$, satisfies the detailed balance property for any state $Μ_$:
\begin{equation}
\label{eq:detail_balance_Z}
\Mu_\Z^\text{eq} \cdot \S + \Mu_\E \cdot \NablaER =\boldsymbol{0}.
\end{equation}
Writing Eq.~\eqref{eq:detail_balance_Z} for the subsets $_i$ ($i=1$ or $2$) yields $Μ_^eq   ·  __i + Μ_  ·  ^__i =0$, hence, the non-singularity of $__1$ implies that the equilibrium state should be:
\begin{equation}
\label{eq:Eq_state_R1}
\Mu_\Z^\text{eq} = - \Mu_\E \cdot \Nabla^\E_{\R_1} \cdot \left( \S_{\R_1} \right)^{-1}.
\end{equation}
However, the presence of cycles in the sub-network may prevent the autocatalytic sub-network from reaching the equilibrium state.
Plugging the expression Eq.~\eqref{eq:Eq_state_R1} into the equation for $_2$ gives the condition for the feasibility of equilibrium.
Namely, the autocatalytic sub-network is unconditionally detailed balance if, and only if:
\begin{equation}
\label{eq:condition_eq}
\widetilde{\Nabla}^\E_{\R_2} = \boldsymbol{0}.
%	\Nabla^\E_{\R_2}  = \Nabla^\E_{\R_1} \cdot \left( \S_{\R_1} \right)^{-1} \cdot \S_{\R_2}.
\end{equation}
In light of the previous discussion, the autocatalytic sub-network is unconditionally detailed balance if, and only if, all the cycles of $$ are also cycles for the full CRN.
From a dynamical standpoint, the criterion for unconditional detailed balance in the autocatalytic CRN is also easily interpreted by looking at Eq.~\eqref{eq:dyn_external_2}.
Indeed, if the autocatalytic sub-network is undriven and $^__2 = 0$, the conservative dynamics, which is fully determined by the state of the autocatalytic species $$, solely describes the influx of external species fueling the autocatalytic reactions.
Furthermore, as $$ is unconstrained from the absence of conservation laws in the autocatalytic sub-networks, nothing prevents the autocatalytic sub-network from reaching the equilibrium state Eq.~\eqref{eq:Eq_state_R1}.

\subsubsection{Decomposition of the fluxes and steady-state}

If the external flux vector of the autocatalytic species, $^$, is non-vanishing, the autocatalytic network will necessary be away from equilibrium, regardless of $^__2$.
Calling $$ (resp. $$) the subset of the autocatalytic species having a non-vanishing (resp. vanishing) entry in external flux vector we have, for all $x ∈$, $I^x ≠0$ (resp. for all $y ∈$, $I^y = 0$).
We can separate the dynamics of the autocatalytic sub-network coupled to the environment to the one that is free:
\begin{gather}
\d_t [\x] = \S^\X \cdot \j + \I^\X ~,\\[0.5em]
\label{eq:dyn_Y}
\d_t [\y] = \S^\Y \cdot \j ~.
\end{gather}
We can decompose the elementary flux vector on the basis defined by the elementary modes of Eq.~\eqref{eq:G_matrix} and the cycles:
\begin{equation}
\label{eq:decomposition_fluxes}
\boldsymbol{j} =
\sum_\varepsilon J^{\varepsilon}(t) \, \boldsymbol{c}_\varepsilon
\sum_{z \in \Z} \J^{z}(t) \,  \boldsymbol{g}_z,
\end{equation}
where $J^ε (t)$ stands for the macroscopic flux along the cycle $c_ε$ and $^z (t)$ is the macroscopic flux along the elementary mode of production $g_z$.
The decomposition Eq.~\eqref{eq:decomposition_fluxes} is valid even if the network is not at steady-state.
The constraint on the right inverse implies that the restriction of the elementary fluxes on $_2$ stems uniquely from the cycles of the sub-network:
\begin{equation}
\boldsymbol{j}^{\R_2} =
\sum_\varepsilon J^{\varepsilon}(t) \, \boldsymbol{c}_\varepsilon^{\R_2}.
\end{equation}
%in particular one can easily express the macroscopic curents along the emergent cycle noticing:
%	content
Steady-state of Eq.~\eqref{eq:dyn_Y} consists of an elementary flux vector, $ȷ$, in the the right nullspace of the sub-matrix $^$.
The latter is spanned by all the cycles of $$ and by the elementary modes of production of the $$ species as, by definition, they leave the $$ species unaffected:
\begin{equation}
\S^\Y \cdot \mathbb{G}_\X = \boldsymbol{0}.
\end{equation}
Hence, the macroscopic fluxes along the elementary mode of production of the $$ species tend to vanish,  $^(t) →0$, while those along the $$ species tend to their steady values, $^(t) →^  $.
From this, we can write the general expression of the steady fluxes as:
\begin{equation}
\label{eq:ss_fluxes}
\overline{\j} =
\sum_\varepsilon \overline{J}^{\, \varepsilon} \, \boldsymbol{c}_\varepsilon
\sum_{x \in \X} \overline{\J}^{\, x}\,  \boldsymbol{g}_x,
\end{equation}

Here it is crucial to stress that, in the decomposition Eqs.~\eqref{eq:decomposition_fluxes} and \eqref{eq:ss_fluxes}, only the second part, originating from the elementary modes of production, tends to produce autocatalytic species.
The other part (consisting of cycles) keeps the number of autocatalytic species constant in time.
More precisely, having $𝒥^z > 0$ implies that species $z$ is effectively produced.
On the contrary, if $𝒥^z < 0$, the corresponding species is instead consumed because its mode of production is performed in the backward sense.
This case corresponds to \emph{recessive} autocatalysis mentioned in \cite{Peng2022}.
As a consequence, Eq.~\eqref{eq:ss_fluxes} makes clear that, to sustain steady production (or consumption) of the autocatalytic species $z$, the latter should be coupled with the environment by the mean of a non-vanishing external flux, $I^z ≠0$; then, it should belong to the subset $$.
At steady-state, when a species is produced (resp. consumed) in the sub-network it should be extracted (resp. injected) to keep its concentration constant, \emph{i.e.}, $^= - I^$.
\vspace{2em}

\paragraph{\textbf{Example 4.}}
In the CRN Eq.~\eqref{eq:CRN_example}, assuming the external species are chemostatted and no external flux couples the autocatalytic species with the  environment ($^= 0$), the long-time behavior of the autocatalytic CRN  is governed by $𝔞$.
When $𝔞 = 1$, $^__2 = 0$ hence, the autocatalytic sub-network converges to an equilibrium state: \\
\begin{minipage}{0.15\textwidth}
\begin{equation*}
\begin{gathered}
\mu_{\sf{Z}_0}^\text{eq} = \mu_{\sf{F}_1} - \mu_{\sf{W}_1}  , \\
\mu_{\sf{Z}_1}^\text{eq} = \mu_{\sf{F}_0} - \mu_{\sf{W}_0} ,
\end{gathered}
\end{equation*}
\end{minipage}
\hspace{1em}
\begin{minipage}{0.3\textwidth}
\begin{equation}
\label{eq:Eq_Example}
\begin{gathered}
\mu_{\sf{Z}_0^{\, \ast}}^\text{eq} = \mu_{\sf{F}_0} + \mu_{\sf{F}_1} - \mu_{\sf{W}_1}, \\
\mu_{\sf{Z}_1^{\, \ast}}^\text{eq} = \mu_{\sf{F}_0} + \mu_{\sf{F}_1} - \mu_{\sf{W}_0} ,
\end{gathered}
\end{equation}
\end{minipage}

\vspace{1em}
\noindent and the elementary fluxes tend to vanish (see Fig. \ref{fig:undriven_CRN}\textsf{\textbf{a}}).
When $𝔞 = 0$, on the contrary,  $^__2 ≠0$ and no solution of Eq.~\eqref{eq:detail_balance_Z} exists in general.
In this case, the autocatalytic CRN reaches a NESS given by a linear combination of its two cycles Eq.~\eqref{eq:emergent_cycles_example}:
\begin{equation}
\overline{\j} = \overline{J}^1 \, \boldsymbol{c}_1 + \overline{J}^2 \, \boldsymbol{c}_2.
\end{equation}
From Eq.~\eqref{eq:c_R1}, the elementary fluxes are completely determined by $ȷ^_2$, which converges to a non-zero value (see Fig. \ref{fig:undriven_CRN}\textsf{\textbf{c}}).
Note that all the details regarding the Figures are specified in the Appendix C. \\
Setting now a non-zero external flux on species $Z_0$ and $Z_1$, we have: $= {Z_0,   Z_1 }$ and $= {Z_0^ ∗,   Z_1^  ∗ }$.
Hence, for both values of $𝔞$, steady-state for the dynamics associated to the $$ species is given by the modes of production of the $$ species and by the cycles of $$:
\begin{equation}
\overline{\j} =
\overline{J}^1 \, \boldsymbol{c}_1 + \overline{J}^2 \, \boldsymbol{c}_2
\overline{\J}^{\sf{Z}_0} \, \boldsymbol{g}_{\sf{Z}_0} + \overline{\J}^{\sf{Z}_1} \, \boldsymbol{g}_{\sf{Z}_1}.
\end{equation}
Enforcing $^Z_i > 0$, implies that $Z_i$ is steadily produced in the sub-network and then released in the environment to keep its concentration constant in the system at steady-state.
This condition can enforced by choosing an appropriate coupling with environment, $I^Z_i$ (see the discussion section and Eq.~\eqref{eq:condition_J_geq_0}).

\begin{figure*}[t]
\centering
\includegraphics{Undriven_CRN_Panel.eps}
\caption{\label{fig:undriven_CRN}
The behavior of the autocatalytic sub-network of Eq.~\eqref{eq:Eq_Example} when $\I^\X = \boldsymbol{0}$ is governed by the value of $a$.
When $a = 1$, \textbf{\textsf{(a)}} \& \textbf{\textsf{(b)}}, the topological condition Eq.~\eqref{eq:condition_eq} is fulfilled, and the autocatalytic CRN reaches equilibrium.
In that case:
\textbf{\textsf{(a)}} the elementary fluxes vanish ;
\textbf{\textsf{(b)}}  the thermodynamic potential (solid orange) is a strictly decreasing function whose time-derivative is given by the opposite of the EPR (solid purple).
When $a = 0$, \textbf{\textsf{(c)}} \& \textbf{\textsf{(d)}}, the autocatalytic CRN aligns with its cycles:
\textbf{\textsf{(c)}}
the elementary fluxes converge to their steady-state solution ;
\textbf{\textsf{(d)}} the thermodynamic potential (solid orange)
results from the contribution of the EPR (solid purple) and the work needed to sustain the cycles (dotted green).
%	 	Once the NESS is established, all the work is dissipated.
The parameters used to generate these plots are specified in Appendix C.
\end{figure*}

\subsection{The case $\boldsymbol{|\Z| = |\R|}$}

Often, for small networks, the number of autocatalytic reactions matches the number of autocatalytic species, \emph{i.e.}, $|| = ||$.
In particular, it is the case for all the elementary cores described in Ref.~\cite{Blokhuis2020}.
%When the number of autocatalytic species is equal to the number of autocatalytic reactions the above discussion simplifies drastically.
In that case, $$ becomes a square non-singular matrix, its right inverse identifies to the usual inverse of $$:
\begin{equation}
\label{eq:inverse_S_square}
\invS \equiv \S^{-1}.
\end{equation}
Note that Eq.~\eqref{eq:inverse_S_square} can be seen as the limiting case of the right-inverse Eq.~\eqref{eq:G_matrix} when the sub-network is square.
Furthermore, $$ has an empty right null-space then it has no cycles.
This implies that $_1$ is the same as $$ and $__1 ≡$, while $_2 = ∅$.
Hence, square sub-networks are always unconditionally detailed balanced:
\begin{equation}
\Mu_\Z^\text{eq} = - \Mu_\E \cdot \NablaER \cdot \S^{-1}.
\end{equation}
In what follows, the case $|| = ||$ is obtained by taking $(__1)^-1 ≡^-1$ and  $^__2 ≡0$.

\section{Thermodynamic of the autocatalytic sub-network}

\subsection{Thermodynamic potential}

The total entropy production rate (EPR), $T Σ̇$, can be written as sum of the EPR of the autocatalytic reactions, $TΣ̇_$, and the one of the additional reactions, $T Σ̇_$:
\begin{align}
T \dot{\Sigma}_{\R}  = - \Mu \cdot \Nabla_{\R} \cdot \j ,
& &
T \dot{\Sigma}_{\nR} = - \Mu_{\E} \cdot \NablaEnR \cdot \v.
\end{align}
Introducing the fundamental chemical forces,
\begin{equation}
\label{eq:chemical_force}
\F = \Mu - \Mu_{\Ep} \cdot \M \, ,
\end{equation}
the two contributions of the EPR can written as \cite{Avanzini2020, Avanzini2021, Avanzini2022}:
\begin{align}
\label{eq:EPR_R}
T \dot{\Sigma}_{\R} = - \F \cdot \Nabla_{\R} \cdot \j,
&  &
%	\label{eq:EPR_nR}
T \dot{\Sigma}_{\nR} = - \F_{\Ef} \cdot \Nabla^{\Ef}_{\nR} \cdot \v.
\end{align}
From Eq.~\eqref{eq:EPR_R} it is clear that, when $= 0$, the full CRN reaches an equilibrium state in which the total EPR vanishes.
Moreover as, by construction, $_ = 1$, the fundamental forces associated to the potential species vanish. \\

\paragraph{\textbf{Example 5.}}
The fundamental forces of network Eq.~\eqref{eq:CRN_example} with $= {F_0,   F_1 }$ are
\begin{equation}
\begin{gathered}
\F_\E = \begin{pNiceArray}{w{c}{.5cm}w{c}{.5cm}w{c}{1.5cm}w{c}{1.5cm}}[last-row]
& 0
& \mu_{\sf{W}_0} -  \mu_{\sf{F}_0}
& \mu_{\sf{W}_1} -  \mu_{\sf{F}_1} \\
\sf{F}_0 & \sf{F}_1 & \sf{W}_0 & \sf{W}_1
\end{pNiceArray}, \\[0.5em]
\F_\Z = \begin{pNiceArray}{cccc}[last-row]
\mu_{\sf{Z}_0}
& \mu_{\sf{Z}_0^{\, \ast}} - \mu_{\sf{F}_0}
& \mu_{\sf{Z}_1}
& \mu_{\sf{Z}_1^{\, \ast}} - \mu_{\sf{F}_1} \\
\sf{Z}_0 & \sf{Z}_0^{\, \ast} & \sf{Z}_1 & \sf{Z}_1^{\, \ast}
\end{pNiceArray}.
\end{gathered}
\end{equation}
As expected, the fundamental forces associated with the two potential species vanish,
$_ = 0$.
One can check that $= 0$ corresponds to equilibrium in the full network.
\vspace{2em}

When all the external fluxes vanish ($I = 0$), the system is closed and the Gibbs free-energy takes its well-known form:
\begin{equation}
\label{eq:close_G}
G(\e, \z) = \sum_{e \in \E} [e] \left( \mu_e - RT \right) + \sum_{z \in \Z} [z] \left( \mu_z - RT \right).
\end{equation}
Now, if some entries in $$ are non-vanishing, the system exchanges \emph{moieties} with the environment:
\begin{equation}
\boldsymbol{m}(\e, \z) = \M \cdot \begin{pmatrix}
[\e] \\[0.5em] [\z]
\end{pmatrix} = \M_\E \cdot [\e] + \M_\Z \cdot [\z].
\end{equation}
Furthermore, the thermodynamic potential,
\begin{equation}
\label{eq:G_def}
\G(\e, \z) = G(\e, \z) - \Mu_{\Ep} \cdot \boldsymbol{m}(\e, \z) + \F_{\E} \cdot \Delta \e (\z),
\end{equation}
captures the behavior of the autocatalytic sub-network when the system is open.
The first two terms in the RHS of Eq.~\eqref{eq:G_def} were previously derived in Refs.~\cite{Rao2016, Avanzini2021, Avanzini2022}: they correspond to the semigrand free energy of an open CRN.
The last term, however, is new: it corresponds to the conservative energy needed to fuel the production of autocatalytic species along their mode of production by the mean of the conservative influx of external species.
Note that, in the non-ideal case (\emph{cf.} Appendix A), the Gibbs free-energy of the closed system has an additional term accounting for the interactions (\emph{cf.} Eq.~\eqref{eq:close_G_ni}).

\subsection{Minimum of $\G$}

The optimum of $$ with respect to the autocatalytic species is:
%	\label{eq:Minimum_G}
%	\partial_{\z} \G = \Mu^\text{eq}_\Z - \Mu_{\Ep} \cdot \M_\Z + \F_\E \cdot  \Nabla^\E_{\R_1} \cdot (\S_{\R_1})^{-1}
%	= \boldsymbol{0}.
\begin{align}
\label{eq:eq_state_1}
\Mu^\text{eq}_\Z
& = \Mu_{\Ep} \cdot \M_\Z - \F_\E \cdot   \Nabla^\E_{\R_1} \cdot (\S_{\R_1})^{-1} \\[0.5em]
\label{eq:eq_state_2}
& = - \Mu_\E \cdot \Nabla^\E_{\R_1} \cdot (\S_{\R_1})^{-1}.
\end{align}
which is in accordance with Eq.~\eqref{eq:Eq_state_R1}.
To derive Eq.~\eqref{eq:eq_state_2} from Eq.~\eqref{eq:eq_state_1}, we used the definition of the chemical forces Eq.~\eqref{eq:chemical_force} along with Eq.~\eqref{eq:Nabla_Ep_R} to express the stoichiometric matrix of the potential species.
Note that the additional term in the potential Eq.~\eqref{eq:G_def} compared to the semigrand free-energy is crucial to obtain the correct expression for the equilibrium state in the sub-network.
Furthermore, one has:
\begin{align}
\label{eq:DeltaG_1}
\Delta \G & = \G(\e, \z) - \G(\e, \z_\text{eq}) \\[0.5em]
&  = RT \, \sum_{z \in \Z} [z] \log \left(  \dfrac{[z]}{[z]_\text{eq}} \right)  - \left( [z] - [z]_\text{eq} \right) \\[0.5em]
\label{eq:DeltaG_2}
& = RT  \, \mathcal{L} \left( \z \Vert \z_\text{eq} \right),
\end{align}
where we introduced the relative entropy between two non-normalized distribution $ℒ ( 𝐚 ‖b ) = ∑_i a_i log(a_i/b_i) - (a_i - b_i)$.
The latter is positive and vanishes if, and only if, the two distributions are the same.
Hence, the equilibrium state is the global minimum of $$ and the convexity of $$ derives from the convexity of the relative entropy.

%\paragraph{Example 4.}
%	\F_{\Ef} =
%	\begin{pNiceArray}{cc}%[last-row]
%		\mu_{\sf{W}_0}-\mu_{\sf{F}_0}, & \mu_{\sf{W}_1}-\mu_{\sf{F}_1} 	 \\
%		%		 \sf{W}_0 & \sf{W}_1
%	\end{pNiceArray},
%	& \quad &
%	\F_{\Z} =
%	\begin{pNiceArray}{cccc}%[last-row]
%		\mu_{\sf{Z}_0}, & \mu_{\sf{Z}_0^{\, \ast}} - \mu_{\sf{F}_0}, &
%		\mu_{\sf{Z}_1}, & \mu_{\sf{Z}_1^{\, \ast}} - \mu_{\sf{F}_1} \\[0.5em]
%		%		\sf{Z}_0 & \sf{Z}_0^{\, \ast} & \sf{Z}_1 & \sf{Z}_1^{\, \ast}
%	\end{pNiceArray}.

\subsection{Decomposition of the EPR}

The time evolution of $$ can be written as:
\begin{widetext}
\begin{equation}
\label{eq:evol_G_1}
\d_t \G = - T \dot{\Sigma}
- \d_t \Mu_{\Ep} \cdot \boldsymbol{m}
+ \d_t \F_{\E} \cdot  \Delta \e
+ \F_{\E} \cdot \left( \I^{\E} + \d_t \Delta \e \right)
+ \F_\Z \cdot \I^\Z.
\end{equation}
\end{widetext}
The terms in the RHS of Eq.~\eqref{eq:evol_G_1} can be gathered as follows:
\begin{gather}
\dot{W}_\text{driv} =  - \d_t \Mu_{\Ep} \cdot \boldsymbol{m} + \d_t \F_{\E} \cdot  \Delta \e, \\[0.5em]
\label{eq:W_nc_E}
\dot{W}_{\text{nc}, ~ \E} = \F_{\E} \cdot \left( \I^{\E} + \d_t \Delta \e \right), \\[0.5em]
\dot{W}_{\text{nc}, ~ \Z} =  \F_{\Z} \cdot \I^\Z,
\end{gather}
such that the time evolution of the free-energy reads
\begin{equation}
\label{eq:second_law_dt_G}
\d_t \G = - T \dot{\Sigma} + \dot{W}_\text{driv} + \dot{W}_{\text{nc}, ~ \E} + \dot{W}_{\text{nc}, ~ \Z}.
\end{equation}
In the decomposition of the EPR: $Ẇ_driv$ represents the driving work rate required to change the equilibrium state  in the autocatalytic sub-network Eq.~\eqref{eq:Eq_state_R1} by varying the chemical potential of the external species.
Because the equilibrium state of the sub-network depends on all the external species, the external force species also contribute to the driving work.
The two non-conservative work rates, $Ẇ_nc$,
encompass the energetic cost to sustain chemical fluxes in the sub-network, putting it away from equilibrium.
Because $_ = 0$, the non-conservative work brought by the potential species vanishes, hence the choice for their name.
With Eq.~\eqref{eq:G_def}, the decomposition of the EPR Eq.~\eqref{eq:second_law_dt_G} constitutes the central result of the present work.
In non-ideal systems, an additional driving work rate pops up when one tunes the interactions by the mean of externally controlled parameters (\emph{cf.} Eq.~\eqref{eq:driving_work_int}).

Importantly, the decomposition of the EPR is valid for an arbitrary coupling with environment $I$.
In particular, it is not the required that the autocatalytic species having a non-vanishing external flux, \emph{i.e.}, the $$ species, are chemostatted.
Subjecting, for example, species $x$ to degradation, its entry in the external flux vector can be written as
\begin{equation}
\label{eq:external_flux_X_species}
I^x = \phi^x_\text{in} (\e, \z) - \phi^x_\text{out} (\e, \z)
\end{equation}
where $ϕ^x_in$ and $ϕ^x_out$ are two positive functions corresponding to the input and output flux on species $x ∈$, respectively.

\begin{figure*}[ht]
\centering
\includegraphics{Efficiency_Panel.eps}
\caption{\label{fig:efficiency}
During the dynamics of an out-of-equilibrium autocatalytic network, driven production compete with energy storage.
\textbf{\textsf{(a)}} Initially, the autocatalytic sub-network discharges, and $\d_t \G < 0$ (gray area), until it eventually reaches the charging phase where free-energy start to accumulate.
The charging progressively stops as the sub-network is entering in the maintenance phase and NESS settles.
\textbf{\textsf{(b)}} The power injected to sustain the cycles of the sub-network that is not dissipated, $\dot{W}_\text{cyc} - T \dot{\Sigma}_\R$, is shared to store free-energy and to drive flows along the cycles.
During the disarging phase, $\dot{W}_\text{cyc} - T \dot{\Sigma}_\R$ is not sufficient to sustain the production of autocatalytic species, as consequence, free-energy must be released from the sub-network.
In the charging phase however, $\dot{W}_\text{cyc} - T \dot{\Sigma}_\R$ drives both the charging and the production of autocatalytic species.
Finally, in the maintenance phase, the power driving the cycles is fully dedicated to sustain the production of autocatalytic species.
\textbf{\textsf{(c)}} The efficiency of the energy storage attains its maximum during the charging phase, where a significant portion of $\dot{W}_\text{cyc} - T \dot{\Sigma}_\R$ serves to increase free-energy, it drops during the maintenance phase as fewer energy is stored.
On the contrary, the efficiency of driven production is increasing and attains its maximum when NESS is established.
\end{figure*}

\subsection{Autonomous autocatalytic network}

In this section, we address the case where the chemical potentials of the external species are kept constant.
Or the case in which there exists a time scale separation such that, in the timescale of the autocatalytic CRN, the chemical potentials the external can be regarded as constant.
In such cases, the driving work rate vanishes:
\begin{equation}
\dot{W}_\text{driv} = 0.
\end{equation}
We refer to this case as \emph{autonomous} autocatalytic sub-network.

\subsubsection{Thermodynamic cost of production}

From Eq.~\eqref{eq:dyn_external_2}, keeping the external species at a constant concentration imposes that the external flux vector is
\begin{equation}
\label{eq:I_E_chemostat}
\I^\E =
- \d_t \Delta \e
- \widetilde{\Nabla}^\E_{\R_2} \cdot \j^{\R_2}
- \NablaEnR \cdot \v,
+ \Nabla^\E_{\R_1} \cdot (\S_{\R_1})^{-1} \cdot \I^\Z.
\end{equation}
Hence, the EPR decomposition Eq.~\eqref{eq:evol_G_1} becomes
\begin{equation}
\label{eq:second_law_autonomous}
\d_t \G =
- T \dot{\Sigma}_{\R}
\dot{W}_\text{cyc}
%	\F_\E \cdot \widetilde{\Nabla}^\E_{\R_2} \cdot \j^{\R_2}
\wprod
%	\left(
%	\F_\E \cdot \Nabla^\E_{\R_1} \cdot (\S_{\R_1})^{-1}  + \F_\Z
%	\right) \cdot \I^\Z
\end{equation}
In the RHS of Eq.~\eqref{eq:second_law_autonomous}, the last term identifies to the work rate needed to sustain the production of the autocatalytic species having a non-vanishing external flux (\emph{i.e.} species $$) by the mean of their elementary modes of production:
\begin{equation}
\label{eq:w_prod}
\wprod =  \left(
\F_\E \cdot \Nabla^\E_{\R_1} \cdot (\S_{\R_1})^{-1}  + \F_\Z
\right) \cdot \I^\Z.
\end{equation}
%It can be interpreted as follow: $\F_\Z \cdot \I^\Z$ corresponds to the force needed to remove the excess of $\Z$ species that are created to maintain their concentration constant ;
Where $_·^__1 ·(__1)^-1 ·^$ quantifies the uptake of external species fueling the elementary modes of production.
Importantly, to fulfill mass conservation, all columns of $^__1 ·(__1)^-1$ must have, at least, one strictly negative entry.
On the other hand,
\begin{equation}
\dot{W}_\text{cyc} = - \F_\E \cdot \widetilde{\Nabla}^\E_{\R_2} \cdot \j^{\R_2},
\end{equation}
accounts for the work rate needed to sustain the cycles of $$.
It vanishes if $^__2 = 0$ because, when this happens, all the cycles of $$ are also cycles for the full network.

As expected, when no external flux couples the autocatalytic species with the environment the work rate associated to the production.
When this occurs, no steady production of the autocatalytic species can be maintained and Eq.~\eqref{eq:second_law_autonomous} becomes:
\begin{equation}
\label{eq:dtG_relaxation_eq}
\d_t \G = - T \dot{\Sigma}_{\R} - \F_\E \cdot \widetilde{\Nabla}^\E_{\R_2} \cdot \j^{\R_2} ~.
\end{equation}
Hence, when $^__2 = 0$, the autocatalytic CRN necessarily evolves to its equilibrium state as $_̣t = - T Σ̇_ < 0$ (see Fig. \ref{fig:undriven_CRN}\textsf{\textbf{b}}).

This finally proves that $Δ$ is a Lyapunov function of the dynamics of the autocatalytic sub-network.
Its minimum is the equilibrium state Eq.~\eqref{eq:Eq_state_R1} and is reached when the sub-network is undriven and when condition Eq.~\eqref{eq:condition_eq} is fulfilled.
Notably, from Eq.~\eqref{eq:dtG_relaxation_eq}, we see that the sub-network can relax locally to equilibrium if $_ = 0$.
%In fact, in that case, the full network converges to its equilibrium state.

\subsubsection{Non-equilibrium steady-state}

Injecting the decomposition of the elementary fluxes Eq.~\eqref{eq:decomposition_fluxes} in the EPR Eq.~\eqref{eq:EPR_R} reads,
\begin{widetext}
\begin{align}
\label{eq:EPR_decomposition}
T \dot{\Sigma}_\R ~
& = ~
- \sum_\varepsilon J^\varepsilon (t) ~ \F_\E \cdot \widetilde{\Nabla}^\E_{\R_2} \cdot \boldsymbol{c}_\varepsilon^{\R_2}
~ - ~
\left. \left. \sum_{z \in \Z} \J^z (t) \right[ \F_\E \cdot \NablaER \cdot \boldsymbol{g}_z + \mathcal{F}_z \right]
\end{align}
\end{widetext}
where the first sum runs over the cycles of $$.
It describes the amount of dissipation along these cycles and can be rewritten as $- _·^__2 ·j^_2$.
Introducing the macroscopic affinities along the  cycles and along the elementary modes of production, respectively defined as
\begin{gather}
A_\varepsilon = - \F_\E \cdot \widetilde{\Nabla}^\E_{\R_2} \cdot \boldsymbol{c}_\varepsilon^{\R_2}, \\[0.5em]
%	& \quad &
\mathcal{A}_z = -\F_\E \cdot \NablaER \cdot \boldsymbol{g}_z - \mathcal{F}_z ~,
\end{gather}
the EPR is rewritten as
\begin{equation}
\label{eq:EPR_affinities}
T \dot{\Sigma}_\R ~ = ~ \sum_\varepsilon J^\varepsilon (t) \, A_{\varepsilon}
~ + ~
\sum_{z \in \Z} \J^z (t) \, \mathcal{A}_{z} .
\end{equation}
Importantly, this decomposition is valid arbitrarily far from equilibrium and not only at steady-state, and corresponds to the well-known Hill-Schnakenberg decomposition \cite{Hill1977, Schnakenberg1976}.
%Once NESS is established, only the macroscopic currents along the modes of the $\X$ species remain:
%	\label{eq:ss_EPR}
%	T \overline{\dot{\Sigma}_\R} ~
%	& = ~
%	- \F_\E \cdot \widetilde{\Nabla}^\E_{\R_2} \cdot \overline{\j}^{\, \R_2}
%	+  \sum_{x \in \X} \overline{\J}^x \, \mathcal{A}_x .
Once the autocatalytic sub-network has settled a NESS, all the work ends up dissipated,
\begin{equation}
T \overline{\dot{\Sigma}_{\R}} =
\overline{\dot{W}}_\text{cyc} +
\overline{\dot{W}}_\text{prod},
\end{equation}
and $_̣t = 0$.

\subsection{Energetics of autonomous autocatalytic networks}

\subsubsection{Energy storage in the autocatalytic sub-network}
%	T \Sigma _\text{cyc} = - \int \d t \hspace{0.5em} \F_\E \cdot \widetilde{\Nabla}^\E_{\R_2} \cdot \overline{\j}^{\R_2}

If from a certain time in the transient dynamics, one finds $_̣t > 0$, free-energy starts to accumulate in the sub-network.
During the first stage of this process a significant portion of the work is converted in free-energy, this phase corresponds to the \emph{charging phase}.
Upon approaching steady state, the sub-network enters in its \emph{maintenance phase}, where the work tends to be totally dissipated to preserve the accumulated free-energy, see Fig.~\ref{fig:efficiency}\textbf{\textsf{a}}-\textbf{\textsf{b}}.
Integrating Eq.~\eqref{eq:second_law_autonomous} where energy storage occurs reads,
\begin{equation}
\label{eq:integrated_second_law}
\Delta \G = - T \Sigma_{\R} + W_\text{cyc} + \Wprod > 0.
\end{equation}
Using chemical systems for \emph{energy storage} was first proposed in Ref.~\cite{Ragazzon2018}.
The efficiency of this process can obtained comparing the gain in free-energy to the amount of work spend \cite{Penocchio2019}:
\begin{equation}
\label{eq:efficiency_ES}
\eta_\text{es} = \dfrac{\Delta \G}{\Wprod + W_\text{cyc}}
= 1 - \dfrac{T \Sigma_{\R}}{\Wprod + W_\text{cyc}}.
\end{equation}
The maximum efficiency is attained during the charging phase and, as the system settles its NESS, the efficiency drops to zero as shown in Fig.~\ref{fig:efficiency}\textbf{\textsf{c}}.
Note that, starting from the equilibrium state Eq.~\eqref{eq:Eq_state_R1} and setting off the external fluxes, the autocatalytic sub-network is ensured to store energy.
In this case, the energy stored is given directly by Eqs.~\eqref{eq:DeltaG_1}-\eqref{eq:DeltaG_2}.

\subsubsection{Thermodynamically favored production}

\begin{figure}[t!]
\advance\leftskip+0em
\includegraphics{Driven_Replication_Panel.eps}
\caption{\label{fig:DR_ss}
During driven production, the cycles of $\S$ drive the flux along the elementary modes of the autocatalytic species that are externally controlled.
As a consequence, the extraction of the autocatalytic species that are produced at steady-state releases energy to the environment.
\textbf{\textsf{(a)}} In that case, the efficiency of driven production at steady-state is given by Eq.~\eqref{eq:ss_efficiency_DR} as function of the forces associated to $\sf{W}_0$ and $\sf{W}_1$.
\textbf{\textsf{(b)}} \& \textbf{\textsf{(c)}} In parallel, it is possible to replicate species $\sf{Z}_0$ and $\sf{Z}_1$, externally controlled by insuring that their macroscopic flux $\mathcal{J}^{\sf{Z}_0}$ and $\mathcal{J}^{\sf{Z}_1}$ are both positive.
\end{figure}

As previously stated, to sustain steady production (or consumption) of autocatalytic species they should be coupled with the environment.
During \emph{driven production}, the power injected to sustain the cycles of $$ at steady state also drives the fluxes along the modes of production of the $$ species,
\begin{align}
\label{eq:DS_cyc_repl}
\overline{\dot{W}}_\text{prod} < 0,
& \quad &
\overline{\dot{W}}_\text{cyc} > 0.
\end{align}
It means that the extractio (or injection) of autocatalytic species through $^$ releases energy to the environment.
Furthermore, if one seeks to increase \emph{all} the $$ species in the environment, we should enforce, in addition, $𝒥^> 0$, element-wise.
As explained in the discussion section, this seemingly difficult condition can be easily enforced at steady-state by choosing an appropriate $^$ (\emph{cf.} Eq.~\eqref{eq:condition_J_geq_0}). \\
When driven production takes place at steady-state, the power performed by the environment to sustain the cycles that is not dissipated serves to drive the steady fluxes along the modes of production, \emph{i.e.},
\begin{equation}
- \overline{\dot{W}}_\text{prod} = \overline{\dot{W}}_\text{cyc} - T \overline{\dot{\Sigma}_\R} .
\end{equation}
From this, the efficiency of driven production is:
\begin{align}
\label{eq:ss_efficiency_DR}
\overline{\eta}_\text{dp}  =
\dfrac{
\overline{\dot{W}}_\text{cyc}
1 - \dfrac{
T \overline{\dot{\Sigma}_\R}
\overline{\dot{W}}_\text{cyc}
\end{align}
\vspace{2em}

\paragraph{\textbf{Example 6.}}
In Fig.~\ref{fig:DR_ss}\textbf{\textsf{a}}, the efficiency of driven production of the network Eq.~\eqref{eq:CRN_example} when $𝔞 = 0$ is represented with respect to the forces associated to its two external force species ($W_0$ and $W_1$) when $Z_0$ and $Z_1$ are the $$ species, coupled with environment (in the figure they are chemostatted).
As shown in Figs.~\ref{fig:DR_ss}\textbf{\textsf{b}}-\textbf{\textsf{c}}, in the region shown in Fig.~\ref{fig:DR_ss}\textbf{\textsf{a}}, the macroscopic fluxes along the modes of the two $$ species are positive.
Hence, $Z_0$ and $Z_1$ are effectively produced in the sub-network and then released in the environment (to keep their concentration constant in the sub-network).
The parameters used to generate the plots are available in Appendix C.
\vspace{2em}

During the transient dynamics, if
\begin{align}
\wprod < 0,
& \quad &
\dot{W}_\text{cyc} > 0,
& \quad &
\d_t \G > 0,
\end{align}
during the transient dynamics, driven production enters in competition with the charging of the sub-network because $Ẇ_cyc   -    T Σ̇_$ serves either to store free-energy in the sub-network, or to sustain the production of autocatalytic species (\emph{cf.} Fig.~\ref{fig:efficiency}\textbf{\textsf{b}}).
In this case, the transient efficiency of driven production reads
\begin{equation}
\label{eq:efficiency_DR}
\eta_\text{dp}  =
\dfrac{- \wprod}{\dot{W}_\text{cyc}}
= 	1 - \dfrac{T \dot{\Sigma}_\R + \d_t \G}{\dot{W}_\text{cyc}},
\end{equation}
which converges to its steady form Eq.~\eqref{eq:ss_efficiency_DR} as the networks converges to NESS.
As a consequence, driven production tends to drastically lower the efficiency of energy storage because the amount of work dedicated to increases the energy is reduced.
Importantly, the presence cycles in the sub-network, whose at least one is not also a cycle for the full network, is required for driven replication.
Indeed, if $$ is square or when the cycles of $$  are also cycles for the full network (\emph{i.e.} when $^__2 = 0$), maintaining steady production of autocatalytic species
is necessary costly for the environment: $Ẇ_prod > 0$.
At best, discharge of the sub-network can provide the sufficient energy to drive production when the sub-network is not in a NESS.
This can be maintained if the sub-network never actually settles a NESS but reaches instead an oscillatory regime.
\vspace{2em}

\paragraph{\textbf{Example 7.}}
When $𝔞 = 0$ in Eq.~\eqref{eq:CRN_example}, the sub-network is able to achieve both charging and driven production during the transient dynamics as represented in Figure \ref{fig:efficiency}.
Initially, $Ẇ_cyc - T Σ̇_< - $ (see Fig.~\ref{fig:efficiency}\textbf{\textsf{b}}), as a consequence, discharge of the sub-network is required to sustain production, \emph{i.e.}, $_̣t < 0$ (see Fig.~\ref{fig:efficiency}\textbf{\textsf{a}}).
During the discharge phase, the efficiency of driven production should be modified to take into account the influx of energy coming from the discharge:
\begin{equation}
\label{eq:efficiency_DR_discharge}
\eta_\text{dp}  =
\dfrac{- \wprod}{\dot{W}_\text{cyc} - \d_t \G}
= 	1 - \dfrac{T \dot{\Sigma}_\R}{\dot{W}_\text{cyc} - \d_t \G}.
\end{equation}
This was used to compute the efficiency of driven production during the discharge phase in Fig.~\ref{fig:efficiency}\textbf{\textsf{c}}.
When charging phase starts,  $_̣t >0$, and the charging competes with driven production because $ Ẇ_cyc - T Σ̇_$ is shared between both processes.
The efficiency of energy storage is maximum during the charging phase when a significant portion of $Ẇ_cyc - T Σ̇_$ is converted into free-energy, it then drops to zero as NESS settles in the sub-network and $_̣t = 0$.
On the contrary, the efficiency of driven production increases and attains its maximum when NESS is established as $Ẇ_cyc - T Σ̇_$ serves only to drive production of autocatalytic species that are then released to the environment. \\
If $𝔞 = 1$ in Eq.~\eqref{eq:CRN_example}, $^__2 = 0$ and $Ẇ_cyc = 0$.
Hence, in that case, maintaining steady production of autocatalytic species is necessarily costly, \emph{i.e.} $Ẇ_prod > 0$.

\subsubsection{Driven synthesis}

Often, the power sustaining the cycles of $$ are not sufficient to solely drives the production.
Hence, in that case, the environment needs to inject further energy to sustain the fluxes along the elementary modes:
\begin{align}
\dot{W}_\text{prod} > 0,
& \quad &
\dot{W}_\text{cyc} > 0.
\end{align}
However, one can promote production of species $_1 ⊂$ to the cost of the other controlled autocatalytic species $_2 = - _1$.
The production of species in these subsets comes with a cost
\begin{gather}
\dot{W}_{\text{prod}, ~ \X_i}
\left(
\F_\E \cdot \Nabla^\E_{\R_1} \cdot {(\S_{\R_1})^{-1}}_{\X_i}  + \F_{\X_i}
\right) \cdot \I^{\X_i} \, ,
\end{gather}
such that:
\begin{align}
\dot{W}_{\text{prod}, ~ \X_1} < 0,
& \quad &
\dot{W}_{\text{prod}, ~ \X_2} > 0.
\end{align}
As before, if one wants species $_1$ to increase the additional condition $𝒥^_1 > 0$ is necessary.
This corresponds to \emph{driven synthesis} of species $_1$ by species $_2$ \cite{Penocchio2019}.
Similarly to Eq.~\eqref{eq:efficiency_DR}, the efficiency of this process at steady state is:
\begin{align}
\overline{\eta}_\text{ds} & = \dfrac{- \overline{\dot{W}}_{\text{prod}, ~ \X_1}}{\overline{\dot{W}}_\text{cyc} + \overline{\dot{W}}_{\text{prod}, ~ \X_2} }.% \\[0.5em]
%	& =
%	1 - \dfrac{T \overline{\dot{\Sigma}_\R}}{\overline{\dot{W}}_\text{prod} + \overline{\dot{W}}_{\text{cyc}, ~ \X_2}}.
\end{align}
%As before, as the sub-network enters in its maintenance phase $\eta_\text{ds}$ increases and it is maximum when the NESS is established.

\section{Summary and Discussion}

In this work we developed a thermodynamic theory of autocatalytic networks.
To do so, we relied on the results obtained in the past decades concerning thermodynamic of chemical reaction networks, making use of the topological constraint fulfilled by autocatalytic networks, \emph{i.e.} their absence of conservation law.
Importantly, our results can be easily transposed to a more specific class of autocatalytic network, such as RAF, as illustrated in the example Eq.~\eqref{eq:CRN_example}.

It is crucial to note that the absence of conservation law in autocatalytic network is possible only because the latter is embedded in a larger network having external species, and possibly external reactions.
Indeed, any reasonable (bio)chemical system is expected to fulfill, at the very least, mass conservation principle.
The physical roles of the external species is then twofold: first, they carry the conservation laws, which are all broken once they are all externally controlled.
Secondly, they fuel both the cycles of the autocatalytic network and the production of autocatalytic species.
Crucially, this last process requires a conservative influx of external species, determined solely by the the autocatalytic species.
This dynamic is central to assess the behavior of the autocatalytic network: when it is the only source of external species in the autocatalytic network, the latter is unconditionally detailed balance.
This stands as a major difference with what was previously known to achieve unconditional detailed balance:
in general, it occurs when there is only potential species \cite{Avanzini2022}.
For autocatalytic networks, however, this property is recovered even in the presence of external force species and originates purely from topology.
The latter dictates how the non-conservative forces brought by the external force species are distributed in the full network: if all the cycles of the autocatalytic network are also cycles of the full network, these forces will not be exerted on the autocatalytic reactions.

%Further, this \emph{displacement} of external species is associated to a conservative work, which appears as an additional term in the thermodynamic potential of the autocatalytic network compared to the usual semigrand free energy of open systems.

Further, we emphasized on the fact that production of autocatalytic species can be distinguished from what occurs along the cycles.
This is reflected in the decomposition we used for the elementary fluxes.
Nevertheless, as both processes need the external species, there exists an interplay between them.
The latter can be exploited to drive production of autocatalytic species by the mean of the cycles of the autocatalytic network, provided they are not also cycles for the full network.
Noticing that cycles in the autocatalytic networks is observed only if the system is rectangular, this suggests that, as far as a thermodynamic is concerned, complexity might viewed be beneficial.

The decomposition of the elementary fluxes makes also clear that, in a autonomous networks, sustaining steady production of autocatalytic species $$ is possible only if they are externally controlled.
The type of coupling with the environment will greatly influences how the exterior controls the production.
Specifically, for ideal solutions, chemostatting species $x$, clamp its concentration.
In this case, the macroscopic flux of its production, $𝒥^x$, should be balance by the external flux at any time:
\begin{equation}
\mathcal{J}^x = - I^x.
\end{equation}
However, in this way, it might be difficult to promote production over consumption at steady-state, \emph{i.e.}, $𝒥^x > 0$.
On the contrary, by imposing a degradation following the general form Eq.~\eqref{eq:external_flux_X_species} one can easily control the sign of the macroscopic flux associated to the production at steady-state.
In particular, taking
\begin{align}
\label{eq:condition_J_geq_0}
\phi^x_\text{in} = 0,
& \quad &
\phi^x_\text{out} \propto [x],
\end{align}
in Eq.~\eqref{eq:external_flux_X_species} ensures that species $x$ is effectively produced at steady-state or, equivalently, $^x > 0$.
This can be done experimentally by putting an autocatalytic system in a CSTR with no influx of species $x$.

Even tough we focused on the case where autocatalytic species grow, the converse point of view, consisting in favoring the production of external species \emph{via} an autocatalytic network is also relevant.
Indeed, in the metabolism one observes that some autocatalytic networks are used to synthesize their external species, at the expense of the autocatalytic species.
The most striking example in cellular metabolism is the tricarboxylic (TCA) cycle, also known as Krebs cycle, in the mitochondria:
the latter runs in its recessive sense to produce GTP/ATP and electrons \cite{Lehninger2024}.
Increasing the level of $CO_2$ and introducing metallic ions tend to drive the Krebs cycle in its \emph{productive} sense for the autocatalytic species (in this case it is usually referred to as the reverse Krebs cycle or rTCA), making it the leading candidate for the emergence of autotrophy \cite{Steffens2021, Muchowska2017}.
The formalism presented here can be easily adapted to study the case where the work expanse to sustain production serves to drive the cycles:
\begin{align}
\dot{W}_\text{cyc} < 0,
& \quad &
\wprod > 0.
\end{align}
An efficiency for this process is also easily derived.

\section{Perspectives}

We surmise that the results obtained can be derived from a more systematic approach addressing the thermodynamics of subsystems.
However, an arbitrary subsystem is expected to have conservation laws, which means that its associated stoichiometric matrix $$ (also called \emph{selection matrix} for general systems \cite{Raux2023}) is not full rank.
In that case, the derivation of a physically relevant decomposition of the elementary fluxes as in Eq.~\eqref{eq:decomposition_fluxes} is not as straightforward, however, linear algebra ensures that there will always exists a subset of the \emph{internal} reactions (or transitions) and species (or state) such that the restriction of $$ on these subsets is square non-singular.
From this restriction one will easily extract a conservative dynamics.

Additionally, autocatalysis occurring in a compartment, such as in cells, will usually be associated to volume growth (or shrink).
In this work, we have neglected this effect, assuming the change in volume can be dropped.
In this regard, Sughiyama \emph{et al.} have recently proposed an innovative framework addressing the chemical thermodynamics of growing system in Ref.~\cite{Sughiyama2022}.
However, it was limited to the case where $$ is square (non-singular), we hope that the results derived here will help to extend their results in the general case, leaving the extension of this framework for future studies.
%We leave for future studies the application of this framework to the case of an autocatalytic network in growing compartment.

\renewcommand{\theequation}{A\arabic{equation}}
\setcounter{equation}{0}

\section{\label{sec:Appendix_A} Appendix A: Extension to non-ideal solutions}

In light of Ref.~\cite{Avanzini2021}, the extension to non-ideal mixtures is straightforward.
Indeed, in that case, the chemical potentials of the species become \cite{Kondepudi2014}:
\begin{equation}
\Mu = \Mu^\text{id} + RT \, \ln  \boldsymbol{\gamma}(\e, \z \, ; \, \boldsymbol{p})
\end{equation}
where the ideal chemical potentials, $Μ^id$, are those of Eq.~\eqref{eq:mu_ideal} ; and $γ$ contains the activity coefficients.
If species $s$ behaves ideally:
\begin{equation}
\gamma_s = 1.
\end{equation}
Deviation from ideality stems from interactions with the species (\emph{e.g.} electrostatic interactions in electrolytes, steric interactions...).
As a consequence, the activity coefficients might depend on the concentrations as well as on some parameters, $p$, tuning the strength of interactions.
The Gibbs free-energy of the closed system has now an an additional term:
\begin{equation}
\label{eq:close_G_ni}
G (\e, \z) = G^\text{id} (\e, \z) + G^\text{int} (\e, \z \, ; \, \boldsymbol{p}),
\end{equation}
where $G^id$ is given by Eq.~\eqref{eq:close_G}.
The activity coefficients are related to the interaction free-energy \emph{via},
\begin{equation}
RT \, \ln  \gamma_s
\dfrac{
\partial G^\text{int}
\partial [s]
\end{equation}
Additionally, local detailed balance Eq.~\eqref{eq:LDB} is affected by the presence of interactions:
\begin{equation}
RT \ln \dfrac{k^{+\rho} (\e, \z \, ; \, \boldsymbol{p})}{k^{-\rho}(\e, \z \, ; \, \boldsymbol{p})}
= -
\biggl[
\boldsymbol{\mu}^\circ + RT \ln \boldsymbol{\gamma}(\e, \z \, ; \, \boldsymbol{p})
\biggr]
\cdot \Nabla_\rho,
\end{equation}
and similarly for the additional reactions $r ∈$.
Now, due to interactions, the kinetic rates constants might now depend on the concentrations, as well as on the external parameters.

Further, for thermodynamic consistency, we assume that the free-energy of the closed system Eq.~\eqref{eq:close_G_ni} is lower bounded and increases superlinearly as the concentrations go to infinity.
As shown in the Ref.~\cite{Avanzini2021} (sections  II.F \& II.G), the thermodynamic potential of the autocatalytic sub-network,
\begin{equation*}
\G(\e, \z) = G(\e, \z) - \Mu_{\Ep} \cdot \boldsymbol{m}(\e, \z) + \F_{\E} \cdot \Delta \e (\z),
\end{equation*}
is lower bounded by the equilibrium state Eq.~\eqref{eq:Eq_state_R1}.
In a nutshell, this comes from the fact that  $- Μ_ ·m(, ) + _ ·Δ()$, can go to minus infinity only linearly with respect to the concentrations while $G$ should grow super-linearly.
Unlike before,
\begin{equation}
\Delta \G = \G(\e, \z) - \G(\e, \z_\text{eq}),
\end{equation}
cannot be written solely in term as a relative entropy.
If $G^in$ is a convex function, $Δ$  can be expressed as a relative entropy and the Bregman divergence \cite{Kobayashi2022} of $G^in$:
\begin{equation}
\Delta \G
RT  \, \mathcal{L} \left( \z \Vert \z_\text{eq} \right)
+ RT \, \mathcal{D}_{G^\text{in}} \left( \z \Vert \z_\text{eq} \right).
\end{equation}
Finally, varying the external parameters tuning the interactions, $p$, bring an additional driving work rate in the decomposition of the EPR Eq.~\eqref{eq:second_law_dt_G}:
\begin{equation}
\label{eq:driving_work_int}
\dot{W}_\text{driv}^\text{int}
\partial_{\boldsymbol{p}} G^\text{int} \cdot \d_t \boldsymbol{p} \,.
\end{equation}

\renewcommand{\theequation}{B\arabic{equation}}
\setcounter{equation}{0}

\section*{Appendix B: Catalytic and non-catalytic pathways}

Considering the following elementary reaction $ρ_0 ∈$:
\begin{equation}
\label{eq:non-catalytic_pathway}
\ce{$R_{\rho_0}$ <=> $P_{\rho_0}$},
\end{equation}
where $R_ρ_0$ and $P_ρ_0$ are the reactant and product complexes, respectively.
Note that, species in $R_ρ_0$ or $P_ρ_0$ can be of both types $$ or $$, provided there is, at least, one autocatalytic species in either of the two (so that $ρ_0 ∈$).
Additionally, the following pathway of reactions in the autocatalytic network:
\begin{equation}
\label{eq:catalytic_pathway}
\ce{
$R_{\rho_0}$ + $C_{\rho_1}$ <=>[\rho_1]
$Z_1^{\, \ast}$ <=>[\rho_2]
%		$Z_2^{\, \ast}$ <=>[\rho_3]
$Z_n^{\, \ast}$
$P_{\rho_0}$ + $C_{\rho_1}$
\end{equation}
corresponds the catalyzed counterpart of reaction $ρ_0$ with $C_ρ_1$ being the complex of species catalyzing the production of $P_ρ_0$ from $R_ρ_0$, and the $Z_i^  ∗$s are linear intermediates.
If one finds both Eqs.~\eqref{eq:non-catalytic_pathway} and \eqref{eq:catalytic_pathway} in the sub-network then, it is straightforward to notice that
\begin{equation}
\label{eq:cycle_pathway}
\boldsymbol{c} = \hspace{0.5em}
\begin{pNiceArray}
\boldsymbol{0} & -1 & 1 & \hdots & 1 & 1 & \boldsymbol{0} \\
\cdots & \rho_0 & \rho_1 & \cdots & \rho_n & \rho_{n+1} & \cdots
\end{pNiceArray}^\top,
\end{equation}
will be a cycle for both the sub-network and the full network.
As a consequence, if an autocatalytic network comprises \emph{only} non-catalytic pathways and their catalytic counterpart (this for example the case in the example CRN Eq.~\eqref{eq:CRN_example} when $𝔞 = 1$) then, necessarily, condition Eq.~\eqref{eq:condition_eq} is fulfilled.
In addition, reaction $ρ_0$ can be taken to be the one upon which the cycle Eq.~\eqref{eq:cycle_pathway} disappear once removed.
Equivalently, $ρ_0 ∈_2$ and $ρ_i ∈_1$ for $i ≥1$.

Note that the same results apply if, instead of one, several catalytic complexes act on the reactions $ρ_i$ ($i ≥1$).

\section*{Appendix C: Parameters to generate the figures}

The plots were generated by numerical integration of the the kinetic equation of Eq.~\eqref{eq:CRN_example}.
For all the figures, the energy are given in dimensionless unit where $RT = 1$, concentrations are in molar unit, and time is given in second.

\paragraph{\textbf{Figure \ref{fig:undriven_CRN}}}

For \textbf{\textsf{a}}, \textbf{\textsf{b}}, \textbf{\textsf{c}} and \textbf{\textsf{d}}, the standard chemical potentials of the autocatalytic species were set to
\begin{align*}
\mu_{\sf{Z}_0}^\circ = 1,
& \quad &
\mu^\circ_{\sf{Z}_1} = 5,
& \quad &
\mu^\circ_{\sf{Z}_0^{\, \ast}} = 10,
& \quad &
\mu^\circ_{\sf{Z}_1^{\, \ast}} = 20;
\end{align*}
and the energetic barrier of reaction in $_1$ to
\begin{align*}
\Delta G^ \circ_\sf{1} = -1,
&  &
\Delta G^ \circ_\sf{2} = 5,
&  &
\Delta G^ \circ_\sf{3} = -0.1,
&  &
\Delta G^ \circ_\sf{4} = 5.
\end{align*}
The standard chemical potential of the external species and the remaining energy barrier are then completely determined from 1) $ΔG^∘_i = - Μ^∘·_i$ for any reaction $i$ and 2) the Wegscheider condition Eq.~\eqref{eq:Wegscheider}.
The kinetic rates are inferred from the energetic barrier with the local detailed balance condition Eq.~\eqref{eq:LDB} with
\begin{align*}
k^{+ \rho} = 1,
& \quad &
\nu^{+ r} = 1,
\end{align*}
for all $ρ∈$ and $r ∈$.
Finally, we impose the following chemostatted concentrations for the external species:
\begin{align*}
[\sf{F}_0] = 10^{-1} \, \mathrm{M},
& \quad &
[\sf{F}_1] = 10^{-2} \, \mathrm{M},
\\[0.5em]
[\sf{W}_0] = 10^{-1} \, \mathrm{M},
& \quad &
[\sf{W}_1] = 10^{-2} \, \mathrm{M}.
\end{align*}

\paragraph{\textbf{Figure \ref{fig:efficiency}}}
For \textbf{\textsf{a}}, \textbf{\textsf{b}} and  \textbf{\textsf{c}}, the standard chemical potentials of the external species were set to
\begin{align*}
\mu_{\sf{F}_0}^\circ = 5,
& \quad &
\mu^\circ_{\sf{F}_1} = 10,
& \quad &
\mu^\circ_{\sf{W}_0} = 10,
& \quad &
\mu^\circ_{\sf{W}_1} = 15;
\end{align*}
and the energetic barrier of reaction in $_1$ to
\begin{align*}
\Delta G^ \circ_\sf{1} = -0.1,
&  &
\Delta G^ \circ_\sf{2} = 2,
&  &
\Delta G^ \circ_\sf{3} = -1,
&  &
\Delta G^ \circ_\sf{4} = 5.
\end{align*}
Again, the remaining standard chemical potentials and the remaining energetic barriers are found with those specified using the constraint.
In particular this high value for the standard chemical potential of the external species favored low energetic autocatalytic species which tend to favor driven production:
\begin{align*}
\mu_{\sf{Z}_0}^\circ = -1,
& \quad &
\mu^\circ_{\sf{Z}_1} \approx -3.1,
& \quad &
\mu^\circ_{\sf{Z}_0^{\, \ast}} \approx 3.9,
& \quad &
\mu^\circ_{\sf{Z}_1^{\, \ast}} = 5.9.
\end{align*}
The chemostats on the external species are,
\begin{align*}
[\sf{W}_0] = 5\cdot 10^{-2} \, \mathrm{M},
& \quad &
[\sf{W}_1] = 10^{-2} \, \mathrm{M},
\\[0.5em]
[\sf{F}_0] = e^\alpha [\sf{W}_0]  ,
& \quad &
[\sf{F}_1] = e^\beta [\sf{W}_1],
\end{align*}
with $α= 3$ and $β= -4$.
Further, the two chemostats on the autocatalytic species read:
\begin{align*}
[\sf{Z}_0] = 5 \cdot 10^{-2} \, \mathrm{M}
& \quad &
[\sf{Z}_1] = 10^{-2} \, \mathrm{M}.
\end{align*}

\paragraph{\textbf{Figure \ref{fig:DR_ss}}}
The parameters of Fig.~\ref{fig:DR_ss} are identical to those of Fig.~\ref{fig:efficiency}, except for the value of $α$ and $β$ for the external chemostats that are now varying.
In that way, the force associated to the waste species are:
\begin{align*}
\F_{\sf{W}_0} = \mu^\circ_{\sf{W}_0} - \mu^\circ_{\sf{F}_0} - \alpha,
& \quad &
\F_{\sf{W}_1} = \mu^\circ_{\sf{W}_1} - \mu^\circ_{\sf{F}_1} - \beta.
\end{align*}

\section*{Acknowledgment}

The author acknowledges fruitful discussions with Paul Raux and David Lacoste and further thanks Yann Sakref and Luis Dinis for their suggestions on an earlier draft.

\begin{thebibliography}{51}%
\makeatletter
\providecommand \@ifxundefined [1]{%
\@ifx{#1\undefined}
\providecommand \@ifnum [1]{%
\ifnum #1\expandafter \@firstoftwo
\else \expandafter \@secondoftwo
\fi
\providecommand \@ifx [1]{%
\ifx #1\expandafter \@firstoftwo
\else \expandafter \@secondoftwo
\fi
\providecommand \natexlab [1]{#1}%
\providecommand \enquote  [1]{``#1''}%
\providecommand \bibnamefont  [1]{#1}%
\providecommand \bibfnamefont [1]{#1}%
\providecommand \citenamefont [1]{#1}%
\providecommand \href@noop [0]{\@secondoftwo}%
\providecommand \href [0]{\begingroup \@sanitize@url \@href}%
\providecommand \@href[1]{\@@startlink{#1}\@@href}%
\providecommand \@@href[1]{\endgroup#1\@@endlink}%
\providecommand \@sanitize@url [0]{\catcode `\\12\catcode `\$12\catcode
`\&12\catcode `\#12\catcode `\^12\catcode `\_12\catcode `\%12\relax}%
\providecommand \@@startlink[1]{}%
\providecommand \@@endlink[0]{}%
\providecommand \url  [0]{\begingroup\@sanitize@url \@url }%
\providecommand \@url [1]{\endgroup\@href {#1}{\urlprefix }}%
\providecommand \urlprefix  [0]{URL }%
\providecommand \Eprint [0]{\href }%
\providecommand \doibase [0]{https://doi.org/}%
\providecommand \selectlanguage [0]{\@gobble}%
\providecommand \bibinfo  [0]{\@secondoftwo}%
\providecommand \bibfield  [0]{\@secondoftwo}%
\providecommand \translation [1]{[#1]}%
\providecommand \BibitemOpen [0]{}%
\providecommand \bibitemStop [0]{}%
\providecommand \bibitemNoStop [0]{.\EOS\space}%
\providecommand \EOS [0]{\spacefactor3000\relax}%
\providecommand \BibitemShut  [1]{\csname bibitem#1\endcsname}%
\let\auto@bib@innerbib\@empty
\bibitem [{\citenamefont {Schuster}(2019)}]{Schuster2019}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.}~\bibnamefont
{Schuster}},\ }\href {https://doi.org/10.1007/s00706-019-02437-z} {\bibfield
{journal} {\bibinfo  {journal} {Monatshefte für Chemie - Chemical Monthly}\
}\textbf {\bibinfo {volume} {150}},\ \bibinfo {pages} {763} (\bibinfo {year}
{2019})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Roy}\ \emph {et~al.}(2021)\citenamefont {Roy},
\citenamefont {Goberman},\ and\ \citenamefont {Pugatch}}]{Roy2021}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
{Roy}}, \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Goberman}},\ and\
\bibinfo {author} {\bibfnamefont {R.}~\bibnamefont {Pugatch}},\ }\href
{https://pnas.org/doi/full/10.1073/pnas.2107829118} {\bibfield  {journal}
{\bibinfo  {journal} {Proceedings of the National Academy of Sciences}\
}\textbf {\bibinfo {volume} {118}} (\bibinfo {year} {2021})}\BibitemShut
\bibitem [{\citenamefont {Lin}\ \emph {et~al.}(2020)\citenamefont {Lin},
\citenamefont {Kussell}, \citenamefont {Young},\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {W.-H.}\ \bibnamefont
{Lin}}, \bibinfo {author} {\bibfnamefont {E.}~\bibnamefont {Kussell}},
\bibinfo {author} {\bibfnamefont {L.-S.}\ \bibnamefont {Young}},\ and\
\bibinfo {author} {\bibfnamefont {C.}~\bibnamefont {Jacobs-Wagner}},\ }\href
{https://doi.org/10.1073/pnas.2013061117} {\bibfield  {journal} {\bibinfo
{journal} {Proceedings of the National Academy of Sciences}\ }\textbf
{\bibinfo {volume} {117}},\ \bibinfo {pages} {27795} (\bibinfo {year}
{2020})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Sakref}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont
{Sakref}}\ and\ \bibinfo {author} {\bibfnamefont {O.}~\bibnamefont
{Rivoire}},\ }\href {https://doi.org/10.1016/j.jtbi.2023.111714} {\bibfield
{journal} {\bibinfo  {journal} {Journal of Theoretical Biology}\ }\textbf
{\bibinfo {volume} {579}},\ \bibinfo {pages} {111714} (\bibinfo {year}
{2024})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Ameta}\ \emph {et~al.}(2021)\citenamefont {Ameta},
\citenamefont {Matsubara}, \citenamefont {Chakraborty}, \citenamefont
{Krishna},\ and\ \citenamefont {Thutupalli}}]{Ameta2021}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
{Ameta}}, \bibinfo {author} {\bibfnamefont {Y.~J.}\ \bibnamefont
{Matsubara}}, \bibinfo {author} {\bibfnamefont {N.}~\bibnamefont
{Chakraborty}}, \bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
{Krishna}},\ and\ \bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
{Thutupalli}},\ }\href {https://doi.org/10.3390/life11040308} {\bibfield
{journal} {\bibinfo  {journal} {Life}\ }\textbf {\bibinfo {volume} {11}},\
\bibinfo {pages} {308} (\bibinfo {year} {2021})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Peng}\ \emph {et~al.}(2022)\citenamefont {Peng},
\citenamefont {Linderoth},\ and\ \citenamefont {Baum}}]{Peng2022}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {Z.}~\bibnamefont
{Peng}}, \bibinfo {author} {\bibfnamefont {J.}~\bibnamefont {Linderoth}},\
and\ \bibinfo {author} {\bibfnamefont {D.~A.}\ \bibnamefont {Baum}},\ }\href
{https://doi.org/10.1371/journal.pcbi.1010498} {\bibfield  {journal}
{\bibinfo  {journal} {PLOS Computational Biology}\ }\textbf {\bibinfo
{volume} {18}},\ \bibinfo {pages} {e1010498} (\bibinfo {year}
{2022})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Hordijk}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {W.}~\bibnamefont
{Hordijk}}\ and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
{Steel}},\ }\href {https://doi.org/10.1007/s11084-014-9374-5} {\bibfield
{journal} {\bibinfo  {journal} {Origins of Life and Evolution of Biospheres}\
}\textbf {\bibinfo {volume} {44}},\ \bibinfo {pages} {111} (\bibinfo {year}
{2014})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Vincent}\ \emph {et~al.}(2021)\citenamefont
{Vincent}, \citenamefont {Colón-Santos}, \citenamefont {Cleaves},
\citenamefont {Baum},\ and\ \citenamefont {Maurer}}]{Vincent2021}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {L.}~\bibnamefont
{Vincent}}, \bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
{Colón-Santos}}, \bibinfo {author} {\bibfnamefont {H.~J.}\ \bibnamefont
{Cleaves}}, \bibinfo {author} {\bibfnamefont {D.~A.}\ \bibnamefont {Baum}},\
and\ \bibinfo {author} {\bibfnamefont {S.~E.}\ \bibnamefont {Maurer}},\
}\href {https://doi.org/10.3390/life11111221} {\bibfield  {journal} {\bibinfo
{journal} {Life}\ }\textbf {\bibinfo {volume} {11}},\ \bibinfo {pages}
{1221} (\bibinfo {year} {2021})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Peng}\ \emph {et~al.}(2020)\citenamefont {Peng},
\citenamefont {Plum}, \citenamefont {Gagrani},\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {Z.}~\bibnamefont
{Peng}}, \bibinfo {author} {\bibfnamefont {A.~M.}\ \bibnamefont {Plum}},
\bibinfo {author} {\bibfnamefont {P.}~\bibnamefont {Gagrani}},\ and\ \bibinfo
{author} {\bibfnamefont {D.~A.}\ \bibnamefont {Baum}},\ }\href
{https://doi.org/10.1016/j.jtbi.2020.110451} {\bibfield  {journal} {\bibinfo
{journal} {Journal of Theoretical Biology}\ }\textbf {\bibinfo {volume}
{507}},\ \bibinfo {pages} {110451} (\bibinfo {year} {2020})}\BibitemShut
\bibitem [{\citenamefont {Xavier}\ \emph {et~al.}(2020)\citenamefont {Xavier},
\citenamefont {Hordijk}, \citenamefont {Kauffman}, \citenamefont {Steel},\
and\ \citenamefont {Martin}}]{Xavier2020}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~C.}\ \bibnamefont
{Xavier}}, \bibinfo {author} {\bibfnamefont {W.}~\bibnamefont {Hordijk}},
\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont {Kauffman}}, \bibinfo
{author} {\bibfnamefont {M.}~\bibnamefont {Steel}},\ and\ \bibinfo {author}
{\bibfnamefont {W.~F.}\ \bibnamefont {Martin}},\ }\href
{https://doi.org/10.1098/rspb.2019.2377} {\bibfield  {journal} {\bibinfo
{journal} {Proceedings of the Royal Society B: Biological Sciences}\ }\textbf
{\bibinfo {volume} {287}},\ \bibinfo {pages} {20192377} (\bibinfo {year}
{2020})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Higgs}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.~G.}\ \bibnamefont
{Higgs}}\ and\ \bibinfo {author} {\bibfnamefont {N.}~\bibnamefont {Lehman}},\
}\href {https://doi.org/10.1038/nrg3841} {\bibfield  {journal} {\bibinfo
{journal} {Nature Reviews Genetics}\ }\textbf {\bibinfo {volume} {16}},\
\bibinfo {pages} {7} (\bibinfo {year} {2015})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Bernhardt}(2012)}]{Bernhardt2012}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {H.~S.}\ \bibnamefont
{Bernhardt}},\ }\href {https://doi.org/10.1186/1745-6150-7-23} {\bibfield
{journal} {\bibinfo  {journal} {Biology Direct}\ }\textbf {\bibinfo {volume}
{7}},\ \bibinfo {pages} {23} (\bibinfo {year} {2012})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Pavlinova}\ \emph {et~al.}(2023)\citenamefont
{Pavlinova}, \citenamefont {Lambert}, \citenamefont {Malaterre},\ and\
\citenamefont {Nghe}}]{Pavlinova2023}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.}~\bibnamefont
{Pavlinova}}, \bibinfo {author} {\bibfnamefont {C.~N.}\ \bibnamefont
{Lambert}}, \bibinfo {author} {\bibfnamefont {C.}~\bibnamefont {Malaterre}},\
and\ \bibinfo {author} {\bibfnamefont {P.}~\bibnamefont {Nghe}},\ }\href
{https://doi.org/10.1002/1873-3468.14507} {\bibfield  {journal} {\bibinfo
{journal} {FEBS Letters}\ }\textbf {\bibinfo {volume} {597}},\ \bibinfo
{pages} {344} (\bibinfo {year} {2023})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Kauffman}(1971)}]{Kauffman1971}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.~A.}\ \bibnamefont
{Kauffman}},\ }\href {https://doi.org/10.1080/01969727108545830} {\bibfield
{journal} {\bibinfo  {journal} {Journal of Cybernetics}\ }\textbf {\bibinfo
{volume} {1}},\ \bibinfo {pages} {71} (\bibinfo {year} {1971})}\BibitemShut
\bibitem [{\citenamefont {Kauffman}(1986)}]{Kauffman1986}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.~A.}\ \bibnamefont
{Kauffman}},\ }\href {https://doi.org/10.1016/S0022-5193(86)80047-9}
{\bibfield  {journal} {\bibinfo  {journal} {Journal of Theoretical Biology}\
}\textbf {\bibinfo {volume} {119}},\ \bibinfo {pages} {1} (\bibinfo {year}
{1986})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Feinberg}(2019)}]{Feinberg2019}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
{Feinberg}},\ }\href {https://doi.org/10.1007/978-3-030-03858-8} {\emph
{\bibinfo {title} {Foundations of Chemical Reaction Network Theory}}}\
(\bibinfo  {publisher} {Springer International Publishing},\ \bibinfo {year}
{2019})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Blokhuis}\ \emph {et~al.}(2020)\citenamefont
{Blokhuis}, \citenamefont {Lacoste},\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
{Blokhuis}}, \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Lacoste}},\
and\ \bibinfo {author} {\bibfnamefont {P.}~\bibnamefont {Nghe}},\ }\href
{https://doi.org/10.1073/pnas.2013527117} {\bibfield  {journal} {\bibinfo
{journal} {Proceedings of the National Academy of Sciences}\ }\textbf
{\bibinfo {volume} {117}},\ \bibinfo {pages} {25230} (\bibinfo {year}
{2020})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Andersen}\ \emph {et~al.}(2021)\citenamefont
{Andersen}, \citenamefont {Flamm}, \citenamefont {Merkle},\ and\
\citenamefont {Stadler}}]{Andersen2021}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~L.}\ \bibnamefont
{Andersen}}, \bibinfo {author} {\bibfnamefont {C.}~\bibnamefont {Flamm}},
\bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Merkle}},\ and\ \bibinfo
{author} {\bibfnamefont {P.~F.}\ \bibnamefont {Stadler}},\ }\href
{https://doi.org/10.48550/arXiv.2107.03086} {\bibfield  {journal} {\bibinfo
{journal} {arXiv [Preprint]}\ } (\bibinfo {year} {2021})}\BibitemShut
\bibitem [{\citenamefont {Gibbs}(1928)}]{Gibbs1928}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~W.}\ \bibnamefont
{Gibbs}},\ }\href
{\emph {\bibinfo {title} {The Collected Works of J. Willard Gibbs}}}\
(\bibinfo  {publisher} {Longmans, Green and Co.},\ \bibinfo {year}
{1928})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {de~Donder}(1927)}]{deDonder1927}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.}~\bibnamefont
{de~Donder}},\ }\href {https://books.google.fr/books?id=MQ7nAAAAMAAJ} {\emph
{\bibinfo {title} {L'affinit{\'e}}}},\ \bibinfo {series} {Memoires de la
Classe des sciences: Collection in-8o}\ No.\ \bibinfo {number} {vol.~1}\
(\bibinfo  {publisher} {Gauthier-Villars},\ \bibinfo {year}
{1927})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {McQuarrie}(1967)}]{McQuarrie1967}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~A.}\ \bibnamefont
{McQuarrie}},\ }\href {https://doi.org/10.2307/3212214} {\bibfield  {journal}
{\bibinfo  {journal} {Journal of Applied Probability}\ }\textbf {\bibinfo
{volume} {4}},\ \bibinfo {pages} {413} (\bibinfo {year} {1967})}\BibitemShut
\bibitem [{\citenamefont {Gillespie}(1992)}]{Gillespie1992}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~T.}\ \bibnamefont
{Gillespie}},\ }\href {https://doi.org/10.1016/0378-4371(92)90283-V}
{\bibfield  {journal} {\bibinfo  {journal} {Physica A: Statistical Mechanics
and its Applications}\ }\textbf {\bibinfo {volume} {188}},\ \bibinfo {pages}
{404} (\bibinfo {year} {1992})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Schnakenberg}(1976)}]{Schnakenberg1976}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
{Schnakenberg}},\ }\href {https://doi.org/10.1103/RevModPhys.48.571}
{\bibfield  {journal} {\bibinfo  {journal} {Reviews of Modern Physics}\
}\textbf {\bibinfo {volume} {48}},\ \bibinfo {pages} {571} (\bibinfo {year}
{1976})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Mou}\ \emph {et~al.}(1986)\citenamefont {Mou},
\citenamefont {li~Luo},\ and\ \citenamefont {Nicolis}}]{Mou1986}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {C.~Y.}\ \bibnamefont
{Mou}}, \bibinfo {author} {\bibfnamefont {J.}~\bibnamefont {li~Luo}},\ and\
\bibinfo {author} {\bibfnamefont {G.}~\bibnamefont {Nicolis}},\ }\href
{https://doi.org/10.1063/1.450623} {\bibfield  {journal} {\bibinfo  {journal}
{The Journal of Chemical Physics}\ }\textbf {\bibinfo {volume} {84}},\
\bibinfo {pages} {7011} (\bibinfo {year} {1986})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Hill}(1983)}]{Hill1983}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.~L.}\ \bibnamefont
{Hill}},\ }\href {https://doi.org/10.1073/pnas.80.10.2922} {\bibfield
{journal} {\bibinfo  {journal} {Proceedings of the National Academy of
Sciences}\ }\textbf {\bibinfo {volume} {80}},\ \bibinfo {pages} {2922}
(\bibinfo {year} {1983})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Rao}\ and\ \citenamefont {Esposito}(2016)}]{Rao2016}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {R.}~\bibnamefont
{Rao}}\ and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Esposito}},\
}\href {https://doi.org/10.1103/PhysRevX.6.041064} {\bibfield  {journal}
{\bibinfo  {journal} {Physical Review X}\ }\textbf {\bibinfo {volume} {6}}
(\bibinfo {year} {2016})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Ge}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {H.}~\bibnamefont
{Ge}}\ and\ \bibinfo {author} {\bibfnamefont {H.}~\bibnamefont {Qian}},\
}\href {https://doi.org/10.1016/j.chemphys.2016.03.026} {\bibfield  {journal}
{\bibinfo  {journal} {Chemical Physics}\ }\textbf {\bibinfo {volume} {472}},\
\bibinfo {pages} {241} (\bibinfo {year} {2016}{\natexlab{a}})}\BibitemShut
\bibitem [{\citenamefont {Avanzini}\ \emph {et~al.}(2021)\citenamefont
{Avanzini}, \citenamefont {Penocchio}, \citenamefont {Falasco},\ and\
\citenamefont {Esposito}}]{Avanzini2021}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {F.}~\bibnamefont
{Avanzini}}, \bibinfo {author} {\bibfnamefont {E.}~\bibnamefont {Penocchio}},
\bibinfo {author} {\bibfnamefont {G.}~\bibnamefont {Falasco}},\ and\ \bibinfo
{author} {\bibfnamefont {M.}~\bibnamefont {Esposito}},\ }\href
{\bibfield  {journal} {\bibinfo  {journal} {The Journal of Chemical Physics}\
}\textbf {\bibinfo {volume} {154}},\ \bibinfo {pages} {94114} (\bibinfo
{year} {2021})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Alberty}(2003)}]{Alberty2003}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {R.~A.}\ \bibnamefont
{Alberty}},\ }\href {https://doi.org/10.1002/0471332607} {\emph {\bibinfo
{title} {Thermodynamics of Biochemical Reactions}}}\ (\bibinfo  {publisher}
{Wiley},\ \bibinfo {year} {2003})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Qian}\ and\ \citenamefont {Beard}(2005)}]{Qian2005}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {H.}~\bibnamefont
{Qian}}\ and\ \bibinfo {author} {\bibfnamefont {D.~A.}\ \bibnamefont
{Beard}},\ }\href {https://doi.org/10.1016/j.bpc.2004.12.001} {\bibfield
{journal} {\bibinfo  {journal} {Biophysical Chemistry}\ }\textbf {\bibinfo
{volume} {114}},\ \bibinfo {pages} {213} (\bibinfo {year}
{2005})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Cengio}\ \emph {et~al.}(2023)\citenamefont {Cengio},
\citenamefont {Lecomte},\ and\ \citenamefont {Polettini}}]{Cengio2022}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.~D.}\ \bibnamefont
{Cengio}}, \bibinfo {author} {\bibfnamefont {V.}~\bibnamefont {Lecomte}},\
and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Polettini}},\ }\href
{https://doi.org/10.1103/PhysRevX.13.021040} {\bibfield  {journal} {\bibinfo
{journal} {Physical Review X}\ }\textbf {\bibinfo {volume} {13}},\ \bibinfo
{pages} {021040} (\bibinfo {year} {2023})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Harunari}\ \emph {et~al.}(2022)\citenamefont
{Harunari}, \citenamefont {Dutta}, \citenamefont {Polettini},\ and\
\citenamefont {Édgar Roldán}}]{Harunari2022}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.~E.}\ \bibnamefont
{Harunari}}, \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont {Dutta}},
\bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Polettini}},\ and\
\bibinfo {author} {\bibnamefont {Édgar Roldán}},\ }\href
{https://doi.org/10.1103/PhysRevX.12.041026} {\bibfield  {journal} {\bibinfo
{journal} {Physical Review X}\ }\textbf {\bibinfo {volume} {12}},\ \bibinfo
{pages} {041026} (\bibinfo {year} {2022})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Despons}\ \emph {et~al.}(2023)\citenamefont
{Despons}, \citenamefont {de~Decker},\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
{Despons}}, \bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont {de~Decker}},\
and\ \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Lacoste}},\ }\href
{http://arxiv.org/abs/2306.02366} {\bibfield  {journal} {\bibinfo  {journal}
{arXiv [Preprint]}\ } (\bibinfo {year} {2023})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Ge}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {H.}~\bibnamefont
{Ge}}\ and\ \bibinfo {author} {\bibfnamefont {H.}~\bibnamefont {Qian}},\
}\href {https://doi.org/10.1103/PhysRevE.94.052150} {\bibfield  {journal}
{\bibinfo  {journal} {Physical Review E}\ }\textbf {\bibinfo {volume} {94}},\
\bibinfo {pages} {3} (\bibinfo {year} {2016}{\natexlab{b}})}\BibitemShut
\bibitem [{\citenamefont {Hordijk}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {W.}~\bibnamefont
{Hordijk}}\ and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
{Steel}},\ }\href {https://doi.org/10.1016/j.jtbi.2003.11.020} {\bibfield
{journal} {\bibinfo  {journal} {Journal of Theoretical Biology}\ }\textbf
{\bibinfo {volume} {227}},\ \bibinfo {pages} {451} (\bibinfo {year}
{2004})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Wachtel}\ \emph {et~al.}(2018)\citenamefont
{Wachtel}, \citenamefont {Rao},\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
{Wachtel}}, \bibinfo {author} {\bibfnamefont {R.}~\bibnamefont {Rao}},\ and\
\bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Esposito}},\ }\href
{https://doi.org/10.1088/1367-2630/aab5c9} {\bibfield  {journal} {\bibinfo
{journal} {New Journal of Physics}\ }\textbf {\bibinfo {volume} {20}},\
\bibinfo {pages} {042002} (\bibinfo {year} {2018})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Avanzini}\ \emph {et~al.}(2020)\citenamefont
{Avanzini}, \citenamefont {Falasco},\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {F.}~\bibnamefont
{Avanzini}}, \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont {Falasco}},\
and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Esposito}},\ }\href
{https://doi.org/10.1088/1367-2630/abafea} {\bibfield  {journal} {\bibinfo
{journal} {New Journal of Physics}\ }\textbf {\bibinfo {volume} {22}},\
\bibinfo {pages} {093040} (\bibinfo {year} {2020})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Rao}\ and\ \citenamefont {Esposito}(2018)}]{Rao2018}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {R.}~\bibnamefont
{Rao}}\ and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Esposito}},\
}\href {https://doi.org/10.1063/1.5042253} {\bibfield  {journal} {\bibinfo
{journal} {Journal of Chemical Physics}\ }\textbf {\bibinfo {volume} {149}},\
\bibinfo {pages} {245101} (\bibinfo {year} {2018})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Avanzini}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {F.}~\bibnamefont
{Avanzini}}\ and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
{Esposito}},\ }\href
{\bibfield  {journal} {\bibinfo  {journal} {The Journal of Chemical Physics}\
}\textbf {\bibinfo {volume} {156}} (\bibinfo {year} {2022})}\BibitemShut
\bibitem [{\citenamefont {Schuster}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
{Schuster}}\ and\ \bibinfo {author} {\bibfnamefont {R.}~\bibnamefont
{Schuster}},\ }\href {https://doi.org/10.1007/BF01171883} {\bibfield
{journal} {\bibinfo  {journal} {Journal of Mathematical Chemistry}\ }\textbf
{\bibinfo {volume} {3}},\ \bibinfo {pages} {25} (\bibinfo {year}
{1989})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Polettini}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
{Polettini}}\ and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
{Esposito}},\ }\href
{\bibfield  {journal} {\bibinfo  {journal} {The Journal of Chemical Physics}\
}\textbf {\bibinfo {volume} {141}} (\bibinfo {year} {2014})}\BibitemShut
\bibitem [{\citenamefont {Hill}(1977)}]{Hill1977}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.~L.}\ \bibnamefont
{Hill}},\ }\href
{\emph {\bibinfo {title} {Free energy transduction in biology : the
steady-state kinetic and thermodynamic formalism}}}\ (\bibinfo  {publisher}
{Academic Press},\ \bibinfo {year} {1977})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Ragazzon}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
{Ragazzon}}\ and\ \bibinfo {author} {\bibfnamefont {L.~J.}\ \bibnamefont
{Prins}},\ }\href {https://doi.org/10.1038/s41565-018-0250-8} {\bibfield
{journal} {\bibinfo  {journal} {Nature Nanotechnology}\ }\textbf {\bibinfo
{volume} {13}},\ \bibinfo {pages} {882} (\bibinfo {year} {2018})}\BibitemShut
\bibitem [{\citenamefont {Penocchio}\ \emph {et~al.}(2019)\citenamefont
{Penocchio}, \citenamefont {Rao},\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {E.}~\bibnamefont
{Penocchio}}, \bibinfo {author} {\bibfnamefont {R.}~\bibnamefont {Rao}},\
and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Esposito}},\ }\href
{https://doi.org/10.1038/s41467-019-11676-x} {\bibfield  {journal} {\bibinfo
{journal} {Nature Communications}\ }\textbf {\bibinfo {volume} {10}},\
\bibinfo {pages} {3865} (\bibinfo {year} {2019})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Nelson}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~L.}\ \bibnamefont
{Nelson}}\ and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Cox}},\
}\href@noop {} {\emph {\bibinfo {title} {Principles of Biochemistry}}}\
(\bibinfo  {publisher} {W.H Freeman and Co.},\ \bibinfo {year} {2024})\
Chap.~\bibinfo {chapter} {16}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Steffens}\ \emph {et~al.}(2021)\citenamefont
{Steffens}, \citenamefont {Pettinato}, \citenamefont {Steiner}, \citenamefont
{Mall}, \citenamefont {König}, \citenamefont {Eisenreich},\ and\
\citenamefont {Berg}}]{Steffens2021}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {L.}~\bibnamefont
{Steffens}}, \bibinfo {author} {\bibfnamefont {E.}~\bibnamefont {Pettinato}},
\bibinfo {author} {\bibfnamefont {T.~M.}\ \bibnamefont {Steiner}}, \bibinfo
{author} {\bibfnamefont {A.}~\bibnamefont {Mall}}, \bibinfo {author}
{\bibfnamefont {S.}~\bibnamefont {König}}, \bibinfo {author} {\bibfnamefont
{W.}~\bibnamefont {Eisenreich}},\ and\ \bibinfo {author} {\bibfnamefont
{I.~A.}\ \bibnamefont {Berg}},\ }\href
{https://doi.org/10.1038/s41586-021-03456-9} {\bibfield  {journal} {\bibinfo
{journal} {Nature}\ }\textbf {\bibinfo {volume} {592}},\ \bibinfo {pages}
{784} (\bibinfo {year} {2021})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Muchowska}\ \emph {et~al.}(2017)\citenamefont
{Muchowska}, \citenamefont {Varma}, \citenamefont {Chevallot-Beroux},
\citenamefont {Lethuillier-Karl}, \citenamefont {Li},\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {K.~B.}\ \bibnamefont
{Muchowska}}, \bibinfo {author} {\bibfnamefont {S.~J.}\ \bibnamefont
{Varma}}, \bibinfo {author} {\bibfnamefont {E.}~\bibnamefont
{Chevallot-Beroux}}, \bibinfo {author} {\bibfnamefont {L.}~\bibnamefont
{Lethuillier-Karl}}, \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
{Li}},\ and\ \bibinfo {author} {\bibfnamefont {J.}~\bibnamefont {Moran}},\
}\href {https://doi.org/10.1038/s41559-017-0311-7} {\bibfield  {journal}
{\bibinfo  {journal} {Nature Ecology \& Evolution}\ }\textbf {\bibinfo
{volume} {1}},\ \bibinfo {pages} {1716} (\bibinfo {year} {2017})}\BibitemShut
\bibitem [{\citenamefont {Raux}\ \emph {et~al.}(2023)\citenamefont {Raux},
\citenamefont {Goupil},\ and\ \citenamefont {Verley}}]{Raux2023}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.}~\bibnamefont
{Raux}}, \bibinfo {author} {\bibfnamefont {C.}~\bibnamefont {Goupil}},\ and\
\bibinfo {author} {\bibfnamefont {G.}~\bibnamefont {Verley}},\ }\href
{https://doi.org/10.48550/arXiv.2309.12922} {\bibfield  {journal} {\bibinfo
{journal} {arXiv [Preprint]}\ } (\bibinfo {year} {2023})}\BibitemShut
\bibitem [{\citenamefont {Sughiyama}\ \emph {et~al.}(2022)\citenamefont
{Sughiyama}, \citenamefont {Kamimura}, \citenamefont {Loutchko},\ and\
\citenamefont {Kobayashi}}]{Sughiyama2022}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont
{Sughiyama}}, \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont {Kamimura}},
\bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Loutchko}},\ and\
\bibinfo {author} {\bibfnamefont {T.~J.}\ \bibnamefont {Kobayashi}},\ }\href
{https://doi.org/10.1103/PhysRevResearch.4.033191} {\bibfield  {journal}
{\bibinfo  {journal} {Physical Review Research}\ }\textbf {\bibinfo {volume}
{4}},\ \bibinfo {pages} {033191} (\bibinfo {year} {2022})}\BibitemShut
\bibitem [{\citenamefont {Kondepudi}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
{Kondepudi}}\ and\ \bibinfo {author} {\bibfnamefont {I.}~\bibnamefont
{Prigogine}},\ }\bibinfo {title} {Modern thermodynamics}\ (\bibinfo
{publisher} {Wiley},\ \bibinfo {year} {2014})\ Chap.~\bibinfo {chapter}
{8}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Kobayashi}\ \emph {et~al.}(2022)\citenamefont
{Kobayashi}, \citenamefont {Loutchko}, \citenamefont {Kamimura},\ and\
\citenamefont {Sughiyama}}]{Kobayashi2022}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.~J.}\ \bibnamefont
{Kobayashi}}, \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Loutchko}},
\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont {Kamimura}},\ and\
\bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont {Sughiyama}},\ }\href
{https://doi.org/10.1103/PhysRevResearch.4.033208} {\bibfield  {journal}
{\bibinfo  {journal} {Physical Review Research}\ }\textbf {\bibinfo {volume}
{4}},\ \bibinfo {pages} {033208} (\bibinfo {year} {2022})}\BibitemShut
\end{thebibliography}%

\end{document}