Conservation laws and work fluctuation relations in chemical reaction networks

Rao, Riccardo; Esposito, Massimiliano

doi:10.1063/1.5042253

We formulate a nonequilibrium thermodynamic description for open chemical reaction networks (CRNs) described by a chemical master equation. The topological properties of the CRN and its conservation laws are shown to play a crucial role. They are used to decompose the entropy production into a potential change and two work contributions, the first due to time dependent changes in the externally controlled chemostats concentrations and the second due to flows maintained across the system by nonconservative forces. These two works jointly satisfy a Jarzynski and Crooks fluctuation theorem. In the absence of work, the potential is minimized by the dynamics as the system relaxes to equilibrium and its equilibrium value coincides with the maximum entropy principle. A generalized Landauer’s principle also holds: the minimal work needed to create a nonequilibrium state is the relative entropy of that state to its equilibrium value reached in the absence of any work.

I. INTRODUCTION

Nonequilibrium thermodynamic descriptions of stochastic (bio-)chemical processes have long since been developed. Among the first, Hill and co-workers studied bio-catalysts as small fluctuating machines operating at steady-state. They introduced the concept of free energy transduction and analyzed how one form of chemical work can drive another one against its spontaneous direction.¹ The importance of decomposing currents into network cycles (i.e., cyclic sets of transitions) was already emphasized. These results were, however, limited to steady-state systems described by linear chemical reaction networks (CRNs). The stochastic as well as the deterministic dynamics of these CRNs is described by the same linear rate equations for, respectively, probabilities or concentrations. They model, for instance, conformational changes of an enzyme or of a membrane transporter. Inspired by these seminal studies, Schnakenberg formulated steady-state thermodynamics for generic Markov jump processes and provided a systematic cycle decomposition for the entropy production (EP) rate.² He considered, in particular, the stochastic description in terms of the Chemical Master Equation (CME)^3,4 of nonlinear chemical reaction networks, i.e., CRNs described at the deterministic level by nonlinear rate equations for concentrations. The Brussels school and Ross and co-workers focused on the connection between the thermodynamic description resulting from the stochastic and the deterministic dynamics.^5–8

With the advent of stochastic thermodynamics,^9–12 the focus moved to the study of fluctuations, rather than focusing on the first two moments. Gaspard first showed that EP fluctuations in nonlinear CRNs at steady state satisfy a fluctuation theorem (FT).¹³ This result was later expressed in terms of currents along Schnakenberg cycles.^14,15 Fluctuations in complex chemical dynamics such as bistability was analyzed, amongst others, by Qian and co-workers.^16–18 A first formulation of stochastic thermodynamics for CRNs beyond steady state was done by Schmiedl and Seifert.¹⁹

Despite this long history, none of these descriptions made use of the specific topology of the CRN encoded in its stoichiometric matrix. Mathematicians know, however, that the CRN topology plays an important role in its deterministic^20,21 as well as stochastic dynamics.^22,23 But the question of how it affects the thermodynamic description was only studied recently: for deterministic dynamics in Refs. 24 and 25 and for stochastic dynamics at steady state in Ref. 26. In this paper, we address this question in full generality for CRNs whose dynamics is stochastic. We will do so by presenting a formulation of stochastic thermodynamics for CRNs which systematically makes use of the conservation laws. Doing so leads to a significantly more informative thermodynamic description. In particular, we decompose the EP into three fundamental dissipative contributions: a newly defined potential change, a driving work contribution due to time dependent changes in the externally controlled chemostats concentrations, and a nonconservative work contribution due to a minimal set of flows maintained across the system by nonconservative forces. In contrast to the traditional chemical work given by minus the free energy change in the chemostats, these two new work contributions are shown to jointly satisfy a finite-time detailed and integral FT when the CRN is initially prepared in an equilibrium state. In turn, the importance of the potential lies in the fact that it is minimized by the relaxation dynamics towards equilibrium in the absence of the first two work contributions, i.e., when the system is detailed-balanced. It can be seen as a Legendre transform with respect to those conservation laws that are broken by the chemostats. At equilibrium, it coincides with the potential obtained from maximizing entropy with broken conservation laws as constrains. We also discuss the connection of our findings to absolute irreversibility,²⁷ to free energy transduction in nonlinear CRNs, and to cycle decompositions of the entropy production. Finally, we derive a nonequilibrium Landauer’s principle for the driving and nonconservative work which generalizes the previous ones to nondetailed-balanced dynamics.^28,29

A. Outline

The paper is organized as follows. In Sec. II, we review the stochastic description of closed and open CRNs and introduce conservation laws and stoichiometric cycles. In Sec. III, the connection with thermodynamics is made. The stochastic reaction rates are expressed in terms of Gibbs potentials via the equilibrium distribution of the closed CRN. Enthalpy balance and entropy balance are defined along stochastic trajectories, and Jarzynski-like FTs for the chemical work are discussed. In Sec. IV, the EP is partitioned into its three contributions. In Sec. V, we analyze open detailed balanced CRNs; more specifically, their relaxation to equilibrium as chemostats are successively introduced. In Sec. VI, finite-time detailed FTs for the driving and nonconservative work are derived. In Sec. VII, the ensemble average description is presented and the nonequilibrium Landauer’s principle is derived. Finally in Sec. VIII, our results are applied on a simple model to show the importance of our formulation for free energy transduction. Throughout the paper, our formalism is illustrated using a simple enzymatic scheme, whereas some technical derivations are given in Appendixes A and B. We also provide a table which lists the symbols used throughout the paper, Table III.

II. STOCHASTIC DYNAMICS AND CRN TOPOLOGY

A. Chemical reaction networks

We consider a homogeneous, isobaric, and isothermal ideal dilute solution made of N_z chemical species, encoded in a vector z. Their integer-valued population n changes due to internal reactions which we label by {ρ_i} for ρ_i = ±1, …, ±N_i

ν_{ρ_{i}} \cdot z ⇌_{k_{- ρ_{i}}}^{k_{ρ_{i}}} ν_{- ρ_{i}} \cdot z .

(1)

In open CRNs, the population of a subset of species, named exchanged species and denoted by y where z ≡ (x, y), varies also due to exchanges with external chemostats denoted by Y. Their effect is modeled by exchange reactions, {ρ_e} for ρ_e = ±1, …, ±N_y (see Fig. 1)

ν_{ρ_{e}}^{y} \cdot y ⇌_{k_{- ρ_{e}}}^{k_{ρ_{e}}} ν_{- ρ_{e}}^{Y} \cdot Y .

(2)

The non-negative integer-valued vectors ${ν_{ρ} \equiv (ν_{ρ}^{x}, ν_{ρ}^{y})}$ for ρ ∈ {ρ_i} ∪ {ρ_e} encode the stoichiometric coefficients of each reaction. Note that each entry of $ν_{ρ_{e}}^{y}$ and $ν_{ρ_{e}}^{Y}$ is nonzero and equal to one only if it corresponds to the species exchanged by ρ_e. Note also that all reactions are assumed elementary and reversible. For any reaction ρ, −ρ denotes its backward counterpart and the sums over ρ includes both + and −. The different types of species are summarized in Table I.

FIG. 1.

View large Download slide

Pictorial representation of an open CRN modeling an enzymatic scheme discussed in example 1.

TABLE I.

In the second column, the symbols used for the various species are listed. The corresponding total number of entries and symbols used to denote their abundance is given in the third and fourth column, respectively. The first column summarizes the name used to refer to these species, while the last one lists the symbol used to collect the abundances of the internal species. Internal species, x and y, are characterized by low populations, n. The population of x can change only because of reactions, whereas that of y is also exchanged with chemostats, which are identified by Y, Eq. (1).

Species		Symbol	Number	Abundance
Internal	$\{) \begin{matrix} Exchanged \end{matrix}$	x	N_x	n_x	$\}) n$
	$\{) \begin{matrix} Exchanged \end{matrix}$	y	N_y	n_y	$\}) n$
	Chemostatted	Y	N_y	[Y]

The topology of the CRN is encoded in its stoichiometric vectors

S_{ρ} ≔ ν_{- ρ} - ν_{ρ} and S_{ρ}^{Y} ≔ ν_{- ρ}^{Y} - ν_{ρ}^{Y} .

(3)

The former quantifies the change of the population induced by a given reaction ρ, whereas the latter quantifies the corresponding amount of chemostatted species that is exchanged. By definition, S_ρ = −S_−ρ and $S_{ρ}^{Y} = - S_{- ρ}^{Y}$ ⁠. Collecting the column vectors S_ρ (respectively, $S_{ρ}^{Y}$ ⁠) corresponding to arbitrarily-chosen forward reactions defines the internal (respectively, external) stoichiometric matrix denoted by $S$ (respectively, $S^{Y}$ ⁠). It is not difficult to see that these can be decomposed as

S \equiv (\begin{matrix} S_{i} & S_{e} \end{matrix}) \equiv (\begin{matrix} S_{i}^{x} & O \\ S_{i}^{y} & S_{e}^{y} \end{matrix})

(4)

and

S^{Y} \equiv (\begin{matrix} S_{i}^{Y} & S_{e}^{Y} \end{matrix}) \equiv (\begin{matrix} O & - S_{e}^{y} \end{matrix}) .

(5)

In closed CRNs, all exchange reactions disappear and the stoichiometric matrix reduces to $S_{i}$ ⁠.

Remark

Previous studies on thermodynamics of CRNs, e.g., Refs. 19, 24, 25, and 30, describe open CRNs by assuming that the exchanged species y are so abundant that they can be regarded as particle reservoirs within the system. As a result, the exchange reactions are disregarded, y are treated as chemostatted, and the stoichiometric matrices read

S_{a l t} = S_{i}^{x} and S_{a l t}^{Y} = S_{i}^{y} .

(6)

In the closed CRNs, the stoichiometric matrix becomes

(S_{a l t}, S_{a l t}^{Y})^{T}

⁠. As we will see, the two approaches are formally very similar, but the former has the advantage of preserving the number of internal species when the CRN is chemostatted. This makes it more suitable for a stochastic description.

Example 1.

For the open CRN in Fig. 1,

x = (E, E^{*}, E^{* *}), y = (A, B), Y = (A_{e}, B_{e})

(7)

and

n = (n_{E}, n_{E^{*}}, n_{E^{* *}}, n_{A}, n_{B}) .

(8)

Internal reactions, ρ_i = ±1, …, ±4, are distinguished from the exchange ones, ρ_e = ±a, ±b, and the stoichiometric matrices read

S = \begin{matrix} E \\ E * \\ E ** \\ A \\ B \end{matrix} \overset{\begin{matrix} + 1 & + 2 & + 3 & + 4 & + a & + b \end{matrix}}{(\begin{matrix} \begin{matrix} - 1 & 1 & - 1 & 1 \\ 1 & - 1 & 0 & 0 \\ 0 & 0 & 1 & - 1 \\ - 1 & 0 & 0 & 1 \\ 0 & 1 & - 1 & 0 \end{matrix} & \begin{matrix} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 1 & 0 \\ 0 & 1 \end{matrix} \end{matrix})}

(9)

and

S^{Y} = \begin{matrix} A_{e} \\ B_{e} \end{matrix} \overset{\begin{matrix} 1 + & 2 + & 3 + & 4 + & a + & b + \end{matrix}}{(\begin{matrix} \begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{matrix} & \begin{matrix} 1 - & 0 \\ 0 & 1 - \end{matrix} \end{matrix})}

(10)

for our arbitrary choice of forward reactions.

Notation

Henceforth, we will use the following notation:

a! = \prod_{i} a_{i}!, a^{\cdot b} = \prod_{i} a_{i}^{b_{i}} and c^{\cdot b} = c^{\sum_{i} b_{i}}

for generic vectors a and b and for a generic constant c. “ln a” must be read as a vector whose entries are the logarithm of the entries of a. 1 denotes a vector whose entries are all equal to 1. Total and partial time derivatives are written as d_t and ∂_t, and the overdot “·” denotes the rates of change of observables which are not state functions.

B. Chemical master equation

In our stochastic description, n is treated as a fluctuating variable and all reactions are regarded as stochastic events. The probability of finding the CRN in the state n at time t is denoted by p_n ≡ p_n(t) and its evolution is ruled by the CME^3,4,31

\begin{matrix} d_{t} p_{n} & = & \sum_{ρ} \{) w_{- ρ} (n + S_{ρ}) p_{n + S_{ρ}} - w_{ρ} (n) p_{n} \}) \\ = & \sum_{m} W_{n m} p_{m}, \end{matrix}

(11)

where the stochastic generator reads

W_{n m} = \sum_{ρ} w_{ρ} (m) \{δ_{n, m + S_{ρ}} - δ_{n, m}\} .

(12)

Since all reactions are assumed elementary, we consider mass-action stochastic reaction rates

w_{ρ} (n) ≔ k_{ρ} \frac{V}{V^{\cdot ν_{ρ}}} {[Y]}^{\cdot ν_{ρ}^{Y}} \frac{n!}{(n - ν_{ρ})!},

(13)

where {k_ρ} denote the rate constants. The dependence on the volume V ensures the correct scaling when taking the large particle limit and guarantees that {k_ρ} are the same as in deterministic descriptions.³² The chemostat concentrations [Y] only appear in exchange reactions ρ_e and quantify the concentration of the exchanged species in the chemostats. Hence, they are real-valued, nonfluctuating, and unaffected by the occurrence of exchange reactions. We assume that [Y] can change over time and their value at each time t is encoded in the driving protocol π_t. This may describe, for instance, the controlled injection of certain molecules across a cell membrane. In such situations, the CRN is said to be subjected to a “driving.” In the absence of driving, the CRNs is instead said to be autonomous.

Equilibrium probability distributions are of crucial importance for our discussion. They satisfy the detailed balance property

w_{ρ} (n) p_{n}^{e q} = w_{- ρ} (n + S_{ρ}) p_{n + S_{ρ}}^{e q}, for all ρ, n .

(14)

This means that the probability current of any reaction ρ occurring from any state n vanishes. Stochastic CRNs which admit a steady-state probability distribution satisfying Eq. (14) are referred to as detailed balanced. Their stochastic thermodynamics will be analyzed in Sec. V.

Example 2.

For the CRN in Fig. 1, the transition rates are

\begin{aligned} w_{+ 1} & = k_{+ 1} n_{A} n_{E}, & w_{- 1} & = k_{- 1} n_{E^{*}}, \\ w_{+ 2} & = k_{+ 2} n_{E^{*}}, & w_{- 2} & = k_{- 2} n_{E} n_{B}, \\ w_{+ 3} & = k_{+ 3} n_{E} n_{B}, & w_{- 3} & = k_{- 3} n_{E^{* *}}, \\ w_{+ 4} & = k_{+ 4} n_{E^{* *}}, & w_{- 4} & = k_{- 4} n_{E} n_{A}, \\ w_{+ a} & = k_{+ a} [A_{e}], & w_{- a} & = k_{- a} n_{A}, \\ w_{+ b} & = k_{+ b} [B_{e}], & w_{- b} & = k_{- b} n_{B} . \end{aligned}

(15)

C. Stochastic trajectories

A stochastic trajectory of duration t, $n_{t}$ ⁠, is defined as a set of reactions {ρ_l} sequentially occurring at times {t_l} starting from n₀ at time t₀. Such trajectories can be generated by a stochastic simulation algorithm.³³ Given the initial state, a trajectory is completely characterized by

j_{ρ} (n, τ) ≔ \sum_{l} δ_{ρ ρ_{l}} δ_{n n_{t_{l}}} δ (τ - t_{l}),

(16)

which encodes the reactions that occur ({ρ_l}), the states from which these occur (⁠ ${n_{t_{l}}}$ ⁠), and the reaction times ({t_l}). The transition index l runs from l = 1 to the last transition prior to time t, N_t. The instantaneous reaction fluxes

J_{ρ} (τ) ≔ \sum_{n} j_{ρ} (n, τ) = \sum_{l} δ_{ρ ρ_{l}} δ (τ - t_{l})

(17)

quantify the instantaneous rate of occurrence of each reaction irrespective of the state from which it occurs. Additionally, we denote the population of the CRN at time τ ∈ [t₀ = 0, t] by n_τ.

The path probability of a trajectory reads

\begin{matrix} P [n_{t}] & = & \prod_{l = 0}^{N_{t}} \exp \{- \int_{t_{l}}^{t_{l + 1}} d τ \sum_{ρ} w_{ρ} (n_{τ}, τ)\} \prod_{l = 1}^{N_{t}} w_{ρ_{l}} (n_{t_{l}}, t_{l}), \end{matrix}

(18)

where $t_{N_{t} + 1} ≔ t$ is the final time of the trajectory. The first term accounts for the probability that the system spends {t_l+1 − t_l} time in the state ${n_{t_{l}}}$ ⁠, while the second term accounts for the probability of transitioning. When averaging Eq. (16) over all stochastic trajectories, we obtain the transition rates, Eq. (13),

⟨ j_{ρ} (n, τ) ⟩ = w_{ρ} (n, τ) p_{n} (τ) .

(19)

Changes of generic observables along trajectories are written as

δ X [n_{t}] = \int_{0}^{t} d τ \{) \dot{X} (n_{τ}, τ) + \sum_{n, ρ} δ X_{ρ} (n, τ) j_{ρ} (n, τ) \}),

(20)

where $\dot{X} (n, π_{τ})$ denotes its change in time while the CRN dwells in the state n (it need not be an exact time derivative) and $δ X_{ρ} (n, π_{τ})$ denotes its finite change along the reaction ρ occurring while in n. By contrast, the changes of state observables $O (n, t)$ can be written as

\begin{matrix} Δ O [n_{t}] & = & O (n_{t}, t) - O (n_{0}, 0) \\ = & \int_{0}^{t} d τ \{) {([\partial_{τ} O (n, τ)]|}_{n_{τ}} + \sum_{n, ρ} Δ_{ρ} O (n, τ) j_{ρ} (n, τ) \}), \end{matrix}

(21)

where $\partial_{τ} O (n, τ)$ is the time derivative of $O (n, τ)$ and

Δ_{ρ} O (n, τ) ≔ O (n + S_{ρ}, τ) - O (n, τ)

(22)

is the difference of $O (n, τ)$ along reactions; see Fig. 2.

FIG. 2.

View large Download slide

Pictorial representation of the change of a state variable observable $O$ along a trajectory. The orange dashed curves represent the changes due to the protocol—the first term in Eq. (21)—while the vertical blue lines represent changes due to reactions—the second term in Eq. (21).

D. Conservation laws

The topological properties of CRNs are encoded in the matrices $S$ and $S^{Y}$ and can be identified via their cokernels and kernels. Conservation laws ℓ are defined as vectors in $coker S$

ℓ \cdot S_{ρ} = 0, for all ρ .

(23)

They identify conserved quantities, called components³⁴

L_{n} ≔ ℓ \cdot n .

(24)

Despite the fact that L_n depends on the stochastic variable n, the probability of observing any specific value L,

P (L) ≔ \sum_{n} p_{n} δ [L_{n}, L],

(25)

is constant over time, i.e., d_tP(L) = 0. δ is a Kronecker delta. More generally, any observable of type $O (L_{n})$ does not fluctuate

d_{t} \sum_{n} p_{n} O (L_{n}) = 0,

(26)

as a direct consequence of the fact that $Δ_{ρ} O (L_{n}) = 0$ ⁠. Clearly, P(L) can be deduced from the initial conditions p_n(0) and only those states for which P(L_n, 0) is nonvanishing have a finite probability of being observed during the subsequent stochastic dynamics.

In closed CRNs, conservation laws (23) follow from

ℓ^{x} \cdot S_{ρ_{i}}^{x} + ℓ^{y} \cdot S_{ρ_{i}}^{y} = 0, for all ρ_{i} .

(27)

We denote a set of linearly independent conservation laws of the closed CRN by $\{ℓ_{λ}\}$ ⁠, and the corresponding components by $\{L_{n}^{λ} ≔ ℓ_{λ} \cdot n\}$ ⁠, for $λ = 1, \dots, N_{λ} ≔ \dim coker S_{i}$ ⁠. The choice of this set is not unique, and different choices have different physical meanings. This set is never empty since the total mass is always conserved. The latter corresponds to a ℓ whose entries are the masses of each species. Physically, the conservation laws of closed CRNs can always be chosen so as to correspond to moieties, which are parts of molecules exchanged between species along reactions or subject to isomerization.³⁵

For open CRNs, the condition identifying conservation laws, Eq. (23), becomes

ℓ^{x} \cdot S_{ρ_{i}}^{x} + ℓ^{y} \cdot S_{ρ_{i}}^{y} = 0, for all ρ_{i},

(28a)

ℓ^{y} \cdot S_{ρ_{e}}^{y} = 0, for all ρ_{e} .

(28b)

We now recall that for all ρ_e, there is one and only one exchanged species for which the corresponding entry of $S_{ρ_{e}}^{y}$ is different from zero. Hence, Eq. (28b) demands that $ℓ^{y}$ = 0 and Eq. (28) become $ℓ^{x} \cdot S_{ρ_{i}}^{x} = 0$ for all ρ_i.

Crucially, any set of independent conservation laws of the open CRN, Eq. (28), denoted by ${ℓ_{λ_{u}}}$ ⁠, for $λ_{u} = 1, \dots, N_{λ_{u}} ≔ \dim coker S < N_{λ}$ ⁠, can be regarded as a subset of the conservation laws of the closed CRN, ${ℓ_{λ}} \equiv {ℓ_{λ_{u}}} \cup {ℓ_{λ_{b}}}$ ⁠, since they satisfy Eq. (27), too. In view of this, we call them unbroken conservation laws. The remaining independent conservation laws, labeled as ${ℓ_{λ_{b}}}$ and referred to as broken, satisfy Eq. (27) while not Eq. (28). They involve exchanged species, $ℓ_{λ_{b}}^{y} \neq 0$ ⁠; hence, $ℓ_{λ_{b}}^{y} \cdot S_{ρ_{e}}^{y} \neq 0$ and the probability distribution of any set ${L_{n}^{λ_{b}} \equiv ℓ_{λ_{b}} \cdot n}$ ⁠,

P ({L_{λ_{b}}}) ≔ \sum_{n} p_{n} \prod_{λ_{b}} δ [) L_{n}^{λ_{b}}, L_{λ_{b}}]),

(29)

changes in time.

Summarizing, in open CRNs, the chemostatting breaks a subset of the conservation laws of the corresponding closed CRN, ${ℓ_{λ_{b}}}$ ⁠. Only the probability distribution of the unbroken components ${L_{n}^{λ_{u}} \equiv ℓ_{λ_{u}} \cdot n}$ ⁠,

P ({L_{λ_{u}}}) ≔ \sum_{n} p_{n} \prod_{λ_{u}} δ [) L_{n}^{λ_{u}}, L_{λ_{u}}]),

(30)

is invariant and completely determined by the initial probability distribution p_n(0). The state space identified by one particular set of values for ${L_{λ_{u}}}$ is called stoichiometric compatibility class.

Example 3.

The CRN in Fig. 1 has two conservation laws

ℓ_{E} = \begin{matrix} \begin{matrix} E & E^{*} & E^{**} & A & B \end{matrix} \\ (\begin{matrix} 1 & 1 & 1 & 0 & 0 \end{matrix}) \end{matrix},

(31a)

ℓ_{b} = \begin{matrix} \begin{matrix} E & E^{*} & E^{**} & A & B \end{matrix} \\ (\begin{matrix} 0 & 1 & 1 & 1 & 1 \end{matrix}) \end{matrix},

(31b)

among which the second is broken. The unbroken conservation law identifies the enzyme moiety and corresponds to the total number of enzyme molecules populating the CRN,

L_{n}^{E} = n_{E} + n_{E^{*}} + n_{E^{* *}}

⁠. Instead, the broken one identifies the moiety A—or equivalently B—

L_{n}^{b} = n_{E^{*}} + n_{E^{* *}} + n_{A} + n_{B}

⁠.

E. Stoichiometric cycles

We can now set the stage for the thermodynamic description based on a stoichiometric cycle decomposition. This section, as well as the other ones discussing cycles (Secs. IV D, VI A, and VII C), may be omitted at a first reading.

Additional information about the CRN topology is provided by the stoichiometric cycles c = {c_ρ} as they are vectors in $\ker S$ ⁠. Equivalently, these satisfy

\sum_{ρ} S_{ρ} c_{ρ} = 0

(32)

and at most one entry for each forward–backward transition pair is nonzero. Since $S$ is integer-valued, any c can always be chosen non-negative-integer-valued. In this way, its entries denote the number of times each transition occurs along a transformation which overall leaves the state n unchanged. Alternatively, a stoichiometric cycle can be seen as a set of reactions ${ρ_{c 1}, ρ_{c 2}, \dots, ρ_{c N_{c}}}$ identifying a closed loop in the state space

n \to n + S_{ρ_{c 1}} \to \dots \to n + \sum_{i = 1}^{N_{c}} S_{ρ_{c i}} = n,

(33)

where $\sum_{i = 1}^{N_{c}} S_{ρ_{c i}} = \sum_{ρ} S_{ρ} c_{ρ} = 0$ ⁠.

We now relate cycles of the closed and open CRNs as previously done for conservation laws. In the closed CRN, the stoichiometric cycles are given by

\sum_{ρ_{i}} S_{ρ_{i}}^{x} c_{ρ_{i}} = 0,

(34a)

\sum_{ρ_{i}} S_{ρ_{i}}^{y} c_{ρ_{i}} = 0 .

(34b)

The entries corresponding to the exchange reactions are taken equal to 0: $c_{ρ_{e}} = 0$ ⁠, for all ρ_e. Let us denote by {c^α}, for $α = 1, \dots, N_{α} ≔ \dim \ker S_{i}$ ⁠, a set of independent stoichiometric cycles of the closed CRN.

In the open CRN, the condition identifying cycles, Eq. (32), reads

\sum_{ρ_{i}} S_{ρ_{i}}^{x} c_{ρ_{i}} = 0,

(35a)

\sum_{ρ_{i}} S_{ρ_{i}}^{y} c_{ρ_{i}} + \sum_{ρ_{e}} S_{ρ_{e}}^{y} c_{ρ_{e}} = 0 .

(35b)

Since the cycles of the closed CRN satisfy Eq. (35), they can be regarded as a subset of an independent set of cycles for the open CRN, {c^α, c^η}. We refer to the additional cycles {c^η}, for $η = 1, \dots, N_{η} ≔ \dim \ker S - \dim \ker S_{i}$ ⁠, as emergent. They are characterized by at least one nonzero entry for {ρ_e}, and the vectors

C_{η}^{Y} ≔ \sum_{ρ} () - S_{ρ}^{Y})) c_{ρ}^{η} = \sum_{ρ_{e}} S_{ρ_{e}}^{y} c_{ρ_{e}}^{η} \neq 0

(36)

quantify the amount of exchanged species flowing in the system from the corresponding chemostats upon completion of c^η. As the concentrations of the chemostats are unaffected by the exchange of particles with the system, the emergent stoichiometric cycles can be thought of as pathways transferring chemicals across chemostats while leaving the internal state of the CRN unchanged.

As first proved in Ref. 24, by applying the rank-nullity theorem to the stoichiometric matrices of the open and closed CRNs, one can show that

N_{y} = N_{λ_{b}} + N_{η} .

(37)

In words, for any exchanged species, either a conservation law is broken or an emergent cycle is created.

Example 4.

The CRN in Fig. 1 has one cycle

c_{int} = \begin{matrix} \begin{matrix} + 1 & + 2 & + 3 & + 4 & + a & + b \end{matrix} \\ (\begin{matrix} 1 & 1 & 1 & 1 & 0 & 0 \end{matrix}) \end{matrix}

(38)

and one emergent cycle

c_{int} = \begin{matrix} \begin{matrix} + 1 & + 2 & + 3 & + 4 & + a & + b \end{matrix} \\ (\begin{matrix} 1 & 1 & 0 & 0 & 1 & - 1 \end{matrix}) \end{matrix}

(39)

Negative entries must be interpreted as reactions occurring in the backward direction. The latter cycle corresponds to the injection of one molecule of A, its conversion into one of B passing via E^*, and its ejection

c_{ext} = \overset{\begin{matrix} A & B \end{matrix}}{(\begin{matrix} 1 & - 1 \end{matrix})}

(40)

We can also check the validity of Eq. (37), as the number of chemostats, 2, equals the number of broken conservation laws, 1, see example 3, plus the number of emergent cycles, 1, Eq. (39).

Remark

Stoichiometric cycles must be distinguished from graph-theoretic cycles, also called loops; see, e.g., Ref. 2. To elucidate this point, we note that the network of transitions of a CRN can be regarded as a semi-infinite graph whose vertices are the accessible states n and whose directed edges are given by the reactions—which are encoded in the stoichiometric matrix $S$ ⁠. Hence, one can see that loops are the recursive appearance of stoichiometric cycles, as in Eq. (33). However, they may not be complete at the boundaries of the graph (low n) due to peculiar topological properties of the CRN; see, e.g., Ref. 26. These observations will be used later to relate different approaches for cycle decomposition of thermodynamic quantities.

III. STOCHASTIC THERMODYNAMICS

We now build a nonequilibrium thermodynamic description on top of the stochastic dynamics. We assume that the solvent acts as a thermal reservoir by keeping the temperature, T, and the pressure constant everywhere. Since particle numbers are low, we can assume that that the time scale in which molecules spatially homogenize is much faster than that of reactions. Therefore, if all reactions could be instantaneously shut down, we would observe an equilibrium mixture of inert species at all times. However, due to reactions, the populations of species and their probability distribution can be far from equilibrium. These hypotheses can be regarded as a special case of local equilibrium³⁶ since temperature, pressure, and density are not only locally well defined, but also constant.

A. Equilibrium of closed CRNs

Equilibrium statistical mechanics requires that the equilibrium distribution of a closed CRN with given values of {L_λ} reads

p^{e q} (n | {L_{λ}}) = \frac{\exp \{- β g_{n}\}}{Z ({L_{λ}})} \prod_{λ} δ [) L_{n}^{λ}, L_{λ}]),

(41)

where

g_{n} = (μ ° - 1 k_{B} T \ln n_{s}) \cdot n + k_{B} T \ln n!

(42)

is the Gibbs free energy of the state n derived in Appendix A. The first term quantifies the energetic contribution of each single molecule: μ° ≡ μ°(T) is the vector of standard-state chemical potentials (see Appendix A), whereas −1k_BT ln n_s is an entropic contribution—constant for all species—since n_s is the population of the solvent. The last term is purely entropic and accounts for the indistinguishability of molecules of the same species. In Eq. (41),

Z ({L_{λ}}) = \sum_{m} \exp \{- β g_{m}\} \prod_{λ} δ [) L_{m}^{λ}, L_{λ}])

(43)

is the partition function, while β = 1/(k_BT). When taking into account an ensemble of components, P({L_λ}), Eq. (41) allows us to write

\begin{matrix} p_{n}^{e q} & = & \sum_{{L_{λ}}} p^{e q} (n | {L_{λ}}) P ({L_{λ}}) \\ = & p^{e q} (n | {L_{n}^{λ}}) P ({L_{n}^{λ}}), \end{matrix}

(44)

which can be regarded as a constrained equilibrium distribution. Hence, $p^{e q} (n | {L_{n}^{λ}})$ is the conditional probability of observing n given the stoichiometric compatibility class it identifies.

Equation (44) can also be written as

p_{n}^{e q} = \exp \{- β [) g_{n} - G_{e q} ({L_{n}^{λ}})])\}

(45)

in terms of the equilibrium Gibbs potential of the CRN

G_{e q} ({L_{λ}}) = k_{B} T \ln P ({L_{λ}}) - k_{B} T \ln Z ({L_{λ}}) .

(46)

It is worth emphasizing that G_eq({L_λ}) is a function solely of the set of components and that $G_{e q} ({L_{n}^{λ}})$ needs to be understood as G_eq evaluated in ${L_{n}^{λ}}$ ⁠. Invoking the hypothesis of local equilibrium, we extend G_eq to arbitrary probability distributions p_n,

G (n) ≔ k_{B} T \ln p_{n} + g_{n},

(47)

and we call it stochastic Gibbs potential, as it is the far-from-equilibrium fluctuating expression of G_eq. In addition to the Gibbs free energy of the state n, g_n, it accounts for the entropic contribution due to the uncertainty of p_n: k_BT ln p_n can indeed be written as −T(−k_B ln p_n), where the term in parentheses is the self-information measured in k_B units.³⁷ For closed CRNs at equilibrium, using Eq. (44), G(n) reduces to G_eq in Eq. (46). Also, its average value, the nonequilibrium Gibbs potential

⟨ G ⟩ = \sum_{n} p_{n} [k_{B} T \ln p_{n} + g_{n}]

(48)

takes its minimum value at equilibrium

\begin{matrix} ⟨G⟩ - {⟨ G_{e q} ⟩}_{L} & = & ⟨G - G_{e q}⟩ \\ = & k_{B} T \sum_{n} p_{n} \ln \frac{p_{n}}{p_{n}^{e q}} \\ \equiv & k_{B} T D (p ‖ p^{e q}) \geq 0 . \end{matrix}

(49)

In the first equality, we used the fact that the equilibrium Gibbs potential depends only on the components

\begin{matrix} {⟨ G_{e q} ⟩}_{L} & \equiv & \sum_{{L_{λ}}} P ({L_{λ}}) G_{e q} ({L_{λ}}) \\ = & \sum_{{L_{λ}}} [\sum_{n} p_{n} \prod_{λ} δ [) L_{n}^{λ}, L_{λ}])] G_{e q} ({L_{λ}}) \\ = & \sum_{n} p_{n} G_{e q} ({L_{n}^{λ}}) . \end{matrix}

(50)

In the last equality of Eq. (49), $D (p ‖ p^{e q})$ is the relative entropy of the transient probability distribution p_n with respect to the equilibrium one $p_{n}^{e q}$ ⁠. It is always positive and vanishes only when $p_{n} = p_{n}^{e q}$ ⁠. We will see later (Sec. VII) that Eq. (49) quantifies exactly the average dissipation of the relaxation to equilibrium.

B. Local detailed balance

The zero-th law of thermodynamics for CRNs requires that closed CRNs relax to equilibrium. To ensure this, the dynamical requirement for detailed balance, Eq. (14), is combined with the equilibrium distribution, Eq. (44). As a result, the local detailed balance ensues

\ln \frac{w_{ρ_{i}} (n)}{w_{- ρ_{i}} (n + S_{ρ_{i}})} = - β Δ_{ρ_{i}} g_{n},

(51)

where $Δ_{ρ_{i}} \cdot$ is defined as in Eq. (22). In agreement with deterministic descriptions, see, e.g., Ref. 25, we recover the relation between the rate constants and the standard-state chemical potentials

\ln \frac{k_{ρ_{i}}}{k_{- ρ_{i}}} = - β (μ ° - k_{B} T 1 \ln [s]) \cdot S_{ρ_{i}},

(52)

in which [s] ≔ n_s/V denotes the concentration of the solvent. The local detailed balance (51) should be regarded as a fundamental property of the stochastic reaction rates of elementary reactions valid beyond closed CRNs. This central concept is well known in stochastic thermodynamics because it provides the connection between stochastic dynamics and nonequilibrium thermodynamics.

In open CRNs,

\ln \frac{w_{ρ} (n)}{w_{- ρ} (n + S_{ρ})} = - β () Δ_{ρ} g_{n} + μ_{Y} \cdot S_{ρ}^{Y}))

(53)

generalizes Eq. (51), where

μ_{Y} = μ_{Y} ° + k_{B} T \ln \{[Y] / [s]\}

(54)

are the chemical potentials of the chemostats. The first contribution accounts for the Gibbs free energy change of the internal species, while the second one accounts for the Gibbs free energy exchanged with the chemostats.

We introduce the transition affinities which quantify the force acting along each transition

A_{ρ} (n) = k_{B} T \ln \frac{w_{ρ} (n) p_{n}}{w_{- ρ} (n + S_{ρ}) p_{n + S_{ρ}}} .

(55)

They measure the distance from detailed balance (14), where they all vanish. Using Eq. (53), they can be rewritten in terms of differences of stochastic Gibbs potential (47) as

A_{ρ} (n) = - Δ_{ρ} G (n) + μ_{Y} \cdot () - S_{ρ}^{Y})) .

(56)

This fundamental relation reveals the thermodynamic nature of the dynamical forces acting along reaction. Its early formulation for deterministic chemical kinetics is due to de Donder.³⁸

We will prove in Sec. VII that our theoretical framework based on Eq. (53) guarantees that closed CRNs described by a CME (11) relax to equilibrium, Eq. (44): the average potential ⟨G⟩ is minimized by the dynamics during the relaxation and hence plays the role of a Lyapunov function.

C. Enthalpy and entropy balance

Starting from the stochastic Gibbs potential (47) and the local detailed balance (53), we now formulate the energy and entropy balance along stochastic trajectories.

The stochastic entropy of the CRNs follows from the derivative of the stochastic Gibbs potential (47) with respect to the temperature

S (n) = - {(\frac{\partial G}{\partial T})}_{n} = - k_{B} \ln p_{n} + s_{n} .

(57)

Similar to G(n), S(n) is the far-from-equilibrium fluctuating expression of the entropy.³⁹ The first term on the rhs is the self-information, while the second is the entropy of the state n

s_{n} = - \frac{\partial g_{n}}{\partial T} = (s ° + k_{B} \ln n_{s}) \cdot n - k_{B} \ln n! .

(58)

It accounts for both the entropic contribution carried by each species, i.e., the standard entropies of formation

s ° = - \frac{\partial μ °}{\partial T},

(59)

and the entropic contribution due to the multiplicity of indistinguishable states. When averaged, we recover the Gibbs–Shannon entropy plus an internal entropy contribution

⟨ S ⟩ = \sum_{n} p_{n} [- k_{B} \ln p_{n} + s_{n}] .

(60)

The enthalpy follows from

H (n) = G (n) + T S (n) = g_{n} + T s_{n} = h \cdot n,

(61)

where

h = μ ° + T s ° = h °

(62)

denotes the vector of standard enthalpies of formation, in agreement with traditional thermodynamics of ideal dilute solutions.³⁴ Likewise, the chemical potentials of the chemostats, Eq. (54), will be decomposed in terms of enthalpic and entropic contributions

μ_{Y} = h_{Y} - T s_{Y},

(63)

where h_Y = $h_{Y} °$ and $s_{Y} = s_{Y} ° - k_{B} \ln \{[Y] / [s]\}$ ⁠.

To recover the enthalpy balance along stochastic trajectories, we write the change of enthalpy as the sum of its changes due to reactions

\begin{aligned} Δ H [n_{t}] & = H (n_{t}) - H (n_{0}) \\ = \int_{0}^{t} d τ \sum_{n, ρ} Δ_{ρ} H (n) j_{ρ} (n, τ), \end{aligned}

(64)

where

Δ_{ρ} H (n) = h \cdot S_{ρ} = \underset{= : Q_{ρ}}{\underset{︸}{\underset{= : Q_{ρ}^{t h r}}{\underset{︸}{h \cdot S_{ρ} + h_{Y} \cdot S_{ρ}^{Y}}} + \underset{= : Q_{ρ}^{c h m}}{\underset{︸}{T s_{Y} \cdot () - S_{ρ}^{Y}))}}}} + \underset{= : W_{ρ}^{c}}{\underset{︸}{μ_{Y} \cdot () - S_{ρ}^{Y}))}}, for all n .

(65)

We used Eqs. (21), (62), and (63). The first two contributions, $Q_{ρ}^{t h r}$ ⁠, account for the heat of reaction, i.e., the heat flowing from the thermal reservoir (the solvent). The third term characterizes the heat flowing from the chemostats, $Q_{ρ}^{c h m}$ ⁠. The first three terms, Q_ρ, integrated along the trajectory quantify the total heat flow

Q [n_{t}] = \int_{0}^{t} d τ \{\sum_{ρ} Q_{ρ}^{t h r} J_{ρ} (τ) + T s_{Y} (τ) \cdot I^{Y} (τ)\},

(66)

where the instantaneous external currents

I^{Y} (τ) ≔ \sum_{ρ} () - S_{ρ}^{Y})) J_{ρ} (τ)

(67)

give the amount of exchanged species injected in the CRN at each time; see Eq. (17).

The last term in Eq. (65), $W_{ρ}^{c}$ ⁠, quantifies the Gibbs free energy exchanged with the chemostats. Once integrated, it gives the chemical work

W_{c} [n_{t}] = \int_{0}^{t} d τ μ_{Y} (τ) \cdot I^{Y} (τ) .

(68)

From Eqs. (64)–(68), the enthalpy balance along a trajectory follows

Δ H [n_{t}] = Q [n_{t}] + W_{c} [n_{t}] .

(69)

This is the expression of the first law of thermodynamics for stochastic CRNs at the trajectory level (cf. Ref. 40, Eq. 2.10).

To recover the entropy balance along stochastic trajectories, we notice that since the entropy is a state function, its change along a trajectory reads

\begin{matrix} Δ S [n_{t}] = \int_{0}^{t} d τ \{) {([- \partial_{τ} k_{B} \ln p_{n} (τ)]|}_{n_{τ}} + \sum_{n, ρ} Δ_{ρ} S (n) j_{ρ} (n, τ) \}), \end{matrix}

(70)

as seen in Eq. (21). The changes along transitions can be recast into

\begin{matrix} T Δ_{ρ} S (n) & = & T Δ_{ρ} s_{n} - k_{B} T \ln \frac{p_{n + S_{ρ}}}{p_{n}} \\ = & \underset{= Q_{ρ}}{\underset{︸}{h \cdot S_{ρ} + h_{Y} \cdot S_{ρ}^{Y} + T s_{Y} \cdot () - S_{ρ}^{Y}))}} \underset{= A_{ρ} (n)}{\underset{︸}{- \underset{= Δ_{ρ} G (n)}{\underset{︸}{[Δ_{ρ} g_{n} + k_{B} T \ln \frac{p_{n + S_{ρ}}}{p_{n}}]}} + \underset{= W_{ρ}^{c}}{\underset{︸}{μ_{Y} \cdot () - S_{ρ}^{Y}))}}}}, \end{matrix}

(71)

where we have used Eq. (61). As highlighted with underbraces, the first three terms are the heat flow along reactions, while the last three terms correspond to the affinity of transition, Eq. (56). When integrating over the whole trajectory, we recover the entropy balance

Δ S [n_{t}] = \frac{1}{T} Q [n_{t}] + Σ [n_{t}],

(72)

where the EP (times the temperature) reads

T Σ [n_{t}] = \int_{0}^{t} d τ \{) {([- \partial_{τ} k_{B} T \ln p_{n} (τ)]|}_{n_{τ}} + \sum_{n, ρ} A_{ρ} (n, τ) j_{ρ} (n, τ) \})

(73a)

\begin{align} = k_{B} T \ln \frac{p_{n_{0}} (0)}{p_{n_{t}} (t)} + \int_{0}^{t} d τ j_{ρ} (n, τ) k_{B} T \ln \frac{w_{ρ} (n, τ)}{w_{- ρ} (n + S_{ρ}, τ)} \end{align}

(73b)

= W_{c} [n_{t}] - Δ G [n_{t}] .

(73c)

The second equality follows from the definition of affinity, Eq. (55), when integrating the changes of the probability distribution. Instead, the third one readily follows from the relationship between affinity and Gibbs potential, Eq. (56). It expresses the overall energy dissipated as the difference between the Gibbs free energy supplied by the chemostats and that changing internally.

Mindful of Eq. (18), the EP can be rewritten as the ratio of the probability of observing the trajectory $n_{t}$ under a forward dynamics driven by a protocol π_t over the probability of observing the backward trajectory $n_{t}^{†}$ under a dynamics driven by the time-reversed protocol π^† such that $π_{τ}^{†} ≔ π_{t - τ}$

T Σ [n_{t}] = k_{B} T \ln \frac{p_{n_{0}} (0) P [n_{t}; π]}{p_{n_{t}} (t) P [n_{t}^{†}; π^{†}]} .

(74)

This central result in stochastic thermodynamics^11,39 was formulated for CRNs in Ref. 19 and clearly shows that the EP measures the statistical asymmetry of a trajectory under time reversal. It implies that the EP satisfies the following integral FT:

⟨\exp \{- Σ / k_{B}\}⟩ = 1,

(75)

where the ensemble average $⟨\cdot⟩$ runs over all trajectories. It represents a refinement of the second law of thermodynamics at the trajectory level. Using Jensen’s inequality, the second law ensues, $⟨Σ⟩ \geq 0$ ⁠.

Remark

Using Eqs. and , the local detailed balance, Eq. , can be rewritten as

k_{B} \ln \frac{w_{ρ} (n)}{w_{- ρ} (n + S_{ρ})} = - \frac{1}{T} Q_{ρ}^{t h r} + s_{Y} \cdot S_{ρ}^{Y} + Δ_{ρ} s_{n} .

(76)

The first term is the entropy change in the thermal bath, the second one is the entropy change in the chemostats, whereas the last one is the internal entropy change of the CRN.

Remark

Chemical work and Gibbs potential are defined up to a gauge, which accounts for the choice of the standard-state chemical potentials. Indeed, let us consider the following transformation:

\begin{aligned} μ ° & \to μ ° + \sum_{λ} a_{λ} ℓ_{λ}, \\ {μ °}_{Y} & \to {μ °}_{Y} + \sum_{λ} a_{λ} ℓ_{λ}^{y}, \end{aligned}

(77)

where the second term is a linear combination of conservation laws. This transformation leaves affinities and EP unchanged, while transforming both the chemical work and the Gibbs potential . The former changes as

W_{c} [n_{t}] \to W_{c} [n_{t}] + \sum_{λ_{b}} a_{λ_{b}} ℓ_{λ_{b}}^{y} \cdot I^{Y} [n_{t}],

(78)

where

I^{Y} [n_{t}] = \int_{0}^{t} d τ I^{Y} (τ)

(79)

are the integrated currents of exchanged species flowing in the system. Likewise, the Gibbs potential becomes

G (n) \to G (n) + \sum_{λ} a_{λ} L_{n}^{λ} .

(80)

Using the properties of conservation laws, Sec. II D, it is easy to verify that

Δ L_{λ_{u}} [n_{t}] = 0, Δ L_{λ_{b}} [n_{t}] = ℓ_{λ_{b}}^{y} \cdot I^{Y} [n_{t}],

(81)

which confirms that the gauge terms cancel in the EP, Eq. .

Alternatively, one can apply the transformation to either (h, h_Y) or (s°,

s_{Y} °

⁠) and investigate how the terms in the entropy balance change. In the former case, one can easily verify that both

Q [n_{t}]

and S(n) are unaltered. In the latter case, instead

\begin{aligned} S (n) & \to S (n) + \sum_{λ} a_{λ} L_{n}^{λ}, \\ Q^{t h r} [n_{t}] & \to Q^{t h r} [n_{t}], and \\ Q^{c h m} [n_{t}] & \to Q^{c h m} [n_{t}] + T \sum_{λ_{b}} a_{λ_{b}} ℓ_{λ_{b}}^{y} \cdot I^{Y} [n_{t}], \end{aligned}

(82)

where we distinguished the thermal and chemical heat contributions.

We thus emphasize that, W_c, G(n), S(n), and Q^chm are not uniquely defined, in contrast to Σ and Q^thr. Despite that, once the gauge is fixed—i.e., the values of the standard-state quantities are chosen—they are useful concepts for characterizing the dissipation of the process. Further discussions on the gauge arising in the work-potential connection will be given in Sec. V C.

Remark

Rather than defining the heat as minus the entropy change in the environment times T, Eqs. (65) and (66), we could have defined it as minus the entropy change in the thermal reservoir times T, Q^thr, thus leaving the chemical part aside. Clearly, this does not affect the EP, but its expression would lose the typical Kelvin–Clausius form, Eq. (72), as it would read $Σ [n_{t}] = Δ S [n_{t}] - \frac{1}{T} Q^{t h r} [n_{t}] - \int_{0}^{t} d τ s_{Y} (τ) \cdot I^{Y} (τ)$ ⁠. These two different but equivalent approaches are not new to nonequilibrium thermodynamics and have been discussed in Ref. 41, Chap. III, Sec. 3, for instance.

D. FT for the chemical work and comparison with previous results

When combining the EP FT (75) with Eq. (73c), we immediately obtain the integral FT for the chemical work

⟨\exp \{- β (W_{c} - Δ G)\}⟩ = 1 .

(83)

However, a Jarzynski-like integral FT^42–45 for the chemical work—i.e., expressions such as $⟨\exp \{- β W_{c}\}⟩ = \exp \{- β Δ G_{e q}\}$ —does not ensue. This relation would require that (i) the process starts and finishes at equilibrium in a closed network, ΔG = ΔG_eq—the condition on the final state can be relaxed, though—and (ii) ΔG_eq is a nonfluctuating quantity along the process so that its exponential can be moved out of the average. However, due to broken conservation laws, G_eq fluctuates along any trajectory of open CRNs.

Let us consider a generic process in which the CRNs is initially closed and at equilibrium, Eq. (44), with a Gibbs free energy $\sum_{{L_{λ}}} P ({L_{λ}}) G_{e q} ({L_{λ}})$ ⁠. The CRN is then open and driven according to some time-dependent protocol, π_τ, for τ ∈ [0, t]. At time t, the CRN is closed again and let to relax to a new equilibrium distribution $p_{n}^{{e q}_{t}}$ ⁠. Since the chemostatting procedure unavoidably breaks some conservation laws, the accessible state space suddenly increases. The final distribution of broken components, $P ({L_{λ_{b}}}; t)$ ⁠, will thus have a support broader than that of the initial distribution, $P ({L_{λ_{b}}}; 0)$ ⁠; see, e.g., Fig. 3. This process is akin to the free expansion of a gas that is initially at equilibrium in a constrained region of space. The crucial point is that the initial state is a constrained, or local, equilibrium with respect to the state space where the dynamics subsequently evolves.

FIG. 3.

View large Download slide

Illustration of the evolution of the probability distribution of the broken components associated with Eq. (31b) in the CRNs in Fig. 1. As the CRN evolves, the state space enlarges as the stochastic dynamics explores states corresponding to different broken components, L_b. The four distributions are obtained by means of 10⁶ trajectories simulated using the stochastic simulation algorithm. All rate constants are equal to 1, whereas the concentrations of the chemostatted species are [A_e] = 17 and [B_e] = 10. The value of the enzyme moiety is L_E = 5.

The stochastic thermodynamics of these processes is characterized by absolute irreversibility.²⁷ Namely, when the EP (74) is integrated over all trajectories to obtain the FT (75), there are some backward trajectories whose corresponding forward probability is vanishing. These are the trajectories leading to values of the broken components not in $s u p p \{P ({L_{λ_{b}}}; 0)\}$ ⁠. Since the EP of these trajectories diverges negatively, see Eq. (74), the expression of the integral FTs (75), as well as (83), is invalidated but can be replaced by $⟨ \exp \{- Σ / k_{B}\} ⟩ = 1 - λ_{S}$ ⁠, where 0 ≤ λ_S ≤ 1 measures the probability of those backward trajectories whose forward one has zero probability.²⁷

Hence, let us assume that $s u p p \{P ({L_{λ_{b}}}; 0)\}$ spans all possible values of ${L_{λ_{b}}}$ so that no absolute irreversibility occurs. By conditioning the average in Eq. (83) upon observation of specific initial and final components (⁠ ${⟨\cdot⟩}_{{L_{λ}}, {L_{λ}^{'}}}$ ⁠), we obtain

\begin{matrix} \sum_{{L_{λ}}} \sum_{{L_{λ}^{'}}} P ({L_{λ}}; 0) P ({L_{λ}^{'}}; t) \\ \times \exp \{β [G_{{e q}_{t}} ({L_{λ}^{'}}) - G_{{e q}_{0}} ({L_{λ}})]\} \\ \times {⟨\exp \{- β W_{c}\}⟩}_{{L_{λ}}, {L_{λ}^{'}}} = 1 . \end{matrix}

(84)

However, this equation cannot be simplified further: since the Gibbs potential depends on the broken components, it fluctuates during the transient dynamics and an average over all components must be taken. As a result, no Jarzynski FT for the chemical work in the Gibbs ensemble can be derived.

In Ref. 19, a Jarzynski relation for the chemical work is derived using the grand canonical ensemble Ref. 19, Eq. (61). Translated into our notation, this result reads

⟨\exp \{- β [W_{c} - Δ (μ^{e q} \cdot n)]\}⟩ = \exp \{- β Δ G_{e q}\},

(85)

where the initial and final equilibrium states are grand canonical

p_{n}^{e q} = \exp \{β [G_{e q} - g_{n} + μ^{e q} \cdot n]\} .

(86)

The grand potential is defined as

G ≔ G - μ^{e q} \cdot n,

(87)

and μ^eq are implicitly defined by

μ_{x}^{e q} \cdot S_{ρ_{i}}^{x} + μ_{y}^{e q} \cdot S_{ρ_{i}}^{y} = 0, for all ρ_{i},

(88)

[Ref. 19, Eq. (27)]. The absence of the exchange transition is due to a different form of chemostatting; see the remark in Sec. II A. The grand potential is naturally suited to describe CRNs in which all species are chemostatted and μ^eq are their chemical potentials. But for most CRNs, where only a subset of species are typically chemostatted, the grand potential is not the most convenient and intuitive potential to work with. The physical interpretation of the contribution −Δ(μ^eq·n) is, for instance, not transparent. In Secs. IV–VIII, we will make use of conservation laws to identify the potential which best describes CRNs where only a subset of species are chemostatted. New work contributions with a transparent physical interpretation will ensue.

IV. CRN-SPECIFIC STOCHASTIC THERMODYNAMICS

We now proceed with our main results. Making use of the conservation laws identified in Sec. II D, we decompose the EP into three fundamental contributions: a potential difference, a contribution due to time-dependent driving, and a minimal set of contributions due to nonconservative chemical forces. To do so, we first decompose the local detailed balance and then proceed with the EP.

A. Entropy production

We start our EP decomposition by partitioning the set of chemostatted species Y into two groups, denoted by Y_p and Y_f. Likewise, the corresponding exchanged species are denoted by y_p and y_f, respectively. The former group is composed by a minimal set of chemostatted species which—when starting from the closed CRN—break all broken conservation laws. In other words, each entry of Y_p breaks exactly one distinct conservation law. The remaining chemostatted species form the latter group. For a given CRN, our partitioning is not unique, but the number of y_p and y_f is uniquely defined: $N_{y_{p}} = N_{λ_{b}}$ and $N_{y_{f}} = N_{y} - N_{λ_{b}}$ ⁠, respectively; see example 5.

We now notice that the linear independence of { $ℓ_{λ}$ } implies that the matrix whose rows are ${ℓ_{λ_{b}}^{y_{p}}}$ is nonsingular. We will denote by ${{\bar{ℓ}}_{λ_{b}}^{y_{p}}}$ the column vectors of the inverse of the latter matrix. By making use of this important property, we can recast the identity

Δ_{ρ} L_{n}^{λ_{b}} \equiv ℓ_{λ_{b}} \cdot S_{ρ} \equiv ℓ_{λ_{b}}^{x} \cdot S_{ρ}^{x} + ℓ_{λ_{b}}^{y_{p}} \cdot S_{ρ}^{y_{p}} + ℓ_{λ_{b}}^{y_{f}} \cdot S_{ρ}^{y_{f}}

(89)

into

S_{ρ}^{y_{p}} = Δ_{ρ} M_{n}^{y_{p}} - \sum_{λ_{b}} {\bar{ℓ}}_{λ_{b}}^{y_{p}} [ℓ_{λ_{b}}^{x} \cdot S_{ρ}^{x} + ℓ_{λ_{b}}^{y_{f}} \cdot S_{ρ}^{y_{f}}],

(90)

where

M_{n}^{y_{p}} ≔ \sum_{λ_{b}} {\bar{ℓ}}_{λ_{b}}^{y_{p}} L_{n}^{λ_{b}} .

(91)

Mindful that $S_{ρ}^{Y} = - S_{ρ}^{y}$ and $ℓ_{λ_{b}}^{x} \cdot S_{ρ_{e}}^{x} = 0$ for all ρ_e, one can use Eq. (90) to rewrite the chemical work along reactions as

- μ_{Y} \cdot S_{ρ}^{Y} = Δ_{ρ} [μ_{Y_{p}} \cdot M_{n}^{y_{p}}] - F_{Y_{f}} \cdot S_{ρ}^{Y_{f}},

(92)

where

F_{Y_{f}} ≔ μ_{Y_{f}} - μ_{Y_{p}} \cdot \sum_{λ_{b}} {\bar{ℓ}}_{λ_{b}}^{y_{p}} ℓ_{λ_{b}}^{y_{f}} .

(93)

A reformulation of the local detailed balance Eq. (53) readily ensues

\ln \frac{w_{ρ} (n)}{w_{- ρ} (n + S_{ρ})} = - β () Δ_{ρ} g_{n} + F_{Y_{f}} \cdot S_{ρ}^{Y_{f}})),

(94)

where

g_{n} ≔ g_{n} - μ_{Y_{p}} \cdot M_{n}^{y_{p}} .

(95)

We now notice that the expression of the potential $g_{n}$ is reminiscent of a Legendre transform of g_n with respect to $M_{n}^{y_{p}}$ ⁠, in which $μ_{Y_{p}}$ are the conjugated intensive fields. To reveal the physical meaning of $M_{n}^{y_{p}}$ ⁠, let us consider the case in which the broken conservation laws correspond to moieties, see Sec. II D, and hence each species can be thought of as a composition of these. Through y_p, some combinations of these moieties are exchanged with the environment. The entries of $M_{n}^{y_{p}}$ quantify the total abundance of these combinations in state n, and hence we refer to $M_{n}^{y_{p}}$ as the moiety population vector. In view of this and the fact that (in general) not all moieties are exchanged, one can interpret $g_{n}$ as the semigrand Gibbs free energy of the state n.³⁴ Note also that, from the definition of broken conservation law, Eq. (27), it follows that $Δ_{ρ_{i}} M_{n}^{y_{p}} = 0$ ⁠, for all ρ_i—viz., internal reactions never create or destroy moieties—whereas for ρ_e only we have that $Δ_{ρ_{e}} M_{n}^{y_{p}} \neq 0$ —viz., exchange reactions introduce or remove moieties. We also mention that an alternative interpretation of $g_{n}$ can be given once we rewrite it as

g_{n} ≔ g_{n} - \sum_{λ_{b}} f_{λ_{b}} L_{n}^{λ_{b}},

(96)

where

f_{λ_{b}} ≔ μ_{Y_{p}} \cdot {\bar{ℓ}}_{λ_{b}}^{y_{p}} .

(97)

In this form, $g_{n}$ is reminiscent of a Legendre transform with respect to the broken components ${L_{n}^{λ_{b}}}$ ⁠, in which ${f_{λ_{b}}}$ are the conjugated intensive fields.

In the second term on the rhs of Eq. (94), $F_{Y_{f}}$ identifies chemical potential gradients imposed by the chemostats on the CRN. Its entries, denoted by ${F_{y_{f}}}$ ⁠, for $y_{f} = 1, \dots, N_{y_{f}}$ ⁠, are a maximal independent set of nonconservative chemical forces: if and only if $F_{Y_{f}} = 0$ ⁠, then the rhs of Eq. (94) is conservative. In this case, the CRN is detailed-balanced since the steady-state probability distribution defined by $p_{n}^{e q} \propto \exp \{- β g_{n}\}$ satisfies the detailed balance property, Eq. (14). Since ${F_{y_{f}}}$ make the CRN non-detailed balanced, we refer to them as fundamental nonconservative chemical forces. Equation (94) is our first major result.

To proceed with our EP decomposition, we combine Eqs. (73b) and (94)

\begin{matrix} T Σ [n_{t}] & = & k_{B} T \ln \frac{p_{n_{0}} (0)}{p_{n_{t}} (t)} - \int_{0}^{t} d τ \sum_{ρ, n} Δ_{ρ} g_{n} (τ) j_{ρ} (n, τ) \\ + \sum_{y_{f}} W_{y_{f}}^{n c} [n_{t}], \end{matrix}

(98)

where

W_{y_{f}}^{n c} [n_{t}] ≔ \int_{0}^{t} d τ F_{y_{f}} (τ) I_{y_{f}} (τ) .

(99)

${I_{y_{f}} (τ)}$ ⁠, for $y_{f} = 1, \dots, N_{y_{f}}$ ⁠, denote the entries of the instantaneous external currents corresponding to Y_f, Eq. (67). We now recall that $g_{n}$ is a state function; hence,

Δ g [n_{t}] = W_{d} [n_{t}] + \int_{0}^{t} d τ \sum_{ρ, n} Δ_{ρ} g_{n} (τ) j_{ρ} (n, τ),

(100)

where

\begin{matrix} W_{d} [n_{t}] ≔ \int_{0}^{t} d τ {([\partial_{τ} g_{n} (τ)]|}_{n_{τ}} = \int_{0}^{t} d τ [- \partial_{τ} μ_{Y_{p}} (τ)] \cdot M_{n_{τ}}^{y_{p}} . \end{matrix}

(101)

Therefore, combining Eqs. (98) and (100), we obtain

T Σ [n_{t}] = - Δ G [n_{t}] + W_{d} [n_{t}] + \sum_{y_{f}} W_{y_{f}}^{n c} [n_{t}],

(102)

where the first term is the difference of stochastic semigrand Gibbs potential

G (n) ≔ k_{B} T \ln p_{n} + g_{n} .

(103)

The EP decomposition in Eq. (102) is a major result of our paper. The first term on the rhs constitutes the conservative force contribution of the EP. It describes the dissipation due to overall changes of thermodynamic state variables: enthalpy, H(n), entropy, S(n), and chemical energy ${μ_{Y_{p}} \cdot M_{n}^{y_{p}}}$ ⁠. The second term, Eq. (101), arises in the presence of time-dependent driving and accounts for the changes caused by manipulations of the chemical potentials $μ_{Y_{p}}$ ⁠. As it is a controlled way of changing the Gibbs free energy landscape of the CRN, we refer to it as driving chemical work. Finally, for each exchanged species Y_f, a nonconservative force contribution (99) arises, ${W_{y_{f}}^{n c}}$ ⁠. All together, they account for the chemical energy flowing between different chemostats across the CRN and we refer to them as nonconservative chemical work contributions. Equation (102) holds for an arbitrary CRN, yet it is CRN-specific, as it is derived using the topological properties of the CRN encoded in the conservation laws. To gain more intuition, we now focus on specific classes of CRNs, whose resulting decomposition is summarized in Table II. In Sec. IV B, we continue our discussion on the work contributions W_d and ${W_{y_{f}}^{n c}}$ ⁠, whereas in example 5 and in Sec. VIII, we evaluate them for specific models. Finally, in Secs. VI and VII, we will further explore the implications of Eq. (102).

TABLE II.

Entropy production for specific processes. “0” (respectively, “✓”) denotes a vanishing (respectively, a finite) contribution.

Dynamics	$- Δ 𝒢$	W_d	W^nc
Autonomous detailed-balanced	✓	0	0
Unconditionally detailed-balanced	✓	✓	0
Autonomous	✓	0	✓
Nonequilibrium steady state	0	0	✓

Dynamics	$- Δ 𝒢$	W_d	W^nc
Autonomous detailed-balanced	✓	0	0
Unconditionally detailed-balanced	✓	✓	0
Autonomous	✓	0	✓
Nonequilibrium steady state	0	0	✓

TABLE III.

List of symbols used throughout the text. The physical quantity that they denote and the equation number in which they are defined are also reported.

Symbol	Physical quantity	Equations
S_ρ	Stoichiometric vectors	(3)
$S$	Stoichiometric matrix	(4) and (5)
$w_{ρ}$ (n)	Stochastic reaction rates	(13)
J_ρ(τ)	Instantaneous reaction fluxes	(17)
$ℓ$	Conservation laws	(23)
L_n	Component	(24)
c	Stoichiometric cycles	(32)
$μ$ _Y	Chemostats chemical potential	(54)
g_n	Gibbs free energy of n	(42)
G(n)	Stochastic Gibbs potential	(47)
⟨G(n)⟩	Nonequilibrium Gibbs potential	(48)
Z	Closed-CRN partition function	(43)
A_ρ(n)	Reaction affinity	(56)
s_n	Entropy of n	(58)
S(n)	Stochastic entropy	(57)
⟨S(n)⟩	Gibbs–Shannon entropy	(60)
H(n)	Enthalpy	(61) and (166)
Q (⁠ $⟨ \dot{Q} ⟩$ ⁠)	Heat flow (rate)	(66) and (167)
W_c (⟨ ${\dot{W}}_{c}$ ⟩)	Chemical work (rate)	(68) and (170)
I^Y	Instantaneous external currents	(67)
Σ (⁠ $⟨ \dot{Σ} ⟩$ ⁠)	Entropy production (rate)	(74) and (171)
$M_{n}^{y_{p}}$	Moiety population vector	(91)
$F_{Y_{f}}$	Fundamental forces	(93)
$I^{Y_{f}}$	Fundamental external currents	(79)
$g_{n}$	Semigrand Gibbs free energy of n	(95)
$G (n)$	Stoch. semigrand Gibbs pot.	(103)
$⟨ G (n) ⟩$	Noneq. semigrand Gibbs pot.	(173)
$Z$	Open-CRN partition function	(119)
W_d (⟨ ${\dot{W}}_{d}$ ⟩)	Driving chem. work (rate)	(101) and (174)
$W_{y_{f}}^{n c}$ (⁠ $⟨ {\dot{W}}_{y_{f}}^{n c} ⟩$ ⁠)	Nonconservative chem. work (rate)	(99) and (175)
$H (n)$	Semigrand enthalpy	(115) and (177)
$A_{η}$	Stoichiometric cycle affinity	(126)
$J_{η}$	Stoichiometric cycle current	(135)
Γ_η	Nonconservative cycle chem. work	(133) and (179)

Symbol	Physical quantity	Equations
S_ρ	Stoichiometric vectors	(3)
$S$	Stoichiometric matrix	(4) and (5)
$w_{ρ}$ (n)	Stochastic reaction rates	(13)
J_ρ(τ)	Instantaneous reaction fluxes	(17)
$ℓ$	Conservation laws	(23)
L_n	Component	(24)
c	Stoichiometric cycles	(32)
$μ$ _Y	Chemostats chemical potential	(54)
g_n	Gibbs free energy of n	(42)
G(n)	Stochastic Gibbs potential	(47)
⟨G(n)⟩	Nonequilibrium Gibbs potential	(48)
Z	Closed-CRN partition function	(43)
A_ρ(n)	Reaction affinity	(56)
s_n	Entropy of n	(58)
S(n)	Stochastic entropy	(57)
⟨S(n)⟩	Gibbs–Shannon entropy	(60)
H(n)	Enthalpy	(61) and (166)
Q (⁠ $⟨ \dot{Q} ⟩$ ⁠)	Heat flow (rate)	(66) and (167)
W_c (⟨ ${\dot{W}}_{c}$ ⟩)	Chemical work (rate)	(68) and (170)
I^Y	Instantaneous external currents	(67)
Σ (⁠ $⟨ \dot{Σ} ⟩$ ⁠)	Entropy production (rate)	(74) and (171)
$M_{n}^{y_{p}}$	Moiety population vector	(91)
$F_{Y_{f}}$	Fundamental forces	(93)
$I^{Y_{f}}$	Fundamental external currents	(79)
$g_{n}$	Semigrand Gibbs free energy of n	(95)
$G (n)$	Stoch. semigrand Gibbs pot.	(103)
$⟨ G (n) ⟩$	Noneq. semigrand Gibbs pot.	(173)
$Z$	Open-CRN partition function	(119)
W_d (⟨ ${\dot{W}}_{d}$ ⟩)	Driving chem. work (rate)	(101) and (174)
$W_{y_{f}}^{n c}$ (⁠ $⟨ {\dot{W}}_{y_{f}}^{n c} ⟩$ ⁠)	Nonconservative chem. work (rate)	(99) and (175)
$H (n)$	Semigrand enthalpy	(115) and (177)
$A_{η}$	Stoichiometric cycle affinity	(126)
$J_{η}$	Stoichiometric cycle current	(135)
Γ_η	Nonconservative cycle chem. work	(133) and (179)

1. Autonomous detailed-balanced CRNs

The CRN is autonomous and all fundamental forces vanish. The trajectory EP becomes minus a potential difference

T Σ [n_{t}] = - Δ G [n_{t}] .

(104)

We will prove in Sec. VII that this is the class of open CRNs which relax to equilibrium and in which the average potential $⟨ G ⟩$ is minimized at equilibrium by the dynamics described by CME (11).

2. Unconditionally detailed-balanced CRNs

The set of species Y_f is empty—i.e., each exchanged species breaks a conservation law—and no fundamental force arises. Hence, these CRNs are detailed-balanced irrespective of the values of μ_Y, but the time-dependent driving prevents them from reaching equilibrium, and their EP reads

T Σ [n_{t}] = - Δ G [n_{t}] + W_{d} [n_{t}] .

(105)

3. Autonomous CRNs

The driving work vanishes and the forces are constant in time. Hence, the EP becomes

T Σ [n_{t}] = - Δ G [n_{t}] + \sum_{y_{f}} F_{y_{f}} I_{y_{f}} [n_{t}] .

(106)

The nonconservative chemical work displays a typical current–force structure. In the long time limit, $Δ G [n_{t}]$ is typically subextensive in time, and we obtain the EP typical of nonequilibrium steady states

T Σ [n_{t}] \overset{t \to \infty}{=} \sum_{y_{f}} F_{y_{f}} I_{y_{f}} [n_{t}]

(107)

[see Eq. (79)]. In other words, $T Σ [n_{t}]$ is dominated by the dissipative flows of chemicals across the CRN.

Remark

For CRN with infinite number of species and reactions—e.g., aggregation–fragmentation and polymerization processes^46–48—the CRN may undergo steady growth regimes in which $Δ G$ is not subextensive in time and cannot be neglected in the long-time limit.

Remark

Our EP decomposition is not unique and different expressions for $g_{n}$ and $F_{Y_{f}}$ correspond to different ways of partitioning Y into Y_p and Y_f.

Example 5.

For the open CRN in Fig. 1, the chemostatted species can be split into Y_p and Y_f in two possible—and trivial—ways: either A is regarded as the species breaking the conservation law , or B. We consider the former choice, y_p = (A) and y_f = (B). Since

ℓ_{A}^{b} = 1

⁠, the only entry of the moiety vector reads

M_{n}^{A} = n_{E^{*}} + n_{E^{* *}} + n_{A} + n_{B} = L_{n}^{b},

(108)

which is equal to the total abundance of the A–B moiety. The intensive variable conjugated to the broken conservation law is equal to the chemical potential of A_e

f_{b} = μ_{A_{e}} .

(109)

The potential thus readily follows from Eq. , or equivalently Eq. ,

g_{n} = g_{n} - μ_{A_{e}} M_{n}^{A} .

(110)

The instantaneous driving work rate associated with any manipulation of the latter potential is

{\dot{W}}_{d} (n) = - \partial_{t} μ_{A_{e}} M_{n}^{A} .

(111)

Once integrated over a trajectory, it gives the driving work, Eq. . Since y_f = (B), the conjugated fundamental chemical force reads

F_{B_{e}} = μ_{B_{e}} - μ_{A_{e}}

(112)

and the instantaneous dissipative contribution due to this force is

{\dot{W}}_{B_{e}}^{n c} = F_{B_{e}} I_{B_{e}},

(113)

where

I_{B_{e}} = J_{+ b} - J_{- b}

⁠. When integrated over a trajectory, it measures the work spent to sustain a current between A_e and B_e across the CRN. A pictorial illustration of the work contributions is given in Fig. 4. The trajectory EP thus reads

\begin{matrix} T Σ [n_{t}] & = & \int_{0}^{t} d τ [) - \partial_{τ} μ_{A_{e}} (τ) M_{n}^{A}]) {|)}_{n_{τ}} - Δ G [n_{t}] \\ + \int_{0}^{t} d τ F_{B_{e}} (τ) I_{B_{e}} (τ) . \end{matrix}

(114)

FIG. 4.

View large Download slide

Pictorial illustration of the work contributions. The driving one arises when the chemical potential of the chemostat A_e changes in time. The nonconservative chemical work, instead, characterizes the sustained conversion of A into B.

B. Energy balance

In Eq. (102), the driving and nonconservative chemical work, W_d and ${W_{y_{f}}^{n c}}$ ⁠, emerge as dissipative contributions. To strengthen their interpretation as work contributions, we now show that they can also be described as part of an energy balance. For this purpose, let us introduce the semigrand enthalpy⁴⁹

H (n) ≔ H (n) - μ_{Y_{p}} \cdot M_{n}^{y_{p}} = G (n) + T S (n) .

(115)

This CRN-specific potential quantifies the portion of energy which is not attributed to volume (−pV, where p is the external pressure) and exchanged moieties, $μ_{Y_{p}} \cdot M_{n}^{y_{p}}$ ⁠. It accounts for the energy stored in its internal chemical composition, i.e., the internal species x and the unbroken components ${L_{λ_{u}}}$ ⁠. When combining its definition with the enthalpy and entropy balances, Eqs. (69), (72), and (102), we obtain

Δ H [n_{t}] = Q [n_{t}] + W_{d} [n_{t}] + \sum_{y_{f}} W_{y_{f}}^{n c} [n_{t}],

(116)

viz., the overall change of semigrand enthalpy is equal to the sum of heat flow, driving, and nonconservative chemicalwork. By analogy with Eq. (69), this can be interpreted as a CRNspecific formulation of the first law.

In Sec. III C, we introduced the chemical work as the Gibbs free energy exchanged with the chemostats, Eq. (68). By comparing Eqs. (61) and (116), we obtain its relationship with W_d and ${W_{y_{f}}^{n c}}$

W_{c} [n_{t}] - Δ [μ_{Y_{p}} \cdot M_{n}^{y_{p}}] = W_{d} [n_{t}] + \sum_{y_{f}} W_{y_{f}}^{n c} [n_{t}] .

(117)

We emphasize that in contrast to the chemical work, the driving one does not account for direct exchanges of Gibbs free energy, but it captures the instantaneous changes of the chemostats Gibbs free energy.

Remark

The driving work is reminiscent of the mechanical work as defined in stochastic thermodynamics. In this framework, $W_{m e c h} [n_{t}] = \int_{0}^{t} d τ {(\partial_{τ} E_{n} (τ)|}_{n_{τ}}$ describes internal energy changes due to external time-dependent control; see, e.g., Refs. 44 and 50. In CRNs, the time-dependent control is exerted via the chemostats, and $W_{d} [n_{t}]$ indeed accounts for this fact.

C. Equilibrium of open CRNs

We have already seen that in the absence of fundamental forces, the rhs of the local detailed balance (94) becomes a state function difference. The steady-state probability distribution

p_{e q} (n | {L_{λ_{u}}}) = \frac{\exp \{- β g_{n}\}}{Z ({L_{λ_{u}}})} \prod_{λ_{u}} δ [) L_{n}^{λ_{u}}, L_{λ_{u}}])

(118)

satisfies the detailed balance property (53) and therefore characterizes the equilibrium of open CRNs. Not accidentally, the relationship between the partition function $Z ({L_{λ_{u}}})$ and that of closed CRNs, Eq. (43),

\begin{matrix} Z ({L_{λ_{u}}}) & = & \sum_{m} \exp \{- β g_{m}\} \prod_{λ_{u}} δ [) L_{m}^{λ_{u}}, L_{λ_{u}}]) \\ = & \sum_{{L_{λ_{b}}}} \exp \{β \sum_{λ_{b}} f_{λ_{b}} L_{λ_{b}}\} Z ({L_{λ}}), \end{matrix}

(119)

is akin to that between canonical and grand canonical partition functions; see, e.g., Ref. 51. With an ensemble of unbroken components, $P ({L_{λ_{u}}})$ ⁠, the constrained equilibrium distribution reads

\begin{matrix} p_{n}^{e q} & = & \sum_{{L_{λ_{b}}}} p_{e q} (n | {L_{λ_{u}}}) P ({L_{λ_{u}}}) \\ = & p_{e q} (n | {L_{n}^{λ_{u}}}) P ({L_{n}^{λ_{u}}}), \end{matrix}

(120)

where $p_{e q} (n | {L_{n}^{λ_{u}}})$ is the probability distribution of observing the state n given its stoichiometric compatibility class. Equation (120) thus generalizes the equilibrium probability distribution (44) to open CRNs.

Importantly, the average semigrand Gibbs potential (103) takes its minimum value at $p_{n}^{e q}$ ⁠, Eq. (120), where it reduces to the equilibrium semigrand Gibbs potential

G_{e q} ({L_{λ_{u}}}) = - k_{B} T \ln Z ({L_{λ_{u}}}) + k_{B} T \ln P ({L_{λ_{u}}}),

(121)

averaged over $P ({L_{λ_{u}}})$ ⁠. Indeed,

⟨ G ⟩ - {⟨ G_{e q} ⟩}_{L_{u}} = ⟨ G - G_{e q} ⟩ = k_{B} T D (p ‖ p_{e q}) \geq 0,

(122)

where

{⟨G_{e q}⟩}_{L_{u}} \equiv \sum_{{L_{λ_{u}}}} P ({L_{λ_{u}}}) G_{e q} ({L_{λ_{u}}}) .

(123)

The first equality follows from the fact that $G_{e q}$ is nonfluctuating, since it depends solely on the unbroken components. As for the Gibbs free energy in closed CRNs, we will show later (Sec. VII) that Eq. (122) quantifies the average dissipation during the relaxation to equilibrium.

D. Dissipation balance along stoichiometric cycles

We can now formulate the EP decomposition in terms of stoichiometric cycle affinities. These are defined as the sum of the transition affinities along stoichiometric cycles ${c \equiv ρ_{c 1}, ρ_{c 1}, \dots, ρ_{c N_{c}}}$

\begin{matrix} A & ≔ & A_{ρ_{c} 1} (n) + A_{ρ_{c} 2} (n + S_{ρ_{c} 1}) + \dots \dots \\ + A_{ρ_{c} N_{c}} (n + \sum_{j = 1}^{N_{c} - 1} S_{ρ_{c} j}) . \end{matrix}

(124)

Using Eq. (56), and the fact that −Δ_ρG(n) vanishes when summed over the loop c, we obtain

A = - μ_{Y} \cdot \sum_{i = 1}^{N_{c}} S_{ρ_{c} i}^{Y} = - μ_{Y} \cdot \sum_{ρ} S_{ρ}^{Y} c_{ρ} .

(125)

Since $\sum_{ρ} S_{ρ}^{Y} c_{ρ}^{α} = 0$ ⁠, those evaluated along the stoichiometric cycles of the closed CRN, {c^α}, always vanish. By contrast, those along the emergent cycles, {c^η}, do not vanish in general

A_{η} = μ_{Y} \cdot C_{η}^{Y}

(126)

[see Eq. (36)]. These affinities can be thus understood as the chemical potential gradient imposed by the chemostats on the cycle.

To rewrite the EP (102) in terms ${A_{η}}$ ⁠, let us highlight their relationship with the fundamental forces

A_{η} = F_{Y_{f}} \cdot C_{η}^{Y_{f}},

(127)

which is obtained when summing the local detailed balance (94) along {c^η} as in Eq. (124). Since the matrix whose columns are ${C_{η}^{Y_{f}}}$ is square and nonsingular—as it can be deduced from the linear independence of the set of emergent cycles—we can invert it and write

F_{Y_{f}} = \sum_{η} {\bar{C}}_{η}^{Y_{f}} A_{η},

(128)

where ${{\bar{C}}_{η}^{Y_{f}}}$ denote the rows of the inverse matrix. This relation clarifies the one-to-one correspondence which lies between ${F_{y_{f}}}$ and ${A_{η}}$ ⁠. Inserting the last expression in the local detailed balance, Eq. (94), we obtain

\ln \frac{w_{ρ} (n)}{w_{- ρ} (n + S_{ρ})} = - β () Δ_{ρ} g_{n} - \sum_{η} A_{η} ζ_{η, ρ})),

(129)

where the coefficients

ζ_{η, ρ} ≔ - {\bar{C}}_{η}^{Y_{f}} \cdot S_{ρ}^{Y_{f}}

(130)

quantify how much each reaction contributes to the emergent cycles. Algebraically, the row vectors {ζ_η} are dual to the cycles, {c^η}

ζ_{η} \cdot c^{η^{'}} = - \sum_{ρ} {\bar{C}}_{η}^{Y_{f}} \cdot S_{ρ}^{Y_{f}} c_{ρ}^{η^{'}} = {\bar{C}}_{η}^{Y_{f}} \cdot C_{η^{'}}^{Y_{f}} = δ_{η, η^{'}} .

(131)

As previously done for Eq. (102), when integrating the trajectory EP (73b) with the local detailed balance (129), we obtain

T Σ [n_{t}] = - Δ G [n_{t}] + W_{d} [n_{t}] + \sum_{η} Γ_{η} [n_{t}] .

(132)

The stochastic semigrand Gibbs potential and the driving work read as in Eqs. (103) and (101), respectively. For each emergent stoichiometric cycle,

Γ_{η} [n_{t}] ≔ \int_{0}^{t} d τ A_{η} (τ) \sum_{ρ} ζ_{η, ρ} J_{ρ} (τ)

(133)

quantifies the chemical work spent to sustain the related cyclic flow of chemicals. For autonomous CRNs,

T Σ [n_{t}] = - Δ G [n_{t}] + \sum_{η} A_{η} J_{η} [n_{t}],

(134)

where

J_{η} [n_{t}] ≔ \int_{0}^{t} d τ \sum_{ρ} ζ_{η, ρ} J_{ρ} (τ)

(135)

quantifies the integrated current along the cycle η. In the long-time limit, in which $Δ G [n_{t}]$ is negligible, we obtain

T Σ [n_{t}] \overset{t \to \infty}{=} \sum_{η} A_{η} J_{η} [n_{t}] .

(136)

When all emergent cycle affinities vanish—as well as when no emergent cycle is created—the CRN becomes detailed-balanced, in agreement with the Kolmogorov–Wegscheider condition.^52–54

We emphasize that the cycle chemical work contributions and currents, Eqs. (133) and (135), can be written as combinations of fundamental external currents, ${I^{Y_{f}}}$ Eq. (67), via Eq. (130). The added value of Eq. (102) over (132) lies in the fact that each force is conjugated to the external current of only one external species.

Remark

An alternative approach that can be used for cycle EP decompositions is the graph-theoretic one based on the identification of the loops appearing in the network of transitions.^2,55 Once these loops are identified, they can be sorted according to the chemostats they are coupled to, as these determine their affinity; see Eq. (124). Equivalently, loops are classified according to the stoichiometric cycle they correspond to. In Ref. 56, a graph-theoretic approach based on loop affinities led to the expression analogous to Eq. (136). By contrast, our cycle EP decomposition is based on a stoichiometric approach: emergent cycles are directly identified by the kernels of $S_{i}$ and $S$ ⁠.

This observation points out the redundancy which is intrinsic in bare graph-theoretic EP decompositions: many loops may be coupled to the same set of reservoirs and thus carry the same affinity, while many others may carry a vanishing affinity—for CRN, these latter are those corresponding to stoichiometric cycles of the closed network, {c^α}. For generic networks, a systematic way of identifying these so-called symmetries was derived in Ref. 57, whereas in Ref. 58, they are used to formulate generic thermodynamic—rather than mere graph-theoretic—EP decompositions.

Example 6.

The emergent cycle affinity corresponding to the emergent stoichiometric cycle reads

A = μ_{B_{e}} - μ_{A_{e}} = F_{B_{e}} .

(137)

The contributions to the corresponding cycle current follow from Eq.

ζ = \begin{matrix} \begin{matrix} + 1 & + 2 & + 3 & + 4 & + a & + b \end{matrix} \\ (\begin{matrix} 0 & 0 & 0 & 0 & 0 & - 1 \end{matrix}) \end{matrix}

(138)

The entries corresponding to the backward reactions are minus those of the forward. Notice that, since the CRN has exactly one emergent cycle, the force and cycle EP decompositions are identical; see Eq. .

V. SEMIGRAND GIBBS POTENTIAL

Here we further elaborate on equilibrium distributions and semigrand Gibbs potentials by addressing three points: (i) the relationship between Eq. (120) and the equilibrium distributions as expressed in chemical reaction network theory; (ii) the role of conservation laws for characterizing the dissipation of CRNs subject to sequential introduction of exchanged species; (iii) the gauge freedom intrinsic to the definition of driving work. This section can be skipped at a first read.

A. Equilibrium distributions in chemical reaction network theory

In Ref. 22 (see also Ref. 59), equilibrium distributions of CRNs are proven to be multi-Poissonian

\begin{matrix} p_{e q} (n | {L_{λ_{u}}}) = \frac{\exp \{n \cdot \ln \{{[z]}_{e q} V\}\}}{n! Z ({L_{λ_{u}}})} \prod_{λ_{u}} δ [) L_{n}^{λ_{u}}, L_{λ_{u}}]), \end{matrix}

(139)

where [z]_eq is the equilibrium concentration distribution of the same CRN described by a set of deterministic rate equations. $Z ({L_{λ_{u}}})$ is again a normalizing factor. To highlight the relationship between this equation and Eqs. (120) and (86), we need to recall that, for deterministic CRNs, thermodynamic equilibrium is defined by the fact that chemical potential differences along all reactions vanish; see Eqs. (88) and (A8). As observed in Ref. 25, this entails that

μ^{e q} = \sum_{λ} f_{λ} ℓ_{λ},

(140)

where {f_λ} are real coefficients depending on μ_Y and ${L_{λ_{u}}}$ ⁠. Those related to the broken components, ${f_{λ_{b}}}$ ⁠, are indeed those appearing in Eq. (97). Using the expression of chemical potential valid in the thermodynamic limit, Eq. (A7), we therefore have

\ln \{{[z]}_{e q} V\} = - β (μ ° - k_{B} T \ln n_{s} - \sum_{λ} f_{λ} ℓ_{λ}),

(141)

from which

\begin{matrix} n \cdot \ln \{{[z]}_{e q} V\} - \ln n! & = & - β (g_{n} - μ^{e q} \cdot n) \\ = & - β () g_{n} - \sum_{λ_{u}} f_{λ_{u}} L_{n}^{λ_{u}})) \end{matrix}

(142)

ensues. At this point, Eqs. (86), (118), and (139) appear identical up to $\sum_{λ_{u}} f_{λ_{u}} L_{n}^{λ_{u}}$ ⁠. However, since this term involves only the unbroken components it vanishes in Eq. (139). This shows the connection between the CRN theoretical and thermodynamic expression of equilibrium distributions.

B. Hierarchies of equilibriums

Here we show that when starting from a closed CRN, a sequential introduction of exchange reactions that keep the CRN detailed balanced drives it down in the semigrand Gibbs potential by equilibrating previously constrained degrees of freedom: the conservation laws; see Fig. 5. Let us imagine a closed CRN whose initial probability distribution is $p_{n} (0) = \sum_{{L_{λ}}} p_{0} (n | {L_{λ}}) P_{0} ({L_{λ}})$ ⁠, where $P_{0} ({L_{λ}}) = \prod_{λ} P_{0}^{λ} (L_{λ})$ ⁠, i.e., different components are independently distributed. As it relaxes to equilibrium, P₀({L_λ}) will not change, while p₀(n|{L_λ}) will relax to Eq. (41). The average dissipation is

\begin{matrix} T ⟨Σ⟩ & = & - Δ ⟨G⟩ \\ = & \sum_{{L_{λ}}} P_{0} ({L_{λ}}) [k_{B} T \sum_{n} p (n | {L_{λ}}) \ln \frac{p (n | {L_{λ}})}{p_{e q} (n | {L_{λ}})}] \\ \equiv & \sum_{{L_{λ}}} P_{0} ({L_{λ}}) [) - Δ ⟨G ({L_{λ}})⟩]) . \end{matrix}

(143)

FIG. 5.

View large Download slide

Pictorial representation of the hierarchy of equilibrium states and the semigrand Gibbs free energy drops following the relaxation to equilibrium when conservation laws are broken.

This expression is obtained when combining the properties of the Gibbs potential, Eq. (49), with the equilibrium distribution of closed CRNs, Eq. (44). It shows that the average drop of Gibbs free energy can be expressed as the weighted average of the drops of Gibbs free energy at given components, $- Δ ⟨G ({L_{λ}})⟩$ ⁠.

We now open the CRN by chemostatting one species. Hence, one conservation law is broken, e.g., the total mass $ℓ_{λ_{1}}$ ⁠, and the CRN relaxes to a new equilibrium, Eq. (120), whose partition function is denoted by $Z_{λ_{1}}$ ⁠, Eq. (119). Clearly, $P_{0}^{λ} (L_{λ})$ ⁠, for λ ≠ λ₁, will not change during the relaxation, and we can rewrite the new equilibrium as

\begin{matrix} p_{e q}^{(λ_{1})} (n) & = & \frac{\exp \{- β g_{n} + β f_{λ_{1}} L_{n}^{λ_{1}}\}}{Z_{λ_{1}} ({L_{n}^{λ}}_{λ \neq λ_{1}})} \prod_{λ \neq λ_{1}} P_{0}^{λ} (L_{n}^{λ}) \\ = & \underset{= p_{e q} (n | {L_{n}^{λ}})}{\underset{︸}{\frac{\exp \{- β g_{n}\}}{Z ({L_{n}^{λ}})}}} \underset{= : P_{e q} (L_{n}^{λ_{1}} | {L_{n}^{λ}}_{λ \neq λ_{1}})}{\underset{︸}{\frac{Z ({L_{n}^{λ}}) \exp \{β f_{λ_{1}} L_{n}^{λ_{1}}\}}{Z_{λ_{1}} ({L_{n}^{λ}}_{λ \neq λ_{1}})}}} \prod_{λ \neq λ_{1}} P_{0}^{λ} (L_{n}^{λ}) . \end{matrix}

(144)

The first term is the equilibrium distribution of the closed CRN, while the second can be interpreted as the equilibrium distribution of the broken component for a given unbroken component. In other words, the final equilibrium can be understood as a closed CRN equilibrium with an equilibrium probability distribution over the broken component. Hence, the average amount of semigrand Gibbs free energy, $G_{λ_{1}} (n) = G (n) - f_{λ_{1}} L_{n}^{λ_{1}}$ ⁠, dissipated during the relaxation can be written as

\begin{matrix} - Δ ⟨G_{λ_{1}}⟩ & = & k_{B} T \sum_{n} p_{e q} (n | {L_{n}^{λ}}) \prod_{λ} P_{0}^{λ} (L_{n}^{λ}) \\ \times \ln \frac{P_{0}^{λ_{1}} (L_{n}^{λ_{1}})}{P_{e q} (L_{n}^{λ_{1}} | {L_{n}^{λ}}_{λ \neq λ_{1}})}, \end{matrix}

(145)

upon application of Eq. (122) with the distributions (44) and (144). When rewriting this expression as a sum over all values of the components and performing the summation over the states of p_eq(n|{L^λ}), we finally obtain

\begin{matrix} - Δ ⟨G_{λ_{1}}⟩ & = & \sum_{{L_{λ}}_{λ \neq λ_{1}}} P_{0}^{λ} (L_{λ}) \\ \times [\sum_{L_{λ_{1}}} P_{0}^{λ_{1}} (L_{λ_{1}}) k_{B} T \ln \frac{P_{0}^{λ_{1}} (L_{λ_{1}})}{P_{e q} (L_{λ_{1}} | {L_{λ}}_{λ \neq λ_{1}})}] \\ = & \sum_{{L_{λ}}_{λ \neq λ_{1}}} P_{0}^{λ} (L_{λ}) [- Δ ⟨G_{λ_{1}} ({L_{λ}}_{λ \neq λ_{1}})⟩] . \end{matrix}

(146)

In the first line, we recognize the relative entropy between the initial probability of the broken component, $P_{0}^{λ_{1}} (L_{λ_{1}})$ ⁠, and the equilibrium one, $P_{e q} (L_{λ_{1}} | {L_{λ}}_{λ \neq λ_{1}})$ ⁠. It is equal to the difference of semigrand Gibbs free energy at a given component, as highlighted in the second line. We thus see that the dissipation following the relaxation from one equilibrium to the other is completely characterized by the equilibration of the initially constrained degrees of freedom.

This procedure can of course be repeated when a further species is chemostatted and it breaks another conservation law. The dissipation is quantified by a difference of semigrand Gibbs free energy, which accounts for the relaxation of the degree of freedom which has been released. When the chemostatting breaks all conservation laws without generating fundamental forces, the CRN finally reaches the global minimum of available semigrand Gibbs free energy, Fig. 5. In this case, the potential becomes the grand potential used in Ref. 19 and discussed in Sec. III D [cf. Eqs. (87), (96), (103), and (140)].

C. $W_{d}$ – $G$ gauge

The driving work and the stochastic semigrand Gibbs potential are defined up to a gauge—distinct from that involving G and W_c—which corresponds to the choice of the components. Let us consider a basis change in the space of conservation laws

ℓ_{λ} \to ℓ_{λ}^{'} = \sum_{λ^{'}} Ω_{λ λ^{'}} ℓ_{λ^{'}},

(147)

with $Ω_{λ_{u} λ_{b}} = 0$ for all λ_u, λ_b so that the unbroken ones preserve their properties. Accordingly, the conjugated intensive variables transform as

f_{λ} \to f_{λ}^{'} = \sum_{λ^{'}} f_{λ^{'}} {\bar{Ω}}_{λ^{'} λ}

(148)

[see Eq. (140)], where $\bar{Ω}$ denotes the inverse of Ω. We now notice that when the sum involves only the broken conservation laws, such a bilinear form becomes

\sum_{λ_{b}} f_{λ_{b}} ℓ_{λ_{b}} \to \sum_{λ_{b}} f_{λ_{b}} ℓ_{λ_{b}} - \sum_{λ_{u}} f_{λ_{u}} ℓ_{λ_{u}},

(149)

where

f_{λ_{u}} ≔ \sum_{λ_{u}^{'} λ_{b}^{'}} f_{λ_{b}^{'}} {\bar{Ω}}_{λ_{b}^{'} λ_{u}^{'}} Ω_{λ_{u}^{'} λ_{u}} .

(150)

Therefore, the instantaneous driving work rate [the integrand of Eq. (101) rewritten with Eq. (97)] and the semigrand potential become

{\dot{W}}_{d} (n) \to {\dot{W}}_{d} (n) + \sum_{λ_{u}} \partial_{t} f_{λ_{u}} L_{n}^{λ_{u}}

(151)

and

G (n) \to G (n) + \sum_{λ_{u}} f_{λ_{u}} L_{n}^{λ_{u}},

(152)

respectively. By contrast, the nonconservative forces—and thus the nonconservative work—is left invariant

F_{Y_{f}} \to F_{Y_{f}} + \sum_{λ_{u}} f_{λ_{u}} ℓ_{λ_{u}}^{y_{f}} = F_{Y_{f}},

(153)

since $ℓ_{λ_{u}}^{y_{f}} = 0$ ⁠. Crucially, the gauge terms in W_d and $- Δ G$ cancel and the EP is unaltered. After all, the physical process is not modified. Notice also that since the gauge term is nonfluctuating, it vanishes for cyclic protocols when integrated over a period.

We thus conclude that driving work and semigrand Gibbs potential are not univocally defined as they are affected by a gauge freedom. The gauge affecting the potential–work connection in stochastic thermodynamics led to debates; see Ref. 60 and references therein. As observed in the latter reference, the problem is rooted in what can be experimentally measured as work, as different experimental setups entail different gauge choices. In our chemical framework, different choices of the broken components involve expressions of the work in which different species appear and whose abundances need to be measured to estimate the work.

Example 7.

To illustrate the potential–work gauge, we use the CRN in Fig. 1. Let us consider the transformation of the set conservation laws, Eq. , identified by the matrix

Ω = (\begin{matrix} 1 & - 1 \\ 0 & 1 \end{matrix}),

(154)

according to which the conservation laws become

ι_{E}^{'} = ι_{E}^{} = \overset{\begin{matrix} E & E * & E ** & A & B \end{matrix}}{(\begin{matrix} 1 & 1 & 1 & 0 & 0 \end{matrix})}

(155a)

ι_{E}^{'} = ι_{b}^{} - ι_{E}^{} = \overset{\begin{matrix} E & E * & E ** & A & B \end{matrix}}{(\begin{matrix} - 1 & 0 & 0 & 1 & 1 \end{matrix})}

(155b)

Using Eq. , the gauge term reads

f_{λ_{u}} (π_{t}) = μ_{A} (π_{t})

(156)

from which we can easily derive the expression for the new driving work rate

{\dot{W}}_{d} (n) = (n_{E} - n_{A} - n_{B}) \partial_{t} μ_{A} .

(157)

The semigrand Gibbs free energy easily follows. We can now highlight the difference between the two definitions of driving work, Eqs. and : while the first entails the measurement of the population of A, B and of the activated complexes E^* and E^**, the latter entails that of A, B and of the free enzyme E. The values of the two expressions will differ except for cyclic protocols integrated over a period.

VI. FLUCTUATION THEOREMS

We now proceed to show that the driving work and the nonconservative chemical work satisfy a finite-time detailed FT. The FT holds for any process, referred to as forward, if the open CRN is initially prepared at equilibrium, Eq. (120). For the sake of simplicity, and without loss of generality, we assume that the initial distribution of unbroken components is $P ({L_{n}^{λ_{u}}}) = \prod_{λ_{u}} δ [) L_{n}^{λ_{u}}, L_{λ_{u}}])$ ⁠. Let π₀ be the initial value of the protocol, which corresponds to equilibrium ruled by $g (π_{0})$ ⁠. At time 0, the driving is activated and the CRN evolves controlled by the protocol π_τ, for τ ∈ [0, t]. The corresponding backward process is again initially prepared at the equilibrium—where $F_{Y_{f}} = 0$ —but the chemical potentials $μ_{Y_{p}}$ must have the same value they have at time t in the forward process. This guarantees that the equilibrium distribution is ruled by $g_{n} (π_{t})$ ⁠. The backward process is driven by the time-reversed protocol, $π_{τ}^{†} ≔ π_{t - τ}$ ⁠, for τ ∈ [0, t] (Fig. 6).

FIG. 6.

View large Download slide

Schematic representation of the forward and backward processes. The relaxation to the equilibrium obtained by shutting down the driving and turning off the forces at time t (respectively, 0) for the forward (respectively, backward) process merely relates the two processes, but it is irrelevant for the FT.

The finite-time detailed FT establishes the relationship between the forward and backward process

\begin{matrix} \frac{P_{t} (W_{d}, {W_{y_{f}}^{n c}})}{P_{t}^{†} (- W_{d}, {- W_{y_{f}}^{n c}})} = \exp \{β () W_{d} + \sum_{y_{f}} W_{y_{f}}^{n c} - Δ G_{e q}))\}, \end{matrix}

(158)

where $P_{t} (W_{d}, {W_{y_{f}}^{n c}})$ is the probability of observing W_d driving work and ${W_{y_{f}}^{n c}}$ nonconservative contributions along the forward process, Eqs. (101) and (99). Instead, $P_{t}^{†} (- W_{d}, {- W_{y_{f}}^{n c}})$ is the probability of observing −W_d driving work and ${- W_{y_{f}}^{n c}}$ nonconservative contributions along the backward process. Finally,

Δ G_{e q} = - k_{B} T \ln \frac{Z (π_{t}, {L_{λ_{u}}})}{Z (π_{0}, {L_{λ_{u}}})}

(159)

is the difference of the equilibrium semigrand Gibbs potential between the backward and forward initial equilibrium states. When integrating this expression over all possible values of W_d and ${W_{y_{f}}^{n c}}$ ⁠, we recover a Jarzynski-like integral FT

⟨\exp \{- β () W_{d} + \sum_{y_{f}} W_{y_{f}}^{n c}))\}⟩ = \exp \{- β Δ G_{e q}\} .

(160)

We emphasize that in contrast to the FT for the chemical work discussed in the first part of Sec. III D, the driving and nonconservative work contributions require that the process starts from the equilibrium state ruled by $G$ ⁠, which is that of open CRNs. As a consequence, there is no break of conservation laws happening during the process and $G_{e q}$ is nonfluctuating. The proof of the FT (158) is given in Appendix B and it hinges on the generating function techniques presented in Ref. 58.

We now discuss some special yet interesting cases of the FT (158). In unconditionally detailed-balance CRNs, the nonconservative work vanishes and we obtain

\frac{P_{t} (W_{d})}{P_{t}^{†} (- W_{d})} = \exp \{β () W_{d} - Δ G_{e q}))\} .

(161)

This is the analog of Crooks’ FT for CRNs^50,61 since solely the work due to external manipulations is involved. By contrast, for autonomous processes, the driving chemical work vanishes and the FT can be formulated as

\frac{P_{t} ({I_{y_{f}}})}{P_{t} ({- I_{y_{f}}})} = \exp \{β \sum_{y_{f}} F_{y_{f}} I_{y_{f}}\},

(162)

which evidences the symmetry that the fluctuations of the fundamental currents [see Eq. (79)] satisfy.

The FT in Eq. (158) is inspired by an analogous result derived in Refs. 58 and 62 in the context of generic Markov jump processes. It is a major result of this paper and its importance is manifold. It holds for processes of finite duration t and it is expressed in terms of measurable chemical quantities. Its only constraint is the initial state, which must be equilibrium. It reveals the most appropriate boundary conditions under which Jarzynski–Crooks-like FTs can be formulated for CRNs: equilibrium distribution of open CRNs. Most importantly, it evidences the merits of our stoichiometric approach based on the identification of conservation laws: it allowed us to characterize the potential describing the equilibrium distribution of open CRNs and to formulate the decomposition of the EP which supports our FTs, Eq. (102).

Remark

A physical interpretation of the argument of the exponential in Eq. (158) follows from the following observation: if, at time t, the driving is stopped and the fundamental forces (93) turned off—viz., set to zero by an appropriate choice of $μ_{Y_{f}}$ ⁠: $μ_{Y_{f}}^{*} ≔ μ_{Y_{p}} \cdot \sum_{λ_{b}} {\bar{ℓ}}_{λ_{b}}^{y_{p}} ℓ_{λ_{b}}^{y_{f}}$ —the CRN relaxes to the initial condition of the backward process. During the relaxation, neither W_d nor ${W_{y_{f}}^{n c}}$ are performed and the related EP is $T Σ_{r e l a x} = G (n, π_{t}) + k_{B} T \ln Z (π_{t}, {L_{λ_{u}}})$ ⁠. The argument of the exponential can thus be interpreted as the EP of the fictitious combined process “forward process + relaxation to the final equilibrium.”

Remark

For autonomous CRNs and arbitrary initial conditions, the steady-state FT follows

\frac{P ({{\dot{I}}_{y_{f}}})}{P ({- {\dot{I}}_{y_{f}}})} \overset{t \to \infty}{=} \exp \{t β \sum_{y_{f}} F_{y_{f}} {\dot{I}}_{y_{f}}\},

(163)

where

P ({{\dot{I}}_{y_{f}}})

is the probability of observing average rates of fundamental external currents

\{) \frac{1}{t} \int_{0}^{t} d τ I_{y_{f}} (τ) \})

equal to

{{\dot{I}}_{y_{f}}}

⁠. Equation can be proved using the large deviation technique used in Ref. 13 in combination with the local detailed balance .

A. FT along stoichiometric cycles

An alternative yet equivalent formulation of the FT (158) is that given in terms of nonconservative contributions along emergent stoichiometric cycles, Eq. (133),

\frac{P_{t} (W_{d}, {Γ_{η}})}{P_{t}^{†} (- W_{d}, {- Γ_{η}})} = \exp \{β () W_{d} + \sum_{η} Γ_{η} - Δ G_{e q}))\},

(164)

where $P_{t} (W_{d}, {Γ_{η}})$ is the probability of observing W_d driving work and {Γ_η} nonconservative contributions along the forward process. We discuss its proof in Appendix B.

Remark

As for the fundamental currents, the local detailed balance can be used to prove a steady-state FT for currents along emergent stoichiometric cycles

\frac{P ({{\dot{J}}_{η}})}{P ({- {\dot{J}}_{η}})} \overset{t \to \infty}{=} \exp \{t β \sum_{η} A_{η} {\dot{J}}_{η}\},

(165)

which is valid for autonomous CRNs and arbitrary initial conditions.

P ({{\dot{J}}_{η}})

is the probability of observing average rates of emergent cycle currents

\{) \frac{1}{t} \int_{0}^{t} d τ \sum_{ρ} ζ_{η, ρ} J_{ρ} (τ) \})

equal to

{{\dot{J}}_{η}}

⁠. In contrast to the analogous FT obtained in Ref. 14, Eq. is achieved using a stoichiometric approach based on the identification of stoichiometric cycles. For this reason, it accounts for the minimal set of nonzero macroscopic affinities.

VII. ENSEMBLE AVERAGE RATES DESCRIPTION

We now summarize our main results for ensemble average rates and discuss the relaxation to equilibrium of detailed-balanced CRNs. We also highlight the difference between an approach that does and does not take into account the topology of the CRN. We do so by recapitulating the procedure to decompose the EP into its fundamental contributions. We end by formulating a nonequilibrium Landauer’s principle.

A. Traditional description

1. Enthalpy balance

The enthalpy balance follows from the time derivative of the average enthalpy, Eq. (61),

d_{t} \sum_{n} p_{n} (h \cdot n) \equiv d_{t} ⟨ H ⟩ = ⟨ \dot{Q} ⟩ + ⟨ {\dot{W}}_{c} ⟩ .

(166)

It characterizes the average rate of change of enthalpy in the same way Eq. (69) characterizes the enthalpy change along stochastic trajectories. The average heat flow rate is given by

⟨ \dot{Q} ⟩ = ⟨ {\dot{Q}}^{t h r} ⟩ + ⟨ {\dot{Q}}^{c h m} ⟩ .

(167)

The first term quantifies the average rate of heat of reaction

⟨ {\dot{Q}}^{t h r} ⟩ = \sum_{ρ} [) h \cdot S_{ρ} + h_{Y} \cdot S_{ρ}^{Y}]) ⟨ J_{ρ} ⟩,

(168)

where ⟨J_ρ⟩ = ∑_n $w_{ρ}$ (n)p_n is the average reaction current. The second term is the average heat flow in the chemostats

⟨ {\dot{Q}}^{c h m} ⟩ = T s_{Y} \cdot ⟨ I^{Y} ⟩,

(169)

where $⟨ I^{Y} ⟩ = \sum_{ρ} (- S_{ρ}^{Y}) ⟨ J_{ρ} ⟩$ are the average external currents, Eq. (19). Instead, the ensemble average chemical work rate

⟨ {\dot{W}}_{c} ⟩ = μ_{Y} \cdot ⟨ I^{Y} ⟩

(170)

quantifies the average rate of exchange of Gibbs free energy with the chemostats.

2. Entropy production rate

At the ensemble average level, the second law of thermodynamics manifests itself in the non-negative average EP rate

\begin{matrix} ⟨ \dot{Σ} ⟩ & = & d_{t} ⟨ S ⟩ - \frac{1}{T} ⟨ \dot{Q} ⟩ \\ = & k_{B} \sum_{n, ρ} w_{ρ} (n) p_{n} \ln \frac{w_{ρ} (n) p_{n}}{w_{ρ} (n + S_{ρ}) p_{n + S_{ρ}}} \geq 0, \end{matrix}

(171)

where ⟨S⟩ = ∑_np_nS(n), Eq. (57). Using the expression for the transition affinity, Eq. (56), it can be recast into

T ⟨ \dot{Σ} ⟩ = ⟨ {\dot{W}}_{c} ⟩ - d_{t} ⟨ G ⟩,

(172)

where the chemical work rate and the average Gibbs potential are given in Eqs. (170) and (48), respectively. Equivalently, Eqs. (166), (171), and (172) can be obtained by directly averaging Eqs. (69), (73a), and (73c), respectively, over all stochastic trajectories.

For closed CRNs, Eq. (172) reduces to $d_{t} ⟨ G ⟩ = - T ⟨ \dot{Σ} ⟩ \leq 0$ ⁠. This relation, together with Eq. (49), shows that: (i) ⟨G⟩ is a Lyapunov function and hence that closed CRNs relax to equilibrium, Eq. (44); (ii) $⟨G⟩ - {⟨ G_{e q} ⟩}_{L} = T ⟨ Σ ⟩$ is the average dissipation during the relaxation to equilibrium.

B. CRN-specific description

1. Entropy production rate

We now summarize the procedure to recover the EP decomposition (102) at the ensemble average level. (i) Identify the broken and unbroken conservation laws, ${ℓ_{λ_{u}}, ℓ_{λ_{b}}}$ ⁠, Sec. II D. (ii) Identify a set of $N_{λ_{b}}$ exchanged species, y_p, for which the matrix whose rows are ${ℓ_{λ_{b}}^{y_{p}}}$ is nonsingular. The columns of its inverse are denoted by ${{\bar{ℓ}}_{λ}^{y_{p}}}$ ⁠. Physically, each species y_p breaks exactly one conservation law. The remaining exchanged species form the set denoted by y_f. (iii) The nonequilibrium semigrand Gibbs potential follows from the average of Eq. (103)

⟨ G ⟩ = \sum_{n} p_{n} [k_{B} T \ln p_{n} + g_{n}] .

(173)

It depends on the vector $⟨ M^{y_{p}} ⟩$ which describes the average population of the combination of moieties whose conservation is broken by the chemostats, Sec. II D and Eq. (91). (iv) The change in time of $⟨ g ⟩$ due to the time-dependent driving describes the average driving work rate, Eq. (101),

⟨ {\dot{W}}_{d} ⟩ = - [) \partial_{t} μ_{Y_{p}}]) \cdot ⟨ M^{y_{p}} ⟩ .

(174)

It quantifies the average amount of work spent to change the chemical potentials of the chemostats Y_p. (v) The second group of exchanged species, y_f, is used to identify the minimal set of fundamental nonconservative forces, $F_{Y_{f}} \equiv {F_{y_{f}}}$ ⁠, Eq. (93). The average nonconservative chemical work rate follows from the product of these forces and their corresponding instantaneous external currents, Eq. (67),

⟨ {\dot{W}}_{y_{f}}^{n c} ⟩ ≔ F_{y_{f}} ⟨ I_{y_{f}} ⟩ .

(175)

They quantify the average work per unit time spent to sustain a net current of species y_f across the CRN. (vi) The average EP rate decomposed as in Eq. (102) finally follows from Eqs. (173)–(175)

T ⟨ \dot{Σ} ⟩ = - d_{t} ⟨ G ⟩ + ⟨ {\dot{W}}_{d} ⟩ + \sum_{y_{f}} ⟨ {\dot{W}}_{y_{f}}^{n c} ⟩ .

(176)

Its three fundamental contributions appear: a conservative force contribution, a time-dependent driving contribution, a minimal set of nonconservative terms.

For open autonomous detailed-balanced CRNs, $F_{Y_{f}} = 0$ ⁠, $\partial_{t} μ_{Y_{p}} = 0$ ⁠, and hence Eq. (176) reduces to $d_{t} ⟨ G ⟩ = - T ⟨ \dot{Σ} ⟩ \leq 0$ ⁠. Recalling Eq. (122), this relation shows that: (i) $⟨ G ⟩$ is a Lyapunov function and hence that these CRNs relax to equilibrium, Eq. (120); (ii) $⟨G⟩ - {⟨ G_{e q} ⟩}_{L_{u}} = T ⟨ Σ ⟩$ is the average dissipation during the relaxation to equilibrium.

2. Enthalpy balance

By averaging Eq. (116), the CRN-specific average enthalpy balance also ensues

d_{t} ⟨ H ⟩ = ⟨ \dot{Q} ⟩ + ⟨ {\dot{W}}_{d} ⟩ + \sum_{y_{f}} ⟨ {\dot{W}}_{y_{f}}^{n c} ⟩,

(177)

which strengthens the interpretation of ⟨ ${\dot{W}}_{d}$ ⟩ and ${⟨ {\dot{W}}_{y_{f}}^{n c} ⟩}$ as average work rate contributions.

C. Average EP along stoichiometric cycles

The average EP decomposition expressed in terms of emergent cycle currents and affinities can be achieved through an analogous recipe. (i) Identify broken and unbroken conservation laws, ${ℓ_{λ_{u}}, ℓ_{λ_{b}}}$ ⁠, as well as stoichiometric and emergent stoichiometric cycles, {c^α, c^η} [Secs. II D and II E]. Follow steps (ii)–(iv) as above. (v) Identify the emergent stoichiometric cycles affinities, Eq. (126), as well as their corresponding average currents ∑_ρζ_η,ρ⟨J_ρ⟩, Eq. (130). (vi) The average EP rate follows from Eqs. (173) and (174) and the emergent stoichiometric cycles currents and affinities

T ⟨ \dot{Σ} ⟩ = - d_{t} ⟨ G ⟩ + ⟨ {\dot{W}}_{d} ⟩ + \sum_{η} ⟨ {\dot{Γ}}_{η} ⟩,

(178)

where

⟨ {\dot{Γ}}_{η} ⟩ = A_{η} \sum_{ρ} ζ_{η, ρ} ⟨ J_{ρ} ⟩

(179)

as in Eqs. (132) and (133).

D. Nonequilibrium Landauer’s principle

We can now formulate the nonequilibrium Landauer’s principle for the driving and nonconservative work. We have already seen that when the driving is stopped and all forces are turned off, the CRN relaxes to equilibrium by minimizing the nonequilibrium semigrand Gibbs potential. Equation (122) can be thus combined with Eq. (176), and by integrating over time, we obtain

\begin{matrix} ⟨ W_{d} ⟩ + \sum_{y_{f}} ⟨ W_{y_{f}}^{n c} ⟩ = Δ {⟨ G_{e q} ⟩}_{L_{u}} + k_{B} T Δ D (p ‖ p_{e q}) + T \underset{\geq 0}{\underset{︸}{⟨ Σ ⟩}} . \end{matrix}

(180)

This fundamental result shows that the minimal cost for transforming a CRN from an arbitrary nonequilibrium state to another is bounded by a relative entropy difference, as depicted in Fig. 7. This entropy is an information-theoretical measure of the dissimilarity between two probability distributions: the actual nonequilibrium one and its corresponding equilibrium, which is used as a reference. For processes starting at equilibrium, $k_{B} T Δ D = k_{B} T D (p_{f} ‖ p_{{e q}_{f}}) \geq 0$ quantifies the minimal cost of producing the final nonequilibrium state. By contrast, for processes relaxing to equilibrium, $k_{B} T Δ D = - k_{B} T D (p_{i} ‖ p_{{e q}_{i}}) \leq 0$ quantifies the maximum amount of work that can be extracted from the initial nonequilibrium state. For transformations in the absence of nonconservative forces (⁠ $F_{Y_{f}} = 0$ ⁠), we obtain the chemical version of the result of Ref. 28. The original Landauer’s principle⁶³ is recovered when considering erasure in a two state system (0 and 1) with identical energies. In this process, the initial equilibrium state (system equally likely to be found in 0 or 1) is transformed into a metastable nonequilibrium one (system found with probability one in 0) via a cyclic protocol. The difference of relative entropy is $Δ D = \ln 2$ and thus ⟨W_d⟩ ≥ k_BT ln 2. Finally, Kelvin’s formulation of the second law is recovered for transformation between equilibrium states in the absence of nonconservative forces, $⟨ W_{d} ⟩ \geq Δ {⟨ G_{e q} ⟩}_{L_{u}}$ ⁠.

FIG. 7.

View large Download slide

Pictorial representation of the transformation between two nonequilibrium probability distributions. The nonequilibrium transformation (blue line) is compared with the equilibrium one (green line). The latter is obtained by shutting down the driving and turning off the forces at each time (dashed gray lines).

Remark

To obtain the Landauer’s principle for ⟨

{\dot{W}}_{d}

⟩ and

{⟨ {\dot{W}}_{y_{f}}^{n c} ⟩}

⁠, the equilibrium states of the open CRN have been used as reference states; see Fig. 7. Alternatively, one could use the equilibrium states of the closed CRN, which are obtained by shutting down all exchange reactions. If one does so and uses Eq. , an analogous Landauer’s principle for the chemical work can be derived

⟨ W_{c} ⟩ = Δ {⟨ G_{e q} ⟩}_{L} + k_{B} T Δ D (p ‖ p_{e q}) + T ⟨Σ⟩ .

(181)

The traditional thermodynamic work relation

⟨ W_{c} ⟩ \geq Δ {⟨ G_{e q} ⟩}_{L}

is recovered for processes whose initial and final states are equilibrium ones.

E. Connection with deterministic descriptions

For CRNs with very abundant populations of species, a deterministic dynamical description in terms of nonlinear rate equations is justified. The corresponding nonequilibrium thermodynamics was analyzed in Ref. 25, where the counterparts of Eqs. (166), (172), and (105) can be found. Following a procedure similar to that described in this paper, one can also formulate the deterministic analog of the EP decomposition (176).

One can also recover the deterministic thermodynamic description from the ensemble average one by performing the thermodynamic limit—n ≫ 1, V ≫ 1, with n/V = : [z] finite, see Appendix A—and assuming that p_n ≃ δ_n,V[z], i.e., the distribution is very peaked around the population that is the solution of the rate equations, V[z].

We conclude with two remarks.

Remark

Not all results valid for stochastic CRNs hold for the deterministic ones. An example is provided by the adiabatic–nonadiabatic EP decomposition introduced in Ref. 64 for generic stochastic processes: it is valid for deterministic CRNs only for complex-balanced CRNs; see Refs. 25 and 65.

Remark

As briefly mentioned in Sec. II A, there is an alternative way of modeling open CRNs in which the exchanged species y are treated as particle reservoirs with very large population. All main results of our paper—i.e., the EP decomposition (102), the finite-time detailed FT (158), and the Landauer’s principle (180)—still hold. The only difference lies in the fact that the different definitions of stoichiometric matrices, Eq. (6), also entail slightly different definitions of broken conservation laws. Besides that, the procedure described in Sec. VII B can be followed in the same way.

VIII. APPLICATION

We now illustrate our EP decompositions (102) and (132) on a CRN displaying more than one fundamental force, which allows us to introduce the phenomenology of free energy transduction. We consider the following active catalytic mechanism:

\begin{matrix} T + E \overset{k_{+ 1}}{⇌} E T \overset{k_{+ 5}}{⇌} E D \overset{k_{+ 4}}{⇌} E + D, \\ E T + S \overset{k_{+ 2}}{⇌} E^{*} \overset{k_{+ 3}}{⇌} E D + P . \end{matrix}

(182)

It describes the T-driven catalysis of S into P, having D as a byproduct; see Fig. 8. All substrates and products are regarded as exchanged species

S \underset{k_{+ s}}{⇌} S_{e}, P \underset{k_{+ p}}{⇌} P_{e}, T \underset{k_{+ t}}{⇌} T_{e}, D \underset{k_{+ d}}{⇌} D_{e} .

(183)

The stoichiometric matrices $S$ and $S^{Y}$ read

\begin{matrix} E \\ ET \\ E * \\ ED \\ S \\ P \\ T \\ D \end{matrix} \overset{\begin{matrix} + 1 & + 2 & + 3 & + 4 & + 5 & + s & + p & + t & + d \end{matrix}}{(\begin{matrix} \begin{matrix} - 1 & 0 & 0 & 1 & 0 \\ 1 & - 1 & 0 & 0 & - 1 \\ 0 & 1 & - 1 & 0 & 0 \\ 0 & 0 & 1 & - 1 & 1 \\ 0 & - 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ - 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \end{matrix} & \begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix} \end{matrix})}

(184)

in which the stoichiometric matrix of the closed CRN is highlighted, and

\begin{matrix} S_{e} \\ P_{e} \\ T_{e} \\ D_{e} \end{matrix} \overset{\begin{matrix} + 1 & + 2 & + 3 & + 4 & + 5 & + s & + p & + t & + d \end{matrix}}{(\begin{matrix} \begin{matrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{matrix} & \begin{matrix} - 1 & 0 & 0 & 0 \\ 0 & - 1 & 0 & 0 \\ 0 & 0 & - 1 & 0 \\ 0 & 0 & 0 & - 1 \end{matrix} \end{matrix})}

(185)

respectively.

FIG. 8.

View large Download slide

Pictorial illustration of the open CRN in Eqs. (182) and (183), from which one can see the active catalytic mechanism more clearly.

We now follow the procedure described in Sec. VII and characterize all terms of Eq. (102). (i) The closed CRN has three independent conservation laws

ι_{E} = \overset{\begin{matrix} E & ET & E * & ED & S & P & D & T & S_{e} & P_{e} & T_{e} & D_{e} \end{matrix}}{(\begin{matrix} 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{matrix})}

(186a)

ι_{E} = \overset{\begin{matrix} E & ET & E * & ED & S & P & D & T & S_{e} & P_{e} & T_{e} & D_{e} \end{matrix}}{(\begin{matrix} 0 & 0 & 1 & 0 & 1 & 1 & 0 & 0 & 1 & 1 & 0 & 0 \end{matrix})}

(186b)

ι_{E} : \overset{\begin{matrix} E & ET & E * & ED & S & P & D & T & S_{e} & P_{e} & T_{e} & D_{e} \end{matrix}}{(\begin{matrix} 0 & 1 & 1 & 1 & 0 & 0 & 1 & 1 & 0 & 0 & 1 & 1 \end{matrix})}

(186c)

The first corresponds to the enzyme moiety and it is unbroken in the open CRN. By contrast, the last two correspond to the moieties S–P and T–D, which are broken in the open CRN. (ii) We choose S_e and T_e as chemostatted species Y_p since the entries of ℓ_S and ℓ_T corresponding to these species identify a nonsingular matrix—it is an identity matrix. (iii) The moiety population vector reads

M_{n}^{y_{p}} = \begin{matrix} \begin{matrix} S \\ T \end{matrix} & (\begin{matrix} n_{E *} + n_{s} + n_{p} \\ n_{ET} + n_{E *} + n_{ED} + n_{T} + n_{D} \end{matrix}) \end{matrix},

(187)

from which the semigrand Gibbs potential $G$ follows, Eqs. (103) and (173). (iv) The driving work rate follows from the scalar product of the vector above and

- \partial_{t} μ_{Y_{p}} = - \partial_{t} μ_{Y_{p}} = \begin{matrix} \begin{matrix} S_{e} \\ T_{e} \end{matrix} & (\begin{matrix} - \partial_{t} μ S_{e} \\ n \partial_{t} μ T_{e} \end{matrix}) \end{matrix},

(188)

Eqs. (101) and (174). (v) The chemostatted species P_e and D_e form the set Y_f and determine the fundamental forces

ℱ_{Y_{f}} = \begin{matrix} (\begin{matrix} ℱ_{P_{e}} \\ ℱ_{D_{e}} \end{matrix}) & = \begin{matrix} P_{e} \\ D_{e} \end{matrix} (\begin{matrix} μ P_{e} - μ S_{e} \\ μ D_{e} - μ T_{e} \end{matrix}) \end{matrix},

(189)

Eq. (93). Together with the instantaneous external currents

I^{Y_{f}} = (\begin{matrix} I_{P_{e}} \\ I_{D_{e}} \end{matrix}) = \begin{matrix} P_{e} \\ D_{e} \end{matrix} (\begin{matrix} J_{p +} - J_{p -} \\ J_{d +} - J_{d -} \end{matrix}),

(190)

they identify the nonconservative contributions, Eq. (99). The first one, $F_{P_{e}} I_{P_{e}}$ ⁠, characterizes the work spent to convert S into P, while the second, $F_{D_{e}} I_{D_{e}}$ ⁠, characterizes that due to the consumption of T. The sum of these terms and the driving work integrated over time contribute to the EP as in Eq. (102).

The similar EP decomposition written in terms of nonconservative contributions along stoichiometric cycles follows when these latter are identified. The kernel of stoichiometric matrix of the closed CRN is empty, while that of the open is spanned by

c_{1} = \overset{\begin{matrix} +1 & +2 & +3 & +4 & +5 & +s & +p & +t & +d \end{matrix}}{(\begin{matrix} 1 & 0 & 0 & 1 & 1 & 0 & 0 & 1 & - 1 \end{matrix})}

(191a)

c_{1} = \overset{\begin{matrix} +1 & +2 & +3 & +4 & +5 & +s & +p & +t & +d \end{matrix}}{(\begin{matrix} 1 & 1 & 1 & 1 & 0 & 1 & - 1 & 1 & - 1 \end{matrix})}

(191b)

which are regarded as emergent stoichiometric cycles. Along the first, the enzyme converts one molecule of T into one of D, while for the second it processes T and S and produces D and P

C_{1}^{Y} = \overset{S_{e} P_{e} T_{e} D_{e}}{(0 0 1 - 1)}

(192a)

C_{2}^{Y} = \overset{Se P_{e} T_{e} D_{e}}{(1 - 1 1 - 1)}

(192b)

At this point, we can proceed from step (v) and determine the affinities

A_{1} = μ_{T_{e}} - μ_{D_{e}},

(193a)

A_{2} = μ_{T_{e}} + μ_{S_{e}} - μ_{D_{e}} - μ_{P_{e}},

(193b)

as well as the related instantaneous currents

J_{1} = J_{+ p} - J_{- p} - J_{+ d} - J_{- d},

(194a)

J_{2} = J_{- p} - J_{+ p} .

(194b)

The nonconservative work follows from the products $A_{1} J_{1}$ and $A_{2} J_{2}$ and the decomposition in Eq. (132) can thus be expressed. The former characterizes the dissipation due to the futile consumption of T, since S is not converted into P. The latter, instead, is the work spent to convert T and S into D and P.

This system can be used to illustrate free energy transduction when one considers the autonomous regime where $F_{D_{e}} < 0$ ⁠, $F_{P_{e}} > 0$ ⁠, but $⟨ {\dot{W}}_{D_{e}}^{n c} ⟩ > - ⟨ {\dot{W}}_{P_{e}}^{n c} ⟩ > 0$ ⁠. Namely, the external current of P_e flows towards the chemostat, $⟨ I_{P_{e}} ⟩ < 0$ (P_e produced), despite the fact that its force is positive, $F_{P_{e}} > 0$ ⁠. This can happen thanks to the free energy provided by the conversion of T_e into D_e, $⟨ {\dot{W}}_{D_{e}}^{n c} ⟩ > 0$ ⁠. In Fig. 9, we illustrate the behavior of the average external currents and work contributions as function of time when the transducer in Fig. 8 is smoothly switched from a nontransducing regime to a transduction one. At early times, $F_{D_{e}} = 0$ ⁠, $F_{P_{e}} > 0$ ⁠, and one observes only a consumption of P_e: $⟨ I_{P_{e}} ⟩ > 0$ and $⟨ I_{D_{e}} ⟩ ≃ 0$ [respectively, orange and blue curves in Fig. 9(b)]. Consequently, the nonconservative work contributions are $⟨ {\dot{W}}_{P_{e}}^{n c} ⟩ > 0$ and $⟨ {\dot{W}}_{D_{e}}^{n c} ⟩ = 0$ [respectively, orange and blue curves in Fig. 9(b)]. By contrast, when the motive force $F_{D_{e}}$ is switched on (at large times), the current $⟨ I_{P_{e}} ⟩$ turns negative, whereas the motive current $⟨ I_{P_{e}} ⟩$ aligns itself with its corresponding force. We thus observe that $⟨ {\dot{W}}_{D_{e}}^{n c} ⟩ > - ⟨ {\dot{W}}_{P_{e}}^{n c} ⟩ > 0$ ⁠. At intermediate times, driving work is extracted following the smooth increase of the motive force [green curve in Fig. 9(b)].

FIG. 9.

View large Download slide

(a) Average external currents and (b) average work rates vs. time for the CRN in Fig. 8. The plots are obtained using 10⁴ trajectories generated via the stochastic simulation algorithm. To simplify the illustration, all substrate and products are treated as chemostatted species. The concentrations of S_e, P_e, and D_e are kept constant ([S_e] = 10, [P_e] = 70, and [D_e] = 10), whereas that of T_e increases according to a logistic function: $[T_{e}] = {[T_{e}]}_{m a x} / (1 + \exp \{- κ (t - t_{0})\})$ (⁠ ${[T_{e}]}_{m a x} = 200$ ⁠, κ = 20, t₀ = 1.5). This mimics the process in which the force that sustains the active catalysis, $F_{D_{e}}$ ⁠, is switched on from 0 to a finite value after t₀. The change of the chemical potential $μ_{T_{e}}$ is plotted in red in the inset. The choice of the rate constants is as follows: k₊₁ = 10³; k₊₂ = 10³; k₊₃ = 10³; k₊₄ = 10³; k₊₅ = 10², whereas the backward rates are obtained by means of Eq. (52) using the following values for the standard-state chemical potentials: $μ_{E} °$ = 1; $μ_{E T} °$ = 3; ${μ °}_{E^{*}} = 4$ ⁠; $μ_{E D} °$ = 2; ${μ °}_{S_{e}} = 1$ ⁠; ${μ °}_{P_{e}} = 2$ ⁠; ${μ °}_{T_{e}} = 10$ ⁠; ${μ °}_{D_{e}} = 1$ ⁠. Since reactions are unimolecular, the constant term −k_BT1 ln[s] is ignored. Finally, k_BT = 1 and the value of the enzyme moiety is L_E = 10.

IX. CONCLUSIONS AND PERSPECTIVES

In this paper, we presented a thorough description of nonequilibrium thermodynamics of stochastic CRNs. The fundamental results of traditional irreversible chemical thermodynamics (viz., enthalpy and entropy balance) are formulated at the level of single trajectories, Eqs. (61) and (72). By making use of the CRN topology and by identifying conservation laws, we decompose the EP into two fundamental work contributions and a semigrand potential difference, Eqs. (102) and (176). The driving work describes the thermodynamic cost of manipulating the CRN by changing the chemical potentials of its chemostats. Instead, the nonconservative work quantifies the cost of sustaining chemical currents through the CRN. These currents prevent the CRN from reaching equilibrium, but when the related fundamental forces vanish (and the chemical potentials of the reservoirs are kept constant in time), the CRN relaxes to equilibrium by minimizing the semigrand Gibbs potential. We elucidate the relationship between this thermodynamic potential and the dynamical potentials used in chemical reaction network theory. Our EP decomposition written in terms of stoichiometric cycle affinities generalizes previous decompositions formulated for linear CRNs or steady-state dynamics.

Two detailed FTs follow from our EP decompositions, Eqs. (158) and (164). They are valid at any time and entirely expressed in terms of physical quantities. Hence, they offer the possibility of validating our findings experimentally and, from a wider perspective, of validating the foundations of stochastic thermodynamics beyond electronic devices or colloidal particles.^66,67 Finally, we derive a nonequilibrium Landauer’s principle for the work contributions, Eq. (180), which quantifies the minimum thermodynamic cost involved in transformations between arbitrary nonequilibrium states. In contrast to early formulations of the latter principle, we consider not only the cost of external manipulations but also that related to sustained currents across the system.

Our EP decomposition identifies the fundamental dissipative contributions in CRNs of arbitrary complexity and thus it can be used to analyze free energy conversion in CRNs beyond single biocatalysts, molecular motors, or sensory systems, which are usually described by linear CRNs.^68–71 The nonconservative work contributions capture Hill’s idea of free energy transduction and extend it to nonlinear CRNs with an arbitrary number of chemical forces. (As illustrated in Sec. VIII, transduction occurs whenever one contribution becomes negative, thus requiring the other ones to be positive and larger than the former in absolute value by virtue of the second law of thermodynamics.) In turn, the driving work contribution allows us to generalize transduction to CRNs with reservoirs externally controlled in time. Hence, our framework can be used to analyze pumping in CRNs,^72,73 namely, mechanisms whose periodic external control sustains a chemical current against its spontaneous direction.

In biochemical information-handling systems^48,70,74,75 and chemical computing,⁷⁶ information is stored and processed at the molecular level. Conservation laws play a crucial role since they enable storing information in the form of nontrivial probability distributions⁷⁷ [see, e.g., Eq. (120)]. Early applications of the nonequilibrium Landauer’s principle proved successful for characterizing the thermodynamic cost of information processing in simple mechanisms.^78,79 Our generalization of this principle could thus be used to analyze biochemical information-handling systems of far greater complexity. This endeavor is important in the light of the current understanding that biological systems evolved by optimizing the gathering and representation of information.^80,81

Noise is known to play an important role in many biochemical processes. Since a complete stochastic description remains both analytically and computationally demanding, developing hybrid stochastic–deterministic descriptions would be of great importance.^25,82,83 Also, many of these processes are regulated by enzymes, thus extending the present theory beyond mass-action kinetics, as already done for deterministic CRNs,⁸⁴ is also necessary.

ACKNOWLEDGMENTS

R.R. warmly thanks A. Wachtel, A. Lazarescu, and M. Polettini for valuable discussions. This work was funded by the National Research Fund of Luxembourg (AFR PhD Grant 2014-2, No. 9114110) and the European Research Council project NanoThermo (ERC-2015-CoG Agreement No. 681456).

APPENDIX A: THERMODYNAMIC POTENTIALS

Using equilibrium statistical mechanics, we derive the equilibrium Gibbs free energy of a CRN in a given state n. Our derivation is similar to that found in Ref. 85, Sec. 3.2, whereas for different approaches, we refer the reader to Refs. 86–89.

We regard the reacting species, labeled by σ = 1, …, N_z, as solutes of an ideal dilute solution in a closed vessel. Since the solvent, s, is much more abundant than the solutes, n_s ≫ ∑_σn_σ. As in ideal solutions, interactions among solutes are negligible and the partition function of the whole solution $Q (T, n, n_{s})$ can be written as the product of single species partition functions, q ≡ {q_σ(T)} and q_s. By idealizing the solution as a lattice gas, in which each site is occupied by one molecule, we obtain

Q (T, n, n_{s}) = \frac{(n_{s} + \sum_{σ} n_{σ})!}{n_{s}! \prod_{σ} n_{σ}!} q_{s} (n_{s}) \prod_{σ} q_{σ}^{n_{σ}} .

(A1)

The combinatoric term accounts for all possible permutations of molecules, in which the overcounting due to the indistinguishability of molecules of the same species is removed. We note that the fact that different molecules might occupy different volumes is neglected.

Since we deal with dilute solutions, q ≡ {q_σ(T)} depends mainly on the temperature and the solutes–solvent interactions, whereas q_s depends on the abundance of solvent as well as the external pressure (which we omit for brevity). Using Stirling’s formula and the high relative abundance of the solvent, the combinatoric term can be approximated as

\frac{(n_{s} + \sum_{σ} n_{σ})!}{n_{s}! \prod_{σ} n_{σ}!} ≃ \prod_{σ} \frac{n_{s}^{n_{σ}}}{n_{σ}!} \equiv \frac{n_{s}^{\cdot n}}{n!} .

(A2)

Using Eq. (A1), the Gibbs free energy of a given state n is thus given by

\begin{matrix} g_{n} & = & - k_{B} T \ln Q (T, n, n_{s}) \\ = & (μ ° - 1 k_{B} T \ln n_{s}) \cdot n + k_{B} T \ln n! + g_{s}, \end{matrix}

(A3)

where

μ ° ≔ - k_{B} T \ln q

(A4)

can be identified as standard chemical potentials. Since the contribution that derives from the solvent, g_s ≔ −k_BT ln q_s(n_s), is constant, it can be set to zero without loss of generality. We emphasize that despite the idealizations that we introduced, Eq. (A3) is consistent with a rigorous approach based on mean-force potentials [cf. Ref. 87, Eq. F.44.a].

The Gibbs free energy changes along internal reactions read

\begin{aligned} Δ_{ρ_{i}} g = g_{n + S_{ρ_{i}}} - g_{n} \\ = (μ ° - 1 k_{B} T \ln n_{s}) \cdot S_{ρ_{i}} + k_{B} T \ln \frac{(n + S_{ρ_{i}})!}{n!} . \end{aligned}

(A5)

1. Thermodynamic limit

For V ≫ 1, n ≫ 1, and finite [z] = n/V, the Gibbs potential (A3) becomes

g_{n} / V ≃ μ \cdot [z] - k_{B} T [z] \cdot 1,

(A6)

where

μ = μ ° + k_{B} T \ln \{[z] / [s]\}

(A7)

are the chemical potentials of solutes in an ideal dilute solution and [s] = n_s/V is the concentration of the solvent. We thus recover the Gibbs free energy density of ideal dilute solutions; see, e.g., Refs. 51 and 90.

When applying the same limit to the Gibbs free energy differences, Eq. (A5), we recover the Gibbs free energies of reaction

Δ_{ρ_{i}} g ≃ μ \cdot S_{ρ_{i}} .

(A8)

This result also justifies the form of the second term in the local detailed balance of exchange reactions, Eq. (53).

Remark

The chemical potentials of ideal dilute solutions obtained in Eq. (A7) are expressed in terms of the concentration of the solvent. By including this term in μ° and introducing a reference concentration for each species [z₀], we recover the common expression for the potential of ideal dilute solutions $μ = \hat{μ} ° + k_{B} T \ln \{[z] / [z_{0}]\}$ ⁠, where the standard-state chemical potential $\hat{μ} ° ≔ μ ° + k_{B} T \ln \{[z_{0}] / [s]\}$ is that measured at the reference concentration.

Summarizing, g_n given in Eq. (A3) characterizes the free energy of each CRN state. In the thermodynamic limit, the traditional potentials of ideal dilute solutions are recovered.

APPENDIX B: PROOFS OF DETAILED FLUCTUATION THEOREMS

To prove the finite time detailed FTs (158), we use moment generating functions and change the notation in favor of one using brackets and operators.

Let $P_{t} (n, W_{d}, {W_{y_{f}}^{n c}})$ be the joint probability of observing a trajectory ending in the state n along which the driving work is W_d while the nonconservative contributions are ${W_{y_{f}}^{n c}}$ ⁠. These probabilities, one for each n, are stacked in the ket $| P_{t} (W_{d}, {W_{y_{f}}^{n c}}) ⟩)$ ⁠. The time evolution of their moment generating function,

\begin{matrix} | Λ_{t} (ξ_{d}, {ξ_{y_{f}}}) ⟩) & ≔ & \int d W_{d} \prod_{y_{f}} d W_{y_{f}}^{n c} \exp \{- ξ_{d} W_{d} - \sum_{y_{f}} ξ_{y_{f}} W_{y_{f}}^{n c}\} \\ \times | P_{t} (W_{d}, {W_{y_{f}}^{n c}}) ⟩), \end{matrix}

(B1)

is ruled by the biased stochastic dynamics

d_{t} | Λ_{t} (ξ_{d}, {ξ_{y_{f}}}) ⟩) = W_{t} (ξ_{d}, {ξ_{y_{f}}}) | Λ_{t} (ξ_{d}, {ξ_{y_{f}}}) ⟩),

(B2)

where the entries of the biased generator are given by

\begin{matrix} W_{m n, t} (ξ_{d}, {ξ_{y_{f}}}) & = & \sum_{ρ} w_{ρ} (n) \{\exp \{- \sum_{y_{f}} ξ_{y_{f}} F_{y_{f}} () - S_{ρ}^{y_{f}}))\}) \\ \times (δ_{m, n + S_{ρ}} - δ_{m, n}\} - ξ_{d} \partial_{t} g_{m} δ_{n, m} . \end{matrix}

(B3)

We denoted the entries of $S_{ρ}^{Y_{f}}$ as ${S_{ρ}^{y_{f}}}$ ⁠. As a consequence of the local detailed balance (94), the stochastic generator satisfies the following symmetry:

W_{t}^{T} (ξ_{d}, {ξ_{y_{f}}}) = B_{t}^{- 1} W_{t} (ξ_{d}, {1 - ξ_{y_{f}}}) B_{t},

(B4)

where the entries of $B_{t}$ are given by

B_{n m, t} ≔ \exp \{- β g_{m} (t)\} δ_{n, m} .

(B5)

Introducing the partition function for the generic equilibrium state identified by the protocol at time τ, $Z_{τ} \equiv Z (π_{τ}, {L_{λ_{u}}}) = \exp {- β G_{{e q}_{τ}}}$ ⁠, the initial condition can be written as

| Λ_{0} (ξ_{d}, {ξ_{y_{f}}}) ⟩) = | p_{{e q}_{0}} ⟩) = B_{0} / Z_{0} | 1 ⟩) .

(B6)

The ket |1⟩ refers to the vector in the state space whose entries are all equal to one.

In order to proceed further, it is convenient to first prove a preliminary result. Let us consider the generic biased dynamics, e.g., Eq. (B2),

d_{t} | Λ_{t} (ξ) ⟩) = W_{t} (ξ) | Λ_{t} (ξ) ⟩),

(B7)

whose initial condition is $|)$ Λ₀(ξ)⟩ = |p(0)⟩. A formal solution of Eq. (B7) is $| Λ_{t} (ξ) ⟩) = U_{t} (ξ) | p (0) ⟩)$ ⁠, where the time-evolution operator reads $U_{t} (ξ) = T_{+} \exp \{\int_{0}^{t} d τ W_{τ} (ξ)\}$ ⁠, with $T_{+}$ being the time-ordering operator. We clearly have $d_{t} U_{t} (ξ) = W_{t} (ξ) U_{t} (ξ)$ ⁠. Let us now consider the following transformed evolution operator:

{\tilde{U}}_{t} (ξ) ≔ X_{t}^{- 1} U_{t} (ξ) X_{0},

(B8)

with $X_{t}$ being a generic invertible operator. Its dynamics is ruled by the following biased stochastic dynamics:

\begin{matrix} d_{t} {\tilde{U}}_{t} (ξ) & = & d_{t} X_{t}^{- 1} U_{t} (ξ) X_{0} + X_{t}^{- 1} d_{t} U_{t} (ξ) X_{0} \\ = & \{d_{t} X_{t}^{- 1} X_{t} + X_{t}^{- 1} W_{t} (ξ) X_{t}\} {\tilde{U}}_{t} (ξ) \\ \equiv & {\tilde{W}}_{t} (ξ) {\tilde{U}}_{t} (ξ), \end{matrix}

(B9)

which allows us to conclude that the transformed time-evolution operator is given by

\tilde{U} (ξ) = T_{+} \exp \{\int_{0}^{t} d τ {\tilde{W}}_{τ} (ξ)\} .

(B10)

From Eqs. (B8)–(B10), we deduce that

\begin{matrix} X_{t}^{- 1} U_{t} (ξ) X_{0} = T_{+} \exp \{\int_{0}^{t} d τ [d_{τ} X_{τ}^{- 1} X_{τ} + X_{τ}^{- 1} W_{τ} (ξ) X_{τ}]\} . \end{matrix}

(B11)

We can now come back to our specific biased stochastic dynamics (B2). The moment generating function of $P_{t} (W_{d}, {W_{y_{f}}^{n c}})$ is given by

\begin{matrix} Λ_{t} (ξ_{d}, {ξ_{y_{f}}}) & = & ⟨ 1 | Λ_{t} (ξ_{d}, {ξ_{y_{f}}}) ⟩ \\ = & ⟨ 1 | U_{t} (ξ_{d}, {ξ_{y_{f}}}) B_{0} / Z_{0} | 1 ⟩ \\ = & ⟨ 1 | \frac{B_{t}}{Z_{t}} B_{t}^{- 1} U_{t} (ξ_{d}, {ξ_{y_{f}}}) B_{0} | 1 ⟩ \frac{Z_{t}}{Z_{0}}, \end{matrix}

(B12)

where $U_{t} (ξ_{d}, {ξ_{y_{f}}})$ is the time-evolution operator of the biased stochastic dynamics (B2). Note that $⟨) 1 |) B_{t} / Z_{t}$ is the equilibrium initial distribution of the backward process $⟨) p_{{e q}_{t}} |)$ ⁠. Using the relation in Eq. (B11), the last term can be rewritten as

\begin{matrix} = & ⟨ p_{{e q}_{t}} | T_{+} \exp \{\int_{0}^{t} d τ [\partial_{τ} B_{τ}^{- 1} B_{τ} + B_{τ}^{- 1} W_{τ} (ξ_{d}, {ξ_{y_{f}}}) B_{τ}]\} | 1 ⟩ \\ \times \exp \{- β Δ G_{e q}\}, \end{matrix}

(B13)

where $Δ G_{e q}$ is defined in Eq. (159). Since $\partial_{τ} B_{τ}^{- 1} B_{τ} = diag \{\partial_{τ} g_{n}\}$ ⁠, the first term in the square bracket can be added to the diagonal entries of the second term, thus giving

\begin{matrix} = & ⟨ p_{{e q}_{t}} | T_{+} \exp \{\int_{0}^{t} d τ [B_{τ}^{- 1} W_{τ} (ξ_{d} - 1, {ξ_{y_{f}}}) B_{τ}]\} | 1 ⟩ \\ \times \exp \{- β Δ G_{e q}\} . \end{matrix}

(B14)

The symmetry (B4) allow us to recast the latter into

\begin{matrix} = & ⟨ p_{{e q}_{t}} | T_{+} \exp \{\int_{0}^{t} d τ W_{τ}^{T} (ξ_{d} - 1, {1 - ξ_{y_{f}}})\} | 1 ⟩ \\ \times \exp \{- β Δ G_{e q}\} . \end{matrix}

(B15)

The crucial step comes as we transform the integration variable from τ to τ^† = t − τ. Accordingly, the time-ordering operator, $T_{+}$ ⁠, becomes an anti-time-ordering one $T_{-}$ ⁠, while the diagonal entries of the biased generator become

W_{m m, t - τ^{†}} (ξ_{d}, {ξ_{y_{f}}}) = \sum_{ρ} w_{ρ} (m, t - τ^{†}) + ξ_{d} \partial_{τ^{†}} g_{m} (t - τ^{†})

(B16)

from which we conclude that

\begin{matrix} W_{n m, t - τ^{†}} (ξ_{d}, {ξ_{y_{f}}}) & = & W_{n m, t - τ^{†}} (- ξ_{d}, {ξ_{y_{f}}}) \\ = : & W_{n m, τ^{†}}^{†} (- ξ_{d}, {ξ_{y_{f}}}) . \end{matrix}

(B17)

$W_{τ^{†}}^{†} (ξ_{d}, {ξ_{y_{f}}})$ is the biased generator of the dynamics subject to the time-reversed protocol, π^†, i.e., the dynamics of the backward process. Equation (B15) thus becomes

\begin{matrix} = & ⟨ p_{{e q}_{t}} | T_{-} \exp \{\int_{0}^{t} d τ^{†} W_{τ^{†}}^{†}^{T} (1 - ξ_{d}, {1 - ξ_{y_{f}}})\} | 1 ⟩ \\ \times \exp \{- β Δ G_{e q}\} . \end{matrix}

(B18)

Upon a global transposition, we can write

\begin{matrix} = & ⟨ 1 | T_{+} \exp \{\int_{0}^{t} d τ^{†} W_{τ^{†}}^{†} (1 - ξ_{d}, {1 - ξ_{y_{f}}})\} | p_{{e q}_{t}} ⟩ \\ \times \exp \{- β Δ G_{e q}\}, \end{matrix}

(B19)

where we also used the relationship between transposition and time-ordering

T_{+} (\prod_{i} A_{t_{l}}^{T}) = (T_{-} \prod_{i} A_{t_{l}})^{T},

(B20)

in which A_t is a generic operator. From the last expression, we readily obtain

\begin{aligned} = ⟨ 1 | U_{t}^{†} (1 - ξ_{d}, {1 - ξ_{y_{f}}}) | p_{{e q}_{t}} ⟩ \exp \{- β Δ G_{e q}\} \\ = Λ_{t}^{†} (1 - ξ_{d}, {1 - ξ_{y_{f}}}) \exp \{- β Δ G_{e q}\}, \end{aligned}

(B21)

where $Λ_{t}^{†} (ξ_{d}, {ξ_{y_{f}}})$ is the moment generating function of $P_{t}^{†} (W_{d}, {W_{y_{f}}^{n c}})$ ⁠. Summarizing, we have the following symmetry:

Λ_{t} (ξ_{d}, {ξ_{y_{f}}}) = Λ_{t}^{†} (1 - ξ_{d}, {1 - ξ_{y_{f}}}) \exp \{- β Δ G_{e q}\},

(B22)

whose inverse Laplace transform gives the FT in Eq. (158).

1. Fluctuation theorem for emergent stoichiometric cycles currents

The finite-time detailed FT for nonconservative contributions along fundamental cycles, Eq. (164), follows the same logic and mathematical steps described above. The moment generating function which now must be taken into account is

\begin{matrix} | Λ_{t} (ξ_{d}, {ξ_{η}}) ⟩) & ≔ & \int d W_{d} \prod_{η} d Γ_{η} \exp \{- ξ_{d} W_{d} - \sum_{η} ξ_{η} Γ_{η}\} \\ \times | P_{t} (W_{d}, {Γ_{η}}) ⟩), \end{matrix}

(B23)

which is ruled by the biased generator whose entries are

\begin{matrix} W_{m n, t} (ξ_{d}, {ξ_{η}}) & = & \sum_{ρ} w_{ρ} (n) \{\exp \{- \sum_{η} ξ_{η} A_{η} ζ_{η, ρ}\}) \\ \times (δ_{m, n + S_{ρ}} - δ_{m, n}\} - ξ_{d} \partial_{t} g_{m} δ_{n, m} . \end{matrix}

(B24)

The symmetry of the latter generator—on top of which the proof is constructed—is based on the expression of the local detailed balance given in Eq. (94)

W_{t}^{T} (ξ_{d}, {ξ_{η}}) = B_{t}^{- 1} W_{t} (ξ_{d}, {1 - ξ_{η}}) B_{t},

(B25)

where the entries of $B_{t}$ are given in Eq. (B5). Following the steps from Eq. (B12) to Eq. (B22), with the definitions and equations in Eqs. (B23)–(B25), proves the FT in Eq. (164).

REFERENCES

1.

T. L.

Hill

,

Free Energy Transduction in Biology

(

Academic Press

,

New York

,

1977

);

Google Scholar

T. L.

Hill

,

Free Energy Transduction and Biochemical Cycle Kinetics

(

Dover

,

2005

).

Google Scholar

2.

J.

Schnakenberg

,

Rev. Mod. Phys.

48

,

571

(

1976

).

https://doi.org/10.1103/revmodphys.48.571

Google Scholar

Crossref

3.

D. A.

McQuarrie

,

J. Appl. Probab.

4

,

413

(

1967

).

https://doi.org/10.2307/3212214

Google Scholar

Crossref

4.

D. T.

Gillespie

,

Physica A

188

,

404

(

1992

).

https://doi.org/10.1016/0378-4371(92)90283-v

Google Scholar

Crossref

5.

J.-l.

Luo

,

C.

Van den Broeck

, and

G.

Nicolis

,

Z. Phys. B: Condens. Matter

56

,

165

(

1984

).

https://doi.org/10.1007/bf01469698

Google Scholar

Crossref

6.

C. Y.

Mou

,

J.-l.

Luo

, and

G.

Nicolis

,

J. Chem. Phys.

84

,

7011

(

1986

).

https://doi.org/10.1063/1.450623

Google Scholar

Crossref

7.

Q.

Zheng

and

J.

Ross

,

J. Chem. Phys.

94

,

3644

(

1991

).

https://doi.org/10.1063/1.459735

Google Scholar

Crossref

8.

M. O.

Vlad

and

J.

Ross

,

J. Chem. Phys.

100

,

7268

(

1994

).

https://doi.org/10.1063/1.466873

Google Scholar

Crossref

9.

K.

Sekimoto

,

Stochastic Energetics

, 1st ed., Volume 799 of Lecture Notes in Physics (

Springer-Verlag Berlin Heidelberg

,

2010

).

Google Scholar

Crossref

10.

C.

Jarzynski

,

Annu. Rev. Condens. Matter Phys.

2

,

329

(

2011

).

https://doi.org/10.1146/annurev-conmatphys-062910-140506

Google Scholar

Crossref

11.

U.

Seifert

,

Rep. Prog. Phys.

75

,

126001

(

2012

).

https://doi.org/10.1088/0034-4885/75/12/126001

Google Scholar

Crossref

PubMed

12.

C.

Van den Broeck

and

M.

Esposito

,

Physica A

418

,

6

(

2015

).

https://doi.org/10.1016/j.physa.2014.04.035

Google Scholar

Crossref

13.

P.

Gaspard

,

J. Chem. Phys.

120

,

8898

(

2004

).

https://doi.org/10.1063/1.1688758

Google Scholar

Crossref

PubMed

14.

D.

Andrieux

and

P.

Gaspard

,

J. Chem. Phys.

121

,

6167

(

2004

).

https://doi.org/10.1063/1.1782391

Google Scholar

Crossref

PubMed

15.

D.

Andrieux

and

P.

Gaspard

,

J. Stat. Phys.

127

,

107

(

2007

).

https://doi.org/10.1007/s10955-006-9233-5

Google Scholar

Crossref

16.

H.

Ge

and

H.

Qian

,

Phys. Rev. Lett.

103

,

148103

(

2009

).

https://doi.org/10.1103/physrevlett.103.148103

Google Scholar

Crossref

PubMed

17.

M.

Vellela

and

H.

Qian

,

J. R. Soc., Interface

6

,

925

(

2009

).

https://doi.org/10.1098/rsif.2008.0476

Google Scholar

Crossref

PubMed

18.

H.

Qian

and

L. M.

Bishop

,

Int. J. Mol. Sci.

11

,

3472

(

2010

).

https://doi.org/10.3390/ijms11093472

Google Scholar

Crossref

PubMed

19.

T.

Schmiedl

and

U.

Seifert

,

J. Chem. Phys.

126

,

044101

(

2007

).

https://doi.org/10.1063/1.2428297

Google Scholar

Crossref

PubMed

20.

F.

Horn

and

R.

Jackson

,

Arch. Ration. Mech. Anal.

47

,

81

(

1972

).

https://doi.org/10.1007/bf00251225

Google Scholar

Crossref

21.

M.

Feinberg

,

Arch. Ration. Mech. Anal.

49

,

187

(

1972

).

https://doi.org/10.1007/bf00255665

Google Scholar

Crossref

22.

D. F.

Anderson

,

G.

Craciun

, and

T. G.

Kurtz

,

Bull. Math. Biol.

72

,

1947

(

2010

).

https://doi.org/10.1007/s11538-010-9517-4

Google Scholar

Crossref

PubMed

23.

D.

Cappelletti

and

C.

Wiuf

,

SIAM J. Appl. Math.

76

,

411

(

2016

).

https://doi.org/10.1137/15m1029916

Google Scholar

Crossref

24.

M.

Polettini

and

M.

Esposito

,

J. Chem. Phys.

141

,

024117

(

2014

).

https://doi.org/10.1063/1.4886396

Google Scholar

Crossref

PubMed

25.

R.

Rao

and

M.

Esposito

,

Phys. Rev. X

6

,

041064

(

2016

).

https://doi.org/10.1103/physrevx.6.041064

Google Scholar

Crossref

26.

M.

Polettini

,

A.

Wachtel

, and

M.

Esposito

,

J. Chem. Phys.

143

,

184103

(

2015

).

https://doi.org/10.1063/1.4935064

Google Scholar

Crossref

PubMed

27.

Y.

Murashita

,

K.

Funo

, and

M.

Ueda

,

Phys. Rev. E

90

,

042110

(

2014

).

https://doi.org/10.1103/physreve.90.042110

Google Scholar

Crossref

28.

M.

Esposito

and

C.

Van den Broeck

,

Europhys. Lett.

95

,

40004

(

2011

).

https://doi.org/10.1209/0295-5075/95/40004

Google Scholar

Crossref

29.

J. M. R.

Parrondo

,

J. M.

Horowitz

, and

T.

Sagawa

,

Nat. Phys.

11

,

131

(

2015

).

https://doi.org/10.1038/nphys3230

Google Scholar

Crossref

30.

H.

Qian

and

D. A.

Beard

,

Biophys. Chem.

114

,

213

(

2005

).

https://doi.org/10.1016/j.bpc.2004.12.001

Google Scholar

Crossref

PubMed

31.

G.

Nicolis

and

I.

Prigogine

,

Self-Organization in Nonequilibrium Systems: From Dissipative Structures to Order Through Fluctuations

(

Wiley-Blackwell

,

1977

).

Google Scholar

32.

T. G.

Kurtz

,

J. Chem. Phys.

57

,

2976

(

1972

).

https://doi.org/10.1063/1.1678692

Google Scholar

Crossref

33.

D. T.

Gillespie

,

Annu. Rev. Phys. Chem.

58

,

35

(

2007

).

https://doi.org/10.1146/annurev.physchem.58.032806.104637

Google Scholar

Crossref

PubMed

34.

R. A.

Alberty

,

Thermodynamics of Biochemical Reactions

(

Wiley-Interscience

,

2003

).

Google Scholar

Crossref

35.

H. S.

Haraldsdóttir

and

R. M. T.

Fleming

,

PLoS Comput. Biol.

12

,

e1004999

(

2016

).

https://doi.org/10.1371/journal.pcbi.1004999

Google Scholar

Crossref

PubMed

36.

D.

Kondepudi

and

I.

Prigogine

,

Modern Thermodynamics: From Heat Engines to Dissipative Structures

, 2nd ed. (

Wiley

,

2014

);

Google Scholar

Crossref

I.

Prigogine

,

Physica

15

,

272

(

1949

).

https://doi.org/10.1016/0031-8914(49)90056-7

Crossref

Google Scholar

37.

T. M.

Cover

and

J. A.

Thomas

,

Elements of Information Theory

, 2nd ed. (

Wiley-Interscience

,

2006

).

Google Scholar

38.

T.

de Donder

,

L’affinité, Mémoires de la Classe des Sciences

(

Gauthier-Villars

,

1927

), Vol. 1.

Google Scholar

39.

U.

Seifert

,

Phys. Rev. Lett.

95

,

040602

(

2005

).

https://doi.org/10.1103/physrevlett.95.040602

Google Scholar

Crossref

PubMed

40.

H.

Callen

,

Thermodynamics and an Introduction to Thermostatistics

(

Wiley

,

1985

).

Google Scholar

41.

S. R.

de Groot

and

P.

Mazur

,

Non-Equilibrium Thermodynamics

(

Dover

,

1984

).

Google Scholar

42.

G.

Bochkov

and

Y.

Kuzovlev

,

Zh. Eksp. Teor. Fiz.

72

,

238

(

1977

).

Google Scholar

43.

G.

Bochkov

and

Y.

Kuzovlev

,

Zh. Eksp. Teor. Fiz.

76

,

1071

(

1979

).

Google Scholar

44.

C.

Jarzynski

,

Phys. Rev. Lett.

78

,

2690

(

1997

).

https://doi.org/10.1103/physrevlett.78.2690

Google Scholar

Crossref

45.

J.

Horowitz

and

C.

Jarzynski

,

J. Stat. Mech.: Theory Exp.

2007

,

P11002

.

https://doi.org/10.1088/1742-5468/2007/11/p11002

Crossref

46.

P.

Krapivsky

,

S.

Redner

, and

E.

Ben-Naim

,

A Kinetic View of Statistical Physics

(

Cambridge University Press

,

2010

).

Google Scholar

Crossref

47.

R.

Rao

,

D.

Lacoste

, and

M.

Esposito

,

J. Chem. Phys.

143

,

244903

(

2015

).

https://doi.org/10.1063/1.4938009

Google Scholar

Crossref

PubMed

48.

D.

Andrieux

and

P.

Gaspard

,

Proc. Natl. Acad. Sci. U. S. A.

105

,

9516

(

2008

).

https://doi.org/10.1073/pnas.0802049105

Google Scholar

Crossref

PubMed

49.

B.

Altaner

,

J. Phys. A: Math. Theor.

50

,

454001

(

2017

).

https://doi.org/10.1088/1751-8121/aa841d

Google Scholar

Crossref

50.

G. E.

Crooks

,

J. Stat. Phys.

90

,

1481

(

1998

).

https://doi.org/10.1023/a:1023208217925

Google Scholar

Crossref

51.

L.

Peliti

,

Statistical Mechanics in a Nutshell

(

Princeton University Press

,

2011

).

Google Scholar

52.

A.

Kolmogoroff

,

Math. Ann.

112

,

155

(

1936

).

https://doi.org/10.1007/bf01565412

Google Scholar

Crossref

53.

F. P.

Kelly

,

Reversibility and Stochastic Networks

(

John Wiley & Sons Ltd.

,

1979

).

Google Scholar

54.

S.

Schuster

and

R.

Schuster

,

J. Math. Chem.

3

,

25

(

1989

).

https://doi.org/10.1007/bf01171883

Google Scholar

Crossref

55.

M.

Polettini

,

Lett. Math. Phys.

105

,

89

(

2014

).

https://doi.org/10.1007/s11005-014-0732-z

Google Scholar

Crossref

56.

D.

Andrieux

and

P.

Gaspard

,

J. Stat. Mech.: Theory Exp.

2007

,

P02006

.

https://doi.org/10.1088/1742-5468/2007/02/p02006

Crossref

57.

M.

Polettini

,

G.

Bulnes Cuetara

, and

M.

Esposito

,

Phys. Rev. E

94

,

052117

(

2016

).

https://doi.org/10.1103/physreve.94.052117

Google Scholar

Crossref

PubMed

58.

R.

Rao

and

M.

Esposito

,

New J. Phys.

20

,

023007

(

2018

).

https://doi.org/10.1088/1367-2630/aaa15f

Google Scholar

Crossref

59.

W. J.

Heuett

and

H.

Qian

,

J. Chem. Phys.

124

,

044110

(

2006

).

https://doi.org/10.1063/1.2165193

Google Scholar

Crossref

PubMed

60.

M.

Campisi

,

P.

Hänggi

, and

P.

Talkner

,

Rev. Mod. Phys.

83

,

771

(

2011

).

https://doi.org/10.1103/revmodphys.83.771

Google Scholar

Crossref

61.

G. E.

Crooks

,

Phys. Rev. E

60

,

2721

(

1999

).

https://doi.org/10.1103/physreve.60.2721

Google Scholar

Crossref

62.

G.

Bulnes Cuetara

,

M.

Esposito

, and

A.

Imparato

,

Phys. Rev. E

89

,

052119

(

2014

).

https://doi.org/10.1103/physreve.89.052119

Google Scholar

Crossref

63.

R.

Landauer

,

IBM J. Res. Dev.

44

,

261

(

2000

).

https://doi.org/10.1147/rd.441.0261

Google Scholar

Crossref

64.

M.

Esposito

,

U.

Harbola

, and

S.

Mukamel

,

Phys. Rev. E

76

,

031132

(

2007

).

https://doi.org/10.1103/physreve.76.031132

Google Scholar

Crossref

65.

H.

Ge

and

H.

Qian

,

Chem. Phys.

472

,

241

(

2016

).

https://doi.org/10.1016/j.chemphys.2016.03.026

Google Scholar

Crossref

66.

K.

Proesmans

,

Y.

Dreher

,

M.

Gavrilov

,

J.

Bechhoefer

, and

C.

Van den Broeck

,

Phys. Rev. X

6

,

041010

(

2016

).

https://doi.org/10.1103/physrevx.6.041010

Google Scholar

Crossref

67.

S.

Ciliberto

,

Phys. Rev. X

7

,

021051

(

2017

).

https://doi.org/10.1103/physrevx.7.021051

Google Scholar

Crossref

68.

U.

Seifert

,

Eur. Phys. J. E

34

,

26

(

2011

).

https://doi.org/10.1140/epje/i2011-11026-7

Google Scholar

Crossref

69.

R.

Rao

and

L.

Peliti

,

J. Stat. Mech.: Theory Exp.

2015

,

P06001

.

https://doi.org/10.1088/1742-5468/2015/06/p06001

Crossref

70.

S.

Bo

,

M.

Del Giudice

, and

A.

Celani

,

J. Stat. Mech.: Theory Exp.

2015

,

P01014

.

https://doi.org/10.1088/1742-5468/2015/01/p01014

Crossref

71.

B.

Altaner

,

A.

Wachtel

, and

J.

Vollmer

,

Phys. Rev. E

92

,

042133

(

2015

).

https://doi.org/10.1103/physreve.92.042133

Google Scholar

Crossref

72.

R. D.

Astumian

,

Annu. Rev. Biophys.

40

,

289

(

2011

).

https://doi.org/10.1146/annurev-biophys-042910-155355

Google Scholar

Crossref

PubMed

73.

M.

Esposito

and

J. M. R.

Parrondo

,

Phys. Rev. E

91

,

052114

(

2015

).

https://doi.org/10.1103/physreve.91.052114

Google Scholar

Crossref

74.

J. M.

Horowitz

and

M.

Esposito

,

Phys. Rev. X

4

,

031015

(

2014

).

https://doi.org/10.1103/physrevx.4.031015

Google Scholar

Crossref

75.

T. E.

Ouldridge

,

C. C.

Govern

, and

P. R.

ten Wolde

,

Phys. Rev. X

7

,

021004

(

2017

).

https://doi.org/10.1103/physrevx.7.021004

Google Scholar

Crossref

76.

D.

Soloveichik

,

M.

Cook

,

E.

Winfree

, and

J.

Bruck

,

Nat. Comput.

7

,

615

(

2008

).

https://doi.org/10.1007/s11047-008-9067-y

Google Scholar

Crossref

77.

W.

Poole

,

A.

Ortiz-Muñoz

,

A.

Behera

,

N. S.

Jones

,

T. E.

Ouldridge

,

E.

Winfree

, and

M.

Gopalkrishnan

, in

DNA Computing and Molecular Programming

, edited by

R.

Brijder

and

L.

Qian

(

Springer

,

Cham

,

2017

), pp.

210

–

231

.

Google Scholar

Crossref

78.

P.

Sartori

and

S.

Pigolotti

,

Phys. Rev. X

5

,

041039

(

2015

).

https://doi.org/10.1103/physrevx.5.041039

Google Scholar

Crossref

79.

T. E.

Ouldridge

and

P. R.

ten Wolde

,

Phys. Rev. Lett.

118

,

158103

(

2017

).

https://doi.org/10.1103/physrevlett.118.158103

Google Scholar

Crossref

PubMed

80.

W.

Bialek

,

Biophysics: Searching for Principles

(

Princeton University Press

,

2012

).

Google Scholar

81.

G.

Tkačik

and

W.

Bialek

,

Annu. Rev. Condens. Matter Phys.

7

(

1

),

89

(

2016

).

https://doi.org/10.1146/annurev-conmatphys-031214-014803

Google Scholar

Crossref

82.

P. C.

Bressloff

and

J. M.

Newby

,

Phys. Rev. E

89

,

042701

(

2014

).

https://doi.org/10.1103/physreve.89.042701

Google Scholar

Crossref

83.

S.

Winkelmann

and

C.

Schütte

,

J. Chem. Phys.

147

,

114115

(

2017

).

https://doi.org/10.1063/1.4986560

Google Scholar

Crossref

PubMed

84.

A.

Wachtel

,

R.

Rao

, and

M.

Esposito

,

New J. Phys.

20

,

042002

(

2018

).

https://doi.org/10.1088/1367-2630/aab5c9

Google Scholar

Crossref

85.

R.

Dirks

,

J.

Bois

,

J.

Schaeffer

,

E.

Winfree

, and

N.

Pierce

,

SIAM Rev.

49

,

65

(

2007

).

https://doi.org/10.1137/060651100

Google Scholar

Crossref

86.

D. A.

McQuarrie

,

Statistical Mechanics

(

Harper & Row

,

1976

).

Google Scholar

87.

B.

Diu

,

C.

Guthmann

,

D.

Lederer

, and

B.

Roulet

,

Éléments de Physique Statistique

(

Hermann

,

1989

).

Google Scholar

88.

D. A.

Beard

and

H.

Qian

,

Chemical Biophysics: Quantitative Analysis of Cellular Systems

(

Cambridge University Press

,

2008

).

Google Scholar

Crossref

89.

T. E.

Ouldridge

,

Nat. Comput.

17

,

3

(

2017

).

https://doi.org/10.1007/s11047-017-9646-x

Google Scholar

Crossref

PubMed

90.

E.

Fermi

,

Thermodynamics

(

Dover

,

1956

).

Google Scholar

2018

Author(s)

Conservation laws and work fluctuation relations in chemical reaction networks

I. INTRODUCTION

A. Outline

II. STOCHASTIC DYNAMICS AND CRN TOPOLOGY

A. Chemical reaction networks

B. Chemical master equation

C. Stochastic trajectories

D. Conservation laws

E. Stoichiometric cycles

III. STOCHASTIC THERMODYNAMICS

A. Equilibrium of closed CRNs

B. Local detailed balance

C. Enthalpy and entropy balance

D. FT for the chemical work and comparison with previous results

IV. CRN-SPECIFIC STOCHASTIC THERMODYNAMICS

A. Entropy production

1. Autonomous detailed-balanced CRNs

2. Unconditionally detailed-balanced CRNs

3. Autonomous CRNs

B. Energy balance

C. Equilibrium of open CRNs

D. Dissipation balance along stoichiometric cycles

V. SEMIGRAND GIBBS POTENTIAL

A. Equilibrium distributions in chemical reaction network theory

B. Hierarchies of equilibriums

C. Wd–G gauge

VI. FLUCTUATION THEOREMS

A. FT along stoichiometric cycles

VII. ENSEMBLE AVERAGE RATES DESCRIPTION

A. Traditional description

1. Enthalpy balance

2. Entropy production rate

B. CRN-specific description

1. Entropy production rate

2. Enthalpy balance

C. Average EP along stoichiometric cycles

D. Nonequilibrium Landauer’s principle

E. Connection with deterministic descriptions

VIII. APPLICATION

IX. CONCLUSIONS AND PERSPECTIVES

ACKNOWLEDGMENTS

APPENDIX A: THERMODYNAMIC POTENTIALS

1. Thermodynamic limit

APPENDIX B: PROOFS OF DETAILED FLUCTUATION THEOREMS

1. Fluctuation theorem for emergent stoichiometric cycles currents

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

C. $W_{d}$ – $G$ gauge