Chemical master equation parameter exploration using DMRG

TL;DR

This work addresses the challenge of obtaining steady-state distributions for well-mixed chemical reaction networks with large state spaces by developing a tensor-network framework based on the Doi-Peliti formalism and matrix product representations. The steady state is computed as the zero-eigenvalue right eigenvector of a non-Hermitian rate operator, using DMRG and TDVP to circumvent sampling, and the method is applied to a seven-species gene toggle switch to map bistability across a parameter grid. The authors introduce a rigorous workflow that yields a kinetic phase diagram, demonstrates that stochastic fluctuations reshape the bistable region relative to mean-field predictions, and shows how excited-state DMRG can extract relaxation timescales for switching. This tensor-network approach enables efficient parameter exploration and has potential to scale to larger CRNs, providing a seed-based strategy that reduces reliance on extensive trajectory sampling for sensitivity analyses and model discovery.

Abstract

Well-mixed chemical reaction networks (CRNs) contain many distinct chemical species with copy numbers that fluctuate in correlated ways. While those correlations are typically monitored via Monte Carlo sampling of stochastic trajectories, there is interest in systematically approximating the joint distribution over the exponentially large number of possible microstates using tensor networks or tensor trains. We exploit the tensor network strategy to determine when the steady state of a seven-species gene toggle switch CRN model supports bistability as a function of two decomposition rates, both parameters of the kinetic model. We highlight how the tensor network solution captures the effects of stochastic fluctuations, going beyond mean field and indeed deviating meaningfully from a mean-field analysis. The work furthermore develops and demonstrates several technical advances that will allow steady-states of broad classes of CRNs to be computed in a manner conducive to parameter exploration. We show that the steady-state distributions can be computed via the ordinary density matrix renormalization group (DMRG) algorithm despite having a non-Hermitian rate operator with a small spectral gap, we illustrate how that steady-state distribution can be efficiently projected to an order parameter that identifies bimodality, and we employ excited-state DMRG to calculate a relaxation timescale for the bistability.

Figures (6)

• Figure 1: Two approaches to stochastic chemical kinetics, illustrated with the gene toggle switch (GTS) model. (Top) The 14 elementary reactions of the model couple together fluctuations in the copy number of the seven chemical species. The model is especially notable for generating a bistable switch. (Middle) Traditionally, the time-evolution for those copy-number fluctuations is sampled by a kinetic Monte Carlo (Gillespie) trajectory that jumps between a set of microstates too numerous to practically enumerate. With sufficiently long trajectories, one can build converged histograms for low-dimensional order parameters like $\Delta$ to sample $P_{\rm ss}(\Delta)$. (Bottom) This manuscript describes an approach to numerically compute those same distributions without sampling. Rather, the elementary reactions are converted into a raising and lowering operator form (Doi-Peliti), allowing the chemical master equation to be expressed in terms of a tensor network operator. The steady-state distribution can be controllably approximated as a matrix product state.
• Figure 2: DMRG returns a tensor network which assigns a probability to each microstate (gray MPS). Summing over those microstates to get the marginal distribution of the order parameter $\Delta$ requires summing over all microstates that yield the particular value of $\Delta$. As $n_{\text{O}}$ does not influence $\Delta$, that coordinate is traced over (contracted with the identity vector $\mathds{1}$). Sums over the occupation numbers of the other species are efficiently realized by building a tensor network of Kronecker delta tensors. Upon contraction, one is left with an order-one tensor (brown circle) whose index ranges over the allowed values of $\Delta$ and whose elements are the marginal steady-state probabilities.
• Figure 3: A kinetic phase diagram shows the GTS model bimodality as a function of protein decomposition rates. Prior mean-field analysis identified a bistable region inside the gray wedge. DMRG calculations, repeated over a grid of $c_{\rm A}^-$ and $c_{\rm B}^-$ values yield the steady-state distribution $\ket{\pi}$ and the marginal $P_{\rm ss}(\Delta)$. Strongly bimodal $P_{\rm ss}(\Delta)$ distributions arise for kinetic parameters within the gray wedge, but stochastic fluctuations cause the bimodal region to shrink relative to the mean-field prediction. Insets highlight two distinct ways the distribution can lose bimodality. The top left inset tracks $P_{\rm ss}(\Delta)$ moving along a diagonal ($c_{\rm A}^- = c_{\rm B}^-$), showing a merging of the two peaks at high decomposition rates, akin to a second order phase transition. Note also that the relevance of the discreteness of the copies of proteins is most apparent in the case of high decomposition rates (i.e., the sea green curve). The bottom right inset tracks $P_{\rm ss}(\Delta)$ moving perpendicular to that diagonal, showing a first-order-like switch between a unimodal A-dominated distribution and a unimodal B-dominated distribution. The kinetic phase diagram is generated with fixed parameters $c_1^+ = c_1^- = c_2^+ = c_2^- = 5$ and $c_\mathrm{A}^+ = c_\mathrm{B}^+ = 1$ using DMRG convergence thresholds described in Appendix \ref{['app:fourth']}.
• Figure 4: Along the diagonal of Fig. \ref{['fig:main_results']}, $c_{\rm A}^- = c_{\rm B}^-$ and $P_{\rm ss}(\Delta)$ supports two peaks. The timescale for switches between positive and negative $\Delta$ spans orders of magnitude. That timescale is estimated from long Gillespie trajectories with or without a maximal occupation number ($M = 20$ and $M = \infty$), showing that the finite truncation at $M = 20$ only barely begins to impact the switching time for the smallest decomposition rates of $c_\mathrm{A}^- = c_\mathrm{B}^- = 0.25$. Excited-state DMRG calculations reproduce the $M = 20$ timescale via Eq. \ref{['eq:dmrgtau']}. Gillespie estimates of switching time require defining two basins with a threshold $d$, one with $\Delta \geq d$ and the other with $\Delta \leq -d$. The mean switching time is the average time between switches, averaged over $10^9$ units of time. The mean switching time is weakly impacted by the chosen threshold $d$, and an appropriately chosen threshold must shift as a function of decomposition rate (see the top left inset of Fig. \ref{['fig:main_results']}). Those shifting basin definitions create the impression that the Gillespie estimates are slightly noisy about the DMRG curve.
• Figure 5: Converging $P_{\rm ss}(\Delta)$ from Gillespie samples requires many switching events, and the Gillespie time needed to generate 100,000 switching events grows dramatically longer as the decomposition rates $c_{\rm A}^- = c_{\rm B}^-$ decrease. An advantage of DRMG is that a converged DMRG calculation can seed the calculation for new decomposition rates. Starting with a converged $\ket{\pi}$ at one value of decomposition rates, we report the number of DMRG sweeps needed to converge and find $\ket{\pi}$ for a new value of decomposition rates which have been decreased by a step size of $\delta c$. By starting at large $c_{\rm A}^- = c_{\rm B}^-$ and decreasing with a small $\delta c$, one can more cheaply construct $P_{\rm ss}(\Delta)$ down through the rare-event regime. Wall time is essentially proportional to the length of the Gillespie time and the number of DMRG sweeps. To guide the eye, an approximate conversion for the single CPU calculations is provided on the right axis, using the fact that a single DMRG sweep took about 0.4 seconds and 100,000 units of Gillespie time could be simulated in about 1 second of wall time.