SIR Peak Infection Challenge

You are the chief epidemiologist of the McMurdo Antartic Research Station, a scientific outpost with a population of 1000 persons. A new and unusual disease, currently known as penguin foot, has struck 10 individuals in the camp, and seems to be spreading.

The disease has several properties. While infected (and contagious), individuals behave as normal but are unable to drive snow machines, the base’s primary mode of evacuation in case of an emergency. (It is universally acknowledged that penguins cannot drive snow machines.) Snow machines can carry one passenger in addition to the driver, and so at any given time, at least half of the population must be able to drive.

The period of infection is on average 10 days, after which individuals return to their normal state and are immune to further infection. The infectivity of the disease is much less well known. The chance of infection from an encounter between an infected and a non-infected person could be anywhere from 0.5% to 10%.

As there are no lasting consequences of penguin foot, the base has continued to operate as normal, with each individual coming into contact with about 25 people per day, with relatively uniform mixing between individuals.

The primary danger posed by the disease is that there will come a point when there aren’t enough people to drive snow machines. Your task is to calculate the maximum number of infections for a range of values of infectivity between .005 and 0.1.

Specifically, you need to know the minimum value of infectivity that is sufficient to incapacitate over half the base at one time, so that as new measurements become available, you can assess the risk of this possibility.

Load libraries

The cell below gets us started by loading the plotting library.

Add a line to this cell to load the pysd library using the standard python import syntax.

%pylab inline
import pysd

Populating the interactive namespace from numpy and matplotlib

Load model

Now, use PySD’s read_vensim function to read the SIR model found at relative location: ‘../../models/Epidemic/SIR.mdl’, and name the resulting object ‘model’.

model = pysd.read_vensim('../../models/Epidemic/SIR.mdl')

Just to practice, run the model, using the .run() function’s params argument to set the value of Infectivity to 0.02.

Save the result in a variable named res.

res = model.run(params={'Infectivity': 0.02})

The cell below contains code to plot the result of your simulation. Execute it and ensure that everything is working properly.

res['Infected'].plot()
plt.xlabel('Days')
plt.ylabel('Number of Individuals Infected');

Identify the peak for this base case run

In the cell below, write an expression to calculate (for your base case) the maximum number of individuals who are infected at any one time. You’ll need to use Pandas syntax to select the column named ‘Infected’ and then use the Pandas functions referenced above.

peak_value = res['Infected'].max()

print('Up to', int(peak_value), 'individuals are infected at one time.')

Up to 4798396 individuals are infected at one time.

Define the range of infectivities we want to sweep over

To identify the worst case scenario, we need to sweep over the plausible values of infectivity, from 0.005 to 0.1, in increments of .005. Our next step is to generate an array of these values. We can use the python package numpy which handles matrix mathematics and array manipulation. (It is common practice to give the numpy module the short handle ‘np’ as I have done below).

You’ll specifically want to use the np.arange(...) function, which extends the python standard range function to handle non-integer values. Consult the numpy documentation to determine the arguments you’ll want to pass in to this function to generate an array that looks like:

[0.005, 0.01, 0.015, ... 0.095, 0.1]

You may have to be creative to ensure that the last value in the array is actually 0.1.

import numpy as np
infectivity_values = np.arange(.005, .105, .005)
print(infectivity_values)

[0.005 0.01  0.015 0.02  0.025 0.03  0.035 0.04  0.045 0.05  0.055 0.06
 0.065 0.07  0.075 0.08  0.085 0.09  0.095 0.1  ]

Evaluate the peak of infections for each value of infectivity

We now want to calculate the peak for the list of possible infectivities, and collect them in a pair of lists. To do this, write a for loop that iterates over each value in the array of infectivites. I’ve started this off for you below.

Within the body of the for loop, run the model with that value of infectivity, saving the result to a variable named res. Then use python’s list append syntax to add the appropriate values to the end of the peak_value_list.

peak_value_list = []

for inf in infectivity_values:
    res = model.run(params={'Infectivity': inf})
    peak_value_list.append(res['Infected'].max())

peak_value_list

[34.86810210467729,
 17376.295842009604,
 3127287.127264763,
 4798396.222053367,
 5490586.578913281,
 6007712.289673515,
 6410082.172490937,
 6733269.976059377,
 6999008.729515837,
 7221809.038370998,
 7412060.684035444,
 7576357.8468962945,
 7719951.693505733,
 7846142.103150383,
 7958727.643037483,
 8060110.328538261,
 8150974.284047906,
 8233434.777704702,
 8308822.446865912,
 8377904.218632583]

Plot the result

Now create a plot showing showing the values of infectivity on the x-axis, and the peak value of the infections on the y axis. Label each axis, and give the plot a title. From this plot we can eyeball the value of infectivity beyond which the peak level of infections rises over 500.

plt.plot(infectivity_values, peak_value_list)
plt.grid()
plt.xlabel('Infectivity')
plt.ylabel('Peak Value of Infections')
plt.title('Peak level of infection as a function of infectivity.');

Bonus Activities

1. Identify the first day that the base might not be able to evacuate, for any value of infectivity.

2. Identify the maximum total number of days that the base might be unable to evacuate, for the range of values of infectivity listed above.

3. Come up with a better story for this example than ‘penguin foot’. =)