Gillespie Algorithm

A Sample of Gillespie's Algorithm (Direct Method) for Autocatalytic Reaction Cycle

Source Code (Ruby)
Source Code (C Language)

Direct Method is one of the exact stochastic simulation algorithms (SSA), which is invented by Gillespie in 1977^[1].
It is not efficient but so simple that we can learn a stochastic method for the first time.
A simple autocatalytic reaction system is implemented in this example by using the Direct Method with Ruby.
Basically, the calculation cost of the DM depends on the total amount of chemicals, namely, O(S) S is the sum of molecules.
In order to improve this sample into one of C language, the style of object oriented coding is avoided as much as possible in the Ruby source code.

Chemical Reactions

  X0 + X1 -- c0 --> X1 + X1
  X1 + X2 -- c1 --> X2 + X2
  X2 + X0 -- c2 --> X0 + X0

The time evolution of chemical species

The correct time evolution by Runge-Kutta method (Ordinary Differencial Equation (ODE) solution) without noise
Source Code (C with GSL)

Sorting Direct Method (SDM)

Source Code (Ruby)

There are some parts which are inefficient for the calculation in the source code above. There are mainly TWO parts that should be improved.

Selecting a reaction
Updating propencities

With the increase of number of reactions, the searching depth will increase. One of the simple ideas for this is to sort the list of reactions in order to search the reaction quickly that has a higher propencity than the others. McCollum et.al has developed Sorting Direct Method (SDM)^[3] in which a reaction is bubbled up in the reaction list when it fires.

Chemical Reactions

 X_(i) + X_(i+1) -- ci --> X_(i+1) + X_(i+1)  (i=0,...,N-2)
 X_(N-1) + X_(0) -- ci --> X_(0) + X_(0)  (i=N-1)
 (N: Number of Chemical Species >= 2)

Seraching Average Times per reaction: DM v.s. SDM

Averagely, the search depth of SDM is lower than the original DM. However, the sorting process becomes overhead. In this example, SDM's speed is mostly the same as DM.

Calculation Time: DM v.s. SDM

In this example, the cost of sorting overhead is slightly bigger than the searching advantage. It is expected that SDM would be effective in the case that a specific reaction often happens with high propencity in a period. However, the order of calculation cost of SDM is still O(N) as well as DM, which means the calculation time increases in proportion to the number of reactions.

Logarithmic Direct Method (LDM)

Source Code (Ruby)

The linear search that is the simplest searching algorithm is used in both DM and SMD. The alternative idea for that is binary search^[6]. The calculation cost for the searching is expected to be O(log(N)).

Calculation Time: DM v.s. LDM

It is no need to sort the reaction list in this algorithm.

Dependency Graph

Source Code (Ruby)

All propencities do not have to be updated when a reaction fires. One of the ideas for this is reducing the updating times of propencity. Gibbson et.al. used a dependency graph of reactions in order to reduce the updatting times^[4]. In this autocatalytic reaction cycle, only three propencities should be updated in every step. For example, if

 X[1] + X[2] -> 2 X[2]

fired, then X[1] and X[2] would change. X[1] and X[2] are related to

 X[0] + X[1] -> 2 X[1] 
 X[1] + X[2] -> 2 X[2] 
 X[2] + X[3] -> 2 X[3]

Then,

 propencity[0], propencity[1], and propencity[2]

should be updated. The effect of this in this example will increase in proportion to the number of reactions.

Calculation Time: DM v.s. Dependency Graph

However, the order of calculation cost in this DM using dependency graph is still O(N) as well as SDM.

Partial-propencity Direct Method (PDM)

Source Code (Ruby)

The original Direct Method and its variants generally calculate propencities for reactions. That is why the order of calculation cost is O(N). Ramaswamy et.al. developed Partial-propencity Direct Method (PDM)^[5] in which propencities for chemical species are calculated in every step. Then, it realized that the cost becomes not O(N) but O(M), where M is the number of chemical species. However, it is limited only for a system of elementary reactions.

chemical reactions

  X0 + X1 -- c0 --> X1 + X1
  X1 + X2 -- c1 --> X2 + X2
  X2 + X0 -- c2 --> X0 + X0

The time evolution of chemical species

Chemical Reactions

 X_(i) + X_(i+1) -- ci --> X_(i+1) + X_(i+1)  (i=0,...,N-2)
 X_(N-1) + X_(0) -- ci --> X_(0) + X_(0)  (i=N-1)
 (N: Number of Chemical Species >= 2)

Calculation Time: PDM

In this example, the number of reactions, N, and the number of chemical species, M, are the same, so that the calculation cost is O(N) and the calculation cost is not so different from the original DM.

Source Code (Ruby)

X_(i+1) + X_(i+1) (i=0,...,N-2) X_(N-1) + X_(0) -- ci --> X_(0) + X_(0) (i=N-1) (N: Number of Chemical Species >= 2)

Calculation Time: DM v.s. PDM v.s. SPDM

-->

Reference: