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

1 Laboratory of Numerical Analysis and Computer Science, Applied Mathematics Section, Gaston Berger University , Saint-Louis, 209-IRD & UMMISCO-UGB, Senegal 2 Mathematics and Applications Laboratory, Mathematics Department, Assane Seck University, Bp: 523, Ziguinchor, 209-IRD & UMMISCO-UGB, Senegal 3 Interdisciplinary Mathematics Institute, Department of Applied Mathematics and Mathematical Analysis, Complutense University of Madrid, Plaza de Ciencias, 3, 28040 Madrid, Spain


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.

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] 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 • Λ ∈ R + is the recruitment rate of persons in state S (person.day −1 ), • µ ∈ [0, 1] is the mortality rate (day −1 ), • β I , β H , β D ∈ R + are the disease effective contact rates (day −1 .person −1 ) of people in compartment I, H and D respectively, • δ, γ ∈ R + denote the transition rates (day −1 ) from compartment E to I and I to H, respectively. • λ ∈ [0, 1] is the disease fatality percentage times the transition rate from compartment H to compartment D. • α ∈ 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.
• θ ∈ R + is the burial rate (day −1 ) of people that died from Ebola.
• τ ∈ [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.

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,t ≥ 0. The controls and their associated parameters are defined as follows: • u 1 (t) ∈ [0, 1] is the control corresponding to the vaccination campaign. First, this control is multiplied by N v N 0 , which corresponds to a reasonable estimation of the maximum capacity of the vaccination campaign, where N v ∈ N is the maximum number of persons that can be vaccinated per day and N 0 ∈ 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 φ ∈ C ∞ (R). We point out that φ(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 φ is a filter used to avoid negative values of S (notice that S (t) = 0 implies that dS dt (t) ≥ 0). Additionally, this function does not alter the maximum capacity 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) ∈ [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 β I and β D by a maximum of c β I β I and c β D β D , respectively. Here, c β I , c β D ∈ [0, 1] are the maximum percentage of reduction for β I and β D respectively. are the maximum percentage of reduction for β H and λ respectively. Furthermore, u 3 is also used to increase the transition rate from compartment H to compartment R, through the term c λ u 3 (t)H(t) (corresponding to the increase of the survival rate). • u 4 (t) ∈ [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)η H I(t). Here, η 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)η R H(t) and −u 4 (t)η D H(t), respectively. Here, η R (day −1 ) and η 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, η H can be estimated taking into account that 1 γ + η 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 η D and η R ). • u 5 (t) ∈ [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 spread-ing the disease outside the considered country. From a modeling point of view, u 5 is used to decrease the value of τ (i.e., the movement of people leaving the country). More precisely, we consider the parameter c τ ∈ [0, 1] which corresponds to the maximum percentage of reduction for τ that can be reached due to the application of those control measures. Then, τ is multiplied by (1 − c τ 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: x with a suitable function f : given by: Considering those controls, we are interested in minimizing the following cost function: where t 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]): Following [19], any pair (x, u), with u(t) ∈ U, for all t ∈ [0, t f ] 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 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 Λ are continuous and admit derivatives relative to x, denoted by D x f (x, u) and D x Λ(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).

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 ϵ > 0 such that, for any other process (x, u) satisfying the constraints of minimization problem (3.4), as well as ∥x − x * ∥ ∞ ≤ ϵ * , we have that J(x * , u * ) ≤ J(x, u). In this terminology, u * is called an optimal control and x * is called an optimal trajectory.
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 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 ≤ t ≤ t 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 ∈ R 7 and x ∈ R 7 , we denote • the transversality condition p(t f ) ∇A(x * (t f )), 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 and 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) (4.7) The transversality condition satisfied by p, as defined in Theorem 1 is actually the terminal condition The function u −→ H(x * (t), p(t), u) from U to 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).
Taking account the restrictions, we derive the following characterization: Which becomes: c τ τ p 1 (t)S * (t) + (p 2 (t) − p 7 (t))E * (t) + p 3 (t)R * (t) 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).

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

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: . Thus, we have: where k 1 = max(c β I β I , c β D β D ). Then, we have: We have: . It follows that: One has: 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 Ω defined above are closed; (d) The running cost Λ(x, u) is continuous in x and u. This implies that Λ is measurable in x and continuous in u. Using Proposition 6.35 in [19] (p 123), we conclude that Λ is LB measurable in x and u; Furthermore, the epigraph of Λ is defined by which is closed. Additionally, as Λ is continuous, it follows that -Λ is lower semicontinuous in (x, u); u −→ Λ(x, u) is convex for each x ∈ Ω, (due to the convexity of u 2 i , i = 1, 2, ..., 5); Indeed, E, D and T are bounded, so there exist m 1 , m 2 , m 3 ≥ 0 such that E ≥ m 1 , D ≥ m 2 and T ≥ m 3 . This implies that Λ(x, u) ≥ λ 0 , with λ 0 = K 4 m 1 + K 5 m 2 + K 6 m 3 .
Thus, the hypotheses of Theorem 23.11 in [19] are satisfied and there exists a solution for problem (5.1).

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: We have that, for Φ 1 Φ 2 : By hypothesis, we have I 1 ≤ 1, S 2 ≤ 1, H 1 ≤ 1 and D 1 ≤ 1. Then, We set and we note that We have that G verifies Thus, 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).
. Now, we can write Then, we obtain that We have that Z verifies: 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.

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.

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 [0, t f ]; 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 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.
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: After these steps, the values of controls u * i , i = 1, ..., 5 are refreshed according to their expressions Eqs (4.10)-(4.18).

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 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 Λ(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.

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.