J. Math. Biol. (1992) 30:693-716

dournal of

mathematical 61ology © Springer-Verlag1992

Dynamic models of infectious diseases as regulators of population sizes Jaime Mena-Lorca I and Herbert W. Hetheote 2'*

1 Instituto de Matemfiticas, Universidad Cat61ica de Valparaiso, Casilla 4059, Valparaiso, Chile 2 Department of Mathematics, Universityof Iowa, Iowa City, IA 52242-1466, USA Received April 4, 1991; receivedin revised form July 31, 1991 Abstract. Five SIRS epidemiological models for populations of varying size are considered. The incidences of infection are given by mass action terms involving the number of infectives and either the number of susceptibles or the fraction of the population which is susceptible. When the population dynamics are immigration and deaths, thresholds are found which determine whether the disease dies out or approaches an endemic equilibrium. When the population dynamics are unbalanced births and deaths proportional to the population size, thresholds are found which determine whether the disease dies out or remains endemic and whether the population declines to zero, remains finite or grows exponentially. In these models the persistence of the disease and disease-related deaths can reduce the asymptotic population size or change the asymptotic behavior from exponential growth to exponential decay or approach to an equilibrium population size. Key words: Epidemiological models - Population dynamics - Thresholds - H o p f

bifurcation

1 Introduction

It is generally accepted in ecology that the sizes of plant and animal populations are influenced and constrained by foraging, predation, competition and limited resources. Infectious diseases can also affect population sizes. Anderson and May (1979) cite examples of the impacts of diseases on laboratory, domestic and wild populations of insects, birds and mammals. Indeed, infectious diseases have also influenced human population sizes and historical events (McNeill, 1976). Even today parasitic diseases involving viruses, bacteria, protozoans, helminths and arthropods combined with low nutritional status are the primary reason for differences in the age-specific survival probabilities in developing and developed countries. Measles still causes significant mortality among malnourished children in developing countries. Recently, human immunodeficiency virus (HIV) infec-

* Research supported by Centers for Disease Control contract 200-87-0515. Support services provided at the Universityof Iowa Center for Advanced Studies

694

J. Mena-Lorcaand H. W. Hethcote

tions and the subsequent development of acquired immunodeficiency syndrome (AIDS) has led to increased death rates in various risk groups throughout the world. Thus it is relevant to analyze models for infectious disease transmission with different demographic structures in order to assess the influences of the disease on the population dynamics. Many epidemiological models are formulated so that the infectious disease spreads in a population which either is a fixed closed population or has a fixed size with balancing inflows and outflows due to births and deaths or migration. Results for the simplest epidemiological models of these types are given in Hethcote (1976, 1989). Surveys of results for models with constant size populations are given in Hethcote et al. (1981) and Hethcote and Levin (1989). These models have been developed more extensively partly because they are easier to analyze than variable population size models and partly because they are often realistic for human diseases. The assumption that the population is closed and fixed is often a reasonable approximation when modeling epidemics where the disease spreads quickly in the population and dies out within a short time (about one year). The assumption that the population has fixed size with deaths balancing births is often reasonable when modeling endemic diseases over many years if the population size is constant or nearly constant and the disease rarely causes deaths. However, if the human or animal population growth or decrease is significant or the disease causes enough deaths to influence the population size, then it is not reasonable to assume that the population size is constant. For epidemiological models with varying population sizes, there are many questions to be answered. Does the disease die out or persist in the population? Is there an epidemic or does the disease remain endemic? How does diseaserelated mortality affect the population size and the disease persistence? How does the disease progress in the population as a function of time? Does the disease incidence oscillate periodically forever or does it approach a steady state? This paper answers these questions for five models of the dynamics of infectious diseases with disease-related deaths in populations with demographic dynamics. When formulating the population dynamics part of an epidemiological model with varying population size, one of the simplest ways to incorporate natural births and deaths is to assume that they are proportional to the population size. If N(t) is the total population size as a function of time t, b is the birth rate constant and d is the death rate constant, then dN/dt = (b - d)N,

N(O) = No

(1.1)

is the initial value problem. The net growth rate is r = b - d so that the population size N(t) grows exponentially for r > 0, is constant if r = 0 and decays exponentially if r < 0. This form of population dynamics is called exponential births and deaths. An approach used by Anderson and May (1979) for laboratory populations of mice is constant immigration at rate A and deaths proportional to the population size with rate constant d. Then the initial value problem for N(t) is dN/dt = A - dN,

m(o) = No.

(1.2)

In this case the population approaches A i d as t goes to infinity. This form is called constant immigration with exponential deaths. The models considered here are of SIRS type, which means that susceptible individuals become infectious, then removed with immunity after recovery from

Dynamic models of infectiousdiseases

695

infection and then susceptible again when the temporary immunity fades away. The number of individuals who are susceptible, infectious and removed at time t are denoted by X(t), Y(t) and Z(t), respectively. All individuals are in one of these classes so their sum adds up to the total population size N(t). The incidence in an epidemiological model is the rate at which susceptibles become infectious. If the time unit is days, then the incidence is the number of new infections per day. The daily contact rate 2 is the average number of adequate contacts per day with other individuals of an infective. Time units of weeks, months or years could also be used. Since X(t)/N(t) is the susceptible fraction of the population, 2X(t)/N(t) is the average number of infection transmissions per infective per day. The standard incidence is 2X(t)Y(t)/N(t), which is the average number of infection transmissions by all infectives Y(t) per day. Note that this same incidence results if 2 is the average number of adequate contacts of a susceptible with other people per day since Y(t)/N(t) is the infectious fraction of the population. An alternate way to formulate the incidence is to use the simple mass action incidence flXY where 13 is called the transmission coefficient. When compared to the standard incidence, the simple mass action incidence implies that fl = 2/N, which implies that the contact rate 2 is proportional to the population size. This would be realistic if X, Y and N were densities per unit land area so that doubling the density N would double the daily contact rate 2. The third incidence is a saturation incidence of the form 2XY/(H + X) where H is a constant. When the number X of susceptibles is large compared to H, this incidence is approximately )~Y. This could correspond to a rare disease in a large population of mostly susceptible individuals so that the transmission depends primarily on the number of infectives. This incidence for a population saturated with susceptibles was proposed by May and Anderson (1979). The recovery rate in the epidemiological models here is vY and the rate of loss of immunity is 6Z. The models include a disease-related death rate of infectives given by ~Y. A disease-related death rate of removed individuals (Busenberg and van den Driessche, 1990) is not included here since this would occur very infrequently. The case-related mortality can influence the population dynamics. The linear terms correspond to residence times with negative exponential distributions in the infectious and removed classes (Hethcote et al., 1981). When all transfer rates are considered, the mean (death-adjusted) infectious period is 1/(~ + a + d) and the mean (death-adjusted) period of temporary immunity is 1/(6 + d). For the standard incidence the contact number a = 2/(7 + a + d) is the average number of adequate contacts of an infective during the infectious period. This contact number is the average number of infection transmissions per infective in a completely susceptible population. The two SIRS models considered in Sects. 2 and 3 have population dynamics given by constant immigration with exponential deaths. The incidences are the simple mass action incidence and the standard incidence. The two SIRS models considered in Sects. 4 and 5 also involve the standard and simple mass action incidences, but the population dynamics are given by exponential births and deaths. For all of these models, thresholds are found which determine whether the disease dies out or remains endemic and whether the population size declines to zero, remains finite or grows exponentially. Specific cases where the persistence of the disease and disease-related deaths cause significant changes in the population dynamics are discussed in Sect. 7. The SIRS model considered in Sect. 6 uses the saturation incidence defined previously and exponential births and deaths. May and Anderson (1979) formu-

696

J. Mena-Lorcaand H. W. Hethcote

lated this model and suggested that it could have periodic solutions. Here the existence of periodic solutions is shown by Hopf bifurcation and the precise parameter region where they occur is determined. The SIRS models analyzed mathematically in Sects. 2 and 5 were formulated and applied to diseases in laboratory populations of mice by Anderson and May (1979). The models in Sects. 3 and 4 are different since they use the standard incidence instead of the simple mass action incidence. Much of the mathematical analysis of the five models is from the unpublished thesis of Mena-Lorca (1988). A more general form of the model in Sect. 4 was formulated and analyzed independently by Busenberg and van den Driessche (1990). Other models with varying population size and disease-related mortality were formulated by Anderson and May (1979, 1981), May and Anderson (1978, 1979) and Anderson et al. (1981). Other models with varying population size have been analyzed by Busenberg and Hadeler (1990), Brauer (1990), Pugliese (1990), and Bremerman and Thieme (1989). Some models for HIV and AIDS have varying population size (Jacquez et al., 1988; Castillo-Chavez et al., 1989; Hyman and Stanley, 1988; Anderson et al., 1988; Thieme and CastilloChavez, 1989). 2 The SIRS model with immigration and simple mass action incidence The model here was formulated by Anderson and May (1979) to model two directly transmitted diseases in a laboratory population of mice. The mice showed some loss of immunity to the bacterium, Pasteurella muris, but their immunity to the pox virus, ectromelia, seemed to be lifelong. In the laboratory experiments of Greenwood et al. and Fenner with these diseases, the mice were added at a constant rate so that the population dynamics are immigration and death. The fit of the model here with simple mass action incidence to the data is quite good and the predictions for the model are consistent with the experiments. See Anderson and May (1979) for details. At the end of Sect. 3, we discuss the implications of the data for the suitability of the two models analyzed. In the transfer diagram, A

~X

ltXY

~Y

yY

~Z

6Z

~X

recall that the number of susceptible, infective and removed individuals as a function of time t are X(t), Y(t) and Z(t), respectively and the total population size is N(t). The parameters in this and other models in this paper are: A= d= //= = 7= =

constant immigration rate natural death rate constant transmission coefficient disease-related death rate constant recovery rate constant loss of immunity rate constant.

Note that Roman letters are used for demographic parameters and Greek letters are used for epidemiologic parameters. Assume that d, u and 6 are nonnegative and that A, p, 7 and 6 + d are positive.

Dynamic models of infectious diseases

697

The differential equations for the SIRS model are: X'(t) = A - dX - flXY + 6Z Y'(O = f l x Y - (7 + ~ + d ) Y

(2.1)

Z ' ( t ) = 7 Y -- (6 + d ) Z N'(t) = A -dN

- ~Y

where one of the equations is redundant since N = X + Y + Z. In the absence of disease, the population size approaches the constant size A / d . For simple mass action incidence the contact number (reproduction number) typically is the product of t , a population size and an average infectious period. Here the threshold quantity which determines the asymptotic behavior is the contact number, = f l ( A / a ) / ( 7 + ~ + d).

(2.2)

For the system (2.1) the first octant in X Y Z space is positively invariant. Because N ' ( t ) < 0 for N > A / d , all paths in the first octant approach, enter or stay inside the subset T defined by X + Y + Z 1 the equilibrium P1 is probably also GAS for the SIRS model with 6 > 0, but this Liapunov function does not work. 3 The SIRS model with immigration and the standard incidence

This SIRS model differs from that in Sect. 2 in only one way. The standard incidence 2 X Y / N is used instead of the simple mass action incidence flXY. The transfer diagram is A

AXY/N

>X

7Y

) Y

6Z

)Z

)X

Dynamic models of infectious diseases

699

The differential equations for this model are X'(t) = A - dX - 2XY/N + 6Z Y'(t) = 2XY/N

- (7 + ct + d ) Y

(3.1)

Z ' ( t ) = 7 Y - (3 + d ) Z N'(t) = A - dN - ~ Y

where the equation for N ( t ) follows from N = X + Y + Z. The arguments in Sect. 2 also show that this model is mathematically and epidemiologically well posed. Here the threshold is the contact number defined in Sect. 1 as a = 2/(7 + ct + d),

(3.2)

which is the average number of adequate contacts of an infective during the infectious period. In the absence of disease, the population size N(t) approaches A/d.

For this model the equilibria in X Y Z

space are Po = ( A / d , O, O) and

Pm = (~e, Ye, Ze) where Xe= me/tr Ye = Ne(1 - l/a)/[1 + 7/(6 + d)]

(3.3)

Ze = y Ye/(6 + d)

N e = A / { d + ~(1

-

1/o-)/[1 + ~,/(3 +

a)]}.

As in Sect. 2, P0 is the only equilibrium in the first octant of X Y Z space when a ~< 1. For a > 1, both equilibria are in the first octant. For this SIRS model the asymptotic behaviors of solution paths of (3.1) are the same as for the model in Sect. 2 and are summarized in Table 1. When the disease persists, the equilibrium population size is reduced from A / d to Ne in (3.3) by the disease-related deaths. Another effect of the disease-related death rate constant ct is to decrease the average infectious period slightly from 1/(7 + d) to 1/(7 + ~ + d). This means that the minimum contact rate 2 so that a > 1 and the disease does not die out must be larger when there are disease-related deaths. The proofs of the results in Table 1 for this SIRS model are quite similar to those for the model in Sect. 2. All paths approach the subset T in the first octant in X Y Z space. For a ~< 1, the Liapunov function L = Y can also be used here to show that Po is GAS. For a > 1, the equilibrium (3.3) can be shown to be locally asymptotically stable by linearization. We remark that for ~ = 0, N ( t ) ~ A / d and in the X + Y + Z = A / d plane, the systems (2.1) and (3.1) reduce to a constant population SIRS model so that global stability results follow from those in Hethcote (1976).

4 The SIRS model with exponential births and deaths and the standard incidence The population dynamics for the models in Sects. 2 and 3 do have a natural death process and disease-related deaths, but the birth process is the rather artificial process of constant immigration. Here the immigration is replaced by a natural birth process in which the birth rate is proportional to the population size. Thus the population dynamics without disease is given by (1.1) and is

700

J. Mena-Lorca and H. W. Hethcote

exponential growth or decline of the population except in the special case of a constant size population when the birth and death rate constants b and d are equal. Using the parameters defined in Sect. 2 and b as the natural birth rate constant, the transfer diagram with the standard incidence is bN

2XY/N

~N

?Y

~ Y

6Z

~Z

~X

The differential equations for the model are: X'(t) = bN - dX - 2XY/N Y'(t) = 2XY/N

-

+ ,~Z

(7 + ~ + d ) r

(4.1)

z ' ( o = ~ Y - (~ + d ) Z N'(O = (b - d)N - ~ Y

where one of the equations is redundant since N = X + Y + Z. By arguments similar to those in Sect. 2, this model is well posed in the first octant of X Y Z space. The net growth rate constant in a disease-free population is r = b - d. In the absence of disease, the population size N ( t ) declines exponentially to 0 if r < 0, remains constant if r = 0 and grows exponentially if r > 0. If disease is present, the population still declines to zero if r < 0. For r > 0 the population can go to zero, remain finite or grow exponentially and the disease can die out or persist depending on the values of several threshold quantities. The contact number a = 2/(7 + ~ + d)

(4.2)

is the average number of adequate contacts of an infective during the mean infectious period 1/(7 + ~ + d ) . A closely related threshold quantity is the modified contact number 0 = 2/(y + ~ + b)

(4.3)

where the death rate constant d in o is replaced by the birth rate constant b. When r > 0 and the disease persists, the net growth threshold ¢

_

ro( 7 ~(a---1) l + ~ - - d

)

(4.4)

determines whether the population size N ( t ) decreases exponentially, remains finite or grows exponentially. The net growth threshold ¢ primarily reflects the combined effects of the disease-related death rate constant ~ and the growth rate constant r. The results for the model (4.l) are summarized in Table 2. Although Table 2 may look complicated, it is not. The modified contact number 0 always determines whether the disease dies out (the fraction I--* 0 if 0 ~< 1) or persists (the fraction I - - . I e if 0 > 1). The population size behavior is usually determined by the net growth rate constant r = b - d except for the three cases where the disease persistence affects the population size behavior. When r = 0, the population would normally remain constant, but for ~ > 0 and 0 = a/> 1, the diseaserelated deaths cause the population to become extinct. In the second and third

Dynamic models of infectious diseases

701

Table 2. Summary of stability results for model (4.1) where GAS means globally asymptotically stable, NS means neutrally stable and NUS means neutrally unstable global attractor* in plane S+I+R= 1

asymptotic behavior of N(t)

asymptotic behavior in X Y Z space

N~0 U~O

r=0, c~=0

O=aO

0=a0, 0>1

0~1

* The point (1, 0, 0) is always the attractor on the side where I = 0

cases, r > 0 so the population would normally grow exponentially, but for 0 > 1, the disease persists and the disease-related deaths cause the population to remain finite when t h = l and to become extinct when q ~ < l . N o t e that ¢ ~ < 1 is equivalent to a(l - 1/a)/[1 + 7/(6 + d)]/> r

(4.5)

which displays more clearly the concept that the disease-related deaths defined by overcome the intrinsic net growth rate r of the population and prevent the population from growing exponentially. We n o w briefly outline the a p p r o a c h used to prove the results in Table 2. Results from the unpublished thesis o f M e n a - L o r c a (1988) and f r o m Busenberg and van den Driessche (1990) are reorganized and extended using different notation. The system (4.1) has a special property which makes it easier to analyze. The right sides o f all equations are h o m o g e n e o u s o f degree 1 so that the different equations below for the susceptible fraction S = X / N , the infectious fraction I = Y / N and the removed fraction R = Z / N do not involve N. S'(t) = b - bS - 2SI + 5R + ~IS I'(t) = 2 S I - (y + ct + b)I + ~ I 2

(4.6)

R ' ( t ) = 71 - (5 + b ) R + a i R N ' ( t ) = (r - a I ) N .

The equations for the fractions S, I and R in the disjoint epidemiological states are redundant since S + 1-4- R = 1. Thus the two dimensional system in I R

702

J. Mena-Lorca and H. W. Hethcote

space (or another two dimensional system) can be analyzed separately and then the population size dynamics can be determined from the differential equation for N. The equilibria for the IR system in (4.6) in the triangle T = {(I, R) [I>~ O, R >~0, I + R 0 and 0 --- tr < 1, there is a line of equilibria along the positive X axis in X Y Z space which are neutrally stable and each solution path in the first octant approaches one of these equilibria. When r = 0, ct > 0 and 0 = a > 1, the differential equation for N is a perturbation of N ' = - ~ I e N so that N ~ 0 as t ~ ~ . When r = 0, a > 0 and 0 = a = 1, numerical calculations show that (X, Y, Z)--+(0, 0, 0); in this case the disease dies out slowly enough that it drives the population to extinction. If the modified contact number satisfies 0 ~< 1, then the infectious fraction 1 decreases to zero. The differential equation for N in (4.6) implies that N(t) grows exponentially for r > 0 and 0 ~< 1. See Coddington and Levinson (1955, pp. 314-317) for results on X ' = (A + B(t))X where B ( t ) ~ 0 . The equation for Y in (4.1) can be written as

Y'(t) = [2 -- (y + ~ + a) + 2 ( X / N - 1)]Y

(4.10)

with X / N ~ 1 since 0 ~< 1 so that Y declines exponentially to 0 for tr < 1 and Y grows exponentially for o- > 1. The solution of Z in (4.1) is Z(t) = Zoe-(cS+d)t-l-e-(6+d)t 7

fo

Y(u)e(a+d)" du

(4.11)

Dynamic models of infectious diseases

703

and the exponential nature of Y shows that Z also declines or grows exponentially depending on the contact number o-. For r > 0 and a = 1 there is a line of equilibria defined by S = 1, Z = ? Y / ( 6 + d ) in S Y Z space which is neutrally stable and all paths in the part of the first octant where So ~< 1 approach an equilibrium on this line. In the case when the modified contact number satisfies 0 > 1, we have ! ~ Ie and R ~ Re so that N ' ( t ) = [r - ale - ct(I - I e ) l N

(4.12)

and N grows exponentially with rate constant Q = r - ale. From Ie = (r -- ~)/~ and the equations from (4.6) for the endemic equilibrium, the R e and Se coordinates satisfy ?(r -- e)

Re -

(6 + d + ~)ct' ?+~+d+~ se2

(4.13)

Substituting these into Se + Ie + Re = 1 leads to the quadratic equation QZ + DQ + E = 0

(4.14)

for ¢ where

0"--

D =__

(4.15) - -

~---r

1+

Using the net growth threshold ~b defined in (4.4), q~ < 1 implies that D > 0, E > 0 and D 2 - 4 E > 0 so that both roots of (4.14) are negative and N decreases exponentially to zero with asymptotic rate constant 0 given by the larger root of (4.14). The condition q~ > 1 implies E < 0 so (4.14) has a positive real root Q < r and N ( t ) grows exponentially with growth rate constant 0- If ~b = 1, then Q = 0 and there is a line of neutrally stable equilibria in X Y Z space given by Xe = N / a , Ye = r N / e . Each path in the first octant approaches one of these equilibria.

5 The SIRS model with exponential births and deaths and the simple mass action incidence The SIRS model here is the same as in Sect. 4 except that the simple mass action incidence B X Y is used. Anderson and May (1979) formulated this model and identified the thresholds. Using the parameters defined previously, the transfer diagram is: bN

flXY

~X

~Y

~Y

6Z

~Z

I

.~X

704

J. Mena-Lorca and H. W. Hethcote

The differential equations for this model are X'(t) = bN - dX - flXY + 6Z Y ' ( t ) = f l X r - (~ + ~ + a ) Y

(5.1)

z ' ( o = 7 r - (6 + d ) z

N ' ( t ) = (b - d ) N - ~ Y .

As in Sect. 2 this model is well posed. In the absence of disease the population size N ( t ) decreases exponentially to zero if r = b - . - d < 0, remains constant if r-= 0 and grows exponentially if r > 0. If disease is present, the population still declines to zero if r < 0 since N ' 1 where a 2 = fl/(7 + ~ + d)

490

2 49>1

~(o + d) /

O, r--ff ~,

Xe = (y + ~ + a)/fl

( (r-----Q)7-~ 8+d+o) ]

U--+oo

(X,Y,Z)--+(X~,oo, oo) x~=(~ +~ + a +e)/[~

If r = 0 and ~ > 0, there is a line of equilibria along the X axis. Local stability analysis shows that the equilibrium (Xe, 0, 0) is neutrally stable if O'2ff~e < 1 and is neutrally unstable if a2Xe > 1 where o 2 = fl/(]) AV(~ "1-d ) . Using N as a Liapunov function, N ' = - a Y e < 0 so the Liapunov-Lasalle theory stated in Sect. 2 implies that solution paths approach the largest positively invariant subset of the Y = 0 plane. In that plane X ' = ( b + 6 ) Z and Z ' = - (b + ~)Z so paths approach the X axis along lines X + Z = C. Thus all solutions with Yo > 0 approach an equilibrium in the line segment where 0"2Xe < 1. For r > 0 the extinction equilibrium (0, 0, 0) is unstable with repulsion along the X axis. An endemic equilibrium is xe = (~, + ~ + g ) / ~ Ye = rX~/(~( 1 -- (a))

(5.6)

Ze = ]~ r e / ( 3 -}- d ) N e = Xe/(1 -- q~)

where the net growth threshold q5 is given by (5.2). This endemic equilibrium is in the first octant of X Y Z space iff q~ < 1. In this case when q5 < 1, the characteristic equation for the endemic equilibrium is W3-k-azW2Jr-alw +%=0

(5.7)

706

J. Mena-Lorca and H. W. Hethcote

where a2=6+d+flYe-r>O al = flYe(flXe -- b) + (flYe - r)(6 + d)

(5.8)

~0 = f l Y e [ ( f l X e -- b)(6 q- d ) - y(r + 6 4- d)]

= r(6 + d ) ( 7 + ~ + d ) > O.

From these equations, we find a2al - ao = (flYe -- r)[(6 + d) 2 + flYe(flXe - b) + (3 + d ) ( f l Y e - r)] + j~YeT(6 + d + r) > 0

since B Y e - r and B X e - b are positive if ~b < 1. Thus the Routh-Hurwitz criteria (Brauer and Nohel, 1969, p. 155) are satisfied if ~b < 1 so the eigenvalues have negative real parts and the endemic equilibrium is locally asymptotically stable. If ~b = 1, then

(s, L R) --, ( 0 , ,r ~(6r? ~ d) ) '

(5.9)

X approaches Are given by (5.6) and Y, Z and N approach infinity. Note that ¢ ~< 1 is equivalent to ~/(1 + ?/(6 + d)) ~> r.

(5.10)

This inequality corresponds to the disease-related deaths being large enough that they prevent the population from growing exponentially with rate constant r as it would in the absence of disease. Thus for O ~< 1, the presence of disease changes the population dynamics significantly. If O > 1, then the disease-related death rate constant ~ is not sufficiently large to overcome the exponential growth, but it does reduce the growth rate constant from r. Let the growth rate constant in the equation for N ' ( t ) be Q = r - ~Ie so that Ie = (r -- 4)/~, (r - Q)y R e -- ~(6 -]- d --}- Q) '

(5.11)

Se=0. The Re expression follows by setting the right side of the R ' ( t ) equation in (5.3) to zero. Also Se = 0 since otherwise S e N would approach infinity so that I(t) would approach infinity by the equation for I'(t) in (5.3). This same equation implies

xe =(~ +~

+a+o)/~.

(5.12)

From Ie + R e = 1, we find that the quadratic equation for the growth rate constant Q is o2+Be+C=0 B=~+6+d+~-r, c = ~(6 + a)(1 - 4J).

(5.13)

707

Dynamic models of infectious diseases

Since B 2 - 4 C > 0 , the quadratic equation (5.13) has one negative and one positive root Q+ if ~b > 1. Also Q+ ~ 0 + as q5 ~ 1+ and Q+ ~ r - as :t ~ 0 +. The root Q+ is the asymptotic growth rate constant when q~ > 1 so N ( t ) has exponential growth asymptotically. In summary for tk > 1, the disease-related deaths reduce the exponential growth rate constant to Q and the limiting susceptible fraction is zero instead of being positive. 6 The SIRS model with saturation incidence

The unique feature of this SIRS model is the saturation incidence of the form 2 X Y / ( H + X ) . When X is much larger than H so the population is saturated with susceptibles, the incidence is approximately 2 Y which depends only on the number of infectives. This saturation incidence was proposed by May and Anderson (1979) as an example which could lead numerically to stable limit cycles. May and Anderson (1978) considered similar models where this saturation incidence reflects the proportion of parasites which are successful in entering a new host. They proposed this saturation incidence in an SIRS model as "the collapsed dynamics of a free living infective state." Using the notation and parameters defined previously, the transfer diagram for the SIRS model with the saturation incidence and exponential births and deaths is bN

2XY/(H + X)

~X

,Y

7Y

6Z

~Z

,X

The differential equations for this model are X'(t) = bN - dX - 2XY/(H + X) + 6Z Y ' ( t ) = 2 X Y / ( H + X ) - (y + ~ + d)Y

(6.1)

z ' ( t ) = ~, Y - (,~ + d ) Z

N ' ( t ) = (b - d ) N - ~ Y.

As in Sect. 2 this model is well posed and the first octant in X Y Z space is positively invariant. The goal is to see how the parameters b, d, ~, c5, ?, )~, and H affect the dynamics of the disease and to find the regions in parameter space where the model has stable limit cycles. Hopf bifurcation is used to determine where periodic solutions appear. In this model the contact number a, modified contact number 0 and net growth threshold ¢ are a = 2/(y + ~ + d),

(6.2)

0 = 2/(y + e + b),

(6.3)

- - -r ( 1 + ~ -?- - ~ ) .

(6.4)

Note that the parameter H could be scaled out of Eqs. (6.1) by changing to new independent variables which are the old variables divided by H. We decided not to do this in order to show explicitly the effects of the parameter H.

708

J. Mena-Lorca and H. W. Hethcote The finite isolated equilibria of (6.1) in X Y Z space are (0, 0, 0) with N = 0

and X e = H / ( a - 1) Ye = rNe/a

(6.5)

Ne=

H (G -

1)(1 - ¢ )

The extinction equilibrium (0, O, 0) is easily found to be locally asymptotically stable if r = b - d < 0 and unstable if r > 0. The endemic equilibrium given by (6.5) is in the first octant only if o- > 1, r > 0, ~ > 0 and ¢ < 1. In this model this endemic equilibrium is locally asymptotically stable for some parameter values and unstable for others. Later in this section we show that when the stability switch occurs, a stable periodic solution arises by H o p f bifurcation. The differential equations for the fractions S = X / N , I = Y / N and R = Z / N which are susceptible, infectious and removed are: S'(t) = b - b S I'(t) -

2SIN H+ SN

2SIN - + r R + ~IS H + SN (7 + ~ + b)I + aI 2

(6.6)

R ' ( t ) = ?I - (6 + b ) R + a i R

N'(t) = (r - aI)N.

The results for this SIRS model (6.1) are summarized in Table 4 and proved below. I f r < 0, then N" ~ r N so that the population always goes to zero and the extinction equilibrium (0, 0, 0) is globally asymptotically stable. Moreover, 2 S N -] I'(t) 0 , Y > ~ O , X + r < . U o } . Theorem 6.1 Consider system (6.1) with r = 0 , a = 0 and the contact number a = 2/(~ + d). I f a 1, then D - {(X, 0) I 0 ~< X ~< No}

Dynamic models Of infectious diseases

709

Table 4. Summary of stability results for model (6.1) where NS means neutrally stable and NUS means neutrally unstable asymptotic behavior of N(t)

N~O

r 1

N = No

line (No, 0, 0) is NUS line (Xe, Ye, Ze) is NS line (X~, 0, 0) is NS if (or - 1)X¢< H line (Xe, 0, 0) is NUS if (tr - 1)X~> H

r = 0, >0 r>0, ~r < 1 (implies 0 < 1) r>0, cr = 1 (implies 0 < 1) r>0,

N ~ Ne

N~

(S, LR)~(1,0, O) (X, Y, Z) ~ ( ~ , 0, 0)

N~

(S, LR)~(1,0, O) (X, Y, Z) ~ ( ~ , Y~,Z~)

~b>l,01,0>1

N~

(X, Y,Z)--*(~, ~ , oo) (S,I,R)-'~(O,[e, Re)

~b 1

(x, Y, z) --, (Xe, ~, ~) (S,I,R)~(S~,I~,Re) (Xe, Ye, Z~) is asy. stable if 2 > 2* (X~, Y~, Ze) is unstable if 2 < 2* with a surrounding stable, periodic solution

is an asymptotic stability region for the endemic equilibrium H

+

These isolated equilibrium points in X Y space where Z = N o - X - Y become lines o f neutrally stable (NS) or neutrally unstable ( N U S ) equilibrium points in the first octant in X Y Z space. I f r = 0 and ~ > 0, then there is a line o f equilibria (Xe, 0, 0) along the X axis which are neutrally stable if (o- - 1)X~ < H and neutrally unstable (with repulsive Y direction) if ( a - 1)Xe > H. As in Sect. 5 the L i a p u n o v function N with N ' = - ~ Y implies that solution paths a p p r o a c h the largest invariant subset o f the plane Y = 0. Since Z" = - ( 6 + d ) Z and X ' + Z ' = 0 in the plane Y = 0, all paths in this plane a p p r o a c h the X axis. Since N must decrease as long as Y is positive and paths starting with Yo > 0 cannot a p p r o a c h the neutrally unstable points on the X axis, they must a p p r o a c h one o f the neutrally stable equilibria on the X axis. I f r > 0 and a < 1, then

Y'(t)=

)a~+X

(7+~+d)

Y 0 and a = 1, then numerical calculations show that (S, L R) 4 ( 1 , 0, 0) and (X, Y, Z ) --, (oo, Y, Z ) where Y and Z depend on the initial conditions. If r > 0, a > 0, and ~b > 1, the model (6.6) is difficult to analyze rigorously because the S, /, R equations do not separate from the N equation and the disease-related deaths do not overcome the natural growth rate constant so the population still grows exponentially. The equation for I in (6.6) can be rewritten as

I'(t) = [2X/(H + X ) - (7 + ~ + b) + ~I]I

(6.10)

which suggests that at an equilibrium, either I = 0 or ]) "Jr O~ -[- b -- ~L .

(6.11)

I'(t) 1, X ~ Xe and S ~ 0 so /e + Re = 1 and Re =y/e/(6 + b - ~L). This implies that /e is the one root in (0, 1) of the quadratic equation ~I2-(6 +b +7 +~)I+6

+b =0.

(6.13)

In this 0 > ! case, N, Y and Z grow asymptotically exponentially with rate constant O = r - ~I~. Numerical calculations have verified that the local results above also give the global behavior of solutions. Finally, we consider the most interesting case when r > 0, ~ > 0, a > 1 and q~ < l so the finite endemic equilibrium point (X~, Ye, Ze) given by (6.5) is in the first octant. The Jacobian at (Xe, Y~, Ze) is [r-LL

d+r-Fo

6+d+r] (6.14)

0

~

-(6 + d)l

where L = (relXe)r(1 -

1/a) =

r (a -- 1)

F

a

( 1 - ~o)'

F=7+~+d.

(6.15)

The characteristic equation of the linearization at (Xe, Ye, Ze) is P(w) = w 3 + a 2 w 2 -}- al w + ao

(6.16)

where a2=L +6 +d-r, ax = L ( F + 6 - r) -- r(6 + d), ao = --L[7(6 + d + r) + (6 + d)(d + r - r)] = rr(6 + d)( 1 - 1/a).

(6.17)

Dynamic models of infectious diseases

711

The R o u t h - H u r w i t z criteria (Brauer and Nohel, 1969) for all roots of (6.16) to have negative real parts are a2 > 0, ao > 0 and a2al - a0 > 0. H o p f bifurcation occurs when one root is negative and two complex conjugate roots have real parts which change sign. If this occurs, then these roots are purely imaginary on a neutral stability surface in parameter space. It is easy to verify that (6.16) has one negative root and two conjugate purely imaginary roots iff a2 > O,

a2al - a0 = 0.

ao > O,

(6.18)

We examine the dependence of these conditions on the parameters in order to determine where H o p f bifurcation might occur and where the equilibrium point is stable and unstable. There are several hypotheses which must be verified in a H o p f bifurcation theorem in order to guarantee the existence of periodic solutions. The H o p f bifurcation theorem in Hassard et al. (1981, Theorem II) is used here because it has associated efficient algorithms in the BIFOR program for computing the form of the bifurcating periodic solutions, their direction, their periods and their stability. The first step is to determine neutral stability surfaces in parameter space where H o p f bifurcation might occur by finding where the conditions (6.18) are satisfied so that two characteristic roots are purely imaginary. Our conditions r > 0, a > 0, a > 1 and 49 < 1 guarantee that ao > 0 and that L is positive so a2 is positive if 6 +d-r

>0.

(6.19)

This condition is assumed for algebraic convenience since a2 > 0 could still be satisfied when 6 + d - r has a small negative value. The third condition in (6.18) can be written as 0 = a2al - ao = P ( L ) = A L 2 + B L + C

A = (6 + d B = (6 + d C = -(6

r) + 7 + ~ 0 2 + (6 + d ) 7 - ra

(6.20)

+ d - r)(6 + d)r.

It is convenient to choose 2 as the bifurcation parameter since L depends on 2, but A, B and C do not. Thus we look for a 2* which depends on the other parameters such that P(L(2*)) = 0. When r > 0, a > 0, a = 2 / F > 1 and 49 < 1, then L(2) given by (6.15) increases monotonically from 0 at 2 = F to Loo = (r / ~ ) F /(1 - 49).

(6.21)

The condition 49 < 1 implies ~ > r so A in (6.20) is positive and the quadratic P ( L ) is concave upward. The assumption (6.19) implies that P ( 0 ) = C < 0. If P ( L ~ ) > 0, then by the intermediate value theorem, there is a unique L* between 0 and L~ such that P ( L * ) = 0. Since L(2) is monotonic increasing from 0 to L ~ , the intermediate value theorem implies that there is a unique 2* in (F, oo) such that L(2*) = L* and hence P(L(2*)) = 0. Thus the characteristic equation (6.16) has purely imaginary roots if P ( L ~ ) is positive. Using (6.21) and D = 6 + d, rD P(Lo~) = [aD -- r(7 + D)] 2

{F2rDA+ x

[BF -

[aD - r(7 + D)] (D - r)(aD - r(7 + D)]}.

(6.22)

712

J. Mena-Lorea and H. W. Hethcote

Since A is positive and ~b < 1, P(Loo) is positive if the last factor in brackets is positive. If > ~,

(6.23)

then (6.24)

B > (D - r)(D - r + ~)

and the last factor is greater than (D - r)[(7 + d ) ( D - r + ~) + ~(~ - r) + r(y + D)]

(6.25)

which is positive since q~ < 1 implies ~ > r. Note that the condition (6.23) is convenient since it greatly simplifies the analysis. Using (6.15) and the quadratic formula on (6.20), the explicit formula for the neutral stability surface 2" in terms of the other five parameters in six parameter space is 2* = 2 A F 2 / { 2 A F - [ - - B + x / B 2 -- 4AC]ct(1 - qS)/r}

(6.26)

where A, B and C are given by (6.20). Let the characteristic roots of (6.16) be # ( 2 ) _ iv(A) and - q ( 2 ) . Since the conditions (6.18) are satisfied, these roots are a pair of purely imaginary roots _+/v(2*) and a negative root -~/(2") on the neutral stability surface defined by (6.26). Moreover, q(2*)vz(2*) = a0(2*) which is positive so the complex part v is nonzero at 2". The only remaining hypothesis for the H o p f bifurcation theorem is the transversality condition #'(2*) # 0. Using/a + / v as a root of (6.16) and separating real and imaginary parts yields ~ 3 _ 3/iV 2

"k-a2/~2 -- a2v 2 + al# + ao = 0

(6.27)

3/x2v -- v 3 + 2a2#v + a i r = 0.

(6.28)

The partial derivatives of the coefficients given by (6.7) are Oa2/O2 = t~L /O2 = FLoo/2 2

(6.29)

8 a , / 0 2 = ( F + 6 - r)FLoo/22

r)]rLoo/2 2.

0ao/82 = - [~,(6 + d - r) + (6 + d)(r + d -

Taking derivatives of (6.27) with respect to 2, we obtain cll 0/z/02 + c12 Ov/O2 = bl where ell =

3~ 2 --

3v 2 + 2a2/.~ --~ al,

c12 = - 6/~v - 2a2 v, bl = [(r + d - F)(6 + d) + y(r + 6 + d)

--

[~2 + V2 _ _

Taking derivatives of (6.28) leads to C21 ~#/02 -[- C22 Or/02 = 62

(F

--

r + 6 ) # ] F L ~ / 2 z.

Dynamic models of infectious diseases

713

where c21 = 6#v + 2a2 v

c= = 3# 2 - 3v 2 + 2a2# + al b2 = [ - 2 # v - (C - r +

(~)v]I'too/2 2.

The denominator matrix in the Cramer's rule solution of the equations above is (a 1 - - 3 v 2 ) 2 q- 4a2 v2 which is positive. Hence the system has a nonzero solution O#/d2 at 2* if the numerator is nonzero. This is

b,1 0 and 7 > c~ are not satisfied. Thus the algebraic calculations are useful since they lead to the explicit formula (6.26) and the numerical calculations are valuable in verifying this formula and establishing its universal applicability. For all of the parameter sets used in the calculations with BIFOR, the bifurcation occurs for 2 < 2* and the bifurcating periodic solutions are stable.

714

J. Mena-Lorcaand H. W. Hethcote

7 Discussion

The SIRS models in Sects. 2 and 3 use the simple mass action incidence and the standard incidence, respectively. Both models have population dynamics corresponding to immigration and exponential deaths so that the populations always persist and remain finite. For these models the effect of the disease on the population dynamics is simple and direct. When the disease remains endemic, the disease-related deaths decrease the equilibrium population size. Although the asymptotic behaviors for both models are summarized in the same way in Table 1, these models are actually quite different. The contact numbers a given in (2.2) and (3.2) are different and the coordinates of the endemic equilibria given by (2.3) and (3.3) are different. As described at the beginning of Sect. 2, Anderson and May (1979) used the SIRS model with simple mass action incidence ~XY to fit data on two diseases in laboratory populations of mice. They found that this model did fit the data and that the predictions of this model were consistent with the experimental data. As described in Sect. 1 the incidence ~XY would seem to be most appropriate for situations where N is the density per unit land area or per unit volume so doubling N would double the density and, consequently, might also double the disease incidence. Thus it might be realistic in a room with mice since doubling the number of mice might double the contact rate and incidence. However, if the mice were infected primarily by their neighbors in nearby cages, then doubling the number of mice in the room might have no effect on the infection rate. Indeed, the description of the experiments states that the space available to the mice was adjusted to keep the population density constant as the absolute levels changed. This suggests that the model in Sect. 3 with the standard incidence 2XY/N might give an equally good or better fit to the data. Could the model with the standard incidence also be used for the data on these diseases in laboratory populations of mice? In both models the graphs of Ne as functions of the disease-related death rate constant ~ show depression for intermediate values of 6. This is reasonable since low values of ~ would have a small effect while high values of ~ would cause the disease to die out since tr decreases as ~ increases. Thus both models are consistent with the observed equilibrium population sizes as a function of ~. For immigration rate A between 0.33 and 6 new mice per day, the observed equilibrium mouse population size data points lie between 20 and 210 and as a function of A, they fall on a straight line with positive intercept of about 15. Expression (2.3) does have Are as a linear function of A with positive intercept, but (3.3) has Are as a linear function of A with zero intercept. This suggests that the model in Sect. 3 with the standard incidence is not as good since it is not consistent with this aspect of the data. Thus the simple mass action incidence appears to be better for the two diseases in laboratory population of mice. Which form of incidence is best for human diseases? If the infection incidence is of the form 2NaXY, then data in Anderson (1982a, p. 157) leads to values of a between -0.93 and -0.97 for five human diseases in community sizes between 1,000 and 400,000. This strongly suggests that the standard incidence corresponding to a = - 1 is much more realistic for human diseases than the simple mass action incidence pXY corresponding to a = 0. In basic SIRS models the contact number is the daily contact rate 2 times the average infectious period. Table 1.4 in Anderson (1982b) shows that for diseases such as measles, whooping cough, chickenpox, diphtheria and scarlet fever, the contact numbers (his

Dynamic models of infectious diseases

715

basic reproductive rates) are quite consistent in various areas and do not vary for a disease by more than a factor of 1.5 from rural to urban areas or from small to large communities. Anderson concludes that the contact numbers (his basic reproductive rates) are only weakly dependent on community size. Schenzle and Dietz (1986) also conclude that contact rates and reproductive numbers are essentially independent of population sizes for human diseases. Thus for human diseases the standard incidence which has the contact rate and the contact number independent of population size seems better than the simple mass action incidence. The SIRS models in Sects. 4, 5 and 6 all have population dynamics given by exponential births and deaths, but they use the standard, simple mass action and saturation incidences. For these models the disease affects the population dynamics dramatically. For each model a net growth threshold ~ is defined which involves r, ~ and other parameters. The inequality q~ < 1 corresponds to the disease-related deaths being sufficiently large to overcome the natural population growth. When r > 0 so that the populations without disease would grow exponentially, the populations have finite equilibrium size or become extinct when ~b < 1 and the disease persists. When r > 0 and ~b > 1, the population still grows exponentially but with a reduced growth rate constant when the disease persists. Details for the three SIRS models are given in Tables 2, 3 and 4. Periodicity in the contact rate can lead to periodic solutions in epidemiological models (Hethcote and Levin, 1989). It is less common for periodic solutions to occur when the coefficients are not periodic, but this does occur for some nonlinear incidences (Liu et al., 1986, 1987; Hethcote and van den Driessche, 1991). The SIRS model in Sect. 6 with the saturation incidence can have asymptotically stable periodic solutions even though there is no periodicity in the contact rate 2. Thus this saturation incidence is another form of nonlinear incidence which leads to periodic solutions. Acknowledgements. We thank Linda Gao and Jinshi Zhou for helpful discussions and suggestions. References Anderson, R. M.: Transmission dynamics and control of infectious disease agents. In: Anderson, R. M., May, R. M. (eds.) Population Biology of Infectious Diseases, pp. 149-176. New York Heidelberg Berlin: Springer 1982a Anderson, R. M.: Directly transmitted viral and bacterial infections of man. In: Anderson, R. M. (ed.) Population Dynamics of Infectious Diseases, pp. 1-37. London New York: Chapman and Hall 1982b Anderson, R. M., Jackson, H. C., May, R. M., Smith, A. D. M.: Population dynamics of fox rabies in Europe. Nature 289, 765-777 (1981) Anderson, R. M., May, R. M.: Population biology of infectious diseases I. Nature 180, 361-367 (1979) Anderson, R. M., May, R. M.: The population dynamics of microparasites and their invertebrates host. Philos. Trans. R. Soc. Lond. B291, 451-524 (1981) Anderson, R. M., May R. M., McLean, A. R.: Possible demographic consequences of AIDS in developing countries. Nature 332, 228-234 (1988) Brauer, F., Nobel, J. A.: Qualitative theory of ordinary differential equations. New York: Benjamin 1969 Brauer, F.: Models for the spread of universally fatal diseases. J. Math. Biol. 28, 451-462 (1990) Bremermann, H. J., Thieme, H. R.: A competitive exclusion principle for pathogen virulence. J. Math. Biol. 27, 179-190 (1989)

716

J. Mena-Lorca and H. W. Hethcote

Busenberg, S. N., van den Driessche, P.: Analysis of a disease transmission model in a population with varying size. J. Math. Biol. 28, 257-270 (1990) Busenberg, S. N., Hadeler, K. P.: Demography and epidemics. Math. Biosci. 101, 41-62 (1990) Castillo-Chavez, C. C., Cooke, K. L., Huang, L., Levin, S. A.: On the role of long periods of infectiousness in the dynamics of AIDS, Part 1, Single population models. J. Math. Biol. 27, 373-398 (1989) Coddington, E. A., Levinson, N.: Theory of ordinary differential equations. New York: McGrawHill 1955 Hale, J. K.: Ordinary differential equations. New York: Wiley-Interscienee 1969 Hassard, B. D., Kazarinoff, N. D., Wan, Y-H.: Theory and applications of Hopf Bifurcation. (Lond. Math. Soc. Lect. Note Ser.) Cambridge London: Cambridge University Press 1981 Hethcote, H. W.: Qualitative analysis for communicable disease models. Math. Biosci. 28, 335-356 (1976) Hethcote, H. W.: Three basic epidemiological models. In: Gross, L., Hallam, T. G., Levin, S. A. (eds.) Applied Mathematical Ecology, pp. 119-144. Berlin Heidelberg New York: Springer 1989 Hethcote, H. W., Levin, S. A.: Periodicity in epidemiological models. In: Gross, L., Hallam, T. G., Levin, S. A. (eds.) Applied Mathematical Ecology, pp. 193-211. Berlin Heidelberg New York: Springer 1989 Hethcote, H. W., Stech, H. W., van den Driessche, P.: Periodicity and stability in epidemic models: a survey. In: Busenberg, S. N., Cooke, K. L. (eds.) Differential Equations and Applications in Ecology, Epidemic and Populations Problems, pp. 65-82. New York: Academic Press 1981 Hethcote, H. W., van den Driessche, P.: Some epidemiological models with nonlinear incidence. J. Math. Biol. 29, 271-287 (1991) Hyman, J. M., Stanley, E. A.: Using mathematical models to understand the AIDS epidemic. Math. Biosci. 90, 415-473 (1988) Jacquez, J. A., Simon, C. P., Koopman, J., Sattenspiel, L., Perry, T.: Modeling and analyzing HIV transmission: The effect of contact patterns. Math. Biosci. 92, 119-199 (1988) Liu, W. M., Hethcote, H. W., Levin, S. A.: Dynamical behavior of epidemiological models with nonlinear incidence rates. J. Math. Biol. 25, 359-380 (1987) Liu, W. M., Levin, S. A., Iwasa, Y.: Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models. J. Math. Biol. 23, 187-204 (1986) May, R. M., Anderson, R. M.: Regulation and stability of host-parasite population interactions. II. Destabilizing processes. J. Anim. Ecol. 47, 219-267 (1978) May, R. M., Anderson, R. M.: Population biology of infectious diseases I1. Nature 280, 455-461 (1979) McNeill, W. H.: Plagues and People. Oxford: Blackwell 1976 Mena-Lorca, J.: Periodicity and stability in epidemiological models with disease-related deaths. Ph.D. Thesis: University Of Iowa 1988 Pugliese, A.: Population models for diseases with no recovery. J. Math. Biol. 28, 65-82 (1990) Schenzle, D., Dietz, K.: Critical population sizes for endemic virus transmission. In: Fricke, W., Hinz, E. (eds.) Ra/imliche persistenz und diffusion yon krankheiten. (Heidelb. Geogr. Arbeiten, vol. 83) Heidelberg: 1986 Thieme, H. R., Castillo-Chavez, C.: On the role of variable infectivity in the dynamics of the human immunodeficiency virus. In: Castillo-Chavez, C. (ed.) Mathematical and statistical approaches to AIDS epidemiology. (Lect. Notes Biomath., vol. 83, pp. 157-176) Berlin Heidelberg New York: Springer 1989