A Modified Stochastic Simulation Algorithm for Time-Dependent Intensity Rates

There are two main approaches in the mathematical modeling of coupled systems of (bio)chemical reactions: continuous, represented by differential equations whose variables are concentrations or discrete, represented by stochastic processes whose variables are numbers of molecules. The latter approach is used mostly for biochemical systems with a low to moderate number of molecules of certain species and this kind of systems are typically modeled as continuous time - discrete state Markov Process. There are exact stochastic algorithms to simulate state trajectories of discrete, stochastic systems and these algorithms are based on methods that are rigorously equivalent to the Master Equation approach. Two of the most widely used methods for simulating the stochastic dynamics of a chemical system are the exact stochastic simulation algorithm (SSA, known also as Gillespie algorithm) and its approximate variant, the tau-leaping algorithm. This paper describes a modified version of SSA - First reaction method - by letting the intensity rates of the reactions to be functions of time. The importance of this adaptation is obvious when considering some classes of biological models (for example, the one involving circadian rhythm). The underlying assumptions are that the system is well stirred such that at any moment, each reactions occur with equal probability at any position, that each reaction, once occurred, completes instantaneously (there are no reactions with delay involved) and that the system is non stiff (there are no different time scales of the reactions involved).