Project Help Notebook by Keith Stroyan

Help NoteBook Copyright 1998 by Academic Press, Inc. - All rights reserved.
Modify this Notebook to prepare your project report.

Project 2:  The Hong Kong Flu Epidemic 1968-69

by A. Student and B. Pupil

with technical help from Da LabMonitor and Help Notebook

Project Introduction

Write this section last.  Begin with the basic question and state your main finding clearly, but with few technical details.  You might even modify the title of this section to make it more interesting.  Following is a rough start.

This report compares the real Hong Kong Flu epidemic of 1968-69 in New York City with the differential equation model of Chapter 2 in the main text, Calculus: The Language of Change.

We work with the best available data, the New York City coroner's estimate of "excess deaths" attributed to the strain of influenza that was first introduced in New York at that time.  These data are NOT a complete record of every person who contracted this flu. Comparison of the data and model takes this into account, as described in detail below.

"Adjusting" the data in order to make a comparison with the model, we find that ___.  But ___.  And ___...

Excess death data from the
New York Hong Kong Flu Epidemic of 1968-69

The data for the Hong Kong flu is derived from the medical reports of coroners in New York City.  During the epidemic, there was an increase in deaths that could be attributed to the Hong Kong flu.  Our data is 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, but the influenza-related pneumonia caused more than expected.

Once we have our data, the question becomes, how to interpret it?  Two factors that we must consider:

First, not everyone who gets the Hong Kong flu is going to die as a result.  We need to use this to scale the excess death data for comparison with the whole population.

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.  People also did not realize immediately that a new strain of influenza was present causing an epidemic.  We use these considerations to "shift the time" between the death observations and the model's predicitions.

After scaling and allowing for time for the epidemic to start, our superposition of model prediction and measured excess deaths looks as follows:  (you may replace this next graph once you have a better one)

Before we go further, we introduce and graph the excess death data.

The Excess Death Data

The following command cells give the excess death data from the New York epidemic.  First in weeks and then in days, so we can more easily compare the data and model.  Enter these commands now.

Excess Death Data in Weeks

Needs["Graphics`Colors`"];

deaths = {14, 31, 50, 66, 156, 190, 156, 108, 68, 77, 33, 65, 24};
weeks = {1,2,3,4,5,6,7,8,9,10,11,12,13};
weekData = Transpose[{weeks,deaths}];
TableForm[Join[{{"Week", "Deaths"}},weekData]]

Excess Death Data in Days, Starting at Zero

deaths = {14, 31, 50, 66, 156, 190, 156, 108, 68, 77, 33, 65, 24};
days = ({1,2,3,4,5,6,7,8,9,10,11,12,13} - 1) 7;
dayData = Transpose[{days,deaths}];
TableForm[Join[{{"Days", "Deaths"}},dayData]]

A Plot of the Daily Data

If you have entered the day data for excess deaths, the next command will plot the data.  (Enter it now and print your report with a graph - after deleting this message.)

deathdays = ListPlot[dayData, AxesLabel->{"Days","Deaths"},PlotJoined -> True,PlotStyle-> Red];

The next section of our report derives the model and uses it to predict the course of a flu epidemic in New York City.

S-I-R Equations

This section of our report derives the "equations of change" for an S-I-R epidemic.

Derivation of the model

Assumptions

We make some basic assumptions in formulating a model for the change of susceptible and infectious people.

ASSUMPTION 1:
S-I-R: We study only diseases that confir immunity.  This means...

ASSUMPTION 2:
Large fixed population:
Although New York City is obviously large, it is not a static population, for the sake of simplicity we will represent it as such - as a result, the only way to get out of the Susceptible group is to become an Infective, after which you must enter the Removed category, which includes both deaths and recoveries.  Since a flu epidemic takes place in __ days, the number of births, deaths, and immigrations bla, bla, blah...

ASSUMPTION 3:
Well mixed bla, bla, blah...  this assumption is used...

Variables

Our model uses the following variables

t = time in days
	s = susceptible fraction  = number susceptible/total population
	i = infectious fraction = number infectious/total population
	r = 1 - s - i = the "removed" fraction, those who are immune and cannot transmit the disease

You may wish to use auxillary variables as follows:

S = number of susceptible people = s*n (see n below in parameters)
	I = number of infectious people = i*n
	R = number of removed (immune) people = r*n

Parameters

b = 1/the number of days a person is infectious
	c = the number of contacts an infectious person makes that could transmit the disease
	a = b c, the daily rate of making contacts
	n = the fixed total population, n = S + I + R (see Assumption 2 above.)

Bla, bla, blah.

We can measure the paremeters  b  and  c for a disease by ....

Daily rate of recovery

If the population consists of  n people total and the fraction infectious is  i,  i*n people are infectious.  The average number of people who recover per day is b* i * n, because bla, bla, blah.

Daily rate of transmission

The daily rate of making "close contacts" is given by the measured parameters b and c,  a = b*c.  If the infectious fraction is  i and the total population is  n,  the total number of contacts is  Bla, bla, blah.

daily rate of close contacts = a * i * n = a I

The fraction  s  of these contacts is made with susceptible people, resulting in new cases.

daily rate of new cases = a * i * n * s = a I (S/n) = a I s

Therefore, the rate of of the fraction of the population

You need some additional explanation here to arrive at the differential equations:

The S-I-R differential equations

(ds/dt) = -a s i

(di/dt) = a s i - b i

Computation with the model or "predicitions"

Initial conditions

An analysis requires a starting condition.  As you can see in the command cell below, the population of New York is around 7,900,000.  Because the Hong Kong flu was a new strain of influenza at the time that the data was taken, the initial condition prior to exposure to the disease was such that all individuals were susceptible.  Due to its great virulency, the disease doesn't require much get started - say 10 people who have the virus move into town, raising the population to 7,900,010.

We used the program from the class notebook SIRsolver to make predicitons with our model differential equations. (You will want to make additional modifications an leave the computations in your report, and then remove this remark.)

Clear[s, i, t, soln, sRate, 
        iRate];
sRate[s_, i_] := -a s i;
\(iRate[s_, i_] := \((a\ s - b)\)\ i; 


      b\  = \ N[1/3.0]\ ; \ \  (*b\  = \ 1/period\ of\ infectiousness*) 

      c\  = \ N[1.4]; 		 (*c\  = \ contact\ number*) 
a\  = \ b\ c\ ;\)		  
startUchange={s[0] == (7900000/7900010),i[0] == (1/790001),s'[t]==sRate[s[t], i[t]],i'[t]==iRate[s[t], i[t]]};
finalTime = 60;
soln = 
        NDSolve[startUchange, {s[t], i[t]}, {t, 0, 
            finalTime}];
Print["Spread of Hong Kong Flu in New \
York."];
Print["Contact number, c = ", 
        c, ", Infectious period 1/b =", (1/b)];
s[t_] = 
        s[t] /. soln[[1]];
i[t_] = 
        i[t] /. soln[[1]];
sPlot  =  
        Plot[s[t], {t, 0, finalTime}, PlotStyle -> {Green}, 
          AxesLabel -> {"t", "s"}];
iPlot  =  
          Plot[i[t], {t, 0, finalTime}, PlotStyle -> {Blue}, 
            AxesLabel -> {"t", "i"}];

Print["s[", finalTime, "] = ", 
        s[finalTime], "    i[", finalTime, "]=  ", 
        i[finalTime]];

The (modified) graph above shows the complete course predicted by out differential equations for a flu epidemic, from a beginning with almost noone infectious to the time when the disease has died out.

The next section compares this prediction with the data.

Matching data and the model predictions

Question on the Data:  When do the deaths peak?

your answer?

Show[deathdays];
TableForm[Join[{{"Days",  "Deaths"}}, dayData]]

Question 2:  When do the model infectives peak?

your answer?  (Note:  You will have to run the model longer to see the peak.)

Show[iPlot];
{i[85], i[86], i[87], i[88], i[89], i[90], i[91], i[92]}

FindMinimum[-i[t], {t, 90}]

The Daily Data, slid over by enough days to match the peak.

Slide the time scale to match the peak. (You can do better than the following sample which we included to help with the computing.  Fix it and remove this statement.)

slide = 40;(*Adjust this so the model peak and data peak match.*)
deaths = {14, 31, 50, 66, 156, 190, 156, 108, 68, 77, 33, 65, 24};
times = ({1,2,3,4,5,6,7,8,9,10,11,12,13} - 1) 7 + slide;
adjustedData = Transpose[{times,deaths}];
deathgraph = ListPlot[adjustedData, AxesLabel->{"Days","Deaths"},PlotJoined -> True,PlotStyle-> Red];
Show[iPlot];

Comparison of the model & the data

Once you have run the model long enough to see where the peak occurs and have adjusted the time scale on the death data so that the peaks match, you still have to line up the veritcal axes.

The proportion of the population that will die from flu is not 100%.  You can either multiply the  model Infectious numbers by a big amount or divide the death data by that amount.  Notice the vertical scale on the next graph and the one afterward: (You can do better than the following sample which we included to help with the computing.  Fix it and remove this statement.)

shrink = 3800. ; (*Adjust this so the model peak and data peak match.*)
deaths={14, 31, 50, 66, 156, 190, 156, 108, 68, 77, 33, 65, 24}/shrink;
slide = 40 ;          (*Adjust as above to match the peaks.*)
times = ({1,2,3,4,5,6,7,8,9,10,11,12,13} - 1) 7 + slide;
adjustedData = Transpose[{times,deaths}]
deathGraph = ListPlot[adjustedData, AxesLabel->{"Days","DeathFraction"},PlotJoined -> True,PlotStyle-> Red];

Show[iPlot,deathGraph];

Justification of the time shift and scaling

Above we found that by shifting the time on the excess death data by __  days and scaling the deaths to be a fraction of ___ we get fairly good agreement between prediction and reality.

The time scaling is reasonable because ___

OR

The time scaling is totally unreasonable because ___

AND

Scaling the deaths to be a fraction of ___ is reasonable because ___

OR

Scaling the deaths to be a fraction of ___ is totally unreasonable because ___

Comparison at the peak

190/shrink
i[89.7]

Bla, bla, blqa

The wiggles at the end

Any way you figure it, there's some explaining to do about the wiggles at the end of the data.

Conclusions

Give your considered conclusion here.