CoCalc Public FilesLabs / Lab6 / Lab6.ipynb
Author: Kimberly Vargas Negrete
Views : 263
Compute Environment: Ubuntu 20.04 (Default)

## Lab 6 Math 119A

#### SIR model for an infectious disease

$S$ is the number (in millions) of suceptible people, $I$ is the number (in millions) of infectious people, and $R$ is the number (in millions) of recovered and immune people. If we ignore births and deaths, the equations that describe the changes in these populations are: $S'=-a SI/P_0\quad I'=a SI/P_0 - b I\qquad R'=b I$ We will assume that our total population is $P_0=330$ million people. We start with 100 infected people, $I(0)=\frac{100}{1,000,000}=0.0001$, $S(0)=P_0-I(0)=329.9999$ million people and $R(0)=0$. (No one has recovered yet.)

Note that if we do not choose to track $R$ the number of recovered people, that the change equations for $S$ and $I$ are defined only in terms of each other.

## Part 1: Trajectories

Throughout this lab, examples will be given from the Romeo and Juliet model. We start by plotting a vector field:

In [6]:
var("J,R")
Jprime(J,R)=R               # formula for J'
Rprime(J,R)=-J              # formula for R'
pvf=plot_vector_field([Jprime,Rprime],(J,-2,2),(R,-2,2),color="green",axes_labels=["J","R"],frame=False,aspect_ratio=1)  # makes a vector field plot
show(pvf)                   # displays vector field plot


Exercise 1. For the SIR model, with $a=1/6$ and $b=1/20$, plot a vector field in the $SI$-plane with $S$ values going from $0$ to $350$ and $I$ values going from $0$ to $200$.

In [5]:
var("S,I")
a=1/6
b=1/20
P0=330
Sprime(S,I)=(-a*S*I)/P0
Iprime(S,I)=(a*S*I)/P0-b*I
vec1=plot_vector_field([Sprime,Iprime], (S,0,350), (I,0,200), color="green", axes_labels=["S","I"], frame=False, aspect_ratio=1)
show(vec1)


Exercise 2. Using the CoCalc functions point() and text() which were introduced in Lab 5, add the following labeled points to your vector field:

1. A point near the beginning of the disease where the number of susceptible people is high and the number of infectious people is low. Label this point "Outbreak Start".

2. Follow the arrows of the vector field and estimate the highest point of the trajectory starting from your first point. Add this high point to the vector field and label it "Max I". This represents the maximum number of people infected by the disease at a single time.

3. Following the trajectory you are imagining down to the $S$-axis and add a point near where you think it would intersect the $S$-axis. Label this point "Outbreak Ends". This point represents the end of the disease, when there are no more infected people.

In [1]:
var("S,I") a=1/6 b=1/20 P0=330 Sprime(S,I)=(-a*S*I)/P0 Iprime(S,I)=(a*S*I)/P0-b*I vec1=plot_vector_field([Sprime,Iprime], (S,0,350), (I,0,200), color="green", axes_labels=["S","I"], frame=False, aspect_ratio=1)+point([300,25], color="blue", size=40)+text("Outbreak Start", [300,25-10])+point([110,130], color="blue", size=40)+text("Max I", [110, 130-10])+point([28,0], color="blue", size=40)+text("Outbreak End", [32, 10]) show(vec1)


From the above exercise you should have gotten an idea of how the number of susceptible and number of infected people change over the course of the disease (i.e. you have an idea of what the trajectory looks like). In the next few exercises we will plot a trajectory to see if it matches your intuition.

Exercise 3. Add the point $I(0)=0.0001$, $S(0)=P_0-I(0)=329.9999$ to your vector field and label it "Start". These correspond to 100 initial infected people and just under 330 million susceptible people.

In [2]:
var("S,I") a=1/6 b=1/20 P0=330 Sprime(S,I)=(-a*S*I)/P0 Iprime(S,I)=(a*S*I)/P0-b*I vec1=plot_vector_field([Sprime,Iprime], (S,0,350), (I,0,200), color="green", axes_labels=["S","I"], frame=False, aspect_ratio=1)+point([300,25], color="blue", size=40)+text("Outbreak Start", [300,25+10])+point([110,130], color="blue", size=40)+text("Max I", [110, 130-10])+point([28,0], color="blue", size=40)+text("Outbreak End", [32, 10])+point([329.999,.0001], color="blue", size=40)+text("Start", [329.999, .0001+10]) show(vec1)

In [ ]:



We are now going to create a trajectory starting from the point in the previous exercise. To do this we will use our Euler2D function from Lab 3:

In [3]:
def Euler2D(Xprime,Yprime,X0,Y0,t_end,dt):
Xlist=[X0]          # starts list for X values
Ylist=[Y0]          # starts list for Tvalues
X=X0                # initial value for X
Y=Y0                # initial value for Y
for t in srange(dt,t_end,dt,include_endpoint=true):  # loops thru time values after tstart
Xp=Xprime(X,Y)  # gives derivative X' at (X,Y)
Yp=Yprime(X,Y)  # gives derivative Y' at (X,Y)
dX=dt*Xp        # computes delta X
dY=dt*Yp        # computes delta Y
X=X+dX          # finds new X value
Y=Y+dY          # finds new Y value
Xlist.append(X) # apends new X value to Xlist
Ylist.append(Y) # apends new Y value to Ylist
XYlist=[Xlist,Ylist]# Combines Xlist and Ylist to make a list of lists
return XYlist #returns a list which is made up of the list of X values and the list of Y values


Euler2D allows us to solve systems of differential equations (such as our Romeo and Juliet model or our SIR model). Here is an example of how we can plot a trajectory using the output of Euler2D:

In [4]:
var("J,R")                  # make J and R symbolic variables
Jprime(J,R)=R               # formula for J'
Rprime(J,R)=-J              # formula for R'
tend=7.
dt=0.001
JRlist=Euler2D(Jprime,Rprime,1.,0.,tend,dt)
Jvals=JRlist[0]             # extracts the J values from JRlist
Rvals=JRlist[1]             # extracts the R values from JRlist
# make a vector field plot
pvf=plot_vector_field([Jprime,Rprime],(J,-2.,2.),(R,-2.,2.),color="green",axes_labels=["J","R"],frame=False,aspect_ratio=1)
## make a trajectory plot
ptr=list_plot(list(zip(Jvals,Rvals)),plotjoined=True,axes_labels=["J","R"],aspect_ratio=1)
## combine vector field and trajectory and show
show(pvf+ptr)


Exercise 4. Run the code for Euler2D above, then use it to plot a trajectory (for our SIR model) starting at the point from the previous exercise. Your time values should start at 0 days and go to 365 days in 1 day increments.

In [7]:
var("S,I") a=1/6 b=1/20 P0=330 I0=.0001 S0=330-I0 Sprime(S,I)=(-a*S*I)/P0 Iprime(S,I)=(a*S*I)/P0-b*I tend=365. dt=1. SIlist=Euler2D(Sprime,Iprime,330,.0001,tend,dt) Svals=SIlist[0] Ivals=SIlist[1] pvf1=plot_vector_field([Sprime,Iprime],(S,0,350),(I,0,200), color="green", axes_labels=["S","I"], frame=False, aspect_ratio=1) ptr1=list_plot(list(zip(Svals,Ivals)),  plotjoined=True, axes_labels=["S","I"], aspect_ratio=1) show(pvf1+ptr1)


We can bring trajectories to life using animations. Here is an example with Romeo and Juliet:

In [8]:
var("J,R")
Jprime(J,R)=R               # formula for J'
Rprime(J,R)=-J              # formula for R'
dt=0.001                    # time step
tend=2.*pi
JRlist=Euler2D(Jprime,Rprime,1.,0.,tend,dt)
tvals=srange(0.,tend,dt,include_endpoint=True)  # list of time values
Jvals=JRlist[0]             # extracts the J values from JRlist
Rvals=JRlist[1]             # extracts the R values from JRlist
pvf=plot_vector_field([Jprime(J,R),Rprime(J,R)],(J,-2.,2.),(R,-2.,2.),color="green",axes_labels=["J","R"],frame=False,aspect_ratio=1)# makes a vector field plot
## trajectory
ptr=list_plot(list(zip(Jvals,Rvals)),plotjoined=True,axes_labels=["J","R"],aspect_ratio=1)
list_frames=[]                     # creates a list for frames to go into
for j in srange(0,len(tvals),200): # loops going from 0 to  the number of time values in tvals skipping by 200 each time.
pt=point([Jvals[j],Rvals[j]], color="red", size=40)  # plots the point at the j-th time value
p=pvf+pt+ptr                   # combine into one plot
list_frames.append(p)          # appends this plot to list_frames
ani=animate(list_frames)           # creates an animation from the list of plots
show(ani,iterations=3)             # displays the animation only repeat 3 times