EBOLA MODEL AND OPTIMAL CONTROL WITH VACCINATION CONSTRAINTS

. The Ebola virus disease is a severe viral haemorrhagic fever syn-drome caused by Ebola virus. This disease is transmitted by direct contact with the body ﬂuids of an infected person and objects contaminated with virus or infected animals, with a death rate close to 90% in humans. Recently, some mathematical models have been presented to analyse the spread of the 2014 Ebola outbreak in West Africa. In this paper, we introduce vaccination of the susceptible population with the aim of controlling the spread of the disease and analyse two optimal control problems related with the transmission of Ebola disease with vaccination. Firstly, we consider the case where the total number of available vaccines in a ﬁxed period of time is limited. Secondly, we analyse the situation where there is a limited supply of vaccines at each instant of time for a ﬁxed interval of time. The optimal control problems have been solved analytically. Finally, we have performed a number of numerical simulations in order to compare the models with vaccination and the model without vacci- nation, which has recently been shown to ﬁt the real data. Three vaccination scenarios have been considered for our numerical simulations, namely: unlim- ited supply of vaccines; limited total number of vaccines; and limited supply of vaccines at each instant of time.


1.
Introduction. Ebola is a lethal virus for humans that is currently under strong research due to the recent outbreak in West Africa and its socioeconomic impact (see, e.g., [5,17,18,21,24,25,30,34,37] and references therein). World Health Organization (WHO) has declared Ebola virus disease epidemic as a public health emergency of international concern with severe global economic burden. At fatal Ebola infection stage, patients usually die before the antibody response. Mainly after the 2014 Ebola outbreak in West Africa, some attempts to obtain a vaccine for Ebola disease have been realized. According to the WHO, results in July 2015 from an interim analysis of the Guinea Phase III efficacy vaccine trial show that VSV-EBOV (Merck, Sharp & Dohme) is highly effective against Ebola [2].
Since 2014, different mathematical models to analyze the spread of the 2014 Ebola outbreak have been presented (see, e.g., [3,4,19,32,33] and references therein). In these models the populations under study are divided into compartments, and the rates of transfer between compartments are expressed mathematically as derivatives with respect to time of the size of the compartments. In a recent work [4], a system of eight nonlinear (fractional) differential equations for a population divided into eight mutually exclusive groups was considered: susceptible, exposed, infected, hospitalized, asymptomatic but still infectious, dead but not buried, died, and completely recovered. By comparing the numerical results of this mathematical model and the real data provided by WHO, the difference in the period of 438 days analyzed is about 7 cases per day. Note that in the day 438 after the beginning of the outbreak, the number of confirmed cases is 15018.
There exist different models for the spreading of Ebola, beginning with the simplest SIR and SEIR models [2,9,10] and later more complex but also more realistic models have been considered [4,24,29]. In [25], a stochastic discrete-time Susceptible-Exposed-Infectious-Recovered (SEIR) model for infectious diseases is developed with the aim of estimating parameters from daily incidence and mortality time series for an outbreak of Ebola in the Democratic Republic of Congo in 1995. In [24], the authors use data from two epidemics (in Democratic Republic of Congo in 1995 and in Uganda in 2000) and built a SEIHFR (Susceptible-Exposed-Infectious-Hospitalized-F(dead but not yet buried)-Removed) mathematical model for the spread of Ebola haemorrhagic fever epidemics taking into account transmission in different epidemiological settings (in the community, in the hospital, during burial ceremonies). In [5], the authors propose a SIRD (Susceptible-Infectious-Recovered-Dead) mathematical model using classical and beta derivatives. In this model, the class of susceptible individuals does not consider new born or immigration. The study shows that, for small portion of infected individuals, the whole country could die out in a very short period of time in case there is no good prevention. In [3], a fractional order SEIR Ebola epidemic model is proposed and the authors show that the model gives a good approximation to real data published by WHO, starting from March 27th, 2014.
Optimal control is a mathematical theory that emerged after the Second World War with the formulation of the celebrated Pontryagin maximum principle, responding to practical needs of engineering, particularly in the field of aeronautics and flight dynamics [31]. In the last decade, optimal control has been largely applied to biomedicine, namely to models of cancer chemotherapy (see, e.g., [23]), and recently to epidemiological models [15,26,36].
In [32], the authors present a comparison between SIR and SEIR mathematical models used in the description of the Ebola virus propagation. They applied optimal control techniques in order to understand how the spread of the virus may be controlled, e.g., through education campaigns, immunization or isolation. In [1], the authors introduce a deterministic SEIR type model with additional hospitalization, quarantine and vaccination components in order to understand the disease dynamics. Optimal control strategies, both in the case of hospitalization (with and without quarantine) and vaccination, are used to predict the possible future outcome in terms of resource utilization for disease control and the effectiveness of vaccination on sick populations. Both in [1] and [32], the authors study optimal control problems with L 2 cost functionals without any state or control constraints. Here, we modify the model analyzed in [4] in order to consider optimal control problems with vaccination constraints. More precisely, we introduce an extra variable for the number of vaccines used, and we compare the hypothetical results if the vaccine were available at the beginning of the outbreak with the results of the model without vaccines. Firstly, we consider an optimal control problem with an end-point state constraint, that is, the total number of available vaccines, in a fixed period of time, is limited. Secondly, we analyze an optimal control problem with a mixed state constraint, in which there is a limited supply of vaccines at each instant of time for a fixed interval of time. Both optimal control problems have been analytically solved. Moreover, we have performed a number of numerical simulations in three different scenarios: unlimited supply of vaccines; limited total number of vaccines to be used; and limited supply of vaccines at each instant of time. From the results obtained in the first two cases, when there is no limit in the supply of vaccines or when the total number of vaccines used is limited, the optimal vaccination strategy implies a vaccination of 100% of the susceptible population in a very short period of time (smaller than one day). In practice, this is a very difficult task because limitations in the number of vaccines and also in the number of humanitarian and medical teams in the affected regions are common. In this direction, the third analyzed case is extremely important since we consider a limited supply of vaccines at each instant of time.
The paper is organized as follows. In Section 2, we recall a mathematical model for Ebola virus. In Section 3, the introduction of effective vaccination for Ebola virus is modeled. An optimal control problem with an end-point state constraint is formulated and solved analytically in Section 4, which models the case where the total number of available vaccines in a fixed period of time is limited. In Section 5, the limited supply of vaccines at each instant of time for a fixed interval of time is mathematically translated into an optimal control problem with a mixed state control constraint. A closed form of the unique optimal control is given. In Section 6, we solve numerically the optimal control problems proposed in Sections 4 and 5. Finally, we end with Section 7 of discussion of the results.
2. Initial mathematical model for Ebola. The total population N under study is subdivided into eight mutually exclusive groups: susceptible (S), exposed (E), infected (I), hospitalized (H), asymptomatic but still infectious (R), dead but not buried (D), buried (B), and completely recovered (C). This model is adapted from [20] and analyzed in [4], where the birth and death rate are assumed to be equal and are denoted by µ, and the contact rate of susceptible individuals with infective, dead, hospitalized and asymptomatic individuals are denoted by β i , β d , β h and β r , respectively. Exposed individuals become infectious at a rate σ. The per capita rate of progression of individuals from the infectious class to the asymptomatic and hospitalized classes are denoted by γ 1 and τ , respectively. Individuals in the dead class progress to the buried class at a rate δ 1 . Hospitalized individuals progress to the buried class and to the asymptomatic class at rates δ 2 and γ 2 , respectively. Asymptomatic individuals become completely recovered at a rate γ 3 . Infectious individuals progress to the dead class at a fatality rate . Dead and buried bodies are incinerated at a rate ξ. We assume that the total population, N = S + E + I + R + H + D + B + C, is constant, that is, the birth and death rates, both denoted by µ, are equal to the incineration rate ξ. The model is mathematically described by the following system of eight nonlinear ordinary differential equations: In Fig. 1, we give a flowchart presentation of model (1). In this flowchart, we identify the compartmental classes as well as the parameters appearing in the model. Moreover, the values of the parameters are given in Table 1. The basic reproduction number (that is, the number of cases one case generates on average over the course of its infectious period, in an otherwise uninfected population) of model (1) can be computed using the associated next-generation matrix method [13]. It is obtained as the spectral radius of the following matrix, known as the next-generation-matrix: a 1 = γ 2 +δ 2 +µ, a 2 = γ 3 +µ, a 3 = δ 1 +ξ, a 4 = σ+µ, and a 5 = γ 1 + +τ +µ.
Therefore, the basic reproduction number R 0 is given by As it is well-known, if the basic reproduction number R 0 < 1, then the infection will stop in the long run; but if R 0 > 1, then the infection will spread in population.
In this section, we have recalled a model for describing the Ebola virus transmission. Now we want to address the question about how to introduce vaccination as a prevention measure. This is analyzed in the next section.
3. Mathematical model for Ebola with vaccination. We now introduce vaccination of the susceptible population with the aim of controlling the spread of the disease. We assume that the vaccine is effective so that all vaccinated susceptible individuals become completely recovered (see, e.g., [6,27] for vaccination in a SEIR model that corresponds to a system of four nonlinear ordinary differential equations). Let us introduce in model (1) a control function u(t), which represents the percentage of susceptible individuals being vaccinated at each instant of time t with t ∈ [0, t f ]. Unfortunately, in many situations, the number of available vaccines does not fulfill the necessities in order to eradicate the disease. In this paper, we consider limitations on the total number of vaccines during a fixed interval of time [0, t f ] and on the number of available vaccines at each instant of time t with t ∈ [0, t f ]. In order to translate this real situation mathematically, we introduce an extra variable W that denotes the number of vaccines used: dW Hence, the model with vaccination is given by the following system of nine nonlinear ordinary differential equations: (2) In model (2), the vaccination parameter is fixed. In the next section, we address the question of how to choose this parameter in an optimal way along time. 4. Optimal control with an end-point state constraint. We start by considering the case where the total number of available vaccines, in a fixed period of time, is limited. We formulate and solve analytically such optimal control problem with end-point state constraint, which will be then solved numerically in Section 6.
We consider the model with vaccination (2) and formulate the optimal control problem with the aim to determine the vaccination strategy u over a fixed interval of time [0, t f ] that minimizes the cost functional where the constants w 1 and w 2 represent the weights associated with the number of infected individuals and on the cost associated with the vaccination program, respectively. We assume that the control function u takes values between 0 and 1. When u(t) = 0, no susceptible individual is vaccinated at time t; if u(t) = 1, then all susceptible individuals are vaccinated at t. Let ϑ denote the total amount of available vaccines in a fixed period of time [0, t f ]. This constraint is represented by Let The optimal control problem consists to find the optimal trajectoryx associated with the optimal controlũ, satisfying the control system (2), the initial conditions x(0) = (18000, 0, 15, 0, 0, 0, 0, 0, 0) (see [4]), the constraint (4), and where the controlũ ∈ Ω minimizes the objective functional (3) with The existence of an optimal controlũ and associated optimal trajectoryx comes from the convexity of the integrand of the cost function (3) with respect to the control u and the Lipschitz property of the state system with respect to state variables x (see, e.g., [8,16] for existence results of optimal solutions). According to the Pontryagin Maximum Principle [31], ifũ ∈ Ω is optimal for the problem (2), (3) with initial conditions (5) and fixed final time t f , then there exists λ ∈ AC([0, t f ]; R 9 ), λ(t) = (λ 1 (t), . . . , λ 9 (t)), called the adjoint vector, such thaṫ where the Hamiltonian H 1 is defined by where b = (−1 0 0 0 0 0 0 1 1) T and Z = 0 with 0 the 8 × 9 null matrix. The minimization condition holds almost everywhere on [0, Solving the minimality condition (7) on the interior of the set of admissible controls Ω gives where the adjoint functions satisfy Since W has initial and terminal conditions, the adjoint functionλ 9 , associated with the state variable W , has no transversality condition. From (8),λ 9 (t) ≡ k, where the constant k must be such that the end point conditions W (0) = 0 and W (t f ) = ϑ are satisfied. As the optimal controlũ can take values on the boundary of the control set [0, 1], the optimal controlũ must satisfỹ The optimal controlũ given by (9) is unique due to the boundedness of the state and adjoint functions and the Lipschitz property of systems (2) and (8).
We would like to note that if we consider the optimal control problem without any restriction on the number of available vaccines, that is, to find the optimal solution (x,ũ), withũ ∈ Ω, which minimizes the cost functional (3) subject to the control system (2), initial conditions (5), and free final conditions (x 1 (t f ), . . . , x 9 (t f )), then the adjoint functions (λ 1 , . . . , λ 9 ) must satisfy transversality conditions λ i (t f ) = 0, i = 1, . . . , 9, and, sinceλ 9 = 0, the optimal control is given bỹ In a concrete situation, the number of available vaccines is always limited. Therefore, it is also important to study the optimal control problem with such kind of constraints. This is done in Section 5. Both problems, with and without constraints, are numerically solved in Section 6.
6. Numerical simulations. We start the numerical simulations by considering an intervention of 90 days, initial conditions given in (5), and the evolution of cumulative confirmed cases based on the data from the World Health Organization (WHO), following all the reports of the disease in the three main affected countries of Western Africa of the 2014 Ebola outbreak, namely, Liberia, Guinea and Sierra Leone. The model (1) with the parameter values from Table 1 fits the real data from WHO, see Fig. 2a. We would like to emphasize that we are just considering the initial period of spreading of the disease, in which the vaccination should be introduced. Looking to a longer period of time (as considered in [4]), then the model fits quite well the real data: in [4, Figure 2] the 2 norm of the difference between the real data and the prediction is 3181, which gives an error of less than 7.3 cases per day, as compared with about 15,000 cases at the end of the outbreak.  Assuming that, in the near future, an effective vaccine against the Ebola virus will be available, as expected by WHO by the end of 2015, we study three different scenarios, which illustrate limitations on the number of vaccines available and on the capacity of administration of the vaccines by the health care services and humanitarian teams working in the affected countries. The three vaccination scenarios are the following: unlimited supply of vaccines; limited total number of vaccines to be used; and limited supply of vaccines at each instant of time.
6.1. Unlimited supply of vaccines. Assume w 1 = w 2 = 1. If we administer from an unlimited supply of vaccines, then the number of total individuals who have an active infection during the 90 days is a decreasing function in time, and is equal to 3.56 individuals at the final time (see Fig. 2b).
If we compare the case where there is no vaccination with the opposite case of unlimited supply of vaccines, we observe that at the end of 90 days the class of completely recovered individuals has approximately 86.5 individuals in the case of no vaccination and 13468 in the case of unlimited vaccination, which represents 74.82 per cent of the total population (see Fig. 3-4). If vaccines are available, then the number of individuals that develop active disease is less than one at the end of 7.5 days and less that 0.1 at the end of 32.4 days. In the case of no vaccination, the class of active infected individuals has 61.7 individuals at the end of 90 days (see Fig. 3). If no vaccination is provided, then the number of deaths, hospitalizations and burials increases from 1.2 to 262.6, when compared to the case of unlimited supply of vaccines (see Fig. 4). The optimal vaccination policy suggested by the solution   Table 1. of the optimal problem, implies a vaccination of 100 per cent of the susceptible population for approximately 1.62 days followed by a fast reduction of the fraction of susceptible population that is vaccinated. This is based on the fact that the vaccine is effective and once all the susceptible population is vaccinated in a short period of time, then the number of susceptible individuals immediately decreases, since they are transferred to the class of completely recovered individuals, as well as the need of vaccination (see Fig. 5a). The previous results show the importance of an effective vaccine for Ebola virus and the very good results that can be attained if the number of available vaccines satisfies the needs of the population. 6.2. Limited total number of vaccines. In Fig. 5b, we observe that at the end of 90 days, 33786 vaccines were used, if the supply of vaccines is unlimited. In this section, we consider the case where the total number of vaccines used in the 90 days period is limited. We consider the case where the total number of vaccines available is lower or equal than the initial number of susceptible individuals (W (90) ≤ 10000, W (90) ≤ 11000, W (90) ≤ 13000, W (90) ≤ 15000, W (90) ≤ 16000, W (90) ≤ 18000) and the case where the total number of vaccines available is bigger than the initial number of susceptible individuals (W (90) ≤ 20000). We first consider w 1 = w 2 = 1. The cumulative confirmed cases (see Fig. 6   the period of 90 days to be equal to 11000, 13000, 15000, 16000 and 18000, then we observe that the optimal control u remains more time at the maximum value 1 when the supply of vaccines is bigger, which means that when the total number of available vaccines is increased there will be resources to vaccinate all susceptible individuals for a longer period of time, which implies a bigger reduction of the number of individuals who get infected by the virus (see Fig. 9a and 9b for the optimal control strategy and respective zoom in the period of vaccination). Consider now the case where the weight constant associated with the cost of implementation of the vaccination strategy, designated by the optimal control u, is bigger than one, for example, consider w 1 = 1 and w 2 = 50, and w 1 = 1 and w 2 = 500. To simplify, consider in both cases W (90) ≤ 10000 and W (90) ≤ 20000. When we increase the weight constant w 2 , the maximum value attained by the optimal control becomes lower than one (see Fig. 10). In the case W (90) ≤ 10000 for w 2 = 50, the optimal control starts with the value u(0) = 0.54 and is a decreasing function with a cost function 344.3. At the end of approximately 3.7 days, the control remains equal to zero. For w 2 = 500, the optimal control starts with the value u(0) = 0.16 and is also a decreasing function, with a cost 399.62. At the end of 13.5 days, it remains equal to zero. The behavior of the optimal state variables S, E, I, R, D, H, B and C are similar.  From previous results, we observe that when there is no limit on the supply of vaccines, or when the total number of used vaccines is limited, the optimal vaccination strategy implies a vaccination of 100 per cent of the susceptible population in a very short period of time, sometimes smaller than one day. But we know that in practice this is a very difficult task, since there are limitations in the number of vaccines available and also in the number of health care workers or humanitarian teams in the regions affected by Ebola virus with capacity to vaccinate such a big number of individuals almost simultaneously. From this point of view, it is important to study the case where there is a limited supply of vaccines at each instant of time. In this section, we consider w 1 = w 2 = 1, a shorter interval of time [0, t f ], with t f = 10, 15, 16, and we assume that at each instant of time there exist only 1000, 1200 and 900 available vaccines, respectively. From our point of view, these numbers of available vaccines at each instant of time and the number of days considered, correspond to possible real scenarios, which are possible to implement in a concrete endemic region and at the same time characterize lack of human and material resources to vaccinate the susceptible population in a short period of time. From the numerical simulations, for such mixed constraints, the number of cumulative confirmed cases increases with time (see Figure 11a). The cost associated with the vaccination campaign, associated with the solution of the optimal control problem with the mixed constraint S(t)u(t) ≤ 1000, t ∈ [0, 10], is equal to 45.8879. Such solution is the less costly of the three considered, followed by the constraint S(t)u(t) ≤ 1200 for t ∈ [0, 15] with a cost of 55.079. The most expensive vaccination strategy is the one associated with the mixed constraint S(t)u(t) ≤ 900, t ∈ [0, 16], with a cost of 59.109. The strategy associated with the constraint S(t)u(t) ≤ 1000 is the one where a lowest number of susceptible individuals completely recover through vaccination, with 7540.9 individuals in the class C at the end of 10 days. If we consider that at each instant of time there are 1200 vaccines available during a period of 15 days, then 12438 completely recover. This is the strategy with more individuals in the class C. If we consider 16 days, but only 900 vaccines available for each instant of time, then only 10839 individuals completely recover (see Fig. 11b). For all three mixed constraint situations, the number of individuals in the classes E, I, R, D, H and B does not change significantly (therefore, the figures with these classes are omitted). As the number of available vaccines represent a small percentage of the susceptible population, in the three cases the optimal vaccination strategies for the constraints S(t)u(t) ≤ 1200 and S(t)u(t) ≤ 900 suggest that the percentage of the susceptible population that is vaccinated is always inferior than 18 percent. In the case of the constraint S(t)u(t) ≤ 1000, this percentage is always inferior to 8 percent (see Fig. 11c). 7. Discussion. We assume that, in a near future, an effective vaccine against the Ebola virus will be available. Under this assumption, three different scenarios have been studied: unlimited supply of vaccines; limited total number of vaccines to be used; and limited supply of vaccines at each instant of time. We have solved the optimal control problems analytically and we have performed a number of numerical simulations in the three aforementioned vaccination scenarios.
Some authors have already considered the optimal control problem with vaccination for Ebola disease, but always with unlimited supply of vaccines [1,32]. It turns out that the solution to this mathematical problem is obvious: the solution consists to vaccinate all susceptible individuals in the beginning of the outbreak. This is a very particular case of our work, investigated in Section 6.1 (see Figure 5). If vaccines are available without any restriction, then one could completely eradicate Ebola in a very short period of time. These results show the importance of an effective vaccine for Ebola virus and the very good results that can be attained if the number of available vaccines satisfy the needs of the population. Unfortunately, such situation is not realistic: in case an effective vaccine for Ebola virus will appear, there always will be restrictions on the number of available vaccines as well as constraints on how to inoculate them in a proper way and in a short period of time; economic problems might also exist.
In our work, for first time in the literature of Ebola, an optimal control problem with state and control constraints has been considered. Mathematically, it represents a health public problem of limited total number of vaccines. The results obtained in Section 6.2 provide useful information on the number of vaccines to be bought, in order to reduce the number of new infections with minimum cost. For example, the results between 10000 and 20000 vaccines (in 90 days) are completely different. With 10000 vaccines, the number of cumulative infected cases continues to increase, while with 20000 vaccines it is already possible to decrease the new infections. The optimal solution, in this case, is similar to the case of unlimited supply of vaccines, that is, it implies a vaccination of 100 per cent of the susceptible population in a very short period of time. In practice, this is an unrealistic task, due to the necessary number of vaccines and humanitarian teams in the regions affected by Ebola. Therefore, we conclude that it is important to study the case where there is a limited supply of vaccines at each instant of time. This was investigated in Section 6.3. This situation is much richer and the optimal control solution is not obvious. For a given number of available vaccines at each instant of time, we have a different solution, which is the optimal rate of susceptible individuals that should be vaccinated. In this case, the optimal control implies the vaccination of a small subset of the susceptible population. It remains the ethical problem of how to choose the individuals to be vaccinated.