In this notebook, we investigate the critical role of random events in shaping the trajectory of an epidemic during its initial stages when the number of infected individuals is limited. We employ the SIR (Susceptible-Infectious-Recovered) model to simulate the spread of an infectious disease over a 100-day period. By introducing randomness into our model, we aim to uncover how chance events early in an outbreak can significantly influence the epidemic’s course. Our findings shed light on the importance of understanding and accounting for stochastic factors in epidemic modeling

Deterministic model

The following contains code for a simple*deterministic*SIR model, which you will be very familiar with .
In[]:=
betaSoln=First@Solve
ℜ
0

β
γ
,β
Out[]=
{βγ
ℜ
0
}
In[]:=
params=<|
ℜ
0
0.75,ββ,γ0.4,1000001,I01|>/.betaSoln
Out[]=

ℜ
0
0.75,βγ
ℜ
0
,γ0.4,1000001,I01
In[]:=
paramsN=Normal@(params//.params)
Out[]=
{
ℜ
0
0.75,β0.3,γ0.4,1000001,I01}
In[]:=
(*Recoveryrate*)​​γ=0.4;​​(*infectionrate*)​​β=0.3;
In[]:=
tMax1=100;
In[]:=
(*FoI*)​​λ[t_]=βi[t]/n;​​(*totalpopulation*)​​n=1000001;
Define the system of differential equations
eqns={​​s'[t]==-λ[t]*s[t],​​i'[t]==λ[t]*s[t]-γ*i[t],​​r'[t]==γ*i[t],​​s[0]==n-1,​​i[0]==1,​​r[0]==0};
In[]:=
soldeter=NDSolve[eqns,{s,i,r},{t,0,tMax1}];
Plot the solutions
Plot[Evaluate[{s[t],i[t],r[t]}/.soldeter],{t,0,tMax1},PlotLegends->{"Susceptible","Infected","Recovered"},AxesLabel->{"Time","Population"},AxesLabel->{"Time","Population"},PlotStyle->{Blue,Red,Green},PlotLegends->{"Susceptible","Infectious","Recovered"},FrameLabel->{"Time (years)","Prevalence (proportion)"},PlotLabel->"Baseline prevalence of susceptible, infected, and recovered animals over time"]
Out[]=
Susceptible
Infected
Recovered
In[]:=
p=Plot[Evaluate[{i[t]}/.soldeter],{t,0,tMax1},PlotLegends->{"Infected"},AxesLabel->{"Time","Population"},PlotLabel->"Prevalence of Infection",PlotStyle->{Red},PlotRange->All]
Out[]=
Infected
Cumulative incidence function
inc[tmax_]:=NIntegrate[λ[t]s[t]/.soldeter,{t,0,tmax}];
In[]:=
cump=Plot[inc[u],{u,0,100},AxesLabel->{"Time","Incidence"},PlotLegends->{"cumulative incidence"},PlotStyle->{Black,Dashed},PlotRange->All]
Out[]=
cumulative incidence
In[]:=
Show[p,cump,PlotRange->All]
Out[]=
Infected
cumulative incidence
For the deterministic model, you should observe that infection goes extinct without causing an epidemic, since
R
0
<1
. As it is a deterministic model, this output will look the same each time we run this code .

Stochastic Model

In comparison to that from a deterministic model, we are going to explore the stochastic model behaviour for the same value of
R
0
.
Let’s first study the prevalence of infection- the number or proportion of people that are infected at a given timepoint .
In[]:=
(*DefinetheSIRmodelreactions*)reactions={s+i->2i,(*Infection*)i->r(*Recovery*)};
Define initial conditions as an Association
initialConditions=<|s->1000000,i->1,r->0|>;
Define rate constants for the reactions
beta=0.3/999999;(*Transmissionrate(β)*)​​gamma=0.4;(*Recoveryrate(γ)*)
Set up the influx (in this case, there is none
influx=<||>;
Define the time parameters
minTime=0;​​maxTime=100;​​timeStep=1;
(the code for GillespieSSA is in the initialization section at the end of the notebook)
Call the GillespieSSA function
sirSimulations=GillespieSSA[reactions,initialConditions,{beta,gamma},influx,{minTime,maxTime,timeStep}];
Extract the results for each population
{Ssim,Isim,Rsim}=sirSimulations;
In[]:=
numSimulations=10;
Perform 10 simulations of the I compartment and store the results
Plot the 10 simulations of the Infectious (I) population
Additionally, we are also interested in the cumulative incidence of infection -in each iteration, how many people have become infected in total over our simulation period of 100 days?
Generate 10 scenarios of cumulative incidence
Plot the cumulative incidence for each scenario

Initialization

CITE THIS NOTEBOOK

Analyzing subcritical pathogen spread with stochastic SIR simulations (R0 < 1)​
by Athanasios Paraskevopoulos​
Wolfram Community, STAFF PICKS, September 6, 2023
​https://community.wolfram.com/groups/-/m/t/3006228