Kermack McKendrick SIR Model

<<Graphics`MultipleListPlot`

Preparations

fn1 = -beta Sx[t] Ix[t];
fn2 = beta Sx[t] Ix[t] - gamma Ix[t];
fn3 = gamma Ix[t];

difeq = {Sx'[t]==fn1,Ix'[t]==fn2,Rx'[t]==fn3,
	Sx[0]==N0-I0,Ix[0]==I0,Rx[0]==0}

N0Init = 200

Dependency of the initial condition

para = {beta->0.001, gamma->0.1, N0 -> N0Init};

intv = 500;
traj1 = NDSolve[difeq/.para/.I0->1,{Sx[t],Ix[t],Rx[t]},{t,0,intv}];
Plot[ Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj1],{t,0,intv}]
phase1 = ParametricPlot[ Evaluate[{Sx[t],Ix[t]}/.traj1],{t,0,intv}]

Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj1]/.t->intv

para = {beta->0.001, gamma->0.1, N0 -> N0Init};

intv = 500;
traj2 = NDSolve[difeq/.para/.I0->50,{Sx[t],Ix[t],Rx[t]},{t,0,intv}];
Plot[ Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj2],{t,0,intv}]
phase2 = ParametricPlot[ Evaluate[{Sx[t],Ix[t]}/.traj2],{t,0,intv}]

Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj2]/.t->intv

para = {beta->0.001, gamma->0.1, N0 -> N0Init};

intv = 500;
traj3 = NDSolve[difeq/.para/.I0->100,{Sx[t],Ix[t],Rx[t]},{t,0,intv}];
Plot[ Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj3],{t,0,intv}]
phase3 = ParametricPlot[ Evaluate[{Sx[t],Ix[t]}/.traj3],{t,0,intv}]

Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj3]/.t->intv

para = {beta->0.001, gamma->0.1, N0 -> N0Init};

intv = 500;
traj4 = NDSolve[difeq/.para/.I0->150,{Sx[t],Ix[t],Rx[t]},{t,0,intv}];
Plot[ Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj4],{t,0,intv}]
phase4 = ParametricPlot[ Evaluate[{Sx[t],Ix[t]}/.traj4],{t,0,intv}]

Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj3]/.t->intv

g1 = Plot[ N0Init - t,{t,0,N0Init}]

Show[g1, phase1,phase2,phase3, phase4, PlotRange->{{0,N0Init},{0,N0Init}},
	AspectRatio->1]

2

What proportion of people remain uninfected?

para = {beta->0.001, gamma->0.1, N0 -> 200};

intv = 400;
traj5 = NDSolve[difeq/.para/.I0->1,{Sx[t],Ix[t],Rx[t]},{t,0,intv}];
Plot[ Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj5],{t,0,intv}]
phase5 = ParametricPlot[ Evaluate[{Sx[t],Ix[t]}/.traj5],{t,0,intv}]

Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj5]/.t->intv

para = {beta->0.001, gamma->0.1, N0 -> 300};

intv = 400;
traj6 = NDSolve[difeq/.para/.I0->1,{Sx[t],Ix[t],Rx[t]},{t,0,intv}];
Plot[ Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj6],{t,0,intv}]
phase6 = ParametricPlot[ Evaluate[{Sx[t],Ix[t]}/.traj6],{t,0,intv}]

Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj6]/.t->intv

para = {beta->0.001, gamma->0.1, N0 -> 400};

intv = 400;
traj7 = NDSolve[difeq/.para/.I0->1,{Sx[t],Ix[t],Rx[t]},{t,0,intv}];
Plot[ Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj7],{t,0,intv}]
phase7 = ParametricPlot[ Evaluate[{Sx[t],Ix[t]}/.traj7],{t,0,intv}]

Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj7]/.t->intv

para = {beta->0.001, gamma->0.1, N0 -> 150};

intv = 400;
traj8 = NDSolve[difeq/.para/.I0->1,{Sx[t],Ix[t],Rx[t]},{t,0,intv}];
Plot[ Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj8],{t,0,intv}]
phase8 = ParametricPlot[ Evaluate[{Sx[t],Ix[t]}/.traj8],{t,0,intv}]

Evaluate[{Sx[t],Ix[t],Rx[t]}/.traj7]/.t->intv

Show[phase5,phase6,phase7,phase8, PlotRange->All,AspectRatio->Automatic]

5

Fraction of those escaping from infection

eq  =  rate  ==  1  -  Exp[-R rate]

fn[x_] := 
      Block[{tmp}, 
tmp  =  
          eq /. R -> x; 
sol  =  
          FindRoot[tmp, {rate, 0.4}]; 
Return[
          rate /. sol];
]

Plot[fn[x], {x, 0, 4}, PlotRange -> {All, {0, 1}}]