[Gillespie algorithm]

---

[Gillespie algorithm]

[Ref] https://en.wikipedia.org/wiki/Gillespie_algorithm

William Feller establishes that the time-to-the-next-jump was exponentially distributed and the probability of the next event is proportional to the rate.

[Example]

Two reaction rules:

(1) A->0 with rate \alpha

(2) A ->A + A with rate \beta

Total reaction rate tot_rate = n*\alpha + n*\beta, where n is the abundance or concentration of particle A.

[step1] Throw dice and determine which event will happen

Reaction one (1) happens n*\alpha/tot_rate, reaction two (2) happens 1-n*\alpha/tot_rate.

[step2] Time update

Time is updated by t+=dt. 

dt is sampled from the exponential function with mean 1/tot_rate, Exp[-dt*tot_rate]*tot_rate.

[step3] Iterate step1-2.

---

[The other way to implement (MC method)]

[Ref] Su-Chan Park, PRE 72, 036111 (2005)

Two reaction rules:

(1) A->0 with rate \alpha

(2) A ->A + A with rate \beta

[step1] Throw dice and determine which interaction will happen

The number of ways to happen reaction one (1) is n, reaction two (2) also n.

We call the total number of interaction ways, n+n = 2n, the number of accessible states M.

Among these two reactions, one is chosen by throwing dice.

In this case reaction (1) will be chosen with half, and the other half reaction (2) will be chosen.

[step2] examine it really happens or not

if reaction (1) is chosen, the reaction occurs with probability \alpha*dt,

if reaction (2) is chosen, the reaction occurs with probability \beta*dt.

The choice of dt independent is from configurations, but it has to satisfy the following condition:

(\alpha*dt<=1 and \beta*dt<=1 ) (since [rate*dt = prob])

[step3] Time update

Time is updated by t+=dt/M. 

[step4] Iterate step1-3.

My interpretation behind the algorithm

Rate = Prob./dt, and here rate is per capita.

During dt, I assume that all possible accessible states M does not change much.

After M updates, t becomes t+dt.

Since all accessible states are once asked during dt on average

that it occurs or not with probability rate*dt, all reactions occur with their rates.