A Model of the 2014 Ebola Epidemic in West Africa with Contact Tracing

A differential equations model is developed for the 2014 Ebola epidemics in Sierra Leone, Liberia, and Guinea. The model describes the dynamic interactions of the susceptible and infected populations of these countries. The model incorporates the principle features of contact tracing, namely, the number of contacts per identified infectious case, the likelihood that a traced contact is either incubating or infectious, and the efficiency of the contact tracing process.The model is first fitted to current cumulative reported case data in each country. The data fitted simulations are then projected forward in time, with varying parameter regimes corresponding to contact tracing efficiencies. These projections quantify the importance of the identification, isolation, and contact tracing processes for containment of the epidemics.

The compartments II(t) and R(t) de-couple from the other compartments, and their values can be obtained from S(t), E(t), I(t), C(t). A schematic diagram of the model is shown in figure 1. The system of differential equations for S(t), E(t), I(t), C(t) iṡ rate of progression to infectiousness − κ(αI(t) + ψC(t))π E ω E E(t) N removal of exposed due to contact tracinġ We note that the transmission terms (involving β and ) and the contact tracing terms (involving κ) are of mass-action form. We note also that the form of all the loss terms in the equations assures that the solutions remain nonnegative for all time.
A major goal of our study is to fit the model to current reported data for Sierra Leone, Liberia, and Guinea. We note that the data available are the cumulative clinical reported cases [15,16], that is, the suspected cases, probable cases and confirmed cases according to the definitions given in WHO [17]. Therefore, if we denote by CU M (t) the cumulative reported cases at time t, then at time t + ∆t we have Hence, in order to fit the data we will use the above equation (2) with ∆t = 1 day and CU M (0) = the initial cumulative reported cases. We note that to simplify the model, we have omitted the infection dynamics of hospital and healthcare workers. Although healthcare worker infections from patients is of great importance and requires major attention, the contribution of healthcare worker infections to new transmissions and to the cumulative reported cases (2) is relatively small.

The Parameters of the Model
The parameters of the model are given in Table 1. The values of the parameters β, , ψ and α are estimated for the three countries using a least square curve fitting algorithm. The parameters σ and ν are taken to be values suggested by references in Table 1. The parameter γ is assumed to be a value such that the case mortality rate (outside hospital) is approximately 80%. In fitting the model to the data, we assume that contact tracing does not occur (κ = 0) during the time period of the data (up until September 23, 2014). The reasons for this assumption are that contact tracing has been insufficient in the three countries and to reduce the number of parameters to estimate. The main goal in incorporating contact tracing is to project forward how effective contact tracing can affect the future number of cases. The basic reproduction number of the model (1) is given by the following formula, computed by the next generation method [13] (Appendix): .

S(t)
The number of susceptible individuals at time t E(t) The number of exposed (incubating and not yet infectious) individuals at time t I(t) The number of infectious individuals at time t C(t) The number of deceased improperly handled at time t Q(t) The number of susceptible individuals under quarantine at time t II(t) The number of infectious individuals under isolation at time t R(t) The number of infected recovered or properly handled deceased at time t N Total national population (assumed to be constant) β Transmission rate excluding improper handling of deceased [1], [9], [11] Transmission rate due to improper handling of deceased [9], [11] κ Average number of contacts traced per identified/isolated infectious individual 1/α Average time from symptoms onset to identification/isolation of infectious individuals independent of contact tracing [1], [7], [9], [11] π I Probability of a contact traced infectious individual is isolated without causing a new case ω I Ratio of probability that contact traced individual is infectious at time of originating case identification to the probability a random individual in the population is infectious π E Probability a contact traced exposed individual is isolated without causing a new case ω E Ratio of probability that contact traced individual is exposed at time of originating case identification to the probability a random individual in the population is exposed 1/γ Average time from symptoms onset to recovery [1], [7], [9], [ The R 0 values we obtain are similar to the R 0 values obtained in [1], [3], [4], [5], [12], [18]. Notice that R 0 does not contain the contact tracing parameters κ, π E , π I , ω E , ω I . We note that R 0 is an incomplete indicator of the epidemic outcome. Even if R 0 > 1, the epidemic will ultimately be eliminated (Appendix), and if R 0 < 1, the number of cases may increase. The value of R 0 is as an approximate measure of the influence of various parameters on the number of secondary cases caused by an infectious individual, which does not account for future infections caused by these secondary cases that may be prevented by contact tracing. Thus R 0 does not provide a complete description of how the parameters affect the epidemic outcome. The removal (isolation) of infectious individuals due to contact tracing can be derived as follows: Let dI(t) = the number of infectious individuals removed due to contact tracing in the time interval (t, t + dt). Then The removal of incubating infected individuals due to contact tracing is derived similarly, and is discussed further below. The underlying assumption is that the probability p E (t) (p I (t)) a contact traced individual is incubating (infectious) at time t is proportional to the probability E(t)/N (I(t)/N ) that a random individual in the population is incubating (infectious) at time t. The proportionality constants ω E and ω I are intrinsic to the infection process of family members, and others closely associated with an originating source case from which the contacts are traced. The contacts' infection likelihood is very different from the general infection likelihood of random individuals in the population. The probabilites p E (t), p I (t) can, in principle, be ascertained from records of contact traced individuals. The probabilities E(t)/N and I(t)/N can be obtained from the solutions of the model. We form the ratios and view ω E and ω I as these proportionality constants. Equivalently, we can view p E (t) as ω E times as likely as E(t)/N and p I (t) as ω I times as likely as I(t)/N . The parameter ω E may be relatively large, since a traced contact will have a much greater chance of being infected than a random individual. For example, suppose that p E (t) is found to be 1/10 in a given time interval and the number of incubating individuals is E(t) = 200. Then, if the population is N = 4 × 10 6 , Further, the ratio ω E /ω I ≈ average incubation time = 1/σ average time from symptoms onset to identification = 1/α which is approximately 2 in our typical parameterizations. We note that the parameter ψ could also be included in this estimation, which would decrease this ratio somewhat. The parameters π E and π I measure the efficiency of the tracking, monitoring, and removing of incubating and infectious contacts. In particular, π E measures how fast public health workers remove incubating individuals upon symptoms onset, and prevent secondary transmissions (so that these individuals are effectively removed at the end of their incubation phase and transition to the II compartment without causing a new case). We note that if t 1 is the time contact tracing begins, then the cumulative number of cases reported between time 0 and time t < t 1 is (αI(s) + ψC(s))ds, and the cumulative number of cases reported at time t > t 1 is Since the total cumulative number of cases, both reported and unreported, at time t is N − S(t), the cumulative number of unreported cases can be determined at any time t > 0. The entire contact tracing process is highly dependent on public health resources, and varies greatly in different locations and epidemic stages. For example, in Sénégal the following policy has been implemented: 1) each identified patient is questioned in order to obtain a complete list of contacts; 2) the contacts are traced; 3) each contact is asked to stay at home; 4) each day, for 21 days, a healthcare worker visits the contacts and verifies whether or not the contacts are showing symptoms.
These protocols are rigorous and have been successful in preventing new cases in Sénégal. On October 17, 2014 the World Health Organization declared the end of the outbreak of the Ebola epidemic in Sénégal (after 42 days with no new cases and with active surveillance demonstrably in place and supported by good diagnostic capacity) [19] .

A More General Model
The model presented is directly applicable to the current epidemics in Sierra Leone, Liberia, and Guinea. We discuss here a more general model that incorporates additional features of Ebola epidemics applicable to other locations. In order to simplify the model presented here, we have neglected several features of contact tracing in the susceptible equation, which would involve a compartment of monitored susceptibles. There are several reasons for this simplification. First of all, strict isolation or monitoring of contact traced individuals not showing symptoms has been problematic in West Africa during this outbreak. Second, accounting for such isolation or monitoring would likely be more complicated than simply removing individuals from the susceptible compartment. Indeed, contact traced susceptibles are in a higher risk group for acquiring infection, since they are likely to be in contact with exposed individuals. For example, family members of an infected individual will continue to be in contact with each other, putting them at higher risk for a "second-order transmission" from another infected family member. Such second-order effects will be considered in a more general model in future work. This more general model involves low risk and high risk susceptibles, S 0 (t) and S 1 (t), respectively, as well as an added monitored compartment M (t). Here M (t) represents contact traced susceptibles sufficiently well-monitored such that they do not cause secondary infections (if infected). Note that the compartment M (t) should include health care workers who have had direct contact with infectious cases. In addition, all contact traced susceptible individuals are assumed to be high-risk. Further, the infection rate of high-risk susceptible individuals who are not well-monitored is assumed to be a linear term in the S 1 equation. A schematic diagram of the more general model is shown in figure 15. The new equations and modified transmission events are as follows: rate of return to low-risk for monitored susceptibleṡ Parameter/Variable Description S 0 (t) The number of low risk susceptible individuals at time t S 1 (t) The number of high risk susceptible individuals at time t M (t) The number of susceptible individuals under effective monitoring at time t ρ Transition rate from low risk to high risk φ Transition rate from high risk to low risk 1/τ Average monitoring period for contact traced susceptible individuals δ Fraction of monitored individuals returned to high risk susceptible θ Infection rate due to high risk (of individuals not well-monitored) π S Probability a contact traced susceptible individual does not cause secondary cases (if infected while monitored) ω S Ratio of the probability a contact traced individual is a high-risk susceptible to the probability a random individual in the population is susceptible ξ Infection rate due to high risk of well-monitored individuals (second order infections) The I(t) and C(t) equations are the same as in (1). Also, the loss term ξM (t) transitions the infected (monitored) individuals into the II(t) compartment, which is decoupled from the equations in (3). Note that these individuals would be removed upon symptoms onset, and because they are well-monitored, they can be considered effectively removed in this transition term. It is possible to generalize the model further by tracking the disease stage of all infected individuals, particularly contact traced individuals. Such a model is best treated with an age of infection variable, which allows tracking of the incubation and infectious disease stages [6]. In future work we will develop these models, with age of infection as a continuum independent variable.

Simulations of the Ebola Epidemic in Sierra Leone
In figure 2 we fit the model without contact tracing to the cumulative reported case data for Sierra Leone from May 27, 2014 to September 23, 2014 (WHO [15,16]). The fit to data can be accomplished with varying combinations of parameters. Here we have used a least-squares algorithm to obtain a choice of parameters with relatively accurate fit. The parameters obtained in the fit yield a basic reproduction number of R 0 = 1.26. The simulation yields the following information about the epidemic on September 23, 2014 (day 119): the ratio of exposed cases to infectious cases is E(119)/I(119) ≈ 2.49; the ratio of improperly handled deceased cases to infectious cases is C(119)/I(119) ≈ 0.57; and the ratio of cumulative reported cases to cumulative unreported cases is ≈ 1.78. These ratios, which are dependent on parameters, are relatively stable at the data end-stage. In figure 3 we add contact tracing to the model and predict the further evolution of the epidemic in Sierra Leone forward from September 23, 2014. The contact tracing parameters α and κ are varied in a sensitivity analysis, while the other contact tracing parameters are held constant. The graphics reveal that a general identification/isolation rate α > 0.3 is required for containing the epidemic. The number κ of contacts traced per identified case is also important if α is smaller. After contact tracing begins and for a short time, the reported cases increase as κ increases, but then the epidemic subsides as contact tracing takes effect, and the reported cases decrease as κ increases.

Simulations of the Ebola Epidemic in Liberia
In figure 4 we fit the model without contact tracing to the cumulative reported case data for Liberia from June 17, 2014 to September 23, 2014 (WHO [15,16]). We have used a least-squares algorithm to obtain a choice of parameters with relatively good fit (other parameter choices will give similar fits). The parameters obtained in the fit yield a basic reproduction number of R 0 = 1.54. The simulation yields the following information about the epidemic on September 23, 2014 (day 98): the ratio of exposed cases to infectious cases is E(98)/I(98) ≈ 3.35; the ratio of improperly handled deceased cases to infectious cases is C(98)/I(98) ≈ 0.58; and the ratio of cumulative reported cases to cumulative unreported cases is ≈ 1.37. These ratios, again depend on parameters, are very stable throughout most of the period of simulation. In figure 5 we add contact tracing to the model and predict the further evolution of the epidemic in Liberia forward from September 23, 2014. The contact tracing parameters α and π E are varied in a sensitivity analysis, while all the other contact tracing parameters are held constant. The graphics again reveal that an identification/isolation rate α > 0.3 is required for containing the epidemic. The role of π E , as the probability of efficiently tracing and monitoring an incubating contacted individual, who is infected, but not yet infectious, is also important for the containment of the epidemic. As in figure 3, the reported cases increase for a short time as π E increases, but then the decrease as π E increases, as contact tracing takes effect. In figure 8 we illustrate the nonlinear effect of varying π E against the cumulative projected cases forward 100 days.

Simulations of the Ebola Epidemic in Guinea
In figure 6 we fit the model without contact tracing to the cumulative reported case data for Guinea from March 25, 2014 to September 23, 2014 (WHO [15,16]). We have used a least-squares algorithm to optimize a choice of parameters for the fitting (other choices will yield a similar fit). The fit is problematic in the middle range of the time period, but the data is very erratic, and the cumulative count even decreases some days. The parameters obtained in the fit yield a basic reproduction number of R 0 = 1.11. The simulation yields the following information about the epidemic on September 23, 2014 (day 182): the ratio of exposed cases to infectious cases is E(182)/I(182) ≈ 2.56; the ratio of improperly handled deceased cases to infectious cases is C(182)/I(182) ≈ 0.66; and the ratio of cumulative reported cases to cumulative unreported cases is ≈ 4.95. These ratios, again depend on parameters, and again are very stable throughout most of the period of simulation. In figure 7 we add contact tracing to the model and predict the further evolution of the epidemic in Guinea forward from September 23, 2014. The contact tracing parameters α and π I are varied in a sensitivity analysis, while all the other contact tracing parameters are held constant. The graphics again reveal that the general identification/isolation rate α > 0.3 is of great importance in containing the epidemic. The efficiency of the contact tracing parameter π I is of lesser importance in these simulations than the parameter π E , especially when α > 0.3 (compare to figure 5). The probability π I of efficiently tracing and isolating an infectious contacted individual is compensated by the general efficiency of isolating and removing an infectious individual (independent of contact tracing). The probability π E of identifying and isolating an incubating contacted individual so that they cause no new infections is dependent solely on the contact tracing process and its efficiency in monitoring and isolation of these incubating individuals upon symptoms onset. After contact tracing begins, the cumulative reported cases and cumulative total cases decrease as π I increases without the switch-over observed in figures 3 and 5. Figures 9,10,11,12,13, and 14 show 100 stochastic simulations compared with the ODE solution for contact tracing in Sierra Leone with different rates of case identification α. The stochastic simulations are generated by simulating a Continuous Time Markov Chain as a continuation of the ODE solution beginning at the last data time point. The contact tracing parameters are assumed to be ω E = 2000, ω I = 1000, π E = 0.5, π I = 0.8 and κ = 20 in Figures 9, 10, and 11. In Figures 12 and 13, we assume κ = 0 (no contact tracing), while κ = 10, π E = 0.1, π I = 0.5 (less effective contact tracing) in Figure 14. All other parameters are taken to be the same as in the previous fit to the Sierra Leone epidemic. The averages of the stochastic model solutions agree with the ODE solutions of (1).

Summary and Conclusions
We have developed a model of Ebola epidemics in West Africa that focuses attention on the elements of public health policies for containment of these epidemics. Our simulations for Sierra Leone, Liberia, and Guinea fit the cumulative reported cases for these countries up to September 23, 2014, and project future epidemic levels forward from September 23, 2014 (based on various parameterizations corresponding to these elements). Our projections indicate that the most important elements for containment of the epidemics within a relatively short time span are that (1) infectious cases (independent of contact tracing) are efficiently reported and isolated, with the average time between the appearance of symptoms and isolation less than 3 days (α > 0.3); (2) contact traced incubating infected cases are efficiently monitored, with average probability of compliance, with isolation upon appearance of symptoms (such that no new cases are caused by individual), greater than 0.5 (π E > 0.5).
Also of importance in mitigation of the epidemics is a reduced rate at which infected deceased are improperly handled (ψ, ν), a sufficient number of contacts traced per identified infectious individual (κ), and an efficient identification and isolation of contact traced infectious individuals (π I ). The model allows quantification of the parameters corresponding to public health controls (α, ψ, κ, π E , π I , ω E , ω I ) for evaluating the impact of public health policies for the evolution of these epidemics.