AN AUTOMATED OPTIMAL VACCINATION CONTROL WITH A MULTI-REGION SIR EPIDEMIC MODEL

1Department of Mathematics and Computer Science, Poly-disciplinary Faculty, Cadi Ayyad University, P.O. Box 4162, Safi, Morocco 2Laboratory of Analysis Modelling and Simulation, Department of Mathematics and Computer Science, Faculty of Sciences Ben M’Sik, Hassan II University of Casablanca, BP 7955, Sidi Othman, Casablanca, Morocco 3Laboratory of Modeling, Analysis, Control and Statistics, Department of Mathematics and Computer Science,


INTRODUCTION
The field of epidemiology with the science of mathematics have been developed to study the transmission laws of epidemics [1]. Mathematical models provide the opportunity to understand how pandemics are spread and transmitted, taking into consideration the fact that these models present a mathematical translation of different hypotheses concerning the process of an epidemic transmission [2,3]. However, in order to define outbreaks of various types of epidemics and to provide insights into disease control and policy formulations, mathematical formulations have been developed [4,5]. Based on this data, effective control and preventive measures are suggested [6]. However, with the effect of spatio-temporal spread of epidemics, mathematical modeling should take into account the geographical criterion to show the spatial spread of an infectious disease within different geographical zones [7].
One of the basic models that was successfully investigated was the Kermack-Mckendrick model. To model an epidemic, the population being studied is divided into three classes labeled S, I, and R [8]: (S) refers to susceptible people who are not infected, but the possibility of transmitting the infection is still existed.
(I) Infected individuals who receive the infection, and able to spread it by contact with other people.
(R) Recovered or removals are individuals who become immune after getting sick, or individuals who are isolated from other members of the group, or ones who die due to the disease [9,10].
The transmission process of an epidemic is described when a population of susceptibles is being introduced into infectious individuals, then the infection is spread in the group through different modes of transmission [11,12].
Recently, the COVID-19 Virus has shown the necessity of taking into consideration the spatial dynamic of epidemics, and described how spatial heterogeneity affects the transmission dynamics of susceptible and infected populations [25,7]. The Corona virus was reported firstly in Wuhan, China, the outbreak was greatly increased and moved to other Chinese cities and multiple countries, moving to other continents [26,27,28,29]. In the history there were also many pandemics that show the spatio-temporal spread of pandemics such as the black plague [30], cholera [31], and others ... [32].
The discrete time multi-regional SIR model is a mathematical modeling of spatial and temporal spread of epidemics, an example of this model is made by [7], the multi-regional model is presented in multiple geographical areas to control the movements of the pandemic, and the infection can be spread from one region to another through travel. However, three main approaches are cited in [18], raise awareness by organizing vaccination campaigns, travel blocking movements coming from infected areas, and treatment. Other models are analyzed in this topic from many researchers in [33,34].
In the history of all these diseases, we can notice their spread from one region to another, and recently the COVID-19 pandemic from its epicenter of Wuhan in China has spread to all parts of the world, which makes taking into account the spatial spread of diseases more important during modeling processes.
The authors in [7] present the first work in the modeling and control of spatio-temporal spread of an epidemic using a multi-region SIR discrete-time model, as a generalization of the concept of classical models and aiming at a description of the evolution of pandemics, Zakary et al proposed a new approach of modeling of the spread of epidemics from one area to another using finite-dimensional models for the Spatio-temporal propagation of epidemics as an alternative of the partial derivatives models which are of infinite dimension. The authors also suggested some control strategies such as awareness-raising, vaccination, and travel-restriction approaches that could prevent specific infectious diseases such as HIV / AIDS, Ebola, or other epidemics in general [7,35,19,18,36], other researchers have shown the power and effectiveness of educational workshops and awareness programs in reducing the number of infected individuals [37,38,39].
In this paper, we propose a new optimal control approach mainly based on a multi-regions discrete-time system and a new form of multi-objective optimization criteria with importance indices and which is subject to multi-points boundary value optimal control problems. With more clarifications and essential details, we devise here a multi-regions discrete model for the study of the spread of an epidemic in M different regions, and analyze the effectiveness of vaccination (or awareness) optimal control strategies when vaccination (or awareness) campaigns are organized in infected zones. Here, we study the case when controls are applied to people who belong to all those regions and which are supposed to be reachable for every agent (nurse, doctor or media) who is responsible for the accomplishment of control strategies followed against the disease.
We consider an area as an infected zone if its number of infected individuals exceeds a threshold defined by the health decision-makers. Therefore, by varying the values of this threshold and then simulating the infection situation for different values of these thresholds shows that it is necessary to think about reducing the time between the first infection and the implementation of the control strategy. Unexpected results that in some situations the neighboring regions infected and its number of infections exceeds the threshold before the number of infections of the region source. This makes the implementation of the control strategies in the neighboring zones more important.
In our modeling approach we divided the studied area Ω into different zones that we call cells. A cell C j ∈ Ω can represent a city, a country or a larger domain. These cells are supposed to be connected by movements of their populations within the domain Ω. We define also a neighboring cells C k of the cell C j all zones connected with C j via every transport mean, thus a cell C j ∈ Ω can have more than one neighboring cell. Here, we suppose that a cell can be infected due to movements of infected people which enter only from its neighboring zones.
We carry out the map of the studied area and then we use different threshold values in the controlled multi-region SIR model to simulate the epidemic spread within the Casablanca-Settat The paper is organized as follows: Section 2. presents the discrete-time multi-region SIR epidemic system. In Section 3., we announce theorems of the existence and characterization of the sought optimal controls functions related to the optimal control approach we propose.
Finally, in section 4., we provide simulations of the numerical results applied to the Casablanca-Settat region and Rabat-Salé-Kénitra region as domain of interest.

MODEL DESCRIPTION AND DEFINITIONS
Based on same modeling assumptions of the reference [7], we assume that there are M geographical regions denoted C j (sub-domains) of studied Ω Where C j can represent a city, a country or a larger domain. We note by V (C j ), the vicinity set, composed by all neighboring cells of C j given by Where C j ∩ C k = / 0 means that there exists at least one mean of transport between C j and C k .
Note that this definition of V (C j ) is more general where it defines a more general form of vicinity regardless the geographical location of zones.
For example, in the Table 1 we can see that the studied area consists of 16 zones.
The multi-regional discrete-time SIR model associated to C j with ε where the disease transmission coefficient β C k j > 0 is the proportion of adequate contacts in domain C j between a susceptible from C j ( j = 1, ..., M) and an infective from another domain C k , d j is the birth and death rate and γ j is the recovery rate ans α C j is the proportion of mortality due to the disease. The biological background requires that all parameters be nonnegative. S

THE MODEL WITH VACCINATION
3.1. Presentation of the model with the control. In this section, we introduce a control variable u C j i that characterizes the effectiveness of the vaccination in the above mentioned model (1)(2)(3). This control in some situations can represent the effect of the awareness and media programs [18,19].
In almost all infectious diseases, the authorities determine the threshold of risk based on many factors, such as availability of medical equipment, budgets, and medical personnel ... Thus, they can wait some time to see the course of events before the intervention. If the number of casualties exceeds this limit, decision-makers have no choice but to start trying to control the situation. This motivate us to define a Boolean function ε associated to domain C j , that will be called the importance function of C j . Where ε C j i is either equaling to 1, in the case when the number of infected of the cell C j at instant i is greater than or equal to the threshold I C j defined by the authorities and health decision-makers, or ε C j i = 0 otherwise. Therefore, we define the importance function ε C j i by the Heaviside step function H as follows Then for a given domain C j ∈ Ω, the model is given by the following equations Our goal is obviously to try to minimize the population of the susceptible group and the cost of vaccination in all affected regions. Our control functions taking values between u C j min and u 3.2. An optimal control approach. We devise in this paper an optimal control approach for each region with different importance functions ε C j i , j = 1, ..., M. We characterize an optimal control that minimize the number of the infected people and maximize the ones in the removed category for all affected regions. Then, we are interested by minimizing the functional R > 0 are the weight constants of control, the infected and the removed in region C j respectively, and u = u C 1 , ...., Here, our goal is to minimize the number of infected people, minimize the systemic costs attempting to increase the number of removed people in each C j (with ε C j i = 1). In other words, we are seeking an optimal control u * such that where U is the control set defined by The sufficient condition for existence of an optimal control for the problem follows from theorem 1. At the same time, by using Pontryagin's Maximum Principle [40] we derive necessary conditions for our optimal control in theorem 2.
For this purpose, we define the Hamiltonian as

Theorem 2. (Necessary Conditions)
Given the optimal control u * and solutions S C j * , I C j * and R C j * , there exists ζ C j k,i , i = 1...N, k = 1, 2, 3, the adjoint variables satisfying the following equations R are the transversality conditions. In addition, Proof. Using Pontryagin's Maximum Principle [40], we obtain the following adjoint equations To obtain the optimality conditions we take the variation with respect to control u C pq i and set it equal to zero and ε Then, we obtain the optimal control  Fig.26, with the initial populations given in Table 1.
By the bounds in U (and U C j ) of the control, it is easy to obtain u C j * i in the following form

NUMERICAL RESULTS
In this section, we present numerical simulations associated to the above mentioned optimal control problem. We write a code in MAT LAB T M and simulated our results for several scenarios. The optimality systems is solved based on an iterative discrete scheme that converges following an appropriate test similar 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 co-state variables obtained at the previous steps. Finally, we execute the previous steps till a tolerance criterion is reached.

Area of interest.
We chose the Casablanca-Settat region and the Rabat-Salé-Kénitra region as the studied area Ω in this paper because we are convinced that we can find some useful data to support our work. They are the most populated and dynamic regions of Morocco, which contain the Rabat city as the capital of Morocco and the seventh largest city in the country with an urban population of around 580,000 inhabitants (2014) and a metropolitan population of more than 1.2 million inhabitants. It is also the capital of the administrative region of Rabat-Salé-Kénitra. They contain also the Casablanca city as the economic and industrial capital of Morocco because with its demographic growth and continuous development of the industrial sector, and the 14 other provinces (see Fig.1), in order to illustrate the objective of our work. zones, therefore, hereafter we note I C j as I min . We suppose as initial states in the area of interest Ω the following values: Susceptible: The real populations given in Table 1.
Infected: 100 infections only in the city of Casablanca, and 0 for the others.
Recovered: We assume that i = 0 represents the first appearance of the epidemic, therefore,  Table 2 for all zones.              C 6 , C 9 , C 11 tends to 0 from the time i = 120. Thus up to the regions C 2 , C 5 and C 10 whose number of susceptible decreases towards 0 at the moment i = 150. We also note that the number of susceptible individuals decreases over time less with the I min = 500 strategy than with I min j = 1000. Fig.17 and Fig.18 show the geographical evolution and graphs of infections in the different regions by applying the vaccination strategy from 500 infected. All regions recognize the same evolution of its infections with the difference in time which begin to grow and the time which registers its maximum value which is counted between 600 and 650 infected. The closest areas to the city of Casablanca C 8 , C 12 , C 13 , C 14 , C 15 and C 16 , begin first from i = 50 and arrives at the peak at time i = 100. Then the regions C 3 , C 4 , C 6 , C 7 , C 9 , C 11 , the least close, its infections rise from the moment i = 75 and reach its maximum value at the instant i = 120, then, remain almost constant until the end of vaccination. After the infected regions C 2 , C 5 and C 10 begins to rise from time i = 100 and reaches its maximum value at the moment i = 150 with 600 infected.   5 We also note that the number of infected is around 600 cases while with the vaccination strategy from 1000 infected the number exceeds 1200 cases. Fig.19 and Fig.20 show the evolution of the recovered in the different regions by applying the vaccination strategy from 500 infected. We find that the number of recovered starts with zero at the beginning then begins to grow but from different periods and exceeds the number 3.10 5 of the recovered, which is more important than the recovered when there are no control strategies, which are at best reached the value than 200 cases. The regions closest to Casablanca start to grow from i = 80 and very quickly reach extreme values which exceed 3.10 5 , then the and reaches the value of 4.5.10 5 . We also note that the number of recovered is the same for the two vaccination strategies after 500 and 1000 with the difference that with 500 the number of recovered increases with a shorter time than with the vaccination strategy from 1000 infected. 4.7. Scenario 3: Application of the Vaccination control from the beginning of the epidemic (I min = 0). In this scenario we assume that the epidemic is well known in other places, thus, we apply the control interventions from the declaration of such epidemic.

CONCLUSION
In this paper, we devised a novel optimization approach that represents an extension of the optimal control approach studied in the work of Zakary et al. in the paper [35]. We applied this new approach to a multi-region discrete epidemic model which has been firstly proposed in [7].
We suggested in this article, a new analysis of infection dynamics in M regions which we supposed to be accessible for health authorities. By defining new importance functions to identify affected zones and then will be treated. We investigated the effectiveness of optimal vaccination control approach, we introduced into the model, control functions associated with appropriate control strategies followed in the targeted regions by mass vaccination campaigns by considering different scenarios. Based on our numerical simulations, we showed the geographical spread of the epidemic and the influence of each region on another and then we deduced the effectiveness of each strategy followed. We concluded that the last scenario of optimal control approach when I min = 0 has given better results than the other cases regarding the maximization of the number of removed individuals and minimization of the spread of infection in all regions studied, but this is clearly the most expensive scenario. Thus, as a result, it is necessary to define small thresholds to control the situation as much as possible.

DATA AVAILABILITY
Data of the actual populations of the Casablanca-Settat region from [43] and for the Rabat-Salé-Kénitra region from [44].

CONFLICT OF INTERESTS
The author(s) declare that there is no conflict of interests.

FUNDING STATEMENT
The author(s) received no financial support for the research, authorship, and/or publication of this article