# Characterizing the conditions for indefinite growth in open chemical
reaction networks

Shesha Gopal Marehalli Srinivas<EMAIL_ADDRESS>Complex Systems and
Statistical Mechanics, Department of Physics and Materials Science, University
of Luxembourg, L-1511 Luxembourg City, Luxembourg  Francesco Avanzini
<EMAIL_ADDRESS>Complex Systems and Statistical Mechanics,
Department of Physics and Materials Science, University of Luxembourg, L-1511
Luxembourg City, Luxembourg  Department of Chemical Sciences, University of
Padova, Via F. Marzolo, 1, I-35131 Padova, Italy  Massimiliano Esposito
<EMAIL_ADDRESS>Complex Systems and Statistical Mechanics,
Department of Physics and Materials Science, University of Luxembourg, L-1511
Luxembourg City, Luxembourg

###### Abstract

The thermodynamic and dynamical conditions necessary to observe indefinite
growth in homogeneous open chemical reaction networks (CRNs) satisfying mass
action kinetics were presented in Srinivas et al. (2023): Unimolecular CRNs
can only accumulate equilibrium concentrations of species while multimolecular
CRNs are needed to produce indefinite growth with nonequilibrium
concentrations. Within multimolecular CRNs, pseudo-unimolecular CRNs produce
nonequilibrium concentrations with zero efficiencies. Nonequilibrium growth
with finite efficiencies requires dynamically nonlinear CRNs. In this paper,
we provide a detailed analysis supporting these results. Mathematical proofs
are provided for growth in unimolecular and pseudo-unimolecular CRNs. For
multimolecular CRNs, four models displaying very distinctive topological
properties are extensively studied, both numerically and partly analytically.

##  I Introduction

Open chemical reaction networks (CRNs) are known to exhibit rich dynamical
behaviors such as nonequilibrium steady states [1], chemical oscillations [2],
chaotic dynamics [3], kinetic non-invertibility [4], patterns [5, 6, 7] and
waves [8, 9]. Recent progress in nonequilibrium thermodynamics of open CRNs
[10, 11, 12, 13, 14, 15, 16] nowadays allows to characterize the energetics of
these complex dynamics. In this work, we study the dynamics and energetics of
indefinite chemical growth in open CRNs: an indefinite increase in the
concentrations of species. The literature on this topic has focused on the
dynamics of CRNs made of irreversible reactions which prevent a consistent
thermodynamics description [17, 18, 19, 20, 21, 22]. Recently, it was proven
numerically that under certain chemosttating conditions, fully reversible open
CRNs can undergo growth [16]. Building on this work, in a companion paper
[23], we argued that growth needs an influx of species into the system at a
constant rate and that growth with the accumulation of nonequilibrium
concentrations of species is only possible in open multimolecular CRNs. In
this paper, we not only provide extensive analytical and numerical support for
these results in a self-contained way, but we also extend their scope by
considering additional models.

The plan of the paper is as follows. In Sec.II we define the dynamics and
thermodynamics of open CRNs, as well as different chemostatting mechanisms
(resp. concentration, flux and mixed control) and indefinite growth. In
Sec.III, we consider unimolecular CRNs. We prove that they only grow under
flux control and that their concentrations always remains close to equilibrium
values, implying a vanishing dissipation and an optimal efficiency of growth
tending to one. In Sec.IV we turn to pseudo-unimolecular CRNs, i.e., a
subclass of multimolecular CRNs displaying linear dynamics. We find that
similarly to unimolecular CRNs, they can only grow under flux control, but
unlike them, concentrations can be far from equilibrium. As a result, the
dissipation scales extensively in time and the efficiency of growth goes to
zero. In Sec. V, we consider four multimolecular CRNs displaying nonlinear
dynamics. Each of them have different topological properties and are studied
(numerically and partly analytically) under the three chemosttating
procedures. We find that growing multimolecular CRNs can show striking
differences in the growth dynamics compared to unimolecular and pseudo-
unimolecular CRNs. Their growing concentrations scale nonlinearly with time
and growth regimes depend on their initial concentrations. Furthermore, we
find that multimolecular CRNs can grow with nonequilibrium concentrations
under flux control and mixed control with an efficiency strictly between zero
and one. Conclusions are drawn in Sec.VI.

##  II Chemical Reaction Networks

###  II.1 Dynamics

We consider CRNs in ideal dilute solutions. The chemical species, labeled
$\alpha\in Z$, are interconverted via elementary [24], reversible, mass-
balanced chemical reactions $\rho$ of the form,

$\bm{\nu}_{+\rho}\cdot\bm{\alpha}\xrightleftharpoons[-\rho]{+\rho}\bm{\nu}_{-\rho}\cdot\bm{\alpha}\,.$
(1)

Here, $\bm{\alpha}=(\dots,\alpha,\dots)^{\intercal}$ and
$\bm{\nu}_{+\rho}=(\dots,\nu_{\alpha,+\rho},\dots)^{\intercal}$ (resp.
$\bm{\nu}_{-\rho}=(\dots,\nu_{\alpha,-\rho},\dots)^{\intercal}$) is the vector
collecting the stoichiometric coefficients of the forward $+\rho$ (resp.
backward $-\rho$) reaction. The topology of CRNs is encoded in the
stoichiometric matrix $\mathbb{S}$ whose columns are given by
$\mathbb{S}_{\rho}=\bm{\nu}_{-\rho}-\bm{\nu}_{+\rho}$.

We consider open CRNs: some species are chemostatted, namely, they are
exchanged with the surroundings. Hence, the concentrations of the chemical
species $\bm{z}=(\dots,[\alpha],\dots)^{\intercal}$ evolve according to the
following rate equation:

$d_{t}\bm{z}=\mathbb{S}\bm{j}(\bm{z})+\bm{I}(\bm{z})\,.$ (2)

The vector $\bm{I}(\bm{z})=(\dots,I_{\alpha},\dots)$ accounts for the
chemostatting procedure by collecting the net fluxes (influx minus outflux) of
each species $\alpha$ between the CRN and the surroundings;
$\mathbb{S}\bm{j}(\bm{z})$ accounts for the change of concentrations due to
the chemical reactions. The vector
$\bm{j}(\bm{z})=(\dots,{j}_{\rho}(\bm{z}),\dots)^{\intercal}$ collects the
reaction currents which are given by the difference between the forward and
the backward fluxes, ${j}_{\rho}(\bm{z})=j_{+\rho}(\bm{z})-j_{-\rho}(\bm{z})$,
satisfying mass action kinetics, i.e.,

$j_{\pm\rho}(\bm{z})={k_{\pm\rho}\prod_{\alpha\in
Z}[\alpha]^{\nu_{\alpha,\pm\rho}}}\equiv
k_{\pm\rho}\bm{z}^{\bm{\nu}_{\pm\rho}}\,.$ (3)

where we introduce the notation $\bm{a}^{\bm{b}}=\Pi_{i}a_{i}^{b_{i}}$. We can
partition the set of species $Z$ into the disjoint subsets of internal $X$ and
chemostatted species $Y$ whose concentrations are given by $\bm{x}$ and
$\bm{y}$, respectively. By definition, the concentrations of the internal
species change only due to the reactions while the concentrations of
chemostatted species change due to the reactions and the chemostatting
procedure. By applying the same splitting to the stoichiometric matrix

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

the rate equation (2) can be rewritten as

$\displaystyle d_{t}\bm{x}$ $\displaystyle=\mathbb{S}^{X}\bm{j}(\bm{z})\,,$
(5) $\displaystyle d_{t}\bm{y}$
$\displaystyle=\mathbb{S}^{Y}\bm{j}(\bm{z})+\bm{I}^{Y}(\bm{z})\,,$ (6)

where the vector $\bm{I}^{Y}(\bm{z})$ is the restriction of the vector
$\bm{I}(\bm{z})$ to the set of chemostatted species $Y$. We consider three
chemostatting procedures: flux control, mixed control, and concentration
control. Under flux control, the influx/outflux of species is constant,

$\bm{I}^{Y}(\bm{z})=\tilde{\bm{I}}\,.$ (7)

Under mixed control, the influx is constant, while the outflux is proportional
to the concentration, i.e.,

$\bm{I}^{Y}(\bm{z})=-\tilde{\mathbb{D}}\bm{y}+\tilde{\bm{I}}\,,$ (8)

where the matrix $\tilde{\mathbb{D}}$, called the extraction matrix, is a
diagonal matrix with entries being the extraction rates $\\{k^{e}_{\alpha}\\}$
with $\alpha\in Y$ and $\tilde{\bm{I}}\geq 0$. Under concentration control,
the concentration of the $Y$ species is held constant, implying,

$\bm{I}^{Y}(\bm{z})=-\mathbb{S}^{Y}\bm{j}(\bm{z}).$ (9)

##### Steady state.

Note that the rate equation (2) can admit two kinds of steady states:
Equilibrium steady states $\bm{z}_{\text{eq}}$ that satisfy
$\bm{j}\left(\bm{z}_{\text{eq}}\right)=0$ and nonequilibrium steady states
$\bm{z}_{\text{ss}}$ that satisfy $\bm{j}\left(\bm{z}_{\text{ss}}\right)\neq
0$, but $\mathbb{S}^{X}\bm{j}(\bm{z}_{\text{ss}})=0$ and
$\mathbb{S}^{Y}\bm{j}(\bm{z}_{\text{ss}})+\bm{I}^{Y}(\bm{z}_{\text{ss}})=0$.
Closed CRNs always relax towards an equilibrium steady state. Open CRNs, if
they relax towards a steady state, generally reach a nonequilibrium steady
state.

##### Coarse-grained dynamics.

When a time scale separation arises between the dynamics of the concentrations
of the internal species $\bm{x}$ and the dynamics of chemostatted species
$\bm{y}$, the former quickly relax to a steady state
$\bm{x}_{\text{ss}}(\bm{y})$ (if it exists) determined by $\bm{y}$. The
reaction currents are then given by

$\bar{\bm{j}}(\bm{y})\equiv{\bm{j}}(\bm{x}_{\text{ss}}(\bm{y}),\bm{y})\,,$
(10)

namely, the steady state current of Eq. (5):
$\mathbb{S}^{X}\bar{\bm{j}}(\bm{y})=0$. Equation (6) can then be coarse-
grained [25, 26] into a closed dynamical equation for the chemostatted
species:

$d_{t}\bm{y}=\mathbb{S}^{Y}\bar{\bm{j}}(\bm{y})+\bm{I}^{Y}(\bm{y})\,.$ (11)

In Sec. V, we will make use of this time scale separation in growing
multimolecular CRNs and compare the solutions of the rate equation (6), i.e.,
$\bm{y}(t)$, with that of the coarse-grained rate equation (11), i.e.,
$\bm{y}_{\text{cg}}(t)$.

###  II.2 Conservation Laws

The linearly-independent left null eigenvectors of the stoichiometric matrix,
$\bm{\ell}^{\lambda}\cdot\mathbb{S}=0$, are called conservation laws: they
identify parts of (or entire) molecules, called moieties, that are preserved
by the reactions. Indeed, their concentrations, defined as
$L^{\lambda}=\bm{\ell}^{\lambda}\cdot\bm{z}$, would be conserved if CRNs were
closed, i.e., $d_{t}L^{\lambda}=\bm{\ell}^{\lambda}\cdot d_{t}\bm{z}=0$ (using
Eq. (2) with $\bm{I}(\bm{z})=0$). Since we consider CRNs with all mass-
balanced chemical reactions, there is always a conservation law, called the
mass conservation law and denoted $\bm{\ell}^{m}$, which involves all the
species: $\bm{\ell}^{m}=(\dots,\ell^{m}_{\alpha},\dots)$ with
$\ell^{m}_{\alpha}\geq 1$. The corresponding concentration,

$L^{m}=\bm{\ell}^{m}\cdot\bm{z}\,,$ (12)

is the mass density. In open CRNs, some moieties are exchanged with the
surroundings and their corresponding concentrations are not conserved anymore:

$d_{t}L^{\lambda}=\sum_{\alpha\in Y}\ell^{\lambda}_{\alpha}~{}{I}_{\alpha}\neq
0\,,$ (13)

where we used Eq. (2) and the definition of conservation laws. This happens
when $\ell^{\lambda}_{\alpha}\neq 0$ for at least one $\alpha\in Y$. The
corresponding conservation laws are said to be broken and labeled
$\\{\bm{\ell}^{\lambda_{b}}\\}$ hereafter. In open CRNs, the mass conservation
law is always broken. On the other hand, the concentrations of some moieties
might still be conserved, i.e., $d_{t}L^{\lambda}=0$, which happens when
$\ell^{\lambda}_{\alpha}=0$ for all $\alpha\in Y$. The corresponding
conservation laws are said to be unbroken and labeled
$\\{\bm{\ell}^{\lambda_{u}}\\}$ hereafter.

Note that the representation (or set) of the conservation laws
$\\{\bm{\ell}^{\lambda}\\}$ is not unique. Different representations identify
different moieties.

###  II.3 Thermodynamics

The theory of nonequilibrium thermodynamics of CRNs is based on two
assumptions. First, all the degrees of freedom, except for the concentrations,
are at equilibrium [13, 15] (the temperature $T$ is set by the solvent and the
solution is assumed incompressible). Thus, all the thermodynamic quantities
have the same form as in equilibrium thermodynamics but are evaluated at
nonequilibrium concentrations. The (vector of) chemical potentials,
$\bm{\mu}=(\dots,\mu_{\alpha},\dots)^{\intercal}$, reads

$\bm{\mu}=\bm{\mu}^{0}+RT~{}\text{ln}\left(\bm{z}\right)\,,$ (14)

where $\bm{\mu}^{0}=(\dots,\mu^{0}_{\alpha},\dots)^{\intercal}$ is the vector
of standard chemical potentials and
$\text{ln}\left(\bm{z}\right)=(\dots,\text{ln}([\alpha]),\dots)^{\intercal}$.
Second, the local detailed balance assumption establishes a correspondence
between dynamics and thermodynamics according to

${RT}~{}\text{ln}\left(\frac{k_{+\rho}}{k_{-\rho}}\right)=-\bm{\mu}^{0}\cdot\mathbb{S}_{\rho}\,.$
(15)

In this framework, the second law reads [13]

${T}{\dot{\Sigma}}={\dot{w}}_{\text{c}}-d_{t}G~{}\geq 0\,.$ (16)

Here, $G(\bm{z})$ is the Gibbs free energy

$G(\bm{z})=\bm{\mu}\cdot\bm{z}-RT||\bm{z}||\,,$ (17)

where $||\bullet||$ denotes the 1-norm of a vector and hence
$||\bm{z}||=\sum_{\alpha\in Z}[\alpha]$. $\dot{\Sigma}$ is the entropy
production rate (EPR),

$\begin{split}T{\dot{\Sigma}}=-\bm{\mu}\cdot\mathbb{S}\bm{j}(\bm{z})=RT\sum_{\rho}j_{\rho}(\bm{z})~{}\text{ln}\left(\frac{j_{+\rho}(\bm{z})}{j_{-\rho}(\bm{z})}\right)\geq
0\,,\end{split}$ (18)

quantifying the dissipated free energy, and ${\dot{w}}_{\text{c}}$ is the
chemical work rate,

${\dot{w}}_{\text{c}}=\sum_{\alpha\in Y}\mu_{\alpha}{I}_{\alpha},$ (19)

accounting for the free energy supplied from the surroundings via the
chemostats.

In this paper, we exploit the decomposition of the chemical work rate into two
contributions: the moiety work rate $\dot{w}_{\text{m}}$ quantifying the
energetic cost supplied to change the concentrations of the exchanged moieties
in the CRN, and the nonconservative work rate $\dot{w}_{\text{nc}}$ accounting
for the energetic cost to drive currents across the CRN. This splitting is
obtained by first recognizing that chemostatting a species does not always
break a conservation law [12, 13, 15]. Thus, we divide the set of chemostatted
species $Y$ into the set of species $Y_{p}\subseteq Y$ that break conservation
laws, called potential species, and the remaining species $Y_{f}=Y\setminus
Y_{p}$, called force species. With this identification, we can associate a
moiety to a single $Y_{p}$ species so that the moiety work rate has the form

$\dot{w}_{\text{m}}=\bm{\mu}_{Y_{p}}\cdot d_{t}\bm{m}\,,$ (20)

where $\bm{\mu}_{Y_{p}}$ is the vector of chemical potentials of the $Y_{p}$
species and $\bm{m}$ is the concentration vector of the corresponding
exchanged moieties. The latter is given by:

$\bm{m}=\big{(}\mathbb{L}^{{b}}_{Y_{p}}\big{)}^{-1}\mathbb{L}^{{b}}\bm{z}\,,$
(21)

where we introduce the matrix $\mathbb{L}^{{b}}$ whose rows are the broken
conservation laws and the matrix $\mathbb{L}_{Y_{p}}^{{b}}$ that is the
(invertible) submatrix of $\mathbb{L}^{{b}}$ with only the columns
corresponding to the $Y_{p}$ species. Indeed, $\bm{m}$ quantifies the
concentration of the moieties defined by a specific representation of the
conservation laws where $\mathbb{L}_{Y_{p}}^{{b}}$ is the identity matrix
namely, a representation where each potential species carries only one
specific moiety [16].

On the other hand, the force species carry moieties that are already carried
by the potential species. This leads to an energetic cost for chemostatting
different species that carry the same moiety [15, 16] which is captured by the
nonconservative work rate
$\dot{w}_{\text{nc}}=\dot{w}_{\text{c}}-\dot{w}_{\text{m}}$,

$\dot{w}_{\text{nc}}=\left(\bm{\mu}_{Y}-\bm{\mu}_{Y_{p}}\left(\mathbb{L}^{b}_{Y_{p}}\right)^{-1}\mathbb{L}^{b}_{Y}\right)\cdot\bm{I}^{Y}\equiv\mathcal{F}_{Y}\cdot\bm{I}^{Y}\,,$
(22)

where we used Eqs. (19), (20) and (21). We emphasize that, although the
summation in Eq. (22) is over all the $Y$ species, since the submatrix of
$\big{(}\mathbb{L}^{b}_{Y_{p}}\big{)}^{-1}\mathbb{L}^{b}_{Y}$ for the $Y_{p}$
species is the identity matrix, the entries of the vector $\mathcal{F}_{Y}$
(also called nonconservative force) are zero for all the potential species
$Y_{p}$. Thus, $\dot{w}_{\text{nc}}=0$ if only the potential species were
chemostatted. However, the entries of the vector $\mathcal{F}_{Y}$ for the
force species $Y_{f}$ are given by the differences between the chemical
potentials of species carrying the same moiety and thus
$\dot{w}_{\text{nc}}\neq 0$ when the force species are chemostatted.

By combining Eq. (20), (22) with (16), the second law becomes:

$T\dot{\Sigma}={\dot{w}}_{\text{nc}}+\dot{w}_{\text{m}}-d_{t}G~{}\geq 0\,.$
(23)

Note that, in closed CRNs,
$\dot{w}_{\text{c}}=\dot{w}_{\text{m}}=\dot{w}_{\text{nc}}=0$, implying
$d_{t}G=-T\dot{\Sigma}\leq 0$. Thus, the Gibbs free energy monotonously
decreases and can be shown to act as a Lyapunov function [13]. At equilibrium,
$d_{t}G=T\dot{\Sigma}=0$. In open CRNs, at a nonequilibrium steady state, the
terms $d_{t}G$ and $\dot{w}_{\text{m}}$ in Eq (23) vanish as they contain
total time derivatives while the constant EPR is balanced by the
nonconservative work rate, i.e., $T\dot{\Sigma}=\dot{w}_{\text{nc}}>0$.

###  II.4 Growth

We define growth as a limiting state of CRNs with unbounded concentrations:
$\lim_{t\to\infty}\rVert\bm{z}(t)-\bm{z}(0)\rVert=\infty$. When CRNs grow,
also the mass density $L^{m}$ in Eq. (12) and the Gibbs free energy
$G(\bm{z})$ in Eq. (17) are unbounded by definition. The converse is also
true. If the mass density $L^{m}$ and the Gibbs free energy $G(\bm{z})$ are
unbounded, then there is at least one species whose concentration is
unbounded. Indeed, if all the concentrations are bounded from above by some
constant $M$, i.e., $[\alpha](t)\leq M~{}\forall~{}\alpha$, then,
$L^{m}(t)\leq|Z|M$ and ${G(\bm{z})}/{RT}\leq CM~{}\text{ln}(M)$, where $|Z|$
is the number of species and $C$ is a constant.

We are now in the position to make three general statements about growth for
any CRN based on the chemostatting procedure.

First, closed CRNs cannot grow since the mass density is conserved.

Second, we consider CRNs under mixed control such that the extraction rates of
all species are finite while the constant influx $\tilde{\bm{I}}$ is arbitrary
and may act only on a subset of species. Continuous-flow Stirred Tank Reactors
(CSTR) constitute a special case where all extraction rates are the same [27,
28]. Then,

$\begin{split}d_{t}L^{m}&=-\bm{\ell}^{m}\cdot\tilde{\mathbb{D}}\bm{z}+\bm{\ell}^{m}\cdot\tilde{\bm{I}}\,,\\\
&=-\sum_{\alpha\in Z}k^{e}_{\alpha}\ell^{m}_{\alpha}[\alpha]+\sum_{\beta\in
Z}\ell^{m}_{\beta}\tilde{I}_{\beta}\,,\end{split}$ (24)

where we used Eqs. (13) and (8). Let ${k}^{e}$ be the smallest extraction
rate, then

$d_{t}L^{m}\leq-{k}^{e}L^{m}+\bm{\ell}^{m}\cdot\tilde{\bm{I}}\,,$ (25)

which implies

$L^{m}(t)\leq
L^{m}(0)e^{-{k}^{e}t}+\left(\bm{\ell}^{m}\cdot\tilde{\bm{I}}\right)\frac{\left(1-e^{-{k}^{e}t}\right)}{{k}^{e}}\,.$
(26)

Thus, the mass density is bounded in time.

Third, flux control always leads to growth when
$\bm{\ell}^{m}\cdot\tilde{\bm{I}}>0$, namely, when the net influx of species
is larger than the net outflux. Indeed, the evolution equation for the mass
density becomes $d_{t}L^{m}=\bm{\ell}^{m}\cdot\tilde{\bm{I}}$ which leads to
the divergence of the mass density linearly in time. Consequently, the
concentrations can grow at most linearly in time.

##  III Unimolecular CRNs

CRNs are said to be unimolecular when they are exclusively composed of
reactions of the form ${\alpha\xrightleftharpoons{}\beta}$. Each vector
$\bm{\nu}_{+\rho}$ (resp. $\bm{\nu}_{-\rho}$) in Eq. (1) has only one non-zero
entry, that is equal to $+1$, corresponding to the reactant (resp. product) of
reaction $\rho$. Unimolecular CRNs can thus be represented as graphs by
mapping species into nodes and reactions into edges. Without loss of
generality, we assume that the graphs have a single connected component.
Physically, this means that for any two species $\alpha$ and $\beta$, there is
a sequence of reactions connecting them. This assumption is not restrictive as
different components correspond to independent unimolecular CRNs. The
stoichiometric matrix $\mathbb{S}$ becomes the incidence matrix of the graph
and admits only one conservation law, or equivalently one moiety, namely, the
mass conservation law $\bm{\ell}^{m}$ with $\ell^{m}_{\alpha}=1$
$\forall\alpha\in Z$.

###  III.1 Dynamics of Closed Unimolecular CRNs

For unimolecular CRNs, mass action kinetics (Eq. (3)) implies that the
reaction currents are linear functions of the concentrations, i.e.,
$\bm{j}(\bm{z})=\Gamma\bm{z}$, where the entries of the matrix $\Gamma$ are of
the form

$\Gamma_{\rho,\alpha}=\begin{cases}k_{+\rho}~{}~{}\text{if}~{}{\nu_{\alpha,+\rho}=1}\\\
-k_{-\rho}~{}\text{if}~{}{\nu_{\alpha,-\rho}=1}\\\ 0~{}\text{else}\end{cases}$
(27)

Thus, for closed unimolecular CRNs, Eq. (2) becomes

$d_{t}\bm{z}=\mathbb{W}\bm{z}\,,$ (28)

with $\mathbb{W}=\mathbb{S}\Gamma$ being an irreducible rate matrix [29].
Indeed, the off-diagonal elements of $\mathbb{W}$ read

$\mathbb{W}_{\alpha,\beta}=\sum_{\rho}\mathbb{W}^{(\rho)}_{\alpha,\beta}=\sum_{\rho}\begin{cases}k_{+\rho}~{}\text{if}~{}\nu_{\alpha,-\rho}=1{~{}\text{and}~{}\nu_{\beta,+\rho}=1}\\\
k_{-\rho}~{}\text{if}~{}\nu_{\alpha,+\rho}=1{~{}\text{and}~{}\nu_{\beta,-\rho}=1}\\\
0~{}\text{else}\end{cases}$ (29)

where we accounted for (possibly) multiple reactions interconverting the same
species, while the diagonal elements are given by

$\mathbb{W}_{\alpha,\alpha}=-\sum_{\beta\neq\alpha}\mathbb{W}_{\beta,\alpha}\,.$
(30)

From the Perron Frobenius theorem, $\mathbb{W}$ admits one zero eigenvalue
while all the other eigenvalues have negative real parts. Because of the local
detailed balance condition Eq. (15), all $\\{\mathbb{W}^{(\rho)}\\}$ matrices,
and consequently $\mathbb{W}$ too, are detailed-balanced: the eigenvector
corresponding to the zero eigenvalue
$\bm{\pi}^{\text{eq}}=(\dots,{\pi}^{\text{eq}}_{\alpha},\dots)^{\intercal}$
satisfies

$\mathbb{W}_{\alpha,\beta}\pi^{\text{eq}}_{\beta}=\mathbb{W}_{\beta,\alpha}\pi^{\text{eq}}_{\alpha}\,,$
(31)

and must be of the form

$\bm{\pi}^{\text{eq}}=\frac{\exp(-\frac{\bm{\mu}^{0}}{RT})}{\sum_{\alpha\in
Z}\exp{\frac{-\mu^{0}_{\alpha}}{RT}}}\,.$ (32)

The equilibrium steady state of Eq. (28) can be written as

$\bm{z}_{\text{eq}}=L^{m}(0)\bm{\pi}^{\text{eq}}=(\dots,[\alpha]_{\text{eq}}(L^{m}(0)),\dots)^{\intercal}\,,$
(33)

since the mass density is a conserved moiety, where we also introduced the
function $[\alpha]_{\text{eq}}(L^{m}(0))\equiv
L^{m}(0){\pi}^{\text{eq}}_{\alpha}$ to stress that the equilibrium
concentration $[\alpha]_{\text{eq}}$ is a function of the mass density
$L^{m}(0)$.

###  III.2 Dynamics of Open Unimolecular CRNs

For open unimolecular CRNs under any chemostatting procedure, Eq. (2) can be
recast in the form

$d_{t}\bm{a}=\mathbb{V}\bm{a}+\bar{\bm{I}}_{a}\,,$ (34)

where $\bm{a}$ is the vector of the dynamical variables, $\mathbb{V}$ is an
appropriately chosen (real) matrix, and $\bar{\bm{I}}_{a}$ is a constant
vector. For the case of flux (resp. mixed) control, Eq. (34) follows from Eq.
(7) (resp. Eq. (8)) along with mapping the vector of concentrations $\bm{z}$
into $\bm{a}$, the detailed balanced rate matrix $\mathbb{W}$ (resp. the
matrix difference $\mathbb{W}-\mathbb{D}$) into $\mathbb{V}$ and
$\tilde{\bm{I}}$ into $\bar{\bm{I}}_{a}$. To show that Eq. (34) can also
represent concentration control (defined in (9)), we write Eqs. (5) and (6)
for the internal species and chemostatted species as

$\displaystyle d_{t}\bm{x}$
$\displaystyle=\mathbb{W}_{XX}\,\bm{x}+\mathbb{W}_{XY}\,\bm{y}\,,$ (35)
$\displaystyle d_{t}\bm{y}$
$\displaystyle=\mathbb{W}_{YX}\,\bm{x}+\mathbb{W}_{YY}\,\bm{y}+\bm{I}^{Y}({\bm{x},\bm{y}})=\bm{0}\,,$
(36)

where the matrices
$\\{\mathbb{W}_{XX},\mathbb{W}_{XY},\mathbb{W}_{YX},\mathbb{W}_{YY}\\}$ result
from applying the splitting $Z=X\cup Y$ to $\mathbb{W}$. Since the
concentrations $\bm{y}$ are constant, only the concentrations $\bm{x}$ are
dynamical variables. By introducing the new detailed-balanced rate matrix
$\hat{\mathbb{W}}$ with off-diagonal entries

$\hat{\mathbb{W}}_{\alpha,\beta}=\mathbb{W}_{\alpha,\beta}\,,$ (37)

with $\alpha,\beta\in X$ and diagonal entries
$\hat{\mathbb{W}}_{\alpha,\alpha}=-\sum_{\beta\in
X}\hat{\mathbb{W}}_{\beta,\alpha}$, the rate equation (35) becomes

$d_{t}\bm{x}=\hat{\mathbb{W}}\bm{x}-\mathbb{D}\bm{x}+\bar{\bm{I}}\,,$ (38)

with the diagonal matrix $\mathbb{D}$ and constant vector $\bar{\bm{I}}$ given
by,

$\displaystyle\mathbb{D}$
$\displaystyle={\hat{\mathbb{W}}-\mathbb{W}_{XX}}\,,$ (39)
$\displaystyle\bar{\bm{I}}$ $\displaystyle=\mathbb{W}_{XY}\,\bm{y}\,.$ (40)

We further assume without loss of generality that $\hat{\mathbb{W}}$ is
irreducible (if not, different blocks of $\hat{\mathbb{W}}$ would decouple in
Eq. (38) and evolve independently). Finally, by mapping $\bm{x}$ into
$\bm{a}$, $\hat{\mathbb{W}}-\mathbb{D}$ into $\mathbb{V}$ and $\bar{\bm{I}}$
into $\bar{\bm{I}}_{a}$, Eq. (38) takes the form of Eq. (34).

###  III.3  No Growth Under Mixed and Concentration Control

The dynamics of open unimolecular CRNs under mixed control (resp.
concentration control) is given by Eq. (34) with
$\mathbb{V}=\mathbb{W}-\mathbb{D}$ (resp. $\hat{\mathbb{W}}-\mathbb{D}$),
$\bm{a}=\bm{z}$ (resp. $\bm{a}=\bm{x}$) and $\bar{\bm{I}}_{a}=\tilde{\bm{I}}$
(resp. $\bar{\bm{I}}_{a}=\bar{\bm{I}}$). In this case, we prove in Appendix A
that all the eigenvalues of the matrix $\mathbb{V}$ have negative real parts
which also implies that $\mathbb{V}$ is invertible. Then, at long times,
$\bm{a}(t)$ relaxes to the steady state $-\mathbb{V}^{-1}\bar{\bm{I}}_{a}$, as
shown in Appendix B.1. We thus conclude that unimolecular CRNs, under mixed or
concentration control, cannot grow.

###  III.4  Growth Under Flux Control

####  III.4.1 Dynamics

The dynamics of open unimolecular CRNs under flux control is given by Eq. (34)
with $\mathbb{V}=\mathbb{W}$, $\bm{a}=\bm{z}$, and
$\bar{\bm{I}}_{a}=\tilde{\bm{I}}$. We solve Eq. (34) in the basis of
eigenvectors of $\mathbb{W}$ (as detailed balanced rate matrices are
diagonalizable [29]) by introducing the invertible matrix of eigenvectors
$\mathbb{P}$ such that

$\mathbb{P}^{-1}\mathbb{W}\mathbb{P}=\mathbb{\Lambda}\,,$ (41)

where $\mathbb{\Lambda}_{\alpha,\beta}=\delta_{\alpha,\beta}\lambda_{\alpha}$
with $\delta_{\alpha,\beta}$ being the Kronecker delta and
$\\{\lambda_{\alpha}\\}$ being the eigenvalues of $\mathbb{W}$. From the
Perron Frobenius theorem, $\mathbb{W}$ has a unique zero eigenvalue, i.e.,
$\lambda_{1}=0$ and the other eigenvalues satisfy $\lambda_{\alpha}<0$ for
$\alpha\neq 1$. Furthermore, $\mathbb{P}_{\alpha,1}=\pi^{\text{eq}}_{\alpha}$
(given in Eq. (32)), and $\mathbb{P}^{-1}_{1,\alpha}=1$. Then, by defining the
vector $\bm{u}=\mathbb{P}^{-1}\bm{z}$, Eq. (34) becomes

$d_{t}\bm{u}=\mathbb{\Lambda}\bm{u}+\mathbb{P}^{-1}\tilde{\bm{I}}\,,$ (42)

whose solution componentwise reads

$\displaystyle u_{1}(t)$
$\displaystyle=L^{m}(t)=L^{m}(0)+\left(\sum_{\alpha\in
Z}\tilde{I}_{\alpha}\right)t\,,$ (43) $\displaystyle u_{\alpha}(t)$
$\displaystyle=u_{\alpha}(0)e^{\lambda_{\alpha}t}+\left(\frac{(\mathbb{P}^{-1}\tilde{\bm{I}})_{\alpha}}{-\lambda_{\alpha}}\right)\left(1-e^{\lambda_{\alpha}t}\right)\,.$
(44)

Hence, the concentrations $\bm{z}=\mathbb{P}\bm{u}$ in the long time limit
read

$\displaystyle[\alpha](t)$
$\displaystyle=\mathbb{P}_{\alpha,1}L^{m}(t)+\sum_{\beta>1}\mathbb{P}_{\alpha,\beta}\left(\frac{(\mathbb{P}^{-1}\tilde{\bm{I}})_{\beta}}{-\lambda_{\beta}}\right)\,,$
(45a)
$\displaystyle=\pi^{\text{eq}}_{\alpha}L^{m}(t)+c_{\alpha}(\tilde{\bm{I}})\,,$
(45b)

where we collected all the time-independent contributions into
$c_{\alpha}(\tilde{\bm{I}})$, and used Eq. (33) to identify
$\pi^{\text{eq}}_{\alpha}L^{m}(t)$ as the equilibrium concentration
$[\alpha]_{\text{eq}}\left(L^{m}(t)\right)$ of the corresponding closed CRN
with mass density $L^{m}(t)$. Equations (43) and (33) show that the mass
density, as well as the equilibrium concentrations
$[\alpha]_{\text{eq}}\left(L^{m}(t)\right)$, grow linearly in time. This,
together with Eq. (45b), implies that also the instantaneous concentrations
$[\alpha](t)$ grow linearly in time. Indeed, the relative (resp. absolute)
difference between $[\alpha](t)$ and
$[\alpha]_{\text{eq}}\left(L^{m}(t)\right)$ decreases in time (resp. is
constant). Thus, we call this type of growth equilibrium growth.

We now examine the reaction currents of equilibrium growth. Let us denote the
current of the reaction $\beta\xrightleftharpoons{+\rho}\alpha$ as $j_{\rho}$.
Then, by using Eqs. (45b) and (31), in the long time limit, the current reads

$\displaystyle j_{\rho}(t)$
$\displaystyle=~{}\mathbb{W}^{(\rho)}_{\alpha,\beta}[\beta](t)-\mathbb{W}^{(\rho)}_{\beta,\alpha}[\alpha](t)\,,$
$\displaystyle=~{}\mathbb{W}^{(\rho)}_{\alpha,\beta}c_{\beta}(\tilde{\bm{I}})-\mathbb{W}^{(\rho)}_{\beta,\alpha}c_{\alpha}(\tilde{\bm{I}})\equiv
J_{\rho}\,.$ (46)

Notice that while the concentrations grow linearly in time close to their
equilibrium value, the reaction currents become constant and do not vanish
because the difference
$c_{\alpha}(\tilde{\bm{I}})=[\alpha](t)-[\alpha]_{\text{eq}}(L^{m}(t))$ is
constant.

####  III.4.2 Thermodynamics

We now analyze the energetics of equilibrium growth at long times using the
two formulations of the second laws given in Eqs. (16) and (23), respectively.
We start by considering the chemical potentials. By using Eqs. (45b) and (32)
in Eq. (14), the chemical potentials read

$\displaystyle\mu_{\alpha}(t)={}$ $\displaystyle
RT\text{ln}\left(L^{m}(t)\right)-RT~{}\text{ln}\left(\sum_{\beta\in
Z}e^{\frac{-\mu^{0}_{\beta}}{{RT}}}\right)$ (47)
$\displaystyle+RT\text{ln}\left(1+\frac{c_{\alpha}}{\pi^{\text{eq}}_{\alpha}L^{m}(t)}\right)\,,$

which grows logarithmically in time. This allows us to write the vector of
chemical potentials as

$\displaystyle\bm{\mu}=$
$\displaystyle\left[RT~{}\text{ln}\left(L^{m}(t)\right)-RT~{}\text{ln}\left(\sum_{\beta\in
Z}e^{\frac{-\mu^{0}_{\beta}}{{RT}}}\right)\right]\bm{\ell}^{m}$ (48)
$\displaystyle+RT\frac{1}{L^{m}(t)}\left(\frac{\bm{c}}{\bm{\pi}^{\text{eq}}}\right)\,,$

by using $\text{ln}(1+x)\approx x$ for $x\ll 1$ and introducing the vector
$\bm{c}/{\bm{\pi}^{\text{eq}}}=\left(\dots,{c_{\alpha}}/{\pi^{\text{eq}}_{\alpha}},\dots\right)^{\intercal}$.
Thus, by plugging Eq. (48) and Eq. (46) in Eq. (18) and using
$\bm{\ell}^{m}\cdot\mathbb{S}=0$, the EPR becomes

$T\dot{\Sigma}=-RT\left(\frac{\bm{c}}{\bm{\pi}^{\text{eq}}}\cdot\mathbb{S}\bm{J}\right)\frac{1}{L^{m}(t)}=\mathcal{O}(t^{-1})\,,$
(49)

with $\bm{J}\equiv\big{(}\dots,J_{\rho},\dots\big{)}^{\intercal}$. Physically,
Eq. (49) shows a monotonic decay of the EPR in time, which results from the
decreasing relative difference between the instantaneous concentrations
$[\alpha](t)$ and the equilibrium concentrations
$[\alpha]_{\text{eq}}\left(L^{m}(t)\right)$ (see Eq. (45b)). On the other
hand, by using Eq. (48) in Eq. (19) and Eq. (17), the chemical work rate and
the rate of change of the Gibbs free energy can be expressed using
$\bm{\ell}^{m}\cdot d_{t}\bm{z}=\bm{\ell}^{m}\cdot{\tilde{\bm{I}}}$ as

$\displaystyle\begin{aligned}
\dot{w}_{c}=&\left[RT~{}\text{ln}\left(L^{m}(t)\right)-RT~{}\text{ln}\left(\sum_{\beta\in
Z}e^{\frac{-\mu^{0}_{\beta}}{{RT}}}\right)\right](\bm{\ell}^{m}\cdot\tilde{\bm{I}})\\\
&+RT\left(\frac{\bm{c}}{\bm{\pi}^{\text{eq}}}\cdot\tilde{\bm{I}}\right)\frac{1}{L^{m}(t)}=\mathcal{O}(\text{ln}(t))\,,\end{aligned}$
(50) $\displaystyle\begin{aligned}
d_{t}G=&\left[RT~{}\text{ln}\left(L^{m}(t)\right)-RT~{}\text{ln}\left(\sum_{\beta\in
Z}e^{\frac{-\mu^{0}_{\beta}}{{RT}}}\right)\right](\bm{\ell}^{m}\cdot\tilde{\bm{I}})\\\
&+RT\left(\frac{\bm{c}}{\bm{\pi}^{\text{eq}}}\cdot
d_{t}\bm{z}\right)\frac{1}{L^{m}(t)}=\mathcal{O}(\text{ln}(t))\,,\end{aligned}$
(51)

showing a logarithmic increase in the long time limit. Note that the EPR in
Eq. (49), the chemical work in Eq. (50) and the time derivative of the Gibbs
free energy in Eq. (51) satisfy the second law (16). Indeed,

$\begin{split}\dot{w}_{c}-d_{t}G=RT\left(\frac{\bm{c}}{\bm{\pi}^{\text{eq}}}\cdot(\tilde{\bm{I}}-d_{t}\bm{z})\right)\frac{1}{L^{m}(t)}=T\dot{\Sigma}\,.\end{split}$
(52)

This physically means that the chemical work is fully converted into the Gibbs
free energy in the long time limit when the EPR vanishes. Thus, the efficiency
of growth, defined in analogy with Ref. [30] as the ratio between the free
energy stored in the system and the chemical work,

$\eta=\frac{d_{t}G}{\dot{w}_{c}}\,,$ (53)

goes to one in the long time limit and therefore equilibrium growth becomes a
thermodynamically reversible process.

We now turn to the formulation of the second law given in Eq. (23). Since
unimolecular CRNs have only one moiety, namely, the mass density, they admit
only one potential species (hereafter labeled $p$). Therefore, the moiety work
rate (20) and the nonconservative work rate (22) become

$\displaystyle\begin{split}\dot{w}_{\text{m}}&=\bigg{[}RT\text{ln}\left(L^{m}(t)\right)-RT~{}\text{ln}\left(\sum_{\beta\in
Z}e^{\frac{-\mu^{0}_{\beta}}{{RT}}}\right)\bigg{]}\left(\bm{\ell}^{m}\cdot\tilde{\bm{I}}\right)\\\
&~{}~{}~{}+RT\frac{c_{p}}{\pi^{\text{eq}}_{p}L^{m}(t)}(\bm{\ell}^{m}\cdot\tilde{\bm{I}})=\mathcal{O}\left(\text{ln}(t)\right)\,,\end{split}$
(54) $\displaystyle\dot{w}_{\text{nc}}$ $\displaystyle=RT\left[\sum_{\alpha\in
Y}\left(\frac{c_{\alpha}}{\pi^{\text{eq}}_{\alpha}}-\frac{c_{p}}{\pi^{\text{eq}}_{p}}\right)\tilde{I}_{\alpha}\right]\frac{1}{L^{m}(t)}=\mathcal{O}(t^{-1})\,,$
(55)

where we used Eq. (48) together with the fact that $\mathbb{L}^{{b}}_{Y_{p}}$
in Eq. (21) is the scalar $1$.

This allows us to stress two main points. First, the difference between
$\dot{w}_{\text{m}}-d_{t}G$ (obtained using Eqs. (54) and (51)) is of order
$\mathcal{O}(t^{-1})$ implying that all the moiety work done to change the
mass density is stored in the Gibbs free energy in the long time limit. In the
same limit, $\dot{w}_{\text{m}}$ and $d_{t}G$ grow logarithmically in time:

$d_{t}{G}\sim\dot{w}_{\text{m}}=\mathcal{O}\left(\text{ln}\left(t\right)\right)\,.$
(56)

Second, the nonconservative work rate (55), namely, the amount of free energy
spent to keep CRNs out of equilibrium, decreases in time at the same rate as
the EPR (49):

$T\dot{\Sigma}\sim\dot{w}_{\text{nc}}=\mathcal{O}\left({t}^{-1}\right)\,.$
(57)

in agreement with the interpretation of this process as equilibrium growth.

###  III.5 Summary

From a dynamical standpoint, growth in open unimolecular CRNs can occur only
under flux control. In this case, the mass density increases linearly in time;
the instantaneous concentrations approach the equilibrium steady-state
concentrations defined by the growing mass density (see Eq. (45b)); all
concentrations scale linearly in time.

From a thermodynamic standpoint, the work done in increasing the mass density
and the change in free energy scale similarly in time (see Eq. (56)), while
the EPR, as well as the nonconservative work, monotonically decay (see Eq.
(57)). Thus, the efficiency of growth (53) goes to one at long times, implying
that growth is a thermodynamically reversible process.

##  IV Pseudo-unimolecular CRNs

Open CRNs are here said to be pseudo-unimolecular when two conditions are
satisfied. First, the set of species $Z$ can be divided into the subsets
$Z=Z_{l}\cup Y_{h}$ such that every reaction is of the form

$\alpha+\sum_{\gamma\in
Y_{h}}\gamma\,\nu_{\gamma,+\rho}\xrightleftharpoons[-\rho]{+\rho}\beta+\sum_{\gamma\in
Y_{h}}\gamma\,\nu_{\gamma,-\rho}\,,$ (58)

for some $\alpha$ and $\beta$ in $Z_{l}$. Second, the species $Y_{h}$ undergo
concentration control while the $Z_{l}$ species can be internal or
chemostatted.

From a topological standpoint, pseudo-unimolecular CRNs can be represented as
graphs (like unimolecular CRNs), by mapping only the $Z_{l}$ species into
nodes and reactions into edges. The corresponding incidence matrix is given by
the substoiochiometric matrix $\mathbb{S}^{l}$ obtained by applying the
splitting $Z=Z_{l}\cup Y_{h}$ to $\mathbb{S}$, i.e.,
$\mathbb{S}=\left(\mathbb{S}^{l},\mathbb{S}^{h}\right)^{\intercal}$. Notice
that this implies that i) pseudo-unimolecular CRNs always admit the
conservation law $\bm{\ell}^{ml}$ with $\ell^{ml}_{\alpha}=1$ for $\alpha\in
Z_{l}$ and $\ell^{ml}_{\alpha}=0$ for $\alpha\in Y_{h}$, representing the mass
conservation law for the $Z_{l}$ species, and ii) $\bm{\ell}^{ml}$ is the only
conservation law with null entries for the $Y_{h}$ species. Thus, the matrix
$\mathbb{L}$ whose rows are the conservation laws can always be written as

$\mathbb{L}=\begin{pmatrix}\bm{\ell}^{ml}\\\
\tilde{\mathbb{L}}\end{pmatrix}=\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\bm{\ell}^{ml}_{l_{\mathrm{}}}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\bm{0}\\\$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\tilde{\mathbb{L}}_{l}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\tilde{\mathbb{L}}_{h}$\hfil\kern
5.0pt\crcr}}}}\right)$}}\;.$ (59)

Note that there is at least one conservation law involving the $Y_{h}$
species, i.e., the mass conservation law of the whole CRN. Furthermore, the
conservation laws involving the $Y_{h}$ species are always broken since the
$Y_{h}$ species undergo concentration control.

From a dynamical standpoint, mass action kinetics (3) implies that the
reaction currents are linear functions of the concentrations
$\bm{z}_{l}=(\dots,[\alpha],\dots)^{\intercal}_{\alpha\in Z_{l}}$ since the
constant concentrations
$\bm{y}_{h}=(\dots,[\alpha],\dots)^{\intercal}_{\alpha\in Y_{h}}$ can be
absorbed (i.e., hidden) in the kinetic constants. Indeed, the reaction
currents $\bm{j}$ can be written as $\bm{j}=\Gamma\bm{z}_{l}$ with the
$\bm{y}_{h}$-dependent entries of the matrix $\Gamma$ reading

$\Gamma_{\rho,\alpha}=\begin{cases}\bar{k}_{+\rho}\equiv
k_{+\rho}{\displaystyle\prod_{\beta\in
Y_{h}}}[\beta]^{\nu_{\beta,+\rho}}~{}~{}\text{if}~{}{\nu_{\alpha,+\rho}=1}\,,\\\
-\bar{k}_{-\rho}\equiv-k_{-\rho}{\displaystyle\prod_{\beta\in
Y_{h}}}[\beta]^{\nu_{\beta,-\rho}}~{}\text{if}~{}{\nu_{\alpha,-\rho}=1}\,,\\\
0~{}\text{else}\,.\end{cases}$ (60)

We thus call the species $Z_{l}$ and $Y_{h}$ dynamically linear and hidden,
respectively.

###  IV.1 Pseudo-Closed Setup

When all the $Z_{l}$ species are internal, the rate equation (2) for the
concentrations $\bm{z}_{l}$ becomes

$d_{t}\bm{z}_{l}=\mathbb{W}\bm{z}_{l}\,,$ (61)

where $\mathbb{W}\equiv\mathbb{S}^{l}\Gamma$ is an irreducible rate matrix
whose diagonal elements are given by Eq. (29) after replacing the kinetic
constants $\\{k_{\pm\rho}\\}$ with $\\{\bar{k}_{\pm\rho}\\}$ defined in Eq.
(60). Since Eq. (61) is similar to the rate equation (28) of closed
unimolecular CRNs, we call this setup pseudo-closed. However, $\mathbb{W}$ in
Eq. (61) is, in general, nondetailed balanced due to the contribution of the
hidden species in $\\{\bar{k}_{\pm\rho}\\}$ unlike $\mathbb{W}$ in Eq. (28).
We consider hereafter only non-detailed balanced rate matrices $\mathbb{W}$ as
the case of detail balanced rate matrices $\mathbb{W}$ is mathematically
equivalent to unimolecular CRNs and is thus already examined in Sec. III.1.

Equation (61) together with the Perron Frobenius theorem implies that
$\bm{z}_{l}$ always reaches a steady state (as $\mathbb{W}$ admits one zero
eigenvalue while all the other eigenvalues have negative real parts). By
denoting the eigenvector of $\mathbb{W}$ corresponding to the zero eigenvalue
by $\bm{\pi}=(\dots,\pi_{\alpha},\dots)_{\alpha\in Z_{l}}^{\intercal}$, the
steady state of Eq. (61) reads

$\bm{z}^{\text{ss}}_{l}=L^{ml}(0)\bm{\pi}=(\dots,[\alpha]_{\text{ss}}({L^{ml}(0))},\dots)^{\intercal}\,,$
(62)

with $L^{ml}(0)=\bm{\ell}^{ml}_{l}\cdot\bm{z}_{l}$. Unlike
$\bm{z}_{\text{eq}}$ in Eq. (33), $\bm{z}^{\text{ss}}_{l}$ is a nonequilibrium
steady state since $\mathbb{W}$ is not detailed balanced. The corresponding
steady-state reaction currents, labeled $j_{\rho}^{\mathrm{ss}}$ for the
$\rho$ reaction interconverting $\beta$ into $\alpha$, are given by

$j_{\rho}^{\mathrm{ss}}=\mathbb{W}^{(\rho)}_{\alpha,\beta}[\beta]_{\text{ss}}-\mathbb{W}^{(\rho)}_{\beta,\alpha}[\alpha]_{\text{ss}}=L^{ml}(0)\bar{j}_{\rho}\,,$
(63)

where we used Eq. (62) and defined
$\bar{j}_{\rho}\equiv\mathbb{W}^{(\rho)}_{\alpha,\beta}\pi_{\beta}-\mathbb{W}^{(\rho)}_{\beta,\alpha}\pi_{\alpha}$.

We now briefly discuss the thermodynamics of the nonequilibrium steady state
$\bm{z}_{\mathrm{ss}}$ in Eq. (62) as it will be also used in Subs. IV.2. By
plugging Eq. (63) in Eq (18) and using
$\mathbb{S}^{l}\bm{j}^{\mathrm{ss}}=\mathbb{W}\bm{z}^{\text{ss}}_{l}=\bm{0}$
(with
$\bm{j}^{\mathrm{ss}}=(\dots,j_{\rho}^{\mathrm{ss}},\dots)^{\intercal}=L^{ml}(0)\bar{\bm{j}}$),
the steady-state EPR reads

$T\dot{\Sigma}^{\text{ss}}=-L^{ml}(0)\bm{\mu}_{h}\mathbb{S}^{h}\bar{\bm{j}}\,,$
(64)

where $\bm{\mu}_{h}$ is the vector of chemical potentials of the hidden
species. By identifying the potential hidden species $Y_{ph}$ and the force
hidden species $Y_{fh}$ (as explained in general in Sec. II.3), and by
applying this splitting to the matrix of broken conservation laws in Eq. (59),
i.e.,
$\tilde{\mathbb{L}}=(\tilde{\mathbb{L}}_{l}\,,\tilde{\mathbb{L}}_{{fh}}\,,\tilde{\mathbb{L}}_{{ph}})$,
and to the vector of chemical potentials, i.e.,
$\bm{\mu}_{{h}}=(\bm{\mu}_{{fh}}\,,\bm{\mu}_{{ph}})^{\intercal}$, the
nonconservative work rate (22) reads

$\begin{split}\dot{w}_{\text{nc}}^{\text{ss}}&=-L^{ml}(0)\big{(}\bm{\mu}_{h}-\bm{\mu}_{{ph}}\big{(}\tilde{\mathbb{L}}_{{ph}}\big{)}^{-1}\tilde{\mathbb{L}}_{h}\big{)}\cdot(\mathbb{S}^{h}\bar{\bm{j}})\\\
&=-L^{ml}(0)~{}(\mathcal{F}_{h}\cdot\mathbb{S}^{h}\bar{\bm{j}})\,,\end{split}$
(65)

where we also used $\bm{I}^{Y_{h}}=-L^{ml}(0)\mathbb{S}^{h}\bar{\bm{j}}$ and
Eq. (63). Equation (65) defines the nonconservative forces $\mathcal{F}_{h}$.
Finally, note that $\dot{w}_{\text{nc}}^{\text{ss}}=T\dot{\Sigma}^{\text{ss}}$
in agreement with the second law Eq. (23) (this can be directly verified by
using
$\tilde{\mathbb{L}}_{h}\cdot\mathbb{S}^{h}\bar{\bm{j}}=-\tilde{\mathbb{L}}_{l}\cdot\mathbb{S}^{l}\bar{\bm{j}}=0$)
which implies that
$\mathcal{F}_{h}\cdot\mathbb{S}^{h}\bar{\bm{j}}=\bm{\mu}_{h}\mathbb{S}^{h}\bar{\bm{j}}$.

###  IV.2 Fully-Open setup

We now consider the case where some of the $Z_{l}$ species are chemostatted.
The set $Z_{l}$ splits into the internal $X_{l}$ and chemostatted $Y_{l}$
dynamically linear species.

From a topological standpoint, chemostatting the $Z_{l}$ species breaks the
conservation law $\bm{\ell}^{ml}$. Thus, the set of potential $Y_{p}$ (resp.
force $Y_{f}$) species includes both dynamically linear and hidden species,
i.e., $Y_{p}={Y_{pl}}\cup Y_{ph}$ (resp. $Y_{f}={Y_{fl}}\cup Y_{fh}$). The
corresponding broken conservation laws and moieties are derived in Appendix C.

####  IV.2.1 Dynamics

The rate equation (61) for fully open pseudo-unimolecular CRNs under any
chemostatting procedure can be written in the form of Eq. (34). For flux
control (resp. mixed control), Eq. (34) follows from the definition in Eq. (7)
(resp. Eq. (8)), and by mapping the vector of dynamically linear
concentrations $\bm{z}_{l}$ into $\bm{a}$, the nondetailed balanced rate
matrix $\mathbb{W}$ (resp. the matrix difference $\mathbb{W}-\mathbb{D}$) into
$\mathbb{V}$, and $\tilde{\bm{I}}$ into $\bar{\bm{I}}_{a}$.

Furthermore, by following the same reasoning as in Sec.III.2, Eq. (34) can
also represent concentration control by mapping $\bm{x}$ into $\bm{a}$,
$\hat{\mathbb{W}}-\mathbb{D}$ into $\mathbb{V}$, and $\bar{\bm{I}}$ into
$\bar{\bm{I}}_{a}$. We emphasize that the key difference between the dynamics
(and thermodynamics) of open unimolecular CRNs and fully open pseudo-
unimolecular CRNs is that the matrices $\mathbb{W}$ and $\hat{\mathbb{W}}$ are
detailed balanced in the former case and nondetailed balanced in the latter
case (see the discussion below Eq. (61)).

####  IV.2.2 No Growth Under Mixed and Concentration Control

The dynamics of fully open pseudo-unimolecular CRNs under mixed control (resp.
concentration control) is given by Eq. (34) with
$\mathbb{V}=\mathbb{W}-\mathbb{D}$ (resp.
${\mathbb{V}=}\hat{\mathbb{W}}-\mathbb{D}$), $\bm{a}=\bm{z}_{l}$ (resp.
$\bm{a}=\bm{x}_{l}$), and $\bar{\bm{I}}_{a}=\tilde{\bm{I}}$ (resp.
$\bar{\bm{I}}_{a}=\bar{\bm{I}}$). According to Appendix A, all the eigenvalues
of $\mathbb{V}$ have negative real parts and $\mathbb{V}$ is invertible. Thus,
in the long time limit, the vector $\bm{a}(t)$ relaxes to the steady state
$-\mathbb{V}^{-1}\bar{\bm{I}}_{a}$ (see Appendix B.1) implying that open
pseudo-unimolecular CRNs do not grow if the dynamically linear species are
chemostatted via mixed control or concentration control, similarly to
unimolecular CRNs under mixed and concentration control (see Sec. III.3).

####  IV.2.3 Growth Under Flux Control

The dynamics of fully open pseudo-unimolecular CRNs under flux control is
given by Eq. (34) with $\mathbb{V}=\mathbb{W}$, $\bm{a}=\bm{z}$ and
$\bar{\bm{I}}_{a}=\tilde{\bm{I}}$. We assume here for simplicity that
$\mathbb{W}$ is diagonalizable but the same dynamic and thermodynamic behavior
emerges even for nondiagonalizable $\mathbb{W}$ as shown in Appendix B.2. By
exactly following the same approach as in Sec.III.4.1, we obtain that the mass
density of the dynamically linear species $L^{ml}(t)$ reads

$L^{ml}(t)=L^{ml}(0)+\left(\sum_{\alpha\in
Z_{l}}\tilde{I}_{\alpha}\right)t\,,$ (66)

which is of the form of Eq. (43) and that the concentrations in the long time
limit read

$\displaystyle[\alpha](t)$
$\displaystyle=\mathbb{P}_{\alpha,1}L^{ml}(t)+\sum_{\beta>1}\mathbb{P}_{\alpha,\beta}\left(\frac{(\mathbb{P}^{-1}\tilde{\bm{I}})_{\beta}}{-\lambda_{\beta}}\right)\,,$
(67a) $\displaystyle=\pi_{\alpha}L^{ml}(t)+c_{\alpha}{(\tilde{\bm{I}})}\,,$
(67b)

where $\mathbb{P}$ is the invertible matrix of eigenvectors of $\mathbb{W}$
written by ordering the eigenvalues as $\lambda_{1}=0$ and
$\text{Re}(\lambda_{\alpha})<0$ for $\alpha>1$, and the convention
$\mathbb{P}_{\alpha,1}=\pi_{\alpha}$ and $\mathbb{P}^{-1}_{\alpha,1}=1$. The
time-independent contributions on the right-hand side of Eq. (67a) are
collected in $c_{\alpha}{(\tilde{\bm{I}})}$ in Eq. (67b). By using Eq. (62),
we identify $\pi_{\alpha}L^{ml}(t)$ as the nonequilibrium steady state
$[\alpha]_{\text{ss}}(L^{ml}(t))$ of the corresponding pseudo-closed CRN with
mass density $L^{ml}(t)$. Equations (66) and (62) show that the mass density
of the dynamically linear species, $L^{ml}(t)$, as well as the steady state
concentrations $[\alpha]_{\text{ss}}\left(L^{ml}(t)\right)$ grow linearly in
time. This together with Eq. (67b) implies that the instantaneous
concentrations $[\alpha](t)$ also grow linearly in time. Since the relative
(resp. absolute) difference between $[\alpha](t)$ and the nonequilibrium
steady state $[\alpha]_{\text{ss}}\left(L^{ml}(t)\right)$ decreases in time
(resp. is constant) (see Eq. (67b),), we call this type of growth
nonequilibrium growth.

We now examine the corresponding reaction currents. By replacing the steady-
state concentrations in Eq. (63) with the concentrations in Eq. (67b), the
reaction current $j_{\rho}(t)$ in the long time limit reads

$\begin{split}j_{\rho}(t)&=\bar{j}_{\rho}L^{ml}(t)+\left(\mathbb{W}^{(\rho)}_{\alpha\beta}c_{\beta}-\mathbb{W}^{(\rho)}_{\beta\alpha}c_{\alpha}\right)\\\
&=\bar{j}_{\rho}L^{ml}(t)+J_{\rho}=\mathcal{O}(t)\,,\end{split}$ (68)

where
$J_{\rho}\equiv\mathbb{W}^{(\rho)}_{\alpha\beta}c_{\beta}-\mathbb{W}^{(\rho)}_{\beta\alpha}c_{\alpha}$.
We note that in nonequilibrium growth both the dynamically linear
concentrations and reaction currents grow linearly in time in contrast to
equilibrium growth where only the former grow in time (compare Eqs.(67b) and
(68) with Eqs (45b) and (46), respectively).

####  IV.2.4 Thermodynamics Under Flux Control

We now analyze the thermodynamics of nonequilibrium growth. In the long time
limit, by using Eq. (67b) in Eq. (14), the chemical potentials of the
dynamically linear species are given by

$\begin{split}\mu_{\alpha}=~{}&RT~{}\text{ln}\left(L^{ml}(t)\right)+RT~{}\text{ln}\left(\pi_{\alpha}e^{\frac{\mu^{0}_{\alpha}}{RT}}\right)\\\
&+RT~{}\text{ln}\left(1+\frac{c_{\alpha}}{\pi_{\alpha}L^{ml}(t)}\right)\,,\end{split}$
(69)

and are logarithmically increasing in time. The corresponding vector of
chemical potentials $\bm{\mu}_{l}$ reads

$\begin{split}\bm{\mu}_{l}=~{}&RT~{}\text{ln}\left(L^{ml}(t)\right)\bm{\ell}^{ml}_{l}+RT~{}\text{ln}\left(\bm{\pi}\exp{\frac{\bm{\mu^{0}}}{RT}}\right)\\\
&+RT\frac{1}{L^{ml}(t)}\left(\frac{\bm{c}}{\bm{\pi}}\right)\,,\end{split}$
(70)

where we used $\text{ln}(1+x)\approx x$ for $x\ll 1$, and introduced
$\text{ln}\left(\bm{\pi}\exp{{\bm{(}{\bm{\mu}^{0}}}/{RT}})\right)=\left(\dots,~{}\text{ln}\left(\pi_{\alpha}\exp({{\mu^{0}_{\alpha}}/{RT}})\right),\dots\right)^{\intercal}$
and ${\bm{c}}/{\bm{\pi}}=(\dots,c_{\alpha}/\pi_{\alpha},\dots)^{\intercal}$.
On the other hand, the chemical potentials of the hidden species
$\bm{\mu}_{h}$ are constant as their concentrations are kept fixed.

By plugging Eqs. (70) and (68) in Eq. (18) and using
$\mathbb{S}^{l}\bm{j}=\mathbb{W}\bm{z}_{l}=\mathbb{W}\bm{c}$ (see Eq. (67b)),
we find that the EPR at long times can be written as

$\begin{split}T\dot{\Sigma}=&\overbrace{L^{ml}(t)(-\bm{\mu}_{h}\mathbb{S}^{h}\bar{\bm{j}})}^{=T\dot{\Sigma}^{\text{ss}}(L^{ml}(t))}-\bm{\mu}_{h}\mathbb{S}^{h}{\bm{J}}\\\
&-RT~{}\text{ln}\left(\bm{\pi}\exp{\frac{\bm{\mu^{0}}}{RT}}\right)\cdot\mathbb{W}\bm{c}=\mathcal{O}(t)\,,\\\
\end{split}$ (71)

where $\bm{J}=\big{(}\dots,J_{\rho},\dots\big{)}^{\intercal}$ and we used Eq.
(64) to identify the leading order term on the right hand side as the steady
state EPR to which the corresponding pseudo-closed CRN with mass density
$L^{ml}(t)$ would relax to. Equation (71) implies that nonequilibrium growth
is a dissipative process with a monotonically increasing EPR resulting from
the growing mass density of the dynamically linear species in Eq. (66).

By putting Eq. (70), $\bm{I}^{Y_{l}}=\tilde{\bm{I}}$ and
$\bm{I}^{Y_{h}}=-\mathbb{S}^{h}\bm{j}$ in Eq. (19), we find that the chemical
work rate increases linearly in the long time limit:

$\displaystyle\dot{w}_{\text{c}}=$
$\displaystyle\overbrace{L^{ml}(t)(-\bm{\mu}_{h}\mathbb{S}^{h}\bar{\bm{j}})-\bm{\mu}_{h}\mathbb{S}^{h}{\bm{J}}\vphantom{RT~{}\text{ln}\Big{(}\bm{\pi}\exp{\frac{\bm{\mu^{0}}}{RT}}\Big{)}\cdot\tilde{\bm{I}}}}^{=\dot{w}_{\mathrm{h}}}$
(72)
$\displaystyle+\overbrace{{RT~{}\text{ln}\big{(}L^{ml}(t)\big{)}(\bm{\ell}^{ml}_{l}\cdot\tilde{\bm{I}})+RT~{}\text{ln}\Big{(}\bm{\pi}\exp{\frac{\bm{\mu^{0}}}{RT}}\Big{)}\cdot\tilde{\bm{I}}}}^{=\dot{w}_{\mathrm{l}}}=\mathcal{O}(t)\,,$

where we also split the chemical work rate into the hidden work rate
$\dot{w}_{\mathrm{h}}\equiv\bm{\mu}_{h}\cdot\bm{I}^{Y_{h}}$ and the linear
work rate $\dot{w}_{\mathrm{l}}\equiv\bm{\mu}_{l}\cdot\tilde{\bm{I}}$. The
hidden (resp. linear) work rate quantifies the contribution of the
chemostatted hidden (resp. linear) species to the total chemical work rate.
Furthermore, by using Eq. (70) and $d_{t}\bm{y}_{h}=\bm{0}$ in Eq. (17), we
find that $d_{t}G$ increases logarithmically in the long time limit:

$\begin{split}d_{t}G&=RT~{}\text{ln}\left(L^{ml}(t)\right)(\bm{\ell}^{ml}_{l}\cdot\tilde{\bm{I}})\\\
&{}~{}~{}~{}+RT~{}\text{ln}\left(\bm{\pi}\exp{\frac{\bm{\mu^{0}}}{RT}}\right)\cdot
d_{t}\bm{z}_{l}=\mathcal{O}(\text{ln}(t))\,,\end{split}$ (73)

where we used $\bm{\ell}^{ml}_{l}\cdot\tilde{\bm{I}}=\bm{\ell}^{ml}_{l}\cdot
d_{t}\bm{z}_{l}$ (see Eqs. (67b) and (66)). Thus, by using Eqs. (73) and (72)
in Eq. (53), we find that the efficiency of nonequilibrium growth Eq.(53) goes
to zero in the long time limit,

$\lim_{t\to\infty}\eta=\lim_{t\to\infty}\frac{RT~{}\text{ln}\left(L^{ml}(t)\right)}{L^{ml}(t)}=0\,,$
(74)

implying that nonequilibrium growth is a thermodynamically irreversible
process in contrast to equilibrium growth (see Subs. III.4.2).

We now turn to the moiety work in Eq. (20) and the nonconservative work in Eq.
(22). Substituting Eqs. (69) and (172) in Eq. (20), the moiety work rate reads

$\displaystyle\dot{w}_{\text{m}}={}$ $\displaystyle
RT~{}\text{ln}\left(L^{ml}(t)\right)(\bm{\ell}^{ml}_{l}\cdot\tilde{\bm{I}})$
(75)
$\displaystyle+RT~{}\text{ln}\left(\pi_{pl}\exp{\frac{\mu^{0}_{pl}}{RT}}\right)(\bm{\ell}^{ml}_{l}\cdot\tilde{\bm{I}})$
$\displaystyle+\bm{\mu}_{ph}\cdot\left(\overline{\mathbb{L}}_{X}\bm{\pi}_{X}+\overline{\mathbb{L}}_{fl}\bm{\pi}_{fl}\right)(\bm{\ell}^{ml}_{l}\cdot\tilde{\bm{I}})=\mathcal{O}(\text{ln}(t))\,,$

where we used
$d_{t}\bm{x}=\bm{\pi}_{X}(\bm{\ell}^{ml}_{l}\cdot\tilde{\bm{I}})$ and
$d_{t}\bm{y}_{fl}=\bm{\pi}_{fl}(\bm{\ell}^{ml}_{l}\cdot\tilde{\bm{I}})$ from
Eq. (67b) and the matrices $\overline{\mathbb{L}}_{X}$,
$\overline{\mathbb{L}}_{fl}$ given in Eq. (172). Note that the leading order
term in Eq. (75) is identical to the leading order terms of both
$\dot{w}_{\mathrm{l}}$ in Eq. (72) and $d_{t}G$ in Eq. (73). Namely,

$\dot{w}_{\mathrm{l}}\sim\dot{w}_{\text{m}}\sim
d_{t}G=\mathcal{O}(\text{ln}(t))\,,$ (76)

which implies that the linear work, as well as the moiety work, is mainly
converted into free energy. This results from the fact that the concentrations
of the $Y_{h}$ species are constant in time and, therefore, they do not
contribute to $d_{t}G$ (see Eq. (17)). Only the concentrations of the
dynamically linear species contribute to $d_{t}G$, implying that the linear
work is mainly converted into free energy. Similarly, as the leading order
contribution to $\dot{w}_{\text{m}}$ is the work done in changing $L^{ml}(t)$
(see Eq. (75)), which does not involve any $Y_{h}$ species, the moiety work
rate and the linear work rate scale similarly in time.

By using Eq. (170) and Eq. (171b), the nonconservative forces defined in Eq.
(22) become

$\mathcal{F}_{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.0pt\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&5.0pt\hfil$\scriptstyle 0~{}~{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\bm{\mu}_{fl}-\mu_{pl}(\ell^{ml}_{pl})^{-1}\bm{\ell}^{ml}_{fl}-\bm{\mu}_{ph}\cdot\overline{\mathbb{L}}_{fl}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle~{}~{}\bm{0}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle~{}~{}\mathcal{F}_{h}$\hfil\kern
5.0pt\crcr}}}}\right)$}}$ (77)

where $\mathcal{F}_{h}$ was defined in Eq. (65). By using the Eqs. (68), (69)
and (77) in Eq. (22), the nonconservative work rate reads

$\displaystyle\dot{w}_{\text{nc}}=$
$\displaystyle\overbrace{-L^{ml}(t)~{}(\mathcal{F}_{h}\cdot\mathbb{S}^{h}\bar{\bm{j}})}^{=\dot{w}_{\text{nc}}^{\text{ss}}(L^{ml}(t))}-\mathcal{F}_{h}\cdot\mathbb{S}^{h}{\bm{J}}$
(78)
$\displaystyle+\left(RT\text{ln}\left(\frac{\bm{\pi}_{fl}}{\pi_{pl}}\exp{\frac{\bm{\mu^{0}_{fl}}-\mu^{0}_{pl}}{RT}}\right)-\bm{\mu}_{ph}\cdot\overline{\mathbb{L}}_{fl}\right)\cdot\tilde{\bm{I}}=\mathcal{O}(t)\,,$

where

$\displaystyle
RT\text{ln}(\bm{\pi}_{fl}/{\pi_{pl}}\exp{((\bm{\mu}^{0}_{fl}-\mu^{0}_{pl})/{RT})})$
(79)
$\displaystyle=(\dots,RT\text{ln}(\pi_{\alpha}/{\pi_{pl}}\exp{((\mu^{0}_{\alpha}-\mu^{0}_{pl})/{RT})}),\dots)^{\intercal}$

for ${\alpha\in Y_{fl}}$ and we used Eq. (65) to identify
$\dot{w}_{\text{nc}}^{\text{ss}}(L^{ml}(t))$, namely, the nonconservative work
at the steady state to which the corresponding pseudo-closed pseudo-
unimolecular CRN with mass density $L^{ml}(t)$ would relax. Note that the
leading order terms of $\dot{w}_{\text{nc}}$ in Eq. (78), the EPR in Eqs. (71)
and $\dot{w}_{\mathrm{h}}$ in Eq. (72) are equal. Namely,

$\dot{w}_{\mathrm{h}}\sim\dot{w}_{\text{nc}}\sim
T\dot{\Sigma}=\mathcal{O}(t)\,,$ (80)

which implies that the hidden work, as well as the nonconservative work rate,
is mainly dissipated.

###  IV.3 Summary

From a dynamical standpoint, growth in fully open pseudo-unimolecular CRNs
resembles growth in unimolecular CRNs. Like in unimolecular CRNs, it can occur
only under flux control; the mass density (of dynamically linear species)
increases linearly in time; the instantaneous concentrations approach the
steady-state concentrations defined by the growing mass density (compare Eqs.
(67b) and (45b)); the concentrations of all (dynamically linear) species scale
linearly in time. However, the steady state is a nonequilibrium (resp.
equilibrium) one for pseudo-unimolecular CRNs (resp. unimolecular CRNs).

From a thermodynamic standpoint, the work done in increasing the mass density
and the change in free energy scale similarly in time like in unimolecular
CRNs (compare Eqs. (76) and (56)). Furthermore, the nonconservative work and
the EPR scale similarly in time like in unimolecular CRNs (compare Eqs. (80)
and (57)). However, pseudo-unimolecular CRNs grow out of equilibrium with the
EPR linearly increasing in time (see Eq. (71)), whereas unimolecular CRNs grow
close to equilibrium with the EPR monotonically decaying in time (see Eq.
(49)).

Unlike for unimolecular CRNs, the efficiency of growth (74) goes to zero for
pseudo-unimolecular CRNs, implying that growth is a thermodynamically
irreversible process.

##  V Multimolecular CRNs

CRNs are said to be multimolecular when at least one chemical reaction
involves two or more species that do not undergo concentration control.
Because of mass action kinetics (3), the corresponding rate equation (2) is
nonlinear and does not admit in genreal analytical solutions. Hence, we will
focus on the four multimolecular CRNs defined in Fig. 1 and use both
analytical methods and numerical simulations to investigate growth under the
three different chemostatting procedures. These CRNs have been chosen because
of their different topological properties. The autocatalytic CRN has one
conservation law corresponding to the mass density,
$L^{m}=[\ch{E}]+2[\ch{E}^{*}]+[\ch{S}]+2[\ch{P}]$. The Michaelis Menten CRN
and the cyclic Michaelis Menten CRN have two conservation laws corresponding
to the total concentration of the enzyme $L^{\ch{E}}=[\ch{E}]+[\ch{ES}]$ and
total concentration of the substrate $L^{\ch{S}}=[\ch{ES}]+[\ch{S}]+[\ch{P}]$.
The minimal metabolic CRN has three conservation laws corresponding to the
total concentration of enzyme,
$L^{\ch{E}}=[\ch{E}]+[\ch{EF}]+[\ch{EW}]+[\ch{E}^{*}]$, the total
concentration of fuel,
$L^{\ch{F}}=[\ch{EF}]+[\ch{EW}]+[\ch{E}^{*}]+[\ch{F}]+[\ch{W}]$, and the total
concentration of substrate, $L^{\ch{S}}=[\ch{E}^{*}]+[\ch{S}]+[\ch{P}]$.

For all numerical simulations in the following, we rescale time,
concentrations, energy/time respectively by i) $1/k_{+1}$, $k_{+1}/k_{+3}$,
$RTk_{+1}^{2}/k_{+3}$ for the autocatalytic CRN, and by ii) $1/k_{-1}$,
$k_{-1}/k_{+1}$, $RTk_{-1}^{2}/k_{+1}$ for the Michaelis Menten, the Cyclic
Michaelis Menten and the Minimal Metabolic CRNs. Furthermore, we disregard
initial transients and focus on long-time dynamics only.

Figure 1:  The CRNs studied in Sec. V, their stoichiometric matrices (with the
horizontal line splitting $\mathbb{S}$ into $\mathbb{S}^{X}$ and
$\mathbb{S}^{Y}$), and their conservation laws. i) The autocatalytic CRN
represents the interconversion of the substrate S into the product P promoted
by the autocatalytic reaction between the enzymes E and $\ch{E}^{*}$; ii) The
Michaelis Menten CRN represents the interconversion of the substrate S into
the product P promoted by the enzymes E and ES; iii) The Cyclic Michaelis
Menten CRN represents the same enzymatic mechanism as the Michaelis Menten
CRN, but also includes a reaction directly interconverting S into P; iv) The
minimal metabolic CRN represents the enzymatic interconversion of the
substrate S into the product P powered by the interconversion of the fuel F
into the waste W and promoted by the enzymes E, EF, EW, and $\ch{E}^{*}$.

###  V.1 No Growth Under Concentration Control

Growth in open CRNs under concentration control has been studied in the
mathematical literature. It has been formally proven that the rate equation
(2) (with fluxes satisfying mass action kinetics) does not admit solutions
with unbounded concentrations for weakly reversible [31] single linkage class
[32] CRNs, as well as strongly endotactic [33] CRNs. This has also been
conjectured to hold for any weakly reversible CRNs [31], but a general proof
is still missing.

We now confirm this conjecture for the autocatalytic CRN in Fig. 1. The other
three CRNs are pseudo-unimolecular CRNs under concentration control and,
consequently, do not grow based on the results of Sec. IV.1.

We start by showing the existence of a steady state for any concentration of
the chemostatted species $\ch{S}$ and $\ch{P}$. According to Eq. (5), the
concentrations of the internal species $[\ch{E}]$ and $[\ch{E}^{*}]$ follow

$\displaystyle d_{t}[\ch{E}]$
$\displaystyle=\bar{k}_{+1}+2k_{+3}[\ch{E}][\ch{E}^{*}]-{k}_{-1}[\ch{E}]-2k_{-3}[\ch{E}]^{3}\,,$
(81a) $\displaystyle d_{t}[\ch{E}^{*}]$
$\displaystyle=\bar{k}_{+2}+k_{-3}[\ch{E}]^{3}-k_{+3}[\ch{E}^{*}][\ch{E}]-k_{-2}[\ch{E}^{*}]\,,$
(81b)

where we introduced the effective kinetic constants $\bar{k}_{+1}\equiv
k_{+1}[\ch{S}]$ and $\bar{k}_{+2}\equiv k_{+2}[\ch{P}]$ since $[\ch{S}]$ and
$[\ch{P}]$ are constant in time. Correspondingly, the mass density
$L^{m}=[\ch{E}]+2[\ch{E}^{*}]+[\ch{S}]+2[\ch{P}]$ follows

$d_{t}L^{m}=\bar{k}_{+1}-k_{-1}[\ch{E}]+2\bar{k}_{+2}-2k_{-2}[\ch{E}^{*}]\,.$
(82)

This implies that the steady-state concentrations $[\ch{E}]_{\text{ss}}$ and
$[\ch{E}^{*}]_{\text{ss}}$ satisfy

$[\ch{E}^{*}]_{\text{ss}}=\frac{\bar{k}_{+1}+2\bar{k}_{+2}}{2{k}_{-2}}-\frac{{k}_{-1}}{2k_{-2}}[\ch{E}]_{\text{ss}}\,.$
(83)

By substituting Eq. (83) into Eq. (81b), we obtain the cubic equation

$\displaystyle h([\ch{E}]_{\text{ss}})\equiv$
$\displaystyle[\ch{E}]_{\text{ss}}^{3}+\left(\frac{k_{-1}k_{+3}}{2k_{-2}k_{-3}}\right)[\ch{E}]_{\text{ss}}^{2}$
(84)
$\displaystyle+\left(\frac{k_{-1}k_{-2}-\bar{k}_{+1}k_{+3}-2\bar{k}_{+2}k_{+3}}{2k_{-2}k_{-3}}\right)[\ch{E}]_{\text{ss}}-\frac{\bar{k}_{+1}}{2k_{-3}}=0\,,$

which admits a unique positive root, i.e., there is a unique steady state
concentration $[\ch{E}]_{\text{ss}}$, because of the Descartes rule of sign
[34]. Furthermore, the steady state concentration $[\ch{E}]_{\text{ss}}$
satisfies

$0<[\ch{E}]_{\text{ss}}<\frac{\bar{k}_{+1}+2\bar{k}_{+2}}{k_{-1}}\,,$ (85)

since $h(0)<0$ and $h((\bar{k}_{+1}+2\bar{k}_{+2})/k_{-1})>0$. This, together
with Eq. (83), implies that $[\ch{E}^{*}]_{\text{ss}}>0$ and, therefore, the
autocatalytic CRN has a well-defined steady state for every value of
$[\ch{S}]$ and $[\ch{P}]$.

We then analyze the stability of the steady state. By using Eqs. (81a) and
(81b), the Jacobian matrix $\mathbb{J}$ of the steady state reads

$\mathbb{J}=\begin{pmatrix}2{k}_{+3}[\ch{E}^{*}]_{\text{ss}}-k_{-1}-6k_{-3}[\ch{E}]_{\text{ss}}^{2}&2{k}_{+3}[\ch{E}]_{\text{ss}}\\\
3k_{-3}[\ch{E}]^{2}_{\text{ss}}-{k}_{+3}[\ch{E}^{*}]_{\text{ss}}&-k_{-2}-k_{+3}[\ch{E}]_{\text{ss}}\end{pmatrix}\,.$
(86)

Its determinant and trace can be written as

$\displaystyle\mathrm{det}(\mathbb{J})=\frac{4k_{-2}k_{-3}[\ch{E}]_{\text{ss}}^{3}+k_{-1}k_{+3}[\ch{E}]_{\text{ss}}^{2}+\bar{k}_{+1}k_{-2}}{[\ch{E}]_{\text{ss}}}>0\,,$
(87a)
$\displaystyle\mathrm{tr}(\mathbb{J})=-\bigg{(}4k_{-3}[\ch{E}]_{\text{ss}}^{2}+k_{+3}[\ch{E}]_{\text{ss}}+k_{-2}+\frac{\bar{k}_{+1}}{[\ch{E}]_{\text{ss}}}\bigg{)}<0\;.$
(87b)

This implies that $\mathbb{J}$ has only negative eigenvalues and, therefore,
the steady state is always locally stable for any values of $[\ch{S}]$ and
$[\ch{P}]$.

Finally, we numerically compute the dynamics of the autocalytic CRN. Figure 2a
shows the typical evolution the concentrations $[\ch{E}]$ and $[\ch{E}^{*}]$:
after a transient, they reach steady state. Figure 2b shows the typical
evolution of the mass density for different initial conditions. Like the
concentrations, the mass density relaxes towards a constant value. We stress
that the same behavior is observed for different values of the kinetic
constants and of the concentrations of the chemostatted species.

Figure 2:  Dynamics of the autocatalytic CRN in Fig. 1 under concentration
control. Evolution of a) Concentrations for the specific initial condition
$[\ch{E}](0)=1.7$ and $[\ch{E}^{*}](0)=1.4$. b) Mass density for thirty
randomly-selected initial conditions. Here, $k_{\pm 1}=k_{\pm 2}=k_{\pm 3}=1$,
$[\ch{S}](0)=10$, and $[\ch{P}](0)=3.1$.

###  V.2 Growth Under Flux Control

Flux control always leads to growth when $\bm{\ell}^{m}\cdot\tilde{\bm{I}}>0$,
as already discussed in Sec. II.4. In this subsection, we analyze the dynamics
and thermodynamics of growth for the CRNs in Fig. 1 and compare them with the
growth of unimolecular and pseudo-unimolecular CRNs examined in Secs. III.4.1,
IV.2.3 and IV.2.4, respectively.

####  V.2.1 Autocatalytic CRN

We consider the autocatalytic CRN in Fig. 1 when $\ch{S}$ is injected with
flux $I_{\ch{S}}>0$ and $\ch{P}$ is extracted with flux $I_{\ch{P}}<0$ such
that the mass density grows at the rate
$d_{t}L^{m}=I_{\ch{S}}+2I_{\text{P}}>0$ (according to Eq. (13)).

Figure 3: Dynamics of the autocatalytic CRN under flux control. Evolution of
a) Concentrations for the initial condition $[\ch{E}](0)=1.3$,
$[\ch{E}^{*}]=0.1$, $[\ch{F}]=1.9$, and $[\ch{W}]=0.4$. b) Relative
difference,
$\lVert\bm{z}(t)-\bm{z}_{\text{eq}}(t)\rVert/\lVert\bm{z}_{\text{eq}}(t)\rVert$
(where the equilibrium is defined as the steady state to which the CRN would
relax if it was closed and with mass density $L^{m}(t)$) and, in the inset,
$[\ch{E}](t)$ and $[\ch{E}]_{\text{eq}}(t)$. Here, $k_{\pm 1}=k_{\pm 2}=k_{\pm
3}=1$, $I_{\ch{S}}=1.0$, and $I_{\ch{P}}=-0.1$.

Figure 3a shows a typical evolution of the concentrations, while Fig. 3b
compares it to the evolution of the equilibrium state to which the CRN would
relax if it was closed and with mass density $L^{m}(t)$, labeled
$\bm{z}_{\text{eq}}(L^{m}(t))$. We first observe that all concentrations grow.
Second, the concentrations grow at different rates, in contrast to what
happens in unimolecular and pseudo-unimolecular CRNs under flux control (see
Eqs. (45b) and (67b)). In particular,
${[\ch{E}](t)\sim[\ch{S}](t)=\mathcal{O}(\sqrt{t})}$ and
${[\ch{E}^{*}](t)\sim[\ch{P}](t)=\mathcal{O}(t)}$. Finally, the concentrations
$\bm{z}(t)$ evolve close to the corresponding equilibrium concentrations
$\bm{z}_{\text{eq}}(L^{m}(t))$ similarly to unimolecular CRNs (see Eq. (45b)).
Indeed, the relative difference
$\lVert\bm{z}(t)-\bm{z}_{\text{eq}}(t)\rVert/\lVert\bm{z}_{\text{eq}}(t)\rVert$
decreases in time. Hence, the autocatalytic CRN under flux control undergoes
equilibrium growth similar to unimolecular CRNs under flux control (see Secs.
III.4).

This specific dynamics can be explained in terms of a time-scale separation
between the slowly evolving mass density and the rapidly evolving
concentrations of the chemical species. On the one hand, the former increases
at the constant rate $d_{t}L^{m}=I_{\ch{S}}+2I_{\ch{P}}>0$. On the other hand,
the latter evolves also because of the chemical reactions whose fluxes,
according to mass-action kinetics (3), increase when the concentrations
increase. Hence, on a time scale in which the mass density $L^{m}$ remains
almost constant, the concentrations $\bm{z}$ can dramatically change and relax
towards the equilibrium $\bm{z}_{\text{eq}}(L^{m}(t))$. This can be captured
mathematically by the ansatz

$\bm{z}(t)=\bm{z}_{\text{eq}}(L^{m}(t))+\bm{c}(t)\,,$ (88)

with $|{c}_{\alpha}(t)|\ll[\alpha]_{\text{eq}}(L^{m}(t))$, resembling Eq.
(45b).

The ansatz (88) allows us to predict the rate of growth of the concentrations
observed in Fig. 3. Indeed, according to mass action kinetics (3), the
equilibrium concentrations $\bm{z}_{\text{eq}}(L^{m}(t))$ must satisfy

$\begin{split}[\ch{E}]_{\text{eq}}(t)&=\frac{k_{+1}}{k_{-1}}[\ch{S}]_{\text{eq}}(t)=\sqrt{\frac{k_{+3}k_{+2}}{k_{-2}k_{-3}}[\ch{P}]_{\text{eq}}(t)}\\\
&=\sqrt{\frac{k_{+3}}{k_{-3}}[\ch{E}^{*}]_{\text{eq}}(t)}\,,\end{split}$ (89)

which, together with
$\bm{\ell}^{m}\cdot\bm{z}_{\text{eq}}(L^{m}(t))=L^{m}(t)$, leads to

$[\text{E}]^{2}_{\text{eq}}(t)+[\text{E}]_{\text{eq}}\underbrace{\left(\frac{1+\frac{k_{-1}}{k_{+1}}}{\frac{2k_{-3}}{k_{+3}}+\frac{2k_{-2}k_{-3}}{k_{+2}k_{+3}}}\right)}_{=a_{1}}-\underbrace{\left(\frac{L^{m}(t)}{\frac{2k_{-3}}{k_{+3}}+\frac{2k_{-2}k_{-3}}{k_{+2}k_{+3}}}\right)}_{=a_{2}(t)}=0\,$
(90)

and therefore,

$[\text{E}]_{\text{eq}}(t)=\sqrt{\left(\frac{a_{1}}{2}\right)^{2}+a_{2}(t)}-\left(\frac{a_{1}}{2}\right)\,.$
(91)

Since $a_{2}(t)\propto L^{m}(t)=\mathcal{O}(t)$, from Eqs. (89) and (91), we
obtain that
$[\ch{E}]_{\text{eq}}(t)\sim[\ch{S}]_{\text{eq}}(t)=\mathcal{O}(\sqrt{t})$
while $[\ch{P}]_{\text{eq}}(t)\sim[\ch{E}^{*}]_{\text{eq}}(t)=\mathcal{O}(t)$.

We now examine the corresponding thermodynamics. We start by considering the
EPR (18), $d_{t}G$ (17) and chemical work rate (19) whose typical evolution is
shown in Fig. 4a.

Figure 4: Thermodynamics of the autocatalytic CRN under flux control during
the dynamics in Fig. 3. Evolution of a) EPR, $d_{t}G$, chemical work rate, and
efficiency of growth. b) EPR, $d_{t}G$, nonconservative work rate, and the
moiety work rate. Here, $\mu^{0}_{\ch{E}}=1$, $\mu^{0}_{\ch{E^{*}}}=2$,
$\mu^{0}_{\ch{S}}=1$, $\mu^{0}_{\ch{P}}=2$.

The EPR continuously decreases until it vanishes. On the other hand, $d_{t}G$
and the chemical work increase in time, and they coincide in the long time
limit. This means that the chemical work is fully converted into Gibbs free
energy at long times and the efficiency of growth, defined as in Eq. (53),
goes to one (see the inset in Fig. 4a). We then turn to the nonconservative
work rate and the moiety work rate, defined in Eqs. (22) and (20),
respectively. Their typical evolution is shown in Fig. 4b (when $\ch{S}$ is
chosen as the potential species): in the long time limit, the nonconservative
work vanishes and the moiety work rate balances the increase of Gibbs free
energy.

This specific behavior of the thermodynamic quantities is similar to
unimolecular CRNs under flux control (see Sec. III.4) and is consistent with
the previous observation of equilibrium growth. It can further be explained
with the ansatz (88). Indeed, by using Eq. (88) with
$||\bm{c}(t)||\ll||\bm{z}_{\text{eq}}(L^{m}(t))||$ in Eq. (14), the chemical
potentials $\bm{\mu}(t)$ read

$\bm{\mu}(t)\approx\bm{\mu}_{\text{eq}}(t)+RT\frac{\bm{c}(t)}{\bm{z}_{\text{eq}}(L^{m}(t))}=\mathcal{O}(\text{ln}(t))\,,$
(92)

where
$\bm{c}(t)/\bm{z}_{\text{eq}}(L^{m}(t))=(\dots,c_{\alpha}(t)/[\alpha]_{\text{eq}}(L^{m}(t)),\dots)^{\intercal}$
and $\bm{\mu}_{\text{eq}}(t)=\bar{\mu}(t)\bm{\ell}^{m}$ with
$\bar{\mu}(t)=\mu^{0}_{\text{E}}+RT\text{ln}\left([\text{E}]_{\text{eq}}(t)\right)=\mathcal{O}(\text{ln}(t))$.
Then, by plugging Eq. (92) in Eq. (18), the EPR becomes

$\begin{split}T\dot{\Sigma}&=-RT\frac{\bm{c}(t)}{\bm{z}_{\text{eq}}(L^{m}(t))}\cdot\mathbb{S}\bm{j}\\\
&\approx-
RT\frac{\bm{c}(t)}{\bm{z}_{\text{eq}}(L^{m}(t))}\cdot\left(d_{t}\bm{z}_{\text{eq}}-\tilde{\bm{I}}\right)\,,\end{split}$
(93)

where we used Eq. (2) and $d_{t}\bm{z}\approx d_{t}\bm{z}_{\text{eq}}$. At
long times, EPR in Eq. (93) vanishes since i) $d_{t}\bm{z}_{\text{eq}}$
approaches a constant value (see Eqs. (89) and (91)), and ii)
$|{c}_{\alpha}(t)|\ll[\alpha]_{\text{eq}}(L^{m}(t))$. Furthermore, by plugging
Eq. (92) in Eqs. (17) and (19) and using
$|{c}_{\alpha}(t)|\ll[\alpha]_{\text{eq}}(L^{m}(t))$, $d_{t}G$ and the
chemical work rate read

$\displaystyle d_{t}G=$
$\displaystyle\bar{\mu}(t)(\bm{\ell}^{m}\cdot\tilde{\bm{I}})+RT\frac{\bm{c}(t)}{\bm{z}_{\text{eq}}(L^{m}(t))}\cdot({\tilde{\bm{I}}}{+}\mathbb{S}\bm{j})=\mathcal{O}(\text{ln}(t))$
(94) $\displaystyle\dot{w}_{c}=$
$\displaystyle\bar{\mu}(t)(\bm{\ell}^{m}\cdot\tilde{\bm{I}})+RT\frac{\bm{c}(t)}{\bm{z}_{\text{eq}}(L^{m}(t))}\cdot{\tilde{\bm{I}}}=\mathcal{O}(\text{ln}(t))\,,$
(95)

implying $d_{t}G\approx\dot{w}_{c}$ and $\eta\to 1$ in agreement with Fig. 4a.

Similarly, by using Eq. (92) in Eqs. (20) and (22), the nonconservative work
rate and the moiety work rate read

$\displaystyle\dot{w}_{\text{m}}$
$\displaystyle=\mu_{\ch{S}}d_{t}L^{m}=\bar{\mu}(t)(\bm{\ell}^{m}\cdot{\tilde{\bm{I}}})+RT\frac{c_{\ch{S}}(t)}{[\ch{S}]_{\text{eq}}(L^{m}(t))}(\bm{\ell}^{m}\cdot{\tilde{\bm{I}}})$
$\displaystyle=\mathcal{O}\left(\text{ln}(t)\right)\,,$ (96)
$\displaystyle\dot{w}_{\text{nc}}$
$\displaystyle=(\mu_{\ch{P}}-2\mu_{\ch{S}})I_{\ch{P}}$
$\displaystyle=RT\left(\frac{c_{\ch{P}}(t)}{[\ch{P}]_{\text{eq}}(L^{m}(t))}-\frac{2c_{\ch{S}}(t)}{[\ch{S}]_{\text{eq}}(L^{m}(t))}\right)I_{\ch{P}}\,,$
(97)

implying that over long times

$\displaystyle T\dot{\Sigma}\approx\dot{w}_{\text{nc}}\approx 0\,,$
$\displaystyle d_{t}G\approx\dot{w}_{\text{m}}=\mathcal{O}(\text{ln}(t))\,.$
(98)

####  V.2.2 Michaelis Menten CRN

We consider the Michaelis Menten CRN in Fig. 1 when the species $\ch{S}$ is
injected with flux $I_{\ch{S}}>0$ and $\ch{P}$ is extracted with flux
$I_{\ch{P}}<0$ such that the moiety concentration $L^{\text{S}}$ grows at rate
$d_{t}L^{\ch{S}}=I_{\text{S}}+I_{\text{P}}>0$.

This chemostatting procedure has three main consequences. First, the
conservation law $\ell^{\ch{E}}$ is unbroken, implying that the concentrations
$[\ch{E}]$ and $[\ch{ES}]$ are bounded. Second, as $L^{\text{S}}$ grows, (at
least one of) the concentrations $[\ch{S}]$ and $[\ch{P}]$ must grow as well.
Consequently, only a subset of species can grow in the Michaelis Menten CRN
under flux control unlike unimolecular CRNs (Sec. III.4), pseudo-unimolecular
CRNs (Sec. IV.2.3), and the autocatalytic CRN (Sec. V.2.1) under flux control.
Third, in the long time limit, when the concentrations $[\ch{S}]$ and
$[\ch{P}]$ are much larger than the concentrations $[\ch{E}]$ and $[\ch{ES}]$,
a time scale separation emerges between the evolution of the chemostatted
concentrations and the internal concentrations which allows us to coarse-grain
the internal concentrations as discussed in Sec. II.1.

Figure 5: Dynamics of the Michaelis Menten CRN under flux control. Evolution
of a) Concentrations. b) Relative difference, given by
$\lVert\bm{y}(t)-\bm{y}_{\text{cg}}(t)\rVert/\lVert\bm{y}_{\text{cg}}(t)\rVert$,
between the concentrations $\bm{y}(t)=([\ch{S}](t),[\ch{P}](t))^{\intercal}$
and the coarse-grained concentrations
$\bm{y}_{\text{cg}}(t)=([\ch{S}]_{\text{cg}}(t),[\ch{P}]_{\text{cg}}(t))^{\intercal}$
(given by Eq. (11)), and comparison between the concentrations and the coarse-
grained concentrations in the inset. Here $[\ch{E}](0)=1.3$,
$[\ch{ES}](0)=0.1$, $[\ch{S}](0)=1.9$, $[\ch{P}](0)=0.4$, $I_{\ch{S}}=1.0$,
$I_{\ch{P}}=-0.1$ and $k_{\pm 1}=k_{\pm 2}=1$.

We show a typical trajectory of the concentrations in Fig. 5a and compare the
dynamics of the actual concentrations $[\ch{S}]$ and $[\ch{P}]$ with the
coarse-grained dynamics $[\ch{S}]_{\text{cg}}$ and $[\ch{P}]_{\text{cg}}$
obtained solving Eq. (11) in Fig. 5b. We note that $[\ch{S}]$ and $[\ch{P}]$
grow linearly in time, $[\ch{ES}]$ saturates to a constant value, while
$[\ch{E}]$ decreases in time. Furthermore, the relative difference between the
vectors $\bm{y}(t)=([\ch{S}](t),[\ch{P}](t))^{\intercal}$ and
$\bm{y}_{\text{cg}}(t)=\left([\ch{S}]_{\text{cg}}(t),[\ch{P}]_{\text{cg}}(t)\right)^{\intercal}$
decreases in time.

We now use the coarse-grained dynamics discussed in Sec. II.1 to determine i)
the growth rate of $[\ch{S}]$ and $[\ch{P}]$, and ii) the evolution of
$[\ch{E}]$ and $[\ch{ES}]$ in the long time limit.

First, by using the diagrammatic method [35] as done in Refs. [36, 25], Eq.
(11) for the Michaelis Menten CRN reads

$\begin{split}d_{t}[\text{S}]_{\text{cg}}-I_{\text{S}}&=-L^{\ch{E}}(0)\left[\frac{k_{+1}k_{+2}[\text{S}]_{\text{cg}}-k_{-1}k_{-2}[\text{P}]_{\text{cg}}}{k_{-1}+k_{+2}+k_{+1}[\text{S}]_{\text{cg}}+k_{-2}[\text{P}]_{\text{cg}}}\right]\\\
&=I_{\text{P}}-d_{t}[\text{P}]_{\text{cg}}\,.\end{split}$ (99)

This, together with the ansatz ${[\ch{S}]_{\text{cg}}=v_{\ch{S}}t+c_{\ch{S}}}$
and ${[\ch{P}]_{\text{cg}}=v_{\ch{P}}t+c_{\ch{P}}}$ (with constants
$c_{\ch{S}}$ and $c_{\ch{P}}$) as well as the approximation
${d_{t}L^{\ch{S}}=I_{\ch{S}}+I_{\ch{P}}\simeq v_{\ch{S}}+v_{\ch{P}}}$ (valid
as long as time scale separation holds, i.e., $d_{t}[\ch{ES}]\simeq 0$), leads
to the quadratic equation

$v_{\ch{S}}=I_{\ch{S}}-L^{\ch{E}}(0)\frac{\left(k_{+1}k_{+2}+k_{-1}k_{-2}\right)v_{\ch{S}}-k_{-1}k_{-2}\left(I_{\ch{S}}+I_{\ch{P}}\right)}{\left(k_{+1}-k_{-2}\right)v_{\ch{S}}+k_{-2}\left(I_{\ch{S}}+I_{\ch{P}}\right)}\,,$
(100)

once the terms of order $\mathcal{O}(1/t)$ are neglected. Equation (100)
admits the following solution

$\begin{split}v_{\ch{S}}&=\frac{k_{-2}I_{\ch P}+(2k_{-2}-k_{+1})I_{\ch
S}+L^{\ch{E}}(0)\left(k_{+1}k_{+2}+k_{-1}k_{-2}\right)}{2(k_{-2}-k_{+1})}\\\
&+\sqrt{\begin{aligned} &\frac{k_{-2}\left(I_{\ch{S}}+I_{\ch{P}}\right)(I_{\ch
S}+k_{-1}L^{\ch{E}}(0))}{k_{+1}-k_{-2}}+\\\ &~{}~{}\Bigg{[}\frac{k_{-2}I_{\ch
P}+(2k_{-2}-k_{+1})I_{\ch
S}+L^{\ch{E}}(0)\left(k_{+1}k_{+2}+k_{-1}k_{-2}\right)}{2(k_{-2}-k_{+1})}\Bigg{]}^{2}\,,\\\
\end{aligned}}\end{split}$ (101)

which implies that the growth rates $v_{\ch{S}}$ and $v_{\ch{P}}$ depend on
the value of $L^{\ch{E}}(0)$. This dependence results from the fact that the
interconversion of $\ch{S}$ into $\ch{P}$ is mediated by the enzymes whose
concentrations $[\ch{E}]$ and $[\ch{ES}]$ are bounded by $L^{\ch{E}}(0)$.

Second, in the coarse-grained dynamics, the concentrations of the internal
species $[\ch{E}]$ and $[\ch{ES}]$ quickly relax to the steady state
determined by $[\ch{S}]$ and $[\ch{P}]$ which, according to the diagrammatic
method [35], read

$\displaystyle[\ch{E}]_{\text{ss}}$
$\displaystyle=\frac{(k_{-1}+k_{+2})L^{\ch{E}}(0)}{k_{-1}+k_{+2}+k_{+1}[\ch{S}]_{\text{cg}}+k_{-2}[\ch{P}]_{\text{cg}}}\,,$
(102a) $\displaystyle[\ch{ES}]_{\text{ss}}$
$\displaystyle=[\ch{E}]_{\text{ss}}\left(\frac{k_{+1}[\ch{S}]_{\text{cg}}+k_{-2}[\ch{P}]_{\text{cg}}}{k_{-1}+k_{+2}}\right)\,.$
(102b)

Equations (102a) and (102b), together with $[\ch{S}]\sim v_{\ch{S}}t$ and
$[\ch{P}]\sim v_{\ch{P}}t$, lead to the scaling $[\ch{E}]=\mathcal{O}(t^{-1})$
and $[\ch{ES}]=\mathcal{O}(1)$ observed in Figs. 5a.

We now turn to the corresponding thermodynamics. We start by considering the
EPR (18), $d_{t}G$ (17) and chemical work rate (19), whose typical time
evolution is plotted in Fig. 6a. On the one hand, the EPR relaxes to a nonzero
constant value and, therefore, the Michaelis Menten CRN under flux control
grows out of equilibrium like pseudo-unimolecular CRNs under flux control (see
Sec. IV.2.4). On the other hand, the chemical work rate and $d_{t}G$ increase
in time approaching each other and, consequently, the efficiency of growth
(53) approaches one unlike in pseudo-unimolecular CRNs (see Sec. IV.2.4). We
then examine the nonconservative work rate (22) and the moiety work rate (20)
whose typical evolution is shown in Fig. 6b (when $\ch{S}$ is chosen as the
potential species). At long times, the nonconservative work saturates to a
constant value while the moiety work rate increases with time and scales like
$d_{t}G$.

Figure 6: Thermodynamics of the Michaelis Menten CRN under flux control during
the dynamics in Fig. 5. Evolution of a) EPR, $d_{t}G$, chemical work rate, and
efficiency of growth (inset) b) EPR, $d_{t}G$, nonconservative work rate, and
moiety work rate. Here, $\mu^{0}_{\ch{E}}=1$, $\mu^{0}_{\ch{ES}}=2$,
$\mu^{0}_{\ch{S}}=1$, $\mu^{0}_{\ch{P}}=1$.

This specific behavior of the thermodynamic quantities can be explained using
the coarse-grained dynamics described in Sec. II.1. By using the ansatz
$[\ch{S}]=v_{\ch{S}}t+c_{\ch{S}}$ and $[\ch{P}]=v_{\ch{P}}t+c_{\ch{P}}$ with
Eq. (14) and neglecting terms of order $\mathcal{O}(1/t)$, we first recognize
that in the long time limit the chemical potentials of S and P read

${\mu}_{\alpha}\approx
RT\text{ln}(t){\ell}^{\ch{S}}_{\alpha}+RT\text{ln}(v_{\alpha})+\mu^{0}_{\alpha}=\mathcal{O}(\text{ln}(t))\,.$
(103)

Then, the EPR can be expressed, according to Ref. [25], as

$T\dot{\Sigma}=-\left(\mu_{\text{P}}-\mu_{\text{S}}\right)\bar{j}_{1}=\mathcal{O}(1)\,,$
(104)

where, using Eq (11), the ansatz $[\ch{S}]=v_{\ch{S}}t+c_{\ch{S}}$ and
$[\ch{P}]=v_{\ch{P}}t+c_{\ch{P}}$, and Eq. (103),

$\bar{j}_{1}=I_{\ch{S}}-v_{\ch{S}}=\mathcal{O}(1)\,,$ (105)

and

$\mu_{\text{P}}-\mu_{\text{S}}=(\mu^{0}_{\ch{P}}-\mu^{0}_{\ch{S}})+RT\text{ln}\left(\frac{v_{\text{P}}}{v_{\text{S}}}\right)=\mathcal{O}(1)\,.$
(106)

Namely, we recover the constant (and positive) entropy production rate in
Figs. 6a and 6b. By using Eq. (103), $d_{t}G$ (17) and the chemical work rate
(19) become

$\displaystyle\begin{split}\dot{w}_{\text{c}}&=RT\text{ln}(t)\left(d_{t}L^{\ch{S}}\right)+RT\sum_{\alpha\in
Y}\left(\mu^{0}_{\alpha}+RT\text{ln}(v_{\alpha})\right){\tilde{I}_{\alpha}}\\\
&=\mathcal{O}(\text{ln}(t))\,,\end{split}$ (107)
$\displaystyle\begin{split}d_{t}G&=\text{RT}\text{ln}(t)\left(d_{t}L^{\ch{S}}\right)+RT\sum_{\alpha\in
Y}\left(\mu^{0}_{\alpha}+RT\text{ln}(v_{\alpha})\right)v_{\alpha}\\\
&=\mathcal{O}(\text{ln}(t))\,,\end{split}$ (108)

namely, they grow logarithmically in time in agreement with Fig. 6a.
Crucially, the EPR in Eq. (104) is subleading with respect to the chemical
work in Eq. (107) and, consequently, the efficiency of growth (53) approaches
one even if the CRN undergoes nonequilibrium growth.

We conclude by examining the moiety work rate and the nonconservative work
rate. Using Eq. (106) in Eqs.(20) and (22), we obtain

$\displaystyle\dot{w}_{\text{m}}$
$\displaystyle=\text{RT}\text{ln}(t)\left(d_{t}L^{\ch{S}}\right)+RT\text{ln}\left(v_{\text{S}}\exp\left(\frac{{\mu}^{0}_{\text{S}}}{RT}\right)\right)d_{t}L^{\ch{S}}$
(109) $\displaystyle=\mathcal{O}(\text{ln}(t))\,,$
$\displaystyle\dot{w}_{\text{nc}}$
$\displaystyle=(\mu_{\text{P}}-\mu_{\text{S}})I_{\text{P}}=\mathcal{O}(1)\,,$
(110)

implying

$\begin{split}d_{t}{G}&\sim\dot{w}_{\text{m}}=\mathcal{O}(\text{ln}(t))\,,\\\
\dot{w}_{\text{nc}}&\sim T\dot{\Sigma}=\mathcal{O}(1)\,,\\\ \end{split}$ (111)

in the long time limit. Note that the values predicted in Eqs. (109) and (104)
are shown in Fig. 6b and are consistent with the numerical simulations.

####  V.2.3 Cyclic Michaelis Menten CRN

We consider the Cyclic Michaelis Menten CRN in Fig. 1 when the species
$\ch{S}$ is injected with flux $I_{\ch{S}}>0$ and $\ch{P}$ is extracted with
flux $I_{\ch{P}}<0$ such that the moiety concentration $L^{\text{S}}$ grows at
rate $d_{t}L^{\ch{S}}=I_{\text{S}}+I_{\text{P}}>0$.

From a dynamic point of view, the Cyclic Michaelis Menten CRN and the
Michaelis Menten CRN (Sec. V.2.2) have the same qualitative behavior. The
concentrations $[\ch{E}]$ and $[\ch{ES}]$ are bounded while $[\ch{S}]$ and
$[\ch{P}]$ grow linearly in time (see their typical evolution in Fig. 7a). At
long times, because of a time scale separation between the chemostatted and
internal species, the evolution of
$\bm{y}(t)=([\ch{S}](t),[\ch{P}](t))^{\intercal}$ converges towards
$\bm{y}_{\text{cg}}(t)=\left([\ch{S}]_{\text{cg}}(t),[\ch{P}]_{\text{cg}}(t)\right)^{\intercal}$
following Eq. (11) (see the typical evolution in Fig. 7b)

Figure 7: Dynamics of the Cyclic Michaelis Menten CRN under flux control.
Evolution of a) Concentrations b) Relative difference
$\lVert\bm{y}(t)-\bm{y}_{\text{cg}}(t)\rVert/\lVert\bm{y}_{\text{cg}}(t)\rVert$,
between the concentrations $\bm{y}(t)=([\ch{S}](t),[\ch{P}](t))^{\intercal}$
and the coarse grained concentrations
$\bm{y}_{\text{cg}}(t)=([\ch{S}]_{\text{cg}}(t),[\ch{P}]_{\text{cg}}(t))^{\intercal}$
(given by Eq. (11)) and comparison between the concentrations and the coarse-
grained concentrations in the inset. Here $[\ch{E}](0)=2.0$,
$[\ch{ES}](0)=1.4$, $[\ch{S}](0)=3.1$, $[\ch{P}](0)=5.2$, $I_{\ch{S}}=0.45$,
$I_{\ch{P}}=-0.10$ and $k_{\pm 1}=1$, $k_{+2}=k_{+3}=e^{0.25}$,
$k_{-2}=k_{-3}=e^{-0.25}$.

We use the coarse-grained dynamics to understand the long-time dynamics of the
concentrations, as done in Sec. V.2.2. Using the diagrammatic method [35], Eq.
(11) for the Cyclic Michaelis Menten CRN reads

$\begin{split}d_{t}[\text{S}]_{\text{cg}}&+k_{+3}[\text{S}]_{\text{cg}}-k_{-3}[\text{P}]_{\text{cg}}-I_{\text{S}}\\\
&=-L^{\ch{E}}(0)\left[\frac{k_{+1}k_{+2}[\text{S}]_{\text{cg}}-k_{-1}k_{-2}[\text{P}]_{\text{cg}}}{k_{-1}+k_{+2}+k_{+1}[\text{S}]_{\text{cg}}+k_{-2}[\text{P}]_{\text{cg}}}\right]\\\
&=I_{\text{P}}+k_{+3}[\text{S}]_{\text{cg}}-k_{-3}[\text{P}]_{\text{cg}}-d_{t}[\text{P}]_{\text{cg}}\,,\end{split}$
(112)

which, together with the ansatz $[\ch{S}]_{\text{cg}}=v_{\ch{S}}t+c_{\ch{S}}$
and $[\ch{P}]_{\text{cg}}=v_{\ch{P}}t+c_{\ch{P}}$ (with $c_{\ch{S}}$ and
$c_{\ch{P}}$ some constant parameters), the approximation
$d_{t}L^{\ch{S}}=I_{\ch{S}}+I_{\ch{P}}\simeq v_{\ch{S}}+v_{\ch{P}}$ (valid as
long as time scale separation holds, i.e., $d_{t}[\ch{ES}]\simeq 0$), leads to

$\begin{split}v_{\ch{S}}=&\left(k_{-3}\left(I_{\ch{S}}+I_{\ch{P}}\right)-\left(k_{+3}+k_{-3}\right)v_{\ch{S}}\right)t\\\
&+I_{\ch{S}}+{\left(k_{-3}c_{\ch{P}}-k_{+3}c_{\ch{S}}\right)}\\\
&-L^{\ch{E}}(0)\frac{\left(k_{+1}k_{+2}+k_{-1}k_{-2}\right)v_{\ch{S}}-k_{-1}k_{-2}\left(I_{\ch{S}}+I_{\ch{P}}\right)}{\left(k_{+1}-k_{-2}\right)v_{\ch{S}}+k_{-2}\left(I_{\ch{S}}+I_{\ch{P}}\right)}\,.\end{split}$
(113)

For consistency, the term proportional to time on the right-hand side of Eq.
(113) must vanish. Hence,

$\displaystyle v_{\text{S}}$
$\displaystyle=\frac{k_{-3}}{k_{+3}+k_{-3}}\left(I_{\text{S}}+I_{\text{P}}\right)\,,$
(114) $\displaystyle v_{\text{P}}$
$\displaystyle=\frac{k_{+3}}{k_{+3}+k_{-3}}\left(I_{\text{S}}+I_{\text{P}}\right)\,.$
(115)

We emphasize that unlike the Michaelis Menten CRN (see Eq. (101)), the growth
rates $v_{\ch{S}}$ and $v_{\ch{P}}$ are independent of $L^{\ch{E}}(0)$ because
of the reaction directly interconverting S into P. Finally, according to the
diagrammatic method [35], $[\ch{E}]_{\text{ss}}$ and $[\ch{ES}]_{\text{ss}}$
still follow Eq. (102) and, therefore, $[\ch{E}]=\mathcal{O}(t^{-1})$ and
$[\ch{ES}]=\mathcal{O}(1)$ (as observed in Figs. 7a).

The main difference between the Cyclic Michaelis Menten CRN and the Michaelis
Menten CRN (see Sec. V.2.2) under flux control concerns the thermodynamics of
growth. We start by considering the EPR (18), $d_{t}G$ (17) and chemical work
rate (19), whose typical time evolution is plotted in Fig. 8a. Like the
Michaelis Menten CRN (see Fig. 6a), the $d_{t}G$ and the chemical work
increase in time and coincide at long times, while the efficiency of growth
(53) goes to one. Unlike the Michaelis Menten CRN, EPR continuously decreases
to zero. This means that the Cyclic Michaelis Menten CRN under flux control
undergoes equilibrium growth like unimolecular CRNs (Sec. III.4) and the
autocatalytic CRN (Sec. V.2.1). We then examine the nonconservative work rate
(22) and the moiety work rate (20) whose typical evolution is shown in Fig. 8b
(when $\ch{S}$ is chosen as the potential species). Like the Michaelis Menten
CRN (see Fig. 6b), the moiety work rate increases with time and scales like
$d_{t}G$. Unlike the Michaelis Menten CRN, the nonconservative work decays
like the EPR.

Figure 8: Thermodynamics of the Cyclic Michaelis Menten CRN under flux control
during the dynamics in Fig. 7. Evolution of a) EPR, $d_{t}G$, chemical work
rate, and efficiency of growth (inset) b) EPR, $d_{t}G$, nonconservative work
rate, and moiety work rate. Here, $\mu^{0}_{\ch{E}}=0.5$,
$\mu^{0}_{\ch{ES}}=2$, $\mu^{0}_{\ch{S}}=1.5$, $\mu^{0}_{\ch{P}}=1$.

This specific behavior of the thermodynamic quantities can be explained as a
consequence of the timescale separation by using the coarse-grained dynamics
discussed in Sec II.1. To do so, we start by recognizing that in the long time
limit, the chemical potentials of the chemostatted species read

$\displaystyle{\mu}_{\alpha}$
$\displaystyle=RT\text{ln}(t){\ell}^{\ch{S}}_{\alpha}+RT\text{ln}\left({v}_{\alpha}\exp\left(\frac{{\mu}^{0}_{\alpha}}{RT}\right)\right)+RT\text{ln}\left(1+\frac{c_{\alpha}}{v_{\alpha}t}\right)$
$\displaystyle=\mathcal{O}(\text{ln}(t))\,,$ (116)

where we used the ansatz $[\ch{S}]=v_{\ch{S}}t+c_{\ch{S}}$ and
$[\ch{P}]=v_{\ch{P}}t+c_{\ch{P}}$ in Eq. (14). Note that in Eq. (116) we did
not neglect the term of order $\mathcal{O}(1/t)$ like in Eq. (103) because it
explains the decrease of the EPR seen in Figs. 8a and 8b. Indeed, according to
Ref. [25], the EPR can be expressed as

$T\dot{\Sigma}=-\left(\mu_{\text{P}}-\mu_{\text{S}}\right)(\bar{j}_{1}+\bar{j}_{3})=\mathcal{O}(1/t),$
(117)

where, by using Eq (11) and the ansatz $[\ch{S}]=v_{\ch{S}}t+c_{\ch{S}}$ and
$[\ch{P}]=v_{\ch{P}}t+c_{\ch{P}}$,

$\bar{j}_{1}+\bar{j}_{3}=I_{\ch{S}}-v_{\ch{S}}=\mathcal{O}(1)\,,$ (118)

while, by using Eqs. (114), (115), (116), the approximation
$\text{ln}(1+x)\approx 1+x$ for $x\ll 1$, together with the local detailed
balance (15) (implying
$v_{\text{S}}\exp(\mu^{0}_{\text{S}}/RT)=v_{\text{P}}\exp(\mu^{0}_{\text{P}}/RT)$),

$\mu_{\ch{P}}-\mu_{\ch{S}}={RT\frac{1}{t}\left(\frac{c_{\ch{P}}}{v_{\text{P}}}-\frac{c_{\ch{S}}}{v_{\text{S}}}\right)}=\mathcal{O}(1/t)\,.$
(119)

We now turn to the other thermodynamic quantities. By using Eq. (116) and
neglecting terms of order $\mathcal{O}(1/t)$, the $d_{t}G$ (17) and the
chemical work rate (19) become

$\displaystyle\begin{split}\dot{w}_{\text{c}}&=RT\text{ln}(t)\left(d_{t}L^{\ch{S}}\right)+RT\sum_{\alpha\in
Y}\left(\mu^{0}_{\alpha}+RT\text{ln}(v_{\alpha})\right){\tilde{I}_{\alpha}}\\\
&=\mathcal{O}(\text{ln}(t))\,,\end{split}$ (120)
$\displaystyle\begin{split}d_{t}G&=\text{RT}\text{ln}(t)\left(d_{t}L^{\ch{S}}\right)+RT\sum_{\alpha\in
Y}\left(\mu^{0}_{\alpha}+RT\text{ln}(v_{\alpha})\right)v_{\alpha}\\\
&=\mathcal{O}(\text{ln}(t))\,,\end{split}$ (121)

namely, they grow logarithmically in time, $d_{t}G\approx\dot{w}_{c}$, and
$\eta\to 1$ in agreement with agreement with Fig. 8b. Furthermore, by using
Eq. (116), the moiety work rate (20) and the nonconservative work rate (22)
read

$\displaystyle\begin{split}\dot{w}_{\text{m}}&=\text{RT}\text{ln}(t)\left(d_{t}L^{\ch{S}}\right){}+RT\text{ln}\left(v_{\text{S}}\exp\left(\frac{{\mu}^{0}_{\text{S}}}{RT}\right)\right)d_{t}L^{\ch{S}}\\\
&=\mathcal{O}(\text{ln}(t))\,,\end{split}$ (122)
$\displaystyle\dot{w}_{\text{nc}}$ $\displaystyle\approx
RT\frac{1}{t}\left(\frac{c_{\ch{P}}}{v_{\text{P}}}-\frac{c_{\ch{S}}}{v_{\text{S}}}\right)I_{\text{P}}=\mathcal{O}(t^{-1})\,.$
(123)

implying

$\begin{split}d_{t}{G}&\sim\dot{w}_{\text{m}}=\mathcal{O}(\text{ln}(t))\,,\\\
\dot{w}_{\text{nc}}&\sim T\dot{\Sigma}=\mathcal{O}(t^{-1})\,,\\\ \end{split}$
(124)

in the long time limit.

Note that direct interconversion of S into P explains why the Cyclic Michaelis
Menten CRN grows close to equilibrium unlike the Michaelis Menten CRN. On the
one hand, the rate of interconversion of S into P (via the enzymatic
mechanism) in the Michaelis Menten CRN is limited by the bounded concentration
$[\ch{E}]$ and $[\ch{ES}]$. This prevents the equilibration between the
chemostatted species, i.e., it establishes a difference between the chemical
potentials of S and P and creates the nonconservative force (106). On the
other hand, the rate of the direct interconversion of S into P in Cyclic
Michaelis Menten CRN increases with the concentrations $[\ch{S}]$ and
$[\ch{P}]$. Their chemical potentials can thus rapidly equilibrate and the
nonconservative force (119) decreases in time.

####  V.2.4 Minimal Metabolic CRN

The minimal metabolic CRN has already been numerically studied in Ref. [16]
when the species S and F are injected at constant rates $I_{\ch{S}}>0$ and
$I_{\ch{F}}>0$ such that that the moiety concentrations $L^{\ch{S}}$ and
$L^{\ch{F}}$ grow with rates $d_{t}L^{\ch{S}}=I_{\text{S}}>0$ and
$d_{t}L^{\ch{F}}=I_{\text{F}}>0$. Correspondingly, the concentrations
$[\ch{S}]$, $[\ch{P}]$, $[\ch{F}]$ and $[\ch{W}]$ grow and the EPR decays in
time, implying that the CRN grows close to equilibrium. This was explained by
noting that $\dot{w}_{\text{nc}}=0$ at all times since the chemostatting
procedure does not create nonconservative forces (see Eq. (22)).

Here, we study the minimal metabolic CRN when the species S and F are injected
at constant rates $I_{\ch{S}}>0$ and $I_{\ch{F}}>0$ while the species P and W
are extracted at constant rates $I_{\ch{P}}<0$ and $I_{\ch{W}}<0$ in such a
way that the moiety concentrations $L^{\text{S}}$ and $L^{\text{F}}$ grow with
rates $d_{t}L^{\ch{S}}=I_{\text{S}}+I_{\text{P}}>0$ and
$d_{t}L^{\ch{F}}=I_{\text{F}}+I_{\text{W}}>0$. This chemostatting procedure
has four main consequences. First, it creates two nonconservative forces (see
Eq. (22)). Second, the conservation law $\ell^{E}$ is unbroken and, therefore,
the concentrations $[\ch{E}]$, $[\ch{ES}]$, $[\ch{EW}]$ and $[\ch{E}^{*}]$ are
bounded. Third, as $L^{\text{S}}$ (resp. $L^{\ch{F}}$) grows, (at least one
of) the concentrations $[\ch{S}]$ and $[\ch{P}]$ (resp. $[\ch{F}]$ and
$[\ch{W}]$) must grow as well. Finally, in the long time limit, when the
concentrations of the chemostatted species are much larger than the
concentrations of the internal species, a time scale separation emerges
between the evolution of the former and the latter implying that the dynamics
of the chemostatted species follow the coarse-grained dynamics discussed in
Sec. II.1

We now show a typical evolution of the concentrations of the chemostatted
species in Fig. 9a and compare it with the coarse-grained dynamics obtained
from Eq. (11) in Fig. 9b. We note that concentrations $[\ch{F}]$ and
$[\ch{S}]$ grow linearly in time while the concentrations $[\ch{P}]$ and
$[\ch{W}]$ saturate to constant values. Furthermore, the relative distance
between the vectors
$\bm{y}(t)=([\ch{F}](t),[\ch{W}](t),[\ch{S}](t),[\ch{P}](t)[\ch{P}](t))^{\intercal}$
and
$\bm{y}_{\text{cg}}(t)=\left([\ch{F}]_{\text{cg}}(t),[\ch{W}]_{\text{cg}}(t),[\ch{S}]_{\text{cg}}(t),[\ch{P}]_{\text{cg}}(t)\right)^{\intercal}$
increases initially and then decreases at long times.

Figure 9:  Dynamics of the minimal metabolic CRN under flux control. Evolution
of a) Concentrations of the chemostatted species $[\ch{F}](t)$, $[\ch{W}](t)$,
$[\ch{S}](t)$ and $[\ch{P}](t)$ b) Relative difference between the vectors
$\bm{y}(t)=([\ch{F}](t),[\ch{W}](t)[\ch{S}](t),[\ch{P}](t))^{\intercal}$ and
$\bm{y}_{\text{cg}}(t)=([\ch{F}]_{\text{cg}}(t),[\ch{W}]_{\text{cg}}(t),[\ch{S}]_{\text{cg}}(t),[\ch{P}]_{\text{cg}}(t))^{\intercal}$
(given by Eq. (11)), and comparison between the concentrations $[\ch{F}](t)$
and $[\ch{P}](t)$ with the coarse-grained concentrations
$[\ch{F}]_{\text{cg}}(t)$ and $[\ch{P}]_{\text{cg}}(t)$ in the inset. Here
$[\ch{E}](0)=1.5$, $[\ch{EF}](0)=[\ch{EW}](0)=[\ch{E}^{*}](0)=0.1$,
$[\ch{F}](0)=[\ch{W}](0)=0.5$, $[\ch{S}](0)=[\ch{P}](0)=1$, $I_{\ch{F}}=2.8$,
$I_{\ch{W}}=-0.1$ , $I_{\ch{S}}=2.0$, $I_{\ch{P}}=-0.05$ and $k_{\pm
1}=k_{-2}=k_{\pm 3}=k_{\pm 4}=k_{\pm 5}=2$ while $k_{+2}=2e$.

We now turn to thermodynamics. The typical time evolution of the EPR (18),
$d_{t}G$ (17) and chemical work rate (19) is plotted in Fig. 10a. The chemical
work rate and $d_{t}G$ increase in time. However, they do not coincide at long
times, unlike in the Autocatalytic CRN in Sec. V.2.1, the Michaelis Menten CRN
in Sec. V.2.2 and the Cyclic Michaelis Menten CRN in Sec. V.2.3, since the EPR
increases in time too. Consequently, the minimal metabolic CRN under flux
control grows out of equilibrium, like the Michaelis Menten CRN (see Sec.
V.2.2), and the efficiency of growth (53), which is initially close to one,
decreases at long times, unlike the Michaelis Menten CRN (see Fig. 6a).

We then examine the nonconservative work rate (22) and the moiety work rate
(20) whose typical evolution is shown in Fig. 10b (when $\ch{S}$ and $\ch{F}$
are chosen as the potential species). At long times, there is a splitting
between the magnitudes of the EPR and the nonconservative work rate on the one
hand and the moiety work rate and $d_{t}G$ on the other hand. Furthermore,

$\begin{split}d_{t}{G}&\sim\dot{w}_{\text{m}}\,,\\\ \dot{w}_{\text{nc}}&\sim
T\dot{\Sigma}>0\,,\\\ \end{split}$ (125)

similar to the Autocatalytic CRN (98), the Michaelis Menten CRN (111), and the
Cyclic Michaelis Menten CRN , (124).

Figure 10:  Thermodynamics of the minimal metabolic CRN under flux control
during the dynamics in Fig. 9. Evolution of a) EPR, $d_{t}G$, chemical work
rate, and efficiency of growth (inset) b) EPR, $d_{t}G$, nonconservative work
rate, and moiety work rate. Here, $\mu^{0}_{\ch{E}}=1$, $\mu^{0}_{\ch{EF}}=3$,
$\mu^{0}_{\ch{EW}}=2$, $\mu^{0}_{\ch{E}^{*}}=4$, $\mu^{0}_{\ch{F}}=2$,
$\mu^{0}_{\ch{W}}=1$, $\mu^{0}_{\ch{S}}=1$, $\mu^{0}_{\ch{P}}=2$.

####  V.2.5 Summary

From a dynamical standpoint, growth in multimolecular CRNs under flux control
can significantly differ from unimolecular and pseudo-unimolecular CRNs. For
instance, i) the concentrations of some species scale as
$\mathcal{O}(\sqrt{t})$ in the Autocatalytic CRN (see Sec. V.2.1), while the
concentrations of all species scale as $\mathcal{O}({t})$ in unimolecular and
pseudo-unimolecular CRNs; ii) the concentrations of only a subset of species
grow in the Michaelis Menten CRN (see Sec. V.2.2).

From a thermodynamic standpoint, some multimolecular CRNs grow close to
equilibrium (e.g., the Autocatalytic CRN and the cyclic Michaelis Menten CRN),
while other multimolecular CRNs grow out ot equilibrium (e.g., the Michaelis
Menten CRN and the minimal metabolic CRN). We noticed that the presence of an
unbroken conservation law is a necessary but not sufficient condition for
nonequilibrium growth. On the one hand, it can prevent the equilibration of
the chemostatted species and create nonconservative forces like in the
Michaelis Menten CRN and the minimal metabolic CRN. On the other hand, if
there are other reaction pathways leading to equilibration, like the direct
interconversion of S into P in the cyclic Michaelis Menten CRN, the CRN grow
at equilibrium (see discussion at the end of Subs. V.2.3).

###  V.3 Growth under Mixed Control

Unimolecular and pseudo-unimolecular CRNs cannot grow under mixed control, as
we proved in Secs. III.3 and Sec. IV.2.1, respectively. On the other hand,
multimolecular CRNs can grow under mixed control as proven in Ref. [16] for
the minimal metabolic CRN.

Here, we first prove in general that the concentrations of extracted
chemostatted species are bounded for CRNs growing under mixed control with a
monotonically increasing mass density, i.e., $d_{t}L^{m}>0$. Indeed, from Eq.
(24) and the definition of mixed control, the mass density $L^{m}$ follows

$d_{t}L^{m}=-\ell^{m}_{Y}\cdot\tilde{\mathbb{D}}\bm{y}+\ell^{m}_{Y}\cdot\tilde{\bm{I}}\,,$
(126)

which, together with $d_{t}L^{m}>0$, implies

$[\alpha]\leq\frac{\ell^{m}_{Y}\cdot\tilde{\bm{I}}}{\ell^{m}_{\alpha}\tilde{\mathbb{D}}_{\alpha,\alpha}}\,,$
(127)

for every chemostatted species $\alpha$ such that
$\tilde{\mathbb{D}}_{\alpha,\alpha}>0$, i.e., every extracted chemostatted
species. Note that this has the obvious implication that Continuous-flow
Stirred Tank Reactors (CSTR) cannot grow as already proven in Sec. II.4.

We now examine the specific dynamics and thermodynamics of the multimolecular
CRNs listed in Fig. 1 under mixed control.

####  V.3.1 Autocatalytic CRN

##### Case I.

We start by considering the autocatalytic CRN when S (resp. P) is injected
with constant rate $I_{\text{S}}>0$ (resp. $I_{\text{P}}>0$) and extracted
with rate ${k}^{e}_{\ch{S}}[\text{S}]$ (resp. ${k}^{e}_{\ch{P}}[\text{P}]$).
We first analytically prove that the autocatalytic CRN cannot grow with a
monotonically increasing mass density, and we then use numerical simulations
to show that the mass density does not grow.

Figure 11: Evolution of the mass density of the Autocatalytic CRN under mixed
control where the species S and P have influx rates $I_{\text{S}}=1.0$ and
$I_{\text{P}}=1.5$ and outflux rates $k^{e}_{\ch{S}}[\text{S}]$ and
$k^{e}_{\ch{P}}[\text{P}]$ for thirty randomly chosen initial conditions and
a) ${I}/{\tilde{k}}=750$, b) ${I}/{\tilde{k}}=75$, c) ${I}/{\tilde{k}}=15$, d)
${I}/{\tilde{k}}=0.75$, where the parameter
${I}/{\tilde{k}}\equiv\text{max}\\{I_{\text{S}},I_{\text{P}}\\}/\text{min}\\{k^{e}_{\ch{S}},k^{e}_{\ch{P}}\\}$.
Here, $k_{\pm 1}=k_{\pm 2}=k_{\pm 3}=1$ and $k^{e}_{\ch{P}}=2.0$.

Let us assume that the autocatalytic CRN grows with $d_{t}L^{m}>0$ which,
according to Eq. (127), implies

$[\text{S}],[\text{P}]\leq\frac{I_{\text{S}}+I_{\text{P}}}{\text{min}({k}^{e}_{\ch
S},{k}^{e}_{\ch P})}\,.$ (128)

Hence, at least one of the concentrations $[\text{E}]$ and $[\text{E}^{*}]$
grow. However, according to Eq. (2) and the definition of mixed control, the
evolution equation of $[\text{S}]$,

$d_{t}[\text{S}]+(k_{+1}+{k}^{e}_{\ch{S}})[\text{S}]=k_{-1}[\text{E}]+I_{\text{S}}\,,$
(129)

has the formal solution

$\displaystyle[\text{S}](t)=$
$\displaystyle[\text{S}](0)e^{-(k_{+1}+k^{e}_{\ch{S}})t}+\frac{I_{\text{S}}}{k_{+1}+k^{e}_{\ch{S}}}\left(1-e^{-(k_{+1}+k^{e}_{\ch{S}})t}\right)$
$\displaystyle+k_{-1}\int_{0}^{t}ds[\text{E}](s)e^{-(k_{+1}+k^{e}_{\ch{S}})(t-s)}\,,$
(130)

implying that $[\text{E}](t)$ must also be bounded. By applying exactly the
same reasoning to the evolution equation of $[\ch{P}]$, we find that
$[\ch{E}^{*}]$ must be bounded too. Thus, the autocatalytic CRN does not grow
with $d_{t}L^{m}>0$.

We now turn to numerical simulations. Figure 11 shows that the mass density
does not grow for different initial conditions and different values of the
parameter
${I}/{\tilde{k}}\equiv{\text{max}\big{(}I_{\text{S}},I_{\text{P}}\big{)}}/\text{min}\big{(}{k}^{e}_{\ch{S}},k^{e}_{\ch{P}}\big{)}$,
representing the ratio between influx and outflux. Note that the autocatalytic
CRN does not grow even for large ${I}/{\tilde{k}}$, namely, large influx
rates.

##### Case II.

We now consider the autocatalytic CRN when S is being injected with rate
$I_{\text{S}}>0$ while $\ch{P}$ is extracted with rate
${k^{e}_{\ch{P}}}[\text{P}]$. As for the previous case, we first analytically
prove that the autocatalytic CRN cannot grow with a monotonically increasing
mass density and then we use numerical simulations to show that the mass
density does not grow.

Let us assume that $d_{t}L^{m}>0$ which, according to Eq. (127), implies

$[\text{P}]\leq\frac{I_{\text{S}}}{k^{e}_{\ch P}}\,.$ (131)

Figure 12:  Evolution of the mass density of autocatalytic CRN under mixed
control when the species S has influx rate $I_{\text{S}}>0$ while $\ch{P}$ is
extracted with a rate $k^{e}_{\ch{P}}[\text{P}]$ for thirty randomly chosen
initial conditions is plotted and a) ${I}/{\tilde{k}}=500$, b)
${I}/{\tilde{k}}=50$, c) ${I}/{\tilde{k}}=10$, d) ${I}/{\tilde{k}}=0.5$, where
${I}/{\tilde{k}}=I_{\text{S}}/k^{e}_{\ch P}$. Here, $k_{\pm 1}=k_{\pm
2}=k_{\pm 3}=1$ and $I_{\text{S}}=1.0$.

Then, the rate equation (2) for P under mixed control

$d_{t}[\text{P}]+(k_{+2}+k^{e}_{\ch P})[\text{P}]=k_{-2}[\text{E}^{*}]\,,$
(132)

has the formal solution,

$\begin{split}[\text{P}](t)=&[\text{P}](0)e^{-(k_{+2}+k^{e}_{\ch P})t}\\\
&+k_{-2}\int_{0}^{t}ds[\text{E}^{*}](s)e^{-(k_{+2}+k^{e}_{\ch
P})(t-s)}\,,\end{split}$ (133)

which implies that $[\text{E}^{*}]$ is bounded (according to the exact same
reasoning as below Eq. (130)). Hence, only $[\ch{S}]$ and $[\ch{E}]$ can grow.
Let us consider the sum $\mathbb{A}=[\ch{S}]+[\ch{E}]$ which is unbounded if
the CRN grows. From Eq. (2) and the definition of mixed control, $\mathbb{A}$
follows

$d_{t}\mathbb{A}=I_{\text{S}}+2k_{+3}[\ch{E}][\ch{E}^{*}]-2k_{-3}[\ch{E}]^{3}\,.$
(134)

By using now that $[\text{E}^{*}]$ is bounded and assuming that
$d_{t}\mathbb{A}>0$, $[\ch{E}](t)$ satisfies
$2k_{-3}[\ch{E}]^{3}<I_{\text{S}}+2k_{+3}[\ch{E}][\ch{E}^{*}]<I_{\text{S}}+2C[\ch{E}]$
(with $C$ being a positive number such that $C>k_{+3}[\ch{E}^{*}]$), or,
equivalently, $[\ch{E}]$ is bounded. Finally, rate equation (2) for $\ch{S}$
satisfies

$d_{t}[\text{S}]+k_{+1}[\text{S}]=k_{-1}[\text{E}]+I_{\text{S}}<C^{\prime}\,,$
(135)

where $C^{\prime}$ is an appropriately chosen positive number, with the formal
solution

$[\text{S}](t)\leq[\text{S}](0)e^{-k_{+1}t}+\frac{C^{\prime}}{k_{+1}}\left(1-e^{-k_{+1}t}\right)<\infty\,,$
(136)

implying that $[\ch{S}]$ is also bounded. Thus, the autocatalytic CRN does not
grow with $d_{t}L^{m}>0$ and $d_{t}\mathbb{A}>0$.

We now turn to numerical simulations. Figure 12 shows that the mass density
does not grow for different initial conditions and different values of the
parameter ${I}/{\tilde{k}}=I_{\text{S}}/k^{e}_{\ch P}$. Note that also in this
case the autocatalytic CRN does not grow even for large ${I}/{\tilde{k}}$,
namely, large influx rates.

####  V.3.2 Michaelis Menten

We consider the Michaelis Menten CRN when the substrate S is injected at the
rate $I_{\ch{S}}>0$ and the product P is extracted at the rate $k^{e}_{\ch
P}[\text{P}]$.

Note that, as in Sec. V.2.2, this chemostatting procedure does not break the
conservation law $\ell^{\text{E}}$ and the concentrations $[\text{E}]$ and
$[\text{ES}]$ are thus bounded. Furthermore, according to Eq. (127), if the
CRN grows with $d_{t}L^{\ch S}>0$, the concentration $[\ch{P}]$ is bounded
above by

$[\text{P}]\leq\frac{I_{\text{S}}}{k^{e}_{\ch P}}\,.$ (137)

Hence, only the concentration $[\ch{S}]$ can grow.

Figure 13: Evolution of $d_{t}L^{\ch S}$ for the Michaelis Menten CRN under
mixed control when S is injected at the constant rate $I_{\text{S}}=1.0$ and P
is extracted at the rate $I_{\text{P}}=-k^{e}_{\ch{P}}[\text{P}]$. a) Initial
condition $[\text{E}](0)=0.3$, $[\text{ES}](0)=0.8$, $[\text{S}](0)=1.9$, and
$[\text{P}](0)=0.4$ b) Initial condition $[\text{E}](0)=0.3$,
$[\text{ES}](0)=0.1$, $[\text{S}](0)=1.9$, and $[\text{P}](0)=0.4$. Here
$k_{\pm 1}=k_{\pm 2}=1$.

Figure 14: Dynamics of the Michaelis Menten CRN under mixed control. Evolution
of a) Concentrations b) Relative difference,
$\lVert\bm{y}(t)-\bm{y}_{\text{cg}}(t)\rVert/\lVert\bm{y}_{\text{cg}}(t)\rVert$,
between the concentrations $\bm{y}(t)=([\ch{S}](t),[\ch{P}](t))^{\intercal}$
and the coarse grained concentrations
$\bm{y}_{\text{cg}}(t)=([\ch{S}]_{\text{cg}}(t),[\ch{P}]_{\text{cg}}(t))^{\intercal}$
(given by Eq. (11)) and comparison between the concentrations and the coarse-
grained concentrations in the inset. Here, $[\text{E}](0)=0.3$,
$[\text{ES}]{0}=0.1$, $[\text{S}](0)=1.9$, $[\text{P}](0)=0.4$,
$I_{\ch{S}}=1.0$, $k^{e}_{\ch{P}}=0.5$, $k_{\pm 1}=k_{\pm 2}=1$.

Strikingly, as shown in Fig. 13, depending on the initial conditions but
independently of the value of $k^{e}_{\ch P}$, $d_{t}L^{\ch{S}}$ can either
vanish or reach a constant positive value. In the latter case, the typical
evolution of the concentrations is shown in Fig. 14a: $[\ch{S}]$ grows
linearly in time, $[\ch{P}]$ and $[\ch{ES}]$ saturate to constant values,
while $[\ch{E}]$ decreases in time. Note that, as shown in Fig. 14b, the
evolution of the concentrations $[\ch{S}]$ and $[\ch{P}]$ is well approximated
by the coarse-grained dynamics obtained solving Eq. (11).

When $[\ch{S}]$ grows, a time scale separation emerges between the evolution
of the concentrations $[\ch{S}]$ and $[\ch{P}]$, on the one hand, and the
concentrations of the internal species, on the other hand. Thus, following
Sec. II.1, we now use the coarse-grained dynamics (11) to determine i) when
the Michaelis Menten CRN grows, ii) the growth rate of $[\ch{S}]$, and iii)
the evolution of $[\ch{P}]$, $[\ch{E}]$ and $[\ch{ES}]$ in the long time
limit. By using the diagrammatic method [35] as done in Refs. [36, 25], Eq.
(11) for the Michaelis Menten CRN under mixed control reads

$\displaystyle d_{t}[\text{S}]$
$\displaystyle=I_{\text{S}}-L^{E}(0)\left[\frac{k_{+1}k_{+2}[\text{S}]-k_{-1}k_{-2}[\text{P}]}{k_{-1}+k_{+2}+k_{+1}[\text{S}]+k_{-2}[\text{P}]}\right]\,,$
(138) $\displaystyle d_{t}[\text{P}]$
$\displaystyle=L^{\ch{E}}(0)\left[\frac{k_{+1}k_{+2}[\text{S}]-k_{-1}k_{-2}[\text{P}]}{k_{-1}+k_{+2}+k_{+1}[\text{S}]+k_{-2}[\text{P}]}\right]-k^{e}_{\ch
P}[\text{P}]\,,$ (139)

which, together with the ansatz $[\text{S}]=v_{\text{S}}t+c_{\ch S}$ (for some
constant parameters $v_{\text{S}}$ and $c_{\ch S}$) leads (in the long time
limit) to

$\displaystyle v_{\text{S}}$ $\displaystyle=I_{\text{S}}-k_{+2}L^{E}(0)\,,$
(140) $\displaystyle d_{t}[\text{P}]$ $\displaystyle=k_{+2}L^{E}(0)-k^{e}_{\ch
P}[\text{P}]\,,$ (141)

by neglecting terms of order $\mathcal{O}(1/t)$. Here, $k_{+2}L^{E}(0)$
represents the rate of interconversion of S into P in the long time limit. On
the one hand, Eq. (140) implies that the concentration $[\ch S]$ (and,
correspondingly, the moiety $L^{\ch S}$) grows, i.e., $v_{\text{S}}>0$, if
$I_{\text{S}}>k_{+2}L^{\ch{E}}(0)$ independently of $k^{e}_{\ch P}$. This
physically means that growth happens when S is injected faster than its
(enzyme-dependent) interconversion into P, which is then extracted. Note that
the numerical simulations shown in Fig. 13a (resp. Fig. 13b) satisfy
$I_{\text{S}}<k_{+2}L^{\ch{E}}(0)$ (resp. $I_{\text{S}}>k_{+2}L^{\ch{E}}(0)$).
On the other hand, Eq. (141) implies that $[\ch{P}]$ reaches the limiting
value $k_{+2}L^{E}(0)/k^{e}_{\ch P}$, consistent with Eq. (137). Finally,
according to the diagrammatic method [35], $[\ch{E}]_{\text{ss}}$ and
$[\ch{ES}]_{\text{ss}}$ still follow Eqs (102a) and, therefore, using Eqs.
(140) and (141), $[\ch{E}]=\mathcal{O}(t^{-1})$ and $[\ch{ES}]=\mathcal{O}(1)$
(as observed in Fig. 14a).

We now examine the thermodynamics of the Michaelis Menten CRN under mixed
control when $I_{\text{S}}>k_{+2}L^{\ch{E}}(0)$. We start by considering the
EPR (18), $d_{t}G$ (17) and chemical work rate (19), whose typical time
evolution is plotted in Fig. 15a. The EPR increases in time implying the
Michaelis Menten CRN under mixed control grows out of equilibrium. The
chemical work rate and $d_{t}G$ increase in time at the same rate as the EPR.
Consequently, the efficiency of growth (53) saturates to a value smaller than
one (see the inset in Fig. 15a). We then examine the nonconservative work rate
(22) and the moiety work rate (20) whose typical evolution is shown in Fig.
15b (when $\ch{S}$ is chosen as the potential species). Both the
nonconservative work rate and the moiety work rate increase with time.
Crucially, there is a splitting between the magnitudes of the EPR and the
nonconservative work rate on the one hand, and the moiety work rate and
$d_{t}G$ on the other hand.

Figure 15: Thermodynamics of the Michaelis Menten CRN under mixed control
under the dynamics in Fig. 14. Evolution of a) EPR, the chemical work rate,
and $d_{t}G$ and efficiency of growth (inset). b) EPR, moiety work rate,
nonconservative work rate, and $d_{t}G$. Here, $\mu^{0}_{\ch{E}}=1$,
$\mu^{0}_{\ch{ES}}=2$, $\mu^{0}_{\ch{S}}=1$, $\mu^{0}_{\ch{P}}=1$.

We explain this behavior of the thermodynamic quantities by using the coarse-
grained dynamics Eqs. (138) and (139) and the ansatz
$[\ch{S}]=v_{\ch{S}}t+c_{\ch{S}}$. We start by recognizing that the chemical
potential (14) of S and P in the long time limit read

$\mu_{\text{S}}=\mu^{0}_{\text{S}}+RT\text{ln}\left(v_{\text{S}}\right)+RT\text{ln}(t)=\mathcal{O}(\text{ln}(t))\,,$
(142a)
$\mu_{\text{P}}=\mu^{0}_{\ch{P}}+RT\text{ln}\left(\frac{k_{+2}L^{E}(0)}{k^{e}_{\ch
P}}\right)=\mathcal{O}(1)\,,$ (142b)

respectively, by neglecting terms of order $\mathcal{O}(1/t)$. Then, the EPR
can be expressed, according to Ref. [25], as

$T\dot{\Sigma}=\left(\mu_{\text{S}}-\mu_{\text{P}}\right)\bar{j}_{1}=\mathcal{O}(\text{ln}(t))\,,$
(143)

where, by using Eq (11), the ansatz $[\ch{S}]=v_{\ch{S}}t+c_{\ch{S}}$, and
Eqs. (142a) and (142b),

$\bar{j}_{1}=k_{+2}L^{E}(0)\,,$ (144)

and

$\begin{split}\mu_{\text{S}}-\mu_{\text{P}}&=(\mu^{0}_{\ch{S}}-\mu^{0}_{\ch{P}})+RT\text{ln}\left(\frac{k^{e}_{\ch
P}v_{\ch{S}}}{k_{+2}L^{E}(0)}\right)+RT\text{ln}(t)\\\
&=\mathcal{O}(\text{ln}(t))\,.\end{split}$ (145)

Namely, we recover the positive increase in the entropy production rate
observed in Figs. 15a and 15b. By using Eqs. (142a) and (142b), $d_{t}G$ (17)
and the chemical work rate (19) become

$\displaystyle\begin{split}\dot{w}_{c}&=RT\text{ln}(t)I_{\ch{S}}+(\mu^{0}_{\ch{S}}I_{\ch{S}}-\mu^{0}_{\ch{P}}k_{+2}L^{E}(0))+RT\text{ln}\left(v_{\text{S}}\right)I_{\ch{S}}\\\
&{}~{}~{}~{}-RT\text{ln}\left(\frac{k_{+2}L^{E}(0)}{k^{e}_{\ch
P}}\right)k_{+2}L^{E}(0)=\mathcal{O}(\text{ln}(t))\,,\end{split}$ (146)
$\displaystyle d_{t}G$
$\displaystyle=RT\text{ln}(t)v_{\ch{S}}+\mu^{0}_{\ch{S}}v_{\text{S}}+RT\text{ln}\left(v_{\text{S}}\right)v_{\ch{S}}=\mathcal{O}(\text{ln}(t))\,,$
(147)

namely, they grow logarithmically in time like the EPR in agreement with Fig.
15a. Furthermore, plugging Eqs. (147) and (146) in Eq. (53), the efficiency
saturates in the long time limit to the value

$\eta=1-\frac{k_{+2}L^{E}(0)}{I_{\text{S}}}=\frac{v_{\ch{S}}}{I_{\ch{S}}}<1\,,$
(148)

as shown in Fig. 15a.

We conclude by examining the moiety work rate and the nonconservative work
rate. By using Eq. (142a) in Eqs.(20) and (22), we obtain

$\displaystyle\dot{w}_{\text{m}}$
$\displaystyle=RT\text{ln}(t)v_{\ch{S}}+RT\text{ln}\left(v_{\text{S}}\right)v_{\ch{S}}+\mu^{0}_{\ch{S}}v_{\text{S}}=\mathcal{O}(\text{ln}(t))\,,$
(149)
$\displaystyle\begin{split}\dot{w}_{\text{nc}}&=RT\text{ln}(t)\bar{j}_{1}+(\mu^{0}_{\ch{S}}-\mu^{0}_{\ch{P}})\bar{j}_{1}\\\
&{}~{}~{}~{}~{}+RT\text{ln}\left(\frac{k^{e}_{\ch
P}v_{\ch{S}}}{k_{+2}L^{E}(0)}\right)\bar{j}_{1}=\mathcal{O}\left(\text{ln}(t)\right)\,,\end{split}$
(150)

implying

$\begin{split}d_{t}{G}&\sim\dot{w}_{\text{m}}=\mathcal{O}(\text{ln}(t))\,,\\\
\dot{w}_{\text{nc}}&\sim T\dot{\Sigma}=\mathcal{O}(\text{ln}(t))\,.\\\
\end{split}$ (151)

in the long time limit.

####  V.3.3 Cyclic Michaelis Menten

We consider the Cyclic Michaelis Menten CRN under mixed control when S is
injected at the rate $I_{\text{S}}>0$ and P is extracted at the rate
$k^{e}_{\ch P}[\ch{P}]$. As in the case of the Michaelis Menten CRN under
mixed control, Sec. V.3.2, the conservation law $\ell^{\text{E}}$ is unbroken
implying that the concentrations $[\text{E}]$ and $[\text{ES}]$ are bounded in
time. Furthermore, if the CRN grows with $d_{t}L^{\ch S}>0$, Eq. (127) becomes

$[\text{P}]\leq\frac{I_{\text{S}}}{k^{e}_{\ch P}}\,,$ (152)

namely, $[\ch{P}]$ is also bounded. Furthermore, by using Eqs. (2), (152) and
$[\ch{ES}]\leq L^{\ch{E}}(0)$, $[\ch{S}]$ satisfies

$d_{t}[\ch{S}]+k_{+3}[\ch{S}]\leq
k_{-1}L^{\ch{E}}(0)+I_{\text{S}}\left(1+\frac{k_{-3}}{k^{e}_{\ch
P}}\right)\,.$ (153)

implying that $[\ch S]$ is bounded too. We conclude that the Cyclic Michaelis
Menten CRN does not grow with $d_{t}L^{\ch S}>0$.

We now verify this result numerically. Figure 16 shows the typical evolution
of the moiety concentration $L^{\ch{S}}$ for thirty different initial
conditions and for four values of the parameter
${I}/{\tilde{k}}=I_{\text{S}}/{k^{e}_{\ch{P}}}$ (which represents the ratio of
influx to outflux rates). We see that the CRN always relaxes to a steady
state, even for large ${I}/{\tilde{k}}$, namely large influx rates.

The direct interconversion of S into P is the reason why the Cyclic Michaelis
Menten CRN does not grow under mixed control unlike the Michaelis Menten CRN
(see Sec. V.3.2). On the one hand, in the Michaelis Menten CRN the rate of
interconversion of S into P (via the enzymatic mechanism) is limited by
$k_{+2}L^{\ch E}(0)$ (see discussion below Eq. (140)). Hence, when $I_{\ch
S}>k_{+2}L^{\ch E}(0)$, $\ch{S}$ accumulates and the CRN grows. On the other
hand, the rate of the direct interconversion of S into P in the Cyclic
Michaelis Menten CRN increases with the concentration $[\ch{S}]$, implying
that S is rapidly interconverted into P which is rapidly extracted as the
concentrations $[\ch S]$ increases. Hence, $\ch{S}$ cannot accumulate and the
CRN does not grow.

Figure 16:  Evolution of $L^{\ch S}$ of the Cyclic Michaelis Menten CRN under
mixed control when the species S is injected at rate $I_{\text{S}}=5.0$ while
$\ch{P}$ is extracted at rate $I_{\text{P}}=-k^{e}_{\ch{P}}[\text{P}]$ for
thirty randomly chosen initial conditions (consistent with $L^{E}(0)=3.0$) and
a) ${I}/{\tilde{k}}=2500$, b) ${I}/{\tilde{k}}=500$, c) ${I}/{\tilde{k}}=50$,
d) ${I}/{\tilde{k}}=0.5$ where
${I}/{\tilde{k}}=I_{\text{S}}/{k^{e}_{\ch{P}}}$. Here, $k_{\pm 1}=k_{\pm
2}=k_{\pm 3}=1$.

####  V.3.4 Minimal Metabolic CRN

We consider the minimal metabolic CRN under mixed control when the species S
and F are injected at rates $I_{\ch{S}}>0$ and $I_{\ch{F}}>0$, respectively,
while the species P and W are extracted at rates ${k^{e}_{\ch P}}[\ch{P}]$ and
${k^{e}_{\ch W}}[\ch{W}]$, respectively.

Figure 17: Dynamics of the minimal metabolic CRN under mixed control when
$\ch{F}$ and $\ch{S}$ are injected at rates $I_{\ch{F}}=1.2$, $I_{\ch{S}}=2.5$
and $\ch{W}$ and $\ch{P}$ are extracted at rates $k^{e}_{\ch{W}}[\ch{W}]$ and
${k^{e}_{\ch{P}}}[\ch{P}]$ with ${k^{e}_{\ch{W}}=k^{e}_{\ch{P}}}=0.5$ .
Evolution of a) Concentrations of the chemostatted species $[\ch{F}](t)$,
$[\ch{W}](t)$, $[\ch{S}](t)$ and $[\ch{P}](t)$ b) Relative difference between
the vectors
$\bm{y}(t)=([\ch{F}](t),[\ch{W}](t)[\ch{S}](t),[\ch{P}](t))^{\intercal}$ and
$\bm{y}_{\text{cg}}(t)=([\ch{F}]_{\text{cg}}(t),[\ch{W}]_{\text{cg}}(t),[\ch{S}]_{\text{cg}}(t),[\ch{P}]_{\text{cg}}(t))^{\intercal}$
(given by Eq. (11)), and comparison between the concentrations $[\ch{F}](t)$
and $[\ch{P}](t)$ with the coarse-grained concentrations
$[\ch{F}]_{\text{cg}}(t)$ and $[\ch{P}]_{\text{cg}}(t)$ in the inset. Here,
$k_{\pm 1}=k_{\pm 2}=k_{\pm 3}=k_{\pm 4}=k_{\pm 5}=1$, $[\ch{E}](0)=0.12$,
$[\ch{ES}](0)=0.68$,$[\ch{EW}](0)=0.6$, $[\ch{E}^{*}](0)=0.04$,
$[\ch{F}](0)=0.12$, $[\ch{W}](0)=0.4$, $[\ch{S}](0)=0.36$,$[\ch{P}](0)=0.72$.

The minimal metabolic CRN under this chemostatting procedure has also been
studied numerically in Ref. [16]. It has been observed that depending on the
values of the parameters $I_{\ch{S}}$, $I_{\ch{F}}$, $k^{e}_{\ch{P}}$ and
$k^{e}_{\ch{W}}$, the CRN either relaxes to a non-equilibrium steady state or
grows. From a dynamical standpoint, only $[\ch{S}]$ and $[\ch{F}]$ grow while
$[\ch{P}]$ and $[\ch{W}]$ saturate to constant values. From a thermodynamic
standpoint, when the CRN grows, a positive increasing EPR balanced by
$\dot{w}_{\text{nc}}$ was observed.

Here, we recreate and explain the above observations when the CRN grows.
Figure 17a shows the typical evolution of the chemostatted species: the
concentrations $[\ch{F}]$ and $[\ch{S}]$ grow while both $[\ch{P}]$ and
$[\ch{W}]$ saturate. Furthermore, as in the case of the Michaelis Menten CRN
(see Sec. V.3.2), the chemostatted concentrations are well approximated by the
coarse-grained dynamics obtained from Eq. (11), as shown by Fig. 17b.

This can be explained by recognizing that the chemostatting procedure does not
break the conservation law $\ell^{E}$ (and the concentrations $[\ch{E}]$,
$[\ch{ES}]$, $[\ch{EW}]$ and $[\ch{E}^{*}]$ are bounded), while it breaks the
conservation laws $\ell^{\ch{F}}$ and $\ell^{\ch{S}}$ with the moiety
concentrations $L^{\text{S}}$ and $L^{\text{F}}$ following as
$d_{t}L^{\ch{S}}=I_{\text{S}}-k^{e}_{\ch{P}}[\ch{P}]$ and
$d_{t}L^{\ch{F}}=I_{\text{F}}-k^{e}_{\ch{W}}[\ch{W}]$. However,

$[\ch{P}],[\ch{W}]\leq\frac{I_{\ch{S}}+I_{\ch{F}}}{\text{min}(k^{e}_{\ch{P}},k^{e}_{\ch{W}})}\,,$
(154)

according to Eq. (127). Hence, only the concentrations $[\ch{F}]$ and
$[\ch{S}]$ can grow, as seen in Fig. 17a. The condition for growth has the
same origin as in the Michaelis Menten CRN under mixed control (see Sec.
V.3.2). The rate of enzyme-mediated interconversion of $\ch{F}$ (resp.
$\ch{S}$) into $\ch{W}$ (resp. $\ch{P}$) is bounded by some function of
$L^{\ch{E}}(0)$. Thus, when $I_{\ch{F}}$ (resp. $I_{\ch{S}}$) is larger than
this rate, $\ch{F}$ (resp. $\ch{S}$) accumulates and the CRN grows.

We now examine the corresponding thermodynamics. The typical evolution of the
EPR (18), $d_{t}G$ (17) and chemical work rate (19) are plotted in Fig. 18a.
All the quantities increase with time and consequently, the efficiency of
growth Eq. (53) saturates to a constant smaller than one (see inset in Fig.
18a), similar to the Michaelis Menten CRN in Fig. 15a.

Figure 18b shows the typical evolution of the nonconservative work rate (22),
the moiety work rate (20) along with the EPR and $d_{t}G$ (when $\ch{S}$ and
$\ch{F}$ are chosen as the potential species): both the nonconservative work
rate and moiety work rate increase with time. Furthermore, we observe a
splitting between the magnitudes of the EPR and the nonconservative work rate,
on the one hand, and the moiety work rate and $d_{t}G$, on the other hand,
implying that

$\begin{split}d_{t}{G}&\sim\dot{w}_{\text{m}}\,,\\\ \dot{w}_{\text{nc}}&\sim
T\dot{\Sigma}>0\,,\\\ \end{split}$ (155)

in the long time limit. We note the similarity between the thermodynamics of
the Metabolic CRN in Figs. 18a, 18b and of the Michaelis Menten CRN in Figs.
15a, 15b.

####  V.3.5 Summary

From the results of Sec. V.3.1 \- V.3.4, we note the following requirements
for growth under mixed control.

Figure 18:  Thermodynamics of the metabolic CRN under mixed control under the
dynamics in Fig. 17. Evolution of a) EPR, the chemical work rate, and $d_{t}G$
and efficiency of growth (inset). b) EPR, moiety work rate, nonconservative
work rate, and $d_{t}G$. Here, $\mu^{0}_{\ch{E}}=1$, $\mu^{0}_{\ch{EF}}=2$,
$\mu^{0}_{\ch{EW}}=2$, $\mu^{0}_{\ch{E^{*}}}=3$, $\mu^{0}_{\ch{F}}=1$,
$\mu^{0}_{\ch{W}}=1$ $\mu^{0}_{\ch{S}}=1$ and $\mu^{0}_{\ch{P}}=1$.

First, from Eq. (127), the CRN needs at least one chemostatted species that is
solely injected (not extracted) to grow. Second, by comparing the
autocatalytic CRN, which does not grow under mixed control (see Sec. V.3.1),
with both the Michaelis Menten CRN and the minimal metabolic CRN, which grow
under mixed control (see Sec. V.3.2 and V.3.4), we expect that CRNs need to
have at least one unbroken conservation law to grow. Finally, given that the
Cyclic Michaelis Menten CRN does not grow under mixed control (see Sec.
V.3.3), we expect that CRNs cannot have any reactions interconverting only the
chemostatted species to grow.

##  VI Conclusions

We analyzed growth in open CRNs under three different chemostatting
procedures. A summary of our results is provided in Fig. 19.

Figure 19: Long time behaviour of the concentrations, EPR and efficiency under
the three different chemostatting procedures for all CRNs. Conditions in
orange show no growth. Conditions in striped blue (resp. solid) show
nonequilibrium (resp. equilibrium) growth. For multimolecular CRNs, scalings
of only a subset of concentrations are depicted. For cases where the exact
scaling of the EPR is not known, nonequilibrium (resp. equilibrium) growth is
indicated by $T\dot{\Sigma}>0$ (resp. $T\dot{\Sigma}\to 0$).

We first showed the importance of the mass conservation law in analyzing
growth. By investigating the dynamics of the mass density, we proved that no
CRN can grow in a Continuous-flow Stirred Tank Reactor setup while any CRN can
grow under flux control

We then analyzed unimolecular and pseudo-unimolecular CRNs. We proved that
these CRNs cannot grow under concentration and mixed control and only grow
under flux control. On one hand, unimolecular CRNs grow close to equilibrium.
We proved that the dissipation scales as $\mathcal{O}(1/t)$ and that the
efficiency of free-energy conversion goes to one. On the other hand, pseudo-
unimolecular CRNs grow far from equilibrium. We proved that the dissipation
scales as $\mathcal{O}(t)$ and that the efficiency goes to zero.

We then studied four multimolecular CRNs with different topological features.
We confirmed numerically and via a fixed point analysis the conjecture that no
growth is possible in the Autocatalytic CRN under concentration control. We
then analyzed growth under flux control. We showed that multimolecular CRNs
display different growth features compared to unimolecular and pseudo-
unimolecular CRNs. For instance, some concentrations scale as
$\mathcal{O}(\sqrt{t})$ in the Autocatalytic CRN and only a subset of
concentrations grow in the Michaelis Menten CRN. Furthermore, the
Autocatalytic CRN and the Cylic Michaelis Menten CRN grow at equilibrium,
while the Michaelis Menten CRN and the Minimal Metabolic CRN grow out of
equilibrium. We finally considered mixed control and found that only the
Michaelis Menten CRN and the Minimal Metabolic CRN can grow. Strikingly, in
both CRNs, some initial conditions relax to steady states while others grow.
When they grow, both CRNs grow out of equilibrium with finite efficiencies.

The analysis of the of the four multimolecular CRNs shows the importance of
unbroken conservation laws for nonequilibrium growth. Indeed, the presence of
unbroken conservation laws is a necessary but not sufficient condition leading
to nonequilibrium growth. On the one hand, they can prevent the equilibration
between the chemostatted species and, consequently, create nonconservative
forces (like in the Michaelis Menten CRN and the Minimal Metabolic CRN). On
the other hand, if there are no unbroken conservation laws (like in the
Autocatalytic CRN) or if there are reaction pathways not involving the
unbroken conservation laws (like in the Cyclic Michaelis Menten CRN), the
chemostatted species can equilibrate leading to equilibrium growth.

Our results on unimolecular and pseudo-unimolecular CRNs rigorously support
the considerations made in the companion paper [23]. We also analyzed in much
greater depth the Autocatalytic CRN and the Michaelis Menten CRN compared to
Ref.[23].

Our work could be extended in various directions. We only considered time-
independent chemostatting procedures but time-dependent ones could produce
interesting results. Sughiyama et al [37] studied the thermodynamics of growth
in detailed balanced CRNs (which always relax to equilibrium), but their
analysis included the effect of volume which we neglected. Adding this effect
in our analysis can in principle be easily done by introducing an equation of
state coupling volume and species concentrations, as done for instance in
Refs. [38, 39]. Autocatalysis did not play an important role in our analysis
of growth. It remains to be seen if autocatalysis may play a more important