Deterministic SIR (Susceptible–Infected–Removed) models applied to varicella outbreaks

Summary

Deterministic SIR models were applied to simulate Susceptible–Infected–Removed and to estimate the threshold condition for varicella outbreaks in children, reported in Medellín, Colombia. The expected numbers of susceptible, infected and removed individuals were compared with observed cases from notification of varicella outbreaks to the local Board of Health and from survey data. The threshold condition was estimated by the basic reproductive ratio and by the relative removal rate, through which measures for preventing and curtailing the outbreaks were identified. The model demonstrated a reasonable fit to the observations, except in two of the six outbreaks which probably reflected under-registration of cases. In order to have prevented these outbreaks, between 4·4% and 52·9% of the susceptible population should have been vaccinated assuming an 85% vaccine effectiveness. Similarly, isolation of affected children should have been increased to between 4·3% and 44·8% per week.

INTRODUCTION

The dynamics of varicella have been studied through mathematical models especially in developed countries, focusing on, e.g.: the underlying process during the epidemic or endemic period for varicella, the morbidity effects of immunization against varicella, the influence of zoster in the epidemiology of varicella, and the seasonal variation in contact rates [1–9]. Nevertheless, few studies have been performed in developing countries like Colombia, whose tropical climate, demographic pattern and socio-economic conditions can influence the dynamical transmission of diseases.

In Colombia, a high level of vaccination against varicella is not provided through The Expanded Programme on Immunization. This disease is one of the most frequently reported (70% of all cases reported to national surveillance system) although notification of cases is not mandatory and outbreaks are only studied through epidemiological field research [10, 11].

The local Board of Health in Medellín, the second largest city in Colombia, does not have sufficient information through surveillance systems to prevent and control varicella outbreaks at school settings where the occurrence of varicella outbreaks in children is a concern of parents and teachers. Varicella is generally a self-limiting illness, but given its high contagiousness and rapid spread, could affect pregnant teachers, immunocompromised children and cause increased absenteeism of parents from work [12, 13].

We simulated the dynamical transmission of varicella outbreaks and illustrated how to identify preventive and curtailing measures in child-care settings, using a mathematical model validated against varicella outbreak data, reported to the local Board of Health by elementary schools and day-care centres for the first term in 2001, and by a survey of parents for each institution that reported a varicella outbreak.

We described the dynamical transmission of varicella outbreaks, using a standard deterministic SIR model whose populations are structured by sub-populations of susceptible or non-immune (S), infected (I) and removed (R) individuals [14, 15] and the forces of infection and removal among them.

Infection occurs through an effective contact between a susceptible population and an infectious individual who transmits the infection according to a specific infection rate. Inversely, removal of the infected population assists in impeding further outbreak dissemination at a specific removal rate through either: isolation or quarantine, immunization, or death [15]. Infected individuals are those who acquire varicella, which for us corresponds to varicella cases during outbreaks, according to a clinical or epidemiological confirmation.

We obtained the threshold to determine the conditions under which such outbreaks could have been started. Using this threshold we identified the number of infected individuals that should have been isolated at the beginning of the outbreak in order to curtail it and the vaccination level that should have been achieved before the outbreak started in order to prevent it altogether [14, 15].

We applied the threshold theorem in the standard manner through the estimation of the basic reproductive ratio, R ₀, which represents the number of secondary infections caused by the introduction of a single infectious case into a completely susceptible population [15–17]. Thus, when R ₀>1, an outbreak occurs [15–17].

We estimated the relative removal rate as an additional criterion of the threshold as described by Bailey [18]. This rate, which is the ratio between the forces of removal and the forces of infection, was obtained by two methods. In method 1, the information required was the number of susceptible individuals at the beginning of the outbreak and the number of cases accumulated when the outbreak ended, which denominated the outbreak size. Method 2, required more information for each unit of time regarding the number of removed individuals and the number of susceptible persons at the onset of outbreaks.

We estimated the relative removal rate as a way of gaining insight about the relationship between all forces and states that influence the epidemic process, including the number of susceptible individuals at the beginning the outbreak at t=0, the outbreak size or number of removed cases when the outbreak ended at t=∞, the number of susceptible individuals, infected and removed persons at each time, and the infection and removal rates.

METHODS

Data compilation

Data were obtained from follow-up investigation of varicella outbreaks in children aged <12 years that were reported to the local Board of Health of Medellín by schools and day-care centres during the first term of 2001.

A survey was designed and conducted on each institution that reported a varicella outbreak, including information regarding age, sex and socio-economic status, and illness history such as the occurrence of previous natural varicella, the date in which the rash first appeared on the skin, and the date that isolation and vaccination against varicella took place. Further monitoring of new cases of varicella was conducted for 1 month after notification of the outbreaks. The information in the survey was provided by the children's parents or caregivers for all cases and school contacts (classmates) of confirmed varicella cases.

The information on socio-economic status was provided by parents according to the official classification. There are six ordinal levels in the socio-economic status used (low: levels 1 and 2; medium: levels 3 and 4; high: levels 5 and 6).

In order to reduce recall biases, the survey included a calendar that responders used to associate important dates such as holidays and periods with key dates related to varicella (e.g. the first onset of the rash or when the child was first isolated). Information was collected within a week after notification of each outbreak. Incomplete surveys were identified and completed with class lists, and, when necessary, by telephone. The information was recorded using Epi-Info v. 6.04d (CDC, Atlanta, GA, USA and WHO, Geneva, Switzerland), and mathematical analyses were performed using Mathcad 2001^® (MathSoft Inc., Cambridge, MA, USA).

Simulation of susceptible and infected individuals and removed cases

A standard deterministic SIR model (Susceptible–Infected–Removed) was used to describe the dynamics of varicella outbreaks [18]. The estimation and validation of parameters were performed based on real data of varicella outbreaks reported in school and day-care centres as previously described.

The dynamics of these outbreaks are usually represented through SIR models, which are constructed with a system of differential equations that reflect the transition rate from the susceptible state to the infected state when they have effective contact, in agreement with the infection rate parameter (β) and the product of the number of susceptible individuals and the number of infected individuals. Similarly, removed state occurs through isolation and recovery of infected persons or through immunization of susceptible individuals according to the removal parameter (γ) or removal rate. The assumptions of the model are as follows [12, 13, 18–20]:

• • Varicella generates long-term immunity.

• • All individuals are equally susceptible.

• • According to the law of mass action, there is a homogeneous mix of susceptible individuals and a single infectious one, which occurs through effective contact among them.

• • The population is closed, so the influence of new susceptible individuals in the range of transmission time of the illness is negligible.

• • The volume of the population is sufficiently large to guarantee a deterministic description.

The system of equations for simulation of the expected numbers of susceptible and infected individuals and removed cases, according to the model, can be reduced to equation (1) for removed individuals:

(1)

where N is the effective population size, S(0) is the number of susceptible individuals at the onset of the outbreak (t=0), R(t) is the number of removed individuals at each time, γ is the removal rate, and ρ is the relative removal rate [18]. This equation was transformed for discrete time with the numeric simulation of the expected SIR curves.

The observed curves of susceptible, infected and removed individuals vs. time were obtained using data collected at each outbreak. The curve of removed cases at each time was obtained from the date of isolation of the children infected with varicella. The observed curve of infected persons by time was calculated by subtracting the accumulated number of removed cases from the amount of accumulated infected individuals up to the given moment at each time. Finally, the observed curve of susceptible individuals by time was calculated by subtracting the sum of instantaneous values of infected persons and the previously obtained removed cases from the effective size of the studied population.

The population's effective size was determined as a result of the difference between the size of the school population and the number of subjects in the immune group that was identified by the previous history of natural varicella and the history of varicella vaccination.

A χ² test for goodness of fit was used to compare the two SIR curves: the observed curve and the one expected using the model. The null hypothesis was that the expected SIR curves fit with observed SIR curves with a 5% level of significance.

Estimating parameters

The outbreak size, R(∞), corresponds to the accumulated number of removed individuals without taking into account the index case at time t(∞) that corresponds to the extinction time of that particular outbreak. R(∞) is estimated using the following equation [18], derived from equation (1):

(2)

The relative removal rate (ρ) was obtained using two methods. Method 1 was based on data obtained in each outbreak regarding the number of removed cases at the end of the outbreak R(∞), the effective population size N, and the number of susceptible individuals at the onset of the outbreak, S(0). With this global information we estimated ρ, derived from equation (2):

(3)

Method 2 required gathering data regarding the outbreak process at each time rather than compiling information at the end of the outbreak, R(∞). It was necessary to discover the number of removed cases at each time R(t), and the number of susceptible individuals at the beginning of the outbreak, S(0).

(4)

The least-squares method was used to estimate ρ, deduced from the equations for susceptible individuals and removed cases of the system of differential equations. The objective was to minimize the difference between the number of expected and observed susceptible individuals, according to a quadratic error. That minimal value of error was obtained using Mathcad Premium v. 10^® (MathSoft Inc.).

The γ parameter was derived from the removal rate definition, where γ=ρβ. The β parameter was derived with the expression β=E(t)/I(t)⋅S(t), where E(t) is the incidence of disease and I(t) is the accumulation of cases or prevalence [18]. The basic reproductive ratio [15], is defined as R ₀=S(0)/ρ.

Applying the threshold condition

An outbreak starts when S(0) is greater than the relative removal rate, or when R ₀>1. Once this threshold for the onset of an outbreak has been identified, two measures are established to prevent it. First, the number of children (X) necessary to vaccinate prior to the outbreak is considered using the following inequality: X>(S(0)−ρ)E, which is applicable when ρ is much greater than R(∞), and E is equal to 85% effectiveness of the vaccine against varicella [20]. Thereafter, the number of children that should be isolated per week is deduced starting from the observed parameter γ, estimating the children's removal rate that should be obtained so that ρ>S(0), which implies that R ₀⩽1.

RESULTS

Description of outbreaks

The data used comes from follow-up investigations of six varicella outbreaks which were reported by schools and child day-care centers from Medellín, Colombia during the first term of 2001. Outbreaks I–V occurred in March 2001 and outbreak VI in May 2001.

Outbreaks were regarded as an unusual increase in the number of cases based on the previous experience of teachers and after confirmation by investigators and public health practitioners using clinical and epidemiological criterion, and according to the number of cases during the same period of time in previous years.

From information of reported outbreaks, 4058 surveys were collected that corresponded to 84·3% of children registered in the schools and day-care centres targeted. In general, the data were collected by surveying more than 80% of all children registered in each institution, except for outbreak V where 62·4% of the total population participated.

The median age of the children was 8 years (average 7·6, s.d.=2·0), the sex ratio was 2:1 with 62·8% of the participants being male and the remaining 30·7%, female.

Outbreaks I and II occurred in child day-care centres with children of low socio-economic status (Table 1). Such settings had a total area between 40 and 50 m² and they had only one recreational space. In these outbreaks the lowest proportion of children vaccinated, a great proportion of susceptibles and a greater outbreak size were observed (Table 2).

Table 1. Population at school and day-care centres (DCC) by socio-economic status. Varicella outbreaks in Medellín, 2001

* Proportion in parentheses.

† Figures within parentheses are percent of studied population in accordance with the total number of children registered.

Table 2. Population effective size, initial susceptible and removed cases of varicella outbreaks in Medellín, 2001

In outbreak I there were 41 reported varicella cases during 5 weeks (Fig.), and a varicella attack rate of 77·4% (Table 2). Varicella cases were children aged between 2 and 5 years (median 3 years, average 3·3, s.d.=0·96), with a sex ratio of 1·7:1 for boys/girls (24/14).

Fig. Observed (–––) and expected (- - -) number of susceptible individuals, infected individuals and removed cases in varicella outbreaks in Medellín, 2001 (outbreaks I–IV).

Outbreak II was studied with a stochastic model because the effective size of the population was lower than that of the other outbreaks [21].

Outbreaks III and IV were reported by schools of intermediate population size, with children from low- to middle socio-economic status (Table 1). In outbreak III 20 varicella cases were reported during 5 weeks (Fig.), with a sex ratio of 2·3:1 (14 boys/6 girls) (median age 7 years, average 6·7, s.d.=1·8). The school had an average of 36 students per classroom, which had a median area of 48·6 m² by group and only one recreational space. In this outbreak, the lowest proportion of susceptibles appeared at the beginning of outbreak.

In outbreak IV, 42 varicella cases were reported during 6 weeks (Fig.), with almost four cases of boys for each girl (33/9) (median age 8 years, average 7·0, s.d.=2·5). There was one recreational space and an average of 35 children in each classroom, and a median area of 25 m², i.e. 1·8 children/m², making this the most densely populated school.

Outbreaks V and VI occurred at schools with a great number of children who were of middle to high socio-economic status (see Table 1). In these schools the greatest proportion of vaccinated children, the smallest size of outbreaks and the lowest attack rates were observed.

In outbreak V, 24 varicella cases (21 boys, 3 girls) were observed among preschool and primary age children (median age 7 years, average 7·0, s.d.=2·0). The school had several recreational spaces, but was used by children of different educational levels. In outbreak VI 73 varicella cases were reported, with a sex ratio of 1·7:1 for boys/girls (43/6) (median age 8 years, average 7·0, s.d.=2·5). This school had schedules and spaces for recreation each educational level.

Estimating parameters and the threshold condition

Table 3 shows relative removal rate (ρ) estimated by two methods: R ₀ values, and the parameters β and γ. According to the data, estimates of ρ showed satisfactory consistency between the two methods, but in all outbreaks, method 1 generated a ρ value lower than the ρ values obtained by method 2. It was observed that the infection rate was greater than the removal rate in all outbreaks, except for outbreak III where the rates were equal.

Table 3. Estimation of parameters (R₀, ρ, β, γ) and control measures for varicella outbreaks in Medellín, 2001

The threshold condition is achieved if S(0) is greater than ρ, or R ₀>1. β is the infection rate and γ is the removal rate.

* The minimum level of vaccination is indicated in parentheses, with regard to S(0).

† Percentage increase in parentheses.

The basic reproduction rate was R ₀>1 for outbreaks I–IV while in outbreaks V and VI it was R ₀⩽1, which means that the threshold condition for the onset of outbreaks was not met. For this reason, the values for parameters β and γ for outbreaks V and VI are not presented.

Table 3 shows that outbreak IV had the highest parameter of infection and removal with nearly six cases per week while five cases were recovered during the same time period. In outbreak III, the infection rate was equal to the removal rate, which is demonstrated through R ₀ being just slightly greater than 1.

Simulating outbreaks

The Figure shows the observed SIR curves for susceptible individuals S(t), infected individuals I(t), and removed cases R(t) compared with expected SIR curves for every outbreak. In general, it was notable that the curve of susceptible individuals decreased as the outbreak progressed over time; removed cases increased from the initial value of zero to the saturation value when the outbreak ended, and the curve of infected individuals initially increased to a maximum value and then decreased to a final value of zero, in line with the typical bell-shaped curve.

The χ² test of goodness of fit for four outbreaks are given at every plot, with α=0·05. An adequate fit between the theoretical model and empirical data was observed, especially for the dynamic of susceptible and removed individuals in all outbreaks, except for outbreak IV.

Preventing outbreaks

Based on the results, these varicella outbreaks could had been prevented if (a) a certain proportion of children had been vaccinated before the outbreaks took place, reducing the susceptible proportion at the beginning of the outbreak S(0), or (b) if the proportion of isolation among the subjects had been greater, increasing the removal rate γ, or decreasing the infection rate β, by reducing the level of crowding in classroom and recreational spaces, as shown in Table 3.

The greatest efforts to prevent varicella outbreaks should have been carried out during outbreak I, where more than half of susceptible children should have been vaccinated before the onset of the outbreak. Additionally, when that outbreak started, the isolation of children should have been increased by 44% with regard to the actual removal rate observed during the outbreak. In outbreaks II and IV, both the vaccination and isolation rates in children should have been increased by 20% in order to prevent the spread of the virus. At the time that outbreak III occurred, it would have been required to vaccinate and/or isolate <5% of the children.

DISCUSSION

We analysed the dynamics of varicella outbreaks notified by schools and child day-care centres in Medellín, Colombia during the first term of 2001. The main purpose was to identify and propose prevention and control measures of outbreaks through the prior vaccination of susceptible individuals and the isolation of infected persons through estimates of the threshold condition. This deterministic threshold for outbreak occurrence is useful because it indicates how the accumulation of susceptible persons at a higher level than the relative removal rate, or R ₀>1, can produce outbreaks. The authors are not aware of any studies where this parameter was estimated for varicella outbreaks under similar conditions.

Research shows that varicella has an infectious period of about 5 days (or 0·714 weeks) [22], then, the removal rate is 1/0·714=1·4 per week per infective, which is the inverse of the infectious period. In the current study estimates of the removal rate (γ) were between 1·62 and 5 per week per infective which represents an observed removal rate greater than that reported in the literature. This may be due to the force of removal instituted by doctors, teachers or family members that had an effect on the isolation of new cases before they completed the infectious period.

In relation to the parameter β, Deguen et al. [8] estimated a parameter β of varicella for large urban populations, with a minimum value of 1·66 per week per infective and a maximum value of 8·9 per week per infective. The estimations of β that were obtained in our study were between 2·5 and 6 per week per infective. Caution must be exercised when comparing the two studies since the current findings are based on outbreaks in schools and day-care centres whose population's effective size is about 1000 children, considered in principle as closed populations, whereas the other study was conducted on a larger, more heterogeneous population.

We estimated the parameters of removal and infection for describing the step-by-step advancement of the dynamical transmission of varicella and we defined the threshold theorem by the relative removal rate corroborated by the basic reproduction rate. With this approach, we gained understanding about the influence these forces possess and what kind of information should be collected through the surveillance of varicella outbreaks by the local Board of Health.

In this study, two methods were implemented to estimate the relative removal rate. Method 1, is of easy application because the information required can be obtained through regular field research of the outbreak, while method 2 requires more detailed information about removed individuals at each time and susceptible individuals before the onset of the outbreak.

In general, method 1 resulted in smaller estimated ρ values compared with the estimates of method 2. Based on this finding, method 1 can be routinely applied to identify the threshold condition for outbreak establishment. This can be done by identifying the number of susceptible children by vaccination or previous varicella before the outbreak occurs (i.e. before the school year begins) and by measuring the number of accumulated cases at the end of the outbreak through notification of all cases of varicella to the local Board of Health.

The results of this approach show that simulation of SIR models have a reasonable fit to observations except in the cases of outbreaks V and VI given that the number of susceptible individuals at the onset of the outbreak was inferior to the relative removal rate and the basic reproductive ratio was R ₀⩽1.

In outbreak VI, the fact that a holiday period began 3 weeks after notification of the index case might have played a part in influencing the interruption of disease transmission. Other probable explanations could be the co-existence of a transmission process inside and outside of the educational institution or the under-registration of varicella cases.

Under-registration of cases can be corroborated in outbreak VI, if equation (3) is solved to deduce R(∞), knowing the number of susceptible individuals in the early stage of the outbreak, the population's effective size and the relative removal rate. This way, between 84 and 100 cases of varicella should have been observed in this outbreak, with an under-registration of 13–27% of cases expected.

From the Figure, it is evident that in outbreak IV, the curve of expected cases reached a maximum value before the curve of the observed cases. The initial recovery from t=0 in the fourth week appears to correspond to the removed cases that do not generate dissemination of varicella inside the school. This phenomenon was possibly due to a large relative removal rate of cases per week per infective or that such a rate was not constant throughout time, a variation that could be incorporated into the model.

The simulations and the χ² tests for goodness of fit show that the SIR model is more adequate to analyse the dynamic of transmission of varicella epidemics for outbreaks I and II where there is a substantial proportion of susceptible individuals at the beginning of the outbreak and closer contact between susceptible and infected individuals.

The design of a model that considers a population's heterogeneity or the non-homogeneous mixture of susceptible and infected individuals, to increase understanding of how socio-economic conditions influence this process and to compare the dynamics of varicella when outbreaks occur in an older population, might be required [24–27].