First Exercise in Epidemiology 606

Department of Epidemiology

Professor Jim Koopman MD MPH

This document serves 3 purposes. It presents a brief introduction to compartmental modeling at a level appropriate for Epidemiology 606 students. It provides a step by step exercise that teaches you many skills in using Stella IIÔ for modeling. This lets you skip doing most of the tutorials that come with Stella IIÔ . Finally it presents three epidemiological models that are cornerstones upon which many epidemiological models are constructed.

For a more thorough introduction to modeling epidemiological processes, you can review Chapters 1, 2, & 3, in the text for the Epid 802 course that is found on the web at http://www.sph.umich.edu/~jkoopman/802Web/Course.htm. But you should only do so if you have a special interest in infectious disease modeling or if you think you would like to do doctoral work with me. Otherwise, the material in this chapter is self-contained. You should be able to do the homework and the exams perfectly well purely on the basis of material presented here.

The computer models of transmission systems we will study are called continuous compartmental models. Continuous refers to the assumption that the entities modeled are quantitatively continuous rather than discrete. The continuous entities we will model are segments of populations. Populations are made up of discrete individuals that cannot be subdivided. But the type of model we use assumes that populations can be broken up into an infinite number of pieces. Making this assumption limits the aspects of a transmission system that can be modeled. But continuous models have provided considerable insight into epidemiological processes despite these limitations.

There are several hidden assumptions that come along with the continuous population assumption. One corollary of the continuous assumption is that all contacts are made instantaneously in time and have no duration. Once a contact is made, there is no ongoing aspect of that contact that increases the likelihood of contacting again the person originally contacted. Actually, there is no contact between persons in models of continuous populations. Continuous population models seek to abstract some essential consequences for a population of contact events between individuals without ever including individuals in the model. By assuming continuity, we avoid several complications and complexities that would occur if we had individuals in our models. This assumption also enables our models to be analyzed using the theory of calculus. The models we will construct and analyze in this course are actually differential equation models. But don't let that scare you off. We will formulate and solve these differential equations using a point and click approach on the computer.

Continuous compartmental models are not the appropriate models for addressing every question regarding the dynamics of epidemiological processes. The continuity of the population unit is a bad assumption if we want to address questions about what happens in small population units, such as families. It is also an inconvenient assumption if we want to model sexual transmission between concurrent partnerships. Continuous compartmental models work best when the questions posed relate to what happens at the broad population level. Our purposes in this chapter are to get insights into the behavior of infections in large populations. Compartmental models are good for such purposes because they are conceptually simple and can be readily solved using numerical methods.

When we say we "solve" a continuous compartmental model of a population process, we mean that the variable values as time progresses are determined by numeric techniques. In Epid 606, we do not learn much about the numerical methods used. We just formulate our models visually by defining compartments and what determines the rate of flow between compartments. Then we decide what compartment values or derived variable values we want to observe over time, put them in a table or a graph, and hit the run button. We let the software worry about the solving. Of course, if you understand how the computer is working, you might avoid a few pitfalls. But in this course we do expect you to learn such details.

Some people call the process just described "simulation". It is more accurate to say that we are numerically solving differential equations. The term simulation is more commonly reserved for models where individual events are simulated. Most simulations are Monte Carlo computer models. In Monte Carlo models chance determines the events that take place in the computer. The models we examine, however, are mathematical models. We use the computer to solve them and by so doing we find the unique values defined by the model structure its parameters, and its initial values. Chance plays no role in the compartmental models we will examine.

Epidemiology 606 students do not need to learn about how to formulate continuous models of transmission systems in terms of differential equations. This class is not intended to prepare you to read the modeling literature. It just provides an introduction that might stimulate you to go on and learn some more modeling.

If you want to read the modeling literature, you will want to learn how to formulate Stella IIÔ models in terms of differential equations. Even more helpful sometimes is being able to construct a Stella IIÔ model from a set of differential equations that appear in an article or book you are reading. Doing that insures a much better understanding of what the article is trying to accomplish. You are quite capable of learning to go back and forth between Stella IIÔ models and ordinary differential equations. There is no real trick to it. But we do not have time to cover that formally in this class. If you decide to acquire this skill, the first four chapters in the 802 course provide a reasonable starting point. They deal with ordinary differential equations. The Anderson and May Text on Infectious Diseases of Humans: Dynamics and Control, starts off with partial differential equations. To keep things simple, we will not address the use of partial differential equations or Stella IIÔ models employing them in this course. Partial differential equations are especially useful for addressing issues related to the age distribution of infection.

Although you will not learn to formulate differential equations, you should understand that the computer models we will construct are intended to approximate smooth, continuous, differential equations. You need to understand this because you must realize that when the graphs you generate are not generating smooth lines, the computer is not accomplishing its intended purpose of solving the differential equations.

Continuous differential equations have an infinitely divisible time dimension. When solved, they always generate smooth curves. At times a line generated by a continuous differential equation may turn around very quickly and look like it is making a sharp angle, but if you blow up the time scale and examine the curve much more finely, you will see that it is truly continuous.

The Stella IIÔ models we construct don't operate continuously overtime like theoretical differential equations do. Stella IIÔ models approximate the performance of differential equations by taking time steps so small that it would make no practical difference if we made them any smaller. The output of the Stella IIÔ models we will examine consist of curves of variable values over time. Using very small time steps makes those curves smooth or continuous. If ever you see sharp spikes or erratic behavior in the curves you generate in the exercises you are asked to complete in this course, you should reduce your time step or "DT". If you are taking a time step of 0.25 years, you may have to reduce it to 0.025 years. If you are at the Stella IIÔ limit for the smallest time step, you have to change time units. For example, you may have to switch from years to days. Getting to short enough steps will make Stella IIÔ output correspond to the values of differential equations.

Before beginning to construct some simple models using the Stella IIÔ software, let us review some elements of transmission models. The list here is not exhaustive.

One broad way to divide the elements of compartmental models is into variables and parameters. Variables are the model elements that the process modeled causes to vary over time. Our variables will be compartments of population. We will represent these compartments by Stella IIÔ "stocks" which are visually presented as rectangles. Parameters are fixed model elements whose values are determined by us before we run the model. We will try to always put our parameters into Stella IIÔ "converters" which are represented visually as circles.

In compartmental models, parameters are model elements that, along with Stocks (compartmental variable levels) determine flows. Parameter values do not depend upon variable values. Variable values do depend upon parameter values.

Using Stella IIÔ it is possible to put parameter values directly into flow regulators without ever having to name the parameters. That is not good practice. By putting all parameters into Stella IIÔ converters and then connecting them to flow regulators with connectors, we more clearly identify our model elements.

The models we will construct are called compartmental models. For the very restricted type of models we will examine, we can think of Stella IIÔ stocks as model compartments. One definition of a compartmental model is a model where the units of the entities modeled never change. The entities modeled can change their status. Indeed that is the crux of the modeling process. We model a change in status as a flow from one compartment to another. But the units of the entities we model never change as they flow.

Many physics and chemistry models are not compartmental models. In physics models, the entities modeled are often transformed in a way that implies unit changes. An example is force getting transformed into acceleration. In Chemistry models, one chemical gets transformed into another chemical. But the human populations we will model will always remain human populations. Stella IIÔ can be used for modes with unit transformations like the Chemistry and Physics models. In Epidemiology 606, however, we do not learn how to construct any model type other than compartmental models.

Population compartments are the "state" variables in our models. All the different values that they can conceivably take on under the conditions we are modeling are called the "state space" of our model. The value of these variables or compartments changes as population flows in and out of them when we run the model over time. In the second model we will examine, we will model birth processes as population flowing into a compartment from a cloud. We will model death processes as population flowing out of a compartment into a cloud. These clouds are connected to the compartments by flow pipes and must be considered an integral part of the compartmental system having the same units as the rest of the system.

We will model all events happening to population as flows from one compartment to another. In other words, the events we will model are changes in status of individuals. In our first two models, we divide everyone (or more properly all the population since our models have no individuals) into three population compartments: Susceptible (S), Infectious (I), and Immune or Removed or Recovered (R). The events are getting infected and recovering from infection.

The art of modeling involves capturing the essence of something while discarding the non-essential. Deciding what are the essential aspects that need to be modeled and determining whether the essence of something has been captured can be viewed as an exercise of scientific judgement. When you construct a model, you are essentially constructing a theory about how a system works. Whether your theory has captured the essence of a system will depend upon how you intend to use your theory. As indicated in the chapter on Transmission System Analysis, models that are quite adequate for some purposes may lack essential elements to accomplish other purposes. We will always want to be complicating and exploring variants of the models we construct in order to better understand the transmission systems that we model. But by beginning extremely simply, we create a base of understanding that enables us to better understand model behavior. Such understanding is needed to avoid mistakes in model construction when we wish to address models capturing more details or more processes.

We begin by dividing the entire population into three classes as indicated above. We disregard whether individuals have symptoms or not in making this classification. We disregard differences between individuals in the degree of their susceptibility. In the real world, exposure to a few organisms might be enough to infect some individuals while other individuals might be able to resist infection by small doses of organisms but still be susceptible to infection when exposed to millions or billions of organisms. But our initial modeling purpose is just to understand how the fact of infection transmission and the acquisition of immunity affect epidemic and endemic infection dynamics. So for now we eliminate all model elements except those that relate to infection transmission and immunity. Thus, we model only one category of susceptible individuals. We do not distinguish individuals by age or by anything else. For many modeling purposes, we will need to distinguish people by age. But we do not have time to get to models with age structure in this course.

Our modeling of the infection process is a gross simplification. In the models we construct, population that flows from the susceptible to the infectious compartment will immediately reach a constant level of contagiousness and stay at that level of contagiousness until they recover. Real infections have an incubation period. Once they start excreting organisms, the level of excretion at first rises exponentially, peaks, and then falls off. Some people have a faster course with less excretion of infectious agent. Some have a long slow infection. Some individuals reach high levels of contagiousness quickly. We disregard all such detail. Everyone is the same in our model and everyone moves from completely susceptible to the same level of contagiousness. Then they move just as suddenly to a state of complete immunity where they stay forever more without the slightest chance of being infected again.

No infections behave exactly like this. Even measles, which we once thought generated lifelong complete immunity, does not generate complete immunity. We now know that the appearance of lifelong immunity was sustained by a degree of susceptibility to infection that did not prevent the initial stages of infection. Rather the role of immunity is to bring infection under control before a stage with symptoms is reached. For second and subsequent measles infections in the pre-vaccine era, immunity controlled infection in the incubation stage before a it progressed to a contagious state. But initial viral replication can often proceed enough to give a boost to immunity. Such boosting gives the appearance of lifelong immunity. But in our model, there will be no immunity boosting infections and still immunity will be permanent.

Our models make another extreme simplification with regard to the time of recovery from infection. We assume that the time to recover is independent of the time one entered the infectious state. In the simple model we will construct with one single infectious state compartment, individuals (more correctly population) who have been infectious for10 years will recover on average at the same time as individuals who have just arrived in the infectious compartment. Such independence of flow from the time of arrival is a characteristic of all single compartments in compartmental models. Within a compartment, there is no distinction of individuals by the time since they arrived in the compartment. Within a compartment, it is assumed that there is complete and instantaneous mixing so that all segments in a compartment have the same chance of leaving at any instant. Stella IIÔ has different classes of stocks. Only the Stella IIÔ "reservoir" stocks that have this thorough and instantaneous mixing characteristic. We do not have time to learn the use of the "conveyor" stock that does keep track of when individuals arrived and puts them out some fixed time later.

A constant flow rate out of a stock generates a distribution that has a particular shape of times that individuals stay in the stock. This distribution is called a negative exponential distribution. I will give here a brief overview of some material in chapter 2 dealing with constant flow rates and negative exponential distributions that is sufficient for the Epidemiology 606 course.

Start Stella IIÔ and get to the constructor page. It should look like the following:

Left click on the "Stock" tool (the square box on the upper left) with your mouse pointer and deposit it on your screen with another left click. Give it the label "I". Next press the "Alt" button and left click deposit another stock to the right of "I". Label it "R". Then click on the flow tool which is just to the right of the stock tool and is colored in the above picture. Now left click on "I" and hold down while you drag the pointer to the middle of "R". R will light up when you have connected the flow to it. Let up on the mouse to complete the connection. Lable the flow "recovery".

Now click on the globe on the left hand border of the Stella IIÔ screen. Question marks will appear where you need to fill in either starting values for the stocks, or flows. Your Stella IIÔ diagram should now look like the following:

How much flow there is from "I" to "R" will depend upon how much "I" there is to flow. We indicate this by connecting "I" to the flow regulator "recovery" with a Stella IIÔ connector. Do this by clicking and releasing on the connector, then clicking down on the "I" compartment and while holding the mouse down, dragging it to the recovery. When this shadows over, let up on the mouse. How much flow there is will also depend upon the recovery rate. The recovery rate is a parameter in our simple model of linear flow from I to R. Let us put a converter in our model which we will label as "recovery rate" and connect it to the recovery flow regulator. Your diagram should now look like the following.

Now we have the outlines of a simple linear flow model from one compartment to another. The next task in building our model is to establish the amount of "I" we start with, the amount of "R" we start with, the value of the recovery rate, and the mathematical relationships in the flow regulator between the various model entities that influence the flow. Let us treat our population as whole so that the number value in any compartment represents the fraction of the total population that is in that compartment. If the total population fluctuated in our model, we couldn't do that. If we had the potential for births not to equal deaths as they do in this model (0 births and 0 deaths), then even if we started out with a total of all population compartments equalling one, that total would change and to get the fraction of the total that is in any compartment we would have to derive a variable in a converter by dividing the compartment value by the sum of all compartment values.

This makes us pause to think what is the meaning of the numbers in our compartments. The numbers never really correspond to a number of individuals because the numbers include fractions of individuals in a way that is really based on the assumption that there are an infinite number of individuals in each compartment. This means that we cannot use continuous compartmental models to determine what is the critical population size to sustain circulation of infectious agents. The nature of the model is just not appropriate for that.

But our present purpose is merely to understand the nature of the distribution of infection times in a compartmental model where individuals flow out of the infectious state at a constant rate that is independent of the time when they flowed into that state. Since our total population will stay constant, we can start with all of a population of size one in the "I" state and none in the "R" state. Double click on these compartments to put these initial population values in them as follows.

Now we have the issue of putting in a value for the rate of recovery. This rate will be per person per unit of time. Using Stella IIÔ we are not forced to choose a unit of time. We could leave it undefined and abstract. But say we are thinking of a real disease that has an average duration of about ten days. To generate an average duration of 10 days, individuals need to flow out of "I" at the rate of 1/10. In epidemiology 606, we will not attempt to give you a deep understanding of why this is. Chapter 2 of Epid 802 attempts to do this. But you can get along for this course just by accepting this fact. A constant flow rate out of a single compartment will generate an average duration of stay in the compartment equal to one divided by the rate. We want to see something about the distribution of times that go into making up this average. So click on the recovery rate converter and put in the value of 1/10 (or 0.1 if you prefer). Now to insure that our time units are labeled as days on any graphs that we wish to generate, click on the "Run" menu at the top of the Stella IIÔ window. Hold down the mouse and go to "time specs". When it darkens, release and you will get a window that looks like the following:

Select days from the "Unit of time:" column. Note that we will be simulating from 0 to 12 units of time in time steps (DT) of 0.25. We will very often want to change these values. But we do not need to change them now. Note also that the integration method used is "Euler's Method". Because we are only modeling continuous differential equations, the Runge-Kutta method 2 or the even more precise method 4 might do us better. But for now leave the setting on the Euler's Method. Click the DT as fraction box. You will see that you can specify the time step either as a decimal or a fraction. Some decimals just won't work so it is always better to use a fraction.

Our last task is to formulate the flow. This is where we have to put in mathematical formulas. For the most part, we will just be using very straightforward formulas in this course. The flow has to be proportional to the number of people in "I" available to flow and to the constrant rate that we established. Thus our flow formula is just "I" times the "recovery rate". Double click on the flow regulator "recovery" and you will see a window that looks like the following.

Stella IIÔ forces you to include all of the units you have connected to the flow regulator and does not allow you to use any other values. In the "Required Inputs" box, click on "I", then since we want to multiply our two inputs, click on the "*" on the keypad to the right of the "Required Inputs" or just enter the "*" from your own keyboard. Finally, click on the "recovery rate". Our formula is complete. It is good practice to explain why we formulated a flow in a particular manner so click on the "Document" box and write in an explanation as follows.

Click on "Hide Document" and go back to the flow window. Note that our flow is only in one direction. One never becomes infectious again after recovery. If we wanted a two-way flow, we could select "BI-FLOW" instead of "UNI-FLOW". The back flow would have to be entered as a negative value after the positive forward flow. But we are just using a one directional flow. Also on this window, note that because we are constructing compartmental models, we do not check the unit conversion box. Click OK.

Now we are ready to run our model. We need to see our output, however. In this course we will adopt the practice of always generating our output both in graphical and in numeric form. From the main Stella IIÔ window, select the graph icon (the first icon in the third set of icons along the top) and deposit it on your diagram. A graph appears. Note that it has units of days because we selected those units earlier. Double click anywhere on that graph. You get a window that looks like the following.

In the "Allowable" window, select "I" and then click on the arrows that move it to the "Selected" window. Now select both the stock "R" and flow "recovery". Move them to the "Selected" window. That gives you a window that looks like the following:

Click OK. This gets you back to the graph window. Go to the run menu and select run. This generates the following graph.

Note that Stella IIÔ automatically scaled each of the variables we plotted. Not having everything on the same scale can be confusing. To put everything on the same scale, double click anywhere on the graph and selected all of the entries in the "Selected" window and then click on the top arrow of the selected entries at the right of the "Selected" window. This gives you the following:

We can set all variables on the same scale of 0 to 1. In this case you don't have to enter the scale you want to set because that is the default. Just press set and then OK. Now your graph looks like the following:

Hand in 1.1:

Hand in a graph like that the one above.

If you run the model for 30 days you get the following:

Curve three is just one tenth the value of curve 1. That is because the flow rate out of "I" which is represented in curve 1, is 0.1. Curve three represents the probability density distribution for time spent in the infectious state. There is more probability that population will flow out of 1 at time zero than at any other time. The density decreases in a negative exponential fashion. The total area under curve 3 equals 1. If you take the average time by summing (actually integrating because this a continuous curve) the value of curve 3 times the corresponding time (and dividing by the total which is 1), you will find that the average time spent is 10.

We said, however, that besides graphic output, we wanted tabular output for our model. Go back to the model diagram in Stella IIÔ . Click on the table icon which is just to the left of the graph icon. Deposit it on the diagram. A table comes up. Double click anywhere on the table. You get the following:

Select the I, R, and recovery rate and put them in the table just as you did for the graph. Instead of reporting every interval, click on the "Every DT" box. The default report interval is 1. You could increase it but 1 is fine for our purposes. Click OK and run the model. Your table should look like the following.

Now go to the run menu, select time specs, and set the DT to 1/8. Open the table and run the model again. You get no difference from the DT = ¼ so your DT setting and integration method are fine. Now, however, set the DT = 1 and run the model again. Now your table should look like the following.

Since choosing a shorter DT changes your values, the DT of 1 is too long using Euler's Method. Let us try using a Runge-Kutta 4, however. Go to the run menu, select time specs, and then select Runge-Kutta 4. Running at a DT=1 now gives you the same values you got with a DT of both ¼ and 1/8. You should always check to see that your integration method and DT are adequate by checking to see that cutting the DT in half does not change your values. For big continuous compartmental models, you will find that a longer DT with the Runge-Kutta 4 will decrease the calculation time.

Handin 1.2

Generate a table like the one above. Explain the calculations that lead "I" to decrease in the pattern that it does and R to increase. This should include an explanation as to where the recovery values come from.

It seems unlikely that the most frequent time spent would be zero. We can change the distribution of the time spent to a more realistic distribution by stringing two or more compartments together. Let us do this.

Generate the following model.

Put one population unit in I 0 and zero in I 1 and R2. Make the duration in each stage 5 days. Figure out the "Stage recovery rate" that you need to use in order to do this. Run the model for 30 days and turn in the following graph to demonstrate that you were able to do this.

Handin 1.3

Answer the following questions: Why does the value of the "Final recovery" flow give the distribution of times spent in the infectious state? What is the most common duration of infection when infection has two sequential recovery stages and each has on average a duration of 5 days? What is the average duration of infectiousness in the two stage infection model just examined?

Until recently, there was a rather cavalier attitude toward formulating contact and transmission processes in transmission models. Contact and transmission were treated as a single process. To establish a science of transmission systems, careful formulation of separate contact and transmission processes will be required. Contact processes need to relate to measurable parameters so that the appropriateness of transmission system models can be checked against real world data. There are two approaches to relating contact processes to measurable events. First, the contact process formulation could include real events that are measurable, such as touching someone, being in the same room with them, sharing toothbrushes, having sex, etc. Second, the formulation may include secondary aspects of the contact processes, such as the location of contact. I believe that there is much promise in this second approach. Formulating contact in a manner that better relates to measurable events is an advanced topic at the forefront of transmission system analysis. We cannot go into detail in Epid 606.

We will formulate contact processes to answer basic questions about transmission dynamics. We will find it possible to gain some insight into transmission dynamics using contact only general concepts regarding the type of contact. The formulation we will use could be of airborne contact, enteric contact, foodborne contact, or sexual contact. We will not, however, use the most common formulation of contact in introductory transmission system models.

The most common formulation of contact and transmission a single parameter equal to the rate at which every pair of individuals in a population makes effective contact capable of transmitting infection. In this formulation, if there are 100 infected individuals and 100 susceptible individuals, there will be 100 times 100 equals 10,000 pairs of individuals who could make a contact where infection would be transmitted. If of these 10,000 pairs of individuals makes effective contact at the rate of 1/1000 per unit time, then new infections would be generated at the rate of 10 per day. If "S" equals the number of Susceptible individuals, if "I" is the number of infectious individuals, and if b is the rate of contact between each pair of individuals that can potentially exist in the population, then the flow of new infections equals SbI.

That formulation breaks down as more structure is added to the populations modeled. Since we will not be adding much structure to our initial models, we could use that formulation. The models we will construct use two parameters instead of one. For that reason, they would be considered by mathematicians to be less virtuous. But I feel it is better to use a formulation with more potential to be elaborated in ways that can relate to data.

In our approach using two parameters to formulate the contact process, we first define a generic contact parameter. The second parameter is then the fraction of that type of contact that will result in transmission when the contact is between a susceptible and an infectious person. We assume that there is only one type of contact, that it is uniform, and that between any pair of individuals, it is symmetric. Air space contact is most likely to be symmetric. Contact where one person leaves contamination on a surface or a food and the other person picks this up is less likely to be symmetric. But even in that case, a symmetric formulation of contact may appropriate, as long as we are not particularly interested in the consequences of the asymmetry.

This approach of using two parameters to formulate contact processes has been thought of mainly as a way of formulating sexual contact processes. To advance a science of transmission systems, I feel we should always define contact in some measurable way. Only then can we estimate transmission probabilities given that contact. I therefore find this two-parameter formulation to be appropriate for all modes of transmission except those mediated by vectors with fluctuating densities or where growth of the organism occurs in the environment. Even for these latter two modes of transmission, however, we can gain important insights into the transmission dynamics with the formulation that we will present.

We label the rate of contact that individuals in the population make with all other individuals as "c". We label the probability of transmission given such a contact as b. We assume that all individuals mix evenly without regard to their infection status. Then the rate of infection experienced by susceptible individuals will be proportional to the number of infectious individuals, the rate at which they make contact, and the probability of transmission given each contact. This rate will not be constant because the number of infected individuals will fluctuate. In the epidemiological literature this rate of infection is called the "force" of infection. In many differential equation models, you will find it labeled as l.

l = cbI.

An important difference between the one parameter and two parameter formulations is what happens when the total size of the population grows. Consider the situation where the number of susceptible individuals is increased but the number of infected individuals stays constant. In the one-parameter formulation, the flow of new infections will increase in proportion to the increasing number of susceptible individuals. In the two-parameter formulation, the flow of newly infected individuals will not increase. Since we will only be using the two-parameter formulation, we will not dwell upon it. In reading the infectious disease modeling literature, however, you will encounter the one parameter formulation very frequently. You may encounter discussions about how that formulation relates to the density of the population being modeled. I find all of that discussion to be misdirected and poorly thought out. With the two parameter formulation, we can model the "c" as a function of density in a meaningful way while in the one parameter formulation we have no such option.

One of the most interesting aspects of infectious diseases is the immune control of infection. The immune system is incredibly complex, effective, resilient, and yet capable of being defeated by resourceful infectious agents in so many ways. We could make our models of immune processes as complex as immunity itself. Some reasonable amount of detail in the immune response may be needed to understand the effects of vaccination and how different strategies of vaccine design might change the effects that vaccination programs have on transmission dynamics. Transmission system models, however, have modeled immune response effects on the infection processes quite simply. Later we will examine a very simple vaccine effect model where the only effects of immunity are to lower transmission probabilities to vaccinated individuals and to shorten the duration and/or contagiousness of infection in individuals with vaccine induced immunity who become infected.

In the simple models we will initially model of transmission systems, however, our modeling of the immune response will be even simpler still. We will assume that recover from infection generates a complete and everlasting state of immunity so that recovered individuals can never be infected again. No infection has these characteristics. But it is worthwhile for us to examine model transmission systems with these characteristics because these model systems begin to give us insights as to how the effects of immunity become manifest at the population level.

To keep things simple in Epid 606, we will only examine transmission system models where the total population stays the same. To insure that happens, infection cannot increase death rates. There are three common ways to keep the population constant. The first is just not to include births or deaths in the model. We will take this approach in the first transmission system we model. The second approach we will take is to make birth rates equal to death rates. Everyone in our population will be subject to the same death rate. People who just came into the population will die at the same rate as people at age 20 or age 100. But since we will mix everyone in homogeneous compartments without keeping their age, we will not be able to determine the age of population that dies.

Our birth processes will be just as brutally simple. Births will not be generated by sex. Every person in the population will generate births at the same rate.

The third approach to balancing birth and death processes is to model the dependence of those processes on population size in a way that always returns the population to the same equilibrium level when something knocks it off that level. This approach is beyond Epid 606 but can be reviewed in Chapter 3 of Epid 802.

John Holland in his book, Hidden Order, demonstrates how adaptive systems such as those we deal with in the biology of infectious agents will inevitably increase their diversity as they evolve. Transmission systems have interacting adaptive agents in the forms of humans (or other animals), and infectious agents. The humans will of course evolve their diversity on a slower time scale than the infectious agents. HLA diversity that has evolved over many millennia is one type of diversity in the host that is well recognized. Antigenic diversity and physiological diversity will both be important for the infectious agents. Later we will examine how transmission systems are affected by diversity in the susceptibility of agents to chemotherapeutics such as antibiotics. In this chapter and the next four, however, we will assume that the infectious agent in the transmission system we are modeling is homogeneous.

Our new molecular abilities to characterize agent diversity offer opportunities to generate data that could serve as a foundation for a science of transmission systems. But in this course, we must deal with the basics before we get into such exciting topics.

We will address very bare bones transmission system models that either generate epidemics or endemic infection levels. We will examine epidemics first in models that have no births or deaths and in which immunity will be complete and everlasting. Then we will examine two different kinds of models that generate endemic levels of infection. We will first examine how adding births and deaths to an infection just like the one we examined for epidemic behavior can cause the emergence of a constant endemic level of infection. Next, we will examine transmission systems where the immune individuals become susceptible again with the passage of time.

The questions we want to address by examining a transmission system that generates epidemics are:

1. What brings an epidemic to an end?
1. Is it the exhaustion of susceptible individuals?
2. Does something have to change with regard to transmission system parameters in order to leave susceptible individuals who have mixed with infected individuals at the same rate as everyone else uninfected at the end of an epidemic?
3. Is there something in the epidemic process that drives the number of infected individuals down to zero even while there are many susceptible individuals left?
2. What determines how fast an epidemic spreads and what fraction of the population will be infected over the course of an epidemic?
3. How many independent parameters are determining the size and shape of an epidemic in a transmission system model with the following characteristics:
1. The population is homogenous
2. The population mixes uniformly
3. Infection is uniform
4. Recovery from infection is not dependent upon how long one has been infected
5. Recovery from infection results in perfect and everlasting immunity

Note that in answering these questions, we are not answering questions about real world transmission systems. But our models capture some aspects of real world transmission systems so that the insights we gain should help us gain insights into real transmission systems.

The first transmission system we will construct will put all individuals into the categories S, I, and R as discussed previously. It will have a contact parameter and a separate transmission probability parameter. It will have one additional parameter that relates to how long infection lasts before it induces immunity.

Construct a model using Stella IIÔ with the following diagram:

Put 0.99999 population units initially in S. Put 0.00001 population units initially in I, and put no population in R. Formulate the "new infection" flow as S*c*B*I. Note that the force of infection would be l = c*B*I so that this formulation is just the rate or force of infection times the number of people experiencing this rate. Put the value of the contact rate c = 4 and the value of the transmission probability B = 1/8.

Formulate the recovery flow as I divided by D. Remember that the average duration is one divided by the recovery so the recovery rate is 1 divided by the average duration. We can enter either the average duration or the recovery rate as the parameter of our model.

Handin 1.4

Answer the question, "At what rate does each infected individual generate infections in other individuals if everyone they encounter is susceptible?" The answer should be in terms of number of infections per infected person per unit of time.

Make sure you can answer the following.

How do the following affect the rate calculated above?

1. Doubling the number of infected individuals a the start of the epidemic
2. Doubling the duration of infection in the first stages of the epidemic
3. Doubling the transmission probability,
4. Doubling the contact rate.

R and R₀

The rate at which each infected individual generates new infections times the average time that they spend generating new infections is the total number of new infections that they generate. This is called the reproductive number. When all contacts are susceptible, it is called the basic reproductive number. Make a derived variable on your model whose value equals the reproductive number. This is the rate at which infected individuals make contact with all other individuals times the fraction of those individuals that are susceptible, times the transmission probability, times the duration of infection. Examine the value of this derived variable on the same curve as "new infection" and as the number of infected individuals. In this case the number is also the fraction of the population that is infected because we set a fixed total population to have the value of one. Scale these two variables so their peaks are both near the top of the graph.

Make sure you can answer the following

At what value of R does the epidemic reach its peak in terms of the total number of infected individuals present? Explain why this is the case. Explain why the curve of new infections peaks earlier than the curve of the total number of infected individuals. (Hint, to help your search for an explanation, examine the recovery flow on the same graph. Make sure you put the recovery flow and the new infection flow on the same scale. Remember, to set the scale, you double click on the graph, select the graph item you want to scale, click on the arrow to the right of that item, enter the values you want to set on the scale, and then click set.

What ends an epidemic?

Run out the model to 120 time units. Examine in a table the values of S, I, R, New infections, and recovery and reproductive number. Rather than just using two decimal points to examine your output, examine them in floating decimal point form. To do this, double click on the header for a column. You will get a window like the following.

Select Free Float in the precision column.

Examine your table to determine the rate at which each infected individual generates new infected individuals at the very beginning of the epidemic and at the end of the epidemic. You may want to examine a graph as well as a table.

Make sure you can answer the following:

Which goes to zero first, S or I? Why does the value of "I" go to zero? Does the reproductive number go to zero? How is it possible for I to go to zero when the reproductive number does not go to zero?

Epidemic Threshold

Reduce the value of "c" in your model to the point where the value of c*B*D equals 1. Examine the pattern of rise of new infections at this value and at values above and below this value. You will see that the threshold value of "c" for there to be an epidemic is at the value where c*B*D (the basic reproduction number) equals 1.

Examination of model behavior at different values of parameter values is a form of "sensitivity analysis". We are examining how sensitive model behavior is to changes in parameter values. In the case of looking to see how model behavior changes qualitatively at different parameter values, we are also looking for thresholds. To facilitate your examination of model behavior as values of parameters change, Stella IIÔ has sensitivity specs. Pull down the Run menu to "Sensi Specs" and you get a window that looks like the following:

Click on "c" in the allowable list and move it to the selected list in the same way you did for graphs and tables. Now click on "c" when it is in the selected list. We will leave the value for the # of runs at three but will select specific values of c for these three runs rather than just picking a bottom level, a top level, and letting Stella IIÔ set the intermediate levels. To set each value, select "Ad hoc" under the variation type. Then enter 1.8 as your first value and set it. Then enter 2.0 as the second value and set it. (This is the value that sets R₀ equal to 1.) Enter 2.2 as your third value and click OK.

Now create a new graph or a new graph page on your old graph icon. To create a new page on the old icon, double click a graph. Then click the new arrow on the bottom right. The kind of graph we want to construct is a comparative graph so at the top left of the window you must select this type of graph. You can only enter one variable into this type of graph. Enter "new infection" and click OK. Now when you select Run under the run menu, you will see S-run. Let it rip. If you are running to 90 time units, you will see that when R₀ = 1, a constant number of newly infected individuals is generated for some time but there is no epidemic. Eventually the number of new infections will go to zero, but the fraction of the population that is susceptible will have to fall first. Set the simulation time to 600 (remember this is under the run menu and then under time specs). Changing the time specs will clean your comparative graphs. Otherwise, comparative graphs, unlike regular graphs, just keep accumulating new curves on top of old curves. Another way to clear a comparative graph is under the "Model" menu to select "restore" and then select "graphs and tables".

To examine model behavior more thoroughly at threshold, turn off the sensi specs. You do this by selecting Sensi-specs under the run menu and then clicking the box where "Sensitivity On" is checked to uncheck it. Set the value of "c" on your diagram to 2 by double clicking on "c". Now go examine a table with values of I, R, and "new infection". Run out to 900 time units. Print out values every 10 time units. Remember that to do this you must have the table set up window. If you are not already in this window, double click on your table to get to it. Uncheck the "Every DT" box and enter the value of 10 for the interval. Examine the values in your table in the "free float" format. Remember, to do this you double click on the header at the top of the column. Note that even if there is no epidemic at this threshold value, transmission will be sustained for some time. When the starting value for I is very low, however, the total number of individuals infected will be small. Check out how many are infected if you start with 0.9 of the population in S and 0.1 in I.

When "c" is 10% above its threshold value, calculate the time it takes the number of infected individuals to double and to quadruple at the start of the epidemic.

Set the value of c back to 4 and examine the behavior of your transmission system model when B is around its threshold value for epidemic take off. Again when B is 10% above its threshold value, calculate the time it takes the number of infected individuals to double and to quadruple at the start of the epidemic. Then set B back to its original value and examine behavior when D is set around its threshold. Again when D is 10% above its threshold value, calculate the time it takes the number of infected individuals to double and to quadruple at the start of the epidemic.

Can the threshold be characterized by any combination of the values c, B, and D whose product equals 1.

Rate of epidemic rise

Handin 1.5

Is the rate of epidemic rise the same for all combinations of values of c, B, and D that give the same value? Explain why or why not.

Final Epidemic Size

At various combinations of c*B*D that equal 2, run the epidemic out to a long time when there are few new infections being generated. Can the final epidemic size be determined just by knowing R₀ and not the values of the individual parameters that make up R₀? Make sure you have varied D and run the transmission for long periods of time to reach your conclusion. Also make sure you are using a DT that will not change your results if you cut it half.

Handin 1.6

Calculate the fraction of the population that is infected at the end of an epidemic at the following values of R₀. 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 8.0, 16.0. Imagine that you have been the health officer assigned to improve hygiene to so that epidemics will have less of an impact. You do a great job and reduce the contact rate to one fourth of its previous level. It so happens, however, that your efforts reduced R₀ from 16 to 4. You are fired for failure. Your replacement is just as good as you are and she reduces R₀ from 4 to 1. Who should get more credit for controlling epidemic risk?

We will now construct two models that have ways to generate new susceptible individuals in the population. The questions we will address with each of these models are:

1. What determines the endemic level of infection?
2. What are some things that cause an endemic infection to have a constant level or to have epidemics that fluctuate around an endemic level?
3. What determines the threshold at which an infection will or will not sustain itself in a population?

What follows assumes that you have learned to manipulate Stella IIÔ from the first part of this exercise.

Construct a model like the following:

Note that the only new element is the birth and death rate "M" and the birth and death flows. The birth flow is the birth and death rate times the total population and the death flows are the death rates times the population dying at that rate. In order to keep the diagram clean, I have used the ghost function at the top right in order to put ghosts in places where I would not have my connectors crossing other things. Just click on the ghost tool, click on the thing you want to ghost, and click where you want to deposit it. The population stays a constant value because the births equal the deaths. The births are always susceptible, however, so that the fraction of the population that is susceptible will build up. The reproductive number is just the basic reproductive number times the fraction of the population that is susceptible.

Set M = to 1/(70*12). Then the average lifetime L=70*12 months or 70 years. This corresponds to an average 70 year survival where our time unit is months. Leave the c=4, D=4, and B=1/8. Note that although D=4, the average duration of infection is now less than 4 because there are two reasons one leaves the I state. One may recover or one may die. The total flow rate out of the I state is M + 1/D or 1/L + 1/D or (D+L)/LD so the average duration is LD/(D+L). Given the fact that L is so much greater than D, this is only very slightly less than D. This average duration, {LD/(D+L)}, is what determines R₀ in this transmission system model.

Make a graph of I and run the simulation for 6000 months. Set the scale of I from 0 to 0.16. You should get a graph that looks like the following:

Let us examine the first epidemic peak more closely. Double click on the graph. In the bottom left corner of the window that appears, choose to display only that part of the graph from 0 to 100 months. Now you see an epidemic that looks very much like the epidemic when we did not have a birth and death process. Go back and choose to examine the second epidemic in the 900 month to 1000 month time period. The second epidemic starts out much slower and does not reach as high of a level.

Before you read on, please think about why this is and give your explanation. One of the major purposes of a modeling tool like Stella IIÔ is to help with your thinking process. There are two practices which will enhance the ability of Stella IIÔ to help your thinking process. The first is to always predict what you will see before you hit the run button. The second is to always pursue an explanation for what you see.

Damped cycles

To help your thinking about why this model generated damped cycles of infection level, examine the values of S and of the reproductive number at the same time as you examine the level of ". Because we are getting big changes and important differences at the very low end of the infection levels may not be very evident on a graph, it may be better to examine a table than a graph to get the best idea of what is going on.

Handin 1.7

Describe the relationships on this graph and provide an explanation as to why the epidemic cycles dampen out.

Equilibrium infection level

To reach an absolute equilibrium, you will have to run this simulation out quite a bit. Equilibrium means that flows into every compartment equal the flows out of every compartment. Examine simulation results on a table to see whether this is the case. In order for the number of infections to stay constant, each infected individual will have to generate exactly one new infection over the course of their infection. That means that the reproduction number must be one. The reproduction number is the fraction of the population that is susceptible times the basic reproduction number. Thus, if we know the basic reproduction number, we should be able to predict what fraction of the population is susceptible at equilibrium. Check that out in this case where we have an R₀ of 2.

In the transmission system just modeled, new susceptible individuals arose through a birth process. Let us examine a different transmission system where there is no immunity gained from infection. Upon recovery from the infected state, one returns to the susceptible state. Since one goes from Susceptible to Infectious to Susceptible, this is called the SIS model of infection.

Most models of gonorrhea are SIS models that assume that infection does not stimulate immunity. Field study data show that people who have recently had a gonorrhea infection have the same rate of new gonorrhea infections as people who have not had an infection. But these observations don't mean that gonorrhea doesn't induce immunity. Every infection induces some immunity.

Two different things could explain why gonorrhea looks like it induces no immunity when in fact it does. First, our detection methods do not distinguish many variants of the gonorrhea organism and the immunity gained is only to specific variants. Although good immunity might be gained to the specific variant, there are so many other variants that the overall rate hardly seems to go down. Another reason that people who have just had an infection have the same rate of infection as people who have not had an infection is that they are probably at higher risk. Many aspects of contact patterns that are hard to measure put one at risk. When someone gets an infection, they have shown that they are in a risky place in the sexual transmission system. Most likely they will still be in this risky place after recovering from their infection. This unmeasured increased risk balances out the immunity that they have gained and makes it look like there was no immunity.

The considerations in the last paragraph indicate that we can't immediately generalize everything that we learn about an SIS transmission system to gonorrhea. But we can still get insights about gonorrhea that will help us think about control programs and endemic levels of infection by examining SIS models. The SIS model we will construct assumes that differences in susceptibility between individuals who have never been infected and individuals who were infected but lost their immunity are unimportant.

To construct an SIS model, we make people flow back and forth between an S and an I compartment. We can do that with a "BiFlow" pipe as follows.

The flow in the open arrow direction is entered as a positive term in the flow regulator labled "NewInfection&Recovery". The flow in the blacked in direction is entered as a negative value. Thus after flow regulator window should be filled in as follows:

Note that the BIFLOW is checked. In the flow regulator, the term before the negative sign (S*c*B*I) is the new infection rate. To reiterate, this can be viewed as the number of potentially transmitting contacts made by the I times the fraction of those contacts that are to susceptible individuals. The term after the negative sign, (I/D) is the number of I times the recovery rate (which is 1/D). When we use the same c, B, and D as in the SIR model (4, 1/8, 4), but now new susceptible individuals are injected at the rate of I/D instead of 1/(70*12), we generate the following curve.

Note that S comes to the same equilibrium level as the SIR transmission system with vital dynamics.

Let us now answer the questions we set out to address. What determines the endemic level of infection? Here, as in the SIR with vital dynamics, the endemic level of susceptibility is related to R₀. When infection is endemic, the reproduction number R must equal 1. {Sorry about using R for two different things but this happens all the time in this literature so I thought I would put on alert for it.}. Therefore, S must equal 1/R₀ since R = R₀S. S will be slightly different for the two models because R₀ will be slightly different.

With the susceptible level of the population determined by R₀, the next issue is what fraction of the remaining population will be in the I state. For the SIS transmission system the answer is easy -- "All of the rest". For the SIR with vital dynamics, can you figure out what fraction of the population will be in the I state. The best way to figure this out is to recognize that the flow into and out of the R compartment must be equal at equilibrium. This provides you with an algebraic equation which you can solve for I. RM = I/D or R/L = I/D {where R is now the population that is immune}. So the fraction of the population that is not S is divided into R and I in the ratio R/I = L/D.

Another question we set out to address was "What are some things that cause an endemic infection to have a constant level or to have epidemics that fluctuate around an endemic level?"

Handin 1.8

Consider the above two models and their behavior and present an explanation as to why there were epidemic cycles in one model before endemic levels were reached but there were no epidemic cycles but only a constant rise to endemicity in the other model.

The last question we were to address was " What determines the threshold at which an infection will or will not sustain itself in a population?" In 605 you learned that R₀ determined endemic thresholds as well as epidemic thresholds. You might want to examine that issue here for the two models presented. You can perform experiments similar to those we used to examine epidemic thresholds.