I'm looking for a decent implementation of Gillespie's Direct Method in Python, as if I code the algorithm myself I'm nigh positive I'll do it inefficiently.
Anyone have a favorite?
(Disclaimer: I haven't used any of the packages or code below.)
Finding any implementation of Gillespie's method in Python was a bit of a challenge; the most fruitful search terms seemed to be "tau leap" or "kinetic Monte Carlo." This blog post implements Gillespie's algorithm, though it's not clear that it's efficient. One of the commenters mentions some other Gillespie/SSA algorithm implementations in Python that solve problems out of a textbook; the relevant problems are 6.3 through 6.6 on this website.
The most promising implementations in Python that I found were StochPy and Python language bindings to COPASI. I also found what looked like promising C++ codes in SPPARKS, STOCKS, and StochKit; they're not Python, but maybe you can wrap around one of them in Python if it looks particularly promising to you.
I would warmly recommend to use StochKit, which is really extremely reliable and efficient. It implements the Gillespie's direct method as well as its most popular optimized variants. Furthermore, it has a built-in mechanism for selecting the appropriate simulation scheme depending on your network.
However, it is written in C++. If you really need to integrate it with a Python, that shouldn't be a problem. I use it extensively in MATLAB scripts (basically, MATLAB writes a descriptor of the simulation, calls StochKit, and finally performs the data analysis). I'm sure you could do pretty much the same thing with NumPy or SciPy.
As part of the lab the developed and maintains StochKit. I am happy to hear that it is highly recommended in the previous answers. However, I wanted to update everyone. There already is a python wrapper for StochKit: GillesPy.
I would also recommend you checkout StochSS, This is an fully functional modeling and simulation IDE that uses StochKit and as one of it's simulations engines. It also includes spatial simulation with the python package PyURDME. Model building, and 3D visualization tools (plus cloud computing integration).
As you may already know if you're asking this question, Gillespie's method is often used in computational systems biology. There are a number of software system that provide implementations; you may want to browse the SBML Software Guide to find some of them (specifically, the ones that are known to support SBML, something that I realize may not be relevant to your needs). However, Geoff Oxberry is correct that Python-based implementations of Gillespie's algorithms seem to be uncommon compared to C++ and Java implementations. In addition to StochPy, I'm aware of the following Python-based packages (in alphabetical order):
If, in your application, you can use software that's coded in something other than Python but can be wrapped with a Python interface, then there is
Finally, if you don't need your simulation to be run locally and can instead use web services, then there are the following options:
These lists are still incomplete, but hopefully they will help.
I wrote a small Python package called OOGillespie whose main feature is that it is rather flexible and embraces object-oriented programming.
Efficiency was not the main priority when writing OOGillespie, but I tried to keep things fast, where it did not lead to a major slowdown. OOGillespie will probably not be able to compete with other implementations for standard tasks, but it may be the only implementation that can handle some non-standard ones. If your bottleneck is computing transition rates from large matrices, OOGillespie may be competitive again by allowing you to use NumPy matrix operations for this purpose. Finally, depending on how experienced you are with classes, OOGillespie may cut down your programming time.
At the end of the day, I do not think that you can simply rank Gillespie implementations without knowing the specific application as a major point is what kind of events you have and how you can calculate tehm.
I have made my own simple implementation and created short tutorial, hope it will be helpfull:
i just converted the gillespie's algorithm into python code:
1.Define total simulation time T and set current time t=0; 2.Initialize the state of the system x – choose randomly or set zero state; 3.Find the sum Q of all the possible transition rates from the current state x; 4.Simulate the time Δt until the next transition by drawing from an exponential distribution with parameter Q; 5.Generate a random number r from a uniform distribution; 6.Set the next state n as follows: If 0 < r < q(x,0)/Q, choose transition 0; if q(x,0)/Q < r < (q(x,0) + q(x,1))/Q choose transition 1, and so on; 7.Update the current time: t=Δt+t; 8.Repeat 3-7 until t≤T.
The source code in python:
additionally, the result of simulation is compared with the exact result using mean squared error (MSE).