Thermodynamic Space of Chemical Reaction Networks

TL;DR

This work develops a universal thermodynamic framework for chemical reaction networks (CRNs) that holds for arbitrary network topologies. By combining local detailed balance with the topology of cycles via elementary flux modes and conservation laws, it derives bounds on stationary reaction affinities and on the accessible space of species concentrations, formalized as the thermodynamic space. Central to the approach are EFMs, conservation laws, and the novel chemical probe, which enable bounds that depend only on global energetic driving and network structure, not on detailed kinetics. The framework is demonstrated on paradigmatic models (e.g., Schlögl bistability, chiral symmetry breaking, self-assembly, reaction-diffusion patterns) and connected to data-driven applications, providing a general tool to predict, constrain, and design non-equilibrium chemical behavior in both natural and artificial systems. TACOS, an open-source package, implements these bounds for any CRN, facilitating practical analysis from thermodynamic properties alone.

Abstract

Living systems operate out of equilibrium, continuously consuming energy to sustain organised, functional states. Their emergent behaviour usually relies on a set of interconnected chemical reaction networks (CRNs) driven by external fluxes that keep some species at fixed concentrations. Hence, uncovering the principles governing the functioning of these CRNs is crucial to understand how living systems generate and regulate complexity. While kinetics plays a key role in shaping detailed dynamical phenomena, the range of operations of a CRN is fundamentally constrained by thermodynamics. Here, we introduce and analytically derive the "thermodynamic space" of a CRN, i.e., the range of accessible stationary concentrations that can be realized under a given energetic budget. We establish analogous bounds for reaction affinities, shedding light on how global thermodynamic properties, such as the total non-equilibrium driving, can limit local non-equilibrium quantities. We illustrate our results in various paradigmatic examples, demonstrating how the onset of complex behaviors is intimately tangled with the presence of non-equilibrium conditions. By providing a general tool for analysing CRNs, the presented framework constitutes a stepping stone to deepen our ability to predict complex out-of-equilibrium phenomena and design artificial chemical systems, starting from the sole knowledge of the underlying thermodynamic properties.

Figures (9)

• Figure 1: Thermodynamic Space of a Chemical Reaction Network (CRN). (a) Schematic representation of an open CRN where internal species transform and assemble via multiple interconnected pathways, driven by an out-of-equilibrium fuel to waste conversion due to an external chemical reservoir. Highlighted are the extremal pathways dissipating the most (red) and the least (blue) free energy. (b) Conceptual illustration of the thermodynamic space (green region) on logarithmic axes for concentrations. It shows that the accessible stationary concentrations are bounded from above by a contribution from the most dissipative pathway (red) and from below by the least dissipative (blue) one. The width of this space (vertical arrows) quantifies the energy budget via the non-equilibrium driving force.
• Figure 2: An example Chemical Reaction Network (CRN). (a) The set of chemical reactions with their associated rates and the corresponding stoichiometric matrix. These quantities define the CRN. $A$ and $B$ denote two different forms of a monomer (e.g., inactive and active), while $C$ is a dimer. $F$ and $W$ are the chemostatted external species. A chemical potential difference $\Delta\mu = \mu_F - \mu_W$ drives the conversion of $A$ into $B$, pushing the system out of equilibrium. (b) The hypergraph representation of the reaction network. Here, the number of edges pointing to each chemical species represents its stoichiometric coefficient in the reaction. (c) The conservation law $\bm{l}^X$ is derived from the stoichiometric matrix $\mathbb{S}^X$ and is orthogonal to the stoichiometric subspace of internal species. As in the main text, concentrations are indicated with lowercase letters. (d) Elementary flux modes (EFMs) involving reaction 4. An EFM represents a cyclic reaction pathway that does not result in a net change of internal species. EFMs are represented by vectors $\bm{e}$, whose entries indicate the number of times each reaction occurs within the cycle, with the sign indicated the direction in which the reaction is performed (Sec. \ref{['sec:cycle']}). The cycle affinity of an EFM, $A_{\bm{e}}$, corresponds to the net chemical potential difference associated with external species consumed or produced along the cycle.
• Figure 3: Example of a chemical probe. (a) A linear CRN is shown with all reactions indicated by black (and gray) lines and the probe indicated by a dashed red line. Arrows denote stationary reaction affinities. The stationary flux associated with the probe vanishes, but its affinity remains finite. In green, we show the EFM conformal to the affinities of the modified network (i.e., original CRN plus the probe), $e^{(F)}_{\hat{\rho}}$. We highlighted nodes and edges involved in this EFM. (b) The EFM of the modified network aligned with the probe affinity, but anti-aligned with the steady-state affinities of the original CRN, $e^{(F)}_{\hat{\rho}}$. Involved edges and nodes are highlighted.
• Figure 4: Thermodynamic space of the Schlögl model. (a) Reaction hypergraph of the Schlögl model with an internal species $X$ and two chemostatted external species $A$ and $B$. (b) The EFM of the Schlögl model, in which the internal species $X$ is created from the environment (indicated as $\emptyset$) and then degraded into it through reactions 1 and 2. (c) Nullcline of the mdoel. Fixed points (zeros of the derivative) are constrained within the thermodynamic space (green-shaded region). (d) Bifurcation diagram with cycle affinity $A_{\bm{e}} = \mu_A - \mu_B$ as the control parameter. Stable (solid lines) and unstable (dashed lines) solutions are bounded by the thermodynamic space (green-shaded region). Red and blue lines in (c) and (d) represent the lower and upper bound of the thermodynamic space, respectively, as derived from Eq. \ref{['eq:con_bound_Schlogl']}. Numerical results in (c)-(d) were obtained with $k_1^{\pm}=1$, $k_2^{\pm} = 8$ and $b = 0.02$, while the driving force was varied by adjusting the chemostatted external concentration $a$. In (c), $a = 0.8$.
• Figure 5: Thermodynamic analysis of chiral symmetry breaking. (a) Thermodynamically consistent Frank model composed of chiral species $R, S$, achiral precursor $A$, and dimer $C$. $A, C$ (green squares) are chemostatted. Dashed red edge indicate the chemical probe. (b) Key Elementary Flux Modes (EFMs) $\hat{\bm{e}}_1, \hat{\bm{e}}_2, \hat{\bm{e}}_3$ involving the probe shown with a vector representation. (c) Molecular transformation pathways, derived from the EFMs in (b), that, as a net effect, convert the species from the left to the right side of the probe ($S\to R$) (d, e) Phase portraits in the $r-s$ concentration plane, displaying nullclines ($dr/dt=0$, orange; $ds/dt=0$, blue), fixed points (stable: $\bullet$; unstable: $\circ$), and system trajectories (streamlines). The green shaded region is the thermodynamically accessible space defined by Eq. \ref{['eq:eq_constant_bound']}. In (d), $\Delta \mu = 0.7RT$, while $\Delta \mu = 2.5RT$ in (e). (f) Pitchfork bifurcation diagram for chiral symmetry breaking: Log-ratio of enantiomer concentrations, $\ln(s^{ss}/r^{ss})$, versus the thermodynamic driving force $\Delta\mu = 2\mu_A - \mu_C$. Solid lines represent stable fixed points; the dashed line indicates the unstable achiral state. The green shaded region, bounded by red and blue lines (theoretical limits at $\ln(s/r) = \pm \Delta\mu/RT$), illustrates the thermodynamic bounds on chiral imbalance. The numerical results were obtained using the parameters $k_0^+ = 1$, $k_0^- = 0.7$, $k_1^- = 16$, $k_1^- = 2$, and $m = 0.5$, while varying the concentration $a$ to control the driving force.