Ebola Model and Optimal Control with Vaccination Constraints

TL;DR

This paper develops and analyzes mathematical models for Ebola transmission incorporating vaccination strategies under various supply constraints, providing insights into optimal control measures to reduce the disease spread.

Contribution

It introduces and analytically solves new optimal control models for Ebola vaccination with supply limitations, enhancing understanding of effective vaccination strategies.

Findings

Vaccination significantly reduces Ebola transmission in simulations.

Limited vaccine supply impacts optimal vaccination timing and coverage.

Models fit real outbreak data well under different supply scenarios.

Abstract

The Ebola virus disease is a severe viral haemorrhagic fever syndrome caused by Ebola virus. This disease is transmitted by direct contact with the body fluids 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 fixed 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 fixed 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 vaccination, which has recently been shown to fit the real data. Three vaccination scenarios have been considered for our numerical simulations, namely: unlimited supply of vaccines; limited total number of vaccines; and limited supply of vaccines at each instant of time.

Tables

Table 1. Table 1. Parameter values for model ( 1 ), corresponding to a basic reproduction number R 0 = 2.287 subscript 𝑅 0 2.287 R_{0}=2.287 . The values of the parameters come from [ 7 , 14 , 22 , 24 , 28 , 34 , 35 ] .

Symbol	Description	Value
$\sigma$	per capita rate at which exposed individuals become infectious	$1/11.4$
$\mu$	death rate	$14/1000$
${\beta }_i$	contact rate of infective and susceptible individuals	$0.14$
${\beta }_d$	contact rate of infective and dead individuals	$0.40$
${\beta }_h$	contact rate of infective and hospitalized individuals	$0.29$
${\beta }_r$	contact rate of infective and asymptomatic individuals	$0.185$
${\gamma }_1$	per capita rate of progression of individuals from the infectious class
	to the asymptomatic class	$1/10$
$ϵ$	fatality rate	$1/9.6$
${\delta }_1$	per capita rate of progression of individuals from the dead class
	to the buried class	$1/2$
${\delta }_2$	per capita rate of progression of individuals from the hospitalized class
	to the buried class	$1/4.6$
${\gamma }_2$	per capita rate of progression of individuals from the hospitalized class
	to the asymptomatic class	$1/5$
$\tau$	per capita rate of progression of individuals from the infectious class
	to the hospitalized class	$1/5$
${\gamma }_3$	per capita rate of progression of individuals from the asymptomatic class
	to the completely recovered class	$1/30$
$\xi$	incineration rate	$14/1000$

Equations

$$\{\begin{array}{l}\frac{dS}{dt}=\mu N-\frac{{\beta }_i}{N}SI-\frac{{\beta }_h}{N}SH-\frac{{\beta }_d}{N}SD-\frac{{\beta }_r}{N}SR-\mu S,\\ \frac{dE}{dt}=\frac{{\beta }_i}{N}SI+\frac{{\beta }_h}{N}SH+\frac{{\beta }_d}{N}SD+\frac{{\beta }_r}{N}SR-\sigma E-\mu E,\\ \frac{dI}{dt}=\sigma E-({\gamma }_1+ϵ+\tau +\mu )I,\\ \frac{dR}{dt}={\gamma }_{1}I+{\gamma }_{2}H-({\gamma }_3+\mu )R,\\ \frac{dD}{dt}=ϵI-{\delta }_{1}D-\xi D,\\ \frac{dH}{dt}=\tau I-({\gamma }_2+{\delta }_2+\mu )H,\\ \frac{dB}{dt}={\delta }_{1}D+{\delta }_{2}H-\xi B\\ \frac{dC}{dt}={\gamma }_{3}R-\mu C.\end{array}$$

$$F{V}^{-1}=\left(\begin{array}{ccccc}{A}_{11}& A_{12}& A_{13}& A_{14}& A_{15}\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right),$$

$$A_{11}=\frac{a_3{\beta }_i\sigma a_1{a}_4+a_3{\beta }_r\sigma \left(a_4{\gamma }_1+\tau {\gamma }_2\right)+{\beta }_dϵ\sigma a_1{a}_4+a_3{\beta }_h\tau \sigma a_1}{a_1{a}_2{a}_3{a}_4{a}_5},$$

$$A_{12}=\frac{{\beta }_r}{a_1}+\frac{{\beta }_r(a_4{\gamma }_1+\tau {\gamma }_2)}{a_1{a}_2{a}_3}+\frac{{\beta }_dϵ}{a_2{a}_3}+\frac{{\beta }_h\tau }{a_2{a}_4},$$

$$A_{13}=\frac{{\beta }_r}{a_1},A_{14}=\frac{{\beta }_d}{a_3},A_{15}=\frac{{\beta }_r{\gamma }_2}{a_1{a}_4}+\frac{{\beta }_h}{a_4}$$

$$a_1={\gamma }_2+{\delta }_2+\mu ,a_2={\gamma }_3+\mu ,a_3={\delta }_1+\xi ,a_4=\sigma +\mu ,\text{and}{a}_5={\gamma }_1+ϵ+\tau +\mu .$$

$$\begin{gathered} R_0=\frac{\sigma }{a_1{a}_2{a}_3{a}_4{a}_5}[a_3\left({\beta }_i({\gamma }_3{a}_1+\mu ({\gamma }_2+{\delta }_2)+{\mu }^2)+{\beta }_r({\gamma }_1{a}_1+{\gamma }_2\tau )+{\beta }_h\tau a_2\right) \\ +{\beta }_dϵ\left(({\gamma }_3{a}_1+\mu ({\gamma }_2+{\delta }_2))+{\mu }^2\right)]. \end{gathered}$$

$$\frac{dW}{dt}(t)=u(t)S(t),\text{subject to the initial condition }W(0)=0.$$

$$\{\begin{array}{l}\frac{dS}{dt}(t)=\mu N-\frac{{\beta }_i}{N}S(t)I(t)-\frac{{\beta }_h}{N}S(t)H(t)-\frac{{\beta }_d}{N}S(t)D(t)\\ -\frac{{\beta }_r}{N}S(t)R(t)-\mu S(t)-S(t)u(t),\\ \frac{dE}{dt}(t)=\frac{{\beta }_i}{N}S(t)I(t)+\frac{{\beta }_h}{N}S(t)H(t)+\frac{{\beta }_d}{N}S(t)D(t)\\ +\frac{{\beta }_r}{N}S(t)R(t)-\sigma E(t)-\mu E(t),\\ \frac{dI}{dt}(t)=\sigma E(t)-({\gamma }_1+ϵ+\tau +\mu )I(t),\\ \frac{dR}{dt}(t)={\gamma }_{1}I(t)+{\gamma }_{2}H(t)-({\gamma }_3+\mu )R(t),\\ \frac{dD}{dt}(t)=ϵI(t)-({\delta }_1+\xi )D(t),\\ \frac{dH}{dt}(t)=\tau I(t)-({\gamma }_2+{\delta }_2+\mu )H(t),\\ \frac{dB}{dt}(t)={\delta }_{1}D(t)+{\delta }_{2}H(t)-\xi B(t),\\ \frac{dC}{dt}(t)={\gamma }_{3}R(t)-\mu C(t)+S(t)u(t),\\ \frac{dW}{dt}(t)=S(t)u(t).\end{array}$$

$$J(u)={\int }_0^{t_f}\left[w_{1}I(t)+w_2{u}^2(t)\right]dt,$$

$$W(t_f)\le \vartheta .$$

$$x(t)=(x_1(t),\dots ,x_9(t))=\left(S(t),E(t),I(t),R(t),D(t),H(t),B(t),C(t),W(t)\right)\in {\mathbb{R}}^9.$$

$$x(0)=(18000,0,15,0,0,0,0,0,0)$$

$$\Omega=\biggl{\{}u(\cdot)\in L^{\infty}(0,t_{f})\,|\,0\leq u(t)\leq 1\biggr{\}}.$$

$$\dot{x}=\frac{\partial {\mathcal{H}}_1}{\partial \lambda }\text{and}\dot{\lambda }=-\frac{\partial {\mathcal{H}}_1}{\partial x},$$

$${\mathcal{H}}_1(x,u,\lambda )=w_1{x}_3+w_2{u}^2(t)+\lambda \left(f(x)+Ax+Bxu\right)$$

$$f=\left(f_1{f}_{2}0000000\right),$$

$$f_1=\mu N-\frac{{\beta }_i}{N}S(t)I(t)-\frac{{\beta }_h}{N}S(t)H(t)-\frac{{\beta }_d}{N}S(t)D(t)-\frac{{\beta }_r}{N}S(t)R(t),$$

$$f_2=\frac{{\beta }_i}{N}S(t)I(t)+\frac{{\beta }_h}{N}S(t)H(t)+\frac{{\beta }_d}{N}S(t)D(t)+\frac{{\beta }_r}{N}S(t)R(t),$$

$$A=\left(\begin{array}{cccccccc}-\mu & 0& 0& 0& 0& 0& 0& 0\\ 0& -\sigma -\mu & 0& 0& 0& 0& 0& 0\\ 0& \sigma & \Lambda & 0& 0& 0& 0& 0\\ 0& 0& {\gamma }_1& -({\gamma }_3+\mu )& 0& {\gamma }_2& 0& 0\\ 0& 0& ϵ& 0& -({\delta }_1+\xi )& 0& 0& 0\\ 0& 0& \tau & 0& 0& -({\gamma }_2+{\delta }_2+\mu )& 0& 0\\ 0& 0& 0& 0& {\delta }_1& {\delta }_2& -\xi & 0\\ 0& 0& 0& {\gamma }_3& 0& 0& 0& -\mu \\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right)$$

$$B=\left(bZ\right),$$

$${\mathcal{H}}_1\left(\tilde{x}(t),\tilde{u}(t),\tilde{\lambda }(t)\right)=\underset{u\in \Omega }{\min }{\mathcal{H}}_1\left(\tilde{x}(t),u,\tilde{\lambda }(t)\right)$$

$$\tilde{u}(t)=-\frac{1}{2}\frac{\left(-{\tilde{\lambda }}_1(t)+{\tilde{\lambda }}_8(t)+{\tilde{\lambda }}_9(t)\right){\tilde{x}}_1(t)}{{\omega }_2},$$

$$\{\begin{array}{l}{\dot{\tilde{\lambda }}}_1=-{\tilde{\lambda }}_1\left(-\frac{{\beta }_i}{N}{\tilde{x}}_3-\frac{{\beta }_h}{N}{\tilde{x}}_6-\frac{{\beta }_d}{N}{\tilde{x}}_5-\frac{{\beta }_r}{N}{\tilde{x}}_4-\mu -\tilde{u}\right)\\ -{\tilde{\lambda }}_2\left(\frac{{\beta }_i}{N}{\tilde{x}}_3+\frac{{\beta }_h}{N}{\tilde{x}}_6+\frac{{\beta }_d}{N}{\tilde{x}}_5+\frac{{\beta }_r}{N}{\tilde{x}}_4\right)-{\tilde{\lambda }}_8\tilde{u}-{\tilde{\lambda }}_9\tilde{u},\\ {\dot{\tilde{\lambda }}}_2=-{\tilde{\lambda }}_2\left(-\sigma -\mu \right)-{\tilde{\lambda }}_3\sigma ,\\ {\dot{\tilde{\lambda }}}_3=-{\omega }_1+{\tilde{\lambda }}_1\frac{{\beta }_i}{N}{\tilde{x}}_1-{\tilde{\lambda }}_2\frac{{\beta }_i}{N}{\tilde{x}}_1-{\tilde{\lambda }}_3\left(-{\gamma }_1-ϵ-\tau -\mu \right)-{\tilde{\lambda }}_4{\gamma }_1-{\tilde{\lambda }}_5ϵ-{\tilde{\lambda }}_6\tau ,\\ {\dot{\tilde{\lambda }}}_4={\tilde{\lambda }}_1\frac{{\beta }_r}{N}{\tilde{x}}_1-{\tilde{\lambda }}_2\frac{{\beta }_r}{N}{\tilde{x}}_1-{\tilde{\lambda }}_4\left(-{\gamma }_3-\mu \right)-{\tilde{\lambda }}_8{\gamma }_3,\\ {\dot{\tilde{\lambda }}}_5={\tilde{\lambda }}_1\frac{{\beta }_d}{N}{\tilde{x}}_1-{\tilde{\lambda }}_2\frac{{\beta }_d}{N}{\tilde{x}}_1-{\tilde{\lambda }}_5\left(-{\delta }_1-\xi \right)-{\tilde{\lambda }}_7{\delta }_1,\\ {\dot{\tilde{\lambda }}}_6={\tilde{\lambda }}_1\frac{{\beta }_h}{N}{\tilde{x}}_1-{\tilde{\lambda }}_2\frac{{\beta }_h}{N}{\tilde{x}}_1-{\tilde{\lambda }}_4{\gamma }_2-{\tilde{\lambda }}_6\left(-{\gamma }_2-{\delta }_2-\mu \right)-{\tilde{\lambda }}_7{\delta }_2,\\ {\dot{\tilde{\lambda }}}_7={\tilde{\lambda }}_7\xi ,\\ {\dot{\tilde{\lambda }}}_8={\tilde{\lambda }}_8\mu ,\\ {\dot{\tilde{\lambda }}}_9=0.\end{array}$$

$$\tilde{u}(t)=\min \left\{1,\max \left\{0,\frac{1}{2}\frac{\left({\tilde{\lambda }}_1(t)-{\tilde{\lambda }}_8(t)-{\tilde{\lambda }}_9(t)\right){\tilde{x}}_1(t)}{{\omega }_2}\right\}\right\}.$$

$$\tilde{u}(t)=\min \left\{1,\max \left\{0,\frac{1}{2}\frac{\left({\tilde{\lambda }}_1(t)-{\tilde{\lambda }}_8(t)\right){\tilde{x}}_1(t)}{{\omega }_2}\right\}\right\}.$$

$$S(t)u(t)\le \vartheta ,\vartheta \ge 0,\text{for almost all}t\in [0,t_f],$$

$$\frac{dW}{dt}(t)=S(t)u(t)$$

$$x(t)=(x_1(t),\dots ,x_8(t))=\left(S(t),E(t),I(t),R(t),D(t),H(t),B(t),C(t)\right)\in {\mathbb{R}}^8.$$

$$\frac{dx(t)}{dt}=f(x(t))+Ax(t)+Bx(t)u(t),$$

Full text

∗Corresponding author: delfim@ua.pt (D. F. M. Torres)

Iván Area

Departamento de Matemática Aplicada II, E. E. Aeronáutica e do Espazo

Campus As Lagoas, Universidade de Vigo, 32004 Ourense, Spain

Faïçal Ndaïrou

African Institute for Mathematical Sciences (AIMS–Cameroon)

P.O. Box 608, Limbe Crystal Gardens, South West Region, Cameroon

Juan J. Nieto

Departamento de Análise Matemática, Estatística e Optimización

Facultade de Matemáticas, Universidade de Santiago de Compostela

15782 Santiago de Compostela, Spain

Cristiana J. Silva and Delfim F. M. Torres∗

Center for Research and Development in Mathematics and Applications (CIDMA)

Department of Mathematics, University of Aveiro, 3810-193 Aveiro, Portugal

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 $\mu$, and the contact rate of susceptible individuals with infective, dead, hospitalized and asymptomatic individuals are denoted by $\beta_{i}$, $\beta_{d}$, $\beta_{h}$ and $\beta_{r}$, respectively. Exposed individuals become infectious at a rate $\sigma$. The per capita rate of progression of individuals from the infectious class to the asymptomatic and hospitalized classes are denoted by $\gamma_{1}$ and $\tau$, respectively. Individuals in the dead class progress to the buried class at a rate $\delta_{1}$. Hospitalized individuals progress to the buried class and to the asymptomatic class at rates $\delta_{2}$ and $\gamma_{2}$, respectively. Asymptomatic individuals become completely recovered at a rate $\gamma_{3}$. Infectious individuals progress to the dead class at a fatality rate $\epsilon$. Dead and buried bodies are incinerated at a rate $\xi$. 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 $\mu$, are equal to the incineration rate $\xi$. The model is mathematically described by the following system of eight nonlinear ordinary differential equations:

[TABLE]

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:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

with

[TABLE]

Therefore, the basic reproduction number $R_{0}$ is given by

[TABLE]

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\in[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\in[0,t_{f}]$. In order to translate this real situation mathematically, we introduce an extra variable $W$ that denotes the number of vaccines used:

[TABLE]

Hence, the model with vaccination is given by the following system of nine nonlinear ordinary differential equations:

[TABLE]

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

[TABLE]

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 $\vartheta$ denote the total amount of available vaccines in a fixed period of time $[0,t_{f}]$. This constraint is represented by

[TABLE]

Let

[TABLE]

The optimal control problem consists to find the optimal trajectory $\tilde{x}$ associated with the optimal control $\tilde{u}$, satisfying the control system (2), the initial conditions

[TABLE]

(see [4]), the constraint (4), and where the control $\tilde{u}\in\Omega$ minimizes the objective functional (3) with

[TABLE]

The existence of an optimal control $\tilde{u}$ and associated optimal trajectory $\tilde{x}$ 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 $\tilde{u}\in\Omega$ is optimal for the problem (2), (3) with initial conditions (5) and fixed final time $t_{f}$, then there exists $\lambda\in AC([0,t_{f}];\mathbb{R}^{9})$, $\lambda(t)=\left(\lambda_{1}(t),\ldots,\lambda_{9}(t)\right)$, called the adjoint vector, such that

[TABLE]

where the Hamiltonian $\mathcal{H}_{1}$ is defined by

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with $\varLambda=-(\gamma_{1}+\epsilon+\tau+\mu)$, and

[TABLE]

where $b=\left(-1\,0\,0\,0\,0\,0\,0\,1\,1\right)^{T}$ and $Z=0$ with [math] the $8\times 9$ null matrix. The minimization condition

[TABLE]

holds almost everywhere on $[0,t_{f}]$. Moreover, transversality conditions $\lambda_{i}(t_{f})=0$, $i=1,\ldots,8$, hold. Solving the minimality condition (7) on the interior of the set of admissible controls $\Omega$ gives

[TABLE]

where the adjoint functions satisfy

[TABLE]

Since $W$ has initial and terminal conditions, the adjoint function $\lambda_{9}$, associated with the state variable $W$, has no transversality condition. From (8), $\tilde{\lambda}_{9}\equiv k$, where the constant $k$ must be such that the end point conditions $W(0)=0$ and $W(t_{f})=\vartheta$ are satisfied. As the optimal control $\tilde{u}$ can take values on the boundary of the control set $[0,1]$, the optimal control $\tilde{u}$ must satisfy

[TABLE]

The optimal control $\tilde{u}$ 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 $(\tilde{x},\tilde{u})$, with $\tilde{u}\in\Omega$, which minimizes the cost functional (3) subject to the control system (2), initial conditions (5), and free final conditions $(x_{1}(t_{f}),\ldots,x_{9}(t_{f}))$, then the adjoint functions $(\lambda_{1},\ldots,\lambda_{9})$ must satisfy transversality conditions $\lambda_{i}(t_{f})=0$, $i=1,\ldots,9$, and, since $\tilde{\lambda}_{9}=0$, the optimal control is given by

[TABLE]

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.

5. Optimal control with a mixed state control constraint

A particularly challenging situation in vaccination programs happens when there is a limited supply of vaccines at each instant of time for a fixed interval of time $[0,t_{f}]$. In order to study this health public problem, from the optimal point of view, we formulate an optimal control problem with a mixed state control constraint (see, e.g., [6]). The cost functional (3) remains (3), the one considered in previous section, as well as the set of admissible controls $\Omega$ (6). The end point state constraint (4) is replaced by the following mixed state control constraint:

[TABLE]

which should be satisfied at almost every instant of time during the whole vaccination program. Analogously to [6], we observe that in our optimal control problem the differential equation

[TABLE]

does not appear neither in the cost and in any other differential equation, nor in the mixed state control constraint. Thus, in this section, the control system does not include the last equation and $x$ is used to denote

[TABLE]

Let us consider the initial conditions (5). The control system can be rewritten in the following way:

[TABLE]

with

[TABLE]

$\varLambda=-(\gamma_{1}+\epsilon+\tau+\mu)$,

[TABLE]

where $b=\left(-1\,0\,0\,0\,0\,0\,0\,1\right)^{T}$ and $Z=0$ with [math] the $7\times 8$ null matrix, and

[TABLE]

with

[TABLE]

and

[TABLE]

It follows from Theorem 23.11 in [11] that problem (2)–(6) has a solution (see also [6]). Let $(\tilde{x},\tilde{u})$ denote such solution. To determine it, we apply the Pontryagin Maximum Principle (see, e.g., Theorem 7.1 in [12]): there exists multipliers $\lambda_{0}\leq 0$, $\lambda\in AC([0,t_{f}];{\mathbb{R}}^{8})$, and $\psi\in L^{1}([0,t_{f}];{\mathbb{R}})$, such that

• •

$\min\{|\lambda(t)|\,:\,t\in[0,t_{f}]\}>\lambda_{0}$ (nontriviality condition);

• •

$\frac{d\lambda(t)}{dt}=-\frac{\partial\mathcal{H}_{2}}{\partial x}(\tilde{x}(t),\tilde{u}(t),\lambda_{0},\lambda(t),\psi(t))$ (adjoint system);

• •

$\lambda(t)B\tilde{x}(t)+\psi(t)\tilde{x}_{1}(t)+\lambda_{0}w_{2}\tilde{u}^{2}(t)\in\mathcal{N}_{[0,1]}(\tilde{u}(t))$ a.e. and

[TABLE]

(minimality condition);

• •

$\psi(t)(\tilde{x}_{1}(t)\tilde{u}(t)-\vartheta)=0$ and $\psi(t)\leq 0$ a.e.;

• •

$\lambda(t_{f})=(0,\ldots,0)$ (transversality conditions);

where the Hamiltonian $\mathcal{H}_{2}$ for problem (2)–(6) is defined by

[TABLE]

and $\mathcal{N}_{[0,1]}(\tilde{u}(t))$ stands for the normal cone from convex analysis to $[0,1]$ at the optimal control $\tilde{u}(t)$ (see, e.g., [11]). The optimal solution $(\tilde{x},\tilde{u})$ is normal (see [6] for details), so we can choose $\lambda_{0}=1$. Analogously to previous section, we obtain a closed form of the unique optimal control $\tilde{u}$:

[TABLE]

The theoretical results obtained in Sections 4 and 5 are illustrated numerically in the next section.

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 $\ell_{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 $I(t)+R(t)+D(t)+H(t)+B(t)-\mu(N-S(t)-E(t)-C(t))$ 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 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)\leq 10000$, $W(90)\leq 11000$, $W(90)\leq 13000$, $W(90)\leq 15000$, $W(90)\leq 16000$, $W(90)\leq 18000$) and the case where the total number of vaccines available is bigger than the initial number of susceptible individuals ($W(90)\leq 20000$). We first consider $w_{1}=w_{2}=1$. The cumulative confirmed cases (see Fig. 6) increases in time in the case $W(90)\leq 10000$ and decreases in the case $W(90)\leq 20000$, $t\in[0,90]$. At the end of the 90 days period, the total number of individuals who got active infection is approximately $76$ and $9.5$ individuals in the case $W(90)\leq 10000$ and $W(90)\leq 20000$, respectively.

In the case $W(90)\leq 10000$, the optimal control takes the maximum value for less than one day (approximately 0.72 days) with a cost equal to $322.74$, and in the case $W(90)\leq 20000$, the optimal control takes the maximum value for approximately 2.2 days with a cost equal to $72.35$. The cost associated to the case $W(90)\leq 20000$ is lower than the one in the case $W(90)\leq 10000$, although more individuals are vaccinated, since the number of individuals in the class $I$ is lower. Namely, in the case $W(90)\leq 10000$, the number of individuals with active infection at the end of 90 days is equal to $I(90)=8.4$ and in the case $W(90)\leq 20000$ the respective number is equal to $I(90)=0.67$. This means that in the case $W(90)\leq 10000$, in a epidemiological scenario corresponding to a basic reproduction number greater than one, 10000 vaccines will not be enough to eradicate the disease. Additionally, if we consider the maximum value for the total number of vaccines used during 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)\leq 10000$ and $W(90)\leq 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)\leq 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.

6.3. Limited supply of vaccines at each instant of time

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)\leq 1000$, $t\in[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)\leq 1200$ for $t\in[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)\leq 900$, $t\in[0,16]$, with a cost of $59.109$. The strategy associated with the constraint $S(t)u(t)\leq 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)\leq 1200$ and $S(t)u(t)\leq 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)\leq 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.

Acknowledgments

Three reviewers deserve special thanks for helpful and constructive comments. The work of Area was partially supported by the Ministerio de Economía y Competitividad of Spain, under grants MTM2012–38794–C02–01 and MTM2016–75140–P, co-financed by the European Community fund FEDER. Ndaïrou acknowledges the AIMS-Cameroon 2014–2015 fellowship, Nieto the partial financial support by the Ministerio de Economía y Competitividad of Spain under grants MTM2010–15314, MTM2013–43014–P and MTM2016–75140–P, and Xunta de Galicia under grants R2014/002 and GRC 2015/004, co-financed by the European Community fund FEDER. Silva was supported through the Portuguese Foundation for Science and Technology (FCT) post-doc fellowship SFRH/BPD/72061/2010. The work of Silva and Torres was partially supported by FCT through CIDMA and project UID/MAT/04106/2013, and by project TOCCATA, reference PTDC/EEI-AUT/2933/2014, funded by Project 3599 – Promover a Produção Científica e Desenvolvimento Tecnológico e a Constituição de Redes Temáticas (3599-PPCDT) and FEDER funds through COMPETE 2020, Programa Operacional Competitividade e Internacionalização (POCI), and by national funds through FCT.