Solving the chemical master equation by a fast adaptive finite state projection based on the stochastic simulation algorithm

Avoiding matrix exponentials for large transition rate matrices

Exact methods for the exponentiation of matrices of dimension N can be computationally expensive in terms of execution time (N3) and memory requirements (N2), not to mention numerical precision issues. A matrix often exponentiated in the natural sciences is the rate matrix. Here, we explore five methods to exponentiate rate matrices, some of which apply more broadly to other matrix types. Three of the methods leverage a mathematical analogy between computing matrix elements of a matrix exponential process and computing transition probabilities of a dynamical process (technically a Markov jump process, MJP, typically simulated using Gillespie). In doing so, we identify a novel MJP-based method relying on restricting the number of "trajectory" jumps that incurs improved computational scaling. We then discuss this method's downstream implications on mixing properties of Monte Carlo posterior samplers. We also benchmark two other methods of matrix exponentiation valid for any matrix (beyond rate matrices and, more generally, positive definite matrices) related to solving differential equations: Runge-Kutta integrators and Krylov subspace methods. Under conditions where both the largest matrix element and the number of non-vanishing elements scale linearly with N-reasonable conditions for rate matrices often exponentiated-computational time scaling with the most competitive methods (Krylov and one of the MJP-based methods) reduces to N2 with total memory requirements of N.

Predicting DNA kinetics with a truncated continuous-time Markov chain method

Predicting the kinetics of reactions involving nucleic acid strands is a fundamental task in biology and biotechnology. Reaction kinetics can be modeled as an elementary step continuous-time Markov chain, where states correspond to secondary structures and transitions correspond to base pair formation and breakage. Since the number of states in the Markov chain could be large, rates are determined by estimating the mean first passage time from sampled trajectories. As a result, the cost of kinetic predictions becomes prohibitively expensive for rare events with extremely long trajectories. Also problematic are scenarios where multiple predictions are needed for the same reaction, e.g., under different environmental conditions, or when calibrating model parameters, because a new set of trajectories is needed multiple times. We propose a new method, called pathway elaboration, to handle these scenarios. Pathway elaboration builds a truncated continuous-time Markov chain through both biased and unbiased sampling. The resulting Markov chain has moderate state space size, so matrix methods can efficiently compute reaction rates, even for rare events. Also the transition rates of the truncated Markov chain can easily be adapted when model or environmental parameters are perturbed, making model calibration feasible. We illustrate the utility of pathway elaboration on toehold-mediated strand displacement reactions, show that it well-approximates trajectory-based predictions of unbiased elementary step models on a wide range of reaction types for which such predictions are feasible, and demonstrate that it performs better than alternative truncation-based approaches that are applicable for mean first passage time estimation. Finally, in a small study, we use pathway elaboration to optimize the Metropolis kinetic model of Multistrand, an elementary step simulator, showing that the optimized parameters greatly improve reaction rate predictions. Our framework and dataset are available at https://github.com/DNA-and-Natural-Algorithms-Group/PathwayElaboration.

An adaptive solution to the chemical master equation using quantized tensor trains with sliding windows

To cope with an extremely large or even infinite state space when solving the chemical master equation in biological problems, a potent strategy is to restrict to a finite state projection (FSP) and represent the transition matrix and probability vector in quantized tensor train (QTT) format, leading to savings in storage while retaining accuracy. In an earlier adaptive FSP-QTT algorithm, the multidimensional state space was downsized and kept in the form of a hyper rectangle that was updated when needed by selectively doubling some of its side dimensions. However, this could result in a much larger state space than necessary, with the effect of hampering both the execution time and stepping scheme. In this work, we improve the algorithm by enabling sliding windows that can dynamically slide, shrink or expand, with updates driven by a number of stochastic simulation algorithm trajectories. The ensuing state space is a considerably reduced hyper rectangle containing only the most probable states at each time step. Three numerical experiments of varying difficulty are performed to compare our approach with the original adaptive FSP-QTT algorithm.

Novel domain expansion methods to improve the computational efficiency of the Chemical Master Equation solution for large biological networks

BACKGROUND

Numerical solutions of the chemical master equation (CME) are important for understanding the stochasticity of biochemical systems. However, solving CMEs is a formidable task. This task is complicated due to the nonlinear nature of the reactions and the size of the networks which result in different realizations. Most importantly, the exponential growth of the size of the state-space, with respect to the number of different species in the system makes this a challenging assignment. When the biochemical system has a large number of variables, the CME solution becomes intractable. We introduce the intelligent state projection (ISP) method to use in the stochastic analysis of these systems. For any biochemical reaction network, it is important to capture more than one moment: this allows one to describe the system's dynamic behaviour. ISP is based on a state-space search and the data structure standards of artificial intelligence (AI). It can be used to explore and update the states of a biochemical system. To support the expansion in ISP, we also develop a Bayesian likelihood node projection (BLNP) function to predict the likelihood of the states.

RESULTS

To demonstrate the acceptability and effectiveness of our method, we apply the ISP method to several biological models discussed in prior literature. The results of our computational experiments reveal that the ISP method is effective both in terms of the speed and accuracy of the expansion, and the accuracy of the solution. This method also provides a better understanding of the state-space of the system in terms of blueprint patterns.

CONCLUSIONS

The ISP is the de-novo method which addresses both accuracy and performance problems for CME solutions. It systematically expands the projection space based on predefined inputs. This ensures accuracy in the approximation and an exact analytical solution for the time of interest. The ISP was more effective both in predicting the behavior of the state-space of the system and in performance management, which is a vital step towards modeling large biochemical systems.

Discrete and continuous models of probability flux of switching dynamics: Uncovering stochastic oscillations in a toggle-switch system

The probability flux and velocity in stochastic reaction networks can help in characterizing dynamic changes in probability landscapes of these networks. Here, we study the behavior of three different models of probability flux, namely, the discrete flux model, the Fokker-Planck model, and a new continuum model of the Liouville flux. We compare these fluxes that are formulated based on, respectively, the chemical master equation, the stochastic differential equation, and the ordinary differential equation. We examine similarities and differences among these models at the nonequilibrium steady state for the toggle switch network under different binding and unbinding conditions. Our results show that at a strong stochastic condition of weak promoter binding, continuum models of Fokker-Planck and Liouville fluxes deviate significantly from the discrete flux model. Furthermore, we report the discovery of stochastic oscillation in the toggle-switch system occurring at weak binding conditions, a phenomenon captured only by the discrete flux model.

Finite state projection for approximating the stationary solution to the chemical master equation using reaction rate equations

When modeling a physical system using a Markov chain, it is often instructive to compute its probability distribution at statistical equilibrium, thereby gaining insight into the stationary, or long-term, behavior of the system. Computing such a distribution directly is problematic when the state space of the system is large. Here, we look at the case of a chemical reaction system that models the dynamics of cellular processes, where it has become popular to constrain the computational burden by using a finite state projection, which aims only to capture the most likely states of the system, rather than every possible state. We propose an efficient method to further narrow these states to those that remain highly probable in the long run, after the transient behavior of the system has dissipated. Our strategy is to quickly estimate the local maxima of the stationary distribution using the reaction rate formulation, which is of considerably smaller size than the full-blown chemical master equation, and from there develop adaptive schemes to profile the distribution around the maxima. We include numerical tests that show the efficiency of our approach.

Bayesian Estimation for Stochastic Gene Expression Using Multifidelity Models

The finite state projection (FSP) approach to solving the chemical master equation has enabled successful inference of discrete stochastic models to predict single-cell gene regulation dynamics. Unfortunately, the FSP approach is highly computationally intensive for all but the simplest models, an issue that is highly problematic when parameter inference and uncertainty quantification takes enormous numbers of parameter evaluations. To address this issue, we propose two new computational methods for the Bayesian inference of stochastic gene expression parameters given single-cell experiments. We formulate and verify an adaptive delayed acceptance Metropolis-Hastings (ADAMH) algorithm to utilize with reduced Krylov-basis projections of the FSP. We then introduce an extension of the ADAMH into a hybrid scheme that consists of an initial phase to construct a reduced model and a faster second phase to sample from the approximate posterior distribution determined by the constructed model. We test and compare both algorithms to an adaptive Metropolis algorithm with full FSP-based likelihood evaluations on three example models and simulated data to show that the new ADAMH variants achieve substantial speedup in comparison to the full FSP approach. By reducing the computational costs of parameter estimation, we expect the ADAMH approach to enable efficient data-driven estimation for more complex gene regulation models.

An adaptive solution to the chemical master equation using tensors

Solving the chemical master equation directly is difficult due to the curse of dimensionality. We tackle that challenge by a numerical scheme based on the quantized tensor train (QTT) format, which enables us to represent the solution in a compressed form that scales linearly with the dimension. We recast the finite state projection in this QTT framework and allow it to expand adaptively based on proven error criteria. The end result is a QTT-formatted matrix exponential that we evaluate through a combination of the inexact uniformization technique and the alternating minimal energy algorithm. Our method can detect when the equilibrium distribution is reached with an inexpensive test that exploits the structure of the tensor format. We successfully perform numerical tests on high-dimensional problems that had been out of reach for classical approaches.

An application of the Krylov-FSP-SSA method to parameter fitting with maximum likelihood

Monte Carlo methods such as the stochastic simulation algorithm (SSA) have traditionally been employed in gene regulation problems. However, there has been increasing interest to directly obtain the probability distribution of the molecules involved by solving the chemical master equation (CME). This requires addressing the curse of dimensionality that is inherent in most gene regulation problems. The finite state projection (FSP) seeks to address the challenge and there have been variants that further reduce the size of the projection or that accelerate the resulting matrix exponential. The Krylov-FSP-SSA variant has proved numerically efficient by combining, on one hand, the SSA to adaptively drive the FSP, and on the other hand, adaptive Krylov techniques to evaluate the matrix exponential. Here we apply this Krylov-FSP-SSA to a mutual inhibitory gene network synthetically engineered in Saccharomyces cerevisiae, in which bimodality arises. We show numerically that the approach can efficiently approximate the transient probability distribution, and this has important implications for parameter fitting, where the CME has to be solved for many different parameter sets. The fitting scheme amounts to an optimization problem of finding the parameter set so that the transient probability distributions fit the observations with maximum likelihood. We compare five optimization schemes for this difficult problem, thereby providing further insights into this approach of parameter estimation that is often applied to models in systems biology where there is a need to calibrate free parameters.

Markov State Models of gene regulatory networks

Background

Gene regulatory networks with dynamics characterized by multiple stable states underlie cell fate-decisions. Quantitative models that can link molecular-level knowledge of gene regulation to a global understanding of network dynamics have the potential to guide cell-reprogramming strategies. Networks are often modeled by the stochastic Chemical Master Equation, but methods for systematic identification of key properties of the global dynamics are currently lacking.

Results

The method identifies the number, phenotypes, and lifetimes of long-lived states for a set of common gene regulatory network models. Application of transition path theory to the constructed Markov State Model decomposes global dynamics into a set of dominant transition paths and associated relative probabilities for stochastic state-switching.

Conclusions

In this proof-of-concept study, we found that the Markov State Model provides a general framework for analyzing and visualizing stochastic multistability and state-transitions in gene networks. Our results suggest that this framework—adopted from the field of atomistic Molecular Dynamics—can be a useful tool for quantitative Systems Biology at the network scale.

Electronic supplementary material

The online version of this article (doi:10.1186/s12918-017-0394-4) contains supplementary material, which is available to authorized users.

Understanding the finite state projection and related methods for solving the chemical master equation

The finite state projection (FSP) method has enabled us to solve the chemical master equation of some biological models that were considered out of reach not long ago. Since the original FSP method, much effort has gone into transforming it into an adaptive time-stepping algorithm as well as studying its accuracy. Some of the improvements include the multiple time interval FSP, the sliding windows, and most notably the Krylov-FSP approach. Our goal in this tutorial is to give the reader an overview of the current methods that build on the FSP.

Computational study of p53 regulation via the chemical master equation

A stochastic model of cellular p53 regulation was established in Leenders, and Tuszynski (2013 Front. Oncol. 3 1-16) to study the interactions of p53 with MDM2 proteins, where the stochastic analysis was done using a Monte Carlo approach. We revisit that model here using an alternative scheme, which is to directly solve the chemical master equation (CME) by an adaptive Krylov-based finite state projection method that combines the stochastic simulation algorithm with other computational strategies, namely Krylov approximation techniques to the matrix exponential, divide and conquer, and aggregation. We report numerical results that demonstrate the extend of tackling the CME with this combination of tools.