# Stochastic Simulation using the Gillespie algorithm

I'm trying to reproduce the simulation with demographic stochasticity in Figure 1 from the paper entitled "Dynamical Resonance Can Account For Seasonality of Influenza Epidemics" (http://www.pnas.org/content/101/48/16915.full).

My code works just fine when $t$ is small $(0<t<0.02)$, but when I increase the value of $t$ (as below), the simulation doesn't finish compiling. Thus, I'm wondering if there's an issue with my code. Ideally, I would like my code to function when $(10<t<20)$ so I can compare my results to the paper. Any input would be greatly appreciated.

Clear[S, Infe, R, No, b0, b1, b, t, Ro, d, L, event, fun, main, list1, list2, data]
No = 500000;
b0 = 500;
b1 = 0.04;
d = 0.02;
L = 4;
list1 = {};
list2 = {};
b[t_] := b0 *(1 + b1* Cos[2 Pi t]);
main[] := Module[{S = 499999, Infe = 1, R = 0, t = 0},
While[t < .03 && Infe > 0, {S, Infe, R, t} = fun[S, Infe, R, t];
AppendTo[list1, Infe]; AppendTo[list2, t] ]];

fun[S0_, Infe0_, R0_, t0_] :=
Module[ {S = S0, Infe = Infe0, R = R0, t = t0},
TransRate = (b[t]*S*Infe)/No;
RecoveryRate = Infe/d;
ImmunityLoss = (No - S - Infe)/L;
TotalRate = TransRate + RecoveryRate + ImmunityLoss;

t = t + RandomReal[ExponentialDistribution[TotalRate]];
event = RandomReal[UniformDistribution[]];

If[event < TransRate/TotalRate,
S = S - 1 ; Infe = Infe + 1,
If[event < (TransRate + RecoveryRate)/TotalRate,
Infe = Infe - 1 ; R = R + 1, R = R - 1 ; S = S + 1]];
{S, Infe, R, t}
]

main[]
data = Transpose@{list2, list1};
ListPlot[data]

-
The value of t approaches your stopping limit of .03 from below slowly. The repeated Append is killing you. Change that operation to a Sow and Reap operation- that will speed things greatly (as in order of magnitude +). Then take a look at the top FAQ questions here to refactor this into a more functional/Mathematica style - the imperative style here is seldom efficient. Lastly, consider just doing it as a Markov chain perhaps - Mathematica has built-in capabilities for these that are quite efficient. –  rasher Jul 21 at 23:19
Also, consider pre-generating the random variates and indexing into them rather than the one-at-a-time way you're doing it now: in cases where there will be many iterations, the set-up/tear-down cost of making individual calls to Random... will get expensive. –  rasher Jul 21 at 23:26
A last thought (I don't have time today to re-write your example incorporating my thoughts above): for your desired range of 10<t<20 you'll on average iterate a quarter to a half billion times. Better to generate the variates in chunks, checking the running sum until your t condition is met, and vectorizing the remaining calculations where possible. –  rasher Jul 21 at 23:45
Thank you Rasher! Your advice helped tremendously. –  Omar Jul 23 at 2:12