Gillespie Stochastic Simulation Algorithm - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-08-19T09:01:59Z https://mathematica.stackexchange.com/feeds/question/119786 http://www.creativecommons.org/licenses/by-sa/3.0/rdf https://mathematica.stackexchange.com/q/119786 23 Gillespie Stochastic Simulation Algorithm István Zachar https://mathematica.stackexchange.com/users/89 2016-07-01T06:58:10Z 2017-02-11T16:46:26Z <p>The <a href="https://en.wikipedia.org/wiki/Gillespie_algorithm" rel="noreferrer">Gillespie SSA</a> is a Monte Carlo stochastic simulation algorithm to find the trajectory of a dynamic system described by a reaction (or interaction) network, e.g. chemical reactions or ecological problems. It was introduced by Dan Gillespie in 1977 (<a href="http://pubs.acs.org/doi/abs/10.1021/j100540a008" rel="noreferrer">see paper here</a>). It is used in case of small molecular numbers (or species abundance) where numerical integration of the related differential equation system is not appropriate due to hard stochastic effects (i.e. the death of a single individual might make a large impact on the population).</p> <p>Can you do it with <em>Mathematica</em>? Is it general enough? Is it fast?</p> https://mathematica.stackexchange.com/questions/119786/-/119788#119788 32 Answer by István Zachar for Gillespie Stochastic Simulation Algorithm István Zachar https://mathematica.stackexchange.com/users/89 2016-07-01T07:15:45Z 2016-07-02T16:45:09Z <p>Yes you can. Below is a fairly general, <em>Mathematica</em>-compiled, fast and robust version.</p> <h2>Examples</h2> <p><strong>1. Michaelis-Menten kinetics</strong></p> <p><a href="https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics">Michaelis-Menten kinetics</a> for enzyme-directed substrate conversion. The enzyme (<code>e</code>) converts the susbtrate (<code>s</code>) through an enzyme-substrate complex (<code>c</code>) to the product (<code>p</code>). For comparison, I've included the deterministic ODE system solved by <code>NDSolve</code>.</p> <pre><code>ClearAll[e, s, c, p, t]; reactions = {e + s -&gt; c, c -&gt; e + s, c -&gt; e + p}; vars = {e, s, c, p}; rates = {1.1, .1, .8}; init = &lt;|e -&gt; 100, s -&gt; 100, c -&gt; 0, p -&gt; 0|&gt;; det = NDSolveValue[{ e'[t] == .9 c[t] - 1.1 e[t] s[t], s'[t] == .1 c[t] - 1.1 e[t] s[t], c'[t] == 1.1 e[t] s[t] - .9 c[t], p'[t] == .8 c[t], e == 100, s == 100, c == 0, p == 0}, vars, {t, 0, 10}]; sto = GillespieSSA[reactions, init, rates, {0, 10}]; op = {PlotStyle -&gt; Thick, PlotTheme -&gt; "Scientific"}; Row@{Plot[Evaluate@Through@det@t, {t, 0, 10}, Evaluate@op, PlotLabel -&gt; "deterministic ODE"], Spacer@10, Plot[Evaluate@Through@sto@t, {t, 0, 10}, Evaluate@op, PlotLabel -&gt; "stochastic SSA"]} </code></pre> <p><img src="https://i.stack.imgur.com/oCSUd.png" alt="Mathematica graphics"></p> <p><strong>2. Lotka-Volterra predator-prey dynamics</strong></p> <p><a href="https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations">Lotka-Volterra dynamics</a>. I omit the conversion to continuous differential equations from now on, leaving it to the "educated reader". <code>x -&gt; Null</code> indicates that the species is removed without producing waste material (at least one that is tracked). Similarly, <code>Null -&gt; x</code> indicates a zero-order reaction where <code>x</code> is generated spontaneously (or is entering from the external environment).</p> <pre><code>ClearAll[y, x]; reactions = {y -&gt; 2 y, y + x -&gt; 2 x, x -&gt; Null}; vars = {y, x}; rates = {1, .005, .6}; init = &lt;|y -&gt; 50, x -&gt; 100|&gt;; sto = GillespieSSA[reactions, init, rates, {0, 100}] </code></pre> <p><img src="https://i.stack.imgur.com/bXm3w.png" alt="Mathematica graphics"></p> <p><strong>3. Circadian cycle</strong></p> <pre><code>ClearAll[a, p, i]; reactions = {a -&gt; 2 a, a -&gt; a + p, p -&gt; i, a + i -&gt; i, i -&gt; Null}; vars = {a, p, i}; init = &lt;|a -&gt; 100, p -&gt; 100, i -&gt; 100|&gt;; rates = {1, .08, .6, .01, .4}; sto = GillespieSSA[reactions, init, rates, {0, 400}]; </code></pre> <p><img src="https://i.stack.imgur.com/QV6fC.png" alt="Mathematica graphics"></p> <p><strong>4. Oregonator</strong></p> <p>The <a href="http://www.scholarpedia.org/article/Oregonator">Oregonator</a> is a model of the <a href="https://en.wikipedia.org/wiki/Belousov%E2%80%93Zhabotinsky_reaction">Belousov-Zhabotinsky</a> reactions. Here I used the <code>resolution</code> argument to speed up evaluation (only every 100th step is stored).</p> <pre><code>ClearAll[a, b, c]; reactions = {b -&gt; a, a + b -&gt; Null, a -&gt; 2 a + c, 2 a -&gt; Null, c -&gt; b}; vars = {a, b, c}; init = &lt;|a -&gt; 1, b -&gt; 2, c -&gt; 3|&gt; rates = {2, .1, 104, .016, 26}; sto = GillespieSSA[reactions, init, rates, {0, 4, 100}]; </code></pre> <p><img src="https://i.stack.imgur.com/IeT6Y.png" alt="Mathematica graphics"></p> <h2>Usage</h2> <p>The function accepts the following arguments:</p> <pre><code> GillespieSSA[ {r1, r2, ...}, (* list of elementary reaction steps as Rules *) &lt;|y1 -&gt; y1, y2 -&gt; y2, ...|&gt;, (* variables with initial values at t==0 *) {c1, c2, ...}, (* reaction rate constants, for each reaction *) &lt;|y1 -&gt; f1, y2 -&gt; f2, ...|&gt;, (* linear in/outflux for each variable; can be left out, in which case 0 will be used*) {mint, maxt, res} (* {start time, end time, step resolution}; resolution can be omitted, defaulting to 1 *) ] </code></pre> <p>It returns an <code>InterpolatingFunction</code> (with <code>InterpolationOrder -&gt; 0</code>) as each variable's solution function, to comply with the result produced by <code>NDSolve</code>. Initial values are taken to be the values at <code>t == mint</code>. The maximal allowed step size (10^7) is hardcoded in <code>iterations</code>, change it for your needs.</p> <p>Note, that the Gillespie method is stochastic: it will convert continuous reaction rates into discrete-valued reaction propensities. It cannot accept reactions where stoichiometric factors are not integer numbers. Also, it gives a different realization every time it is run (if the random generator is not reseeded) due to its stochastic nature. Moreover, at small molecular/species amounts, stochastic effects could cause extinction and produce different results as in the deterministic, continuous case (<code>NDSolve</code>).</p> <h2>Code</h2> <pre><code>ClearAll[GillespieSSA]; GillespieSSA[res : {__Rule}, in_Association, rateconst_?VectorQ, influx_Association: &lt;||&gt;, {mint_?NumberQ, maxt_?NumberQ, dstep_Integer: 1}] := Module[ {vars, reactant, product, balance, propensities, initialValues, fluxRates, symRates, rep, stepList, iterations = 10^7, step, compiled, times, rest}, (* Pre-generating a list is much faster than iteratively calling one-by-one. *) stepList = N@RandomVariate[ExponentialDistribution@1, iterations + 1]; {vars, initialValues} = {Keys@in, Values@in}; {reactant, product} = Outer[Coefficient[#1, #2] &amp;, #, vars, 1] &amp; /@ Transpose@(List @@@ res); balance = product - reactant; propensities = Inner[Binomial[#2, #1] &amp;, reactant, vars, Times]* PadRight[rateconst, Length@res, 1]; fluxRates = If[influx === &lt;||&gt;, 0 &amp; /@ vars, vars /. influx]; Block[{count}, rep = Thread[vars -&gt; Table[Indexed[count, i], {i, Length@vars}]]; symRates = propensities /. rep; compiled = ReleaseHold[ Hold@Compile[{ {init, _Integer, 1}, {flux, _Integer, 1}, {bal, _Integer, 2}, {dtList, _Real, 1}, {min, _Real}, {max, _Real}, {iter, _Integer}, {resol, _Integer}}, Module[{ count = init, rates, i = 1, c = 1, t = min, dt, ff, f, range, r, fReal = 0. &amp; /@ flux, rateSum, data = Internal`Bag[]}, rates = "SymbolicRates"; rateSum = N@Total@rates; range = Range@Length@rates; Internal`StuffBag[data, Internal`Bag[Join[{t}, N@count]]]; While[Total@count &gt; 0. &amp;&amp; rateSum &gt; 0. &amp;&amp; t &lt;= max &amp;&amp; i &lt;= iter, i++; dt = dtList[[i]]/rateSum; t = t + dt; r = RandomChoice[rates -&gt; range]; (* Fractional part is carried over to minimize undersampling error *) ff = (flux*dt) + fReal; f = IntegerPart@ff; fReal = ff - f; (* `count` is maintained as an integer not to loose precision. *) count = Max[0, #] &amp; /@ (count + bal[[r]] + f); If[Mod[i, resol] == 0, c++; Internal`StuffBag[data, Internal`Bag[Join[{t}, N@count]]];]; rates = "SymbolicRates"; rateSum = N@Total@rates; ]; If[t &lt; max &amp;&amp; i &lt; iter, c++; Internal`StuffBag[data, Internal`Bag[Join[{max}, N@count]]]]; Table[Internal`BagPart[Internal`BagPart[data, j], All], {j, c}] ], Parallelization -&gt; True , RuntimeAttributes -&gt; Listable, RuntimeOptions -&gt; "Speed", CompilationOptions -&gt; {"InlineExternalDefinitions" -&gt; True, "InlineCompiledFunctions" -&gt; True} ] /. "SymbolicRates" -&gt; symRates]; {times, rest} = {First@#, Rest@#} &amp;@Transpose@compiled[ Round@initialValues, N@fluxRates, Round@balance, N@stepList, N@mint, N@maxt, Round@iterations, dstep]; Interpolation[Transpose@{times, #}, InterpolationOrder -&gt; 0] &amp; /@ rest ]]; </code></pre>