AN OPTIMAL CHEMOPROPHYLAXIS AND TREATMENT CONTROL FOR A SPATIOTEMPORAL TUBERCULOSIS MODEL

unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract. In this work, we propose a Susceptible-Exposed-Infected-Recovered (SEIR) spatiotemporal model that characterizes the dynamics of tuberculosis disease by taking into consideration the spatial heterogeneity; in order to provide a realistic description of this disease. Based on an existing model, we add the Laplacian term in each class to describe the spatial mobility of its individuals, which led us to a SEIR reaction-diffusion system. Then, controls with treatment and chemoprophylaxis are incorporated to reduce the latently infected (exposed) and actively infected individual populations to fight against the spread of the disease. Theoretically the existence, positivity and boundness of state systems have been proved, also the existence of controls has been shown, and a characterization of the controls in terms of the state functions and the adjoint functions has been provided. To illustrate the effectiveness of our theoretical results, we give numerical simulations for several scenarios. Our results indicate that the control effect is effective if both control strategies controls on treatment and chemoprophylaxis strategies are used simultaneously.


INTRODUCTION
Tuberculosis (abbreviated TB) is a common infectious disease caused by Mycobacterium tuberculosis bacteria (Mtb). It mainly attacks the lungs (for pulmonary tuberculosis). Moreover, tuberculosis could affect the central nervous system, the circulatory system, the urinary-genital system, the bones, the joints, and even the skin. Tuberculosis could be spread through coughing, sneezing, kissing, through spitting of people who have pulmonary tuberculosis, as well as the use of non-sterilized utensils (plates, water glasses) of an infected person. In rare cases, a pregnant woman with active tuberculosis can infect the fetus ( [22]). On the one hand, this disease can only be spread through people who are infected with the active tuberculosis virus; on the other hand, the ones with an inactive tuberculosis virus cannot infect others. The spread of the virus from an infected person to another depending on the number of droplets contaminated and expelled by the patient, the aeration of the environment, the exposure time to the contaminating virus and the virulence of the Mycobacterium tuberculosis (Mtb) strain. Stopping the spreading of this disease could be done by isolating (putting in quarantine) the infected patients and by the effective commencement of the anti-tuberculosis treatment ( [4,11,29,22]). Lately, about 95% of 8 million new (TB) cases each year are recorded in the developing countries, in which 80 of the cases are recorded on people aged 19 to 59 years old ( [22]).
As a matter of fact, mathematical models, especially compartmental models, have played an essential role in fighting against infectious diseases since their birth by Kermack 1927. Several infectious diseases have been modeled using compartmental mathematical models, particularly tuberculosis ( [5,14,15,30,6,13,8]), these works have modeled the dynamics of this disease in several regions of the world. In the work ( [1]) authors give a SEIR model based on some previous works, in their model they consider the aspects of exogenous reinfection, disease relapse as well as primary active (TB) cases and natural recovery, and they incorporate chemoprophylaxis.
The results already mentioned neglected the spatial effect on the spread of the disease, statistics in the world ( [16,27]) and particularly in Morocco ( [25]) have shown that there are high-risk and low-risk regions, from this remark and based on previous works ( [32,9,17,18,26]), where authors have introduced models that take into consideration the spatial effect on the spread of diseases using partial differential equations, we added to the model of ( [3]) terms of spatial diffusion in each equation in order to express the mobility of individuals in high-risk region.
Several types of interventions exist, vaccines and treatments, one treatment given to latent people (chemoprophylaxis) and another concerns infectious individuals. In this work, we consider (spatiotemporal) optimal control strategy associated with chemoprophylaxis and treatment of latently and actively infected individuals with (TB). Also, many researches have focused on this topic and other related topics [33,34]. The theory of optimal control is a branch of mathematics used to act and control the dynamical systems, in the literature several works applied this theory to fight against tuberculosis using ODE models ( [5,12,14,19,20,15,8]). In this work, we will extend this method to a PDE's model using the variational approach. Optimal control techniques are used to characterize a pair of optimal controls in terms of state and adjoint functions. The optimality system is solved numerically based on a discrete iterative scheme that converges following an appropriate test related to the forward-backward sweep method (FBSM) ( [21]). Numerical simulations show that the control effect is effective if chemoprophylaxis and treatment strategies are used simultaneously.
The structure of the present paper is as follows : Section 2 is devoted to the basic mathematical model and the associated optimal control problem. In Section 3, we prove the existence of a strong global solution for our system. In Section 4, we prove the existence of an optimal solution. Necessary optimality conditions are established in Section 5. As an application, the numerical results associated with our control problem are given in section 6. Finally, we conclude the paper in Section 7.

THE BASIC MATHEMATICAL MODEL
2.1. The model without controls. We shall consider the model of tuberculosis (TB) proposed by C.P. Bhunu et al. ([3]), in this model the human population is divided into four compartments : susceptible individuals S T , those exposed to (TB) E T (latently infected), those infected and displaying symptoms of (TB) I T , and those who have recovered from sickness, R T . It is assumed that susceptible humans are recruited into the population at per capita rate Λ. The total variable population size at time t is given by, N(t) = S T (t) + E T (t) + I T (t) + R T (t). Susceptible individuals acquire (TB) infection following contact with an active infectious individual at a rate λ = β cI T N , β is the probability that one susceptible individual becomes infected by an infectious individual, and c is the per capita contact rate. Susceptible individuals infected are moved to latently infected class at a rate f λ , where f is the probability that an infected enters the latent stage. The latently infected individuals move to active (TB) at rates k for endogenous reactivation and δ 1 λ for exogenous re-infection respectively. Susceptible individuals infected with (Mtb) progress to the infective class at a rate of (1 − f )λ and represent the main active cases (TB). When an individual is in the active stage of the disease, he may recover normally at rate p and then moves in the recovered class R T (although they may contain certain live bacilli).
Individuals in R T are not fully immunized against (Mtb) infection and become infected at a rate δ 2 λ and progress to E T , as primary infection provides certain immunity. Certain individuals from R T become infected again at a rate of q. µ is the normal mortality rate in every class supposed to be identical, there is a death because of disease, therefore individuals infected have died at a rate d. Rates of treatment for individuals with latent infections and for the infectious are supposed to be r 1 and r 2 , respectively. The basic model of (TB) takes the following form: As mentioned above in the introduction, statistics showed that in several regions in the world, the spatial factor plays a major role in the transmission and the propagation of the disease, this lead us to think to introduce in our model (2.1) terms that express the spatial effect. Based on several previous works who have adopted this approach in several similar cases ([32, 9, 17, 18, 26]) using PDE's, we choose to integrate the ordinary spatial diffusion term in each compartments S T , E T , I T , R T , in order to describe the mobility of the population. Let Ω be a fixed and bounded domain in R 2 with smooth boundary ∂ Ω and η is the outward unit normal vector on the boundary. The time belongs to a finite interval [0, T ], while x varies in a bounded domain Ω, the system (2.1) becomes: The positive constants α 1 , α 2 , α 3 , and α 4 denote the corresponding diffusion for the susceptible population, exposed population, infected population and recovered population, respectively.
Initial conditions and no-flux boundary conditions are given by With the aim of demonstrating the impact of spatial factor, and mobility contribution to tuberculosis transmission, we present a simulation of our model over a 5-year period ( [15,10,28]), for further information on the values of the parameters and numerical approach see paragraph 6 named (Numerical results). Figures 1, 2, 3, and 4 present the numerical results for the numbers of susceptible, exposed, infected, and recovered people, respectively. We consider two situations, in the first one the disease starts from the corner (1), and second one, the disease starts from the middle (2). As can be seen, in both situations, individuals who are susceptible will be exposed and after a period of incubation get infected, therefore spread the disease to reach the whole population, however, the only difference between these two cases (1) and (2) is that the disease spreads rapidly in the second case, which indicates the danger of the disease when it starts in the middle, which illustrates the significance of applying a spatial approach. Regarding the class of recovered, there is a small number of recovered persons ( approximately 5 individuals ).
The observations obtained in these simulations lead us to propose an adapted control strategy taking these observations into account. The strategy proposed here is to introduce two controls, the first being chemoprophylaxis to minimize the individual exposed (latently infected) and the second being incorporated treatment to minimize the individual actively infected. 2.2. The model with controls. As mentioned in the last paragraph, the number of latently infectious and infected people increases considerably, motivated by the work of ([1]), a strategy of control is being introduced that involves two types of treatments, with the first ζ 1 (t, x) representing the r 1 chemoprophylaxis effort of individuals who are latently infected, in order to minimize the number of individuals who are infected. Whereas the control ζ 2 (t, x) represents the treatment effort r 2 of individuals actively infected, to augment more recovered individuals.
Under these new hypotheses, we will obtain r 1 ζ 1 individuals who will move from the class E T  (latently infected) to the recovered class R T , and r 2 ζ 2 will move to the class R T from the class of infected I T . Using the newly added changes, then our system with controls becomes as follows.
Initial conditions and no-flux boundary conditions are given by The boundary conditions imply that the population does not diffuse across the boundary and we define our objective functional as Subject to the state system given by (2.5-2.7). The objective of our work is to minimize the infected population and the cost of implementing the control by using possible minimal control variables ζ i for i = 1, 2. In the objective functional the quantity ρ 1 and ρ 2 represents the weight constant of the latently infected population and actively infected population, respectively. η 1 and η 2 are weight constants for mechanisms on chemoprophylaxis control and treatment control, respectively. The terms describe, the costs associated to the mechanisms on chemoprophylaxis control and treatment control, respectively. The square of the controls variables reflects the severity of the side effects of the mechanisms on chemoprophylaxis and treatment. Our objective is to find control functions such that , where the control set is defined as For biological reasons, the following are assumed to hold:

EXISTENCE OF SOLUTION
As the model (2.5-2.7) describes the population for biological reasons, the population S T , E T , I T , and R T should remain non-negative and bounded. We study in this section the existence of a global strong solution, positivity, and boundedness of solutions of problem for (2.5-2.7). Let y = A denotes the linear operator defined as follows With the domain of A is defined by Let Ω be a bounded domain from T R 2 , with the boundary smooth enough, y 0 i ≥ 0 on Ω (for i = 1, 2, 3, 4), the problem (2.5-2.7 ) has a unique (global) strong solution y ∈ W 1,2 ([0, T ] : H (Ω )) such that y i ∈ L (T, Ω )∩L ∞ (Q) for i = 1, 2, 3, 4 . In addition y 1 , y 2 , y 3 and y 4 are non-negative. Furthermore there exists C > 0 (independent of (ζ 1 , ζ 2 )) for all t ∈ [0, T ] Proof. To prove the existence of a (global) strong solution for system (2.5-2.7), now we write system (2.5-2.7) as shown in (((7.1)) see Appendix). Let The system (3.4) represent the nonlinear term of (2.5) and we consider the function h (y (t)) = (h 1 (y (t)) , h 2 (y (t)) , h 3 (y (t)) , h 4 (y (t))), then we can rewrite the system (2.5-2.7) in H(Ω) as follows Notice that the function h(t, y) defined in (3.5) are not Lipschitz continuous with respect to y, uniformly for t ∈ [0, T ]. Therefore, we cannot apply Theorem (7.1) (see appendix) for our problem directly.
Step 1: This step studies the local existence of positive solutions to system (2.5)-(2.7) in view of Theorem (7.1) (see appendix). We use a truncation procedure for h. For a fixed positive integer k > 0, let us define the function sets and consider the following auxiliary problem: where h k t, y k = h k 1 t, y k , h k 2 t, y k , h k 3 t, y k , h k 4 t, y k . Here, for each index i, h k i t, y k are defined as follows where [y i ] Ds i means that y i ∈ Ds i , and As the operator A defined in (3.1-3.2) is dissipating, self-adjoint and generates a C 0 -semigroup of contractions on H (Ω )( [31]), it is clear that function h k t, y k becomes Lipschitz continuous in y k uniformly with respect to t ∈ [0, T ]. Therefore, theorem (7.1) (see appendix) assures problem (2.5-2.7) admits a unique strong solution In order to show that Note that this system has a solution given by and analogously, we have Thus we have proved that By the first equation of (2.5), we obtain Using the regularity of y k 1 and the Green's formula, we can write Since y k i L ∞ (Q) for i = 1, 2, 3, 4 are bounded independently of (ζ 1 , ζ 2 ) and y 0 1 ∈ H 2 (Ω ), we deduce that We make use of (3.6), (3.9), and (3.10), in order to get and concluding that the inequality in (3.3) holds for i = 1, similarly for y k 2 , y k 3 and y k 4 . In ordre to show the positiveness of y k i for i = 1, 2, 3, 4, we write the system (2.5) in the form It is obvious to see that the functions F 1 y k 1 , y k 2 , y k 3 , y k 4 , F 2 y k 1 , y k 2 , y k 3 , y k 4 , F 3 y k 1 , y k 2 , y k 3 , y k 4 , and F 4 y k 1 , y k 2 , y k 3 , y k 4 are continuously differentiable satisfying for all y k 1 , y k 2 , y k 3 , y k 4 ≥ 0. Since initial values of system (2.5) are non-negative, we deduce that y k 1 (t, x) ≥ 0, y k 2 (t, x) ≥ 0, y k 3 (t, x) ≥ 0, and y k 1 (t, x) ≥ 0, ∀ (t, x) ∈ Q (see ( [24])). Now we particularize k > 0 large enough such that For example, we can take k > 2max y 0 i L ∞ (Ω ) , i = 1, 2, 3, 4 . Let θ ∈ (0, T ) be maximal with property (3.12). By (3.8)-(3.12), it is clear that y k i (t, x) < k , for (t, x) ∈ [0, θ ] × Ω and i = 1, 2, 3, 4. So, h k (t, y 1 , y 2 , y 3 , y 4 ) coincides with h(t, y 1 , y 2 , y 3 , To this end, we first introduce P = y 1 + y 2 + y 3 + y 4 then . This leads to the estimate 0 < P (t, x) ≤

THE EXISTENCE OF THE OPTIMAL SOLUTION
In this section, we will prove the existence of an optimal control for the problem (2.8) subject to reaction diffusion system (2.5-2.7) and (ζ 1 , ζ 2 ) ∈ U ad . The main result of this section is the following theorem. admits an optimal solution (y * , (ζ * 1 , ζ * 2 )).
Analogously, we have for y n i −→ y * i in L 2 (Ω) for i = 2 , 3 , 4 uniformly with respect to t. From the boundedness of ∆y n i in L 2 (Q), which implies it is weakly convergent in L 2 (Q) on a subsequence denoted again y n i then for all distribution ϕ y n i → y * i weakly in L 2 0, T ; H 2 (Ω) , i = 1, 2, 3, 4 y n i → y * i weakly star in L ∞ 0, T ; H 1 (Ω) , i = 1, 2, 3, 4 We now show that y n i y n j → y * i y * j for i = 1, 2, 3, 4 and j = 1, 2, 3, 4 strongly in L 2 (Q), we write y n i y n j − y * i y * j = (y n i − y * i ) y n j + y * i y n j − y * j and we make use of the convergences y n i −→ y * i strongly in L 2 (Q), i = 1, 3, 3, 4, y n j −→ y * j strongly in L 2 (Q), j = 1, 3, 3, 4 and of the boundedness of y * i , y n j in L ∞ (Q), then y n i y n j → y * i y * j strongly in L 2 (Q). We use 0 < N n and 0 < N * , and of the boundedness of N * , N n in L ∞ (Q), we deduce that y n i y n j N n → y * i y * j N * for i = 1, 2, 3, 4 and j = 1, 2, 3, 4.
Since ζ n 1 and ζ n 2 are bounded, we can assume that ζ n 1 → ζ * 1 and ζ n 2 → ζ * 2 weakly in L 2 (Q) on a subsequence denoted again ζ n 1 and ζ n 2 . Since U ad is a closed and convex set in L 2 (Q), it is weakly closed, so (ζ * 1 , ζ * 2 ) ∈ U ad We now show that and making use of the convergences y n i −→ y * i strongly in L 2 (Q) f or i = 1, 2, i = 1, 2 and ζ n i −→ ζ * i weakly in L 2 (Q), f or i = 1, 2, one obtains that ζ n i y n i+1 → ζ * i y * i+1 weakly in L 2 (Q) f or i = 1, 2.

NECESSARY OPTIMALITY CONDITIONS
 ∈ U ad , in this section, we show the optimality conditions to problem (2.5-2.7), and we find the characterization of optimal control. First , we need the Gateaux differentiability of the mapping ζ → y(ζ ). For this reason, denoting by y ε = y ε 1 , y ε 2 , y ε 3 , y ε 4 = (y 1 , y 2 , y 3 , y 4 ) (ζ ε ) and y * = y * 1 , y * 2 , y * 3 , y * 4 = (y 1 , y 2 , y 3 , y 4 ) (ζ * ) the solution of (2.5-2.7) corresponding to ζ ε and ζ * respectively, where (y * , ζ * ) is an optimal pair such that ζ ε = ζ * + εζ ∈ U ad (for ε > 0 small) and ζ ∈ L 2 (Q) 2 . We put θ = β c N , , we have the following theorem for i = 1, 2, 3, 4 and we denote S ε the system (2.5-2.7) corresponding to ζ ε and S * the system (2.5-2.7) corresponding to ζ * , subtracting system S ε from S * . As in ( [3]), authors supposed in their stability study of this model that the total population N does not vary in terms of S T , E T , I T , and R T . Thus, we take into account the same condition to obtain: with the homogeneous Neumann boundary conditions We prove that Y ε i are bounded in L 2 (Q) uniformly with respect to ε. For this, we put then, we can rewrite system (5.2) as is the semi-group generated by A, then the solution of (5.5) can be expressed as On the other hand the coefficients of the matrix H ε are bounded uniformly with respect to ε, using Gronwall's inequality, we have where β > 0 (i = 1, 2, 3, 4). Then Hence . Hence, then system (5. 2-5.4) can be written in the form and its solution can be expressed as By (5.6) and (5.10) one deduces that Thus all the coefficients of the matrix H ε tend to the corresponding coefficients of the matrix H in L 2 (Q). An application of Gronwall's inequality yields to Y ε i → Y i in L 2 (Q) as ε → 0, for i = 1, 2, 3, 4. The proof is complete.
Proof. Like in theorem (3.1), by making the change of variable s = T − t and the change of x) , (t, x) ∈ Q, i = 1, 2, 3, 4. We can easily prove the existence of the solution to this lemma .
To obtain the necessary conditions for the optimal control problem, applying standard optimality techniques, analyzing the objective functional and using relationships between the state and adjoint equations, we obtain a characterization of the optimal control. Theorem 5.2. Let ζ * be an optimal control of (2.5-2.7) and let y * ∈ W 1,2 ([0, T ] ; H (Ω )) with y * i ∈ G (T, Ω ) for i = 1, 2, 3, 4 be the optimal state, that is y * is the solution to (2.5-2.7) with the control ζ * . Then, and ζ * 2 = min ζ max 2 , max 0, Proof. We suppose ζ * is an optimal control and y * = y * 1 , y * 2 , y * 3 , y * 4 = (y 1 , y 2 , y 3 , y 4 ) (ζ * ) are the corresponding state variables. Consider ζ ε = ζ * + εh ∈ U ad and corresponding state solution y ε = y ε 1 , y ε 2 , y ε 3 , y ε 4 = (y 1 , y 2 , y 3 , y 4 ) (ζ ε ), we have (5.14)  . We use (5.1) and (5.13), we have  and U ad is convex, as the minimum of the objective functional is attained at ζ * it is seen that We take h = v − ζ * and we use (5.14-5.15) then (Ω)) 2 dt ≥ 0 for all v ∈ U ad By standard arguments varying v, we obtain and ζ * 2 = min ζ max 2 , max 0, 6. NUMERICAL RESULTS 6.1. Numerical method. In this section, the numerical results that reinforce and demonstrate the effectiveness of our strategy of control are provided. This strategy involves applying two types of treatment to exposed and to infected individuals respectively, to control the spread of the disease of tuberculosis. Concerning the numerical method, we give numerical simulations to our optimality system which is formulated by state equations with initial and boundary conditions (2.5-2.7), adjoint equations with transversality conditions (5.12), and optimal control characterization (5.13). We apply the forward-backward sweep method (FBSM) ( [21]) to solve • Test 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, go to Step 1.

Numerical results.
The authors in ( [1]) estimated that the initial susceptible, exposed, where the disease begins from the middle of Ω ) where we put 5 exposed to (TB) people infected people, 2 infected people, 0 removed people and 38 susceptibles for our numerical simulation.
To explain and demonstrate the influence of every control and its effect on spreading the disease, we have chosen to follow three different scenarios. We have simultaneously optimized the chemoprophylaxis (ζ 1 ) and treatment (ζ 2 ) controls in the first scenario. But in the second scenario, the treatment control (ζ 2 ) is maintained constant rather than optimized, whereas the chemoprophylaxis control (ζ 1 ) is optimized, and then in the third scenario, optimal treatment and constant chemoprophylaxis given to the infected persons are used.
6.2.1. Optimal chemoprophylaxis and treatment. In figures 5, 6, 7, and 8, when we apply our strategy of spatiotemporal control with two treatments. We assume that optimal treatments start on day t = 1 being the same day when the infection is identified in Ω. The effect of spatiotemporal controls treatment is very significant in reducing the propagation of infection.
Indeed, in figure 6, the infected population density after 4 years decreases from 35 infected in cases where there is no control, to 2 infected in cases where there is optimal control. In figure   6, we can see that there is a maximum number of individuals who are eliminated at about 43 people against less than 4 in case of no control, which is extremely beneficial and indicates the value of our strategy of control. 6.2.2. Constant treatment. In figures 5 − 11, the numerical results are studied with a single chemoprophylaxis of the population exposed and a constant treatment from the population infected. We find that in this scenario, the number of exposed persons decreases but less significantly than the first scenario, although the number of infected persons stays high and the     optimized in order to see its impact on the evolution of the system. According to the figures 12, 13 and 14, we notice that the number of exposed individuals remains high and they spend a To wrap this up, if the mechanisms on chemoprophylaxis and treatment are both are optimized, this will allow us to cure the infected population and at the same time to block the new infected cases within the exposed population. If we optimize only the treatment, we are going to cure the infected population at the present time, but the sickness won't be cured. Unlike the treatment, the mechanisms of chemoprophylaxis allows us to stop and block the infection in the future, but doesn't help in to cure infected population in the present.

Chemoprophylaxis and treatment.
In order to shed light on the optimal drug dosages, cost of treatment and chemotherapy, used to eradicate the transmission of tuberculosis in the population, we numerically present the optimal controls in figures (15)(16). This indicates that in the objective functional, the goal is not only to minimize the number of infected and exposed but also to reduce the costs of drugs, We can clearly see from figures (15-16) that the distribution  of drugs, treatment and chemotherapy, increases until reaching 1 after 3 years then it will cancel itself after 5 years Which gives a precise and optimal medication regimen

CONCLUSION
In this paper, we considered an existing model of tuberculosis that describes the dynamics of this disease with ODE's. Our approach is to extend this model using PDE's in order to give a more realistic description of the tuberculosis (TB). The application of a distributed optimal control pair for a spatiotemporal tuberculosis model with chemoprophylaxis and treatment are incorporated to reduce the actively infected individuals. We showed the existence of an optimal control pair and existence of solutions to the state. We used techniques of optimal control theory to characterize the controls, and derived the optimality system. The optimality system, which is composed by the system state, the dual system and the characteristic of the control, is solved numerically based on the forward-backward sweep method (FBSM). Numerical simulations of the resulting optimality system showed that the optimal strategy is more effective if we apply the combination of chemoprophylaxis and treatment together, so we can, very effectively, reduce the spread of infection as well as to block the new infected cases within the exposed individuals.
where A is a linear operator defined on a Banach space X, with the domain D(A) and g : [0, T ] × X → X is a given function. If X is a Hilbert space endowed with the scalar product (·, ·), then the linear operator A is called dissipative if (Az, z) ≤ 0, (∀z ∈ D(A)). Proof. or the proof of the theorem, we shall recall the following well-known results [33, Proposition 1.2 of chapter IV] (see also [34,35]). First recall a general existence result which we use in the sequel (Proposition 1.2, p. 175, [2]; see also [23,31]).

ACKNOWLEDGMENTS
Research reported in this paper was supported by the Moroccan Systems Theory Network.