The stochastic simulation algorithm (SSA) is a really common method to model stochastic systems with a continous index (usually time) and a discrete state space (originally, number of molecules interacting in a gas phase). It is a workhorse in the field of systems biology, where biomolecular reactions with discrete molecule numbers are important and common. There are tons of great software packages to run the SSA in whatever language you want, but the algorithm itself is quite simple, and worth taking a minute to code up. I recently started playing a bit with Julia, and thought I’d code up a simple SSA.

Algorithm

The codes below follow the direct method to the SSA. The idea (without going into details about well-mixedness, system size, etc.) is to randomly pick the reaction that will happen next in the system, based on the probability, or propensity that the reaction occurs in the infinitesimal time step \(dt\).

For a simple example, let’s start with a model of birth and death. In this model, molecules, cells, prairie dogs, or whatever you want can do two things: be born and die. They enter the world with a propensity \(k_r dt\) and they exit the world with a propensity \(\gamma x dt\), where \(x\) is the number of things in the world. Here, propensity means “probability that the event happens”. When propensities are affine linear in \(x\) (you can write it as \(y=mx+b\)), it’s convenient to write them in terms of reactions that multiply \(x\), which we’ll call \(\mathbf{W0}\) and terms that don’t, which we’ll call \(\mathbf{W1}\). Therefore, we can write the propensity for such systems as

\[\mathbf{w} = \mathbf{W0}\mathbf{x}+\mathbf{W1}.\]

Or in the simple birth-death model, as

\[\mathbf{w} = \begin{pmatrix} \gamma \\ 0 \end{pmatrix}x + \begin{pmatrix} 0 \\ k_r \end{pmatrix}.\]

The next thing we need to run the SSA is the stoichiometry matrix, \(\mathbf{S}\), which tells us how species change for each reaction. For the codes below, I have each row of \(\mathbf{S}\) be a different species, and each column be a different reaction. For our simple birth and death model, this amounts to an \(\mathbf{S}\) matrix which looks like

\[\mathbf{S} = \begin{pmatrix} 1 & -1 \end{pmatrix}\]

because we have just 1 species, \(x\), and 2 reactions, one of which increases \(x\) by 1, and one of which decreases \(x\) by 1.

As Gillespie put in his original paper, we now need to answer two questions:

  • When does the next reaction happen?
  • What reaction happens?

Without going into exhaustive detail, the time of the next reaction is an exponential random variable with a mean given by the sum of the propensities \(\mathbf{w}\).
We can simply draw a random number from the exponential distribution and get the time of the next reaction. Keep in mind that the propensities \(\mathbf{w}\) will change at each time step, so the exponential distribution that determines the time until the next reaction is also different at each time step.

To answer the second question, we just remember that the propensities are like probabilities - if a reaction has a high propensity, it is more likely to happen, and vice versa for low propensity reactions. We want to simply sample from this distribution, which we can do using the CDF trick.

After selecting the reaction that happens next and when that reaction happens, we simply advance the time index, and update the state index accordingly. Do this over and over and you’ve generated a sample path or trajectory for this process.

Below are simple codes for the SSA in python, MATLAB, and Julia, with a link to the source. In terms of performance (for a simple birth and death model), Julia and MATLAB perform about equally well, and python is quite a bit slower.

Code

MATLAB

function end_state = ssa(t_f, x0, stoich_mat, W1,W0)

  x = x0;
  t = 0;
  while (t < t_f)

    a = W1*x+W0;
    a0 = sum(a);

    u1 = rand();
    u2 = rand();

    % Update time
    tau = min(log(1/u1)/a0, t_f - t);
    t = t + tau;

    % Update state
    r = find(cumsum(a) > u2*a0, 1, 'first');
    x = x + stoich_mat(r, :);
  end
  end_state = x;
end

source

Python

def ssa(W0,W1,S,tf,x0):
    '''
    a simple SSA. 
    '''
    x = x0 
    t = 0
    while t<tf:
        a = W1.dot(x)+W0
        a0 = np.sum(a.ravel()) 
        # get some random numbers
        u1 = np.random.rand()
        u2 = np.random.rand()

        # update the time vector 
        tau = np.min([np.log(1/u1)/float(a0), tf-t])
        t+=tau

        # update the state
        r = np.where(np.cumsum(a)>u2*a0)[0][0]
        x += S[r,:].reshape((1,1))
    return x

source

Julia

function ssa(W0,W1,S,tf,x0)
    # Simple ssa function in julia. 
    x = x0
    t = 0 
    while t<tf
        a = W1*x+W0
        a0 = sum(a)
        
        # Get some random numbers. 
        u1 = rand()
        u2 = rand()

        # update the time
        tau = minimum([log(1/u1)/a0, tf - t]);
         t = t+tau

         # update state
         r = find(cumsum(a).>u2*a0)[1]
         x = x + S[r,:]
     end
     return x
end

source

Thanks to Huy Vo for the initial MATLAB code and kind suggestions.