# A Multiregion Discrete-Time Epidemic Model ofMycobacterium tuberculosis Infections: An Optimal Control Approach

Authors: Rachid Bouajaji; Hassan Laarabi; Mostafa Rachik; Abdelhadi Abta
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

The main goal of this article is to devise the spatial-temporal spread of TB, in multiple neighboring domains, taking into account the epidemiological diversity of their populations. However, since both the environment and any population are spatially heterogeneous, it is obviously desirable to include spatial structure into an epidemic model. Individuals with tuberculosis can spread the disease by moving from one area to another. In addition, people travel by air between cities, so diseases can be spread quickly between very distant places (as was the case with the COVID-19). In our model, each region’s studied population is divided into five compartments S, L1, I, L2, and R. Further, we introduce in our discrete systems three control variables which represent the effectiveness rates of vaccination, travel-blocking operation, and treatment. We focus in our study to control the outbreaks of an epidemic that affects a hypothetical population belonging to a specific region. Firstly, we analyze the epidemic model when the control strategy is based on the vaccination control only, and secondly, when the travel-blocking control is added, we finish with the introduction of the treatment control. The optimal control theory, based on Pontryagin’s maximum principle, is applied thrice in this paper, for the characterizations of the vaccination, travel-blocking, and treatment controls. The numerical results associated with the multipoint boundary value problems are obtained based on the forward-backward sweep method combined with progressive-regressive Runge–Kutta fourth-order schemes.


---

## Body

## 1. Introduction

Conventional mathematical models have focused only on the temporal spread of infection but ignored or neglected spatial dynamics, although the spatial spread of the epidemic was actually observed and many infectious diseases were transmitted from one region to another. Since an infectious case is first found at one location and then the disease spreads to other areas, the spatial variation of the population plays an important role in the control and eradication of this disease. However, due to the large mobility of people within a country or even worldwide, models that ignore the effect of space are not sufficient to give a realistic picture of the spread of the diseases. The epidemic could spread over large and remote areas. It can reach continents, and among the examples in human history there is the black death (plague) which spread in Europe in the 13th century, then measles and smallpox between the 16th and 17th centuries which spread in the New World [1, 2], HIV/AIDS in the year 1981, the West Nile virus that emerged in North America in the late 1990s [3], SARS that appeared in Asia in 2003 [4], and recently the new corona virus (COVID-19) that appeared at the end of 2019 [5]. It continues to spread still in the world and jumps from one continent to another, especially due to the explosive development of means of air, sea, and land transport, which contributed to a wide spread. The overlap of many different components makes modeling the spatial spread of infectious diseases a complex task. And spatial heterogeneity requires formulating a small model for each region (city, country, and so on). One possible approach is to consider people traveling between distinct geographic regions (subfields of the field of study), assuming that contagion does not occur while traveling. With regard to diseases which affect animals and insects (vectors) and which can also be transmitted to humans [6], it is very important to define a precise and clear pattern of not only the contact/mixing between humans and different species but also the movement patterns of each species involved, including humans. The importance of modeling the spatial spread of infectious diseases can be highlighted by the outbreaks of corona disease in many countries around the world. Before this disease was diagnosed at the end of 2019, infected people were moving from city to city, and after diagnosing some infected cases and identifying the problem, it took some time for the political decision to prevent all movement, especially from the affected areas, and the regions were actually divided according to the number of infected cases it contains (from one affected continent to another, or from one affected country to another, or from one affected city to another …) [7]. Thus, in order to control spatial spread of the disease, it was necessary to follow all precautions and take all necessary measures to reduce displacement to and from the affected areas.The mathematical modeling approach is one of the useful tools and can be used to study the dynamical behavior of a disease in a better way and to develop some useful strategies for possible eradication and control. In the last few years, several transmission models based on different approaches have been suggested to analyze the dynamics of TB disease in various regions. Most of the researchers who have studied mathematical models of tuberculosis have focused on the continuous-time models described by ordinary differential equations (see [8–15]). We observe that most of those researchers focused on the continuous-time models described by the differential equations. It is noted that, in recent years, more and more attention has been given to discrete-time models (see [16–19] and the references cited therein). The reasons for adopting discrete modeling are as follows. Firstly, the statistical data are collected at discrete moments (day, week, month, or year). So, it is more direct and more accurate and timely to describe the disease using discrete-time models than continuous-time models. Secondly, the use of discrete-time models can avoid some mathematical complexities such as choosing a function space and regularity of the solution. Thirdly, the numerical simulations of continuous-time models are obtained by the way of discretization. The deterministic, spatial disease models have the advantage that they allow for the investigation, either analytically or numerically, of various spatiotemporal control strategies. The oldest control method for the spatial spread of infectious diseases is quarantine through travel restrictions.A study of the impact on disease dynamics of this type of quarantine was undertaken by Sattenspiel and Herring [20], as applied to the model of Sattenspiel and Dietz [21]. Numerically, this intervention method was shown to have little impact on the ultimate course of the epidemic.There is increasing interest in studying the spatial spread of infectious diseases [22–25]. Most of the models studied have been partly continuous because of their mathematical tractability. In Allen et al. [26], a discrete-time, age-independent SIR-type epidemic model with n subgroups is formulated and analyzed and then applied to a measles epidemic model on a university campus. In similar pattern, Zakary et al. [7] devised a discrete-time SIR model that describes the propagation of a disease in a population of individuals who travel between p regions (domains). In the same way, in our article, we will study the propagation of tuberculosis in discrete-time case between several regions. Tuberculosis is a contagious disease, secondary to infection with “bacillus of Koch” (Mycobacterium tuberculosis). This bacterial agent is transmitted by air via the droplets contaminated by the bacterium, which are suspended in the exhaled air by patients, especially during coughing. Inhaling a small number of contaminated droplets is enough to infect an individual. The displacement of populations (travelers and refugees) has largely contributed recently in the spread of the disease in the world. People who are more likely to acquire TB infection are the following:(1)
People recently exposed to someone who has symptomatic TB disease.(2)
People who live in congregate settings with high-risk persons.(3)
People who live or have lived in countries where TB is common.(4)
People who are healthcare workers who are in contact with TB patients when proper infection control procedures are not followed.Many people who acquire TB infection do not have symptoms and may never develop TB disease. These people have latent TB infections (LTBIs) [27]. After exposure to the tuberculosis bacillus, some people developed a primary infection, which is controlled by the immune system in 90% of cases: tuberculosis is then said to be “latent.” The bacillus remains in the body, but the immune system prevents it from multiplying. In 10% of those infected, the bacillus is not sufficiently controlled by the immune system, and these people develop a form of so-called “active” tuberculosis, which will cause disease and complications.The Moroccan Ministry of Health registers nearly 30,000 cases of tuberculosis each year, including new cases and cases of relapse. The incidence rate is around 87 cases per 100,000 inhabitants, of which pulmonary tuberculosis represents half. They also noted that this disease affects young people aged 15 to 45, adding that among the causes of the spread of this disease are socio-economic factors, especially housing conditions, poverty, and malnutrition [28, 29].In our model, the population total is divided into five categories:(i)
(S): susceptible, who have never encountered the Mycobacterium.(ii)
(L1): early latent, that is, individuals recently infected (less than two years) but not infectious.(iii)
(I): infected, that is, individuals who have active tuberculosis and are infectious.(iv)
(L2): persistent latent, that is, individuals who were infected and remain latent.(v)
(R): recovered, that is, individuals who were previously infected and treated.We are interested in studying the dynamics of Koch bacillus spread at several region in a discrete-time model. In addition, in order to find the best strategy to reduce the number of susceptible, infected who have active MTB or recently and persistent infected latent, and to rise the number of recovered, we will use three control, namely, vaccination, travel-blocking, and treatment programs. The paper is organized as follows: in Section2, the model. S, L1, I, L2, and R is described for a multiregion discrete model. In Sections 3 and 4, we present the optimal control problem for the proposed model where we give some results concerning the existence of the optimal controls and we characterize these optimal controls using Pontryagin’s maximum principle in discrete time; illustrative numerical simulations are given at the end. Finally, we conclude the paper in Section 5.


## 2. Formulation of the Multiregion Discrete Mathematical Model

We assume that there arep geographical regions (domains) of the domain studied Ω (in Figure 1, Ω is the region of oriental with p = 8).Figure 1
Region of Oriental in Morocco. This region is divided into eight subregions (or domains)Ωj1⩽j⩽8: Oujda Angad Ω1, Berkane Ω2, Nador Ω3, Driouch Ω4, Guercif Ω5, Taourirt Ω6, Jerada Ω7, and Figuig Ω8.LetΩ=∪j=1pΩj, and let NiΩj be the population of domain Ωj at time i, i.e., the number of individuals who are physically present in Ωj, both residents and travelers.According to the disease transmission mechanism, the host population ofΩj is grouped into five epidemiological compartments, and let SiΩj,L1,iΩj,IiΩj,L2,iΩj, and RiΩj be the number of individuals in the susceptible, early latent, infective, persistent latent, and removed compartments of Ωj at time i, respectively. In addition to the death and recruitment, there are population movements among those five epidemiological compartments from time unit i to time i + 1. We assume that the recruited individuals (by birth and immigration) are constant and enter the susceptible compartment. To make the model a little more realistic, but in order to work with a constant overall population, we suppose that birth and death occur with the same rate. In addition, we suppose that individuals who are out of their domain do not give birth, and so birth occurs in the home domain at a per capita rate dj>0. And death takes place anywhere with a per capita rate dj. After one unit time, the susceptible individuals may stay in the susceptible compartment, or get infected and move to the infectious compartment, or die. Disease transmission is assumed to occur between individuals present in a given domain Ωj. The individuals in the infectious compartment may continue to be infected, or recover and be transferred to the recovered compartment, or die. Moreover, here we have assumed that there is no mortality due to infection. For instance, since the average life span of Moroccan people is 76.45 years, we have d=1/76.45yr−1 [30]. Individuals in the early latent compartment L1 can progress either to active disease I with rate ϕδ or to a persistent latent infection L2 with rate 1−Φδ, following the approach in [31]. Parameter ϕ reflects that only 5% of infected individuals will ever develop active TB [15]. We choose d such that the progression rate from early infections to active disease is ϕδ=0.6yr−1, which roughly approximates the data by Styblo [32], describing the proportions of disease development after conversion. For the rates of reactivation we adopt ω=0.0003yr−1 for untreated latent infections [33] and for those who have undergone a therapeutic intervention. As in Gomes et al. [34], the partial susceptibility factor affecting the rate of exogenous reinfection of untreated individuals, s, is fixed at 0.25, in accordance with the highest estimates of protection conferred by BCG vaccination (see [35]). In treated patients, this factor becomes σR, for which several exploratory values are adopted. Treatment of different infection stages is implemented at specific rates: τ0 applies to active TB and represents the rate of recovery (typically as a result of treatment, though here it also accounts for the infrequent natural recovery); τ1 and τ2 apply, respectively, to the latent classes L1 and L2 as the rates at which chemotherapy or a postexposure vaccine is applied. The rate τ0 is fixed at 2yr−1, corresponding to an average duration of infectiousness of 6 months, while τ1 and τ2 are considered at different exploratory values (see [15]). The values of the model parameters presented in control system (2) are given in Table 1. The disease transmission in a given domain Ωj at time i is modeled using standard incidence, given by(1)∑k=1pβjkIiΩkNiΩjSiΩj,where the disease transmission coefficient βjk>0 is the proportion of adequate contacts in domain Ωj between a susceptible from Ωj and an infective from another domain Ωk. The multiregion discrete-time TB model associated with Ωj is written as follows:(2)Si+1Ωj=SiΩj+djNiΩj−∑k=1pβjkIiΩkNiΩjSiΩj−djSiΩj,L1,i+1Ωj=L1,iΩj+∑k=1pβjkIiΩkNiΩjSiΩj+∑k=1pβjkIiΩkNiΩjσjL2,iΩj+σRjRiΩj−δj+τ1j+djL1,iΩj,Ii+1Ωj=IiΩj+ϕjδjL1,iΩj+ωjL2,iΩj+ωRjRiΩj−τ0j+djIiΩj,L2,i+1Ωj=L2,iΩj+1−ϕjδjL1,iΩj−∑k=1pβjkIiΩkNiΩjσjL2,iΩj−ωj+τ2j+djL2,iΩj,Ri+1Ωj=RiΩj+τ0jIiΩj+τ1jL1,iΩj+τ2jL2,iΩj−∑k=1pβjkIiΩkNiΩjσRjRiΩj−ωRj+djRiΩj,where NiΩj=SiΩj+L1,iΩj+IiΩj+L2,iΩj+RiΩj is the population size corresponding to the domain Ωj at time i. It is clear that the population size remains constant for all i≥0; in fact(3)Ni+1Ωj=Si+1Ωj+L1,i+1Ωj+Ii+1Ωj+L2,i+1Ωj+Ri+1Ωj.=SiΩj+djNiΩj−djSiΩj+L1,iΩj−djL1,iΩj+IiΩj−djIiΩj+L2,iΩj−djL2,iΩj+RiΩj−djRiΩj.=SiΩj+L1,iΩj+IiΩj+L2,iΩj+RiΩj+djNiΩj−djSiΩj+L1,iΩj+IiΩj+L2,iΩj+RiΩj.=NiΩj.Table 1
Parameter values.
SymbolDescriptionValueSourceDNatural death rate1/76.45[30]δRate at which individuals leaveL112yr−1[12]ϕProportion of individuals going to compartment I0.05[33]ωRate of endogenous reactivation for persistent latent infections0.0003yr−1[36]ωRRate of endogenous reactivation for treated individuals0.00002yr−1[12]σFactor reducing the risk of infection as a result of acquired immunity to a previous infection for persistent latent individuals0.25[12]σRFactor reducing the risk of infection as a result of acquired immunity to a previous infection for treated individuals0.25[12]τ0Rate of recovery under treatment of active TB2yr−1[12]τ1Rate of recovery of early latent individuals under postexposure interventions2yr−1[12]τ2Rate of recovery of persistent latent individuals under postexposure interventions1yr−1[12]NTotal population10000; 20000; 30000EstimatedTTotal simulation duration20yr


## 3. The Model with Vaccination and Travel Blocking

### 3.1. Presentation of the Model

LetP=1,2,…,p and let us denote by PH⊂P the set of indexes of domain at high risk; we introduce a control variable νjΩk that characterizes the travel-blocking operation, in order to restrict movements from domains Ωk,k∈PH to Ωj, where(4)νijΩk<1,∀k∈PHk≠j,νijΩk=1,elsewhere,∀i≥0.We also introduce a control variableuiΩj that characterizes the effectiveness of vaccination in the aforementioned model (2). Then, for a given region Ωj the model is given as follows:(5)Si+1Ωj=SiΩj+djNiΩj−∑k=1pνijΩkβjkIiΩkNiΩjSiΩj−djSiΩj−uiΩjSiΩj,L1,i+1Ωj=L1,iΩj+∑k=1pνijΩkβjkIiΩkSiΩjNiΩj+∑k=1pβjkIiΩkNiΩjσjL2,iΩj+σRjRiΩj−δj+τ1j+djL1,iΩj,Ii+1Ωj=IiΩj+ϕjδjL1,iΩj+ωjL2,iΩj+ωRjRiΩj−τ0j+djIiΩj,L2,i+1Ωj=L2,iΩj+1−ϕjδjL1,iΩj−∑k=1pβjkIiΩkNiΩjσjL2,iΩj−ωj+τ2j+djL2,iΩj,Ri+1Ωj=RiΩj+τ0jIiΩj+τ1jL1,iΩj+τ2jL2,iΩj−∑k=1pβjkIiΩkNiΩjσRjRiΩj−ωRj+djRiΩj+uiΩjSiΩj.


#### 3.1.1. An Optimal Control Approach

We are interested in controlling the population of regionΩj. The problem is to minimize the objective functional given by(6)JjνjΩk,uΩj=AjIITΩj−AjRRTΩj+∑k∈PH∑i=1T−1AjIIiΩj−AjRRiΩj+Bk2νijΩk2+Cj2uiΩj2,where Bk > 0, Cj > 0, AjI > 0, and AjR > 0 are the weight constants of controls, the infected, and the removed in region Ωj, respectively. νjΩk=νjΩkk∈PH where νjΩk=ν0jΩk,…,νT−1jΩk and uΩj=u0Ωj,…,uT−1Ωj. Our goal is to minimize the infected and the cost of applying controls and increase the number of recovered in Ωj.In other words, we are seeking optimal controlνjΩk∗ and uΩj∗ such that(7)JjνjΩk∗,uΩj∗=minνjΩk∈Vj,uΩj∈UjJjνjΩk,uΩj,where Vj and Uj are the control sets defined by Vj=νjΩlmesurable0<νminjΩl≤νijΩl≤νmaxjΩl<1,i=0,1,…,T−1,l∈PH. Uj=uΩjmesurable0<uminΩj≤uiΩj≤umaxΩj<1,i=0,1,…,T−1.By using Pontryagin’s maximum principle [37], we derive necessary conditions for our optimal controls. For this purpose, we defined the Hamiltonian as(8)HiΩj=∑k∈PHAjIIiΩj−AjRRiΩj+Bk2νijΩk2+Ck2uiΩj2+λ1,i+1jSiΩj+djNiΩj−∑k=1pνijΩkβjkIiΩkNiΩjSiΩj−djSiΩj−uiΩjSiΩj+λ2,i+1jL1,iΩj+∑k=1pνijΩkβjkIiΩkSiΩjNiΩj+∑k=1pβjkIiΩkNiΩjσjL2,iΩj+σRjRiΩj−δj+τ1j+djL1,iΩj+λ3,i+1jIiΩj+ϕjδjL1,iΩj+ωjL2,iΩj+ωRjRiΩj−τ0j+djIiΩj+λ4,i+1jL2,iΩj+1−ϕjδjL1,iΩj−∑k=1pβjkIiΩkNiΩjσjL2,iΩj−ωj+τ2j+djL2,iΩj+λ5,i+1jRiΩj+τ0jIiΩj+τ1jL1,iΩj+τ2jL2,iΩj−∑k=1pβjkIiΩkNiΩjσRjRiΩj−ω1j+djRiΩj+uiΩjSiΩj.Theorem 1.
Given an optimal controlνjΩk∗,k∈PH and solutions SΩj∗,L1Ωj∗,IΩj∗,L2Ωj∗ and RΩj∗, there exist λs,ij,i=1,2,…,T,j=1,2,…,p,s=1,2,3,4,5. The adjoint variables satisfy the following equations:(9)λ1,ij=1−djλ1,i+1j+∑k=1pνijΩkβjkIiΩkNiΩjλ2,i+1j−λ1,i+1j+uiΩjλ5,i+1j−λ1,i+1j,λ2,ij=1−djλ2,i+1j+δjλ4,i+1j−λ2,i+1j+ϕjδjλ3,i+1j−λ4,i+1j+τ1jλ5,i+1j−λ2,i+1j,λ3,ij=AjI+1−djλ3,i+1j+βjjSiΩjNiΩjλ2,i+1j−λ1,i+1j+σjβjjL2,iΩjNiΩjλ2,i+1j−λ4,i+1j+σRjβjjRiΩjNiΩjλ2,i+1j−λ5,i+1j+τ0jλ5,i+1j−λ3,i+1j,λ4,ij=1−djλ4,i+1j+∑k=1pσjβjkIiΩkNiΩjλ2,i+1j−λ4,i+1j+ωjλ3,i+1j−λ4,i+1j+τ2jλ5,i+1j−λ4,i+1j,λ5,ij=−AjR+1−djλ5,i+1j+∑k=1pσRjβjkIiΩkNiΩjλ2,i+1j−λ5,i+1j+ωRjλ3,i+1j−λ5,i+1j,where λ1,Tj=0;λ2,Tj=0;λ3,Tj=AjI;λ4,Tj=0; and λ5,Tj=−AjR are the transversality conditions.
In addition, we have(10)νijΩl∗=minmaxνminjΩl,βjlBlIiΩl∗SiΩl∗NiΩlλ1,i+1l−λ2,i+1l,νmaxjΩl,i=0,1,…,T−1,l∈PH,(11)uiΩj∗=minmaxuminΩj,SiΩj∗λ1,i+1j−λ5,i+1jCj,umaxΩj,i=0,1,…,T−1.Proof.
Using Pontryagin’s maximum principle [37, 38], we obtain the following adjoint equations:(12)λ1,ij=∂Hi∂SiΩj=1−∑k=1pνijΩkβjkIiΩkNiΩj−dj−uiΩjλ1,i+1j+∑k=1pνijΩkβjkIiΩkNiΩjλ2,i+1j+uiΩjλ5,i+1j,λ2,ij=∂Hi∂L1,iΩj=1−δj−τ1j−djλ2,i+1j+ϕjδjλ3,i+1j+1−ϕjδjλ4,i+1j+τ1jλ5,i+1j,λ3,ij=∂Hi∂IiΩj=AjI−λ1,i+1jβjjSiΩjNiΩj+λ2,i+1jβjj1NiΩj+βjjσjL2,iΩj+σRjRiΩjNiΩj+λ3,i+1j1−τ0j−dj−λ4,i+1jβjjσjL2,iΩjNiΩj+λ5,i+1jτ0j−βjkσRjRiΩjNiΩj,λ4,ij=∂Hi∂L2,iΩj=λ2,i+1j∑k=1pβjkσjIiΩkNiΩj+λ3,i+1jωj+λ4,i+1j1−τ2j−ωj−dj−∑k=1pβjkσjIiΩkNiΩj+λ5,i+1jτ2j,λ5,ij=∂Hi∂RiΩj=−AjR+λ2,i+1j∑k=1pσRjβjkIiΩkNiΩj+λ3,i+1jωRj+λ5,i+1j1−ωRj−dj−∑k=1pβjkσRjIiΩkNiΩj,with λ1,Tj=0,λ2,Tj=0,λ3,Tj=AjI,λ4,Tj=0,λ5,Tj=−AjR,j=1,2,…,p.
To obtain the optimality conditions, we take the variation with respect to controlsνijΩl and uiΩj and set it equal to zero:(13)∂Hi∂νijΩl=BlνijΩl−λ1,i+1jβjlIiΩlNiΩjSiΩj+λ2,i+1jβjlIiΩlNiΩjSiΩj=0,∂Hi∂uiΩj=CjuiΩj−λ1,i+1jSiΩj+λ5,i+1jSiΩj=0,and then we obtain the optimal control pair(14)νijΩl=βjlBlIΩlSiΩjNiΩjλ1,i+1j−λ2,i+1j,uiΩj=SiΩjλ1,i+1j−λ5,i+1jCj.
By the bounds inVj and Uj of the controls, it is easy to obtain νijΩl∗ and uiΩj∗ in the following form:(15)νijΩl∗=minmaxνminjΩl,βjlBlIiΩlSiΩlNiΩlλ1,i+1l−λ2,i+1l,νmaxjΩl,i=0,1,…,T−1,l∈PH,uiΩj∗=minmaxuminΩj,SiΩjλ1,i+1j−λ5,i+1jCj,umaxΩj,i=0,1,…,T−1.


#### 3.1.2. Numerical Results

We now present numerical simulations associated with the aforementioned optimal control problem. We write a code inMATLABTM and simulate our results using different data. The optimality systems are solved based on an iterative discrete scheme that converges following an appropriate test similar to the one related to the forward-backward sweep method (FBSM). The state system with an initial guess is solved forward in time and then the adjoint system is solved backward in time because of the transversality conditions. Afterwards, we update the optimal control values using the values of state and costate variables obtained at the previous steps. Finally, we execute the previous steps till a tolerance criterion is reached.The multiregion TB epidemic model we suggest here is applicable for any number of regionsp. For example, in order to show the importance of our work, we choose p=3, i.e., we consider three regions (domains) Ω1,Ω2, and Ω3 with different parameters based on [12] and given in Table 1. Other epidemiological and numerical parameters are presented in Tables 1 and 2, respectively. The initial value of each category, SΩi0=76/120NΩi;L1Ωi0=36/120NΩi;IΩi0=4/120NΩi;L2Ωi0=2/120NΩi; and RΩi0=1/120NΩi, is given in Table 2, which is based on [39] (Figures 2–13).Table 2
Parameter values ofβ utilized for the resolution of multiregion discrete systems (2), (5), and (16), leading to simulations obtained from Figures 2–13, with the initial conditions S0,L10,I0,L20, and R0 associated with the three regions Ω1,Ω2, and Ω3.
NS (0)L10I (0)L20R (0)βΩ11000063343083334166830.05Ω2200001266661686663341660.21Ω33000019000925010005002500.075Figure 2
The state variables in three regions without control.
(a)(b)(c)Figure 3
Comparison for each state variable for the three regions without control.
(a)(b)(c)(d)(e)Figure 4
The state variables in three regions with BCG vaccination and travel-blocking control.
(a)(b)(c)Figure 5
Comparison for each state variable for the three regions with BCG vaccination and travel-blocking control.
(a)(b)(c)(d)(e)Figure 6
The BCG vaccination and travel-blocking control in three regions.
(a)(b)Figure 7
Comparison for each state variable for the three regions with BCG vaccination, travel blocking, and treatment control.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 8
The state variables in three regions with controls vaccination, travel blocking, and treatment.
(a)(b)(c)Figure 9
Comparison of the susceptible in each region without control and with control in two cases.
(a)(b)(c)Figure 10
Comparison of the early latent in each region without control and with control in two cases.
(a)(b)(c)Figure 11
Comparison of the infected in each region without control and with control in two cases.
(a)(b)(c)Figure 12
Comparison of the persistent latent in each region without control and with control in two cases.
(a)(b)(c)Figure 13
Comparison of the recovered in each region without control and with control in two cases.
(a)(b)(c)Since the BCG vaccine in Morocco is mandatory, we considered the BCG vaccination control in the three regions given by (11). We are also interested in controlling Ω1 and Ω3, by assuming that Ω2 is at high risk of infection (see Figure 2); then, PH=2. We introduce the second control ν1Ω2,ν3Ω2 given by (10), in order to reduce entry of infected people from region Ω2.In Figure2, we can observe that in the absence of control and in the presence of an epidemic which spreads in three regions characterized by different parameters, the number of susceptible individuals varies from SΩ10=6334,SΩ20=12666, and SΩ30=19000 as initial conditions to 5088, 12094, and 19230, respectively, in the three regions; on the other hand, when the BCG vaccination controls uΩ1uΩ2anduΩ3, given by (11), are introduced in the three regions and we blocked the movement of infected individuals to the Ω2 region by introducing the controls v1Ω2,v3Ω2 given by (10), we observe that the number of susceptible subjects decreases from the initial state in the three regions to 703, 545, and 3556, respectively (see Figure 4).To understand well, we compared each state variable in the three regions without control (see Figure3), and with the introduction of vaccination and blocking the movement (see Figure 5), we can see clearly in Figure 5 that the curves of the susceptible decrease in all regions especially in second region after 2 years. On the other hand, the recovered curves increase markedly after the implementation of the first strategy which proves the impact of these controls.In Figure4, we can clearly see that after two years, there is effectiveness of vaccination on all the regions and especially on the second region controlled by the blockage of displacement, which is seen through the decrease in the curves related to the number of susceptible. At the same time, we see that we have a rapid increase in the number of recovered, especially in the second region.In Figure6(a), it can be seen that the higher the value of the vaccination control (uΩ2), the more the number of susceptible is reduced, and this is clearly seen in Figure 5. On the other hand, in Figure 6(b), we find that when the travel-blocking control takes small values, i.e., when the contact between susceptible and infected individuals is too low, we obtain good results. This can be seen by comparing regions one and three for the number of infected (see Figure 5).


## 3.1. Presentation of the Model

LetP=1,2,…,p and let us denote by PH⊂P the set of indexes of domain at high risk; we introduce a control variable νjΩk that characterizes the travel-blocking operation, in order to restrict movements from domains Ωk,k∈PH to Ωj, where(4)νijΩk<1,∀k∈PHk≠j,νijΩk=1,elsewhere,∀i≥0.We also introduce a control variableuiΩj that characterizes the effectiveness of vaccination in the aforementioned model (2). Then, for a given region Ωj the model is given as follows:(5)Si+1Ωj=SiΩj+djNiΩj−∑k=1pνijΩkβjkIiΩkNiΩjSiΩj−djSiΩj−uiΩjSiΩj,L1,i+1Ωj=L1,iΩj+∑k=1pνijΩkβjkIiΩkSiΩjNiΩj+∑k=1pβjkIiΩkNiΩjσjL2,iΩj+σRjRiΩj−δj+τ1j+djL1,iΩj,Ii+1Ωj=IiΩj+ϕjδjL1,iΩj+ωjL2,iΩj+ωRjRiΩj−τ0j+djIiΩj,L2,i+1Ωj=L2,iΩj+1−ϕjδjL1,iΩj−∑k=1pβjkIiΩkNiΩjσjL2,iΩj−ωj+τ2j+djL2,iΩj,Ri+1Ωj=RiΩj+τ0jIiΩj+τ1jL1,iΩj+τ2jL2,iΩj−∑k=1pβjkIiΩkNiΩjσRjRiΩj−ωRj+djRiΩj+uiΩjSiΩj.


### 3.1.1. An Optimal Control Approach

We are interested in controlling the population of regionΩj. The problem is to minimize the objective functional given by(6)JjνjΩk,uΩj=AjIITΩj−AjRRTΩj+∑k∈PH∑i=1T−1AjIIiΩj−AjRRiΩj+Bk2νijΩk2+Cj2uiΩj2,where Bk > 0, Cj > 0, AjI > 0, and AjR > 0 are the weight constants of controls, the infected, and the removed in region Ωj, respectively. νjΩk=νjΩkk∈PH where νjΩk=ν0jΩk,…,νT−1jΩk and uΩj=u0Ωj,…,uT−1Ωj. Our goal is to minimize the infected and the cost of applying controls and increase the number of recovered in Ωj.In other words, we are seeking optimal controlνjΩk∗ and uΩj∗ such that(7)JjνjΩk∗,uΩj∗=minνjΩk∈Vj,uΩj∈UjJjνjΩk,uΩj,where Vj and Uj are the control sets defined by Vj=νjΩlmesurable0<νminjΩl≤νijΩl≤νmaxjΩl<1,i=0,1,…,T−1,l∈PH. Uj=uΩjmesurable0<uminΩj≤uiΩj≤umaxΩj<1,i=0,1,…,T−1.By using Pontryagin’s maximum principle [37], we derive necessary conditions for our optimal controls. For this purpose, we defined the Hamiltonian as(8)HiΩj=∑k∈PHAjIIiΩj−AjRRiΩj+Bk2νijΩk2+Ck2uiΩj2+λ1,i+1jSiΩj+djNiΩj−∑k=1pνijΩkβjkIiΩkNiΩjSiΩj−djSiΩj−uiΩjSiΩj+λ2,i+1jL1,iΩj+∑k=1pνijΩkβjkIiΩkSiΩjNiΩj+∑k=1pβjkIiΩkNiΩjσjL2,iΩj+σRjRiΩj−δj+τ1j+djL1,iΩj+λ3,i+1jIiΩj+ϕjδjL1,iΩj+ωjL2,iΩj+ωRjRiΩj−τ0j+djIiΩj+λ4,i+1jL2,iΩj+1−ϕjδjL1,iΩj−∑k=1pβjkIiΩkNiΩjσjL2,iΩj−ωj+τ2j+djL2,iΩj+λ5,i+1jRiΩj+τ0jIiΩj+τ1jL1,iΩj+τ2jL2,iΩj−∑k=1pβjkIiΩkNiΩjσRjRiΩj−ω1j+djRiΩj+uiΩjSiΩj.Theorem 1.
Given an optimal controlνjΩk∗,k∈PH and solutions SΩj∗,L1Ωj∗,IΩj∗,L2Ωj∗ and RΩj∗, there exist λs,ij,i=1,2,…,T,j=1,2,…,p,s=1,2,3,4,5. The adjoint variables satisfy the following equations:(9)λ1,ij=1−djλ1,i+1j+∑k=1pνijΩkβjkIiΩkNiΩjλ2,i+1j−λ1,i+1j+uiΩjλ5,i+1j−λ1,i+1j,λ2,ij=1−djλ2,i+1j+δjλ4,i+1j−λ2,i+1j+ϕjδjλ3,i+1j−λ4,i+1j+τ1jλ5,i+1j−λ2,i+1j,λ3,ij=AjI+1−djλ3,i+1j+βjjSiΩjNiΩjλ2,i+1j−λ1,i+1j+σjβjjL2,iΩjNiΩjλ2,i+1j−λ4,i+1j+σRjβjjRiΩjNiΩjλ2,i+1j−λ5,i+1j+τ0jλ5,i+1j−λ3,i+1j,λ4,ij=1−djλ4,i+1j+∑k=1pσjβjkIiΩkNiΩjλ2,i+1j−λ4,i+1j+ωjλ3,i+1j−λ4,i+1j+τ2jλ5,i+1j−λ4,i+1j,λ5,ij=−AjR+1−djλ5,i+1j+∑k=1pσRjβjkIiΩkNiΩjλ2,i+1j−λ5,i+1j+ωRjλ3,i+1j−λ5,i+1j,where λ1,Tj=0;λ2,Tj=0;λ3,Tj=AjI;λ4,Tj=0; and λ5,Tj=−AjR are the transversality conditions.
In addition, we have(10)νijΩl∗=minmaxνminjΩl,βjlBlIiΩl∗SiΩl∗NiΩlλ1,i+1l−λ2,i+1l,νmaxjΩl,i=0,1,…,T−1,l∈PH,(11)uiΩj∗=minmaxuminΩj,SiΩj∗λ1,i+1j−λ5,i+1jCj,umaxΩj,i=0,1,…,T−1.Proof.
Using Pontryagin’s maximum principle [37, 38], we obtain the following adjoint equations:(12)λ1,ij=∂Hi∂SiΩj=1−∑k=1pνijΩkβjkIiΩkNiΩj−dj−uiΩjλ1,i+1j+∑k=1pνijΩkβjkIiΩkNiΩjλ2,i+1j+uiΩjλ5,i+1j,λ2,ij=∂Hi∂L1,iΩj=1−δj−τ1j−djλ2,i+1j+ϕjδjλ3,i+1j+1−ϕjδjλ4,i+1j+τ1jλ5,i+1j,λ3,ij=∂Hi∂IiΩj=AjI−λ1,i+1jβjjSiΩjNiΩj+λ2,i+1jβjj1NiΩj+βjjσjL2,iΩj+σRjRiΩjNiΩj+λ3,i+1j1−τ0j−dj−λ4,i+1jβjjσjL2,iΩjNiΩj+λ5,i+1jτ0j−βjkσRjRiΩjNiΩj,λ4,ij=∂Hi∂L2,iΩj=λ2,i+1j∑k=1pβjkσjIiΩkNiΩj+λ3,i+1jωj+λ4,i+1j1−τ2j−ωj−dj−∑k=1pβjkσjIiΩkNiΩj+λ5,i+1jτ2j,λ5,ij=∂Hi∂RiΩj=−AjR+λ2,i+1j∑k=1pσRjβjkIiΩkNiΩj+λ3,i+1jωRj+λ5,i+1j1−ωRj−dj−∑k=1pβjkσRjIiΩkNiΩj,with λ1,Tj=0,λ2,Tj=0,λ3,Tj=AjI,λ4,Tj=0,λ5,Tj=−AjR,j=1,2,…,p.
To obtain the optimality conditions, we take the variation with respect to controlsνijΩl and uiΩj and set it equal to zero:(13)∂Hi∂νijΩl=BlνijΩl−λ1,i+1jβjlIiΩlNiΩjSiΩj+λ2,i+1jβjlIiΩlNiΩjSiΩj=0,∂Hi∂uiΩj=CjuiΩj−λ1,i+1jSiΩj+λ5,i+1jSiΩj=0,and then we obtain the optimal control pair(14)νijΩl=βjlBlIΩlSiΩjNiΩjλ1,i+1j−λ2,i+1j,uiΩj=SiΩjλ1,i+1j−λ5,i+1jCj.
By the bounds inVj and Uj of the controls, it is easy to obtain νijΩl∗ and uiΩj∗ in the following form:(15)νijΩl∗=minmaxνminjΩl,βjlBlIiΩlSiΩlNiΩlλ1,i+1l−λ2,i+1l,νmaxjΩl,i=0,1,…,T−1,l∈PH,uiΩj∗=minmaxuminΩj,SiΩjλ1,i+1j−λ5,i+1jCj,umaxΩj,i=0,1,…,T−1.


### 3.1.2. Numerical Results

We now present numerical simulations associated with the aforementioned optimal control problem. We write a code inMATLABTM and simulate our results using different data. The optimality systems are solved based on an iterative discrete scheme that converges following an appropriate test similar to the one related to the forward-backward sweep method (FBSM). The state system with an initial guess is solved forward in time and then the adjoint system is solved backward in time because of the transversality conditions. Afterwards, we update the optimal control values using the values of state and costate variables obtained at the previous steps. Finally, we execute the previous steps till a tolerance criterion is reached.The multiregion TB epidemic model we suggest here is applicable for any number of regionsp. For example, in order to show the importance of our work, we choose p=3, i.e., we consider three regions (domains) Ω1,Ω2, and Ω3 with different parameters based on [12] and given in Table 1. Other epidemiological and numerical parameters are presented in Tables 1 and 2, respectively. The initial value of each category, SΩi0=76/120NΩi;L1Ωi0=36/120NΩi;IΩi0=4/120NΩi;L2Ωi0=2/120NΩi; and RΩi0=1/120NΩi, is given in Table 2, which is based on [39] (Figures 2–13).Table 2
Parameter values ofβ utilized for the resolution of multiregion discrete systems (2), (5), and (16), leading to simulations obtained from Figures 2–13, with the initial conditions S0,L10,I0,L20, and R0 associated with the three regions Ω1,Ω2, and Ω3.
NS (0)L10I (0)L20R (0)βΩ11000063343083334166830.05Ω2200001266661686663341660.21Ω33000019000925010005002500.075Figure 2
The state variables in three regions without control.
(a)(b)(c)Figure 3
Comparison for each state variable for the three regions without control.
(a)(b)(c)(d)(e)Figure 4
The state variables in three regions with BCG vaccination and travel-blocking control.
(a)(b)(c)Figure 5
Comparison for each state variable for the three regions with BCG vaccination and travel-blocking control.
(a)(b)(c)(d)(e)Figure 6
The BCG vaccination and travel-blocking control in three regions.
(a)(b)Figure 7
Comparison for each state variable for the three regions with BCG vaccination, travel blocking, and treatment control.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 8
The state variables in three regions with controls vaccination, travel blocking, and treatment.
(a)(b)(c)Figure 9
Comparison of the susceptible in each region without control and with control in two cases.
(a)(b)(c)Figure 10
Comparison of the early latent in each region without control and with control in two cases.
(a)(b)(c)Figure 11
Comparison of the infected in each region without control and with control in two cases.
(a)(b)(c)Figure 12
Comparison of the persistent latent in each region without control and with control in two cases.
(a)(b)(c)Figure 13
Comparison of the recovered in each region without control and with control in two cases.
(a)(b)(c)Since the BCG vaccine in Morocco is mandatory, we considered the BCG vaccination control in the three regions given by (11). We are also interested in controlling Ω1 and Ω3, by assuming that Ω2 is at high risk of infection (see Figure 2); then, PH=2. We introduce the second control ν1Ω2,ν3Ω2 given by (10), in order to reduce entry of infected people from region Ω2.In Figure2, we can observe that in the absence of control and in the presence of an epidemic which spreads in three regions characterized by different parameters, the number of susceptible individuals varies from SΩ10=6334,SΩ20=12666, and SΩ30=19000 as initial conditions to 5088, 12094, and 19230, respectively, in the three regions; on the other hand, when the BCG vaccination controls uΩ1uΩ2anduΩ3, given by (11), are introduced in the three regions and we blocked the movement of infected individuals to the Ω2 region by introducing the controls v1Ω2,v3Ω2 given by (10), we observe that the number of susceptible subjects decreases from the initial state in the three regions to 703, 545, and 3556, respectively (see Figure 4).To understand well, we compared each state variable in the three regions without control (see Figure3), and with the introduction of vaccination and blocking the movement (see Figure 5), we can see clearly in Figure 5 that the curves of the susceptible decrease in all regions especially in second region after 2 years. On the other hand, the recovered curves increase markedly after the implementation of the first strategy which proves the impact of these controls.In Figure4, we can clearly see that after two years, there is effectiveness of vaccination on all the regions and especially on the second region controlled by the blockage of displacement, which is seen through the decrease in the curves related to the number of susceptible. At the same time, we see that we have a rapid increase in the number of recovered, especially in the second region.In Figure6(a), it can be seen that the higher the value of the vaccination control (uΩ2), the more the number of susceptible is reduced, and this is clearly seen in Figure 5. On the other hand, in Figure 6(b), we find that when the travel-blocking control takes small values, i.e., when the contact between susceptible and infected individuals is too low, we obtain good results. This can be seen by comparing regions one and three for the number of infected (see Figure 5).


## 3.1.1. An Optimal Control Approach

We are interested in controlling the population of regionΩj. The problem is to minimize the objective functional given by(6)JjνjΩk,uΩj=AjIITΩj−AjRRTΩj+∑k∈PH∑i=1T−1AjIIiΩj−AjRRiΩj+Bk2νijΩk2+Cj2uiΩj2,where Bk > 0, Cj > 0, AjI > 0, and AjR > 0 are the weight constants of controls, the infected, and the removed in region Ωj, respectively. νjΩk=νjΩkk∈PH where νjΩk=ν0jΩk,…,νT−1jΩk and uΩj=u0Ωj,…,uT−1Ωj. Our goal is to minimize the infected and the cost of applying controls and increase the number of recovered in Ωj.In other words, we are seeking optimal controlνjΩk∗ and uΩj∗ such that(7)JjνjΩk∗,uΩj∗=minνjΩk∈Vj,uΩj∈UjJjνjΩk,uΩj,where Vj and Uj are the control sets defined by Vj=νjΩlmesurable0<νminjΩl≤νijΩl≤νmaxjΩl<1,i=0,1,…,T−1,l∈PH. Uj=uΩjmesurable0<uminΩj≤uiΩj≤umaxΩj<1,i=0,1,…,T−1.By using Pontryagin’s maximum principle [37], we derive necessary conditions for our optimal controls. For this purpose, we defined the Hamiltonian as(8)HiΩj=∑k∈PHAjIIiΩj−AjRRiΩj+Bk2νijΩk2+Ck2uiΩj2+λ1,i+1jSiΩj+djNiΩj−∑k=1pνijΩkβjkIiΩkNiΩjSiΩj−djSiΩj−uiΩjSiΩj+λ2,i+1jL1,iΩj+∑k=1pνijΩkβjkIiΩkSiΩjNiΩj+∑k=1pβjkIiΩkNiΩjσjL2,iΩj+σRjRiΩj−δj+τ1j+djL1,iΩj+λ3,i+1jIiΩj+ϕjδjL1,iΩj+ωjL2,iΩj+ωRjRiΩj−τ0j+djIiΩj+λ4,i+1jL2,iΩj+1−ϕjδjL1,iΩj−∑k=1pβjkIiΩkNiΩjσjL2,iΩj−ωj+τ2j+djL2,iΩj+λ5,i+1jRiΩj+τ0jIiΩj+τ1jL1,iΩj+τ2jL2,iΩj−∑k=1pβjkIiΩkNiΩjσRjRiΩj−ω1j+djRiΩj+uiΩjSiΩj.Theorem 1.
Given an optimal controlνjΩk∗,k∈PH and solutions SΩj∗,L1Ωj∗,IΩj∗,L2Ωj∗ and RΩj∗, there exist λs,ij,i=1,2,…,T,j=1,2,…,p,s=1,2,3,4,5. The adjoint variables satisfy the following equations:(9)λ1,ij=1−djλ1,i+1j+∑k=1pνijΩkβjkIiΩkNiΩjλ2,i+1j−λ1,i+1j+uiΩjλ5,i+1j−λ1,i+1j,λ2,ij=1−djλ2,i+1j+δjλ4,i+1j−λ2,i+1j+ϕjδjλ3,i+1j−λ4,i+1j+τ1jλ5,i+1j−λ2,i+1j,λ3,ij=AjI+1−djλ3,i+1j+βjjSiΩjNiΩjλ2,i+1j−λ1,i+1j+σjβjjL2,iΩjNiΩjλ2,i+1j−λ4,i+1j+σRjβjjRiΩjNiΩjλ2,i+1j−λ5,i+1j+τ0jλ5,i+1j−λ3,i+1j,λ4,ij=1−djλ4,i+1j+∑k=1pσjβjkIiΩkNiΩjλ2,i+1j−λ4,i+1j+ωjλ3,i+1j−λ4,i+1j+τ2jλ5,i+1j−λ4,i+1j,λ5,ij=−AjR+1−djλ5,i+1j+∑k=1pσRjβjkIiΩkNiΩjλ2,i+1j−λ5,i+1j+ωRjλ3,i+1j−λ5,i+1j,where λ1,Tj=0;λ2,Tj=0;λ3,Tj=AjI;λ4,Tj=0; and λ5,Tj=−AjR are the transversality conditions.
In addition, we have(10)νijΩl∗=minmaxνminjΩl,βjlBlIiΩl∗SiΩl∗NiΩlλ1,i+1l−λ2,i+1l,νmaxjΩl,i=0,1,…,T−1,l∈PH,(11)uiΩj∗=minmaxuminΩj,SiΩj∗λ1,i+1j−λ5,i+1jCj,umaxΩj,i=0,1,…,T−1.Proof.
Using Pontryagin’s maximum principle [37, 38], we obtain the following adjoint equations:(12)λ1,ij=∂Hi∂SiΩj=1−∑k=1pνijΩkβjkIiΩkNiΩj−dj−uiΩjλ1,i+1j+∑k=1pνijΩkβjkIiΩkNiΩjλ2,i+1j+uiΩjλ5,i+1j,λ2,ij=∂Hi∂L1,iΩj=1−δj−τ1j−djλ2,i+1j+ϕjδjλ3,i+1j+1−ϕjδjλ4,i+1j+τ1jλ5,i+1j,λ3,ij=∂Hi∂IiΩj=AjI−λ1,i+1jβjjSiΩjNiΩj+λ2,i+1jβjj1NiΩj+βjjσjL2,iΩj+σRjRiΩjNiΩj+λ3,i+1j1−τ0j−dj−λ4,i+1jβjjσjL2,iΩjNiΩj+λ5,i+1jτ0j−βjkσRjRiΩjNiΩj,λ4,ij=∂Hi∂L2,iΩj=λ2,i+1j∑k=1pβjkσjIiΩkNiΩj+λ3,i+1jωj+λ4,i+1j1−τ2j−ωj−dj−∑k=1pβjkσjIiΩkNiΩj+λ5,i+1jτ2j,λ5,ij=∂Hi∂RiΩj=−AjR+λ2,i+1j∑k=1pσRjβjkIiΩkNiΩj+λ3,i+1jωRj+λ5,i+1j1−ωRj−dj−∑k=1pβjkσRjIiΩkNiΩj,with λ1,Tj=0,λ2,Tj=0,λ3,Tj=AjI,λ4,Tj=0,λ5,Tj=−AjR,j=1,2,…,p.
To obtain the optimality conditions, we take the variation with respect to controlsνijΩl and uiΩj and set it equal to zero:(13)∂Hi∂νijΩl=BlνijΩl−λ1,i+1jβjlIiΩlNiΩjSiΩj+λ2,i+1jβjlIiΩlNiΩjSiΩj=0,∂Hi∂uiΩj=CjuiΩj−λ1,i+1jSiΩj+λ5,i+1jSiΩj=0,and then we obtain the optimal control pair(14)νijΩl=βjlBlIΩlSiΩjNiΩjλ1,i+1j−λ2,i+1j,uiΩj=SiΩjλ1,i+1j−λ5,i+1jCj.
By the bounds inVj and Uj of the controls, it is easy to obtain νijΩl∗ and uiΩj∗ in the following form:(15)νijΩl∗=minmaxνminjΩl,βjlBlIiΩlSiΩlNiΩlλ1,i+1l−λ2,i+1l,νmaxjΩl,i=0,1,…,T−1,l∈PH,uiΩj∗=minmaxuminΩj,SiΩjλ1,i+1j−λ5,i+1jCj,umaxΩj,i=0,1,…,T−1.


## 3.1.2. Numerical Results

We now present numerical simulations associated with the aforementioned optimal control problem. We write a code inMATLABTM and simulate our results using different data. The optimality systems are solved based on an iterative discrete scheme that converges following an appropriate test similar to the one related to the forward-backward sweep method (FBSM). The state system with an initial guess is solved forward in time and then the adjoint system is solved backward in time because of the transversality conditions. Afterwards, we update the optimal control values using the values of state and costate variables obtained at the previous steps. Finally, we execute the previous steps till a tolerance criterion is reached.The multiregion TB epidemic model we suggest here is applicable for any number of regionsp. For example, in order to show the importance of our work, we choose p=3, i.e., we consider three regions (domains) Ω1,Ω2, and Ω3 with different parameters based on [12] and given in Table 1. Other epidemiological and numerical parameters are presented in Tables 1 and 2, respectively. The initial value of each category, SΩi0=76/120NΩi;L1Ωi0=36/120NΩi;IΩi0=4/120NΩi;L2Ωi0=2/120NΩi; and RΩi0=1/120NΩi, is given in Table 2, which is based on [39] (Figures 2–13).Table 2
Parameter values ofβ utilized for the resolution of multiregion discrete systems (2), (5), and (16), leading to simulations obtained from Figures 2–13, with the initial conditions S0,L10,I0,L20, and R0 associated with the three regions Ω1,Ω2, and Ω3.
NS (0)L10I (0)L20R (0)βΩ11000063343083334166830.05Ω2200001266661686663341660.21Ω33000019000925010005002500.075Figure 2
The state variables in three regions without control.
(a)(b)(c)Figure 3
Comparison for each state variable for the three regions without control.
(a)(b)(c)(d)(e)Figure 4
The state variables in three regions with BCG vaccination and travel-blocking control.
(a)(b)(c)Figure 5
Comparison for each state variable for the three regions with BCG vaccination and travel-blocking control.
(a)(b)(c)(d)(e)Figure 6
The BCG vaccination and travel-blocking control in three regions.
(a)(b)Figure 7
Comparison for each state variable for the three regions with BCG vaccination, travel blocking, and treatment control.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 8
The state variables in three regions with controls vaccination, travel blocking, and treatment.
(a)(b)(c)Figure 9
Comparison of the susceptible in each region without control and with control in two cases.
(a)(b)(c)Figure 10
Comparison of the early latent in each region without control and with control in two cases.
(a)(b)(c)Figure 11
Comparison of the infected in each region without control and with control in two cases.
(a)(b)(c)Figure 12
Comparison of the persistent latent in each region without control and with control in two cases.
(a)(b)(c)Figure 13
Comparison of the recovered in each region without control and with control in two cases.
(a)(b)(c)Since the BCG vaccine in Morocco is mandatory, we considered the BCG vaccination control in the three regions given by (11). We are also interested in controlling Ω1 and Ω3, by assuming that Ω2 is at high risk of infection (see Figure 2); then, PH=2. We introduce the second control ν1Ω2,ν3Ω2 given by (10), in order to reduce entry of infected people from region Ω2.In Figure2, we can observe that in the absence of control and in the presence of an epidemic which spreads in three regions characterized by different parameters, the number of susceptible individuals varies from SΩ10=6334,SΩ20=12666, and SΩ30=19000 as initial conditions to 5088, 12094, and 19230, respectively, in the three regions; on the other hand, when the BCG vaccination controls uΩ1uΩ2anduΩ3, given by (11), are introduced in the three regions and we blocked the movement of infected individuals to the Ω2 region by introducing the controls v1Ω2,v3Ω2 given by (10), we observe that the number of susceptible subjects decreases from the initial state in the three regions to 703, 545, and 3556, respectively (see Figure 4).To understand well, we compared each state variable in the three regions without control (see Figure3), and with the introduction of vaccination and blocking the movement (see Figure 5), we can see clearly in Figure 5 that the curves of the susceptible decrease in all regions especially in second region after 2 years. On the other hand, the recovered curves increase markedly after the implementation of the first strategy which proves the impact of these controls.In Figure4, we can clearly see that after two years, there is effectiveness of vaccination on all the regions and especially on the second region controlled by the blockage of displacement, which is seen through the decrease in the curves related to the number of susceptible. At the same time, we see that we have a rapid increase in the number of recovered, especially in the second region.In Figure6(a), it can be seen that the higher the value of the vaccination control (uΩ2), the more the number of susceptible is reduced, and this is clearly seen in Figure 5. On the other hand, in Figure 6(b), we find that when the travel-blocking control takes small values, i.e., when the contact between susceptible and infected individuals is too low, we obtain good results. This can be seen by comparing regions one and three for the number of infected (see Figure 5).


## 4. The Model with Travel Blocking, Vaccination, and Treatment

### 4.1. Presentation of the Model

We are now going to introduce a third control variablewiΩj that characterizes the treatment of the infected in order to reduce them.Then, for a given regionΩj, the model is given as follows:(16)Si+1Ωj=SiΩj+djNiΩj−∑k=1pνijΩkβjkIiΩkNiΩjSiΩj−djSiΩj−uiΩjSiΩj,L1,i+1Ωj=L1,iΩj+∑k=1pνijΩkβjkIiΩkSiΩjNiΩj+∑k=1pβjkIiΩkNiΩjσjL2,iΩj+σRjRiΩj−δj+τ1j+djL1,iΩj,Ii+1Ωj=IiΩj+ϕjδjL1,iΩj+ωjL2,iΩj+ωRjRiΩj−τ0j+djIiΩj−wiΩjIiΩj,L2,i+1Ωj=L2,iΩj+1−ϕjδjL1,iΩj−∑k=1pβjkIiΩkNiΩjσjL2,iΩj−ωj+τ2j+djL2,iΩj,Ri+1Ωj=RiΩj+τ0jIiΩj+τ1jL1,iΩj+τ2jL2,iΩj−∑k=1pβjkIiΩkNiΩjσRjRiΩj−ωRj+djRiΩj+uiΩjSiΩj+wiΩjIiΩj.


#### 4.1.1. An Optimal Control Approach

We are interested in controlling the population of regionΩj. The problem is to minimize the objective functional given by(17)JjνjΩk,uΩj=AjIITΩj−AjRRTΩj+∑k∈PH∑i=1T−1AjIIiΩj−AjRRiΩj+Bk2νijΩk2+Cj2uiΩj2+Dj2wiΩj2, where Bk > 0, Cj > 0, Dj > 0, AJI>0, and AJR>0 are the weight constants of controls, the infected, and the removed in region νijΩk=ν0jΩk,…,νT−1jΩk, uΩj=u0Ωj,…,uT−1Ωj, and wΩj=w0Ωj,…,wT−1Ωj.Our goal is to minimize the infected and the cost of applying controls and increase the number of recovered inΩj.In other words, we are seeking optimal controlνjΩk∗ and uΩj∗ and wΩj∗ such that(18)JjνjΩk∗,uΩj∗,wΩj∗=minνjΩk∈Vj,uΩj∈Uj,wΩj∈WjJjνjΩk,uΩj,wΩj,where Vj,Uj, and Wj are the control sets defined by(19)Vj=νjΩlmesurable /0<νminjΩl≤νijΩl≤νmaxjΩl<1,i=0,1,…,T−1,l∈PH,Uj=uΩjmesurable /0<uminΩj≤uiΩj≤umaxΩj<1,i=0,1,…,T−1,Wj=wΩjmesurable /0<wminΩj≤wiΩj≤wmaxΩj<1,i=0,1,…,T−1.By using Pontryagin’s maximum principle [37], we derive necessary conditions for our optimal controls. For this purpose, we defined the Hamiltonian as(20)HiΩj=∑k∈PHAjIIiΩj−AjRRiΩj+Bk2νijΩk2+Cj2uiΩj2+Dj2wiΩj2+λ1,i+1jSiΩj+djNiΩj−∑k=1pνijΩkβjkIiΩkNiΩjSiΩj−djSiΩj−uiΩjSiΩj+λ2,i+1jL1,iΩj+∑k=1pνijΩkβjkIiΩkSiΩjNiΩj+∑k=1pβjkIiΩkNiΩjσjL2,iΩj+σRjRiΩj−δj+τ1j+djL1,iΩj+λ3,i+1jIiΩj+ϕjδjL1,iΩj+ωjL2,iΩj+ωRjRiΩj−τ0j+dj−wiΩjIiΩj+λ4,i+1jL2,iΩj+1−ϕjδjL1,iΩj−∑k=1pβjkIiΩkNiΩjσjL2,iΩj−ωj+τ2j+djL2,iΩj+λ5,i+1jRiΩj+τ0jIiΩj+τ1jL1,iΩj+τ2jL2,iΩj−∑k=1pβjkIiΩkNiΩjσRjRiΩj−ω1j+djRiΩj+uiΩjSiΩj+wiΩjIiΩj.Theorem 2.
Given optimal controlsνjΩk∗,k∈PH,uΩj∗,wΩj∗ and solutions SΩj∗,L1Ωj∗,IΩj∗,L2Ωj∗ and RΩj∗, there exist λs,ij,fori=1,2,…,T,j=1,2,…,pands=1,2,3,4,5. The adjoint variables satisfy the following equations:(21)λ1,ij=1−djλ1,i+1j+∑k=1pνijΩkβjkIiΩkNiΩjλ2,i+1j−λ1,i+1j+uiΩjλ5,i+1j−λ1,i+1j,λ2,ij=1−djλ2,i+1j+δjλ4,i+1j−λ2,i+1j+ϕjδjλ3,i+1j−λ4,i+1j+τ1jλ5,i+1j−λ2,i+1j,λ3,ij=AjI+1−djλ3,i+1j+βjjIiΩjNiΩjλ2,i+1j−λ1,i+1j+σjβjjL2,iΩjNiΩjλ2,i+1j−λ4,i+1j+σRjβjjRiΩjNiΩjλ2,i+1j−λ5,i+1j+wiΩjλ5,i+1j−λ3,i+1j+τ0jλ5,i+1j−λ3,i+1j,λ4,ij=1−djλ4,i+1j+∑k=1pσjβjkIiΩkNiΩjλ2,i+1j−λ4,i+1j+ωjλ3,i+1j−λ4,i+1j+τ2jλ5,i+1j−λ4,i+1j,λ5,ij=−AjR+1−djλ5,i+1j+∑k=1pσRjβjkIiΩkNiΩjλ2,i+1j−λ5,i+1j+ωRjλ3,i+1j−λ5,i+1j,where λ1,Tj=0;λ2,Tj=0;λ3,Tj=AjI;λ4,Tj=0; and λ5,Tj=−AjR are the transversality conditions.
In addition, we have(22)νijΩl∗=minmaxνminjΩl,βjlBlIiΩl∗SiΩl∗NiΩlλ1,i+1l−λ2,i+1l,νmaxjΩl,i=0,1,…,T−1,l∈PH,uiΩj∗=minmaxuminΩj,SiΩj∗λ1,i+1j−λ5,i+1jCj,umaxΩj,i=0,1,…,T−1,wiΩj∗=minmaxwminΩj,IiΩj∗λ3,i+1j−λ5,i+1jDj,wmaxΩj,i=0,1,…,T−1.Proof.
Using Pontryagin’s maximum principle [37], we obtain the following adjoint equations:(23)λ1,ij=∂Hi∂SiΩj=1−∑k=1pνijΩkβjkIiΩkNiΩj−dj−uiΩjλ1,i+1j+∑k=1pνijΩkβjkIiΩkNiΩjλ2,i+1j+uiΩjλ5,i+1j,λ2,ij=∂Hi∂L1,iΩj=1−δj−τ1j−djλ2,i+1j+ϕjδjλ3,i+1j+1−ϕjδjλ4,i+1j+τ1jλ5,i+1j,λ3,ij=∂Hi∂IiΩj=AjI−λ1,i+1jβjjSiΩjNiΩj+λ2,i+1jβjjSiΩjNiΩj+βjjσjL2,iΩj+σRjRiΩjNiΩj+λ3,i+1j1−τ0j−dj−wiΩj−λ4,i+1jβjjσjL2,iΩjNiΩj+λ5,i+1jτ0j−βjjσRjRiΩjNiΩj+wiΩj,λ4,ij=∂Hi∂L2,iΩj=λ2,i+1j∑k=1pβjkσjIiΩkNiΩj+λ3,i+1jωj+λ4,i+1j1−τ2j−ωj−dj−∑k=1pβjkσjIiΩkNiΩj+λ5,i+1jτ2j,λ5,ij=∂Hi∂RiΩj=−AjR+λ2,i+1j∑k=1pσRjβjkIiΩkNiΩj+λ3,i+1jωRj+λ5,i+1j1−ωRj−dj−∑k=1pβjkσRjIiΩkNiΩj,with λ1,Tj=0,λ2,Tj=0,λ3,Tj=AjI,λ4,Tj=0,λ5,Tj=−AjR,j=1,2,…,p.
To obtain the optimality conditions, we take the variation with respect to controlsνijΩl,uiΩj, and wΩj and set it equal to zero:(24)∂Hi∂νijΩl=BlνijΩl−λ1,i+1jβjlIiΩlNiΩjSiΩj+λ2,i+1jβjlIiΩlNiΩjSiΩj=0,∂Hi∂uiΩj=CjuiΩj−λ1,i+1jSiΩj+λ5,i+1jSiΩj=0,∂Hi∂wiΩj=DjwiΩj−λ3,i+1jIiΩj+λ5,i+1jIiΩj=0,and then we obtain the optimal controls(25)νijΩl=βjlBlIΩlSiΩjNiΩjλ1,i+1j−λ2,i+1j,uiΩl=SiΩjλ1,i+1j−λ5,i+1jCj,wiΩl=IiΩjλ3,i+1j−λ5,i+1jDj.
By the bounds inVj,Uj, and Wj of the controls, it is easy to obtain νijΩl∗,uiΩj∗, and wΩj∗ in the following form:(26)νijΩl∗=minmaxνminjΩl,βjlBlIiΩlSiΩlNiΩlλ1,i+1l−λ2,i+1l,νmaxjΩl,i=0,1,…,T−1,l∈PH,uiΩj∗=minmaxuminΩj,SiΩjλ1,i+1j−λ5,i+1jCj,umaxΩj,i=0,1,…,T−1,wiΩj∗=minmaxwminΩj,IiΩjλ3,i+1j−λ5,i+1jDj,wmaxΩj,i=0,1,…,T−1.


#### 4.1.2. Numerical Results

In this section, we will also introduce the third control (a very important control); it is the treatment control, and of course we force all the infected to take anti-tuberculosis drugs; it is observed, from these subfigures, that using this control strategy, there is a significant decrease in the number of infected individuals treated compared to those untreated, especially in the second region (Figure7).In Figure7, we can clearly see that after first year of disease spread, there is effectiveness of treatment in all the regions and especially in the second region controlled by the blockage of displacement. We observe that the number of infected individuals increases from IΩ10=334,IΩ20=666, and IΩ30=1000 as initial conditions to 347, 691, and 1036 after two years, respectively, in the three regions, in the absence of controls; on the other hand, when we introduced the treatment control in the three regions, the number of infected individuals decreases to 124, 114, and 269 after 2 years, respectively. This is a direct result of the introduction of the treatment in the three regions, and for that we notice in the second region that the red curve reaches its maximum before two years.In Figure8, we notice that the number of susceptible decreases from SΩ10=6334,SΩ20=12666, and SΩ30=19000 as initial conditions to 215, 296, and 833 in the three regions, respectively; it is because of the control effect especially in the second high-risk region. Similarly, the number of individuals cured in the second region increases from 166 people without control to 15267 with the effect of vaccination and blockage of movement and finally treatment. The impact of the introduction of the treatment has given good results with the decrease in the number of infected individuals.In Figure9, we present the populations of susceptible in each region in three cases without control (blue curves) and with control in two strategies: with vaccination and travel blocking (red curves) and with vaccination, travel blocking, and treatment (green curves); we see that the number of susceptible decreases quickly after implementing the two strategies.In Figure10, we present the populations of early latent in each region in three cases without control (blue curves) and with control in two strategies: with vaccination and travel blocking (red curves) and with vaccination, travel blocking, and treatment (green curves); we see that the number of early latent decreases, and of course, the best strategy is the second.In Figure11, when we compare the infected in each region without control and with control in two strategies, we clearly see that the green curves decrease rapidly to zero. This is due to the effect of the treatment on the infected individuals.In Figure12, the comparison of persistent latent shows us the efficiency of the introduction of controls; it is therefore seen that the infected with long controlled duration does not increase like those uncontrolled.In Figure13, by comparing the curves corresponding to the recovered for each region, we can see that there is a significant difference between the individuals recovered with and without control; this is due to the impact of the control strategies used, and for the second region, we have a large increase due to the blocking of movement towards this last region and the majority of those infected are treated.In conclusion, when we compare the curves obtained in the numerical simulations, we see that the strategy which includes vaccination, treatment, and blocking of trips is more effective than the first strategy (vaccination and blocking of trips).


## 4.1. Presentation of the Model

We are now going to introduce a third control variablewiΩj that characterizes the treatment of the infected in order to reduce them.Then, for a given regionΩj, the model is given as follows:(16)Si+1Ωj=SiΩj+djNiΩj−∑k=1pνijΩkβjkIiΩkNiΩjSiΩj−djSiΩj−uiΩjSiΩj,L1,i+1Ωj=L1,iΩj+∑k=1pνijΩkβjkIiΩkSiΩjNiΩj+∑k=1pβjkIiΩkNiΩjσjL2,iΩj+σRjRiΩj−δj+τ1j+djL1,iΩj,Ii+1Ωj=IiΩj+ϕjδjL1,iΩj+ωjL2,iΩj+ωRjRiΩj−τ0j+djIiΩj−wiΩjIiΩj,L2,i+1Ωj=L2,iΩj+1−ϕjδjL1,iΩj−∑k=1pβjkIiΩkNiΩjσjL2,iΩj−ωj+τ2j+djL2,iΩj,Ri+1Ωj=RiΩj+τ0jIiΩj+τ1jL1,iΩj+τ2jL2,iΩj−∑k=1pβjkIiΩkNiΩjσRjRiΩj−ωRj+djRiΩj+uiΩjSiΩj+wiΩjIiΩj.


### 4.1.1. An Optimal Control Approach

We are interested in controlling the population of regionΩj. The problem is to minimize the objective functional given by(17)JjνjΩk,uΩj=AjIITΩj−AjRRTΩj+∑k∈PH∑i=1T−1AjIIiΩj−AjRRiΩj+Bk2νijΩk2+Cj2uiΩj2+Dj2wiΩj2, where Bk > 0, Cj > 0, Dj > 0, AJI>0, and AJR>0 are the weight constants of controls, the infected, and the removed in region νijΩk=ν0jΩk,…,νT−1jΩk, uΩj=u0Ωj,…,uT−1Ωj, and wΩj=w0Ωj,…,wT−1Ωj.Our goal is to minimize the infected and the cost of applying controls and increase the number of recovered inΩj.In other words, we are seeking optimal controlνjΩk∗ and uΩj∗ and wΩj∗ such that(18)JjνjΩk∗,uΩj∗,wΩj∗=minνjΩk∈Vj,uΩj∈Uj,wΩj∈WjJjνjΩk,uΩj,wΩj,where Vj,Uj, and Wj are the control sets defined by(19)Vj=νjΩlmesurable /0<νminjΩl≤νijΩl≤νmaxjΩl<1,i=0,1,…,T−1,l∈PH,Uj=uΩjmesurable /0<uminΩj≤uiΩj≤umaxΩj<1,i=0,1,…,T−1,Wj=wΩjmesurable /0<wminΩj≤wiΩj≤wmaxΩj<1,i=0,1,…,T−1.By using Pontryagin’s maximum principle [37], we derive necessary conditions for our optimal controls. For this purpose, we defined the Hamiltonian as(20)HiΩj=∑k∈PHAjIIiΩj−AjRRiΩj+Bk2νijΩk2+Cj2uiΩj2+Dj2wiΩj2+λ1,i+1jSiΩj+djNiΩj−∑k=1pνijΩkβjkIiΩkNiΩjSiΩj−djSiΩj−uiΩjSiΩj+λ2,i+1jL1,iΩj+∑k=1pνijΩkβjkIiΩkSiΩjNiΩj+∑k=1pβjkIiΩkNiΩjσjL2,iΩj+σRjRiΩj−δj+τ1j+djL1,iΩj+λ3,i+1jIiΩj+ϕjδjL1,iΩj+ωjL2,iΩj+ωRjRiΩj−τ0j+dj−wiΩjIiΩj+λ4,i+1jL2,iΩj+1−ϕjδjL1,iΩj−∑k=1pβjkIiΩkNiΩjσjL2,iΩj−ωj+τ2j+djL2,iΩj+λ5,i+1jRiΩj+τ0jIiΩj+τ1jL1,iΩj+τ2jL2,iΩj−∑k=1pβjkIiΩkNiΩjσRjRiΩj−ω1j+djRiΩj+uiΩjSiΩj+wiΩjIiΩj.Theorem 2.
Given optimal controlsνjΩk∗,k∈PH,uΩj∗,wΩj∗ and solutions SΩj∗,L1Ωj∗,IΩj∗,L2Ωj∗ and RΩj∗, there exist λs,ij,fori=1,2,…,T,j=1,2,…,pands=1,2,3,4,5. The adjoint variables satisfy the following equations:(21)λ1,ij=1−djλ1,i+1j+∑k=1pνijΩkβjkIiΩkNiΩjλ2,i+1j−λ1,i+1j+uiΩjλ5,i+1j−λ1,i+1j,λ2,ij=1−djλ2,i+1j+δjλ4,i+1j−λ2,i+1j+ϕjδjλ3,i+1j−λ4,i+1j+τ1jλ5,i+1j−λ2,i+1j,λ3,ij=AjI+1−djλ3,i+1j+βjjIiΩjNiΩjλ2,i+1j−λ1,i+1j+σjβjjL2,iΩjNiΩjλ2,i+1j−λ4,i+1j+σRjβjjRiΩjNiΩjλ2,i+1j−λ5,i+1j+wiΩjλ5,i+1j−λ3,i+1j+τ0jλ5,i+1j−λ3,i+1j,λ4,ij=1−djλ4,i+1j+∑k=1pσjβjkIiΩkNiΩjλ2,i+1j−λ4,i+1j+ωjλ3,i+1j−λ4,i+1j+τ2jλ5,i+1j−λ4,i+1j,λ5,ij=−AjR+1−djλ5,i+1j+∑k=1pσRjβjkIiΩkNiΩjλ2,i+1j−λ5,i+1j+ωRjλ3,i+1j−λ5,i+1j,where λ1,Tj=0;λ2,Tj=0;λ3,Tj=AjI;λ4,Tj=0; and λ5,Tj=−AjR are the transversality conditions.
In addition, we have(22)νijΩl∗=minmaxνminjΩl,βjlBlIiΩl∗SiΩl∗NiΩlλ1,i+1l−λ2,i+1l,νmaxjΩl,i=0,1,…,T−1,l∈PH,uiΩj∗=minmaxuminΩj,SiΩj∗λ1,i+1j−λ5,i+1jCj,umaxΩj,i=0,1,…,T−1,wiΩj∗=minmaxwminΩj,IiΩj∗λ3,i+1j−λ5,i+1jDj,wmaxΩj,i=0,1,…,T−1.Proof.
Using Pontryagin’s maximum principle [37], we obtain the following adjoint equations:(23)λ1,ij=∂Hi∂SiΩj=1−∑k=1pνijΩkβjkIiΩkNiΩj−dj−uiΩjλ1,i+1j+∑k=1pνijΩkβjkIiΩkNiΩjλ2,i+1j+uiΩjλ5,i+1j,λ2,ij=∂Hi∂L1,iΩj=1−δj−τ1j−djλ2,i+1j+ϕjδjλ3,i+1j+1−ϕjδjλ4,i+1j+τ1jλ5,i+1j,λ3,ij=∂Hi∂IiΩj=AjI−λ1,i+1jβjjSiΩjNiΩj+λ2,i+1jβjjSiΩjNiΩj+βjjσjL2,iΩj+σRjRiΩjNiΩj+λ3,i+1j1−τ0j−dj−wiΩj−λ4,i+1jβjjσjL2,iΩjNiΩj+λ5,i+1jτ0j−βjjσRjRiΩjNiΩj+wiΩj,λ4,ij=∂Hi∂L2,iΩj=λ2,i+1j∑k=1pβjkσjIiΩkNiΩj+λ3,i+1jωj+λ4,i+1j1−τ2j−ωj−dj−∑k=1pβjkσjIiΩkNiΩj+λ5,i+1jτ2j,λ5,ij=∂Hi∂RiΩj=−AjR+λ2,i+1j∑k=1pσRjβjkIiΩkNiΩj+λ3,i+1jωRj+λ5,i+1j1−ωRj−dj−∑k=1pβjkσRjIiΩkNiΩj,with λ1,Tj=0,λ2,Tj=0,λ3,Tj=AjI,λ4,Tj=0,λ5,Tj=−AjR,j=1,2,…,p.
To obtain the optimality conditions, we take the variation with respect to controlsνijΩl,uiΩj, and wΩj and set it equal to zero:(24)∂Hi∂νijΩl=BlνijΩl−λ1,i+1jβjlIiΩlNiΩjSiΩj+λ2,i+1jβjlIiΩlNiΩjSiΩj=0,∂Hi∂uiΩj=CjuiΩj−λ1,i+1jSiΩj+λ5,i+1jSiΩj=0,∂Hi∂wiΩj=DjwiΩj−λ3,i+1jIiΩj+λ5,i+1jIiΩj=0,and then we obtain the optimal controls(25)νijΩl=βjlBlIΩlSiΩjNiΩjλ1,i+1j−λ2,i+1j,uiΩl=SiΩjλ1,i+1j−λ5,i+1jCj,wiΩl=IiΩjλ3,i+1j−λ5,i+1jDj.
By the bounds inVj,Uj, and Wj of the controls, it is easy to obtain νijΩl∗,uiΩj∗, and wΩj∗ in the following form:(26)νijΩl∗=minmaxνminjΩl,βjlBlIiΩlSiΩlNiΩlλ1,i+1l−λ2,i+1l,νmaxjΩl,i=0,1,…,T−1,l∈PH,uiΩj∗=minmaxuminΩj,SiΩjλ1,i+1j−λ5,i+1jCj,umaxΩj,i=0,1,…,T−1,wiΩj∗=minmaxwminΩj,IiΩjλ3,i+1j−λ5,i+1jDj,wmaxΩj,i=0,1,…,T−1.


### 4.1.2. Numerical Results

In this section, we will also introduce the third control (a very important control); it is the treatment control, and of course we force all the infected to take anti-tuberculosis drugs; it is observed, from these subfigures, that using this control strategy, there is a significant decrease in the number of infected individuals treated compared to those untreated, especially in the second region (Figure7).In Figure7, we can clearly see that after first year of disease spread, there is effectiveness of treatment in all the regions and especially in the second region controlled by the blockage of displacement. We observe that the number of infected individuals increases from IΩ10=334,IΩ20=666, and IΩ30=1000 as initial conditions to 347, 691, and 1036 after two years, respectively, in the three regions, in the absence of controls; on the other hand, when we introduced the treatment control in the three regions, the number of infected individuals decreases to 124, 114, and 269 after 2 years, respectively. This is a direct result of the introduction of the treatment in the three regions, and for that we notice in the second region that the red curve reaches its maximum before two years.In Figure8, we notice that the number of susceptible decreases from SΩ10=6334,SΩ20=12666, and SΩ30=19000 as initial conditions to 215, 296, and 833 in the three regions, respectively; it is because of the control effect especially in the second high-risk region. Similarly, the number of individuals cured in the second region increases from 166 people without control to 15267 with the effect of vaccination and blockage of movement and finally treatment. The impact of the introduction of the treatment has given good results with the decrease in the number of infected individuals.In Figure9, we present the populations of susceptible in each region in three cases without control (blue curves) and with control in two strategies: with vaccination and travel blocking (red curves) and with vaccination, travel blocking, and treatment (green curves); we see that the number of susceptible decreases quickly after implementing the two strategies.In Figure10, we present the populations of early latent in each region in three cases without control (blue curves) and with control in two strategies: with vaccination and travel blocking (red curves) and with vaccination, travel blocking, and treatment (green curves); we see that the number of early latent decreases, and of course, the best strategy is the second.In Figure11, when we compare the infected in each region without control and with control in two strategies, we clearly see that the green curves decrease rapidly to zero. This is due to the effect of the treatment on the infected individuals.In Figure12, the comparison of persistent latent shows us the efficiency of the introduction of controls; it is therefore seen that the infected with long controlled duration does not increase like those uncontrolled.In Figure13, by comparing the curves corresponding to the recovered for each region, we can see that there is a significant difference between the individuals recovered with and without control; this is due to the impact of the control strategies used, and for the second region, we have a large increase due to the blocking of movement towards this last region and the majority of those infected are treated.In conclusion, when we compare the curves obtained in the numerical simulations, we see that the strategy which includes vaccination, treatment, and blocking of trips is more effective than the first strategy (vaccination and blocking of trips).


## 4.1.1. An Optimal Control Approach

We are interested in controlling the population of regionΩj. The problem is to minimize the objective functional given by(17)JjνjΩk,uΩj=AjIITΩj−AjRRTΩj+∑k∈PH∑i=1T−1AjIIiΩj−AjRRiΩj+Bk2νijΩk2+Cj2uiΩj2+Dj2wiΩj2, where Bk > 0, Cj > 0, Dj > 0, AJI>0, and AJR>0 are the weight constants of controls, the infected, and the removed in region νijΩk=ν0jΩk,…,νT−1jΩk, uΩj=u0Ωj,…,uT−1Ωj, and wΩj=w0Ωj,…,wT−1Ωj.Our goal is to minimize the infected and the cost of applying controls and increase the number of recovered inΩj.In other words, we are seeking optimal controlνjΩk∗ and uΩj∗ and wΩj∗ such that(18)JjνjΩk∗,uΩj∗,wΩj∗=minνjΩk∈Vj,uΩj∈Uj,wΩj∈WjJjνjΩk,uΩj,wΩj,where Vj,Uj, and Wj are the control sets defined by(19)Vj=νjΩlmesurable /0<νminjΩl≤νijΩl≤νmaxjΩl<1,i=0,1,…,T−1,l∈PH,Uj=uΩjmesurable /0<uminΩj≤uiΩj≤umaxΩj<1,i=0,1,…,T−1,Wj=wΩjmesurable /0<wminΩj≤wiΩj≤wmaxΩj<1,i=0,1,…,T−1.By using Pontryagin’s maximum principle [37], we derive necessary conditions for our optimal controls. For this purpose, we defined the Hamiltonian as(20)HiΩj=∑k∈PHAjIIiΩj−AjRRiΩj+Bk2νijΩk2+Cj2uiΩj2+Dj2wiΩj2+λ1,i+1jSiΩj+djNiΩj−∑k=1pνijΩkβjkIiΩkNiΩjSiΩj−djSiΩj−uiΩjSiΩj+λ2,i+1jL1,iΩj+∑k=1pνijΩkβjkIiΩkSiΩjNiΩj+∑k=1pβjkIiΩkNiΩjσjL2,iΩj+σRjRiΩj−δj+τ1j+djL1,iΩj+λ3,i+1jIiΩj+ϕjδjL1,iΩj+ωjL2,iΩj+ωRjRiΩj−τ0j+dj−wiΩjIiΩj+λ4,i+1jL2,iΩj+1−ϕjδjL1,iΩj−∑k=1pβjkIiΩkNiΩjσjL2,iΩj−ωj+τ2j+djL2,iΩj+λ5,i+1jRiΩj+τ0jIiΩj+τ1jL1,iΩj+τ2jL2,iΩj−∑k=1pβjkIiΩkNiΩjσRjRiΩj−ω1j+djRiΩj+uiΩjSiΩj+wiΩjIiΩj.Theorem 2.
Given optimal controlsνjΩk∗,k∈PH,uΩj∗,wΩj∗ and solutions SΩj∗,L1Ωj∗,IΩj∗,L2Ωj∗ and RΩj∗, there exist λs,ij,fori=1,2,…,T,j=1,2,…,pands=1,2,3,4,5. The adjoint variables satisfy the following equations:(21)λ1,ij=1−djλ1,i+1j+∑k=1pνijΩkβjkIiΩkNiΩjλ2,i+1j−λ1,i+1j+uiΩjλ5,i+1j−λ1,i+1j,λ2,ij=1−djλ2,i+1j+δjλ4,i+1j−λ2,i+1j+ϕjδjλ3,i+1j−λ4,i+1j+τ1jλ5,i+1j−λ2,i+1j,λ3,ij=AjI+1−djλ3,i+1j+βjjIiΩjNiΩjλ2,i+1j−λ1,i+1j+σjβjjL2,iΩjNiΩjλ2,i+1j−λ4,i+1j+σRjβjjRiΩjNiΩjλ2,i+1j−λ5,i+1j+wiΩjλ5,i+1j−λ3,i+1j+τ0jλ5,i+1j−λ3,i+1j,λ4,ij=1−djλ4,i+1j+∑k=1pσjβjkIiΩkNiΩjλ2,i+1j−λ4,i+1j+ωjλ3,i+1j−λ4,i+1j+τ2jλ5,i+1j−λ4,i+1j,λ5,ij=−AjR+1−djλ5,i+1j+∑k=1pσRjβjkIiΩkNiΩjλ2,i+1j−λ5,i+1j+ωRjλ3,i+1j−λ5,i+1j,where λ1,Tj=0;λ2,Tj=0;λ3,Tj=AjI;λ4,Tj=0; and λ5,Tj=−AjR are the transversality conditions.
In addition, we have(22)νijΩl∗=minmaxνminjΩl,βjlBlIiΩl∗SiΩl∗NiΩlλ1,i+1l−λ2,i+1l,νmaxjΩl,i=0,1,…,T−1,l∈PH,uiΩj∗=minmaxuminΩj,SiΩj∗λ1,i+1j−λ5,i+1jCj,umaxΩj,i=0,1,…,T−1,wiΩj∗=minmaxwminΩj,IiΩj∗λ3,i+1j−λ5,i+1jDj,wmaxΩj,i=0,1,…,T−1.Proof.
Using Pontryagin’s maximum principle [37], we obtain the following adjoint equations:(23)λ1,ij=∂Hi∂SiΩj=1−∑k=1pνijΩkβjkIiΩkNiΩj−dj−uiΩjλ1,i+1j+∑k=1pνijΩkβjkIiΩkNiΩjλ2,i+1j+uiΩjλ5,i+1j,λ2,ij=∂Hi∂L1,iΩj=1−δj−τ1j−djλ2,i+1j+ϕjδjλ3,i+1j+1−ϕjδjλ4,i+1j+τ1jλ5,i+1j,λ3,ij=∂Hi∂IiΩj=AjI−λ1,i+1jβjjSiΩjNiΩj+λ2,i+1jβjjSiΩjNiΩj+βjjσjL2,iΩj+σRjRiΩjNiΩj+λ3,i+1j1−τ0j−dj−wiΩj−λ4,i+1jβjjσjL2,iΩjNiΩj+λ5,i+1jτ0j−βjjσRjRiΩjNiΩj+wiΩj,λ4,ij=∂Hi∂L2,iΩj=λ2,i+1j∑k=1pβjkσjIiΩkNiΩj+λ3,i+1jωj+λ4,i+1j1−τ2j−ωj−dj−∑k=1pβjkσjIiΩkNiΩj+λ5,i+1jτ2j,λ5,ij=∂Hi∂RiΩj=−AjR+λ2,i+1j∑k=1pσRjβjkIiΩkNiΩj+λ3,i+1jωRj+λ5,i+1j1−ωRj−dj−∑k=1pβjkσRjIiΩkNiΩj,with λ1,Tj=0,λ2,Tj=0,λ3,Tj=AjI,λ4,Tj=0,λ5,Tj=−AjR,j=1,2,…,p.
To obtain the optimality conditions, we take the variation with respect to controlsνijΩl,uiΩj, and wΩj and set it equal to zero:(24)∂Hi∂νijΩl=BlνijΩl−λ1,i+1jβjlIiΩlNiΩjSiΩj+λ2,i+1jβjlIiΩlNiΩjSiΩj=0,∂Hi∂uiΩj=CjuiΩj−λ1,i+1jSiΩj+λ5,i+1jSiΩj=0,∂Hi∂wiΩj=DjwiΩj−λ3,i+1jIiΩj+λ5,i+1jIiΩj=0,and then we obtain the optimal controls(25)νijΩl=βjlBlIΩlSiΩjNiΩjλ1,i+1j−λ2,i+1j,uiΩl=SiΩjλ1,i+1j−λ5,i+1jCj,wiΩl=IiΩjλ3,i+1j−λ5,i+1jDj.
By the bounds inVj,Uj, and Wj of the controls, it is easy to obtain νijΩl∗,uiΩj∗, and wΩj∗ in the following form:(26)νijΩl∗=minmaxνminjΩl,βjlBlIiΩlSiΩlNiΩlλ1,i+1l−λ2,i+1l,νmaxjΩl,i=0,1,…,T−1,l∈PH,uiΩj∗=minmaxuminΩj,SiΩjλ1,i+1j−λ5,i+1jCj,umaxΩj,i=0,1,…,T−1,wiΩj∗=minmaxwminΩj,IiΩjλ3,i+1j−λ5,i+1jDj,wmaxΩj,i=0,1,…,T−1.


## 4.1.2. Numerical Results

In this section, we will also introduce the third control (a very important control); it is the treatment control, and of course we force all the infected to take anti-tuberculosis drugs; it is observed, from these subfigures, that using this control strategy, there is a significant decrease in the number of infected individuals treated compared to those untreated, especially in the second region (Figure7).In Figure7, we can clearly see that after first year of disease spread, there is effectiveness of treatment in all the regions and especially in the second region controlled by the blockage of displacement. We observe that the number of infected individuals increases from IΩ10=334,IΩ20=666, and IΩ30=1000 as initial conditions to 347, 691, and 1036 after two years, respectively, in the three regions, in the absence of controls; on the other hand, when we introduced the treatment control in the three regions, the number of infected individuals decreases to 124, 114, and 269 after 2 years, respectively. This is a direct result of the introduction of the treatment in the three regions, and for that we notice in the second region that the red curve reaches its maximum before two years.In Figure8, we notice that the number of susceptible decreases from SΩ10=6334,SΩ20=12666, and SΩ30=19000 as initial conditions to 215, 296, and 833 in the three regions, respectively; it is because of the control effect especially in the second high-risk region. Similarly, the number of individuals cured in the second region increases from 166 people without control to 15267 with the effect of vaccination and blockage of movement and finally treatment. The impact of the introduction of the treatment has given good results with the decrease in the number of infected individuals.In Figure9, we present the populations of susceptible in each region in three cases without control (blue curves) and with control in two strategies: with vaccination and travel blocking (red curves) and with vaccination, travel blocking, and treatment (green curves); we see that the number of susceptible decreases quickly after implementing the two strategies.In Figure10, we present the populations of early latent in each region in three cases without control (blue curves) and with control in two strategies: with vaccination and travel blocking (red curves) and with vaccination, travel blocking, and treatment (green curves); we see that the number of early latent decreases, and of course, the best strategy is the second.In Figure11, when we compare the infected in each region without control and with control in two strategies, we clearly see that the green curves decrease rapidly to zero. This is due to the effect of the treatment on the infected individuals.In Figure12, the comparison of persistent latent shows us the efficiency of the introduction of controls; it is therefore seen that the infected with long controlled duration does not increase like those uncontrolled.In Figure13, by comparing the curves corresponding to the recovered for each region, we can see that there is a significant difference between the individuals recovered with and without control; this is due to the impact of the control strategies used, and for the second region, we have a large increase due to the blocking of movement towards this last region and the majority of those infected are treated.In conclusion, when we compare the curves obtained in the numerical simulations, we see that the strategy which includes vaccination, treatment, and blocking of trips is more effective than the first strategy (vaccination and blocking of trips).


## 5. Conclusion

Tuberculosis represents a considerable burden on health services all over the world because the process is long and very expensive, so it is impossible to immunize the entire population, especially for developing countries. In this case, it is necessary to target the regions most likely to catch tuberculosis, i.e., “regions of fragility”; in this way, we will reduce the cost and the effort deployed by public authorities. This is the reason why we chose to work on multiregion models. This paper analyzes two different optimal control issues for a TB multiregion discrete-time model. In the first part, we introduce two controls that characterize the effectiveness of vaccination and blocking trips, where we vaccinated some of the susceptible people in all regions; except that in the high-risk region, we vaccinated as many susceptible people as possible, in order to reduce the number of susceptible and increase the number of recovered; at the same time, the displacement to high-risk areas was prohibited by controlling contact between susceptible people and infected individuals at high risk of infection. The numerical results associated with this first part assured us of the effectiveness of the vaccination and the impact of the travel blocking on reducing the spread of the disease from one region to another. We extended the first part to a second part by adding a processing control for the infected individuals, and above all, we concentrated the efforts on the high-risk region and we treated most of the individuals in this region. Our numerical results show that this strategy (the strategy that combines vaccination, treatment, and blocking trips) is more effective than the first strategy (see Figure11). In conclusion, we have shown the interest of each case studied by a comparison between the numerical simulations provided. We have demonstrated the importance and effectiveness of our approach in controlling the spread of infection. Control programs based on these strategies can effectively reduce the number of infected cases and increase the number of individuals removed from a specific region. As a perspective of this work and in order to make our model more realistic, we will introduce the notion of delay to model the latency period.



---

*Source: 5594778-2021-06-07.xml*