On the analysis of a multi-regions discrete SIR epidemic model: an optimal control approach.

In this paper, we devise a discrete time SIR model depicting the spread of infectious diseases in various geographical regions that are connected by any kind of anthropological movement, which suggests disease-affected people can propagate the disease from one region to another via travel. In fact, health policy-makers could manage the problem of the regional spread of an epidemic, by organizing many vaccination campaigns, or by suggesting other defensive strategies such as blocking movement of people coming from borders of regions at high-risk of infection and entering very controlled regions or with insignificant infection rate. Further, we introduce in the discrete SIR systems, two control variables which represent the effectiveness rates of vaccination and travel-blocking operation. 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. The multi-points boundary value problems, associated to the optimal control problems studied here, are obtained based on a discrete version of Pontryagin's maximum principle, and resolved numerically using a progressive-regressive discrete scheme that converges following an appropriate test related to the Forward-Backward Sweep Method on optimal control. © Springer-Verlag Berlin Heidelberg 2016.

Introduction

Mathematical modeling has played important role in describing and analyzing the evolution of infectious diseases. In 1927, Kermack and McKendrick [1] were the first researchers on mathematical epidemiology to propose the SIR model which is one of the classical epidemic models that has a compartment structure. In SIR systems, the host population is divided into three epidemiological groups: the susceptibles (S) (individuals not yet infected with the disease), the infectives (I) (individuals who have been infected with the disease and are capable of spreading the disease to those in the susceptible category), and the removed (R) (individuals who have been infected and then removed from the disease. Those in this category are not able to be infected again or transmit the infection to others). The transmission dynamics of an infectious disease is described by modeling the population movements among those epidemiological compartments. In these types of models, individuals develop an immunity to some diseases such as chicken pox and SARS. In addition, contacts between infected and susceptible individuals can take part in the spread of the epidemic.

On a large geographical scale, the disease becomes spatially mobile to different regions due to the movement of people from a region to another. The infection transport network becomes more complex when people come via airplane for instance, from further places.

The mathematical models conventional of the spread of the disease, have usually ignored or overlooked the spatial dynamics, while the spatial spread of infection has been observed many times [2]. In fact, there exist in the history of infectious diseases, many examples of infection that was observed spreading with respect to space. Sometimes, epidemics spread over large areas, and can even reach continents. Such cases include the Black Death (plague) which emerged in Europe in the 1300s, measles and smallpox in the New World between the 1500s and 1600s, HIV/AIDS in 1981, West Nile virus in North America in the late 1990s, and SARS in Asia in 2003 [3].

The spatial spread of infectious diseases is a phenomenon that involves many different components, which makes its modeling a complex task. Spatial heterogeneity requires a formulation of a small model for each region (city, or state, etc). One possible approach is to consider the travel of individuals between discrete geographical regions (subdomains of domain of study), assuming that the transmission does not take place during travel. For diseases that affect animals and insects (vectors), and can also be passed on to humans [4], it is important to accurately model not only the contact/mixing between humans and the various species, but also the movement patterns of each of the species involved, including humans.

The importance of modelling spatial spread of infectious disease can be highlighted by considering the Ebola outbreaks disease in several sub-Saharan African countries [5]. Before the diagnosis of that disease, infected people moved from one city to another. After some infected cases were diagnosed, and the problem was recognized, it took some time until the political decision was made to ban all movements. At that time, the infected cases had already been divided among a large number of regions. The problem was exacerbated. Then, in order to control the spatial spread of the disease, it has became necessary to consider all parts that people can visit. Hence, this study showcases the significance of integration of spatial dynamics in mathematical models of the spread of infectious diseases.

There is an increasing interest in the study and application of spatial spread [3, 6, 7]. Most of the models studied, have been partly continuous because of their mathematical tractability. In Allen et al. [8], a discrete-time, age-independent SIR type epidemic model with n subgroups is formulated and analyzed, then applied to a measles epidemic on a university campus. In similar pattern, we devise a discrete-time SIR model that describes the propagation of a disease in a population of individuals who travel between p regions (domains).

Figure 1 illustrates an example of discrete geographical domains of region of Casablanca-Settat (Morocco) where , that image was originally made based on information from [9]. Births and deaths are included but the population size of each domain remains constant. One reason for the upsurge of discrete epidemic models is that discrete models have advantages in describing an infectious disease since epidemic data are usually collected in discrete time units, which would make it more convenient to use discrete-time models [10] and the ease with which data can be compared to the simulated results [8].

Fig. 1

Region of Casablanca-Settat in Morocco. This region is divided into nine subregions (or domains) : Casablanca (), Mohammedia (), El jadida (), Nouaceur (), Mediouna (), Ben Slimane (), Berrechid (), Settat () and Sidi Bennour () [9]

In the literature, there have been many studies of discrete epidemic models [11-15] and references therein. In Reference Allen et al. [11], explained the advantage of SIS and SIR discrete-time models in approximating the more well-known continuous-time epidemic models. As regards to the positivity of solutions, discrete-time equations often provide positive solutions. However, differential equations give similar behaviors of solutions when time step is approaching zero [12].

There are different approaches to model the evolution of infectious diseases in discrete time. The recurrent difference equations from the discretisation of continuous differential equation models is one of the direct modeling approaches and frequently used, since this type of model can be well understood in the application under reasonable assumptions knowing that there are limits on the range of parameters [13, 14]. We will use this approach to formulate the model considered in this paper. Optimal control theory is well used as an available and effective option for diseases control, mainly Tuberculosis [16-20], Malaria [21, 22], HIV [23-31], Hepatitis [32, 33], Vector borne diseases [34], Cancer [35-40], and other diseases [41-47], but there are very few applications with both space and time as discrete variables. The aforementioned studies, especially those related to infectious diseases, have mainly focused on the optimization of the intervention for single region, rather than the intervention for multiple regions, hence, this work has a novel control application in addition to providing a framework to analyze spatial control strategies. The aim of this work is not to consider a special disease but to set up an optimal control problem relative to multi-regional discrete-time SIR model, applicable for any type of infectious diseases. Firstly, we represent the percentage of vaccinated susceptible populations as a control function of time in the SIR model corresponding to the targeted region aiming to control. Hence, the purpose of this optimal control (vaccination) strategy is to minimize the infected and susceptible individuals, and to maximize the total number of recovered individuals in that region by using the minimum possible cost of applying this control and simultaneously, to investigate the sensitivity of the susceptible individuals by the control (vaccination). We illustrate how the optimal control theory and the percentage of the vaccination (control) variables can be applied to minimize the susceptible and infected individuals and increase the removed ones.

The importance of taking into account the spatial spread of epidemic, is manifested in showing the influence of infection rates in regions at high-risk on theirs neighbors. Secondly, another control variable is added, characterizing the travel-blocking operation which attempts to block contacts between susceptibles of the controlled region (by vaccination) and infected individuals of domains at high-risk. We derive the optimality system for the multi-regional SIR model with the percentage of vaccinated individuals and travel-blocking operation. Then we find an optimal control strategy for this model and solve numerically this system by using an iterative procedure.

The rest of this paper is organized as follows: We describe the model without controls in Sect. 2. Objective functional and the analysis of optimal control with numerical simulations, for the case when only the vaccination control is used, are given in Sect. 3. Section 4 includes the analysis of optimal control problem and numerical simulations of the case when the travel-blocking control is also added. Finally, we conclude in Sect. 5.

The multi-regions discrete SIR model

We assume that there are p geographical regions (domains) of the domain studied (In Fig. 1, is the region of Casabanca-Settat with ).

Let , and let be the population of domain at time i, i.e., the number of individuals who are physically present in , both residents and travelers. According to the disease transmission mechanism, the host population of is grouped into three epidemiological compartments, let ,  and be the number of individuals in the susceptible, infective, and removed compartments of at time i, respectively. In addition to the death and recruitment, there are population movements among those three epidemiological compartments from time unit i to time . 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 . And death takes place anywhere with a per capita rate . 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 . The individuals in the infective compartment can keep being the infective, or get recovery and be transferred to the recovered compartment, or die, assuming that there is no mortality due to infection. The individuals in the recovered compartment never leave the compartment unless they die.

The disease transmission in a given domain at time i is modelled using standard incidence, given bywhere the disease transmission coefficient is the proportion of adequate contacts in domain between a susceptible from and an infective from another domain .

The multi-regions discrete-time SIR model associated to is written as follows where is the birth and death rate and is the recovery rate. The biological background requires that all parameters be non-negative.

is the population size corresponding to domain at time i. It is clear that the population size remains constant for all , in fact

The model with vaccination only

Presentation of the model

We introduce a control variable that characterizes the effectiveness of treatment (vaccination) in the above mentioned model (1–3). Then, for a given region 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 treatment, while increasing the population in the removed group. Our control function is assumed taking values between and , where and , .

An optimal control approach

We are interested in controlling the population of region . Then, the problem is to minimize the objective functional given bywhere are the weight constants of control, the infected and the removed group respectively, and . Our goal is to minimize the infected group, minimize the systemic costs attempting to increase the number of the removed individuals in . In other words, we are seeking an optimal control such thatwhere is the control set defined by

The sufficient condition for existence of an optimal control for the problem follows from Theorem 1 in [48], and at the same time, by using the Pontryagin’s Maximum Principle [49] we derive necessary conditions for our optimal control. For this purpose, we define the Hamiltonian as:

Theorem 3.2.1

(Sufficient conditions) For the optimal control problem given by (7) along with the state Eqs. (4)–(6), there exists a control such that

Proof

See Dabbs ([48], Theorem 1).

Theorem 3.2.2

(Necessary Conditions) Given an optimal control and solutions , there exist , the adjoint variables satisfying the following equations: where , are the transversality conditions. In addition

Using Pontryagin’s Maximum Principle [49], and setting , , and we obtain the following adjoint equations:with . To obtain the optimality conditions we take the variation with respect to control and set it equal to zeroThen we obtain the optimal controlBy the bounds in of the control, it is easy to obtain in the following form

Numerical results

We now present numerical simulations associated to the above mentioned optimal control problem. We write a code in and simulated our results using different data. 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 costate variables obtained at the previous steps. Finally, we execute the previous steps till a tolerance criterion is reached.

The multi-regions SIR epidemic model we suggest here, is applicable for any number of regions p. For an example, and in order to show the importance of our work, we choose i.e, we consider three regions (domains) , and with different parameters cited in Table 1. We try to control by the vaccination term given by (12).

Table 1

Parameters values of and utilized for the resolution of all multi-regions discrete systems (1)–(3), (4)–(6) and (15)–(17), and then leading to simulations obtained from Figs. 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 and 12, with the initial conditions and associated to the three regions and

	S_0	I_0	R_0	β	d	γ
Ω_1	10,000	900	0	0.051	0.16	0.0025
Ω_2	8000	900	0	0.2481	0.1346	0.003
Ω_3	7000	50	0	0.11	0.219	0.025

Fig. 2

Shapes of state functions of the discrete system (1)–(3) where there is yet no control. a State variables S, I and R associated to the region . b State variables S, I and R associated to the region . c State variables S, I and R associated to the region

Fig. 3

Shapes of state functions of the discrete system (4)–(6) when the control is now introduced. given by (12). a Optimal state variables and associated to the region . b Optimal state variables and associated to the region . c Optimal state variables and associated to the region

Fig. 4

Shape of the optimal control function showing the values taken by the effectiveness rate of the vaccination in the first region during one month

Fig. 5

Impact of the control severity weight on the number of the removed people of (a), and also on the values of the control (b)

Fig. 6

Impact of different values of the infection rate on the shape of the state functions S(t), I(t) and R(t) in region

Fig. 7

Simulations of the optimal state functions and associated to the multi-regions discrete system (14)–(16) with controls and given by (24) and (25) respectively. a Dynamics of the population in region . b Dynamics of the population in region . c Dynamics of the population in region

Fig. 8

The optimal effectiveness of the vaccination control strategy: (see a), the optimal effectiveness of the travel-blocking control strategy: (see b)

Fig. 9

Comparison of control with and without use of control

Fig. 10

Comparison of control with and without use of control with different choices of . a in the presence of control . b in the absence of control

Fig. 11

Impact of the control severity weight , with , on the number of the infected people in (simulation a), and also on the control (simulation b)

Fig. 12

Impact of the control severity weight , with , on the number of the infected people in (simulation a), and also on the control (simulation b)

In Fig. 2, we can observe that in the absence of a control and in the presence of an epidemic that spreads in three regions characterized by different parameters, the number of infected individuals rise from , and as initial conditions to 5730, 5430 and 3200 in the three regions respectively. Once a vaccination control is introduced in the system (4)–(6), particularly in the equations that describes dynamics of , and functions associated to the first region, we can deduce its effect on decreasing the number of infected people in Fig. 3, from 5730 when there was yet no vaccination strategy, to 2750 when there is the control . One of the major benefits of that control, is to increase the number of the removed people, and this can be observed in the case (a) of Fig. 3, where the number of the removed people becomes approximately equal to 5500, and that can obviously prove the effectiveness of the vaccination strategy in the first region with a rate that varies in Fig. 4, from a value equal to 0.15 towards a value equal to 0.3, and which also proves that by a control taking only a nonzero value close to 0, we can reach our goal with a significant number of the removed people.

We observe in Fig. 5 that the more the control severity weight is small, the control is more important and then there is an increasing of the number of removed people in (a), because has obviously an impact on the values of the control (see simulation (b)) from the control characterization (12), and also these results are obtained from the fact that the most important role of the vaccination control is to increase the number of the removed people .

Figure 6 is added here, to show the importance of another control associated to a travel-blocking operation, instead of following only a vaccination strategy. In fact, the number of infected people in is increased by hosting people coming from , and we can see the impact of the infection rate on the system (1)–(3) of in Fig. 6, and we deduce that the more is big, the more the number of infected people in the first region becomes huge.

The model with vaccination and travel-blocking

Let , and denote by , the set of indexes of domains at high-risk. We introduce a control variable that characterizes the travel-blocking operation, in order to restrict movements from domains to . WhereThen for a given region the model is given as follows:

The problem is to minimize the objective functional given bywhere are the weight constants of controls, the infected and the recovered respectively, , and .

Our goal is to minimize the infected group, minimize the cost of applying controls and increase the number of recovered in .

In other words, we are seeking optimal controls and such thatwhere and are the control sets defined by where and , .

By using Pontryagin’s Maximum Principle, [49] we derive necessary conditions for our optimal controls. For this purpose, we define the Hamiltonian as:

Theorem 4.2.1

Given optimal controls , and solutions , there exist , the adjoint variables satisfying the following equations: where are the transversality conditions.

In addition, we have

Using Pontryagin’s Maximum Principle [49] and (13) we obtain the following adjoint equations:with . To obtain the optimality conditions we take the variation with respect to controls ( and ) and set it equal to zeroThen we obtain the optimal control pairBy the bounds in and of the controls, it is easy to obtain and in the following form

Remark 4.2.1

If h is the cardinality of then there are controls.

In this subsection, we keep the same data used above (Table 1) to simplify the comparison of results obtained in both cases. , and we are also interested in controlling , by assuming that is at high-risk of infection (see Fig. 2), then . We introduce the second control given by (25), in order to reduce entry of infected people from region .

Because of the impact of the infection rate associated to only one region on neighbor regions, by increasing the number of its infected people via travel as it was shown in Fig. 6, and in order to show the advantage of the travel-blocking control , in decreasing the number of infected people who travel to the first region from the second region, we deduce from Fig. 7 that the number of infected people in the first region decreases from 5730 when there was no controls yet in Fig. 2, and from 2750 when the control is introduced alone in Fig. 3, towards a smaller number equal to 575 when the control is added, and which can prove that the travel-blocking operation was successful to prevent the disease from spread. We should note that the more is small, the more we can obtain good results regarding the evolution of infection rate. In fact, it is easy to deduce from the discrete control system (14)–(16), that when is closer to 0, the term “” converges to 0 and when is far from 0, the term “” participates to an increase of the function. In contrast, the more is bigger and far from 0, the more we can obtain satisfactory results regarding the evolution of the number of the removed people.

Figure 8 presents the simulations of and , and we can see that varies from 0.09 ( in Fig. 4) to 0.31 ( in Fig. 4) and varies from 0.0275 to 0.05 which proves that with this small close to 0, the decrease of the number of infected people in Fig. 7 is more important than the case where there was no travel-blocking yet. In addition, the effect of the lacking amount of the function in Fig. 8 compared to the simulation in Fig. 4, could be replaced by the effect of control , and by introducing both controls and , we would obtain better results that prove the importance of the travel-blocking control. In general, it does not occur a very big change of the control function from the case when there is no control term yet to the case when is introduced to (14)–(16) as we can see in Fig. 9, but with changing the values of the severity weight associated to the control in cases (a) and (b) of Fig. 10, we can observe more the impact of on , and the lacking effect of could be seen replaced by the effect of , because for different values of , the control becomes smaller in the presence of in case (a), and bigger in case (b) when is not introduced yet to (4)–(6).

Figure 11 also shows the impact of on the values of , and we deduce that more is important, the more is also important. In fact, by analogous conclusions from Fig. 5, we deduce that in the case of a travel-blocking control that the more the control severity weight is big, the more the number of infected people in (a) is big, because has obviously an impact on the values of the control (see simulation Fig. 11b) as we explained above (the more is small the more is big), and also these results are obtained from the fact that the most important role of the travel-blocking control is to decrease the number of the infected people without forgetting that the more is closer to zero, the more the number infected becomes insignificant. Opposite results are obtained in Fig. 12 by changing the values of the control severity weight , and we obtain small (see simulation Fig. 12b) whenever is taken small values because has obviously an impact on the values of the control from the control characterization (25).

For both Figs. 10 and 11, we can explain those numerical results by the fact that more is big, the more is small (could obviously be deduced by the location of the severity weight in the denominator of the characterization formula (12) of in Theorem 3.2.2, and then the effect of becomes important and valuable. We can obtain analogous results for when we change the values of the severity weight associated to the control , as we see in Fig. 12.

Conclusion

We analyzed in this work two different optimal control problems for a discrete-time multi-regions SIR epidemic model. In the first part, only one control that characterizes the effectiveness of vaccination is used, and where the optimal control problem was subject of an optimization criterion represented by the minimization of an objective function aiming to minimize the number of infected people and the cost of vaccination while maximizing the number of removed people in the targeted region aiming to control. Numerical results, associated to this first part, showed the effectiveness of this vaccination while the targeted region was still sensitive to infection rates in its neighboring domains at high-risk. Then, we extended the first part to a second part by adding a travel-blocking control, in which we have controlled contacts between susceptible people of the controlled region (by vaccination) and the infected of domains at high-risk of infection. We showed the advantage of each case studied with a comparison between the numerical simulations they provided. We have identified optimal control strategies for several values of the controls severity weights to show the importance and the effectiveness of our approach in controlling the infection spread. Control programs that follow those strategies can effectively reduce the number of infected cases and increase the number of the removed individuals in a specific region.