# SIR model¶

#### Please see section 2.6.3 in the textbook for description of all the parameters and variables used below.¶

Let us work with nondimensionalized form of the SIR model:

Defining $x=\frac{S}{N}$, $y=\frac{I}{N}$ and $z=\frac{R}{N}$

The charcateristic time scale in the problem is $t_c=\frac{1}{\gamma}$. Defining new dimesnionaless quantity $\tau=\gamma t$

The dimensionless equations are:

\begin{eqnarray} x' & =& -R_0xy \\ y' & =& R_0xy-y \\ z' & =& y \end{eqnarray}

Where $R_0=\displaystyle\frac{\alpha N}{\gamma}$ (the basic reporduction number) and the primes in the equations stand for derivative w.r.t $\tau$. The equation for $z$ is redudant as $z$ can be found by $z(t)=1-x(t)-y(t)$. We know if $x(0)=x_0$, $y(0)=y_0$ and $z(0)=0$ then $x_0+y_0=1$. Also, note that $x+y \leq 1$.

An explicit solution for the $x$ and $y$ is also possible, if we divide $y'$ by $x'$ i.e., evaluating $y'/x'$ by using the above initial conditions, we obtain $y=1-x+\frac{1}{R_0}\ln \frac{x}{x_0}$.

Below we plot the phase plot and the time series for the two cases $R_0>1$ and $R_0<1$.

In [1]:
# import all the packages required
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

In [2]:
# Solve the differential equations and plot the phase portrait
from scipy.integrate import odeint

# the value of R0=3.5
def fun(X,t):
R0=3.5
x=X[0];y=X[1];
dx=-R0*x*y; dy=R0*x*y-y
return [dx,dy]

plt.figure(figsize=(10,6))

# Times at which the solution is to be computed.
t = np.linspace(0, 120, 2000)

# Initial x and y
xo1=np.linspace(0, 1, 150)
yo1=1-xo1

# plot x+y<=1 line
plt.plot(xo1,yo1,'k--',lw=2)

#let's plot some orbits/trajectories

Xinit=np.array([xo1[148],yo1[148]])
Xsolu = odeint(fun,Xinit, t)
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=4)

Xinit=np.array([xo1[110],yo1[110]])
Xsolu = odeint(fun,Xinit, t)
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=4)

Xinit=np.array([xo1[70],yo1[70]])
Xsolu = odeint(fun,Xinit, t)
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=4)

Xinit=np.array([xo1[50],yo1[50]])
Xsolu = odeint(fun,Xinit, t)
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=4)

Xinit=np.array([xo1[10],yo1[10]])
Xsolu = odeint(fun,Xinit, t)
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=4)

# let's plot the critical point as well
#plt.plot(0.0,0.0,'or',ms=12)

# stuff required for plotting arrows/vectors on the plot
xmax = plt.xlim(xmin=0)[1];ymax = plt.ylim(ymin=0)[1]

xl = np.linspace(0.0,xmax,20)
yl = np.linspace(0.0,ymax,20)
x1 , y1  = np.meshgrid(xl, yl)
[u,v]= fun([x1,y1],0);

# normalize the lengths of the vectors
uvn=np.sqrt(u**2 + v**2); uvn[uvn==0]=1.
u=u/uvn; v=v/uvn;

plt.quiver(x1[x1+y1<=1],y1[x1+y1<=1],u[x1+y1<=1],v[x1+y1<=1],uvn[x1+y1<=1],scale=40,cmap=plt.cm.gray)
plt.xlim(0,1); plt.ylim(0,1)
plt.xlabel('$x(t)$',fontsize=30); plt.ylabel('$y(t)$',fontsize=30)
plt.text(0.8,0.8,'$R_0>1$',fontsize=30)
plt.tick_params(axis='both', which='major', labelsize=15)

plt.show;

In [3]:
# Solve the differential equations and plot the time series

plt.figure(figsize=(10,8))

t = np.linspace(0, 10, 2000)
Xinit=np.array([xo1[140],yo1[140]])
Xsolu = odeint(fun,Xinit, t)
plt.plot(t,Xsolu[:,0],'b',lw=2,label='Susceptible')
plt.plot(t,Xsolu[:,1],'r',lw=2,label='Infected')
plt.plot(t,1-Xsolu[:,1]-Xsolu[:,0],'g',lw=2,label='Recovered')

plt.xlabel('$t$',fontsize=30); #plt.ylabel('$x(t)$',fontsize=20)
plt.text(8,0.8,'$R_0>1$',fontsize=30)
plt.xlim(0,10); plt.ylim(0,1)
plt.tick_params(axis='both', which='major', labelsize=15)

plt.show;

In [4]:
# Solve the differential equations and plot the phase portrait
from scipy.integrate import odeint

# the value of R0=3.5
def fun(X,t):
R0=0.75
x=X[0];y=X[1];
dx=-R0*x*y; dy=R0*x*y-y
return [dx,dy]

plt.figure(figsize=(10,6))

# Times at which the solution is to be computed.
t = np.linspace(0, 120, 2000)

# Initial x and y
xo1=np.linspace(0, 1, 150)
yo1=1-xo1

# plot x+y<=1 line
plt.plot(xo1,yo1,'k--',lw=2)

#let's plot some orbits/trajectories
Xinit=np.array([xo1[148],yo1[148]])
Xsolu = odeint(fun,Xinit, t)
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=4)

Xinit=np.array([xo1[110],yo1[110]])
Xsolu = odeint(fun,Xinit, t)
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=4)

Xinit=np.array([xo1[70],yo1[70]])
Xsolu = odeint(fun,Xinit, t)
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=4)

Xinit=np.array([xo1[50],yo1[50]])
Xsolu = odeint(fun,Xinit, t)
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=4)

Xinit=np.array([xo1[10],yo1[10]])
Xsolu = odeint(fun,Xinit, t)
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=4)

# let's plot the critical point as well
#plt.plot(0.0,0.0,'or',ms=12)

# stuff required for plotting arrows/vectors on the plot
xmax = plt.xlim(xmin=0)[1];ymax = plt.ylim(ymin=0)[1]

xl = np.linspace(0.0,xmax,20)
yl = np.linspace(0.0,ymax,20)
x1 , y1  = np.meshgrid(xl, yl)
[u,v]= fun([x1,y1],0);

# normalize the lengths of the vectors
uvn=np.sqrt(u**2 + v**2); uvn[uvn==0]=1.
u=u/uvn; v=v/uvn;

plt.quiver(x1[x1+y1<=1],y1[x1+y1<=1],u[x1+y1<=1],v[x1+y1<=1],uvn[x1+y1<=1],scale=40,cmap=plt.cm.gray)
plt.xlim(0,1); plt.ylim(0,1)
plt.xlabel('$x(t)$',fontsize=30); plt.ylabel('$y(t)$',fontsize=30)
plt.text(0.8,0.8,'$R_0<1$',fontsize=30)
plt.tick_params(axis='both', which='major', labelsize=15)

plt.show;

In [5]:
plt.figure(figsize=(12,8))

t = np.linspace(0, 10, 2000)
Xinit=np.array([xo1[60],yo1[60]])
Xsolu = odeint(fun,Xinit, t)
plt.plot(t,Xsolu[:,0],'b',lw=2,label='Susceptible')
plt.plot(t,Xsolu[:,1],'r',lw=2,label='Infected')
plt.plot(t,1-Xsolu[:,1]-Xsolu[:,0],'g',lw=2,label='Recovered')

plt.xlabel('$t$',fontsize=30); #plt.ylabel('$x(t)$',fontsize=20)
plt.text(8,0.6,'$R_0<1$',fontsize=30)