1. Introduction

Since its discovery in 1976 in the country of Zaire (now the Democratic Republic of the Congo), four primary strains of the Ebola virus disease (EVD) have been identified: Zaire ebolavirus, Sudan ebolavirus, Reston ebolavirus, and Côte d'Ivoire ebolavirus [19]. The Zaire ebolavirus is the strain of the 2014 West Africa outbreak and the strain that has caused the deadliest recorded EVD outbreaks [41]. Figure 1 classifies the 19 outbreaks of EVD with more than 5 reported cases (excluding the 2014 West Africa Outbreak) by mortality rate [12]. Zaire ebolavirus is the strain responsible for all of the outbreaks in the highest mortality rate category (80-100%) and four out of five of the outbreaks in the second highest mortality rate category (60-79%).

Ebolavirus can only be transmitted through direct contact with the bodily fluids of an infectious individual or body [22]. An individual infected with ebolavirus will not become infectious or show symptoms for several days (range: 2-21 days, mean $\pm$ std: $6.3\pm 3.31$ days, median: 5.5 days) [20,45,40,25]. After the incubation period the individual becomes infectious and starts to show symptoms, which may include fever, vomiting, hemorrhagic bleeding, and multiple organ dysfunction which can lead to death [20,50,11]. During the 2014 West Africa EVD outbreak, it was discovered that even after a male recovers, the virus can remain in his semen anywhere from three to nine months after the onset of symptoms and can be transmitted to sexual partners [18,15]. Additionally, an individual who has died due to EVD remains infectious for a period of up to seven days [51] and, if not handled properly prior to burial or cremation, can provide another route of virus transmission [33,50,51].

---

1.1. Timeline of the 2014 West Africa Outbreak

While the 2014 West Africa EVD outbreak originated in the country of Guinea during December of 2013 [8], Sierra Leone and Liberia were also greatly affected [20,9,67]. A timeline of major outbreak response plan implementations during the 2014 West Africa EVD outbreak is shown in Figure 2.

The timeline in Figure 2 divides the outbreak into five disjoint segments according to the response plans implemented over the course of the epidemic. At the beginning of the outbreak, there was no response plan in place because the World Health Organization (WHO) did not declare the outbreak severe enough to require immediate support until mid July of 2014. The initial response plan (WHO Phase 1), implemented July 31,2014, by the WHO, was intended to stop the transmission of EVD in the most affected countries (Sierra Leone, Liberia, and Guinea), as well as provide neighboring countries with support to prevent further spread of the disease [66,65,63]. Following this first-response plan, the United Nations (UN) created the UN Mission for Ebola Emergency Response (UNMEER) plan on September 19,2014 [58,63]. The UNMEER plan added to the efforts of the WHO plan already under way and was intended to isolate EVD cases, ensure safe burials for as many EVD deaths as possible, assist in facilitating case management, and educate individuals in identifying EVD symptoms and going to hospitals when infected [58]. By mid November 2014, the implementation of UNMEER was underway [58,63]. In January of 2015, the WHO shifted the focus of their plan to increasing their capacity for identifying and tracing cases and their contacts; this shift in focus has been termed ''Phase 2'' of the WHO Ebola Response Plan [63]. The UNMEER response plan continued in conjunction with Phase 2 of the WHO response plan. In July 2015, the UNMEER response plan ended [58], and the WHO Ebola Response Plan shifted to Phase 3. This phase focused on stopping all remaining sources of EVD transmission by working in conjunction with international, national, and community partners to improve the speed of case identification, establish and maintain health facilities, and improve Ebola survivor support [63].

---

1.2. Previous EVD models

Several mathematical and computational models simulating the spread of EVD through a population have been proposed. These models have been used to explore a variety of questions about EVD outbreaks and to estimate rates of transmission and the basic reproduction number of particular outbreaks. Previously proposed models have noted the importance of EVD outbreak models including features like modes of transmission from deceased infectious [59,52,50,21,53], explicitly defining distinct susceptible groups for both high and low risk of exposure to EVD [35], differentiating the infectious in hospitals from the infectious not in hospitals [7,52,50,60], and defining healthcare workers as distinct from other susceptible individuals [50]. However, none of these previous models have included all four of these features. Several models simulate the effects of possible intervention strategies and control measures such as isolation, quarantine, and social intervention to reduce transmission rates during an outbreak [52,35,2,50,40,21,53]. However, the 2014 EVD outbreak hastened the development of a prophylactic vaccine [32] the impact of which has yet to be fully explored.

The clinical trials of the recently developed EVD vaccine used a ring vaccination strategy [31]. The ring vaccination strategy identifies and vaccinates contacts of infectious individuals, as well as contacts of contacts. Recent studies have used a network-based transmission agent-based model [14] and spatially explicit agent-based models [1,46] to examine the impact of ring vaccination of EVD outbreaks using data from the 2014 West Africa Outbreak. The explicit inclusion of space within the model follows the modeling approach of ring vaccination models for foot and mouth disease [47,37].

Many of the previous EVD models are used to estimate the basic reproduction number ($\mathcal{R}_0$) of the 1995 Congo and 2000 Uganda EVD outbreaks [13,40] and the 2014 West Africa EVD outbreak [59,2,50,35,26,53,60]. While the basic reproduction number is a useful estimate of the severity of an outbreak, the estimated value of $\mathcal{R}_0$ is only valid using parameter values representing the beginning of an outbreak. If parameter values like transmission rates change over the course of an outbreak due to the implementation of control measures, then a time-dependent reproduction number is more appropriate. Two possibilities are the actual reproduction number ($\mathcal{R}_A(t)$) which is calculated using incidence and prevalence data [3], and the effective reproduction number ($\mathcal{R}_E(t)$) which is calculated using the size of the susceptible population and the value of the basic reproduction number ($\mathcal{R}_0$) [17]. The effective reproduction number has been calculated for the 1995 Congo and 2000 Uganda EVD outbreaks by Chowell, et al. [13], and for the 2014 EVD outbreak by Browne, et al. [7]. However, the lack of key EVD outbreak features in the Chowell and Brown models, specifically transmission from deceased infectious individuals, leaves room for improvement.

---

2. Model description and assumptions

We have developed an ordinary differential equations model that simulates an Ebola epidemic within a population and takes into account: (1) transmission through contact with infectious EVD individuals and deceased EVD bodies, (2) the heterogeneity of the risk of becoming infected with EVD, and (3) the increased survival rate of infected EVD patients due to greater access to trained healthcare providers via an increase in healthcare workers within the population. We apply this model to the population of Sierra Leone during the 2014 EVD outbreak.

The total population is divided into 12 states (see Figure 3). There are three classes of susceptible individuals: healthcare workers who have a high risk of exposure ($S_W$), members of the general population (i.e., not healthcare workers) with a high risk of exposure ($S_H$), and members of the general population with a low risk of exposure ($S_L$). Low risk general population members become high risk based on the rate at which an average individual's contacts become infected, which is determined by the density of infected EVD patients and dead bodies. High risk members of the general population return to low risk by the rate at which an infectious contact ceases to be infectious either by recovery or the burial of the deceased individual's body.

Once an individual becomes infected, they move to the exposed class ($E$) where they are asymptomatic and not infectious to others. When an exposed individual becomes symptomatic, they become infectious to susceptible individuals and move from the exposed class to either the infectious class ($I$) or infectious in hospitals class ($\hat{I}$). We assume individuals in the $\hat{I}$ class can only transmit the virus to healthcare workers (i.e., the $S_W$ class) and not the general susceptible population (i.e., the $S_L$ and $S_H$ classes), thus effectively quarantining the $\hat{I}$ from a large portion of the population. An infectious individual will either die and move to the deceased class ($D$, or $\hat{D}$ for those $\hat{I}$ individuals who died in hospitals), or they will recover and move to recovering but still infectious class ($R_I$). The bodies of individuals who have died remain infectious until they are buried ($B$). Individuals in the $R_I$ class no longer show symptoms of EVD, and are recovering from the disease but are still infectious (with a lower transmission probability) until they move to the recovered non-infectious class ($R$) after several months [15,18]. We assume that individuals in the $R_I$ class will not die from the ebolavirus. Furthermore, since individuals in the $R_I$ class no longer show EVD symptoms and ''appear'' to have recovered, we assume that their close contacts are no longer considered highly susceptible. Thus, individuals in the $R_I$ class only infect individuals in low susceptible class ($S_L$). Individuals in the recovered and non-infectious class ($R$) are assumed to have immunity from further ebolavirus infections. This assumption is made based on evidence that recovery from an EVD infection provides immunity for at least 10 years [35,10], and our model is designed to simulate the 2014 EVD outbreak, which lasted about 2 years.

Given the availability of a prophylactic vaccine, our model allows for the vaccination of susceptible individuals, effectively moving them to the recovered and non-infectious class ($R$). Vaccinated healthcare workers are moved into a separate class ($V_W$) where they continue their role as healthcare workers in the model.

Background death due to natural causes occurs in all three susceptible classes ($S_L$, $S_H$, $S_W$), the $R$ and $R_I$ classes, and the vaccinated healthcare worker class ($V_W$). Births into the population are assumed to only occur in the low susceptible class ($S_L$). Susceptible healthcare workers are the only individuals allowed to migrate into and out of the population.

The model equations are given in System (1). See Table 1 for a description of each of the model's parameters.

$\begin{align} S_L'(t) &= \Omega -\frac{\beta_{1}^{L}I + \beta_{2}^{L}D + \beta_{3}R_{I}}{N}S_{L} -\varepsilon \rho_L S_{L} + \sigma_{L}S_{H} -\sigma_{H}\frac{I+\hat{I}+D+\hat{D}}{N}S_{L} -\mu S_{L} \label{eq1} \end{align}$

(1.1)

$\begin{align} S_H'(t) &= -\frac{\beta_{1}^{H}I + \beta_{2}^{H}D}{N}S_{H} -\varepsilon \rho_H S_{H} + \sigma_{H}\frac{I+\hat{I}+D+\hat{D}}{N}S_{L} -\sigma_{L}S_{H} -\mu S_{H} \label{eq2} \end{align}$

(1.2)

$\begin{align} S_W'(t) &= \psi (\phi + I+\hat{I}+D+\hat{D}) -\frac{\beta_{1}^{W}\hat{I} + \beta_{2}^{W}\hat{D}}{N}S_{W} -\varepsilon \hat{\rho}S_{W} -\mu S_{W} \label{eq3} \end{align}$

(1.3)

$\begin{align} E'(t) &= \frac{(\beta_{1}^{H} I + \beta_{2}^{H}D)S_{H} + (\beta_{1}^{L}I + \beta_{2}^{L}D + \beta_{3}R_{I}) S_{L} + (\beta_{1}^{W}\hat{I} + \beta_{2}^{W}\hat{D})S_{W}}{N} -\kappa E \label{eq4} \end{align}$

(1.4)

$\begin{align} I'(t) &= \kappa(1-\alpha)E -\gamma I \label{eq5} \end{align}$

(1.5)

$\begin{align} \hat{I}'(t) &= \kappa \alpha E -\hat{\gamma} \hat{I} \label{eq6} \end{align}$

(1.6)

$\begin{align} R_I'(t) &= \gamma(1 -\delta) I + \hat{\gamma} \left(1 -\delta \nu \frac{\phi}{\phi + S_{W} + V_W}\right) \hat{I} -\eta R_{I} -\mu R_{I} \label{eq7} \end{align}$

(1.7)

$\begin{align} R'(t) &= \varepsilon (\rho_H S_{H} + \rho_L S_{L}) + \eta R_{I} -\mu R \label{eq8} \end{align}$

(1.8)

$\begin{align} D'(t) &= \gamma \delta I -\lambda D \label{eq10} \end{align}$

(1.9)

$\begin{align} \hat{D}'(t) &= \hat{\gamma}\left(\delta \nu \frac{\phi}{\phi + S_{W} + V_W}\right) \hat{I} -\hat{\lambda}\hat{D} \label{eq11} \end{align}$

(1.10)

$\begin{align} B'(t) &= \lambda D + \hat{\lambda}\hat{D} \label{eq12} \end{align}$

(1.11)

$\begin{align} V_W'(t) &= \varepsilon \hat{\rho} S_{W} -\mu V_{W} \label{eq9} \end{align}$

(1.12)

In System (1), $N(t)$ represents the sum of all 12 states at time $t$, that is

$\begin{align*} N(t)= & S_H(t)+S_L(t)+S_W(t)+E(t)+I(t)+\hat{I}(t)+R_I(t)+R(t)+D(t)+\hat{D}(t)\\ &+B(t)+V_W(t). \end{align*}$

Since transmission of the Ebola virus requires direct contact with bodily fluids, we assume transmission is frequency dependent; see Equations (1.1, 1.2, 1.3, and 1.4). In response to the increasing severity of the epidemic in West Africa, many countries around the world sent doctors and healthcare workers to help treat EVD patients [48]. Additionally, during the WHO's initial response, one of the primary objectives was to train individuals to care for EVD patients and attend to those who had died of EVD [63]. Thus, our model assumes that the migration rate of new healthcare workers into the population depends on the number of currently infectious live individuals and dead bodies at time $t$ during the outbreak ($I(t)+\hat{I}(t)+D(t)+\hat{D}(t)$) such that greater numbers of infectious correspond to higher rates of new healthcare workers entering the population [62]; see Equation (1.3). However, our model also assumes there would still be new healthcare workers entering the population even in the absence of an outbreak, albeit at a low rate. We additionally assume that a greater abundance of healthcare workers present in the population will correlate to a greater number of hospitalized patients recovering and thus inversely scale the death rate due to EVD; see Equations (1.7 and 1.10).

---

3. Initial conditions and model parameters

The initial time $t=0$ corresponds to April 27,2014 (see Figure 2), which is 30 days prior to the initial date the CDC began recording outbreak data. The initial population size, $N(0)$, is the total population of Sierra Leone (5,743,725 as of July 2014) [16]. This population is then distributed among the different classes of the model. The healthcare worker class ($S_W$) begins as 0.039% of the total population (an estimation based on the physician and nurse/midwife densities in Sierra Leone) [16]. The infectious, non-hospitalized class ($I$) contains 10 individuals, with five times that many individuals placed in the high susceptible class ($S_H$). The rest of the population begins in the low susceptible ($S_L$) class. All other classes contain no individuals at time $t=0$.

The values for parameters $\mu$, $\Omega$, $\kappa$, $\delta$, $\gamma$, $\hat{\lambda}$, and $\eta$ were taken from clinical and demographic studies of EVD [15,18,45,50], the CDC [12], and the CIA World Factbook [16]. See Table 1 for the source of each parameter value. The burial rate of hospitalized EVD patients (${\raise0.5ex\hbox{$\scriptstyle 1$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle \hat{\lambda}$}}$) removes infectious dead bodies within the seven days in which they are still infectious [51]. The value of $\delta$ is estimated as the average of the death rates of Zaire ebolavirus outbreaks that had more than 5 infected individuals (as reported by the CDC [12]) not including the 2014 outbreaks.

The infectious period of individuals who go to hospitals with Ebola symptoms (${\raise0.5ex\hbox{$\scriptstyle 1$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle \hat{\gamma}$}}$) was estimated to be slightly faster than the known infectious period of those individuals who did not go to hospitals (${\raise0.5ex\hbox{$\scriptstyle 1$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle \gamma$}}$). Additionally, the burial rate of individuals who do not go to hospitals ($\lambda$), was estimated to be slightly lower than for hospitalized EVD patients ($\hat{\lambda}$) because there were no regulations affecting the burials of those who died of EVD outside of hospitals. However, we assume that the burial rate of deceased individuals who were not hospitalized still removes infectious dead bodies within the seven days in which they are still infectious [51].

The parameters $\sigma_H$ and $\sigma_L$, the rates at which individuals move between the low and high susceptible classes, were estimated using parameters that already had known literature values. We assume the rate at which individuals move from low to high exposure ($\sigma_H$) is the same rate at which individuals become symptomatic and infectious ($\kappa$). We assume the rate at which individuals move from high to low exposure ($\sigma_L$) is equal to the inverse of the sum of the number of days it takes individuals who did not go to hospitals to recover or die (${\raise0.5ex\hbox{$\scriptstyle 1$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle \gamma$}}$) and the number of days it takes for those individuals who die from EVD not in hospitals to be buried (${\raise0.5ex\hbox{$\scriptstyle 1$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle \lambda$}}$); that is, $\sigma_L = 1/\left({\raise0.5ex\hbox{$\scriptstyle 1$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle \gamma$}} + {\raise0.5ex\hbox{$\scriptstyle 1$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle \lambda$}}\right)$.

The migration rate of healthcare workers into the population due to the outbreak ($\psi$) varied by response plan. It is assumed that when no response plan had been implemented, there was less influx of healthcare workers ($\psi = {\raise0.5ex\hbox{$\scriptstyle 0.0005$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 7$}}$, an increase of 0.0005 healthcare workers per infectious individual per week). However, when response plans were put into place, the rate at which healthcare workers entered the population increased ($\psi = {\raise0.5ex\hbox{$\scriptstyle 0.002$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 7$}}$, an increase of 0.002 healthcare workers per infectious individual per week).

The transmission rates within our model are a product of the rate at which susceptible individuals come into contact with infectious individuals and the probability that a single infectious contact will transmit the virus to a susceptible individual. Unique transmission rates are used for each susceptible class and for each mode of transmission (i.e., from infectious individuals and infectious dead bodies). See Table 1 for the definition of each transmission rate and Table 3 for transmission rate values used. The implementation of the response plans resulted in behavioral changes that altered contact rates and thus transmission rates. Using incidence data from the CDC [9] and System (1), we have estimated the transmission rates $\beta_1^H$, $\beta_2^H$, $\beta_1^L$, $\beta_2^L$, $\beta_1^W$, $\beta_2^W$, and $\beta_3$ for each response period. Additionally, the proportion of individuals who go to hospitals once infected ($\alpha$), as well as the effectiveness of the healthcare workers in reducing deaths due to EVD ($1-\nu$), changed over the different response periods. Thus we have also fitted the values of $\alpha$ and $\nu$ for each response period. Note that the smaller the value of $\nu$, the greater the reduction in deaths due to EVD in hospitals.

---

3.1. Fitting $\beta$, $\alpha$, and $\nu$ values

For each response period, $n=5000$ possible parameter sets were generated using Latin Hypercube sampling (LHS) over the uniform distributions for the $\beta$, $\alpha$, and $\nu$ parameters given in Table 3 (each column represents a single response period). See [6,5] for a detailed description and example of LHS. For each parameter set for the first response period, $t\in[0,95]$, System (1) is used to simulate the EVD outbreak in Sierra Leone from April 27,2014 to July 31,2014. For a given parameter set $i$, where $i = 1, \dots, n$, we define $CI_i(t)$ and $CD_i(t)$ as the cumulative number of infections and deaths, respectively, up to time $t$ as estimated by System (1), where

$CI_i(t) = \int_0^t \kappa E(u) du \text{and} CD_i(t) = \int_0^t \gamma\delta I(u) + \hat{\gamma}\delta\nu\frac{\phi}{\phi+S_W(u)+V_W(u)} du.$

The parameter set that best describes the outbreak over this time period minimizes the weighted error

$h_i=\sum\limits_{t\in\tau} \sqrt{w(t)\left[\left(CD_i(t) - CD^*_t\right)^2 + \left(CI_i(t) - CI^*_t\right)^2\right]}, \label{eqn:errorfnc}$

(2)

where $CI^*_t$ and $CD^*_t$ are the cumulative number of infections and deaths, respectively, as reported by the CDC [9] at time $t$, and $\tau$ is the set of times at which the CDC provides data for the cumulative number of infections and deaths for a given time period. The weight function $w(t)$ in Equation (2) prioritizes the model more closely matching the data at the end of a time period, and is given by

$w(t) = 1-0.5e^{-a(t-b)}, \label{eqn:weightfnc} \nonumber$

where the values of $a$ and $b$ for each intervention period are given in Table 2.

Next, for each parameter set for the second response period, $t\in[95,202]$, System (1) is used to simulate the EVD outbreak in Sierra Leone from July 31,2014 to November 15,2014. Each simulation uses as its initial conditions the terminal conditions from the first response period using the parameter set that minimized the weighted error (Equation (2)). The cumulative number of infections and deaths are calculated over the time period for each simulation (i.e., each parameter set). The parameter set that best describes the outbreak over this second response period is the one that minimizes the weighted error given the CDC cumulative infection and death data for July 31,2014, through November 15,2014 [9].

The process is repeated for response periods $t\in[202,248]$, $t\in[248,460]$, and $t\in[460,647]$, each time using as the initial conditions the terminal conditions of the previous response period generated from the parameter set minimizing the weighted error (Equation (2)). The $\beta$, $\alpha$, and $\nu$ values found using this method are shown in Table 3. The number of data points for each time period (21 in $[0,95]$; 33 in $[95,202]$; 9 in [202,248]; 31 in $[248,460]$; 26 in $[460,647]$) is relatively low when compared to the number of parameters being fit (9 parameters), and thus this method cannot be guaranteed to uniquely identify the fitted parameters. However, in the absence of more data, we accept these parameter estimations as sufficient for predicting the cumulative infections and deaths of the 2014 EVD outbreak according to our model.

The results in Table 3 show that the proportion of infected individuals who seek medical treatment ($\alpha$) was higher in the time periods in which intervention measures were in place compared to the time period where there was no response. Furthermore, the value of $\alpha$ was substantially higher ($>65\%$) by the time both the WHO Phase 1 and UNMEER plans were in place (i.e., $t>202$). Table 3 also reveals that the transmission rates ($\beta$ parameters) did not decrease for all subpopulation as the response strategies were implemented. In fact, some $\beta$ values decreased during initial response periods and then increased for later response periods. Lastly, Table 3 shows that, during the response periods with intervention measures, heathcare workers were more effective at reducing deaths due to EVD than during the no-reponse period ($t\in[0,95]$). In fact, the parameter fitting suggests that the healthcare workers were most effective in reducing deaths due to EVD during the WHO Phase 1 intervention period ($t\in[95,248]$).

Figure 4 shows the cumulative number of infections and deaths over time as reported by the CDC (points) and as predicted by the model (solid curve) simulated over $t\in[0,647]$ using the parameter sets for each response period that minimizes the weighted error (Equation (2)).

---

4. Reproduction numbers

Epidemic reproduction numbers provide a measure of the severity of an outbreak of an infectious disease. The most commonly calculated reproduction number is the basic reproduction number, $\mathcal{R}_0$, which measures the number of secondary infections that a single infected individual causes in a completely susceptible population [29]. Additionally, the derivation of $R_0$ using the next-generation method (see [23]) provides a threshold condition determining the stability of the disease-free equilibrium. Specifically, the introduction of an initial infection with conditions under which $\mathcal{R}_0 > 1$ leads to an epidemic, while outbreaks with $\mathcal{R}_0 < 1$ will die out. In Section 4.1, we derive and evaluate $\mathcal{R}_0$. However, since the estimated value of $\mathcal{R}_0$ is only valid using parameter values representing the beginning of an outbreak, and the $\beta$, $\alpha$, and $\nu$ parameters change with each response period, we calculate the time-dependent actual reproduction number, $\mathcal{R}_A(t)$ (see Section 4.2). In Section 4.3 we define the effective reproduction number and discuss why it is an inappropriate measure for this model of the 2014 EVD outbreak.

---

4.1. Basic reproduction number, $\mathcal{R}_0$

An expression for $\mathcal{R}_0$ in terms of model parameters is derived using the next-generation method [23], which determines the stability of the disease-free equilibrium. For the model described in System (1), the disease-free equilibrium is

$\begin{equation}S_H^* =E^* = I^* = \hat{I}^* = R_I^* = D^* = \hat{D}^*= B^* = 0\nonumber\end{equation}$

$S_L^* = \frac{\Omega}{\epsilon \rho_L + \mu}, S_W^* = \frac{\psi \phi}{\epsilon \hat{\rho} + \mu}, R^* = \frac{\Omega \epsilon \rho_L }{\mu (\epsilon \rho_L + \mu)}, V_W^* = \frac{\epsilon \hat{\rho} \psi \phi}{\mu (\epsilon \hat{\rho} + \mu)},\nonumber$

and the expression for $\mathcal{R}_0$ is

$\begin{array}{l}\mathcal{R}_0 &=&[\mu (\hat{\gamma} \hat{\lambda} (\beta_2^L \gamma \delta (\eta + \mu) + \lambda (-\beta_3 \gamma (\delta - 1) + \beta_1^L (\eta + \mu))) (\mu + \epsilon \rho_W) (\mu + \psi) \Omega + \nonumber \\ && \alpha (\gamma \lambda (\eta + \mu) (\mu + \epsilon \rho_L) \phi \psi (\beta_2^W \hat{\gamma} \delta \mu \nu + \beta_1^W \hat{\lambda} (\mu + \psi)) -\hat{\gamma} \hat{\lambda} (\mu + \epsilon \rho_W) \nonumber \\ && (\beta_2^L \gamma \delta (\eta + \mu) (\mu + \psi) + \lambda (\beta_3 \gamma \delta (\mu (\nu - 1) - \psi) + \beta_1^L (\eta + \mu) (\mu + \psi))) \Omega))] / \nonumber \\ &&[\gamma \hat{\gamma} \lambda \hat{\lambda} (\eta + \mu) (\mu + \epsilon \rho_L) (\mu + \epsilon \rho_W) (\mu + \psi) (\phi \psi + \Omega)]. \nonumber\end{array}$

Evaluating $\mathcal{R}_0$ using the parameter values found in Table 1 and the first column of Table 3, we find $\mathcal{R}_0 = 1.33$. This is consistent with estimates of $\mathcal{R}_0$ for the 2014 EVD outbreak in Sierra Leone using other models that estimated $1.26 < \mathcal{R}_0<2.53$ [2,35,59]. Since $\mathcal{R}_0>1$, the disease-free equilibrium is unstable and the initial outbreak will start to spread through the population, which is what occurred in Sierra Leone.

---

4.2. Actual reproduction number, $\mathcal{R}_A(t)$

The actual reproduction number for an outbreak, $\mathcal{R}_A(t)$, is defined as the average number of secondary cases caused by an individual who is infectious at time $t$ [3]. It is calculated using the prevalence and incidence of infections, and measures the extent of the disease in a population at any time during the outbreak [3]. The prevalence at time $t$, denoted $p(t)$, is the total number of infected individuals in the population at time $t$. For the EVD model described by System (1),

$p(t) = I(t) + \hat{I}(t) + R_I(t) + D(t) + \hat{D}(t). \label{eq:prevalence} \nonumber$

The incidence at time $t\geq1$, denoted $j(t)$ with $t$ measured in days, is the number of new infections generated over the course of one day, $[t-1,t]$. For System (1),

$j(t) = \int_{t-1}^t \kappa E(s) ds. \label{eq:incidence} \nonumber$

Give the incidence and prevalence, the actual reproduction number is defined as

$\mathcal{R}_A(t) = \xi\frac{j(t)}{p(t)},\label{eqn:RA}\nonumber$

where $\xi$ is the average duration of infectiousness [3]. Since System (1) contains multiple infectious classes and each class has a different length of infectiousness, $\xi$ is calculated as a weighted average of the infection length of each infectious class; $\xi$ is weighted by the proportion of individuals who enter each infectious class (i.e., the proportion who go to hospitals ($\alpha$) and the proportion who die ($\delta$)). Thus

$\xi =(1-\alpha)\frac{1}{\gamma}+\alpha\frac{1}{\hat{\gamma}}+\delta \left((1-\alpha)\frac{1}{\gamma}+\alpha\frac{1}{\hat{\gamma}} \right) +(1-\delta)\frac{1}{\eta}. \nonumber$

Figure 5 shows $\mathcal{R}_A(t)$ over the entire time horizon, where the fitted value of $\alpha$ for each time period (see Table 3) is used. Within the first response period we see a spike in the actual reproduction number. This is due to the fact that at the beginning of an epidemic the prevalence level will be low, but the incidence level can be relatively large, as was the case in the first several weeks of the 2014 EVD outbreak. It is not until the third response period that $\mathcal{R}_A(t)$ drops below 1. Thus it is not until the combination of the WHO Phase I response and the implementation of UNMEER that the average infectious individual is infecting less than one susceptible individual. The value of $\mathcal{R}_A(t)$ falling below 1 indicates that the epidemic would begin to die out, and, indeed, the points of inflection in the cumulative infection and death curves (see Figure 4) also occur within the third time period.

The average actual reproduction number for the $i^{th}$ time period $[t_i,t_{i+1}]$ is

$\bar{\mathcal{R}}_A^i = \frac{1}{t_{i+1}-t_i}\int_{t=t_i}^{t_{i+1}} \mathcal{R}_A(s) ds,$

where $t_0 = 0$, $t_1 = 95$, $t_2 = 202$, $t_3 = 248$, $t_4 = 460$, and $t_5 = 647$. The average actual reproduction number for each time period is shown in Table 4.

---

4.3. Effective reproduction number, $\mathcal{R}_E(t)$

The effective reproduction number, $\mathcal{R}_E(t)$, measures the number of secondary cases caused by an infectious individual in a population that is not completely naïve to the disease. Specifically, $\mathcal{R}_E(t)$ accounts for a population's reduced susceptibility due to previous exposure or vaccination [17]. In general, the effective reproduction number is defined as the product of the basic reproduction number and the proportion of the population that is susceptible [17,13,54]. For the model described by System (1),

$\mathcal{R}_E(t) = \mathcal{R}_0 \left( \frac{S_H(t)+S_L(t)+S_W(t)}{N(t)}\right). \label{eq:RE}$

(3)

It should be noted that $\mathcal{R}_E(t)$ is not constructed to define a threshold describing when an epidemic will occur, as is the case with $\mathcal{R}_0$. Instead, $\mathcal{R}_E(t)$ scales the value of $\mathcal{R}_0$. Furthermore, since the effective reproduction number is defined in terms of $\mathcal{R}_0$ and since $\mathcal{R}_0$ is only valid using parameter values representing the beginning of an outbreak, $\mathcal{R}_E(t)$ would only be a valid reproduction number over the first response period, $t\in[0,95]$. Thus, we do not use the effective reproduction number in our analysis of the effectiveness of the implemented responses to the 2014 outbreak or our analysis of the potential impact of an Ebola vaccine.

---

5. Use of a prophylactic vaccine

The magnitude of the 2014 EVD outbreak in West Africa hastened the development of a vaccine against Ebola. In July 2015, it was reported that the rVSV-ZEBOV vaccine showed 100% (95% CI 74.7-100.0; $p =0.0036$) efficacy in those vaccinated immediately after confirmation of an infected contact [32,31]. The vaccine was less effective for individuals for whom vaccination was delayed to 21 days after confirmation of an infected contact [32,31]. To examine the potential impact of the rVSV-ZEBOV vaccine, we (1) evaluated the impact of administering the vaccine during each phase of the outbreak (see Section 5.1), and (2) determined the proportion of the population that would need to be vaccinated prior to an initial infection in order to prevent an epidemic (see Section 5.2).

---

5.1. Distribution of a prophylactic vaccine during the outbreak

Though a prophylactic vaccine is designed to be administered to individuals prior to exposure in order to prevent an outbreak, if a prophylactic vaccine against ebolavirus had existed during the beginning of the 2014 EVD outbreak in West Africa, it likely would have been administered to individuals not yet exposed to the virus (the susceptible population). To determine what impact this would have had on the outbreak, using System (1), we calculated the cumulative number of infections and deaths that would have occurred in Sierra Leone if a vaccination campaign had started at $t=95$ with the beginning of the WHO Phase I response, at $t=202$ with the beginning of the UNMEER intervention strategy, at $t=248$ with the beginning of the WHO Phase 2 response, or at $t=460$ with the beginning of the WHO Phase 3 response. These calculations were compared with the cumulative infections and deaths predicted by System (1) in the absence of a vaccination campaign, thus yielding the number of infections and deaths that would have been prevented by a vaccination campaign.

We consider four different vaccination strategies. In each case, we assume that the vaccine has 100% efficacy ($\epsilon=1$) [32,31]. The first two strategies vaccinate 90% of healthcare workers ($S_W$) over one year ($\hat{\rho} ={\raise0.5ex\hbox{$\scriptstyle 0.9$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 365$}}$) and either 30% of the general susceptible population ($S_H+S_L$) over one year ($\rho_H = \rho_L = {\raise0.5ex\hbox{$\scriptstyle 0.3$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 365$}}$, a moderate vaccination campaign) or 90% of the general susceptible population ($\rho_H = \rho_L = {\raise0.5ex\hbox{$\scriptstyle 0.9$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 365$}}$, an ambitious vaccination campaign). We also consider two ring vaccination strategies. Ring vaccination is an ambitious vaccination strategy which requires a list of contacts, as well as contacts of contacts, of an infectious individual to be traced and vaccinated as quickly as possible after their possible exposure to the virus [31]. In each ring vaccination strategy we assume 99% of healthcare workers ($S_W$) are vaccinated each week ($\hat{\rho} ={\raise0.5ex\hbox{$\scriptstyle 0.99$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 7$}}$), 95% of the highly susceptible individuals are vaccinated every three weeks ($\rho_H = {\raise0.5ex\hbox{$\scriptstyle 0.95$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 21$}}$), and either 1% of low susceptible individuals are vaccinated over one year ($\rho_L = {\raise0.5ex\hbox{$\scriptstyle 0.01$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 365$}}$) or 30% of low susceptible individuals are vaccinated over one year ($\rho_L = {\raise0.5ex\hbox{$\scriptstyle 0.3$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 365$}}$).

Table 5 shows the cumulative number of infections and deaths prevented if the vaccination campaign had been introduced at the beginning of each of the intervention periods. If the moderate vaccination campaign (the first strategy discussed above and in Table 5) had begun with the first intervention response (WHO Phase 1 starting at $t=95$), then our model predicts 5,573 infections and 1,678 deaths would have been prevented in Sierra Leone. Though this is a significant improvement over the baseline of the no-vaccination campaign, it still fails to prevent over 60% of the predicted cumulative infections and over 60% of the predicted cumulative deaths. Furthermore, the later the beginning of the vaccination campaign occurs (as $t$ increases), the fewer infections and deaths there are prevented. Given that the primary efficacy results of the rVSV-ZEBOV vaccine were published in July 2015, if the vaccination campaign had started shortly thereafter with the implementation of the WHO Phase 3 response plan ($t=460$; July 31,2015), only 111 infections and 24 deaths would have been prevented. The case of the more ambitious vaccination campaign (the second strategy discussed above and in Table 5) results in the prevention of greater numbers of infections and deaths. However, even if the more ambitious vaccination campaign had started with the beginning of the WHO Phase 1 response ($t=95$), there still would have been over 5,100 infections and over 1,600 deaths, which is considerably larger than any previous EVD outbreak (the largest previous outbreak having occurred in Uganda in 2000-2001 with 425 infections and 224 deaths [12]). The two ring vaccination campaigns show very little to no improvement over the moderate vaccination campaign. This means that it would be a waste of time and resources to utilize a ring vaccination campaign during an outbreak with similar transmission rates to those that were present in the 2014 West Africa outbreak.

---

5.2. Calculating the vaccination threshold

A prophylactic vaccine is designed to be administered prior to exposure in order to prevent an outbreak. Thus we determined the proportion of the population that would have needed to have been vaccinated prior to the 2014 EVD outbreak in Sierra Leone in order to prevent that outbreak from occurring; we refer to this proportion as the vaccination threshold.

The vaccination threshold is determined by simulating System (1) assuming no intervention and with the same initial conditions as described in Section 3 (but with only one initially infected individual). Additionally, at time $t=0$ a proportion $x$ of the general population ($S_H+S_L$) and a proportion $y$ of the healthcare workers ($S_W$) have been vaccinated and are thus moved to the $R$ and $V_W$ classes, respectively. Specifically, $S_H(0) = 5(1-x)I(0)$, $I(0) = 1$, $S_W(0) = 0.00039(1-y)N(0)$, $V_W(0) = 0.00039yN(0)$, $R(0) = xN(0)$, $S_L(0) = N(0) - S_H(0) - I(0) - S_W(0) - V_W(0) - R(0)$, and $E(0) = \hat{I}(0) = R_I(0) = D(0) = \hat{D}(0) = B(0) = 0$, where $N(0)$ is the population size of Sierra Leone. We refer to $x$ as the vaccination coverage.

Given that healthcare workers are at a greater risk of exposure during an outbreak, we assume that any efforts to vaccinate a population against ebolavirus prior to an outbreak would preferentially vaccinate healthcare workers. Thus, if $x>0$, then we assume 95% of the healthcare workers will have been vaccinated as well (i.e., $y=0.95$). The day of the outbreak's peak is defined as $t^*$ such that $I(t^*)+\hat{I}(t^*) = \displaystyle\max_{t}\left[I(t)+\hat{I}(t)\right]$. We determine the minimum $x$ for which $t^*=0$, i.e., the outbreak never ''takes off''.

Table 6 shows the effect of vaccination prior to an initial outbreak given various levels of vaccination coverage ($x$). For each proportion $x$, System (1) is simulated using parameters representing four vaccination scenarios. All four scenarios assume no intervention other than vaccination, and thus the model is simulated with parameters from Table 1 and the first column of Table 3. The first scenario assumes vaccination only occurs prior to $t=0$ (i.e., $\rho_L=\rho_H=\hat{\rho} = 0$); the second scenario assumes there is a continuing vaccination campaign using the parameters of the moderate vaccination campaign from Section 5.1 (i.e., $\rho_L = \rho_H = {\raise0.5ex\hbox{$\scriptstyle 0.3$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 365$}}$, $\hat{\rho}={\raise0.5ex\hbox{$\scriptstyle 0.9$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 365$}}$); the third scenario assumes a ring vaccination strategy with ambitious vaccination levels for healthcare workers ($S_W$) and the highly susceptible population ($S_H$) (i.e., $\rho_L={\raise0.5ex\hbox{$\scriptstyle 0.01$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 365$}}$, $\rho_H={\raise0.5ex\hbox{$\scriptstyle 0.95$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 21$}}$, $\hat{\rho}={\raise0.5ex\hbox{$\scriptstyle 0.99$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 7$}}$), and lastly the fourth scenario assumes an even more ambitious ring vaccination strategy with ambitious levels of vaccination for healthcare workers ($S_W$), the highly susceptible population ($S_H$), and also the low susceptible population ($S_L$) (i.e., $\rho_L={\raise0.5ex\hbox{$\scriptstyle 0.3$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 365$}}$, $\rho_H={\raise0.5ex\hbox{$\scriptstyle 0.95$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 21$}}$, $\hat{\rho}={\raise0.5ex\hbox{$\scriptstyle 0.99$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 7$}}$). Each simulation is run until $I(t) + \hat{I}(t) < 1$; i.e., no new cases. In addition to calculating the day of the outbreak's peak ($t^*$), the cumulative number of cases is calculated for each simulation.

In the first scenario, where vaccination only occurs prior to $t=0$, if as little as 36% of the general susceptible population is vaccinated prior to an initial infection (1 initial case), then the outbreak is unable to spread through the population (i.e., $t^*=0$). If there is a continuing vaccination campaign beyond $t=0$, then as little as 30% of the general susceptible population needs to be vaccinated prior to an initial infection. With both of the ring vaccination strategies, as little as 35% of the general susceptible population needs to be vaccinated prior to an initial infection.

In all four scenarios in Table 6, the proportion of the individuals in the $E$ class entering the hospital is 40.7% ($\alpha = 0.407$) and the effectiveness of hospital workers in preventing deaths due to EVD is 30.9% ($1-\nu = 0.309$). Since the parameters $\alpha$ and $\nu$ both represent critical elements in controlling and EVD outbreak, using System (1) we calculated the cumulative number of cases for varying values of $\alpha$ and $\nu$, along with varying values in vaccination coverage ($x$) and the proportion of the low susceptible population which is vaccinated per year ($\rho_L\times 365$). All other parameter values were taken from Table 1 and the first column of Table 3 (the same used in generating Table 6) except that we assumed $\rho_H={\raise0.5ex\hbox{$\scriptstyle 0.95$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 21$}}$ and $\hat{\rho} = {\raise0.5ex\hbox{$\scriptstyle 0.99$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 7$}}$. The matrix plots in Figure 6 show the relationships between $\alpha$, $\nu$, $x$, and $\rho_L$. Taken together, the graphs in Figure 6 show that the cumulative number of infections (CI) is most sensitive to changes in $\alpha$ and $x$, and least sensitive to changes in $\nu$.

---

6. Discussion

We have proposed a model that simulates the temporal dynamics of an EVD outbreak. We used this model, along with cumulative infections and death data from the CDC, to determine subpopulation-specific transmission rates ($\beta$'s), the proportion of symptomatic individuals who sought medical treatment ($\alpha$), and the effectiveness of healthcare workers to reduce deaths due to EVD ($\nu$) for each response period over the course of the outbreak. The parameter fitting showed that the proportion of symptomatic individuals who sought medical treatment was substantially higher in the time periods where the WHO and UN were providing intervention measures as compared to the time period where there was no response. One of the means by which the WHO and UN planned to slow transmission was to get more symptomatic individuals to hospitals and clinics. Our model indicates that the WHO and UN were successful in their efforts. Additionally, the parameter fitting revealed that healthcare workers were most effective in reducing deaths due to EVD during the WHO Phase 1 response period. Lastly, the parameter fitting also revealed that the transmission rates ($\beta$ parameters) did not decrease for all subpopulations as the WHO and UN response strategies were implemented, demonstrating that the effects of controls strategies on the transmission to specific subpopulations can be non-intuitive and are thus aided by the examination of mathematical models.

We also used our model to show that, even if there had been an ambitious rVSV-ZEBOV vaccination campaign during the 2014 EVD outbreak, it would not have been enough to halt the epidemic. This is due to the fact that rVSV-ZEBOV is a prophylactic vaccine designed to prevent infection and does not treat those already infected. To determine the effectiveness of the rVSV-ZEBOV vaccine in preventing an outbreak we calculated the vaccination threshold. Based on our findings we recommend that the Sierra Leone government aim to vaccinate at least 95% of their healthcare workers and 40% of the general population. Compared to vaccination thresholds needed to achieve herd immunity [30] in other diseases ($\geq70\%$ for diseases like polio [42], smallpox [4,24], and measles [28]), this is a relatively low vaccination threshold. Of course, higher vaccination coverage with rVSV-ZEBOV of general population will provide greater protection. By the end of the 2014 EVD outbreak in Sierra Leone roughly 17,000 individuals had survived the disease and would thus have immunity. Though this is a sizable number, it is still less than 1% of the total population of Sierra Leone and is thus not sufficient for preventing another outbreak. Additionally, if another outbreak does occur, our findings show that efforts spent on increasing the proportion of infected who seek treatment from healthcare professionals will provide the greatest reduction in cumulative infections if it can be assume that individuals, once in a hospital or healthcare facility, are no longer able to transmit the disease to the general population, thus acting as a type of quarantine.

Since the waning of the 2014 EVD outbreak, there have been flare-ups. For example, in late February and early March 2016, there were eight EVD cases in Guinea after the outbreak had been declared over on December 29,2015 [61]. These flare-ups are due to the virus remaining present in individuals who are recovering and show the importance of the inclusion of the $R_I$ class in the model. Though our model cannot predict flare-ups due to its lack of stochasticity, the inclusion of the $R_I$ class does allow for an outbreak to persist longer in a population before self-extinguishing. In response to the March 2016 flare-up in Guinea, hundreds of individuals who were in contact with the eight infected were identified and vaccinated using a ring-vaccination strategy [61]. Though the model we propose is not a network model and thus cannot explicitly simulate ring vaccination, the division of the general susceptible population into high and low risk of exposure provides a means for examining the effects of a ring vaccination strategy through differential vaccination rates of the high and low risk susceptible individuals. Our results show that given an outbreak, a greater reduction in cumulative infections will be achieved through a moderate vaccination campaign that does not distinguish between high and low susceptible individuals than with a ring vaccination campaign. This is due to the fact that the $R_I$ class remains infectious with a lower transmission probability for a prolonged period and can infect individuals in the low susceptible class.

In light of the 2014 EVD outbreak in West Africa, ebolavirus has become a disease of much greater concern than it historically has been. To aid in understanding the dynamics of the disease, many mathematical models have been proposed. The novel model we propose adds another tool for understanding EVD dynamics, and our analysis of this novel model shows the impact of the implemented response strategies on transmission rates and other parameters in Sierra Leone. Lastly, and most importantly, our analysis provides a vaccination threshold for the new rVSV-ZEBOV vaccine in Sierra Leone. Furthermore, the ability of our model in predicting the outcome of an array of vaccination coverages makes our model a useful tool for determining future vaccination policy.

---