A COMPUTER SIMULATION MODEL RUBELLA EPIDEMIC
OF A
JOAN S. HORWITZ Program in Biostatistics, School of Public Health and Tropical Medicine, Tulane University. 1430 Tulane Avenue. New Orleans. Louisiana 70112. USA and DOUGLAS C. MONTGOMERI School of Industrial and Systems Engineering. Georgia Institute of Technology. Atlanta. Georgia 30337. U.S.A.
Abstract ~Thispaper presents a deterministic computer simulation model of the rubella epidemic process. The model consists of four difference equations which relate the number of infectant individuals in a closed population to the number of susceptible individuals. the number of immune individuals. the contact rate, the birth and death rates, and the proportion of vaccinated susceptibles. After computer simulation, the model is validated using data from the East North Central section of the United States. Satisfactory agreement between predicted and reported cases is obtained. It is demonstrated that under reporting of actual cases is a significant factor in modeling the epidemic process. The efficiency of various vaccination programs in eradicating rubella is investigated. Simulation
Mathematical model
Rubella
Preventive medicine
Epidemics
INTRODUCTION
The modeling of a disease process provides better insight into the epidemiology of the disease and is a valuable tool in evaluating effective control programs. This paper deals specifically with a computer simulation model of rubella. Better known as German measles, rubella is of both social and economic importance. In the 196465 epidemic alone, there were approx 12,5OO,OOOrubellacases of which 20,000 produced children with congenital rubella syndrome. a disease found in the offspring of pregnant rubella victims which can be characterized by hearing loss, cardiovascular defects, eye defects, developmental and neurological defects and often cases of spontaneous abortion. It is estimated that this one epidemic cost the American public $1,461,274,000 in physician and hospital services. loss of potential earning and special provision for those with congenital rubella syndrome [l]. The frequency of rubella is characterized by a period of 69 yr between epidemics with a strong seasonal pattern within periods. The seasonal peak appears in the spring with the low point in the late summer. This paper presents a deterministic model of the rubella epidemic process. The model has two objectives. First. the model is constructed so as to both represent rubella epidemiology and provide a reasonable *‘fit”to real data conforming to marked patterns in frequency. This deterministic model includes three classes of rubella: susceptibles, infectants and immunes. Second. the mode1 projects the course of rubella under alternative vaccination programs. Although the data used to develop this model was taken from the East North Central section of the United States, the mode1 has general applicability. It would be particularly IX9
useful in developing nations without especially when vaccination programs DEVELOPMENT
the funds to sponsor mass vaccination programs. on a smaller scale could eradicate the disease. OF’
A RUBELLA
MODEL
The rubella process can be described in terms of the network in Fig. I. Individuals in the population arc classified according to their state of disease. The susccptibles are those not presently infected with rubella but arc capable of becoming infected v4:hen exposed to the disease. The infectants have the disease (with or without overt symptoms) and are capable of transmitting the disease to others. The immune individuals have at one time had the disease and arc not capable of being reinfected.
I
I
In addition. births and deaths taking place in the population arc constantly altering the total population size. In the case of rubella. all newborns become susceptible (an infant actually retains his mother’s immunity for approx h months: however. for the purpose of the model. this period of temporary immunity will not bc included). Therefore. all births will bc added to the susceptible population. On the other hand. deaths are from all groups, or the total population. In many diseases the death rate among infectants may be higher than among the well population. Although one may die from rubella (as with measles encephalitis). the probability is very small and for the purpose of this model will be considered to be zero. In other words, it is assumed that rubella does not increase the risk of death and the death rate will be constant among all groups. Some additional assumptions and parameters have been considered in the development of this model. Chemotherapy, the treatment of disease by chemical reagents. does not alter the course of rubella in terms of a more rapid recovery and is only for the purpose of treating symptoms. Since it does not modify the length of disease nor the infectivity of the disease. chemotherapy will not be a parameter in this model. Prophylaxis does exist in the form of a relatively new rubella vaccine. Although antibody levels (a measure indicating immunity) arc generally lower in vaccinated individuals as compared to those observed in response to natural rubella infection. this vaccine does protect against illness on natural esposurc. or adequate contact. According to the Public Health Service Advisory Committee on Immunization Practices. [2] such antibody levels in vaccinated individuals have declined slightly over a 5yr observation period (although immunity persists), but only continued observation can absolutely determine whether longterm protection can bc cupcctcd. The vaccine and natural infection will be assumed to provide longterm protection or immunity against rubella. Since the model is a shortrange prcdictivc device, the authors t’ecl this is a valid assumption.
A computer
simulation
model of a rubella
epidemic
191
As a last consideration, immigration and emigration will not be included because the exact nature of these would be impossible to define. Therefore, the model will be restricted to a closed population. The rubella epidemic process can now be described by the following difference equations : Cf +
1, +
P, +
1 =
w,s,/pt
1
I, +
=
1 =
p,
+

YC,
(1)
c,  yz,+ vs,
(3)
P,  1’Pt
(4)
where S, = the C, = the I, = the P, = the
/?, = ,u =
and
y= r =
number number number number
of of of of
individuals in the susceptible category at time t, individuals in the infectant (or case) category at time t, individuals in the immune category at time t, individuals in the total population (P = C + S + I) at time t, a contact rate between susceptibles and infectants leading to infection of the susceptible individual at time t. the birth rate, the death rate. the percentage of susceptibles being vaccinated.
The term /lJ,S,/P, represents the number of new infectants for time period t + 1. Assuming complete recovery each period. the total infectants at time t + 1 equals the new infectants @‘,C,SJP,) minus the deaths (yC,). The susceptibles at time t + 1 equals the susceptibles at time t (S,) minus the new infectants, plus the births (pP,), minus the deaths (yS,), minus the vaccinated susceptibles (vS,). Similarly, the immunes at time t + I equals the immunes at time t plus the recovered infectants (C’,),minus the deaths (~1,) plus the newly vaccinated susceptibles (vS,). The total population, thus, at time t + 1 equals the population at time t (P,), plus the births minus the deaths (yP,). A deterministic simulation was chosen so that large populations could be analyzed. The model is deterministic because given the initial numbers of susceptibles, infectants and immunes, along with the birth, death, and contactrates, the future state of the epidemic process can be uniquely determined. Solving equations (lH4) at time t = 1, the results can be used to resolve these equations at time t = 2, and similarly up to time t = 26. The computer program is written in FORTRAN IV and final runs were made on an IBM 7044. A program listing is given in Horwitz [3]. SOURCE OF DATA AND INITIAL CONDITIONS The basis for most of the rubella data used in this study is the Center for Disease Control (CDC), Atlanta, Georgia. In 1966, the Conference of State and Territorial Epidemiologists placed rubella on the list of notifiable diseases. i.e. each state is required to report the number of infected cases at the end of each week. Prior to 1966, many states voluntarily reported rubella cases, but this did not necessarily mean that the entire state was reporting.
Although data are now submitted in a weekly telegraphic report. the completeness in reporting as well as the type and accuracy of the information varies both between and within states. It is wellknown that data from certain sections of the country are more rcliable than others. Particulariy reliable is the East North Central section of the United States. On this basis the East North C’cntral section was chosen as the test population for study. This section includes the states of Ohio. Indiana. Illinois. Michigan and Wisconsin. Due to data availability. the simulation study begins in 1966. Based upon fitting the model to 1966 data. expected data for future years will he predicted. It ia important to note that most previous models of this type have been tested on hypothetical populations. In few cases has the predictive value hen wdunttxi on the basis of comparing the predictions to real data. The model is applied in interials of ?ueck periods. the approximate length of time ;I patient is infectious. Since rubella is characterized by a fairly constant period of infection. the ~1st’ of discrete time periods is reasonable. Another reason for Zweek intervals involves reporting. Since rubella is reported weekly. Zweek reporting smooths the data. The model is dependent upon certain initial conditions and the rates involved. These conditions arc the numbers of persons in each of the three categories (infectants. immunes and susccptibles). The sum of the three catcgorics is the total population. The total study population is equal to the total population of the East North Central section of the United States at the beginning of 1966. According to L!.S. census data the population at that time was 3x.4x0.000. The initial number of infectants is based on the number of reported casts during the first 2 weeks of 1966. For the East North Central section. this number is 416 infectants, IO3 in Week 1 and 3 I? in Week 2. The initial number of immunes is much more difficult to define. as this statistic is not record&. It can not bc counted on the basis of a person’s rcmembcring having a historq of the disease. Several stratified random serosurvcys [l] have been conducted in different arcas of the country to define agespecific immunity on the basis of persons with rubella hemagglutination inhibition (HI) antibody. Most of these studies have placed the per cent of seroimmuncs in the total population. in 1966. somewhere between 75 and 80 per cent. It is important to remember that since rubella reporting is “incomplete and diagnostic accuracy variable, and since a significant proportion of rubella infections are subclinical” [I]. serological data may be the only way to determine the information regarding susceptibility. On this basis, the initial number of immunes will be assumed to be 30.000,OOO (approx 78 per cent). The remaining initial population parameter is the number of susceptibles. Since S = P  C  I. the initial number of susceptibles is X.479.584 (approx 22 per cent). The initial rates to be estimated are ,D. ;’ and /j. The birth rate. ,L’and the death rate. 7. were based on 1970 census data [4J which give vearly birth and death rates for different sections of the country. Since there is little variation in these rates, let ,L!~= 19. I per thousand where /lr is the yearly birth rate, and let ;tr = 9.7 per thousand where ;I~ is the ycarlq death rate. However. since this simulation 1s based on 26 twoweek period. the birth and death rates need to be expressed as biweekly rates. In this case ,LL= 0.000728. where 11 is the biweekly rate and ;’ = 0.000375, where 7 is the biweekly death rate. In choosiag p. the contact rate. we observe that it is a measure of the susceptibility of the host. the infectivity of the rubella virus. and the social conditions in the population
A computer
simulation
model of a rubella
Table 1. Initial population Population
category
epidemic
I’)?
values Number
individuals
X.479.584 416 30,1)00.000 38,4xn,ooo
Susceptibles Cases (infectants) Immunes Total population
at risk. It is obvious that these are subject to change over time and may be influenced by such factors as children being in and out of school. Therefore. if the reported rubella cases are used to calculate this contact rate. it seems reasonable to assume that /l varies over time within each yearly cycle. From the epidemiology of rubella, it is known that there is a strong seasonal pattern in reported rubella cases per year. Based on this cyclic variation. /j must be timedependent and will be denoted by [j,. t = 1. . . ...26. The &‘s will be calculatcd from 1966 rubella data and then substituted into the model to predict the following years. Referring back to equation (l), we express the number of infectants at time period t + 1. sayC,, ,. in terms of the number of infectants at time t as C I + , = P,S,C,IP,  c,.
(5)
Solving equations (5) for 0, we obtain P, = VJM(C, Table 2. Contact
f
+ i/C,) + Y). rates. /j,‘s
P, 1 5.02I 6.07
6.769 5.517 4.9, I
4.313 4.131 J.X1;6 4.671 4.86 I
3.282 I.506 ‘I.156 7.7I7 4.902 3.140 4.6?0 3.732 4.66 I
?.hN 23 74 25 26
4.069 5.530 4,897 5.062 4.504 3.845
(6)
194
JOAU
S.
HOKWITL and D. c‘. MONTGOML~1
Substituting initial values for susceptibles. immunes, total population size. deathrates, and the reported number of infected cases per week for Iyr. all the needed information is known. Thus. the pt’s can be computed for I = 1,2....,26 from equation (6). The initial population values and the values of/j, arc shown in Tables 1 and 3, respectively. IMPLEMENTATION When employing the model using these conditions, it was found that the epidemic process did not continue very far into the first year of predictions (1967). In attempting to find a cause for the failure of the model, the possibility of under reporting was reconsidered. Under reporting can be attributed to the fact that due to the lack of side effects associated with rubella, a victim often does not see a physician, especially when it is the second or third infection in a family. In addition. even the medically attended patients are usually not reported to the local or state health department. According to Witte et al. [S]; It is estimated that only one case of clinical rubella in 20 is reported. Although these surveillance data cannot be considered quantitatively accurate, it is useful to depict the trends and patterns of rubella occurrence in the United States. If this is the case, it is reasonable that a model based on such under reporting should quickly diminish since such a small number of cases would not keep the epidemic going. To test this hypothesis. a correction factor was introduced increasing the reported cases by a multiplicative constant. The results of a simulation with an arbitrarily chosen multiplicative factor of 22 applied to the infectant portion of the population are displayed in Fig. 2. The cases have been accumulated also by 4week periods in Fig. 3 and by Xweek periods in Fig. 4 in order to depict what is occurring by months and by seasons. The total number of infectants are displayed as reported. Note that although the multiplicative constant has been incorporated into the model. the predicted data has been redivided by this constant before being recorded in Figs. 3 4. We have found that the model is relatively insensitive to the choice of correction factor. From inspection of Figs. 2~~3, we ohserve that the model appears to adequately characterize the rubella epidemic process.
A computer
simulation
model of a rubella
epidemic
The actual cases (those reported to CDC) in 1967 and 1968 were 8,311 and 11,232, respectively. The predicted cases, produced by the model for these 2 years were 8,132 and 13.923. It is interesting to note from these graphs that the predicted epidemic peaks later than actual. a characteristic also found in Abbey’s [6] ReedFrost application to German measles. In both 1967 and 1968. the rise in predicted cases occurs at a slower pace than actual, but upon reaching a peak, the predicted cases decline just as rapidly (and at basically the same point in time) as the actual cases. Based on the results obtained from this simulation study, we conclude that under reporting is a significant factor which must be taken into account if a viable mathematical model is to be developed for rubella. Under reporting most likely plays an important role in other
of dealing with this phenomenon must be devcloped so that future studies may bc more accurate. The largest discrepancy in 196X predictions occurs at the latter part of the year, the last 4 to 6 months primarily. The predict4 cases for this time period arc consistantly greater than the actual reported cases. This is probably due to vaccination programs (although small in scope) that were initially implemented at approximately this time. Through 30 June 1970. 3,i25.925 doses of rubella vaccine wrrc administered through state and local health departments. [7] There arc two difficulties in incorporating this knowledge into the model. First. there is no accurate data regarding the susceptibility of the people receiving the vaccine. For csamplc. the procedure for administering in the schools is to sond a note home to the parcnts requesting permission to vaccinate their children. If the note is signed. the child reccivcs the vaccine. Whether or IIOL the child has had a histor! of rubella may or may not bc known to the parent. and could cvcn have no bearing on the parent’s granting permission. In other words, the population receiving prophylaxis may not all be susceptibles. In fact. CDC’ estimates that 30 per cent of the children receiving immunization arc already immune [7]. Second. these 3 million reflect only vaccinations given through the health dcpartmcnts. Accurate accounting of vaccine dispensed through private physicians and hospitals is not :kvailablc. epidemic processes as well. Effective means
EXTENT
AND
TIME
OF
VACCINATION
PROGRAMS
The advent of rubella vaccine in the late 1960’s has provided a solution to the rubella control problem. At first it seemed that this program would not be undertaken due to lack of funds. However, by late 1970 enough money was available to complete a mass immunization program consisting of 50 70 million vaccines. Although the United States was able to produce such funds. it is not always easy for other nations to do so. With limited resources. it is advantageous to have a method for testing the effects of alternative control programs. The simulation model presented in this paper provides the means for such a test. It will bc assumed that all those receiving the vaccine ;IIc susceptibles. In addition. it is understood that vaccinating on a random basis will not iaccinatc comparable numbers of immuncs in each age group. but rather a number proportional to the number of susccptibles in each group. As characterized h) equations (7) and (3). vaccinated susceptibles are removed from the suscoptiblc population and transferred to the immune population. Since all indications point towards longterm immunity. vaccinated individuals will remain immune.
Per cent susccpt~hlrs 01’\,
I II”,,
vaccinated
75” . ,,
in I%9 5w,,
75”,,
A computer
simulation
model of a rubella
197
epidemic
The control programs under consideration include no vaccination (results of which are displayed in Fig. 2) 10 per cent, 25 per cent, 50 per cent and 75 per cent vaccination of the susceptible population. Table 3 summarizes the results of these different vaccination programs, all taking place entirely within year 4, or 1969. As can be seen from Table 3, it is only necessary to vaccinate 25 per cent of the susceptible population in order to stop the spread of rubella infection in the entire population. It is important to remember. however, that in order to be sure of vaccinating 25 per cent of the susceptible population, it is usually necessary to vaccinate more individuals because vaccination programs are usually conducted without being sure of each individual’s status. Also, infection may reoccur if it is introduced from outside the population under study. A necessary assumption of this model was that the population be considered closed, i.e. no immigration or emigration. If infection dies out in a closed population. it cannot start again. It can also be seen in Table 3 that a twofold increase in the per cent vaccinated (in the case of increasing from 25 to 50 per cent) produces only a onethird decrease in the number of cases in the epidemic. It is necessary to triple the 25 per cent program to 75 per cent in order to halve the number of cases, from 24,693 to 12,592. Table 4. Timing
study with 25 per cent susceptibles
vaccinated
Year
1969
Year vaccination takes place 1970
1971
1969 1970 1971 1972 1973
24. Ill 516 0 0 0
42.015 35.380 I10 0 0
42,015 53,760 10.598 30 0
Total size
24,693
77,525
105,393
Timing also plays an important part in the effect vaccination programs have on the course of rubella. Considering a total vaccination program of 25 per cent, Table 4 shows that if the 25 per cent program is delayed just 1 yr. the number of total cases increases from 24,693 to 77,525 and increases to 105,393 if the same program is delayed 2 yr. SUMMARY The primary purpose of this research was to relate both the mathematical theory of epidemics and the epidemiological characteristics of rubella in the form of a mathematical model. The method of solution chosen was a discretetime computer simulation model in the form of a set of deterministic equations that would describe such an epidemic process in a large closed population. The population under study was the entire East North Central section of the United States beginning in 1966, an area considered very reliable in reporting rubella cases. The population was considered to be unstructured, there being no difference between individuals other than state of disease. The initial conditions were chosen from census data for population size, birth and deathrate for this area in 1966. The population was divided into three disjoint groups (susceptibles, infectants and immunes) on the basis of the results
19x
JOAN S. HORWITZ and D. C. MONTGOMFR\
of serological testing. The contact rates between these individuals were based on 1 yr reported cases. It is shown that the hypothesized model yields an adequate representation of the rubella epidemic process. The cases predicted by the model are reasonably close to the number actually reported. The largest discrepancy occurred in the latter part of 1968. and is likely due to the implementation of vaccination programs. It is also demonstrated that under reporting of actual cases is a significant factor which must be accounted for. Several alternative vaccination programs were compared beginning in 1969 and it was found that only 25 per cent of the susceptible population need be vaccinated in order to completely eliminate the rubella infection from the population. Vaccinating any larger fraction produced diminishing returns. The timing of the immunization program is also an important factor. REFERFNCES I. Rubella survtiliancc. N.C.D.C‘.. Atlanta. No. I (IYhY). 2. Rubella virus vaccine. Morhidirr ctrld Mortalrr~ Wcrkl~ Rtzp. 20, 304305 (197 I ). 3. JOAV S. HORWITZ. A mathematical model of an epidemic process. Masters thesis, Georgia Institute of Technology ( 1972). 4. United States Dcpartmcnt of Health. F.ducatlon and Welfare. Public Health Servw. vital statistics of the United States. (1967). Washington, Government printing office (1969). 5, J. J. WITTE. Epidemiology of rubella. .4rn. ./. Dis. C’hiltl. 118, 107 (1952). 6. H. AIXBES. Examination of ReedFrost theory, Hww~ Bid. 24, 201 2.53 (1957). 7. R. WALLACL. N.C.D.C.. Atlanta. Immunization branch. personal communication (1970).