Epidemiological models
(an introductory notebook by prof. Paolo Rossi, Mathematics Department, U. of Padova)

SIS model (with demographics)

Epidemiological Parameters

T=365;(*final time*)
T0=0;(*initial time*)
pop=60*(10^6);(*total population*)
Infected0=1600;(*infected at t=0*)
Ro0=1.36;(*approximate basic reproduction number (Subscript[R, 0])*)
ItoS=7;(*disease duration*)
μ=0*9.40*(10^(-3))*(365^(-1));(*normal death and birth rate*)
S0=1-I0; 
I0=Infected0/pop;
γ=1/ItoS;
λ=Ro0*γ;

s=NDSolve[{II'[t]==-λ II[t]^2+(λ-γ-μ)II[t],II[0]==I0},II[t],{t,0,T}];
Plot[pop*Evaluate[II[t]/.s],{t,0,T}, PlotRange->{0,pop},PlotLegends->{"Infected"},PlotLabel->"SIS MODEL\n Basic reproduction number (Subscript[R, 0]) = "<>ToString[N[λ/γ]]<>"\n Infection duration: "<>ToString[ Evaluate[(1/γ)]]<> " days\n Contact rate: "<>ToString[ λ]<>" per day\n",ImageSize->800]

Case of yearly periodic variation in contact rate: λ=λ(t)

Epidemiological Parameters

T=5*365;(*final time*)
T0=0;(*initial time*)
pop=60*(10^6);(*total population*)
Infected0=1600;(*infected at t=0*)
Ro0=1.36;(*approximate basic reproduction number (Subscript[R, 0])*)
ItoS=7;(*disease duration*)
A=0.2; (*intensity of seasonal variation*)
μ=0*9.40*(10^(-3))*(365^(-1));(*normal death and birth rate*)
S0=1-I0; 
I0=Infected0/pop;
γ=1/ItoS;
λ=Ro0*γ;

s=NDSolve[{II'[t]==-λ(1-A Sin[(π/365)t]^2)II[t]^2+(λ-γ-μ)II[t],II[0]==I0},II[t],{t,0,T}];
Plot[pop*Evaluate[II[t]/.s],{t,0,T}, PlotRange->{0,pop},PlotLegends->{"Infected"},PlotLabel->"SIS MODEL\n Contact rate: "<>ToString[ λ]<>" per day\n Average infection duration: "<>ToString[ Evaluate[(1/γ)]]<> " days",ImageSize->800]

SIR model

SIR model with demographics

Epidemiological Parameters

T=3/12*365;(*final time*)
T0=0;(*initial time*)
pop=60*(10^6);(*total population*)
Infected0=1600;(*infected at t=0*)
Ro01=10.;(*approximate basic reproduction number (Subscript[R, 0]) option 1*)
Ro02=3.;(*approximate basic reproduction number (Subscript[R, 0]) option 2*)
ItoR=7;(*disease duration*)
μ=0*9.40*(10^(-3))*(365^(-1));(*normal death and birth rate*)
S0=1-I0; 
I0=Infected0/pop;
γ=1/ItoR;
λ1=Ro01*γ;
λ2=Ro02*γ;

s1=NDSolve[{S'[t]==-λ1 S[t]II[t]+μ-μ S[t],II'[t]==λ1 S[t]II[t]+(-γ-μ)II[t],S[0]==S0,II[0]==I0},{S[t],II[t]},{t,0,T}];
s2=NDSolve[{S'[t]==-λ2 S[t]II[t]+μ-μ S[t],II'[t]==λ2 S[t]II[t]+(-γ-μ)II[t],S[0]==S0,II[0]==I0},{S[t],II[t]},{t,0,T}];
Plot[Evaluate[pop*{II[t]/.s1,(1-S[t]-II[t])/.s1,II[t]/.s2,(1-S[t]-II[t])/.s2}],{t,0,T},PlotRange->{0,pop},PlotLegends->{"Infected      ((Subscript[R, 0]) = "<>ToString[N[Ro01]]<>")", "Immunized (\!\(\*SubscriptBox[\(R\), \(\(\\ \)\(0\)\)]\) = "<> ToString[N[Ro01]]<>")","Infected      ((Subscript[R, 0]) = "<>ToString[N[λ2/γ]]<>")", "Immunized ((Subscript[R, 0]) = "<> ToString[N[λ2/γ]]<>")"},PlotLabel->"SIR MODEL\n Average infection duration: "<>ToString[ Evaluate[(1/γ)]]<> " days",ImageSize->800]

Effects of social distancing

Epidemiological Parameters

T=10/12*365;(*final time*)
T0=0;(*initial time*)
pop=60*(10^6);(*total population*)
Infected0=60;(*infected at t=0*)
Ro0=10.;(*approximate basic reproduction number (Subscript[R, 0])*)
ItoR=20;(*disease duration*)
tsd1=23; (*social distancing starts option 1*)
tsdf1=100; (*social distancing ends option 1*)
B1=0.; (*intensity of social distancing option 1*)
tsd2=23; (*social distancing starts option 2*)
tsdf2=100; (*social distancing ends option 2*)
B2=0.91; (*intensity of social distancing option 2*)
μ=0*9.40*(10^(-3))*(365^(-1));(*normal death and birth rate*)
S0=1-I0; 
I0=Infected0/pop;
γ=1/ItoR;
λ=Ro0*γ;

s1=NDSolve[{S'[t]==-λ (1-B1(UnitStep[t-tsd1]-UnitStep[t-tsdf1]))S[t]II[t]+μ-μ S[t],II'[t]==λ (1-B1 (UnitStep[t-tsd1]-UnitStep[t-tsdf1])) S[t]II[t]+(-γ-μ)II[t],S[0]==S0,II[0]==I0},{S[t],II[t]},{t,0,T}];
s2=NDSolve[{S'[t]==-λ (1-B2 (UnitStep[t-tsd2]-UnitStep[t-tsdf2]))S[t]II[t]+μ-μ S[t],II'[t]==λ (1-B2 (UnitStep[t-tsd2]-UnitStep[t-tsdf2])) S[t]II[t]+(-γ-μ)II[t],S[0]==S0,II[0]==I0},{S[t],II[t]},{t,0,T}];
Plot[Evaluate[pop*{II[t]/.s1,II[t]/.s2}],{t,0,T},PlotRange->{0,pop},PlotLegends->{"Infected (" <> ToString[Evaluate[B1*100]]<>"% drop in contact rate from day "<> ToString[tsd1]<>" to day "<>ToString[Evaluate[tsdf1]] <>")","Infected (" <> ToString[Evaluate[B2*100]]<>"% drop in contact rate from day "<> ToString[tsd2]<>" to day "<>ToString[Evaluate[tsdf2]] <>")"},PlotLabel->"SIR MODEL\n (Subscript[R, 0]) = "<>ToString[N[λ/γ]]<>"\n Infection duration: "<>ToString[ (1/γ)]<> " days\n Contact rate: "<>ToString[N[λ]]<>" per day",ImageSize->800]

SIRS model

SIRS model with demographics

Epidemiological Parameters

T=6*365;(*final time*)
T0=0;(*initial time*)
pop=60*(10^6);(*total population*)
Infected0=60;(*infected at t=0*)
Ro0=4.;(*approximate basic reproduction number (Subscript[R, 0])*)
ItoR=20.;(*disease duration*)
RtoS=600; (*immunization duration*)
μ=0*9.40*(10^(-3))*(365^(-1));(*normal death and birth rate*)
S0=1-I0; 
I0=Infected0/pop;
γ=1/ItoR;
λ=Ro0*γ;
ξ=1/RtoS;

s1=NDSolve[{S'[t]==-λ S[t]II[t]+μ-μ S[t]+ξ(1-S[t]-II[t]),II'[t]==λ S[t]II[t]+(-γ-μ)II[t],S[0]==S0,II[0]==I0},{S[t],II[t]},{t,0,T}];
Plot[Evaluate[pop*{II[t]/.s1,(1-S[t]-II[t])/.s1}],{t,0,T},PlotRange->{0,pop},PlotLegends->{"Infected", "Immunized"},PlotLabel->"SIRS MODEL\n (Subscript[R, 0]) = "<>ToString[N[λ/γ]]<>"\n Infection duration: "<>ToString[ (1/γ)]<> " days\n Contact rate: "<>ToString[N[λ]]<>" per day\n immunization duration: "<>ToString[ Evaluate[(1/ξ)]]<> " days",ImageSize->800]

Case of yearly periodic variation in contact rate: λ=λ(t) (the good season makes it better, but also causes a recrudescence during the next winter)

Epidemiological Parameters

T=6*365;(*final time*)
T0=0;(*initial time*)
pop=60*(10^6);(*total population*)
Infected0=60;(*infected at t=0*)
Ro0=4.;(*approximate basic reproduction number (Subscript[R, 0])*)
ItoR=20.;(*disease duration*)
RtoS=280; (*immunization duration*)
A=0.4; (*intensity of seasonal variation*)
μ=0*9.40*(10^(-3))*(365^(-1));(*normal death and birth rate*)
S0=1-I0; 
I0=Infected0/pop;
γ=1/ItoR;
λ=Ro0*γ;
ξ=1/RtoS;

s1=NDSolve[{S'[t]==-(λ(1-A Sin[(π/365)t]^2))S[t]II[t]+μ-μ S[t]+ξ(1-S[t]-II[t]),II'[t]==(λ(1-A Sin[(π/365)t]^2)) S[t]II[t]+(-γ-μ)II[t],S[0]==S0,II[0]==I0},{S[t],II[t]},{t,0,T}];
Plot[Evaluate[pop{II[t]/.s1,(1-S[t]-II[t])/.s1}],{t,0,T},PlotRange->{0,pop},PlotLegends->{"Infected", "Immunized","Seasonal modulation of contact rate"},PlotLabel->"SIRS MODEL\n (Subscript[R, 0]) = "<>ToString[N[λ/γ]]<>"\n Infection duration: "<>ToString[ (1/γ)]<> " days\n Contact rate: "<>ToString[N[λ]]<>" per day\n immunization duration: "<>ToString[ Evaluate[(1/ξ)]]<> " days",ImageSize->800]

SEIRS model with increased mortality for infected, seasonal variation and social distancing (activated at given time or after a given number of deaths)

Epidemiological Parameters

T=24/12*365;(*final time*)
T0=0;(*initial time*)
pop=60*(10^6);(*total population*)
A=0.; (*intesity of seasonal modulation*)
B1=0.85; (*intesity of social distancing by deaths*)
maxDeaths=1000;(*deaths at which social distancing is triggered*)
tsdfD=T; (*social distancing triggered by deaths ends*)
B2=0.85;  (*intesity of social distancing by time*)
tsd=70; (*social distancing starts*)
tsdf=T; (*social distancing ends*)
Exposed0=60;(*exposed at t=0*)
Ro0=5.;(*approximate basic reproduction number (Subscript[R, 0])*)
CFR=0.034;(*case fatality ratio*)
EtoI=5;(*incubation time*)
ItoR=21;(*time from advanced stage to recovery*)
RtoS=(10^10);(*immunization time*)
μ=0*9.40*(10^(-3))*(365^(-1)); (*normal death rate*)
η=0*9.40*(10^(-3))*(365^(-1)); (*normal birth rate*)
Dsd=maxDeaths/pop;
R0=0;
E0=Exposed0/pop;
S0=1-E0;
I0=0;
γ=1/ItoR;
λ=Ro0*γ;
a=1/EtoI;
ξ=1/RtoS;
ν=-((CFR γ)/(-1+CFR));

s1=NDSolve[{S'[t]==-λ (1-A Sin[(π/365)t]^2)(1-B1 (UnitStep[(1-EE[t]-S[t]-II[t]-R[t])-Dsd]-UnitStep[t-tsdfD]))S[t]II[t]+η(S[t]+EE[t]+II[t]+R[t])-μ S[t]+ξ R[t],EE'[t]==λ (1-A Sin[(π/365)t]^2)(1-B1 (UnitStep[(1-EE[t]-S[t]-II[t]-R[t])-Dsd]-UnitStep[t-tsdfD]))S[t]II[t]-μ EE[t]-a EE[t],II'[t]==a EE[t]+(-γ-(ν+μ))II[t],R'[t]==-ξ R[t]+γ II[t]-μ R[t],EE[0]==E0,R[0]==R0,S[0]==S0,II[0]==I0},{S[t],II[t],R[t],EE[t]},{t,0,T}];
s2=NDSolve[{S'[t]==-λ (1-A Sin[(π/365)t]^2)(1-B2(UnitStep[t-tsd]-UnitStep[t-tsdf]))S[t]II[t]+η(S[t]+EE[t]+II[t]+R[t])-μ S[t]+ξ R[t],EE'[t]==λ (1-A Sin[(π/365)t]^2)(1-B2(UnitStep[t-tsd]-UnitStep[t-tsdf]))S[t]II[t]-μ EE[t]-a EE[t],II'[t]==a EE[t]+(-γ-(ν+μ))II[t],R'[t]==-ξ R[t]+γ II[t]-μ R[t],EE[0]==E0,R[0]==R0,S[0]==S0,II[0]==I0},{S[t],II[t],R[t],EE[t]},{t,0,T}];Plot[Evaluate[pop*{II[t]/.s1,II[t]/.s2}],{t,0,T},PlotRange->All,PlotLegends->{"Infected (" <> ToString[Evaluate[B1*100]]<>"% drop in contact rate after "<> ToString[N[Dsd*pop]]<>" deaths)","Infected ("<> ToString[Evaluate[B2*100]]<>"% drop in contact rate from day "<> ToString[N[tsd]]<>" to day "<>ToString[N[tsdf]] <>")"},PlotLabel->"SEIRS MODEL\n Contact rate: "<>ToString[λ]<>" per day\n Average incubation period: "<>ToString[ Evaluate[(1/a)]]<> " days\n Average infectious period: "<>ToString[ Evaluate[(1/γ)]]<> " days\n Average immunization period: "<>ToString[ Evaluate[(1/ξ)]]<> " days\n Case fatality ratio: "<>ToString[N[100 CFR]]<>"%\n Maximum summer drop in contact rate: "<> ToString[N[100 A]]<> "%",ImageSize->800]
Plot[Evaluate[pop*{(1-EE[t]-S[t]-II[t]-R[t])/.s1,(1-EE[t]-S[t]-II[t]-R[t])/.s2}],{t,0,T},PlotRange->Automatic,PlotLegends->{"Dead (" <> ToString[Evaluate[B1*100]]<>"% drop in contact rate after "<> ToString[N[Dsd*pop]]<>" deaths)","Dead (" <> ToString[Evaluate[B2*100]]<>"% drop in contact rate from day "<> ToString[N[tsd]]<>" to day "<>ToString[N[tsdf]] <>")"},PlotLabel->"SEIRS MODEL\n Contact rate: "<>ToString[λ]<>" per day\n Average incubation period: "<>ToString[ Evaluate[(1/a)]]<> " days\n Average infectious period: "<>ToString[ Evaluate[(1/γ)]]<> " days\n Average immunization period: "<>ToString[ Evaluate[(1/ξ)]]<> " days\n Case fatality ratio: "<>ToString[N[100*CFR]]<>"%\n Maximum summer drop in contact rate: "<> ToString[N[100 A]]<> "%",ImageSize->800]