Pade SSA: A frequency domain method for estimating the dynamics of stochastic reaction networks

Dynamic analysis and control of living cells relies on mathematical representations of cellular processes that are themselves modelled as biomolecular reaction networks. Stochastic models for biomolecular reaction networks are commonly employed for analysing intracellular networks having constituent species with low-copy numbers. In such models, the main object of interest is the probability distribution of the state vector of molecular counts which evolves according to a set of ordinary differential equations (ODEs) called the Chemical Master Equation (CME). Typically this set is very large or even infinite, making the CME practically unsolvable in most cases. Hence the outputs based on the CME solution, like the statistical moments of various state components, are generally estimated with Monte Carlo (MC) procedures by simulating the underlying Markov chain with Gillespie's Stochastic Simulation Algorithm (SSA). However to obtain statistical reliability of the MC estimators, often a large number of simulated trajectories are required, which imposes an exorbitant computational burden. The aim of this paper is to present a frequency domain method for mitigating this burden by exploiting a small number of simulated trajectories to robustly estimate certain intrinsic eigenvalues of the stochastic dynamics. This method enables reliable estimation of time-varying outputs of interest from a small number of sampled trajectories and this estimation can be carried out for several initial states without requiring additional simulations. We demonstrate our method with a couple of examples.