I am trying to implement Michaelis–Menten kinetics using Gillespie Algorithm.
$$\text{Reaction 1:}\quad E+S\rightarrow C$$
$$\text{Reaction 2:}\quad C\rightarrow E+S$$ $$\text{Reaction 1:}\quad C\rightarrow E+P$$
Each reaction has a rate with a weight
$$\text{Reaction rate 1:}\quad k_1 ES$$
$$\text{Reaction rate 2:}\quad k_2 C$$ $$\text{Reaction rate 3:}\quad k_3C$$
and time is distributed exponentially. $$t\sim Exp(rate)$$
Minimun t means that reaction occures.
It kinda works for two reactions. How can I modify the code for 3 or more reactions. I have seen this code but I would like to learn coding and algorithm.
{k1 = 1.1, k2 = 0.1, k3 = 0.8, e = 100, s = 100, c = 1, p = 1,
numOfReaction = 100, numOfsim = 10};
sim = NestList[(
Δt1 =
RandomVariate@ExponentialDistribution[k1 #[[2]] #[[3]]];
Δt2 =
RandomVariate@ExponentialDistribution[k2 #[[4]]];
Δt = Min[Δt1, Δt2];
If[Δt1 < Δt2, {#[[
1]] + Δt, #[[2]] - 1, #[[3]] - 1, #[[4]] +
1}, {#[[1]] + Δt, #[[4]] - 1, #[[2]] +
1, #[[3]] + 1}]) &, {0, e, s, c}, numOfReaction];
ListStepPlot[{Transpose@{sim[[All, 1]], sim[[All, 2]]},
Transpose@{sim[[All, 1]], sim[[All, 4]]}},
PlotLegends -> {"Simulation"},
PlotStyle -> Directive[AbsoluteThickness[0.2]], Frame -> True,
PlotTheme -> "Detailed", FrameLabel -> {"Time", "Population"},
ImageSize -> Large]
I have tried this code but it does not work.
NestList[(
Δt1 =
RandomVariate[ExponentialDistribution[k1 #[[2]] #[[3]]]];
Δt2 =
RandomVariate[ExponentialDistribution[k2 #[[4]]]];
Δt3 =
RandomVariate[ExponentialDistribution[k3 #[[4]]]];
Δt =
Min[{Δt1, Δt2, Δt3}];
If[Δt1 < Δt2, {#[[
1]] + Δt, #[[2]] - 1, #[[3]] - 1, #[[4]] +
1}, {#[[4]] - 1, #[[2]] + 1, #[[3]] + 1}];
If[Δt1 < Δt3, {#[[
1]] + Δt, #[[2]] - 1, #[[3]] - 1, #[[4]] +
1}, {#[[4]] - 1, #[[2]] + 1, #[[5]] + 1}];
If[Δt2 < Δt3, {#[[
1]] + Δt, #[[4]] - 1, #[[2]] + 1, #[[3]] +
1}, {#[[4]] - 1, #[[2]] + 1, #[[4]] + 1}]) &, {0, e, s,
c}, numOfReaction]