A new epidemics–logistics model: Insights into controlling the Ebola virus disease in West Africa

European Journal of Operational Research 265 (2018) 1046–1063

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Innovative Applications of O.R.

A new epidemics–logistics model: Insights into controlling the Ebola virus disease in West Africa I.˙ Esra Büyüktahtakın a,∗, Emmanuel des-Bordes b, Eyyüb Y. Kıbıs¸ c a b c

Department of Mechanical and Industrial Engineering, New Jersey Institute of Technology, Newark, NJ, USA Department of Industrial and Manufacturing Engineering, Wichita State University, Wichita, KS, USA Huether School of Business, The College of Saint Rose, Albany, NY, USA

a r t i c l e

i n f o

Article history: Received 4 January 2017 Accepted 20 August 2017 Available online 24 August 2017 Keywords: (S) Decision support systems Infectious disease Spatially explicit optimization Epidemic control Ebola virus disease

a b s t r a c t Compartmental models have been a phenomenon of studying epidemics. However, existing compartmental models do not explicitly consider the spatial spread of an epidemic and logistics issues simultaneously. In this study, we address this limitation by introducing a new epidemics–logistics mixed-integer programming (MIP) model that determines the optimal amount, timing and location of resources that are allocated for controlling an infectious disease outbreak while accounting for its spatial spread dynamics. The objective of this proposed model is to minimize the total number of infections and fatalities under a limited budget over a multi-period planning horizon. The present study is the first spatially explicit optimization approach that considers geographically varying rates for disease transmission, migration of infected individuals over different regions, and varying treatment rates due to the limited capacity of treatment centers. We illustrate the performance of the MIP model using the case of the 2014–2015 Ebola outbreak in Guinea, Liberia, and Sierra Leone. Our results provide explicit information on intervention timing and intensity for each specific region of these most affected countries. Our model predictions closely fit the real outbreak data and suggest that large upfront investments in treatment and isolation result in the most efficient use of resources to minimize infections. The proposed modeling framework can be adopted to study other infectious diseases and provide tangible policy recommendations for controlling an infectious disease outbreak over large spatial and temporal scales. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Infectious disease outbreaks have been an immense challenge for humanity. The recent example of an infectious disease outbreak is the 2014–2015 Ebola outbreak (World Health Organization, 2014c). Despite recent signs of slowing down, the deadly Ebola virus disease (EVD) has spread rapidly across West Africa including the three most-affected countries–Guinea, Liberia, and Sierra Leone–since early March 2014. As of June 12, 2016, the Centers for Disease Control and Prevention (CDC) reported that the total number of Ebola cases and deaths in West Africa has reached 28,652 and 11,325, respectively (Centers for Disease Control & Prevention (CDC), 2015). The Ebola virus is considered one of the deadliest viruses ever known, with a case fatality proportion varying from 25 to 90% in past outbreaks (World Health Organization, 2014c). Ebola symptoms include sudden onset of fever, weakness, vomiting, diarrhea,

∗

Corresponding author. ˙ E-mail addresses: [email protected], [email protected] (I.E. Büyüktahtakın).

http://dx.doi.org/10.1016/j.ejor.2017.08.037 0377-2217/© 2017 Elsevier B.V. All rights reserved.

headache, and sore throat. Severe and fatal stages of the disease are accompanied by internal and external bleeding and multiple organ failures (Chowell & Nishiura, 2014; World Health Organization, 2014c). EVD is transmitted among humans via direct contact with bodily fluids from an infected person or indirect contact with contaminated surfaces (World Health Organization, 2014c). EVD transmission is further aggravated by traditional West African funeral practices such as washing, touching, and kissing the body and unsafe burial of Ebola-infected bodies (WHO Ebola Response Team, 2014; World Health Organization, 2014c). Ebola has no cure, vaccine or specific treatment for infected individuals. Therefore, shortterm intervention strategies to control transmission include medical treatment in the form of supportive therapy, such as maintenance of fluids (World Health Organization, 2014c), quarantine, isolation of cases, contact tracing, and safe burial practices. Safe burial consists of sanitizing and placing the dead body in a disinfected body bag, and burying it in a grave at least 2 meters deep (Nielsen et al., 2015). On the other hand, Ebola treatment centers (ETCs) play a significant role in isolating and treating infected individuals, thereby helping to curb the Ebola outbreak (World Health Organization, 2014b).

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

Controlling epidemics such as Ebola over a spatial scale and in a prompt manner is a challenging optimization problem. First, infectious diseases spread through the air or other mechanisms that need to be accounted for when modeling an epidemic (Cohen, 20 0 0). Second, capturing spatial considerations is essential because critical disease parameters such as transmission and fatality rates among the population may change in various geographical regions. In addition, infectious diseases such as Ebola can be easily transmitted from one region (e.g., city, country, continent) to another region (other cities, countries, and continents) due to the mobility and migration of the population (Wang & Zhao, 2004). Third, the epidemic intervention resources are limited and insufficient in most circumstances, forcing policymakers and healthcare providers to determine the most effective policies to allocate resources strategically and rapidly to the right locations for controlling an infectious disease outbreak. Compartmental models in epidemiology decompose the members of a population into different subgroups (compartments) according to the stages of the infectious disease and characterize the rates of transiting from one subgroup into another, duration of the infection, and contact rates (see, e.g., Bailey, 1975; Frauenthal, 2012). For example, the susceptible-infected-recovered (SIR) model and its variations have been widely used to predict the growth of epidemics and to inform intervention strategies (Anderson & May, 1991; Kermack & McKendrick, 1932; Tebbens & Thompson, 2009). Epidemic models have been developed to analyze a number of infectious diseases, including Ebola (Althaus, Gsteiger, & Low, 2014; Chowell, Hengartner, Castillo-Chavez, Fenimore, & Hyman, 2004; Chowell & Nishiura, 2014; Legrand, Grais, Boelle, Valleron, & Flahault, 2007; Meltzer et al., 2014; Pandey et al., 2014; Zaman, Kang, & Jung, 2009), tuberculosis (Long, Vaidya, & Brandeau, 2008), measles (Babad et al., 1995), influenza (Rachaniotis, Dasaklis, & Pappis, 2012), plague (Allen, Brauer, Van den Driessche, & Wu, 2008), human immunodeficiency virus (HIV) (Lasry, Zaric, & Carter, 2007), smallpox (Ferguson et al., 2003), and anthrax (Craft, Wein, & Wilkins, 2005). Some studies focus on discrete-time epidemic models (Allen, 1994) in order to alleviate the difficulties pertaining to differential equations. While the spatial characteristic is a central feature of many epidemic control problems, where infections show migration or dispersal patterns, most previous work has oversimplified the problem by ignoring the spread of disease over space. To our knowledge, none of the previous studies on Ebola modeling has incorporated the logistics of treatment using a spatially explicit optimization model in order to allocate limited resources to control the spread of the disease. Most studies on Ebola have considered continuous-time differential equations and stochastic simulation models, which are difficult to include in an optimization model. Furthermore, the current epidemic modeling literature does not consider treatment capacity and thus assumes a constant treatment rate, which leads to suboptimal solutions for capacity allocation over time and space. The objective of this paper is to develop a general spatiotemporal optimization modeling framework to help decisionmakers minimize infections and death due to an epidemic. The model provides information on the spread dynamics of infections, and where and when to allocate limited resources. The strength of the model is demonstrated by a case of determining effective strategies to suppress the 2014–2015 Ebola outbreak. To the best of our knowledge, there have not been any known spatially explicit optimization studies that integrates both epidemic growth and optimal resource allocation modeling in one optimization model for controlling an infectious disease outbreak over a finite planning horizon. This paper is organized as follows. Section 2 reviews the relevant literature on epidemic disease modeling and logistics for

1047

controlling disease outbreaks. In Section 3, we present the epidemics–logistics optimization model including the susceptibleinfected-treated (and quarantined)-recovered-funeral-buried (SITRFB) compartment model. Section 4 presents a real case regarding the 2014–2015 Ebola outbreak in Guinea, Liberia, and Sierra Leone and provides related economic, geographic, population, and disease transmission data. More detailed information on the data and data sources is provided in Appendix A. In Section 5, we present the validation of the proposed model, detailed computational experiments, and numerical results to explain the Ebola spatial spread and disease dynamics in each of the three most-affected nations. Here we also provide insights into optimal intensity, location, and timing of intervention for controlling the Ebola epidemic. Finally, in Section 6, we summarize our findings and present policy implications and recommendations. Section 7 presents conclusions and directions for future research. 2. Literature review and paper contributions Operations Research (OR) approaches have been implemented to incorporate treatment capacity and evaluate the performance of different control measures for Ebola and other epidemic diseases. Previous OR or management science (MS) approaches on epidemics control involve compartmental model-based simulations (Legrand et al., 2007; Meltzer et al., 2014; Pandey et al., 2014), differential equations (Craft et al., 2005; Kaplan, Craft, & Wein, 2003), costeffectiveness analysis of resource allocation (Tebbens & Thompson, 2009; Zaric, Brandeau, & Barnett, 2000), network models (Ancel, Newman, Martin, & Schrag, 2003; Eubank et al., 2004), simulation optimization (Lee, Pietz, Benecke, Mason, & Burel, 2013; Longini, Halloran, Nizam, & Yang, 2004; Patel, Longini, & Halloran, 2005), and mathematical programming (Ren, Ordóñez, & Wu, 2013; Tanner, Sattenspiel, & Ntaimo, 2008; Yarmand, Ivy, Denton, & Lloyd, 2014). Various type of facility location models, including maximal covering, location-allocation, and P-median models, have also been studied in response to large-scale emergencies and disaster management (Jia, Ordóñez, & Dessouky, 2007a,b). A number of surveys regarding OR/MS contribution to epidemics and disaster control can be found in (Altay & Green, 2006; Dasaklis, Pappis, & Rachaniotis, 2012; Dimitrov & Meyers, 2010). The majority of mathematical models that study the logistics of controlling epidemics use simulation methods (e.g., Arinaminpathy & McLean, 2009; Berman & Gavious, 2007; Dasaklis, Rachaniotis, & Pappis, 2017; Liu & Zhao, 2012; Longini et al., 2007; Porco et al., 2004; Riley & Ferguson, 2006). Some of these studies model simulations using network relations among disease compartments to study the spatio-temporal aspects of the problem (e.g., Berman & Gavious, 2007; Longini et al., 2007; Porco et al., 2004; Riley & Ferguson, 2006). Those simulation models usually evaluate various policies (e.g., vaccination of contacts, limited response capacity, heterogeneity in symptoms and infectiousness, more rapid diagnosis due to public awareness, and isolation of cases) on intervention to control the disease spread (Porco et al., 2004). A number of studies have used agent-based simulation models and social networks to explicitly capture interactions and contacts among sub-populations, such as susceptible, infectious, recovered, or immuned of the Ebola virus disease (Ajelli et al., 2016; Kurahashi & Terano, 2015; Siettos, Anastassopoulou, Russo, Grigoras, & Mylonakis, 2011; 2016; Wells et al., 2015). These simulations are used to assess the impact of various intervention strategies, such as ring vaccination (Kurahashi & Terano, 2015), case isolation (Wells et al., 2015), and burial practices (Siettos, Anastassopoulou, Russo, Grigoras, & Mylonakis, 2016). Some of the dynamic models are driven by the spatial correlation of individuals in the population (e.g., Wells et al., 2015). They study spatial heterogeneity in

1048

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

the epidemic but do not provide explicit information on intervention timing and intensity for each specific region of the most affected countries (Ajelli et al., 2016). Sophisticated compartmental model-based simulations can incorporate important features of an epidemic to capture the realism of the disease dynamics. Yet, they are limited by only evaluating a relatively small set of control strategies rather than considering the full spectrum of control policy options (Dimitrov & Meyers, 2010). Simulation optimization also requires a large amount of computer time per simulation run and uses heuristics, which do not guarantee optimality (Tanner et al., 2008). Recently, researchers study mathematical optimization such as linear and integer programming, non-linear optimization, and stochastic programming to analyze all possible alternatives and determine the optimal one for controlling an epidemic (Ren et al., 2013; Tanner et al., 2008; Yarmand et al., 2014). Previous mathematical programming models that study the integration of epidemics control and logistics/operations management generally use simulations to generate different scenarios of disease growth as an input into the optimization model. For example, Dasaklis et al. (2017) study the smallpox outbreak problem in two modules, one for simulating the disease’s progression and the other for formulating a linear programming model to optimally distribute a set of medical and ancillary supplies to affected populations. On the other hand, Liu and Zhao (2012) present two models including an SEIR epidemic and emergency distribution model to study a bioterrorism problem. Simulating the epidemic model, they generate scenarios to analyze the logistics of the emergency distribution problem. Mathematical programming models of epidemic control mostly focus on determining optimal vaccination strategies. Tanner et al. (2008) develop a stochastic programming framework to find the optimal vaccination policy to minimize the total cost of the intervention policy for controlling infectious disease epidemics under parameter uncertainty. Using a chance constraint, the authors target reducing the basic reproductive number to less than one. Yarmand et al. (2014) present a compartmentalbased simulation model of the disease spread to estimate the probability of epidemic containment in each region for different vaccination levels. Then they use vaccination outcome probabilities obtained from the simulation model to formulate a two-stage stochastic program in order to allocate the vaccine doses more efficiently. Ren et al. (2013) develop a Taylor approximation representation of smallpox disease propagation to define disease transmission rates with or without mass vaccination. Then they input those disease transmission rates in a multi-city resource allocation model in which a limited amount of vaccination is distributed with the objective of minimizing the total number of fatalities due to a smallpox outbreak. The authors use a gravity model to determine the fraction of newly infected cases in one city that appear in another city. They employ heuristics to solve the resulting large-scale mixed-integer program. None of these former studies take into account the capacity of hospitals, and they do not address logistics issues such as the number and size of treatment centers to locate in order to intervene in an infectious disease outbreak. Our approach also differs from previous work by explicitly formulating disease transition among compartments, spatial spread dynamics (e.g., movement of infected and susceptible individuals among neighboring regions), and logistics in an optimization model. As opposed to different mathematical programming approaches that use simulation to generate scenarios into the optimization models, we capture the epidemic disease growth and resource allocation problems in one single optimization model. To our knowledge, our paper is also the first mathematical programming study on optimally allocating resources for controlling the Ebola outbreak.

In this paper, we develop a general spatio-temporal optimization modeling framework to help decision-makers minimize infections and death due to an infectious disease outbreak. Our approach contributes to the epidemiology and OR literature, in the following ways: (i) We develop a new epidemics–logistics mixed-integer programming (MIP) model of the epidemic control problem. The model minimizes the total number of people infected and deceased (without treatment) due to the epidemic over a finite planning horizon, subject to availability constraints on resources. Unlike previous models, our model takes into account the spread of the epidemic among multiple regions and the logistics of epidemic control simultaneously in a spatio-temporal setting. A regionalized multi-country landscape setting is considered to represent the spatially heterogeneous growth of the disease, spread, infections, and control costs. This model incorporates different transmission, progression, and recovery rates for different population compartments and each specific location defining a discrete-time SITR-FB model. (ii) This optimization model accounts for fixed and variable costs of treatment resources that are bounded by the budget constraint when determining where ETCs including isolation units should be allocated over space and time for optimal control intervention. The total costs for establishing treatment and isolation units, providing medical treatment, tracing contacts of infected people, and logistics as well as safe burial cannot exceed the available budget allocated for them. (iii) As opposed to previous epidemic compartment models, transmission rates between infected (I) and treated (T) compartments are not constant but instead depend on treatment capacity and the number of infected people receiving medical service. (iv) The migration of individuals plays a major role in the spread of infectious diseases (Soto, 2009). To explain how these movement patterns contribute to the further spread of an epidemic, we incorporate the migration of susceptible and infected individuals between regions of a country in our MIP model. (v) We validate the predictions of the MIP model by demonstrating the impact of actual interventions. We graphically show that the proposed model provides a good approximation of the outbreak data in all three countries over a 60week period. We also use a paired-t-test, which proves that there is no statistically significant difference between the weekly new infections and cases predicted by the proposed model and the original outbreak data. (vi) The computational difficulty of the model is magnified with biological and mathematical complexity. To tackle the complex nature of the model, the non-linear capacity constraints are linearized leading to an MIP model. The MIP model is preferable than non-linear and continuous epidemic models because the implementation of discrete linear models is more straightforward for non-technical users, thus providing an advantage in the public health world (Brauer, Feng, & Castillo-Chavez, 2010). (vii) Although the MIP model is general and can be adopted to study other infectious disease outbreaks, we illustrate the results of the model on the real case of the 2014–2015 Ebola outbreak in Guinea, Liberia, and Sierra Leone, where both epidemic features and logistic issues are addressed. Our results provide insight into the future growth and spread of the epidemic under various resource allocation scenarios as well as the impact of critical problem parameters such as disease transmission rates and intervention budget.

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

1049

Fig. 1. Network representation of transmission dynamics of Ebola and SITR-FB model for each location (k, l). Note that, for clarity, only flow rates are shown on the arcs of the figure.

Furthermore, our in-depth spatial and temporal analysis provides decision-makers with policies and strategies for the optimal distribution of intervention budget among various operations, including treatment, isolation, contact tracing, and safe burial, and the allocation of resources, such as ETCs, to the right locations and at the right time. 3. Epidemics–logistics model This section presents the proposed epidemics–logistics MIP optimization model that describes the transmission dynamics of the Ebola virus. We note that although the compartmental model is structured to capture the disease dynamics of Ebola, it can be modified for other types of infectious disease outbreaks using disease-specific compartments and rates. The compartment model is described for a geographical area that is represented by a twodimensional grid with K¯ rows and L¯ columns. For example, each location (k, l) represents a region of the three most-affected West African nations, which is defined by row k and column l. The period index j ∈ J captures the dynamic nature of the disease transmission as well as resource allocation decisions in the model. In the application of the model, each period j corresponds to one week of a 60-week treatment planning time horizon. 3.1. Epidemic compartmental model The biological part of the MIP model (SITR-FB) describes the Ebola disease transmission among six mutually exclusive compartments in each spatial location (k, l), as illustrated in Fig. 1. In the SITR-FB model, compartments (nodes) in the network correspond to the health status of individuals with respect to the EVD in a population, while edges describe the transition between given compartments. As the disease progresses through a population, susceptible individuals (S) are infected and thus move to the infectious (I) compartment via either person-to-person contact in the community at a periodic rate of σ1(k,l ) or touching Ebola-related dead bodies that are not yet buried during traditional funerals at a periodic rate of σ2(k,l ) . At compartment I, a fraction μ1(k,l ) of the infected individuals will die and move to the funeral (F) compartment in 1/λ1(k,l ) weeks, while a proportion μ3(k,l ) will recover (R) with immunity in 1/λ2(k,l ) weeks without treatment.

Therefore, both μ1(k,l ) λ1(k,l ) and μ3(k,l ) λ2(k,l ) represent weekly rates of transmission. Additionally, a subset of the infected individuals will be hospitalized for treatment (T) based on the value of the treatment capacity variable C (j k,l ) , i.e., availability of beds in the ETCs, determined by solving the optimization model. Note that

the diagram in Fig. 1 shows no constant transition rate from I to T since all I individuals seeking treatment are hospitalized in an ETC based on the availability of space. We also assume that individuals who are not hospitalized will remain in the community and spread the disease further. Also, a portion (1 − μ2(k,l ) ) of individuals in the treatment compartment will recover at a periodic rate of γ2(k,l ) , while a fraction (μ2(k,l ) ) of them will die at a periodic

rate of γ1(k,l ) . The buried (B) compartment ensures that deceased individuals in the funeral compartment are safely buried at a rate of γ3(k,l ) . Finally, in order to demonstrate the impact of movement between regions of a given country on the disease spread, we consider the migration of susceptible and infected individuals into and out of S and I compartments with rates of (α (h, q) , ν (k, l) ) and (φ (h, q) , ρ (k, l) ), respectively, for location (k, l). Each parameter is defined in detail in the notation section, and corresponding case study values are presented in Table 3.

3.2. Model notations Model notations that will be used throughout the rest of this paper follow: Sets and indices: J N K L j (k, l ) M(k,l )

Set of time periods, J = {0, . . . , J¯} Set of ETC types, N = {1, . . . , N¯ } Set of rows, K = {1, . . . , K¯ } Set of columns, L = {1, . . . , L¯ } Index for time period where j ∈ J Index for location where k ∈ K and l ∈ L Set of all surrounding regions of location (k, l )

Transition parameters used to describe the rate of movement between compartments: μ1(k,l ) μ2(k,l ) μ3(k,l ) σ1(k,l ) σ2(k,l ) 1

λ1(k,l ) 1

λ2(k,l ) 1

γ1(k,l ) 1

γ2(k,l ) (k,l )

γ3

Disease fatality proportion without treatment in location (k, l ) Disease fatality proportion while receiving treatment in location (k, l ) Disease survival proportion without treatment in location (k, l ) Transmission rate per person due to community interaction in location (k, l ) Transmission rate per person during traditional funeral ceremony in location (k, l ) Average time (in periods) from infection to death without treatment in location (k, l ) Average time (in periods) from infection to recovery without treatment in location (k, l ) Average time (in periods) from treatment to death in location (k, l ) Average time (in periods) from treatment to recovery in location (k, l ) Safe burial rate of Ebola-related dead bodies in location (k, l )

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

1050

Other parameters: (k,l )

b1 j

(k,l )

b2 j

) g(nk,l j

κn

1

2 π (k,l ) (k,l ) θ (k,l ) ϑ (k,l ) υ (k,l ) τ (k,l ) ς (k,l ) α (h,q)→(k,l ) φ (h,q)→(k,l ) ν (k,l )→(h,q) ρ (k,l )→(h,q)

Unit cost of treatment for an infected individual in location (k, l ) at end of period j Unit cost of safe burial for a dead body in location (k, l ) at end of period j Fixed cost of establishing type n ETC in location (k, l ) at end of period j Capacity (number of beds) of type n ETC Total available budget for treatment Total available budget for safe burial Initial number of susceptible individuals in location (k, l ) Initial number of infected individuals in location (k, l ) Initial number of treated individuals in location (k, l ) Initial number of recovered individuals in location (k, l ) Initial number of unburied dead bodies (funerals) in location (k, l ) Initial number of buried dead bodies in location (k, l ) Initial treatment capacity in location (k, l ) Migration rate of susceptible individuals from surrounding locations (h, q ) to location (k, l ) Migration rate of infected individuals from surrounding locations (h, q ) to location (k, l ) Migration rate of susceptible individuals from location (k, l ) to surrounding locations (h, q ) Migration rate of infected individuals from location (k, l ) to surrounding locations (h, q )

State variables: S(jk,l )

I (j k,l )

T j(k,l ) R(jk,l )

Fj(k,l )

B(jk,l ) Sˆ(k,l ) j

S˜(jk,l ) Iˆ(j k,l ) I˜(j k,l )

Number of susceptible individuals in location (k, l ) at end of period j Number of infected individuals in location (k, l ) at end of period j Number of individuals receiving treatment in location (k, l ) at end of period j Number of recovered individuals in location (k, l ) at end of period j Number of deceased individuals due to the epidemic in location (k, l ) at end of period j Number of buried individuals in location (k, l ) at end of period j Number of of period j Number of of period j Number of period j Number of period j

susceptible individuals migrating into location (k, l ) at end susceptible individuals migrating from location (k, l ) at end infected individuals migrating into location (k, l ) at end of infected individuals migrating from location (k, l ) at end of

Decision variables: C (j k,l ) I¯(j k,l ) ) yn(k,l j

Total capacity (number of beds) of established ETCs in location (k, l ) at end of period j Number of infected individuals hospitalized (and quarantined) in location (k, l ) at end of period j Number of type n ETCs established in location (k, l ) at end of period j

3.3. Model assumptions The epidemic compartmental model is constructed based on Ebola characteristics. In our SITR-FB model, we assume that individuals in the recovered compartment receive lifelong immunity (Griggs, 2014). In addition, we do not account for natural births and deaths. This is a valid assumption because the considered planning horizon for the epidemic is relatively short. We consider the migration of susceptible and infected individuals among regions of a given country for treatment or business. We do not consider the migration of individuals under treatment because treated people are assumed to be quarantined in an isolation unit at the same time. Here, we neglected the emigration of funerals and recovered individuals. However, emigration of other compartments could be easily incorporated into the model using a similar approach to ours. The infection is transmitted primarily either by person-to-person contact through interaction, or during traditional funerals and burial ceremony via washing and touching Ebolarelated dead bodies. We note that the number of new infections in our model is a function of the transmission rate and the number of

infected people and funerals. However, the transmission rate is assumed not to be a function of susceptible people because the number of susceptible populations will not decline enough to make a noticeable difference in the infection rate for Ebola or other similar infectious disease outbreak on a large spatial scale. This assumption also alleviates the unnecessary non-linearity due to defining the transmission rate as a function of susceptible people. Furthermore, because the latency (exposed) period of Ebola can be as short as two days (World Health Organization, 2014a), we assume that once infected, susceptible individuals become infectious immediately and move into the infected compartment, as in the work of Zaman et al. (2009). This could eliminate unnecessary complications, while an exposed compartment could be easily incorporated into the model. We assume that healthcare workers are protected from the Ebola virus, and all treated patients are quarantined at the same time in all treatment centers. We also assume that once a treatment facility is established, it stays open until the end of the planning horizon. Once moved to the buried compartment, dead bodies can no longer infect new people. We assume that dead bodies (funerals) are safely buried, and the related safe burial costs is reflected in the budget. In our analysis, the treatment budget is restricted, while the burial budget is kept ample in order to analyze situations, which could lead to a high number of dead individuals due to tight limitations on treatment resources. In such cases, if both treatment and burial costs are restricted under the same budget, then the high cost of buried people may lead to infeasibility. It is also more likely that different decision-makers decide on the allocation of treatment and burial budget for controlling the epidemic. 3.4. Epidemics–logistics MIP model formulation In this section, we propose an optimization model that handles the epidemiological and logistical considerations of an infectious disease simultaneously. Using the above assumptions, notations, and transitions among Ebola epidemic compartments, as shown in Fig. 1, we formulate the epidemics–logistics model as follows. 3.4.1. Objective function The objective function of the proposed model is to minimize the total number of infected individuals and deaths of infected people who do not receive treatment, in all regions over the finite planning horizon. Therefore, the objective function (1) is formulated as

  Min 1 + μ1(k,l ) λ1(k,l ) I (j k,l )

(1)

j∈J k∈K l∈L

3.4.2. Epidemic (Disease transmission) constraints Here, we describe the disease transmission constraints (2)– (12) to describe the population dynamics presented in the SITR-FB model as demonstrated in Fig. 1. Initial condition constraint We first define the initial conditions for the number of individuals in susceptible, infected, treated, recovered, funeral, and buried compartments, and the total ETC capacity as

S0(k,l ) = π (k,l ) ,

I0(k,l ) = (k,l ) ,

T0(k,l ) = θ (k,l ) ,

R0(k,l ) = ϑ (k,l ) ,

F0(k,l ) = υ (k,l ) ,

B0(k,l ) = τ (k,l ) ,

C0(k,l ) = ς (k,l ) ,

k ∈ K, l ∈ L.

(2)

Migration constraint The immigration of susceptible and infected individuals into location (k,l) from neighboring locations (h, q) ∈ M(k, l) with the corresponding rates of immigration is defined in Eqs. (3) and (4), respectively, as follows:

Sˆ(jk,l ) =



(h,q )∈M(k,l )

α (h,q)→(k,l ) S(jh,q) ,

j ∈ J, k ∈ K, l ∈ L,

(3)

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

Iˆ(j k,l ) =



φ (h,q)→(k,l ) I (j h,q) ,

(h,q )∈M(k,l )

j ∈ J, k ∈ K, l ∈ L.

(4)

The out-migration of susceptible and infected individuals from location (k, l) to neighboring locations (h, q) ∈ M(k, l) with the corresponding rates of emigration is defined in Eqs. (5) and (6), respectively, as follows:

S˜(jk,l ) = I˜(j k,l ) =



(h,q )∈M(k,l )

 (h,q )∈M(k,l )

ν (k,l )→(h,q) S(jk,l ) ,

j ∈ J, k ∈ K, l ∈ L,

(5)

ρ (k,l )→(h,q) I (j k,l ) ,

j ∈ J, k ∈ K, l ∈ L.

(6)

) S((k,l j+1 )

(k,l )

= Sj

ˆ(k,l )

+ Sj

˜(k,l )

− Sj

(k,l ) (k,l )

− σ1

Ij

(k,l ) (k,l )

− σ2

Fj

j ∈ J\J¯, k ∈ K, l ∈ L.

(7)

Infected individuals. According to Constraint (8), the number of infected individuals at the end of period j + 1 is equal to the number of infected individuals in the previous period, plus immigration of infected individuals into location (k, l), minus the emigration of infected from location (k, l), plus newly infected individuals due to contacting infected and funerals, and minus infected individuals, who recovered, died, or accepted for treatment at the end of period j. Constraint (8) defining the number of infected people for all periods j ∈ J\J¯ and over all locations k ∈ K, l ∈ L is given as

 ) I((k,l = I (j k,l ) + Iˆ(j k,l ) − I˜(j k,l ) + σ1(k,l ) I (j k,l ) + σ2(k,l ) Fj(k,l ) − μ1(k,l ) λ1(k,l ) j+1 ) (k,l )

+ μ3

 (k,l ) (k,l )

λ2

Ij

¯(k,l )

− Ij

(8)

Treated individuals. Constraint (9) represents the total number of individuals undergoing treatment in location (k, l) at the end of period j + 1, which is equal to the treated people in period j plus infected individuals accepted for treatment based on the availability of beds (I¯) minus treated individuals who die, or recover while receiving treatment. k,l ) (k,l ) T((j+1 + I¯(j k,l ) − ) = Tj



   μ2(k,l ) γ1(k,l ) + 1 − μ2(k,l ) γ2(k,l ) Tj(k,l )

j ∈ J\J¯, k ∈ K, l ∈ L.

(9)

Recovered individuals. Constraint (10) provides the cumulative number of individuals in location (k, l) who recover with or without treatment. According to Constraint (10), the total recovered individuals at the end of period j + 1 is equal to the recovered individuals in the previous period plus the newly recovered individuals coming from compartments T or I.



) R((k,l = R(jk,l ) + 1 − μ2(k,l ) j+1 )

 (k,l ) (k,l ) γ2 Tj + μ3(k,l ) λ2(k,l ) I (j k,l )

j ∈ J\J¯, k ∈ K, l ∈ L.

(10)

Dead individuals. Constraint (11) represents the total number of unburied funerals at the end of period j + 1 resulting from infected and treated individuals minus the buried dead bodies. k,l ) F((j+1 )

(k,l )

= Fj

3.4.3. Logistics and operations management constraints Here, we describe constraints (13)–(15) regarding logistics and operations management. Constraints (13)–(15) are simultaneously solved with constraints defining the disease dynamics (2)–(12) in the optimization model. Budget constraint. Constraint (13) imposes a budget limitation on the sum of the fixed costs for locating different-sized treatment facilities (e.g., ETCs) and the variable cost of treating infected individuals over all locations (k, l) in all periods j.





k∈K l∈L

Population dynamics constraints Susceptible individuals. Constraint (7) defines the number of susceptible individuals at the end of period j + 1 to be equal to the number of susceptible individuals in the previous period, plus the immigration into location (k, l), minus the emigration of susceptible individuals from location (k, l), and minus newly infected individuals as a function of I- and F-variables at the end of period j.

(k,l ) (k,l ) (k,l )

+ μ1

λ1 I j

(k,l )

+ μ2

(k,l ) (k,l )

γ1

Tj

Fj

(11)

Buried individuals. Constraint (12) determines the cumulative number of buried dead bodies resulting from the outbreak. ) B((k,l = B(jk,l ) + γ3(k,l ) Fj(k,l ) j+1 )

j ∈ J\J¯, k ∈ K, l ∈ L.

  (k,l ) (k,l )  (k,l ) (k,l ) gn j yn j + b1 j T j j∈J\{0,J¯} n∈N



≤ 1

(13)

j∈J

Treatment capacity constraint. Constraint (14) determines the total capacity (cumulative number of beds) in location (k, l) until the end of period j, including the initial capacity C0(k,l ) .

C (j k,l ) =

j  

(k,l ) κn ynr + C0(k,l )

j ∈ J\{0, J¯}, k ∈ K, l ∈ L.

(14)

r=1 n∈N

Constraint (15) provides the number of infected individuals who become hospitalized (and quarantined) based on the number of available beds in ETCs in location (k, l) at the end of period j. According to constraint (15), if the number of infected individuals, I (j k,l ) , is more than the remaining capacity of the ETCs in location (k, l), then only C (j k,l ) − T j(k,l ) number of infected individuals can be accommodated in the ETCs. Otherwise, if there are enough available beds in the ETCs, all infected individuals in location (k, l) will be admitted to the ETCs. Therefore, constraint (15) implies that the number of people who are hospitalized (and quarantined) for treatment (I¯) will be equal to the minimum between the number of infected people that seek treatment and the available number of beds.

I¯(j k,l ) = min{I (j k,l ) , C (j k,l ) − T j(k,l ) }



j ∈ J\J¯, k ∈ K, l ∈ L.

Note that the available capacity C (j k,l ) − T j(k,l )

(15)



is non-negative, as

implied by the non-negativity of I and I¯. Because constraint (15) is non-linear, it is linearized by defining additional binary variables and inequalities to represent the multiplication of continuous and binary variables (see Appendix B for details of linearization). 3.4.4. Non-negativity and integrality constraints Constraint (16) describes non-negativity restrictions on the number of susceptible, infected, treated, funeral, buried, and recovered individuals.

S(jk,l ) , I (j k,l ) , I¯(j k,l ) , T j(k,l ) , R(jk,l ) , Fj(k,l ) , B(jk,l ) ≥ 0

j ∈ J, k ∈ K, l ∈L. (16)

Constraint (17) imposes integer requirements on the number of type n ETCs to be opened in location (k, l) at the end of period j. Furthermore, if the number of infected individuals is less than 1, ) then the value of the integer variable yn(k,l is forced to be zero, i.e., j no ETC will be opened.

yn (jk,l ) ∈ {0, 1, 2, . . .};

) yn(k,l ≤ I (j k,l ) j

n ∈ N, j ∈ J\{0, J¯}, k ∈ K, l ∈L. (17)

(k,l ) (k,l )

− γ3

j ∈ J\J¯.k ∈ K, l ∈ L,

1051

(12)

3.4.5. Solving epidemics–logistics model The epidemics–logistics model presented in Eqs. (1)–(17) above is a mixed-integer nonlinear programming (MINLP) formulation due to the non-linear treatment capacity Constraint (15). We obtain an equivalent MIP formulation by replacing the non-linear Constraint (15) with equivalent linear Constraints (B-18b)–(B-18l)

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

1052

Table 1 Summary of Ebola treatment cost components (World Health Organization, 2014a). Fixed costs are given for a 50 (100)-bed ETC. Cost description

Fixed cost

Variable cost (per week)

Ebola treatment center Isolation unit center (IUC) Laboratory diagnosis Subnational technical services Contact tracing Safe burial Total

$ 386,0 0 0 (694,80 0)

$ 4405

$112,500

$566.41

$10 0,0 0 0

$270

Safe burial cost (per week)

$1125 $564 $ 598,500 ($1,077,300)

$6929.41

$563.63 $563.63

(see Appendix B). Thus, the epidemics–logistics model can be represented by the MIP model defined in Eqs. (1)–(14), (B-18b)–(B18l), (16), and (17). This MIP model integrates the disease dynamics and logistic considerations into one optimization model. The MIP formulation above is directly solved by CPLEX 12.6. Different than other models that simulate a number of logistics allocation scenarios in compartmental disease models, the MIP model envisions current and future growth of the disease as well as impacts of all possible resource-allocation scenarios on disease progression simultaneously, and determine the best strategy to minimize the number of infected and dead people as a result of the disease. In the next section, we demonstrate the use of the proposed MIP model in a case study involving control of the 2014– 2015 Ebola outbreak in the three most-affected West African countries. 4. Case study The optimization model introduced in Section 3 is employed on a real case study concerning the 2014–2015 Ebola outbreak in the three most affected West African nations. Data is gathered from the World Health Organization (WHO) and the literature regarding recent Ebola outbreaks in Guinea, Liberia, and Sierra Leone. These Sub-Saharan West African countries share a common border and have similar ethnic and cultural affinity. The broken and poor healthcare system in these three most-affected nations makes it difficult for infected and hospitalized individuals with Ebola to get the proper medical care they need, hence, contributing to the continuous growth of the Ebola epidemic. For example, before the Ebola crisis occurred, Liberia had only 0.1 physicians, 2.7 nurses, and 8 hospital beds per 10,0 0 0 people, while Sierra Leone had 0.2 physicians and 1.7 nurses per 10,0 0 0 people (United Nations, 2014a). In addition, when the outbreak started, no Ebola isolation units or ETCs existed in Guinea, Liberia, and Sierra Leone (World Health Organization, 2015). Due to the weak health infrastructure in those countries, we focus our attention on new ETCs to be built for providing service to Ebola patients. In our case study, we consider two types (50- and 100-bed) of ETCs to provide supportive care to infected people (United Nations, 2014a). The data set used for our numerical analysis includes intervention cost data (Table 1), geographic information, and population data for each region of Guinea, Liberia, and Sierra Leone (Table 2), as well as estimated disease-transition parameters (Table 3) (see Appendix A for a detailed explanation of the data, data sources, and assumptions). The periods in our MIP model are set to weeks, and thus both the daily transition and monthly cost data parameters are converted into weekly parameters. Each of the considered countries – namely Guinea, Liberia, and Sierra Leone – are divided into regions, as previously provided in Table 2 (see Appendix A for details). These six regions are repre-

sented in the form of a two-by-three matrix as input to the MIP formulation. We also consider subnational-level technical services and related logistics costs (Table 1) to ensure that infected people are transported to any of the ETCs in their region if the ETC capacity allows. We assume that country boundaries are closed for protection from Ebola epidemic, and thus migration only takes place within the regions of a county. Furthermore, we use the population ratio for each region of a country to compute the initial number of cases in that region. The population ratios used in this case study are also given in Table 2. The geographic population data presented in Table 2 is used as an estimate for the initial number of susceptible individuals, which provides an upper bound for the number of infected individuals. Migration rates among regions are also estimated by using the migration data presented in (Bay, 2003). Base case data includes economic, geographical, population, and epidemiological information, which is summarized in Tables 1–3. Initial conditions for the number of infected individuals in Guinea, Liberia, and Sierra Leone are obtained using the data reported in the WHO database (World Health Organization, 2014b; World Health Organization, Regional Office for Africa, 2014), and presented for different intervention start dates in Table 4. All cases with different intervention start dates (note each date is 4-weeks, approximately one month, apart) are solved for a 60-week planning horizon. 5. Results In this section, we validate the predictions of the MIP model with respect to the actual data and intervention strategies. We solve the full MIP model (1)–(14), (B-18b)–(B-18l), (16) and (17) over a 60-week period, from August 30, 2014, to October 24, 2015, using the case study data provided in Tables 1–4. All computational experiments are performed on a desktop computer equipped with Windows 7 with 3.4 GHz and 16 gigabyte memory, using CPLEX 12.6. We then present numerical results and analysis obtained by solving the proposed model for the three West African nations most affected by Ebola. Results imply that our MIP model provides a very good fit of Ebola outbreak empirical data and thus could be used to gain insights into the transmission dynamics of EVD. 5.1. Model validation The model parameter values are obtained from various sources in the literature (Bay, 2003; Camacho et al., 2014; Rivers, Lofgren, Marathe, Eubank, & Lewis, 2014; WHO, 2016b; WHO Ebola Response Team, 2014; World Health Organization, 2014a) as summarized in Tables 1–3. The proposed MIP model is then validated against the current outbreak data (World Health Organization, 2014b) in terms of the cumulative number of cases and funerals from August 30, 2014, to October 24, 2015. The actual number of ETCs opened in each considered region of the three countries and opening dates are used and fed into ) the model by fixing the value of the variable yn(k,l to the actual j number of ETCs with type n in period j and location (k, l). For example, five 100-bed ETCs were opened in Sierra Leone between December 13 and 20, 2014, which corresponds to the 15th week ) in our model. Consequently, the value of y2(k,l is set at 5, where ,15 (k, l) corresponds to the coordinate of Sierra Leone with respect to Guinea and Liberia. Therefore, the proposed model is solved ) with fixed yn(k,l values based on the actual ETCs opened (see supj plementary document providing information about actual number of ETCs opened, their location, and ETC opening times). Thus, the model is validated by adding additional parameters to the model, i.e., by fixing the number of treatment facilities and their capacities

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

1053

Table 2 Regions and population sizes in Guinea, Liberia, and Sierra Leone (Bay, 2003). Guinea

Population

Ratio

Liberia

Population

Ratio

Sierra Leone

Population

Ratio

Upper (G1 ) Middle (G2 ) Lower (G3 ) Total

4,308,494 2,670,567 3,649,911 10, 628, 972

0.41 0.25 0.34 1.00

Northern (L1 ) Southern (L2 )

2,233,091 1,243,517

0.64 0.36

Sierra Leone (SL)

4,976,871

1.00

3, 476, 608

1.00

4, 976, 871

1.00

Table 3 Summary of model parameters obtained from literature for Ebola outbreak in Guinea, Liberia, and Sierra Leone. Parameter

Description

Data

Reference

Guinea

Liberia

Sierra Leone

0.40 0.35 0.60 0.93

0.39 0.29 0.61 2.22

0.31 0.24 0.69 2.50

Estimated from WHO (WHO, 2016b) Estimated from WHO (WHO, 2016b) Estimated WHO Ebola Response Team (2014)

2.50 1.00

2.63 2.27

2.85 2.50

WHO Ebola Response Team (2014) WHO Ebola Response Team (2014)

γ3(k,l ) σ1(k,l )

Fatality proportion without treatment Fatality proportion with treatment Survival proportion without treatment Average time (in weeks) from symptoms onset to funeral Average time (in weeks) from symptoms onset to recovery Average time (in weeks) from treatment to funeral (death) Average time (in weeks) from treatment to recovery Safe burial rate Transmission rate in community

1.56 0.73 0.33

2.27 0.74 0.28

2.37 0.71 0.32

σ2(k,l )

Transmission rate at traditional funeral

0.68

0.68

0.73

WHO Ebola Response Team (2014) Rivers et al. (2014) Estimated from Camacho et al. (2014); Rivers et al. (2014) Estimated from Camacho et al. (2014)

μ1(k,l ) μ2(k,l ) μ3(k,l ) 1/λ1(k,l ) 1/λ2(k,l ) 1/γ1(k,l ) 1/γ2(k,l )

Table 4 Number of initial infections for different intervention start dates (World Health Organization, 2014b; World Health Organization, Regional Office for Africa, 2014). Country

Number of initial infections (in 2014) August 30

September 27

October 25

November 22

Guinea Liberia Sierra Leone

218 684 604

425 1627 1416

597 2217 2396

879 4152 5201

Table 5 Statistical analysis to compare model predictions of weekly cases and funerals and actual outbreak data. Country

Mean

Two-tailed paired-t-test

Outbreak

Predicted

t-stat

t-critical

p-value

2.002

0.478 0.970 0.589

Cases

Guinea Sierra Leone Liberia

52.6 217.7 154.9

54.2 216.1 167.3

0.713 0.038 −0.543

Funerals

Guinea Sierra Leone Liberia

35.1 65.4 68.6

34.6 71.7 72.5

0.290 −0.611 −0.905

0.770 0.544 0.369

based on their real opening dates. Consequently, the predicted outbreak data, which is the outcome of the model, is compared with the real outbreak data given in the WHO database (World Health Organization, 2014b). Fig. 2 provides a visual comparison of the predicted and outbreak data for Ebola cases and fatalities in Guinea, Liberia, and Sierra Leone between August 30, 2014 and October 24, 2015. Fig. 2 demonstrates that the proposed model provides a good-fit for the cumulative number of cases and funerals in all three countries. Although the cumulative predicted cases and fatalities are slightly overestimated in Liberia and Sierra Leone, paired-t-test results prove that the proposed model provides statistically similar results with respect to the outbreak data for the time period between August 30, 2014, and October 24, 2015 (Table 5). As opposed to comparing the cumulative data points, the paired t-test assists in analyzing the differences between the pairs of weekly predicted

and the actual data. Table 5 shows that there is no statistical difference between the predicted weekly number of new cases and new fatalities, and the actual outbreak data, because all p-values are greater than 0.05. 5.2. Impact of intervention budget on Ebola epidemic in Guinea, Liberia, and Sierra Leone In this subsection, the impact of imposing different budget levels for treatment intervention on the total number of infections and deaths is analyzed. The columns of Table 6 report the total budget available, model solution time, optimality gap, objective function value, countries, their corresponding regions, optimal treatment budget allocations determined by the MIP model, number of 50- and 100-bed ETCs, and cumulative number of treated, infected, deceased, and recovered individuals in each of the six regions. All those values in Table 6 are obtained by solving the MIP model for different budget limitations with a time limit of 3600 CPU seconds. As shown here, the model is difficult to solve under a tight budget limitation ($30 million and $45 million). This is because the increase in budget allows for the allocation of more ETCs, thus making it more difficult for the solver to decide on the optimal ETC assignment among multiple possible locations. However, when the budget is ample, resources are provided to all locations that are in need of treatment, and thus the spatial allocations of ETCs and treatment resources are optimally decided within a second. In the absence of control intervention, the cumulative number of infected individuals and deaths increases exponentially to high levels in the considered countries (Table 6). Specifically, the number of cases is estimated to peak at 25,338 for Guinea, 39,002 for Liberia, and 163,853 for Sierra Leone, respectively, under no control intervention. However, when control intervention is introduced in the form of an available budget to set up ETCs, the impacts of the outbreak are reduced considerably (Table 6) due to the rapid isolation and treatment of the cases, thus substantially reducing the total number of infected and deceased individuals. If a budget limitation of $30 million is given, each region of Guinea as well as L2

1054

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

Fig. 2. Comparison of model predictions with actual outbreak data for cumulative number of infected cases (I) and funerals (F) in Guinea, Liberia, and Sierra Leone when opened ETCs are employed in the model.

region of Liberia receives a 50-bed treatment facility, while Sierra Leone is assigned four 100-bed treatment facilities, thus reducing the total number of fatalities significantly. It can also be observed that the ETC allocation scheme changes in all three countries and their regions for different budget limitations, thus implying that the utilization of an optimization model is critical for efficient resource allocation during an infectious disease outbreak. Finally, out of a total budget of $60 million, only $53.9 million has been optimally allocated to all six regions in order to open two 50 and sixteen 100-bed ETCs for isolating and treating infected individuals. From this budget, $3.3 million, $1.6 million, and $2.5 million are allocated to regions G1 , G2 , and G3 , respectively,in order to open and operate two 50-bed and two 100-bed ETCs in Guinea. Interestingly, a similar amount of overall budget ($23.5 million) is allocated to both Liberia and Sierra Leone, to open seven 100bed ETCs in each country until the end of the 60-week planning horizon. We observe that as the available budget increases, the model tends to allocate more 100-bed ETCs in order to treat more infected people while also saving from the fixed cost of opening ETCs. Therefore, the cumulative number of treated individuals increases proportional to the available budget, reaching to 5129, with an ample treatment budget of $60 million. In addition, the cumulative number of infected individuals and deaths decrease by approximately 97% of the total number of infected and dead individuals under no treatment. This is because control interventions, such as quarantine and treatment, effectively reduce contact with infected people, thus decreasing disease transmission among the population and eventually slowing the spread of the epidemic. 5.3. Impact of different intervention starting dates In this section, we solve the MIP model presented in Section 3 for different intervention starting dates (August 30,

September 27, October 25, and November 22) with an ample budget level. We then compare the impact of optimal intervention strategies in each case with respect to the resulting number of infected individuals and deaths in Guinea, Liberia, and Sierra Leone during a 60-week planning horizon. Original outbreak data is also provided in order to evaluate the differences between optimal intervention results and real data. In all cases, we observe that treatment considerably lessens the disease’s impact, compared to the no-treatment case, regardless of the starting date (Fig. 3). However, the earlier the starting date of the intervention, the fewer the infections and funerals. In particular, we observe a considerable reduction in infected and deceased individuals if the intervention starts on August 30 in all three countries. Fig. 3 shows that model predictions for different starting dates clearly deviate from each other for each of the three countries. All intervention strategies reduce the number of infections and lessen the exponential growth of the disease. However, starting treatment intervention on November 22 is too late for all countries compared to the original outbreak data. Despite a sufficient budget allocation for ETCs and treatments, intervention starting on November 22 leads to a notably higher number of infected and deceased individuals compared to intervention starting in August, September, and October in all countries. Therefore, as intervention is postponed further, the number of cases and fatalities can increase significantly, even with intense intervention efforts. For all intervention starting dates except November 22, Liberia exhibits the highest number of infections and deaths, while Sierra Leone exhibits the number of infections when intervention is delayed to November 22, due to its higher funeral transmission rates and longer infection-to-funeral period. This result implies that applying an early intervention strategy is particularly important for

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

1055

Table 6 Optimal allocation of budget and ETCs, and corresponding number of cases in Guinea, Liberia, and Sierra Leone for 60 periods when treatment starts on August 30, 2014. Note that burial costs are not included in the OptBudget. Budgeta

Timeb

Gapc

Objective

Country

Regiond

OptBudgete

50 / 100-bed

$0

2

0

609,917

Guinea

G1 G2 G3 L1 L2 SL − G1 G2 G3 L1 L2 SL − G1 G2 G3 L1 L2 SL − G1 G2 G3 L1 L2 SL −

$0 $0 $0 $0 $0 $0 $0 $2,281,357 $1,624,893 $2,008,093 $0 $3,078,268 $20,783,428 $29,776,039 $0 $945,002 $2,498,961 $12,976,401 $7,761,174 $20,783,428 $44,964,966 $3,333,545 $1,624,932 $2,481,711 $14,553,584 $8,931,471 $23,650,003 $53,990,246

0/0 0/0 0/0 0/0 0/0 0/0 0/0 1/0 1/0 1/0 0/0 1/0 0/4 4/4 0/0 1/0 0/1 0/2 1/1 0/4 2/8 1/1 1/0 0/1 0/4 0/3 0/7 2/16

Liberia

$30 million

3600

0.06

107,764

Sierra Leone Total Guinea

Liberia

$45 million

3600

0.38

40,386

Sierra Leone Total Guinea

Liberia

$60 million

2

0

4694

Sierra Leone Total Guinea

Liberia Sierra Leone Total a b c d e f g h i j

f

Treatedg

Infectedh

Funerali

Recoveredj

0 0 0 0 0 0 0 242 148 203 0 357 2377 3085 0 50 205 1561 878 2378 5072 239 148 202 1478 822 2240 5129

10,294 6320 8724 25,150 13,852 163,853 228, 193 559 333 460 25,149 13,580 2202 41, 724 10,282 6333 458 1897 1098 2202 22, 270 538 333 456 1539 857 1936 5659

6257 3841 5302 10,258 5694 47,129 78, 481 387 232 319 10,257 5559 914 17, 668 6249 3849 318 975 560 914 12, 865 374 231 317 825 459 829 3038

3411 2093 2890 12,959 7134 90,966 119, 453 172 101 141 12,959 7061 1289 21, 723 3406 2097 140 923 538 1289 8393 164 101 139 714 397 1106 2621

Overall available budget for treatment, Solution time in CPU seconds, |Best integer solution-best lower bound|/best lower bound, Region represents divisions within each country, OptBudget represents optimal total treatment budget allocation by the MIP model over all regions and periods, Total number of 50- and 100-bed ETCs established, Cumulative number of treated individuals, Cumulative number of infected individuals, Cumulative number of dead individuals, Cumulative number of recovered individuals.

locations where disease transmission is high and the duration is long. Our results have a one-to-one correspondence with the basic reproduction number, R0 , which is defined as the expected number of secondary cases of a typical single infected individual during his or her entire infectious period, in a population which is entirely susceptible (Heffernan, Smith, & Wahl, 2005). If R0 < 1, disease transmission is not self-sustaining and is unable to generate a major epidemic. On the other hand, an epidemic is likely to occur whenever R0 > 1. If R0 = 1, each existing infection is causing one new infection, and the disease will stay alive without epidemic (Chowell et al., 2004). By observing the cumulative growth function of cases in Fig. 3, we can infer about the value of R0 . As shown in Fig. 3, R0 is constantly changing for each country, time period, as well as intervention starting date, and thus is based on the capacity allocated for disease intervention. For example, based on the original outbreak case data in Sierra Leone in Fig. 3, when the disease is first exponentially growing, we can infer that R0 is larger than 1; however R0 starts to decrease at the inflection point of the cumulative growth function, and becomes less than or equal to 1 when the growth function reaches an asymptote, i.e., the epidemic stops growing. Based on the results of the intervention that starts on November 22, 2014 in Sierra Leone, R0 is larger than 1 at the beginning of the intervention period, and then it reduces to 1 by March 28, 2015. Similar interpretations apply to model predictions on the cumulative number of cases and fatalities under varying intervention starting dates for Liberia and Guinea.

5.4. Budget allocation based on optimal solution In this section, we present the optimal budget allocation for treatment components (fixed and variable costs of treatment) and safe burial for Guinea, Liberia, and Sierra Leone for different intervention starting dates (August 30, September 27, October 25, and November 22), as shown in Fig. 4. Again, for this experiment, we set the budget limitation for both the treatment and safe burials to ample values in the MIP formulation in order to determine the total budget amount required to minimize total infections and deaths under each case with different intervention starting dates. Here, while we do not minimize the allocated budget, we compute the total budget (e.g., cost of capacity allocations) by fixing the optimal solution in the left-hand side of the budget constraint (13). Similarly, letting Fj∗(k,l ) represent the number of buried people based on the MIP model solution, we compute the total allocated budget (k,l ) (k,l ) ∗(k,l ) for funerals using b2 j γ3 Fj . j∈J k∈K l∈L

Fig. 4 suggests that the major Ebola cost driver for all four different starting dates of control intervention is the variable treatment cost (e.g., operational and treatment costs in ETC and isolation units, laboratory diagnosis, subnational technical services, and contact tracing), followed by the fixed cost of treatment, and finally burial cost. When control intervention is introduced on August 30, the model solution allocates two 50-bed and two 100bed ETCs to Guinea, and seven 100-bed ETCs to both Liberia and Sierra Leone, with an optimal budget allocation of $54.83 million including burial costs. In this case, we observe that the variable

1056

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

Fig. 3. Original outbreak data and model predictions for cumulative number of the cases (I) and funerals (F) in presence of control strategies in Guinea, Liberia, and Sierra Leone when treatment starts on August 30, September 27, October 25, and November 22, 2014 (note ten-week intervals on x-axis) with optimal treatment budget allocation.

treatment cost accounts for approximately 64.8% ($35.5 million), while the fixed cost of treatment is 33.6% ($18.4 million) and the burial cost accounts for 1.6% ($0.84 million) of the optimal budget allocated. If control intervention starts late in the course of the epidemic, then the cost of treatment is expected to increase significantly. For example, we observe that if treatment starts on September 27, then a minimum budget of $123.51 million is required to reduce the initial 3468 infections over the planning horizon. In addition, if we introduce control interventions on October 25 (November 22), then the model solution allocates an optimal budget of $180.9 million ($357.1 million) with an initial 5210 (10,232) cases in order to minimize the total number of infections and deaths. This result shows that for each month of delay in intervention, the re-

quired optimal amount of budget to control the epidemic nearly doubles. Although a significantly larger treatment budget is allocated in delayed-intervention cases, total infections and deaths by the end of the planning horizon are considerably larger compared to early treatment efforts (Fig. 3). Our results also suggest that the MIP model allocates the majority of the budget in the form of treatment cost and cost of allocating ETCs to countries by taking the initial number of infections as well as the future potential number of infections and deaths in each country into account, as shown in Fig. 4. For instance, the August and September interventions reveal that the model allocates most of the budget to Liberia, where the initial number of infections is higher than the number of infections in Guinea and Sierra Leone. On the other hand, when treatment starts on

Fig. 4. Cost analysis ($US in millions [M]), optimal budget allocation, and total required budget (TRB) for different starting dates of control intervention in Guinea, Liberia, and Sierra Leone, where control interventions start on August 30, September 27, October 25, and November 22, with ample treatment budget.

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

1057

Fig. 5. Optimal number and spatial allocation of ETCs represented by 50- and 100-bed ETC symbols at end of periods (weeks) one and two when intervention starts on August 30, September 27, and October 25 in 2014. Total number of open ETCs remains same after period 2 until end of period 60. The number next to “x” represents the number of ETCs opened for the corresponding type of ETC.

October or November, most of the budget is allocated to Sierra Leone, due to the highest initial number of infections as well as high disease transmission rates that aggravate the number of infections in this country.

5.5. Optimal spatial allocation of ETCs Fig. 5 illustrates the optimal number and location of 50- and 100-bed ETCs for August 30, September 27, and October 25 intervention starting dates. We demonstrate the optimal size and location of ETCs at the end of the first and second periods over a 60-period horizon, because the number of ETCs does not change after the second period. We note that the MIP model envisions current and future growth simultaneously, thus optimally allocating resources over time and various regions based on the severity of the outbreak. In Fig. 5, we observe that 50- and/or 100-bed facilities are established in each region of Guinea, Liberia, and Sierra Leone at the end of both period one (September 5) and two (September 12), right after intervention begins on August 30, while the total number of open facilities remains the same until the end of the planning horizon. When intervention is postponed to September 27 and October 25, the number of allocated ETCs increases in order to suppress the disease and prevent further spread. In particular, as the number of initial infections increases with the delay of intervention, more 100-bed ETCs are allocated than 50-bed facilities in an effort to treat more cases and save from the fixed cost of opening ETCs. Fig. 5 also demonstrates that the number of ETCs increases in Sierra Leone and Liberia from periods one (October 6) to two (October 13) when intervention starts on September 30, while it remains the same after ETCs are allocated at the end of period one (November 2) in the October 25 intervention case. The results of this experiment imply that most ETC allocations are made early in the planning horizon in order to immediately reduce the size of the epidemic and most efficiently use treatment facilities throughout the planning horizon.

Table 7 Number of initial infections for different intervention start dates. Region

Population

Number of initial infections (in 2014) June 16

August 30

November 22

G1 G2 G3 SL L1 L2

4,308,494 2,670,567 3,649,911 4,976,871 2,233,091 1,243,517

30 14 22 49 15 0

89 55 75 604 439 245

367 210 297 5201 2657 1495

Table 8 Sensitivity analysis for migration for early border closure. With migration

June 16

G1 G2 G3 SL L1 L2

Without migration

Infected

Funeral

Infected

Funeral

254.65 143.17 205.84 171.7 30.98 18.7

174.4 97.46 140.32 81.03 20.8 12.19

254.64 143.2 205.47 171.7 30.98 0

174.39 97.49 140.31 81.03 20.8 0

5.6. Impact of intervention with and without migration Computational experiments demonstrate that isolating and treating infected individuals in ETCs early in the planning horizon eventually reduces the number of infections to 0 by the end of May 9, 2015 (period 36), when intervention efforts begin on August 30 and migration is allowed among regions within the same country. All of the regions are completely cleared of the disease which demonstrates that EVD can be eradicated with the timely establishment of optimal number of ETCs as well as continuing treatment and isolation efforts. Table 7 provides population information and number of initial infections with different dates. Table 8 shows the impact of migration and border closure early in the epidemic. The column “With

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

1058

Table 9 Sensitivity analysis for migration among regions of countries. With migration Migration rate

Without migration Base

x10

x50

x100

0

Infected

Funeral

Infected

Funeral

Infected

Funeral

Infected

Funeral

Infected

Funeral

Aug 30

G1 G2 G3 SL L1 L2

539 334 457 1936 1540 857

374 232 317 829 826 460

540 333 458 1936 1541 856

375 232 318 829 827 459

542 329 459 1936 1547 855

376 228 320 829 830 457

544 323 462 1936 1555 846

378 224 321 829 835 451

539 335 457 1936 1540 857

374 232 317 829 826 460

Nov 22

G1 G2 G3 SL L1 L2

2040 1166 1651 14420 8541 4803

1357 776 1098 4902 4225 2376

2041 1162 1652 14, 420 8549 4795

1359 773 1100 4902 4231 2370

2050 1145 1661 14, 420 8585 4777

1365 760 1107 4902 4254 2354

2060 1124 1672 14, 420 8630 4723

1373 743 1115 4902 4284 2320

2039 1166 1650 14, 420 8540 4804

1357 776 1098 4902 4225 2376

Migration” represents the original model results without any borders between regions, while “Without Migration” column represents the results when borders between the regions of Guinea and Liberia are closed, i.e., no movement among regions within a country is allowed. We also perform sensitivity analysis on migration rates between the regions of Guinea and Liberia by multiplying the base migration rate by 10, 50, and 100, and present the resulting number of infected and funerals under migration as well as without migration in Table 9. Results suggest that border closure may increase or decrease the epidemic’s burden in a region, depending on the net migration rate, which is the difference between the number of immigrants and emigrants. In particular, border closure reduces the number of infected individuals in highly populated regions with a positive net migration rate, i.e., the number of people coming into the area is larger than the number of people leaving the area (e.g., G1, G3, and L1). On the other hand, border closure may further increase the number of infected individuals in low-population regions with a negative net migration rate, i.e., the number of people coming into the area is smaller than the number of people leaving the area (e.g., G2 and L2). This is because border closure hinders the movement from the poor to the developed regions, in turn preventing infected individuals from benefiting from treatment centers in the developed regions, thus preventing a more effective utilization of already-opened ETCs. Closing borders only helps if a region receives more migration from other parts of the country than departure from it. Furthermore, closing borders among regions within a country is beneficial when the infection is newly observed in a given region and has not spread to other locations, such is the case of L2 (Table 8). As shown in Table 8, the aforementioned impacts of migration compared to border closure escalate as the movement among regions of Guinea and Liberia intensifies.

5.7. Parameter sensitivity analysis We also perform a sensitivity analysis to identify key parameters influencing the growth of the disease and to provide insights into effective control strategies. This analysis includes community and funeral transmission rates, transition rate from I to F, and burial rate under no treatment. We also analyze parameters under treatment, such as fatality rate as well as transition rates from T to F and from T to R. In this analysis, either each considered parameter or each rate in the base case is varied, while other parameters are held constant to assess their impact on the total infections and deaths in each experiment (see supplementary material for detailed results for sensitivity analysis on model parameters).

Our analysis suggests that community transmission and funeral transmission rates are the most sensitive parameters in the model, followed by burial rates. Under control measures where ample treatment resources are applied, total infections and deaths considerably decrease, and infected individuals are kept in isolation. Therefore, transition rates under treatment are observed to be less significant compared to the community and funeral transmission parameters. Our analysis suggests that in order to stop the spread of Ebola and consequently reduce fatalities, policymakers must focus primarily on reducing transmission and contact rates aside from the rapid introduction of treatment efforts. Reducing community and funeral transmission rates by isolating infected people and changing human behavior through community education can help to bring the epidemic under control. Further increasing resources for faster burial of Ebola-related dead bodies and public information campaigns to seek rapid intensive care would have a positive effect on reducing the impact of the epidemic.

6. Discussion In this paper, we introduced a new epidemics–logistics MIP model that accounts for the spatial spread of an infectious disease outbreak while providing optimal resource allocation strategies for intervening an epidemic. This model considers realistic aspects of an epidemic disease, such as geographically varying parameters, migration of infected individuals over different regions, impact of limited resources on disease dynamics, and thus a varying treatment rate based on the capacity of treatment facilities. The objective of the proposed MIP model is to minimize the total number of infections and fatalities under a limited budget over a finite planning horizon. The proposed model optimally allocates the budget as well as the ETCs over time among the three most-affected countries simultaneously, which would only be possible by solving the optimization model. The model is also used to evaluate the impact of a large number of intervention and budget allocation strategies on controlling an epidemic. We illustrate the model performance on a case study of the 2014–2015 Ebola outbreak in Guinea, Liberia, and Sierra Leone. We analyze the spatial and temporal spread of the disease under no-control and control interventions, provide the optimal number and location of treatment centers, and analyze both the cumulative number of cases and deaths when intervention efforts start on different dates and under different budgets. Computational results help identify key problem parameters, such as the total treatment budget, intervention start dates, disease transmission rates, and burial rates, and their impacts on the total number of infected individuals and deaths due to the Ebola epidemic.

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

Furthermore, we validate our MIP model using actual ETC opening dates, and locations and number of ETCs opened over the planning horizon. We assess the model performance by applying a paired-t-test, which analyzes the mean differences of outbreak and predicted data by comparing the paired data points. The statistical analysis as well as visual comparison of outbreak and predicted data demonstrate that the model predictions in terms of cumulative number of cases and deaths closely fit the outbreak data from the WHO Ebola response roadmap situation reports. Our validation results further confirms that our model is a correct representation of the growth of the disease under actual treatment intervention and capacity allocations. Thus we validate both the correctness of disease transmission and logistics constraints presented in the epidemics–logistics model. Our model differs from other agent-based simulation models in a number of ways. For example, agent-based models analyze the interactions among individuals within a relatively small population, while our model focuses on interactions among larger populations. Agent-based models can capture individual-level interactions as well as heterogeneity among disease compartments, such as age-specific transmission rates. However, the complexity of a simulation model also escalates with population size. Therefore, while agent-based models are helpful for analyzing smaller populations and the fine details of disease transmission, optimization models could help with capturing the disease dynamics over larger-size populations and bigger landscapes. Furthermore, agent-based simulations evaluate the effect of interventions one at a time, while our optimization model oversees the impact of all possible interventions and budget allocation scenarios on the growth of the disease simultaneously and determines the optimal intervention strategy to minimize resulting infections and fatalities. Our computational study provides a number of important results and insights for controlling the Ebola outbreak. Without intervention, the number of cases and deaths is shown to grow exponentially. This result suggests that, if left uncontrolled, most of the susceptible population in Guinea, Liberia, and Sierra Leone will be impacted by the EVD, and the Ebola epidemic will possibly move to other countries in Africa and the world, thus leading to a global disaster, as also speculated in the reports and publications of WHO (Team, 2016; WHO, 2016a). We observe that the actual number of infections and deaths due to the EVD exhibits a logistic growth form in all countries. That is, first the disease shows an exponential growth and then, over time, slows down, reaching an asymptote. In particular, the implementation of optimal control intervention at the initial stage of the epidemic, August 30, 2014, eradicates the infestation in Guinea, Sierra Leone, and Liberia by May 9, 2015. All results from the case study consistently show that early control reduces infection and fatalities considerably more compared to delayed intervention. While a late implementation of control intervention requires intense response efforts and a significantly larger treatment budget, it still leads to a substantially greater cumulative number of infections and deaths compared to the case when intervention is early. Hence, the immediate isolation and hospitalization of infectious individuals at the onset of the disease are key strategies in reducing the number of infections and controlling the Ebola epidemic. The necessity of immediately intervening the EVD at the onset of the epidemic is also declared as one of the most important lessons learned from the 2014–2015 Ebola outbreak by WHO, CDC, and United Nations Educational, Scientific and Cultural Organization (UNESCO) (Team, 2016; WHO, 2016c). Our results regarding the significance of rapid intervention also support the study of Legrand et al. (2007), who shows that the control of past Ebola epidemics in Democratic Republic of Congo (DRC) in 1995 and in Uganda in 20 0 0 was heavily dependent on the timing of the intervention. Those results are also confirmed by the case in which the 2014 Ebola outbreak in Nigeria was successfully contained

1059

because of early intervention and the quick isolation of cases (Fasina et al., 2014). Contact tracing and increasing healthcare and community services to ensure prompt diagnosis are useful strategies that could help with the rapid hospitalization of cases and thus prevent further spread of the disease. Our experimental results also highlight the fact that taking an early action of control is particularly important for countries where the disease progression is fast, such as in the Liberia case. This result is consistent with the UN’s report showing that most of the funds for Ebola response efforts is allocated to Liberia in order to curb the epidemic (United Nations, 2015b). The UNESCO also emphasizes that the scope and budget of activities foreseen in Liberia are expected to be proportionally larger than that of Sierra Leone and Guinea, despite the fact that the overall population is smaller (UNESCO, 2014). The organization explains this result by the continuous highest incidence rate of Ebola in Liberia, which is compatible with our results. Results also imply that the number of infected individuals depends not only on the timing of the intervention but also available budget determining the treatment capacity. Among all intervention options, a major portion of the intervention budget is optimally allocated to treatment (medical treatment, isolation, diagnosis, and contact tracing), while the fixed cost of establishing ETCs and cost of safe burials are significantly less compared to variable treatment expenses. When the budget is limited, available funds are distributed among countries based on the initial number of infections, population density, disease transmission rates, and geographical extent of the country. Computational results also show that the optimal solution can change considerably for varying budget limitations, thus highlighting the importance of utilizing an optimization model over intuitive decisions, which could lead to suboptimal strategies. It is interesting to note that as of May 31, 2015, according to the United Nations (UN) report, $409.04 million, $1,013.12 million, and $653.2 million have been deployed to Guinea, Liberia, and Sierra Leone as funds for Ebola response efforts to support these affected countries (United Nations, 2015b). In Fig. 4, we observe that the actual amount of allocated funding is far beyond the suggested budget allocation for the intervention starting on November, 22, 2014, implying that a large amount of intervention funds was provided to these countries after the end of November 2014 (United Nations, 2015a). This result is also supported by the actual infection growth, which shows an exponential trend for the first 15 weeks after August 30, 2014, until December 13, 2014. Therefore, total funding provided for the most affected countries arrived too late for controlling the epidemic, leading to a significant number of Ebola cases and deaths as a result of the Ebola outbreak. Our finding has also been supported in many reports written about lessons from the Ebola epidemic (Grépin, 2015; Moon et al., 2015). The analysis of critical problem parameters provides insights into effective intervention strategies. In particular, disease transmission rates, burial rate, and recovery rate are found to be vital parameters affecting the number of infections and deaths. Therefore, the immediate isolation of cases to reduce community contact rates, advances in medication that improve recovery rates, precautions taken during traditional funeral practices to avoid direct contact with the dead body, safe burials, and increased educational campaigns to improve Ebola disease awareness can effectively reduce disease transmission rates and the number of secondary infections, thus helping to take control of the epidemic. Our model could be used as a decision support tool to aid policymakers in understanding disease dynamics and making the most effective decisions to fight and eradicate epidemics in a timely manner. One specific policy implication of the case study is that large upfront investments in terms of treatment facilities with isolation units result in the most efficient use of resources to

1060

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

minimize infections and take control of the Ebola epidemic. This result holds, even when intervention efforts are delayed. Lessons learned from Ebola (Team, 2016; WHO, 2016c) highlights the importance of capacity allocation as a quick response strategy to intervening in the epidemic, as shown by our model results. In addition to providing isolation units in the ETCs, confinement of the disease can also be achieved by establishing communitybased clinics to ensure that sick individuals are locally separated from the public, and also by closing borders to neighboring countries in order to reduce the spread of the disease to other nations (Chen et al., 2014). On the other hand, our results suggest that limiting movement within regions of a country may further increase the number of infected individuals because borders within a country could also limit the effective utilization of established ETCs. This result is consistent with the UN’s call for lifting general travel and border restrictions that could possibly undermine the efforts to respond to the Ebola outbreak (United Nations, 2014b). On the other hand, our results suggest that enforcing borders between regions within a country would only help to reduce the spread of the disease when the epidemic is at its initial stages and has not spread to other regions. This result could support the action of Liberia’s government to encourage residents in three quarantined areas of the country to limit their movements to help prevent the spread of the deadly virus (Radio France Internationale, 2014). Furthermore, larger regions, which include capitals and attract many people from outside, may also benefit from closing their borders, rather than smaller regions where more people emigrate than immigrate, as consistent with the CDC’s initial response to the Ebola epidemic in West Africa by taking immediate border health measures (Cohen, 2016). 7. Conclusion and future work In this paper, we present a new epidemics–logistics optimization model. Our MIP approach addresses limitations of the current epidemiology and OR literature by doing the following: (1) integrating logistics issues into a novel spatially explicit epidemic model; (2) considering the capacity of treatment facilities in an epidemic model, thus addressing the limiting assumption of constant treatment rates used in traditional compartmental models; (3) using spatially varying rates for disease transmission; (4) incorporating the migration of susceptible and infectious individuals among regions of a given country; and (5) providing a computational advantage over continuous epidemic models that use differential equations. We illustrate the use of the model in determining optimal intervention strategies for controlling the 2014–2015 Ebola epidemic. However, our epidemics–logistics MIP approach is general, and given the available data, it could be applied to a wide range of diseases, including smallpox, influenza, and tuberculosis, in order to identify the most efficient control strategy for controlling various epidemic diseases. Future research could address some of the limitations in both epidemiological and logistics parts of the proposed model. For example, our model assumes a homogenous population in each region of the considered countries. In a future extension of the model, the heterogeneity of contact rates within each region could be considered. In particular, age-specific characteristics, such as transition rates among different age groups, could be incorporated into the model (Büyüktahtakın, Kıbıs¸ , Cobuloglu, Houseman, & Lampe, 2015). Disease transmission within hospitals could also be considered separately by incorporating a new compartment to represent hospital infections. Although we have not considered this in our study, our epidemic model could easily be extended by adding a compartment class of exposed (infected but not yet infectious) individuals.

As proven by our validation approach, our model is a valuable tool for decision making. It is also important to note that the predictions of even the best optimization models are as good as the data provided to them. Data may vary over time and spatial location, and those variations could possibly impact some of the model results. One possible future approach to handle the potential pitfalls regarding the volatility in data is to incorporate uncertainty explicitly in the model parameters. This could be achieved by extending the model into a scenario-based stochastic optimization formulation, which includes various scenarios regarding the uncertain parameters with their corresponding probabilities. Some of the assumptions made in the logistics part of the model could also be handled in future work. For example, treatment and intervention funds are rarely available at the onset of the disease, and even if funds are immediately available, their arrival time to each region of the affected countries is uncertain. A stochastic version of the model could incorporate the potential delay in the arrival time of the intervention budget. Future studies could also investigate minimizing both total number of infected individuals and the cost of intervention efforts simultaneously. Nevertheless, proposed future considerations will further increase the complexity of the model, necessitating specialized algorithms, such as cutting planes and decomposition approaches, to solve the problem optimally over large spatial and temporal scales. Acknowledgments We gratefully acknowledge the support of the National Science Foundation CAREER Award under Grant # CBET-1554018. We are also grateful for the helpful comments of three anonymous referees and the editor. Appendix A A. Details about case study data and parameter estimation A.1. Economic data This subsection provides details regarding cost data used in the case study. All cost data is estimated from the WHO Ebola response roadmap report (World Health Organization, 2014a). A summary of the weekly cost data including fixed cost, variable treatment cost for infected individuals undergoing treatments, unit burial cost for dead individuals, and unit tracing cost for infected individuals under treatment is provided in Table 1. Below, we describe each of these cost components in detail. Fixed cost data: The primary drivers for fixed cost in our optimization model are opening 50- and 100-bed ETCs, isolation unit centers (IUCs) within each ETC, and laboratory diagnosis centers. The fixed cost for establishing a 50-bed ETC to accommodate infected individuals for treatment is estimated to be $386,0 0 0. Another relevant component is the cost of establishing an IUC where infected individuals are placed in isolation, which is estimated to be $112,500. Finally, setting up the laboratory diagnosis for sample testing has a onetime fee of $10 0,0 0 0. We further assume that the fixed cost for opening a 100-bed ETC including IUC and laboratory diagnosis to be 80% more than the fixed cost of opening a 50bed ETC, as provided in Table 1. Variable treatment cost data (operational data): Our variable treatment-cost data mainly consists of the operation and maintenance (O&M) cost for ETCs, IUCs, laboratory diagnosis, subnational technical coordination, and contact tracing per infected person and per week. The results of our survey indicated a range of O&M costs for these services. First, these costs include $881,0 0 0 per month for running a 50bed ETC where hospitalized individuals receive treatment.

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

Second, the O&M cost for the IUC is also estimated to be $2266 per month per infected individual. Third, we assume that $129,600 would be needed to operate and maintain the laboratory diagnosis, estimated at 600 procedures per month. Based on this cost scheme, laboratory diagnosis costs $270 for an estimated five procedures for each patient per week. Fourth, the variable cost for a subnational-level technical services (information management, technical expertise, social mobilization, and logistics management) is assumed to be $450,0 0 0 per month to cover a 100-bed ETC facility based on the WHO situation report estimation. Fifth, according to their report, approximately $225,450 is needed each month to track and process one thousand contacts using a team of one hundred people. We assume an estimated contact-tracing cost of $564 for tracing ten suspected contacts of an infected individual. Safe burial cost data: For each Ebola-related fatality, it is necessary that safe burial and household disinfection procedures are thoroughly carried out by a team of trained personnel. The WHO report states that approximately $56,363 per week is needed to perform one hundred Ebola-related safe-burial procedure. A.2. Population and geographical data Table 2 provides the most updated population size in Guinea, Liberia, and Sierra Leone obtained from the work of Bay (2003). For this case study, we partition the three countries into six regions based on the ratio of their square mile area relative to Sierra Leone. This leads to Guinea being split into three regions, Liberia into two regions, and Sierra Leone being classified as only one region. The given data set shows that Guinea has the largest population size of 10,628,972, with an area of 94,926 square miles and consisting of three geographic regions: upper, middle, and lower. Upper Guinea consists of Boké, Conakry, and Kindia, with a total population of 4,308,494, while Middle Guinea consists of Labé, Mamou, and Faranah, with a combined population of 2,670,567. Lastly, Lower Guinea is classified as Kankan and N’zérékoré, with a total population of 3,649,911. Liberia, which has a population of 3,476,608 and an area of 43,0 0 0 square miles, is classified into two geographic regions: northern and southern. The northern region, including Bomi, Bong, Gbarpolu, Grand Cape Mount, Lofa, Margibi, and Montserrado, has a total population of 2,233,091. The southern region, including Grand Bassa, Grand Gedeh, Grand Kru, Maryland, Nimba, River Cess, River Gee, and Senoe, has a population of 1,243,517. Finally, Sierra Leone, with a population of 4,976,871 and an area of 27,699 square miles, is categorized as its own region. A.3. Epidemiological data Table 3 presents the disease transition parameters along with their symbols, description, case study values, and reference resources that are used to estimate data values. Case fatality proportion without treatment: In this paper, we define case fatality proportion without treatment as the percentage of Ebola-related deaths among cases that do not receive medical treatment from an ETC facility over the course of the disease. The case fatality proportion from past Ebola outbreaks is known to vary between 25 and 90% (WHO, 2016b). We estimate the fatality proportion for Guinea, Liberia, and Sierra Leone to be 40%, 39%, and 31%, respectively. Case fatality proportion with treatment: Similarly, the case fatality proportion with treatment is defined as the percentage of Ebola-related deaths among cases who receive medical treatment from an ETC facility over the course of the disease. Based on WHO (2016b), the case fatality proportion

1061

with treatment is estimated to be 0.35, 0.29, and 0.24 in Guinea, Liberia and Sierra Leone, respectively. Average time from onset of symptoms to funeral: Based on the work of the WHO Ebola Response Team (2014), we estimated the average time from the onset of symptoms to funeral as 0.93 weeks (6.4 ± 5.3 days), 2.22 weeks (7.9 ± 8.0 days), and 2.50 weeks (8.6 ± 6.9 days) in Guinea, Liberia, and Sierra Leone, respectively. Average time from onset of symptoms to recovery: In addition, this study defines the mean length of recovery as the length of infection. Based on information from the WHO Ebola Response Team (2014), we estimated the fitted values for the length of infection for infected individuals as 2.50 (16.3 ± 6.1 days), 2.63 (15.4 ± 8.2 days), and 2.85 (17.2 ± 6.2 days) weeks in Guinea, Liberia, and Sierra Leone, respectively. Average time from treatment to funeral (death): Based on data from the WHO Ebola Response Team (2014), we estimate the average time from treatment to death for infected individuals as 1.0 week in Guinea. Furthermore, we assume 2.27 and 2.50 weeks as the average time from treatment to death for infected individuals in Liberia and Sierra Leone, respectively. Average time from treatment to recovery (in weeks) Furthermore, the average time it takes an infected individual receiving medical care to recover is estimated to be 1.56 weeks (11 ± 5.4 days) in Guinea, 2.27 weeks (12.8 ± 8.1 days) in Liberia, and 2.37 weeks (12.4 ± 5.8 days) in Sierra Leone (WHO Ebola Response Team, 2014). Rate of safe burial: We assume that it takes approximately 4.5 ± 2 days for an Ebola-related dead body to be safely buried (Rivers et al., 2014). Therefore, we estimate a safe burial rate of 0.73, 0.71, and 0.74 for Guinea, Liberia, and Sierra Leone, respectively. Transmission rate in the community: Community transmission rate is estimated based on the study of (Camacho et al., 2014) and Rivers et al. (2014), considering both person-toperson contact during community interaction and fraction of the individuals who are infected during the hospitalization period (0.198 − 0.397 at 95% confidence interval). Therefore, transmission rate in the community is estimated to be 0.33 in Guinea, 0.32 in Sierra Leone, and 0.28 in Liberia. Transmission rate at traditional funerals: Transmission rates via person-to-person contact during traditional funerals are found to change within the interval 0.08–2.00 because cultural practices pertaining to burial might vary among African countries (Camacho et al., 2014). Therefore, we assume the best estimates for transmission rate during funerals are 0.68 in Guinea and Liberia, and 0.73 in Sierra Leone. B. Linearization Here we linearize the logical constraint that describes the number of hospitalized individuals in Eq. (15) (Kıbıs¸ & Büyüktahtakın, 2017). The logical term in Eq. (15) can be rewritten as







I¯(j k,l ) = C (j k,l ) − T j(k,l ) z(jk,l ) + I (j k,l ) 1 − z(jk,l )



j ∈ J\J¯, k∈K, l ∈L, (B-18a)

(k,l )

where z j represents a binary variable, which takes a value of 1 when the number of infected individuals to be hospitalized is restricted by the number of available beds in the ETCs, and zero when all infected individuals are hospitalized for treatment in the ETCs due to available capacity. In order to imply that I¯(j k,l ) chooses the minimum0 value between (C (j k,l ) − T j(k,l ) ) and I (j k,l ) in Eq. (B-18a), the following constraints are added to the model:

I¯(j k,l ) ≤ C (j k,l ) − T j(k,l )

j ∈ J\J¯, k ∈ K, l ∈ L,

(B-18b)

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E.

1062

References

I¯(j k,l ) ≤ I (j k,l )

j ∈ J\J¯, k ∈ K, l ∈ L.

(B-18c)

Eqs. (B-18a) that depending on the minimum value  –(B-18c) ensure  between C (j k,l ) − T j(k,l )

and I (j k,l ) , the binary variable z(jk,l ) will re-

ceive a value of 0 or 1. Note that Equation (B-18a) is still not linear due to quadratic terms, both of which are the product of continuous and binary variables. Therefore, Equation (B-18a) is linearized as follows: Let U j(k,l ) be an auxiliary variable, and substitute the new vari-





able U j(k,l ) with C (j k,l ) − T j(k,l ) z(jk,l ) . Furthermore, let HLB ≤ C (j k,l ) − T j(k,l ) ≤ HUB , where HLB and HUB are lower and upper bounds for

C (j k,l ) − T j(k,l ) , respectively. In order to ensure that U j(k,l ) is equal



to



C (j k,l ) − T j(k,l ) z(jk,l ) , the following constraints are added to the

model:

U j(k,l ) ≤ HUB z(jk,l ) ,

j ∈ J\J¯, k ∈ K, l ∈ L,

(B-18d)

U j(k,l ) ≥ HLB z(jk,l ) ,

j ∈ J\J¯, k ∈ K, l ∈ L,

(B-18e)









U j(k,l ) ≤ C (j k,l ) − T j(k,l ) − HLB 1 − z(jk,l ) ,

j ∈ J\J¯, k ∈ K, l ∈ L, (B-18f)









U j(k,l ) ≥ C (j k,l ) − T j(k,l ) − HUB 1 − z(jk,l ) ,

j ∈ J\J¯, k ∈ K, l ∈ L. (B-18g)

Examining successively the two possible values of the binary variable z(jk,l ) , we observe that Eqs. (B-18d)–(B-18g) ensure that the

value of U j(k,l ) is equal to C (j k,l ) − T j(k,l ) . A similar procedure is followed in order to linearize the second quadratic term I (j k,l ) (1 − z(jk,l ) ) in Eq. (B-18a). Let W j(k,l ) be a continuous variable, and replace it with the quadratic term I (j k,l ) (1 −

z(jk,l ) ). In addition, let ILB ≤ I (j k,l ) ≤ IUB , where ILB and IUB represent

the lower and upper bounds for I (j k,l ) , respectively. Therefore, the following four constraints are added to the model to ensure that the value of W j(k,l ) is equal to I (j k,l ) (1 − z(jk,l ) ):









W j(k,l ) ≤ IUB 1 − z(jk,l ) , W j(k,l ) ≥ ILB 1 − z(jk,l ) ,

j ∈ J\J¯, k ∈ K, l ∈ L,

(B-18h)

j ∈ J\J¯, k ∈ K, l ∈ L,

(B-18i)

W j(k,l ) ≤ I (j k,l ) − ILB z(jk,l ) ,

j ∈ J\J¯, k ∈ K, l ∈ L,

(B-18j)

W j(k,l ) ≥ I (j k,l ) − IUB z(jk,l ) ,

j ∈ J\J¯, k ∈ K, l ∈ L.

(B-18k)

Thus, Eq. (B-18a) becomes equivalent to the linear constraint

I¯(j k,l ) = U j(k,l ) + W j(k,l ) ,

j ∈ J\J¯, k ∈ K, l ∈ L,

(B-18l)

where U j(k,l ) and W j(k,l ) are restricted by Eqs. (B-18d)–(B-18k). We obtain an MIP model by replacing the non-linear Eq. (15) with equivalent Constraints (B-18b)–(B-18l). Supplementary material Supplementary material associated with this article can be found, in the online version, at 10.1016/j.ejor.2017.08.037.

Ajelli, M., Merler, S., Fumanelli, L., y Piontti, A. P., Dean, N. E., Longini, I. M., et al. (2016). Spatiotemporal dynamics of the Ebola epidemic in Guinea and implications for vaccination and disease elimination: A computational modeling analysis. BMC Medicine, 14(1), 130. Allen, L. J. (1994). Some discrete-time SI, SIR, and SIS epidemic models. Mathematical Biosciences, 124(1), 83–105. Allen, L. J., Brauer, F., Van den Driessche, P., & Wu, J. (2008). Mathematical epidemiology: 1945. Springer. Altay, N., & Green, W. G. (2006). OR/MS research in disaster operations management. European Journal of Operational Research, 175(1), 475–493. Althaus, C. L., Gsteiger, S., & Low, N. (2014). Ebola virus disease outbreak in Nigeria: Lessons to learn. (Accessed January 12, 2015). https://peerj.com/preprints/569v1. pdf. Ancel, M. L., Newman, M., Martin, M., & Schrag, S. (2003). Applying network theory to epidemics: Control measures for mycoplasma pneumoniae outbreaks. Emerging Infectious Diseases, 9(2), 204–210. Anderson, R. M., & May, R. M. (1991). Infectious diseases of humans: 1. Oxford: Oxford University Press. Arinaminpathy, N., & McLean, A. (2009). Logistics of control for an influenza pandemic. Epidemics, 1(2), 83–88. Babad, H., Nokes, D., Gay, N., Miller, E., Morgan-Capner, P., & Anderson, R. (1995). Predicting the impact of measles vaccination in England and Wales: Model validation and analysis of policy options. Epidemiology and Infection, 114(2), 319–344. Bailey, N. T. (1975). The mathematical theory of infectious diseases and its applications. New York: Hafner Press. Bay, M. (2003). GeoHive global statistics. Reference Reviews, 17(6), 26–27. Berman, O., & Gavious, A. (2007). Location of terror response facilities: A game between state and terrorist. European Journal of Operational Research, 177(2), 1113–1133. Brauer, F., Feng, Z., & Castillo-Chavez, C. (2010). Discrete epidemic models. Mathematical Biosciences, 7(1), 1–33. Büyüktahtakın, I.˙ E., Kıbıs¸ , E. Y., Cobuloglu, H. I., Houseman, G. R., & Lampe, J. T. (2015). An age-structured bio-economic model of invasive species management: insights and strategies for optimal control. Biological invasions, 17(9), 2545–2563. Camacho, A., Kucharski, A., Funk, S., Breman, J., Piot, P., & Edmunds, W. (2014). Potential for large outbreaks of Ebola virus disease. Epidemics, 9, 70–78. Centers for Disease Control and Prevention (CDC) (2015). 2014 Ebola outbreak in West Africa – Case counts. (Accessed October 12, 2015) http://www.cdc.gov/vhf/ ebola/outbreaks/2014- west- africa/case- counts.html. Chen, T., Leung, R. K.-K., Liu, R., Chen, F., Zhang, X., Zhao, J., & Chen, S. (2014). Risk of imported Ebola virus disease in China. Travel Medicine and Infectious Disease, 12(6), 650–658. Chowell, G., Hengartner, N. W., Castillo-Chavez, C., Fenimore, P. W., & Hyman, J. (2004). The basic reproductive number of Ebola and the effects of public health measures: The cases of Congo and Uganda. Journal of Theoretical Biology, 229(1), 119–126. Chowell, G., & Nishiura, H. (2014). Transmission dynamics and control of Ebola virus disease (EVD): A review. BMC Medicine, 12(1), 196. Cohen, M. L. (20 0 0). Changing patterns of infectious disease. Nature, 406(6797), 762–767. Cohen, N. J. (2016). Travel and border health measures to prevent the international spread of Ebola. MMWR Supplements, 65. Craft, D. L., Wein, L. M., & Wilkins, A. H. (2005). Analyzing bioterror response logistics: The case of anthrax. Management Science, 51(5), 679–694. Dasaklis, T. K., Pappis, C. P., & Rachaniotis, N. P. (2012). Epidemics control and logistics operations: A review. International Journal of Production Economics, 139(2), 393–410. Dasaklis, T. K., Rachaniotis, N., & Pappis, C. (2017). Emergency supply chain management for controlling a smallpox outbreak: The case for regional mass vaccination. International Journal of Systems Science: Operations & Logistics, 4(1), 27– 40. Dimitrov, N. B., & Meyers, L. A. (2010). Mathematical approaches to infectious disease prediction and control. INFORMS TutORials in Operations Research, 7, 1– 25. Eubank, S., Guclu, H., Kumar, V. A., Marathe, M. V., Srinivasan, A., Toroczkai, Z., & Wang, N. (2004). Modelling disease outbreaks in realistic urban social networks. Nature, 429(6988), 180–184. Fasina, F., Shittu, A., Lazarus, D., Tomori, O., Simonsen, L., Viboud, C., & Chowell, G. (2014). Transmission dynamics and control of Ebola virus disease outbreak in Nigeria, July to September 2014. European Surveillance, 19(40), 20920. Ferguson, N. M., Keeling, M. J., Edmunds, W. J., Gani, R., Grenfell, B. T., Anderson, R. M., & Leach, S. (2003). Planning for smallpox outbreaks. Nature, 425(6959), 681–685. Frauenthal, J. C. (2012). Mathematical modeling in epidemiology. Springer Science & Business Media. Grépin, K. A. (2015). International donations to the Ebola virus outbreak: Too little, too late? BMJ, 350, h376. Griggs, M. B. (2014). What happens after someone survives Ebola?(Accessed March, 03, 2015). http://www.popsci.com/article/science/what-happensafter- someone- survives- ebola. Heffernan, J., Smith, R., & Wahl, L. (2005). Perspectives on the basic reproductive ratio. Journal of the Royal Society Interface, 2(4), 281–293.

˙ Büyüktahtakın et al. / European Journal of Operational Research 265 (2018) 1046–1063 I.E. Jia, H., Ordóñez, F., & Dessouky, M. (2007a). A modeling framework for facility location of medical services for large-scale emergencies. IIE Transactions, 39(1), 41–55. Jia, H., Ordóñez, F., & Dessouky, M. M. (2007b). Solution approaches for facility location of medical supplies for large-scale emergencies. Computers & Industrial Engineering, 52(2), 257–276. Kaplan, E. H., Craft, D. L., & Wein, L. M. (2003). Analyzing bioterror response logistics: The case of smallpox. Mathematical Biosciences, 185(1), 33–72. Kermack, W. O., & McKendrick, A. G. (1932). Contributions to the mathematical theory of epidemics. II. The problem of endemicity. Proceedings of the Royal Society of London. Series A, 138(834), 55–83. Kıbıs¸ , E. Y., & Büyüktahtakın, I.˙ E. (2017). Optimizing invasive species management: A mixed-integer linear programming approach. European Journal of Operational Research, 259(1), 308–321. Kurahashi, S., & Terano, T. (2015). A health policy simulation model of smallpox and Ebola haemorrhagic fever. In Agent and multi-agent systems: Technologies and applications (pp. 405–415). Springer. Lasry, A., Zaric, G. S., & Carter, M. W. (2007). Multi-level resource allocation for HIV prevention: A model for developing countries. European Journal of Operational Research, 180(2), 786–799. Lee, E. K., Pietz, F., Benecke, B., Mason, J., & Burel, G. (2013). Advancing public health and medical preparedness with operations research. Interfaces, 43(1), 79–98. Legrand, J., Grais, R., Boelle, P., Valleron, A., & Flahault, A. (2007). Understanding the dynamics of Ebola epidemics. Epidemiology and Infection, 135(04), 610–621. Liu, M., & Zhao, L. (2012). An integrated and dynamic optimisation model for the multi-level emergency logistics network in anti-bioterrorism system. International Journal of Systems Science, 43(8), 1464–1478. Long, E. F., Vaidya, N. K., & Brandeau, M. L. (2008). Controlling co-epidemics: Analysis of HIV and tuberculosis infection dynamics. Operations Research, 56(6), 1366–1381. Longini, I. M., Halloran, M. E., Nizam, A., & Yang, Y. (2004). Containing pandemic influenza with antiviral agents. American Journal of Epidemiology, 159(7), 623–633. Longini, I. M., Halloran, M. E., Nizam, A., Yang, Y., Xu, S., Burke, D. S., et al. (2007). Containing a large bioterrorist smallpox attack: A computer simulation approach. International Journal of Infectious Diseases, 11(2), 98–108. Meltzer, M. I., Atkins, C. Y., Santibanez, S., Knust, B., Petersen, B. W., Ervin, E. D., et al. (2014). Estimating the future number of cases in the Ebola epidemic: Liberia and Sierra Leone, 2014–2015. Morbidity and Mortality Weekly Response Surveillance Summaries, 63(Supplement 3), 1–14. Moon, S., Sridhar, D., Pate, M. A., Jha, A. K., Clinton, C., Delaunay, S., et al. (2015). Will Ebola change the game? Ten essential reforms before the next pandemic. The report of the Harvard-LSHTM Independent Panel on the Global Response to Ebola. The Lancet, 386(10 0 09), 2204–2221. Nielsen, C. F., Kidd, S., Sillah, A., Davis, E., Mermin, J., Kilmarx, P. H., et al. (2015). Improving burial practices and cemetery management during an Ebola virus disease epidemic–Sierra Leone, 2014. Morbidity and Mortality Weekly Report, 64(1), 20–27. Pandey, A., Atkins, K. E., Medlock, J., Wenzel, N., Townsend, J. P., Childs, J. E., et al. (2014). Strategies for containing Ebola in West Africa. Science, 346(6212), 991–995. Patel, R., Longini, I. M., & Halloran, M. E. (2005). Finding optimal vaccination strategies for pandemic influenza using genetic algorithms. Journal of Theoretical Biology, 234(2), 201–212. Porco, T. C., Holbrook, K. A., Fernyak, S. E., Portnoy, D. L., Reiter, R., & Aragón, T. J. (2004). Logistics of community smallpox control through contact tracing and ring vaccination: a stochastic network model. BMC Public Health, 4(1), 34. Rachaniotis, N. P., Dasaklis, T. K., & Pappis, C. P. (2012). A deterministic resource scheduling model in epidemic control: A case study. European Journal of Operational Research, 216(1), 225–231. Radio France Internationale (2014). Liberia taking difficult decisions to limit movement in Ebola areas. (Accessed April 1, 2017). http://en.rfi.fr/africa/ 20140805- liberia- taking- difficult- decisions- limit- movement- ebola- areas. Ren, Y., Ordóñez, F., & Wu, S. (2013). Optimal resource allocation response to a smallpox outbreak. Computers & Industrial Engineering, 66(2), 325–337. Riley, S., & Ferguson, N. M. (2006). Smallpox transmission and control: Spatial dynamics in great britain. Proceedings of the National Academy of Sciences, 103(33), 12637–12642. Rivers, C. M., Lofgren, E. T., Marathe, M., Eubank, S., & Lewis, B. L. (2014). Modeling the impact of interventions on an epidemic of Ebola in Sierra Leone and Liberia. PLoS Currents, 6, 1–16.

1063

Siettos, C., Anastassopoulou, C., Russo, L., Grigoras, C., & Mylonakis, E. (2011). Modeling the 2014 Ebola virus epidemic–agent-based simulations, temporal analysis and futurepredictions for Liberia and Sierra Leone. PLoS Currents, 7, 1–22. Siettos, C. I., Anastassopoulou, C., Russo, L., Grigoras, C., & Mylonakis, E. (2016). Forecasting and control policy assessment for the Ebola virus disease (EVD) epidemic in Sierra Leone using small-world networked model simulations. BMJ Open, 6(1), e008649. Soto, S. (2009). Human migration and infectious diseases. Clinical Microbiology and Infection, 15(s1), 26–28. Tanner, M. W., Sattenspiel, L., & Ntaimo, L. (2008). Finding optimal vaccination strategies under parameter uncertainty using stochastic programming. Mathematical Biosciences, 215(2), 144–151. Team, W. E. R. (2016). After Ebola in West Africa: Unpredictable risks, preventable epidemics. The New England Journal of Medicine, 2016(375), 587–596. Tebbens, R. J. D., & Thompson, K. M. (2009). Priority shifting and the dynamics of managing eradicable infectious diseases. Management Science, 55(4), 650– 663. UNESCO (2014). UNESCO’s response to Ebola. (Accessed January, 09, 2016). http:// unesdoc.unesco.org/images/0 023/0 02311/231158e.pdf. United Nations (2014a). Backgrounder for ECOSOC special session on Ebola. (Accessed November 15, 2014). http://www.un.org. United Nations (2014b). With spread of Ebola outpacing response, security council adopts resolution 2177 (2014) urging immediate action, end to isolation of affected states. (Accessed April 1, 2017). https://www.un.org/press/en/2014/ sc11566.doc.htm. United Nations (2015a). Resources for Results I, Office of the UN Special Envoy on Ebola. (Accessed October 17, 2015). https://ebolaresponse.un.org/sites/default/ files/20141124_rri_0.pdf. United Nations (2015b). Resources for Results IV, Office of the UN Special Envoy on Ebola. (Accessed October 17, 2015). http://ebolaresponse.un.org/sites/default/ files/resources_for_results_iv.pdf. Wang, W., & Zhao, X.-Q. (2004). An epidemic model in a patchy environment. Mathematical Biosciences, 190(1), 97–112. Wells, C., Yamin, D., Ndeffo-Mbah, M. L., Wenzel, N., Gaffney, S. G., Townsend, J. P., et al. (2015). Harnessing case isolation and ring vaccination to control Ebola. PLoS Neglected Tropical Diseases, 9(5), e0 0 03794. WHO (2016a). Ebola in West Africa: Heading for catastrophe?(Accessed March, 30, 2017). http://www.who.int/csr/disease/ebola/ebola- 6- months/west- africa/en/. WHO (2016b). Ebola virus disease. (Accessed May, 09, 2016). http://www.who.int/ mediacentre/factsheets/fs103/en/. WHO (2016c). WHO leadership statement on the Ebola response and WHO reforms. (Accessed March, 30, 2017). http://www.who.int/csr/disease/ebola/ joint- statement- ebola/en/. WHO Ebola Response Team (2014). Ebola virus disease in West Africa: The first 9 months of the epidemic and forward projections. The New England Journal of Medicine, 371(16), 1481–1495. World Health Organization (2014a). Ebola response roadmap. (Accessed September 28, 2014). http://www.who.int/csr/resources/publications/ebola/ response-roadmap/en/. World Health Organization (2014b). Ebola response roadmap situation reports. (Accessed December 10, 2014). http://apps.who.int/ebola/en/current-situation/ ebola- situation- report. World Health Organization (2014c). Ebola virus disease. (Accessed June 11, 2015). http://www.who.int/mediacentre/factsheets/fs103/en/. World Health Organization (2015). Factors that contributed to undetected spread of the Ebola virus and impeded rapid containment. (Accessed January, 09, 2016). http://www.who.int/csr/disease/ebola/one- year- report/factors/en/. World Health Organization, Regional Office for Africa (2014). Disease outbreak news. (Accessed October 20, 2015). http://www.afro.who.int/en/ clusters- a- programmes/dpc/epidemic- a- pandemic- alert- and- response/ outbreak-news.html. Yarmand, H., Ivy, J. S., Denton, B., & Lloyd, A. L. (2014). Optimal two-phase vaccine allocation to geographically different regions under uncertainty. European Journal of Operational Research, 233(1), 208–219. Zaman, G., Kang, Y. H., & Jung, I. H. (2009). Optimal treatment of an SIR epidemic model with time delay. BioSystems, 98(1), 43–50. Zaric, G. S., Brandeau, M. L., & Barnett, P. G. (20 0 0). Methadone maintenance and HIV prevention: A cost-effectiveness analysis. Management Science, 46(8), 1013–1031.