Mathematical Modelling, Simulation, and Optimal Control of the 2014 Ebola Outbreak in West Africa

The Ebola virus is currently one of the most virulent pathogens for humans. The latest major outbreak occurred in Guinea, Sierra Leonne and Liberia in 2014. With the aim of understanding the spread of infection in the affected countries, it is crucial to modelize the virus and simulate it. In this paper, we begin by studying a simple mathematical model that describes the 2014 Ebola outbreak in Liberia. Then, we use numerical simulations and available data provided by the World Health Organization to validate the obtained mathematical model. Moreover, we develop a new mathematical model including vaccination of individuals. We discuss different cases of vaccination in order to predict the effect of vaccination on the infected individuals over time. Finally, we apply optimal control to study the impact of vaccination on the spread of the Ebola virus. The optimal control problem is solved numerically by using a direct multi-shooting method.


Introduction
Ebola is a lethal virus for humans. It was previously confined to Central Africa but recently was also identified in West Africa [1]. As of October 8, 2014, the World Health Organization (WHO) reported 4656 cases of Ebola virus deaths, with most cases occurring in Liberia [2]. The extremely rapid increase of the disease and the high mortality rate make this virus a major problem for public health. Typically, patients present fever, malaise, abdominal, headache, and asthenia. One week after the onset of symptoms, a rash often appears followed by haemorrhagic complications, leading to death after an average of 10 days in 50%-90% of infections [3][4][5][6]. Ebola is transmitted through direct contact with blood, bodily secretions and tissues of infected ill or dead humans and non-human primates [7][8][9].
The main aim of this work is to understand the dynamics of the Liberian population infected by Ebola virus in 2014, by using an appropriate mathematical model. More precisely, we begin by considering a system of ordinary differential equations to describe the 2014 Ebola outbreak, which is nothing else than an epidemic SIR model, that is, a model based on the division of the population into three groups: the Susceptible, the Infected, and the Recovered [10][11][12]. We simulate the model obtained by parameters estimated on November 4, 2014, by Kaurov [13]. Our objectives are to better understand the outbreak and to predict the effect of vaccination on the infected individuals over time.
In order to deal with the epidemic of Ebola, governments have decided to implement tough measures. For example, some have applied quarantine procedures while others have opted for mass vaccination plans as a precaution to the epidemic [14][15][16][17]. In this context, we consider an optimal control problem to study the effect of vaccination on the spread of virus.
The text is organized as follows. In Section 2 we introduce a basic mathematical model to describe the dynamics of the Ebola virus during the 2014 outbreak in Liberia. After the mathematical modelling, we use the obtained model to simulate it in Section 3 with the parameters estimated from recent statistical data based on the WHO report of the 2014 Ebola outbreak [18]. Following [13], in Section 4 we study different cases of vaccination of individuals by adding a vaccination term to the model. Section 5 presents an optimal control problem, which we use to study the impact of a vaccination campaign on the spread of the virus. We end with Section 6 of conclusions and future work, where the results are summarized and some research perspectives presented.

The Basic Mathematical Model
The objective of this section is to describe mathematically the dynamics of the population infected by the Ebola virus. The dynamics is described by a system of differential equations. This system is based on the common SIR (Susceptible-Infectious-Recovery) epidemic model, where the population is divided into three groups: the susceptible group, denoted by S(t), the infected group, denoted by I(t), and the recovered group, denoted by R(t). The total population, assumed constant during the short period of time under study, is given by N = S(t) + I(t) + R(t).
To create the equations that describe the population of each group along time, let us start with the fact that the population of the susceptible group will be reduced as the infected come into contact with them with a rate of infection β. This means that the change in the population of susceptible is equal to the negative product of β with S(t) and I(t): Now, let us create the equation that describes the infected group over time, knowing that the population of this group changes in two ways: (i) people leave the susceptible group and join the infected group, thus adding to the total population of infected a term βS(t)I(t); (ii) people leave the infected group and join the recovered group, reducing the infected population by −µI(t).
This is written as Finally, the equation that describes the recovery population is based on the individuals recovered from the virus at rate µ. This means that the recovery group is increased by µ multiplied by I(t): (3) Figure 1 shows the relationship between the three variables of our SIR model. As shown in the next section, this simple model describes well the 2014 Ebola outbreak in Liberia for suitable chosen parameters β and µ and appropriate initial conditions S(0), I(0) and R(0).

Numerical Simulations of the Basic Model and Discussion
With the aim to approximate the 2014 Ebola outbreak in Liberia, we now simulate the complete set of equations (1)-(3) that describe the SIR model: where β is the rate of infection and µ the rate of recovery. The data for the outbreak in Liberia is based on the cumulative numbers of total reported cases (confirmed and probable) taken from the official data reported by the Ministry of Health of Liberia up to October 26, 2014, as provided by the World Health Organisation (WHO) [18] (in particular, the Ebola virus disease cases reported each week from Liberia are given in Figure 2 of [19]). The model fits well the reported data of confirmed (means infected) cases in Liberia. Note that the exact population size N does not need to be known to estimate the model parameters as long as the number of cases is small compared to the total population size [20]. To estimate the parameters of the model, we adapted the initialisation of S and I with the reported data of [19] by fitting the actual data of confirmed cases in Liberia. The data that we are studying is the number of individuals with confirmed infection between July 14 and October 12, 2014. The result of fitting is shown in Figure 2. The comparison between the curve of infected obtained by our simulation and the reported data of confirmed cases by WHO shows that the mathematical model (4) fits well the real data by using β = 0.000318 as the rate of infection and µ = 0.0175 as the rate of recovery. The initial numbers of susceptible, infected, and recovered groups, are given by S(0) = 460, I(0) = 12, and R(0) = 0, respectively. The choice of S(0) is in agreement with the data shown in Figure 2 of [19], as the biggest value in the histogram is 460. To interpret the simulation results, we describe the evolution of each group of individuals over time and the link between them. Figure 3 shows that the susceptible group begins to plummet and, simultaneously, the number of infected begins to rise. This is due to how infectious the Ebola virus is. When we compare our results with those of Astacio et al. [21], when they modeled the 1976 outbreak of Ebola in Yambuku, a small village in Mongala District in northern Democratic Republic of Congo (previously Zaire), and the 1995 outbreak of Ebola in Kikwit, in the southwestern part of the Democratic Republic of Congo, we note that the evolution of each group obtained in our simulations is reasonable. The main difference between the    2014 outbreak of Ebola in Liberia and the previous outbreaks of Ebola in Yambuku and Kikwit studied in [21] is at the end of the evolution of infected: while in the previous 1976 and 1995 outbreaks of Ebola in Yambuku and Kikwit the infected group goes to zero, now in the 2014 outbreak of Ebola in Liberia the infected group remains important with a high number of infected individuals by the virus, which describes the reality of what is currently happening in Liberia. One can conclude that the 2014 Ebola outbreak of Liberia is much more dangerous than the previous Ebola outbreaks.

Modelling and Simulation of Ebola Outbreak with Vaccination
In this section, we simulate the model by using parameters estimated on November 4, 2014, by Kaurov [13], who studied statistically the recent data of the World Health Organisation (WHO) [18]. In his statistical study, Kaurov studied the outbreak by modeling it with Wolfram's Mathematica language. Here we use his SIR model for modelling and time-discrete equations for resolution of the set of equations through Matlab.

The Basic Model with Kaurov's Estimation of Parameters
Let us recall the complete set of equations considered by [13]: where, as before, β is the rate of infection and µ is the rate of recovery. The parameters obtained by Kaurov, knowing that 95% of population is susceptible and 5% of population is infected, are β = 0.2 and µ = 0.1 [13]. In agreement, the initial susceptible, infected, and recovered populations are given, respectively, by S(0) = 0.95, I(0) = 0.05, and R(0) = 0 (note that in contrast with Sections 2 and 3, the variables are now dimensionless). Our numerical simulations were done with the ODE solver of Matlab. The results are shown in Figure 4. To interpret the results, we describe the evolution of each group over time and the links between them. Figure 4 shows that the susceptible group begins to plummet due to how infectious the virus is and, simultaneously, the infected group's number begins to rise. An important dimensionless quantity that plays a key role in any SIR model is the well-known basic reproduction number R 0 [22]. This number tells us how fast the disease will spread at the epidemic. For model (5), the basic reproduction number is simply given by β/µ. Let us now discuss the comparison between the different values of R 0 obtained recently. On September 2, 2014, Althaus studied the case of the 2014 Ebola outbreak occurred in Guinea, in Sierra Leone, and in Liberia [20]. In his study, Althaus used the WHO data occurred in the three countries. He found different values for the basic reproduction number, depending on the country: in Guinea R 0 = 1.51, in Sierra Leone the value of R 0 is given by 2.53, while in Liberia R 0 is equal to 1.59. Let us now compare between these values and the value for the basic reproduction number studied by Kaurov, knowing that his study, which was done on November 4, 2014, is more recent than the one of Althaus, done September 2, 2014. As the parameters estimated statistically by Kaurov are for the three countries, our comparison will be between R 0,A , the mean of the three basic reproduction numbers found by Althaus, and R 0,K , the number found more recently by Kaurov. The value of R 0,A is given by 1.876, where the more recent R 0,K is given by 2. We claim that the basic reproduction number increased between September and November 2014.

Dynamics in Case of Vaccination
Infectious diseases have tremendous influence on human life. In last decades, controlling infectious diseases has been an increasingly complex issue. A strategy to control infectious diseases is through vaccination. Now our idea is to study the effect of vaccination in practical Ebola situations. The case of the French nurse cured of Ebola with the help of an experimental vacine is a proof of the possibility of treatment by vaccinating infected individuals [23].
In this section, we improve the SIR model described in Section 4.1. The system of equations that describe the Susceptible-Infectious-Recovery (SIR) model with vaccination is: where β is the rate of infection, µ is the rate of recovery, and υ is the percentage of individuals vaccinated every day [24]. Figure 5 shows the relationship between the variables of the SIR model with vaccination (6).

Numerical Simulation in Case of Vaccination
We simulate the SIR model with vaccination (6)      to zero in about 90 days (see Figure 6), where in the presence of vaccination the curve of the infected group goes to zero in about 45 days (see Figure 7). This shows the efficiency of using vaccination in controlling Ebola. Table 1 presents the maximum percentage of infected individuals, corresponding to different rates of vaccinations, and the respective number of days for the peak of infected.

The Optimal Control Problem
Optimal control has been recently used with success in a series of epidemiological problems: see [17,[25][26][27] and references therein. In this section, we present an optimal control problem by introducing into the model (5) a control u(t) representing the vaccination rate at time t. The control u(t) is the fraction of susceptible individuals being vaccinated per unit of time. We assume that all susceptible individuals that are vaccinated are transferred directly to the removed class. Then, the mathematical model with control is given by the following system of nonlinear differential equations: where S(0) ≥ 0, I(0) ≥ 0, and R(0) ≥ 0 are given. Our goal is to reduce the infected individuals and the cost of vaccination. Precisely, our optimal control problem consists of minimizing the objective functional subject to 0 ≤ u(t) ≤ 0.9 for t ∈ [0, t f ], where the parameters A ≥ 0 and t f denote, respectively, the weight on cost and the duration of the vaccination program. We present the numerical solution to this optimal control problem, comparing the results with the numerical simulation results obtained in Section 4 via Matlab. As in Section 4, the rate of infection is given by β = 0.2, the recovered rate is given by µ = 0.1, and we use S(0) = 0.95, I(0) = 0.05, and R(0) = 0 for the initial number of susceptible, infected, and recovered populations, respectively (i.e., at the beginning 95% of population is susceptible and 5% is already infected with Ebola). Note that when u(t) ≡ 0 model (7) reduces to (5); and when u(t) ≡ υ model (7) reduces to (6).
The numerical simulation of the system without control (u(t) ≡ 0) was done with Matlab in Section 4.1. For our numerical solution of the optimal control problem, we have used the ACADO solver [28], which is based on a multiple shooting method, including automatic differentiation and based ultimately on the semidirect multiple shooting algorithm of Bock and Pitt [29]. ACADO is a self-contained public domain software environment written in C++ for automatic control and dynamic optimization. Figures 8 and 9 show the significant difference in the number of infected individuals with and without control. In fact, in case of optimal control, the function I decreases greatly: the maximum value on the infected curve I under optimal control is 5.6% against 17.9% without any control (see Figure 9).
In Figure 10, we see that the number of susceptible S decreases more rapidly during the vaccination campaign. It reaches 11% at the end of the      campaign against 19% in the absence of optimal control. Figure 11 shows that the number of recovered individuals increases rapidly. The number R(t f ) at the end of the optimal control vaccination period is 88.7% instead of 80.5% without control. Figure 12 gives a representation of the optimal control u(t).

Conclusion
We investigated a mathematical model that provides a good description of the 2014 Ebola outbreak in Liberia: the model fits well the data of confirmed cases in Liberia provided by the World Health Organization. The evolution of infected individuals shows that the infected population remains important, describing the actual evolution of the virus in Liberia. This evolution is different than Ebola's outbreak in the last decades. Our study of the SIR model with vaccination shows that vaccination is a very efficient factor in reducing the number of infected individuals in a short period of time and increasing the number of recovered individuals. We also discussed an optimal control problem related to the impact of a vaccination strategy on the spread of the 2014 outbreak of Ebola in West Africa. The numerical simulation of both systems, with and without control, shows that an optimal control strategy greatly helps to reduce the number of infected and susceptible individuals and increases the number of recovered individuals.
As future work, we plan to include in our study other factors. In particular, we intend to include in the mathematical model the treatment of infected individuals with a quarantine procedure. Another interesting line of research is to investigate the effect of impulsive vaccination.