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 5:  Three Peaks in an S-I-R Epidemic

by A. Student and B. Pupil

with technical help from Da LabMonitor and Help Notebook

My Project Introduction

Write this 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.

HINT 5.1: 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:
Fixed population:  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...

ASSUMPTION 4:
Large population 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

HINTs 5.2 & 3: Review of Calculus & Max-min Theory

Increasing, decreasing, and the derivative

What does the derivative have to do with increasing and decreasing.  How do you answer HINT 5.2?

The Maximum of a Function

One way to find the maximum of a function is to find the place there its graph changes from increasing to decreasing. In terms of calculus, this would be ...

What if the function changes from up to down more than once??

The Maximum Number Infectious

The quantity  i  is largest when  s = 1/c, because ... (see the project after HINT 5.4)

...

A computer example

We modified the program SIRsolver to show an example of the peak.

Needs["Graphics`Colors`"];

Clear[s, i, t, soln, sRate, iRate]; sRate[s_, i_] := -a s i; iRate[s_, i_] := (a s - b) i; c = 6.8; b = (1/11.0); a = b c; startUchange={s[0] == (1/3.) ,i[0] == (1/30.) ,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}]; 
      
s[t_]\  = s[t] /. soln[\([1]\)]; 
i[t_] = i[t] /. soln[\([1]\)]; 

      Print["Spread of rubella in a large fixed population."]; 

      Print["Contact number, c = ", c, ", Infectious period 1/b =", 
        1\/b]; 
sPlot\  = \ 
        Plot[{1/c, s[t]}, {t, 0, finalTime}, PlotRange -> {0, 1/3. }, 
          PlotStyle -> {Purple, Blue}]; 

      iPlot\  = \ Plot[i[t], {t, 0, finalTime}, PlotStyle -> {Red}]; \)

It looks like the peak occurs approximately where  t=18.  Let's plot that:

\(Show[sPlot, iPlot, Graphics[Line[{{18, 0}, {18, 0.3}}]], 
      AxesLabel -> {"t", "i peak"}]; 
lastS = s[t] /. soln; 

    Print["s[18.0] = ", s[18.0], ",  1/c = ", 1/c, \ 
      ", s[18]-1/c = ", s[18.0] - 1/c]; \)

The Maximum Number of New Cases

The rate of producing new cases  is largest when  s-i = 1/c, because ... (see the project after HINT 5.5)

...

A computer example of our general result

Clear[s, i, t, soln, sRate, iRate]; sRate[s_, i_] := -a s i; iRate[s_, i_] := (a s - b) i; c = 6.8; b = (1/11.0); a = b c; startUchange={s[0] == (1/3.) ,i[0] == (1/30.) ,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}]; 

    s[t_]\  = s[t] /. soln[\([1]\)]; 
i[t_] = i[t] /. soln[\([1]\)]; 

    f[t_] = a\ s[t]\ i[t]; 

    Print["Spread of rubella in a large fixed population."]; 

    Print["Contact number, c = ", c, ", Infectious period 1/b =", 
      1\/b]; 
iPlot\  = \ Plot[i[t], {t, 0, finalTime}, PlotStyle -> {Red}]; 
    
fPlot = \ Plot[f[t], {t, 0, finalTime}, PlotStyle -> {DarkGreen}]; 

    criticalPlot = \ 
      Plot[s[t] - i[t] - 1/c, {t, 0, finalTime}, PlotStyle -> {Blue}, 
        PlotRange -> {\(- .1\),  .1}]; \)

So we look at the graph and then locate the critical value accurately:

FindRoot[s[t] - i[t] == 1/c,  {t, 10}]

So why would we use this result to plot:

\(Show[iPlot, fPlot, Graphics[Line[{{10.5, 0}, {10.5, 0.1}}]], 
      AxesLabel -> {"t", "New Case Peak"}]; 
lastS = s[t] /. soln; 

    Print["s[10.5]-i[10.5] = ", s[10.5] - i[10.5], ",  1/c = ", 1/c]; 
    \)

We could also check our analytical work by asking Mathematica to do the following computaions, because ...

Plot[f'[t], {t, 0, 20}];

FindRoot[f'[t] == 0,  {t, 10}]

Solution of the critical equation using the invariant

We know in general that the quantity

s + i - (1/c)Log[s]

is constant for all time in an S-I-R epidemic.  We can compute this constant from the initial values and then solve the system of equations

s + i - (1/c)Log[s]=k
s - i = (1/c)

... you provide the reasons

c(2s-k) = 1+Log[s]

Now we solve this equation for the case of rubella:

\(c\  = \ 6.8; 
invar[s_, i_] := s + i - Log[s]/c; 

    k\  = \ invar[s[0], i[0]]; 

    Plot[{c\ \((2  s - k)\), \ 1 + Log[s]}, {s, 0, 1}, 
      PlotStyle -> {Red, Blue}]; 

    FindRoot[c \((2\ s - k)\) == 1 + Log[s], {s,  .5}]\)

\(plot1 = Plot[{s[t], 0.229}, {t, 0, 20}, PlotStyle -> {Purple, Yellow}]; 
    
plot2 = 
      Plot[{s[t] - i[t], 1/c}, {t, 0, 20}, PlotStyle -> {Brown, Green}]; 

    Show[plot1, plot2]; \)

So we see ...

The Maximum Number Net New Cases

The rate of increase in infectious people is largest when

(s^2) - s i - ((2 s)/c) + (1/(c^2))= 0

because ...
... see HINT 5.6

A computer example

You provide the code using ideas above

Solution of the critical equation using the invariant

Solve simultaneously

(s^2) - s i - ((2 s)/c) + (1/(c^2))= 0
s + i - (1/c)Log[s]=k

... remembering that we could plot both as i = function of s, for example, at least in specific cases of  c  and  k...

Comparison of the Times of the Peaks

We did a test case of one particular disease and initial condition and found the three peaks as follows:

s = 1/c

s - i = 1/c

\(\(																
    s\^2 - s\ i - \(2\ s\)\/c + 1\/c\^2 = 0\)\)

Alltogether:

If we hope to intervene in an epidemic, which peak is most important to head off?  Which happens first?