Trigger warning¶
When Dr. Capehart does this in a live session, you get to hear the full (and sometimes uncensored) rant about when a certain PhD student brought their child to sit (in his good chair) with a very distinctive cough to find out by an email the day after that he had been exposed to Whopping Cough.
Get your Tdap!
Introduction¶
This a Python program is an example of a "Kermack-McKendrick Model" in epidemic modeling.
This is also a good example of
- Scripting in Python
- Solving Differential Equations (in this case Brute-Force Finite Differencing) in Python
- Creating a Policy-Relevant Model
- Illustrating how to code coupled systems of equations.
- Why you should be getting your flu shot.
The basics of a Kermack-McKendrick Model is that you have [typically] three sub-populations in a community where something icky is going around. This is normally called an "SIR" Model:
- Susceptibles : people who in principle can get sick (in this case we will bundle a block of people who are immunized into this block)
- Infected : people who have the disease and can spread it into the Susceptable population
- Recovered : people who have had the disease and are now immune.
- To this category we're adding a "Fatalities" category because we are monsters. "We also don't use Salk"
Other versions of this system can distribute the Infected category into subgroups that symptomatic (changing their behavior in contrast to the other classes [not including the dead], contagious [at different levels of symptomatic or asymptomatic]).
You can break the "Susceptible" group into those with temporary immunity from their mothers, permanent immunity by random chance, and immunity from immunization.
You can have recovered people who have been infected lose their immunity.
In short, you can let your Pestilence Horseman run wild because “it’s only a model”
Remember that this is a simple model. The pro's working in operational epidemiology have far more sophisticated models.
Sketching the Model Environment and the Kermack-McKendrick Equations¶
It's always nice to sketch and doodle our models. The black letters will be our variables. (We're not doing the "cordon sanitaire" option (where healthy people wall themselves off from the sick people) in this example.
This creates the following system of equations with the labeled parameters above
$$ \begin{aligned} \frac{dS}{dt} &= -\beta SI \\ \\ \frac{dI}{dt} &= +\beta SI - \gamma I - \omega I \\ \\ \frac{dR}{dt} &= +\gamma I \\ \\ \frac{dF}{dt} &= +\omega I \\ \\ \end{aligned} $$The above $\beta$ term is a combination of $\kappa$, $p$, and $i$. We will show how they are constructed as we move through this exercise.
################################################
#
# Libraries
#
import numpy as np
import scipy as sci
import scipy.integrate as sciint
import sympy as sym
import matplotlib.pyplot as plt
# Notice here that this is using a very different
# method to pull in functions.
#
# I normally don't import libraries this way
#
# This is simply (as the name implies) to bring
# in a youtube video.
from IPython.display import YouTubeVideo
#
################################################
Setting up the Code¶
To do this we will break this process down into several sub-sections
- set up the time extent of our simulation, "Time Control"
- set up the initial state of our players, "Initial or Boundary Conditions"
- set up the rules for the simulation, "Initialize Parameters"
First let’s set up the time control of our simulation. There are variants of this system where you can do truly horrific things to the players to this system such as lose their immunity after a while, add new babies who can get sick easily or have temporary maternal immunity, recovered populations become sterile, or weaker with a higher natural death rate. AES/CEE 615 (Environmental Systems Modeling) explores some of these scenarios at different levels depending on how ~disturbed~ creative the students are for that given semester.
But in this case we are playing a short-game lasting about a month. We’d like some high temporal resolution. So let’s start a simulation of 30 days with a time granularity (time step, $\Delta t$) of about 0.125 days (3-hrs).
################################################
#
# Time Control
#
#
# Total Duration and Time Interval
#
TimePeriod = 30 # days
DeltaT = 3./24. # hours -> days
#
# Number of timesteps from 0 to TimePeriod
#
nTime = int(TimePeriod / DeltaT + 1)
#
# Time Array from 0 to 30
#
# This is different from some array generators
# in python since linspace acutally gives you
# the stopping value you *asked* for! (No
# need to add an extra time step.)
#
Time = np.linspace(start = 0,
stop = TimePeriod, # inclusive must be an integer
num = nTime)
#
################################################
Setting up Initial/Boundary Conditions¶
The above equations need a starting point. In this case, it's our cast of characters.
- Starting Susceptables, $S_0$ (in this run this includes people who are immunized, so technically this just all healthy people).
- Initial Infected Cases, $I_0$ (This/These are our "Patient(s) Zero. Pro-Tip: Don't be these guys.) Immunization Rant will be down further in this notebook!
- Recovered, $R_0$ (not to be confused with another $R_o$ that we will discussed shortly)
- Fatalites, $F_0$ (because as mentioned above, we are monsters.)
################################################
#
# Initial Conditions
#
S_0 = 998.0 # Uninfected People (includes Vaccinated)
I_0 = 1.0 # Patients Zero: Arg! We *HATE* that guy!
R_0 = 0.0 # Recovered: 0 since we're just starting
F_0 = 0.0 # Fatalities: Again, 0 since we're just starting
#
################################################
Creating our Arrays.¶
Here we are going borrow from a certain British TV cooking presenter and repeat our mantra for this class, "maximum satisfaction with minimal effort."
We already have an array as long as we want, our Time[] array. All we have to do is just borrow it, multuply it by zero and we have our arrays for our 4 population groups.
################################################
#
# Setting up our arrays
#
S = Time*0 # Uninfected People (multiplying time by zero)
I = Time*0 # Infected (multiplying time by zero)
R = Time*0 # Recovered (multiplying time by zero)
F = Time*0 # Fatalities (multiplying time by zero)
#
# And we can initialize our first (or zero'th since it's python)
# time step from our above initial conditions
#
k = 0
S[k] = S_0 # Uninfected People (includes Vaccinated)
I[k] = I_0 # Infected
R[k] = R_0 # Recovered
F[k] = F_0 # Recovered
print(" Time[", k, "]=", Time[k],
" S[", k, "]=", S[k],
" I[", k, "]=", I[k],
" R[", k, "]=", R[k],
" F[", k, "]=", F[k])
#
################################################
Setting up Basic Disease ~Rules~ Parameters¶
Every model needs a set of quantified rules. For diseases this is difficult to explicitly nail down. As such, we often rely on a set of expected probabilities. Here are the traits we'll be working with
Now onto the disease: We will need to have the
- contact rate: the likely number of people who an infected person comes into contact per unit time ($\kappa$)
- Units: probability of infection per risky contact
- infection rate: the probability of infection when a truly vulerable person is exposed per unit time ($\kappa$)
- Units: person-to-person contacts (of $S$'s, $I$'s and $R$'s) per sick person per unit time
- recovery rate: the probability of recovery per unit time for a given infected person ($\gamma$)
- Units: probability of a given sick person recovering per unit time
- fatality rate: the probability per unit time of an infected person dying from the disesase ($\omega$)
- Units: probability of a given sick person to take the eternal dirt nap it per unit time.
- $\beta$ from the above equation system is a little more complex than the original case and is a combination of several of these terms. We'll address this later.
- Units: new infections per contact between sick and susceptable people per unit time
- immunization fraction: If we want to include factors like artificial or natural immunity ($p$) we can include that, too. Again, we aren't ready for that yet but we'll assume it's zero for now.
- Units: Fraction of people ~smart enough~ not stupid enough to listen to Andrew Wakefield and Jenny McCarthy.
- Education Outreach: my vaccinated child doesn't have autism, he just doesn't like you.
You can get this through some unit analysis. If you follow the above diagram you can see how the parts come together. Beta is officially a time dependant field but it is calculated as a function of values form a single time step. This makes this a "diagnostic" variable. $S\left(t\right)$, $I\left(t\right)$, $R\left(t\right)$, and $F\left(t\right)$ since they are represented in the above equations as time-derivatives.
$$\beta\left(t\right)=\frac{i \kappa (1-p)}{S\left(t\right)+I\left(t\right)+R\left(t\right)}$$################################################
#
# Disease Traits
#
contact_rate = 6.00 # typical number of contacts (exposure) per case per day
infection_rate = 0.25 # probabilty of infection on exposure per case per day
recovery_rate = 0.50 # mean chance of recovery per day per case per day
fatality_rate = 0.10 # mean chance of fatality per case per day
#
# Public Health Traits
#
immunization_fraction = 0.0 # fraction of artificial immunity
#
################################################
Value-added Parameters¶
Basic Reproductioin Number, $R_o$¶
If you've seen the medical thriller "Contagion" the CDC Epidemiologist played by Kate Winslet gives a description of something called the "Basic Reproductive Number," or "R-nought" ($R_o$), factor while plague-splaining the spread of the disease to city planners.
[
# embed youtube video
YouTubeVideo('VrATMF_FB9M')
The $R_o$ describes how many people one infected person wil infect.
Typical values for diseases we may encouter in these scenarios include
| Disease | Means of Contagion | ${R_0}$ |
|---|---|---|
| Measles | Airborne | 12-18 |
| Diphtheria | Saliva | 6-7 |
| Smallpox | Airborne droplet | 5–7 |
| Polio | Fecal-oral route | 5–7 |
| Rubella | Airborne droplet | 5–7 |
| Pertussis | Airborne droplet from coughing in the prof's office | 5-6 |
| Mumps | Airborne droplet | 4–7 |
| (Conjunctivitis) "Pink Eye" | Dirty Lab Goggles | ~4 |
| HIV/AIDS | Sexual contact | 2–5 |
| SARS | Airborne droplet | 2–5 |
| COVID-19 (ests.) | Airborne droplet | 1.4–3.8 |
| Influenza (1918 pandemic strain) | Airborne droplet | 2–3 |
| Ebola (2014 Ebola outbreak) | Bodily fluids | 1.5-2.5 |
| MERS | Airborne droplet | 0.3-0.8 |
$R_o$'s greater than 1 indicate a disease that will spread through the population. The higher the $R_o$, the faster the spread. $R_o$'s below 1 will likely self-contain due to low probability if infection (not every case will result in a new case but an influx of sick people may create a critical mass of "players" for some "Susceptables" to loose to the odds). Using our parameters we can derive $R_o$ via unit analysis as
$$R_o =\frac{i \kappa}{\gamma +\omega}$$################################################
#
# Basic Reproductive Number
#
Ro = infection_rate * contact_rate / (recovery_rate+fatality_rate)
print("Basic Reproductive Number, Ro =", Ro)
#
################################################
#We'reDoingShots!¶
In this case we can play with a fraction of the population that enter the game with either natural immunity or have been immunized. This is the $p$ parameter above.
This also allows us to bring up the idea of herd immunity which is one of the major battlegrounds in the vaccine "controversy."
Let's consider
Assumung that $R=0$ and $S\gg I$,
$$\begin{aligned} \frac{dI}{dt} &= \frac{ i \kappa \left( 1-p \right) }{S+I+ \left( R=0 \right)} S I - \left( \gamma + \omega \right) I \\ 0 &= \frac{i \kappa\left(1-p\right)}{\left[(S+I) \approx S\right]}S\;I-\left(\gamma +\omega \right)I \\ 0 &= \frac{i \kappa\left( 1-p \right) }{S} S I - \left( \gamma +\omega \right) I \\ 0 &= I \left[ \frac{i \kappa \left( 1-p \right) }{S}S-\left( \gamma +\omega \right) \right] \\ 0 &= I \left[ i \kappa \left(1-p\right)- \left( \gamma +\omega \right) \right] \\ 0 &= i \kappa \left( 1-p \right)- \left( \gamma +\omega \right) \\ i \kappa\left(1-p\right) &= \gamma +\omega \\ 1-p &= \frac{\gamma +\omega }{i \kappa} \\ 1-p &= \frac{1}{R_o} \\ p &= 1 - \frac{1}{R_o} \end{aligned} $$This "critical" value of public immunity fraction, $p$, or $p_c$, is the "herd immunity." Thus for a $R_o$ of 12 (the low end of the communicability of Measles), to force the functional $R_o$ for Measles down to zero or less, over 90% of the populus must have artificial immunity.
################################################
#
# Herd Immunity
#
herd_immunity = 1 - 1/Ro
print("Herd Immunity Threshold, pc =",
herd_immunity)
#
################################################
Drafting our Prognostic and Diagnostic Equations¶
In this scenario we are going to show off SciPy's integrating features. Most of them will fit into the following setup.
We will need a single function that contains all of our prognostic equations, $\frac{dS}{dt}$, $\frac{dI}{dt}$, $\frac{dR}{dt}$, and $\frac{dF}{dt}$. Since two of our equations need a value for $\beta$, that gets put in there too.
The Functions for our finite differencing approach¶
Here we will create a single function that bundles all of our equations. In the Differential Equations Hardcore Version of this exercise, we will use a single function with all four equations.
Also we will have an equation for $\beta$ since we will have an opportunity to get fancy with that function as we play.
################################################################################################
#
# Our Four Beta, S, I, R, F function
#
# We pass a vector of [S, I, R, F] through our system.
#
#
# The Beta Function (the rate of aquisition of disease per
# susceptable-infected contacts.)
#
def beta(S,I,R):
beta_out = (1-immunization_fraction)*(infection_rate*contact_rate)/(S+I+R)
return(beta_out)
#
# The Susceptable Equation
#
def dSdt(S,I,R):
dSdt_out = -beta(S,I,R) * I * S
return(dSdt_out)
#
# The Infected Equation
#
def dIdt(S,I,R):
dIdt_out = beta(S,I,R) * I * S - (recovery_rate + fatality_rate) * I
return(dIdt_out)
#
# The Recovered Equation
#
def dRdt(I):
dRdt_out = recovery_rate * I
return(dRdt_out)
#
# The Fatalities Equation
#
def dFdt(I):
dFdt_out = fatality_rate * I
return(dFdt_out)
#
################################################
Time to Bunny-Hop through Time! Euler's Method¶
We are doing a finite differencing approach similar to our earlier adventures where we learned programming with Mathcad.
There are fancier ways to do it but the goal here is more to show you some traditional structures in programming, such as loops and if-then-blocks.
Here is the basic math of any generic prognostic (time-dependant) function.
- Take a differential equations that's a function of time. That $f_x(x)$ function represents all the forcings that are being inflicted on our mystery value x that change its value over time. $$\frac{dx}{dt}=f(x)$$
- Whistle in the dark through your Calc 1 class and hope that $f_x(x)$ doesn't change much over a small discrete change of your independant integrating value (for us that's time), a.k.a., $\Delta t$ $$\frac{\Delta x}{\Delta t}=f(x)$$
- Let's now treat time like a numberline with distinct intervals of space one $\Delta t$ with our understanding of what's going on only resolvable at those discrete points on the hopscotch board that is now our modeled universe. To move forward we base our change on $x$ on our jump to the current time-step for which we're calculating, $t_{\tau}$, on our values that we just solved on our previous timestep, $t_{\tau-1}$. (Hence it's important to chose timesteps where our valus of x don't change MUCH as we move in our steps forward in time. This often means small values of $\Delta t$. Much as in our scripts were we played with Pi using Reimann Sums)
- And from there it's just a little Algebra-Fu™️ to solve for the value of x on our landing point (a.k.a., $x_{\tau}$)
Applying this to our scenario we have four equations moving together through time.
$$\begin{aligned} S_{\tau} &= S_{\tau-1}+f_S(x_{\tau-1})\Delta t \\ I_{\tau} &= I_{\tau-1}+f_I(x_{\tau-1})\Delta t \\ R_{\tau} &= R_{\tau-1}+f_R(x_{\tau-1})\Delta t \\ F_{\tau} &= F_{\tau-1}+f_F(x_{\tau-1})\Delta t \end{aligned}$$Our $f_S$, $f_I$, $f_R$, and $f_F$ are the functions we created above!
Overall this method is called finite differencing. The specific method of marching forward-in-time shown above is call Euler's Method.
To solve this we just need to put things into a loop!
For demonstration purposes, we're going to set this up one step at a time but copying our previous work each step where we add a feature or two to show how we do it specifically in the Python language.
Just remember that above we already have the following "victory laps" already done:
- Created our Time and Time Step Variables
- Created our Modeling Variable Arrays & Initial Conditions
- Created our Model Parameters (e.g., infection rate, recovery rate...)
- Created our Functions for the Kermack-McKendrick Equations.
So onward...
Victory Lap #5. Build the time loop.¶
This is going to be a very simple loop that walks through our time. Here we can loop through our time indices which we'll call $k$ which will go from 0 to nTime.
Remember though that we have to do that wierd thing that Python does where the last falue requested in an array or series is the one AFTER the last one. This is supposedly because this makes it prettier in languages that start with zero (though those of us who have worked with languages that start their indexing at 0 are used to just subtracting one to the number of values you want...).
We'll use numpy.arange() for this
################################################
#
# Step 5: Build the Time Loop
# We are also going to
#
k = 0
print(" Time[", k, "]=", Time[k])
for k in np.arange(start = 1, stop = nTime):
print(" Time[", k, "]=", Time[k])
#
################################################
Victory Lap #6. Insert the equations...¶
Now let's our equations into our loop.
################################################
#
# Step 6: Plumbing Test!
#
k = 0
print(" Time[", k, "]=", Time[k],
" S[", k, "]=", S[k],
" I[", k, "]=", I[k],
" R[", k, "]=", R[k],
" F[", k, "]=", F[k])
for k in np.arange(start = 1, stop = nTime):
S[k] = S[k-1] + dSdt(S[k-1], I[k-1], R[k-1]) * DeltaT
I[k] = I[k-1] + dIdt(S[k-1], I[k-1], R[k-1]) * DeltaT
R[k] = R[k-1] + dRdt(I[k-1]) * DeltaT
F[k] = F[k-1] + dFdt(I[k-1]) * DeltaT
# Because you can have negative personalities
S[k]
print(" Time[", k, "]=", Time[k],
" S[", k, "]=", S[k],
" I[", k, "]=", I[k],
" R[", k, "]=", R[k],
" F[", k, "]=", F[k])
#
################################################
Just looking at numbers is boring.
Plotting the Results!¶
Let's Plot it! Here we are going to pull each value in sol.T so we can specifically customize each color.
################################################
#
# Plot Results
#
plt.plot(Time, S, "g", # susceptables
Time, I, "r", # infected
Time, R, "b", # recovered
Time, F, "k") # fatalities
plt.title("SIRF Model System, $R_o$="+
str(Ro) +
", $p_c$=" +
str(herd_immunity) +
", $p$=" +
str(immunization_fraction))
plt.xlim(0, TimePeriod)
plt.ylim(0, (S_0 + I_0 + R_0 + F_0))
plt.xlabel('Time (Days)')
plt.ylabel("Populations")
plt.legend(["Susceptables",
"Infected",
"Recovered",
"Fatalities"])
plt.show()
#
################################################
Exploring Herd Immunity¶
Let's now explore three values of Herd Immunity. We will run three cases with 25%, 50% and 75% of the populus being vaccinated.
################################################
#
# Plot Results
#
# We're also going to go for a longer game (one year)
TimePeriod = 365.
nTime = int(TimePeriod / DeltaT + 1)
Time = np.linspace(start = 0,
stop = TimePeriod, # this is *inclusive*
num = nTime)
#
# Setting up our arrays
#
S = Time*0 # Uninfected People (multiplying time by zero)
I = Time*0 # Infected (multiplying time by zero)
R = Time*0 # Recovered (multiplying time by zero)
F = Time*0 # Fatalities (multiplying time by zero)
#
# And we can initialize our first (or zero'th since it's python)
# time step from our above initial conditions
#
k = 0
S[k] = S_0 # Uninfected People (includes Vaccinated)
I[k] = I_0 # Infected
R[k] = R_0 # Recovered
F[k] = F_0 # Recovered
# This is the only time we need to do this above initialization
# until we change the size of our arrays.
# remember X[k=0] doesn't change in our loops
for immunization_fraction in [0.00,0.25, 0.50, 0.75, 1.00]:
############################################
#
# Run the Solver
#
for k in np.arange(start = 1, stop = nTime):
S[k] = S[k-1] + dSdt(S[k-1], I[k-1], R[k-1]) * DeltaT
I[k] = I[k-1] + dIdt(S[k-1], I[k-1], R[k-1]) * DeltaT
R[k] = R[k-1] + dRdt(I[k-1]) * DeltaT
F[k] = F[k-1] + dFdt(I[k-1]) * DeltaT
#
# For this graph we will also want to sum up the total
# number of people who have contracted the disease
#
Total_Impacted_Cases = I + R + F
#
# Plot Each Graph
#
plt.plot(Time, S, "g", # susceptables
Time, I, "r", # infected
Time, R, "b", # recovered
Time, F, "k", # fatalities
Time, Total_Impacted_Cases, "m") # impacted cases
plt.xlim(0, TimePeriod)
plt.ylim(0,(S_0 + I_0 + R_0 + F_0))
plt.title("SIRF Model System, $R_o$="+
str(Ro) +
", $p_c$=" +
str(herd_immunity) +
", $p$=" +
str(immunization_fraction))
plt.xlabel('Time (Days)')
plt.ylabel("Populations")
plt.legend(["Susceptables",
"Infected",
"Recovered",
"Fatalities",
"Impacted"])
plt.show()
#
############################################
#
################################################
Since this is more about learning programming tricks in Python rather than learning we can do this again... not with a loop but with a 4-panel.
I'm going to "steal" the 4-panel section from our plotting tutorial for a lot of this material because I forget this operation way too often.
#############################################
#
# Make a 2x2 plot with differennt logs on the
# axes to demonstrate leveraging the ax-and-fig
# method of plotting.
#
#
# We'll start with the subplot command which
# will let us make 2 rows (first argument)
# and 2 columns.
#
# The second modifier pair instruct the matplotlib
# "share" the axis attributes for each plot.
#
# When we do this, FIG will represent the atributes
# of the whole plotting area.
#
# The acutal "subplots" are typically called "ax"
# which you can consider to be the
#
fig, ax = plt.subplots(2, # 2 rows in the y
2, # 2 cols in the x
sharex='col', # set the colum (x axes to match)
sharey='row') # set the rows (y axes to match)
# # # # # # # # # # # # # # # # # # # # # # # # #
# Let's have our Top Left (0,0) plot represent the p=0 case
# # # # # # # # # # # # # # # # # # # # # # # # #
immunization_fraction = 0.00
for k in np.arange(start = 1, stop = nTime):
S[k] = S[k-1] + dSdt(S[k-1], I[k-1], R[k-1]) * DeltaT
I[k] = I[k-1] + dIdt(S[k-1], I[k-1], R[k-1]) * DeltaT
R[k] = R[k-1] + dRdt(I[k-1]) * DeltaT
F[k] = F[k-1] + dFdt(I[k-1]) * DeltaT
#
# For this graph we will also want to sum up the total
# number of people who have contracted the disease
#
Total_Impacted_Cases = I + R + F
ax[0, 0].plot(Time, S, "g", # susceptables
Time, I, "r", # infected
Time, R, "b", # recovered
Time, F, "k", # fatalities
Time, Total_Impacted_Cases, "m") # impacted cases
ax[0, 0].set_title("$p$="+str(immunization_fraction),
loc='left')
# Let's have our Top Right (0, 1) plot represent the p=0.25 case
# we are only going to put the legend on one pannel
# (this Lower Right one)
immunization_fraction = 0.25
for k in np.arange(start = 1, stop = nTime):
S[k] = S[k-1] + dSdt(S[k-1], I[k-1], R[k-1]) * DeltaT
I[k] = I[k-1] + dIdt(S[k-1], I[k-1], R[k-1]) * DeltaT
R[k] = R[k-1] + dRdt(I[k-1]) * DeltaT
F[k] = F[k-1] + dFdt(I[k-1]) * DeltaT
#
# For this graph we will also want to sum up the total
# number of people who have contracted the disease
#
Total_Impacted_Cases = I + R + F
ax[0, 1].plot(Time, S, "g", # susceptables
Time, I, "r", # infected
Time, R, "b", # recovered
Time, F, "k", # fatalities
Time, Total_Impacted_Cases, "m") # impacted cases
ax[0, 1].set_title("$p$="+str(immunization_fraction),
loc='right')
# # # # # # # # # # # # # # # # # # # # # # # # #
# Let's have our Lower Left (1, 0) plot represent the p=0.25 case
# # # # # # # # # # # # # # # # # # # # # # # # #
immunization_fraction = 0.50
for k in np.arange(start = 1, stop = nTime):
S[k] = S[k-1] + dSdt(S[k-1], I[k-1], R[k-1]) * DeltaT
I[k] = I[k-1] + dIdt(S[k-1], I[k-1], R[k-1]) * DeltaT
R[k] = R[k-1] + dRdt(I[k-1]) * DeltaT
F[k] = F[k-1] + dFdt(I[k-1]) * DeltaT
#
# For this graph we will also want to sum up the total
# number of people who have contracted the disease
#
Total_Impacted_Cases = I + R + F
ax[1, 0].plot(Time, S, "g", # susceptables
Time, I, "r", # infected
Time, R, "b", # recovered
Time, F, "k", # fatalities
Time, Total_Impacted_Cases, "m") # impacted cases
ax[1, 0].set_title("$p$="+str(immunization_fraction),
loc='left')
# # # # # # # # # # # # # # # # # # # # # # # # #
# Let's have our Lower Right (1, 1) plot represent the p=0.75 case
# which for us is on the safer side of our needed herd immunity
# the legend goes here because it's rather big.
# # # # # # # # # # # # # # # # # # # # # # # # #
immunization_fraction = 0.75
for k in np.arange(start = 1, stop = nTime):
S[k] = S[k-1] + dSdt(S[k-1], I[k-1], R[k-1]) * DeltaT
I[k] = I[k-1] + dIdt(S[k-1], I[k-1], R[k-1]) * DeltaT
R[k] = R[k-1] + dRdt(I[k-1]) * DeltaT
F[k] = F[k-1] + dFdt(I[k-1]) * DeltaT
#
# For this graph we will also want to sum up the total
# number of people who have contracted the disease
#
Total_Impacted_Cases = I + R + F
ax[1, 1].plot(Time, S, "g", # susceptables
Time, I, "r", # infected
Time, R, "b", # recovered
Time, F, "k", # fatalities
Time, Total_Impacted_Cases, "m") # impacted cases
ax[1, 1].set_title("$p$="+str(immunization_fraction),
loc='right')
ax[1, 1].legend(["Susceptables",
"Infected",
"Recovered",
"Fatalities",
"Total Impacted"])
# add the big title
fig.suptitle("SIRF Model System, $R_o$="+
str(Ro) +
", $p_c$=" +
str(herd_immunity))
# you have to kludge a common axis here which is
# rather disapointing (maybe the next version of
# matplotlib will have it.)
#
# we use the fig.text function
# the first two arguments are the x and y centers of
# the plot in the larger plotter space
#
fig.text(0.50, 0.00,
'Time (days)',
ha = 'center')
fig.text(0.00, 0.50,
'Populations',
va = 'center',
rotation = 'vertical')
plt.show()
#
############################################
Exploring Herd Immunity¶
Let's look more closely at how the herd immunity chagnes the overall impact of the diseases (counting results)
This also shows how to loop through the system and load the results into an array
################################################
#
# Plot Results
#
#
# let's first make an array to carry our
# immune factor
#
# Make the immunity (dp = 0.01 so we have 101 slots)
immune_array = np.linspace(0, 1,101)
# clone the immune array as zeros for the impact array
impact_array = np.zeros(immune_array.size)
for i in range(0, immune_array.size):
####################################
#
# Load Immunization Fraction
#
immunization_fraction = immune_array[i]
#
# Run the ODE Solver
#
for k in np.arange(start = 1, stop = nTime):
S[k] = S[k-1] + dSdt(S[k-1], I[k-1], R[k-1]) * DeltaT
I[k] = I[k-1] + dIdt(S[k-1], I[k-1], R[k-1]) * DeltaT
R[k] = R[k-1] + dRdt(I[k-1]) * DeltaT
F[k] = F[k-1] + dFdt(I[k-1]) * DeltaT
#
# For this graph we will also want to sum up the total
# number of people who have contracted the disease
#
Total_Impacted_Cases = I + R + F
#
#
impact_array[i] = Total_Impacted_Cases[nTime-1]
#
####################################
plt.plot(immune_array,
impact_array,
"m")
plt.vlines(herd_immunity,
0, (S_0 + I_0 + R_0 + F_0),
'b')
plt.legend(["Total Impacted",
"Herd Immunity"])
plt.xlim(0, 1)
plt.ylim(0, (S_0 + I_0 + R_0 + F_0))
plt.title("SIRF Model System, $R_o$="+
str(Ro) +
", $p_c$=" +
str(herd_immunity))
plt.xlabel("Fraction of Population Immunized, $p$")
plt.ylabel("Total Impacted by Disease")
plt.show()
#
################################################
So if you don't have a legitimate medical reason not to take a given immunization (e.g., alergies to the flu vaccines, immunocompromised, etc.), GET YOUR FRIGGIN' SHOTS!!!!!
(and if your curious about the title of the exercise... click below.)
#Tom Lehrer Sing's "I got it from Agness"
YouTubeVideo('R6qFG0uop9k')
Version Information¶
################################################################
#
# Loading Version Information
#
%load_ext version_information
%version_information version_information, numpy, scipy, matplotlib
#
################################################################