Lab 3 - The Spread of Epidemics Solutions

Names:

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.  Then move the cursor to the right of the second red line, and have your partner type his or her 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.

Report

No formal report is required for this laboratory.  You will need to provide answers to the questions asked in the gray cell boxes provided.  (The box will expand to fit your answer.)  Please write answers as short paragraphs, using complete sentences.

Name the notebook with your email addresses (separated by a period), followed by .lab3.  For instance, students whose email addresses were 1student and 2 student should name the notebook

1student.2student.lab3

You will need to save the notebook to one of the following places:
1)	A floppy disk you provide.
2)	Temporary storage of some sort on the hard drive of the machine you are using (if it is not your own computer), or to your computer's hard drive.
3)	Temporary storage on the Academic Computing Server.

You should save your notebook frequently as you work through the laboratory.  If your computer crashes, you will at least be able to get back to where you were when you last saved the lab.

Once you have finished your laboratory, drag the icon for the notebook from the place you saved it to the Homework Drop folder for our course on the CS/Math Server.  You will get a message that you do not have enough access priviledges to see what you drop in the folder.  Say OK to this message, and the file will be copied.

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:  As β increases, the epidemic occurs earlier, runs for a shorter period of time, and reaches a greater peak value.  If β decreases below λ, there is a major change in the behavior.  The graph of the infecteds then becomes strictly decreasing.

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

ANS:  Increasing β makes the first (positive) term in I'(t) larger.  This means that the growth of the infecteds should occur more rapidly.  Since the second (negative) term then gets larger faster, and the susceptibles in the positive term decreases rapidly, we then get a rapid decrease in the number of infecteds and so the duration of the epidemic should shorten.  This is exactly what we see when we plot the numerical solutions.

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:  The effect of decreasing λ turns out to be essentially the same as increasing β, and vice-versa.  If λ decreases, the epidemic occurs earlier, becomes shorter, and has a greater maximum.  Increasing λ does the opposite, until it passes β.  Once larger than β we get only a strictly decreasing I(t).

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

ANS:  Decreasing λ reduces the size of the negative term in I'(t), allowing for more rapid growth.  Again, this makes S(t) decrease more rapidly as well, so we then get a shorter epidemic, but with a greater peak.  This is exactly what we saw in the numerical solution.

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:  This was answered above.  The value where the change in the character of the graph occurs is when λ=β.

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

ANS:  Since at t = 0 we have S(t) ≈ n, the positive term in I'(t) is approximately β I(t), and the negative term is λ I(t).  These essentially cancel if β=λ.  (They don't exactly cancel since S(0) < n, which is why I(t) still goes to zero rather than staying constant.)

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:  As i0 increases, the epidemic occurs earlier.  The peak value does not change much until i0 gets sufficiently large, at which point it increases as well.

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

ANS:  The first (positive) term in I'(t) at t = 0 is essentially concave down quadratic in I(0), since S(0) = n - I(0).  This means that increasing I(0) will make the first term grow a bit faster than the second term, at least at first.  This shifts the epidemic earlier, and causes a greater peak.  Eventually, when I(0) is about n/2, we would then see a decrease in the first term as I(0) increases, meaning the graph of the infecteds should begin to look strictly decreasing.  This actually occurs at a higher i0 than I asked you to use, so you might not have seen this effect.  The model does indicate what we saw in the numerical solution.

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=1800;
hshift=50;
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:  The model seems to do a good job, provided our assumption about excess deaths is correct.  It does not do quite as good a job after the 118 day or so mark, however, and even is a bit off at the 80 day mark.  The match is good enough for me to believe that the model does approximate this flu epidemic, and might be useful in studing other epidemics.

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:  One assumption I am not entirely happy with is the idea that once you have recovered, you can no longer get the disease.  We might need to factor in a "relapse" term in the model to be a bit more accurate.  Another thought is that the coefficients β and λ (especially β) might be time-dependent, as people begin to take precautions to avoid exposure to the virus.  This might be especially true if we consider the Aids epidemic.  For this virus, the value of λ is quite tiny (people do not recover naturally, and only recently has there been progress in reducing the effects of the virus in an infected individual).  The main obstacle to the spread of this disease is the effort people take to avoid being exposed to the virus.  As a first approximation, however, the model seems to do quite well.

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.277;
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:  The fraction is between 0.276 and 0.277.  We need to have about 27.7% of the population immune to prevent an influenza epidemic.

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?

.277/.95

ANS:  About 30% (29.2%) of the population would have to be vaccinated in order to prevent the epidemic.

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:  While it is feasible to vaccinate 30% of the population against influenza, most people do not regard this disease as important enough to get vaccinations.  In addition, there are different strains of the virus, which require different vaccines.  As to the situation in 1968, I doubt that there was enough lead time to have prevented the epidemic from happening.

Measles

For measles, the value of β is about 2.0, and the value of λ is about 0.13.  (The higher β indicates a more virulent disease, while the lower λ indicates a longer infection time.)  You will need to change the value for fraction in the next cell.  (We will continue to use the 1968 NYC population size while doing this, although the actual population size is essentially irrevelant.)

Clear[s, i, r, t, β, n, λ, i0, r0, tend, 
      solution, fraction]
fraction=0.935;
n = 7900000;
β = 
        2.00;
λ = 0.13;
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 a measles epidemic?

ANS:  The fraction is between 0.934 and 0.935.  This means that about 93.5% of the population needs to be immune to prevent a measles epidemic.

Assuming there is a vaccine that is 95% effective against measles (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?

.935/ .95

ANS:  We need to have about 98.5% of the population vaccinated to prevent an epidemic.  This figure, although high, has almost been reached in elementary school students.

There is a program that innoculates elementary school students against measles.  Current estimates are that about 97% of students entering elementary school are innoculated against measles.  Is this sufficient to prevent epidemics?  (By the way, there have been some outbreaks at colleges and universities in the past few years.  This particular program was not in effect at the time most current college age students were in elementary school.)

ANS:  This is almost enough, but slighly short of what we need.   As compliance with the program improves, the 98.5% figure might be reached soon.

Rubella

For rubella ("German measles"), the value of β is about 0.66, and the value of λ is about 0.09.  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.864;
n = 7900000;
β = 
        0.66;
λ = 0.09;
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 a rubella epidemic?

ANS:  The fraction is between 0.863 and 0.864.  This means that about 86.4% of the population needs to be immune to prevent a rubella epidemic.

Assuming there is a vaccine that is 95% effective against rubella (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?

.864/.95

ANS:  We would need to vaccinate about 91% of the population to prevent a rubella epidemic.

Which disease appears easier to control; measles or rubella?

ANS:  Of the two, rubella seems easier to control.  Neither is anywhere as easy to control as a particular strain of influenza, but the stakes are somewhat higher with these diseases.

If our elementary school program also provides innoculations against rubella, and the same 97% of the incoming students have been innoculated, should this be effective at preventing future rubella epidemics?

ANS:  That percentage would do the trick for rubella.

More about finding the fraction.

You should be able to use the system of differential equations to find the value of ((Subscript[r, 0])/n)(fraction of population) that you estimated graphically for each of the diseases above.  How should this fraction relate to β and λ?

{1 - .34/.47, 1 - .13/2, 1 - .09/.66}

ANS:  In order to have a strictly decreasing graph for I(t), we need I'(0) to be negative.  This will happen if S(0)/n β = λ.  Since S(0) = n - I(0) - R(0), S(0)/n is approximately 1 - R(0)/n.  The fraction R(0)/n should then be about 1 - λ/β.  The values of this for each of the diseases are in the cell above this one.  These agree almost perfectly with the estimated values from the graphs.