The S-I-R Model: 
A Second Computation with Lists

by
K. D.  Stroyan
University of Iowa

copyright 2007 by Keith Stroyan - All rights reserved.

Notebook Overview

The purpose of this notebook is to get you acquainted with how Mathematica is used to solve problems with recursive computation and "lists."

The main drawback of the First S-I-R NoteBook is that we must enter a computation for each time step.  If we compute hourly, this would take 24 steps each day.  The program below uses an command to calculate a whole list of values in any size step.  The results are stored in a "list" and the list can be printed or graphed.

References to the Text

In Section 2.2  we derive the epidemic model that up-dates the susceptible and infectious population in shorter time steps, dt, say.

S[t+dt] = S[t] - a (S[t]/n) I[t] dt
I[t+dt] = I[t] + ( a (S[t]/n)  - b) I[t] ) dt

Starting from S[0] and I[0], we compute S[dt] and I[dt], then use that to compute S[dt + dt] = S[ 2 dt ] and I[ 2 dt], etc.  This recursive process is a function,

Input  {S[t], I[t]} -> Output {S[t+dt], I[t+dt]}

Let's call this function G:

G[{S[t],I[t]}] = {S[t+dt],I[t+dt]}

Recursive computation means we do the "G" operation over and over,

{S[0],I[0]}
{S[dt],I[dt]} = G[{S[0],I[0]}]
{S[dt+dt],I[dt+dt]} = G[{S[dt],I[dt]}] = G[ G[{S[0],I[0]}] ]
{S[3 dt],I[3 dt]} = G[{S[2 dt],I[2 dt]}] = G[ G[{S[dt],I[dt]}] ] = G[G[G[{S[0],I[0]}]]]
	etc.

Mathematica's NestList[.] command performs exactly this computation, so all we need to do to automate FirstSIR is understand how to define the big function G and apply NestList[.].

NestList Explained

The NestList[.] command computes the list {x0, G[x0], G[G[x0]], ... } the number of times you specify.

Clear[G,x];
NestList[G,x,5]

Enter the following to see Mathematica's own explanation:

??NestList

Some Simple  NestList[.]  Examples

Enter the following computations to see what happens.  In this case, the process function  G  adds  dt  to  t  and  2 dt  to  x.  We start at  t = 0, so get times 0, 0+0.1, 0.1 + 0.1, ... and we start at  x = 1, getting  positions of  1, 1 + 0.2, 1.2 + 0.2, 1.4 + 0.2, ...

G[{t_,x_}]:={t+dt,x+2 dt};
t0=0;
x0=1;
dt=0.1`;
finalTime=2;
solnList=NestList[G,{t0,x0},Floor[(finalTime/dt)]]
ListPlot[solnList,Joined->True]

Enter the following computations to see what happens when we take smaller steps.

G[{t_,x_}]:={t+dt,x+2 dt};
t0=0;
x0=1;
dt=0.01`;
finalTime=2;
solnList=NestList[G,{t0,x0},Floor[(finalTime/dt)]]
ListPlot[solnList,Joined->True]

Practice Exercise

Now you try.  Recursively compute the pairs, {0,1}, {0+1/2, 1/2}, {1/2+1/2,(1/2)*(1/2)},...{t+dt,t*dt}, ... , for  t = 0 to 10 in steps of 1/2.

A Step-by-step Computation for S-I-R

The various steps in computing an epidemic are given in the separate numbered cells below.  They are combined later in the three month computation cell.

1)  The Disease-Specific Parameters

First we define our rate constants a and b, and give the initial numbers of susceptibles, infectives and recovered as before.  Initially, time t0 = 0.

Type the correct expressions where we have put ?? and <Enter> this computation.

c=N[6.8];
b=N[(1/11)];
a=??;
Print["Contact number = ",c,", Infectious period = ",((" days.")/b)];
Print["Daily rate of contacts, a = ",a];
n=30000;
Print["Total population, n = ",n," so  rNr = ",n," - sNr - iNr"]

2)  The Program Parameters, dt and tfinal

We will write a general purpose program that can work with any time step size,  dt, and compute the values of  sNr, iNr  in steps of  dt  up to time  tfinal.  The next cell assigns the values of these program parameters.  (<Enter> this computation now.)

dt=N[(1/2)]
finalTime=4

3)  The Initial Values

Mathematica is designed to work well with lists, so we will store our computations in a list.  Each entry in the list will be a list consisting of the independent variable time,  t, followed by the dependent variable, sNr[t], iNr[t]  or  rNr[t].  Our main loop below appends each new set of values to these lists.  We head the list with the initial values assigned a few cells above:

(<Enter> this list assignment now.)

t0=0.0;
sNr0=10000.0;
iNr0=1000.0;
initialValues={t0,sNr0,iNr0}

4)  The Main Recursion

To update the values of  sNr, iNr, and rNr, we will use the Do loop below.	The syntax in general is in this closed cell with simple examples:  (Enter the next cell to do the main computation with the previously assigned parameters.)

sRate[sNr_,iNr_]:=-((a sNr iNr)/n);
iRate[sNr_,iNr_]:=(((a sNr)/n)-b) iNr;
G[{t_,sNr_,iNr_}]:={t+dt,sNr+sRate[sNr,iNr] dt,iNr+iRate[sNr,iNr] dt}
solnList=NestList[G,initialValues,Floor[(finalTime/dt)]]

5)  Numerical Output

The recursive process function G  produces a list of triples and we would like a list of pairs  {time, susceptibles} and a list of pairs {times, infectives}.  This is technical, so look at your computed list of susceptible data by entering the next few cells and watching how the lists are minipulated

Transpose[solnList]

Take[Transpose[solnList],{1}]

Take[Transpose[solnList],{2}]

Take[Transpose[solnList],-1]

Flatten[Take[Transpose[solnList],-1],1]

suscs=Flatten[Take[Transpose[solnList],{2}],1];
times=Flatten[Take[Transpose[solnList],{1}],1];
infes=Flatten[Take[Transpose[solnList],-1],1];
sPairs=Transpose[{times,suscs}]
iPairs=Transpose[{times,infes}]

Print[TableForm[sPairs]]

6)  Graphical Output

Use the lists sData, iData and rData to draw graphs of the solutions to the S-I-R model by Entering the following:

ListPlot[sPairs,Joined->True,AxesLabel->{time,susceptibles},PlotRange->{0,n}]
ListPlot[iPairs,Joined->True,AxesLabel->{time,infected},PlotRange->{0,n}]

Three Months of Rubella in One Computation:

The next cell combines all the computations you did above for rubella, but changes the duration to 90 days.  <Enter> the computation.

c=N[6.8`];b=N[(1/11)];
a=b c;
n=30000;Print["Spread of rubella in a population of ",n];
Print["Contact number, c = ",c,", Infectious period 1/b =",(1/b)];
Print["Daily contact rate, a =",a];
t0=0.`;sNr0=10000.`;
iNr0=1000;
initialValues={t0,sNr0,iNr0};
dt=N[(1/2)];finalTime=N[90];
steps=Floor[(finalTime/dt)];
Print["Increments of dt = ",dt," for ",steps," recursive steps."];
sRate[sNr_,iNr_]:=-((a sNr iNr)/n);
iRate[sNr_,iNr_]:=(((a sNr)/n)-b) iNr;
G[{t_,sNr_,iNr_}]:={t+dt,sNr+sRate[sNr,iNr] dt,iNr+dt iRate[sNr,iNr]};
solnList=NestList[G,initialValues,steps];
suscs=Flatten[Take[Transpose[solnList],{2}],1];
times=Flatten[Take[Transpose[solnList],{1}],1];
infes=Flatten[Take[Transpose[solnList],-1],1];
sPairs=Transpose[{times,suscs}];
iPairs=Transpose[{times,infes}];
ListPlot[sPairs,Joined->True,AxesLabel->{time,susceptibles}]
ListPlot[iPairs,Joined->True,AxesLabel->{time,infected}]
Print["The last  5  {time , susceptible}  pairs:"];
Print[Take[sPairs,-5]];

Exercises

Exercise 2.2.3

1)   Complete a S-I-R table for four days  with a time step of 1/2 day.  Graph S, I, and R vs  t.  This can be done by entering the above cells  in the step-by-step explanation for a time step of 1/2 and a final time of four days.

2)  Compute and graph a 3 month rubella epidemic using the combined computation.