A multi-region discrete time mathematical modeling of the dynamics of Covid-19 virus propagation using optimal control.

We study in this work a discrete mathematical model that describes the dynamics of transmission of the Corona virus between humans on the one hand and animals on the other hand in a region or in different regions. Also, we propose an optimal strategy to implement the optimal campaigns through the use of awareness campaigns in region j that aims at protecting individuals from being infected by the virus, security campaigns and health measures to prevent the movement of individuals from one region to another, encouraging the individuals to join quarantine centers and the disposal of infected animals. The aim is to maximize the number of individuals subjected to quarantine and trying to reduce the number of the infected individuals and the infected animals. Pontryagin's maximum principle in discrete time is used to characterize the optimal controls and the optimality system is solved by an iterative method. The numerical simulation is carried out using Matlab. The Incremental Cost-Effectiveness Ratio was calculated to investigate the cost-effectiveness of all possible combinations of the four control measures. Using cost-effectiveness analysis, we show that control of protecting susceptible individuals, preventing their contact with the infected individuals and encouraging the exposed individuals to join quarantine centers provides the most cost-effective strategy to control the disease. © Korean Society for Informatics and Computational Applied Mathematics 2020.

Introduction

The new Corona Virus (Covid-19) appeared in Wuhan, China on September 31, 2019, where the World Health Organization has received multiple cases of pneumonia in this city. The characteristics of the virus did not match any known viruses. This new virus comes from a group of viruses that include those that cause common cold, SARS syndrome, and respiratory syndrome in the Middle East. This virus was temporarily called the emerging corona virus (Covid-19) [1].

This virus took several forms and names, among them MERS-COV. It first appeared in 2012 in the Kingdom of Saudi Arabia and then in 2015 in South Korea and the number of people infected with this disease was high. These countries used the centers of treatment to control the spread of this disease along with other preventive measures. Some studies conducted in 2015 showed that 1046 patients were infected with the virus MERS-COV and cases of deaths reached 298 among which 67 cases were in the United Arab Emirates. In 2018, more than 500 camels died due to this disease in Saudi Arabia, but the number of deaths due to the new CORONA virus in humans has increased by up to 50 according to initial reports in mid-December 2019 in China [2].

The registered cases continued to increase rapidly in early 2020, with a total of 4593 confirmed cases of the new virus and 106 deaths according to the WHO World Health Organization report published on 28 January 2020. Chinese authorities announced that the first cases of the new virus came from the Chinese city of Wuhan, which is home to more than 10 million people, and has since spread to some Chinese provinces and cities, and the number of infections outside China has reached 56 infections in 14 other countries including: Japan, Korea, Thailand, Vietnam, Australia, Singapore, France, Germany, Nepal, Malaysia, USA, Canada, Cambodia and Sri Lanka. Corona virus symptoms include respiratory symptoms, fever, cough, shortness of breath, and difficulty of breathing. In most severe cases, the infection may cause pneumonia, severe acute respiratory syndrome, kidney failure and even death [3].

There are several mathematical modeling studies that have emerged which aim at understanding the corona virus and describing its dynamics. For example, [4] formulated a mathematical model for MERS transmission dynamics and estimating transmission rates. They estimated the basic reproduction number using the estimates of the transmission rates in the first two periods. Drosten et al [5] provided a description of a fatal case of MERS-CoV infection and associated phylogenetic analyses. Guery et al [6] analyzed the clinical features of infected cases. Tahir. M et al [7] used a non-linear mathematical model to study the dynamics of the transmission of MERS-CoV in human population through an agent known as camel. They divided their population into two classes: humans including (susceptible humans S(t) , infected humans I(t), recovered individuals R(t)), and camels including (Susceptible camels and Infected camels . They have found the basic reproduction number and discussed the stability analysis of their model using the basic reproduction number. They found that the model is locally asymptotically stable at disease free equilibrium when and when endemic equilibrium exists and becomes stable. In their research, equilibrium points are also globally asymptotically stable under certain conditions. Besides these, many mathematical models were formulated to investigate the dynamics of infectious diseases with control (For example: [8-14]...).

Additionally, most of the previous research has focused on continuous modeling. In this research, we will adopt discrete-time modeling where statistical data are collected at a discrete-time (day, week, month, and year). Treatment and vaccination of some patients are also given at a discrete-time. Therefore, it is more direct, more convenient, and more accurate to describe a phenomena using discrete time modeling compared to continuous time modeling. The use of discrete-time models may avoid some mathematical complications such as the choice of the area of function and uniformity of the solution. Thus, difference equations appear as a more natural way to describe epidemiological models. Moreover, numerical solutions are used for discretionary differential equations and this encourages us to employ the difference equations directly. The numerical exploration of discrete time models is easy, and therefore can be easily implemented by mathematicians [7, 15–18].

Besides the aforementioned works, we will study the dynamics of a mathematical model of the novel Corona virus which contains the following additions:

A discrete-time mathematical modeling,

Dividing space into multi regions,

A compartment that represents the number of the exposed individuals in region j who are carriers of the disease and have no symptoms,

A compartment that represents the number of the quarantined individuals who receive treatment in medical centers in region j.

The transmission of the disease from animals to humans,

The death rate induced by the infected humans in region j (),

The death rate induced by the individuals exposed to quarantine in region j( ),

The death rate induced by the infected animals in region j (),

The Corona virus propagation phenomenon differs from one region to another (city, village, neighborhood...). This difference is apparent between villages and cities. It basically takes place due to the educational level, human resources, public services, the interaction between humans and animals or the interaction that occurs among humans or animals according to each region. For example, most inhabitants of the city have means of transportation, roads, hospitals equipped with medical tools and pharmacies unlike rural regions which lack many of these public services. This topic of multi region or multi groups models has been widely studied by many researchers for example [19-22].

Our model consists of two groups, the human population and the animal population in the region j. The human population is divided into five different classes, namely: Susceptible humans Exposed humans Infected humans Quarantined humans and Recovered humans The animal population is divided into two different classes, namely: Susceptible animals and Infected animals Throughout this research, we seek to find the optimal strategies to minimize the number of the exposed humans, the infected humans and the infected animals and also to maximize the number of the exposed humans who will join quarantine centers. In order to achieve this purpose, we use optimal control strategies associated with four types of controls: the first represents awareness campaigns in region j that aim at introducing health measures to protect individuals from being infected by the virus. The second represents security campaigns and health measures to prevent the movement of individuals from one region to another. The third represents the effort to encourage the exposed humans to join quarantine centers. The forth control represents the disposal of the infected animals.

The paper is organized as follows. In Sect. 2, we present Corona disease, its spread manner and the proposed preventions protocols. In Sect. 3, we present our multi-regions discrete mathematical model that describes the dynamics of the Corona disease propagation. In Sect. 4, we present the optimal control problem for the proposed model and we characterize these optimal controls using Pontryagin’s maximum principle in discrete time. Numerical Simulations and Cost-effectiveness Analysis are given in Sect. 5. Finally, we conclude the paper in Sect. 6.

Corona disease, methods of its spread and the proposed safety protocols

The spread within the same region

Corona virus is one of the most harmful and pathogenic viruses that humans and animals have ever known. It specially attacks the human respiratory system. Corona virus is an infectious virus, which can be transmitted from an animal to a human and also from a human to another. The virus can be transmitted from an infected person to other people through the air, due to sneezes and coughs, and through direct contact with an infected individual via touching or greeting by hand. Contaminated places with the disease can also transmit the virus when touched.

However, WHO is trying to stop this disease by insisting on some recommendations. These recommendations aime at reducing exposure and transmission of the disease, including personal hygiene, respiratory hygiene and safe nutritional practices. The proposed safety measures contain the following:

Washing hands with soap and water, or rubbing hands with an alcohol cleanser, covering the mouth and nose with a medical mask, or bending the elbow when coughing or sneezing,

Avoiding contact with anyone with symptoms of a cold or influenza and seeking medical attention in the event of a fever, cough and difficulty breathing,

When visiting open markets, direct unprotected contact with live animals and surfaces that come into contact with animals should be avoided,

Cooking well, especially meat,

Preventing and reducing contact with susceptible individuals, infected individuals or infected animals within the same region,

Preventing the movement of susceptible individuals, infected individuals, and infected animals from one region to another.

Encouraging the affected individuals to join medical centers or subject themselves to quarantine inside each region,

Creating hospitals and quarantine areas [1].

The spread from one region to another

Corona virus first appeared in the Chinese city of Wuhan, then it moved to some neighboring cities, like: Huanggang, Yingcheng, Xiantao, Tianmen, Wuxue, Suizhou, Xinyang, Jiujiang, Yueyang due to of the movement of travelers, workers students and animals from one region to another. It then moved outside China to other countries like the United States of America, Japan, South Korea, and Thailand (see Fig. 1) [23].

Fig. 1

Map of Covid-19 virus propagation in the world

The multi-regions discrete mathematical model

Our model is divided into two mathematical systems, the system of human population and the system of animals population taking into consideration the specificity of each region and the movements of the individuals between different regions. We propose the global area of interest represented by the geographic union of the regions which are the fields of this study. These regions are supposed to be connected by the movements of their populations. The area can be divided into sub-domains of j() regions which represent a country, a city, a town, or a small domain. Let and p is the total number of the region.

The and are the numbers of the infectious classes associated with the domain j (). We note that the transmission of the Corona disease takes place through the contact between human–human or human–animal. The unit k can correspond to days, months or years.

For each region j , we write the number of the human population: susceptible humans , exposed humans , infected humans , quarantined humans and recovered humans as functions of discrete time k. The total human population is denoted by with . Similarly, we write the number of animals in each region j: susceptible animals and infected animals , respectively as functions of discrete-time k. The total animal class is denoted by with .

The following system describes the propagation of the disease in a given population for different regions:where and are the given initial states.

The susceptible human population represents the susceptible individuals who are infected including all population groups (children, young, adults, and elderly) in region j. The susceptible human population is increased by the recruitment of individuals at a rate in region j and decreased by an effective contact with the exposed humans at a rate ( represents the transmission rate from susceptible individuals to exposed individuals) and ( represents the transmission rate from susceptible individuals to infected individuals) and due to infection from the infected animals. Finally, the susceptible human population suffers natural mortality (at a rate ) in the same region and in a different region.

The exposed humans are the individuals carrying the disease and have no symptoms in region j. This compartment is increased at the rates and . The number of some exposed individuals decreases when the exposed individuals become infected at a rate represents the transmission rate from exposed individuals to infected individuals). It is decreased by natural death (at a rate ).

The infected humans are the individuals infected by the Corona disease. The infected population is increased by the exposed individuals at a rate and decreased by infected individuals who are subjected to quarantine at a rate . It is decreased by natural death (at the rate ) and is the death rate of infected human population.

The quarantined humans are the individuals who are subjected to quarantine in region . The quarantined population is increased by infected individuals (at a rate ). It is decreased by natural death (at the rate ) and the death rate of the infected humans.

The recovered humans represent the recovered individuals in region j . It is also increased when the quarantined individuals become recovered at a rate . The number of the recovered individuals of this population decreases by natural death (at the rate ).

The susceptible animals represents the susceptible infectious animals (chicken, serpent, bat, cat,...) in region j . This compartment of susceptible animals is increased by the recruitment of animals into the compartment at a rate in region j and decreased by an effective contact with the infected animals at a rate (which represents the transmission rate from susceptible animals to infected animals) and (which represents the transmission rate from susceptible animals in region j to infected animals of region r ). Finally, the susceptible animals suffer natural mortality (at a rate ) in region j.

The infected animals refer to the infected animals in region j. The infected animal’s compartment is increased by the rate in region j and . Finally, the infected animals suffer natural mortality (at a rate ) and infectious death rate of infected animals in region j.

The variables , , , , and are the numbers of the individuals in the seven classes at time k respectively. The unit k can correspond to periods, phases or years. It depends on the frequency of the survey studies as needed. The graphical representation of the proposed model is shown in Fig. 2.

Fig. 2

Schematic diagram of the seven infectious classes in the model

The total population size at time k in region j is denoted by with

The optimal control problem

At the present time, there is no treatment and no vaccine for the novel Corona disease 2019. As a result, many countries adopted some strategies such as: reducing contact with the susceptible individuals, infected individuals or infected animals within the same region, imposing tough measures to restrict the movement of individuals and animals from one region to another, putting the affected and susceptible individuals in quarantine inside each region and creating hospitals and quarantined areas.

The strategy of control that we proposed aims at minimizing the number of exposed humans infected humans , infected animals and maximize the number of the quarantined individuals during the time step to T in the domain j and also to minimize the cost spent in on awareness campaigns and an on security campaigns.

In the model (1), there are four controls , , . The first control represents the effort of awareness campaigns in regionj that aim at introducing health measures to protect individuals from being infected by the virus (cover all cuts, shoes and long sleeve shirts when handling animals, washing hands thoroughly on a regular basis, taking shower and cleaning up work place and home...). Therefore, the term is used to reduce the force of infections. The second control represents the security campaigns and health measures to prevent the movement of individuals from region j to region r (where, and p is the total number of the regions) and the term is used to reduce the force of infections. The control represents the effort to encourage the exposed individuals to join quarantine centers in the same region. Finally, measures the effort to dispose of the infected animals. So, the controlled mathematical system is given by the following system of difference equations.where and are the given initial states.

Then, the problem is to minimize the objective functionalWhere the parameters and for are the cost coefficients at time k in region j. T is the final time.

In other words, we seek the optimal controls , , and such thatWhere is the set of admissible controls defined byThe sufficient condition for the existence of an optimal control for problem (2-3) comes from the following theorem.

Theorem 1

There exists the optimal control and such thatsubject to the control system (3) with initial conditions.

Proof

See Dabbs ([32], Theorem 1).“Dabbs. K (2010) Optimal control in discrete pest control models. University of Tennessee Honors Thesis Projects”

We apply the discrete version of Pontryagin’s Maximum Principle [17, 24–28]. The key idea is introducing the adjoint function to attach the system of difference equations to the objective functional resulting in the formation of a function called the Hamiltonian. This principle converts the problem of finding the control to optimize the objective functional subject to the state difference equation with initial condition to find the control to optimize the Hamiltonian pointwise (with respect to the control).

Now, we have the Hamiltonian at time step k, defined bywhere is the right side of the system of difference equations (3) concerning the region j of the i th state variable at time step ..

Theorem 2

Given an optimal control , and solutions , , , , and of corresponding state system (3), there exist adjoint functions, , , , , ,and satisfyingwith the transversality conditions at time TFurthermore, for , the optimal controls ,and are given by

The Hamiltonian at time step k in region j is given byFor , the adjoint equations and transversality conditions can be obtained by using Pontryagin’s Maximum Principle, in discrete time, given in [10, 11, 18, 23, 27] such thatFor , the optimal controls , and can be solved from the optimality condition,that isWe haveHowever, the control attached to this case will be eliminated and removed. By the bounds in of the controls, it is easy to obtain , and in the form of (11-12-13-14).

Numerical Simulations and Cost-effectiveness Analysis

Numerical Simulations

In this formulation, there were initial conditions for the state variables and terminal conditions for the adjoints i.e. the optimality system is a two-point boundary value problem with separated boundary conditions at times step and . We solve the optimality system by an iterative method with forward solving of the state system followed by backward solving of the adjoint system. We start with an initial guess for the controls at the first iteration and then before the next iteration we update the controls by using the characterization. We continue until convergence of successive iterates is achieved.

We present some numerical simulations in order to illustrate our theoretical results. We consider system (3) with the following parameter values in region1 and region 2 (see Tables 1 and 2).

Table 1

The description of the parameters used for the definition of discrete time systems of region (1)

Sh,01	Eh,01	Ih,01	Qh,01	Rh,01	Sa,01	Ia,01	Nh1	bh1
1000	400	200	200	40	1000	600	1440	100
μh1	α1,h1	α2,h1	α3,h1	β1,h1	β2,h1	β3,h1	β12,h1	β22,h1
0, 065	0, 05	0, 001	0, 001	0, 4	0, 5	0, 6	0, 3	0, 4
δ1,h1	δ2,h1	δ3,a1	μa1	ba1	β4,a1	β52,a1	Na1
0, 0001	0, 0001	0, 0001	0, 025	40	0, 15	0, 4	1600

Table 2

The description of the parameters used for the definition of discrete time systems of regions (2)

Sh,02	Eh,02	Ih,02	Qh,02	Rh,02	Sa,02	Ia,02	Nh2	bh2
1200	500	300	300	60	1200	500	2360	100
μh2	α1,h2	α2,h2	α3,h2	β1,h2	β2,h2	β3,h2	β11,h2	β21,h2
0, 065	0, 04	0, 002	0, 002	0, 5	0, 4	0, 65	0, 3	0, 4
δ1,h2	δ2,h2	δ3,a2	μa2	ba2	β4,a2	β51,a2	Na2
0, 0002	0, 0002	0, 0002	0, 035	50	0, 2	0, 4	1700

This model is valid for all regions of the world, but the simulation will be limited to the study of a model that is composed of two regions: and represent neighboring regions of the region1.

We begin by presenting the evolution solution of our model (1) in region1 with contact and without controls with region 2 as shown in Figs. 3, 4 and 5.

Fig. 3

Evolution of the exposed humans as a function of time in region 1

Fig. 4

Evolution of the infected humans as a function of time in region 1

Fig. 5

Evolution of the infected animals as a function of time in region 1

Figures 3, 4 and 5 shows that the number of the exposed individuals, infected individuals and infected animals in region 1 has increased due to their contact with individuals or animals from region 2 compared to the absence of contact in region 1.

The aim is to highlight the specificities of each region when formulating control strategies which seek to adopt the following strategies:

Strategy 1: protecting susceptible individuals from contacting the infected individuals in the same region 1

To realize this strategy, we apply only the control i.e. the implementation a healthy protocol programs and awareness campaigns in region 1 to protect susceptible individuals from contacting infected individuals and infected animals in region1. (see the Figs. 6–7).

Fig. 6

Evolution of the exposed humans as a function of time in region 1

Fig. 7

Evolution of the infected humans as a function of time in region 1

Figure 6 shows that the number of the exposed individuals in region1 decreased from 868.52 (without controls) to 804.34 (with controls) at the end of the implementation of the proposed strategy. Figure 7 demonstrates that the number of the infected individuals in region 1 decreases from 657.01 (without controls) to 608.64 (with controls) at the end of the implementation of the proposed strategy.

Strategy 2: protecting and preventing susceptible individuals from contacting the infected individuals in the same region or in other regions.

To realize this strategy, we apply only the controls and i.e. the implementation of awareness campaigns in region 1 that aim at protecting individuals from being infected by the virus, health measures and security campaigns to prevent the movement of individuals from region 1 to another region. (see the figures below).

Figure 8 shows that the number of the exposed individuals in region 1 decreases from 868.52 (without controls) to 624.03 (with controls) at the end of the implementation of the proposed strategy. Figure 9 demonstrates that the number of the infected humans in region 1 decreases from 657.01 (without controls) to 472.32 (with controls) at the end of the implementation of the proposed strategy.

Fig. 8

Evolution of the exposed humans as a function of time in region 1

Fig. 9

Evolution of the infected humans as a function of time in region 1

Strategy 3: protecting susceptible individuals, preventing their contact with the infected individuals and encouraging the exposed individuals to join quarantine centers.

To realize this strategy, we apply only the control , and i.e. the implementation of awareness campaigns in region 1 that aim at protecting individuals from being infected by the virus, in addition to health measures and security campaigns inorder to prevent the movement of individuals from region 1 to another region. Also, we encourage the exposed individuals to join quarantine centers. (see the figures below).

Figure 10 shows that the number of the exposed individuals in region1 decreases from 868.52 (without controls) to 482.05 (with controls) at the end of the implementation of the proposed strategy. Figure 11 demonstrates that the number of the infected individuals in region 1 decreases from 657.01 (without controls) to 364.95 (with controls) at the end of the implementation of the proposed strategy. Also, the number of the quarantined individuals increases significantly from 10.15 (without controls) to 224.57 (with controls). (see Fig. 12)). Through the use of those controls, the above mentioned strategy has been proved effective.

Fig. 10

Evolution of the exposed humans as a function of time in region 1

Fig. 11

Evolution of the infected humans as a function of time in region 1

Fig. 12

Evolution of the quarantined humans as a function of time in region 1

Strategy 4: protecting susceptible individuals, preventing their contact with the infected individuals, encouraging the exposed individuals to join quarantine centers and the disposal of the infected animals. .

To realize this strategy, we apply only the controls , , and i.e. the implementation of awareness campaigns in region 1 which aim at protecting individuals from being infected by the virus, in addition to health measures, security campaigns to prevent the movement of individuals from region 1 to another region, encouraging the exposed individuals to join the quarantine centers and the disposal of the infected animals. (see the figures below).

Figure 13 shows that the number of the exposed individuals in region1 decreases from 868.52 (without controls) to 432.47 (with controls) at the end of the implementation of the proposed strategy. Figure 14 demonstrates that the number of the infected individuals in region1 decreases from 657.01 (without controls) to 328.69 (with controls) at the end of the implementation of the proposed strategy. Also, the number of the quarantined individuals increases significantly from 10.15 (without controls) to 202.31 (with controls) (see Fig. 15). Finally, the number of the infected animals decreases significantly from 1511 (without controls) to 986.72 (with controls) (see Fig. 16). Through the use of those controls, the above mentioned strategy has been proved effective.

Fig. 13

Evolution of the exposed humans as a function of time in region 1

Fig. 14

Evolution of the infected humans as a function of time in region 1

Fig. 15

Evolution of the quarantined humans as a function of time in region 1

Fig. 16

Evolution of the infected animals as a function of time in region 1

Cost-effectiveness analysis

In this section, we analyse the cost-effectiveness of the previous four senarios and strategies by comparing between these four control strategies to determine the most cost-effective strategy. Following the method as applied in several studies [29-31], we evaluate the costs using the incremental cost-effectiveness ratio (ICER). This ratio used to compare the differences between the costs and health outcomes of two competing intervention strategies.

The ICER is defined as the quotient of the difference in costs in strategies i and j, by the difference in infected averted in strategies i and j (i, j .

Given two competing strategies E and F, where strategy F has higher effectiveness than strategy E (), the ICER values are calculated as follow:Where the total costs (TC) and the total cases averted (TA) are defined, during a given period for strategy i for by:where and corresponds to the person unit cost of the four possible interventions, while (, )is the optimal solution associated to the optimal control (, , ). Using the simulation results and we ranked, in the Table 3 our control strategies in order of increased numbers of averted infections.

Table 3

Total costs and total averted infections for strategies 1–4,

Strategy	Total averted infections (TA)	Total cost (TC)
4	-10,632.10^4	11,331.10^4
3	-6,3735.10^4	0,32775.10^4
2	-4,1167.10^4	9,1125.10^4
1	-1,016.10^4	1,1253.10^4

Strategy 4 is compared with strategy 3 with respect to increased effectiveness, in reference to Table 3. So:Since , then strategy 4 is less effective than strategy 3. Therefore, strategy 4 is excluded from the set of alternatives.

Next, strategy 3 is compared with strategy 2. The ICER values for strategy 3 and strategy 2 are calculated below:Since , then strategy 2 is less effective than strategy 3. Therefore, strategy 2 is excluded from the set of alternatives.

Finally, strategy 3 is compared with strategy 1. The ICER values for strategy 3 and strategy 1 are calculated below:Since , then strategy 1 is less effective than strategy 3. Therefore, strategy 1 is excluded from the set of alternatives.

Thus, the conclusion is that strategy 3 (protecting susceptible individuals, preventing their contact with the infected individuals and encouraging the exposed individuals to join quarantine centers) is the most effective strategies.

Conclusion

In this work, we formulated a multi-region discrete mathematical model that describes the dynamics of transmission of the novel Corona virus between humans on the one hand and animals on the other hand in a region or in different regions. The strategies that we introduced include four controls: the first represents awareness programs, the second represents health measures and security campaigns, the third represents the effort to encourage the exposed individuals and the infected individuals to join quarantine centers and the forth control represents the disposal of the infected animals. Pontryagin’s maximum principle in discrete time was used to characterize the optimal controls and the optimality system was solved by an iterative method. The numerical simulation was carried out using Matlab. Numerically, the optimal control strategies that include protectiing individuals from being infected by the virus, security campaigns and health measures to prevent the movement of individuals from one region to another, encouraging the individuals to join quarantine centers and the disposal of infected animals (Strategies 1, 2, 3 and 4) showed significant differences in the number of the exposed humans, infected humans, quarantined humans and infected animals compared to the case where controls are not used. Using cost-effectiveness analysis, we showed that strategy 3 (protecting susceptible individuals, preventing their contact with the infected individuals and encouraging the exposed individuals to join quarantine centers) is the optimal and the most effective strategies.