Team:Tsinghua-A/Modeling/GILLESPIE/part1

A more accurate model using Gillespie algorithm

To make the model more convincing, we use Gillespie algorithm, which is a kind of stochastic simulation algorithm to generate a statistically correct trajectory (possible solution) of a stochastic equation.

Summary of the steps to run the algorithm

“Step 0 (Initialization). Input the desired values for the M reaction constants c_1,. . ., c_Mand the N initial molecular population numbers X_1,. . ., X_M. Set the time variable t and the reaction counter n both to zero. Initialize the unit-interval uniform random number generator (URN).
Step 1. Calculate and store the M quantities a_1 = h_1 c_1,, . , , a_M = h_M c_M for the current molecular population numbers, where h_v is that function of X_1,. . ., X_M defined in (1). Also calculate and store as a_0 the sum of the M a_v values.
Step 2. Generate two random numbers r1 and r2 using the unit-interval uniform random number generator, and calculateτandμaccording to (2) and (3).
Step 3. Using theτandμvalues obtained in step 2, increase t byτ, and adjust the molecular population levels to reflect the occurrence of one R_μ, reaction; e.g., if R_μ is the reaction in (4), then increase X_1 by 1 and decrease X_2 by 1. Then increase the reaction counter n by 1 and return to step 1. “(Danlel T. Gillespie, 1977, p. 2345)[7]

We define for each reaction R_μ and h_μ according to the equations below.

Schematic of the stochastic simulation al algorithm (Danlel T. Gillespie, 1977)[7]

Simplified summary of the steps to run the algorithm

1. Initialization:
Initialize the number of molecules in the system, reactions constants, and random number generators.
2. Monte Carlo step:
Generate random numbers to determine the next reaction to occur as well as the time interval. The probability of a given reaction to be chosen is proportional to the number of substrate molecules.
3. Update:
Increase the time step by the randomly generated time in Step 2. Update the molecule count based on the reaction that occurred.
4. Iterate:
Go back to Step 2 unless the number of reactants is zero or the simulation time has been exceeded.