Thermodynamics of large-scale chemical reaction networks

Abstract

Chemical and biological networks can describe a wide variety of processes, from gene regulatory networks to biochemical oscillations. Modeled by chemical master equations, these processes are inherently stochastic, as fluctuations dominate deterministic order at mesoscopic scales. These classic many-body processes suffer from the so-called curse of high dimensionality, which makes exact mathematical descriptions exponentially expensive to compute. The exponential cost renders the study of the thermodynamic properties of such systems out of equilibrium intractable and forces approximations of system noise or assumptions of continuous particle numbers. Here, we use tensor networks to numerically explore the thermodynamics of chemical processes by directly solving the ensemble solution of the chemical master equation with efficient (sub-exponential) computational cost. We provide accurate estimates of the entropy production rate, heat flux, chemical work, and nonequilibrium thermodynamic potentials, free from sampling errors or mean-field approximations. We illustrate our results through a dissipative self-assembly model. In this way, we show how tensor networks can inform the design of efficient chemical processes in previously unattainable regimes.

Figures (6)

• Figure 1: (a) Calculating the rate of heat flux requires the change in probability $\ket{\dot{p}_t}$. As a TN, this comes from contracting the MPO $H$ with the MPS $\ket{\tilde{p}_t}$. The energy operator only acts on one the $j$th term. In practice, it is more efficiennot toot actually calculate $\ket{\dot{p}_t}$, but to build the environment tensors (grey rectangles) from the terms $k < j$ and $k > j$, as shown by the second and third terms. (b) The chemical work results from the sum of MPO, MPS contractions with the independent $\ket{\boldsymbol{1}}$ tensors, each multiplied by the chemical potential and stochiometric coefficient associated with each reaction. (c) The Renyi-$2$ entropy is advantageous due to being the quotient of two efficient tensor network contractions. The numerator is $\bra{p_t}H\ket{p_t}$, while the denominator is the inner product $\bra{p_t}\ket{p_t}$.
• Figure 2: (a) Dashed blue lines represent the entropy production rate, while the solid lines approximate the entropy production rate using $\dot{\Sigma}_2$ obtained from the Renyi-$2$ entropy. Colors correspond to volume in liters, with the largest volume being thirty times larger than the smallest, which corresponds to a state space $\approx 4e^5$ times larger. Differences in the dissipation measures stems from using the Renyi $2$-entropy rate in place $\dot{S}$. The difference between $\dot{S}$ and $\dot{S}_2$ is large for small systems (b), but decreases as the system size grows. (c) shows the decomposition of $\dot{\Sigma}$ into different energy rates for a volume $V = 6e^{-23}$L. Summing each contribution to the dissipation according to Eq. \ref{['eq:2ndlaw']} recovers the exact dissipation rate shown by the solid black line.
• Figure 3: (a) The EPR density for thirty volumes between $V=1e^{-23}$L and $V=3e^{-22}$L. For small systems $\sigma$ behaves unpredictably, but as the volume grows, each EPR density overlaps with each other, showing that $\sigma$ is becoming intensive and making $\dot\Sigma$ extensive. (b) Shows $\sigma$ at each time step as a function of volume. Again, as the volume increases the percent error is decreasing, but the progress is not constant. Different time steps show non-monotonic errors highlighted with a handful of representative error functions in black. Unlike asymptotic extensivity arguments of traditional thermodynamics, we see how the EPR is becoming extensive over the course of both the system's evolution and size.
• Figure 4: (a) $\dot{G}_2$ plotted for different volumes, indicated by the colorbar. As the volume increases, the magnitude of the free-energy rate increases. The background colors demarcate different thermodynamic functions undertaken by the system at $V=3e-22$L. Each period of functioning correlates to the expected number of molecules (b). Looking at the regions from left to right, when the system builds free-energy (green region), the system is growing the number of activated monomers M$^*$ with efficiency $\eta_{\dot{G}^+}$ (c). Increasing $\eta_{\dot{G}^+}$ means less of the chemical work is being dissipated, more is being converted to free-energy. The system briefly begins using internal free-energy, blue region as $\dot{W}_{chem} > 0$, with efficiency $\eta_{\dot{\Sigma}}$ d. For the majority of time, (red region) the system is using free-energy, $\dot{G}_2 < 0$, and doing work against the environment $\dot{W}_{chem} < 0$, which corresponds to the increase in assembled structure, A$_2$. The efficiency (e) says what percent is being dissipated versus being used externally. The final phase (blue region) the, system is harnessing chemical work $\dot{W}_{chem} > 0$ and dissipating free-energy $\dot{G}_2 < 0$ in order to maintain the number of assembled molecules. An increasing efficiency (f) says that more free-energy is being dissipated with increasing volume.
• Figure 5: (a) The percent difference in $\dot\Sigma_2$ using a bond dimension of $D=30$ and $D= 50$. We see that, as the volume grows, the difference in dissipation rates grows due to the larger system requiring more singular values to accurately capture the EPR. The modest percent difference at larger volume, despite using significantly less singular values, shows that the truncation is not removing significant singular values. (b) Conservation of probability is another vital measure of how well the TN is approximating the actual probability distribution. The solid lines illustrate how as the system size grows, eventually the BD must be increased to conserve probability. The lines with symbols show the total probability of the TN in time when $D=50$.