I am trying to model a simple birth and death process in Mathematica using the Gillespie Algorithm.(adsbygoogle = window.adsbygoogle || []).push({});

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 (Text):

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 Forums | Science Articles, Homework Help, Discussion**

Dismiss Notice

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

Loading...

Similar Threads for Gillespie Algorithm Monte |
---|

MATLAB Sample Code for Monte Carlo estimation of pi |

**Physics Forums | Science Articles, Homework Help, Discussion**