Deterministic Modeling of Infectious Diseases: Measles Cycles and the Role of Births and Vaccination
H Trottier, P Philippe
Keywords
chaos, damped oscillations, deterministic modeling, dynamics, infectious diseases, random variations
Citation
H Trottier, P Philippe. Deterministic Modeling of Infectious Diseases: Measles Cycles and the Role of Births and Vaccination. The Internet Journal of Infectious Diseases. 2002 Volume 2 Number 2.
Abstract
This paper is third in our series on deterministic modeling (DM) of infectious diseases. The first paper dealt with theories and methods of DM. The second provided sensitivity analyses based on the SEIR model and showed how changes in the transition rate from susceptibles to infectious (force of infection) and variations in the latent and infectiousness periods can alter the dynamics of a population with no births and deaths. The current paper aims at providing a more realistic portrait of the infection by introducing births and deaths in the SEIR model, vaccination, and finally time variation in the contact rate. This paper shows how the recruitment of susceptibles by births yields the well-known cycling of epidemics. We also show how vaccination can alter the population pattern of infections in the same way that the birth rate does it. Varying the birth or vaccination rates is however insufficient to avoid tapering off of epidemics with time. Since the tapering off is at variance with real data, seasonality in the contact/transmission rate was added to the equations. Including seasonal variations introduced chaos in the pattern of population infection. The resulting sensitivity to demographic data is sufficient to explain most epidemic patterns at least from a qualitative standpoint. We conclude this series of papers on population patterns of infections by urging for efforts to understand how complex work and evolve since this provides a key to appropriate public-health interventions. The population dynamic approach can help control many types of population disease patterns.
Introduction
This paper is an application of the theory and methods of deterministic modeling (DM) to common childhood diseases. It is the last in a series of three papers dealing with DM. The first paper (1) was concerned with the theory and methods of DM. The second paper (2) looked into model sensitivity and showed how changing the transition rate from susceptibles to infectious or varying key parameters of the SEIR model could alter the profile of an epidemic for a population of constant size with no new entries from births (new susceptibles) and no mortality. Albeit simple, this model showed that the reproduction number (R) of the infection is the rate-limiting step of epidemics; in particular, R must be greater than 1 for epidemics to thrive. In the particular case of a constant population size, no additional disease case arises in the population after a certain time because no new susceptible individual enters the population. It follows that once individuals have had the disease, they become immune and the epidemic dies out.
The current paper aims at providing a better understanding of the dynamics of epidemics by allowing for entries of new susceptibles (by births) and mortality in the course of time. Altering the differential equations to deal with births and deaths changes the population dynamics because new individuals can get diseased continuously. In this case, multiple epidemics are possible. We, therefore, show that recruiting new susceptibles by births generate the well-known cycling (oscillations) of most epidemics. The frequency of cycling over time will thus depend on the birth rate; elevated birth rates will entice tight cycling over time, and low birth rates loose cycling. Further, if repetitive epidemics separated by more or less short inter-epidemic intervals can be elicited by adding births to the susceptibles compartment of the model, then vaccination will have similar bearing on epidemics through its depleting of the number of susceptibles. Indeed, effective vaccination acts as though the number of susceptibles were reduced; this will therefore affect the contact rate and the reproduction number R. Therefore, natality and vaccination have common effects (3).
Epidemics in large unvaccinated populations often allows for a very clear pattern of recurrence every 2 to 3 years (4). For example, glancing over the reported cases of measles in New York from 1928 to 1971 allows one to clearly see a regular pattern (Chart 1). From 1945 until the widespread use of the vaccine in 1963, outbreaks of measles occurred more or less every two years (5).
This regular recurrence of the infection can be rescued by deterministic models by introducing new susceptibles (birth rate) in the SEIR model. In the second article of this series, we saw that if R is higher than 1, the number of cases increases to reach a peak and then decreases. Consequently, no further case arises out of the population after a certain time because all susceptible individuals have been processed to the recovereds compartment. There is, therefore, only one epidemic. The situation is different if we allow for entries of new susceptibles in the population in the course of time; new epidemics are possible because the new entries can turn over to disease cases. We show in the following that damped oscillations over time arise in this context. Once the first few epidemic peaks have happened the occurrence of disease tends to decrease until an endemic situation finally settles in. As we shall see later, the dampening of oscillations is at variance with real data sets; therefore, the equations of the population dynamic will need to be reparametrized to insure the regularity of epidemics; the parameter at stake to embed into equations is the seasonal variation of the contact rate.
A Model Of Measles Including Birth And Death Rates
An appropriate model for measles is the SEIR model (2). The latter model can be adapted to deal with births and deaths in the following way:
where:
S (t) = number of susceptibles at time t
E (t) = number of exposeds (infecteds) but not yet infectious at time t
I (t) = number of infectious at time t
R(t) = number of recovereds at time t
λ = the rate (force) of infection per unit time (the time unit is an interval from t to t+δ where δ is a very small)
f = the rate at which an infected individual becomes infectious per unit time
r = the rate at which an infectious individual recovers per unit time.
b= birth rate per unit time
m= death rate per unit time
Assuming a starting population of size N and a birth rate b per unit time, the number of new susceptibles per unit time is given by b*N. Moreover, m*S(t), m*E(t), m*I(t) and m*R(t) reflect the number of susceptibles, infecteds, infectious, and immunizeds, respectively, who die per unit time. Thus, in this model the rate of change per unit time in the number of susceptibles, infecteds, infectious, and immunizeds is given by the following differential equations:
dS/dt = λ*S(t) + b*N – m*S(t)
dE/dt = λ*S(t) – f*E(t) – m*E(t)
dI/dt = f*E(t) – r*I(t) – m*I(t)
dR/dt = r*I(t) – m*R(t)
The transition rate per unit time (day) for measles (see second paper (2) for details) is therefore (for a population 1 000 000 individuals):
The daily death rate can be estimated by the life expectancy, which is around 78 years for Canada. The daily death rate can be estimated by the inverse of the number of days that an individual can expect to live (i.e., 365 days * 78 years). We suppose that the population is stationary, i.e., the birth rate is equal to the death rate. The daily death and birth rates are therefore:
The simulations below were undertaken with software Model Maker (Family Genetix). Only simulated series of infected cases are presented.
Results
Birth rate
Let us now peruse the dynamics of measles after the introduction of one infectious individual in a healthy unvaccinated stationary population of 1M individuals (R=15, latent period= 8 days, period of infectivity=7 days, life expectancy=78 years, and birth rate=death rate as above). To better fix ideas, let us recall that R=15 is the measles reproduction number for a completely unvaccinated population. Chart 2 shows that the introduction of only one infectious immigrant at time zero produces epidemic cycles for many years to come. This shows that the signal characterizing the introduction of an infected in the population has resonance that dies away only slowly in time. Further, once an infectious case has entered the population (given R>1), further recruitment of new susceptibles by births is sufficient to elicit oscillations over time. Chart 2 time scale displays periods of 10 years; thus, one can note an increase in the number of 10-year-period cycles from two to three. The inter-epidemic period is therefore of about three years. The size of the epidemics also reduces as time goes by.
Figure 5
Let us now hypothesize a second situation with a population birth rate higher than the one above. Let us posit that the birth rate = death rate = 0.0000548 (60 % higher). All other parameters remain constant. This situation is in line with that observed in underdeveloped countries (high birth rate, high death rate; R=15). Now, the number of epidemics in each 10-year period (Chart 3) has increased over that with a lower birth and death rates. Not only is the inter-epidemic period now about two years at the beginning of the observation period, but tapering off (damped oscillations) occurs earlier along with a more constant number of cases. Therefore, the larger the recruitment of susceptibles by births, the larger the number of cycles per time period and the less the variance. This is because the larger number of susceptibles available converts more regularly to infecteds and infectious, thereby providing less opportunity for more isolated but larger (and unexpected) epidemics. By the way, a large vaccination rate depleting the number of susceptibles may be thought of as the opposite of a large recruitment rate by births. It is well known that vaccination translates to an increase in the time lapse between epidemics with a concomitant reduction in epidemic size. Thus, vaccination has the same effect as lowering the birth rate.
Figure 6
Now, let us examine what would occur with lowering the birth rate by 60% (in comparison with Chart 2), but with the death rate maintained to Chart 2 value. Chart 4 shows that the number of epidemics then decreases to one per 10-year period. Simply decreasing the birth rate at 60% of its original value thus increases the variance of epidemics size.
Vaccination
Now, what is the impact of vaccination on the dynamic. For example, what would be the outcome if 60% of individuals were vaccinated. Vaccination reduces the disease transmission rate because there is fewer susceptibles (effect similar to lowering the birth rate). If R=15 in the unvaccinated case, this means that one infectious will infect 15 individuals. But if 60% of the population is vaccinated, 9 individuals (60% of 15) will escape disease (by going directly in immunes compartment) and only 6 would be exposed to effective contacts. Let us therefore examine the dynamics of measles by considering R=6 with a birth and mortality rate of 0.0000351 (stationary population). To allow for the relevant changes, one must modify the susceptibles and immunes differential equations (in red in the text) that now include one more component such that:
dS/dt = λ*S(t) + b*N*
where b is the birth rate, N the population size, and p the vaccination rate. Thus, p represents the proportion of new births that go directly to the recovereds compartment with no processing through the exposed and infectious states, and (1-p) is the proportion of new susceptibles by births (60% of births go directly to the immunes compartment because of vaccination and 40% to the susceptibles). Further, to take care of the fact that 60% of individuals are already immunized in the population, we make the simulation with a population of 400 000 susceptibles and 600 000 immunes. The force (???of infection is then????*I(t) = (R/N*D)* I(t) = (6/400 000*7)*I(t) = 0.00000214*I(t). All other parameters are constant and identical to those of Chart 2 (birth rate= death rate= 0.000351). The number of new cases that occur with a vaccination effectiveness of 60% are presented in Chart 5. Clearly, the dynamic pattern is identical to that of Chart 4. This means that reducing the birth rate by 60% has effects on the dynamic similar to vaccinating with 60% effectiveness. In particular, the epidemics are more dispersed and more unpredictable. The resonance effect also continues to taper off with time.
Chaos
The tapering-off effect is at variance with real observations (see the NY City case in Chart 1). The tapering-off effect is due to the limitations of the SEIR model as we use it here. Most modeling studies of measles are predicated on seasonal variations of the annual contact rate, high during school months and low in-between. Let us attempt to take care of this refinement. One may tackle seasonal variation in various ways. More often than not, authors have been using a sine function to introduce variation in the contact rate, but a sine function does not truly capture reality. More recently, term-time variations have been used (3).
Term-time variation in North America involves an annual transmission rate low during the summer holidays (say, from June 24th to August 28th) and high for the remaining days of the year. This variation in the annual contact rate is more in line with the school year risks than a sinusoidal variation. Things therefore go like this for R:
From day 175 to day 240 (June 24th to August 28th), i.e., 65 days at low contact rate:
(R/ infectiousness period) - y)/population
and from day 1 to day 175 and day 240 to day 365, i.e., 300 days at high contact rate:
(R /infectiousness period) + x)/population
where x and y are unknowns to be estimated from the data. Many different combinations of x and y are possible, but the solution to the equations is unique if account is taken of the average beta that must be equal to {(R/ infectiousness period)/population}. The solution is therefore given by:
{((R/ infectiousness period)+x)*300 +((R/ infectiousness period)-y)*65 }/{365*population} ={(R/ infectiousness period)/population},
which simplifies to: 300 x = 65 y
If we arbitrarily posit that x=0.2, then y=0.923; this determines beta to be : 0.000002343 {((15/7)+0.2)/1000000} for school days (300 days), and 0.00000122 {((15/7)-0.923)/1000000} for the remaining 65 days (summer holidays). This solution fits an average beta of {(15/7)/1000000}. Chart 6 displays the step function of the annual variation in the contact rate to be taken on by the differential equations (Chart 6 scale is arbitrary).
The introduction of this simple parameter change in the transmission rate allows for a definite improvement in the epidemic pattern. Once the transient values have died away (the first half of Chart 7 diagram), the dynamic mimics reality much more closely.
Figure 10
Tapering-off has withdrawn, and a regular pattern of epidemic occurrence has now settled down. One can also note chance-like variations from one epidemic to the next similar to those observed in real populations. This is chaos, i.e., order embedded in disorder. Order (epidemic regularity) is ascribed to the deterministic differential equations, and disorder (chance-like variations) to the effect of the nonlinearities that allow for sensitivity to initial conditions. This means that any small variation in the epidemic outset conditions (demographic data) will determine a pattern-specific epidemic, i.e., an epidemic, the details of which are reducible to no other. Therefore, using nonlinear deterministic equations with few refinements (birth and death rates, seasonal variation in the contact rate) has allowed us to recreate the pattern of true epidemics.
Discussion
The SEIR model including birth and death rates aimed at providing a more realistic picture of the dynamic of epidemics. We have shown that the nonlinear equations of the SEIR model while providing for the dynamic recruitment of susceptibles by birth could rescue the oscillations observed in real populations. We have also shown how vaccination could alter the population pattern of infections in the same way that the birth rate does it. The SEIR model nevertheless departed from reality as oscillations were tapering off with time. We concluded that this pattern was at variance with real data as unvaccinated populations exhibit continuous biennial or triennial cycles. This suggested that other factors ought to be taken into account by the model to more closely mimic true epidemics. Seasonality in the contact/transmission rate could explain this inconsistency.
The picture of population infections presented in this article nevertheless remains static. True epidemics have patterns of infection that change with time. The SEIR model including birth and death rates, vaccination, and seasonal variation in the transmission rate however possesses all the basic characteristics to comply with any observed dataset (6). This means that the key parameters of the populations dynamics are embedded in the model. For example, it would suffice to slowly decrease the birth rate in time to reduce the pool of susceptibles, entice lower epidemic peaks, introduce more dispersion and, finally, more (apparently-) random fluctuations. An appropriate combination of the key parameters would thus explain most epidemic patterns. Be this as it may, it would be impossible to mimic more than qualitatively the details of epidemics; this is because of sensitivity to initial conditions that can alter patterns of incidence in tricky ways that would necessarily come up as apparent randomness.
Furthermore, the SEIR model remains crude as some other exogenous factors have not been considered in this paper; this is the case for immigration and the mixing pattern by age (age-dependant transmission). Also, we have considered the introduction of only one infected case in the population while in reality several could have been introduced over time. We have also assumed that all individuals are born susceptible. This is a reasonable assumption as long as small infants comprise a relatively small proportion of the total population. It would nevertheless be possible to include another compartment reflecting the maternally-derived immunized individuals who would become susceptible to infection at a given age and at a specified rate. Such refinements have not been embedded in the models because it is surmised they do not contribute much to the overall picture of the transmission dynamic. At any rate, they would have been responsible for additional variation.
This series of papers has shown that the dynamics of population infection hinge on both exogenous and endogenous factors. In particular, small variations in the demographic of a population can add to the random appearance of epidemics. Lately, Turchin and Taylor (7) have developed a method (response-surface methodology) to tell apart the role of endogenous and exogenous factors. Knowing what is up to the environment is crucial because the environment is amenable to intervention. In principle, it would suffice to understand the system that processes the environmental inputs to pinpoint how to intervene and to what extent intervention would yield the expected outcome (8, 9). Efforts to understand how complex systems learn (keep endogenous memory of the past) and evolve (respond to exogenous challenging factors) are needed to better master the key parameters of public-health intervention (10). The population dynamic approach used here can help control many types of population disease patterns.
Acknowledgments
We would like to thank the Social Sciences and Humanities Research Council of Canada for financial support to HT.
Correspondence to
University of Montreal
Social and preventive Medicine
2375 Chemin de la côte Ste-Catherine
CP 6128, Succ A
Montreal, Qc
H3C 3J7
Canada