Calculation of the Average Duration of an Illness and Computation of the Reproduction Number in the SIR Model....

by Douglas E. Lewit, mathematics student at Northeastern Illinois University and secondary education major at National-Louis University

I prepared this Maple worksheet as part of a presentation to Professor Mubayi's lab group at Northeastern Illinois University.  Every member of the research group explores a different aspect of how mathematics is used to study public health, which is Professor Mubayi's primary speciality or research focus.  During this presentation, I explore two different SIR models.  I begin with a simpler SIR model that has only two parameters, and β.  α, in this context, is essentiatlly an infection rate, while β, within this context, refers to the rate of recovery.  Of great importance in the study of mathematical epidemiology is a value called the reproduction number, usually denoted as The exact numerical value of  varies from model to model and problem to problem, but in general the reproduction number reveals the "spreadability" of the illness being studied.  If  is less than 1, the illness will die out before reaching epidemic proportions.  If  is greater than 1, then government officials should be alerted as soon as possible because the illness has the potential to spread to epidemic proportions.

For my research I studied two articles.

The title of the first article that I used for my research is S-I-R Model of Epidemics, Part 1, Basic Model and Examples, Revised September 22, 2005.  The article is posted on the internet (apparently without an author's name) and can be found at http://www.me.rochester.edu/courses/ME406/webexamp6/sir1.pdf.  Although the file is a .pdf, much of the mathematical work for the article was done using Wolfram Mathematica.  As I was reading through the article, I had to translate some of the Mathematica code into Maple code.    

The title of the second article is Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, dated June 27th, 2001, written and researched by P. van den Driessche from the Department of Mathematics and Statistics, University of Victoria, Victoria, BC, Canada, V8W 3P4, and James Watmough from the Department of Mathematics and Statistics, University of New Brunswick, Fredericton, NB, Canada, E3B 5A3.

First off, thanks for attending my presentation on Saturday morning.  I realize that most people just want to sleep late on Saturday mornings!  So thanks again for making the trip.  Okay, let's get started.  To begin with, I want to use the following system of differential equations to show why the average duration of an illness is

where S stands for the number of susceptible individuals in the population, I is the number of infected individuals, and finally R is the number of removed (i.e. dead) or recovered individuals.  The underlying assumption here is that the population is constant.  (The sick people aren't having kids during their period of sickness.  Makes sense to me.  When you're sick, do you really want to take care of children and screaming babies?)

Consider the second differential equation.  The underlying assumption here is that α is usually a lot smaller than β, and so for the sake of generalizing and simplifying the problem the author of the article that I read temporarily discarded the parameter, α.  Thus the second differential equation becomes, with some initial condition, I(0)=_C1, where _C1 is some constant greater than 0.

Therefore, the average value of I(t) on the interval (0, ∞) is approximately but remember that our conclusion is based upon the assumption that α is a lot smaller than β.  As α gets larger, well, our assumption (and hence our result) falls to pieces, okay?  Now let's turn our attention to finding the reproduction number for a particular case of an infectious model.

Just checking to see whether or not I can use the variable 'e'.  It's possible, but now I no longer have access to the variable 'e' for entering numbers in scientific notation.  For example, 0.01 can be entered in Maple as 1e-2.  Now that I'm using the variable 'e' for another purpose, I no longer have that option.

Unassigning a variable is simple!  Try it for yourself.

In the article, the functions that constitute the right-hand sides of the differential equations are split into two parts: 1) infection coming in, and 2) infection going out.  In the second case, you have to factor out a negative sign in order to implement the author's method.  Let's see how to do this with Maple.

A couple comments are in order here.  First of all, S/N is assumed to be 1, but at the disease-free equilibrium point (DFE) the number of susceptibles is assumed to be the same as N, the total population.  At the DFE the treatment group has a population of 0, hence T=0 in the matrix above.  Simple, right?  :-)   Well maybe not, but if you think about it for a couple of minutes it starts to make sense.  Also, when using Maple's Jacobian command, there are 2 required arguments.  (There's an optional third argument as well, but let's save that one for a later date.)  The first argument must be a vector while the second argument must be a list.  Data types are extremely important in Maple.  If you don't use the correct data types, your Maple commands won't work, plain and simple.  Be careful!  Maple's built-in Help menu is very helpful.  Spend some time looking it over before you begin your next Maple project.

This agrees with the computation in the article written by P. van den Driessche and James Watmough.  (Look at the bottom of page 7 of their article to see the final result of their calculation.  The only difference is that they factored a part of the denominator of their final expression.  But my Maple result really does agree with their result.)