Kermack McKendrick SIR Model

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}}]

SIR model with demography

deq = {Sx'[t] == mu K - beta Sx[t] Ix[t] - mu Sx[t],
	Ix'[t] == beta Sx[t]Ix[t] - mu Ix[t] - alpha Ix[t] - gamma Ix[t],
Rx'[t] == gamma Ix[t] - mu Rx[t],
Sx[0]==S0, Ix[0]==I0, Rx[0]==R0}

Slight disease

para = {mu->0.1,beta->0.001, gamma->0.1,alpha->0.05, K->1000}

deq/.para

deqN = deq/.para/.{S0->1000,I0->0.01, R0->0}

sol = NDSolve[deqN,{Sx[t],Ix[t],Rx[t]},{t,0,100}]

Plot[Evaluate[{Sx[t],Ix[t],Rx[t]}/.sol],{t,0,100},PlotRange->All,PlotStyle->{RGBColor[0,0,1],RGBColor[1,0,0],RGBColor[0,1,0]}]

Plot[Evaluate[Sx[t]+Ix[t]+Rx[t]/.sol],{t,0,80},PlotRange->{0,1200}]

Mildly virulent disease

para = {mu->0.1,beta->0.001, gamma->0.1,alpha->0.3, K->1000}

deq/.para

deqN = deq/.para/.{S0->1000,I0->0.01, R0->0}

sol = NDSolve[deqN,{Sx[t],Ix[t],Rx[t]},{t,0,100}]

Plot[Evaluate[{Sx[t],Ix[t],Rx[t]}/.sol],{t,0,100},PlotStyle->{RGBColor[0,0,1],RGBColor[1,0,0],RGBColor[0,1,0]}]

Plot[Evaluate[Sx[t]+Ix[t]+Rx[t]/.sol],{t,0,100},PlotRange->{0,1200}]

Severely virulent disease

para = {mu->0.1,beta->0.001, gamma->0.1,alpha->0.7, K->1000}

deq/.para

deqN = deq/.para/.{S0->1000,I0->0.01, R0->0}

sol = NDSolve[deqN,{Sx[t],Ix[t],Rx[t]},{t,0,300}]

Plot[Evaluate[{Sx[t],Ix[t],Rx[t]}/.sol],{t,0,300},PlotRange->All,PlotStyle->{RGBColor[0,0,1],RGBColor[1,0,0],RGBColor[0,1,0]}]

Plot[ Evaluate[Ix[t]/.sol],{t,0,300},PlotStyle->RGBColor[1,0,0],PlotRange->All]

Plot[Evaluate[Sx[t]+Ix[t]+Rx[t]/.sol],{t,0,300},PlotRange->{0,1200}]

eqmI = mu/beta(K-eqmS)/eqmS/.eqmS->(mu + alpha + gamma)/beta

para2 = {mu->0.1,beta->0.001, gamma->0.1, K->1000}

eqmI/.para2

lossByDisease = alpha eqmI/mu

Plot[lossByDisease/.para2,{alpha,0,0.8}]