An optimal control model to design strategies for reducing the spread of the Ebola virus disease

1.   Introduction

Describing a phenomenon using mathematical models is useful for understanding the dynamics of human and animal diseases. Thus, modeling and simulation are important decision tools that can be used for this purpose ^[1]. They can be adapted to different diseases according to their characteristics, and, thus, can be used to handle real situations such as those arising from recent pandemics: Ebola virus disease in 2014–2016 and SARS-CoV-2 since 2019 ^[2,3,4,5]. From a mathematical point of view, there exist several studies proposing models to estimate and control the spread of Ebola virus disease, such as the ones proposed in ^[3,6,7] and ^[8]. In February 2021, an outbreak of EVD was declared by the World Health Organization (WHO) both in Guinea and the Democratic Republic of Congo. Shortly after the infections were detected, national health authorities, with support of WHO and partners, mounted a swift response in fighting these outbreaks. On June 2021, the outbreak was declared over with a total of 16 confirmed and 7 probable cases reported. Of whom, 12 people died. Recently, in August 17, 2021, WHO confirmed the first EVD case in Ivory Coast. Additionally, another suspect case and 9 contact cases have been identified and followed up. So, it becomes urgent to take imminent measures to fight against EVD. Many authors focus on this question and propose a large variety of mathematical models to study the propagation of diseases in human and animal populations ^[1,9,10,11].

Optimal control theory is another part of the mathematics that is used to manage the spread of a disease allowing to take decisions in complex biological situations ^[12]. For example, Zaman et al.^[13] applied optimal control to a vector-borne disease (such as, malaria, dengue fever or West Nile virus) with direct transmission in host population. They proved the existence of the optimal problem and established some numerical simulations to support theoretical results. Additionally, Yusuf et al. ^[12] proposed an optimal control problem based on a SIR model with vaccination and treatment as possible controls. They studied an optimal combination of vaccination and treatment strategies to minimize the cost of those control measures.

In a previous work ^[3], the authors established a deterministic spatial-temporal epidemiological model called Be-CoDiS (between-countries disease spread) to simulate the spread of human diseases in a considered area. Be-CoDiS was validated by considering the 2014–2016 West-African Ebola virus disease (EVD). Ebola is a human and primate virus disease that causes a high mortality rate (between 50% and 90%) ^[14,15]. During the period from December 2013 up to March 2016, several important outbreaks were reported in West Africa (Guinea, Liberia, Sierra Leone and Nigeria). Furthermore, some isolated cases were detected in other countries such as Mali, Senegal, the USA, the United Kingdom, Italy and Spain ^[16]. In another recent paper ^[17], the authors performed a stability and sensitivity analysis of Be-CoDiS. They first studied the equilibrium states of simplified versions of this model, limited to the cases of one or two countries and then determined their basic Reproduction ratios. Then, they established the global stability of the disease-free equilibrium (DFE) for the two simplified versions and illustrated the theoretical results by considering numerical simulations based on data from the 2014–16 West African Ebola virus epidemic.

In this work, we study an optimal control problem based on a simplified version of Be-CoDiS (limited to the case of one country). To this aim, we first formulate the optimal control problem by explaining each control variable. The controls represent the measures that can be applied in order to prevent and treat the disease (such as, prevention campaigns, vaccination, detection, hospitalization or quarantine). Those measures aim to reduce and eradicate Ebola in the population.Then, we show the existence and uniqueness of optimal solution and characterize it. We note that, with respect to existing literature focusing on models for controlling Ebola outbreaks ^[7,18], here, we consider at the same time a large variety of control measures, by using 5 controls terms.

This work is organized as follows. In Section 2, we recall the formulation of a simplified version of Be-CoDiS limited to one country. In Section 3, we first formulate the optimal control problem, then explain each of the control variable and establish some assumptions regarding the model, known as the classical regularity hypotheses. In Section 4, we characterize the optimal controls using the Pontryagin maximum principle ^[19]. The existence and uniqueness of the optimal solution are discussed in Sections 5 and 6. In Section 7, we present some numerical experiments that illustrate the pertinence of the controls on the disease spread.

2.   Considered epidemiological model

The model we are going to study is a compartmental model whose compartments are described below (see ^{[3,15,20,21,22]}):

● Susceptible (denoted by $S$): This compartment contains people who are not infected by the disease pathogen.

● Infected (denoted by $E$): People in this compartment are infected by the disease pathogen but do not have visible clinical signs (e.g., fever, hemorrhages, etc.). They cannot infect other people yet. This period is called the incubation period. At the end of it, people move to the Infectious compartment.

● Infectious (denoted by $I$): People start developing clinical signs and can infect other people. After the infectious period, which is the mean duration of a person in this compartment, infectious people are placed under sanitary cares by authorities and then classified as Hospitalized.

● Hospitalized (denoted by $H$): People in this compartment are under treatment, but can still infect other people. After the hospitalization period, people can either recover (and then pass to the Recovered compartment) or die (Dead compartment). We precise that compartment $H$ does not contain hospitalized people who cannot infect other people any more. They are included in the Recovered compartment explained below.

● Dead (denoted by $D$): People that died from the disease, who are not buried yet and can still infect other people due to contacts with their body. After a fixed average time, the body is buried.

● Recovered (denoted by $R$): Contains people who have survived the disease. People in this compartment develop a natural immunity to the disease pathogen and are no longer infectious.

After the hospitalization of infected people, various control measures may be applied by the authorities in order to control the propagation of the disease (see ^[7,8,23,24] and recently ^[18]).

Here, we restrict our study to the case of the evolution of the epidemic inside a single country. This assumption is reasonable according to existing literature studying the evolution of epidemics within a country and considering similar models (see, e.g., for COVID-19 studies ^[2,5,25]). Of course other alternatives are possible, such as, adding the interaction with more countries (as done in ^[3]) or considering smaller areas. For the sake of simplicity, we consider that $S$, $E$, $I$, $H$, $R$ and $D$ denote the ratio of people in each compartment in the considered country (rather than the total number of people).

Finally, we assume that the model coefficients are constant.

A diagram of this model for one country is shown in Figure 1.

Under these assumptions, the evolution of the epidemic is modeled by

where

● $\Lambda \in \mathbb{R}^+$ is the recruitment rate of persons in state $S$ (person.day$^{-1}$),

● $\mu \in [0, 1]$ is the mortality rate (day$^{-1}$),

● $\beta_I, \beta_H, \beta_D \in \mathbb{R}^+$ are the disease effective contact rates (day$^{-1}$.person$^{-1}$) of people in compartment $I$, $H$ and $D$ respectively,

● $\delta, \gamma \in \mathbb{R}^+$ denote the transition rates (day$^{-1}$) from compartment $E$ to $I$ and $I$ to $H$, respectively.

● $\lambda \in [0, 1]$ is the disease fatality percentage times the transition rate from compartment $H$ to compartment $D$.

● $\alpha \in \mathbb{R}^+$ is the disease survival percentage (1 minus the disease fatality percentage) times the transition rate from state $H$ to compartment $R$. For the sake of simplicity, we assume that this transition rate is the same as the one from $H$ to $D$.

● $\theta \in \mathbb{R}^+$ is the burial rate (day$^{-1}$) of people that died from Ebola.

● $\tau \in [0, 1]$ is the daily rate (%) of the movement of people in states $S$, $E$ and $R$ (people in other compartments are not supposed to travel due to their health situation) leaving the country.

As shown in ^[17], the solutions of Eq (2.1) satisfy that $(S(t), E(t), I(t), H(t), R(t), D(t))\in [0, 1]^6$.

3.   Control problem formulation

Now, we consider the controlled version of model (2.1) given by

where $T$ measures the cumulative number of infected persons (i.e., in state $E$) leaving the system during the simulated time interval (and, thus, the considered geographical area). It is used as a measure of the risk of spreading the disease outside the studied country. We point out that the last equation of system (3.1) is not coupled with the other equations. Thus, we can solve the six first equations of that system and then, the solution of the last one can be computed as follows:

for any $t, \hat{t} \geq 0$.

The controls and their associated parameters are defined as follows:

● $u_1(t) \in [0, 1]$ is the control corresponding to the vaccination campaign. First, this control is multiplied by $\dfrac{N_v}{N_0}$, which corresponds to a reasonable estimation of the maximum capacity of the vaccination campaign, where $N_v \in \mathbb{N}$ is the maximum number of persons that can be vaccinated per day and $N_0 \in \mathbb{N}$ denotes the initial size of the population at the beginning of the simulation. Secondly, $u_1(t)$ is also multiplied by the function

It is easy to prove that $\varphi\in {\mathcal C}^\infty (\mathbb{R})$. We point out that $\varphi (S(t)) = 1$ when $S(t) > 1/(2N_0)$ (in particular, when there is, at least, one person in the Susceptible compartment) and tends to 0 when $S(t)$ tends to 0. Function $\varphi$ is a filter used to avoid negative values of $S$ (notice that $S(t) = 0$ implies that $\dfrac{{\rm d} S}{{\rm d}t}(t) \geq 0$). Additionally, this function does not alter the maximum capacity $\frac{N_v}{N_0}$ of the vaccination campaign $u_1(t)$, until reaching the case of less than one person in the Susceptible compartment.

● $u_2(t) \in [0, 1]$ is the control corresponding to prevention campaigns. Those campaigns (such as educational campaign) aim to reduce effective contacts with infected persons and people that died from Ebola. From a modeling point of view, $u_2$ may reduce the values $\beta_I$ and $\beta_D$ by a maximum of $c_{\beta_I }\beta_{I}$ and $c_{\beta_D }\beta_{D}$, respectively. Here, $c_{\beta_I }, c_{\beta_D } \in \left[0, 1 \right]$ are the maximum percentage of reduction for $\beta_{I}$ and $\beta_{D}$ respectively.

● $u_3(t) \in [0, 1]$ is the control corresponding to the increase of sanitary measures in hospital (such as better treatment and better sanitary conditions). Those measures aim to reduce effective contacts with infected patients in hospital and to reduce the disease mortality. From a modeling point of view, $u_3$ may reduce the values $\beta_H$ and $\lambda$ by a maximum of $c_{\beta_H }\beta_{H}$ and $c_{\lambda }\lambda$, where $c_{\beta_H }, c_{\lambda }\in \left[0, 1 \right]$ are the maximum percentage of reduction for $\beta_H$ and $\lambda$ respectively. Furthermore, $u_3$ is also used to increase the transition rate from compartment $H$ to compartment $R$, through the term $c_{\lambda} u_3(t)H(t)$ (corresponding to the increase of the survival rate).

● $u_4(t) \in [0, 1]$ is the control corresponding to early detection campaigns. Those campaigns aim reducing the time between the apparition of clinical signs and the hospitalization. From a modeling point of view, $u_4$ is used to increase the transition rate from compartment $I$ to compartment $H$ (and thus, decreasing the time of a person in state $I$) through the term $u_4(t)\eta_H I(t)$. Here, $\eta_H$ (day$^{-1}$) corresponds to the maximum value that the transition rate from $I$ to $H$ can be increased. As a consequence, the time of a person in state $H$ can be increased. Thus, the transition rates from compartment $H$ to compartment $R$ and from compartment $H$ to compartment $D$ are reduced through the terms $- u_4(t)\eta_R H(t)$ and $- u_4(t)\eta_D H(t)$, respectively. Here, $\eta_R$ (day$^{-1}$) and $\eta_D$ (day$^{-1}$) correspond to the maximum value that the transition rate from $H$ to $R$ and from $H$ to $D$ can be decreased, respectively. Notice that, $\eta_H$ can be estimated taking into account that $\dfrac{1}{\gamma+\eta_H}$ is the minimum average number of transition days that one can expect between compartments $I$ and $H$ (a similar idea can be applied to estimate $\eta_D$ and $\eta_R$).

● $u_5(t) \in [0, 1]$ is the control corresponding to the application of quarantine measures (movement between geographical areas are limited). Those control measures aim reducing the risk of spreading the disease outside the considered country. From a modeling point of view, $u_5$ is used to decrease the value of $\tau$ (i.e., the movement of people leaving the country). More precisely, we consider the parameter $c_{\tau} \in [0, 1]$ which corresponds to the maximum percentage of reduction for $\tau$ that can be reached due to the application of those control measures. Then, $\tau$ is multiplied by $(1-c_{\tau} u_5(t))$.

A diagram of this controlled model for one country is shown in Figure 2.

Let us denote $u = (u_1, u_2, u_3, u_4, u_5)$ and $x = (S, E, I, H, R, D, T)$. Then system (3.1) can be written as:

with a suitable function

given by:

Considering those controls, we are interested in minimizing the following cost function:

where $t_{\rm f} > 0$ is the considered final time and $c_1, ..., c_5$, $K_1, ..., K_6$ are weight coefficients. Parameters $c_1, ..., c_5$ correspond to economical costs related to the implementation of the control measures. We note that we could remove one of the constants $c_i$ or $K_j$ and get an equivalent problem, but we keep all of them for simplicity.

Therefore, we consider the following standard control problem of the form (see p 437 in ^[19]):

where $U = [0, 1]^5.$

and $\Omega = \{(S, E, I, H, R, D, T)\in [0, 1]^7 \; : \; S+E+I+H+R+D = 1 \}$.

Following ^[19], any pair $(x, u)$, with $u(t)\in U$, for all $t \in \left[0, t_{\rm f} \right]$ and $x$ solution of system (3.2), is called a process of the underlying control system (3.2).

Remark 1. We could consider $J$ as a function depending only on $u, S, E, D, T$ (i.e., $J(u, S, E, D, T)$) and even depending only on $u$ (i.e., $J(u)$), since once we set $x(0)$, for every control $u$ we assume that there exists a unique solution $x$ of system (3.2).

We can also rewrite the cost function (see p 436 in ^[19]), as

where $A(x) = K_1 E+ K_2 D+ K_3 T$ is the endpoint cost corresponding to the minimization of $E, D$ and $T$ at the final time (see also ^[26,27] and ^[28] for more cost functions of this form) and

is the running cost.

It is easy to see that the function $A$ is continuously differentiable. The functions $f$ and $\Lambda$ are continuous and admit derivatives relative to $x$, denoted by $D_xf(x, u)$ and $D_x\Lambda(x, u)$ respectively, which are also continuous. Therefore, this problem satisfies the classical regularity hypotheses (see, e.g., p 437, Section 22 in ^[19]). Those assumptions imply that the cost $J(x, u)$ is well defined for any process $(x, u)$ (see, again ^[19]).

In the following section, we aim to characterize $(x^*, u^*)$ minimizing the cost function $J$. Then, in Sections 5 and 6, we show the existence and uniqueness of solution of the minimization problem (3.4).

4.   Characterization of the optimal control

Here, we apply a version of Pontryagin Maximum Principle given in ^[19], p 438 to characterize the solutions of control problem (3.4). See ^[29] for other versions.

Definition 1. (see p 437 in ^[19]) Let $(x^*, u^*)$ be a given process satisfying the constraints of the minimization problem (3.4). This process is called a local minimizer provided that, there exists some $\epsilon > 0$ such that, for any other process $(x, u)$ satisfying the constraints of minimization problem (3.4), as well as $\|x-x^*\|_\infty\leq \epsilon$^*, we have that $J(x^*, u^*)\leq J(x, u)$. In this terminology, $u^*$ is called an optimal control and $x^*$ is called an optimal trajectory.

^*Given a continuous function $x:[0, t_{\rm f}]\rightarrow \mathbb{R}^7$, we denote by $\left\| x\right\| _\infty$ the norm defined by $\| x\| _\infty = \max_{t\in [0, t_{\rm f}]} \{ \left| x_i(t)\right|, i \in \{ 1, ..., 7\} \}.$

Before formulating the Pontryagin maximum principle, we recall some definitions and tools related to our problem (see p 436–439 in ^[19]).

In addition to system (3.2), which can be written as:

we consider the following system of equations in the auxiliary variables

$p = (p_1, p_2, p_3, p_4, p_5, p_6, p_7)$:

Since the system of Eq (4.2) is linear and homogeneous with respect to $p$, given the functions $x$ and $u$, for any final condition $p_i(t_{\rm f})$, it admits a unique solution, denoted by $p = (p_1, p_2, p_3, p_4, p_5, p_6, p_7)$, which is defined on the entire interval $0\leq t\leq t_{\rm f}$ on which $x$ and $u$ are defined.

In order to combine Eqs (4.1) and (4.2), we consider the following function called the Hamiltonian of the system:

Thus, systems (4.1) and (4.2) can be rewritten as

For $p\in \mathbb{R}^7$ and $x\in \mathbb{R}^7$, we denote

From the Corollary 22.3 of Theorem 22.2 in ^[19] (p 438–439), we state the following theorem which gives necessary conditions for the solutions of the minimization problem (3.4):

Theorem 1. [Pontryagin maximum principle]

Let the process $(x^*, u^*)$ be a local minimizer for problem (3.4) under the classical regularity hypotheses mentioned at the end of Section 3, with $U$ bounded. Then, there exists a function $p:\left[0, t_{\rm f}\right] \longrightarrow \mathbb{R}^7$ satisfying

● the transversality condition $p(t_{\rm f})\ \nabla A(x^*(t_{\rm f})),$

● the adjoint equation $-p'(t) = D_x \mathcal{H}(x^*(t), p(t), u^*(t))$,

● the maximum condition $\mathcal{H}(x^*(t), p(t), u^*(t)) = M(x^*(t), p(t))$ = $h$(constant), $a.e.$ for $t\in \left[0, t_{\rm f}\right]$.

Proof. See ^[19], p 514–519.

Now, we apply Pontryagin maximum principle (see Theorem 1) to give necessary conditions for any solution of the minimization problem (3.4).

Theorem 2. Suppose that $x^*$ and $u^*$ are optimal for the minimization problem (3.4). Then, $u^*$ satisfies

where

and

Here, $p$ is the so-called adjoint vector solution to the adjoint equation

and the terminal condition

Proof 1. According to Eq (4.3), we have the following Hamiltonian for system (3.4):

Let $p$ be the function defined in Theorem 1. It satisfies the adjoint equation (4.6)

The transversality condition satisfied by $p$, as defined in Theorem 1 is actually the terminal condition

The function $u \longmapsto \mathcal{H}(x^*(t), p(t), u)$ from $U$ to $\mathbb{R}$ attains its maximum at the point $u = u^*(t)$ (see Theorem 1).

In order to characterize the optimal control, we follow the technique proposed in ^[30] (see p 12). We differentiate the Hamiltonian with respect to $u_1, u_2, u_3, u_4$, $u_5$ and we solve $\dfrac{\partial \mathcal{H}}{\partial u_1} = 0$, $\dfrac{\partial \mathcal{H}}{\partial u_2} = 0$, $\dfrac{\partial \mathcal{H}}{\partial u_3} = 0$, $\dfrac{\partial \mathcal{H}}{\partial u_4} = 0$ and $\dfrac{\partial \mathcal{H}}{\partial u_5} = 0$ on the interior of the control set and then we use the property (boundness) of the control set $U$.

We obtain that:

$\textbf{a)}$ $\dfrac{\partial \mathcal{H}}{\partial{u_1}} = 0$ gives

When $0\leq \dfrac{\frac{N_v}{N_0}(p_5(t)-p_1(t))}{c_1 \varphi(S^*(t))} \leq 1$ the maximizing value is satisfied by taking $u^*_1(t) = \dfrac{\frac{N_v}{N_0}(p_5(t)-p_1(t))}{c_1 \varphi(S^*(t))}\in \left[0, 1\right]$. If $\dfrac{\frac{N_v}{N_0}(p_5(t)-p_1(t))}{c_1 \varphi(S^*(t))} > 1$, then $u^*_1(t) = 1$ (by projection). Similarly, $\dfrac{\frac{N_v}{N_0}(p_5(t)-p_1(t))}{c_1 \varphi(S^*(t))} < 0$ implies that $u^*_1(t) = 0$.

We note that $\varphi(S^*(t)) = 0$ if $S^*(t) = 0$, and in this case, there is no susceptible people to vaccinate, then $u^*_1(t) = 0$ and $u^{*}_{1}(t)\notin \left[0, 1\right]$. If $S^*(t) > 0$, then $\varphi(S^*(t)) > 0.$

Now, we can summarize (see again ^[19], p 444):

Which can be rewritten as:

We procede similarly for the other controls:

$\textbf{b)}$ $\dfrac{\partial \mathcal{H}}{\partial{u_2}} = 0$ implies that

Let us set $g(t) = \dfrac{(p_1(t)-p_2(t))S^*(t)(c_{\beta_I}\beta_II^*(t)+c_{\beta_D}\beta_DD^*(t))}{c_2}$

It follows the following characterization:

And then

$\textbf{c)}$ $\dfrac{\partial \mathcal{H}}{\partial{u_3}} = 0$ gives

Let us denote $C(t) = \dfrac{(p_1(t)-p_2(t))S^*(t)c_{\beta_H}\beta_HH^*(t)+(p_5(t)-p_6(t))c_\lambda \lambda H^*(t)}{c_3}$.Then, the restrictions on the control give:

and thus the following characterization holds:

$\textbf{d)}$ $\dfrac{\partial \mathcal{H}}{\partial{u_4}} = 0$ gives

Let us denote

Taking account the restrictions, we derive the following characterization:

Which becomes:

$\textbf{e)}$ Finally, $\dfrac{\partial \mathcal{H}}{\partial{u_5}} = 0$ gives

The characterization is given by:

which can be rewritten as:

The optimal control and the state are approximated by solving the optimality system, which is a combination of the state system (3.1), the adjoint system (4.7), the boundary conditions and the characterization of the optimal control Eqs (4.10)–(4.12)–(4.14)–(4.16)–(4.18).

By substituting the values of $u^*_1, u^*_2, u^*_3, u^*_4, u^*_5$ in the control system (3.1), we obtain the following optimality system (4.19)–(4.20).

5.   Existence of optimal solution

In this section, we prove that, under suitable assumptions, there exists an optimal solution for system (3.4).

To do so, we rewrite system (3.4) as (see ^[19], p 480):

where

and

$Q = \left[0, t_{\rm f} \right] \times \Omega.$

Then, we prove the following result:

Theorem 3. There exists an optimal solution to the problem (5.1).

Proof 2. To prove this result, we apply the Theorem 23.11 proposed in ^[19] (see p 481):

(a) We can see that each $G_i, i = 1, 2, ..., 5$ is continuous in $x$. Furthermore:

– $G_1(x) = \bigg(-\frac{N_v}{N_0}\varphi(S), 0, 0, 0, \frac{N_v}{N_0}\varphi(S), 0, 0\bigg)^T$. Thus, we have:

where $M_1 = 2\frac{N_v}{N_0}$

– $G_2(x) = \bigg(S(c_{\beta_I}\beta_II+c_{\beta_D}\beta_DD), -S(c_{\beta_I}\beta_II+c_{\beta_D}\beta_DD), 0, 0, 0, 0, 0 \bigg)^T$

where $k_1 = \max (c_{\beta_I}\beta_I, c_{\beta_D} \beta_D)$. Then, we have:

where $M_2 = 2k_1$

– $G_3(x) = \bigg (c_{\beta_H}\beta_HHS, -c_{\beta_H}\beta_HHS, 0, 0, c_\lambda \lambda H, -c_\lambda \lambda H, 0\bigg)$.

where $k_2 = \max (c_{\beta_H}\beta_H, c_\lambda \lambda)$. Thus

With $M_3 = 2k_2$

– $G_4(x) = \bigg (0, 0, -\eta_HI, \eta_HI+(\eta_R+\eta_D)H, -\eta_R H, -\eta_DH, 0 \bigg)^T$. We have:

where $k_3 = \max (2 \eta_H, 2\eta_R+2\eta_D)$. It follows that:

where $M_4 = 2k_3$

– $G_5(x) = \bigg (c_\tau \tau S, c_\tau \tau E, 0, 0, c_{\tau}\tau R, 0, -c_\tau \tau E)^T$. One has:

where $M_5 = 2c_\tau \tau$

Thus, we conclude that each $G_i(x), i = 1, 2, 3, 4, 5$ has linear growth.

(b) For almost every $t$, the set $U = [0, 1]^5$ is closed and convex by definition;

(c) The sets $Q$ and $\Omega$ defined above are closed;

(d) The running cost $\Lambda(x, u)$ is continuous in $x$ and $u$. This implies that $\Lambda$ is measurable in $x$ and continuous in $u$. Using Proposition 6.35 in ^[19] (p 123), we conclude that $\Lambda$ is LB measurable in $x$ and $u$;

Furthermore, the epigraph of $\Lambda$ is defined by

which is closed. Additionally, as $\Lambda$ is continuous, it follows that

– $\Lambda$ is lower semicontinuous in $(x, u)$;

– $u\longmapsto \Lambda(x, u)$ is convex for each $x \in \Omega$, (due to the convexity of $u_i^2, i = 1, 2, ..., 5$);

– $\Lambda(x, u) \geq K_4E+K_5D+K_6T, \forall x \in \Omega$. Indeed, $E$, $D$ and $T$ are bounded, so there exist $m_1, m_2, m_3 \geq 0$ such that $E \geq m_1$, $D \geq m_2$ and $T \geq m_3$. This implies that $\Lambda(x, u) \geq \lambda_0$, with $\lambda_0 = K_4m_1+K_5m_2+K_6m_3$.

(e) The projection $\{\alpha \in \Omega :(\alpha, \beta)\in \Omega\times \Omega \;{\mbox{for some}} \;\beta \in \Omega\}$ of $\Omega$ is bounded by definition.

(f) Let $r > 1$. By hypothesis, $\forall t\in \left[0, t_{\rm f}\right], u(t) \in U \Longrightarrow \left| u(t) \right| \leq 1$.

Let $k(t) = 1, \forall t$, we obtain that $\left| u(t) \right| \leq k(t)$ for a.e. $t$.

Thus, the hypotheses of Theorem 23.11 in ^[19] are satisfied and there exists a solution for problem (5.1).

6.   Uniqueness of the optimality system

Now, we aim to prove that the optimal controls given by Eqs (4.10), (4.12), (4.14), (4.16) and (4.18) are unique. To do so, following the ideas introduced by ^[31] (see p 435), we verify that the state and adjoint functions $f_i$ and $p_i$, $i = 1, 2, ..., 7$ are bounded and are Lipschitz functions. First, we rewrite the state equation on the following form:

where

where $x = (\mu+\lambda+\alpha-u_4(t)(\eta_H+\eta_D))$, $y = \alpha -u_4(t)\eta_R+c_\lambda u_3(t)\lambda$, $z = (1-c_\lambda u_3(t))\lambda -u_4(t)\eta_D$, $A_1 = -(\mu+(1-c_\tau u_5(t))\tau)$, $B_1 = -(\mu+\delta+(1-c_\tau u_5(t))\tau)$

$\phi = \bigg(S, E, I, H, R, D, T\bigg)^T$ and $\Phi_{t} = \bigg (\frac{{\rm d} S}{{\rm d}t}, \frac{{\rm d} E}{{\rm d}t},$ $\frac{{\rm d} I}{{\rm d}t}, \frac{{\rm d} H}{{\rm d}t}, \frac{{\rm d} R}{{\rm d}t}, \frac{{\rm d} D}{{\rm d}t}, \frac{{\rm d} T}{{\rm d}t} \bigg)^T,$

where $A_2 = S((1-c_{\beta_I}u_2(t))\beta_II+(1-c_{\beta_H}u_3(t))\beta_HH+(1-c_{\beta_D}u_2(t))\beta_DD)$.

One has

where $\Phi_1 = (S_1, E_1, I_1, H_1, R_1, D_1, T_1)^T$ and $\Phi_2 = (S_2, E_2, I_2, H_2, R_2, D_2, T_2)^T$.

Let $k_1 = \max (\beta_{I_1}, \beta_{I_2}), k_2 = \max (\beta_{H_1}, \beta_{H_2}), k_3 = \max (\beta_{D_1}, \beta_{D_2})$,

$k_4 = \max (c_{\beta_{I_1}}\beta_{I_1}\left| u_2\right|, c_{\beta_{I_2}}\beta_{I_2}\left| u_2\right|)$, $k_5 = \max (c_{\beta_{H_1}}\beta_{H_1}\left| u_3\right|, c_{\beta_{H_2}}\beta_{H_2}\left| u_3\right|)$ and

$k_6 = \max (c_{\beta_{D_1}}\beta_{D_1}\left| u_2\right|, c_{\beta_{D_2}}\beta_{D_2}\left| u_2\right|)$.

We have that, for $\Phi_1 \ne \Phi_2$:

Set $m_1 = \max (2k_1, 2k_4), m_2 = \max(2k_2, 2k_5), m_3 = \max(2k_3, 2k_6)$

By hypothesis, we have $I_1 \leq 1$, $S_2 \leq 1$, $H_1 \leq 1$ and $D_1 \leq 1$. Then,

with $M = 2m_1 +2m_2 +2m_3+2\dfrac{N_v}{N_0}$.

We set

and we note that

It follows that

Now, let

We have that $G$ verifies

Thus,

where

We conclude that $G$ is Lipschitz. Additionally, as the control variables $u_i(t), i = 1, 2, 3, 4, 5$ are bounded, we deduce that $G$ is bounded.

We follow the same techniques considering the adjoint system (4.7).

Let $Z(\Phi) = B \Phi + K(\Phi)$, where

where $b = (1-c_\tau u_5(t))\tau$, $e_1 = (\mu+(1-c_\tau u_5(t))\tau)$, $e_2 = \mu+\delta+(1-c_\tau u_5(t))\tau$, $e_3 = (\gamma +u_4(t)\eta_H$, $e_4 = \mu-c_{\lambda} u_3 \lambda+(\alpha +u_4(t))$, $z_1 = ((1-c_{\beta_I}u_2(t))\beta_II+(1-c_{\beta_H}u_3(t))\beta_HH+(1-c_{\beta_D}u_2(t))\beta_DD)$, $z_2 = \frac{N_v}{N_0}\varphi'(S)u_1(t)+c_1u_1(t)^2\varphi(S)\varphi'(S)$, $y = \mu+\lambda+\alpha -u_4(t)(\eta_R+\eta_D)$, $w = \mu+\gamma +u_4(t)\eta_H$, $v = (\lambda -u_4(t)\eta_D+c_\lambda u_3(t)\lambda)$, $(1-c_\tau u_5(t))$.

Now, we can write

From $\Phi_1 \ne \Phi_2$, we have that $\| \Phi_{1} - \Phi_{2} \| \ne 0$. Then, we obtain that

where

We have that $Z$ verifies:

where

Thus, the function $Z$ is Lipschitz. Since controls $u_i(t), i = 1, 2, 3, 4, 5$ are bounded, then we have the boundedness of $Z$.

We have shown that the state and adjoint functions are bounded and are Lipschitz functions, so the optimal solution is unique.

7.   Numerical simulations

In this section, we first propose a numerical approach to approximate the optimal solutions characterized in Theorem 2. Then, we propose numerical experiments to illustrate the interest of the proposed implementation.

7.1.   Numerical scheme

Here, we use the Runge Kutta fourth order method developed in ^[32] (see p 49–52) to solve the optimality system found in the previous sections. This method can be summarized as follows: (see again ^[32], p 50)

Step 1: Choose an arbitrary value of $u^*$ (generally $u^* = 0$) over $\left[0, t _{\rm f} \right]$;

Step 2: Using initial conditions $x(t_0) = x_0$ and the value of $u^*$ above, solve the state system with explicite scheme;

Step 3: Including the transversality conditions $p(t_{\rm f})$ and considering the expressions of $u^*$ and $x^*$ estimated previously, solve the adjoint system $p(t)$ with implicit scheme;

Step 4: Refresh the expression of $u^*$ by replacing $x(t)$ and $p(t)$ by their expressions.

Step 5: Check convergence. If values of the variables in the current and previous iterations are close enough, return the actual values as solutions. Else, return to Step 2.

We use the following notations for the state and adjoint equations:

Given a step size $h$, the approximation of each state variable $x_i, i = 1, 2, ..., 7$ is given by:

where

For the adjoint vector, the approximation is given backward in time and is of the form:

where

After these steps, the values of controls $u_i^*, i = 1, ..., 5$ are refreshed according to their expressions Eqs (4.10)–(4.18).

7.2.   Numerical experiments

In order to illustrate the pertinence of the controls, we implement the algorithm proposed previously. We consider several numerical results detailed below. To this aim, we considered a particular Scilab implementation of the numerical schemes presented in Section 7.1.

Based on the values proposed in ^[17] for the case of the 2014-2016 West-African Ebola virus disease epidemic, we consider the following data: $\beta_I = 0.2500; \beta_H = 0.0195; \beta_D = 0.2400; \mu = 0.021; \tau = 0.0000024;$ $\theta = 0.7500; \delta = 0.022; \gamma = 0.15; \lambda = 0.1177; \alpha = 0.1040; c_1 = 10;$ $c_2 = 20; c_3 = 30; c_4 = 40; c_5 = 50; c_\tau = 0.1;$ $c_\lambda = 0.2; c_{\beta_I} = 0.3; c_{\beta_D} = 0.3; c_{\beta_H} = 0.4; \eta_H = 0.6; \eta_R = 0.7;$ $\eta_D = 0.8; K_1 = 1000;$ $K_2 = 2000; K_3 = 8000; K_4 = 1200; K_5 = 500; K_6 = 1500; N_0 = 10000000;$ $N_v = 5000$ and $\Lambda = 0.0217$.

Furthermore, to obtain results with a relevant graphical interpretation, we set $S_0 = 0.999; E_0 = 0.001; I_0 = 0; H_0 = 0;$ $R_0 = 0; D_0 = 0; T_0 = 0,$.

We solve the state system with and without controls and compare the obtained results below.

In Figures 3–5, we show the evolution of susceptible ($S$) and exposed ($E$) populations, the evolution of infectious ($I+H$) and recovered ($R$) populations and the evolution of deaths ($D$) and cumulative number of infected persons leaving the system ($T$), respectively. On those graphs, the red curves represent the absence of controls in the considered population and the blue ones, the presence of controls.

We can see in Figure 3 that, in absence of controls, the number of susceptible people increases due te the effect of the recruitment term $\Lambda(t)$. On the other hand, if controls are applied, this number remains quite stable. Indeed, the vaccination campaign move people from this state to the recovered state, and thus allow to reduce the number of person that can be affected by the disease. Focusing of exposed people, we remark that, with controls, the number of infected persons per day decreases slowly, whereas without control measure it explodes dramatically. This show that the obtained optimal control measures seem to be efficient in reducing the impact of the epidemic.

Focusing on Figures 4 and 5, we observe that the number of infectious people ($I+H$) and deaths ($D$) decreases and remains quite low in presence of control measures. But if no controls are applied, those number increase rapidly. As expected, in case of control measures, the number of recovered people increases due to the sanitary measures in hospital and vaccination campaigns. When there is no control, this number tends to decrease with time.

Finally, we observe on Figure 5, that, when controls are applied, the cumulative number of infected persons leaving the country ($T$) is restrained, limiting the spread of the disease outside the affected territory. When no control measures are applied, as expected, this number growth quickly.

In Figures 6–8, we present the plots of the optimal controls. We observe an increasing behavior for all control measures. The model seems to indicate that focusing on detection campaigns ($u_4$) seems to be primordial as its slope increases faster that other measures $u_1$, $u_2$ and $u_3$. Then, vaccination campaign ($u_1$) and increasing sanitary measures in hospital ($u_3$) should start since the beginning of the hazard. The start of prevention campaigns seems to be delayed. Regarding the quarantine measures ($u_5$), they should be applied strongly since the beginning of the epidemic.

8.   Conclusions

Since the recent outbreaks of Ebola virus disease, that occurred in some Western African countries, several measures have been taken by sanitary authorities in order to control the disease. Additionally, with the latest cases detected in Guinea, Democratic Republic of Congo and Ivory Coast, more vigilance is required to avoid future outbreaks, and so, mathematical studies are useful tools to tackle this situation.

In this work, we established an optimal control problem of a simplified version of Be-CoDiS limited to the case of one country. The control variables corresponded to a wide range of possible control measures such as vaccination campaign, educational campaign, treatment of hospitalized persons, early detection campaign, and quarantine measures. We presented the existence and uniqueness of an optimal control. Using Pontryagin Maximal Principal, we gave necessary conditions for this optimality and obtained the characterization of optimal controls.

Finally, considering a representative numerical experiment, we implemented theoretical results and studied the pertinence of those controls. We see that the combined effect of the considered control measures seem to have a strong impact on the disease magnitudes.

Indeed, vaccination campaigns reduced the reservoir of possible affected people. The prevention campaigns allowed to decrease the contacts with infected individuals and those that died from Ebola and thus, to limit the propagation of the disease. The early detection campaigns and the increase of sanitary measures helped to reduce the number of deaths in a short time interval. Finally, the application of quarantine measure allowed to limit the propagation of the disease outside the affected area.

We note that we have considered in this work the propagation of the disease inside a single country. It might be interesting in future works to consider the interactions between two or more countries related by the migration flow.

Confilct of interest

The authors declare that they have no conflict of interest.

Acknowledgment

This work was carried out thanks to the financial support from the "Institut de Recherche et de Développement" through projects UMI 209-IRD and UMMISCO UGB and the Spanish "Ministry of Science and Innovation" under Project PID2019-106337GB-I00.