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.legend(loc=5, borderaxespad=0.,fontsize=20)
    
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.legend(loc=2, borderaxespad=0.,fontsize=20)
    
plt.xlabel('$t$',fontsize=30); #plt.ylabel('$x(t)$',fontsize=20)
plt.xlim(0,10); plt.ylim(0,1)
plt.text(8,0.6,'$R_0<1$',fontsize=30)
plt.tick_params(axis='both', which='major', labelsize=15)

plt.show;