Autocatalytic Reaction System by Gillespie's Algorithm (Direct Method)

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 here by using the Direct Method on Ruby and C language. Basically, the calculation cost of the DM depends on the total amount of chemicals, namely, O(S), S is the sum of molecules. The style of object oriented coding in the Ruby code is avoided as much as possible in order to improve this sample into one of C language from Ruby code.

Source Code

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

An example of time evolution of chemical species


The correct time evolution by Ordinary Differencial Equation (ODE) solution (Eular method)


Source Code
Eular method is the simplest numerical solution for ordinary differencial equations, and it is so simple method to lean a numerical calculation for differential equations for the first time as well as the Direct Method. But actually the calculation error is larger than the other methods such as Runge-Kutta method, so that it is not a practical method.

Properties of Autocatalytic Reaction System

This Autocatalytic Reaction system is so-called an Harmonic Oscillator. A typical example of oscilator is Spring, which is sometimes expressed by sin (cos) function. Another famous oscillator model is Lotka-Volterra model, which is famous as a predetor-preyer model. These are expressed mathematically by Replicator Equation, in other words, these are mathematically equally to the same differencial equations of this Autocatalytic Reaction system. Although it oscilates regularly until 4 elements system, it becomes chaotic oscilation over 5 elements. It is, for example, a system of 10 chemicals below.


Another important oscilator property is that the amplitude depends on the initial value. In other words, the system is unstable for external noise. For example, the following graph is the result where the x[0] value is added 100.0 only once at the time 0.5.

You can understand this behavior more clearly in the phase diagram where the x axis is x[0] and y axis is x[1] below. The volume of the phase space changes at the time 0.5.

The Autocatalytic Reaction system calculated by DM method becomes actually an oscillator with noise. So then it will stop suddenly because of the randomness (noise). In more detail to say, two of three values becomes 0 and all the reaction probabilities become 0, then the system stops.

Dynamic Link Library (libdm.so)

Sometimes a useful function should be separated into a library file. This is a sample below. The common calculation part of Direct Method is coded in dm.h and dm.c. You can download all the files and just type 
  make
on your terminal (Windows may need, for example, Cygwin environment to compile them).

Source Code
  1. dm.h
  2. dm.c
  3. main.c
  4. Makefile

Also, you can use this dm library for your own simulation. Change parameters, init function, and update_p function in main.c, then make. That's it. If you want to use the library in the other system, include dm.h header file and put dm.h and libdm.so on the current directory as well as your main .c file. Then
  gcc (your main .c) -I. -L. -ldm
. You do not have to recompile libdm.so file, because it is a dynamic link library. To recompile libdm.so again, 
  gcc -shared dm.c -o libdm.so


Ruby Extented Library

Ruby gives us a fun because of its very flexible and easy grammer. But one of disadvantages is its speed. In particular, this sort of numerical calculation, a single process and including many iteration processes, is very slow in comparison to a compiler language code such as C language. In order to speed it up, we can use C language library from Ruby code, which is called Ruby Extended Library. Here is an example of the Ruby Extended Library to use dm dynamic link library above.

Source Code
  1. dm.c
  2. extconf.rb
  3. main.rb
You should make a binary file from the source code. Put the files above and dm.h, libdm.so in the same directory (Do NOT put the other files). Then type 
  ruby extconf.rb
  make
extconf.rb creates Makefile from the files in the directory. Then, you can get a Ruby Extended Library dm.so (or dm.bundle). You can call the library from Ruby code as well as Ruby library. main.rb is a sample to call the library. It calculate the same Autocatalytic Reaction system. 

  ruby main.rb
will output the same result as well. The speed is about 4 times faster than pure Ruby code, but about 15 times slower than the original C code. One of useful points is that we do not have to recompile source code again. Also we can change the update function during the program running. The followings are the time estimate in each code on Mac OS X 10.6.7 (Intel Core 2 Duo, 1.6 GHz). 
  (Pure C)
  $ time ./main.exe 

  real  0m0.051s
  user  0m0.031s
  sys 0m0.005s

  (Ruby with Extended Library)
  $ time ruby main.rb 

  real  0m0.744s
  user  0m0.720s
  sys 0m0.008s

  (Pure Ruby)
  $ time ruby dm.rb 

  real  0m3.062s
  user  0m3.045s
  sys 0m0.010s


Reference:
  1. EGillespie, D.T., Exact stochastic simulation of coupled chemical reactions, The journal of physical chemistry, 81(25), 2340--2361, 1977
  2. GSL-GNU Scientific Library