Limiting Values for the s-i-r Model

by
K. D.  Stroyan
University of Iowa

copyright 2007 by Keith Stroyan - All rights reserved.

Notebook Overview

At all times during the course of a model s-i-r epidemic, the  quantity

i + s - Log[s]/c        Log[s]
i + s - ------
          c

is constant.  The parameter  c  is the constant contact number,  c = a/b, but  s = s[t]  and  i = i[t]  vary in such a way that the expression above does not change.  (Text Example 2.6.2, repeated in this NoteBook below, explains how calculus is used to compute this expression.)

If we know the initial values of  s  and  i, as well as  c, we can first compute the constant

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

and then compute the limiting value of  s by solving

k = s[∞] - ((Log[s[∞]])/c)

If we do not know  c, this invariant quantity can be used to measure it.  If we know the susceptible fraction before and after an epidemic, we can solve the two equations in two unknowns,  c  and  k:

k ≈  s[0] - ((Log[s[0]])/c)

k = s[∞] - ((Log[s[∞]])/c)

In this case it is not hard to show that

c ≈ ((Log[s[0]]-Log[s[∞]])/(s[0] - s[∞]))

References to the Text

Example 2.6.2.  How Calculus finds the Time-Invariant Equation

Without justifying why this works, we compute the ratio

{di/dt}/{ds/dt} = (a s i - b i)/(-a s i) di
{--}
 dt    a s i - b i
---- = -----------
 ds     -(a s i)
{--}
 dt

cancel dt and do some algebra

di=((1/(c s))-1) ds=(1/(c s))ds-ds

In Chapter 13, we will learn that if y = s + k for any constant, then dy/ds = 1 and that if dg/ds = (1/c)(1/s) for any constant c, then g = (1/c) Log[s], the natural logarithm.

These facts of calculus let us antidifferentiate as follows:

di=(1/(c s))ds-ds

∫di=∫(1/(c s))ds-∫ds
∫di=(1/c)∫(1/s)ds-∫ds

i = (1/c)Log[s] - s + k , for some constant k

This is the important time invariant equation for the epidemic model.  We study this is in detail at the end of Chapter 23.  It is needed in the project on the maximum in an epidemic.

Step-by-step Computation

This group of cells shows you how to use Mathematica to approximate the limiting fraction of susceptible people left after an epidemic.

We first calculate the quantity at the start of the epidemic and denote the value by  k.

The contact number for a rubella epidemic is given in the next cell.  The initial fractions of infected and susceptible people are then used to calculate the constant  k  when you press the <Enter> key with the insert cursor in the active cell.  (Do this.)

c=6.8;
s=N[(1/3)];
i=N[(1/30)];
k=i+s-((Log[s])/c)

Since the limiting value of  i[t]  is 0, the limiting value

s = sL = lim s[t]

will satisfy the equation

i + s - Log[s]/c == k

with  i = 0, where  k  is the same constant value you just computed.

A value of  sL  that satisfies this equation will be the limiting proportion of susceptibles.  We cannot solve this equation symbolically, but Mathematica has general "root finding" procedures that find numerical approximations to equations.

First, let's plot the curve i = k - s - Log[s]/c  and look for the graphical solution of the  s-value where  i = 0:

Clear[s,i];
i[s_]:=k-s+((Log[s])/c);
Plot[i[s],{s,0.01`,0.5`},PlotRange->{0,1},AxesLabel->{"s","i"},PlotStyle->{Orange}]

The graph you just computed crosses the  s-axis (where i = 0) somewhere around  0.05  and again around  0.4.  We are looking for the smaller root as the limiting percentage of susceptibles.  (The larger one corresponds to an epidemic that starts with only a tiny fraction of infecteds and around  40% susceptibles.)  Once we have an idea of where the root is, we can use the Mathematica function, "FindRoot" to get an accurate numerical answer.

In FindRoot below we use == instead of = or :=.   Mathematica uses  =  to assign values,  :=  to make definitions and  ==  to denote equations.  FindRoot needs to be provided with a number near the root to start looking for it.  We give it  0.05  by inspecting our graph.

FindRoot[i[s]==0,{s,0.05}]

From this output, we can conclude that after the rubella epidemic has run its course, approximately 3.5% of the population will still be susceptible.

Use the above procedure to find roots with other initial guesses such as  0.4, 0.3, 0.2, 0.1.  Which of these lead to the small root and which lead to the large initial susceptible root?

FindRoot[i[s]==0,{s,0.4}]

Complete Sample Computation

The following is a summary of the limiting susceptible computation:

c=6.8`;
s=N[(1/3)];
i=N[(1/30)];
k=i+s-((Log[s])/c);
Clear[s,i];
i[s_]:=k-s+((Log[s])/c);
Plot[i[s],{s,0.01`,0.5`},PlotStyle->{Thickness[0.01`],Orange}]

initialGuess = ??;
FindRoot[i[s]==0,{s,initialGuess}]

Exercises

Exercise 2.6.1: Math Nerd Fever

At the beginning of semester last year, 96% of the students were not interested in calculus, but word of the chapter on epidemics spread around campus infecting the population with mythical Math Nerd Fever.  Vitims of the disease find that they become permanently interested in calculus, but are only successful in spreading this interest for a brief period.  As a result of this infection, at the end of the academic year 44% of the population was interested in calculus.  On the average, how many new math nerds does each infected student create?

For Math Nerd Fever, we know  s[0] = 0.96 and  s[ Infinty ] = ??.  What can we assume about  i[0]  and  i[∞]?  Use the invariant equation to find the contact number  c.

Exercise 2.6.3: Analytical Solution to s[ Infinity ]

1)      Use the sample computations above to find k = i + s - Log[ s]/c  when
	c=6.8, s[0]=1/3, and i[0]=1/30.

2)      Then find all solutions in s of the equation

k = s - Log[s]/c        Log[s]
k = s - ------
          c

using your k value from part one.  Notice from the graph that there are two roots.  The smaller one is s[∞].