Thermodynamic Efficiency in Dissipative Chemistry

TL;DR

This paper develops a thermodynamic framework for analyzing the efficiency of dissipative chemical systems operating under nonequilibrium conditions, enabling better control and optimization of artificial chemical processes.

Contribution

It introduces a theoretical approach to quantify energy, work, and dissipation in open chemical systems, advancing the understanding of their thermodynamic efficiency.

Findings

Established thermodynamic notions for dissipative chemistry.

Provided a method to quantify efficiency of nonequilibrium chemical processes.

Laid the groundwork for performance analysis of artificial chemical systems.

Abstract

Chemical processes in closed systems are poorly controllable since they always relax to equilibrium. Living systems avoid this fate and give rise to a much richer diversity of phenomena by operating under nonequilibrium conditions. Recent experiments in dissipative self-assembly also demonstrated that by opening reaction vessels and steering certain concentrations, an ocean of opportunities for artificial synthesis and energy storage emerges. To navigate it, thermodynamic notions of energy, work and dissipation must be established for these open chemical systems. Here, we do so by building upon recent theoretical advances in nonequilibrium statistical physics. As a central outcome, we show how to quantify the efficiency of such chemical operations and lay the foundation for performance analysis of any dissipative chemical process.

Figures

Tables

Table 1. Table A.1: Parameters used for the energy storage model depicted in figure  1 of the main text. Values of the backward kinetic constants were obtain through Eq. ( 9 ) in order to assure thermodynamic consistency. Note that [ W ] delimited-[] W [\mathrm{W}] is kept fixed, while we used [ F ] delimited-[] F [\mathrm{F}] to tune ℱ fuel subscript ℱ fuel \mathcal{F}_{\mathrm{fuel}} in the various discussions.

$k_{+1F}$	$5{\text{M}}^{-1}{\text{s}}^{-1}$
$k_{+1W}$	$1\cdot {10}^{-3}{\text{M}}^{-1}{\text{s}}^{-1}$
$k_{+2}$	$1{\text{M}}^{-1}{\text{s}}^{-1}$
$k_{+3F}$	$1\cdot {10}^{-6}{\text{s}}^{-1}$
$k_{+3W}$	$5{\text{s}}^{-1}$
$k_{+4}$	$1\cdot {10}^{-1}{\text{s}}^{-1}$

Equations

$${\mathcal{W}}_{\mathrm{fuel}}=\text{upDelta}\mathcal{G}+T\Sigma ,$$

$${\eta }_{\mathrm{es}}=\frac{\text{upDelta}\mathcal{G}}{{\mathcal{W}}_{\mathrm{fuel}}}=1-\frac{T\Sigma }{{\mathcal{W}}_{\mathrm{fuel}}}.$$

$$\text{upDelta}\mathcal{G}=RT\sum _{\begin{array}{c}X=\mathrm{M},\\ {\mathrm{M}}^{\ast },{\mathrm{A}}_2^{\ast },{\mathrm{A}}_2\end{array}}\left[[X]\ln \frac{[X]}{[X{]}_{\mathrm{eq}}}-[X]+[X{]}_{\mathrm{eq}}\right]\ge 0,$$

$${\dot{\mathcal{W}}}_{\mathrm{fuel}}=-{\dot{\mathcal{W}}}_{\mathrm{ext}}+T\dot{\Sigma }.$$

$${\eta }_{\mathrm{ds}}=-\frac{{\dot{\mathcal{W}}}_{\mathrm{ext}}}{{\dot{\mathcal{W}}}_{\mathrm{fuel}}}=1-\frac{T\dot{\Sigma }}{{\dot{\mathcal{W}}}_{\mathrm{fuel}}},$$

$${\mathrm{d}}_t\underset{[\mathbf{X}]}{\underbrace{\left(\begin{array}{c}[\mathrm{M}]\\ [{\mathrm{M}}^{\ast }]\\ [{\mathrm{A}}_2^{\ast }]\\ [{\mathrm{A}}_2]\end{array}\right)}}=\underset{{\mathbb{S}}^{\mathrm{X}}}{\underbrace{\left(\begin{array}{cccccc}-1& -1& 0& 0& 0& 2\\ 1& 1& -2& 0& 0& 0\\ 0& 0& 1& -1& -1& 0\\ 0& 0& 0& 1& 1& -1\end{array}\right)}}\cdot \underset{\mathbf{J}={\mathbf{J}}_{+}-{\mathbf{J}}_{-}}{\underbrace{\left(\begin{array}{c}{k}_{+1\mathrm{F}}[\mathrm{F}][\mathrm{M}]-k_{-1\mathrm{F}}[{\mathrm{M}}^{\ast }]\\ k_{+1\mathrm{W}}[\mathrm{W}][\mathrm{M}]-k_{-1\mathrm{W}}[{\mathrm{M}}^{\ast }]\\ k_{+2}[{\mathrm{M}}^{\ast }{]}^2-k_{-2}[{\mathrm{A}}_2^{\ast }]\\ k_{+3\mathrm{F}}[{\mathrm{A}}_2^{\ast }]-k_{-1\mathrm{F}}[{\mathrm{A}}_2][\mathrm{F}]\\ k_{+3\mathrm{W}}[{\mathrm{A}}_2^{\ast }]-k_{-1\mathrm{W}}[{\mathrm{A}}_2][\mathrm{W}]\\ k_{+4}[{\mathrm{A}}_2]-k_{-2}[\mathrm{M}{]}^2\end{array}\right)}},$$

$$0={\mathrm{d}}_t\underset{[\mathbf{Y}]}{\underbrace{\left(\begin{array}{c}[\mathrm{F}]\\ [\mathrm{W}]\end{array}\right)}}=\underset{{\mathbb{S}}^{\mathrm{Y}}}{\underbrace{\left(\begin{array}{cccccc}-1& 0& 0& 2& 0& 0\\ 0& -1& 0& 0& 2& 0\end{array}\right)}}\cdot \underset{\mathbf{J}={\mathbf{J}}_{+}-{\mathbf{J}}_{-}}{\underbrace{\left(\begin{array}{c}{k}_{+1\mathrm{F}}[\mathrm{F}][\mathrm{M}]-k_{-1\mathrm{F}}[{\mathrm{M}}^{\ast }]\\ k_{+1\mathrm{W}}[\mathrm{W}][\mathrm{M}]-k_{-1\mathrm{W}}[{\mathrm{M}}^{\ast }]\\ k_{+2}[{\mathrm{M}}^{\ast }{]}^2-k_{-2}[{\mathrm{A}}_2^{\ast }]\\ k_{+3\mathrm{F}}[{\mathrm{A}}_2^{\ast }]-k_{-1\mathrm{F}}[{\mathrm{A}}_2][\mathrm{F}]\\ k_{+3\mathrm{W}}[{\mathrm{A}}_2^{\ast }]-k_{-1\mathrm{W}}[{\mathrm{A}}_2][\mathrm{W}]\\ k_{+4}[{\mathrm{A}}_2]-k_{-2}[\mathrm{M}{]}^2\end{array}\right)}}+\underset{\mathbf{I}}{\underbrace{\left(\begin{array}{c}{I}_{\mathrm{F}}\\ I_{\mathrm{W}}\end{array}\right)}},$$

$$\begin{split}{\mu}_{X}&={\mu}^{\circ}_{X}+RT\ln\frac{[{X}]}{[0]}\,,&{\mu}_{Y}&={\mu}^{\circ}_{Y}+RT\ln\frac{[{Y}]}{[0]}\,,\end{split}$$

$$RT\ln \frac{k_{+\rho }}{k_{-\rho }}=-\sum _X{\mu }_X^{\circ }{\mathbb{S}}_{\rho }^X-\sum _Y{\mu }_Y^{\circ }{\mathbb{S}}_{\rho }^Y.$$

$$A_{\rho }^{\mathrm{eq}}=-{\sum }_X{\mu }_X^{\mathrm{eq}}{\mathbb{S}}_{\rho }^X-{\sum }_Y{\mu }_Y^{\mathrm{eq}}{\mathbb{S}}_{\rho }^Y=0,$$

$$J_{\rho }^{\mathrm{eq}}=J_{+\rho }^{\mathrm{eq}}-J_{-\rho }^{\mathrm{eq}}=0.$$

$$T\dot{\Sigma }=RT\sum _{\rho }{J}_{\rho }\ln \frac{J_{+\rho }}{J_{-\rho }}\ge 0.$$

$$T\dot{\Sigma }=-{\mathrm{d}}_{t}G+{\dot{\mathcal{W}}}_{\mathrm{chem}},$$

$$G={\sum }_X[X]\left({\mu }_X-RT\right)+{\sum }_Y[Y]\left({\mu }_Y-RT\right)$$

$${\dot{\mathcal{W}}}_{\mathrm{chem}}={\sum }_Y{\mu }_Y{I}_Y={\mu }_{\mathrm{F}}{I}_{\mathrm{F}}+{\mu }_{\mathrm{W}}{I}_{\mathrm{W}}$$

$${\mathbf{\ell }}_{\mathrm{M}}$$

$$L_{\mathrm{M}}={\mathbf{\ell }}_{\mathrm{L}}\cdot \left(\begin{array}{c}[\mathbf{X}]\\ [\mathbf{Y}]\end{array}\right)=[\mathrm{M}]+[{\mathrm{M}}^{\ast }]+2[{\mathrm{A}}_2^{\ast }]+2[{\mathrm{A}}_2],$$

$${\mathrm{d}}_t{L}_{\mathrm{M}}=0$$

$$L_{\mathrm{W}}={\mathbf{\ell }}_{\mathrm{W}}\cdot \left(\begin{array}{c}[\mathbf{X}]\\ [\mathbf{Y}]\end{array}\right)=[{\mathrm{M}}^{\ast }]+2[{\mathrm{A}}_2^{\ast }]+[\mathrm{F}]+[\mathrm{W}].$$

$${\mathrm{d}}_t{L}_{\mathrm{W}}:=I_{\mathrm{F}}+I_{\mathrm{W}}.$$

$$T\dot{\Sigma }=-{\mathrm{d}}_t\mathcal{G}+{\dot{\mathcal{W}}}_{\mathrm{fuel}},$$

$$\begin{array}{rl}\mathcal{G}& ={\sum }_X[X]\left({\mu }_X-RT\right)+{\sum }_Y[Y]\left({\mu }_Y-RT\right)-{\mu }_{\mathrm{W}}{L}_{\mathrm{W}}\\ =& [\mathrm{M}]{\mu }_{\mathrm{M}}+[{\mathrm{A}}_2]{\mu }_{{\mathrm{A}}_2}+[{\mathrm{M}}^{\ast }]\left({\mu }_{{\mathrm{M}}^{\ast }}-{\mu }_{\mathrm{W}}\right)+[{\mathrm{A}}_2^{\ast }]\left({\mu }_{{\mathrm{A}}_2^{\ast }}-2{\mu }_{\mathrm{W}}\right)+[\mathrm{F}]\left({\mu }_{\mathrm{F}}-{\mu }_{\mathrm{W}}\right)+\\ & -RT\left([\mathrm{M}]+[{\mathrm{A}}_2]+[{\mathrm{M}}^{\ast }]+[{\mathrm{A}}_2^{\ast }]+[\mathrm{F}]+[\mathrm{W}]\right)\end{array}$$

$${\dot{\mathcal{W}}}_{\mathrm{fuel}}:=I_{\mathrm{F}}({\mu }_{\mathrm{F}}-{\mu }_{\mathrm{W}}).$$

$$0={\frac{\mathrm{d}\Lambda }{\mathrm{d}[\mathrm{M}]}|}_{\mathrm{eq}}$$

$$0={\frac{\mathrm{d}\Lambda }{\mathrm{d}[{\mathrm{A}}_2]}|}_{\mathrm{eq}}$$

$$0={\frac{\mathrm{d}\Lambda }{\mathrm{d}[{\mathrm{M}}^{\ast }]}|}_{\mathrm{eq}}$$

$$0={\frac{\mathrm{d}\Lambda }{\mathrm{d}[{\mathrm{A}}_2^{\ast }]}|}_{\mathrm{eq}}$$

$$\begin{array}{rl}{\mathcal{G}}_{\mathrm{eq}}& =\lambda L_{\mathrm{M}}-RT\left([\mathrm{M}{]}_{\mathrm{eq}}+[{\mathrm{A}}_2{]}_{\mathrm{eq}}+[{\mathrm{M}}^{\ast }{]}_{\mathrm{eq}}+[{\mathrm{A}}_2^{\ast }{]}_{\mathrm{eq}}+[\mathrm{F}{]}_{\mathrm{eq}}+[\mathrm{W}{]}_{\mathrm{eq}}\right)\\ & =[\mathrm{M}]{\mu }_{\mathrm{M}}^{\mathrm{eq}}+[{\mathrm{A}}_2]{\mu }_{{\mathrm{A}}_2}^{\mathrm{eq}}+[{\mathrm{M}}^{\ast }]\left({\mu }_{{\mathrm{M}}^{\ast }}^{\mathrm{eq}}-{\mu }_{\mathrm{W}}\right)+[{\mathrm{A}}_2^{\ast }]\left({\mu }_{{\mathrm{A}}_2^{\ast }}^{\mathrm{eq}}-2{\mu }_{\mathrm{W}}\right)+\\ & -RT\left([\mathrm{M}{]}_{\mathrm{eq}}+[{\mathrm{A}}_2{]}_{\mathrm{eq}}+[{\mathrm{M}}^{\ast }{]}_{\mathrm{eq}}+[{\mathrm{A}}_2^{\ast }{]}_{\mathrm{eq}}+[\mathrm{F}{]}_{\mathrm{eq}}+[\mathrm{W}{]}_{\mathrm{eq}}\right),\end{array}$$

$$\mathcal{G}-{\mathcal{G}}_{\mathrm{eq}}=RT\sum _X\left[[X]\ln \frac{[X]}{[X{]}_{\mathrm{eq}}}-[X]+[X{]}_{\mathrm{eq}}\right]\ge 0.$$

$${\mathrm{d}}_t\underset{[\mathbf{X}]}{\underbrace{\left(\begin{array}{c}[\mathrm{M}]\\ [{\mathrm{M}}^{\ast }]\\ [{\mathrm{A}}_2^{\ast }]\\ [{\mathrm{A}}_2]\end{array}\right)}}=\underset{{\mathbb{S}}^{\mathrm{X}}}{\underbrace{\left(\begin{array}{cccccc}-1& -1& 0& 0& 0& 2\\ 1& 1& -2& 0& 0& 0\\ 0& 0& 1& -1& -1& 0\\ 0& 0& 0& 1& 1& -1\end{array}\right)}}\cdot \underset{\mathbf{J}={\mathbf{J}}_{+}-{\mathbf{J}}_{-}}{\underbrace{\left(\begin{array}{c}{k}_{+1\mathrm{F}}[\mathrm{F}][\mathrm{M}]-k_{-1\mathrm{F}}[{\mathrm{M}}^{\ast }]\\ k_{+1\mathrm{W}}[\mathrm{W}][\mathrm{M}]-k_{-1\mathrm{W}}[{\mathrm{M}}^{\ast }]\\ k_{+2}[{\mathrm{M}}^{\ast }{]}^2-k_{-2}[{\mathrm{A}}_2^{\ast }]\\ k_{+3\mathrm{F}}[{\mathrm{A}}_2^{\ast }]-k_{-3\mathrm{F}}[{\mathrm{A}}_2][\mathrm{F}{]}^2\\ k_{+3\mathrm{W}}[{\mathrm{A}}_2^{\ast }]-k_{-3\mathrm{W}}[{\mathrm{A}}_2][\mathrm{W}{]}^2\\ k_{+4}[{\mathrm{A}}_2]-k_{-4}[\mathrm{M}{]}^2\end{array}\right)}}+\left(\begin{array}{c}2I_{\mathrm{ext}}\\ 0\\ 0\\ -I_{\mathrm{ext}}\end{array}\right),$$

$$T\dot{\Sigma }={\dot{\mathcal{W}}}_{\mathrm{chem}},$$

$${\dot{\mathcal{W}}}_{\mathrm{chem}}={\mu }_{\mathrm{F}}{I}_{\mathrm{F}}+{\mu }_{\mathrm{W}}{I}_{\mathrm{W}}+2{\mu }_{\mathrm{M}}{I}_{\mathrm{ext}}-{\mu }_{{\mathrm{A}}_2}{I}_{\mathrm{ext}}.$$

Full text

Traditional chemical thermodynamics deals with closed systems, which always evolve towards equilibrium. At equilibrium, all reaction currents — defined as the forward reaction fluxes minus the backwards ($J_{\rho}=J_{+\rho}-J_{-\rho}$, where $\rho$ labels the reactions) — eventually vanish. The first thermodynamic description of nonequilibrium chemical processes was achieved by the Brussels school founded by de Donder and perpetuated by Prigogine Prigogine and Defay (1954); Prigogine (1967), but they focused on few reactions close to equilibrium in the so-called linear regime. However, processes such as dissipative self-assembly are open chemical reaction networks (CRN) involving many reactions operating far away from equilibrium. The openness arises from the presence of one or more chemostats, i.e. particle reservoirs coupled with the system which externally control the concentrations of some species — just like thermostats control temperatures — and allow for matter exchanges. Open CRN can then be thought of as thermodynamic machines powered by chemostats. Two operating regimes may be distinguished, reminiscent of stroke and steady-state engines. In the first, work is used to induce a time-dependent change in the species abundances that could never be reached at equilibrium. An example could be the accumulation of a large amount of molecules with a high free energy content as in dissipative self-assembly, or the depletion of some undesired species as in metabolite repair Linster et al. (2013). In the second, work is used to maintain the system in a nonequilibrium stationary state which continuously transduces an input work into useful output work. Currently no framework exists to asses how efficient and powerful such chemical engines can be. We provide one grounded in the recently established nonequilibrium thermodynamics of CRN Rao and Esposito (2016); Falasco et al. (2018), which was born from the combination of state-of-the-art statistical mechanics Jarzynski (2011); Seifert (2012); Zhang et al. (2012); Ge et al. (2012); Van den Broeck and Esposito (2015) and mathematical CRN theory Horn and Jackson (1972); Feinberg (1972). Rigorous concepts of free energy, chemical work and dissipation valid far from equilibrium reveal crucial. They provide the basis for thermodynamically meaningful definitions of efficiencies and optimal performance in the different operating regimes. In the following, energy storage (ES) and dissipative synthesis (DS) will be analyzed as models epitomizing the first and the second operating regime, respectively, but our findings apply to any dissipative chemical process.

In energy storage, an open CRN initially at equilibrium with high concentrations of low-energy molecules and low concentrations of high-energy ones, is brought out of equilibrium with the aim to increase the concentrations of the high-energy species. This process is reminiscent of charging a capacitor via the coupling to a voltage generator. In the context of supramolecular chemistry, the concept of ES was proposed by Ragazzon and Prins Ragazzon and Prins (2018). An insightful model capturing its main features is described in Fig. 1. Its thermodynamic analysis, detailed in the Supplementary Information, will now be outlined. The accumulation of the high energy species $\mathrm{A}_{2}$ is enabled when chemostats set a positive chemical potential difference of fuel and waste, i.e. $\mathcal{F}_{\mathrm{fuel}}=\mu_{\mathrm{F}}-\mu_{\mathrm{W}}>0$. This implies the injection of $\mathrm{F}$ molecules at a rate $I_{\mathrm{F}}$ and the extraction of $\mathrm{W}$ at rate $I_{\mathrm{W}}$. The resulting power (i.e., work per unit of time) performed on the system by the fueling mechanism is $\dot{\mathcal{W}}_{\mathrm{fuel}}=I_{\mathrm{F}}\,\mathcal{F}_{\mathrm{fuel}}$ Rao and Esposito (2016); Falasco et al. (2018). The proper way to quantify the energy content of an open CRN is via its nonequilibrium free energy $\mathcal{G}$. During the charging process, only part of the work, namely $\upDelta\mathcal{G}$, is dedicated to accumulate the high energy species $\mathrm{A}_{2}$ Ragazzon and Prins (2018) and is stored as free energy in the system. The remaining fraction, namely $T\Sigma$, is dissipated according to the second law of thermodynamics

[TABLE]

where $T$ is temperature and $\Sigma\geq 0$ the entropy production which only vanishes at equilibrium. The thermodynamic efficiency of the ES process is thus the ratio

[TABLE]

Eq. 1 has been used to derive the second inequality. We emphasize that each of these contributions has an explicit expression in terms of concentrations and rate constants, see Supplementary Information. For instance, the energy stored at any time with respect to equilibrium is given by the expression

[TABLE]

which can be recognized as the information measure called relative entropy Cover and Thomas (2006). Since $\upDelta\mathcal{G}$ has to be positive in ES, the second law implies that work has to be positive (done on the system). It also ensures that $\eta_{\mathrm{es}}$ is bounded between zero and one.

We simulated an ES process and plotted the dynamics of concentrations as well as efficiency and its contributions in Fig. 2. The process can be divided into a charging and a maintenance phase. During the former, the system energy grows ($\mathrm{d}_{t}\mathcal{G}>0$) in a way which correlates with the accumulation of the high energy species $\mathrm{A}_{2}$.

The process can be quite efficient as a significant portion of the work is converted into free energy. However, in the maintenance phase, the system reaches a nonequilibrium steady state. The efficiency drops towards zero (proportional to the inverse time) as the entire work is being spent to preserve the energy previously accumulated ($\mathrm{d}_{t}\mathcal{G}\simeq 0$). Figure 3 focuses on the maintenance phase for different values of $\mathcal{F}_{\mathrm{fuel}}$. It shows that by driving the system away from equilibrium, one can reach species abundances which are very different with respect to the equilibrium ones. It also shows that the accumulation of free energy does not necessarily coincide with an increase in concentration of the most energetic species $\mathrm{A}_{2}$. Indeed, while at low values of $\mathcal{F}_{\mathrm{fuel}}$ the accumulation of $\mathcal{G}$ correlates with $[\mathrm{A}_{2}]$, beyond a critical point, $\mathrm{A}_{2}$ starts to be depleted while energy continues getting stored by further moving away the concentration distribution from equilibrium.

As we have seen, the crucial part of energy storage is the charging phase, as the maintenance phase is purely dissipative and burns chemical work without any energy gain. In order to make use of the energy accumulated during the charging phase, a mechanism extracting the energetic species from the system must be introduced. This complementary but distinct working regime of an open CRN will now be considered.

In dissipative synthesis, an energetic species which accumulates thanks to a fueling process is continuously extracted from a system in a nonequilibrium steady state. One may consider for instance processes where the product either evaporates, precipitates or undergoes other fast transformations while being rapidly replaced by reactants. By building upon the above ES scheme, a simple way to model DS is to add an ideal extraction/injection mechanism to the CRN (orange dashed arrows in Fig. 1). This mechanism removes the assembled molecule $\mathrm{A}_{2}$ and renews two $\mathrm{M}$ molecules at a rate $I_{\mathrm{ext}}=k_{\mathrm{ext}}[\mathrm{A}_{2}]$. In doing so, we model the synthesis of molecules that are strongly unfavored at equilibrium, a strategy used by Nature Desai and Mitchison (1997); Howard (2001); Hess and Ross (2017) and which may be within reach of supramolecular chemists Boekhoven et al. (2010, 2015); Sorrenti et al. (2017).

From the thermodynamic standpoint detailed in the Supplementary Information, the input power spent by the fueling mechanism, $\dot{\mathcal{W}}_{\mathrm{fuel}}=I_{\mathrm{F}}\,\mathcal{F}_{\mathrm{fuel}}=I_{\mathrm{F}}\,(\mu_{\mathrm{F}}-\mu_{\mathrm{W}})$, is now not just dissipated as $T\dot{\Sigma}$, but part of it is used to sustain the production of $\mathrm{A}_{2}$:

[TABLE]

The output power released by the extraction mechanism, $\dot{\mathcal{W}}_{\mathrm{ext}}=I_{\mathrm{ext}}(2\mu_{\mathrm{M}}-\mu_{\mathrm{A_{2}}})$, is negative when DS occurs. In this context the thermodynamic efficiency is thus given by

[TABLE]

where Eq. 4 has been used to derive the second equality. It is bounded between zero and one when DS occurs.

In Fig. 4, we simulated DS for various working conditions by varying $k_{\mathrm{ext}}$ and $\mathcal{F}_{\mathrm{fuel}}$. We start our analysis by considering a given value of $\mathcal{F}_{\mathrm{fuel}}$. As $k_{\mathrm{ext}}$ is increased, $\eta_{\mathrm{ds}}$ first grows to an optimal value before decreasing towards negative values where the DS regime ends (see Fig. 4a). At the same time $I_{\mathrm{ext}}$ increases until it reaches a plateau. This happens when $k_{\mathrm{ext}}$ overcomes the ability of the system to sustain high values of $[\mathrm{A}_{2}]$ (Fig. 4b). Eventually the drop in $[\mathrm{A}_{2}]$ is such that $2\mu_{\mathrm{M}}-\mu_{\mathrm{A_{2}}}>0$, thus resulting in the loss of the DS regime. We now fix $k_{\mathrm{ext}}$ and increase $\mathcal{F}_{\mathrm{fuel}}$ (Fig. 4a and 4c). The DS regime starts at a threshold value, when $[\mathrm{A}_{2}]$ becomes high enough. After that, both $[\mathrm{A}_{2}]$ and the efficiency grow to an optimal value before decreasing again. This time however, the efficiency remains positive as $[\mathrm{M}]$ drops together with $[\mathrm{A}_{2}]$. Fig. 4c shows another important feature. As $\mathcal{F}_{\mathrm{fuel}}$ is increased, $I_{\mathrm{ext}}$ first increases too, but eventually reaches a maximum after which it decreases. This phenomenon is an hallmark of far-from-equilibrium physics which could not happen in a linear regime, namely when $k_{\mathrm{ext}}$ and $\mathcal{F}_{\mathrm{fuel}}$ are small. Remarkably, the global maximum of the efficiency in Fig. 4a is reached in a region far from equilibrium. We note that it corresponds to values of $\mathcal{F}_{\mathrm{fuel}}$ close to the one maximizing $[\mathrm{A}_{2}]$ in the maintenance phase of ES (see Fig. 3) and to values of $k_{\mathrm{ext}}$ of order one resulting in $I_{\mathrm{ext}}$ which do not overly deplete $[\mathrm{A}_{2}]$. We finally turn to the lines of maximum efficiency and efficiency at maximum power, where the maximization is done with respect to $k_{\mathrm{ext}}$ at a given $\mathcal{F}_{\mathrm{fuel}}$ (Fig. 4a). Since these two lines typically do not coincide, the study of the tradeoffs is the object of a rich field called finite-time thermodynamics Benenti et al. (2017). Interestingly, while these two lines cannot coincide in the linear regime (see Supplementary Information), we see that they do intersect far-from-equilibrium, not far from the global maximum of the efficiency. Our analysis thus allowed us to identify a region of good tradeoff between power and efficiency for the model of DS we introduced. In order to emphasize the fact that all the interesting features that we identified in DS occur far-from-equilibrium, we analyze in detail in the Supplementary Information the linear regime of DS. After identifying the Onsager matrix, we are able to reproduce the entire bottom-left part of Fig. 4a analytically as well as the behavior of the maximum efficiency and efficiency at maximum power in that region.

Thermodynamics was born from the effort to systematize the performance of steam engines. Open chemical reaction networks, which are at the core of recent efforts in artificial synthesis and ubiquitous in living systems, can be seen as chemical engines. In the spirit of this analogy, in this letter we built a chemical thermodynamic framework which enables us to systematically analyze the performance of two fundamental dissipative chemical processes. The first, energy storage, is concerned with the time dependent accumulation of high energy species far from equilibrium and is currently raising significant attention from supramolecular chemists. The second, dissipative synthesis, aims at continuously extracting the previously obtained high energy species and provides a simple and insightful instance of energy transduction. In doing so, we identified their optimal regimes of operation. Crucially they lie far-from-equilibrium in regions unreachable using conventional linear regime thermodynamics. We emphasize that the methods developed in this paper can in principle be applied to any open chemical reaction network and thus provide the basis for future performance studies and optimal design of dissipative chemistry. They are thus destined to play a major role in bioengineering and nanotechnologies.

Appendix A Energy Storage

Aa Dynamics

The evolution in time of the concentrations of the species $\mathrm{M}$, $\mathrm{M}^{\ast}$, $\mathrm{A}_{2}^{\ast}$, and $\mathrm{A}_{2}$ is ruled by the rate equations

[TABLE]

where $[\mathrm{F}]$ and $[\mathrm{W}]$ are the concentrations of fuel and waste species. Since these latter are externally kept constant by the chemostats, the balance equations for their concentrations read

[TABLE]

with $I_{\mathrm{F}}$ and $I_{\mathrm{W}}$ denoting the external currents of fuel and waste flowing from the chemostats. We denote by $X=\mathrm{M},\mathrm{M}^{\ast},\mathrm{A}_{2},\mathrm{A}_{2}^{\ast}$ the internal species, by $Y=\mathrm{F},\mathrm{W}$ the chemostatted ones, and label by $\rho=1\mathrm{F},1\mathrm{W},2,3\mathrm{F},3\mathrm{W},4$ all reactions.

Ab Thermodynamics

We imagine an isothermal, isobaric, and well-stirred ideal dilute solution. Then, each species is thermodynamically characterized by chemical potentials of the form

[TABLE]

where ${\mu}^{\circ}_{X}$ and ${\mu}^{\circ}_{Y}$ are standard-state chemical potentials and $[0]$ is the standard-state concentration.

Dynamics and thermodynamics are related via the hypothesis of local detailed balance, which relates the ratio of rate constants to the differences of standard-state chemical potentials along reactions

[TABLE]

At equilibrium, the thermodynamic forces driving each reaction, also called affinities, vanish

[TABLE]

as well as all reaction currents

[TABLE]

The dissipation of the process is captured by the entropy production (EP) rate, also vanishing at equilibrium

[TABLE]

Using the rate equations and the local detailed balance, Eq. (9), one can rewrite this quantity as

[TABLE]

where

[TABLE]

is the Gibbs free energy, while

[TABLE]

is the chemical work per unit time exchanged with the chemostats.

One can also show that if the CRN were closed (fuel and waste not chemostatted) it would relax to equilibrium by minimizing $G$ Rao and Esposito (2016). Fuel and waste are however chemostatted and we need to identify the conditions for equilibrium in the open CRN. To do so we preliminary identify the topological properties of the network.

The stoichiometric matrix $\mathbb{S}\equiv(\mathbb{S}^{\mathrm{X}},\mathbb{S}^{\mathrm{Y}})^{\mathsf{T}}$ (see Eqs. (6) and (7)) encodes the topological properties of the CRN. We can access these properties by determining its cokernel, which is spanned by

[TABLE]

The first of these vectors identifies a conserved quantity

[TABLE]

which is proved using the rate equations Eq. (6) and Eq. (7).

The second vector identifies what we call a broken conserved quantity

[TABLE]

Using again the rate equations, it can be shown that

[TABLE]

Namely, $L_{\mathrm{W}}$ changes only due to the exchange of fuel and waste with the chemostats. If the CRN were closed, $L_{\mathrm{W}}$ would be constant. Using Eq. (21), we can rewrite the entropy production in Eq. (13) as

[TABLE]

where

[TABLE]

is a semigrand Gibbs potential, and

[TABLE]

is the fueling chemical work per unit of time (i.e., the fueling power). The derivation of Eq. (24) for an arbitrary CRN is discussed in Refs. Rao and Esposito (2016); Falasco et al. (2018); Rao and Esposito (2018).

If $\mu_{\mathrm{F}}=\mu_{\mathrm{W}}$, Eq. (22) shows that $\mathcal{G}$ is a monotonically decreasing function in time, given that $T\dot{\Sigma}\geq 0$ by virtue of the second law of thermodynamics. Its minimum value — i.e., the equilibrium value — under the constraint given by the conservation law (Eq. (19)) is found by minimizing the function $\Lambda=\mathcal{G}-\lambda L_{\mathrm{M}}$, where $\lambda$ is the Lagrange multiplier corresponding to $L_{\mathrm{M}}$. The equilibrium concentrations thus satisfy the following conditions

[TABLE]

The equilibrium semigrand Gibbs potential reads

[TABLE]

which leads by direct calculation to equation 3 in the main text:

[TABLE]

Therefore, when $\mu_{\mathrm{F}}=\mu_{\mathrm{W}}$, the quantity $\mathcal{G}-\mathcal{G}_{\mathrm{eq}}$ is a Lyapunov function for the open network relaxing to equilibrium. When $\mathcal{F}_{\mathrm{fuel}}=\mu_{\mathrm{F}}-\mu_{\mathrm{W}}\neq 0$, the fueling chemical work in Eq. (22) does not vanish, and the system is prevented from reaching equilibrium.

Equation 1 in the main text is obtained by integrating Eq. (22) from time $t=0$ to a generic time $t$. In our simulation of energy storage, we focused on the special case in which the system at time $t=0$ is at equilibrium ($\mathcal{F}_{\mathrm{fuel}}=0$).

Ac Parameters

With reference to the model in figure 1 of the main text, the following parameters were used for all the simulations:

Appendix B Dissipative Synthesis

Ba Dynamics

With the addition of the extraction mechanism, the evolution in time of the concentrations of the species $\mathrm{M}$, $\mathrm{M}^{\ast}$, $\mathrm{A}_{2}$, and $\mathrm{A}_{2}^{\ast}$ is ruled by the following rate equations

[TABLE]

where the current of extraction reads $I_{\mathrm{ext}}=k_{\mathrm{ext}}[\mathrm{A}_{2}]$.

We examine this system at the steady state, in which all concentrations are stationary: $\mathrm{d}_{t}[X]_{\mathrm{ss}}=0$ for all $X$. Their expressions are not analytical, but can be easily obtained numerically, thus showing that the steady state state is unique within a broad range of values that we examined.

Bb Thermodynamics

For the dissipative synthesis model at the steady state, Eq. (13) becomes

[TABLE]

where the chemical work per unit of time now reads

[TABLE]

In order to construct the entropy balance as in equation 4 of the main text, we once again need to consider conservation vectors (18) and (18), i.e. a basis of the cokernel of $\mathbb{S}$.

[TABLE]

Now, both of these vectors identify broken conserved quantities. The former corresponds to the conserved quantity relative to the monomer

[TABLE]

In the framework of Ref. Rao and Esposito (2016), this is a broken conservation law because of the presence of the extraction mechanism. Here its value does not change by construction of the model, since every $\mathrm{A}_{2}$ which is exchanged is readily replaced by 2 $\mathrm{M}$ molecules

[TABLE]

The latter represents the $\mathrm{F}$/$\mathrm{W}$ conservation law

[TABLE]

which is broken by the fueling mechanism

[TABLE]

At the steady state all time derivative vanish, and we can use Eq. (37) to recast the chemical work per unit of time in Eq. (30) into

[TABLE]

where

[TABLE]

is the input power, and

[TABLE]

is the output power. Combining Eq. (38) with Eq. (29), we obtain equation 4 of the main text.

Bc Supplementary plots for $-\dot{\mathcal{W}}_{\mathrm{ext}}$ and $\dot{\mathcal{W}}_{\mathrm{fuel}}$

Bd Linear Regime

For $k_{\mathrm{ext}}=0$ and $\mathcal{F}_{\mathrm{fuel}}=\mu_{\mathrm{F}}-\mu_{\mathrm{W}}=0$, the entropy production at the steady state vanishes, and hence the steady state is an equilibrium one ($[X]_{\mathrm{eq}}$). For

[TABLE]

the entropy production is close to zero and hence the system is in a linear regime close to equilibrium. In this regime, we can write the steady-state concentrations as $[X]_{\mathrm{ss}}=[X]_{\mathrm{eq}}(1+f_{X}/RT)$, where $f_{X}\ll RT$ for all $X$ encode the linear shifts from equilibrium. Regarding the chemostatted ones, without loss of generality, we write $[\mathrm{F}]=[\mathrm{F}]_{\mathrm{eq}}(1+\mathcal{F}_{\mathrm{fuel}}/RT)$ and $[\mathrm{W}]=[\mathrm{W}]_{\mathrm{eq}}$, where $\mu^{\circ}_{\mathrm{F}}+RT\ln[\mathrm{F}]_{\mathrm{eq}}=\mu^{\circ}_{\mathrm{W}}+RT\ln[\mathrm{W}]_{\mathrm{eq}}$. In this way, when approximating the chemical potentials of the chemostats using the fact that $\mathcal{F}_{\mathrm{fuel}}\ll RT$, the equality in Eq. 41b is recovered.

By insterting the above expressions into the rate equations, Eq. (28) and (7), one obtains the analytical solution of the dissipative synthesis model at the steady-state in the linear regime. Indeed, by discarding second order terms and exploiting the properties of the equilibrium state ($J_{+\rho}^{\mathrm{eq}}=J_{-\rho}^{\mathrm{eq}}$), the rate equations read

[TABLE]

for the internal and chemostatted species, respectively. The extraction current is given by $I_{\mathrm{ext}}=k_{\mathrm{ext}}[\mathrm{A}_{2}]_{\mathrm{eq}}$, while the matrix $\mathbb{M}$ is a 6 by 6 matrix which encodes both the topology and the kinetics of the linear regime dynamics

[TABLE]

where

[TABLE]

are the equilibrium forward fluxes. The labels $\mathrm{X}$ and $\mathrm{F}$ in Eq. (42) select blocks of $\mathbb{M}$ corresponding to internal and fuel species, respectively, as sketched in the scheme below.

Let us now introduce the index “$\mathrm{a}$” to denote the activated species which are neither exchanged nor extracted ($\mathrm{M}^{\ast}$ and $\mathrm{A}_{2}^{\ast}$), whereas the index “$\mathrm{e}$” denotes the extracted/injected species ($\mathrm{A}_{2}$ and $\mathrm{M}$). The rate equations can thus be further split into

[TABLE]

We now observe that from the definition of conservation law, the following constraint holds

[TABLE]

which implies that

[TABLE]

This allows us to rewrite Eqs. (45) as

[TABLE]

We now solve the first of the three equations above for the vector in parenthesis, using the fact that $\mathbb{M}^{\mathrm{a}}_{\mathrm{a}}$ is nonsingular.

[TABLE]

This follows from the fact that $\mathbb{M}^{\mathrm{a}}_{\mathrm{a}}$ is Gramian Horn and Johnson (1985), and $\mathbb{S}^{\mathrm{a}}$ contains linearly independent vectors. Therefore, the last two equations in (48) can be recast into

[TABLE]

Changing signs conveniently, we can rewrite the above equations in terms of the Onsager matrix $\mathbb{L}$, which expresses the linear dependence of currents from forces when the system is close to equilibrium:

[TABLE]

Indeed, in the linear regime the chemical force associated to the extraction currents is $2\mu_{\mathrm{M}}-\mu_{\mathrm{A}_{2}}=2f_{\mathrm{M}}-f_{\mathrm{A}_{2}}$. The entries of the Onsager matrix are given by

[TABLE]

We can use Eq. (55) to analytically evaluate the efficiency $\eta_{\mathrm{ds}}$ introduced in equation 5 of the main text, as well as the output power $\dot{\mathcal{W}}_{\mathrm{ext}}$, in terms of $k_{\mathrm{ext}}$ and $\mathcal{F}_{\mathrm{fuel}}$, namely the control parameters in the model:

[TABLE]

When $\mathcal{F}_{\mathrm{fuel}}$ is kept fixed, the values of $k_{\mathrm{ext}}$ which extremise $\eta_{\mathrm{ds}}$ and $-\dot{\mathcal{W}}_{\mathrm{ext}}$ are readily found by deriving the previous expressions and look for the unique stable points:

[TABLE]

The above equations define the sets of points of maximum efficiency and efficiency at maximum power for any value of $\mathcal{F}_{\mathrm{fuel}}$ within the linear regime. By equating the right hand sides of Eq. (58) and Eq. (59), one obtains that these two expressions coicide if and only if $\mathbb{L}_{12}=\mathbb{L}_{21}=0$, which is never the case for coupled currents.

When evaluated using the parameters in SI_Table A.1, Eq. (56) reads

[TABLE]

where the cross coefficients are equal according to the Onsager reciprocal relations.

When the analytical solution is plotted against $k_{\mathrm{ext}}$ and $\mathcal{F}_{\mathrm{fuel}}$, we obtain the plot in SI_Fig. B.6b, where both maximum efficiency and efficiency at maximum power are highlighted as in figure 4 of the main text. An enlargement of the linear region of figure 4 of the main text is shown in SI_Fig. B.6a.