Gillespie Algorithm (Monte Carlo Simulation) for simple process in Mathematica

In summary: It wasn't until later that mathematicians realized that you could solve the problem by taking derivatives.In summary, the Gillespie algorithm does not converge to a steady state value for RNA using small numbers for the initial conditions. However, if you use larger numbers for the initial conditions, the deterministic and stochastic solutions match perfectly.
  • #1
Lucid Dreamer
25
0
I am trying to model a simple birth and death process in Mathematica using the Gillespie Algorithm.

I am using 1 DNA molecule that is transcribed to mRNA with rate k1,
[tex] \mbox{DNA} \longrightarrow \mbox{DNA + RNA} [/tex]
and the transcribed RNA are subject to degradation with rate k2,
[tex] \mbox{RNA} \longrightarrow \mbox{0} [/tex]
I expect my DNA numbers to remain the same (since DNA is not being produced or degraded), and the steady state RNA number occurs when [tex] \frac{\partial (RNA)}{\partial t}= 0 [/tex]
for,
[tex] \frac{\partial (RNA)}{\partial t} = \mbox{(# of DNA)k1 - (# of RNA)k2}[/tex]
Thus, deterministically, the steady state RNA number is given by [tex] \mbox{(# of DNA)}\frac{\mbox{k1}}{\mbox{k2}} [/tex]
Using the initial conditions # of DNA = 1, k1 = 0.01, k2 = 0.00577, I expect to get a steady state mRNA number of 1.73

I model this process stochastically in Mathematica, but my steady state numbers for the deterministic and stochastic solution do not match. Using the gillespie algorithm, my steady state value is close to 2.3.

The funny thing is that if I start my code with the initial condition # of DNA = 50, k1= 0.01, k2= 0.00577, my deterministic and stochastic solution mach perfectly (in this regime, I expect steady mRNA to be 86.6)

So when I choose small numbers for my initial condition, the deterministic and stochastic solution do not match, but if I choose larger initial conditions, then there is a perfect match between the deterministic and stochastic solution.

Here is the code I used. Any help is appreciated (both in solving this problem and in making the code more efficient)

Code:
Clear[dna0, mRNA0, protein0, time0, dna, mRNA, protein, time, ndna, 
nmRNA, nprotein, ntime, k1, k2, k3, k4, b, n, m, concTime, 
concTimeAvg, rand, randexp]
(*dna -> dna + nMRNA [k1]*)
(*mRNA -> 0 [k2]*)
time0 = 0(*Initial Time*);
dna0 = 1(*Initial number of DNA molecules*);
mRNA0 = 0(*Initial number of mRNA molecules*);
n = 100(*Number of iterations*);
m = 1000(*Number of runs*);
k1 = 0.01(*mRNA production rate*);
k2 = 0.00577(*mRNA degredation rate*);

(*Single iteration*)
concentration[nd_, nm_, t0_, n_] := NestList[
  (
    a1 = #[[1]]*k1;
    a2 = #[[2]]*k2;
    a0 = a1 + a2;
    randexp = 
     a0 // ExponentialDistribution // 
      RandomVariate;(*step time forward*)
    rand = RandomReal[];
    If[rand <= a1/a0,
     {#[[1]], #[[2]] + 1, #[[3]] + 
       randexp},(*add one mRNA molecule*)
     {#[[1]], #[[2]] - 1, #[[3]] + 
       randexp}(*subtract one mRNA molecule*)]
    ) &,
  {nd, nm, t0}, n]

{dna, mRNA, time} = Transpose[concentration[dna0, mRNA0, time0, n]];

(*Average over m runs*)
concentrationAverage[nd_, nm_, t0_, n_, m_] := 
 Table[concentration[nd, nm, t0, n], {m}]
{dnaAverage, mRNAAverage, timeAverage} = 
  concentrationAverage[dna0, mRNA0, time0, n, m] // N // Mean // 
   Transpose;
ListLinePlot[Transpose[{timeAverage, dnaAverage}]]
ListLinePlot[Transpose[{timeAverage, mRNAAverage}], PlotRange -> All]
 
Physics news on Phys.org
  • #2
LucidDreamer,

try to weight your averages by Delta t at each simulation step.
 
  • #3
The reason you aren't getting a match when using small numbers is because the system will then require more sampling, say... infinity? Well, not infinity... but in theory you always need an infinite number of calculations to get true convergence, so you should consider running more calculations to see if you get convergence. I also do monte carlo simulations and the general rule is always do more computing if you aren't sure of the results.
 
  • #4
Hi Lucid Dreamer, welcome to PF!

I ran your code and reproduced your results. It's not a matter of more sampling. It's a matter of applying (continuous) derivatives to try to reproduce a discrete process. When you have small DNA values, the system keeps bumping against mRNA=0, which it can't go below. As a result, the distribution of final values is strongly right-skewed and the mean is larger than one would expect by solving the continuous equation. As the DNA value increases, the approximation of a continuous system gets better and better.

The problem of small numbers is why Gillespie developed his algorithm in the first place!
 
  • #5



The Gillespie algorithm, also known as Monte Carlo simulation, is a powerful tool for modeling stochastic processes in a deterministic system. It allows for the simulation of individual events in a system, taking into account the randomness and variability of these events. In this case, the Gillespie algorithm is being used to model a simple birth and death process in Mathematica.

The process being modeled involves the transcription of DNA into mRNA and the subsequent degradation of the mRNA. The rates of these events, k1 and k2 respectively, are known and the initial conditions for the number of DNA molecules and the production and degradation rates are provided.

The deterministic solution for this process predicts a steady state number of mRNA molecules, but when using the Gillespie algorithm, the steady state value obtained is different. This difference is likely due to the inherent randomness in the simulation and the small initial condition chosen.

To improve the accuracy of the simulation, it may be helpful to increase the number of iterations and runs, as well as to optimize the code for efficiency. Additionally, it may be useful to compare the results of the stochastic simulation to other methods of modeling the same process to ensure the accuracy of the Gillespie algorithm.
 

What is the Gillespie Algorithm?

The Gillespie Algorithm, also known as the Monte Carlo Simulation, is a computational method used to simulate and analyze complex systems that involve random events over time. It is commonly used in fields such as physics, chemistry, biology, and economics.

What is the purpose of using the Gillespie Algorithm?

The Gillespie Algorithm is used to model and analyze systems that involve stochastic processes. It allows scientists to generate probabilistic predictions and understand the behavior of these systems over time.

How does the Gillespie Algorithm work?

The Gillespie Algorithm works by simulating the behavior of a system by generating random numbers to represent the occurrence of events. These random numbers are compared to the probabilities of different events happening, and the system is updated accordingly. This process is repeated until a desired time has been reached or a specific outcome has been achieved.

What are the advantages of using the Gillespie Algorithm?

One of the main advantages of using the Gillespie Algorithm is its ability to model and analyze non-linear and complex systems. It also allows for the incorporation of stochasticity, which is important in many real-world systems. Additionally, the algorithm can easily be implemented in computer programs, making it a widely accessible tool for scientists.

What are some common applications of the Gillespie Algorithm?

The Gillespie Algorithm is commonly used in many fields of science, such as ecology, epidemiology, and systems biology. It has also been used in economics to model financial markets and in engineering for reliability analysis. Additionally, the algorithm has been applied to study chemical reactions, genetic networks, and evolution.

Back
Top