# Nonequilibrium Thermodynamics of Non-Ideal Chemical Reaction Networks

Francesco Avanzini<EMAIL_ADDRESS>Complex Systems and Statistical
Mechanics, Department of Physics and Materials Science, University of
Luxembourg, L-1511 Luxembourg  Emanuele Penocchio<EMAIL_ADDRESS>Complex Systems and Statistical Mechanics, Department of Physics and Materials
Science, University of Luxembourg, L-1511 Luxembourg  Gianmaria Falasco
<EMAIL_ADDRESS>Complex Systems and Statistical Mechanics, Department
of Physics and Materials Science, University of Luxembourg, L-1511 Luxembourg
Massimiliano Esposito<EMAIL_ADDRESS>Complex Systems and
Statistical Mechanics, Department of Physics and Materials Science, University
of Luxembourg, L-1511 Luxembourg

###### Abstract

All current formulations of nonequilibrium thermodynamics of open chemical
reaction networks rely on the assumption of non-interacting species. We
develop a general theory which accounts for interactions between chemical
species within a mean-field approach using activity coefficients.
Thermodynamic consistency requires that rate equations do not obey to standard
mass-action kinetics, but account for the interactions with concentration
dependent kinetic constants. Many features of the ideal formulations are
recovered. Crucially, the thermodynamic potential and the forces driving non-
ideal chemical systems out of equilibrium are identified. Our theory is
general and holds for any mean-field expression of the interactions leading to
lower bounded free energies.

##  I Introduction

Thermodynamics studies the interconversions of energy. It was originally
formulated as an equilibrium theory in the 19th century. Phenomenological
extensions close to equilibrium, in the so-called linear regime, were
introduced in the first half of the 20th century Prigogine (1961); de Groot
and Mazur (1984). The earliest formulations beyond the linear regime were
developed for stochastic chemical reactions in the second half of the 20th
century Hill (1977); Schnakenberg (1976); Jiu-li  _et al._ (1984); Mou  _et
al._ (1986); Ross  _et al._ (1988). But it is only in the 21st century that
stochastic thermodynamics was systematized to characterize dissipative
processes occurring arbitrarily far from equilibrium and their fluctuations
Jarzynski (2011); Seifert (2012); den Broeck and Esposito (2015). In recent
years, building on the work of Prigogine and co-workers Kondepudi and
Prigogine (2015), stochastic thermodynamics was formulated for stochastic
Gaspard (2004); Schmiedl and Seifert (2007); Rao and Esposito (2018a) and
extended to deterministic Qian and Beard (2005); Polettini and Esposito
(2014); Rao and Esposito (2016) chemical reaction networks (CRNs) maintained
in a nonequilibrium regime with particle reservoirs called chemostats.

Nonequilibrium thermodynamics of CRNs is a powerful tool to analyze chemical
complexity. It has been used for reaction diffusion systems to quantify the
energetic cost of creating patterns Falasco  _et al._ (2018), sustaining
chemical waves Avanzini  _et al._ (2019) and powering a chemical cloaking
device Avanzini  _et al._ (2020a). The growth process of macromolecules like
copolymers Andrieux and Gaspard (2008); Blokhuis and Lacoste (2017); Gaspard
(2020a) and biomolecules Rao  _et al._ (2015) has been characterized. Also
applications to chaotic CRNs have been considered Gaspard (2020b). Recently,
the connections of the theory to information geometry has been investigated
Yoshimura and Ito (2020). Furthermore, different strategies to maintain
nonequilibrium regimes, like the use of finite chemostats Fritz  _et al._
(2020) and serial transfers in closed reactors Blokhuis  _et al._ (2018), have
been explored.

All current formulations of nonequilibrium thermodynamics of CRNs are however
based on the assumption that the chemical species do not interact, except via
chemical reactions. The purpose of this paper is to generalize nonequilibrium
thermodynamics to non-ideal CRNs described by deterministic rate equations.
Interactions are treated with a mean-field approach using the same activity
coefficients introduced in equilibrium thermodynamics Hołyst and Poniewierski
(2012), but expressed in terms of nonequilibrium concentrations. Exploiting
the local detailed balance assumption, we impose thermodynamic consistency on
the dynamics. The standard mass-action kinetics has to be modified introducing
concentration dependent kinetic constants to account for the effects of
interactions as already recognized in the literature Boudart (1968); Laidler
(1987); Butt (1999).

Crucially, we generalize the decomposition of the entropy production developed
in Ref. Rao and Esposito, 2018a to non-ideal CRNs. This allows us to determine
the thermodynamic potential and the forces driving non-ideal chemical systems
out of equilibrium. The expression of the forces is remarkably similar to the
corresponding one for ideal CRNs since all the effects due to the interactions
are hidden in the activity coefficients. The thermodynamic potential still
acts as a Lyapunov function in detailed balanced systems, and is always lower
bounded by its equilibrium value. However, the difference between the
nonequilibrium and the equilibrium value of the potential cannot be solely
expressed in terms of a relative entropy as for ideal CRNs. Furthermore, the
exchange currents controlling the chemostatted species cannot be defined only
dynamically as in ideal CRNs, but they must account for the interactions via
the activity coefficients.

To develop our theory in a self-contained way, we proceed as follows. In Sec.
II we discuss the network theory of chemical reactions. After introducing the
basic notation in Subs II.1, we examine the dynamics in Subs. II.2. We then
define conservation laws and cycles in Subs. II.3 and II.4, respectively. In
Subs. II.5, we identify under which conditions the existence of equilibrium
states is granted, namely, CRNs are detailed balanced. We develop our
thermodynamic theory in Sec. III. After introducing the general setup in Subs.
III.1, we connect the dynamics to the thermodynamics using the local detailed
balance assumption in Subs. III.2. The expressions for the exchange currents
controlling the chemostatted species are derived in Subs. III.3. We then make
use of conservation laws to decompose the entropy production rate in Subs.
III.4. At steady state, the entropy production has a second physically
meaningful decomposition based on the cycles as we show in Subs. III.5. The
properties of the thermodynamic potential are discussed in Subs. III.6, where
we derive its lower bound, and in Subs. III.7, where we show that it cannot be
written as its equilibrium value plus a relative entropy as in ideal systems.

Throughout the manuscript, we use the CRN (9) to illustrate our results. This
example represents a catalytic process power by a proton transfer in a system
where the relevant interactions are the electrostatic interactions between
charged species. We summarize our results and discuss their implications in
Sec. IV.

This work should be particularly relevant for energetic considerations in
biosystems, where many metabolic processes involve storing energy in ion
gradients across membranes Voet and Voet (2010), and electrochemical systems,
where CRNs are coupled to electronic circuits Bard and Faulkner (2001).

##  II Non-ideal Chemical Reaction Networks

###  II.1 Setup

We consider systems composed of interacting chemical species, identified by
the label $\alpha\in Z$, homogeneously distributed in a constant volume $V$.
Our description can treat gas phases and dilute solutions where the volume of
the solution is overwhelmingly dominated by the solvent. The chemical species
undergo elementary chemical reactions G. (1993), identified by the index
$\rho\in\mathcal{R}$,

$\boldsymbol{\alpha}\cdot\boldsymbol{\nu}_{+\rho}\ch{<=>[+\rho][-\rho]}\boldsymbol{\alpha}\cdot\boldsymbol{\nu}_{-\rho}\,,$
(1)

with $\boldsymbol{\alpha}=(\dots,\alpha,\dots)^{\intercal}$ the vector of
chemical species and $\boldsymbol{\nu}_{\pm\rho}$ the vector of stoichiometric
coefficients of the forward/backward reaction $\pm\rho$. The set of chemical
reactions (1) defines the CRN. Our framework considers open CRNs. We split the
set of all the species $Z$ into two disjoint subsets: the internal species $X$
and the chemostatted species $Y$. The former undergo only the chemical
reactions (1). The latter too, but they are also externally exchanged. We will
discuss the dynamical implications and the thermodynamic meaning of the
chemostatting procedure in Subs. II.2 and Subs. III.3, respectively.

###  II.2 Dynamics

The state of deterministic CRNs with constant volume is specified by the
concentration vector
${\boldsymbol{z}}(t)=(\dots,[\alpha](t),\dots)^{\intercal}$. Its dynamics
follows the rate equation

$\mathrm{d}_{t}{\boldsymbol{z}}(t)=\mathbb{S}\boldsymbol{j}({\boldsymbol{z}}(t))+\boldsymbol{I}(t)\,,$
(2)

where we introduced the stoichiometric matrix $\mathbb{S}$, the current vector
$\boldsymbol{j}({\boldsymbol{z}})$ and the exchange current vector
$\boldsymbol{I}(t)$. The first term on the r.h.s. of Eq. (2), i.e.,
$\mathbb{S}\boldsymbol{j}({\boldsymbol{z}}(t))$, accounts for the
concentration changes due to the chemical reactions (1). The second term,
i.e., $\boldsymbol{I}(t)$, accounts for the external matter flows.

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

$j^{\rho}({\boldsymbol{z}})=j^{+\rho}({\boldsymbol{z}})-j^{-\rho}({\boldsymbol{z}})\,.$
(3)

In case of ideal CRNs, the fluxes $j^{\pm\rho}({\boldsymbol{z}})$ are
expressed in terms of mass-action kinetics de Groot and Mazur (1984); Pekař
(2005); Laidler (1987):

$j^{\pm\rho}({\boldsymbol{z}})={k}^{\text{id}}_{\pm\rho}{\boldsymbol{z}}^{\boldsymbol{\nu}_{\pm\rho}}\,,$
(4)

where ${k}^{\text{id}}_{\pm\rho}$ are the kinetic constants of the
forward/backward reaction $\pm\rho$ and we used the following notation
$\boldsymbol{a}^{\boldsymbol{b}}=\prod_{i}a_{i}^{b_{i}}$. In case of non-ideal
CRNs, mass-action kinetics is not thermodynamically consistent as we show in
Subs. III.2. To take into account the effects of the interactions, we assume
that the kinetic constants may depend on the concentrations
$\big{\\{}{k}_{\pm\rho}({\boldsymbol{z}})\big{\\}}$ as widely acknowledged in
the literature Boudart (1968); Laidler (1987); Butt (1999). Thus, the reaction
fluxes are given by the general expressions

$j^{\pm\rho}({\boldsymbol{z}})={k}_{\pm\rho}({\boldsymbol{z}}){\boldsymbol{z}}^{\boldsymbol{\nu}_{\pm\rho}}\,.$
(5)

The specific ${\boldsymbol{z}}$-dependence of
$\\{{k}_{\pm\rho}({\boldsymbol{z}})\\}$ derives from the particular model used
to describe the interactions. We develop our theory only assuming that
${k}_{\pm\rho}({\boldsymbol{z}})$ satisfy the conditions introduced in Subs.
II.5 and Subs. III.2 to be thermodynamically consistent.

The exchange current vector $\boldsymbol{I}(t)$ specifies the external matter
flows. It has null entries for the internal species, i.e., $I^{\alpha}(t)=0$
for $\alpha\in X$: the concentration of the internal species changes only
because of the chemical reactions (1) by definition. The entries of
$\boldsymbol{I}(t)$ for the chemostatted species, i.e., $I^{\alpha}(t)$ for
$\alpha\in Y$, are derived in Subs. III.3 because the thermodynamic meaning of
the chemostatting procedure must be taken into account.

We can apply the splitting of the chemical species into internal and
chemostatted ones to the stoichiometric matrix,

$\mathbb{S}=\begin{pmatrix}\mathbb{S}{}^{X}\\\ \mathbb{S}{}^{Y}\\\
\end{pmatrix}\,,$ (6)

and the concentration vector
${\boldsymbol{z}}=(\boldsymbol{x},\boldsymbol{y})$. Analogously, the rate
equation (2) becomes

$\displaystyle\mathrm{d}_{t}\boldsymbol{x}(t)=\mathbb{S}{}^{X}\boldsymbol{j}(\boldsymbol{x}(t),\boldsymbol{y}(t))\,,$
(7)
$\displaystyle\mathrm{d}_{t}\boldsymbol{y}(t)=\mathbb{S}{}^{Y}\boldsymbol{j}(\boldsymbol{x}(t),\boldsymbol{y}(t))+\boldsymbol{{I}^{Y}}(t)\,,$
(8)

with
$\boldsymbol{{I}^{Y}}(t)=(\dots,I^{\alpha}(t),\dots)^{\intercal}_{\alpha\in
Y}$ collecting the not null entries of $\boldsymbol{I}(t)$. Note that the Eqs.
(7) and (8) are only a reformulation of Eq. (2).

######  Example 1.

We consider the following CRN (see also Fig. 1):

Figure 1: Drawing of the CRN (9). A membrane (green dashed line) divides the
system into compartments. A membrane enzyme $\ch{E}$ binds the substrate
$\ch{S}$. The interconversion of the substrate $\ch{S}$ into the product
$\ch{P}$ is coupled to the proton $\ch{H^{+}}$ transfer from compartment
${{(1)}}$ to compartment ${{(2)}}$. In parallel, the anions A^- move from
compartment ${{(1)}}$ to compartment ${{(2)}}$. The $\ch{S}$, $\ch{P}$,
$\ch{H^{+}}$ are chemostatted (yellow arrows). Other cations C^+ and anions
D^- are present in compartment ${{(1)}}$ and ${{(2)}}$.

$\begin{split}\ch{A^{-}_{{{(1)}}}&<=>[${+1}$][${-1}$]A^{-}_{{{(2)}}}}\\\
\ch{S_{{{(2)}}}+E_{{{(3)}}}&<=>[${+2}$][${-2}$]ES_{{{(3)}}}}\\\
\ch{H^{+}_{{{(1)}}}+ES_{{{(3)}}}&<=>[${+3}$][${-3}$]EP_{{{(3)}}}+H^{+}_{{{(2)}}}}\\\
\ch{EP_{{{(3)}}}&<=>[${+4}$][${-4}$]E_{{{(3)}}}+P_{{{(2)}}}}\end{split}$ (9)

where the superscripts indicate the electric charge and the subscripts
indicate the compartment where a chemical species is located. Here, $\alpha\in
Z=\\{\ch{E_{{{(3)}}}},\ch{ES_{{{(3)}}}},\ch{EP_{{{(3)}}}},\ch{A^{-}_{{{(1)}}}},\ch{A^{-}_{{{(2)}}}},\ch{H^{+}_{{{(1)}}}},\ch{H^{+}_{{{(2)}}}},\ch{S_{{{(2)}}}},\ch{P_{{{(2)}}}}\\}$.
The CRN (9) represents the transformation of the substrate S into the product
P inside compartment ${{(2)}}$ catalysed by the membrane enzyme E. The
interconversion of the complex ES into EP on the membrane (compartment
${{(3)}}$) is coupled to the transfer of protons H^+ from compartment
${{(1)}}$ to compartment ${{(2)}}$. Independently, anions A^- are free to move
from compartment ${{(1)}}$ to compartment ${{(2)}}$. Other cations C^+ and
anions D^- are present in compartment ${{(1)}}$ and ${{(2)}}$, but they are
not involved in any chemical reaction. For this reason, the species C^+ and
D^- are not included in $Z$. Because of the coupling between the catalytic
reaction and the proton transfer, the CRN (9) resembles the working mechanism
of the membrane enzyme ATP synthase Huxley  _et al._ (2000). All the protons,
the substrate and the product are chemostatted, i.e.,
$Y=\\{\ch{H^{+}_{{{(1)}}}},\ch{H^{+}_{{{(2)}}}},\ch{S_{{{(2)}}}},\ch{P_{{{(2)}}}}\\}$
and
$X=\\{\ch{E_{{{(3)}}}},\ch{ES_{{{(3)}}}},\ch{EP_{{{(3)}}}},\ch{A^{-}_{{{(1)}}}},\ch{A^{-}_{{{(2)}}}}\\}$.
For this kind of systems the strongest interactions are the electrostatic
interactions between charged species inside the same compartment.

The stoichiometric matrix of the CRN (9) is specified as

${\mathbb{S}}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{E_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$-1$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{ES_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{EP_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$-1$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt\cr}}}}\right)$}}\,.$ (10)

∎

###  II.3 Conservation Laws

The linearly independent vectors $\\{\boldsymbol{\ell}^{\lambda}\\}$ in the
cokernel of the stoichiometric matrix

$\boldsymbol{\ell}^{\lambda}\cdot\mathbb{S}=0$ (11)

are the so-called conservation laws Polettini and Esposito (2014); Rao and
Esposito (2016). All the scalars
$L^{\lambda}({\boldsymbol{z}})\equiv\boldsymbol{\ell}^{\lambda}\cdot{\boldsymbol{z}}$
would be conserved quantities if the CRN was closed, namely,
$\boldsymbol{I}(t)=0$. Indeed,
$\mathrm{d}_{t}L^{\lambda}({\boldsymbol{z}}(t))=\boldsymbol{\ell}^{\lambda}\cdot\mathbb{S}\boldsymbol{j}({\boldsymbol{z}}(t))=0$.
In closed CRNs the total mass is conserved and, consequently, the set
$\\{\boldsymbol{\ell}^{\lambda}\\}$ is never empty.

When CRNs are open, some conservation laws do not correspond anymore to
conserved quantities. Hence, we split the set of conservation laws
$\\{\boldsymbol{\ell}^{\lambda}\\}$ into two disjoint subsets: the unbroken
conservation laws $\\{\boldsymbol{\ell}^{\lambda_{u}}\\}$ and the broken
conservation laws $\\{\boldsymbol{\ell}^{\lambda_{b}}\\}$. The unbroken
conservation laws are the largest subset of conservation laws that can be
written with null entries for the chemostatted species, i.e.,
$\ell^{\lambda_{u}}_{\alpha}=0$ for $\alpha\in Y$. Their corresponding scalar
quantities
$L^{\lambda_{u}}({\boldsymbol{z}})=\boldsymbol{\ell}^{\lambda_{u}}\cdot{\boldsymbol{z}}$
are still conserved:

$\mathrm{d}_{t}L^{\lambda_{u}}({\boldsymbol{z}}(t))=\underbrace{\boldsymbol{\ell}^{\lambda_{u}}\cdot\mathbb{S}\boldsymbol{j}({\boldsymbol{z}}(t))}_{=0}+\sum_{\alpha\in
Y}\underbrace{\ell^{\lambda_{u}}_{\alpha}}_{=0}I^{\alpha}(t)=0\,.$ (12)

The broken conservation laws are the other conservation laws
$\\{\boldsymbol{\ell}^{\lambda_{b}}\\}=\\{\boldsymbol{\ell}^{\lambda}\\}\setminus\\{\boldsymbol{\ell}^{\lambda_{u}}\\}$.
Their corresponding scalar quantities
$L^{\lambda_{b}}({\boldsymbol{z}})=\boldsymbol{\ell}^{\lambda_{b}}\cdot{\boldsymbol{z}}$
are, in general, not conserved and they give rise to balance equations:

$\mathrm{d}_{t}L^{\lambda_{b}}({\boldsymbol{z}}(t))=\underbrace{\boldsymbol{\ell}^{\lambda_{b}}\cdot\mathbb{S}\boldsymbol{j}({\boldsymbol{z}}(t))}_{=0}+\sum_{\alpha\in
Y}\underbrace{\ell^{\lambda_{b}}_{\alpha}}_{\neq 0}I^{\alpha}(t)\neq 0\,,$
(13)

where the variations of $L^{\lambda_{b}}({\boldsymbol{z}}(t))$ are due to the
matter exchanged through the currents $\\{I^{\alpha}(t)\\}_{\alpha\in Y}$. In
open CRNs the total mass is not conserved and, consequently, the set
$\\{\boldsymbol{\ell}^{\lambda_{b}}\\}$ is never empty.

Chemostatting a species does not always break a conservation law Rao and
Esposito (2016, 2018a); Falasco  _et al._ (2018); Rao and Esposito (2018b). We
thus distinguish the species ${Y_{p}}\subseteq Y$ that break the conservation
laws from the others ${Y_{f}}=Y\setminus{Y_{p}}$. This classification is not
unique and different choices have different physical meanings. It should be
motivated by the physical role played by each chemostatted species as it will
become clear in Subs. III.4. We apply the same splitting to the stoichiometric
matrix $\mathbb{S}{}^{Y}$ and the matrix $\mathbb{L}^{{b}}$ whose rows are the
broken conservation laws $\\{\boldsymbol{\ell}^{\lambda_{b}}\\}$:

$\displaystyle\mathbb{S}{}^{Y}=\begin{pmatrix}\mathbb{S}{}^{Y_{f}}\\\
\mathbb{S}{}^{Y_{p}}\end{pmatrix}\,,$
$\displaystyle\mathbb{L}^{{b}}=\begin{pmatrix}\mathbb{L}^{{b}}_{X},\mathbb{L}^{{b}}_{{Y_{f}}},\mathbb{L}^{{b}}_{{Y_{p}}}\end{pmatrix}\,,$
(14)

where $\mathbb{L}^{{b}}_{X}$, $\mathbb{L}^{{b}}_{{Y_{f}}}$, and
$\mathbb{L}^{{b}}_{{Y_{p}}}$ have
$\\{{\ell}^{\lambda_{b}}_{\alpha}\\}_{\alpha\in X}$,
$\\{{\ell}^{\lambda_{b}}_{\alpha}\\}_{\alpha\in{Y_{f}}}$, and
$\\{{\ell}^{\lambda_{b}}_{\alpha}\\}_{\alpha\in{Y_{p}}}$ as entries,
respectively. Note that the number of ${Y_{p}}$ species is equal to the number
of broken conservation laws, by definition. The matrix
$\mathbb{L}^{{b}}_{{Y_{p}}}$ is then square and nonsingular, and so it can be
inverted. We will exploit this properties in Subs. III.4 where we decompose
the entropy production rate.

######  Example 2.

Given the stoichiometric matrix (10) of the CRN (9), there are five
conservation laws,

$\begin{split}\boldsymbol{\ell}^{\ch{E}}=\hbox{}\vbox{\kern
0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{E_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{ES_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{EP_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr}}}}\right)$}}\,,\text{ }\text{ }\text{
}\text{ }\boldsymbol{\ell}^{\ch{A}}=&\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern
0.0pt\kern 2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{E_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{ES_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{EP_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr}}}}\right)$}}\,,\text{ }\text{ }\text{
}\text{ }\boldsymbol{\ell}^{\ch{H}}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern
0.0pt\kern 2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{E_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{ES_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{EP_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr}}}}\right)$}}\,,\\\
\boldsymbol{\ell}^{\ch{HS}}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern
0.0pt\kern 2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{E_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{ES_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{EP_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr}}}}\right)$}}\,,&\text{ }\text{ }\text{
}\text{ }\boldsymbol{\ell}^{\ch{S}}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern
0.0pt\kern 2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{E_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{ES_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{EP_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr}}}}\right)$}}\,.\end{split}$ (15)

Four of them have a clear physical meaning: the total concentration of the
enzyme is given by
$L^{\ch{E}}=\boldsymbol{\ell}^{\ch{E}}\cdot{\boldsymbol{z}}=[\ch{E_{{(3)}}}]+[\ch{ES_{{(3)}}}]+[\ch{EP_{{(3)}}}]$,
of the anions by
$L^{\ch{A}}=\boldsymbol{\ell}^{\ch{A}}\cdot{\boldsymbol{z}}=[\ch{A^{-}_{{(1)}}}]+[\ch{A^{-}_{{(2)}}}]$,
of the protons by
$L^{\ch{H}}=\boldsymbol{\ell}^{\ch{H}}\cdot{\boldsymbol{z}}=[\ch{H^{+}_{{(1)}}}]+[\ch{H^{+}_{{(2)}}}]$,
and of the substrate by
$L^{\ch{S}}=\boldsymbol{\ell}^{\ch{S}}\cdot{\boldsymbol{z}}=[\ch{ES_{{(3)}}}]+[\ch{EP_{{(3)}}}]+[\ch{S_{{(2)}}}]+[\ch{P_{{(2)}}}]$.
The conserved quantity
$L^{\ch{HS}}=\boldsymbol{\ell}^{\ch{HS}}\cdot{\boldsymbol{z}}=[\ch{ES_{{(3)}}}]+[\ch{S_{{(2)}}}]+[\ch{H^{+}_{{(2)}}}]$
comes from the coupling between the consumption of the substrate and the
proton transfer.

When the species $\ch{H^{+}_{{{(1)}}}}$, $\ch{H^{+}_{{{(2)}}}}$,
$\ch{S_{{{(2)}}}}$, and $\ch{P_{{{(2)}}}}$ are chemostatted, the three
conservation laws $\boldsymbol{\ell}^{\ch{H}}$, $\boldsymbol{\ell}^{\ch{HS}}$,
and $\boldsymbol{\ell}^{\ch{S}}$ are broken. We identify the set ${Y_{p}}$ as
${Y_{p}}=\\{\ch{H^{+}_{{{(2)}}}},\ch{S_{{{(2)}}}},\ch{P_{{{(2)}}}}\\}$, but
different choices are possible (see the discussion in Appendix A of Ref.
Avanzini  _et al._ , 2019). Thus,
${Y_{f}}=Y\setminus{Y_{p}}=\\{\ch{H^{+}_{{{(1)}}}}\\}$. The matrix
$\mathbb{L}^{{b}}$ whose rows are the broken conservation laws is given by

$\mathbb{L}^{{b}}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{E_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{ES_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{EP_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr}}}}\right)$}}\,,$ (16)

where the grey vertical lines mark the split of $\mathbb{L}^{{b}}$ into
$\mathbb{L}^{{b}}_{X}$, $\mathbb{L}^{{b}}_{{Y_{f}}}$, and
$\mathbb{L}^{{b}}_{{Y_{p}}}$. ∎

###  II.4 Cycles

The linearly independent vectors $\\{\boldsymbol{c}_{\iota}\\}$ in the kernel
of the stoichiometric matrix

$\mathbb{S}\boldsymbol{c}_{\iota}=0$ (17)

are the so-called internal cycles Polettini and Esposito (2014); Rao and
Esposito (2016). They represent sequences of reactions that upon completion
leave the all concentrations unchanged. Any linear combination of internal
cycles gives a steady-state current vector of closed CRNs,
$\mathrm{d}_{t}{\boldsymbol{z}}(t)=\mathbb{S}\overline{\boldsymbol{j}}=0$ with
$\overline{\boldsymbol{j}}=\sum_{\iota}\boldsymbol{c}_{\iota}\psi^{\iota}$.

When CRNs are open, the steady-state current vector must i) leave the
concentrations of the internal species unchanged, i.e.,
$\mathbb{S}{}^{X}\overline{\boldsymbol{j}}=0$, and ii) be balanced by the
exchanged current vector, i.e.,
$\mathbb{S}{}^{Y}\overline{\boldsymbol{j}}+\boldsymbol{{I}^{Y}}=0$. Thus,
$\overline{\boldsymbol{j}}\in\ker(\mathbb{S}{}^{X})$ and, consequently, can be
written as a linear combination of the right-null eigenvectors

$\mathbb{S}{}^{X}\boldsymbol{c}_{\xi}=0\,.$ (18)

The set $\\{\boldsymbol{c}_{\xi}\\}$ includes the internal cycles
$\\{\boldsymbol{c}_{\iota}\\}$ and, in general, other vectors called emergent
cycles $\\{\boldsymbol{c}_{\epsilon}\\}$. The steady-state current vector can
therefore be written as

$\overline{\boldsymbol{j}}=\sum_{\iota}\boldsymbol{c}_{\iota}\psi^{\iota}+\sum_{\epsilon}\boldsymbol{c}_{\epsilon}\psi^{\epsilon}\,.$
(19)

Because of the rank-nullity theorem for the stoichiometric matrix $\mathbb{S}$
and $\mathbb{S}{}^{X}$, every time a species is chemostatted either a
conservation law is broken or an emergent cycle arises, namely,

$|Y|=|{\lambda_{b}}|+|\epsilon|\,,$ (20)

with $|Y|$, $|{\lambda_{b}}|$ and $|\epsilon|$ the number of chemostatted
species, broken conservation laws and emergent cycles, respectively. This
implies that the number of ${Y_{f}}$ species is equal to the number of
emergent cycles, i.e., $|{Y_{f}}|=|\epsilon|$.

######  Example 3.

Given the stoichiometric matrix (10) of the CRN (9), there is only one
emergent cycle and no internal cycles,

$\boldsymbol{c}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr}}}}\right)$}}\,.$ (21)

∎

###  II.5 Equilibrium

The equilibrium steady-state ${\boldsymbol{z}}_{\text{eq}}$ of the rate
equation (2), if it exists, is characterized by vanishing reaction currents,

$\boldsymbol{j}({\boldsymbol{z}}_{\text{eq}})=0\,.$ (22)

This, together with Eq. (5), implies

$\frac{{k}_{+\rho}({\boldsymbol{z}}_{\text{eq}})}{{k}_{-\rho}({\boldsymbol{z}}_{\text{eq}})}={\boldsymbol{z}}_{\text{eq}}^{\mathbb{S}_{\rho}}\,.$
(23)

If a solution ${\boldsymbol{z}}_{\text{eq}}$ of Eq. (22) (or equivalently of
Eq. (23)) exists, the CRN is said to be detailed balanced. Closed CRNs must be
detailed balanced for thermodynamic consistency. It is not granted that open
CRNs are detailed balanced. It depends on the chemostatting procedure.

Here, we show that open CRNs are necessarily detailed balanced if the number
of chemostatted species is equal to the number of broken conservation laws,
i.e., $Y={Y_{p}}$, or, equivalently, there are no emergent cycles. In Subs.
III.4, we show that equilibrium states can exist also when $Y\neq{Y_{p}}$ and
there are emergent cycles. For thermodynamic consistency, namely the local
detailed balance condition that will be introduced in Eq. (40), the product of
the forward rate constants along every internal cycle (17) must equal that of
the backward rate constants:

$\prod_{\rho}\left(\frac{{k}_{+\rho}({\boldsymbol{z}})}{{k}_{-\rho}({\boldsymbol{z}})}\right)^{c^{\rho}_{\iota}}=1\,,$
(24)

where $c^{\rho}_{\iota}$ is the $\rho$ entry of the $\iota$ internal cycle.
This is know as Wegscheider’s condition in the framework of ideal CRNs
Schuster and Schuster (1989). Hence, only $|\rho|-|\iota|$ ratios
${k}_{+\rho}({\boldsymbol{z}})/{k}_{-\rho}({\boldsymbol{z}})$ are independent
(with $|\rho|$ the number of reactions and $|\iota|$ the number of internal
cycles). This means that Eq. (23) imposes $|\rho|-|\iota|$ constraints. This
is independent on whether CRNs are closed or open.

In closed CRNs, only $|\alpha|-|\lambda|$ of the equilibrium concentrations
are independent (with $|\alpha|$ the number of all the species in $Z$ and
$|\lambda|$ the number of conservation laws). The others are fixed by the
initial conditions through the conserved quantities:

$L^{\lambda}=\boldsymbol{\ell}^{\lambda}\cdot{\boldsymbol{z}}(0)=\boldsymbol{\ell}^{\lambda}\cdot{\boldsymbol{z}}_{\text{eq}}\,.$
(25)

Therefore, Eq. (23) imposes $|\rho|-|\iota|$ constraints for
$|\alpha|-|\lambda|$ independent variables. Because of the rank-nullity
theorem for the stoichiometric matrix $\mathbb{S}$,

$|\rho|-|\iota|=|\alpha|-|\lambda|\,,$ (26)

the existence of a solution to Eq. (23) is granted. Closed CRNs must be
detailed balanced and hence each solution of Eq. (23) must be physically
meaningful, i.e., $[\alpha]_{\text{eq}}\in\mathbb{R}_{\geq 0}$
$\forall\alpha\in Z$. This imposes constraints on the possible expressions of
$\big{\\{}{k}_{\pm\rho}({\boldsymbol{z}})\big{\\}}$.

In open CRNs, $|\lambda|-|{\lambda_{b}}|$ concentrations are fixed by the
initial conditions through the unbroken conserved quantities:

$L^{\lambda_{u}}=\boldsymbol{\ell}^{{\lambda_{u}}}\cdot{\boldsymbol{z}}(0)=\boldsymbol{\ell}^{{\lambda_{u}}}\cdot{\boldsymbol{z}}_{\text{eq}}\,.$
(27)

Other $|Y|$ concentrations are constrained by the chemostatting procedure. The
final number of independent concentrations is given by

$|\alpha|-(|\lambda|-|{\lambda_{b}}|)-|Y|=|\alpha|-|\lambda|-|{Y_{f}}|=|\rho|-|\iota|-|{Y_{f}}|\,,$
(28)

where we used $|Y|-|{\lambda_{b}}|=|{Y_{f}}|$ and Eq. (26). Therefore, Eq.
(23) imposes $|\rho|-|\iota|$ constraints for $|\rho|-|\iota|-|{Y_{f}}|$
independent variables. The existence of a solution is always granted if and
only if $|{Y_{f}}|=0$, namely, all the chemostatted species break a
conservation law (i.e., $Y={Y_{p}}$). Equivalently, the existence of a
solution is always granted if and only if there are no emergent cycles,
$|\epsilon|=0$. The solution is physical meaningful because it must also be a
solution of the corresponding closed CRN.

These are sufficient conditions for the existence of an equilibrium state
${\boldsymbol{z}}_{\text{eq}}$. For ideal detailed balanced CRNs, i.e.,
evolving according to mass-action kinetics (4), there is a unique equilibrium
state for every stoichiometric compatibility class
$\Omega\left(\big{\\{}L^{\lambda_{u}}\big{\\}}\right)$, namely, for every
manifold in the concentration space characterized by specific and unique
values of all the conserved quantities Horn and Jackson (1972); Schuster and
Schuster (1989). For non-ideal CRNs, this property is not granted anymore and
there might be more than one equilibrium state for every stoichiometric
compatibility class. The existence of multiple equilibria in the same
stoichiometric compatibility class will become relevant in Subs. III.6.

######  Example 4.

It is not granted that the open CRN (9) is detailed balanced because the set
${Y_{f}}$ is not empty, i.e., ${Y_{f}}=\\{\ch{H^{+}_{{{(1)}}}}\\}$, and there
is the emergent cycle (21). ∎

##  III Thermodynamics

###  III.1 Setup

Nonequilibrium thermodynamics of CRNs presumes that all degrees of freedom
other than concentrations are equilibrated. The temperature $T$ is set by a
thermal reservoir (e.g., the solvent in dilute solutions), diffusion processes
are fast enough to keep the chemical species homogeneously distributed, and
the fields responsible for the interactions, e.g., electrostatic fields, relax
instantaneously to their mean-field values. For dilute solutions, the pressure
$p$ is set by the environment of the solution. In absence of reactions, the
CRN would be an equilibrated mixture. In this way, thermodynamic state
functions can be specified by their equilibrium form but expressed in terms of
nonequilibrium concentrations.

We assume that the free energy of non-ideal CRNs, for a given temperature $T$
and volume $V$, is a mean-field function
$G({\boldsymbol{z}},{\boldsymbol{e}}{})$ of the concentrations
${\boldsymbol{z}}$ and of some externally controlled parameters
${\boldsymbol{e}}{}=(\dots,{e}{}_{i},\dots)^{\intercal}$ that tune the
interactions. We chose the symbol $G$ for the free energy because it
corresponds to the Gibbs free energy when dealing with dilute solutions. If we
considered gas phases in a constant volume, $G$ would be the Helmholtz free
energy. For instance, Debye-Hückel theory McQuarrie (1975) derives
$G({\boldsymbol{z}},{\boldsymbol{e}}{})$ for electrolyte solutions, and van
der Waals theory Kondepudi and Prigogine (2015) provides
$G({\boldsymbol{z}},{\boldsymbol{e}}{})$ for a gas of interacting and finite-
size molecules. Standard computational approaches are also available to
compute the free energy of generic non-ideal systems Chipot and Pohorille
(2007); Lelièvre  _et al._ (2010). Without loss of generality,
$G({\boldsymbol{z}},{\boldsymbol{e}}{})$ can be written as the sum of an ideal
term $G^{\text{id}}({\boldsymbol{z}})$ and a term accounting for the effects
of the interactions $G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})$:

$G({\boldsymbol{z}},{\boldsymbol{e}}{})=G^{\text{id}}({\boldsymbol{z}})+G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})\,.$
(29)

The free energy contribution carried by each species $\alpha$ is given by the
chemical potential Holyst and Poniewierski (2012):

${\mu}_{\alpha}({\boldsymbol{z}})=\frac{\partial
G({\boldsymbol{z}},{\boldsymbol{e}}{})}{\partial[\alpha]}=\mu^{\text{id}}_{\alpha}([\alpha])+RT\ln\gamma_{\alpha}({\boldsymbol{z}},{\boldsymbol{e}}{})\,.$
(30)

By assumption, $\mu^{\text{id}}_{\alpha}([\alpha])$ is the ideal chemical
potential in either gas phases at constant volume $V$ or in dilute solutions,

$\mu^{\text{id}}_{\alpha}([\alpha])=\frac{\partial
G^{\text{id}}({\boldsymbol{z}})}{\partial[\alpha]}=\mu^{\circ}_{\alpha}+RT\ln[\alpha]\,,$
(31)

with $R$ the gas constant and $\mu^{\circ}_{\alpha}$ the standard chemical
potential. The so-called activity coefficient
$\gamma_{\alpha}({\boldsymbol{z}},{\boldsymbol{e}}{})$ accounts for the
effects of the interactions

$RT\ln\gamma_{\alpha}({\boldsymbol{z}},{\boldsymbol{e}}{})=\frac{\partial
G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})}{\partial[\alpha]}\,.$ (32)

In equilibrium thermodynamics of non-ideal CRNs,
$\gamma_{\alpha}({\boldsymbol{z}},{\boldsymbol{e}}{})$ is often treated as a
constant parameter, while in our framework it evolves in time since the
concentrations ${\boldsymbol{z}}$ are dynamical variables. The specific
expression of $\gamma_{\alpha}({\boldsymbol{z}},{\boldsymbol{e}}{})$, or
equivalently of $G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})$, depends
on the particular model used to describe the interactions. We develop our
theory only assuming that $G({\boldsymbol{z}},{\boldsymbol{e}}{})$ is lower
bounded and increases superlinearly as concentrations go to infinity. These
constraints are necessary to guarantee thermodynamic consistency, i.e.,
detailed balanced CRNs relax towards an equilibrium state, as we show in Subs.
III.6. The free energy provided by the Debye-Hückel theory McQuarrie (1975)
satisfies these conditions: as concentrations become infinitely large,
$G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})$ goes to minus infinity
slower than $G^{\text{id}}({\boldsymbol{z}})$ goes to infinity.

We collect all the chemical potentials in the vector

$\boldsymbol{{\mu}}({\boldsymbol{z}})=\big{(}\dots,{\mu}_{\alpha}({\boldsymbol{z}}),\dots\big{)}^{\intercal}\,,$
(33)

that can be written as

$\begin{split}\boldsymbol{{\mu}}({\boldsymbol{z}})&=\boldsymbol{\mu^{\circ}}+RT\ln{\boldsymbol{z}}+RT\ln{\boldsymbol{\gamma}}({\boldsymbol{z}},{\boldsymbol{e}}{})\\\
&=\boldsymbol{\mu^{\circ}}+RT\ln{\boldsymbol{z}}+\boldsymbol{\nabla}_{\boldsymbol{z}}G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})\,,\end{split}$
(34)

with
$\boldsymbol{\mu^{\circ}}=(\dots,\mu^{\circ}_{\alpha},\dots)^{\intercal}$,
${\boldsymbol{\gamma}}({\boldsymbol{z}},{\boldsymbol{e}}{})=(\dots,\gamma_{\alpha}({\boldsymbol{z}},{\boldsymbol{e}}{}),\dots)^{\intercal}$
and
$\boldsymbol{\nabla}_{\boldsymbol{z}}=(\dots,\partial/\partial[\alpha],\dots)^{\intercal}$.

######  Example 5.

We assume that the only relevant interactions in the CRN (9) are the
electrostatic interactions between charged species inside the same
compartment. For illustrative reasons, we describe the electrostatic potential
in each compartment as a linear function of the net charge density instead of
using the Debye-Hückel theory McQuarrie (1975). This is consistent with a
standard model of the transmembrane potential in mitochondria Magnus and
Keizer (1997, 1998); Fall and Keizer (2001). The interaction free energy
$G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})$ is thus given by the
electrostatic potential energy. It can be written as

$G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})=G^{\text{in}}_{{{(1)}}}({\boldsymbol{z}},{\boldsymbol{e}}{})+G^{\text{in}}_{{(2)}}({\boldsymbol{z}},{\boldsymbol{e}}{})\,,$
(35)

where $G^{\text{in}}_{(i)}({\boldsymbol{z}},{\boldsymbol{e}}{})$, i.e., the
electrostatic potential energy in the $i$ compartment, is a quadratic function
of the net charge density in the $i$ compartment:

$G^{\text{in}}_{(i)}({\boldsymbol{z}},{\boldsymbol{e}}{})=\frac{F\kappa}{2}\big{(}\underbrace{-[\ch{A^{-}_{(i)}}]+[\ch{H^{+}_{(i)}}]+[\ch{C^{+}_{(i)}}]-[\ch{D^{-}_{(i)}}]}_{\equiv
Q_{(i)}({\boldsymbol{z}},{\boldsymbol{e}}{})}\big{)}^{2}\,,$ (36)

with $F$ the Faraday constant and $\kappa$ a generic proportionality constant.
On the membrane (compartment ${{(3)}}$), $G^{\text{in}}_{{(3)}}=0$ since the
chemical species are not charged and do not interact. Here, the parameters
${\boldsymbol{e}}{}$ are the concentrations of the species that are not
involved in the chemical reactions, i.e.,
${\boldsymbol{e}}{}=\big{(}[\ch{C^{+}_{{(1)}}}],[\ch{C^{+}_{{(2)}}}],[\ch{D^{-}_{{(1)}}}],[\ch{D^{-}_{{(2)}}}]\big{)}$,
but interact with species $\boldsymbol{\alpha}$. These concentrations may be
controlled with some external processes.

The chemical potentials of the species in the CRN (9) are specified as

$\begin{split}&{\mu}_{\ch{E_{{(3)}}}}({\boldsymbol{z}})=\mu^{\text{id}}_{\ch{E_{{(3)}}}}([\ch{E_{{(3)}}}])=\mu^{\circ}_{\ch{E}}+RT\ln[\ch{E_{{(3)}}}]\,,\\\
&{\mu}_{\ch{ES_{{(3)}}}}({\boldsymbol{z}})=\mu^{\text{id}}_{\ch{ES_{{(3)}}}}([\ch{ES_{{(3)}}}])=\mu^{\circ}_{\ch{ES}}+RT\ln[\ch{ES_{{(3)}}}]\,,\\\
&{\mu}_{\ch{EP_{{(3)}}}}({\boldsymbol{z}})=\mu^{\text{id}}_{\ch{EP_{{(3)}}}}([\ch{EP_{{(3)}}}])=\mu^{\circ}_{\ch{EP}}+RT\ln[\ch{EP_{{(3)}}}]\,,\\\
&{\mu}_{\ch{A^{-}_{{(1)}}}}({\boldsymbol{z}})=\mu^{\circ}_{\ch{A^{-}}}+RT\ln[\ch{A^{-}_{{(1)}}}]-F\kappa
Q_{{(1)}}({\boldsymbol{z}},{\boldsymbol{e}}{})\,,\\\
&{\mu}_{\ch{A^{-}_{{(2)}}}}({\boldsymbol{z}})=\mu^{\circ}_{\ch{A^{-}}}+RT\ln[\ch{A^{-}_{{(2)}}}]-F\kappa
Q_{{(2)}}({\boldsymbol{z}},{\boldsymbol{e}}{})\,,\\\
&{\mu}_{\ch{H^{+}_{{(1)}}}}({\boldsymbol{z}})=\mu^{\circ}_{\ch{H^{+}}}+RT\ln[\ch{H^{+}_{{(1)}}}]+F\kappa
Q_{{(1)}}({\boldsymbol{z}},{\boldsymbol{e}}{})\,,\\\
&{\mu}_{\ch{H^{+}_{{(2)}}}}({\boldsymbol{z}})=\mu^{\circ}_{\ch{H^{+}}}+RT\ln[\ch{H^{+}_{{(2)}}}]+F\kappa
Q_{{(2)}}({\boldsymbol{z}},{\boldsymbol{e}}{})\,,\\\
&{\mu}_{\ch{S_{{(2)}}}}({\boldsymbol{z}})=\mu^{\text{id}}_{\ch{S_{{(2)}}}}([\ch{S_{{(2)}}}])=\mu^{\circ}_{\ch{S}}+RT\ln[\ch{S_{{(2)}}}]\,,\\\
&{\mu}_{\ch{P_{{(2)}}}}({\boldsymbol{z}})=\mu^{\text{id}}_{\ch{P_{{(2)}}}}([\ch{P_{{(2)}}}])=\mu^{\circ}_{\ch{P}}+RT\ln[\ch{P_{{(2)}}}]\,,\end{split}$
(37)

where we assumed that the standard chemical potentials do not depend on the
compartment, and the net charge densities
$\big{\\{}Q_{(i)}({\boldsymbol{z}},{\boldsymbol{e}}{})\big{\\}}$ are specified
in Eq. (36). Note that the chemical potentials of the species
$\ch{E_{{(3)}}}$, $\ch{ES_{{(3)}}}$, $\ch{EP_{{(3)}}}$, $\ch{S_{{(2)}}}$ and
$\ch{P_{{(2)}}}$ are the same as for ideal CRNs. This is a consequence of
accounting only for the interactions between charge species. ∎

###  III.2 Local Detailed Balance

Here, we build the connection between dynamics and nonequilibrium
thermodynamics. We use the local detailed balance condition that binds the
ratios between the forward and backward currents
$\\{j^{+\rho}({\boldsymbol{z}})/j^{-\rho}({\boldsymbol{z}})\\}$ to the free
energies of reaction $\\{\Delta_{\rho}G({\boldsymbol{z}})\\}$:

$RT\ln\frac{j^{+\rho}({\boldsymbol{z}})}{j^{-\rho}({\boldsymbol{z}})}=-\Delta_{\rho}G({\boldsymbol{z}})\,.$
(38)

The latter is given by

$\Delta_{\rho}G({\boldsymbol{z}})=\boldsymbol{{\mu}}({\boldsymbol{z}})\cdot\mathbb{S}_{\rho}\,.$
(39)

This, together with Eq. (5), implies a constraint for the kinetic constants
${k}_{\pm\rho}({\boldsymbol{z}})$ of the chemical reactions:

$\begin{split}RT\ln\frac{{k}_{+\rho}({\boldsymbol{z}})}{{k}_{-\rho}({\boldsymbol{z}})}&=-(\boldsymbol{\mu^{\circ}}+RT\ln{\boldsymbol{\gamma}}({\boldsymbol{z}},{\boldsymbol{e}}{}))\cdot{\mathbb{S}_{\rho}}\\\
&=-(\boldsymbol{\mu^{\circ}}+\boldsymbol{\nabla}_{\boldsymbol{z}}G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{}))\cdot\mathbb{S}_{\rho}\,.\end{split}$
(40)

Equation (40) shows that the ${\boldsymbol{z}}$ dependence of the kinetic
constants, which represents the breakdown of mass-action kinetics for non-
ideal CRNs, comes from the interactions. These can affect the relative
stability of reagents and products and, consequently, the relative rate of the
forward and backward reactions. However, the local detailed balance condition
highlights only the ${\boldsymbol{z}}$ dependence of the antisymmetric part
$\sqrt{{k}_{+\rho}({\boldsymbol{z}})/{k}_{-\rho}({\boldsymbol{z}})}$ of the
kinetic constants. This does not exclude that also the symmetric part
$\sqrt{{k}_{+\rho}({\boldsymbol{z}}){k}_{-\rho}({\boldsymbol{z}})}$ of the
kinetic constants is ${\boldsymbol{z}}$ dependent as already discussed in the
literature Butt (1999). In the limit
${\boldsymbol{\gamma}}({\boldsymbol{z}},{\boldsymbol{e}}{})\to 1$, or
equivalently
$\boldsymbol{\nabla}_{\boldsymbol{z}}G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})\to
0$, the local detailed balance condition for ideal CRNs is recovered,

$RT\ln\frac{{k}^{\text{id}}_{+\rho}}{{k}^{\text{id}}_{-\rho}}=-\boldsymbol{\mu^{\circ}}\cdot\mathbb{S}_{\rho}\,,$
(41)

and, consequently, the dynamics satisfies mass-action kinetics (4).

######  Example 6.

We specify Eq. (40) for the CRN (9) assuming that the standard chemical
potentials of A^- and H^+ do not depend on the compartment (namely
$\mu^{\circ}_{\ch{A^{-}_{{{(1)}}}}}=\mu^{\circ}_{\ch{A^{-}_{{{(2)}}}}}$ and
$\mu^{\circ}_{\ch{H^{+}_{{{(1)}}}}}=\mu^{\circ}_{\ch{H^{+}_{{{(2)}}}}}$):

$\displaystyle
RT\ln\frac{{k}_{+1}({\boldsymbol{z}})}{{k}_{-1}({\boldsymbol{z}})}=F\kappa\big{(}Q_{{(2)}}({\boldsymbol{z}},{\boldsymbol{e}}{})-Q_{{(1)}}({\boldsymbol{z}},{\boldsymbol{e}}{})\big{)}\,,$
(42) $\displaystyle
RT\ln\frac{{k}_{+2}({\boldsymbol{z}})}{{k}_{-2}({\boldsymbol{z}})}=RT\ln\frac{{k}^{\text{id}}_{+2}}{{k}^{\text{id}}_{-2}}=-(\mu^{\circ}_{\ch{ES_{{{(3)}}}}}-\mu^{\circ}_{\ch{E_{{{(3)}}}}}-\mu^{\circ}_{\ch{S_{{{(2)}}}}})\,,$
(43) $\displaystyle
RT\ln\frac{{k}_{+3}({\boldsymbol{z}})}{{k}_{-3}({\boldsymbol{z}})}=-\big{(}\mu^{\circ}_{\ch{EP_{{{(3)}}}}}-\mu^{\circ}_{\ch{ES_{{{(3)}}}}}+F\kappa\big{(}Q_{{(2)}}({\boldsymbol{z}},{\boldsymbol{e}}{})-Q_{{(1)}}({\boldsymbol{z}},{\boldsymbol{e}}{})\big{)}\big{)}\,,$
(44) $\displaystyle
RT\ln\frac{{k}_{+4}({\boldsymbol{z}})}{{k}_{-4}({\boldsymbol{z}})}=RT\ln\frac{{k}^{\text{id}}_{+4}}{{k}^{\text{id}}_{-4}}=-(\mu^{\circ}_{\ch{P_{{{(2)}}}}}+\mu^{\circ}_{\ch{E_{{{(3)}}}}}-\mu^{\circ}_{\ch{EP_{{{(3)}}}}})\,.$
(45)

∎

###  III.3 Chemostatting

Thermodynamically, the chemostatting procedure directly determines the
chemical potentials $\\{{\mu}_{\alpha}\\}$ of the $Y$ species. For ideal CRNs,
where ${\mu}_{\alpha}({\boldsymbol{z}})=\mu^{\text{id}}_{\alpha}([\alpha])$,
this means that the chemostatting procedure directly determines the
concentrations, $[\alpha]=\exp(({\mu}_{\alpha}-\mu^{\circ}_{\alpha})/RT)$. For
non-ideal CRNs, the equivalence

${\mu}_{\alpha}=\mu^{\circ}_{\alpha}+RT\ln[\alpha]+RT\ln\gamma_{\alpha}({\boldsymbol{z}},{\boldsymbol{e}}{})\text{
}\text{ }\text{ }\forall\alpha\in Y\,,$ (46)

does not directly determine the concentrations $\boldsymbol{y}$ due to the
dependence on the interactions, i.e.,
$[\alpha]=\exp(({\mu}_{\alpha}-\mu^{\circ}_{\alpha})/RT)/\gamma_{\alpha}({\boldsymbol{z}},{\boldsymbol{e}}{})$
$\forall\alpha\in Y$. Hence, $\boldsymbol{I}(t)$ in Eq. (2) (or
$\boldsymbol{{I}^{Y}}(t)$ in Eq. (8)) describes the exchange currents such
that the chemical potentials of the $Y$ species are externally determined.

We now derive the explicit form of $\boldsymbol{{I}^{Y}}(t)$. The readers who
may want to skip the details of the derivation can go directly to Eq. (51).
The constraints imposed by Eq. (46) must be satisfied for every time $t$ of
the dynamics. Thus, also their time derivative must vanish:

$\begin{split}\frac{RT}{[\alpha](t)}\mathrm{d}_{t}[\alpha](t)+\boldsymbol{\nabla}_{\boldsymbol{z}}g^{\text{in}}_{\alpha}\big{(}{\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t)\big{)}\cdot\mathrm{d}_{t}{\boldsymbol{z}}(t)+&\\\
+\boldsymbol{\nabla}_{\boldsymbol{e}}{}g^{\text{in}}_{\alpha}\big{(}{\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t)\big{)}\cdot\mathrm{d}_{t}{\boldsymbol{e}}{}(t)-\mathrm{d}_{t}{\mu}_{\alpha}(t)&=0\,.\end{split}$
(47)

Here, we used the following compact notation
$g^{\text{in}}_{\alpha}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))=RT\ln\gamma_{\alpha}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))$
and we accounted for the possible time-dependent control of the chemical
potentials $\\{{\mu}_{\alpha}\\}$ and of the parameters ${\boldsymbol{e}}{}$.
The set of Eqs. (47) for every $\alpha\in Y$ can be rewritten as

$\mathbb{W}(t)(\mathbb{S}\boldsymbol{j}({\boldsymbol{z}}(t))+\boldsymbol{I}(t))+\boldsymbol{{\dot{g}}_{Y}^{\text{in}}}(t)-\mathrm{d}_{t}\boldsymbol{{\mu}}_{Y}(t)=0\,,$
(48)

by introducing the vectors

$\boldsymbol{{\dot{g}}_{Y}^{\text{in}}}(t)=\left(\dots,\boldsymbol{\nabla}_{\boldsymbol{e}}{}g^{\text{in}}_{\alpha}\big{(}{\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t)\big{)}\cdot\mathrm{d}_{t}{\boldsymbol{e}}{}(t),\dots\right)^{\intercal}_{\alpha\in
Y}$ (49)

and
$\boldsymbol{{\mu}}_{Y}(t)=(\dots,{\mu}_{\alpha}(t),\dots)^{\intercal}_{\alpha\in
Y}$, and the matrix $\mathbb{W}(t)$ whose entries are

$\mathbb{W}^{\alpha}_{\beta}(t)=\frac{RT}{[\alpha](t)}\delta_{\alpha,\beta}+\frac{\partial
g^{\text{in}}_{\alpha}\big{(}{\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t)\big{)}}{\partial[\beta]}\,,$
(50)

with $\alpha\in Y$ the row index and $\beta\in Z$ the column index. We notice
that the submatrix $\mathbb{W}^{Y}$, whose entries are
$\\{\mathbb{W}^{\alpha}_{\beta}\\}_{\alpha,\beta\in Y}$, is square and not
nonsingular so that it can be inverted. We thus obtain an explicit expression
for $\boldsymbol{{I}^{Y}}(t)$:

$\boldsymbol{{I}^{Y}}(t)=-\Big{(}\mathbb{W}^{Y}(t)\Big{)}^{-1}\Big{(}\mathbb{W}(t)\mathbb{S}\boldsymbol{j}({\boldsymbol{z}}(t))+\boldsymbol{{\dot{g}}_{Y}^{\text{in}}}(t)-\mathrm{d}_{t}\boldsymbol{{\mu}}_{Y}(t)\Big{)}\,.$
(51)

Systems are said to be autonomous when the chemostatting procedure maintains
the chemical potentials $\boldsymbol{{\mu}}_{Y}$ and the parameters
${\boldsymbol{e}}{}$ constant, i.e.,
$\mathrm{d}_{t}\boldsymbol{{\mu}}_{Y}(t)=0$ and
$\boldsymbol{{\dot{g}}_{Y}^{\text{in}}}(t)=0$. Otherwise, systems are said to
be nonautonomous. The matrices $\mathbb{W}(t)$ and $\mathbb{W}^{Y}(t)$ link
the value of the external currents to the concentrations of the internal
species that interact with the chemostatted ones. In the limit
${\boldsymbol{\gamma}}({\boldsymbol{z}},{\boldsymbol{e}}{})\to 1$, or
equivalently
$\boldsymbol{\nabla}_{\boldsymbol{z}}G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})\to
0$, $\boldsymbol{{I}^{Y}}(t)$ for ideal CRNs is recovered

$\begin{split}\boldsymbol{{I}^{Y}}(t)&=-\mathbb{S}{}^{Y}\boldsymbol{j}({\boldsymbol{z}}(t))+\Big{(}\mathbb{W}^{Y}(t)\Big{)}^{-1}\mathrm{d}_{t}\boldsymbol{\mu}^{\text{id}}_{Y}(t)\\\
&=-\mathbb{S}{}^{Y}\boldsymbol{j}({\boldsymbol{z}}(t))+\mathrm{d}_{t}\boldsymbol{y}(t)\,,\end{split}$
(52)

controlling directly the concentrations Rao and Esposito (2016) according to
the protocol $\mathrm{d}_{t}\boldsymbol{y}(t)$.

######  Example 7.

We specify the different terms in Eq. (51) for the CRN (9). According to Eq.
(50), $\mathbb{W}(t)$ is given by

$\scaleto{\mathbb{W}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{E_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{ES_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{EP_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil${-F\kappa}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil${\frac{RT}{[\ch{H^{+}_{{{{(1)}}}}}]}+F\kappa}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil${-F\kappa}$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil${\frac{RT}{[\ch{H^{+}_{{{{(2)}}}}}]}+F\kappa}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil${\frac{RT}{[\ch{S_{{{{(2)}}}}}]}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil${\frac{RT}{[\ch{P_{{{{(2)}}}}}]}}$\hfil\kern
5.0pt\cr}}}}\right)$}}}{60pt}\,,$ (53)

and hence $(\mathbb{W}^{Y}(t))^{-1}$ is given by

$(\mathbb{W}^{Y})^{-1}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil${\frac{[\ch{H^{+}_{{{{(1)}}}}}]}{RT+F\kappa[\ch{H^{+}_{{{{(1)}}}}}]}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil${\frac{[\ch{H^{+}_{{{{(2)}}}}}]}{RT+F\kappa[\ch{H^{+}_{{{{(2)}}}}}]}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil${\frac{[\ch{S_{{{{(2)}}}}}]}{RT}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil${\frac{[\ch{P_{{{{(2)}}}}}]}{RT}}$\hfil\kern
5.0pt\cr}}}}\right)$}}\,.$ (54)

Here, we did not write the explicit time dependence of the concentrations and
of the matrices for the sake of compactness. The vector
$\boldsymbol{{\dot{g}}_{Y}^{\text{in}}}(t)$ is specified as

$\boldsymbol{{\dot{g}}_{Y}^{\text{in}}}(t)=\hbox{}\vbox{\kern
0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$F\kappa\big{(}\mathrm{d}_{t}[\ch{C^{+}_{{(1)}}}](t)-\mathrm{d}_{t}[\ch{D^{-}_{{(1)}}}](t)\big{)}$\hfil\kern
5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$F\kappa\big{(}\mathrm{d}_{t}[\ch{C^{+}_{{(2)}}}](t)-\mathrm{d}_{t}[\ch{D^{-}_{{(1)}}}](t)\big{)}$\hfil\kern
5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr}}}}\right)$}}\,.$ (55)

∎

###  III.4 Decomposition of the Entropy Production Rate

We split the entropy production rate into different contributions. To do so we
adapt the decomposition of the entropy production for stochastic ideal CRNs
developed in Ref. Rao and Esposito, 2018a to deterministic non-ideal CRNs.
This is not straightforward and represents a major result of this work. The
readers who may want to skip the details of the derivation can go directly to
Eq. (68).

We start by expressing the entropy production rate of deterministic CRNs
Polettini and Esposito (2014); Rao and Esposito (2016),

$\dot{\Sigma}(t)=R\sum_{\rho\in\mathcal{R}}\left(j^{+\rho}({\boldsymbol{z}}(t))-j^{-\rho}({\boldsymbol{z}}(t))\right)\ln\frac{j^{+\rho}({\boldsymbol{z}}(t))}{j^{-\rho}({\boldsymbol{z}}(t))}\geq
0\,.$ (56)

Using the local detailed balance (38) and the specific expressions for the
free energies of reaction (39), we get

$T\dot{\Sigma}(t)=-\boldsymbol{{\mu}}({\boldsymbol{z}}(t))\cdot\mathbb{S}\boldsymbol{j}({\boldsymbol{z}}(t))\,.$
(57)

By explicitly considering the contributions of the $X$, ${Y_{f}}$ and
${Y_{p}}$ species, we obtain

$\small
T\dot{\Sigma}(t)=-\left(\boldsymbol{{\mu}}_{X}({\boldsymbol{z}}(t))\cdot\mathbb{S}{}^{X}+\boldsymbol{{\mu}}_{Y_{f}}(t)\cdot\mathbb{S}{}^{Y_{f}}+\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot\mathbb{S}{}^{Y_{p}}\right)\boldsymbol{j}({\boldsymbol{z}}(t))\,,$
(58)

with
$\boldsymbol{{\mu}}_{i}(\bullet)=(\dots,{\mu}_{\alpha}(\bullet),\dots)^{\intercal}_{\alpha\in
i}$ and $i=\\{X,{Y_{f}},{Y_{p}}\\}$. Note that we do not explicitly write the
${\boldsymbol{z}}(t)$ dependence of the chemical potentials of the $Y$ species
since they are externally controlled.

A part of the free energy externally exchanged through the currents
$\\{I^{\alpha}(t)\\}$ modifies the free energy of the system. Another part is
transferred through the system. To proceed with the entropy production
decomposition, we need to distinguish these two parts. To do so we exploit the
conservation laws. We consider

$\mathbb{L}^{{b}}\mathbb{S}=0=\mathbb{L}^{{b}}_{X}\mathbb{S}{}^{X}+\mathbb{L}^{{b}}_{{Y_{f}}}\mathbb{S}{}^{Y_{f}}+\mathbb{L}^{{b}}_{{Y_{p}}}\mathbb{S}{}^{Y_{p}}\,,$
(59)

and, recalling that the matrix $\mathbb{L}^{{b}}_{{Y_{p}}}$ can be inverted
(see Subs. II.3), we write

$\mathbb{S}{}^{Y_{p}}=-(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}(\mathbb{L}^{{b}}_{X}\mathbb{S}{}^{X}+\mathbb{L}^{{b}}_{{Y_{f}}}\mathbb{S}{}^{Y_{f}})\,.$
(60)

Equation (60) binds the net variation of the number of molecules for the
${Y_{p}}$ species to that of the $X$ and ${Y_{f}}$ species. We substitute this
result in Eq. (58):

$\begin{split}T\dot{\Sigma}(t)=&-\Big{(}\boldsymbol{{\mu}}_{X}({\boldsymbol{z}}(t))\cdot\mathbb{S}{}^{X}+\boldsymbol{{\mu}}_{Y_{f}}(t)\cdot\mathbb{S}{}^{Y_{f}}+\\\
&-\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}(\mathbb{L}^{{b}}_{X}\mathbb{S}{}^{X}+\mathbb{L}^{{b}}_{{Y_{f}}}\mathbb{S}{}^{Y_{f}})\Big{)}\boldsymbol{j}({\boldsymbol{z}}(t))\,.\end{split}$
(61)

By adding and subtracting
$\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}_{{Y_{p}}}\mathbb{S}{}^{Y_{p}}$,
the entropy production rate becomes

$T\dot{\Sigma}(t)=-\boldsymbol{\mathcal{F}}({\boldsymbol{z}}(t))\cdot\mathbb{S}\boldsymbol{j}({\boldsymbol{z}}(t))\,,$
(62)

where

$\boldsymbol{\mathcal{F}}({\boldsymbol{z}}(t))=\big{(}\boldsymbol{{\mu}}({\boldsymbol{z}}(t))\cdot\mathbb{1}-\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}\big{)}^{\intercal}$
(63)

with $\mathbb{1}$ the identity matrix. Note that we just introduced a term in
the expression of the entropy production rate, i.e.,
$\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}$,
whose net contribution vanishes since $\mathbb{L}^{{b}}\mathbb{S}=0$.
Introducing this term allows us to recognize the proper thermodynamic
potential of open non-ideal CRNs.

To do so, we define the following state function

$\Psi({\boldsymbol{z}}(t))=\boldsymbol{\mathcal{F}}({\boldsymbol{z}}(t))\cdot{\boldsymbol{z}}(t)\,,$
(64)

whose time derivative according to the rate equation Eq. (2) reads

$\begin{split}\mathrm{d}_{t}\Psi({\boldsymbol{z}}(t))&=\boldsymbol{\mathcal{F}}({\boldsymbol{z}}(t))\cdot\mathbb{S}\boldsymbol{j}({\boldsymbol{z}}(t))+\boldsymbol{\mathcal{F}}_{Y}(t)\cdot\boldsymbol{{I}^{Y}}(t)+\\\
&+RT\mathrm{d}_{t}\lVert{\boldsymbol{z}}(t)\rVert-\mathrm{d}_{t}\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}{\boldsymbol{z}}(t)+\\\
&+\mathrm{d}_{t}\big{(}\boldsymbol{\nabla}_{{\boldsymbol{z}}}G^{\text{in}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))\big{)}\cdot{\boldsymbol{z}}(t)\,,\end{split}$
(65)

where
$\boldsymbol{\mathcal{F}}_{Y}(t)=(\dots,\mathcal{F}_{\alpha}(t),\dots)_{\alpha\in
Y}$, $\lVert{\boldsymbol{z}}(t)\rVert=\sum_{\alpha}[\alpha](t)$. The time
dependence of ${\boldsymbol{e}}{}(t)$ accounts for the (possible) time
dependent manipulation of the parameters.

To procede with the decomposition, we first recognize that the term
$\boldsymbol{\mathcal{F}}_{Y}(t)\cdot\boldsymbol{{I}^{Y}}(t)$ simplifies to
$\boldsymbol{\mathcal{F}}_{Y_{f}}(t)\cdot\boldsymbol{{I}^{{Y_{f}}}}(t)$ since
$\mathcal{F}_{\alpha}=0$ for $\alpha\in{Y_{p}}$. We then notice that the term
$\mathrm{d}_{t}\big{(}\boldsymbol{\nabla}_{{\boldsymbol{z}}}G^{\text{in}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))\big{)}\cdot{\boldsymbol{z}}(t)$
can be written as

$\displaystyle\mathrm{d}_{t}\big{(}$
$\displaystyle\boldsymbol{\nabla}_{{\boldsymbol{z}}}G^{\text{in}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))\big{)}\cdot{\boldsymbol{z}}(t)=$
$\displaystyle=$
$\displaystyle\mathrm{d}_{t}\big{(}\boldsymbol{\nabla}_{{\boldsymbol{z}}}G^{\text{in}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))\cdot{\boldsymbol{z}}(t)\big{)}-\boldsymbol{\nabla}_{{\boldsymbol{z}}}G^{\text{in}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))\cdot\mathrm{d}_{t}{\boldsymbol{z}}(t)$
$\displaystyle=$
$\displaystyle\mathrm{d}_{t}\big{(}\boldsymbol{\nabla}_{{\boldsymbol{z}}}G^{\text{in}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))\cdot{\boldsymbol{z}}(t)\big{)}-\mathrm{d}_{t}G^{\text{in}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))+$
$\displaystyle+\boldsymbol{\nabla}_{{\boldsymbol{e}}{}}G^{\text{in}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))\cdot\mathrm{d}_{t}{\boldsymbol{e}}{}(t)$
(66)

by using the time derivative of the interaction free energy

$\begin{split}\mathrm{d}_{t}G^{\text{in}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))=&\boldsymbol{\nabla}_{{\boldsymbol{z}}}G^{\text{in}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))\cdot\mathrm{d}_{t}{\boldsymbol{z}}(t)+\\\
&+\boldsymbol{\nabla}_{{\boldsymbol{e}}{}}G^{\text{in}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))\cdot\mathrm{d}_{t}{\boldsymbol{e}}{}(t)\,,\end{split}$
(67)

with
$\boldsymbol{\nabla}_{{\boldsymbol{e}}{}}=(\dots,\partial/\partial{e}{}_{i},\dots)^{\intercal}$.

We now exploit Eq. (65) to express
$\boldsymbol{\mathcal{F}}({\boldsymbol{z}}(t))\cdot\mathbb{S}\boldsymbol{j}({\boldsymbol{z}}(t))$
as a function of the other terms and we substitute it in Eq. (62). By doing so
and collecting all the full time derivatives in a single term denoted
$\mathrm{d}_{t}\mathcal{G}$, we obtain the following decomposition of the
entropy production rate, which constitutes a central result of this work,

$T\dot{\Sigma}(t)=-\mathrm{d}_{t}\mathcal{G}({\boldsymbol{z}}(t))+\dot{w}_{\text{nc}}(t)+\dot{w}_{\text{driv}}(t)\geq
0\,.$ (68)

Here, we introduced the semigrand free energy of non-ideal systems as

$\begin{split}\mathcal{G}({\boldsymbol{z}})=&\boldsymbol{{\mu}}({\boldsymbol{z}})\cdot{\boldsymbol{z}}-RT\lVert{\boldsymbol{z}}\rVert-\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot\boldsymbol{m}({\boldsymbol{z}})+\\\
&+G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})-\boldsymbol{\nabla}_{{\boldsymbol{z}}}G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})\cdot{\boldsymbol{z}}\end{split}$
(69)

with the concentrations of the so-called moieties:

$\boldsymbol{m}({\boldsymbol{z}})=(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}{\boldsymbol{z}}\,.$
(70)

They represent parts of (or entire) molecules which are exchanged with the
environment. The semigrand free energy (69) can be rewritten in a more
appealing ways as

$\begin{split}\mathcal{G}({\boldsymbol{z}})&=\boldsymbol{\mu^{\text{id}}}({\boldsymbol{z}})\cdot{\boldsymbol{z}}-RT\lVert{\boldsymbol{z}}\rVert-\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot\boldsymbol{m}({\boldsymbol{z}})+G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})\\\
&=G({\boldsymbol{z}},{\boldsymbol{e}}{})-\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot\boldsymbol{m}({\boldsymbol{z}})\end{split}$
(71)

by using Eq. (34). The term
$\boldsymbol{\mu^{\text{id}}}({\boldsymbol{z}})\cdot{\boldsymbol{z}}-RT\lVert{\boldsymbol{z}}\rVert$
is the free energy of closed ideal CRNs, i.e.,
$G^{\text{id}}({\boldsymbol{z}})$, accounting for the chemical energy of every
species. The term $G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})$ is the
interaction free energy. Their sum gives the free energy
$G({\boldsymbol{z}},{\boldsymbol{e}}{})$ in Eq. (29). The term
$\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot\boldsymbol{m}({\boldsymbol{z}})$ is the
energetic contribution of the matter exchanged with the environment. It
vanishes for closed systems. Since $\mathcal{G}({\boldsymbol{z}})$ is a state
function, its time derivative vanishes at steady-state.

The driving work rate $\dot{w}_{\text{driv}}(t)$ can be split into the sum of
the chemical and the interaction driving work rates
$\dot{w}_{\text{driv}}(t)=\dot{w}_{\text{driv}}^{\text{ch}}(t)+\dot{w}_{\text{driv}}^{\text{in}}(t)$
that are are specified as

$\dot{w}_{\text{driv}}^{\text{ch}}(t)=-\mathrm{d}_{t}\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot\boldsymbol{m}({\boldsymbol{z}}(t))$
(72)

and

$\dot{w}_{\text{driv}}^{\text{in}}(t)=\boldsymbol{\nabla}_{{\boldsymbol{e}}{}}G^{\text{in}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))\cdot\mathrm{d}_{t}{\boldsymbol{e}}{}(t)\,,$
(73)

respectively. The chemical driving work (72) accounts for the time dependent
control of only the chemical potentials $\boldsymbol{{\mu}}_{Y_{p}}(t)$. The
interaction driving work (73) accounts for the time dependent control of the
interaction free energy through the time evolution of the parameters
${\boldsymbol{e}}{}(t)$. They represent the energetic cost of modifying the
equilibrium that is defined if $Y={Y_{p}}$ (as shown in Subs. II.5, open CRNs
with $Y={Y_{p}}$ have a well defined equilibrium state for every value of
${\boldsymbol{e}}{}$). These two terms vanish in autonomous systems.

The nonconservative work rate is given by

$\dot{w}_{\text{nc}}(t)=\boldsymbol{\mathcal{F}}_{{Y_{f}}}(t)\cdot\boldsymbol{{I}^{{Y_{f}}}}(t)\,.$
(74)

It quantifies the energetic cost to sustain fluxes of chemical species through
the CRN by means of the fundamental nonconservative forces

$\boldsymbol{\mathcal{F}}_{{Y_{f}}}(t)=\left(\boldsymbol{{\mu}}_{Y_{f}}(t)\cdot\mathbb{1}-\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}_{{Y_{f}}}\right)^{\intercal}\,.$
(75)

These are the forces keeping the system out of equilibrium. Indeed, they
vanish when the CRN is detailed balanced. This always happens when all the
chemostatted species break a conservation law $Y={Y_{p}}$ and, for this
reason, the CRNs are said unconditionally detailed balanced. Therefore, the
${Y_{f}}$ species are named force species. The ${Y_{p}}$ species are named
potential species since they fix the equilibrium conditions. Thus, the
splitting $Y={Y_{f}}\cup{Y_{p}}$ must be motivated by the different
thermodynamic role of the ${Y_{f}}$ species and ${Y_{p}}$ species. The
nonconservative work accounts for autonomous mechanisms keeping the system out
of equilibrium unlike the chemical driving work (72) and the interaction
driving work (73) . We stress that the expression of the fundamental force
(75) for non-ideal CRNs is analogous to that found for ideal CRNs Rao and
Esposito (2018a). The effects of the interactions that modify the value of the
nonconservative work (74) are hidden in the expressions of the exchanged
currents (51). Note also that the force (75) can vanish, and hence the CRN
becomes detailed balanced, even if $Y\neq{Y_{p}}$. Indeed,
$\boldsymbol{\mathcal{F}}_{{Y_{f}}}=0$ when the chemical potentials of the
${Y_{f}}$ species satisfy
$\boldsymbol{{\mu}}_{Y_{f}}^{\intercal}=\boldsymbol{{\mu}}_{Y_{p}}\cdot(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}_{{Y_{f}}}$.
Thus, the fundamental nonconservative forces are generated by the difference
between the actual values of chemical potentials of the ${Y_{f}}$ species,
i.e., $\boldsymbol{{\mu}}_{Y_{f}}$, and their equilibrium values if only the
${Y_{p}}$ species were chemostatted, i.e.,
$\boldsymbol{{\mu}}_{Y_{p}}\cdot(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}_{{Y_{f}}}$.

Finally, we notice that for autonomous
($\dot{w}_{\text{driv}}^{\text{ch}}(t)=\dot{w}_{\text{driv}}^{\text{in}}(t)=0$)
detailed balanced ($\dot{w}_{\text{nc}}(t)=0$) CRNs, the semigrand free energy
decreases monotonously in time, i.e.,
$\mathrm{d}_{t}\mathcal{G}(t)=-T\dot{\Sigma}(t)\leq 0$. The semigrand free
energy is dissipated.

######  Example 8.

We specify the moieties (70), the interaction driving work (73), the
nonconservative forces (75) and the nonconservative work (74) for the open CRN
(9). We start by giving an explicit expression for the matrix
$(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}$ (with $\mathbb{L}^{{b}}$
given in Eq. (16)):

$(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}=\hbox{}\vbox{\kern
0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{E_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{ES_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{EP_{{{(3)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{A^{-}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$-1$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr}}}}\right)$}}\,.$ (76)

The grey vertical lines mark the split of
$(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}$ into
$(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}_{X}$,
$(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}_{{Y_{f}}}$, and
$(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}_{{Y_{p}}}=\mathbb{1}$. By
using Eq. (76) in the definition of the concentrations of the moieties (70),
we obtain that

$\boldsymbol{m}({\boldsymbol{z}})=\begin{pmatrix}[\ch{H^{+}_{{{(1)}}}}]+[\ch{H^{+}_{{{(2)}}}}]\\\
[\ch{ES_{{{(3)}}}}]+[\ch{S_{{{(2)}}}}]-[\ch{H^{+}_{{{(1)}}}}]\\\
[\ch{EP_{{{(3)}}}}]+[\ch{P_{{{(2)}}}}]+[\ch{H^{+}_{{{(1)}}}}]\end{pmatrix}\,.$
(77)

Only the first moiety has a straightforward interpretation as the protons. The
others are linear combination of broken conserved quantities accounting for
the fragments of molecules exchange through the CRN.

We now consider the interaction driving work (73). By using the expression of
the interaction free energy for this example (35), we obtain

$\small\begin{split}\dot{w}_{\text{driv}}^{\text{in}}(t)=F\kappa\Big{(}&Q_{{(1)}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))\big{(}\mathrm{d}_{t}[\ch{C^{+}_{{{(1)}}}}](t)-\mathrm{d}_{t}[\ch{D^{-}_{{{(1)}}}}](t)\big{)}+\\\
+&Q_{{(2)}}({\boldsymbol{z}}(t),{\boldsymbol{e}}{}(t))\big{(}\mathrm{d}_{t}[\ch{C^{+}_{{{(2)}}}}](t)-\mathrm{d}_{t}[\ch{D^{-}_{{{(2)}}}}](t)\big{)}\Big{)}\,,\end{split}$
(78)

accounting for the variation of the semigrand free energy due to possible
nonautonomous processes changing the concentrations of the cations
$\ch{C^{+}}$ and the anions D^- in the two compartments.

Finally, we turn to the nonconservative forces (75). There is only one
nonconservative force for the open CRN (9),

$\begin{split}\mathcal{F}(t)&={\mu}_{\ch{H^{+}_{{{(1)}}}}}(t)-\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}\mathbb{L}^{{b}}_{{Y_{f}}}\\\
&=-(\mu^{\text{id}}_{\ch{P_{{{(2)}}}}}(t)+{\mu}_{\ch{H^{+}_{{{(2)}}}}}(t)-{\mu}_{\ch{H^{+}_{{{(1)}}}}}(t)-\mu^{\text{id}}_{\ch{S_{{{(2)}}}}}(t))\,,\end{split}$
(79)

since
$\boldsymbol{{\mu}}_{Y_{p}}(t)=({\mu}_{\ch{H^{+}_{{{(2)}}}}}(t),\mu^{\text{id}}_{\ch{S_{{{(2)}}}}}(t),\mu^{\text{id}}_{\ch{P_{{{(2)}}}}}(t))^{\intercal}$
(recall that the chemical potential of $\ch{S_{{{(2)}}}}$ and
$\ch{P_{{{(2)}}}}$ is the same as for ideal CRNs). It represents the
thermodynamic force (i.e., the free energy of reaction) of the following
effective chemical reaction

$\ch{S_{{{(2)}}}+H^{+}_{{{(1)}}}<=>H^{+}_{{{(2)}}}+P_{{{(2)}}}}$ (80)

between chemostatted species. Note that the role of the chemostatted species
is precisely that of keeping this effective reaction out of equilibrium by
imposing the thermodynamic force (79). The corresponding nonconservative work
reads

$\dot{w}_{\text{nc}}(t)=-(\mu^{\text{id}}_{\ch{P_{{{(2)}}}}}(t)+{\mu}_{\ch{H^{+}_{{{(2)}}}}}(t)-{\mu}_{\ch{H^{+}_{{{(1)}}}}}(t)-\mu^{\text{id}}_{\ch{S_{{{(2)}}}}}(t))I^{\ch{H^{+}_{{{(1)}}}}}(t)\,.$
(81)

We stress that, despite Eq. (81) resembling the corresponding expression for
ideal CRNs, the interactions affect the value the exchange current
$I^{\ch{H^{+}_{{{(1)}}}}}(t)$. ∎

###  III.5 Steady-State Entropy Production

When autonomous CRNs reach a nonequilibrium steady state
$\overline{{\boldsymbol{z}}}$, the entropy production rate (68) simplifies to

$T\overline{\dot{\Sigma}}=\overline{\dot{w}}_{\text{nc}}=\boldsymbol{\mathcal{F}}_{{Y_{f}}}\cdot\overline{\boldsymbol{I}}{}^{{Y_{f}}}\,,$
(82)

highlighting the physical role of the nonconservative work: it provides the
energy to balance dissipation and keep CRNs out of equilibrium.

An equivalent decomposition can be achieved using the cycle splitting of the
current vector (19) in Eq. (57):

$\begin{split}T\overline{\dot{\Sigma}}&=-\boldsymbol{{\mu}}\cdot\mathbb{S}\overline{\boldsymbol{j}}\\\
&=-\boldsymbol{{\mu}}\cdot\mathbb{S}\bigg{\\{}\sum_{\iota}\boldsymbol{c}_{\iota}\psi^{\iota}+\sum_{\epsilon}\boldsymbol{c}_{\epsilon}\psi^{\epsilon}\bigg{\\}}\\\
&=-\boldsymbol{{\mu}}_{Y}\cdot\hat{\mathbb{S}}{}^{Y}\boldsymbol{\psi}\end{split}$
(83)

where we used Eq. (17) and (18). The vector $\boldsymbol{\psi}$ collects the
coefficients $\\{\psi^{\epsilon}\\}$, i.e.,
$\boldsymbol{\psi}=(\dots,\psi^{\epsilon},\dots)^{\intercal}$. Each $\epsilon$
column of the effective stoichiometric matrix
$\hat{\mathbb{S}}{}^{Y}_{\epsilon}=\mathbb{S}{}^{Y}\boldsymbol{c}_{\epsilon}$
codifies the net stoichiometry of the chemostatted species along the
$\epsilon$ emergent cycle.

This alternative decomposition of the entropy production has two advantages.
First, it emphasizes the role of the emergent cycles in the steady-state
dissipation. They are out of equilibrium effective reactions interconverting
chemostatted species according to the stoichiometry codified in
$\hat{\mathbb{S}}{}^{Y}$. Second, this decomposition can be used to develop
thermodynamically consistent coarse-graining theories as it has been done for
ideal CRNs Wachtel  _et al._ (2018); Avanzini  _et al._ (2020b).

######  Example 9.

We specify the effective stoichiometric matrix $\hat{\mathbb{S}}{}^{Y}$ for
the CRN (9) given the emergent cycle (21):

${\hat{\mathbb{S}}{}^{Y}}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(1)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{H^{+}_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{S_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color@gray@fill{0.50}\ch{P_{{{(2)}}}}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr}}}}\right)$}}\,.$ (84)

Notice that the corresponding effective reaction is equal to the one we
inferred from the nonconservative force in Eq. (80). This equivalence does not
hold in general. ∎

###  III.6 Lower Bound of the Semigrand Free Energy

We now examine how the assumptions on the free energy
$G({\boldsymbol{z}},{\boldsymbol{e}}{})$ introduced in Subs. III.1 imply that
the semigrand free energy $\mathcal{G}({\boldsymbol{z}})$ (71) is lower
bounded and thus guarantee thermodynamic consistency.

Since $G({\boldsymbol{z}},{\boldsymbol{e}}{})$ is lower bounded and increases
superlinearly as concentrations go to infinity by assumption, also
$\mathcal{G}({\boldsymbol{z}})$ (see Eq. (71)) is lower bounded for finite
values of the chemical potential of the chemostatted species. Indeed, as
concentrations become infinitely large,
$-\boldsymbol{{\mu}}_{Y_{p}}\cdot\boldsymbol{m}$ can go to minus infinity
slower (linearly with ${\boldsymbol{z}}$) than
$G({\boldsymbol{z}},{\boldsymbol{e}}{})$ goes to infinity (superlinearly with
${\boldsymbol{z}}$). Thus, $\exists\mathcal{G}_{\text{min}}$ such that

$\mathcal{G}({\boldsymbol{z}})\geq\mathcal{G}_{\text{min}}\,,\text{ }\text{
}\text{ }\forall{\boldsymbol{z}}\,.$ (85)

This ensures thermodynamic consistency: autonomous detailed balanced systems
relax towards an equilibrium state. Given that i)
$\mathrm{d}_{t}\mathcal{G}({\boldsymbol{z}}(t))=-T\dot{\Sigma}(t)\leq 0$ (see
Subs. III.4), ii) $\dot{\Sigma}(t)=0$ only at equilibrium and iii)
$\mathcal{G}({\boldsymbol{z}}(t))$ cannot decrease indefinitely because of Eq.
(85), the system must reach an equilibrium state,
$\lim_{t\to+\infty}{\boldsymbol{z}}(t)={\boldsymbol{z}}_{\text{eq}}$, and
$\mathcal{G}({\boldsymbol{z}}(t))\geq\mathcal{G}({\boldsymbol{z}}_{\text{eq}})$,
$\forall t$.

The specific equilibrium state ${\boldsymbol{z}}_{\text{eq}}$ to which the
system relaxes depends on the initial conditions. In ideal CRNs, there is a
unique equilibrium state for each stoichiometric compatibility class, here
labeled $\Omega\left(\big{\\{}L^{\lambda_{u}}\big{\\}}\right)$ (see Subs.
II.5). In non-ideal CRNs, there might be multiple equilibria for every
stoichiometric compatibility class. Thus, we introduce
$\Omega\left({\boldsymbol{z}}_{\text{eq}}\right)\subseteq\Omega\left(\big{\\{}L^{\lambda_{u}}\big{\\}}\right)$
as the manifold in the concentration space such that any state in
$\Omega\left({\boldsymbol{z}}_{\text{eq}}\right)$ will relax towards the
specific equilibrium ${\boldsymbol{z}}_{\text{eq}}$. This implies that

$\mathcal{G}({\boldsymbol{z}})\geq\mathcal{G}({\boldsymbol{z}}_{\text{eq}})\,,\text{
}\text{ }\text{
}\forall{\boldsymbol{z}}\in\Omega\left({\boldsymbol{z}}_{\text{eq}}\right)\,.$
(86)

Let us stress two more points. First, the bound (86) does not provide any
constraint for the values of the semigrand free energy evaluated on different
manifolds. Consequently, it is possible that
$\exists{\boldsymbol{z}}\notin\Omega\left({\boldsymbol{z}}_{\text{eq}}\right)$
such that
$\mathcal{G}({\boldsymbol{z}})<\mathcal{G}({\boldsymbol{z}}_{\text{eq}})$.
Second, $\mathcal{G}_{\text{min}}$ must be an equilibrium value of the
semigrand free energy. Indeed, if there was a nonequilibrium state
${\boldsymbol{z}}$ such that
$\mathcal{G}({\boldsymbol{z}})=\mathcal{G}_{\text{min}}$, the relaxation of
${\boldsymbol{z}}$ would necessarily correspond to a decrease of the semigrand
free energy which is not possible because of Eq. (85).

We now turn to nonautonomous and/or nondetailed balanced CRNs. The expression
of the semigrand free energy (71) is exactly the same as for the corresponding
autonomous detailed balanced CRNs. Thus, Eq. (86) still holds with
${\boldsymbol{z}}_{\text{eq}}$ now the equilibrium state of the corresponding
detailed balanced CRNs with $Y={Y_{p}}$. However, the evolution of
${\boldsymbol{z}}(t)$ is constrained on a specific stoichiometric
compatibility class $\Omega\left(\big{\\{}L^{\lambda_{u}}\big{\\}}\right)$,
but not (necessarily) on a specific manifold
$\Omega\left({\boldsymbol{z}}_{\text{eq}}\right)$. This means that different
bounds (86) can be identified during the same nondetailed balanced dynamics.

In conclusion, the semigrand free energy is always lower bounded by
$\mathcal{G}({\boldsymbol{z}}_{\text{eq}})$ where
${\boldsymbol{z}}_{\text{eq}}$ is the equilibrium state to which the system
would relax if it was detailed balanced. This, together with
$\mathrm{d}_{t}\mathcal{G}({\boldsymbol{z}}(t))=-T\dot{\Sigma}(t)\leq 0$ for
autonomous detailed balanced CRNs, makes the semigrand free energy defined in
Eq. (69) the proper thermodynamic potential of open non-ideal CRNs. In other
words,
$\mathcal{G}({\boldsymbol{z}}(t))-\mathcal{G}({\boldsymbol{z}}_{\text{eq}})$
is a Lyapunov function.

###  III.7 Relative Entropy

We show here that, because of the interactions, the difference
$\mathcal{G}({\boldsymbol{z}}(t))-\mathcal{G}({\boldsymbol{z}}_{\text{eq}})$
cannot be solely expressed as a relative entropy. This is a significant
difference with respect to previous studies on the thermodynamics of ideal
CRNs Rao and Esposito (2016, 2018a) and ideal reaction-diffusion systems
Falasco  _et al._ (2018); Avanzini  _et al._ (2019).

To this aim, we notice that, because of the local detailed balance (38), the
equilibrium chemical potentials satisfy

$\boldsymbol{{\mu}}({\boldsymbol{z}}_{\text{eq}})\cdot\mathbb{S}=0\,.$ (87)

Thus, $\boldsymbol{{\mu}}({\boldsymbol{z}}_{\text{eq}})$ can be written as a
linear combination of conservation laws (11),
$\boldsymbol{{\mu}}({\boldsymbol{z}}_{\text{eq}})=\sum_{\lambda}\mathpzc
f_{\lambda}\boldsymbol{\ell}^{\lambda}$. Recalling that the conservation laws
split into unbroken and broken conservation laws in open CRNs, we obtain

$\boldsymbol{{\mu}}({\boldsymbol{z}}_{\text{eq}})\cdot{\boldsymbol{z}}_{\text{eq}}=\boldsymbol{{\mu}}({\boldsymbol{z}}_{\text{eq}})\cdot{\boldsymbol{z}}(t)-\sum_{\lambda_{b}}\mathpzc
f_{\lambda_{b}}\boldsymbol{\ell}^{\lambda_{b}}\cdot{\boldsymbol{z}}(t)+\sum_{\lambda_{b}}\mathpzc
f_{\lambda_{b}}\boldsymbol{\ell}^{\lambda_{b}}\cdot{\boldsymbol{z}}_{\text{eq}}\,.$
(88)

By recalling that the unbroken conservation laws are those with null entries
for the chemostatted species, hence
${\mu}_{\alpha}({\boldsymbol{z}}_{\text{eq}})=\sum_{\lambda_{b}}\mathpzc
f_{\lambda_{b}}\ell^{\lambda_{b}}_{\alpha}$ for $\alpha\in Y$, and the
expressions for the concentrations of the moieties (70), we verify that

$\displaystyle\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot\boldsymbol{m}({\boldsymbol{z}}(t))=\sum_{\lambda_{b}}\mathpzc
f_{\lambda_{b}}\boldsymbol{\ell}^{\lambda_{b}}\cdot{\boldsymbol{z}}(t)\,,$
(89)
$\displaystyle\boldsymbol{{\mu}}_{Y_{p}}(t)\cdot\boldsymbol{m}({\boldsymbol{z}}_{\text{eq}})=\sum_{\lambda_{b}}\mathpzc
f_{\lambda_{b}}\boldsymbol{\ell}^{\lambda_{b}}\cdot{\boldsymbol{z}}_{\text{eq}}\,.$
(90)

By using Eq. (88), Eq. (89), and Eq. (90), we find that

$\begin{split}\mathcal{G}({\boldsymbol{z}}(t))-\mathcal{G}({\boldsymbol{z}}_{\text{eq}})=&\left(\boldsymbol{{\mu}}({\boldsymbol{z}}(t))-\boldsymbol{{\mu}}({\boldsymbol{z}}_{\text{eq}})\right)\cdot{\boldsymbol{z}}(t)+\\\
&-RT(\lVert{\boldsymbol{z}}(t)\rVert-\lVert{\boldsymbol{z}}_{\text{eq}}\rVert)+\\\
&+G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})-\boldsymbol{\nabla}_{{\boldsymbol{z}}}G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})\cdot{\boldsymbol{z}}+\\\
&-G^{\text{in}}({\boldsymbol{z}}_{\text{eq}},{\boldsymbol{e}}{})+\boldsymbol{\nabla}_{{\boldsymbol{z}}}G^{\text{in}}({\boldsymbol{z}}_{\text{eq}},{\boldsymbol{e}}{})\cdot{\boldsymbol{z}}_{\text{eq}}\,.\\\
\end{split}$ (91)

Finally, by using Eq. (34) to split the ideal and the non-ideal contribution
of the chemical potential, we obtain

${\begin{split}\mathcal{G}({\boldsymbol{z}}(t))-\mathcal{G}({\boldsymbol{z}}_{\text{eq}})=&RT\mathcal{L}({\boldsymbol{z}}(t)\parallel{\boldsymbol{z}}_{\text{eq}})+\\\
&+G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})-G^{\text{in}}({\boldsymbol{z}}_{\text{eq}},{\boldsymbol{e}}{})+\\\
&-\boldsymbol{\nabla}_{{\boldsymbol{z}}}G^{\text{in}}({\boldsymbol{z}}_{\text{eq}},{\boldsymbol{e}}{})\cdot({\boldsymbol{z}}-{\boldsymbol{z}}_{\text{eq}})\,,\end{split}}$
(92)

where we introduced the relative entropy for non-normalized concentration
distributions

$\mathcal{L}(\boldsymbol{a}\parallel\boldsymbol{b})=\sum_{i}a_{i}\ln\left(\frac{a_{i}}{b_{i}}\right)-(a_{i}-b_{i})\geq
0\,.$ (93)

Equation (92) proves that
$\mathcal{G}({\boldsymbol{z}}(t))-\mathcal{G}({\boldsymbol{z}}_{\text{eq}})$
cannot be written only as a relative entropy, but additional terms emerge from
the interactions. It is not granted that these additional terms are always
positive unless $G^{\text{in}}({\boldsymbol{z}},{\boldsymbol{e}}{})$ is a
convex function of ${\boldsymbol{z}}$. Nevertheless, if
${\boldsymbol{z}}(t)\in\Omega\left({\boldsymbol{z}}_{\text{eq}}\right)$,
$\mathcal{G}({\boldsymbol{z}}(t))-\mathcal{G}({\boldsymbol{z}}_{\text{eq}})\geq
0$ as proved in Subs. III.6.

######  Example 10.

The difference
$\mathcal{G}({\boldsymbol{z}})-\mathcal{G}({\boldsymbol{z}}_{\text{eq}})$ in
the CRN (9), where ${\boldsymbol{z}}_{\text{eq}}$ is the equilibrium to which
the system would relax if only the ${Y_{p}}$ species where chemostatted, can
be written as

$\begin{split}\mathcal{G}({\boldsymbol{z}})-\mathcal{G}({\boldsymbol{z}}_{\text{eq}})=&RT\mathcal{L}({\boldsymbol{z}}(t)\parallel{\boldsymbol{z}}_{\text{eq}})+\\\
&+\frac{F\kappa}{2}\big{(}Q_{{(1)}}({\boldsymbol{z}},{\boldsymbol{e}}{})-Q_{{(1)}}({\boldsymbol{z}}_{\text{eq}},{\boldsymbol{e}}{})\big{)}^{2}+\\\
&+\frac{F\kappa}{2}\big{(}Q_{{(2)}}({\boldsymbol{z}},{\boldsymbol{e}}{})-Q_{{(2)}}({\boldsymbol{z}}_{\text{eq}},{\boldsymbol{e}}{})\big{)}^{2}\,.\end{split}$
(94)

Since the interaction free energy (35) is a convex function, the terms in Eq.
(94) emerging from the electrostatic interactions are always positive, i.e.,
$\big{(}Q_{(i)}({\boldsymbol{z}},{\boldsymbol{e}}{})-Q_{(i)}({\boldsymbol{z}}_{\text{eq}},{\boldsymbol{e}}{})\big{)}^{2}\geq
0$. ∎

######  Example 11.

We examine the difference
$\mathcal{G}({\boldsymbol{z}})-\mathcal{G}({\boldsymbol{z}}_{\text{eq}})$ for
the following closed CRN:

$\ch{B<=>[{+1}][{-1}]A^{-}+C^{+}}\,,$ (95)

when the electrostatic interactions are described within Debye-Hückel theory
McQuarrie (1975):

$G^{\text{in}}({\boldsymbol{z}})=-\frac{4F\kappa}{h^{3}}\bigg{\\{}\frac{h^{2}}{2}\mathcal{I}({\boldsymbol{z}})-h\sqrt{\mathcal{I}({\boldsymbol{z}})}+\ln\Big{[}h\sqrt{\mathcal{I}({\boldsymbol{z}})}+1\Big{]}\bigg{\\}}\,,$
(96)

where $\mathcal{I}({\boldsymbol{z}})$ is the ionic strength

$\mathcal{I}({\boldsymbol{z}})=\frac{1}{2}([\ch{A^{-}}]+[\ch{C^{+}}])\geq
0\,,$ (97)

and $\kappa$ and $h$ are positive constants depending on the specific system
(see Ref. McQuarrie, 1975 for their explicit expressions). In this case
$\mathcal{G}({\boldsymbol{z}})-\mathcal{G}({\boldsymbol{z}}_{\text{eq}})$
becomes

$\begin{split}\mathcal{G}({\boldsymbol{z}})-\mathcal{G}({\boldsymbol{z}}_{\text{eq}})=&RT\mathcal{L}({\boldsymbol{z}}(t)\parallel{\boldsymbol{z}}_{\text{eq}})+\\\
&+\big{(}G^{\text{in}}({\boldsymbol{z}})-G^{\text{in}}({\boldsymbol{z}}_{\text{eq}})\big{)}+\\\
&+2F\kappa\big{[}\mathcal{I}({\boldsymbol{z}})-\mathcal{I}({\boldsymbol{z}}_{\text{eq}})\big{]}\frac{\sqrt{\mathcal{I}({\boldsymbol{z}}_{\text{eq}})}}{1+h\sqrt{\mathcal{I}({\boldsymbol{z}}_{\text{eq}})}}\,.\end{split}$
(98)

Here the terms emerging from the electrostatic interactions may be negative
because the interaction free energy provided by the Debye-Hückel theory is not
a convex function. ∎

##  IV Conclusions

In this work, we developed a thermodynamic theory of non-ideal CRNs. We used
activity coefficients to account for interactions within a mean-field
approach. Our framework is general and holds for any (thermodynamic
consistent) expressions of the interactions. Exploiting the local detailed
balance assumption, we showed that the dynamics of non-ideal CRNs does not
follow mass-action kinetics unlike for ideal-CRNs. Generalizing some results
developed for ideal-CRNs, we used conservation laws and cycles to derive
physically meaningful decompositions of the entropy production rate. The first
(68) is derived exploiting the conservation laws and allowed us to determine
the proper thermodynamic potential and the forces driving non-ideal CRNs out
of equilibrium. The second (83) is derived exploiting the cycles and allowed
us to express the steady-state entropy production using only the chemical
potential of the chemostatted species and the currents along the emergent
cycles. We proved that the thermodynamic potential acts as a Lyapunov function
in detailed balanced CRNs.

Crucially, we showed that the thermodynamic structure of non-ideal CRNs is the
same as in ideal CRNs. In particular, the conservation laws and cycles
determine whether systems are detailed balance or not and define the forces
(75) maintaining nonequilibrium regimes in the same way for ideal and non-
ideal CRNs.

These results are fundamental to characterize processes of energy transduction
in biosystems and electrochemical systems where interacting ions are
ubiquitous. As a perspective, one could exploit the same strategy we followed
in this work to specialize the stochastic thermodynamics of CRNs to
interacting systems. We leave these points to future investigations.

##  V Acknowledgments

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

## References

  * Prigogine (1961) I. Prigogine,  _Introduction to Thermodynamics of Irreversible Processes_ , 2nd ed. (Interscience Publishers, New York, 1961).
  * de Groot and Mazur (1984) S. R. de Groot and P. Mazur,  _Non-Equilibrium Thermodynamics_ (Dover, 1984).
  * Hill (1977) Terrell L. Hill,  _Free Energy Transduction in Biology_ (Academic Press, New York, 1977).
  * Schnakenberg (1976) J. Schnakenberg, “Network theory of microscopic and macroscopic behavior of master equation systems,” Rev. Mod. Phys. 48, 571–585 (1976).
  * Jiu-li  _et al._ (1984) Luo Jiu-li, C. Van den Broeck, and G. Nicolis, “Stability criteria and fluctuations around nonequilibrium states,” Z. Phys. B Condensed Matter 56, 165–170 (1984).
  * Mou  _et al._ (1986) Chung Yuan Mou, Jiu-li Luo, and Gregoire Nicolis, “Stochastic thermodynamics of nonequilibrium steady states in chemical reaction systems,” J. Chem. Phys. 84, 7011–7017 (1986).
  * Ross  _et al._ (1988) John Ross, Katharine L. C. Hunt, and Paul M. Hunt, “Thermodynamics far from equilibrium: Reactions with multiple stationary states,” J. Chem. Phys. 88, 2719–2729 (1988).
  * Jarzynski (2011) Christopher Jarzynski, “Equalities and inequalities: Irreversibility and the second law of thermodynamics at the nanoscale,” Annu. Rev. Condens. Matter Phys. 2, 329–351 (2011).
  * Seifert (2012) Udo Seifert, “Stochastic thermodynamics, fluctuation theorems and molecular machines,” Rep. Prog. Phys. 75, 126001 (2012).
  * den Broeck and Esposito (2015) C. Van den Broeck and M. Esposito, “Ensemble and trajectory thermodynamics: A brief introduction,” Physica A 418, 6 – 16 (2015).
  * Kondepudi and Prigogine (2015) Dilip Kondepudi and Ilya Prigogine,  _Modern Thermodynamics. From heat engines to dissipative structures_ , 2nd ed. (Wiley, 2015).
  * Gaspard (2004) Pierre Gaspard, “Fluctuation theorem for nonequilibrium reactions,” The Journal of Chemical Physics 120, 8898–8905 (2004).
  * Schmiedl and Seifert (2007) Tim Schmiedl and Udo Seifert, “Stochastic thermodynamics of chemical reaction networks,” J. Chem. Phys. 126, 044101 (2007).
  * Rao and Esposito (2018a) Riccardo Rao and Massimiliano Esposito, “Conservation laws and work fluctuation relations in chemical reaction networks,” J. Chem. Phys. 149, 245101 (2018a).
  * Qian and Beard (2005) Hong Qian and Daniel A. Beard, “Thermodynamics of stoichiometric biochemical networks in living systems far from equilibrium,” Biophys. Chem. 114, 213 – 220 (2005).
  * Polettini and Esposito (2014) Matteo Polettini and Massimiliano Esposito, “Irreversible thermodynamics of open chemical networks. I. emergent cycles and broken conservation laws,” J. Chem. Phys. 141, 024117 (2014).
  * Rao and Esposito (2016) Riccardo Rao and Massimiliano Esposito, “Nonequilibrium thermodynamics of chemical reaction networks: Wisdom from stochastic thermodynamics,” Phys. Rev. X 6, 041064 (2016).
  * Falasco  _et al._ (2018) Gianmaria Falasco, Riccardo Rao, and Massimiliano Esposito, “Information thermodynamics of turing patterns,” Phys. Rev. Lett. 121, 108301 (2018).
  * Avanzini  _et al._ (2019) Francesco Avanzini, Gianmaria Falasco, and Massimiliano Esposito, “Thermodynamics of chemical waves,” J. Chem. Phys. 151, 234103 (2019).
  * Avanzini  _et al._ (2020a) Francesco Avanzini, Gianmaria Falasco, and Massimiliano Esposito, “Chemical cloaking,” Phys. Rev. E 101, 060102 (2020a).
  * Andrieux and Gaspard (2008) David Andrieux and Pierre Gaspard, “Nonequilibrium generation of information in copolymerization processes,” Proc. Natl. Acad. Sci. USA 105, 9516–9521 (2008).
  * Blokhuis and Lacoste (2017) Alex Blokhuis and David Lacoste, “Length and sequence relaxation of copolymers under recombination reactions,” The Journal of Chemical Physics 147, 094905 (2017).
  * Gaspard (2020a) Pierre Gaspard, “Template-directed growth of copolymers,” Chaos 30, 043114 (2020a).
  * Rao  _et al._ (2015) Riccardo Rao, David Lacoste, and Massimiliano Esposito, “Glucans monomer-exchange dynamics as an open chemical network,” J. Chem. Phys. 143, 244903 (2015).
  * Gaspard (2020b) Pierre Gaspard, “Stochastic approach to entropy production in chemical chaos,” Chaos 30, 113103 (2020b).
  * Yoshimura and Ito (2020) Kohei Yoshimura and Sosuke Ito, “Information geometric inequalities of chemical thermodynamics,” (2020), arXiv:2005.08444, arXiv:2005.08444 [cond-mat.stat-mech] .
  * Fritz  _et al._ (2020) Jonas H. Fritz, Basile Nguyen, and Udo Seifert, “Stochastic thermodynamics of chemical reactions coupled to finite reservoirs: A case study for the brusselator,” J. Chem. Phys. 152, 235101 (2020).
  * Blokhuis  _et al._ (2018) Alex Blokhuis, David Lacoste, and Pierre Gaspard, “Reaction kinetics in open reactors and serial transfers between closed reactors,” J. Chem. Phys. 148, 144902 (2018).
  * Hołyst and Poniewierski (2012) Robert Hołyst and Andrzej Poniewierski,  _Thermodynamics for Chemists, Physicists and Engineers_ (Springer, 2012).
  * Boudart (1968) Michel Boudart,  _Kinetics of Chemical Processes_ (Prentice-Hall, 1968).
  * Laidler (1987) K. J. Laidler,  _Chemical Kinetics_ (Harper Collins Publishers, 1987).
  * Butt (1999) John B. Butt,  _Reaction kinetics and reactor design_ (Marcel Dekker, Inc., 1999).
  * Voet and Voet (2010) D. Voet and J. G. Voet,  _Biochemistry_ (J. Wiley & Sons, 2010).
  * Bard and Faulkner (2001) Allen J. Bard and Larry R. Faulkner,  _Electrochemical Methods: Fundamentals and Applications_ , 2nd ed. (Wiley, 2001).
  * G. (1993) Svehla G., “Nomenclature of kinetic methods of analysis (iupac recommendations 1993),”  _Pure and Applied Chemistry_ ,  65, 2291 (1993).
  * Pekař (2005) Miloslav Pekař, “Thermodynamics and foundations of mass-action kinetics,” Prog. React. Kinet. Mech. 30, 3–113 (2005).
  * Huxley  _et al._ (2000) A. F. Huxley, R. M. Simmons, Kazuhiko Kinosita, Ryohei Yasuda, Hiroyuki Noji, and Kengo Adachi, “A rotary molecular motor that can work at near 100% efficiency,” Philos. Trans. R Soc. Lond. B Biol. Sci. 355, 473–489 (2000).
  * Rao and Esposito (2018b) Riccardo Rao and Massimiliano Esposito, “Conservation laws shape dissipation,” New J. Phys. 20, 023007 (2018b).
  * Schuster and Schuster (1989) Stefan Schuster and Ronny Schuster, “A generalization of wegscheider’s condition. implications for properties of steady states and for quasi-steady-state approximation,” J. Math. Chem. 3, 25–42 (1989).
  * Horn and Jackson (1972) F. Horn and R. Jackson, “General mass action kinetics,” Arch. Ration. Mech. Anal. 47, 81–116 (1972).
  * McQuarrie (1975) D.A. McQuarrie,  _Statistical Mechanics_ (Harper & Row, 1975).
  * Chipot and Pohorille (2007) Christophe Chipot and Andrew Pohorille,  _Free Energy Calculations_ (Springer, Berlin, 2007).
  * Lelièvre  _et al._ (2010) Tony Lelièvre, Mathias Rousset, and Gabriel Stoltz,  _Free Energy Computations: A Mathematical Perspective_ (Imperial College Press, 2010).
  * Holyst and Poniewierski (2012) Robert Holyst and Andrzej Poniewierski,  _Thermodynamics for Chemists, Physicists and Engineers_ (Springer, Dordrecht (NL), 2012).
  * Magnus and Keizer (1997) G. Magnus and J. Keizer, “Minimal model of beta-cell mitochondrial $\text{{C}a}^{2+}$ handling,” Am. J. Physiol. Cell Physiol. 273, C717–C733 (1997).
  * Magnus and Keizer (1998) Gerhard Magnus and Joel Keizer, “Model of $\beta$-cell mitochondrial calcium handling and electrical activity. I. cytoplasmic variables,” Am. J. Physiol. Cell Physiol. 274, C1158–C1173 (1998).
  * Fall and Keizer (2001) Christopher P. Fall and Joel E. Keizer, “Mitochondrial modulation of intracellular $\text{{C}a}^{2+}$ signaling,” J. Theor. Biol. 210, 151 – 165 (2001).
  * Wachtel  _et al._ (2018) Artur Wachtel, Riccardo Rao, and Massimiliano Esposito, “Thermodynamically consistent coarse graining of biocatalysts beyond michaelis–menten,” New J. Phys. 20, 042002 (2018).
  * Avanzini  _et al._ (2020b) Francesco Avanzini, Gianmaria Falasco, and Massimiliano Esposito, “Thermodynamics of non-elementary chemical reaction networks,” New J. Phys. 22, 093040 (2020b).