The OOGillespie Module¶

This module provides a class that allows to simulate a stochastic process with discrete, Poissonian events with the Gillespie algorithm.

In contrast to other implementations it embraces object-oriented programming, which has its advantages and disadvantages:

It is easy to use if you are already familiar with object-oriented programming in Python.
It may be good for you if its way of doing things comes to you naturally.
It is very flexible.
It is not fast – however, your bottleneck may be elsewhere anyway, e.g. computing the rates.

A Simple Example: Birth and Death¶

We want to implement a simple birth–death process where births happen with a constant rate of \(α=0.1\) and each individual dies with a rate of \(ω=0.001\), starting at a population of \(N=0\).

We start with importing the Gillespie class as well as setting the control parameters:

from oogillespie import Gillespie, Event

α = 0.1
ω = 0.001

We now start defining the class for our specific process. It inherits from Gillespie and most of the remainder of this example will happen in the class definition.

class birth_death(Gillespie):

The first of the method we need is initialise, which is automatically called upon initialisation and sets the initial state of the system. It gets passed parameters from the constructor calls. Here we only need to set \(N\), which we want to default to 0:

	def initialise(self,N=0):
		self.N = N

Next we implement the birth event. For this, we need a method that is called upon such an event and increases \(N\) by 1. To register this method as an event and specify its rate, we use the decorator Event and provide the constant rate \(α\) as an argument:

	@Event(α)
	def birth(self):
		self.N += 1

Now, we implement the death event. The main difference to birth events is that the rate is not constant. Therefore we first define a method death_rate which calculates the rate of the event. We then pass this method as an argument to the Event decorator.

	def death_rate(self):
		return ω*self.N
	
	@Event(death_rate)
	def death(self):
		self.N -= 1

Finally, we implement a function state that returns whatever we consider relevant information about the state of the system. Here it is the time and \(N\):

	def state(self):
		return self.time,self.N

This concludes the definition of the class. We can now simulate a realisation of such a process by iterating over an instance of this class. To make the loop terminate, we have to specify either the maximum number of steps or the maximum time as an argument. We here want the maximum number of steps to be 10000. Finally, to keep this example simple, we print the result, but we could as well further process or plot it:

for time,N in birth_death(max_steps=10000):
	print(time,N)

Taken together, our script looks like this:

from oogillespie import Gillespie, Event

α = 0.1
ω = 0.001

class birth_death(Gillespie):
	def initialise(self,N=0):
		self.N = N
	
	@Event(α)
	def birth(self):
		self.N += 1
	
	def death_rate(self):
		return ω*self.N
	
	@Event(death_rate)
	def death(self):
		self.N -= 1
	
	def state(self):
		return self.time,self.N

for time,N in birth_death(max_steps=10000):
	print(time,N)

An Advanced Example: Rock, Paper, Scissors¶

We want to implement a process with three populations \(N_0, N_1, N_2\) whose individuals transform into each other. Individuals of population \(N_i\) are born with a rate \(α_i\). Transformations from \(i\) to \(j\) happen with a rate \(N_i·A_{ij}·N_j\).

More specifically, we want \(α = (0.2, 0.01, 0.01)\) and the transformations to be cyclical (rock–paper–scissors):

\[\begin{split}A = \pmatrix{ 0 & 0.5 & 0 \\ 0 & 0 & 0.3 \\ 0.1 & 0 & 0 \\ }.\end{split}\]

We start with importing the Gillespie class as well as setting the control parameters as NumPy arrays:

from oogillespie import Gillespie, Event
import numpy

α = numpy.array([ 0.2, 0.01, 0.01 ])

A = numpy.array([
		[  0  , 0.5 ,  0  ],
		[  0  ,  0  , 0.3 ],
		[ 0.1 ,  0  ,  0  ],
	])

We then continue with defining the class for our process. Nothing new happens in the initialise and state functions, except that N is an array – and we return a copy of it to avoid changes to the internal state carrying over:

class RPS(Gillespie):
	def initialise(self,N=None):
		self.N = numpy.asarray(N or [0,0,0])
	
	def state(self):
		return self.time,self.N.copy()

In contrast to the previous example, there are now three different birth events, one for each population. However, being healthily lazy, we do not want to write three nearly identical functions for this but want to have a parametrised event, where the parameter specifies which population increases.

To achieve this, we pass a sequence (α) instead of a number as an argument to the Event decorator. This way, we specify that the birth event comes in different versions with different rates, one for each element of α. Whenever such an event happens, the respective index of α is passed to birth as an additional argument.

	@Event(α)
	def birth(self,i):
		self.N[i] += 1

For the transformation event, two things change:

The rate depends on the current state. This is handled like before, i.e., by passing a method instead of a number to the Event decorator.
The event depends on two parameters. This is reflected by the parameters i and j of transform and transform_rate returning a two-dimensional array of rates instead of a number.

	def transform_rate(self):
		return self.N*A*self.N[:,None]
	
	@Event(transform_rate)
	def transform(self,i,j):
		self.N[i] -= 1
		self.N[j] += 1

Altogether our code would look like this:

from oogillespie import Gillespie, Event
import numpy

α = numpy.array([ 0.2, 0.01, 0.01 ])

A = numpy.array([
		[  0  , 0.5 ,  0  ],
		[  0  ,  0  , 0.3 ],
		[ 0.1 ,  0  ,  0  ],
	])

class RPS(Gillespie):
	def initialise(self,N=None):
		self.N = numpy.asarray(N or [0,0,0])
	
	def state(self):
		return self.time,self.N.copy()
	
	@Event(α)
	def birth(self,i):
		self.N[i] += 1
	
	def transform_rate(self):
		return self.N*A*self.N[:,None]
	
	@Event(transform_rate)
	def transform(self,i,j):
		self.N[i] -= 1
		self.N[j] += 1

for time,N in RPS(max_steps=100000):
	print(time,*N)

Command reference¶

class Event(rates)¶

Decorator that marks a method (of a subclass of Gillespie) as an event.

There are several valid use cases of the argument rates of the decorator and the number of arguments of the decorated method:

The event method has no arguments other than self, and rates is a non-negative number specifying the rate of the event.
The event has one argument other than self, and rates is a sequence of non-negative numbers. In this case, the numbers specify the rates of different variants of the event. If an event happens, the location of the respective rate in the sequence is passed as an argument to the event method.
A generalisation of the above: The event has k arguments other than self, and rates is a nested sequence of non-negative numbers as an argument, with k levels of nesting. All sequences on the same level must have the same length. For example, for k=2, rates is a matrix. If an event happens, the location of the respective rate in the sequence is passed as arguments to the event method.
rates is a method returning a number or (nested) sequence of numbers as above. This method must take self and only self as an argument.

In any case the number, sequence, or method providing the rates must be defined before the event method.

set_parent(parent)¶: Sets the instance of the class (usually a subclass of Gillespie) on which the event is to be executed. Without this, the event is useless. Gillespie calls this when registering Events.

actions()¶: Yields zero-argument functions (actions) for each variant of the event.

class Gillespie(max_steps=1000, max_t=numpy.inf, time=0, seed=None, **kwargs)¶

Base class for Gillespie models. This class only works if inherited from and if the methods initialise and state are replaced. Also, at least one method has to be marked as an event with the respective decorator.

When an instance of the class is iterated over, the actual simulation is performed and the iterator returns the output of state after each step.

Parameters:

time: The starting time.
max_steps
max_t: The maximum number of steps or time, respectively. Before either of the two is exceeded (or all event rates become zero), the simulation is aborted.
seed: Seed for random number generation. If None, system time or similar is used (like random.seed).
kwargs: Keyword arguments to be forwarded to initialise.

abstract initialise()¶: You must overwrite this method. It is called when initialising the integrator. Use it to set internal parameters to initial values. It gets passed surplus arguments from the constructor.

abstract state()¶: You must overwrite this method. It is called to determine the return value when iterating/simulating. Use it to return whatever properties you are interested in.

The OOGillespie Module¶

A Simple Example: Birth and Death¶

An Advanced Example: Rock, Paper, Scissors¶

Command reference¶

Table of Contents

This Page