key: cord-0037149-o4s1evca
authors: Brauer, Fred; Castillo-Chavez, Carlos
title: Models for Endemic Diseases
date: 2011-10-03
journal: Mathematical Models in Population Biology and Epidemiology
DOI: 10.1007/978-1-4614-1686-9_10
sha: 42cd440c5cfca0b0901f9d2b7c423e145cda8d32
doc_id: 37149
cord_uid: o4s1evca

We have been studying SIR models, in which the transitions are from susceptible to infective to removed, with the removal coming through recovery with full immunity (as in measles) or through death from the disease (as in plague, rabies, and many other animal diseases). Another type of model is an SIS model in which infectives return to the susceptible class on recovery because the disease confers no immunity against reinfection. Such models are appropriate for most diseases transmitted by bacterial or helminth agents, and most sexually transmitted diseases (including gonorrhea, but not such diseases as AIDS, from which there is no recovery). One important way in which SIS models differ from SIR models is that in the former there is a continuing flow of new susceptibles, namely recovered infectives. Later in this chapter we will study models that include demographic effects, namely births and deaths, another way in which a continuing flow of new susceptibles may arise.

with r = β N −γ and with K = N −γ/β , our qualitative result tells us that if β N −γ < 0 or β N/γ < 1, then all solutions of the model (10.2) with nonnegative initial values except the constant solution I = K − β /γ approach the limit zero as t → ∞, while if β K/γ > 1, then all solutions with nonnegative initial values except the constant solution I = 0 approach the limit K −γ/β > 0 as t → ∞. Thus there is always a single limiting value for I, but the value of the quantity β K/γ determines which limiting value is approached, regardless of the initial state of the disease. In epidemiological terms this says that if the quantity β K/γ is less than one, the infection dies out in the sense that the number of infectives approaches zero. For this reason the equilibrium I = 0, which corresponds to S = K, is called the disease-free equilibrium. On the other hand, if the quantity β K/γ exceeds one, the infection persists. The equilibrium I = K − γ/β , which corresponds to S = γ/β , is called an endemic equilibrium.

As we have seen in epidemic models, the dimensionless quantity β K/γ is called the basic reproduction number or contact number for the disease, and it is usually denoted by R 0 . In studying an infectious disease, the determination of the basic reproduction number is invariably a vital first step. The value one for the basic reproduction number defines a threshold at which the course of the infection changes between disappearance and persistence. Since β K is the number of contacts made by an average infective per unit time and 1/γ is the mean infective period, R 0 represents the average number of secondary infections caused by each infective over the course of the infection. Thus, it is intuitively clear that if R 0 < 1, the infection should die out, while if R 0 > 1, the infection should establish itself. In more highly structured models than the simple one we have developed here, the calculation of the basic reproduction number may be much more complicated, but the essential concept remains, that of the basic reproduction number as the number of secondary infections caused by an average infective over the course of the disease. However, there is a difference from the behavior of epidemic models. Here, the basic reproduction number determines whether the infection establishes itself or dies out, whereas in the SIR epidemic model the basic reproduction number determines whether there will be an epidemic.

We were able to reduce the system of two differential equations (10.1) to the single equation (10. 2) because of the assumption that the total population S + I is constant. If there are deaths due to the disease, this assumption is violated, and it would be necessary to use a two-dimensional system as a model. We shall consider this in a more general context in the next section.

A model for a disease from which infectives recover with no immunity against reinfection and that includes births and deaths is

3)

describing a population with a density-dependent birth rate Λ (N) per unit time, a proportional death rate μ in each class, and a rate α of departure from the infective class through recovery or disease death and with a fraction f of infectives recovering with no immunity against reinfection. In this model, if f < 1, the total population size is not constant, and K represents a carrying capacity, or maximum possible population size, rather than a constant population size. It is easy to verify that

If we add the two equations of (10.11) and use N = S + I, we obtain

We will carry out the analysis of the SIS model only in the special case f = 1, so that N is the constant K. The system (10.11) is asymptotically autonomous and its asymptotic behavior is the same as that of the single differential equation

where S has been replaced by K − I. But (10.4) is a logistic equation that is easily solved analytically by separation of variables or qualitatively by an equilibrium analysis. We find that I → 0 if Kβ (K) < (μ + α), or R 0 < 1 and I → I ∞ > 0 with

The endemic equilibrium, which exists if R 0 > 1, is always asymptotically stable. If R 0 < 1, the system has only the disease-free equilibrium and this equilibrium is asymptotically stable. The verification of these properties remains valid if there are no births and deaths. This suggests that a requirement for the existence of an endemic equilibrium is a flow of new susceptibles either through recovery without immunity against reinfection or through births. Exercises 1. Modify the SIS model (10.1) to the situation in which there are two competing strains of the same disease, generating two infective classes I 1 , I 2 under the assumption that coinfections are not possible. Does the model predict coexistence of the two strains or competitive exclusion? 2. * A communicable disease from which infectives do not recover may be modeled by the pair of differential equations S = −β SI, I = β SI.

Show that in a population of fixed size K such a disease will eventually spread to the entire population. 3. * Consider a disease spread by carriers who transmit the disease without exhibiting symptoms themselves. Let C(t) be the number of carriers and suppose that carriers are identified and isolated from contact with others at a constant per capita rate α, so that C = −αC. The rate at which susceptibles become infected is proportional to the number of carriers and to the number of susceptibles, so that S = −β SC. Let C 0 and S 0 be the numbers of carriers and susceptibles, respectively, at time t = 0. 4. * Consider a population of fixed size K in which a rumor is being spread by word of mouth. Let y(t) be the number of people who have heard the rumor at time t and assume that everyone who has heard the rumor passes it on to r others in unit time. Thus, from time t to time (t + h). the rumor is passed on hry(t) times, but a fraction y(t)/K of the people who hear it have already heard it, and thus there are only hry(t) K−y(t) K people who hear the rumor for the first time. Use these assumptions to obtain an expression for y(t + h) − y(t), divide by h, and take the limit as h → 0 to obtain a differential equation satisfied by y(t). 5. At 9 AM one person in a village of 100 inhabitants starts a rumor. Suppose that everyone who hears the rumor tells one other person per hour. Using the model of the previous exercise, determine how long it will take until half the village has heard the rumor. 6. * If a fraction λ of the population susceptible to a disease that provides immunity against reinfection moves out of the region of an epidemic, the situation may be modeled by a system S = −β SI − λ S, I = β SI − αI.

Show that both S and I approach zero as t → ∞.

Epidemics that sweep through a population attract much attention and arouse a great deal of concern. We have omitted births and deaths in our description of epidemic models because the time scale of an epidemic is generally much shorter than the demographic time scale. In effect, we have used a time scale on which the number of births and deaths in unit time is negligible.

However, there are diseases that are endemic in many parts of the world and that cause millions of deaths each year. To model a disease that may be endemic we need to think on a longer time scale and include births and deaths. A reference describing the properties of many endemic diseases is Anderson and May (1991) . For diseases that are endemic in some region, public health physicians would like to be able to estimate the number of infectives at a given time as well as the rate at which new infections arise. The effects of quarantine or vaccination in reducing the number of victims are of importance, just as in the treatment of epidemics. In addition, the possibility of defeating the endemic nature of the disease and thus controlling or even eradicating the disease in a population is worthy of study.

A model of Kermack and McKendrick (1932) includes births in the susceptible class proportional to total population size and a death rate in each class proportional to the number of members in the class. This model allows the total population size to grow exponentially or die out exponentially if the birth and death rates are unequal. It is applicable to such questions as whether a disease will control the size of a population that would otherwise grow exponentially. We shall return to this topic, which is important in the study of many diseases in less-developed countries with high birth rates. To formulate a model in which total population size remains bounded we could follow the approach suggested by Hethcote (1976) in which the total population size K is held constant by making birth and death rates equal. Such a model is

Because S + I + R = K, we can view R as determined when S and I are known and consider the two-dimensional system

We shall examine a slightly more general SIR model with births and deaths for a disease that may be fatal to some infectives. For such a disease, the class R of removed members should contain only recovered members, not members removed by death from the disease. It is not possible to assume that the total population size remains constant if there are deaths due to disease; a plausible model for a disease that may be fatal to some infectives must allow the total population to vary in time. The simplest assumption to allow this is a constant birth rate Λ , but in fact the analysis is quite similar if the birth rate is a function Λ (N) of total population size N.

Let us analyze the model

where N = S + I + R, with a mass action contact rate, a constant number of births Λ per unit time, a proportional natural death rate μ in each class, and a rate of recovery or disease death α of infectives with a fraction f of infectives recovering with immunity against reinfection. In this model, if f = 1 the total population size approaches a limit K = Λ /μ. Then K is the carrying capacity of the population. If f < 1, the total population size is not constant, and K again represents a carrying capacity or maximum possible population size, rather than a population size. We view the first two equations as determining S and I, and then consider the third equation as determining N once S and I are known. This is possible because N does not enter into the first two equations. Instead of using N as the third variable in this model we could have used R, and the same reduction would have been possible. If the birth or recruitment rate Λ (N) is a function of total population size, then in the absence of disease the total population size N satisfies the differential equation

The carrying capacity is the limiting population size K, satisfying

The condition Λ (K) < μ ensures the asymptotic stability of the equilibrium population size K. It is reasonable to assume that K is the only positive equilibrium, so that Λ (N) > μN for 0 ≤ N ≤ K. For most population models,

However, if Λ (N) represents recruitment into a behavioral class, as would be natural for models of sexually transmitted diseases, it would be plausible to have Λ (0) > 0, or even to consider Λ (N) to be a constant function. If Λ (0) = 0, we require Λ (0) > μ, because if this requirement is not satisfied, there is no positive equilibrium and the population would die out even in the absence of disease. We have used a mass action contact rate for simplicity, even though a more general contact rate would give a more accurate model, just as in the epidemics considered in the preceding section. With a general contact rate and a density-dependent birth rate we would have a model

If f = 1, so that there are no disease deaths, the equation for N is

so that N(t) approaches a limiting population size K. The theory of asymptotically autonomous systems [Castillo-Chavez and Thieme (1993), Markus (1956), Thieme (1994), Thieme and Castillo-Chavez (1993)] implies that if N has a constant limit then the system is equivalent to the system in which N is replaced by this limit. Then the system (10.6) is the same as the system (10.5) with β replaced by the constant β (K) and N by K, and Λ (N) replaced by the constant Λ (K) = μK.

We shall analyze the model (10.5) qualitatively. In view of the remark above, our analysis will also apply to the more general model (10.6) if there are no disease deaths. Analysis of the system (10.6) with f < 1 is much more difficult. We will confine our study of (10.6) to a description without details.

The first stage of the analysis is to note that the model (10.5) is a properly posed problem. That is, since S ≥ 0 if S = 0 and I ≥ 0 if I = 0, we have S ≥ 0, I ≥ 0 for t ≥ 0, and since N ≤ 0 if N = K, we have N ≤ K for t ≥ 0. Thus the solution always remains in the biologically realistic region S ≥ 0, I ≥ 0, 0 ≤ N ≤ K if it starts in this region. By rights, we should verify such conditions whenever we analyze a mathematical model, but in practice this step is frequently overlooked.

Our approach will be to identify equilibria (constant solutions) and then to determine the asymptotic stability of each equilibrium. Asymptotic stability of an equilibrium means that a solution starting sufficiently close to the equilibrium remains close to the equilibrium and approaches the equilibrium as t → ∞, while instability of the equilibrium means that there are solutions starting arbitrarily close to the equilibrium that do not approach it. To find equilibria (S ∞ , I ∞ ) we set the right side of each of the two equations equal to zero. The second of the resulting algebraic equations factors, giving two alternatives. The first alternative is I ∞ = 0, which will give a disease-free equilibrium, and the second alternative is β S ∞ = μ + α, which will give an endemic equilibrium, provided that β S ∞ = μ + α < β K. If I ∞ = 0 the other equation gives S ∞ = K = Λ /μ. For the endemic equilibrium the first equation gives

We linearize about an equilibrium (S ∞ , I ∞ ) by letting y = S − S ∞ , z = I − I ∞ , writing the system in terms of the new variables y and z and retaining only the linear terms in a Taylor expansion. We obtain a system of two linear differential equations,

The coefficient matrix of this linear system is

We then look for solutions whose components are constant multiples of e λt ; this means that λ must be an eigenvalue of the coefficient matrix. The condition that all solutions of the linearization at an equilibrium tend to zero as t → ∞ is that the real part of every eigenvalue of this coefficient matrix be negative. At the disease-free equilibrium the matrix is

which has eigenvalues −μ and β K − μ − α. Thus, the disease-free equilibrium is asymptotically stable if β K < μ + α and unstable if β K > μ + α. Note that this condition for instability of the disease-free equilibrium is the same as the condition for the existence of an endemic equilibrium.

In general, the condition that the eigenvalues of a 2 × 2 matrix have negative real part is that the determinant be positive and the trace (the sum of the diagonal elements) be negative. Since β S ∞ = μ + α at an endemic equilibrium, the matrix of the linearization at an endemic equilibrium is

and this matrix has positive determinant and negative trace. Thus, the endemic equilibrium, if there is one, is always asymptotically stable. If the quantity

is less than one, the system has only the disease-free equilibrium, and this equilibrium is asymptotically stable. In fact, it is not difficult to prove that this asymptotic stability is global, that is, that every solution approaches the disease-free equilibrium. If the quantity R 0 is greater than one then the disease-free equilibrium is unstable, but there is an endemic equilibrium that is asymptotically stable. Again, the quantity R 0 is the basic reproduction number. It depends on the particular disease (determining the parameter α) and on the rate of contacts, which may depend on the population density in the community being studied. The disease model exhibits a threshold behavior: If the basic reproduction number is less than one, the disease will die out, but if the basic reproduction number is greater than one, the disease will be endemic. Just as for the epidemic models of Chapter 1, the basic reproduction number is the number of secondary infections caused by a single infective introduced into a wholly susceptible population, because the number of contacts per infective in unit time is β K, and the mean infective period (corrected for natural mortality) is 1/(μ + α).

There are two aspects of the analysis of the model (10.6) that are more complicated than the analysis of (10.5). The first is in the study of equilibria. Because of the dependence of Λ (N) and β (N) on N, it is necessary to use two of the equilibrium conditions to solve for S and I in terms of N and then substitute into the third condition to obtain an equation for N. Then by comparing the two sides of this equation for N = 0 and N = K it is possible to show that there must be an endemic equilibrium value of N between 0 and K.

The second complication is in the stability analysis. Since (10.6) is a threedimensional system that cannot be reduced to a two-dimensional system, the coefficient matrix of its linearization at an equilibrium is a 3 × 3 matrix, and the resulting characteristic equation is a cubic polynomial equation of the form

The Routh-Hurwitz conditions

are necessary and sufficient conditions for all roots of the characteristic equation to have negative real part. A technically complicated calculation is needed to verify that these conditions are satisfied at an endemic equilibrium for the model (10.6).

The asymptotic stability of the endemic equilibrium means that the compartment sizes approach a steady state. If the equilibrium had been unstable, there would have been a possibility of sustained oscillations. Oscillations in a disease model mean fluctuations in the number of cases to be expected, and if the oscillations have long period, that could also mean that experimental data for a short period would be quite unreliable as a predictor of the future. Epidemiological models that incorporate additional factors may exhibit oscillations. A variety of such situations is described in [Hethcote and Levin (1989), Hethcote, Stech, and van den Driessche (1981)].

The epidemic models of the previous chapter also exhibited a threshold behavior of a slightly different kind. For these models, which were SIR models without births or natural deaths, the threshold distinguished between a dying out of the disease and an epidemic, or short-term spread of disease.

From the third equation of (10.5) we obtain

where N = S + I + R. From this we see that at the endemic equilibrium, N = K − (1 − f )αI/μ, and the reduction in the population size from the carrying capacity K is

The parameter α in the SIR model may be considered as describing the pathogenicity of the disease. If α is large, it is less likely that R 0 > 1. If α is small, then the total population size at the endemic equilibrium is close to the carrying capacity K of the population. Thus, the maximum population decrease caused by disease will be for diseases of intermediate pathogenicity.

Numerical simulations indicate that the approach to endemic equilibrium for an SIR model is like a rapid and severe epidemic if the epidemiological and demographic time scales are very different. The same happens in the SIS model. If there are few disease deaths, the number of infectives at endemic equilibrium may be substantial, and there may be damped oscillations of large amplitude about the endemic equilibrium. For both the SIR and SIS models we may write the differential equation for I as

which implies that whenever S exceeds its endemic equilibrium value S ∞ , I is increasing, and epidemic-like behavior is possible. If R 0 < 1 and S < K, it follows that I < 0, and thus I is decreasing. Thus, if R 0 < 1, I cannot increase and no epidemic can occur.

Next, we will turn to some applications of SIR and SIS models.

In order to prevent a disease from becoming endemic it is necessary to reduce the basic reproduction number R 0 below one. This may sometimes be achieved by immunization. If a fraction p of the Λ newborn members per unit time of the population is successfully immunized, the effect is to replace K by K(1 − p), and thus to reduce the basic reproduction number to R 0 (1 − p). The requirement R 0 (1 − p) < 1 gives 1 − p < 1/R 0 , or

A population is said to have herd immunity if a large enough fraction has been immunized to ensure that the disease cannot become endemic. The only disease for which this has actually been achieved worldwide is smallpox, for which R 0 is approximately 5, so that 80 percent immunization does provide herd immunity. For measles, epidemiological data in the United States indicate that R 0 for rural populations ranges from 5.4 to 6.3, requiring vaccination of 81.5 percent to 84.1 percent of the population. In urban areas R 0 ranges from 8.3 to 13.0, requiring vaccination of 88.0 percent to 92.3 percent of the population. In Great Britain, R 0 ranges from 12.5 to 16.3, requiring vaccination of 92 percent to 94 percent of the population. The measles vaccine is not always effective, and vaccination campaigns are never able to reach everyone. As a result, herd immunity against measles has not been achieved (and probably never can be). Since smallpox is viewed as more serious and requires a lower percentage of the population be immunized, herd immunity was attainable for smallpox. In fact, smallpox has been eliminated; the last known case was in Somalia in 1977, and the virus is maintained now only in laboratories. The eradication of smallpox was actually more difficult than expected, because high vaccination rates were achieved in some countries but not everywhere, and the disease persisted in some countries. The eradication of smallpox was possible only after an intensive campaign for worldwide vaccination [Hethcote (1978)].

In order to calculate the basic reproduction number R 0 for a disease, we need to know the values of the contact rate β and the parameters μ, K, and α. The parameters μ, K, and α can usually be measured experimentally, but the contact rate β is difficult to determine directly. There is an indirect means of estimating R 0 in terms of the life expectancy and the mean age at infection that enables us to avoid having to estimate the contact rate. In this calculation, we will assume that β is constant, but we will also indicate the modifications needed when β is a function of total population size N. The calculation assumes exponentially distributed life spans and infective periods. In fact, the result is valid so long as the life span is exponentially distributed, but if the life span is not exponentially distributed, the result could be quite different.

Consider the "age cohort" of members of a population born at some time t 0 and let a be the age of members of this cohort. If y(a) represents the fraction of members of the cohort who survive to age (at least) a, then the assumption that a fraction μ of the population dies per unit time means that y (a) = −μy(a). Since y(0) = 1, we may solve this first order initial value problem to obtain y(a) = e −μa . The fraction dying at (exactly) age a is −y (a) = μy(a). The mean life span is the average age at death, which is ∞ 0 a[−y (a)]da, and if we integrate by parts, we find that this life expectancy is

Since y(a) = e −μa , this reduces to 1/μ. The life expectancy is often denoted by L, so that we may write

The rate at which surviving susceptible members of the population become infected at age a and time t 0 + a is β I(t 0 + a). Thus, if z(a) is the fraction of the age cohort alive and still susceptible at age a, z (a) = −[μ + β I(t 0 + a)]z(a). Solution of this linear first-order differential equation gives

The mean length of time in the susceptible class for members who may become infected, as opposed to dying while still susceptible, is

and this is the mean age at which members become infected. If the system is at an equilibrium I ∞ , this integral may be evaluated, and the mean age at infection , denoted by A, is given by

For our model the endemic equilibrium is

and this implies

This relation is very useful in estimating basic reproduction numbers. For example, in some urban communities in England and Wales between 1956 and 1969 the average age of contracting measles was 4.8 years. If life expectancy is assumed to be 70 years, this indicates R 0 = 15.6. The derivation of A = 1/β I ∞ is obtained from considering surviving susceptible members at each age. This is the value that would be obtained from data giving the fraction of susceptibles at each age. However, if average age at infection has the normal meaning of average age at which those people who become infected do become infected, then the calculation would be different. The susceptible population at age a is a fraction e −(μ+β I ∞ ) of the number of newborn members, and the incidence of new infections is β I ∞ e −(μ+β I ∞ ) . This would lead to an average age at infection

If β is a function β (N) of total population size, the relation (10.10) becomes

If disease mortality does not have a large effect on total population size, in particular if there is no disease mortality, this relation is very close to (10.10). The relation between age at infection and basic reproduction number indicates that measures such as inoculations, which reduce R 0 , will increase the average age at infection. For diseases such as rubella (German measles), whose effects may be much more serious in adults than in children, this indicates a danger that must be taken into account: While inoculation of children will decrease the number of cases of illness, it will tend to increase the danger to those who are not inoculated or for whom the inoculation is not successful. Nevertheless, the number of infections in older people will be reduced, although the fraction of cases that are in older people will increase.

Many common childhood diseases, such as measles, whooping cough, chicken pox, diphtheria, and rubella, exhibit variations from year to year in the number of cases.

These fluctuations are frequently regular oscillations, suggesting that the solutions of a model might be periodic. This does not agree with the predictions of the model we have been using in this section; however, it would not be inconsistent with solutions of the characteristic equation, which are complex conjugate with small negative real part corresponding to lightly damped oscillations approaching the endemic equilibrium. Such behavior would look like recurring epidemics. If the eigenvalues of the matrix of the linearization at an endemic equilibrium are −u ± iv, where i 2 = −1, then the solutions of the linearization are of the form Be −ut cos(vt + c), with decreasing "amplitude" Be −ut and "period" 2π v . For the model (10.5) we recall from (10.7) that at the endemic equilibrium we have

and from (10.8), the matrix of the linearization is

The eigenvalues are the roots of the quadratic equation

If the mean infective period 1/(γ + α) is much shorter than the mean life span 1/μ, we may neglect the terms that are quadratic in μ. Thus, the eigenvalues are approximately

and these are complex with imaginary part μ(R 0 − 1)(γ + α). This indicates oscillations with period approximately

We use the relation μ(R 0 − 1) = μL/A and the mean infective period τ = 1/(γ + α) to see that the interepidemic period T is approximately 2π √ Aτ. Thus, for example, for recurring outbreaks of measles with an infective period of 2 weeks or 1/26 year, in a population with a life expectancy of 70 years with R 0 estimated as 15, we would expect outbreaks spaced 2.76 years apart. Also, since the "amplitude" at time t is e −μR 0 t/2 , the maximum displacement from equilibrium is multiplied by a factor e −(15)(2.76)/140 = 0.744 over each cycle. In fact, many observations of measles outbreaks indicate less damping of the oscillations, suggesting that there may be additional influences that are not included in our simple model. To explain oscillations about the endemic equilibrium, a more complicated model is needed. One possible generalization would be to assume seasonal variations in the contact rate. This is a reasonable supposition for a childhood disease most commonly transmitted through school contacts, especially in winter in cold climates. Note, however, that data from observations are never as smooth as model predictions and models are inevitably gross simplifications of reality that cannot account for random variations in the variables. It may be difficult to judge from experimental data whether an oscillation is damped or persistent.

In the model (10.5), the demographic time scale described by the birth and natural death rates μK and μ and the epidemiological time scale described by the rate (α + γ) of departure from the infective class may differ substantially. Think, for example, of a natural death rate μ = 1/75, corresponding to a human life expectancy of 75 years, and epidemiological parameters α = 0 and γ = 25, describing a disease from which all infectives recover after a mean infective period of 1/25 year, or two weeks. Suppose we consider a carrying capacity K = 1000 and take β = 0.1, indicating that an average infective makes (0.1)(1000) = 100 contacts per year. Then R 0 = 4.00, and at the endemic equilibrium we have S ∞ = 250.13, I ∞ = 0.40, R ∞ = 749.47. This equilibrium is globally asymptotically stable and is approached from every initial state.

However, if we take S(0) = 999, I(0) = 1, R(0) = 0, simulating the introduction of a single infective into a susceptible population, and solve the system numerically, we find that the number of infectives rises sharply to a maximum of 400 and then decreases to almost zero in a period of 0.4 year, or about 5 months. In this time interval the susceptible population decreases to 22 and then begins to increase, while the removed (recovered and immune against reinfection) population increases to almost 1000 and then begins a gradual decrease. The size of this initial "epidemic" could not have been predicted from our qualitative analysis of the system (10.5). On the other hand, since μ is small compared to the other parameters of the model, we might consider neglecting μ, replacing it by zero in the model. If we do this, the model reduces to the simple Kermack-McKendrick epidemic model (without births and deaths) of Section 9.2.

If we follow the model (10.5) over a longer time interval, we find that the susceptible population grows to 450 after 46 years, then drops to 120 during a small epidemic with a maximum of 18 infectives, and exhibits widely spaced epidemics decreasing in size. It takes a very long time before the system comes close to the endemic equilibrium and remains close to it. The large initial epidemic conforms to what has often been observed in practice when an infection is introduced into a population with no immunity, such as the smallpox inflicted on the Aztecs by the invasion of Cortez.

If we use the model (10.5) with the same values of β , K, and μ, but take α = 25, γ = 0 to describe a disease fatal to all infectives, we obtain very similar results. Now the total population is S + I, which decreases from an initial size of 1000 to a minimum of 22 and then gradually increases and eventually approaches its equilibrium size of 250.53. Thus, the disease reduces the total population size to one-fourth of its original value, suggesting that infectious diseases may have large effects on population size. This is true even for populations that would grow rapidly in the absence of infection, as we shall see in the next section.

In order to describe a model for a disease from which infectives recover with immunity against reinfection and that includes births and deaths as in the model (10.5), we may modify the model (10.5) by removing the equation for R and moving the term αI describing the rate of recovery from infection to the equation for S . This gives the model

describing a population with a constant number of births μK per unit time, a proportional death rate μ in each class, and a fraction γ of infectives dying from infection and a fraction α of infectives recovering with no immunity against reinfection. In this model, if γ > 0, the total population size is not constant and K represents a carrying capacity, or maximum possible population size, rather than a constant population size. The analysis of the model (10.11) is very similar to that of the SIR model (10.5), except that there is no equation for R to be eliminated. The only difference is the additional term αI in the equation for S , and this does not change any of the qualitative results. As in the SIR model we have a basic reproductive number

and I ∞ given by (10.7), which is asymptotically stable. There are, however, differences that are not disclosed by the qualitative analysis. If the epidemiological and demographic time scales are very different, for the SIR model we observed that the approach to endemic equilibrium is like a rapid and severe epidemic. The same happens in the SIS model, especially if there is a significant number of deaths due to disease. If there are few disease deaths, the number of infectives at endemic equilibrium may be substantial, and there may be oscillations of large amplitude about the endemic equilibrium. For both the SIR and SIS models we may write the differential equation for I as

which implies that whenever S exceeds its endemic equilibrium value, I is increasing, and epidemic-like behavior is possible. If R 0 < 1 and S < K, it follows that I < 0, and thus I is decreasing. Thus, if R 0 < 1 I cannot increase and no epidemic can occur.

1. Recurrent outbreaks of measles and other childhood diseases have previously been explained by an interaction between intrinsic epidemiological forces generating dampened oscillations and seasonal and/or stochastic excitation. The following model shows that isolation or quarantine (i.e., sick individuals stay at home and have a reduced infective impact) can create self-sustained oscillations.

In the model considered here the population is divided into susceptibles (S), infectives (I), isolated or quarantined individuals (Q), and recovered individuals (R), for whom permanent immunity is assumed. Let N denote the total population, and let A = S + I + R denote the active (nonisolated) individuals. The model takes the form

All newborns are assumed to be susceptible. μ is the per capita mortality rate, σ is the per capita infection rate of an average susceptible individual provided everybody else is infected, γ is the rate at which individuals leave the infective class, and ξ is the rate at which individuals leave the isolated class; all are positive constants.

(i) Show that the total population size N is constant.

(ii) Give the meanings of 1/μ, 1/γ, 1/ξ and their units. 

Express the new parameters in terms of the old parameters. Check that all the new parameters and variables are dimensionless. (iv) Study the stability of the equilibrium point (0, 0, 0) and derive a basic reproductive number R 0 . (v) Use a computer algebra system to demonstrate that (10.13) has periodic trajectories. Use the parameter values ν = 0.0002, θ = 0.0156, and ξ close to θ 2 (1 − θ ). You also need to choose proper initial values.

In the SIR models that we have studied, it has been assumed that the immunity received by recovery from the disease is permanent. This is not always true, since there may be a gradual loss of immunity with time. In addition, there are often mutations in a virus, and as a result the active disease strain is sufficiently different from the strain from which an individual has recovered, that the immunity received may wane. Temporary immunity may be described by an SIRS model in which a rate of transfer from R to S is added to an SIR model. For simplicity, we confine our attention to epidemic models, without including births, natural deaths, and disease deaths, but the analysis of models including births and deaths would lead to the same conclusions. Thus we begin with a model

with a proportional rate θ of loss of immunity.

Since N = (S + I + R) = 0, the total population size N is constant, and we may replace R by N − S − I and reduce the model to a two-dimensional system

(10.14)

Equilibria are solutions of the system

and there is a disease-free equilibrium S = α/β ,

The matrix of the linearization of (10.14) at an equilibrium (S, I) is

At the disease-free equilibrium, A has the sign structure

This matrix has negative trace and positive determinant if and only if β N < α, or R 0 < 1. At an endemic equilibrium, the matrix has sign structure

and thus always has negative trace and positive determinant. We see from this that as in other models studied in this chapter, the disease-free equilibrium is asymptotically stable if and only if the basic reproducton number is less than 1 and the endemic equilibrium, which exists if and only if the basic reproduction number exceeds 1, is always asymptotically stable. However, it is possible for a different SIRS model to have quite different behavior.

We consider an SIRS model, that assumes a constant period of temporary immunity following recovery from the infection in place of an exponentially distributed period of temporary immunity. It will turn out that the endemic equilibrium for this model may be unstable, thus giving an example of a generalization that leads to new possibilities for the behavior of a model.

We add the assumption that there is a temporary immunity period of fixed length ω, after which recovered infectives revert to the susceptible class. The resulting model is described by the system

Since N = S + I + R is constant, we may discard the equation for R and use a two-dimensional model 

Equilibria are given by I = 0 or β S = α. There is a disease-free equilibrium S = N, I = 0. There is also an endemic equilibrium for which β S = α. However, the two equations for S and I give only a single equilibrium condition. To determine the endemic equilibrium (S ∞ , I ∞ ) we must write the equation for R in the integrated form

The characteristic equation at an equilibrium is the condition that the linearization at the equilibrium have a solution whose components are constant multiples of e λt . In the ordinary differential equation case this is just the equation that determines the eigenvalues of the coefficient matrix, a polynomial equation, but in the general case it is a transcendental equation. The result on which our analysis depends is that an equilibrium is asymptotically stable if all roots of the characteristic equation have negative real part, or equivalently that the characteristic equation have no roots with real part greater than or equal to zero.

To linearize about an equilibrium (S ∞ , I ∞ ) of (10.15) we substitute

and neglect the quadratic term, giving the linearization

The characteristic equation is the condition on λ that this linearization have a solution

and this is

This reduces to

At the disease-free equilibrium S ∞ = N, I ∞ = 0, this reduces to a linear equation with a single root λ = β N − α, which is negative if and only if R 0 = β N/α < 1. We think of ω and N as fixed and consider β and α as parameters. If α = 0, the equation (10.16) is linear and its only root is −β I ∞ < 0. Thus, there is a region in the (β , α) parameter space containing the β -axis in which all roots of (10.16) have negative real part. In order to find how large this stability region is we make use of the fact that the roots of (10.16) depend continuously on β and α. A root can move into the right half-plane only by passing through the value zero or by crossing the imaginary axis as β and α vary. Thus, the stability region contains the β -axis and extends into the plane until there is a root λ = 0 or until there is a pair of pure imaginary roots λ = ±iy with y > 0. Since the left side of (10.16) is positive and the right side of (10.16) is negative for real λ ≥ 0, there can not be a root λ = 0.

The condition that there is a root λ = iy is To satisfy the first condition it is necessary to have ωα > 1, since | sin ωy| ≤ |ωy| for all y. This implies, in particular, that the endemic equilibrium is asymptotically stable if ωα < 1. In addition, it is necessary to have sin ωy < 0. There is an infinite sequence of intervals on which sin ωy < 0, the first being π < ωy < 2π. For each of these intervals the equations (10.18) define a curve in the (β , α) plane parametrically with y as parameter. The region in the plane below the first of these curves is the region of asymptotic stability, that is, the set of values of β and α for which the endemic equilibrium is asymptotically stable. This curve is shown for ω = 1, N = 1 in Figure 10 .1. Since R 0 = β /α > 1, only the portion of the (β ,α) plane below the line γ = β is relevant. The new feature of the model of this section is that the endemic equilibrium is not asymptotically stable for all parameter values. What is the behavior of the model if the parameters are such that the endemic equilibrium is unstable? A plausible suggestion is that since the loss of stability corresponds to a root λ = iy of the characteristic equation, there are solutions of the model behaving like the real part of e iyt , that is, there are periodic solutions. This is exactly what does happen according to a very general result called the Hopf bifurcation theorem, which says that when roots of the characteristic equation cross the imaginary axis a stable periodic orbit arises.

From an epidemiological point of view, periodic behavior is unpleasant. It implies fluctuations in the number of infectives, which makes it difficult to allocate resources for treatment. It is also possible for oscillations to have a long period. This means that if data are measured over only a small time interval the actual behavior may not be displayed. Thus, the identification of situations in which an endemic equilibrium is unstable is an important problem. 

Determine the basic reproduction number of the model (10.19 ) and find all endemic equilibria.

Many parts of the world experienced very rapid population growth in the eighteenth century. The population of Europe increased from 118 million in 1700 to 187 million in 1800. In the same time period, the population of Great Britain increased from 5.8 million to 9.15 million, and the population of China increased from 150 million to 313 million [McNeill (1976) ]. The population of the English colonies in North America grew much more rapidly than this, aided by substantial immigration from England, but the native population, which had been reduced to one-tenth of their previous size by disease following the early encounters with Europeans and Euro-pean diseases, grew even more rapidly. While some of these population increases may be explained by improvements in agriculture and food production, it appears that an even more important factor was the decrease in the death rate due to diseases. Disease death rates dropped sharply in the eighteenth century, partly from better understanding of the links between illness and sanitation and partly because the recurring invasions of bubonic plague subsided, perhaps due to reduced susceptibility. One plausible explanation for these population increases is that the bubonic plague invasions served to control the population size, and when this control was removed, the population size increased rapidly.

In developing countries it is quite common to have high birth rates and high disease death rates. In fact, when disease death rates are reduced by improvements in health care and sanitation, it is common for birth rates to decline as well, as families no longer need to have as many children to ensure that enough children survive to take care of the older generations. Again, it is plausible to assume that population size would grow exponentially in the absence of disease but is controlled by disease mortality.

The SIR model with births and deaths of Kermack and McKendrick (1932) includes births in the susceptible class proportional to population size and a natural death rate in each class proportional to the size of the class. Let us analyze a model of this type with birth rate r and a natural death rate μ < r. For simplicity we assume that the disease is fatal to all infectives with disease death rate α, so that there is no removed class and the total population size is N = S + I. Our model is

From the second equation we see that equilibria are given by either I = 0 or β S = μ + α. If I = 0, the first equilibrium equation is rS = μS, which implies S = 0, since r > μ. It is easy to see that the equilibrium (0, 0) is unstable. What actually would happen if I = 0 is that the susceptible population would grow exponentially with exponent r − μ > 0. If β S = μ + α, the first equilibrium condition gives

Thus, there is an endemic equilibrium, provided r < α + μ, and it is possible to show by linearizing about this equilibrium that it is asymptotically stable. On the other hand, if r > α + μ, there is no positive equilibrium value for I. In this case we may add the two differential equations of the model to give

and from this we may deduce that N grows exponentially. For this model either we have an asymptotically stable endemic equilibrium or population size grows exponentially. In the case of exponential population growth, we may have either vanishing of the infection or an exponentially growing number of infectives. If only susceptibles contribute to the birth rate, as may be expected if the disease is sufficiently debilitating, the behaviour of the model is quite different. Let us consider the model 

which has the same form as the Lotka-Volterra predator-prey model of population dynamics. This system has two equilibria, obtained by setting the right sides of each of the equations equal to zero, namely (0, 0) and an endemic equilibrium ((μ + α)/β , (r − μ)/β ). It turns out that the qualitative analysis approach we have been using is not helpful, since the equilibrium (0, 0) is unstable and the eigenvalues of the coefficient matrix at the endemic equilibrium have real part zero. In this case, the behavior of the linearization does not necessarily carry over to the full system. However, we can obtain information about the behaviour of the system by a method that begins with the elementary approach of separation of variables for firstorder differential equations. We begin by taking the quotient of the two differential equations and using the relation I S = dI dS to obtain the separable first-order differential equation

Separation of variables gives

Integration gives the relation

where c is a constant of integration. This relation shows that the quantity

is constant on each orbit (path of a solution in the (S, I) plane). Each of these orbits is a closed curve corresponding to a periodic solution. This model is the same as the simple epidemic model of Section 9.2 except for the birth and death terms, and in many examples the time scale of the disease is much faster than the time scale of the demographic process. We may view the model as describing an epidemic initially, leaving a susceptible population small enough that infection cannot establish itself. Then there is a steady population growth until the number of susceptibles is large enough for an epidemic to recur. During this growth stage the infective population is very small, and random effects may wipe out the infection, but the immigration of a small number of infectives will eventually restart the process. As a result, we would expect recurrent epidemics. In fact, bubonic plague epidemics did recur in Europe for several hundred years. If we modify the demographic part of the model to assume limited population growth rather than exponential growth in the absence of disease, the effect would be to give behavior like that of the model studied in the previous section, with an endemic equilibrium that is approached slowly in an oscillatory manner if R 0 > 1. (i) Show that there are always two equilibria, an extinction equilibrium (0, 0) and a coexistence equilibrium with

(ii) Show that both equilibria are unstable, in fact saddle points.

This book is concerned primarily with theoretical models for natural phenomena. Such models necessarily contain parameters whose values must be estimated in order to make it possible to compare model predictions with real-life data. The challenges of connecting models and data and the validation of models are critically important in science. In fact, one of the first recorded modeling contributions in the field of epidemiology was that of Daniel Bernoulli (1766), which focused on the increase in the average life expectancy generated by elimination of a lethal disease, a study related directly to data on a single smallpox outbreak, an intellectual contribution that has been analyzed in detail and expanded by Dietz and Heesterbeeck (2002).

In the current debate concerning whether the United States population should be vaccinated against smallpox (in order to prepare for a possible terrorist attack) critics of mass vacci-nation have argued that the risks associated with widespread usage of the current smallpox vaccine would outweigh the potential benefits. Over 200 years ago, opponents of smallpox variolation in the 18th century used the same arguments. They argued that: (i) inoculation was risky because artificial smallpox could cause mortality, and (ii) inoculation programs could increase the transmission of smallpox (because individuals inoculated with artificial smallpox could transmit smallpox [ibid, 286-287]. Blower concludes that:

It is not clear how influential Bernoulli's paper was in influencing public health policy, but it remains a classic paper as it was the first known mathematical analysis that was used to try to influence public health policy.

[ibid, 287]. Efforts to connect models to data increased dramatically with the onset of the HIV epidemic, with emphasis on the estimation of the incubation period distribu 

The formulation of the standard ordinary least squares formulation usually involves two models: mathematical and statistical. The mathematical model used here is a compartmental mathematical model, which is given in terms of a nonlinear system of ordinary differential equations involving a parameter vector θ . Specifically, we have

where x(t; θ ) ∈ R n denotes the state variable vector at time t and θ ∈ R p denotes the parameter vector. The statistical model linked to the process generated by the dynamical system is formulated under the assumption that the model output and associated random deviations (measurement error) are captured by the random variables

where z(t j ; θ 0 ) denotes the output of the mathematical model. Usually the model output is a functional of the state variable vector, that is, z(t; θ ) = F (x(t; θ )), and in equation (10.23) the model output is evaluated at θ = θ 0 , the "true" parameter vector. The random variables E j model the random deviations away from z(t, θ 0 ) and are assumed to satisfy 1. E j 's are independent and identically distributed random variables; 2. E[E j ] = 0 for every j; 3. var(E j ) = σ 2 0 < ∞, for every j. The ordinary least squares problem arises from efforts to minimize [Y j − z(t j ; θ )] 2 over the set of parameter vectors θ constrained by a prespecified feasible region, here denoted by Θ . The minimizer is a random variable, called the estimator θ OLS and given in this context by

(10.24)

The usefulness of the estimator θ OLS derives from its statistical properties. A classical asymptotic result, in the spirit of a central limit theorem, establishes that for n sufficiently large, this estimator has a p-multivariate normal sampling distribution. In other words,

where Σ 0 = n −1 σ 2 0 Ω −1 0 and Ω 0 = lim n→∞ 1 n χ(θ 0 ) T χ(θ 0 ), provided this limit exists and the matrix Ω 0 is nonsingular. The p× p matrix Σ 0 is the covariance matrix of the θ OLS estimator, while the matrix χ(θ 0 ) is n × p and is called the sensitivity matrix, with its (i, j) entry defined by

[Banks and Tran (2009)]. The theoretical quantities θ 0 , σ 2 0 , and Σ 0 are in general unknown. In practice, one has only the data associated with a single realization y i of the observation process Y i (for i = 1,..., n) and has no option but to solve the minimization problem, that is, the computation of an estimate forθ OLS under these conditions. That is, we carry out the following minimization process:

The estimate obtained is then used to approximate the error assuming constant variance σ 2 0 via the approximation

The covariance matrix Σ 0 can also be approximated usingθ OLS andσ 2 OLS as follows:

We proceed to apply the above OLS methodology to a well-known compartmental model and an influenza data set. In 1978, an outbreak of influenza was reported in a boarding school for boys in the United Kingdom [Communicable Disease Surveillance Centre (1978) ]. The single outbreak is modeled using the classical SIR compartmental model given by the following set of nonlinear differential equations:

The dataset reported in [Communicable Disease Surveillance Centre (1978)] corresponds to the prevalence of influenza, and therefore, the output in this case is modeled by

from which the statistical model for the observation process becomes Y j = I(t j ; θ 0 ) + E j for j = 1,..., n. (10.25) In this case, it was reported that N = 763, and so the conditions at the start of this outbreak can be assumed to be S 0 = 762 and I 0 = 1. The parameter vector to be estimated involves the transmission coefficient β and the recovery rate γ, that is, θ = (β , γ) . The data set in [Communicable Disease Surveillance Centre (1978)] is denoted by y i for i = 1,..., 12, and it is used to computê

The minimization leads to an estimate of θ = (β , γ) . The minimization can be carried out in multiple ways. There are in fact several optimization algorithms (for example, Nelder-Mead simplex), and, for example, the computing software Matlab (Mathworks, Inc) includes the following appropriate optimization routines fminsearch, fmincon, lsqcurvefit, lsqnonlin. 

Numerical values are computed by solving

where

In this case (influenza data), the OLS estimation of the covariance matrix iŝ

with the standard errors calculated by taking the square root of the diagonal entries inΣ OLS . The estimates, within one standard error, are therefore given bŷ

Estimation can go further. For example, it is often important to estimate the effective reproductive number R(t), defined by the expression

In order to estimate R(t), we rewrite the equation for the infective population (10.25) as follows:

Here, R(t) = S N R 0 is the effective reproduction number and R 0 = β /γ is the basic reproductive number). We have R(t) ≤ R 0 , with the upper bound, the basic reproductive number, being achieved only when the entire population is susceptible. Here, R(t) is defined as the product of the transmission rate and the average infectious period, that is, the effective reproductive number gives the average number of secondary infections caused by a single infective individual, at a given susceptible fraction. The prevalence of infection increases or decreases according to whether R 0 (t) is greater than or less than one. Since there is no replenishment of the susceptible population in this model, R 0 (t) decreases over the course of an outbreak as susceptible individuals become infected.

Usingθ OLS we can compute the relevant point-estimate (without uncertainty bounds) curve given by 

As we have seen, in a large variety of models, the behavior when R 0 < 1 is different from the behavior when R 0 > 1. More precisely, as R 0 increases through 1 there is an exchange of stability between the disease-free equilibrium and the endemic equilibrium (which is negative as well as unstable and thus biologically meaningless if R 0 < 1). There is a bifurcation, or change in equilibrium behavior, at R 0 = 1, but the equilibrium infective population size depends continuously on R 0 . Such a transition is called a forward, or transcritical, bifurcation. It is also possible, as we have seen in an SIRS model, that the endemic equilibrium for R 0 > 1 may be unstable, depending on the distribution of infective periods. However, it would be a serious mistake to assume that this normal situation is universal. It has been noted [Dushoff, Huang and Castillo-Chavez(1998), Hadeler and Castillo-Chavez (1995), Hadeler and van den Driessche (1997), Kribs-Zaleta and Velasco-Hernandez (2000) ] that in epidemic models with multiple groups and asymmetry between groups or multiple interaction mechanisms, it is possible to have a very different bifurcation behavior at R 0 = 1. There may be multiple positive endemic equilibria for values of R 0 < 1. Typically, there is an interval of values of R 0 on which there are two asymptotically stable equilibria, one disease-free and one endemic, and an unstable endemic equilibrium between them. Such a situation is called a backward bifurcation at R 0 = 1.

The qualitative behavior of an epidemic system with a backward bifurcation differs from that of a system with a forward bifurcation in at least three important ways.

If there is a forward bifurcation at R 0 = 1, it is not possible for a disease to invade a population if R 0 < 1, because the system will return to the disease-free equilibrium I = 0 if some infectives are introduced into the population. On the other hand, if there is a backward bifurcation at R 0 = 1 and enough infectives are introduced into the population to put the initial state of the system above the unstable endemic equilibrium with R 0 < 1, the system will approach the asymptotically stable endemic equilibrium.

Other differences are observed if the parameters of the system change to produce a change in R 0 . With a forward bifurcation at R 0 = 1 the equilibrium infective population remains zero so long as R 0 < 1 and then increases continuously as R 0 increases. With a backward bifurcation at R 0 = 1, the equilibrium infective population size also remains zero so long as R 0 < 1 but then jumps to the positive endemic equilibrium as R 0 increases through 1. In the other direction, if a disease is being controlled by means that decrease R 0 , it is sufficient to decrease R 0 to 1 if there is a forward bifurcation at R 0 = 1, but it is necessary to bring R 0 well below 1 if there is a backward bifurcation.

We have been assuming homogeneous mixing of members of the population being studied, and this is certainly unrealistically simple. Members of the population may differ, for example, in rate of contact or in location. In the study of sexually transmitted diseases, differences in activity levels are important aspects. Contact rates may be age-dependent, and this would suggest the use of age-structured models. Spatial dependence may take two forms, the local diffusion of members of the population, which would lead to partial differential equations of reaction-diffusion types, and travel between communities, which would lead to patch or metapopulation models.

Models incorporating one or more of these kinds of heterogeneity can be developed and analyzed. Inevitably, their analysis involves more structure, equations, and parameters, as well as more sophisticated mathematical methods.

There are other modes of transmission of communicable diseases that can be described by compartmental models. Some infections can be transmitted vertically, for example, from mother to daughter prior to birth [Busenberg and Cooke (1993) ]. Another form is transmission by a vector. For example, malaria is transmitted back and forth between humans and mosquitoes. Thus an infected mosquito may bite a human and thus infect the human. An uninfected mosquito may bite an infected human and become infected, but infection is not transmitted directly from human to human or from mosquito to mosquito [Ross (1911) . Sexually transmitted diseases that are transmitted by heterosexual contact are also examples of vector transmission, with male and females playing the roles of the two species.

Consider an SIR model eqrefeqsec8a41. For measles, typical parameter choices are μ = 0.02, β = 1800, α = 100, K = 1 (to normalize carrying capacity to 1) [Engbert and Drepper (1984) ].

Show that for these parameter choices R 0 ≈ 18 and to achieve herd immunity would require vaccination of about 95 percent of the susceptible population.

In practice, it is not possible to vaccinate 95 percent of a population because not all members of the population would come to be vaccinated and not all vaccinations are successful. One way to avoid recurring outbreaks of disease is "pulse vaccination" [Agur, Mazor, Anderson, and Danon (1993), Shulgin, Stone, and Agur (1998), Stone, Shulgin, and Agur (2000)]. The basic idea behind pulse vaccination is to vaccinate a given fraction p of the susceptible population at intervals of time T with T (depending on p) chosen to ensure that the number of infectives remains small and approaches zero. In this project we will give two approaches to the calculation of a suitable function T (p).

The first approach depends on the observation that I decreases so long as S < Γ < (μ + γ)/β . We begin by vaccinating pΓ members, beginning with

Then S(t) is greater than the solution of the initial value problem

Solve this initial value problem and show that the solution obeys

Thus a suitable choice for T (p) is

Calculate T (p) for p = m/10 (m = 1, 2,..., 10).

The second approach is more sophisticated. Start with I = 0, S = μK − μS. We let t n = nT (n = 0, 1, 2,...) and run the system for 0 ≤ t ≤ t 1 = T . Then we let S 1 = (1 − p)S(t 1 ). We then repeat, i.e., for t 1 ≤ t ≤ t 2 , S(t) is the solution of S = μK − μS, S(t 1 ) = S 1 , and S 2 = (1 − p)S 1 . We obtain a sequence S n in this way.

and for t n ≤ t ≤ t n+1 ,

Show that the solution is periodic if

and that the periodic solution is

It is possible to show by linearizing about this periodic solution that the periodic solution is asymptotically stable if

If this condition is satisfied, the infective population will remain close to zero.

Show that this stability condition reduces to

Question 6. Use a computer algebra system to graph T (p), where T is defined implicitly by

Compare this expression for T with the one obtained earlier in Question 2 in this project. A larger estimate for a safe value of T would save money by allowing less frequent vaccination pulses.

We model a general discrete-time SIS model with two competing strains in a population with discrete and nonoverlapping generations. This model arises from a particular discretization in time of the corresponding SIS continuous-time stochastic model for two competing strains.

S n population of susceptible individuals in generation n I 1 n population of infected individuals with strain 1 in generation n I 2 n population of infected individuals with strain 2 in generation n T n total population in generation n f recruitment function Parameters μ per capita natural death rate γ i per capita recovery rate for strain i α i per capita infection rate for strain i Construction of the model equations: The model assumes that (i) the disease is not fatal; (ii) all recruits are susceptible and the recruitment function depends only on T n ; (iii) there are no coinfections; (iv) death, infections, and recoveries are modeled as Poisson processes with rates μ, α i , γ i (i = 1, 2); (v) the time step is measured in generations; (vi) the populations change only because of "births" (given by the recruitment function), deaths, recovery, and infection of a susceptible individual for each strain; (vii) individuals recover but do not develop permanent or temporary immunity, that is, they immediately become susceptible again. By assumption we have that the probability of k successful encounters is a Poisson distribution, which in general has the form p(k) = e −β β k /k!, where β is the parameter of the Poisson distribution. In our context, only one success is necessary. Therefore, when there are no successful encounters, the expression p(0) = e −β represents the probability that a given event does not occur. For example, the probability that a susceptible individual does not become infective is Prob(not being infected by strain i) = e −α i I i n , and, Prob(not recovering from strain i) = e −γ i I i n . Hence, Prob(not being infected)= Prob(not being infected by strain 1)Prob(not being infected by strain 2) = e −α 1 I 1 n e −α 2 I 2 n .

Now the probability that a susceptible does become infected is given by 1 − e −α i I i n . Then, Prob(infected by strain i) = Prob(infected). Prob(infected by strain i | infected) =(1 − e −(α 1 I 1 n +α 2 I 2 n ) ) α i I i n α 1 I 1 n +α 2 I 2 n .

(a) Using the above discussion, show that the dynamics are governed by the system S n+1 = f (T n ) + S n e −µ e −(α 1 I 1 n +α 2 I 2 n ) + I 1 n e −µ (1 − e −γ 1 ) + I 2 n e −µ 1 − e −γ 2 , (10.26) I 1 n+1 = α 1 S n I 1 n α 1 I 1 n + α 2 I 2 n e −µ (1 − e −(α 1 I 1 n +α 2 I 2 n ) ) + I 1 n e −µ e −γ 1 ,

where T n = S n + I 1 n + I 2 n . 

Check that this is the case. (d) Study the disease dynamics at a demographic equilibrium, that is, at a point where T ∞ = T ∞ e −µ + f (T ∞ ). Substitute S n = T ∞ − I n 1 − I 2 n where T ∞ is a stable demographic equilibrium, that is, assume T 0 = T ∞ to get the following equations:

System (10.28) describes the dynamics of a population infected with the two strains at a demographic equilibrium. Show that in system (10.28), if R 1 = e −µ T ∞ α 1 1−e −(µ+γ 1 ) < 1 and R 2 = e −µ T ∞ α 2 1−e −(µ+γ 2 ) < 1, then the equilibrium point (0, 0) is asymptotically stable. (e) Interpret biologically the numbers R i , i = 1, 2.

(f) Consider f (T n ) = Λ , where Λ is a constant. Show that T n+1 = Λ + T n e −µ and that

(g) Consider f (T n ) = rT n (1−T n )/k), and show that in this case the total population dynamic is given by,

and that the fixed points are

whenever r + e −μ > 1. 

Consider the following SIS model with dispersion between two patches, Patch 1 and Patch 2, where in Patch i ∈ {1, 2} at generation t, S i (t) denotes the population of susceptible individuals; I i (t) denotes the population of infecteds assumed infectious; T i (t) ≡ S i (t) + I i (t) denotes the total population size. The constant dispersion coefficients D S and D I measure the probability of dispersion by the susceptible and infective individuals, respectively. Observe that we are using a different notation from what we have used elsewhere, writing variables as a function of t rather than using a subscript for the independent variable in order to avoid needing double subscripts:

where r i is a positive constant.

(a)Using computer explorations, determine whether it is possible to have a globally stable disease-free equilibrium on a patch (without dispersal) where the full system with dispersal has a stable endemic equilibrium. Do you have a conjecture? (b) Using computer explorations determine whether it is possible to have a globally stable endemic equilibrium on a patch (without dispersal) where the full system with dispersal has a stable disease-free equilibrium. Do you have a conjecture? 

When one tries to fit epidemiological data over a long time interval to a model, it is necessary to include births and deaths in the population. Throughout the book we have considered population models with birth and death rates that are constant in time. However, population growth often may be fit better by assuming a linear population model with a time-dependent growth rate, even though this does not have a model-based interpretation. There could be many reasons for variations in birth and death rates; we could not quantify the variations even if we knew all of the reasons. Let r(t) = dN dt /N denote the time-dependent per capita growth rate. To estimate r(t) from linear interpolation of census data, proceed as follows:

1. Let N i and N i+1 be the consecutive census measurements of population size taken at times t i and t i+1 respectively.

Show that in this case,

Use the data of Table 10 .2 to estimate the growth rate r(t) for the population of the USA. Figure 10 .4 shows the time evolution of the USA mortality rate. This mortality rate is fit well by

with μ 0 = 0.01948, μ f = 0.008771, t 1/2 = 1912, and Δ = 16.61. Then the "effective birth rate" b(t) is defined as the real birth rate plus the immigration rate.

Estimate b(t) using r(t) = b(t) − μ(t), with r(t) found in Question 1.

Consider an SEIR disease transmission model. We assume that:

(a) An average infective individual produces β new infections per unit of time when all contacts are with susceptibles but that otherwise, this rate is reduced by the ratio S/N. (b) Individuals in the exposed class E progress to the infective class at the per capita rate k. (c) There is no disease-induced mortality or permanent immunity, and there is a mean infective period of 1/γ.

We define γ = r + μ. The model becomes:

(a) Show that the mean number of secondary infections (belonging to the exposed class) produced by one infective individual in a population of susceptibles is Q 0 = β /γ. (b) Assuming that k and μ are time-independent, show that R 0 is given by Q 0 f , where f = k/(k + μ). What is the epidemiological interpretation of Q 0 f ?

The usual measure of the severity of an epidemic is the incidence of infective cases. The incidence of infective cases is defined as the number of new infective individuals per year. If we take one year as the unit of time, the incidence of infective cases is given approximately by kE. The incidence rate of infective cases per 100,000 population is given approximately by 10 5 kE/N. Tuberculosis (TB) is an example of a disease with an exposed (noninfective) stage. Infective individuals are called active TB cases. Estimated incidence of active TB in the USA was in a growing phase until around 1900 and then experienced a subsequent decline. The incidence rate of active TB exhibited a declining trend from 1850 (See Table 10 .3 and Figure 10 .5). The proportion of exposed individuals who survive the latency period and become infective is f = k k+μ . The number f will be used as a measure of the risk of developing active TB by exposed individuals. (Table 10. 3)?

to active TB depends strongly on the standard of living. An indirect measure of the standard of living can be obtained from the life expectancy at birth. The observed life expectancy for the USA is approximated well by the sigmoid shape function Assume that the risk f varies exactly like life expectancy, that is, assume that f is given by

We refine the model (10.30) by replacing the parameter k by the variable expression μ f (t)/(1 − f (t) and k + μ by μ/(1 − f (t), obtained from the relation f = k/(k + μ). Since the time scale of the disease is much faster than the demographic time scale, the recovery rate r is approximately equal to γ. This gives the model

Simulate TB epidemics starting in 1700 using the model (10.33) with γ = 1 years −1 and β = 10 years −1 , both constant, and with μ(t) given by (10. 29) and f (t) given by (10.32) . Find values of f 0 and f f for which an accurate reproduction of the observed TB trends ( 

Leishmaniasis is a vector-borne disease caused by a protozoan parasites and transmitted by the bite of certain species of sand fly (referred as "vector"). Leishmaniasis is found in many tropical and subtropical countries. The most serious and potentially fatal, if left untreated, form of leishmaniasis is a "visceral" form. The Indian state of Bihar is one of the major foci of visceral leishmaniasis (VL) in the world. Two mathematical models of the spread of VL in Bihar are shown in Figure 10 .7. The models (variables defined in Table 10 .12) incorporate the possibility that individuals seek treatment at private (T ) or public (G) health facilities. The treatment of individuals at private health facilities results in underreporting of cases and deaths in the state, since private practitioners are not required by law to report cases and deaths. Reported-incidence and -mortality data for 2003 can be used to obtain estimates of model-dependent underreporting levels (1− p). Moreover, the transmission coefficients (β in Model I and λ h and λ v in Model II) are also unknown for this disease.

Model I can be written as

The reproduction number for Model I is 

Here, m is the per capita average number of sand flies (assumed constant), C is the mean rate of bites per sand fly, β hv is the transmission "probability" per bite from an infectious sand fly, and Z N v is the proportion of infectious sand flies in the vector population. The proportion of bites of susceptible sand flies on infectious humans is modeled by I N h , that is, sand flies bite the host population at random. The infection rate of susceptible sand flies is

where β vh is the transmission "probability" per infectious human bite by a susceptible sand fly.

Model II is given by

In Model II, the rate of infection for susceptible humans can be modeled by Table 10 .12 for definitions) with λ h ≡ mCβ hv . Here, m is the per capita average number of sand flies (assumed constant), C is the mean rate of bites per sand fly, β hv is the transmission "probability" per bite from an infectious sand fly, and Z N v is the proportion of infectious sand flies in the vector population. The proportion of bites of susceptible sand flies on infectious humans is modeled by I N h , that is, sand flies bite the host population at random. The infection rate of susceptible sand flies is F v (t) = λ v I(t) N h (t) with λ v = Cβ vh where β vh is the transmis-sion "probability" per infectious human bite by a susceptible sand fly. Even though we consider a systematic sub-representation of λ h and λ v , we do not estimate their underrepresentation parameters. That is, we estimate transmission coefficients λ h and λ v as a lumped parameters and do not estimate m, C β hv and β vh explicitly. The reproduction number for Model II is

Estimates of the controlled (since treatment modifies the infectious period) reproduction number (R c ) and the proportion of reported cases (p) are generated from the models. Berkeley Madonna or Matlab (built-in routine lsqcurvefit) curve-fitting tools can be used to estimate model parameters. As an example, parameters estimates for Model I are shown in Table 10 .12. The "best" fit of the model to the cumulative number of reported cases (official data; Table 10 .12) using a least squares fitting procedure is used to generate the unknown parameter estimates. However, for Model II, multiple data (cases and deaths) can be used to estimate parameters. In the models, the cumulative reported cases from time t 0 to t are given by

where k ∈ {1, 2,..., 12}, representing 12 months. Corresponding estimates of C(t k ) (i.e.,C k ) can be found from the data given in Table 10 .12. Hence, the estimation problem, for example in the case of Model I, is then to find optimal values (β and p) of β and p such that Similarly, if mortality data (D k , k ∈ {1, 2,..., 12}) are to be used then the cumulative number of deaths from the model is given by

couple of weeks (it lasts longer in older, younger, and other immunocompromised individuals). A successful infection occurs only in the event that pneumococci are able to spread to and colonize another part of the body, such as the lungs, ear, or blood, etc.

While there are some vaccines routinely used against these infections, it is a lively research area with many open questions. Of the vaccines available, the older one protect against infection only, leaving the colonization phase unaffected, and targets the 23 most common of the over 90 serotypes (analogous to strains of a virus). The somewhat inefficient protection it provides to juvenile and older individuals is not permanent, and children, one of the most affected groups, are unresponsive to the vaccine. Another vaccine, targeting seven common serotypes, is very effective against invasive infection, particularly in children, and may provide protection against colonization. However, whether this aspect of protection is beneficial is under investigation, since it may also provide selective pressure for previously uncommon, yet more invasive, serotypes not targeted by the vaccine. A model of these dynamics was developed in [Sutton, Banks, and Castillo-Chavez (2008)], and is discussed here. See this paper and references therein for more information on the epidemiology of these infections and vaccines developed to target them.

A schematic of pneumococcal infection dynamics including vaccination is shown in Figure 10 .8. The model equations are given by

where the infection rate κ(t) is the oscillatory function κ(t) = κ 0 {1 + κ 1 cos[ω(t − τ)]}.

1. Give interpretations for the parameters φ , ρ, ε, and δ . The simplified model excluding vaccination is given by Many model parameters are not usually available in typical sources such as census data, as is usually the case with epidemiological models, and further are not directly measurable from surveillance data. However, to study prevention strategies, such as vaccination, in a specific population, we need to know all parameter values. In [Sutton, Banks, and Castillo-Chavez (2008) Below is a sample Matlab program to estimate parameters θ = (β , γ) T from "data" (data were actually generated from a forward solution of the model run with the true parameters) in which the number of infections during a week is reported for 14 weeks. The code below fits the sample case notification data using a basic SIR model S = −β SI,

using the initial guess (0.00202, 0.0022) T for parameters (β , γ) T . The data here are without any noise, or error. That is, these data were calculated from solutions of the above SIR model exactly. Of course, there is always some error in actual data, since even a very good model is only an approximation of the epidemiological/physical processes occurring that generate the observations. is the case, suggest what type of data may be more informative for that parameter.

Case notification data to be used to estimate parameters using a simple SIR model, with noise added as denoted in the column headings below: 8. In 2005, the studied population was prescribed a new vaccine nationally, the long-term effects of which were not clear at the time. Vaccine companies publish estimates of the efficacy of released vaccines, but these may vary depending on the ability of the immunized population to mount a response to them. The effect of this vaccine on the colonization susceptibility was unclear. Estimate the vaccine efficacies ε, δ from the data below, assuming unchanged loss of protection ρ and an increased vaccination rate of φ = 0.007984. Data of type d 1 are monthly reported cases of (total) infections, and type d 2 are vaccinated cases. Use the state variable values from the end (last time point) of the "best fit" solution to the answer to the previous question as initial conditions. 9. The reproductive number for the full model with vaccination R φ is given by R φ = β s 0 (ρ + α + δ κ + μ + φ ) + εβ s 0 v (ρ + φ + α + κ μ ) (φ + α + κ + μ)(ρ + α + δ κ + μ) − φ ρ + β s 0 κ(ρ + α + δ κ + μ + δ φ) + εβ s 0 v (κρ + δ κ(φ + α + κ + μ)) (γ + η + μ)[(φ + α + κ + μ)(ρ + α + δ κ + μ) − φ ρ] , (10.37) where κ ≈ κ 0 and s 0 and s 0 v are the proportions of the population unvaccinated and vaccinated, respectively, at the disease-free equilibrium. With the vaccine efficacy parameters just obtained, what can be said about the possibility of effectively vaccinating this population against IPD with the given vaccines?

Use Model II and incidence data in Table 10.12 to estimate λ h

Compute underreporting incidence levels and estimate R II c using model estimates from step

Repeat the first two steps using both data sets (i.e., incidence and mortality data sets in Table 10

Census of India

These most notably include pneumonia, meningitis, and bacteremia. S. pneumoniae, or the pneumococcus, is commonly part of the normal flora of healthy individuals, and is spread through casual contacts via respiratory droplets. Colonization of the nasopharyngeal region is asymptomatic and typically reversed in a References

It is not possible to obtain a good fit of the data of Table 10 .3 to the model (10.30). It is necessary to use a refinement of the model that includes time-dependence in the parameters, and the next step is to describe such a model. 1:n,4) ; t data = 1:numel(data); figure(1); plot(t data,model,t data,data,'.') title('Model fit to weekly cases reported') xlabel('t (weeks)') ylabel('Cases') fprintf('The estimated parameters are beta = %d and gamma = %d',theta hat(1),theta hat (2)) function f=simple SIR(t,y,theta) f=zeros(4,1); f(1)=-theta(1) * y(1) * y(2); f(2)=theta(1) * y(1) * y(2) -theta(2) * y(2); f(3)=theta(2) * y(2); f(4) = theta(1) * y(1) * y(2); function F=obj fcn(p,t,n,IC) [t,z] =ode45(@simple SIR,t,IC,[],p); F = zeros(n,1); F(1:n) = z(2:n+1,4)-z(1:n,4); 6. After running the sample code to estimate the parameter values from the "noiseless" data, use it to estimate the same parameters, starting from the same initial guess (0.00202, 0.0022) T , using the data sets below with an increasing amount of error (1%, 5%, 10%) in the observations. Interpret your results for the parameter estimates as the observational error in the data increases. It may be that one parameter is more reliably estimated from these data than the other. If that