# Free-Energy Transduction in Chemical Reaction Networks: from Enzymes to
Metabolism

Artur Wachtel  Department of Molecular, Cellular and Developmental Biology,
Yale University, New Haven (CT), U.S.A.  Riccardo Rao  Simons Center for
Systems Biology, School of Natural Sciences, Institute for Advanced Study,
08540 Princeton (NJ), U.S.A.  Massimiliano Esposito
<EMAIL_ADDRESS>Complex Systems and Statistical Mechanics,
Department of Physics and Materials Science, University of Luxembourg, 162a,
Avenue de la Faïencerie, 1511 Luxembourg, G. D. Luxembourg

###### Abstract

We provide a rigorous definition of free-energy transduction and its
efficiency in arbitrary—linear or nonlinear—open chemical reaction networks
(CRNs) operating at steady state. Our method is based on the knowledge of the
stoichiometric matrix and of the chemostatted species (i.e. the species
maintained at constant concentration by the environment) to identify the
fundamental currents and forces contributing to the entropy production.
Transduction occurs when the current of a stoichiometrically balanced process
is driven against its spontaneous direction (set by its force) thanks to other
processes flowing along their spontaneous direction. In these regimes, open
CRNs operate as thermodynamic machines. After exemplifying these general ideas
using toy models, we analyze central energy metabolism. We relate the
fundamental currents to metabolic pathways and discuss the efficiency with
which they are able to transduce free energy.

##  I Introduction

Living organisms are out-of-equilibrium systems which harness chemical free
energy from their surrounding to fuel their own function [1]. To do so, free
energy from external energy-rich molecules is transferred into energetic
mediator molecules—primarily adenosine triphosphate (ATP)—which are then used
to power the otherwise-unfavorable reactions that are essential for life. The
resulting low-energy molecules and heat are released back into the environment
as waste. At the core of this process lies energy metabolism: a finely-
regulated network of coupled biochemical reactions that—through enzymatic
catalysis—enable transformations on the environment that would hardly occur
spontaneously [2, 3].

In order to study the thermodynamics of living systems, traditional chemical
thermodynamics is of limited use. Indeed it considers transformations between
initial and final equilibrium distributions of chemicals (typically the
equilibrium state before and after reactant and product are brought in
contact) and time is absent from the description. In contrast, energy
metabolism is a dynamical process that never reaches equilibrium, except upon
death of the living organism. The thermodynamics of living systems is thus
much more akin to that of thermodynamic machines, which transduce free energy
thanks to reservoirs at different temperatures, pressures or chemical
potentials. These chemical machines must be studied using modern formulations
of nonequilibrium thermodynamics, as these formulations combine
phenomenological irreversible thermodynamics with dynamical rate laws. To
analyze their performance, identifying the conjugated pairs of thermodynamic
force and current which determine the entropy production (or dissipation) is a
crucial step. While this is easily done for macroscopic machines [4], as well
as for simple biological or synthetic molecular machines [5, 6, 7, 8, 9, 10,
11, 12], slightly more complex machines described by stochastic dynamics on
finite state graphs already require some graph theoretical methods, like those
pioneered by T.L. Hill in his work on free-energy transduction [5] as well as
by J. Schnakenberg [13, 14]. Furthermore, for chemical reaction networks
(CRNs), these theories are restricted to chemical systems with elementary
first or pseudo-first order reactions, i.e. systems whose deterministic
dynamics in terms of concentrations is linear and equivalent to the dynamics
for the probability to find a molecule of a given type. These restrictions can
be justified when studying the energetics of single enzymes, but become
completely inadequate to describe even simple metabolic networks, where many
reactions contain multimolecular steps and are described by enzymatic rate
laws (i.e. for non-elementary reactions).

In this paper, we generalize these ideas, leaving these limitations behind. We
provide a general theory of free-energy transduction for any open CRN. To do
so, we make use of the topology of the networks (encoded in the stoichiometric
matrix) and of the chemostatted species to determine the fundamental
thermodynamic forces and currents that are necessary to define free-energy
transduction and its efficiency. We then demonstrate the value of the theory
by analyzing free-energy transduction in central metabolism, highlighting the
relation between metabolic pathways and fundamental thermodynamic currents.

In section II we review the dynamics and nonequilibrium thermodynamics of open
CRNs. In particular, we show how the stoichiometric matrix can be used to
identify the fundamental thermodynamic forces and currents controlling the
dissipation. In section III we show how to make use of these results to
analyze simple and illustrative CRNs, thus paving the way for a systematic
study of free-energy transduction and its efficiency. In section IV we show
that our method is not limited to the analysis of simple systems by applying
it to central energy metabolism. We first consider glycolysis, the
tricarboxylic acid cycle and the electron transport chain, separately and
combined. Ethanol fermentation and gluconeogenesis are subsequently discussed.
We highlight the relation between metabolic pathways and fundamental currents
and analyze their efficiency. In section V we formalize the systematic
procedure to analyze free-energy transduction and its efficiency in arbitrary
CRNs. Section VI contains a discussion and the conclusions.

##  II Chemical Reaction Networks

After defining the stoichiometry and kinetics, we review how to use
thermodynamics to analyze the dissipation caused by exchanges of free energy
with the environment [15, 16, 6, 17, 18].

###  II.1 Stoichiometry

A CRN describes the interconversion of a set of chemical species $\ce{Z1}$, …,
$\ce{Z_{n}}$ by means of a set of chemical reactions
$\rho_{1},\dots,\rho_{m}$. Each reaction $\rho_{j}$ is characterized by an
equation

$\displaystyle\nu_{1,j}\ce{Z_{1}}+\dots+\nu_{n,j}\ce{Z_{n}}$
$\displaystyle\rightleftharpoons\bar{\nu}_{1,j}\ce{Z_{1}}+\dots+\bar{\nu}_{n,j}\ce{Z_{n}}\,,$
(1)

which encodes that each time reaction $\rho_{j}$ occurs, it will destroy
$\nu_{i,j}\geq 0$ reactant molecules of type $\ce{Z_{i}}$ and produce
$\bar{\nu}_{k,j}\geq 0$ product molecules of type $\ce{Z_{k}}$. Their
differences $\mathbb{S}_{i,j}=\bar{\nu}_{i,j}-\nu_{i,j}$ are the
_stoichiometric coefficients_ that represent the net production of each
species in each reaction. We combine the latter into an $n\times m$ matrix
$\mathbb{S}$ called _stoichiometric matrix_.

As imposed by the conservation laws of physical chemistry, we always assume
that any such reaction is balanced in terms of electrons and atomic nuclei,
and consequently also balanced in terms of mass and charge: every proton needs
to be explicitly included in the stoichiometry. Furthermore every reaction is
reversible: each reaction $\rho_{+j}$ has a naturally associated reverse
reaction $\rho_{-j}$ where reactants and products are interchanged.

###  II.2 Dynamics

In addition to the stoichiometry encoded in the stoichiometric matrix
$\mathbb{S}$, each reaction $\rho_{j}$ comes with a _reaction current_
$J_{j}(\bm{z})$. These currents characterize the net rate of reaction when the
species concentrations are $\bm{z}$. The precise expression of
$\bm{J}(\bm{z})$ is determined by the reaction _kinetics_ , which depends on
both the type of reaction (elementary, enzyme-catalyzed, etc.), and on the
type of solution (ideal dilute solution, non-ideal dilute solution, etc.).

The reaction currents and the stoichiometric coefficients govern the dynamics
of the concentrations $\bm{z}$ via the deterministic differential equations,

$\displaystyle\dot{\bm{z}}=\mathbb{S}\bm{J}(\bm{z})+\bm{I}(\bm{z})\,,$ (2)

known as the _rate equations_ of chemical kinetics. In open CRNs,
$\bm{I}(\bm{z})$ accounts for the external sources and sinks. $\bm{I}(\bm{z})$
is absent in closed CRNs. Both here and in the following it is sufficient for
our purposes to limit our discussion to steady states, where concentrations
are stationary: $\bm{0}=\dot{\bm{z}}=\mathbb{S}\bm{J}+\bm{I}$ (the dependence
on $\bm{z}$ is sometimes omitted for brevity). We will briefly discuss
extensions beyond stationarity in the discussion. In case of stationarity,
reaction currents and external currents have to balance out:
$\mathbb{S}\bm{J}=-\bm{I}$. These steady states are in general out of
thermodynamic equilibrium.

###  II.3 Thermodynamics

We consider dilute solutions kept at constant temperature $T$. For each
chemical species $i$, the chemical potential $\mu_{i}(\bm{z})$ characterizes
the free-energy content of that species when the concentrations are $\bm{z}$
111For isobaric solutions, $\mu_{i}(\bm{z})$ must be regarded as the Gibbs
free-energy content of species $i$. As for the kinetics, the precise
expression of $\mu_{i}(\bm{z})$ depends on the type of solution. For instance,
for ideal dilute solutions, $\mu_{i}(\bm{z})$ solely depends on the
concentration of $i$,

$\displaystyle\mu_{i}(\bm{z})=\mu^{\circ}_{i}+RT\ln
z_{i}=RT\ln\frac{z_{i}}{z_{i,\mathrm{eq}}}\,,$ (3)

where $\mu^{\circ}_{i}$ is the standard-state chemical potential of species
$i$, $R$ is the gas constant; see e.g. [20]. For non-ideal dilute solutions
$\mu_{i}(\bm{z})$ additionally depends on the concentrations of all chemicals
via an activity coefficient $\gamma_{i}(\bm{z})$, which effectively accounts
for the non-negligible interactions between chemicals,

$\mu_{i}(\bm{z})=\mu^{\circ}_{i}+RT\ln z_{i}+RT\ln\gamma_{i}(\bm{z})\,.$ (4)

The stoichiometric difference of chemical potentials serves as a thermodynamic
force for the reactions:

$\displaystyle-\Delta_{j}G(\bm{z})=-\mathbb{S}_{j}^{\mathsf{T}}\bm{\mu}(\bm{z})\,,$
(5)

where $\mathbb{S}_{j}$ is the $j$th colum of the stoichiometric matrix. The
overall dissipation in the reaction network is given by the product of forces
and currents for all reactions:

$\displaystyle T\sigma=-\sum_{j}J_{j}\,\Delta_{j}G\geq 0\,.$ (6)

The _entropy production rate_ (EPR) $\sigma$ quantifies the rate of entropy
change in the system and in the reservoirs (i.e. the solution and the
chemostats) [16]. The second law of thermodynamics guarantees that $\sigma$ is
non-negative, and that it vanishes solely at thermodynamic equilibrium, where
the thermodynamic forces of all reactions vanish.

In general, the _dissipative contribution_ of each reaction,
$-J_{j}\,\Delta_{j}G$, need not be non-negative. As we will discuss soon, the
presence of negative contributions is indicative of transduction mechanisms.
At the same time, there are important types of reaction for which
$-J_{j}\,\Delta_{j}G$ is always non-negative: elementary reactions in ideal
and non-ideal solutions [16, 18], as well as cytosolic enzymes [6, 17]. For
these types of reaction, the kinetics and thermodynamics are linked by the
_local detailed balance conditions_

$\displaystyle J_{j}(\bm{z})$
$\displaystyle=\phi_{+j}(\bm{z})-\phi_{-j}(\bm{z})$ (7a)
$\displaystyle-\Delta_{j}G(\bm{z})$
$\displaystyle=RT\ln\frac{\phi_{+j}(\bm{z})}{\phi_{-j}({\bm{z}})}\,,$ (7b)

where $\phi_{+j}(\bm{z})$ (resp. $\phi_{-j}(\bm{z})$) is the precise rate at
which the forward $+j$ (resp. reverse $-j$) reaction occurs. For example, for
elementary reactions in ideal solutions, the reaction rates follow the
mass–action kinetics: $\phi_{+j}(\bm{z})=k_{+j}\prod_{i}z_{i}^{\nu_{i,j}}$ and
$\phi_{-j}(\bm{z})=k_{-j}\prod_{i}z_{i}^{\bar{\nu}_{i,j}}$, where $k_{+j}$ and
$k_{-j}$ are the rate constants. When Eqs. (7) are satisfied, it is easy to
see that $J_{j}$ and $-\Delta_{j}G$ are aligned for any $\bm{z}$, and hence
$-J_{j}\,\Delta_{j}G\geq 0$. Since any chemical process can be described as a
network of elementary chemical reactions, the aforementioned fact constitutes
a proof that $\sigma\geq 0$.

Any thermodynamically consistent description ensures that when the system is
closed, $\bm{I}=0$, the steady state corresponds to thermodynamic equilibrium
where all reaction currents vanish, $\bm{J}(\bm{z})=\bm{0}$, as well as the
dissipation, $\sigma=0$. Indeed, the local detailed conditions (7) do so. In
the special case of mass–action kinetics, they also imply Wegscheider’s
conditions [21, 22]).

###  II.4 Network topology and dissipation

We now make use of the topology of the CRN encoded in its stoichiometric
matrix. To that end, we split the chemical concentrations into two groups:
$\bm{z}=(\bm{x},\bm{y})$. The species $\ce{X}$ are entirely _internal_ to the
reaction network and cannot be exchanged with the outside world, i.e.
$\bm{I}^{\ce{X}}=0$. Consequently, their dynamics are entirely due to
reactions:

$\displaystyle\dot{\bm{x}}$
$\displaystyle=\mathbb{S}^{\ce{X}}\,\bm{J}(\bm{x},\bm{y})\,.$ (8)

In contrast, the _chemostatted_ species $\ce{Y}$ are controlled by the
environment:

$\displaystyle\dot{\bm{y}}$
$\displaystyle=\mathbb{S}^{\ce{Y}}\,\bm{J}(\bm{x},\bm{y})+\bm{I}^{\ce{Y}}\,.$
(9)

Such control can result from external reservoirs fixing their chemical
potentials (or equivalently their concentrations in an ideal solution) [18] or
from mechanisms fixing their in- or out-fluxes $\bm{I}^{\ce{Y}}$ [23]. Note
that the splitting of the species results in a corresponding splitting of the
stoichiometric matrix
$\mathbb{S}=(\mathbb{S}^{\ce{X}},\mathbb{S}^{\ce{Y}})^{\mathsf{T}}$.

At steady state, these two chemostatting procedures coincide and the reaction
currents satisfy

$\displaystyle 0$ $\displaystyle=\mathbb{S}^{\ce{X}}\,\bm{J}^{\ast}\,,$ (10)
$\displaystyle 0$
$\displaystyle=\mathbb{S}^{\ce{Y}}\,\bm{J}^{\ast}+\bm{I}^{\ce{Y}}\,.$ (11)

The steady-state currents thus need to be in the null-space of the reduced
stoichiometric matrix $\mathbb{S}^{\ce{X}}$:
$\bm{J}^{\ast}\in\ker\mathbb{S}^{\ce{X}}$ at steady state. Importantly, this
space contains the null-space of $\mathbb{S}$,
$\ker\mathbb{S}\subseteq\ker\mathbb{S}^{\ce{X}}$, whose elements represent
combinations of reaction currents that when taken together neither produce nor
consume any molecule: they return the system back to its initial state. We
call these _internal cycles_ and let $\set{\bm{c}^{\ell}}$ denote a maximal
independent set of these—i.e. a basis of $\ker\mathbb{S}$. The null-space
$\ker\mathbb{S}^{\ce{X}}$ is spanned by all combinations of reactions that
leave the concentrations of internal species unvaried. This encompasses all
internal cycles, but also accounts for additional cycles: The steady state
must exchange chemostatted species in order to faithfully represent a
nonequilibrium state. We call these _additional_ cycles _emergent_ , and
denote a maximal independent set of them by $\set{\bm{C}^{\varepsilon}}$.
Their stoichiometry is given by
$\mathbb{S}\bm{C}=(\bm{0},\mathbb{S}^{\ce{Y}}\bm{C})$. The internal species do
not appear and only the chemostatted species remain. The positive terms are
products and the negative terms are reactants. Every emergent cycle thus
defines an _effective reaction_ amongst the chemostatted species, which is
fully balanced since each reaction in the network it was derived from was
fully balanced. Altogether they characterize the thermodynamic transformations
that the CRN performs on the environment.

The thermodynamic importance of internal and emergent cycles becomes apparent
when decomposing the steady-state currents
$\bm{J}^{\ast}=\sum_{\ell}J_{\ell}\bm{c}^{\ell}+\sum_{\varepsilon}J_{\varepsilon}\bm{C}^{\varepsilon}$
according to internal and emergent cycles. Using this decomposition, the
dissipation simplifies to

$\displaystyle
T\sigma=-\sum_{j}J^{*}_{j}\Delta_{j}G=-\sum_{\varepsilon}J_{\varepsilon}\Delta_{\varepsilon}G\,.$
(12)

In the last equality, we exploited the fact that the internal cycles forces
vanish

$\displaystyle-\sum_{j}c_{j}\,\Delta_{j}G=-\bm{c}\mathbb{S}^{\mathsf{T}}\bm{\mu}=-\left(\mathbb{S}\bm{c}\right)^{\mathsf{T}}\bm{\mu}=0\,.$
(13)

Therefore, while the internal cycle currents $J_{\ell}$ do not vanish in
general, they do not produce dissipation. Instead, the emergent cycle forces

$\displaystyle-\Delta_{\varepsilon}G=-\sum_{j}\bm{C}^{\varepsilon}_{j}\Delta_{j}G=-\sum_{j}\bm{C}_{j}^{\varepsilon}\mathbb{S}^{\mathsf{T}}_{j}\bm{\mu}=(\mathbb{S}^{\ce{Y}}\bm{C}^{\varepsilon})^{\mathsf{T}}\bm{\mu}^{\ce{Y}}\,,$
(14)

typically do not vanish. They only involve the chemical potentials of the
chemostatted species since the contributions of internal species, by
construction, cancel in the summation.

The dissipation in a nonequilibrium steady state is thus entirely captured by
the currents and forces along the emergent cycles [15]. From a thermodynamic
standpoint they thus constitute the fundamental set of currents and forces in
the entire open CRN.

We emphasize that our framework holds far-from-equilibrium in regimes where
the fluxes are nonlinear functions of the forces. Linear irreversible
thermodynamics is recovered close to equilibrium, where the fluxes become
linear in the forces
$J_{\varepsilon}\approx\sum_{\varepsilon^{\prime}}\mathcal{L}_{\varepsilon\varepsilon^{\prime}}\Delta_{\varepsilon^{\prime}}G$,
where $\mathcal{L}$ is the Onsager matrix of the open CRN.

##  III Free-energy transduction in simple models

Before proceeding, we define a _chemical process_ as a generic stoichiometric
relation between reactants and products, like Eq. (1). Such a stoichiometric
relation must be balanced with respect to electrons and atomic nuclei—and thus
charge and mass—but the path leading from reactants to products can be
arbitrary. In this way, a chemical process is defined independently from the
kinetics that characterises its rate of occurrence. The stoichiometry of the
processes is however sufficient to determine its thermodynamic force via Eq.
(5), since chemical potentials are state functions.

An example of a chemical process is the oxidation of glucose:

$\displaystyle\ce{Glu+6O2}\rightleftharpoons\ce{6CO2+6H2O}\,.$ (15)

This process can be the result of combustion or of cellular respiration. While
the kinetics is different in these two cases, we still consider it to be the
same chemical process with the same thermodynamic force. An emergent cycle
$\bm{C}^{\varepsilon}$ also defines a chemical process.

We now have all the ingredients at hand to characterize the free-energy
transduction of chemical processes. In this section we start by discussing
simple models that may describe enzymatic mechanisms or membrane transport
processes. In the next section, we will scale up our analysis to metabolic
processes.

###  III.1 Independent processes

Figure 1:  (Left) Closed CRN made of two _independent_ reactions. (Right) The
two emergent cycles that arise when the CRN is opened by chemostatting the
species A, B, C, and D.

Let us consider the CRN

$\displaystyle\ce{A}\rightleftharpoons\ce{D}\ \ \ \ \ \ \ \
\ce{B}\rightleftharpoons\ce{C}$ (16)

depicted in Fig. 1 (left). Each reaction could consist of many intermediate
steps, as long as these are _independent_ from the other reaction.
Furthermore, since each reaction must preserve mass, they can be thought of as
different conformations of a molecule. The CRN is opened by chemostatting the
species A, B, C, and D as depicted in Fig. 1 (right). The two emergent cycles
in this system are identical to the two reactions (16).

We now consider the two independent chemical processes:
$p_{1}\colon\ce{A}\rightleftharpoons\ce{D}$ and
$p_{2}\colon\ce{B}\rightleftharpoons\ce{C}$. We note that in this trivial
example, reactions, emergent cycles and processes all coincide. These two
processes happen with net currents $J^{\text{ind}}_{p_{1}}$ and
$J^{\text{ind}}_{p_{2}}$, and their direction is determined by the external
forces,

$\displaystyle-\Delta_{p_{1}}G=\mu_{\ce{A}}-\mu_{\ce{D}}\,,\quad-\Delta_{p_{2}}G=\mu_{\ce{B}}-\mu_{\ce{C}}\,.$
(17)

Both processes reach thermodynamic equilibrium if and only if
$\mu_{\ce{A}}=\mu_{\ce{D}}$ and $\mu_{\ce{B}}=\mu_{\ce{C}}$.

The total entropy-production rate is the sum of the rates associated with each
of the two emergent cycles, which coincide with the two independent processes:

$\displaystyle
T\sigma^{\text{ind}}=-J^{\text{ind}}_{p_{1}}\Delta_{p_{1}}G-J^{\text{ind}}_{p_{2}}\Delta_{p_{2}}G\,.$
(18)

The assumption of independence between the two reactions allows us to state
that each of the processes needs to satisfy the second law independently.
Hence,

$\displaystyle
T\sigma^{\text{ind}}_{p_{1}}\coloneqq-J^{\text{ind}}_{p_{1}}\Delta_{p_{1}}G\geq
0\,,\quad
T\sigma^{\text{ind}}_{p_{2}}\coloneqq-J^{\text{ind}}_{p_{2}}\Delta_{p_{2}}G\geq
0\,.$ (19)

This makes the two independent processes purely dissipative, and implies that
each reaction current is following the direction of its thermodynamic force.
Free-energy transduction cannot happen in this case.

We will now consider different ways in which these two processes might be
coupled resulting in transduction.

###  III.2 Tightly coupled processes

Figure 2:  (Left) Closed CRN made of three reactions. (Right) The single
emergent cycle that arises when A, B, C, and D are chemostatted. This CRN
implements the same two processes as in Fig. 1, but with intermediate steps
that couple the two processes _tightly_ , i.e. via a fixed stoichiometry.

Let us assume that the same four species A, B, C, and D we introduced in
§III.1 are the only chemostats in the following CRN with arbitrary kinetics:

$\displaystyle\ce{A+B}\rightleftharpoons\ce{AB}\rightleftharpoons\ce{CD}\rightleftharpoons\ce{C+D}\,.$
(20)

We consider the same two processes as before:
$p_{1}\colon\ce{A}\rightleftharpoons\ce{D}$ and
$p_{2}\colon\ce{B}\rightleftharpoons{C}$. This time the CRN has only one
emergent cycle with stoichiometry

$\displaystyle\ce{A+B}\rightleftharpoons\ce{C+D}\,,$ (21)

as illustrated in Fig. 2 to highlight the difference with respect to the
previous case.

The two individual processes $p_{1}\colon\ce{A}\rightleftharpoons\ce{D}$ and
$p_{2}\colon\ce{B}\rightleftharpoons{C}$ come now bound with a fixed
stoichiometric ratio of unity: a molecule of type D can only be created when
also a molecule of type C is created. Therefore, we say that the two processes
are _tightly coupled_. The force of the emergent cycle is the sum of the
individual forces:

$\displaystyle-\Delta_{\text{tight}}G=\mu_{\ce{A}}+\mu_{\ce{B}}-\mu_{\ce{C}}-\mu_{\ce{D}}=-\Delta_{p_{1}}G-\Delta_{p_{2}}G\,.$

Let $J^{\text{tight}}$ be the current of the emergent cycle in the steady
state, then the entropy-production rate is

$\displaystyle
T\sigma=-J^{\text{tight}}\Delta_{\text{tight}}G=-J^{\text{tight}}\Delta_{p_{1}}G-J^{\text{tight}}\Delta_{p_{2}}G\,.$
(22)

This expression is, in total, always non-negative: $T\sigma\geq 0$. Hence, we
can deduce that either

  1. 1.

$-\Delta_{p_{1}}G=\Delta_{p_{2}}G$ and the system is at equilibrium,
$J^{\text{tight}}=0$ .

  2. 2.

$-\Delta_{p_{1}}G,-\Delta_{p_{2}}G<0$ or $-\Delta_{p_{1}}G,-\Delta_{p_{2}}G>0$
which makes the system purely dissipative.

  3. 3.

$-\Delta_{p_{1}}G$ and $-\Delta_{p_{2}}G$ have different sign. Then the system
can be understood to transduce free energy. In the following we assume,
without restriction, that $-\Delta_{p_{1}}G\,,J^{\text{tight}}>0$ while
$-\Delta_{p_{1}}G>\Delta_{p_{2}}G>0$.

In the latter situation, we can say that the tight coupling mechanism uses the
free energy released by the process $p_{1}$ as input to drive the (output)
process $p_{2}$ uphill in free energy:

$\displaystyle
T\sigma^{\text{tight}}_{p_{1}}\coloneqq-J^{\text{tight}}\Delta_{p_{1}}G>0\,,\quad
T\sigma^{\text{tight}}_{p_{2}}\coloneqq-J^{\text{tight}}\Delta_{p_{2}}G<0\;.$
(23)

This example represents the simplest form of a chemical machine. The machine
itself is given by the internal species AB and CD and operates on an
environment given by the chemostats A, B, C, and D. The efficiency of free-
energy transduction for this machine is the ratio of the output over the input
free energy,

$\displaystyle\eta^{\text{tight}}\coloneqq-\frac{T\sigma^{\text{tight}}_{p_{2}}}{T\sigma^{\text{tight}}_{p_{1}}}=\frac{\Delta_{p_{2}}G}{-\Delta_{p_{1}}G}\leq
1\,.$ (24)

This efficiency is bounded from above by unity, since
$T\sigma^{\text{tight}}=T\sigma^{\text{tight}}_{p_{1}}+T\sigma^{\text{tight}}_{p_{2}}\geq
0$. Maximal efficiency is achieved only at equilibrium, where both
$J^{\text{tight}}=0$ and $-\Delta_{\text{tight}}G=0$ wich implies
$T\sigma^{\text{tight}}_{p_{1}}=T\sigma^{\text{tight}}_{p_{2}}=0$ while
$-\Delta_{p_{1}}G=\Delta_{p_{2}}G$ may be arbitrarily large.

Also note that $\eta^{\text{tight}}$ is entirely expressed in terms of the
thermodynamic forces. The kinetics of the coupling process only affect the
current, which cancels out. This is the distinctive feature of tightly coupled
processes.

###  III.3 Loosely coupled processes

Figure 3:  Closed CRN describing two loosely coupled processes mediated by an
enzymatic mechanism. The open CRN obtained by chemostatting A, B, C and D
(Fig. 4) is characterized by two linearly independent emergent cycles. One
choice of basis for the emergent cycles is depicted in Fig. 4.

A different form of coupling can be achieved by considering some form of
catalyst, $\ce{E}$, that couples the conversion of $\ce{A}$ and $\ce{B}$ into
$\ce{D}$ and $\ce{C}$ using the CRN:

$\displaystyle\ce{A+E}\rightleftharpoons\ce{EA}$
$\displaystyle\rightleftharpoons\ce{ED}\rightleftharpoons\ce{E+D}$ (25)
$\displaystyle\ce{B+EA}\rightleftharpoons\ce{BEA}$
$\displaystyle\rightleftharpoons\ce{CED}\rightleftharpoons\ce{C+ED}$ (26)

also depicted in Fig. 3.

With the species A, B, C, and D chemostatted, we have two emergent cycles and
we have to choose a basis. The three possible choices are:

  * •

$\ce{A}\rightleftharpoons\ce{D}$ and $\ce{B}\rightleftharpoons\ce{C}$

  * •

$\ce{A+B}\rightleftharpoons\ce{C+D}$ and $\ce{B}\rightleftharpoons\ce{C}$

  * •

$\ce{A+B}\rightleftharpoons\ce{C+D}$ and $\ce{A}\rightleftharpoons\ce{D}$

Figure 4: The two coupled emergent cycles $\ce{A}\rightleftharpoons\ce{D}$ and
$\ce{B}\rightleftharpoons\ce{C}$ of the catalytic reaction network shown in
Fig. 3.

Since we are still interested in the chemical processes
$p_{1}\colon\ce{A}\rightleftharpoons\ce{D}$ and
$p_{2}\colon\ce{B}\rightleftharpoons\ce{C}$—without loss of generality—we
choose the first choice of emergent cycles, which are depicted in Fig. 4. In
this way each emergent cycle carries a process. These cycles have the same
forces as the independent ones in section III.1. The fact that there is a
kinetic coupling does not affect the thermodynamic equilibrium, nor the
thermodynamic forces. However, two currents are needed now to characterize
dissipation:

$\displaystyle T\sigma^{\text{cat}}$
$\displaystyle=-J^{\text{cat}}_{p_{1}}\Delta_{p_{1}}G-J^{\text{cat}}_{p_{2}}\Delta_{p_{2}}G$
(27)
$\displaystyle=T\sigma^{\text{cat}}_{p_{1}}+T\sigma^{\text{cat}}_{p_{2}}\,.$
(28)

Note that the two currents are _not_ the currents coming from the kinetics of
the independent processes used in §III.1, but from the catalytic coupling and
thus depend on the total amount of available catalyst.

Again, since the overall dissipation is non-negative, we have three cases: (1)
thermodynamic equilibrium, (2) pure dissipation, or (3) transduction. In the
latter case we assume, without restriction, $T\sigma^{\text{cat}}_{p_{1}}>0$
while $T\sigma^{\text{cat}}_{p_{2}}<0$. Then the efficiency of free-energy
transduction is

$\displaystyle\eta^{\text{cat}}=-\frac{T\sigma^{\text{cat}}_{p_{2}}}{T\sigma^{\text{cat}}_{p_{1}}}=\frac{J^{\text{cat}}_{p_{2}}\Delta_{p_{2}}G}{-J^{\text{cat}}_{p_{1}}\Delta_{p_{1}}G}=\eta^{\text{tight}}\frac{J^{\text{cat}}_{p_{2}}}{J^{\text{cat}}_{p_{1}}}\leq
1\,.$ (29)

This efficiency is different from the tightly coupled case: It depends on the
kinetics via the currents. In general, there is no way to tell which of the
two efficiencies is larger.

In this example, the internal species all contain the catalyst E in different
states of binding to other molecules. The catalyst can thus be seen as the
machine acting on an environment given by the chemostatted species.

###  III.4 Hybrid cases

Figure 5:  A closed CRN representing a hybrid between the two CRNs given in
Figs. 2 and 3. The intermediate steps of one conversion process acts as a
catalyst for another. Chemostatting the six base species results in the two
emergent cycles shown in Fig. 6.

We now consider the CRN depicted in Fig. 5 that is a hybrid between the
tightly and loosely coupled CRNs introduced earlier

$\displaystyle\ce{A+B}\rightleftharpoons\ce{AB}$
$\displaystyle\rightleftharpoons\ce{CD}\rightleftharpoons\ce{C+D}$ (30)
$\displaystyle\ce{AB+U}\rightleftharpoons\ce{ABU}$
$\displaystyle\rightleftharpoons\ce{CDV}\rightleftharpoons\ce{CD+V}$ (31)

With the six chemostats $\ce{A}$, $\ce{B}$, $\ce{C}$, $\ce{D}$, $\ce{U}$,
$\ce{V}$, we have two emergent cycles

  * •

$\varepsilon_{1}\colon\ce{A+B}\rightleftharpoons\ce{C+D}$

  * •

$\varepsilon_{2}\colon\ce{U}\rightleftharpoons\ce{V}$

as depicted in Fig. 6, and the dissipation reads

$\displaystyle
T\sigma^{\text{hyb}}=T\sigma^{\text{hyb}}_{\varepsilon_{1}}+T\sigma^{\text{hyb}}_{\varepsilon_{2}}\,.$
(32)

Figure 6: The two emergent cycles of the chemostatted network introduced in
Fig. 5. (left) Emergent cycle $\varepsilon_{1}$ which tightly couples the
processes $p_{1}$ and $p_{2}$. (right) Emergent cycle $\varepsilon_{2}$ which
transforms U into V.

As before, we are interested in the chemical processes:
$p_{1}\colon\ce{A}\rightleftharpoons\ce{D}$ and
$p_{2}\colon\ce{B}\rightleftharpoons\ce{C}$ which are both tightly coupled by
the emergent cycle $\varepsilon_{1}$. But contrary to the tight coupling case
we now have an additional dissipative process transforming U into V on the
emergent cycle $\varepsilon_{2}$.

The dissipation can thus be expressed in terms of _three_ contributions:

$\displaystyle
T\sigma^{\text{hyb}}=T\sigma^{\text{hyb}}_{p_{1}}+T\sigma^{\text{hyb}}_{p_{2}}+T\sigma^{\text{hyb}}_{\varepsilon_{2}}\,.$
(33)

Transduction arises when $T\sigma^{\text{hyb}}_{p_{2}}<0$, but we need to
distinguish two cases. If the current of the additional process follows its
force, $T\sigma^{\text{hyb}}_{\varepsilon_{2}}>0$, it should be treated as an
input process, and the efficiency is

$\displaystyle\eta^{\text{hyb}}$
$\displaystyle=-\frac{T\sigma^{\text{hyb}}_{p_{2}}}{T\sigma^{\text{hyb}}_{p_{1}}+T\sigma^{\text{hyb}}_{\varepsilon_{2}}}=\frac{J^{\text{hyb}}_{\varepsilon_{1}}\Delta_{p_{2}}G}{-J^{\text{hyb}}_{\varepsilon_{1}}\Delta_{p_{1}}G-J^{\text{hyb}}_{\varepsilon_{2}}\Delta_{\varepsilon_{2}}G}$
(34)
$\displaystyle=\eta^{\text{tight}}\left(1+\frac{J^{\text{hyb}}_{\varepsilon_{2}}\Delta_{\varepsilon_{2}}G}{J^{\text{hyb}}_{\varepsilon_{1}}\Delta_{p_{1}}G}\right)^{-1}=\eta^{\text{tight}}\left(1+\frac{\sigma^{\text{hyb}}_{\varepsilon_{2}}}{\sigma^{\text{hyb}}_{p_{1}}}\right)^{-1}\,,$
(35)

where we used the fact that since the coupling between $p_{1}$ and $p_{2}$ is
tight, their currents will be equal in strength. Since
$\sigma_{\varepsilon_{2}}>0$, the efficiency of transduction in this hybrid
case is smaller than if we only had $p_{1}$ tightly coupled to $p_{2}$. If
instead $T\sigma_{p_{2}}<0$ and $T\sigma_{\varepsilon_{2}}<0$, then
$\varepsilon_{2}$ should be treated as an additional output process, and the
efficiency should be defined as

$\displaystyle\eta^{\text{hyb'}}=-\frac{T\sigma_{p_{2}}+T\sigma_{\varepsilon_{2}}}{T\sigma_{p_{1}}}=\eta^{\text{tight}}-\frac{\sigma_{\varepsilon_{2}}}{\sigma_{p_{1}}}\,.$
(36)

In this case, since multiple output processes are obtained from the same
input, the efficiency increases with respect to the tightly coupled case.

###  III.5 Counter-Example to Transduction: Mass transfer

Figure 7: Mass transfer is a counter example to transduction.

In order to define transduction, it is crucial to to make sure that the two
processes—inputs and outputs—are fully balanced in terms of electrons and
atomic nuclei. Here we highlight this fact.

Let us re-consider the example of _tight coupling_ as given in Eq. (20) or
Fig. 2. Only looking at this reaction network or its emergent cycle
stoichiometry Eq. (21) cannot guarantee that the species $\ce{A}$ and $\ce{D}$
have the same charge and atomic composition, nor that there is an actual
process that could transform $\ce{A}$ into $\ce{D}$. Let us assume for
instance that $\ce{A}=\ce{DU}$ and $\ce{C}=\ce{BU}$, meaning that the net
stoichiometry of the emergent cycle would be

$\displaystyle\ce{DU}+\ce{B}\rightleftharpoons\ce{BU}+\ce{D}\,,$ (37)

as depicted in Fig. 7. This represents the net transfer of the moiety $\ce{U}$
from $\ce{D}$ to $\ce{B}$. In this case, there is no way to decompose this
process into two independent chemical processes and to naturally decompose the
entropy-production rate as we did in Eq. (22). Mass transfer does not
correspond to free-energy transduction.

###  III.6 Nonlinear Case Study

Figure 8:  Nonlinear closed CRN given by the reactions (38) and (39).

The models discussed so far are pseudo-unimolecular: when considering the
chemostatted species constant in time, all reactions are unimolecular and
dynamically linear in the internal species. But the theory in section II, as
well as our coming main result in section V, do not rely on this fact.

To illustrate this point, we consider the non-linear CRN

$\displaystyle\ce{A}\rightleftharpoons\ce{X}_{1}\ \;\
\ce{B}\rightleftharpoons\ce{X3}\ \;$ $\displaystyle\
\ce{C}\rightleftharpoons\ce{X}_{1}\ \;\ \ce{D}\rightleftharpoons\ce{X}_{2}$
(38) $\displaystyle\ce{X}_{1}+\ce{X}_{2}$
$\displaystyle\rightleftharpoons\ce{X}_{3}+\ce{X}_{4}$ (39)

depicted in Fig. 8. Its central reaction (39) takes two internal species into
two different internal species, making it dynamically nonlinear even when the
chemostats A, B, C, and D have fixed concentrations. Indeed, the chemostatting
produces the single emergent cycle
$\ce{A}+\ce{B}\rightleftharpoons\ce{C}+\ce{D}$ depicted in Fig. 9.

Figure 9: Emergent cycle $\ce{A}+\ce{B}\rightleftharpoons\ce{C}+\ce{D}$
obtained when chemostatting (A, B, C, D) in the CRN of Fig. 8.

This example provides an alternative way to tightly couple the processes
$p_{1}$ and $p_{2}$. From the perspective of the reservoirs, it is essentially
identical to the system in Fig. 2. However, the nonlinear nature of its
internal dynamics will have an impact on the global kinetics of the system:
The total current through this system has a different dependence on the
chemostats’ concentrations than the linear system discussed earlier. Thus, the
global kinetics can be a way to differentiate the two different
implementations of tight coupling.

##  IV Free-energy transduction in central metabolism

We now apply the framework introduced in the previous sections to a
biochemically relevant case: central energy metabolism. The function of
central energy metabolism is the transfer of energy from energy-rich
molecules—most importantly glucose—into $\ce{ATP}$ [2, 3].

We first consider cellular respiration and its three components: glycolysis,
the tricarboxylic acid (TCA) cycle, and the electron transport chain with ATP
synthesis. We show that the corresponding metabolic pathways are recovered as
emergent cycles, and hence they can be interpreted as chemical processes. From
the perspective of free-energy transduction, we find that all these energy-
conversion pathways are tightly coupled. This allows us to estimate the
efficiency of these chemical machines from the $\Delta G$s alone.

We then conclude our analysis by considering two extensions of glycolysis:
fermentation and gluconeogenesis. While the former can be interpreted as a
tightly-coupled process, the latter cannot, being composed of several separate
processes.

We refer to App. A for a detailed discussion of all individual reaction steps.

###  IV.1 Cellular Respiration

With its high yield of ATP, cellular respiration plays a central role in
energy metabolism, and is shared among most cellular life forms with
essentially the same structure. The three components of respiration are
glycolysis in which glucose is converted into pyruvate while synthesizing
ATPs, TCA cycle in which pyruvate is used to reduce the coenzymes CoQ and
NAD+, and electron transport chain in which the reducing power of NADH is used
to synthesize ATP. We will be discussing these pathways individually in the
following sections.

Importantly, when all the reactions of these three pathways (see (80a–j), (81)
and (82)) are taken together as part of one single network and when H2O, H+,
ATP, ADP, P_i, O2, CO2, Glu are regarded as chemostatted, then one observes
the single emergent cycle

$\displaystyle\ce{Glu+6O2+26ADP+26P_{\mathrm{i}}+26H+}$

$\rightleftharpoons$

(40) $\displaystyle\ce{6CO2+32H2O+26ATP}$

depicted in Fig. 10.

Figure 10:  The single emergent cycle of cellular respiration is formed by a
combination of glycolysis, the TCA cycle, as well as the electron transport
chain. Each reaction has indicated multiplicities, so as to generate the net
balance given in (40). The coloring highlights our interpretation of two
coupled processes in the net reaction: red/blue is the synthesis of ATP and
green/brown is the breakdown of glucose (Glu) into two molecules of pyruvate
(Pyr).

The choice of the chemostatted species is suggested by viewing cellular
respiration as a cellular process operating on the environment. Glu, O2, CO2,
H2O, and H+ are chemicals that cells exchange with the environment, whereas
the synthesis of ATP from ADP and P_i is the desired cellular chemical output.

Following this reasoning and using the information about the balance of
charges and atomic nuclei, we recognize that the net stoichiometric balance of
cellular respiration is the result of two chemical processes tightly coupled
together: the burning of glucose into carbon dioxide and the synthesis of ATP,

$\displaystyle\ce{Glu+6O2}\rightleftharpoons\ce{6CO2+6H2O}$ (41)
$\displaystyle\ce{ADP+P_{\mathrm{i}}+H+}\rightleftharpoons\ce{ATP+H2O}\,,$
(42)

The stoichiometric coupling of 26 ATP per Glu in central energy metabolism
depends on the species [3]. Here we consider the bacterium _E. coli_. For
eukaryotes, the organization of central energy metabolism is slightly
different: for them the TCA cycle and the electron transport chain are
performed in the mitochondria. The transport into and out of mitochondria is
partially active and thus itself consumes some $\ce{ATP}$. Bacteria do not
have specialized organelles for this function. This difference, together with
sightly different enzymes for the electron transport chain, results in a
different global stoichiometry.

Under standard physiological conditions (see App. A), the full oxidation of
glucose (41) releases about
$2910\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ ($1180\,RT$), while the
free energy of the $\ce{ATP}$ hydrolysis reaction is
$46\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ ($18.4\,RT$). Thus, the
theoretical (physical) limit to the conversion rate is roughly 64 $\ce{ATP}$
per 1 glucose. Considering that cellular respiration produces 26 $\ce{ATP}$
per glucose, this makes for an efficiency of about 41%:

$\eta_{\mathrm{CR}}=\frac{26\,\Delta_{\mathrm{ATPsynth}}G}{-\Delta_{\mathrm{GluOxidat}}G}=0.41\,,$
(43)

where

$\displaystyle-\Delta_{\mathrm{GluOxidat}}G$
$\displaystyle=\mu_{\ce{Glu}}+6\mu_{\ce{O2}}-6\mu_{\ce{CO2}}-6\mu_{\ce{H2O}}$
(44) $\displaystyle=$2910\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$$
$\displaystyle-\Delta_{\mathrm{ATPsynth}}G$
$\displaystyle=\mu_{\ce{ADP}}+\mu_{\ce{Pi}}+\mu_{\ce{H+}}-\mu_{\ce{ATP}}-\mu_{\ce{H2O}}$
(45) $\displaystyle=-$46\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$\,.$

We now turn to the individual analysis of the three pathways constituting
cellular respiration.

####  IV.1.1 Glycolysis

Glycolysis is the first part of cellular respiration, and its individual steps
are the enzymatic reactions depicted in Fig. 11, see also Eqs. (80a–j). The
chemostatting of Glu, Pyr, ATP, ADP, P_i, NAD+, NADH, H+ and H2O, produces a
single emergent cycle, which re-traces the traditional understanding of the
glycolysis pathway. It begins with one molecule of glucose (Glu) and ends with
two molecules of pyruvate (Pyr), producing a net gain of 2 ATP per glucose. In
order to realize this process, the glucose is phosphorylated twice before it
is broken into two 3-carbon bodies, DHAP and GAP. These can then be easily
converted into each other. In GAP form, the two pieces are phosphorylated
again via the help of NAD+, which is reduced in the process. Finally, each of
the carbon bodies is de-phosphorylated twice, resulting in the total
production of 4 ATP, leaving a net gain of 2 ATP per glucose.

Figure 11:  Single emergent cycle of glycolysis which couples the synthesis of
ATP (red/blue) and the breakdown of glucose into two molecules of pyruvate
(green/brown).

The net balance of glycolysis coincides with the effective reaction that the
emergent cycle performs on the chemostats

$\displaystyle\ce{Glu+2ADP+2P_{\mathrm{i}}+2NAD+}$

$\rightleftharpoons$

(46) $\displaystyle\ce{2Pyr+2ATP+2NADH+2H++2H2O}\,.$

We can easily identify the synthesis of ATP from ADP as one of the processes
involved, Eq. (42). The other process is instead the conversion of glucose
into pyruvate:

$\displaystyle\ce{Glu+2NAD+}$
$\displaystyle\rightleftharpoons\ce{2Pyr+2NADH+4H+}$ (47)

Under standard physiological conditions (details in App. A) the conversion of
glucose into pyruvate, Eq. (47), releases roughly
$170\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ ($68\,RT$). Recalling
that (i) the free energy required for ATP synthesis, Eq. (45), is
$46\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ ($18.4\,RT$), and (ii)
two ATPs are hydrolyzed during glycolysis, one can estimate the efficiency of
this pathway,

$\eta_{\mathrm{glyco}}=\frac{2\,\Delta_{\mathrm{ATPsynth}}G}{-\Delta_{\ce{Glu}\rightarrow\ce{Pyr}}G}=0.54\,,$
(48)

with

$\displaystyle-\Delta_{\ce{Glu}\rightarrow\ce{Pyr}}G$
$\displaystyle=\mu_{\ce{Glu}}+2\mu_{\ce{NAD+}}-2\mu_{\ce{Pyr}}-2\mu_{\ce{NADH}}-4\mu_{\ce{H+}}$
(49) $\displaystyle=$170\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$\,.$

Note that the production of NADH in Eq. (47) is a side-effect. This buildup is
generally considered the main reason why glycolysis needs to be coupled with
additional reactions in order to achieve a proper cyclic process. We will
discuss small additions to glycolysis in sections IV.2 and IV.3, alcoholic
fermentation and gluconeogenesis, respectively.

####  IV.1.2 Tricarboxylic Acid Cycle

The tricarboxylic acid (TCA) cycle is a series of enzymatic reactions (81)
that gradually decompose pyruvate into CO2 while synthesizing ATP and reducing
NAD+ to NADH, as well as reducing ubiquinone (CoQ) to ubiquinole (CoQH2). The
CRN of the TCA cycle is depicted in Fig. 12 together with the single emergent
cycle resulting from chemostatting Pyr, ATP, ADP, P_i, NADH, NAD+, CO2, CoQ,
CoQH2, H+, and H2O given by

$\displaystyle\ce{ADP+CoQ+2H2O+4NAD++Pyr+P_{\mathrm{i}}}$

$\rightleftharpoons$

(50) $\displaystyle\ce{ATP+3CO2+CoQH2+2H++4NADH}$

This cycle tightly couples two processes: the synthesis of ATP, Eq. (42), and
the anaerobic decomposition of pyruvate

$\displaystyle\ce{Pyr+3H2O+4NAD++CoQ}$

$\rightleftharpoons$

(51) $\displaystyle\ce{3CO2+4NADH+CoQH2+3H+}$

Figure 12:  Single emergent cycle of the tricarboxylic acid (TCA) cycle which
couples the synthesis of ATP (red/blue) to the decomposition of Pyruvate
(green/brown).

The free energy released by the decomposition of pyruvate (standard
physiological conditions) is about
$150\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ ($60\,RT$), see App. A.
Since the whole process is tightly coupled, and just one ATP is hydrolized for
any pyruvate decomposed, one finds a modest efficiency:

$\eta_{\mathrm{TCA}}=\frac{\Delta_{\mathrm{ATPsynth}}G}{-\Delta_{\mathrm{PyrDec}}G}=$0.30\text{\,}\,$,$
(52)

with

$\displaystyle-\Delta_{\mathrm{PyrDec}}G$
$\displaystyle=\mu_{\ce{Pyruvate}}+3\mu_{\ce{H2O}}+\mu_{\ce{CoQ}}+4\mu_{\ce{NAD+}}$
$\displaystyle\;-3\mu_{\ce{CO2}}-\mu_{\ce{CoQH2}}-4\mu_{\ce{NADH}}-3\mu_{\ce{H+}}$
(53) $\displaystyle=$150\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$\,.$

####  IV.1.3 Electron Transport Chain

Figure 13:  Emergent cycle of the electron transport chain couples the
oxidation of NADH (green/brown) to the synthesis of ATP (red/blue) via a
proton motive force (concentration gradient of protons across the membrane).
Note that CoQ and CoQH2 remain exactly balanced in this cycle. For that
reason, if they are regarded as internal species, the electron transport chain
is characterized solely by this single emergent cycle. Instead, if they are
chemostatted, two emergent cycles arise, this one and a second one depicted in
Fig. 14.  Figure 14:  With ubiquinol and ubiquinone chemostatted, a second
emergent cycle appears in the electron transport chain. It oxidizes ubiquiol
(CoQH2) and thus restores ubiquinone (CoQ) for the TCA cycle.

The last part of cellular respiration is the electron transport chain,
depicted in Figs. 13 and 14, see also (82). It uses the coenzyme Q in its
reduced (CoQH2, ubiquinol) and oxidized (CoQ, ubiquinone) forms as an
intermediate carrier of hydrogen to oxidize NADH and thus to restore NAD+. In
doing so, the electron transport chain pumps protons from the cytoplasm (H+)
into the periplasm ($\ce{H}^{+}_{\mathrm{p}}$). The resulting proton motive
force is then used by ATP synthase to synthesize ATP from ADP.

We note that, when considering the electron transport chain in isolation, it
makes most sense to consider the two forms of coenzyme Q as internal species.
In this case the network gives rise to a single emergent cycle (Fig. 13). The
net effect of this emergent cycle is the synthesis of two ATP for each
oxidized NADH:

$\displaystyle\ce{2ADP+3H++NADH+\tfrac{1}{2}O2+2P_{\mathrm{i}}}$

$\rightleftharpoons$

(54) $\displaystyle\ce{2ATP+3H2O+NAD+}\,.$

This tightly couples two processes: the synthesis of ATP, Eq. (42), and the
oxidation of NADH,

$\displaystyle\ce{NADH+\tfrac{1}{2}O2+H+}$
$\displaystyle\rightleftharpoons\ce{NAD++H2O}\,.$ (55)

The free energy released upon oxidation of NADH is about
$217\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ ($88\,RT$), and the
efficiency is estimated to be 42%:

$\eta_{\mathrm{ETC}}=\frac{2\Delta_{\mathrm{ATPsynth}}G}{-\Delta_{\mathrm{NADHoxidat}}G}=$0.42\text{\,}\,$,$
(56)

with

$\displaystyle-\Delta_{\mathrm{NADHoxidat}}G$
$\displaystyle=\mu_{\ce{NADH}}+\tfrac{1}{2}\mu_{\ce{O2}}+\mu_{\ce{H+}}-\mu_{\ce{NAD+}}-\mu_{\ce{H2O}}$
(57) $\displaystyle=$217\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$\,.$

The emergent cycle of cellular respiration (40) is composed of reactions from
all of the three pathways, and different reactions need to be performed with
different multiplicity in order to balance the internal species. It is however
not possible to deduce its net balance solely from the emergent cycles of the
three isolated pathways (46), (50), and (54). The reason is that to properly
connect the electron transport chain to the TCA cycle it is necessary to also
treat both forms of coenzyme Q as chemostatted in the former. Doing so gives
rise to a second emergent cycle in the electron transport chain, depicted in
Fig. 14, which generates one ATP per CoQH2:

$\displaystyle\ce{ADP+H++CoQH2+\tfrac{1}{2}O2+P}_{\mathrm{i}}$

$\rightleftharpoons$

(58) $\displaystyle\ce{ATP+2H2O+CoQ}\,.$

As a result, the emergent cycle of cellular respiration (40) is now reproduced
by composing: $1\times$ glycolysis (46), $2\times$ the TCA cycle (50),
$10\times$ the first emergent cycle of the electron transport chain (54), and
$2\times$ the second one (58).

###  IV.2 Ethanol Fermentation

Fermentation is a way to re-cycle NADH when operating under anaerobic
conditions. This allows the cell to use glycolysis for the synthesis of ATP in
a cyclic fashion. In Fig. 15 we show the reactions of alcoholic fermentation
as an emergent cycle.

Figure 15:  Single emergent cycle of alcoholic fermentation extending
glycolysis which is the result of a tight coupling of two processes: synthesis
of ATP (red/blue) and decomposition of glucose into ethanol and carbon dioxide
(green/brown).

It shares most of its reactions with regular glycolysis, but extends them by
two steps: cleavage of CO2 from pyruvate and subsequent addition of hydrogen
via oxidation of NADH, (83). Thus, when considering the total net balance,
ethanol fermentation reads

$\displaystyle\ce{Glu+2ADP+2P_{\mathrm{i}}+2H+}$

$\rightleftharpoons$

(59) $\displaystyle\ce{2Ethanol+2CO2+2ATP+2H2O}$

We see that NAD+ and NADH need not be regarded as chemostatted: they are now
internal to the reaction network. In contrast, CO2 and Ethanol take their role
as chemostats. In this form it is easy to see how fermentation tightly couples
the two processes

$\displaystyle\ce{ADP+P_{\mathrm{i}}+H+}$
$\displaystyle\rightleftharpoons\ce{ATP+H2O}$ $\displaystyle\ce{Glu}$
$\displaystyle\rightleftharpoons\ce{2Ethanol+2CO2}$ (60)

with a stoichiometric ratio of 2:1.

Since the conversion of glucose into ethanol and CO2, Eq. (60), releases
roughly $267\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ ($108\,RT$), the
overall efficiency of such fermentative ATP production is

$\eta_{\mathrm{EF}}=\frac{2\Delta_{\mathrm{ATPsynth}}G}{-\Delta_{\ce{Glu}\rightarrow\ce{Eth}}G}=$0.34\text{\,}\,$,$
(61)

with

$\displaystyle-\Delta_{\ce{Glu}\rightarrow\ce{Eth}}G$
$\displaystyle=\mu_{\ce{Glu}}-2\mu_{\ce{Eth}}-2\mu_{\ce{CO2}}-4\mu_{\ce{H+}}$
(62) $\displaystyle=$267\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$\,.$

###  IV.3 Gluconeogenesis and Futile Cycles

The pathway of gluconeogenesis enables cells to synthesize glucose from
pyruvate, see (80). The pathway of gluconeogenesis is however not the exact
reverse of glycolysis as it involves some different enzymatic reactions and
couples to different molecules. These modifications lead to a different
emergent cycle, depicted in Fig. 16, with net balance

$\displaystyle\ce{2Pyr+2NADH+2HCO3+2ATP+2GTP+2H++2H2O}$

$\rightleftharpoons$

(63) $\displaystyle\ce{Glu+2NAD++2CO2+2ADP+2GDP+4P_{\mathrm{i}}}\,.$

Gluconeogenesis can be split into two processes: the synthesis of glucose from
pyruvate

$\displaystyle\ce{2Pyr+2NADH+4H+}$
$\displaystyle\rightleftharpoons\ce{Glu+2NAD+}\,,$ (64)

as well as

$\displaystyle\ce{ATP+GTP+HCO3+H2O}$

$\rightleftharpoons$

(65) $\displaystyle\ce{ADP+GDP+2P_{\mathrm{i}}+CO2+H+}\,.$

Here, (64) serves as the output process, while (65) is the input process. Note
that the output process of gluconeogenesis (i.e. the synthesis of glucose) is
the reverse of the input process of glycolysis (47). However, the input
process of gluconeogenesis (65) involves different molecules than the output
process of glycolysis. Therefore, gluconeogenesis and glycolysis are not exact
reversals of each other.

We note that the input process (65) could in fact be decomposed further into
three separate processes: the hydrolysis of ATP, the hydrolysis of GTP, as
well as the dehydration of carbonic acid. But since all three are input
processes, separating them all only clutters the notation and adds nothing to
the present discussion of gluconeogenesis.

Figure 16:  Emergent cycle of the combined glycolysis and gluconeogenesis
network which represents the gluconeogenesis pathway. The additional reactions
of glycolysis are faded out.  Figure 17:  The internal cycle of the combined
glycolysis and gluconeogenesis network.  Figure 18:  The emergent cycle of the
combined glycolysis and gluconeogenesis network that is commonly referred to
as futile.

When the CRNs of glycolysis and gluconeogenesis are combined, one obtains
three emergent cycles and an internal one. The three emergent cycles are given
by glycolysis (Fig. 11), gluconeogenesis (Fig. 16), and a _futile cycle_ that
just hydrolyzes ATP to ADP without coupling to additional conversions (Fig.
18)

$\displaystyle\ce{ATP+H2O}$
$\displaystyle\rightleftharpoons\ce{ADP+P_{\mathrm{i}}+H+}\,.$ (66)

The internal cycle (which is also present if the CRN is closed) is depicted in
Fig. 17. One can easily verify that contrary to emergent cycles, upon
completing every reaction in the internal cycle, not only internal species but
also chemostatted ones are left unchanged.

Whether one considers the performance of ATP synthesis (42) or of glucose
synthesis (64), tight coupling is broken when the reactions of glycolysis and
gluconeogenesis are combined. Indeed, the dissipation is now carried by three
emergent cycles:

$T\sigma=-J_{\mathrm{glyco}}\Delta_{\mathrm{glyco}}G-J_{\mathrm{futile}}\Delta_{\mathrm{futile}}G-J_{\mathrm{gluco}}\Delta_{\mathrm{gluco}}G\,,$
(67)

where $J_{\mathrm{glyco}}$, $J_{\mathrm{futile}}$, and $J_{\mathrm{gluco}}$
are the currents associated with each emergent cycle. The corresponding forces
read

$\displaystyle-\Delta_{\mathrm{glyco}}G$
$\displaystyle=-\Delta_{\mathrm{Glu}\rightarrow\mathrm{Pyr}}G-2\Delta_{\mathrm{ATPsynth}}G$
(68) $\displaystyle-\Delta_{\mathrm{futile}}G$
$\displaystyle=\Delta_{\mathrm{ATPsynth}}G$ (69)
$\displaystyle-\Delta_{\mathrm{gluco}}G$
$\displaystyle=-\Delta_{\mathrm{Pyr}\rightarrow\mathrm{Glu}}G+2\Delta_{\mathrm{ATPsynth}}G\,,$
(70)

with (49), (45), and

$\displaystyle-\Delta_{\mathrm{Pyr}\rightarrow\mathrm{Glu}}G$
$\displaystyle=\Delta_{\mathrm{Glu}\rightarrow\mathrm{Pyr}}G$ $\displaystyle\
\
+2(\mu_{\ce{GTP}}+\mu_{\ce{HCO3}}-\mu_{\ce{CO2}}-\mu_{\ce{GDP}}-\mu_{\ce{P_{\mathrm{i}}}})\,.$
(71)

If we now want to characterize the efficiency of ATP synthesis, we need to
account for the fact that the output process is now carried by three emergent
cycles. Therefore,

$\eta_{\mathrm{ATPsynth}}=\frac{(2J_{\mathrm{glyco}}-2J_{\mathrm{gluco}}-J_{\mathrm{futile}})\Delta_{\mathrm{ATPsynth}}G}{-J_{\mathrm{glyco}}\Delta_{\ce{Glu}\rightarrow\ce{Pyr}}G-J_{\mathrm{gluco}}\Delta_{\ce{Pyr}\rightarrow\ce{Glu}}G}\,.$
(72)

If $J_{\mathrm{futile}}$ and $J_{\mathrm{gluco}}$ vanish, we recover the tight
coupling efficiency of glycolysis alone, Eq. (48). Under normal physiological
conditions, glycolysis and gluconeogenesis are tightly regulated so that they
will not occur at the same time [3]. Futile cycles are also known to be
regulated. Characterizing the effect of these regulations on the efficiency
would be very interesting but is beyond the scope of this paper.

We note that most literature on cellular respiration considers the combination
of the internal and futile cycle (Figs. 17 and 18) as being two different
futile cycles, which both have the net stoichiometry of ATP hydrolysis. Seeing
these two cycles as two futile cycles driven by the same force, or as an
internal cycle and a futile cycle, is very much a question of taste. None of
the dynamical characterizations will depend on this choice. Even the total
dissipation rate will be invariant under this choice. However, the values of
the individual cycle currents and cycle forces will change.

##  V Defining free-energy transduction and its efficiency

Our work results in the following prescription for the identification of free-
energy transduction in CRN described by a stoichiometric matrix $\mathbb{S}$:

  1. 1.

Find the internal cycles in the network—i.e. combinations of reactions that
leave the concentrations of all species unchanged—by finding a basis spanning
the kernel of $\mathbb{S}$: $\set{\bm{c}^{\ell}}$ such that
$\mathbb{S}\bm{c}^{\ell}=\bm{0}$.

  2. 2.

Identify which chemical species are exchanged with the environment—i.e. the
chemostatted species $\ce{Y}$—and split
$\mathbb{S}=(\mathbb{S}^{\ce{X}},\mathbb{S}^{\ce{Y}})^{\mathsf{T}}$
accordingly: $\mathbb{S}^{\ce{X}}$ identifies the block corresponding to
internal species and $\mathbb{S}^{\ce{Y}}$ to chemostatted species,
respectively.

  3. 3.

Find the emergent cycles of the network—i.e. combinations of reactions that
change the concentrations of the chemostatted species but not those of the
internal species—by finding the _additional_ vectors to $\set{\bm{c}^{\ell}}$
needed to span the kernel of $\mathbb{S}^{\ce{X}}$: $\bm{C}$ such that
$\mathbb{S}\bm{C}=(\bm{0},\mathbb{S}^{\ce{Y}}\bm{C})\neq\bm{0}$. The net
stoichiometries of the emergent cycles are given by
$\mathbb{S}^{\ce{Y}}\bm{C}$.

  4. 4.

Compute, measure or estimate the currents along the emergent cycles.

  5. 5.

Use your knowledge about the charge and atomic composition of the involved
molecules to decompose the emergent cycles further: some emergent cycles may
be the result of several tightly-coupled processes.

  6. 6.

Write the total dissipation as the product of forces and currents of the
emergent cycles

$\displaystyle
T\sigma=-\sum_{\text{cycles}}\Delta_{\text{cyc}}G\,J_{\text{cyc}}\,,$ (73)

and apply the additional decomposition of cycles into processes on the forces:

$\displaystyle
T\sigma=-\sum_{\text{proc}}\Delta_{\text{proc}}G\,J_{\text{proc}}\,.$ (74)

In this decomposition all the process currents $J_{\text{proc}}$ that belong
to the same emergent cycle are identical to the emergent cycle current
$J_{\text{cyc}}$.

  7. 7.

Identify the sign of the dissipative contributions
$\set{-\Delta_{\text{proc}}G\,J_{\text{proc}}}$. Three cases are possible.

    1. (i)

All contributions are positive. The CRN is purely dissipative.

    2. (ii)

Some contributions are positive, some are negative. The CRN is transducing
free energy. The positive contributions are the input, the negative
contributions are the output,

$\displaystyle T\sigma=T\sigma_{\text{input}}+T\sigma_{\text{output}}\geq
0\,,$ (75)

and the efficiency of free-energy transduction is

$\displaystyle\eta=-\frac{T\sigma_{\text{output}}}{T\sigma_{\text{input}}}\leq
1\,.$ (76)

    3. (iii)

All contributions vanishing individually is highly unlikely as it would
require immense fine tuning of external concentrations to equilibrium values.

In our considerations we have found one special case for Eq. (76): Tight
coupling. In this case only one emergent cycle appears, and its current will
cancel in the ratio—making the efficiency depend on the ratio of chemical
potential differences and thus energetics alone. In the general case, the
efficiency of free-energy transduction depends both on energetics via chemical
potentials as well as on the kinetics via the currents of the involved
processes.

##  VI Discussion and conclusions

We provided the conceptual tools to analyze complex open CRNs as free-energy-
transducing chemical machines. Our approach is based on the knowledge of the
stoichiometric matrix and of the species exchanged with the environment (i.e.
the chemostatted species). At steady state, only the forces corresponding to
the emergent cycles (combinations of chemical potentials of chemostatted
species) and their corresponding currents (combinations of currents of
chemostatted species) determine the dissipation, see Eq. (12). These emergent
currents and forces are crucial to define transduction and its efficiency:
they are needed to identify the chemical processes that serve as input and
output in the dissipation.

An important outcome of this work is that transduction is generally an
emergent phenomenon that occurs at the network level when multiple elementary
reactions are coupled to each other. It typically cannot happen at the level
of single elementary reactions, which very often involve mass transfer (see
Sec. III.5). Two coupled isomerizations may however be seen as an exception to
the rule.

We focused our study on steady states, but generalizations to transient and
nonautonomous regimes (e.g. situations in which the concentrations of the
chemostats changes in time) are possible [16, 24, 18]. These generalization
have been used for coarse-graining schemes [17], as well as in the context of
energy storage [25].

The framework that we propose in this article provides a general theory of
free-energy transduction valid for _nonlinear_ CRNs. Theories of free-energy
transduction have been previously established [5], but they were limited to
linear CRNs. These theories have been successful in studying molecular motors
because the enzymatic dynamics is well described by a linear network in which
every transition corresponds to a pseudo-unimolecular reaction [6, 7, 8, 9,
10, 11, 12]. By formulating a theory for nonlinear CRN, we paved the way to
analyze much more complex CRNs, such as metabolic networks. These networks
contain numerous bi-molecular reactions that could not be described by linear
networks.

When applied to metabolic CRNs, our analysis shows that the emergent cycles
coincide with known metabolic pathways. We found that, when taken either
separately or together, the three main parts of central energy
metabolism—glycolysis, the TCA cycle, and the electron transport
chain—correspond to a single emergent cycle, which is distinctive of tightly-
coupled free-energy transduction. We could therefore study their transduction
efficiency solely based on stoichiometric and thermodynamic considerations.
However, tight coupling is broken when additional pathways are added—such as
gluconeogenesis—as they give rise to additional and possibly futile emergent
cycles. In this case, the transduction efficiencies depend on the currents,
and assessing them requires some information about the kinetics. This goes
beyond the scope of this study but can be envisaged in the future.

Of course, metabolic regulation plays a key role in metabolism and properly
describing this aspect remains an open challenge. Our current thermodynamics
framework remains applicable to study any post-translational regulation
because such a regulation only affects the kinetics—and not the forces. In
contrast, translational regulations have their own energetics, which would
need to be accounted for. However, the latter operate on longer times scales,
and our method could still provide good descriptions for short time scales.

Our method could also be of value for metabolic network reconstruction
approaches such as flux or energy balance analysis [26, 27, 28], which seek to
determine the reaction currents based on constrained optimizations. In this
context, thermodynamics has been recently recognized as an important factor to
constrain reconstruction methods [29, 30, 31, 32] but also to study metabolism
in general [33, 34]. Further work is required to elucidate these connections.

The largest network we considered in this paper is cellular metabolism. We now
briefly comment on some of the questions that may arise when considering CRNs
of larger size and complexity, such as genome-wide metabolic networks or
possibly atmospheric and biogeochemical CRNs. Up to small modifications, the
basic principle of our approach should continue to apply. The question of
which species to treat as chemostatted is largely contingent. It is a special
instance of the classical question in thermodynamics of what to treat as the
system and what to treat as the environment. For example, if the changes in
concentration of an internal species in large abundance can be neglected over
the relevant time scales of the experiment, one can treat it as chemostatted.
The same would be true if external control mechanisms may be at work to
control its concentration. For a given choice of chemostats, the procedure to
find the emergent cycles is simple and can be implemented using a simple
linear algebra algorithm. The task is indeed to find a set of independent
integer vectors that belong to the null space of the open-system
stoichiometric matrix but not to the null space of the closed-system
stoichiometric matrix. As already discussed, this basis is not unique, and
some bases may provide clearer insights about how the CRN transduces free
energy than other bases. When composing cellular respiration in terms of
glycolysis, TCA cycle and electron transport chain, we have seen that our
methodology can be used to study larger CRNs in terms of smaller open CRNs by
chemostatting suitable internal species. This suggest that modular approaches
to complex CRN may be developed in the future.

We end by noting that the present work also provides insights on the validity
of the second law (6) and of the local detailed balance (7) for non-elementary
reactions. Indeed, since the latter should be viewed as processes (including
emergent cycles), the validity of the second law (6) is still ensured by (12),
but not the local detailed balance which imposes the reaction to flow along
its force (see Eq. (7) and below). Instead, as we have seen, processes may
undergo transduction and flow against their force if they are coupled to other
processes. To treat a process as if it were an elementary reaction satisfying
local detailed balance, one must make sure that when chemostatting its
products and reactants, only a single emergent cycle arises which does not
share any internal species with any other process in the CRN. This happens for
many enzymatic reactions for instance. A more systematic discussion on this
can be found in Refs. [6, 17].

## Acknowledgements

This work was supported by grants from the Simons Foundation (691552, RR), and
the European Research Council, project NanoThermo (ERC-2015-CoG Agreement No.
681456, ME).

##  Appendix A Reactions of the Central Energy Metabolism

We expand on the thermodynamic data and the metabolic chemical reactions used
in section IV.

###  A.1 Thermodynamics data under physiological conditions

reactant |  $\mu$ ($\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$) |  | reactant |  $\mu$ ($\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$)
---|---|---|---|---
ATP | $-2281\pm 3$ |  |  CO2 | $-403\pm 6$
ADP | $-1406\pm 2$ |  |  O2 | $-1\pm 6$
AMP | $-531\pm 2$ |  |  CoQ | $790\pm 20$
P_i | $-1072.6\pm 1.5$ |  |  CoQH2 | $770\pm 20$
NAD+ | $-1146\pm 13$ |  | Glucose | $-409.4\pm 1.3$
NADH | $-1080\pm 13$ |  | Pyruvate | $-355.2\pm 1.5$
H2O | $-151.9\pm 1.5$ |  | Ethanol | $65\pm 3$
Table 1:  Values of chemical potentials used to estimate free energy changes and efficiencies and collected from _eQuilibrator_ [35].  Metabolic Process |  $-\Delta G$ ($\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$) | $\eta$
---|---|---
Glucose oxidation | 2910 $\pm$ 50 ($=1180\pm 20\,RT$) | –
ATP synthesis | -46 $\pm$ 4 ($=18.4\pm 1.7\,RT$) | –
Cellular respiration | 1730 $\pm$ 120 ($=700\pm 50\,RT$) | 0.41 $\pm$ 0.04
Glu$\rightarrow$Pyr | 170 $\pm$ 40 ($=68\pm 15\,RT$) | –
Glycolysis | 80 $\pm$ 40 ($=31\pm 15\,RT$) | 0.54 $\pm$ 0.12
Pyruvate decomposition | 150 $\pm$ 80 ($=60\pm 30\,RT$) | –
TCA cycle | 110 $\pm$ 80 ($=40\pm 30\,RT$) | 0.30 $\pm$ 0.16
NADH oxidation | 217 $\pm$ 18 ($=88\pm 8\,RT$) | –
Electron transport chain | 130 $\pm$ 20 ($=51\pm 8\,RT$) | 0.42 $\pm$ 0.05
Glu$\rightarrow$Eth | 267 $\pm$ 13 ($=108\pm 5\,RT$) | –
Ethanol fermentation | 175 $\pm$ 16 ($=71\pm 6\,RT$) | 0.34 $\pm$ 0.03
Table 2:  Main metabolic processes discussed in section IV and estimates of
the corresponding free energy changes under standard physiological conditions.
For tightly-coupled transduction, we also report the estimated efficiency. All
errors are evaluated using the standard propagation of errors (square root of
the sum of squares).

In Tab. 1, we report the chemical potentials of the chemostatted species used
to estimate thermodynamic forces and efficiencies. These values capture
typical physiological conditions—which we referred to as “standard
physiological conditions”—and are collected from _eQuilibrator_ [35]. To be
specific, the activity coefficient $\gamma_{i}(\bm{z})$ in the chemical
potential for non-ideal dilute solutions, Eq. (4), is characterized using the
phenomenological _extended Debye–Hückel relation_

$\ln\gamma_{i}(\bm{z})=-\frac{aq_{i}^{2}\sqrt{I}}{1+b\sqrt{I}}\,,$ (77)

where $q_{i}$ is the charge of the molecule $i$, $I$ is the _ionic strength_
of the solution, and $a$ and $b$ are phenomenological coefficients, see e.g.
Ref. [20, §3.6]. Since protons are often regarded as chemostatted—pH fixed—it
is convenient to report the values of the chemical potentials using pH as a
reference [20, §3.6] using the so-called transformed chemical potentials [20,
§4.1]

$\mu_{i}(\bm{z})=\mu_{i}(\bm{z})-N_{i}\mu_{\ce{H+}}(\bm{z})\,,$ (78)

where $N_{i}$ is the number of proton moieties that the molecule $i$ contains,
and

$\mu_{\ce{H+}}(\bm{z})=\mu^{\circ}_{\ce{H+}}+RT\mathrm{pH}\ln
10+RT\ln\gamma_{\ce{H+}}(\bm{z})$ (79)

is the chemical potential of free protons written in terms of pH. A similar
discussion holds for magnesium ions Mg++, which bind to $\ce{ATP}$ ions and
are also regarded as chemostatted. The values of chemical potential reported
in Tab. 1 correspond to the choice pH $=7.5$, pMg $=3.0$, $I=0.25$,
$z_{i}=10^{-3}$ M (millimolar concentrations, except for water,
$z_{\ce{H_{2}O}}\simeq 55$ M), and $T=298.15$K ($RT\simeq 2.5$), see Tab. 1.

In Tab. 2 we report the overall free energy change, Eq. (5), and the
thermodynamic efficiency, Eq. (24), for the different tightly-coupled chemical
processes discussed in the main text.

###  A.2 Glycolysis and Gluconeogenesis

The enzymatic reactions constituting glycolysis are

$\displaystyle\ce{ATP+Glu}$ $\displaystyle\rightleftharpoons\ce{ADP+G{6}P+H+}$
(80a) $\displaystyle\ce{G{6}P}$ $\displaystyle\rightleftharpoons\ce{F{6}P}$
(80b) $\displaystyle\ce{ATP+F{6}P}$
$\displaystyle\rightleftharpoons\ce{ADP+F{1,\\!6}BP+H+}$ (80c)
$\displaystyle\ce{F{1,\\!6}BP}$ $\displaystyle\rightleftharpoons\ce{DHAP+GAP}$
(80d) $\displaystyle\ce{DHAP}$ $\displaystyle\rightleftharpoons\ce{GAP}$ (80e)
$\displaystyle\ce{GAP+NAD++P_{\mathrm{i}}}$
$\displaystyle\rightleftharpoons\ce{{1,\\!3}BPG+H++NADH}$ (80f)
$\displaystyle\ce{{1,\\!3}BPG+ADP}$
$\displaystyle\rightleftharpoons\ce{{3}PG+ATP}$ (80g)
$\displaystyle\ce{{3}PG}$ $\displaystyle\rightleftharpoons\ce{{2}PG}$ (80h)
$\displaystyle\ce{{2}PG}$ $\displaystyle\rightleftharpoons\ce{H2O+PEP}$ (80i)
$\displaystyle\ce{ADP+H++PEP}$ $\displaystyle\rightleftharpoons\ce{ATP+Pyr}$
(80j) Additionally, gluconeogenesis comprises the reactions
$\displaystyle\ce{ATP+HCO3+Pyr}$
$\displaystyle\rightleftharpoons\ce{ADP+H++OAA+P_{\mathrm{i}}}$ (80k)
$\displaystyle\ce{GTP+OAA}$ $\displaystyle\rightleftharpoons\ce{CO2+GDP+PEP}$
(80l) $\displaystyle\ce{F{1,\\!6}BP+H2O}$
$\displaystyle\rightleftharpoons\ce{F{6}P+P_{\mathrm{i}}}$ (80m)
$\displaystyle\ce{G{6}P+H2O}$
$\displaystyle\rightleftharpoons\ce{Glu+P_{\mathrm{i}}}$ (80n)

When the species H2O, H+, ATP, ADP, P_i, NAD+, NADH, GTP, GDP, CO2, HCO3, Glu,
Pyr are chemostatted, we find three emergent cycles in the glycolysis &
gluconeogenesis network:

The first cycle is a _futile cycle_ with net balance corresponding to
$\ce{ATP}$ hydrolysis. The second cycle is the glycolysis pathway decomposing
glucose ($\ce{Glu}$) into two pyruvate ($\ce{Pyr}$) while phosphorylating two
$\ce{ADP}$ into two $\ce{ATP}$. The last cycle is the gluconeogenesis pathway
producing glucose from pyruvate. In addition, this system has a purely
internal cycle composed of the reactions
$\text{(A4a)}+\text{(A4n)}-\text{(A4m)}-\text{(A4c)}$.

###  A.3 Tricarboxylic Acid (TCA) Cycle

The TCA cycle is composed of the Pyruvate dehydrogenase reaction

$\displaystyle\ce{CoASH+NAD++Pyr}$
$\displaystyle\rightleftharpoons\ce{AcetylCoA+CO2+NADH}$ (81a) and the
Tricarboxylic Acid Cycle reactions $\displaystyle\ce{AcetylCoA+H2O+OAA}$
$\displaystyle\rightleftharpoons\ce{Citrate+CoASH+H+}$ (81b)
$\displaystyle\ce{Citrate}$
$\displaystyle\rightleftharpoons\ce{Aconitate+H2O}$ (81c)
$\displaystyle\ce{Aconitate+H2O}$
$\displaystyle\rightleftharpoons\ce{Isocitrate}$ (81d)
$\displaystyle\ce{Isocitrate+NAD+}$
$\displaystyle\rightleftharpoons\ce{CO2+KG+NADH}$ (81e)
$\displaystyle\ce{CoASH+KG+NAD+}$
$\displaystyle\rightleftharpoons\ce{CO2+NADH+SucCoA}$ (81f)
$\displaystyle\ce{ADP+P_{\mathrm{i}}+SucCoA}$
$\displaystyle\rightleftharpoons\ce{ATP+CoASH+Succ}$ (81g)
$\displaystyle\ce{CoQ+Succ}$
$\displaystyle\rightleftharpoons\ce{CoQH2+Fumarate}$ (81h)
$\displaystyle\ce{Fumarate+H2O}$ $\displaystyle\rightleftharpoons\ce{Malate}$
(81i) $\displaystyle\ce{Malate+NAD+}$
$\displaystyle\rightleftharpoons\ce{H++NADH+OAA}$ (81j)

When chemostatting H+, ATP, ADP, P_i, H2O, CO2, NAD+, NADH, CoQ, CoQH2, Pyr
and when each reaction happens once, we arrive at the single emergent cycle
with the net stoichiometry given in Eq. (50). This cycle represents the
oxidation of pyruvate to CO2 and the simultaneous reduction of ubiquinone
(CoQ) and NAD+.

###  A.4 Electron Transport Chain and ATP Synthesis

The reactions composing the electron transport chain are

$\displaystyle\ce{CoQ+5H++NADH}$
$\displaystyle\rightleftharpoons\ce{CoQH2+4H^{+}_{P}+NAD+}$ (82a)
$\displaystyle\ce{CoQH2+4H++{\textstyle\frac{1}{2}}O2}$
$\displaystyle\rightleftharpoons\ce{CoQ+4H^{+}_{P}+H2O}$ (82b)
$\displaystyle\ce{ADP+4H^{+}_{P}+P_{\mathrm{i}}}$
$\displaystyle\rightleftharpoons\ce{ATP+3H++H2O}$ (82c)

When chemostatting ATP, ADP, P_i, H+, H2O, O2, NAD+, NADH, CoQ, and CoQH2, we
observe two emergent cycles. The first originates from the reactions
$\text{(A6a)}+\text{(A6b)}+2\>\text{(A6c)}$, and has a net stoichiometry given
in Eq. (54). This cycle corresponds to the oxidation of NADH via molecular
oxygen in order to synthesize ATP. The second cycle originates from the
reactions $\text{(A6b)}+\text{(A6c)}$, and has a net stoichiometry given in
Eq. (58). This cycle corresponds to the synthesis of ATP via the oxidation of
CoQH2. As discussed in the main text, when CoQ and CoQH2 are regarded as
internal species, only the first of these two cycles emerges.

###  A.5 Ethanol Fermentation

In ethanol fermentation, the reactions composing glycolysis, Eq. (80a–j), are
complemented by

$\displaystyle\ce{Pyr+H+}$ $\displaystyle\rightleftharpoons\ce{Ethanal+CO2}$
(83a) $\displaystyle\ce{Ethanal+NADH+H+}$
$\displaystyle\rightleftharpoons\ce{Ethanol+NAD+}$ (83b)

When replacing Pyruvate, NAD+, and NADH for Ethanol as chemostats, one finds
the emergent cycle described by the chain of reactions $\text{(A4a--
e)}+2\>\text{(A4f--j)}+2\>\text{(A7a,b)}$ and whose net balance is given by
Eq. (59).

###  A.6 Cellular Respiration

Cellular respiration comprise all reactions of glycolysis (without
gluconeogenesis) (80), the TCA cycle (81), and the electron transport chain
(82). By identifying H2O, H+, ATP, ADP, P_i, Glucose, CO2, O2 as chemostatted
species, cellular respiration can be seen as the chemical process written in
Eq. (40). This balance corresponds to the emergent cycle depicted in figure
10. It is composed of the reactions (A4a–e) \+ 2 (A4f–j) \+ 2 (A5a–j) \+ 10
(A6a) \+ 12 (A6b) \+ 22 (A6c) .

## References

  * Yang  _et al._ [2021] X. Yang, M. Heinemann, J. Howard, G. Huber, S. Iyer-Biswas, G. L. Treut, M. Lynch, K. L. Montooth, D. J. Needleman, S. Pigolotti, J. Rodenfels, P. Ronceray, S. Shankar, I. Tavassoly, S. Thutupalli, D. V. Titov, J. Wang, and P. J. Foster, Physical bioenergetics: Energy fluxes, budgets, and constraints in cells, Proceedings of the National Academy of Sciences 118, e2026786118 (2021).
  * Nelson and Cox [2008] D. Nelson and M. Cox,  _Lehninger Principles of Biochemistry_ (W. H. Freeman, 2008).
  * Voet  _et al._ [2016] D. Voet, J. G. Voet, and C. W. Pratt,  _Fundamentals of biochemistry_ (Wiley, 2016).
  * Bejan [2016] A. Bejan,  _Advanced engineering thermodynamics_ (John Wiley & Sons, 2016).
  * Hill [1977] T. L. Hill,  _Free energy transduction in biology: the steady-state kinetic and thermodynamic formalism_ (Academic Press, New York (NY), 1977).
  * Wachtel  _et al._ [2018] A. Wachtel, R. Rao, and M. Esposito, Thermodynamically consistent coarse graining of biocatalysts beyond Michaelis–Menten, New Journal of Physics 20, 042002 (2018).
  * Jülicher  _et al._ [1997] F. Jülicher, A. Ajdari, and J. Prost, Modeling molecular motors, Rev. Mod. Phys. 69, 1269 (1997).
  * Parmeggiani  _et al._ [1999] A. Parmeggiani, F. Jülicher, A. Ajdari, and J. Prost, Energy transduction of isothermal ratchets: Generic aspects and specific examples close to and far from equilibrium, Phys. Rev. E 60, 2127 (1999).
  * Parrondo and de Cisneros [2002] J. M. R. Parrondo and B. J. de Cisneros, Energetics of brownian motors: a review, Applied Physics A 75, 179 (2002).
  * Seifert [2012] U. Seifert, Stochastic thermodynamics, fluctuation theorems and molecular machines, Reports on Progress in Physics 75, 126001 (2012).
  * Brown and Sivak [2020] A. I. Brown and D. A. Sivak, Theory of nonequilibrium free energy transduction by molecular machines, Chem. Rev. 120, 434 (2020).
  * Amano  _et al._ [2022] S. Amano, M. Esposito, E. Kreidt, D. A. Leigh, E. Penocchio, and B. M. W. Roberts, Insights from an information thermodynamics analysis of a synthetic molecular motor, Nature Chemistry 14, 530 (2022).
  * Schnakenberg [1976] J. Schnakenberg, Network theory of microscopic and macroscopic behavior of master equation systems, Rev. Mod. Phys. 48, 571 (1976).
  * Andrieux and Gaspard [2007] D. Andrieux and P. Gaspard, Fluctuation theorem for currents and schnakenberg network theory, J. Stat. Phys. 127, 107 (2007).
  * Polettini and Esposito [2014] M. Polettini and M. Esposito, Irreversible thermodynamics of open chemical networks. i. emergent cycles and broken conservation laws, Journal of Chemical Physics 141, 024117 (2014).
  * Rao and Esposito [2016] R. Rao and M. Esposito, Nonequilibrium thermodynamics of chemical reaction networks: Wisdom from stochastic thermodynamics, Physical Review X 6, 041064 (2016).
  * Avanzini  _et al._ [2020] F. Avanzini, G. Falasco, and M. Esposito, Thermodynamics of non-elementary chemical reaction networks, New J. Phys. 22, 093040 (2020).
  * Avanzini  _et al._ [2021] F. Avanzini, E. Penocchio, G. Falasco, and M. Esposito, Nonequilibrium thermodynamics of non-ideal chemical reaction networks, J. Chem. Phys. 154, 094114 (2021).
  * Note [1] For isobaric solutions, $\mu_{i}(\bm{z})$ must be regarded as the Gibbs free-energy content of species $i$.
  * Alberty [2003] R. A. Alberty,  _Thermodynamics of biochemical reactions_ (Wiley-Interscience, 2003).
  * Horn and Jackson [1972] F. Horn and R. Jackson, General mass action kinetics, Archive for Rational Mechanics and Analysis 47, 81 (1972).
  * Schuster and Schuster [1989] S. Schuster and R. Schuster, A generalization of wegscheider’s condition. implications for properties of steady states and for quasi-steady-state approximation, Journal of Mathematical Chemistry 3, 25 (1989).
  * Avanzini and Esposito [2022] F. Avanzini and M. Esposito, Thermodynamics of concentration vs flux control in chemical reaction networks, The Journal of Chemical Physics 156, 014116 (2022).
  * Falasco  _et al._ [2018] G. Falasco, R. Rao, and M. Esposito, Information thermodynamics of turing patterns, Phys. Rev. Lett. 121, 108301 (2018).
  * Penocchio  _et al._ [2019] E. Penocchio, R. Rao, and M. Esposito, Thermodynamic efficiency in dissipative chemistry, Nat. Commun. 10, 3865 (2019).
  * Beard  _et al._ [2002] D. A. Beard, S.-d. Liang, and H. Qian, Energy balance for analysis of complex metabolic networks, Biophys. J. 83, 79 (2002).
  * Yang  _et al._ [2005] F. Yang, H. Qian, and D. A. Beard, Ab initio prediction of thermodynamically feasible reaction directions from biochemical network stoichiometry, Metab. Eng. 7, 251 (2005).
  * Palsson [2015] B. Ø. Palsson,  _Systems Biology: Constraint-based Network Reconstruction and Analysis_ (Cambridge University Press, 2015).
  * Niebel  _et al._ [2019] B. Niebel, S. Leupold, and M. Heinemann, An upper limit on gibbs energy dissipation governs cellular metabolism, Nature Metabolism 1, 125 (2019).
  * Saldida  _et al._ [2020] J. Saldida, A. P. Muntoni, D. de Martino, G. Hubmann, B. Niebel, A. M. Schmidt, A. Braunstein, A. Milias-Argeitis, and M. Heinemann, Unbiased metabolic flux inference through combined thermodynamic and 13c flux analysis, bioRxiv 10.1101/2020.06.29.177063 (2020).
  * Oftadeh  _et al._ [2021] O. Oftadeh, P. Salvy, M. Masid, M. Curvat, L. Miskovic, and V. Hatzimanikatis, A genome-scale metabolic model of saccharomyces cerevisiae that integrates expression constraints and reaction thermodynamics, Nature Communications 12, 10.1038/s41467-021-25158-6 (2021).
  * Akbari  _et al._ [2021] A. Akbari, J. T. Yurkovich, D. C. Zielinski, and B. O. Palsson, The quantitative metabolome is shaped by abiotic constraints, Nature Communications 12, 10.1038/s41467-021-23214-9 (2021).
  * Bar-Even  _et al._ [2012] A. Bar-Even, A. Flamholz, E. Noor, and R. Milo, Thermodynamic constraints shape the structure of carbon fixation pathways, Biochim. Biophys. Acta (BBA) - Bioenergetics 1817, 1646 (2012).
  * Noor  _et al._ [2014] E. Noor, A. Bar-Even, A. Flamholz, E. Reznik, W. Liebermeister, and R. Milo, Pathway thermodynamics highlights kinetic obstacles in central metabolism, PLoS Comput Biol 10, 1 (2014).
  * Noor  _et al._ [2013] E. Noor, H. S. Haraldsdóttir, R. Milo, and R. M. T. Fleming, Consistent estimation of gibbs energy using component contributions, PLOS Computational Biology 9, 1 (2013).