The Spread of Epidemics

Name:

Place your names in the space below.  Use the mouse to put the cursor just to the right of the red line, and type your name.

>>>|

Purpose

In this laboratory, we will look at the S-I-R Model for the spread of an epidemic.  We will see the effects of varying the epidemic parameters, and then see if the model matches historical data from the 1968 Hong Kong Flu epidemic for some choice of these parameters.  Finally, we will use the model to predict the success of various types of innoculation programs on three diseases.

Preliminaries

The S-I-R Model for the spread of epidemics uses a system of three differential equations to predict the sizes of three different classifications of people in the population, namely susceptibles, infecteds, and recovereds.  We make some simplifying assumptions about the epidemic.  First, we assume that initially all members of the population are either susceptible to the disease, or actually have the disease.  We assume the disease is passed by contacts between susceptibles and infecteds at a rate proportional to the product of their numbers.  A fixed fraction of the infecteds recover from the disease each day.  Once an individual recovers from the disease, we assume he or she cannot get it again, at least during this round of the epidemic.

If we let s(t) be the number of susceptibles, i(t) the number of infecteds, r(t) be the number of recovereds, n be the population size, β be the number of contacts per day in which an infected person actually gives the disease to a susceptible person , and λ be the percentage of people who recover each day, then the system of differential equations for this model is

s' (t) = -(β/n) s (t) i (t)
i'(t) = (β/n) s (t) i (t)- λ i(t)
r'(t) = λ i(t).

We have an n in the denominator because ((s(t))/n) is the fraction of people an infected individual meets in that are susceptible to the disease.  The value of λ is roughly the reciprocal of the amount of time it takes the average infected individual to recover from the disease.

Plotting solutions of systems of equations uses the Mathematica command NDSolve to create the functions to be plotted. Its syntax is very similar to that of DSolve.  The next cell shows how we can use this command to reproduce a plot of the susceptibles (in green), the infecteds (in red) and the recovereds (in blue).

Clear[s, i, r, t, β, n, γ, tend, 
      solution]
n = 
        7900000;
β =  .6;
λ =  .34;
i0 = 10;
tend=135;
solution = 
        NDSolve[{s'[t] == -(β/n) s[t] i[t], 
            s[0] == 
              n - i0, i'[
                t] == (β/n) s[t] i[t] - λ i[t], 
            i[0] == i0, r'[t] == λ i[t], 
            r[0] == 0}, {s[t], i[t], r[t]}, {t, 0,  
            tend}];
Plot[
        Evaluate[{s[t], i[t], r[t]}  /.  solution], {t, 0, tend}, 
        PlotStyle -> {RGBColor[0, 1, 0], RGBColor[1, 0, 0], 
            RGBColor[0, 0, 1]}];

Project

In the input cells that follow, you will need to adjust any values that appear in magenta.  You also may wish to adjust any values in blue if doing so produces a better graph for your analysis.  Since all graphs are self-scaling vertically, be sure to keep an eye on the vertical scale when you analyze the effect of changes you make to the parameters.

What happens when you change β?

Keep λ fixed, but see what happens as you vary β between 0.1 and 2.0.

Clear[s, i, r, t, β, n, λ, i0, tend, 
      solution]
n = 7900000;
β=.47;
λ =  .34;
i0 = 10;
tend=135;
solution = 
        NDSolve[{s'[t] == -(β/n) s[t] i[t], 
            s[0] == 
              n - i0, i'[
                t] == (β/n) s[t] i[t] - λ i[t], 
            i[0] == i0, r'[t] == λ i[t], 
            r[0] == 0}, {s[t], i[t], r[t]}, {t, 0,  
            tend}];
Plot[
        Evaluate[i[t]  /.  solution], {t, 0, tend}, 
        PlotStyle -> RGBColor[1, 0, 0]];

What happens as you change β?  You may need to increase the time interval for some of the graphs.  Be sure to stay alert for automatic changes in the scale on the vertical axis.

ANS:

Do these changes make sense from an intuitive understanding of the model?  Explain.

ANS:

What happens when you change λ?

This time, keep β fixed, but see what happens as you vary λ between 0.1 and .65.

Clear[s, i, r, t, β, n, λ, i0, tend, 
      solution]
n = 7900000;
β =  .47;
λ=.34;
i0 = 10;
tend=135;
solution = 
        NDSolve[{s'[t] == -(β/n) s[t] i[t], 
            s[0] == 
              n - i0, i'[
                t] == (β/n) s[t] i[t] - λ i[t], 
            i[0] == i0, r'[t] == λ i[t], 
            r[0] == 0}, {s[t], i[t], r[t]}, {t, 0,  
            tend}];
Plot[
        Evaluate[i[t]  /.  solution], {t, 0, tend}, 
        PlotStyle -> RGBColor[1, 0, 0]];

What happens as we change λ?  You may need to increase the time interval for some of the graphs.  Be sure to stay alert for automatic changes in the scale on the vertical axis.

ANS:

Do these changes make sense from an intuitive understanding of the model?  Explain.

ANS:

There should have been a major change in the character of the graph of i(t) at a particular value of λ.  What was that change, and at what λ value did it occur?

ANS:

Explain how you can find this λ value from the system of differential equations.

ANS:

What happens when we change the number of initial infecteds?

Finally, we will change the number of initial infecteds, but leave β and λ alone.  Vary i0 between 1 and 1,000,000, using powers of 10.

Clear[s, i, r, t, β, n, λ, i0, tend, 
      solution]
n = 
        7900000;
β =  .47;
λ =  .34;
i0=10;
tend=135;
solution = 
        NDSolve[{s'[t] == -(β/n) s[t] i[t], 
            s[0] == 
              n - i0, i'[
                t] == (β/n) s[t] i[t] - λ i[t], 
            i[0] == i0, r'[t] == λ i[t], 
            r[0] == 0}, {s[t], i[t], r[t]}, {t, 0,  
            tend}];
Plot[
        Evaluate[i[t]  /.  solution], {t, 0, tend}, 
        PlotStyle -> RGBColor[1, 0, 0]];

What happens as we change i0?  Be sure to stay alert for automatic changes in the scale on the vertical axis.

ANS:

Do these changes make sense from an intuitive understanding of the model?  Explain.

ANS:

Matching the data from the 1968 NYC Hong Kong Flu epidemic.

In 1968, New York City was hit by a strain of flu nicknamed "Hong Kong."  This was a new strain (at the time) to which people had no prior immunity.  The disease rapidly reached epidemic proportions.

We start this section with an excerpt from Keith D. Stroyan, Calculus using Mathematica, Scientific Projects and Mathematical Background, Project 3.

In the study of virulent diseases, one of the biggest problems facing investigators is the availability of sufficient data for analysis.  In order to model a disease, we need to know the proportions of the population that are susceptible, immune, and removed.  However, it is often very hard to determine these proportions because a large number of cases of a given disease go unreported (consider how many times you’ve had the flu, and out of those times, how many times you actually saw a doctor and so were “counted” as an infected).  For this reason, certain diseases lend themselves to study because, in one way or another, they make their presence known to the medical community.  For example, in the case of polio, the disease was severe enough to require medical attention, ensuring that there would be records of all cases.  In the case of the Hong Kong flu, which you will be studying here, the data was derived from a slightly different source.

The data for the Hong Kong flu is derived from the medical reports of coroners in New York City.  During the period of the epidemic, the course of the disease can be traced as a function of the excess of individuals who died due to influenza-related pneumonia in the New York City area.  The term “excess” in the above sentence refers to the fact that, in New York City, you expect to get a certain number deaths due to influenza-related pneumonia.   However, during the epidemic, there was a large increase in influenza-related deaths which could be attributed to the Hong Kong flu.

Once we have our data, the question becomes, how to interpret it?  There are two factors that we must consider.  First, not everyone who gets the Hong Kong flu is going to die as a result.  And second, those who do die as a result probably will not do so right away, but rather will have the disease for an extended period of time before expiring.  On the basis of these two facts, then, we can hypothesize that the number of excess deaths is proportional to the number of infected individuals from a couple of weeks before [the emphasis is ours].

The next cell enters the "excess deaths" data into Mathematica.  I set it to be an input cell, so it has already been evaluated, provided you answered "Yes" to the initialization request.

Clear[days, deaths]\(days = 7\ Range[13]; 

    deaths = {14, 28, 50, 66, 156, 190, 156, 108, 68, 77, 33, 65, 24}; \)

It is important to realize that there is no significance to the starting time t=0 days in the “Excess Flu Deaths” plot.  In particular, there is no reason to assume that t=0 in this plot corresponds to the “start” of the epidemic.  Keep this in mind as we compare the S-I-R model against “Excess Flu Deaths.”

Here is a plot of infecteds from the model.  The value .47 for β  and .34 for λ are empirically derived values for influenza.

Clear[s, i, r, t, β, n, λ, i0, tend, 
      solution]
n = 7900000;
β=0.47;
λ=0.34;
i0 = 10;
tend=135;
solution = 
        NDSolve[{s'[t] == -(β/n) s[t] i[t], 
            s[0] == 
              n - i0, i'[
                t] == (β/n) s[t] i[t] - λ i[t], 
            i[0] == i0, r'[t] == λ i[t], 
            r[0] == 0}, {s[t], i[t], r[t]}, {t, 0,  
            tend}];
infectedPlot = 
        Plot[Evaluate[i[t]  /.  solution], {t, 0, tend}, 
          PlotStyle -> RGBColor[1, 0, 0]];

The next cell superimposes a plot of the "excess deaths" data on the plot we got above.  You will need to adjust the vscale and hshift to get the graphs to (roughly) match.  Recall that the "excess death" data does not necessarily begin at the starting time of the epidemic, so we may need a positive shift.  The red plot is i(t), and the blue plot is the scaled and shifted "excess deaths" plot.

vscale=1000;
hshift=0;
deathPlot = 
        ListPlot[Transpose[{days + hshift, vscale deaths}], 
          PlotJoined -> True, PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
Show[infectedPlot, deathPlot, 
        DisplayFunction -> $DisplayFunction];

You may wish to go back and change β and/or λ to see if you get a better match.  If you do, be sure to evaluate both input cells that produce the plots (namely, the two input cells above this cell).

Do you think the model approximates the actual 1968 epidemic?  Why or why not?  Keep in mind how the graph above was produced.

ANS:

Your thoughts on the model.

Use the cell below to write a few sentences about the model.  Which assumptions seem reasonable, and which need to be improved?  Describe at least two additional features of real epidemics you might want to try and include in an improved model.

ANS:

Using vaccination programs to prevent epidemics.

The basic idea behind a vaccination program is to move people directly into the "recovered" category without having them get the disease.  An epidemic can be prevented if we innoculate a large enough percentage of the population to make i(t) a strictly decreasing function of t.   We can model the effect of innoculating a percentage of the population by simply changing the initial conditions.  The number of people who have been made immune to the disease are subtracted from the initial number of susceptibles, and added to the initial number of recovereds.

For each of the following diseases, we want to estimate the fraction of the population that has to be innoculated in order to prevent an epidemic from happening.  Since not all innoculations work, we will first find the fraction of the population that has to be immune to the disease, and then adjust for the effectiveness of the vaccine.  To find the smallest fraction that must be immune, we can do a "bisection" type approach.  If the fraction is 0, there will be an epidemic, while if the fraction is 1 there will not be one.  The smallest fraction that has to be immune is then between 0 and 1.  We can test what happens if the fraction is 0.5, and find out if that is sufficient.  We can continue to bisect to get a decent approximation.  If the fraction is sufficient, i(t) will be strictly decreasing.  Remember, we want to find the smallest fraction we would need.

Influenza

For this, we will assume β is 0.47 and λ is 0.34.  You will need to change the value for fraction in the next cell.

Clear[s, i, r, t, β, n, λ, i0, r0, tend, 
      solution, fraction]
fraction=0.5;
n = 7900000;
β = 
        0.47;
λ = 0.34;
i0 = 
        10;
r0 = n fraction;
tend=135;
solution = 
        NDSolve[{s'[t] == -(β/n) s[t] i[t], 
            s[0] == 
              n - i0 - r0, i'[
                t] == (β/n) s[t] i[t] - λ i[t], 
            i[0] == i0, r'[t] == λ i[t], 
            r[0] == r0}, {s[t], i[t], r[t]}, {t, 0,  
            tend}];
Plot[
        Evaluate[i[t]  /.  solution], {t, 0, tend}, 
        PlotStyle -> RGBColor[1, 0, 0]];

What fraction of the population must be immune in order to prevent an influenza epidemic?

ANS:

Assuming there is a vaccine that is 95% effective against influenza (meaning that 95% of the people who get the vaccine become immune to the disease), what percentage of the population would need to be innoculated in order to prevent an epidemic from beginning?

ANS:

If this vaccine existed in 1968, do you think it would have been possible to prevent the epidemic?  You should include whether or not this would be feasible if there was only a short time period to do the innoculations.

ANS: