Modelling with Python

21/05/2015

Jens Hahn

2. Discrete stochastic modelling

Gillespie Algorithm (SSA)

This time we want to have a look on a stochastic model. These models are not deterministc, so every simulation will give you a slightly different result. Furthermore, these models are discrete, so you're looking on the behaviour of single molecules. Since these models describe single molecules, they are quite slow, but surely more accurate for the simulation of systems with only few particles, e.g. gene regulation.

Let's start

We will use the same reaction as in the ODE part before, a molecule A and a molecule B react to a molecule C in a reversible reaction: $$\textrm{A} + \textrm{B} \rightleftharpoons \textrm{C}$$

Reaction propensities

In this formalism we have to define reaction propensities for every reaction. They consist basically of two parts:

Firstly, the actual probability for the reaction to occur when 2 molecules meet Let's call them $\textrm{k}_1$ and $\textrm{k}_2$.

Secondly, we need to think about the fact that there are several molecules per species. Quite similar to the mass action rate law, only here its about possibilities of pairs (permutations), h.

The actual reaction propensity is the product of the probability and the permutation term.

Here:
R1: second order
$\textrm{A} + \textrm{B} \Rightarrow \textrm{h}_1 = \textrm{N}_A \times \textrm{N}_B$

R2: first order
$\textrm{C} \phantom{ + \textrm{B}} \Rightarrow \textrm{h}_2 = \textrm{N}_C$

Propensity a = h k

$$\textrm{k}_1 = 0.2$$$$\textrm{k}_2 = 0.1$$$$\textrm{a}_1 = \textrm{k}_1 * \textrm{h}_1$$$$\textrm{a}_2 = \textrm{k}_2 * \textrm{h}_2$$

These propensities will give as the information which reaction actually occurs.

Time step

We need also the information at what time a reaction occurs. Therefore, we will use this formula: $$\tau = - \frac{1}{a}\ln(\textrm{ran_num})$$ In this formula, a is the sum of the reaction propensities and ran_num is a random number between 0 and 1.
The formula basically describes, that with increasing time without any reaction, the probability for any reaction to occur also increases (exponentially).

Simulation algoritm

Gillespie SSA:
This algorithm consists of two random numbers between 0 and 1. One will tell you which reaction occurs (use normalised propensities!) and the other one will tell you when the reaction occurs.

Python

Let's start the simulation. We begin with defining the initial states and the paramters:

Now it's your turn! What do we need to simulate our model for... let's say 100 seconds?

Here some help:

1. The simulation time
2. Lists to save the simulation results
3. A while loop with two random numbers
4. Some matplotlib to visuali