Bounding the stationary distributions of the chemical master equation via mathematical programming

Abstract

The stochastic dynamics of biochemical networks are usually modelled with the chemical master equation (CME). The stationary distributions of CMEs are seldom solvable analytically, and numerical methods typically produce estimates with uncontrolled errors. Here, we introduce mathematical programming approaches that yield approximations of these distributions with computable error bounds which enable the verification of their accuracy. First, we use semidefinite programming to compute increasingly tighter upper and lower bounds on the moments of the stationary distributions for networks with rational propensities. Second, we use these moment bounds to formulate linear programs that yield convergent upper and lower bounds on the stationary distributions themselves, their marginals and stationary averages. The bounds obtained also provide a computational test for the uniqueness of the distribution. In the unique case, the bounds form an approximation of the stationary distribution with a computable bound on its error. In the non-unique case, our approach yields converging approximations of the ergodic distributions. We illustrate our methodology through several biochemical examples taken from the literature: Schlögl's model for a chemical bifurcation, a two-dimensional toggle switch, and a model for bursty gene expression.

Comprehensive Review of Models and Methods for Inferences in Bio-Chemical Reaction Networks

“…The PDE approximation can be obtained assuming Taylor expansion of CME (Schnoerr et al, 2017). The error bounds for numerically computed stationary distributions of CME are obtained in (Kuntz et al, 2017). CME for hierarchical BRNs consisting of dependent and independent sub-networks is solved analytically in (Reis et al, 2018).…”

Control Variates for Stochastic Simulation of Chemical Reaction Networks

“…For the equilibrium distribution, for example, the time-derivative of all moments is equal to zero. This directly yields constraints that have been used for parameter estimation at steady-state [5] and bounding moments of the equilibrium distribution using semi-definite programming [14,15,21]. The latter technique of bounding moments has been successfully adapted in the context of transient analysis [12,29,30].…”

“…Here, such approximations are not necessary, since we will apply the moment dynamics in the context of stochastic sampling instead of trying to integrate (6). Apart from integration strategies, setting (6) to zero has been used as a constraint for parameter estimation at steady-state [5] and bounding moments at steady-state [11,15,21]. The extension of the latter has recently lead to the adaption of these constraints to a transient setting [12,30].…”

Optimization-based synthesis of stochastic biocircuits with statistical specifications

“…In current synthetic biology, the majority of studies uses Monte Carlo based stimulations of single cell trajectories to approximately evaluate these statistics of population distributions [20]. To complement the time-consuming nature of the Monte Carlo approach, other computational tools are available to directly quantify population distributions [21,22] and raw moments [23,24,25,26] without running simulations. However, these methods are not designed to directly compute descriptive statistics of biocircuits, which makes it difficult to further develop systematic design tools that can handle statistical biocircuit specifications.…”

“…1). Building upon a moment computation approach [26,27,28], we formulate the biocircuit design problems in the form of convex optimization programs [29], which enable efficient evaluation of the statistics and its sensitivity without running time-consuming simulations of single-cell trajectories. Our optimization based synthesis approach is capable of characterizing feasible design space that satisfies multiple and possibly incompatible performance specifications with mathematical rigor.…”