Optimal control problem for a tuberculosis model with multiple infectious compartments and time delays

a MAEGE Laboratory, FSJES Ain Sebaa, Hassan II University, Casablanca, Morocco b Laboratory of dynamical systems, Mathematical Engineering Team, Chouaib Doukkali University, El Jadida, Morocco c MACS Laboratory, FS Ain Chock, Hassan II University, Casablanca, Morocco d LAMS Laboratory, FS Ben M’Sik, Hassan II University, Casablanca, Morocco mohamed.elhia@gmail.com, balatif.maths@gmail.com, boujallal@gmail.com, m rachik@yahoo.fr


Introduction
Infectious diseases, the scourge of humanity, have marked the history of human societies. Over the centuries and throughout the world, they have always been the major cause of death. They generate considerable and sometimes unbearable socio-economic, demographic, cultural, health and safety costs.
In order to be able to deal with these devastating diseases, a relevant tool should be made available to facilitate public health decision-making.
In this context, mathematicians and epidemiologists have long worked together to create mathematical models that allow competent authorities to prepare in advance to react quickly and effectively if an epidemic breaks out. Since Daniel Bernoulli's famous work [1], mathematical models have become one of the most important tools used in fighting against epidemics. They help to understand the dynamics of infectious diseases and to estimate the effect of different control and prevention strategies. An interesting overview of the use of mathematical models *Corresponding Author in epidemiology can be found in [2][3][4][5][6]. The first modern model in mathematical epidemiology is the simple SIR model proposed by Kermack and McKendrick [7] to model plague and cholera epidemics. In this model, a population is divided into three compartments: susceptible (S) (those who are not infected but may become infected), infectious (I) (those who are infected and can transmit the disease) and recovered or cured individuals (R). As such, this model was not appropriate for many infectious diseases with a latency period such as influenza or tuberculosis, where a person is infected but not infectious. Therefore, an extension of the standard SIR model was proposed, including a fourth compartment (E) of exposed persons (those who are infected but not infectious). This extension has been called the SEIR model. (for more details on the compartmental models in epidemiology, we refer to [8] and [9]).
Tuberculosis or TB, as it is commonly called, is an old disease that has affected humans for thousands of years. The disease is caused by the Bacterium Mycobacterium Tuberculosis (MTB) that generally affects the lungs, but can also affect other parts of the body like the brain and spine. The TB bacteria are spread through the air from a person to another. When people with lung TB cough, spit, sneeze, speak, or sing, they spread the TB germs into the air [10]. People who have intense contact with a TB patient in poorly ventilated areas are the most likely to become infected. Tuberculosis is curable and preventable but may be fatal if not treated properly [11]. (for more details on tuberculosis, we refer the interested reader to [12][13][14][15]).
For several years TB represents a major global health problem. Millions of people continue to fall sick with TB each year worldwide. It is one of the top 10 causes of death worldwide, ranking above HIV/AIDS as one of the leading causes of death from an infectious disease. According to report by WHO [16], there were, in 2017, 10 million new TB cases and 1.6 million people died from the disease (including 300 000 deaths among people with HIV). There were cases in all countries and age groups, but overall 90% were adults in their most productive years (aged ≥ 15 years) and 87% of the world's cases were in developing countries.
The first mathematical model for TB dynamics was proposed by Waaler et al. [17]. The authors divided the population into three classes and constructed the model according to the epidemiological characteristics of the TB transmission. Recently, several researchers have proposed various dynamical models and developed a theoretical framework to understand TB transmission and control its spread (see, e.g., [18][19][20][21][22][23] and references cited therein).
Moualeu et al. [24] presented a deterministic model for the transmission dynamics of Tuberculosis in the context of weak diagnosis capacity. Optimal control theory is used to obtain a cost-effective balance of two different intervention methods. Huo and Zou [25] constructed a TB model with treatment at home and treatment in hospital. It lies emphasis on the modelling effects of treatment at home.
In [26] Li and Ma investigated the global dynamics of a TB model that considers the prevention effect and latent delay. Silva et al. [27] introduced delays in a TB model, representing the time delay on the diagnosis and commencement of treatment of individuals with active TB infection.
Kar and Mondal [28] presented a basic Tuberculosis model including exogenous re-infection, endogenous reactivation, and the re-infection among the treated individuals. The local stability analysis of the equilibrium is shown and an optimal control based on treatment strategy is solved using Pontryagin's maximum principle.
Yang et al. [29] considers in their TB model new and relapse infections. Using Lyapunov functions, it showed that the global dynamic is completely determined by the basic reproduction number R 0 .
Altaf Khan et al. [30] proposed a fractional order model for TB dynamics with relapse using Atangana Baleanu derivative. Kim et al. [31] proposes optimal control strategies for reducing the number of high-risk latent and infectious TB patients with minimum intervention implementation costs.
A comparison of some control strategies for a model with Caputo time fractional derivative is proposed by Yildiz [32]. Purwati et al. [33], in their model take into account a discrete agestructured population. The existence and stability of the model equilibrium are discussed based on the basic reproduction ratio. Then the optimal control strategy is applied for controlling the transmission of TB in child and adult populations. Using a basic SEI model with saturated incidence rate, Baba et al. [34] studied the effect of optimal controller and awareness. It is noted that a considerable number of studies that deal with TB modelling use the basic SEI or SEIR model with only a single compartment of infectious people (see [35][36][37][38][39][40][41]). However, considering a single infectious compartment is no longer convenient, as it fails to take into account that some infected individuals can be detected and treated whereas others remain undetected and untreated. Furthermore, according to the study in [42], patients diagnosed with tuberculosis who do not initiate treatment represent an important failing in the provision of care. For instance, the proportion of patients diagnosed with tuberculosis who experienced pre-treatment loss to follow-up in Africa ranged from 6 to 38%. It is also indicated, in this study, that some reasons for the loss to follow-up were Health-systemrelated obstacles, which included dissatisfaction with long waiting times in health services, the need for repeated visits and delays in receiving the results of Sputum Smears. Some reasons for not starting treatment for tuberculosis were patientrelated (e.g. difficulty getting time off from work or a lack of understanding of tuberculosis and its severity or the potential benefits of treatment). Other reasons were disease-related (e.g. weakness and fatigue). That is the reason why including the loss to follow-up phenomenon, detection and treatment into a TB model is too important.
Taking into consideration the works mentionned above and as an extension of a previous work [43] in which we have discussed the global dynamic of an SIR model with two stages of infection, we propose here an optimal control problem based on a more realistic SEIR model for TB that includes two infectious levels, multiple time delays and control terms. In our model, the first infectious level contains undetected infectious individuals, i.e. individuals who are infectious but have not yet performed the TB test. The second level involves individuals with Tuberculosis who accessed tests. As only a part of individuals who carried out a TB test are notified and start treatment while others are lost during the diagnostic process or did not initiate treatment, we consider a second infectious level which includes two subclasses of diagnosis patients, namely treated and lost to follow up individuals. Furthermore, it is well known that there is an incubation period for TB [44] and the treatment should follow a certain process to be efficient [45] and this matter justifies the use of two time delays representing the incubation period and the time needed for treatment in our model. To control the spread of TB, we consider three control terms which represent prevention, detection efforts and the efforts that prevent the failure of treatment. The aim of this work is to investigate the impact of prevention, detection and treatment strategies on reducing the number of exposed, undiagnosed and lost to follow-up individuals when the dynamics of the TB transmission is governed by a delayed model. To our knowledge, this work is the first to apply optimal control theory on a TB model which considers the loss to follow-up phenomenon, detection of infectious individuals and multiple time delays.
The paper is organised as follows: in Section 2, we present our mathematical model with time delays and control terms. In Section 3, we present some proprieties of the model. We formulate the optimal control problem in Section 4, prove the existence of a solution and put forward the control expression. In Section 5, we propose numerical simulations and discussion based on different control strategies. Finally, the conclusions are given in Section 6.

Mathematical model
The mathematical model we consider here is a deterministic compartmental model, composed of a delayed differential equation and three control terms.
At any time, an individual is in one of the six following compartments: Susceptible (S): healthy individuals who are not yet exposed to TB; Exposed (E): individuals who are in the latent period; Undiagnosed infectious (I): people who have active TB but have not been confirmed by a test; Treated infectious (T ): people who have been diagnosed as having active TB and follow a therapeutic program; Lost individuals (L): they accessed TB testing but are not under treatment, either because they were lost during the diagnostic process, did not receive any treatment or they did not complete their treatment; Recovered (R): they are individuals who were previously infected and then recover from active TB through natural recovery or after completing successfully their treatment.
We assume that all recruitment is into the susceptible compartment and occur at a constant rate Λ. The natural death rate, denoted µ, is constant across all compartments. TB is assumed to be fatal for infectious individuals, that is why we define additional death rates d I , d T and d L . The transmission of TB occurs following an adequate contact between a susceptible and infectious in (I), (T ) or (L). Due to the non-linear contact dynamics in the population, we use the incidence function β 1 I S N , β 2 T S N and β 3 L S N to indicate successful transmission of TB, where β i , i = 1, 2, 3, denote the effective contact rate with infectious individuals in compartment I, T and L respectively. Thus, new infections are produced at the rate All newly infected individuals develop latent TB and enter the compartment E. We also assume that the latent period of the disease is represented by the time delay τ 1 . People infected at time t−τ 1 become infectious at time t, they join the undetected compartment I at rate k (k is the rate at which individuals leave the latent class by becoming infectious). Among the undiagnosed infectious some of them are diagnosed and start their treatment at a rate δ 1 , while others are diagnosed but for some reason they do not follow any treatment, they enter the compartment L with a rate δ 2 . After a period of time suffering from TB, some individuals in the compartment L decide to go to the hospital with a rate γ 1 . Also among the infectious who had begun their treatment some of them will give up treatment and will enter the class L at a rate γ 2 . The patient under treatment might need time to recover, so we consider an other time delay τ 2 that represents the treatment duration. A patient who starts his treatment at time t − τ 2 will recover at time t. Lost and treated individuals progress to the recovered class with constant rates α 1 and α 2 respectively.
The control strategy we adopt consists of introducing three control parameters v 1 , v 2 and v 3 representing the following: v 1 : the efforts of preventing susceptible individuals from becoming infectious individuals. These efforts include awareness program, isolation and any other distancing measurement that can limit contact between susceptible and infectious people; v 2 : the efforts made to detect undiagnosed patients through screening tests and putting them under a therapeutic program. v 3 : the efforts that prevents the failure of treatment by providing financial aid to patients, consolidating the number of medical staff and increasing the monitoring of people at risk until they complete successfully their treatment.
The controls v 1 (t), v 2 (t) and v 3 (t), which are function of time t, are assumed to be bounded with 0 ≤ v i (t) ≤ 1 for i = 1, 2, 3. A flow chart of our model is given in Figure 1. The dynamics of the model is governed by the following system of delayed differential equations: where The initial conditions for system (1) take the form 3. Some proprieties of the model

Boundedness of trajectories
Let N the total population (N = S + E + I + T + L + R). The trajectories of the system (1) are bounded. Indeed, by adding all equations in (1), one has Thus, So, all possible solutions of the system (1) enter the region

Existence of solutions
Let The system (1) can be writing as follows The function F satisfies where M 1 , M 2 and M 3 are positive constants, independent of state variables and Moreover, one has Thus, it follows that the function G is uniformly Lipschitz continuous. From the boundedness of the controls v i , i = 1, 2, 3 and the restriction on the state variables, we conclude that there exist a solution of the system (1). (see [46])

The basic reproduction number
The basic reproduction number is the average number of secondary cases produced by a single infective individual which is introduced into an entirely susceptible population. Proposition 1. The basic reproduction number of the system (1) is given by where Proof. The system (1) has a disease-free equlibrium (DFE) D * = Λ µ , 0, 0, 0, 0 . In order to calculate the basic reproduction number, we use the Van Den Driesseche and Watmough next generation approach [47] and techniques reported in [48,49]. The next generation method consists of considering only the infected classes E, I, T and L.
We put where F is the rate of new infections in each class, and V describes the other flows in the infected classes.
The Jacobian matrix of F and V at the diseasefree equilibrium D * are given by The basic reproduction number is defined, according to Van den Driessche and Watmough [47], as the spectral radius of the next generation matrix, F V −1 . Since F V −1 is a rank one matrix, the only non-zero eigenvalue is given by (3).
is known explicitly, we can examine the sensitivity of R 0 with respect to the control terms and predict the relative change of R 0 with respect to each control. The sensitivity of R 0 with respect to a parameter ω is measured by the so called sensitivity index ∂R 0 ∂ω [50]. In this study, we get the following results: Proposition 2. The control v 1 decreases R 0 for any value of the system parameters, while the controls v 2 and v 3 can decrease or increase R 0 depending on the system parameters values.
Proof. Using the expression of R 0 given by (3) we have 1, 2, 3) and the terms A, B, C, D and E are positive one can easily check, from (4), that ∂R 0 ∂v 1 is negative for any value of the system parameters. So, we conclude that the increase in v 1 leads to the decrease in R 0 (v 1 , v 2 , v 3 ). Therefore, the control v 1 plays an important role in controlling the disease regardless of other system parameters.
In the expressions (5) and (6) we can see that the signs of ∂R 0 ∂v 2 and ∂R 0 ∂v 3 maybe positive or negative according to the values taken by the system parameters. So the impact of v 2 and v 3 on increasing or decreasing R 0 (v 1 , v 2 , v 3 ) is depending on the system parameters.

Optimal control problem
In this section we present our optimal control problem, discuss the existence of the optimal control and then give a characterization of optimal control terms. Note the state variable R does not appear in the first five equations of (1). So, the other variables do not depend on R, and we can limit our study to the following system The control strategy aims at minimizing the number of exposed (E), undetected (I) and lost individuals (L) as well as minimizing the cost of this strategy. Mathematically, for a fixed terminal time t f , the problem is to minimize the objective functional where A i ≥ 0( for i = 1, 2, 3), denote weights that balance the size of the terms. In other words, we seek the optimal values v * 1 , v * 2 and v * 3 of the controls v 1 , v 2 and v 3 , such that with U is the set of admissible controls defined by Remark 1. In the following, to avoid some mathematical complexities, we consider that the total population N remains constant during the control period.
In order to derive the necessary conditions for the optimal control, the Pontryagin's Maximum Principale with delay given in [51] is used. This principal converts the problem (7)-(9) into a problem of minimizing a Hamiltonian, H, defined by By applying Pontryagin's Maximum Principale with delay [51] and the existence result for optimal control and corresponding optimal states from the study [52], we obtain the following theorem: Theorem 1. Consider the optimal control problem (7)- (9). There exists an optimal control (v * 1 , v * 2 , v * 3 ) ∈ U and corresponding solutions S * , E * , I * , T * and L * such that Furthermore, there exists adjoint functions λ 1 ,λ 2 ,λ 3 , λ 4 and λ 5 , such that with the transversality conditions . . , 5 and j = 1, 2. Moreover, the optimal controls v * i ( for i = 1, 2, 3) are given by v * 1 (t) = min 1, max 0, (λ2 − λ1) λS N A1 v * 2 (t) = min 1, max 0, Proof. The existence of the optimal control is obtained from Fleming and Rishel [52] (see Corollary 4.1). The adjoint equations and transversality conditions can be obtained by using Pontryagin's Maximum Principle with delay in the state and control variables [51] such that The optimal controls v * i ( for i = 1, 2, 3) can be solved from the optimality conditions ∂H ∂v i = 0 ( for i = 1, 2, 3) (13) with further simplification of (13) and special attention on the bounds of controls as defined in U, we obtain (12).

Numerical simulation
In this section, we present the numerical solution of the optimality system which is a twopoint boundary value problem, with separated boundary conditions at times t 0 = 0 and t f . We use an iterative method based on a combination of forward and backward difference approximation, which converges when a tolerance criterion is reached. Taking into account the nature of the optimal control problem with time delays, we extend the iterative Forward-Backward Sweep Method (FBSM) in [53,54] to our delayed optimality system. First, we consider a step size h > 0, τ = max(τ 1 , τ 2 ) and (m, n) ∈ N 2 with τ = mh and t f − t 0 = nh. Then, we consider m knots to left of t 0 and right of t f and we obtain the following partition Next, approximations in term of nodal points of the state variables, the adjoint variables and the controls are given by the vectors X = respectively. The numerical resolution process is summarized in the following algorithm: Step 1 • Initialization on the left for state variables. For i = −m, . . . , 0, do • Initialization on the right for adjoint variables. For i = n, . . . , n + m, do Step 2 For i = 0, . . . , n − 1, • Solving the state system forward in time using the stored value for X and V . • Solving the adjoint system backward in time using the stored value for the controls and the state variable.
Step 3 Updating the controls using the formula (12).
Step 4 Testing the convergence : If the difference of values of these variables in this iteration and the last iteration is sufficiently small, output the obtained current values as solutions. If the difference is not considerably small, return to Step 2.
Note that as reported in [55], the risk of developing the disease subsequently is much higher in the 5 years following infection. So, for our simulations we fixed the total control period with t f = 5 years.
The parameter values are given in Table 1 and are expressed per year. These values are chosen to ensure that the reproduction number R 0 is greater than 1. Therefore, the disease will spread in the population, and we can show the effectiveness of our control strategies. Progression rate from I to T 0.3 assumed δ2 Progression rate from I to L 0.5 assumed γ1 Progression rate from L to T 0.2 [57] γ2 Progression rate from T to L 0.298 assumed α1 Recovery rate of L 0.25 [58] α2 Recovery rate of T 0.53 [59] Firstly, let us explore the effect of time delays on the dynamics of the system (7) without controls. It is well known that a proportion of people infected with TB bacteria develops active tuberculosis within a finite time, the latent period τ 1 is the range of 1 − 3 years [44]. Also, the current standard treatment for active TB consist of taking antibiotics for at least six to nine months [45], so we have considered two values of τ 2 (τ 2 = 0.5 or 0.75 year).
In Figure 2, we present the evolution of the compartments E, I and L using different values of the time delay τ 1 and fixing the value of τ 2 at 0.5 year. We observe that introducing time delays into our model can have a profound effect on the number of the infected; the smaller the time delay is, the higher the number of infected individuals in E, I and L. Note that for the case where τ 2 is set to 0.75 year and τ 1 is changing we obtained the same results as in Figure 2, which can be explained by the small difference between τ 2 = 0.5year and τ 2 = 0.75 year.
In the case where τ 1 is fixed (e.g τ 1 = 1 year) and τ 2 varies, we see in Figure 3 that we obtain the same solution for E, I and L, due to the small variation in τ 2 . However, the difference when varying τ 2 is seen in the number of the treated people and a slight difference is observed at the end of the period (see Figure 3 (d)).
Moreover, it is shown from the Figures 2 and 3 that the number of exposed, infectious and lost individuals increased significantly with and without time delays. So, there is a real need to find an effective strategy to control the spread of TB. (a) Exposed.      In the rest of this section, let us chose τ 1 = 1 year and τ 2 = 0.5 year, and investigate the effectiveness of control strategies to achieve several objectives. Objective 1: Preventing susceptible individuals from becoming infectious individuals.
Given the major role of contact in transmitting TB between susceptible and infectious people and the importance of prevention programs in limiting the number of new cases, we propose an optimal strategy based on the control v 1 . Recall that v 1 represents any measurement that can reduce contact between susceptible and infectious individuals, such as awareness programs, distancing, or isolation. In Figure 4, we depict the evolution of the infected individuals in E, I and L compartments over time in the uncontrolled case and when the control v 1 is implemented alone. We observe that the number of infected in E, I and L rose dramatically over the given period when no control is exerted. While, in the presence of the preventing control v 1 , we can see that the number of individuals decreases sharply, and shows a decline of 74%, 66% and 47% for E, I and L respectively in the final year, compared to no-control case.  Objective 2: Encouraging detection and therapeutic programs.
To achieve this goal, we propose the control v 2 which includes screening and any other effort that facilitates access to therapeutic programs. Figure  5 displays the optimal solution of T , E, I and L when only the control v 2 is considered. Compared to the uncontrolled case, there is a steady increase in the number of patients under treatment. Over the control period, the average number of T increases by 14%. It is interesting to note that the use of v 2 alone also has an impact on the improvement of the results obtained for E, I and L, there is a reduction by 6%, 23% and 16% at the end of the period, respectively.   Objective 3: Offering support and followup for patients under treatment.
It should be remembered that in some countries, especially in Sub-Saharan Africa, some infectious cases that were detected are lost during the diagnosis process or after having started their treatment for financial, cultural and psychological reasons. To deal with this situation, we propose an optimal control strategy using the control v 3 that represents measurements which can prevent the loss phenomenon (like financial support and monitoring). Figure 6 shows that the control v 3 reduces the number of L and I by 16% and 10% respectively at the end of the period. While its impact on reducing the number of exposed individuals is minimal (3%).    Figure 7 represents the optimal controls when each one is used alone. In this figure we can see that full effort on v 1 is applied during almost the entire control period. The control v 2 is also fully used but only for the first 2.5 years and then decreased smoothly. We also observe that the control v 3 is the least used.   At this point, it can be concluded that in terms of reducing E, I and L simultaneously, the prevention control v 1 is the most efficient when only one control is used. However, combining controls is potentially more effective than using any one alone. Different combinations of optimal controls and strategies (v 1 with v 2 , v 2 with v 3 , . . .) can be used to achieve other objectives depending on the particularity of the disease in each country.
Let us now find the most effective combinations among the three controls. In Table 2 we present the number of infected individuals E +I +L at the end of the given period, when different interventions are combined. According to these results, the most efficient choices are those including the prevention control v 1 , namely is less efficient than the other strategies.
As an illustrative example, we depict in Figure  8 the evolution of E, I and L when the three controls are implemented. There is a clear improvement in the results obtained. In the final year, the number of individuals in E, I and L are reduced by 75%, 71% and 60%, respectively.  Overall, the numerical results under various interventions show that the prevention control v 1 is the most useful control. If we have to choose a strategy based on a combination of controls, then this one must include v 1 which contributes positively to the improvement of the final results.

Conclusion
The work presented here contributes to the growing literature on applying control techniques to epidemiology. Our contribution consists of formulating and solving an optimal control based on a more realistic SEIR model for TB, which includes time delays, three infectious compartments and three controls. It is well known that TB's symptoms appear awhile after infection and patients must follow a therapeutic protocol for a certain time. Thus, involving delays representing the latent period and the treatment duration in our model is of particular importance in TB's dynamics modelling; it allows having more complicated dynamics and a more consistent model with the real situation. Although some patients are diagnosed and treated, it turns out that others, despite being screened, are lost to follow-up for cultural, socio-economic and health system-related reasons. Thus, instead of considering a single compartment which groups all the infectious patients, it is more appropriate to divide the infectious people into three categories, namely, undiagnosed infectious, diagnosed patients who are under treatment, and finally, diagnosed patients who are lost to follow-up. Such a classification best reflects the role of screening and treatment and the importance of the loss to follow-up phenomena and provides a better understanding of TB transmission.
To control the spread of TB, three different control strategies were considered, namely, prevention efforts (v 1 ) (like awareness program and isolation), detection efforts (v 2 ) and efforts that prevent the failure of treatment (v 3 ) (e.g. financial support and monitoring). By introducing these controls in our model, we aim to achieve three objectives: • Preventing susceptible individuals from becoming infectious individuals. • Encouraging detection and therapeutic programs. • Offering support and follow-up for patients under treatment.
It should be noted here that the control system for TB with delays, which we propose in this paper, can easily be adapted for other diseases that were previously modelled using the standard SEIR model. Using the next-generation matrix method, we got the expression of the basic reproduction number R 0 with controls. The sensitivity analysis of R 0 with respect to the terms of the control shows that the v 1 plays an important role in TB control. In other words, we found that an increase of v 1 resulted in a decrease in R 0 regardless of system parameters; while the impact of v 2 and v 3 on the R 0 increase or decrease is depending on the values taken by the other system parameters.
To better understand the impact of delays on our system without control, we proposed numerical simulations by setting τ 1 and changing τ 2 and vice versa. We found that τ 1 has more influence on the dynamics of our model whereas the impact of τ 2 is limited to a slight change in the number of people who are under treatment, this is due to the little difference that exists between the two values of τ 2 that we considered (0.5 year and 0.75 year). Also, we found that the smaller the τ 1 and τ 2 values are, the higher the number of infected people is.
The main concern of this study is to investigate the impact of the proposed strategies on reducing the number of infected people. For this purpose, we have designed an optimal control problem in which we seek to minimize the number of exposed E, undetected I and lost individuals L as well as the cost associated to the implementation of the control measures. To investigate the effectiveness of the control strategies, several cases associated with the three control strategies were considered. Single control and different combinations were compared in this paper.
Numerical results show that when only one control is to be applied, the best choice is the prevention control. This control that reduces contact between susceptible and infectious people has a great effect on the number of individuals in E, I and L compartment. In contrast, if we use v 2 or v 3 alone, we observe that despite there is a decline in the number of people in I and L compartment, their effect on reducing the number of exposed individuals is minimal. When a combination of controls is to be employed, our results suggest that it should include the control v 1 . In this case, the reduction of E + I + L at the final time t f range from 72 to 73%. For future work, it would be interesting to extend the structure of our model by incorporating age effect and/or spatial diffusion, and observe how this can affect the optimal dynamics of our model. In addition, we intend to propose a direct method to find optimal control based on the viability theory.