Discrete Stochastic Simulations

We consider discrete stochastic simulations that are modelled with a set of species and a set of reactions that transform the species' amounts. Instead of using a continuum approximation and dealing with species mass or concentration, the amount of each species is a non-negative integer which is the population. Depending on the species, this could be the number of molecules or the number of organisms, etc. Reactions transform a set reactants into a set of products, each being a linear combination of species with integer coefficients.

Consider a system of N species represented by the state vector X(t) = (X1(t), ... XN(t)). Xn(t) is the population of the nth species at time t. There are M reaction channels which change the state of the system. Each reaction is characterized by a propensity function am and a state change vector Vm = (Vm1, ..., VmN). am dt is the probability that the mth reaction will occur in the infinitesimal time interval [t .. t + dt). The state change vector is the difference between the state after the reaction and before the reaction.

To generate a trajectory (a possible realization of the evolution of the system) one starts with an initial state and then repeatedly fires reactions. To fire a reaction, one must answer the two questions:

  1. When will the next reaction fire?
  2. Which reaction will fire next?
Let the next reaction have index μ and fire at time t + τ. Let α be the sum of the propensities. The time to the next reaction is an exponentially distributed random variable with mean 1 / α ; the probability density function is P(τ = x) = α e- α x. The index of the next reaction to fire is a discrete random variable with probability mass function P(μ = m) = am / α.