INTRODUCTION
Since its discovery in 1976, 14 confirmed outbreaks of Ebola haemorrhagic fever (EHF) have been reported [1–10], seven of which have occurred since 2000 (see Table 1) [1–Reference Okware8]. Very little is known about the virus: natural reservoirs are poorly identified but may include fruit bats; vaccine and therapeutic strategies are under development [Reference Rouquet11–Reference Leroy16]. The typical natural history of the disease begins with an average incubation period of 1–2 weeks. Patients present most frequently with fever, asthenia, diarrhoea, abdominal pain, headache, arthralgia, myalgia, sore throat, dysphagia, and conjunctivitis [1, 2, Reference Okware8, Reference Bray17–Reference Ndambi19]. One week after the onset of symptoms a rash often appears followed by haemorrhagic complications, leading to death after an average of 10 days in 50–90% of infections. Survivors may experience severe asthenia, hearing loss, ocular signs and recovery usually occurs in 2 weeks to 2 months after the onset of symptoms. Most individuals acquire infection after direct contact with blood, bodily secretions and tissues of infected ill or dead humans and non-human primates [2, Reference Dowell20–Reference Peters and LeDuc22]. There is evidence that individuals (health-care workers, relatives) may become infected following contacts with patients' body fluids or direct contact with patients during a visit at the hospital or participation in traditional burial ceremonies [Reference Dowell20, Reference Roels23, Reference Hewlett and Amola24]. Ebola is unlikely be transmitted during the incubation period and transmissibility increases with duration of disease and direct contact with infected individuals during the late stages of illness [Reference Dowell20, Reference Borio21]. During the 1976 outbreak in Democratic Republic of Congo (DRC, formerly Zaire), 86 (26·7%) of the 318 cases were infected either from a contaminated syringe or through needle-stick injury [2]. However, for epidemics that occurred after 1990, infection from contaminated syringes or through needle-stick injuries has not been documented. Although there is evidence of asymptomatic carriers, the very low levels of virus detected in these individuals suggest they do not pose a significant source of transmission [Reference Leroy25, Reference Baxter26].
Table 1. Confirmed outbreaks of Ebola (excluding isolated cases)
CFR, Case-fatality ratio; DRC, Democratic Republic of Congo (formerly Zaire).
Bold text corresponds to the epidemics analysed in this paper.
Two studies have proposed estimates of the average number of secondary infections generated by one primary case of Ebola in an entirely susceptible population [Reference Chowell27, Reference Ferrari, Bjornstad and Dobson28]; this quantity is called the basic reproduction number and is denoted R 0. The first study proposed a compartmental model and fitted it to historical data to estimate the R 0. The second study proposed estimates of R 0 based on the chain binomial model. These two works did not study the contribution of the different settings for transmission (in the community, in the hospital, during burial ceremonies) in the estimation of R 0.
Here, we analyse previously published EHF data with a stochastic compartmental model which incorporates explicitly the settings of the transmission in the community, in the hospital and during burial ceremonies. Our goal is to better understand and to provide insight into where control interventions should be targeted in the future. We subdivided the infectious phases into three stages to account for transmission in the community, in the hospital (including isolation wards), and after death during traditional burial. We provide maximum-likelihood estimates of R 0 for EHF and the portion of each infectious phase represented in this value. We propose an analysis of the potential effects of control interventions on the dynamics of an Ebola epidemic.
METHODS
Data
We analysed data from two recent epidemics briefly described below (Table 2).
Table 2. The stochastic compartmental model
S, Number of susceptible individuals; E, number of exposed individuals; I, number of infectious cases in the community; H, number of hospitalized cases; F, number of cases who are dead but not yet buried; R, number of individuals removed from the chain of transmission; βI, transmission coefficient in the community; βH, transmission coefficient at the hospital; βF, transmission coefficient during funerals. θ 1 is computed in order that θ% of infectious cases are hospitalized. δ1, δ2 are computed in order that the overall case-fatality ratio is δ. The inverse of the mean duration of the incubation period is α. The mean duration from symptom onset to hospitalization is γh−1, γdh−1 is the mean duration from hospitalization to death, and γi−1 denotes the mean duration of the infectious period for survivors. The mean duration from hospitalization to end of infectiousness for survivors is γih−1 and γf−1 is the mean duration from death to burial. Values presented in days in Tables 3 and 5 were converted to weeks for computation. Transmission coefficients are expressed in weeks−1.
Kikwit, DRC, 1995 [Reference Khan3, Reference Muyembe-Tamfum29–Reference Kerstiens and Matthys31]
The epidemic took place in Kikwit (≈200 000 inhabitants) and its surroundings. A total of 315 cases were identified with an 80% hospitalization rate and an 81% case-fatality ratio (CFR). First interventions were implemented on 4 May 1995: cases were placed in an isolation ward at hospital, body burial was done by the International Committee of the Red Cross, gloves and personal protection were distributed in households and community education concerning risk of transmission was implemented. EHF was confirmed on 10 May 1995. We fitted the model to the dates of onset of first symptoms which were available (291 Ebola cases out of 315 cases). These data are illustrated in Figure 1.
Fig. 1. Observed data (□) and fitting curves (—) for (a) the 1995 DRC epidemic and (b) 2000 Uganda epidemic.
Gulu District, Uganda, 2000 [4, 5]
Most of the 425 presumptive cases (confirmed and clinical) occurred in the district of Gulu (≈470 000 inhabitants) and the CFR was ∼53%. EHF was confirmed on 15 October 2000 but suspicion was strong enough to set up an isolation ward on 10 October 2000. Follow-up of contacts, community education and cessation of traditional burial were implemented. We fitted the model to the dates of onset of symptoms of 418 cases with available data. These data are illustrated in Figure 1.
The model
We developed a stochastic compartmental model where individuals are classified as: (1) susceptible individuals (S) who can be infected by Ebola virus following a contact with infectious cases; (2) exposed individuals (E) who have been infected by Ebola virus but are not yet infectious or symptomatic; (3) symptomatic and infectious individuals in the community (I); (4) hospitalized Ebola cases (H) who are infectious; (5) dead Ebola cases (F) who may transmit the disease during funerals; and (6) individuals removed from the chain of transmission (R, cured or dead and buried). The model structure, disease phases and parameters definitions are described in Table 2. Simulations of the model were performed using Gillespie's first reaction method [Reference Gillepsie32]. A transition rate, depending only on the present state of the population, is allocated to each transition λi (see Table 2). At each iteration of the algorithm, a time τi is drawn from an exponential distribution with parameter λi for each transition. The next transition μ is the transition that has the minimum time to occurence (τμ). Counts in each compartment are updated accordingly.
Understanding the dynamics of the disease
We fitted the dynamic model to morbidity data from the two epidemics described above. Parameter estimates for the model were drawn from the literature where available and estimated otherwise. Table 3 provides epidemic specific parameters values used for simulations. For both epidemics, we made the following assumptions:
(a) The entire population was considered initially susceptible.
(b) Interventions were completely efficient after the date indicated in Table 3 and not efficient at all before this date.
(c) Before interventions, the population was exposed to cases within the community, hospitalized and dead cases, as hospitals were open to the general community prior to interventions.
(d) After interventions, no transmission occurred at hospital or during burial and transmission in the community decreased. Then, the transmission coefficients at hospital and during burial are set to 0 and the transmission coefficient in the community is decreased by a factor (1−z).
(e) After developing symptoms, the mean infectious period for cases who survived was 10 days and dead patients remained infectious for an average of 2 days after their death.
(f) The model was initialized with the number of index cases indicated in Table 3.
(g) All observed cases (except index cases) were assumed to be related to human-to-human transmission.
We estimated the transmission rate in the community before interventions (βI), the transmission rate at hospital (βH), the transmission rate during traditional funerals (βF) and the efficacy of interventions in the community (z) fitting the model to morbidity data from the two aforementioned epidemics. We estimated these four parameters using approximate maximum likelihood. To evaluate the approximate likelihood of one set of parameters, we assumed that the weekly incidences were Poisson distributed with parameters equal to the average weekly incidences over 700 runs of the model (we verified that average incidences were stable when we performed 700 runs of the model with a given set of parameters). To determine the maximum-likelihood estimates and their 95% confidence intervals, we computed the likelihood of sets of parameters generated by Latin Hypercube Sampling (LHS) and we assumed that twice the difference of log-likelihood values was χ 2 distributed with the degrees of freedom equal to the number of estimated parameters [Reference Riley34, Reference Iman and Conover35]. LHS ensures that input data for the parameter value identification simulations cover the sampling space. The expression for R 0 (see Appendix) was determined following the method of Diekmann & Heesterbeek [Reference Diekmann and Heesterbeek36, Reference van den Driessche and Watmough37]. We found that R 0 can be written as the sum of three terms: a first term relative to transmission in the community, a second term to transmission during hospitalization and a third term to transmission during traditional burial.
Table 3. Epidemiological features of two outbreaks
* This is approximately the average duration from death to burial.
† No data was available for Uganda outbreak; the value used for DRC in 1995 was applied. This value corresponds to the most infectious period.
‡ No data was available for Uganda outbreak; the rate reported in DRC in 1995 was applied.
Simulation of scenarios of control measures
To evaluate the impact of control interventions on the two epidemic profiles, we performed two multivariate uncertainty and sensitivity analysis. We studied the effect of varying the following parameters: (1) the time to intervention T, (2) the hospitalization rate of Ebola cases after intervention (t>T), (3) the efficacy of the isolation ward and barrier nursing within the isolation ward (decrease of βH for t>T), (4) the efficacy of interventions during body burial for t>T, (5) the mean duration between onset and hospitalization for t >T. These parameters are called ‘intervention parameters’.
For the first multivariate sensitivity analysis, the parameters (except the intervention parameters, the size of the population and the number of index cases) were set to the values observed or estimated for the 1995 DRC epidemic (see Tables 4 and 5). For the second multivariate sensitivity analysis, these parameters were fixed to the values observed or estimated for the 2000 Uganda epidemic (see Tables 4 and 5). For both multivariate sensitivity analyses, epidemics were simulated in a population of 100 000 inhabitants starting with one index case (see Table 5).
Table 4. Parameter estimates of the model
CI, Confidence interval.
Table 5. Values of parameters for the multivariate sensitivity analysis
Bold text corresponds to the intervention parameters studied in the multivariate sensitivity analysis.
For each multivariate sensitivity analysis, 500 different sets of intervention parameters were generated by LHS assuming that the five intervention parameters were uniformly distributed [Reference Blower38]. The time to intervention was distributed between 4 and 10 weeks, the hospitalization rate after interventions varied between 0% and 100%, the mean time between onset and hospitalization varied between 1 and 5 days. The efficacy of interventions at hospital (1−z H) varied between 50% and 100%. The efficacy of interventions after death (1−z F) varied between 75% and 100%. These ranges are also reported in Table 5. For each set of parameters generated with the LHS (a scenario), we simulated 700 epidemics and computed the mean size of these epidemics between 1 and 51 weeks after the onset of symptoms of the index case. We computed the partial rank correlation coefficients (PRCCs) between each varying parameter and the mean size of the epidemic at weeks 1–51 [Reference Blower38]. PRCCs quantify the linear relationship between the ranks of one input variable (each intervention parameter) and the output variable (the epidemic size), after the linear influence of the ranks of the other variables has been eliminated.
RESULTS
Understanding the dynamics of the disease
For each epidemic, the best-fit model-based epidemic curve is plotted with the observed data in Figure 1. Figure 2 represents distributions of the peak of the incidence and distributions of the final size of epidemics for 1000 simulated epidemics with the parameters set to their best estimates. We estimated the R 0 at 2·7 (95% CI 1·9–2·8) for the 1995 DRC outbreak and 2·7 (95% CI 2·5–4·1) for the 2000 Uganda outbreak (see Table 4). For the 1995 DRC epidemic, the community component of R 0 accounted for 0·5 (95% CI 0·4–1·9), the hospitalization component for 0·4 (95% CI 0·0–2·2) and the burial component for 1·8 (95% CI 0·0–2·3). The transmission coefficient in the community (excluding hospital and burial) was reduced to 88% (95% CI 22–100) of its initial value after introduction of control measures. For the 2000 Uganda epidemic, transmission in the community, hospitals and traditional burial respectively accounted for 2·6 (95% CI 0·3–2·8), 0·01 (95% CI 0·0–3·5) and 0·1 (95% CI 0·0–3·2) in the value of R 0. The transmission coefficient in the community (excluding hospital and burial) was reduced to 12% (95% CI 8–99) of its initial value after introduction of control measures. After implementation of interventions, the effective reproduction number (neglecting the depletion of susceptible individuals) is 0·4 (95% CI 0·3–0·6) for the DRC epidemic and 0·3 (95% CI 0·2–0·4) for the Uganda epidemic.
Fig. 2. Distribution of the peak of the weekly incidences (left panels) and the final size (right panels) of the epidemic obtained with 1000 runs. The histograms represent the distribution of the peak and the final size of the simulated epidemics when parameters are set to their maximum-likelihood estimates. Black crosses represent observed data.
Simulation of scenarios of control measures
Figure 3 represents the PRCC between the mean size of the epidemic at each week and the five parameters studied: the time to interventions, the hospitalization rate of Ebola cases after interventions, the efficacy of interventions at hospital, the efficacy of interventions after death and the mean duration between onset and hospitalization of Ebola case after interventions.
Fig. 3. Partial rank correlation coefficients (PRCCs) between the cumulative incidences and the five studied intervention parameters. Epidemics were simulated in a population of 100 000 inhabitants with one index case and with values of parameters (except the five intervention parameters) estimated with data from (a) the 1995 DRC epidemic and (b) data from the 2000 Uganda epidemic. These figures represent the PRCC between each varying parameter and the epidemic size x weeks after the onset of symptoms of the index case (x varying between 1 and 51).
Impact of interventions when parameters are set to the values observed or estimated from the 1995 DRC epidemic
When the parameters of the model (except the five intervention parameters) were set to the values relating to the 1995 DRC epidemic, for 50% (5%, 95% respectively) of the set of parameters generated with the LHS, the mean size of the epidemic at week 51 was under 190 cases (35, 960 respectively). The PRCCs show that the most important intervention parameter for the control of the epidemic is the time to intervention (0·99 at week 15, 0·93 at week 50). The mean size of the epidemic is also linked to the efficacy of interventions after death (−0·88 at week 51), the hospitalization rate after interventions (−0·71 at week 51) and the efficacy of interventions at hospital (−0·53 at week 51). The PRCCs for the mean time between the onset of symptoms and hospitalization are smaller than the PRCCs for other parameters but larger than 0·25 (0·33 at week 20 and 0·25 at week 51).
Impact of interventions when parameters are set to the values observed or estimated from the 2000 Uganda epidemic
When the parameters of the model (except the five intervention parameters) were set to the values relating to the 2000 Uganda epidemic, for 50% (5%, 95% respectively) of the set of parameters generated with the LHS, the mean size of the epidemic at week 51 was under 135 cases (35, 590 respectively).
The PRCCs show that the most important intervention parameter for the control of the epidemic is the time before interventions are instituted (over 0·99 after week 10). The mean size of the epidemic is also linked to the hospitalization rate after interventions (−0·88 at week 50) and the mean time between the onset of symptoms and hospitalization (0·41 at week 50). The PRCCs for efficacy of interventions at hospital are about −0·06 after week 17. The PRCC for the efficacy of interventions after death is about −0·17 at week 50.
In both scenarios, the time to intervention was identified as a key parameter for the control of the epidemic size. Furthermore, PRCCs for the hospitalization rate were negative and larger than 0·5 in absolute values. Although less important, the mean time between onset of symptoms and hospitalization has an impact on the predicted epidemic size. For the transmission pattern identified in DRC in 1995, the size of the epidemic is strongly related to the efficacy of interventions at hospital and after death. On the contrary, the efficacy of interventions at hospital and after death was not identified as a key parameter when transmission at hospital and during traditional burial was small. Thus, when the transmission at hospital decreases by 50% at least after interventions, the rapid hospitalization of cases after the onset of symptoms may reduce the size of the epidemic.
DISCUSSION
We presented the results of a dynamic model for the spread of EHF and fitted it to two historical epidemics. First, we estimated the R 0 at 2·7 (95% CI 1·9–2·8) for the 1995 DRC epidemic and at 2·7 (95% CI 2·5–4·1) for the 2000 Uganda epidemic. Our study allowed quantifying transmission in different settings during the two epidemics. According to our estimates, the term of R 0 concerning the transmission during traditional burial was estimated at 1·8 (95% CI 0·0–2·3) for the DRC epidemic and at 0·1 (95% CI 0·0–3·2) for the Uganda epidemic. For the Uganda epidemic, transmission in the community seems to have played an important role. Although confidence intervals are wide, these results suggest different roles of community, hospital and burial-related transmission in the two epidemics studied here. The term of R 0 associated with traditional burial increases with the CFR (see Appendix). This may explain a slightly more important role of funerals in the spread of the disease during the 1995 DRC epidemic, as the CFR was greater in this epidemic. However, a higher reproduction rate during burial may also indicate less precautions or increased contacts with cadavers at this time.
After interventions, neglecting depletion of susceptible individuals, we found that the effective reproduction number dropped to 0·4 (95% CI 0·3–0·6) for the DRC epidemic and 0·3 (95% CI 0·2–0·4) for the Uganda epidemic. This meant a large decrease in the point estimate of the transmission rate in the community in Uganda, where point estimates of hospital and burial transmission rates were small from the beginning. For the DRC epidemic, setting the transmission rate during burial to zero after interventions allowed an important decrease of effective reproduction number. Why intervention in the community appeared much more efficacious in Uganda remains speculative.
Second, we performed a multivariate sensitivity analysis of the model in order to identify the most important parameters for the control of the epidemic. It appears that for both pattern of transmission identified (DRC, 1995 and Uganda, 2000), the time to intervention, the hospitalization rate and the mean time between onset and hospitalization after institution of control interventions were related to epidemic size. Thus, the size of the epidemic could be reduced further by reinforcing interventions like contact tracing which could allow rapid hospitalization of cases after they develop the first symptoms. These results were obtained assuming that the transmission coefficient during hospitalization and burial would decrease by at least 50% for hospitalization and 75% for burial. In case of a non-identified nosocomial source of transmission, this assumption of the efficacy of the interventions at hospital would probably be too optimistic. For the transmission pattern identified with data from the 1995 DRC epidemic, the efficacy of interventions at hospital and after death was also a key parameter for the control of the epidemic size.
When we assumed that the transmission coefficient during burial would decrease at least by 50% instead of 75% the mean size of the epidemic obtained with the uncertainty analysis could reach 40 000 cases with the transmission pattern of the 1995 DRC epidemic. However, this does not take into account behavioural changes in the contact process and reinforcement of control interventions that would probably occur in such a case.
Since we fitted the model to available data on reported symptom onset dates (291 on 315 cases for the DRC epidemic and 418 on 425 cases for the Uganda epidemic) this may have led to an underestimation of R 0. At the beginning of the epidemic, there may have been additional cases that went unreported.
We applied our model to the two largest epidemics with sufficient data and did not study the previous smaller epidemics. However, using a stochastic version of the model, we illustrate the fact that a unique set of parameters can lead to epidemics of various sizes and notably small epidemics. Thus, our estimates of R 0 are also compatible with the occurrence of smaller outbreaks.
We assumed that all cases were related to human-to-human transmission except the first cases. Estimations of contact rates or efficacy of interventions in the community setting could be biased if a large proportion of cases were the result of routes or modes of transmission, other than the human-to-human transmission. Moreover, we did not differentiate between transmission through contact with patients or their body fluids and inoculation by contaminated instruments, but needle-stick contamination was not documented for the two epidemics we studied. We assumed that after the control interventions were put in place, there was no transmission at hospital or after death of patients. During the DRC epidemic, only three health-care workers developed EHF after barrier-nursing procedures were initiated. It is possible that two of them became infected prior to the arrival of the intervention team and the introduction of safety precautions [Reference Guimard30].
The type of model we used assumes homogeneous mixing within the population (even during hospitalization at the beginning of the epidemic and for the traditional burial), which may be too simplistic, notably in countries, where the structure of the community favours infections in households. Considering the effect of social networks and exploring the potential occurrence of super-spreading events could be an interesting issue for future research. In a recent paper, Lloyd-Smith et al. show the importance of super-spreading events in the spread of epidemics, notably for SARS and measles [Reference Lloyd-Smith39]. This may be an important issue for Ebola, but lack of available data prevents any clear conclusion.
Our estimates of R 0 are within the range (1·34–3·65) of the two other works in the field [Reference Chowell27, Reference Ferrari, Bjornstad and Dobson28]. In addition, two risk assessment studies have examined the risk related to contacts with Ebola cases during different states of illness [Reference Dowell20, Reference Francesconi40]. A risk factor study on modes of transmission prior to the institution of barrier precautions and other public health measures during the 1995 epidemic in Kikwit, DRC showed a pattern of increasing risk with exposures to patients in the later phases of illness. This study also showed significant association between contracting the disease and touching a cadaver [Reference Dowell20]. A retrospective risk factor assessment of cases and their contacts during the outbreak in Uganda showed significant association between the disease and contact with patients during illness. On the contrary, association between the disease and contact during body burial was not significant and point estimates of the prevalence proportion ratios were close to 1 [Reference Francesconi40]. Our estimates of the role of the community, hospital and traditional burial are consistent with these risk factor studies. Traditional burial ceremonies may have played a more important role in DRC and less so in Uganda where community transmission may have been a more significant source of infection. As our confidence intervals are wide, the results presented here present general trends in the relative components of transmission.
Our results show that rapid implementation of interventions and, when barrier nursing and isolation wards are efficient, rapid hospitalization of cases are key factors for the control of Ebola epidemics.
APPENDIX. Formulae for R 0
The formula for R 0 is the spectral radius of the next generation matrix of the following system:
Following the method described in van den Driessche [Reference van den Driessche and Watmough37] we determined the expression for R 0.
where
γi−1 is the mean duration of the infectious period for patients who survived to their illness, γd−1 is the mean duration of the infectious period for patients who died, γf−1 is the mean duration of the infectious period between death and burial, γh−1 is the mean duration between onset of symptoms and hospitalization, θ 1 is computed in order that θ% of infectious cases are hospitalized (see Table 3), δ1 and δ2 are computed to obtain a case fatality ratio at δ (see Table 3), N is the size of the population, βI is the transmission rate in the community, βH is the transmission rate after hospitalization, βF is the transmission rate during traditional burial.
ACKNOWLEDGEMENTS
This programme is partially funded by the French Armament Procurement Agency.
DECLARATION OF INTEREST
None.
