0
$\begingroup$

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]

enter image description here

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]

enter image description here

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

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

1 Answer 1

0
$\begingroup$
(*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]

enter image description here

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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