# Exploring Variability: 100 Stochastic Simulations of Cumulative Incidence

We know $$R_0=0.75$$ and $$\gamma=0.4$$. So with the following code I calculate $$\beta$$

betaSoln =
First@Solve[Subscript[\[GothicCapitalR], 0] == \!$$TraditionalForm\ \*FractionBox[\(\(\$$$$\[Beta]$$\), $$\[Gamma]$$]\), \[Beta]]

params = <|
Subscript[\[GothicCapitalR], 0] ->
0.75, \[Beta] -> \[Beta], \[Gamma] -> 0.4, \[ScriptCapitalN] ->
1000000, I0 -> 1|> /. betaSoln

paramsN = Normal@(params //. params)


Deterministic Model

(*Given parameters*)\[Gamma] = 0.4;
\[Beta] = 0.3;
tMax = 100;
n = 1000000;

(*Define the force of infection (\[Lambda])*)
\[Lambda][t_] := \[Beta] i[t]/n;

(*Define the system of differential equations*)
eqns = {s'[t] == -\[Lambda][t] s[t],
i'[t] == \[Lambda][t] s[t] - \[Gamma] i[t], r'[t] == \[Gamma] i[t],
s[0] == 1000, i[0] == 1, r[0] == 0};

(*Solve the system*)
sol = NDSolve[eqns, {s, i, r}, {t, 0, tMax}][[1]];

(*Cumulative incidence function*)
inc[tmax_] := NIntegrate[\[Lambda][t] s[t] /. sol, {t, 0, tmax}];

cump = Plot[inc[u], {u, 0, 100}, AxesLabel -> {"Time", "Incidence"},
PlotLegends -> {"cumulative incidence"},
PlotStyle -> {Black, Dashed}, PlotRange -> All];

p = Plot[Evaluate[{i[t]} /. sol], {t, 0, tMax},
PlotLegends -> {"Infected"}, AxesLabel -> {"Time", "Population"},
PlotStyle -> {Red}, PlotRange -> All];

Show[p, cump, PlotRange -> All]



So far we don't see anything unusual.

Now I am trying to define the stochastic Model for 100 simulations

(*Define the SIR model reactions*)
reactions = {s + i -> 2 i,(*Infection*)i -> r (*Recovery*)};

(*Define updated initial conditions as an Association*)
initialConditions = <|s -> 999999, i -> 1, r -> 0, cumulativeinc -> 0|>;

(*Define rate constants for the reactions*)
beta = 0.44/999999; (*Transmission rate (\[Beta])*)
gamma = 0.4; (*Recovery rate (\[Gamma])*)

(*Set up the influx (in this case,there is none)*)
influx = <||>;

(*Define the time parameters*)
minTime = 0;
maxTime = 100;
timeStep = 1;

(*Call the GillespieSSA function*)
sirSimulations =
GillespieSSA[reactions, initialConditions, {beta, gamma},
influx, {minTime, maxTime, timeStep}];

(*Extract the results for each population*)
{Ssim, Isim, Rsim} = sirSimulations[[1 ;; 3]];

(*Extract the numerical values from the InterpolatingFunction*)
IsimValues = Isim["ValuesOnGrid"];

(*Calculate the cumulative incidence*)
CumulativeIncidence = Accumulate[IsimValues];

ListLinePlot[CumulativeIncidence, PlotRange -> All,
FrameLabel -> {"Time", "Cumulative Incidence"},
PlotLabel -> "Cumulative Incidence Plot",
AxesStyle -> Directive[Black, 12], ImageSize -> 500]


I am trying to plot 100 simulations but it gives only one

• It is not clear how your stochastic model looks like? Commented Sep 5, 2023 at 14:44
• @AlexTrounev I am trying to solve it with Gillespie Algorith Commented Sep 5, 2023 at 15:30
• Well, we know algorithm. But how stochastic model looks like? Commented Sep 6, 2023 at 2:04

(*Generate 100 scenarios of cumulative incidence*)
numScenarios = 100;
cumulativeIncidenceData =
Table[sirSimulations =
GillespieSSA[reactions, initialConditions, {beta, gamma},
influx, {minTime, maxTime, timeStep}];
{Ssim, Isim, Rsim} = sirSimulations[[1 ;; 3]];
IsimValues = Isim["ValuesOnGrid"];
Accumulate[IsimValues], {numScenarios}];

(*Plot the cumulative incidence for each scenario*)
ListLinePlot[cumulativeIncidenceData, PlotRange -> All,
PlotLegends -> Range[numScenarios],
AxesLabel -> {"Time", "Cumulative Incidence"},
PlotLabel -> "Cumulative Incidence Plot (100 Scenarios)",
AxesStyle -> Directive[Black, 12], ImageSize -> 600]
`

If anyone else has any other ideas I would be interested in sharing them with me