5
$\begingroup$

I am trying to implement Michaelis–Menten kinetics using Tau-Leaping. It works fine as expected. Only issue here is when substrate reaches zero the algorithm stop since PoissonDistribution[0] gives error. What I want is if #[[3]] = 0, then set #[[3]] = 1 and run until #[[4]]>0. I don't know how to use if condition here. Any suggestion? Thanks.

SeedRandom[1234];
With[{τ = 0.001, k1 = 1.1, k2 = 0.1, k3 = 0.8, e = 200, s = 200, c = 1, p = 1},
  sim = 
    NestWhileList[
      ({Δp1, Δp2, Δp3} = 
          {RandomVariate @ PoissonDistribution[τ k1  #[[2]] #[[3]]], 
           RandomVariate @ PoissonDistribution[τ k2  #[[4]]], 
           RandomVariate @ PoissonDistribution[τ k3  #[[4]]]};
        {#[[1]] + τ, 
         #[[2]] - Δp1 + Δp2 + Δp3, 
         #[[3]] - Δp1 + Δp2, 
         #[[4]] + Δp1 - Δp2 - Δp3, #[[5]] +Δp3}) &, 
    {0, e, s, c, p}, 
    #[[3]] > 0 &]];

{t, e, s, c, p} = Transpose @ sim;

ListStepPlot[
  {Transpose @ {t, e}, Transpose @ {t, s}, Transpose @ {t, c}, Transpose@{t, p}}, 
  Frame -> True, PlotTheme -> "Detailed", 
  FrameLabel -> {"Time", "Population"}, ImageSize -> Large, 
  PlotLegends -> {"Enzyme", "Substrate", "Complex", " Product"}]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Can you provide a SeedRandom that exhibits the failure? I can't reproduce the problem. (Not that I could solve the issue if I could reproduce the problem....) $\endgroup$ – bobthechemist Nov 16 '17 at 2:18
  • $\begingroup$ Thanks for the comments. I added SeedRandom. The code produce the image that I have attached but it stop in a short time before complex reach to zero. $\endgroup$ – OkkesDulgerci Nov 16 '17 at 3:20
7
$\begingroup$

Hack #1: add 10^-100 to the argument of PoissonDistribution. This won't appreciably affect the rates, but will prevent PoissonDistribution[0] from occurring.

This reveals another potential problem: reactions can use more than the number of molecules that exist, leading to negative numbers, which then leads to problems with negative arguments to PoissonDistribution.

Hack #2: use Min to make sure reactions don't consume more reactants than exist. (Just thought of this one, please comment if you see any problems with it).

Together:

With[{τ = 0.001, k1 = 1.1, k2 = 0.1, k3 = 0.8, e = 200, s = 200, 
  c = 1, p = 1, eps = 10^-100}, 
  sim = NestWhileList[({Δp1, Δp2, Δp3} = {
    Min[{RandomVariate@PoissonDistribution[eps + τ k1 #[[2]] #[[3]]],
      #[[2]], #[[3]]}],
    Min[{RandomVariate@PoissonDistribution[eps + τ k2 #[[4]]],
      #[[4]]}],
    Min[{RandomVariate@PoissonDistribution[eps + τ k3 #[[4]]],
      #[[4]]}]};
    {#[[1]] + τ, #[[2]] - Δp1 + Δp2 + Δp3, #[[3]] - Δp1 + Δp2,
     #[[4]] + Δp1 - Δp2 - Δp3, #[[5]] + Δp3}) &,
    {0, e, s, c, p}, #[[4]] > 0 &]
];

{t, e, s, c, p} = Transpose@sim;

ListStepPlot[{Transpose@{t, e}, Transpose@{t, s}, Transpose@{t, c}, 
  Transpose@{t, p}}, Frame -> True, PlotTheme -> "Detailed", 
  FrameLabel -> {"Time", "Population"}, ImageSize -> Large, 
  PlotLegends -> {"Enzyme", "Substrate", "Complex", " Product"}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ It works like a charm!!! Thanks.. $\endgroup$ – OkkesDulgerci Nov 16 '17 at 3:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.