Lesson 11: Python Code for Gillespie Algorithm

NOTE - I made a slight mistake in the if / elif statement checks in this video. The "less than" conditions should actually be "less than or equal to". Failing to check for the equality can lead to errors in rare cases where the random product exactly equals the sum of the rates. So I have updated the code to be correct on Github. Sorry about that!

This video is about the Python code we can use to run the Gillespie Algorithm for a simulation the transcription of mRNA from a gene.

Here is a link to the original Gillespie paper:

Here is a video of Jeff Gore, an MIT professor, explaining the Gillespie Algorithm in more depth:

All the code from my videos is available on my Github:

If you're are trying to catch up on the basics of Python, differential equations, or cell biology, here are some good resources:

Python Basics-

Differential Equations

Cell Biology

Central Dogma

Thanks for watching and let me know if you have any questions!

[MikeSaintAntoine]

Hi everyone, thanks for watching! Here's the link to the code on my Github in case you want to try it for yourself (also feel free to add me as a friend on Github if you want 🙂) :

[qingqiwen]

it's much more easy and clear for me to understand! THANKS

[nikhilbhatia]

What if I want to use Gillespie algorithm to find density profiles at steady state for a lattice of size L where particle enter at a specific rate move throughout the lattice at different rates and exit the last lattice site with some different rate? I am doing simulations using Gillespie for master equation obtained for each site. Am I doing this right?

[PinakiSwain]

Thanks for uploading this. Very clear. Just a comment. It looks to me that the first condition in if block and the second condition in the elif block are always going to be true. Can't we replace them with
if rand * rate_sum <= 0:
Production statement
else:
Degradation statement

[mr.z]

Thanks for your awesome video! I have a question related to rate_sum. rate_sum = k+gamma*current_X, as current_X has a chance to be negative then rate_sum may be negative as well although it's very unlikely but still possible. But rate_sum will be used as exponential distribution lambda which should be always positive. May I ask how do you constrain the rate_sum to make sure it's always positive?

[deyomash]

Heyoh! I constantly get the error that my t and (variable) is of different dimension when i want to plot. I did an SIR model but essentially the same, i just added I and R to +1 or -1, shouldnt be much difference right?