Solving the chemical master equation using sliding windows

The chemical master equation (CME) is a system of ordinary differential equations that describes the evolution of a network of chemical reactions as a stochastic process. Its solution yields the probability density vector of the system at each point in time. Solving the CME numerically is in many cases computationally expensive or even infeasible as the number of reachable states can be very large or infinite. We introduce the sliding window method, which computes an approximate solution of the CME by performing a sequence of local analysis steps. In each step, only a manageable subset of states is considered, representing a ``window'' into the state space. In subsequent steps, the window follows the direction in which the probability mass moves, until the time period of interest has elapsed. We construct the window based on a deterministic approximation of the future behavior of the system by estimating upper and lower bounds on the populations of the chemical species. In order to show the effectiveness of our approach, we apply it to several examples previously described in the literature. The experimental results show that the proposed method speeds up the analysis considerably, compared to a global analysis, while still providing high accuracy. The sliding window method is a novel approach to address the performance problems of numerical algorithms for the solution of the chemical master equation. The method efficiently approximates the probability distributions at the time points of interest for a variety of chemically reacting systems, including systems for which no upper bound on the population sizes of the chemical species is known a priori.

Propagation Models for Biochemical Reaction Networks

Maria-Emanuela-Canini Mateescu

In this thesis we investigate different ways of approximating the solution of the chemical master equation (CME). The CME is a system of differential equations that models the stochastic transient behaviour of biochemical reaction networks. It does so by describing the time evolution of probability distribution over the states of a Markov chain that represents a biological network, and thus its stochasticity is only implicit. The transient solution of a CME is the vector of probabilities over the states of the corresponding Markov chain at a certain time point t, and it has traditionally been obtained by applying methods that are general to continuous-time Markov chains: uniformization, Krylov subspace methods, and general ordinary differential equation (ODE) solvers such as the fourth order Runge-Kutta method. Even though biochemical reaction networks are the main application of our work, some of our results are presented in the more general framework of propagation models (PM), a computational formalism that we introduce in the first part of this thesis. Each propagation model N has two associated propagation processes, one in discrete-time and a second one in continuous-time. These propagation processes propagate a generic mass through a discrete state space. For example, in order to model a CME, N propagates probability mass. In the discrete-time case the propagation is done step-wise, while in the continuous-time case it is done in a continuous flow defined by a differential equation. Again, in the case of the chemical master equation, this differential equation is the equivalent of the chemical master equation itself where probability mass is propagated through a discrete state space. Discrete-time propagation processes can encode methods such as the uniformization method and the fourth order Runge-Kutta integration method that we have mentioned above, and thus by optimizing propagation algorithms we optimize both of these methods simultaneously. In the second part of our thesis, we define stochastic hybrid models that approximate the stochastic behaviour of biochemical reaction networks by treating some variables of the system deterministically. This deterministic approximation is done for species with large populations, for which stochasticity does not play an important role. We propose three such hybrid models, which we introduce from the coarsest to the most refined one: (i) the first one replaces some variables of the system with their overall expectations, (ii) the second one replaces some variables of the system with their expectations conditioned on the values of the stochastic variables, (iii) and finally, the third one, splits each variable into a stochastic part (for low valuations) and a deterministic part (for high valuations), while tracking the conditional expectation of the deterministic part. For each of these algorithms we give the corresponding propagation models that propagate not only probabilities but also the respective continuous approximations for the deterministic variables.

EPFL2011

Hybrid Numerical Solution of the Chemical Master Equation

Maria-Emanuela-Canini Mateescu, Verena Wolf

We present a numerical approximation technique for the analysis of continuous-time Markov chains that describe networks of biochemical reactions and play an important role in the stochastic modeling of biological systems. Our approach is based on the construction of a stochastic hybrid model in which certain discrete random variables of the original Markov chain are approximated by continuous deterministic variables. We compute the solution of the stochastic hybrid model using a numerical algorithm that discretizes time and in each step performs a mutual update of the transient probability distribution of the discrete stochastic variables and the values of the continuous deterministic variables. We implemented the algorithm and we demonstrate its usefulness and efficiency on several case studies from systems biology.

2010

Stabilized Numerical Methods for Stochastic Differential Equations driven by Diffusion and Jump-Diffusion Processes

Adrian Blumenthal

Stochastic models that account for sudden, unforeseeable events play a crucial role in many different fields such as finance, economics, biology, chemistry, physics and so on. That kind of stochastic problems can be modeled by stochastic differential equations driven by jump-diffusion processes. In addition, there are situations, where a stochastic model is based on stochastic differential equations with multiple scales. Such stochastic problems are called stiff and lead for classical explicit integrators such as the Euler-Maruyama method to time stepsize restrictions due to stability issues. This opens the door for stabilized explicit numerical methods to efficiently tackle such situations. In this thesis we introduce first a stabilized multilevel Monte Carlo method for stiff stochastic differential equations. Using S-ROCK methods we show that this approach is very efficient for stochastic problems with multiple scales, but also for nonstiff problems with a significant noise part. Further, we present an improved version of the stabilized multilevel Monte Carlo method by considering S-ROCK methods with a higher weak order of convergence. Then we extend the S-ROCK methods to jump-diffusion processes. We study in detail the strong order of convergence of the newly introduced methods and we discuss the corresponding mean square stability domains. In the next part we present the multilevel Monte Carlo method for jump-diffusion processes. We state and prove a theorem that indicates the computational cost required to achieve a certain mean square accuracy. In the numerical section we compare the multilevel Monte Carlo approach to two variance reduction techniques, the antithetic and the control variates. We also show how the S-ROCK method for jump-diffusion processes, introduced in this thesis, can be used to create a stabilized multilevel Monte Carlo method for jump-diffusions that handles stiffness and considers the inclusion of jumps at the same time. Finally, we propose in this thesis a variable time stepping algorithm that uses S-ROCK methods to approximate weak solutions of stiff stochastic differential equations. A rigorous analytical study is carried out to derive a computable leading term of the time discretization error and an adaptive algorithm is suggested that adapts the time grid and adjusts the number of stages of the S-ROCK method simultaneously.