# Thermodynamics of Concentration vs Flux Control in Chemical Reaction
Networks

Francesco Avanzini<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

We investigate the thermodynamic implications of two control mechanisms of
open chemical reaction networks. The first controls the concentrations of the
species that are exchanged with the surroundings, while the other controls the
exchange fluxes. We show that the two mechanisms can be mapped one into the
other and that the thermodynamic theories usually developed in the framework
of concentration control can be applied to flux control as well. This implies
that the thermodynamic potential and the fundamental forces driving chemical
reaction networks out of equilibrium can be identified in the same way for
both mechanisms. By analyzing the dynamics and thermodynamics of a simple
enzymatic model we also show that, while the two mechanisms are equivalent at
steady state, the flux control may lead to fundamentally different regimes
where systems achieve stationary growth.

##  I Introduction

Stochastic thermodynamics Sekimoto (2010); Jarzynski (2011); Seifert (2012);
den Broeck and Esposito (2015) provides a rigorous ground to characterize the
energetics of systems driven arbitrarily far from equilibrium. It has been
formulated for open chemical reaction networks (CRNs) undergoing a stochastic
Gaspard (2004); Schmiedl and Seifert (2007); Rao and Esposito (2018a) or a
deterministic dynamics Qian and Beard (2005); Polettini and Esposito (2014);
Rao and Esposito (2016). Crucially, the latter formulation is equivalent to
the former in macroscopic limit Ge and Qian (2016a, b); Rao and Esposito
(2016, 2018a), where fluctuations become negligible and the stochastic
dynamics converges to the deterministic one Kurtz (1970, 1971, 1972), and to
Gibbs’ chemical thermodynamics at equilibrium. In these frameworks, CRNs are
driven out of equilibrium by the exchanges of some chemical species with the
surroundings which are in general represented in terms of infinitely large
(but also finite Fritz  _et al._ (2020)) reservoirs, called chemostats,
controlling the concentrations of the exchanged species. This theory has been
used to quantify the energetic cost for creating patterns Falasco  _et al._
(2018) and waves Avanzini  _et al._ (2019), and sustaining the growth process
of macromolecules, like copolymers Andrieux and Gaspard (2008); Blokhuis and
Lacoste (2017); Gaspard (2020) and biomolecules Rao  _et al._ (2015). Its
connections to information geometry Yoshimura and Ito (2021a), as well as to
thermodynamic uncertainty relations and speed limits Yoshimura and Ito
(2021b), have been recently investigated. Furthermore, the topological
properties of CRNs, and in particular the conservation laws, have been
exploited to decompose the entropy production rate in such a way as to
identify the proper thermodynamic potential of open CRNs as well as the
fundamental forces and chemostats breaking the detailed balance condition and
driving CRNs out of equilibrium Rao and Esposito (2018a); Avanzini  _et al._
(2021).

Controlling the concentrations of some species via chemostats is not the only
way the surroundings can affect CRNs. They can also control the exchange
fluxes. For instance, if the mitochondrial metabolism is analyzed as an open
CRN, its surroundings correspond to the cytosolic processes which do not
control the concentrations of some species inside the mitochondrion, but
instead they determine its exchange fluxes Magnus and Keizer (1997, 1998);
Fall and Keizer (2001). Furthermore, in industrial continuous-flow stirred
tank reactors Aris (1989) as well as in biotechnological continuous cultures
von Stockar (2013), systems are fed by reactants at a constant rate while some
other species are continuously extracted with a controlled outflow. Hence,
recent studies started considering flux control mechanisms to drive CRNs out
of equilibrium rather than concentration control via chemostats. Some of these
studies focused on continuous-flow stirred tank reactors which have been
analyzed from a kinetic Craciun and Feinberg (2005) as well as a thermodynamic
Blokhuis  _et al._ (2018) standpoint. This formulation has for example been
used to study the emergence of homochirality Laurent  _et al._ (2021a, b) and
reformulate the thermodynamic evolution theorem in modern terms Hochberg and
Ribó (2020). Other studies considered CRNs under both concentration and flux
control, and investigated how the sensitivity of CRNs to changes in the
extraction fluxes can be adjusted by nonequilibrium futile cycles Qian and
Beard (2006).

In this work, we provide a systematic comparison between open CRNs under
concentration control and flux control from a thermodynamic standpoint. We
focus on deterministic CRNs whose dynamical description is summarized in Sec.
II. In Sec. III, we discuss concentration control and its thermodynamic
characterization, in particular the decomposition of the entropy production
rate in terms of thermodynamic potential and fundamental forces. We turn to
flux control in Sec. IV where we show that it can be mapped into a specific
protocol of concentration control. This proves that the same thermodynamic
theories developed for CRNs under concentration control can be used for CRNs
under flux control, including the decomposition of the entropy production
rate. In Sec. V, we investigate the different thermodynamic properties of the
two control mechanisms by considering a model system of an enzymatic reaction.
We focus on how the entropy production rate, the thermodynamic potential and
the work done by the fundamental forces evolve depending on the control
mechanism. General conclusions are drawn in Sec. VI where we also discuss the
implications of our results which should be particularly relevant to
characterize the energetics of different growth regimes in autocatalytic
systems.

##  II Dynamics of Open CRNs

CRNs are considered here as ideal dilute solutions composed of chemical
species ${\boldsymbol{\alpha}}=(\dots,\alpha,\dots)^{\intercal}$ which
interconvert via chemical reactions. In the macroscopic limit, the evolution
of the abundances of the chemical species, expressed in terms of
concentrations ${\boldsymbol{z}}=(\dots,[\alpha],\dots)^{\intercal}$, follows
a deterministic rate equation:

$\frac{\mathrm{d}}{\mathrm{d}t}{\boldsymbol{z}}={\mathbb{S}}{\boldsymbol{j}}+{\boldsymbol{I}}\,.$
(1)

The first term on the r.h.s. of Eq. (1), i.e., ${\mathbb{S}}{\boldsymbol{j}}$,
accounts for the variation of the concentrations due to the chemical
reactions. Each column ${\boldsymbol{S}}_{\rho}$ of the stoichiometric matrix
${\mathbb{S}}$ encodes the net variations in terms of number of molecules of
each species undergoing (the elementary G. (1993)) reaction $\rho$,

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

where $\boldsymbol{\nu}_{+\rho}$ (resp. $\boldsymbol{\nu}_{-\rho}$) is the
vector of the stoichiometric coefficients of the forward (resp. backward)
reaction. Thus,
${\boldsymbol{S}}_{\rho}=\boldsymbol{\nu}_{-\rho}-\boldsymbol{\nu}_{+\rho}$.
Each entry $j_{\rho}$ of the vector
${\boldsymbol{j}}=(\dots,j_{\rho},\dots)^{\intercal}$ specifies the net
current of every reaction $\rho$ as the difference between the forward and
backward one, i.e., $j_{\rho}=j_{+\rho}-j_{-\rho}$, which satisfy mass-action
kinetics:

$j_{\pm\rho}=k_{\pm\rho}{\boldsymbol{z}}^{\boldsymbol{\nu}_{\pm\rho}}\,,$ (3)

with $k_{\pm\rho}$ the kinetic constants of the forward/backward reaction and
${\boldsymbol{a}}^{\boldsymbol{b}}=\prod_{i}a_{i}^{b_{i}}$.

The second term on the r.h.s. of Eq. (1) is the exchange vector
${\boldsymbol{I}}=(\dots,{I}_{\alpha},\dots)^{\intercal}$ accounting for net
matter flux of each species between the CRN and the surroundings independently
of the control mechanism. It vanishes if CRNs are closed. In open CRNs, we
split the chemical species into the internal species $\alpha\in X$, namely
those that are not exchanged (i.e., ${I}_{\alpha}=0$), and the exchanged
species $\alpha\in Y$. By applying the same splitting to the concentration
vector ${\boldsymbol{z}}=(\boldsymbol{x},\boldsymbol{y})$ and the
stoichiometric matrix,

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

the rate equation (1) can be specialized as

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{x}={\mathbb{S}}^{X}{\boldsymbol{j}}\,,$
(5)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{y}={\mathbb{S}}^{Y}{\boldsymbol{j}}+{\boldsymbol{I}}^{Y}\,,$
(6)

with ${\boldsymbol{I}}^{Y}=(\dots,{I}_{\alpha},\dots)^{\intercal}_{\alpha\in
Y}$ collecting all the nonw null entries of ${\boldsymbol{I}}$.

Closed CRNs must be detailed balanced, namely an equilibrium steady state
${\boldsymbol{z}}_{\text{eq}}$ such that
${\boldsymbol{j}}({\boldsymbol{z}}_{\text{eq}})=0$ always exists. In open
CRNs, the detailed balanced condition can be broken because of the exchanges
of the chemical species.

##  III Concentration Control

A first way in which the matter exchanges with the surroundings can affect the
dynamics of CRNs is by controlling the concentrations of the $Y$ species.
Thus, $\boldsymbol{y}$ cease to be dynamical variables, unlike the
concentrations of the internal species $\boldsymbol{x}$, and become externally
controlled parameters, while Eq. (6) becomes merely a definition of the
exchange fluxes ${\boldsymbol{I}}^{Y}$.

We specialize this concentration control in two classes. The first one
maintains the concentrations $\boldsymbol{y}$ constant, i.e.,
${\mathrm{d}\boldsymbol{y}}/{\mathrm{d}t}=0$, and the corresponding exchange
fluxes are given by

${\boldsymbol{I}}^{Y}=-{\mathbb{S}}^{Y}{\boldsymbol{j}}\,.$ (7)

In this case, CRNs are said to be autonomous since the control parameters
$\boldsymbol{y}$ have fixed values. The second one forces the concentrations
$\boldsymbol{y}$ to follow a specific protocol, i.e.,
${\boldsymbol{y}}={\boldsymbol{\pi}}$, and the corresponding exchange fluxes
must satisfy

${\boldsymbol{I}}^{Y}=\frac{\mathrm{d}}{\mathrm{d}t}{\boldsymbol{\pi}}-{\mathbb{S}}^{Y}{\boldsymbol{j}}\,.$
(8)

In this case, CRNs are said to be nonautonomous.

The concentration control is the most common one in the framework of
nonequilibrium thermodynamics. Physically, it corresponds to coupling a CRN
with external reservoirs called chemostats as in Fig. 1.

Figure 1: Pictorial illustration of the CRN (24) when subjected to
(autonomous) concentration control via chemostats fixing the concentrations of
F, S, W, P.

Each chemostat exchanges a specific $Y$ species according to the process

$\alpha\ch{<=>[\end{equation}][}$

] α_c that is assumed to always be at equilibrium. From a thermodynamic
standpoint, this implies that the chemical potential $\mu_{\alpha}$ of each
$\alpha\in Y$ species, namely its contribution to the free energy, is
controlled by the chemostat, i.e, $\mu_{\alpha}=\mu_{\alpha_{c}}$. The
chemical potentials are given by

$\mu_{\alpha}=\mu^{\circ}_{\alpha}+RT\ln[\alpha]\,,$ (9)

where $\mu^{\circ}_{\alpha}$ is the standard chemical potential of the species
$\alpha$, $R$ is the gas constant, and $T$ is the temperature fixed by the
solvent. Consequently, the chemostats control the concentrations of the $Y$
species,

$[\alpha]=\exp\Big{(}\frac{\mu_{\alpha_{c}}-\mu^{\circ}_{\alpha}}{RT}\Big{)}\text{
}\text{ }\text{ }\forall\alpha\in Y\,,$ (10)

that are said to be chemostatted. Depending on whether the chemical potentials
$\mu_{\alpha_{c}}$ are constant in time or not, the concentrations
$\boldsymbol{y}$ have fixed values or not. Note that in previous studies
Schmiedl and Seifert (2007); Polettini and Esposito (2014); Falasco  _et al._
(2019) the $Y$ species are directly regarded as chemostats within the CRN.
These two interpretations are formally equivalent, but the former has the
advantage of being easier to compare to the flux control as we will see.
Within this framework, the entropy production rate of elementary reactions,

$\dot{\Sigma}=R\sum_{\rho}\big{(}j_{+\rho}-j_{-\rho})\ln\frac{j_{+\rho}}{j_{-\rho}},$
(11)

has been decomposed Rao and Esposito (2018a); Avanzini  _et al._ (2021) as

$T\dot{\Sigma}=-\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{G}+\dot{w}_{\text{nc}}+\dot{w}_{\text{driv}}\,,$
(12)

where the semi-grand free energy $\mathcal{G}$ is the proper thermodynamic
potential of open CRNs; $\dot{w}_{\text{nc}}$ is the nonconservative work rate
accounting for the energetic cost of breaking detailed balance;
$\dot{w}_{\text{driv}}$ is the driving work rate quantifying the energetic
cost of changing the equilibrium state to which CRNs would relax if they were
detailed balance. The derivation of Eq. (12) is discussed in detail in Ref.
Rao and Esposito (2018a) for stochastic and in Ref. Avanzini  _et al._ (2021)
for deterministic CRNs. Here, we recap the general strategy and we provide the
explicit expressions of $\mathcal{G}$, $\dot{w}_{\text{nc}}$, and
$\dot{w}_{\text{driv}}$. To do so, we have to introduce the conservation laws
and the splitting of the chemostatted species $Y$ into potential ${Y_{p}}$ and
force species ${Y_{f}}$. The conservation laws Polettini and Esposito (2014);
Rao and Esposito (2016) are linearly independent vectors
$\\{{\boldsymbol{\ell}}^{\lambda}\\}$ in the cokernel of the stoichiometric
matrix,

${\boldsymbol{\ell}}^{\lambda}\cdot{\mathbb{S}}=0\,.$ (13)

Indeed, all the scalars
$L^{\lambda}\equiv{\boldsymbol{\ell}}^{\lambda}\cdot{\boldsymbol{z}}$ would be
conserved quantities if CRNs were closed:
${\mathrm{d}L}/{\mathrm{d}t}={\boldsymbol{\ell}}^{\lambda}\cdot{\mathbb{S}}{\boldsymbol{j}}=0$.
These scalars (or their linear combinations) quantify the concentrations of
the so-called moieties, i.e., parts of (or entire) molecules that remain
intact in all the reactions. Note that the total mass must be conserved in
closed CRNs and thus the set $\\{{\boldsymbol{\ell}}^{\lambda}\\}$ is never
empty. In open CRNs, the set of conservation laws can be split into the
unbroken conservation laws $\\{{\boldsymbol{\ell}}^{\lambda_{u}}\\}$ and the
broken conservation laws $\\{{\boldsymbol{\ell}}^{\lambda_{b}}\\}$. The former
are the largest subset of conservation laws that can be written with null
entries for the $Y$ species, i.e, $\ell^{\lambda_{u}}_{\alpha}=0$
$\forall\alpha\in Y$. The scalar quantities
$L^{\lambda_{u}}={\boldsymbol{\ell}}^{\lambda_{u}}\cdot{\boldsymbol{z}}$ are
thus conserved even if CRNs are opened:

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

The later are the other conservation laws and the scalar quantities
$L^{\lambda_{b}}={\boldsymbol{\ell}}^{\lambda_{b}}\cdot{\boldsymbol{z}}$ are
in general not conserved when CRNs are opened, i.e.,

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

because the corresponding moieties are exchanged with the chemostats. Note
that the total mass is not conserved in open CRNs and thus the set
$\\{{\boldsymbol{\ell}}^{{\lambda_{b}}}\\}$ is never empty. Chemostatting a
species does not always break a conservation law Rao and Esposito (2016,
2018a, 2018b); Falasco  _et al._ (2018). The potential species ${Y_{p}}$ are
those that break the conservation laws when chemostatted, while the force
species are the other species ${Y_{f}}=Y\setminus{Y_{p}}$. On the one hand,
every time a ${Y_{p}}$ species is chemostatted a new moiety is exchanged
between the CRN and the chemostats. The concentrations of the exchanged
moieties are expressed in terms of linear combinations 111 The choice of the
conservation laws (13) is not unique, and different choices identify different
moieties. The linear combination in Eq. (16) selects a particular set of
conservation laws and, correspondingly, of moieties such that each moiety is
exchanged with only one chemostat as long as no ${Y_{f}}$ species are
chemostatted: $\mathrm{d}_{t}\boldsymbol{m}={\boldsymbol{I}}_{Y_{p}}$ of the
broken conserved quantities (15) according to

$\boldsymbol{m}=(\mathbb{L}^{{b}}_{{Y_{p}}})^{-1}{\boldsymbol{L}}^{b}$ (16)

where ${\boldsymbol{L}}^{b}=(\dots,L^{\lambda_{b}},\dots)^{\intercal}$ and
$\mathbb{L}^{{b}}_{{Y_{p}}}$ is the (invertible) submatrix for the ${Y_{p}}$
species of the matrix $\mathbb{L}^{{b}}$ whose rows are the broken
conservation laws (i.e., $\mathbb{L}^{{b}}$ and $\mathbb{L}^{{b}}_{{Y_{p}}}$
have $\\{\ell^{\lambda_{b}}_{\alpha}\\}$ and
$\\{\ell^{\lambda_{b}}_{\alpha}\\}_{\alpha\in{Y_{p}}}$ as entries,
respectively). On the other hand, when a ${Y_{f}}$ species is chemostatted,
the corresponding chemostat exchanges a moiety with the CRN that is already
exchanged with another one thus establishing a flux of a moiety between two
(or more) chemostats 222 When also the ${Y_{f}}$ species are chemostatted, the
concentration of each moiety in Eq. (16) depends on the flux with the
corresponding ${Y_{p}}$ chemostat and on the fluxes with the ${Y_{f}}$
chemostats that exchange the same moiety according to
$\mathrm{d}_{t}\boldsymbol{m}={\boldsymbol{I}}_{Y_{p}}+\big{(}\mathbb{L}^{{b}}_{{Y_{p}}}\big{)}^{-1}\mathbb{L}^{{b}}_{{Y_{f}}}{\boldsymbol{I}}_{Y_{f}}$
(with $\mathbb{L}^{{b}}_{{Y_{f}}}$ the submatrix for the ${Y_{f}}$ species of
$\mathbb{L}^{{b}}$). The different role played by the ${Y_{p}}$ and ${Y_{f}}$
specie provides a rigorous ground to identify the thermodynamic potential of
open CRNs and to decompose the entropy production rate (12). As in equilibrium
thermodynamics when passing from the canonical to the grand canonical
ensemble, the thermodynamic potential $\mathcal{G}$ is obtained from the Gibbs
free energy $G$ by eliminating the energetic contributions of the matter
exchanged with the particle reservoir that in this context corresponds to the
chemostats. The latter accounts for the concentrations of the moieties (16)
times the chemical potential of the ${Y_{p}}$ species. Thus,

$\mathcal{G}=\underbrace{{\boldsymbol{\mu}}\cdot{\boldsymbol{z}}-RT\lVert{\boldsymbol{z}}\rVert}_{=\text{
}G}-{\boldsymbol{\mu}}_{Y_{p}}\cdot\boldsymbol{m}\,,$ (17)

with ${\boldsymbol{\mu}}=(\dots,\mu_{\alpha},\dots)^{\intercal}$,
${\boldsymbol{\mu}}_{Y_{p}}=(\dots,\mu_{\alpha},\dots)^{\intercal}_{\alpha\in{Y_{p}}}$,
and $\lVert{\boldsymbol{z}}\rVert=\sum_{\alpha}[\alpha]$. Exploiting the
conservation law, one can further verify that $\mathcal{G}$ is lower bounded
by its equilibrium value, i.e. $\mathcal{G}\geq\mathcal{G}_{\text{eq}}$. The
nonconservative work rate,

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

quantifies the energetic cost of sustaining fluxes of moieties between
chemostats by means of the nonconservative force
${\boldsymbol{\mathcal{F}}}{}_{Y_{f}}=\big{(}{\boldsymbol{\mu}}_{Y_{f}}\cdot\mathbb{1}-{\boldsymbol{\mu}}_{Y_{p}}\cdot\big{(}\mathbb{L}^{{b}}_{{Y_{p}}}\big{)}^{-1}\mathbb{L}^{{b}}_{{Y_{f}}}\big{)}^{\intercal}$
(with $\mathbb{L}^{{b}}_{{Y_{f}}}$ the submatrix for the ${Y_{f}}$ species of
$\mathbb{L}^{{b}}$, i.e., the matrix with entries
$\\{\ell^{\lambda_{b}}_{\alpha}\\}_{\alpha\in{Y_{f}}}$, and $\mathbb{1}$ the
identity matrix). When there are no ${Y_{f}}$ species, each moiety is
exchanged with a single ${Y_{p}}$ chemostat and there are no fluxes between
chemostats. Thus, $\dot{w}_{\text{nc}}$ vanishes. If there are ${Y_{f}}$
species, $\dot{w}_{\text{nc}}$ can still vanish if

${\boldsymbol{\mu}}_{Y_{f}}=\big{(}{\boldsymbol{\mu}}_{Y_{p}}\cdot\big{(}\mathbb{L}^{{b}}_{{Y_{p}}}\big{)}^{-1}\mathbb{L}^{{b}}_{{Y_{f}}}\big{)}^{\intercal}\,.$
(19)

The driving work rate,

$\dot{w}_{\text{driv}}=-\Big{(}\frac{\mathrm{d}}{\mathrm{d}t}{\boldsymbol{\mu}}_{Y_{p}}\Big{)}\cdot\boldsymbol{m}\,,$
(20)

accounts for the time dependent manipulation of the chemical potentials of the
moieties, and so of their concentrations, via the ${Y_{p}}$ chemostats. Since
the chemical potential of a species (9) is univocally determined by its
concentration, the driving work rate obviously vanishes if
${\mathrm{d}\boldsymbol{y}}/{\mathrm{d}t}=0$, namely under an autonomous
concentration control. We consider now three simple cases to illustrate the
physical meaning of the thermodynamic quantities given in Eq. (17), (18) and
(20). First, when either no species are exchanged or only the concentrations
of the ${Y_{p}}$ species are maintained constant by the chemostats, Eq. (12)
simplifies to ${\mathrm{d}\mathcal{G}}/{\mathrm{d}t}=-T\dot{\Sigma}\leq 0$.
This, together with the lower bound $\mathcal{G}\geq\mathcal{G}_{\text{eq}}$,
implies that the CRN eventually reaches an equilibrium state after a
nonequilibrium transient dynamics in which the semigrand free energy
$\mathcal{G}$ is dissipated, i.e., it is detailed balanced. The particular
equilibrium state to which the CRN relaxes depends on the concentrations of
the ${Y_{p}}$ species. Second, when the concentrations of the ${Y_{p}}$ are
driven, the CRN is still detailed balanced, but the equilibrium state to which
it would relax changes in time. If the driving is fast enough, the CRN never
reaches equilibrium and is thus maintained in a nonequilibrium transient
dynamics. In this case, Eq. (12) becomes
${\mathrm{d}\mathcal{G}}/{\mathrm{d}t}=-T\dot{\Sigma}+\dot{w}_{\text{driv}}$
and $\dot{w}_{\text{driv}}$ accounts for the energetic cost of driving the
equilibrium state. For this reason $\dot{w}_{\text{driv}}$ depends only on the
${Y_{p}}$ species. Third, when the concentrations of the ${Y_{p}}$ species are
maintained constant while also the ${Y_{f}}$ species are exchanged, detailed
balance is broken: the CRN is maintained out of equilibrium both in the
transient dynamics and at steady state. Indeed, Eq. (12) becomes
${\mathrm{d}\mathcal{G}}/{\mathrm{d}t}=-T\dot{\Sigma}+\dot{w}_{\text{nc}}$
which implies that the dissipation cannot vanish at steady state since
$T\overline{\dot{\Sigma}}={\overline{\dot{w}}_{\text{nc}}}$. Thus,
$\dot{w}_{\text{nc}}$ quantifies the energetic cost of breaking detailed
balance and ${\boldsymbol{\mathcal{F}}}{}_{Y_{f}}$ is the nonconservative
force preventing the CRN from reaching an equilibrium state.

##  IV Flux Control

A second way in which the matter exchanges with the surroundings can affect
the dynamics of CRNs is by controlling the exchange fluxes
${\boldsymbol{I}}^{Y}$. The concentrations $\boldsymbol{y}$ are not controlled
parameters, but they are still dynamical variables whose evolution follows Eq.
(6). We can specialize this flux control in two classes. In the first one, the
exchange fluxes ${\boldsymbol{I}}^{Y}$ are either constant or depend on the
concentrations ${\boldsymbol{z}}$. We refer to this regime as autonomous since
the control parameters ${\boldsymbol{I}}^{Y}$ are not externally manipulated.
In the second one, the exchange fluxes ${\boldsymbol{I}}^{Y}$ can be forced to
follow an externally imposed protocol. We refer to this regime as
nonautonomous. The flux control physically corresponds to coupling a CRN with
other external processes whose net effects on the concentrations of the $Y$
species are quantified by the exchange fluxes ${\boldsymbol{I}}^{Y}$:

$\alpha\ch{<->[{I}_{\alpha}][\end{equation}]}$

… , where we used the symbol $\ch{<->}$ to represent the fact that external
processes can either inject or extract the species $\alpha\in Y$. For example,
some species can be injected in the CRN at a certain constant rate and
extracted with an outflow that is proportional to their concentrations as in
Fig. 2.

Figure 2: Pictorial illustration of the CRN (24) when subjected to
(autonomous) flux control. The species S and F are continuously injected at
the constant rates ${I}_{\ch{S}}$ and ${I}_{\ch{F}}$, respectively. The
species P and W are extracted at concentration-dependent rates
${I}_{\ch{P}}=-k_{\ch{P}}[\ch{P}]$ and ${I}_{\ch{W}}=-k_{\ch{W}}[\ch{W}]$,
respectively.

This particular strategy is equivalent to a continuous-flow stirred tank
reactor when applied to all the species Blokhuis  _et al._ (2018).
Furthermore, chemostating can be considered as particular external processes
(III) that fix the concentrations $\boldsymbol{y}$. It is natural to wonder
whether the decomposition of the entropy production rate (12) derived for CRNs
coupled to chemostats can still be applied to this different framework or not.
We argue here that it can be by mapping the flux control into a nonautonomous
concentration control. Consider as first case a generic CRN where the exchange
fluxes ${\boldsymbol{I}}^{Y}_{1}$ are controlled (here we use the subscript 1
to stress that these are the exchange fluxes of the first case). We do not
make any assumption on the expressions of the entries of
${\boldsymbol{I}}^{Y}_{1}$: they might be constant, or depend on the
concentrations ${\boldsymbol{z}}$ or follow an externally imposed protocol.
What matters is that the concentrations of the chemical species at a generic
time $t$,

${\boldsymbol{z}}_{1}=(\boldsymbol{x}{}_{1},\boldsymbol{y}{}_{1})\,,$ (21)

are univocally determined by solving Eq. (1) (or equivalently Eq. (5) and (6))
for given expressions of ${\boldsymbol{I}}^{Y}_{1}$ and a given initial
condition ${\boldsymbol{z}}_{0}=(\boldsymbol{x}{}_{0},\boldsymbol{y}{}_{0})$.
Consider as second case the same CRN with the same initial condition, but
where the concentrations of the $Y$ species are controlled. If the
concentrations $\boldsymbol{y}{}_{2}$ follow an external protocol
${\boldsymbol{\pi}}_{2}$ which is equivalent to their evolution in the first
case, that is

${\boldsymbol{\pi}}_{2}=\boldsymbol{y}{}_{1}\text{ }\text{ }\text{ }\forall
t\,,$ (22)

then also the evolution of the exchange fluxes ${\boldsymbol{I}}^{Y}_{2}$ and
the internal concentrations will necessarily be the same for the whole
dynamics:

${\boldsymbol{I}}^{Y}_{2}={\boldsymbol{I}}^{Y}_{1}\,,\text{ }\text{ }\text{
}\boldsymbol{x}{}_{2}=\boldsymbol{x}{}_{1}\,,\text{ }\text{ }\text{ }\forall
t\,.$ (23)

This implies that the CRN cannot know whether it is subjected to a
concentration control or a flux control: it only experiences the net effects
of the coupling with the surroundings, i.e., the exchange fluxes
${\boldsymbol{I}}^{Y}$. Therefore, the free energy balance between the CRN and
the surroundings must be exactly the same in the two cases and, consequently,
the decomposition of the entropy production in Eq. (12) still holds. The
difference between the two ways of controlling CRNs regards only their
surroundings. In particular, we emphasize that while chemostats are external
reservoirs at equilibrium, the external processes can be out of equilibrium.
Furthermore, the derivation Rao and Esposito (2018a); Avanzini  _et al._
(2021) of the semigrand free energy $\mathcal{G}$, the nonconservative work
rate $\dot{w}_{\text{nc}}$ and the driving work rate $\dot{w}_{\text{driv}}$
(summarized in Sec. III) is based only on which moieties (broken conserved
quantities) are exchanged with the surroundings. It does not require any
assumption on the kind of processes breaking the conservation laws. Thus,
$\mathcal{G}$, $\dot{w}_{\text{nc}}$ and $\dot{w}_{\text{driv}}$ are given by
Eq. (17), (18), and (20), respectively, whether CNRs are subjected to
concentration control or flux control, and their physical meaning is exactly
the same in both cases.

##  V Model System

We now consider the following CRN,

$\begin{split}\ch{F+E&<=>[${+1}$][${-1}$]EF}\\\
\ch{EF&<=>[${+2}$][${-2}$]EW}\\\ \ch{EW&<=>[${+3}$][${-3}$]E+W}\\\
\ch{S+EF&<=>[${+4}$][${-4}$]E^{*}}\\\
\ch{E^{*}&<=>[${+5}$][${-5}$]EW+P}\end{split}$ (24)

representing the active interconversion of the substrate S into the product P
catalyzed by the enzyme E and powered by the interconversion of the fuel F
into the waste W. The species EF (resp. EW) is the complex enzyme-fuel (resp.
enzyme-waste) while E^* is the complex binding both the fuel and the
substrate. Given the stoichiometric matrix of the CRN (24),

${{\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}1$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}2$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}3$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}4$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}5$\hfil\kern
5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\ch{E}$\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&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}\ch{EF}$\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$-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}\ch{EW}$\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&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\ch{E^{*}}$\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&5.0pt\hfil$-1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\ch{P}$\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$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\ch{W}$\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&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\ch{S}$\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&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\ch{F}$\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$0$\hfil\kern 5.0pt\cr}}}}\right)$}}\,,$ (25)

we identify the conservation laws

${\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}\ch{E}$\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}\ch{EF}$\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}\ch{EW}$\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}\ch{E^{*}}$\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}\ch{P}$\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}\ch{W}$\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}\ch{S}$\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}\ch{F}$\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}\ch{E}$\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}\ch{EF}$\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}\ch{EW}$\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}\ch{E^{*}}$\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}\ch{P}$\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}\ch{W}$\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}\ch{S}$\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}\ch{F}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr}}}}\right)$}}\,,\text{ }\text{ }\text{
}\text{ }{\boldsymbol{\ell}}^{\ch{F}}=\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}\ch{E}$\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}\ch{EF}$\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}\ch{EW}$\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}\ch{E^{*}}$\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}\ch{P}$\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}\ch{W}$\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}\ch{S}$\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}\ch{F}$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr}}}}\right)$}}\,,$ (26)

and the concentrations of the corresponding moieties, i.e., the enzyme, the
substrate and the fuel moiety:
$L^{\ch{E}}=[\ch{E}]+[\ch{EF}]+[\ch{EW}]+[\ch{E^{*}}]$,
$L^{\ch{S}}=[\ch{E^{*}}]+[\ch{P}]+[\ch{S}]$,
$L^{\ch{F}}=[\ch{EF}]+[\ch{EW}]+[\ch{E^{*}}]+[\ch{W}]+[\ch{F}]$. In the
following, we examine how the different thermodynamic quantities introduced in
the decomposition of the entropy production rate (12) change according to the
control mechanism in some prototypical cases. We use $k_{+3}$,
$k_{+3}/k_{+1}$, and $RT(k_{+3})^{2}/k_{+1}$ as units of measure for time,
concentration, and thermodynamic quantities, respectively. We assume that
$k_{\pm 1}=k_{-2}=k_{\pm 3}=k_{\pm 4}=k_{\pm 5}=2$, $k_{+2}=2e$,
$\mu^{\circ}_{\ch{E}}=\mu^{\circ}_{\ch{W}}=\mu^{\circ}_{\ch{S}}=1$,
$\mu^{\circ}_{\ch{EW}}=\mu^{\circ}_{\ch{P}}=\mu^{\circ}_{\ch{F}}=2$,
$\mu^{\circ}_{\ch{EF}}=3$, $\mu^{\circ}_{\ch{E^{*}}}=4$, and the initial
condition is given by $[\ch{E}]=0.15$, $[\ch{EF}]=0.1$, $[\ch{EW}]=0.1$,
$[\ch{E^{*}}]=0.1$, $[\ch{P}]=0.5$, $[\ch{W}]=0.5$, $[\ch{S}]=1$,
$[\ch{F}]=1$. The solution of the rate equation (1), where the currents
satisfy mass-action kinetics (3), has been numerically computed using the
standard SciPy ODE solver and then employed to determine the thermodynamic
quantities using their definitions.

###  V.1 Detailed Balanced Case

We start by considering the case where only the species
$Y=\\{\ch{S},\ch{F}\\}$ are exchanged with the surroundings. The conservation
laws ${\boldsymbol{\ell}}^{\ch{S}}$ and ${\boldsymbol{\ell}}^{\ch{F}}$ are
broken and, consequently, i) both the species $\ch{S}$ and $\ch{F}$ are
potential species ${Y_{p}}$ and ii) the CRN is (unconditionally) detailed
balanced. According to Eq. (16), the concentrations of the exchanged moieties
are given by the concentrations of the substrate and the fuel moiety,

$\boldsymbol{m}=\begin{pmatrix}L^{\ch{S}}\\\ L^{\ch{F}}\\\ \end{pmatrix}$ (27)

since $\mathbb{L}^{{b}}_{{Y_{p}}}$ is the identity matrix in this particular
case where $\mathbb{L}^{{b}}$, whose rows are the broken conservation laws, is
speficied 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}\ch{E}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\ch{EF}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\ch{EW}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\ch{E^{*}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\ch{P}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\ch{W}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\ch{S}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\ch{F}$\hfil\kern
5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\boldsymbol{{\boldsymbol{\ell}}}^{\ch{S}}$\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&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\vrule$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\vrule$\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}\boldsymbol{{\boldsymbol{\ell}}}^{\ch{F}}$\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$1$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\vrule$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$\color[rgb]{0.50,0.50,0.50}\vrule$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt\cr}}}}\right)$}}\,.$ (28)

Thus, the semigrand free energy (17) is given by

$\mathcal{G}=G-\mu_{\ch{S}}L^{\ch{S}}-\mu_{\ch{F}}L^{\ch{F}}\,,$ (29)

the driving work (20) is specified as

$\dot{w}_{\text{driv}}=-\Big{(}\frac{\mathrm{d}}{\mathrm{d}t}\mu_{\ch{S}}\Big{)}L^{\ch{S}}-\Big{(}\frac{\mathrm{d}}{\mathrm{d}t}\mu_{\ch{F}}\Big{)}L^{\ch{F}}\,.$
(30)

while the nonconservative work (18) vanishes because there are no ${Y_{f}}$
species.

Figure 3:  Exchange fluxes (a) and thermodynamic quantities (b) of the CRN
(24) when the species S and F are subjected to an autonomous concentration
control. Concentrations of P, W, S, and F (c), and thermodynamic quantities
(d) of the CRN (24) when the species S and F are subjected to an autonomous
flux control with ${I}_{\ch{S}}=0.025$ and ${I}_{\ch{P}}=0.05$.

By applying an autonomous concentration control on the $Y$ species, the CRN
will eventually reach an equilibrium steady state. In the transient, the
exchange fluxes ${I}_{\ch{S}}$ and ${I}_{\ch{F}}$ balance the variation of the
concentration $[\ch{S}]$ and $[\ch{F}]$ due to the chemical reactions (see
Fig. 3a). When equilibrium is reached, all the reaction currents vanish and
hence also the exchange fluxes (see Eq. (7) and Fig. 3a). From a thermodynamic
stand point, the dynamics corresponds to a pure dissipation of the semigrand
free energy as shown in Fig. 3b. Indeed, no driving work is performed and Eq.
(12) simplifies to ${\mathrm{d}\mathcal{G}}/{\mathrm{d}t}=-T\dot{\Sigma}\leq
0$. On the other hand, by imposing an autonomous flux control that
continuously injects $\ch{S}$ and $\ch{F}$ in the CRN at constant rate, namely
${I}_{\ch{S}}$ and ${I}_{\ch{F}}$ have constant positive values, the
concentrations $[\ch{P}]$, $[\ch{W}]$, $[\ch{S}]$, and $[\ch{F}]$ continuously
grow as shown in Fig. 3c. The concentrations $[\ch{E}]$, $[\ch{EF}]$,
$[\ch{EW}]$, and $[\ch{E^{*}}]$ are instead bounded since the concentration of
the enzyme moiety $L^{\ch{E}}$ is still conserved (in other words,
${\boldsymbol{\ell}}^{\ch{E}}$ is an unbroken conservation law). In this case,
the dynamics never reaches a steady state. From a thermodynamic standpoint,
the control mechanism performs a driving work which maintains the CRN out of
equilibrium (see. Fig. 3d). In the long time limit, the driving work
continuously extracts free energy at (almost) constant rate
($\dot{w}_{\text{driv}}\simeq-0.08$) which is provided by the constant
decrease of the semigrand free energy
(${\mathrm{d}\mathcal{G}}/{\mathrm{d}t}<\dot{w}_{\text{driv}}<0$). In
parallel, the dissipation is very small
($0<\dot{\Sigma}\ll|{\mathrm{d}\mathcal{G}}/{\mathrm{d}t}|$) which implies
that the CRN is constantly close to an equilibrium state that is changing in
time because of the driving work. Note that this is independent of the
magnitude of ${I}_{\ch{S}}$ and ${I}_{\ch{F}}$. In the long time limit, the
concentrations of $\ch{S}$ and $\ch{F}$ are large enough that the exchange
processes act as a slow perturbation thus allowing the CRN to stay close to
equilibrium.

###  V.2 Nondetailed Balanced Case

We consider now the case where also the species P and W are exchanged, i.e.,
$Y=\\{\ch{P},\ch{W},\ch{S},\ch{F}\\}$. Since no further conservation laws are
broken, namely ${\boldsymbol{\ell}}^{\ch{E}}$ is still an unbroken
conservation law, the species P and W are force species ${Y_{f}}$. Thus, the
expressions of $\mathcal{G}$ and $\dot{w}_{\text{driv}}$ are still the same as
in Eq. (29) and (30), respectively, while the nonconservative work is now
given by

$\dot{w}_{\text{nc}}=(\mu_{\ch{P}}-\mu_{\ch{S}}){I}_{\ch{P}}+(\mu_{\ch{W}}-\mu_{\ch{F}}){I}_{\ch{W}}\,,$
(31)

since also the matrix $\mathbb{L}^{{b}}_{{Y_{f}}}$ is the identity matrix in
this particular case given $\mathbb{L}^{{b}}$ in Eq. (28).

Figure 4: Exchange fluxes (a) and thermodynamic quantities (b) of the CRN (24)
when the species P, W, S and F are subjected to an autonomous concentration
control. Exchange fluxes (c) and thermodynamic quantities (d) of the CRN (24)
when the species P, W, S and F are subjected to an autonomous flux control
with ${I}_{\ch{P}}=-0.05[\ch{P}]$, ${I}_{\ch{W}}=-0.1[\ch{W}]$,
${I}_{\ch{S}}=0.025$, and ${I}_{\ch{F}}=0.05$. Concentrations of P, W, S, and
F (e), and thermodynamic quantities (f) of the CRN (24) when the species P, W,
S and F are subjected to an autonomous flux control with
${I}_{\ch{P}}=-1[\ch{P}]$, ${I}_{\ch{W}}=-2[\ch{W}]$, ${I}_{\ch{S}}=1$, and
${I}_{\ch{F}}=2$.

When an autonomous concentration control maintains the concentrations of the
species F, S, W, and P constant (like in Fig. 1), the CRN relaxes towards a
nonequilibrium steady state. Indeed, the CRN continuously exchanges the $Y$
species with the surroundings (see the fluxes ${\boldsymbol{I}}^{Y}$ in Fig.
4a). At steady state, the species F (resp. S) is constantly injected into the
CRN, i.e., $\overline{{I}}_{\ch{F}}>0$ (resp. $\overline{{I}}_{\ch{S}}>0$),
while the species W (resp. P) is constantly extracted at the same rate, i.e.,
$\overline{{I}}_{\ch{W}}=-\overline{{I}}_{\ch{F}}<0$ (resp.
$\overline{{I}}_{\ch{P}}=-\overline{{I}}_{\ch{S}}<0$). This physically means
that the CRN operates as a steady state engine interconverting $\ch{F}$ into
$\ch{W}$ and $\ch{S}$ into $\ch{P}$. From a thermodynamic standpoint, the
dynamics continuously dissipates (free) energy (see Fig. 4b). In the
transient, the dissipated energy $\dot{\Sigma}$ is provided by the decrease of
the semigrand free energy ${\mathrm{d}\mathcal{G}}/{\mathrm{d}t}<0$, and by
the nonconservative work $\dot{w}_{\text{nc}}>0$. At steady state, the
dissipation due to the continuous interconversion of $\ch{F}$ into $\ch{W}$
and $\ch{S}$ into $\ch{P}$ is solely balanced by the nonconservative work,
i.e., $\overline{\dot{\Sigma}}=\overline{\dot{w}}_{\text{nc}}$, since the
semigrand free energy is constant, i.e.,
${\mathrm{d}\mathcal{G}}/{\mathrm{d}t}=0$. We now turn to the situation where
an autonomous flux control is applied to the CRN. As in Fig. 2, we consider
that the species S and F are continuously injected in the CRN at constant
rates ${I}_{\ch{S}}$ and ${I}_{\ch{F}}$, while the species P and W are
extracted at rates that are proportional to their concentrations, i.e.,
${I}_{\ch{P}}=-k_{\ch{P}}[\ch{P}]$ and ${I}_{\ch{W}}=-k_{\ch{W}}[\ch{W}]$.
Depending on the injection fluxes ${I}_{\ch{S}}$ and ${I}_{\ch{F}}$ and the
extraction rate constants $k_{\ch{P}}$ and $k_{\ch{W}}$, the dynamics can be
significantly different as shown in Fig. 4. When for example
${I}_{\ch{P}}=-0.05[\ch{P}]$, ${I}_{\ch{W}}=-0.1[\ch{W}]$,
${I}_{\ch{S}}=0.025$, and ${I}_{\ch{F}}=0.05$, the CRN relaxes towards a
steady state like when the concentration control is applied. At steady state
the extraction fluxes of the species ${\ch{P}}$ and ${\ch{W}}$ balance exactly
the injection fluxes of the species ${\ch{S}}$ and ${\ch{F}}$ (Fig. 4c):
$\overline{{I}}_{\ch{W}}=-\overline{{I}}_{\ch{F}}$ and
$\overline{{I}}_{\ch{P}}=-\overline{{I}}_{\ch{S}}$. The CRN operates again as
a steady state engine interconverting $\ch{F}$ into $\ch{W}$ and $\ch{S}$ into
$\ch{P}$. From a thermodynamic standpoint (see Fig. 4d), we identify three
different regimes. For $t<1.7$, the energy dissipated by the dynamics
$\dot{\Sigma}$ is provided by the decrease of the semigrand free energy
${\mathrm{d}\mathcal{G}}/{\mathrm{d}t}<0$, by the nonconservative work
$\dot{w}_{\text{nc}}>0$ and also by the driving work
$\dot{w}_{\text{driv}}>0$. For $1.7<t<26$, the energy provided by the
nonconservative work and the driving work is greater than the dissipation,
$\dot{w}_{\text{nc}}+\dot{w}_{\text{driv}}-\dot{\Sigma}>0$, leading to an
increase of the semigrand free energy
${\mathrm{d}\mathcal{G}}/{\mathrm{d}t}>0$. For $t>26$, steady state is
reached: the driving work as well as the time derivative of the semigrand free
energy vanish and only the nonconservative work balances the dissipation,
i.e., $\overline{\dot{\Sigma}}=\overline{\dot{w}}_{\text{nc}}$, like in the
case of concentration control. Thus, the steady state dynamics and
thermodynamics under flux control (unlike the transient ones) are
qualitatively the same as under concentration control. However, the specific
steady state to which the CRN relaxes depends control mechanism. This can be
easily verified by comparing the values of the exchange fluxes in Fig. 4a and
4c as well as of the thermodynamic quantities in Fig. 4b and 4d. When for
example ${I}_{\ch{P}}=-1[\ch{P}]$, ${I}_{\ch{W}}=-2[\ch{W}]$,
${I}_{\ch{S}}=1$, and ${I}_{\ch{F}}=2$, the CRN never reaches steady state.
The concentrations of S and F continuously grow (Fig. 4e), while the
concentrations of P and W become constant and very small. This happens because
the species S and F are injected too quickly into the CRN to be effectively
interconverted into P and W and extracted: ${I}_{\ch{F}}>-{I}_{\ch{W}}$ and
${I}_{\ch{S}}>-{I}_{\ch{P}}$. Thus, they accumulate. From a thermodynamic
standpoint (see Fig. 4f), the dissipation continuously increases and it is
mostly balanced by the nonconservative work, i.e.,
$\dot{\Sigma}\simeq\dot{w}_{\text{nc}}$. On the other hand, the driving work
continuously extracts energy from the CRN which is provided by the decrease of
the semigrand free energy, i.e.,
$\dot{w}_{\text{driv}}\simeq-{\mathrm{d}\mathcal{G}}/{\mathrm{d}t}$. This
thermodynamic behavior is not specific to the model: it is a direct
consequence of the CRN being close to a nonequilibrium steady state that
changes in time. Indeed, in the long time limit the concentrations of
$\ch{S}$, $\ch{F}$ are large enough that the exchange processes act as a
perturbation, while the concentrations P and W are almost constant. This
implies that the entropy production rate is (mostly) balanced by the non
conservative work (as if the CRN was at steady state), while the energetic
cost of changing the steady state accounted by
${\mathrm{d}\mathcal{G}}/{\mathrm{d}t}$ is provided by the driving work. In
this regime, the dynamics and thermodynamics under flux control are
significantly different from those under concentration control.

##  VI Discussion

We discuss now the general differences between concentration and flux control
focusing in particular on the various contributions to the entropy production
rate (12). In the initial transient of the dynamics, the flux control always
leads to a time dependent variation of the concentrations of the exchanged
species which can be obtained only with an nonautonomous concentration
control. This happens because the exchange fluxes do not balance the variation
to the concentrations due to the chemical reactions. Thus, part of the free
energy exchange between CRNs and the surroundings is always given by a driving
work contribution if the flux control is applied. If CRNs are nondetailed
balanced, another contribution to the free energy exchange comes from the
nonconservative work independently of whether concentration control or flux
control is imposed. In the long time limit, CRNs show (qualitatively) the same
phenomenology with both control mechanisms if a steady state is reached: the
dissipation is balanced by the nonconservative work, i.e.,
$\overline{\dot{\Sigma}}={\overline{\dot{w}}_{\text{nc}}}$. On the other hand,
the relaxation towards a steady state is not always granted under flux
control. When this is the case, we found that the variation to the semigrand
free energy of CRNs is (almost) balanced by the driving work, i.e.,
${\mathrm{d}\mathcal{G}}/{\mathrm{d}t}\simeq\dot{w}_{\text{driv}}$, while the
dissipated free energy (in case of nondetailed balanced CRNs) is (mostly)
provided by the nonconservative work, i.e.,
$\dot{\Sigma}\simeq\dot{w}_{\text{nc}}$. This is a consequence of the fact
that CRNs evolve close to a steady state which changes in time. If CRNs are
detailed balanced, the entropy production rate becomes very small compared to
the absolute value of the variation of the semigrand free energy, i.e.,
$\dot{\Sigma}\ll|{\mathrm{d}\mathcal{G}}/{\mathrm{d}t}|$, because CRNs evolve
close to an equilibrium steady state. Our results have been obtained for ideal
dilute solutions, but they can be easily generalized to non-ideal solutions by
recognizing that the chemostats directly control the chemical potentials of
the exchanged species via the chemical reaction (III) instead of their
concentrations Avanzini  _et al._ (2021). Our work opens the way to the
investigation of the energetic cost of growth in autocatalytic processes.
Using kinetic arguments, mathematicians have proven that for a specific type
of CRNs (called single linkage class CRNs) continuous growth is prohibited
under concentration control Anderson (2011a, b). According to preliminary
numerical observations, our thermodynamic approach indicates that the
concentration control cannot provide enough energy to balance both the
dissipation and the increase of the semigrand free energy due to continuous
autocatalytic growth. However, we find that flux control might maintain the
CRN close to a nonequilibrium steady state where the concentrations of the
autocatalytic species continuously increase so powering continuous growth. We
leave the rigorous characterization of these numerical observations to future
studies.

##  VII Acknowledgments

This research was funded by the European Research Council project NanoThermo
(ERC-2015-CoG Agreement No. 681456). We thank the reviewers for their valuable
comments.

## References

  * Sekimoto (2010) K. Sekimoto,  _Stochastic Energetics_ , Lecture Notes in Physics (Springer-Verlag Berlin Heidelberg, 2010).
  * 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).
  * Gaspard (2004) Pierre Gaspard, “Fluctuation theorem for nonequilibrium reactions,” J. Chem. Phys. 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).
  * Ge and Qian (2016a) Hao Ge and Hong Qian, “Nonequilibrium thermodynamic formalism of nonlinear chemical reaction systems with waage–guldberg’s law of mass action,” Chem. Phys. 472, 241 – 248 (2016a).
  * Ge and Qian (2016b) Hao Ge and Hong Qian, “Mesoscopic kinetic basis of macroscopic chemical thermodynamics: A mathematical theory,” Phys. Rev. E 94, 052150 (2016b).
  * Kurtz (1970) Thomas G. Kurtz, “Solutions of ordinary differential equations as limits of pure jump markov processes,” J. Appl. Prob. 7, 49–58 (1970).
  * Kurtz (1971) T. G. Kurtz, “Limit theorems for sequences of jump markov processes approximating ordinary differential processes,” J. Appl. Prob. 8, 344–356 (1971).
  * Kurtz (1972) Thomas G. Kurtz, “The relationship between stochastic and deterministic models for chemical reactions,” J. Chem. Phys. 57, 2976–2978 (1972).
  * 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).
  * 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).
  * 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 (2020) Pierre Gaspard, “Template-directed growth of copolymers,” Chaos 30, 043114 (2020).
  * 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).
  * Yoshimura and Ito (2021a) Kohei Yoshimura and Sosuke Ito, “Information geometric inequalities of chemical thermodynamics,” Phys. Rev. Research 3, 013175 (2021a).
  * Yoshimura and Ito (2021b) Kohei Yoshimura and Sosuke Ito, “Thermodynamic uncertainty relation and thermodynamic speed limit in deterministic chemical reaction networks,” Phys. Rev. Lett. 127, 160601 (2021b).
  * Avanzini  _et al._ (2021) Francesco Avanzini, Emanuele Penocchio, Gianmaria Falasco, and Massimiliano Esposito, “Nonequilibrium thermodynamics of non-ideal chemical reaction networks,” J. Chem. Phys. 154, 094114 (2021).
  * 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).
  * Aris (1989) Rutherford Aris,  _Elementary Chemical Reactor Analysis_ (Dover, 1989).
  * von Stockar (2013) Urs von Stockar,  _Biothermodynamics: The Role of Thermodynamics in Biochemical Engineering_ (EPFL Press, 2013).
  * Craciun and Feinberg (2005) Gheorghe Craciun and Martin Feinberg, “Multiple equilibria in complex chemical reaction networks: I. the injectivity property,” SIAM J. Appl. Math. 65, 1526–1546 (2005).
  * 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).
  * Laurent  _et al._ (2021a) Gabin Laurent, David Lacoste, and Pierre Gaspard, “Emergence of homochirality in large molecular systems,” Proc. Natl. Acad. Sci. USA 118 (2021a), 10.1073/pnas.2012741118.
  * Laurent  _et al._ (2021b) Gabin Laurent, Pierre Gaspard, and David Lacoste, “A robust transition to homochirality in complex chemical reaction networks,”  (2021b), arXiv:2107.14634 [cond-mat.stat-mech] .
  * Hochberg and Ribó (2020) David Hochberg and Josep M. Ribó, “Thermodynamic evolution theorem for chemical reactions,” Phys. Rev. Research 2, 043367 (2020).
  * Qian and Beard (2006) Hong Qian and Daniel A. Beard, “Metabolic futile cycles and their functions: a systems analysis of energy and control,” IEE Proceedings - Systems Biology 153, 192–200 (2006).
  * G. (1993) Svehla G., “Nomenclature of kinetic methods of analysis (iupac recommendations 1993),”  _Pure and Applied Chemistry_ ,  65, 2291 (1993).
  * Falasco  _et al._ (2019) Gianmaria Falasco, Tommaso Cossetto, Emanuele Penocchio, and Massimiliano Esposito, “Negative differential response in chemical reactions,” New J. Phys. 21, 073005 (2019).
  * Rao and Esposito (2018b) Riccardo Rao and Massimiliano Esposito, “Conservation laws shape dissipation,” New J. Phys. 20, 023007 (2018b).
  * Note (1) The choice of the conservation laws (13) is not unique, and different choices identify different moieties. The linear combination in Eq. (16) selects a particular set of conservation laws and, correspondingly, of moieties such that each moiety is exchanged with only one chemostat as long as no ${Y_{f}}$ species are chemostatted: $\mathrm{d}_{t}\boldsymbol{m}={\boldsymbol{I}}_{Y_{p}}$.
  * Note (2) When also the ${Y_{f}}$ species are chemostatted, the concentration of each moiety in Eq. (16) depends on the flux with the corresponding ${Y_{p}}$ chemostat and on the fluxes with the ${Y_{f}}$ chemostats that exchange the same moiety according to $\mathrm{d}_{t}\boldsymbol{m}={\boldsymbol{I}}_{Y_{p}}+\leavevmode@ifvmode{\hbox{$\left(\vbox to1.0pt{}\right.$}}\vvmathbb L^{{b}}_{{Y_{p}}}\leavevmode@ifvmode{\hbox{$\left)\vbox to1.0pt{}\right.$}}^{-1}\vvmathbb L^{{b}}_{{Y_{f}}}{\boldsymbol{I}}_{Y_{f}}$ (with $\vvmathbb L^{{b}}_{{Y_{f}}}$ the submatrix for the ${Y_{f}}$ species of $\vvmathbb L^{{b}}$).
  * Anderson (2011a) David F. Anderson, “Boundedness of trajectories for weakly reversible, single linkage class reaction systems,” J. Math. Chem. 49, 2275 (2011a).
  * Anderson (2011b) David F. Anderson, “A proof of the global attractor conjecture in the single linkage class case,” SIAM J. Appl. Math. 71, 1487–1508 (2011b).


# Deceiving Semantic Shortcuts on Reasoning Chains:
How Far Can Models Go without Hallucination?

Bangzheng Li1, Ben Zhou2, Fei Wang1, Xingyu Fu2, Dan Roth2 Work in progress.
Muhao Chen3
1University of Southern California
2University of Pennsylvania
3University of California, Davis
<EMAIL_ADDRESS>
###### Abstract

Despite the recent advancement in large language models (LLMs) and their high
performances across numerous benchmarks, recent research has unveiled that
LLMs suffer from hallucinations and unfaithful reasoning. This work studies a
specific type of hallucination induced by semantic associations. Specifically,
we investigate to what extent LLMs take shortcuts from certain keyword/entity
biases in the prompt instead of following the correct reasoning path. To
quantify this phenomenon, we propose a novel probing method and benchmark
called EureQA. We start from questions that LLMs will answer correctly with
utmost certainty, and mask the important entity with evidence sentence
recursively, asking models to find masked entities according to a chain of
evidence before answering the question. During the construction of the
evidence, we purposefully replace semantic clues (entities) that may lead to
the correct answer with distractor clues (evidence) that will not directly
lead to the correct answer but require a chain-like reasoning process. We
evaluate if models can follow the correct reasoning chain instead of short-
cutting through distractor clues. We find that existing LLMs lack the
necessary capabilities to follow correct reasoning paths and resist the
attempt of greedy shortcuts. We show that the distractor semantic associations
often lead to model hallucination, which is strong evidence that questions the
validity of current LLM reasoning.

##  1 Introduction

Figure 1:  Surface-level semantic cues can be eliminated by recursively
replacing the entities with their descriptions.

Significant progress of large language models (LLMs) has showcased their
emergent abilities in a wide range of reasoning tasks, encompassing competence
in program synthesis Kuznia et al. (2022), mathematical reasoning Mishra et
al. (2022), symbolic logic reasoning Gaur and Saunshi (2023), and common-sense
reasoning Feng et al. (2023). Innovations in inference methods have further
improved language models’ reasoning capabilities, either by generating
intermediate steps toward the ultimate solution Wei et al. (2022) or by
decomposing complex inquiries into manageable sub-problems Zhou et al.
(2022a); Yao et al. (2023). Research has further demonstrated that when
provided with these techniques, LLMs have achieved exceptional results on
multi-hop QA benchmarks Geva et al. (2021).

Figure 2:  EureQA

Nonetheless, it remains under-investigated whether the effectiveness of LLMs’
reasoning pattern OpenAI (2023) is genuinely based on a reasoning path or if
these models predominantly depend on semantic associations to generate a
plausible answer. Figure 1 shows how ChatGPT111In this paper, we consistently
utilize gpt-3.5-turbo-0301 snapshot for ChatGPT. fails, as the reasoning depth
grows. Initially, the original question, requesting the identification of a
single individual subject to two distinct conditions, can be easily resolved
by aggregating the producers associated with the film “Gardener of Eden” and
the cast members affiliated with “J. Edgar,” thereafter discerning the
intersecting name in both sets. However, the correct answer given by ChatGPT
could alternatively be influenced by the high occurrence of “Leonardo
DiCaprio” in close association with “Gardener of Eden” and “J. Edgar” in the
model’s training data. In such a scenario, ChatGPT may generate the correct
response simply through this semantic association.

To investigate whether LLMs can follow a correct reasoning path and resist the
attempt to take “greedy shortcuts” during inference, we introduce EureQA
(Extending Underlying reasoning Chains in QA), a sophisticated multi-hop
question answering dataset, meticulously designed for diminishing semantic
associations and gauging the capability of LLMs in undertaking extensively
chained reasoning processes. In EureQA, we design a method that removes
semantic shortcuts that will lead to the correct answer while adding
“false/noisy semantic cues” that are distracting and irrelevant to the correct
answer. We also guarantee a simple and consistent reasoning path and only use
common facts that should be memorized by the models. In this way, if an LLM
fail on such tasks, it directly suggests that the LLM is taking incorrect
shortcuts during inference. For instance, we modify the original question in
Figure 1 by substituting “J. Edgar” with a placeholder dubbed “Film A” and
supplementing a descriptive sentence about the film. Our intent with this
strategy is to lessen the model’s reliance on direct information from the film
name. We can automatically extend the depth of a reasoning chain by further
replacing entities in the descriptive sentence about “Film A” with
corresponding abstract names. In particular, this generation methodology is
developed on a knowledge base, ensuring that the information employed predates
the training data of the LLMs under examination.

We then subject LLMs to a comprehensive evaluation, employing a variety of
prominent prompting techniques, to provide insights into their effectiveness
in answering questions requiring extended procedures of reasoning. Our
findings suggest that an assessed LLM 222In this context, ChatGPT., as
illustrated in Figure 1 with a paraphrased response, is incapable of
conducting the reasoning procedure for a question with a depth of three.
Although the model have memorized ample information necessary to answer the
question in its parametric knowledge, the model cannot generate the final
response through a logical reasoning trajectory.

To summarize, the contribution of our work are three-fold. First, we propose a
novel method for generating question answering data with extended reasoning
chains, allowing manipulation of the depth of reasoning. Second, we propose
EureQA, a QA dataset specifically designed for evaluating LLMs with reduced
semantic associations. With EureQA, we can examine the reasoning trajectory of
LLMs towards the final answer. Finally, our extensive experimental findings
suggest that contemporary LLMs predominantly depend on semantic associations
in question answering. As entities in questions replaced with alternate
inquires, LLMs are unable to adhere to an accurate reasoning process necessary
for problem resolution.

##  2 EureQA

In this section, we introduce the dataset EureQA for evaluating LLMs on
extended reasoning chains. We start with the definition and generation
approach of a reasoning chain (Section 2.1), followed by the processes for
question generation (Section 2.2) and refinement (Section 2.3).

###  2.1 Reasoning Chain

In EureQA, every question is constructed through an implicit reasoning chain,
as depicted in Figure 2, part a. This chain is structured in layers, each of
which comprises three components: an entity $e_{i}$, an associated fact
$f_{i}$ about $e_{i}$, and a relation $r_{i}$ which connects $e_{i}$ and the
entity in the succeeding layer of the chain, referred to as $e_{i+1}$. The
identity of $e_{i}$ is constrained by $f_{i}$, $e_{i+1}$ and $r_{i}$ together.
The existence and uniqueness of $e_{i}$ are subject to selecting an
appropriate $f_{i}$. By stacking up these layers, one can form an extended
reasoning chain. By identifying the final entity of the last layer, a viable
reasoning chain can be rendered.

To automatically collect such data, we use a structured knowledge base DBpedia
Auer et al. (2007), because DBpedia provides accurate facts and relations
between entities, and all these facts and relations should have been memorized
by most LLMs since Wikipedia is seen during pre-training. We establish on the
_dbpedia-snapshot-2021-09_ data dump to guarantee that the information of
entities is accessible to LLMs at their training phase. The criteria for
choosing valid $e_{i}$, $r_{i}$ and $f_{i}$ to form chains are the
following:333Being “valid” here refers to satisfying the following criteria
for data extraction.

Entity. A valid entity $e_{i}$ must be accurately acknowledged within the
knowledge base as an entity instead of any other type. In our case, it should
be prefixed with http://dbpedia.org/resource/.

Relation. A valid relation $r_{i}$ should connect $e_{i}$ with a plurality of
potential candidate entities $e_{i+1}$. This ensures no direct semantic
shortcut between $e_{i}$ and $e_{i+1}$. For example, the relation “award” is
deemed valid for the entity “Tiger Woods” because multiple triplets, such as
(_Tiger Woods_ , award, _FedDx Cup_) and (_Tiger Woods_ , award, _Vardon
Trophy_), are documented within DBpedia. In contrast, the relation “college”
is invalid as it corresponds to the only instance, (_Tiger Woods_ , college,
_Stanford University_), in our knowledge base.

Fact. A valid fact $f_{i}$ consists of an relation $r_{i}^{f}$ as well as an
entity $e_{i}^{f}$. Distinct from the criteria for relation selection, the
candidate $e_{i}^{f}$, which is connected by $r_{i}^{f}$ via $e_{i}$, must be
unique to ensure that the model can effectively narrow down the candidates of
$e_{i}$ based on $f_{i}$. In this case, “college” is a valid $r_{i}^{f}$ for “
_Tiger Woods_ ” as the entity $e_{i}$ because _Stanford University_ is the
only $e_{i}^{f}$ that fits in the triplet.

Following the previously established criteria, we employ a random-walk
algorithm to fabricate chains of reasoning. The procedure initiates with a
designated seed entity $e_{1}$ in the knowledge base. Subsequently, a relation
$r_{1}$, an associated fact $f_{1}$, and a subsequent entity $e_{2}$ are
chosen randomly to construct a valid layer of the reasoning chain. This
operation is continuously executed until no more valid layer is found or the
cumulative number of layers equals $N_{max}$. We set $N_{max}=5$ in this work.
We perform at most 100 random walks starting from each seed entity to optimize
efficiency and remove duplicated chains. Seed entities are sourced via two
methods, either through prompting ChatGPT or by data extraction from
https://today.yougov.com/.

This approach to reasoning chain generation for EureQA guarantees the
following three properties of each chains:

Reasoning-dependent. Every intermediate layer relies on information from the
subsequent layer in the chain for resolution. Consequently, the model must
commence from the terminal layer and engage in sequential reasoning throughout
the chain.

Length-flexible. The extent of the reasoning chain can be conveniently
adjusted by adding or removing layers, thereby enabling an evaluation of the
depth of reasoning.

Determinism-adjustable. The determinism of the reasoning chain can be altered
by omitting fact $f_{i}$ in a layer, facilitating an evaluation of the model’s
performance in addressing questions with multiple potential answers.

###  2.2 Question Generation

In EureQA, each question, denoted as $q$, is articulated in natural language
and adheres to the above reasoning chain. A layer-by-layer procedure is
employed to translate this structured chain into a human-readable text. Here
we denote $q=q_{0}+q_{1}+...+q_{n}+m_{n+1}$ where $\\{q_{i}|1<=i<=n\\}$
indicates the sub-question for the $i$-th layer and “$+$” stands for
concatenation. $q_{0}$ and $m_{n+1}$ will be introduced later in this
subsection. Figure 2’s part b provides a tangible example of this translation
process. Every layer comprises an entity $e_{i}$, a relation $r_{i}$, a
succeeding entity $e_{i+1}$, and an associated fact $f_{i}$, which includes
the corresponding $r_{i}^{f}$ and $e_{i}^{f}$. The initial step involves
translating the triplets, ($e_{i}$, $r_{i}$ $e_{i+1}$) and ($e_{i}$,
$r_{i}^{f}$, $e_{i}^{f}$) into two statements, $s_{i}^{relation}$ and
$s_{i}^{fact}$, respectively. For this, we generate these statements using
ChatGPT through in-context learning, coupled with two exemplars crafted by
humans.

The subsequent step involves substituting each entity $e$ in every statement
with a placeholder, as illustrated in bottom right of part b in Figure 2. In
particular, each layer now associates with a pair of statements, namely
$s_{i}^{relation}$ and $s_{i}^{fact}$, which derive from the translation of
the relational triplet ($e_{i}$, $r_{i}$ $e_{i+1}$) and the factual triplet
($e_{i}$, $r_{i}^{f}$, $e_{i}^{f}$) correspondingly. In this context, we
proceed to solely obfuscate $e_{i}$ and $e_{i+1}$ in $s_{i}^{relation}$ and
$s_{i}^{fact}$ by substituting them with their respective hypernyms, $h_{i}$
and $h_{i+1}$, whilst preserving the identity of $e_{i}^{f}$. Hypernym details
are extracted from DBpedia if available and sought from LLM otherwise. ChatGPT
is employed in our approach with two exemplars obtained from DBpedia, ensuring
a consistent level of granularity of hypernym. Each hypernym is subsequently
appended with a specific label, which is assigned in an alphabetical sequence,
starting from the letter “A”. For example, the third layer of a reasoning
chain could hold “Actor C stars in TV Series D” as $s_{3}^{relation}$ while
$s_{3}^{fact}$ could be “Actor C is the artist for Along on Christmas Day.”
Here, $e_{3}$ is substituted by “Actor C” and $e_{4}$, which also appears in
the fourth layer, is replaced by “TV Series D.” The resulting statements with
masked entities are denoted as $m_{i}^{relation}$ and $m_{i}^{fact}$
respectively.

The final step of question generation involves generating an interrogative
question, represented as $q_{0}$, and an entity information statement, denoted
as $m_{n+1}$, pertaining to $q_{n}$. As evidenced in bottom left of part b in
Figure 2, the selection of the initial interrogative question is contingent
upon the hypernym of the entity $e_{1}$ in the first layer. To elaborate,
$q_{0}$ starts with “Who is …" if the hypernym $h_{1}$ corresponds to a person
or “What is …" in all other cases. In the example shown in Figure 2, as “Wes
Anderson” has the hypernym “director”. The generation process for $m_{n+1}$ is
essentially uncomplicated. Given that $m_{i}^{relation}$ in the terminal layer
exhibits the relationship between $e_{n}$ and $e_{n+1}$ and no additional
layers exist to identify $e_{n+1}$, we generate $m_{n+1}$ utilizing the
template “$h_{n+1}$ is $e_{n+1}$.”.

###  2.3 Question Refinement

To ascertain the viability of a reasoning chain and its derivative query, we
use a systematic and algorithmic approach to conduct a viability filtering on
generated chains. As mentioned above, each layer of a specific chain specifies
$e_{i}$ by both a relational triplet ($e_{i}$, $r_{i}$ $e_{i+1}$) and a
factual triplet ($e_{i}$, $r_{i}^{f}$, $e_{i}^{f}$). By excluding $e_{i}$ from
these triplets, we query our knowledge base with both ($?variable$, $r_{i}$
$e_{i+1}$) and ($?variable$, $r_{i}^{f}$, $e_{i}^{f}$). The layer is deemed
viable if $e_{i}$ is the solitary solution to the $?variable$ according to the
knowledge base. Accordingly, we only retain those reasoning chains where each
layer has been judged viable as per the filtration outcomes.

In our research, our central objective is to analyze reasoning capabilities
without relying on memorizing factual information. As such, we have also
conducted a knowledge filtering procedure to ensure that questions are devoid
of any knowledge requirement. To illustrate this, we consider a question,
$q=q_{0}+q_{1}+...+q_{n}+m_{n+1}$, and verify both $s_{i}^{relation}$ and
$s_{i}^{relation}$, where $1<=i<=n$. In more explicable terms, we evaluate
knowledge by presenting ChatGPT with the prompt “ _Is this statement correct:
[ $s_{i}$] Yes or No?_” A statement is considered as memorized if the LLM
response includes a “Yes.”444we applied self-consistency strategy Wang et al.
(2022) with three runs for a majority vote. A $q_{i}$ is deemed to have
satisfied the examination criteria if and only if both of $s_{i}^{relation}$
and $s_{i}^{relation}$ are categorized as memorized.

With these two filtering procedures applied, we have collected 1,600 unique
questions for EureQA.

| depth of 5 | depth of 4 | depth of 3 | depth of 2 | depth of 1
---|---|---|---|---|---
| direct | icl | direct | icl | direct | icl | direct | icl | direct | icl
Llama-2-7b-Chat | 25.0 | 22.2 | 36.9 | 30.3 | 43.3 | 43.3 | 44.7 | 44.0 | 57.6 | 50.7
Llama-2-13b-Chat | 25.5 | 21.3 | 41.5 | 31.2 | 43.8 | 37.1 | 47.5 | 46.5 | 59.8 | 58.0
Llama-2-70b-Chat | 37.5 | 57.2 | 50.3 | 60.7 | 46.4 | 62.4 | 50.9 | 59.8 | 60.3 | 72.6
ChatGPT | 47.5 | 57.3 | 51.8 | 62.8 | 53.2 | 62.1 | 60.6 | 64.7 | 61.4 | 73.5
GPT-4 | 68.1 | 77.8 | 67.6 | 81.4 | 64.4 | 80.9 | 62.7 | 75.5 | 79.5 | 87.7
Table 1: Accuracy of various LLMs across different depths of reasoning. We
evaluate two prompt strategies: direct zero-shot prompt and icl with two
examplers.

##  3 Experiment Setup

This section presents the experiment configurations employed for assessing the
long-chain reasoning capabilities of LLMs through EureQA. The discourse begins
with a discussion on the model configurations (Section 3.1), subsequently
touching upon the adopted prompting methodologies (Section 3.2), and
ultimately, it scrutinizes the evaluation configurations (Section 3.3).

###  3.1 Model Configuration

We consider state-of-the-art LLMs including both open- and closed-source ones,
ranked in the number of parameters: Llama-2-7b-Chat, Llama-2-13b-Chat,
Llama-2-70b-Chat Touvron et al. (2023),555Based on our empirical findings,
Chat versions of Llama-2 exhibit more natural and reasonable responses to our
reasoning-based questions than other versions. ChatGPT (gpt-3.5-turbo-0301)
and GPT-4 (gpt-4-0314). All Llama series of models run with a temperature
$\tau$ of 0.6 while all GPT series of models run with a temperature $\tau$ of
0.7.

###  3.2 Prompting Methods

We employ two prompting methods to assess LLMs on EureQA: namely, direct and
icl. The direct technique presents raw questions to the model without any
supplementary context. On the other hand, icl method utilizes a few-shot
approach, preceding each query with two context-specific exemplars. Each
exemplar is characterized by a sample problem, succeeded by the statement
"Let’s solve this question step by step" and a step-by-step solution written
by humans. In these solutions, the problem is analyzed right from the final
layer through the initial layer to give the correct answer, thereby obviating
any necessity for backtracking during the problem-solving procedure.

###  3.3 Evaluation Protocol

We evaluate the accuracy of LLMs towards the correct answer. Considering that
EureQA functions as a free-form QA benchmark and our observation indicates
that LLMs typically respond with comprehensive reasoning processes, we employ
a specific criterion for string matching: If the correct answer, identifiable
as an entity, is present in the LLM response, we deem such a response as
correct.

We evaluate LLMs on questions with different layers of underlying reasoning
chain. Specifically, we have 1,600 5-layer questions, 200 4-layer questions,
200 3-layer questions, 200 2-layer questions and 200 1-layer questions.
Experiments are conducted on these sets respectively.

##  4 Results

The results of the evaluated LLMs on EureQA, as displayed in Table 1,
illustrate the variations in performance across different depths of reasoning.
In general, with the entities recursively substituted by the descriptions of
reasoning chaining layers, and therefore eliminating surface-level semantic
cues, these models generate more incorrect answers. When the reasoning depth
increases from one to five, there is a notable decline in performance for all
models, with an average accuracy decrease of 21.3% when using the icl prompt.

To identify the errors in these LLM responses, we randomly sampled 20 data
points where GPT-4 with icl, the best performing method in this benchmark, has
generated incorrect answers. The most frequent error pattern lies in
hallucination during reasoning process, encompassing 80% of the cases. In such
cases, an intermediate reasoning stage erroneously conceived an inaccurate,
albeit plausible answer, culminating in accumulated errors in the final
answer. For instance, the “founder of a rocket company” was misconstrued as “
_Elon Musk_ , the founder of SpaceX”, rather than the correct answer of “
_Jeff Bezos_ , the founder of Blue Origin.” This phenomenon proves our claim
that LLM relies on semantic cues for problem-solving instead of cognitive
reasoning.

It can also be concluded from the results that the demonstrated human-written
exemplars do help improve the performance. For Llama-2-70b-Chat, ChatGPT, and
GPT-4, using the icl prompt consistently leads to better performance than the
direct prompt. This suggests that these models can benefit from examples that
provide explicit and well-crafted human-written reasoning processes. These
exemplars encourage the models to resolve the questions with the correct
reason instead of guessing with semantic cues. Nonetheless, even with clear
reasoning processes provided, models still tend to fall for incorrect
shortcuts.

A different behavior was observed in Llama-2-7b-Chat and Llama-2-13b-Chat,
wherein direct prompting gained better results over icl prompting. To
understand this occurrence of these smaller models, we randomly sampled 20
data points, where direct prompting provided accurate answers, in contrast to
the incorrect responses from icl. The primary identified errors can be
recognized as the weakness of these small models in processing long texts,
which constitute approximately one-fifth of the cases. In these cases, models
confused information included in the exemplars with the details in the queried
questions. In this case of icl, models follow the exemplars to elaborate each
step of the deductive process, which results in generating significantly more
tokens before the prediction of the subsequent token. These extra tokens,
although helpful for the human reasoning process, distract these models from
gathering semantic cues and lead to an incorrect answer.

##  5 Related Work

Hallucination. Existing research on analyzing LLMs hallucination mainly focus
on three aspects: input-conflicting hallucination, context-conflicting
hallucination, and fact-conflicting hallucination Zhang et al. (2023). Input-
conflicting hallucination happens when the LLMs outputs are divergent from
prompts Maynez et al. (2020). Self-contradictory outputs can occur in long-
form or multi-turn answers, implying that LLMs lose track of the core input
information during generation Liu et al. (2023). Factual based errors – mis-
aligns with established world knowledge – are most widely studied since they
can happen in various LLMs to affect multiple tasks’ performance Min et al.
(2023); Li et al. (2023). However, none of these existing works look into
hallucination on reasoning-required questions and test whether the models are
following the correct reasoning path versus the semantic token bias.

Reasoning capability of LLMs. The success of in-context learning Brown et al.
(2020) and instruction tuning Wei et al. (2021) have inspired various works to
solve reasoning tasks by prompting LLMs. More advanced prompting strategies
have been proposed to enhance reasoning capabilities of LLMs. These
investigations take advantage of the heuristic nature of human problem-solving
procedures and incorporate them into textual prompts as guidance for the
models. An example of this can be seen in chain-of-thought prompting Wei et
al. (2022), a method in which intermediate steps are generated leading up to a
final solution. Yao et al. (2023) and Zhou et al. (2022b) proposed to
decompose the questions into simpler, manageable sub-problems to as a method
to facilitate complex reasoning. Nonetheless, this work questions the
exhibited “reasoning capabilities” by distinguishing it with the semantic
shortcuts in the questions.

##  6 Conclusion

This paper proposes a novel QA benchmark for probing LLM hallucinations
induced by semantic associations. We introduced a systematic method for
generating question answering data with extended reasoning chains. This data
enable us to examine the reasoning trajectories of LLMs and our current
experimental results indicate that LLMs predominantly depend on semantic
shortcuts for reasoning and such behavior contributes to its failure. This
ongoing research will produce more comprehensive analyses and resources in the
future.

## References

  * Auer et al. (2007) Sören Auer, Christian Bizer, Georgi Kobilarov, Jens Lehmann, Richard Cyganiak, and Zachary Ives. 2007.  Dbpedia: A nucleus for a web of open data.  In _Proceedings of the 6th International The Semantic Web and 2nd Asian Conference on Asian Semantic Web Conference_ , ISWC’07/ASWC’07, page 722–735, Berlin, Heidelberg. Springer-Verlag.
  * Brown et al. (2020) Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al. 2020.  Language models are few-shot learners.  _Advances in neural information processing systems_ , 33:1877–1901.
  * Feng et al. (2023) Yu Feng, Ben Zhou, Haoyu Wang, Helen Jin, and Dan Roth. 2023.  Generic temporal reasoning with differential analysis and explanation.  In _Proceedings of the 61st Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 12013–12029, Toronto, Canada. Association for Computational Linguistics.
  * Gaur and Saunshi (2023) Vedant Gaur and Nikunj Saunshi. 2023.  Reasoning in large language models through symbolic math word problems.  In _Findings of the Association for Computational Linguistics: ACL 2023_ , pages 5889–5903, Toronto, Canada. Association for Computational Linguistics.
  * Geva et al. (2021) Mor Geva, Daniel Khashabi, Elad Segal, Tushar Khot, Dan Roth, and Jonathan Berant. 2021.  Did aristotle use a laptop? a question answering benchmark with implicit reasoning strategies.  _Transactions of the Association for Computational Linguistics_ , 9:346–361.
  * Kuznia et al. (2022) Kirby Kuznia, Swaroop Mishra, Mihir Parmar, and Chitta Baral. 2022.  Less is more: Summary of long instructions is better for program synthesis.  In _Proceedings of the 2022 Conference on Empirical Methods in Natural Language Processing_ , pages 4532–4552, Abu Dhabi, United Arab Emirates. Association for Computational Linguistics.
  * Li et al. (2023) Junyi Li, Xiaoxue Cheng, Wayne Xin Zhao, Jian-Yun Nie, and Ji-Rong Wen. 2023.  Halueval: A large-scale hallucination evaluation benchmark for large language models.  _arXiv e-prints_ , pages arXiv–2305.
  * Liu et al. (2023) Nelson F Liu, Kevin Lin, John Hewitt, Ashwin Paranjape, Michele Bevilacqua, Fabio Petroni, and Percy Liang. 2023.  Lost in the middle: How language models use long contexts.  _arXiv preprint arXiv:2307.03172_.
  * Maynez et al. (2020) Joshua Maynez, Shashi Narayan, Bernd Bohnet, and Ryan McDonald. 2020.  On faithfulness and factuality in abstractive summarization.  _arXiv preprint arXiv:2005.00661_.
  * Min et al. (2023) Sewon Min, Kalpesh Krishna, Xinxi Lyu, Mike Lewis, Wen-tau Yih, Pang Wei Koh, Mohit Iyyer, Luke Zettlemoyer, and Hannaneh Hajishirzi. 2023.  Factscore: Fine-grained atomic evaluation of factual precision in long form text generation.  _arXiv preprint arXiv:2305.14251_.
  * Mishra et al. (2022) Swaroop Mishra, Matthew Finlayson, Pan Lu, Leonard Tang, Sean Welleck, Chitta Baral, Tanmay Rajpurohit, Oyvind Tafjord, Ashish Sabharwal, Peter Clark, and Ashwin Kalyan. 2022.  LILA: A unified benchmark for mathematical reasoning.  In _Proceedings of the 2022 Conference on Empirical Methods in Natural Language Processing_ , pages 5807–5832, Abu Dhabi, United Arab Emirates. Association for Computational Linguistics.
  * OpenAI (2023) OpenAI. 2023.  Gpt-4 technical report.
  * Touvron et al. (2023) Hugo Touvron, Louis Martin, Kevin Stone, Peter Albert, Amjad Almahairi, Yasmine Babaei, Nikolay Bashlykov, Soumya Batra, Prajjwal Bhargava, Shruti Bhosale, et al. 2023.  Llama 2: Open foundation and fine-tuned chat models.  _arXiv preprint arXiv:2307.09288_.
  * Wang et al. (2022) Xuezhi Wang, Jason Wei, Dale Schuurmans, Quoc Le, Ed Chi, Sharan Narang, Aakanksha Chowdhery, and Denny Zhou. 2022.  Self-consistency improves chain of thought reasoning in language models.  _arXiv preprint arXiv:2203.11171_.
  * Wei et al. (2021) Jason Wei, Maarten Bosma, Vincent Y Zhao, Kelvin Guu, Adams Wei Yu, Brian Lester, Nan Du, Andrew M Dai, and Quoc V Le. 2021.  Finetuned language models are zero-shot learners.  _arXiv preprint arXiv:2109.01652_.
  * Wei et al. (2022) Jason Wei, Xuezhi Wang, Dale Schuurmans, Maarten Bosma, Fei Xia, Ed Chi, Quoc V Le, Denny Zhou, et al. 2022.  Chain-of-thought prompting elicits reasoning in large language models.  _Advances in Neural Information Processing Systems_ , 35:24824–24837.
  * Yao et al. (2023) Shunyu Yao, Dian Yu, Jeffrey Zhao, Izhak Shafran, Thomas L. Griffiths, Yuan Cao, and Karthik Narasimhan. 2023.  Tree of Thoughts: Deliberate problem solving with large language models.
  * Zhang et al. (2023) Yue Zhang, Yafu Li, Leyang Cui, Deng Cai, Lemao Liu, Tingchen Fu, Xinting Huang, Enbo Zhao, Yu Zhang, Yulong Chen, et al. 2023.  Siren’s song in the ai ocean: A survey on hallucination in large language models.  _arXiv preprint arXiv:2309.01219_.
  * Zhou et al. (2022a) Ben Zhou, Kyle Richardson, Xiaodong Yu, and Dan Roth. 2022a.  Learning to decompose: Hypothetical question decomposition based on comparable texts.  _EMNLP_.
  * Zhou et al. (2022b) Denny Zhou, Nathanael Schärli, Le Hou, Jason Wei, Nathan Scales, Xuezhi Wang, Dale Schuurmans, Claire Cui, Olivier Bousquet, Quoc Le, et al. 2022b.  Least-to-most prompting enables complex reasoning in large language models.  _arXiv preprint arXiv:2205.10625_.