key: cord-0771154-wlncrdv4 authors: Khajji, Bouchaib; Kada, Driss; Balatif, Omar; Rachik, Mostafa title: A multi-region discrete time mathematical modeling of the dynamics of Covid-19 virus propagation using optimal control date: 2020-05-08 journal: J Appl Math Comput DOI: 10.1007/s12190-020-01354-3 sha: 9f1421f795084d05cda18dcd08dc9bec99fac178 doc_id: 771154 cord_uid: wlncrdv4 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. 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 Extended author information available on the last page of the article 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 [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 nonlinear 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 (S c ) and Infected camels (I c )). They have found the basic reproduction number R 0 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 E 0 when R 0 < 1 and when R 0 > 1, 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] [9] [10] [11] [12] [13] [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] [16] [17] [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 E j h . that represents the number of the exposed individuals in region j who are carriers of the disease and have no symptoms, * A compartment Q j h . 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 (δ j 1 ), * The death rate induced by the individuals exposed to quarantine in region j( δ j 2 ), * The death rate induced by the infected animals in region j (δ j 3 ), 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] [20] [21] [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: The animal population is divided into two different classes, namely: Susceptible animals S j a (t) and Infected animals I j a (t). 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 multiregions 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. 2 Corona disease, methods of its spread and the proposed safety protocols 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]. 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] . Our model is divided into two mathematical systems, the system S Similarly, we write the number of animals in each region j: susceptible animals S j a,k and infected animals I j a,k , respectively as functions of discrete-time k. The total animal class is denoted by N j a,k with N j a,k = S j a,k + I j a,k . The following system describes the propagation of the disease in a given population for different regions: ,0 ≥ 0 are the given initial states. The susceptible human population (S j h ) 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 b j h in region j and decreased by an effective contact with the exposed humans at a rate β j 1 (β j 1 represents the transmission rate from susceptible individuals to exposed individuals) and β j 2 ( β j 2 represents the transmission rate from susceptible individuals to infected individuals) and β j 3 due to infection from the infected animals. Finally, the susceptible human population suffers natural mortality (at a rate μ j h ) in the same region and in a different region. The exposed humans (E j h ) are the individuals carrying the disease and have no symptoms in region j. This compartment is increased at the rates β j 1 , β j 2 and β j 3 . The number of some exposed individuals decreases when the exposed individuals become infected at a rate α j 1 (α j 1 represents the transmission rate from exposed individuals to infected individuals). It is decreased by natural death (at a rate μ j h ). The infected humans (I j h ) are the individuals infected by the Corona disease. The infected population is increased by the exposed individuals at a rate α j 1 and decreased by infected individuals who are subjected to quarantine at a rate α j 2 . It is decreased by natural death (at the rate μ j h ) and δ j 1 is the death rate of infected human population. The quarantined humans (Q j h ) are the individuals who are subjected to quarantine in region Ω j . The quarantined population is increased by infected individuals (at a rate α j 2 ). It is decreased by natural death (at the rate μ j h ) and the death rate δ j 2 of the infected humans. The recovered humans (R j h ) represent the recovered individuals in region j . It is also increased when the quarantined individuals become recovered at a rate α j 3 . The number of the recovered individuals of this population decreases by natural death (at the rate μ j h ). The susceptible animals (S j a ) 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 S j a at a rate b j a in region j and decreased by an effective contact with the infected animals at a rate β j 4 (which represents the transmission rate from susceptible animals to infected animals) and β j 5,r (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 μ j a ) in region j. The infected animals (I j a ) refer to the infected animals in region j. The infected animal's compartment is increased by the rate β j 4 in region j and β j 5,r . Finally, the infected animals suffer natural mortality (at a rate μ j a ) and infectious death rate δ j 3 of infected animals in region j. a,k and I j a,k 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 . The total population size at time k in region j is denoted by The strategy of control that we proposed aims at minimizing the number of exposed humans E j h,k , infected humans I j h,k , infected animals I j a,k and maximize the number of the quarantined individuals Q j h,k during the time step k = 0 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 u j represents the effort of awareness campaigns in region j 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 (1−u j ) is used to reduce the force of infections. The second control v j represents the security campaigns and health measures to prevent the movement of individuals from region j to region r (where, r = j, 1 ≤ r ≤ p and p is the total number of the regions) and the term (1 − v j ) is used to reduce the force of infections. The control w j represents the effort to encourage the exposed individuals to join quarantine centers in the same region. Finally, θ j measures the effort to dispose of the infected animals. So, the controlled mathematical system is given by the following system of difference equations. a,0 ≥ 0, and I j a,0 ≥ 0 are the given initial states. Then, the problem is to minimize the objective functional Where the parameters A In other words, we seek the optimal controls u j , v j ,w j and θ j such that Where U ad is the set of admissible controls defined by The sufficient condition for the existence of an optimal control (u j * k , v j * k , w j * k , θ j * k ) for problem (2-3) comes from the following theorem. Theorem 1 There exists the optimal control u j * k , v j * k , w j * k and θ j * k ) such that subject to the control system (3) We apply the discrete version of Pontryagin's Maximum Principle [17, [24] [25] [26] [27] [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 H j k at time step k, defined by where f j i,k+1 is the right side of the system of difference equations (3) concerning the region j of the ith state variable at time step k + 1.. Proof The Hamiltonian at time step k in region j is given by For k = 0, 1, . . . , T − 1, 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 that For k = 0, 1, . . . , T − 1, the optimal controls u * j However, the control attached to this case will be eliminated and removed. By the bounds in U ad of the controls, it is easy to obtain u * j k , v * j k , w * j k and θ * j k in the form of (11-12-13-14). 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 k = 0 and k = T . 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) . 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: j = 1 and j = 2 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. 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 u 1 i.e. the implementation a healthy protocol programs and awareness campaigns in region 1 to protect susceptible . 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 To realize this strategy, we apply only the controls u 1 and v 1 i.e. the implementation of awareness campaigns in region 1 that aim at protecting individuals from being 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. 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 u 1 , v 1 and w 1 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. 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 u 1 , v 1 ,w 1 and θ 1 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. 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] [30] [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 ∈ {1, 2, 3, 4}). Given two competing strategies E and F, where strategy F has higher effectiveness than strategy E (T A(F)>T A(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 i = 1, 2, 3, 4 by: T A(i) = where A 1 1 , A 1 2 , A 1 3 and A 1 4 corresponds to the person unit cost of the four possible interventions, while (E 1 * h,k , I 1 * h,k , I 1 * a,k )is the optimal solution associated to the optimal control (u 1 * , v 1 * , w 1 * , θ 1 * ). Using the simulation results and A 1 1 = A 1 2 = A 1 3 = A 1 4 = 1, we ranked, in the Table 3 our control strategies in order of increased numbers of averted infections. Strategy 4 is compared with strategy 3 with respect to increased effectiveness, in reference to Table 3 Since I C E R(3)