SIRS model analysis for Math 5525

by Nasser M. Abbasi. April 10, 2020. Using Mathematica 12.1

SetDirectory[NotebookDirectory[]]

Simply evaluate this notebook to run the program. Or just evaluate the CELL below. The whole program is in one cell below.

Manipulate[
tick;

Module[{ic,eq1,eq2,eq3,S,R,I0,t,initialR,sol},
If[initialS+initialI>1000,initialI=1000-initialS];
initialR=1000-(initialS+initialI);
ic={S[0]==initialS,I0[0]==initialI,R[0]==initialR};
eq1=S'[t]==-beta S[t] *I0[t]+mu*R[t];
eq2=I0'[t]==beta*S[t]*I0[t]-v I0[t];
eq3=R'[t]==v*I0[t]-mu*R[t];
sol=NDSolve[{eq1,eq2,eq3,ic},{S,I0,R},{t,0,60}];
Grid[{
{Grid[{
{"Transmission rate (beta)  ",padIt2[beta,{5,4}]},
{"Recovery rate (v) ",padIt2[v,{3,3}]},
{"reinfection rate (mu) ",padIt2[mu,{3,3}]},
{"v/beta ",padIt2[v/beta,{6,3}]},
{},
{If[1000≥v/beta,
Row[{Style["case (a)",Bold]," population >= v/beta"}],
Row[{Style["case (b)",Bold]," population <= v/beta"}]
]}},Frame->True
]
},
{Plot[{
Evaluate[S[t]/.sol],
Evaluate[I0[t]/.sol],
Evaluate[R[t]/.sol],
If[addSI,
Evaluate[S[t]/.sol]+Evaluate[I0[t]/.sol],Nothing],
If[addSI,1000]

},
{t,0,60},
PlotLegends->If[addSI,{"susceptible (S)","Infected (I)","Recovered (R)","S+I","Total population"},{"susceptible (S)","Infected (I)","Recovered (R)"}],
PlotStyle->If[addSI,{Blue,Red,Black,Dashed,Green},{Blue,Red,Black}],
GridLines->Automatic,GridLinesStyle->LightGray,
AxesLabel->{"Time (days)","Population size"},ImageSize->500,
ImagePadding->30,
PlotRange->{{0,60},{0,1000}}
]
}}]
],
{{beta,0.0005,"Transmission rate (beta)"},0.0001,0.003,0.0001,Appearance->"Labeled"},
{{v,0.05,"Recovery rate (v)"},0.001,0.4,0.001,Appearance->"Labeled"},
{{mu,0.05,"Recovery to Susceptible rate (mu)"},0.001,0.4,0.001,Appearance->"Labeled"},
{{initialI,5,"Initial infections size I(0)"},0,1000,1,Appearance->"Labeled"},
{{initialS,995,"Initial susceptible size S(0)"},0,1000,1,Appearance->"Labeled"},
Grid[{{
Button[Text[Style["First critical point",11]],{initialI=0,initialS=1000,tick=Not[tick]},ImageSize->{100,35}],
Button[Text[Style["second critical point",11]],
tick=Not[tick];
If[1000≥ v/beta,
{initialI=mu/(v+mu)*(1000-v/beta),initialS=v/beta}
,Nothing
],ImageSize->{130,35}],
Button[Text[Style["Add S+I curve",11]],{addSI=True,tick=Not[tick]},ImageSize->{100,35}],
Button[Text[Style["Remove S+I curve",11]],{addSI=False,tick=Not[tick]},ImageSize->{110,35}]

}
}
],
{{tick,False},None},
{{addSI,True},None},
TrackedSymbols:>{beta,v,mu,initialI,initialS,tick},
SaveDefinitions->False,
Initialization:>(padIt2[v_,f_List]:=AccountingForm[v,f,NumberSigns->{"",""},
NumberPadding->{"0","0"},SignPadding->True])
]