Adaptive Discrete Galerkin Methods Applied to the Chemical Master Equation

In systems biology, the stochastic description of biochemical reaction kinetics is increasingly being employed to model gene regulatory networks and signaling pathways. Mathematically speaking, such models require the numerical solution of the underlying evolution equation, known as the chemical master equation (CME). Until now, the CME has primarily been treated by Monte Carlo techniques, the most prominent of which is the stochastic simulation algorithm [D. T. Gillespie, J. Comput. Phys., 22 (1976), pp. 403-434]. The paper presents an alternative, which focuses on the discrete partial differential equation (PDE) structure of the CME. This allows us to adopt ideas from adaptive discrete Galerkin methods as first suggested by Deuflhard and Wulkow [IMPACT Comput. Sci. Engrg., 1 (1989), pp. 269-301] for polyreaction kinetics and independently developed by Engblom. From the two different options for discretizing the CME as a discrete PDE, Engblom chose the method of lines approach (first space, then time), whereas we strongly advocate use of the Rothe method (first time, then space) for clear theoretical and algorithmic reasons. Numerical findings at two rather challenging problems illustrate the promising features of the proposed method and, at the same time, indicate lines of necessary further improvement of the method worked out here.