www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Mathematical modeling and analysis of Covid-19 infection spreads with restricted optimal treatment of disease incidence D. Pal1,∗, D. Ghosh2, P.K. Santra3, G.S. Mahapatra2 1Chandrahati Dilip Kumar High School (H.S.) Chandrahati 712504, West Bengal, India 2Department of Mathematics, National Institute of Technology Puducherry Karaikal-609609, India 3Maulana Abul Kalam Azad University of Technology Kolkata-700064, India ∗Correspondence email: pal.debkumar@gmail.com Received: 13 August 2020, accepted: 14 June 2021, published: 17 July 2021 Abstract— This paper presents the current situa- tion and how to minimize its effect in India through a mathematical model of infectious Coronavirus disease (COVID-19). This model consists of six compartments to population classes consisting of susceptible, exposed, home quarantined, government quarantined, infected individuals in treatment, and recovered class. The basic reproduction number is calculated, and the stabilities of the proposed model at the disease-free equilibrium and endemic equilib- rium are observed. The next crucial treatment con- trol of the Covid-19 epidemic model is presented in India’s situation. An objective function is considered by incorporating the optimal infected individuals and the cost of necessary treatment. Finally, optimal control is achieved that minimizes our anticipated objective function. Numerical observations are pre- sented utilizing MATLAB software to demonstrate the consistency of present-day representation from a realistic standpoint. Keywords-Novel coronavirus; SEHGIR model; Basic Reproduction number; Stability; Optimal con- trol I. INTRODUCTION Recently, the coronavirus disease has turned out to be a pandemic over almost the whole world. The basic indication of this infection is ordinary fever, cough, and breathing problems. This virus also showed the capability to produce serious health problems among a specific group of individuals, including the aged populace as well as patients with cardiovascular disease and diabetes [1]. However, a clear picture of the nature of this Copyright: © 2021 Pal et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 infection spreads with restricted optimal treatment of disease incidence, Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 1 of 20 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... epidemiology is still being explained [2]. This virus transmits from human to human as a result, Covid-19 disease spread across the globe, and the total number of active Covid-19 cases increases day by day ([3]-[10]). In particular, India has become the stiffest affected country with Covid-19 endemic [11] due to its very high pop- ulation densities. The figure for positive Covid-19 infection started to increase from 4th March 2020. As of 8th May 2020, a total of 59690 confirmed COVID cases, together with 17887 recovered and 1986 deaths in India [12]. Different precaution measures ([13]-[16]) have been taken by the In- dian Government to maintain social distance [17] among the huge numbers of the population in India. There is no specific medicine for Covid-19 infection to date. Therefore, doctors recommended different treatments via medication to COVID-19 patients depending on their symptoms. The therapy and vaccine yet to get, spread of Coved 19 diseases can be restricted via appropriate precautionary measures like quarantined mechanisms ([18]-[20]), individual safeguard from the infected individ- ual by using social distancing [21], etc. As the Covid-19 virus spread very quickly throughout the world, so various mathematical models depending on the pandemic outbreak ([22]-[32], [33], [34], [35]) have been performed already. Wu et al. [36] studied the dynamics behind the spread of Covid- 19 virus world-wise using SEIR model. Read et al. [37] developed a Covid-19 SEIR model based on Poisson-distributed daily time augmentations. Paul et al. [38] presented a mathematical model on Covid-19 incorporating the different safety strate- gies to protect the citizens from the virus. Sardar et al. [39] proposed a mathematical model to identify the lockdown effect of the spreading of Covid-19 disease in India. Pal et al. [40] explored a Covid-19 based SEQIR model to understand India’s disease situation. This paper introduces a six-compartmental Covid-19 infection model by separating the total populace into six classes, purposely susceptible, exposed, home quarantined, government quaran- tined, infected individuals in treatment as well as recovered class. We introduce treatment control in the model to assimilate realistic and biologically significant in the pandemic situation. A brief de- scription of the necessary and sufficient conditions for the existence of multi-objective optimal control is provided in Section 2. The model derivation and preliminaries are explained in Section 3. The basic properties of our proposed model structure are discussed in Section 4. In section 5, we introduce the concept of the basic reproduction number (R0) [41]. Next, we deal with disease-free equilibrium (DFE) (E0) and endemic equilibrium (E1) points of the system. It is clear that Covid- 19 infection is not only community health trouble [42] but also a tremendous societal and monetary shock on the developing countries. Therefore, it is an essential concern to control ([43]-[46]) the spread of Covid-19 infection in India by adopting several optimal control policies. In Section 6, we have formulated the Covid-19 epidemic model with control treatment. This section provides us a procedure to find optimal control [47] u(t) that increases the recovery rate as well as minimizes the cost associated with the treatment. Analytical results are obtained in the previous sections are numerically verified in Section 7 with the help of realistic values of the model parameters using MATLAB. Lastly, a general conclusion about our proposed model structure is provided in Section 8. II. MULTI-OBJECTIVE OPTIMAL CONTROL Suppose x(t) ∈ X ⊂ Rn represents the state variables of a system and u(t) ∈ U ⊂ Rm represents the control variables at time t, with t0 ≤ t ≤ tf . An optimal control problem consists of finding a piecewise continuous control u(t) and the associated state x (t) that optimizes a cost function J [u(t),x (t)]. The majority of mathemat- ical models that use the optimal control theory rely on Pontryagin’s Maximum Principle, a first-order condition for finding the optimal solution. Theorem 1. (Pontryagin’s Maximum Principle [48]) If u∗ (t) and x∗(t) are optimal for the problem max u J [u(t),x (t)] , (1) Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 2 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... where J [u(t),x (t)] = max u ∫ tf t0 f (t,x (t) ,u (t)) dt, subject to { dx dt = g (t,x (t) ,u (t)) , x (t0) = x0, then there exists a piecewise differentiable adjoint variable λ (t) such that H(t,x∗(t),u(t),λ(t))≤H(t,x∗(t),u∗(t),λ(t))) for all controls u at each time t, where the Hamiltonian H is given by H(t,x(t),u(t),λ(t)) = f (t,x (t) ,u (t)) + λg (t,x (t) ,u (t)) (2) and{ λ ′ (t) = − ∂H (t,x∗ (t) ,u∗ (t) ,λ (t)) ∂x , λ (tf ) = 0. While Pontryagin’s Maximum Principle gives the necessary conditions for the existence of an optimal solution, the following theorem provides sufficient conditions. Theorem 2. (Arrow Sufficiency Theorem [49]) For the optimal control problem (1), the conditions of the maximum principle are sufficient for the global minimization of J [u(t),x (t)], if the minimized Hamiltonian function H, defined in (2), is convex in the variable x for all t in the time interval [t0, tf ] for a given λ. One of the major side effects of vaccina- tion/treatment is the creation of drug resis- tant virus/bacteria which eventually leads to drug failure (due to ineffectiveness of the vac- cine/treatment). Optimal control has been used to curb the creation of drug resistant virus/bacteria or drug failure (at the same time reducing the cost of treatment or vaccination) by imposing a condition that monitors the global effect of the vaccina- tion/treatment program. Hence if x(t) represents the group of individuals to be vaccinated/treated and u(t) ∈ U represents the control on vaccina- tion/treatment, where the control set U is given by U ={u(t) :v0≤ u(t)≤v1, Lebesgue measurable}, then, the following objective functions are to be minimized simultaneously: I1 (u) = ∫ tf t0 x (t) dt and I2 (u) = ∫ tf t0 um (t) dt, for m > 0, and the optimal solution can be represented as min u {I1 (u) ,I2 (u)} . In general, there does not exist a feasible solution that minimizes both objective functions simulta- neously. Hence, the Pareto optimality concept is used to simultaneously find optimal control u∗ that minimizes both objective functions. III. DERIVATION AND PRELIMINARIES OF COVID-19 MODEL This section develops a mathematical model of COVID-19 transmission with the subsequent suppositions: The underlying human population is split up into six mutually exclusive compartments, namely, susceptible (S), exposed (infected but not yet infectious) (E), home quarantined population (population were exposed to the virus but viewing light symptoms of coronavirus disease and stay at home isolation) (H), government quarantined pop- ulation (population was infective in symptomatic phase, i.e., showing symptoms of coronavirus dis- ease and stay at Government observation places for isolation) (G), infected (I), and recovered class (R) (infectious people who have cleared or recovered from coronavirus infection). Therefore, the total human population N(t) = S(t) + E(t) + H(t) + G(t) + I(t) + R(t). This model involves certain assumptions which consist of the following: (i) The susceptible population (S) comprises individuals who have not yet been infected by Covid-19, but may be infected through Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 3 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... contact with both types of home quarantined (H) and government quarantined (G) peo- ple. (ii) The exposed population (E) comprises indi- viduals infected with Covid-19 infection but not infectious. (iii) The infective population in home quaran- tined phase (H) comprises individuals who have Covid-19 infection with light symp- toms (but capable of infecting) and quaran- tined at home for isolation. (iv) The infective population in the symptomatic phase (G) comprises individuals who have developed Covid-19 infection with compli- cations and various symptoms but their test report yet not come positive and are quar- antined by the Government facility for iso- lation. (v) The infected population (I), whose COVID 19 test is positive clinically and stayed at hospital for treatment (incapable of infect- ing others). The infected individuals com- ing from home and Government quarantined compartments if their test report comes pos- itive. (vi) The recovered class (R) consists of those who become healed from the disease by treatment or quarantined program. (vii) The susceptible individuals become infected by adequate contact with infective individu- als in the asymptomatic phase (home quar- antined) and symptomatic phase (Govern- ment quarantined), and enter into the ex- posed class. The susceptible population is also decreased due to natural death. (viii) The exposed population is entered into the home quarantined, government quarantined, and infected population, respectively. The said population is also diminished due to natural death. (ix) One part of home quarantined individuals enters into the infected population, and the other becomes recovered. This population is also decreased by natural death. (x) One part of the government quarantined in- dividuals enters into the infected population, and the other becomes recovered. This indi- vidual is also decreased by natural death. (xi) One part of the infected population enters into the recovered class. Other individuals are decreased due to infection and natural death. (xii) Home quarantined (asymptomatic), govern- ment quarantined (symptomatic), and the infected population recover from the coron- avirus disease and enters into the recovered class. The recovered population diminishes by natural death. The parameters of the Covid-19 model are pre- sented as follows: Λ : The recruitment rate of susceptible from embedding population. α1 : The coefficient of transmission rate from home quarantined to susceptible individuals, and the expression gives the transmission rate: α1H(t)S(t). α2 : The coefficient of transmission rate from Government quarantined popula- tion to susceptible individuals, and the transmission rate is: α2G(t)S(t). β1 : The fraction of exposed individuals that will start to show light symptoms of Covid-19 (but remains capable of infect- ing others) and move to the class H. β2 : The rate at which the exposed individuals become infected by Covid 19 infection and move to the class I. β3 : The fraction of the exposed individ- uals that will start to show symptoms of infection and move to the class G. γ2 : The rate of home quarantined in- dividuals eventually show disease symp- toms and move to class I. γ1 : The recovery rate of the home quarantined population H. σ2 : The rate at which government quar- antined individuals eventually show dis- ease symptoms and move to class I. σ1 : The recovery rate of the Government quarantined population G. Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 4 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... d2 : The disease-related death rate of infective population in the infected phase I � : The recovery rate of the infected population I. d1 : The natural death rate of all human epidemiological classes. In our proposed Covid-19 model, S(t), E(t), H(t), G(t), I(t), and R(t) denote the numbers of susceptible, exposed, home quarantined, gov- ernment quarantined, infected, and recovered, re- spectively. Through the contact between suscep- tible and home quarantined populations, a part of the susceptible population, i.e., α1H(t)S(t), becomes infected and enters into exposed class. Similarly, through the contact between susceptible and government quarantined populations, a part of the susceptible people, i.e., α2G(t)S(t), becomes infected and enters into the exposed category. The fraction of the home quarantined population γ2 will start to show symptoms of and move to the class I. Another portion of the home quarantined population γ1 is recovered from infection due to treatment or quarantined process and move to the recovered class R. Similarly, a fraction of the government quarantined community σ2 will start to show symptoms of Covid 19 and move to the class I. Other portion of the home quarantined population σ1 is recovered from infection due to treatment or quarantined process and move to the recovered class R. A fraction of infected individuals � is recovered from infection through treatment in Hospital and move to recovered class R. Another fraction d2 of the infected individuals is diminished due to the disease-related death rate of the infective population. From every class, a part of the inhabitants is reduced at the natural death rate d1. We diagrammatically represent the flow of in- dividuals from one class to another in Fig. 1. Therefore, our proposed mathematical model of the Covid-19 infection is presented through the Fig. 1. Pictorial representation of proposed Covid 19 model for Indian scenario following set of non-linear differential equation dS dt = Λ − (α1H + α2G) S −d1S, (3) dE dt = (α1H + α2G) S − (β1 + β2 + β3 + d1) E, dH dt = β1E −γ1H −γ2H −d1H, dG dt = β3E −σ1G−σ2G−d1G, dI dt = β2E + γ2H + σ2G−d1I −d2I − �I, dR dt = �I −d1R + γ1H + σ1G; with initial conditions: S(0) > 0,E(0) ≥ 0,H(0) ≥ 0, G(0) ≥ 0,I(0) ≥ 0,R(0) ≥ 0. (4) The SEHGIR model formulation (3) can be Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 5 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... rewritten as: dS dt = Λ − (α1H + α2G) S −d1S, (5) dE dt = (α1H + α2G) S −AE, dH dt = β1E −BH, dG dt = β3E −CG, dI dt = β2E + γ2H + σ2G−DI, dR dt = γ1H + σ1G + �I −d1R; where A = β1 + β2 + β3 + d1, B = γ1 + γ2 + d1, C = σ1 + σ2 + d1, D = d1 + d2 + �, with initial conditions (4). IV. FUNDAMENTAL PROPERTIES A. Positivity of the solutions Theorem 3. Each solution of the proposed system (5) under conditions (4) satisfy S(t) > 0, E(t) ≥ 0, H(t) ≥ 0, G(t) ≥ 0, I(t) ≥ 0, R(t) ≥ 0 for all values of t ≥ 0. Proof: The first equation of the system (5), can be written dS dt = Λ − (α1H + α2G) S −d1S = Λ −ψS; where ψ = (α1H + α2G) − d1. Thereafter by integration, we obtain the following expression S(t)=S(0) exp ( − ∫ t 0 ψ(s)ds ) +Λ exp ( − ∫ t 0 ψ(s)ds )∫ t 0 e ∫ t 0 ψ(v)dvds>0. Hence S(t) is non-negative for all t. From the next equation of (5), we get, dE dt ≥−AE. This equation provides E(t) ≥ E(0) exp(−At) ≥ 0. Also, from the remaining equations and with the help of initial conditions, we obtain H(t) ≥ H(0) exp(−Bt) ≥ 0, G(t) ≥ G(0) exp(−Ct) ≥ 0, I(t) ≥ I(0) exp(−Dt) ≥ 0, as well as R(t) ≥ R(0) exp(−d1t) ≥ 0. So, it is observed that S(t) > 0, E(t) ≥ 0, H(t) ≥ 0, G(t) ≥ 0, I(t) ≥ 0, R(t) ≥ 0 for all values of t ≥ 0. Hence the theorem. B. Invariant region Theorem 4. The feasible region Γ defined by Γ = { (S,E,H,G,I,R) ∈ R6+ : 0 < N ≤ Λ η } , where η = min{d1,d1 + d2} is positively invari- ant for the system (3). Proof: Let ((S(0),E(0),H(0),G(0),I(0),R(0)) ∈ Γ. Adding the equations of the system (3) we obtain dN dt = Λ−d1S−d1E−d1H−d1H−d1G−(d1+d2)I−d1R. Therefore, dN dt +ηN = Λ−(d1−η)S−(d1−η)E−d1H −(d1−η)H−(d1−η)G−(d1 +d2−η)I (6) −(d1−η)R ≤ Λ, where η = min{d1,d1 + d2}. The solution N(t) of the differential equation (6) has the following property, 0 < N(t) ≤ N(0) exp(−ηt) + Λ η (1 − exp(−ηt)), where N(0) represents the sum of the initial values of the variables. As t → ∞, we have 0 < N(t) ≤ Λ η . Also, if N(0) ≤ Λ η then N(t) ≤ Λ η Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 6 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... for all t. This means Λ η is the upper bound of N. On the other hand, if N(0) > Λ η implies N(t) will decrease to Λ η . This means that if N(0) > Λ η , the solution (S(t),E(t),H(t),G(t),I(t),R(t)) enters Γ or approaches it asymptotically. Hence it is positively invariant under the flow induced by the systems (3) and (4). Thus in Γ the mathematical model (3) with ini- tial conditions (4) is well-posed epidemiologically. Hence it is sufficient to study the dynamics of the model in Γ. V. EXISTENCE OF EQUILIBRIUM AND STABILITY ANALYSIS In this section, we will study the existence and stability behavior of the system (3) at various equilibrium points. The equilibrium points of the system (2.1) are: (i) Disease-free equilibrium (DFE): E0 ( Λ d1 , 0, 0, 0, 0, 0 ) , (ii) Endemic equilibrium: E1(S ∗,E∗,H∗,G∗,I∗,R∗). A. The basic reproduction number The basic reproduction number (BRN) ([50]- [53]) of the system (3) will be obtained by the next-generation matrix method [54]. Let z = (E(t),H(t),G(t),I(t),S(t),R(t))T , the proposed Covid-19 system (3) can be written in the following form: dz dt = z(z) −υ(z); where z(z) =   (α1H + α2G) S 0 0 0 0 0   , υ(z) =   AE −β1E + BH −β3E + CG −(β2E + γ2H + σ2G) + DI −(γ1H + σ1G + �I) + d1R −Λ + (α1H + α2G) S + d1S   . The Jacobian matrices of z(z) and υ(z) at the DFE E0 are as follows, respectively: Dz(E0) =  F4×4 0 00 0 0 0 0 0   , Dυ(E0) =  V4×4 0 00 0 0 0 0 0   , where F =   0 Λα1 d1 Λα2 d1 0 0 0 0 0 0 0 0 0 0 0 0 0   , V =   A 0 0 0 −β1 B 0 0 −β3 0 C 0 −β2 −γ2 −σ2 D   . Following [54], R0 = ρ ( FV −1 ) where ρ is the spectral radius of the next-generation matrix (FV −1). Thus, from the model (3), we have the following expression for BRN R0 : R0 = Λ d1 1 ABC [α1β1C + α2β3B]. Notice that Λ d1 is the number of susceptibles at the DFE. B. Existence of endemic equilibrium E1(S ∗,E∗,H∗,G∗,I∗,R∗) In this section, we will analyze the ex- istence of a non-trivial endemic equilibrium E1(S ∗,E∗,H∗,G∗,I∗,R∗) of the system (3). To find the endemic equilibrium of the system (3), we consider the following: S(t) > 0,E(t) > 0,H(t) > 0, G(t) > 0,I(t) > 0,R(t) > 0 Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 7 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... and dS dt = 0, dE dt = 0, dH dt = 0, dG dt = 0, dR dt = 0. (7) From the third, fourth, fifth, and sixth equation of (7) we obtain H∗ = β1E ∗ B , G∗ = β3E ∗ C , I∗ = E∗ D [β2 + γ2β1 B + β3σ2 C ], R∗ = E∗ d1 [ � D {β2 + γ2β1 B + β3σ2 C }+ γ1β1 B + β3σ1 C ] . Now from dE dt = 0 and using the values of H∗and G∗, we get, S∗ = Λ d1R0 > 0. Again, putting the value of S∗ in the first equation of (7) we gain, E∗ = Λ A [ 1 − 1 R0 ] . Hence, E∗ has a unique positive solution iff R0 > 1. Summarizing the above discussions, we arrive at the following result. Theorem 5. The system (3) has a DFE E0( Λ d1 , 0, 0, 0, 0, 0), which exists for all pa- rameter values. If R0 > 1 the system (3) also admits a unique endemic equilibrium E1(S ∗,E∗,H∗,G∗,I∗,R∗). C. Asymptotic behavior For the stability of DFE E0( Λd1 , 0, 0, 0, 0, 0) we consider the theorems given below Theorem 6. The DFE E0 of the system (3) is locally asymptotically stable if R0 < 1. Proof: See Appendix A. Theorem 7. The diseases free equilibrium (DFE) E0( Λ d1 , 0, 0, 0, 0, 0)is globally asymptotically sta- ble (GAS) in R6+ for the system (3) if R0 < 1 and becomes unstable if R0 > 1. Proof: We rewrite the system (3) as dX dt = F(X,V ), dV dt = G(X,V ), G(X, 0) = 0, where X = (S,R) ∈ R2 (the number of uninfected individuals compartments), V = (E,H,G,I) ∈ R4 (the number of infected individ- uals compartments), and E0( Λd1 , 0, 0, 0, 0, 0) is the DFE of the system (3). The global stability of the DFE is guaranteed if the following two conditions are satisfied: (i) For dX dt = F(X, 0), X∗ is globally asymptot- ically stable in R2. (ii) G(X,V ) = BV − Ĝ(X,V ), Ĝ(X,V ) ≥ 0 for (X,V ) ∈ Ω, where B = DV G(X∗, 0) is a Metzler matrix, and Ω is the positively invariant set to the model (3). Following Castillo-Chavez et al. [55], we check for aforementioned conditions. For system (3), F(X, 0) = [ Λ −d1S 0 ] , B =   −A Λα1 d1 Λα2 d1 0 β1 −B 0 0 β3 0 −C 0 β2 γ2 σ2 −D   and Ĝ(X,V ) =   ( Λ d1 −S)(α1H + α2G) 0 0 0   . Clearly, Ĝ(X,V ) ≥ 0 (using Theorem 2), when- ever the state variables are inside Ω (the positively invariant set of the model (3)). Again, it is clear that X∗ = ( Λ d1 , 0)T is a globally asymptotically stable equilibrium of the system dX dt = F(X, 0). Hence, the theorem follows. Theorem 8. The endemic equilibrium point E1(S ∗,E∗,H∗,G∗,I∗,R∗) of the system (3) is locally asymptotically stable if R0 > 1, B1B2 − B3 > 0 and B1B2B3 −B21B4 −B 2 3 > 0 Proof: See Appendix B. Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 8 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... VI. PROPOSED COVID-19 MODEL WITH CONTROL In this section, the primary focus is to set up an optimal control problem of the epidemic model (3). In the present situation of the Covid- 19 outbreak, it is highly essential to construct an optimal control problem so that the total amount of drug is minimized. Here we take one control variable u(t) on the recovery rate of the infectious individuals in the infected phase with treatment in the hospital. Therefore, our epidemic model with one control and the same initial conditions (4) becomes: dS dt = Λ − (α1H(t) + α2G(t)) S(t) −d1S(t), dE dt = (α1H + α2G) S − (β1 +β2 +β3 +d1) E, dH dt =β1E −γ1H −γ2H −d1H, dG dt =β3E −σ1G−σ2G−d1G, dI dt =β2E(t) + γ2H(t) + σ2G(t) (8) −(d1 + d2)I(t) −u(t)I(t), dR dt =u(t)I(t) + γ1H(t) + σ1G(t) −d1R(t). The control function u(t), 0 ≤ u(t) ≤ 1 represents the fraction of the infected individuals who are identified and will be treated. When u(t) is close to 1 then the treatment failure is low, but the implementation cost is high. For the model (8), the single-objective cost functional to be minimized is given by the objective functional ([56]-[59]) J(u(t)) = ∫ tf 0 [G1I + 1 2 G2u 2]dt; (9) with G1 > 0 and G2 > 0, where we want to minimize the infectious group I while also keeping the cost of treatment u(t) low. The term G1I represents the cost of infection, while the term 1 2 G2u 2 represents the cost of treatment. The goal is to find an optimal control, u∗, such that J (u∗) = min{J(u) : u ∈ U} , (10) where U ={u : u is Lebesgue measurable, 0 ≤ u ≤ 1, t ∈ [0, tf ] } (11) Applying the Pontryagins Maximum Principle, we have the following result( S ∗ ,E ∗ ,H ∗ ,G ∗ ,I ∗ ,R ∗ ) of the system (8), that minimizes J(u) over U. Theorem 9. There exists an optimal control u∗and corresponding solutions ( S ∗ ,E ∗ ,H ∗ ,G ∗ ,I ∗ ,R ∗ ) of the system (8), that minimizes J(u) over U. Furthermore, there exist adjoint functions λi(t), i = 1, 2, 3, 4, 5, 6, such that dλ1 dt = (λ1−λ2)(α1H+α2G)+λ1d1, dλ2 dt = (λ2−λ3)β1 +(λ2−λ5)β2 +(λ2−λ4)β3 +λ2d1, dλ3 dt = (λ1−λ2)α1S+(λ3−λ6)γ1 +(λ3−λ5)γ2 +λ3d1, dλ4 dt = (λ1−λ2)α2S+(λ4−λ6)σ1 +(λ4−λ5)σ2 +λ4d1, dλ5 dt = (λ5−λ6)u+(d1 + d2)λ5−G1, dλ6 dt =d1λ6; with transversality conditions λi(tf ) = 0, i = 1, 2, 3, 4, 5, 6, and the control u∗ satisfies the optimality condition u∗ = min{max{0, (λ5 −λ6)I ∗ G2 }, 1}. Proof: The Hamiltonian is defined as follows: Ĥ =G1I+ 1 2 G2u 2 +λ1[Λ−(α1H+α2G) S−d1S] +λ2[(α1H + α2G) S −AE] (12) +λ3[β1E −BH] + λ4[β3E −CG] +λ5[β2E + γ2H + σ2G− (d1 + d2)I −uI] +λ6[uI + γ1H + σ1G−d1R], Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 9 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... where λi (i = 1, 2, 3, 4, 5, 6) are the adjoint functions to be determined. The form of the adjoint equations and transver- sality conditions are expected results from Pon- tryagin’s Maximum Principle [60]. The adjoint system can be obtained as follows: dλ1 dt = − ∂Ĥ ∂S = (λ1−λ2) (α1H+α2G)+λ1d1, dλ2 dt = − ∂Ĥ ∂E = (λ2 −λ3)β1 + (λ2 −λ5)β2 +(λ2 −λ4)β3 + λ2d1, dλ3 dt = − ∂Ĥ ∂H = (λ1 −λ2)α1S + (λ3 −λ6)γ1 +(λ3 −λ5)γ2 + λ3d1, dλ4 dt = − ∂Ĥ ∂G = (λ1 −λ2)α2S + (λ4 −λ6)σ1 +(λ4 −λ5)σ2 + λ4d1, dλ5 dt = − ∂Ĥ ∂I = (λ5−λ6)u+(d1 +d2)λ5−G1, dλ6 dt = − ∂Ĥ ∂R = d1λ6. (13) The transversality conditions (or boundary condi- tions) are λi(tf ) = 0, i = 1, 2, 3, 4, 5, 6. (14) By the optimality condition, at u = u∗(t) we have ∂Ĥ ∂u = G2u ∗ − (λ5 −λ6)I ∗ = 0 ⇒ u∗(t) = (λ5−λ6)I ∗ G2 . (15) By using the bounds for the control u(t), we get u∗ =   (λ5−λ6)I ∗ G2 , if 0 ≤ (λ5−λ6)I ∗ G2 ≤ 1. 0, if (λ5−λ6)I ∗ G2 ≤ 0. 1, if (λ5−λ6)I ∗ G2 ≥ 1. In compact notation: u∗ = min { max { 0, (λ5 −λ6)I ∗ G2 } , 1 } . (16) Using (16), we obtain the following optimality system: dS dt = Λ − (α1H + α2G) S −d1S, (17) dE dt = (α1H + α2G) S −AE, dH dt =β1E −BH, dG dt =β3E −CG, dI dt =β2E + γ2H + σ2G− (d1 + d2)I −min { max { 0, (λ5 −λ6)I ∗ G2 } , 1 } I, dR dt = min { max { 0, (λ5 −λ6)I ∗ G2 } , 1 } I +γ1H + σ1G−d1R, dλ1 dt = (λ1 −λ2) (α1H + α2G) + λ1d1, dλ2 dt = (λ2 −λ3)β1 + (λ2 −λ5)β2 +(λ2 −λ4)β3 + λ2d1, dλ3 dt = (λ1 −λ2)α1S + (λ3 −λ6)γ1 +(λ3 −λ5)γ2 + λ3d1, dλ4 dt = (λ1 −λ2)α2S + (λ4 −λ6)σ1 +(λ4 −λ5)σ2 + λ4d1, dλ5 dt = (λ5−λ6) min{max{0, (λ5−λ6)I ∗ G2 }, 1} +(d1 + d2)λ5 −G1, dλ6 dt =d1λ6; subject to the following conditions: S(0) > 0, E(0) ≥ 0, H(0) ≥ 0, G(0) ≥ 0, I(0) ≥ 0, R(0) ≥ 0 and λ1(tf ) = 0, λ2(tf ) = 0, λ3(tf ) = 0, λ4(tf ) = 0, λ5(tf ) = 0, λ6(tf ) = 0. This completes the proof. Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 10 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 −10 0 0.5 1 1.5 x 10 −10 0 0.5 1 1.5 2 2.5 3 3.5 α 1 α 2 R e p ro d u c ti o n N u m b e r( R 0 ) Fig. 2. Sensitivity of R0 to α1 and α2, rest of the parameters are based on Table 1 VII. NUMERICAL SIMULATIONS The current section presents some computer simulations to assess the proposed model’s ap- plicability for the Covid-19 scenario. The simu- lation is carried out based on available data of pandemic infection in India. Also, these numerical simulation is very much crucial from a practical viewpoint. Estimating the parameters of the model for India, we have studied the proposed Covid-19 system. The main objective is to study the effects of two quarantined population parameters α1 and α2, to show the impact of these parameters on the pandemic curve through the graphical presen- tation. By changing the values of the mentioned parameters, we observe the infected population’s behavior for 60 days from 2nd April for their particular base values. Table 1 and Table 2 give the values of the model parameters and initial population density, respectively. Based on Table 1, the BRN is R0 = 3.0909, which is much greater than one. Hence the in- fection spread so quickly in India. Therefore, it needs to take the right policy to reduce the value of R0 much less to 1. For the proposed model, a graphical presentation of R0 to α1and α2 is given in Figure 2. TABLE I MODEL PARAMETERS FOR COVID-19 SYSTEM Parameters Values (Unit) Data Source Λ 50000 day−1 [61] α1 2 × 10−10 day−1 Estimated α2 1 × 10−10 day−1 Estimated β1 0.4 day−1 Assumed β2 1 × 10−6 day−1 Assumed β3 0.05 day−1 Assumed γ1 0.15 day−1 Estimated γ2 0.0028 day−1 Estimated σ1 0.15 day−1 Estimated σ2 0.002 day−1 Estimated � 0.06 day−1 Estimated d1 2 × 10−5 day−1 [61] d2 0.001 day−1 Estimated TABLE II PRELIMINARY POPULATION DENSITY FOR COVID-19 MODEL S(0) E(0) H(0) G(0) I(0) R(0) 12×108 2×105 2×105 5×104 1649 5×104 In Figure 3, the ’red’ curve presents an infected individual for this proposed model, and the bar diagram is the actual infected individual as per our available data. Figure 3 depicts that the ac- tual infected individual almost coincides with our Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 11 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... 1/4 10/4 20 /4 30/4 10/5 20/5 29/5 0 2 4 6 8 10 x 10 4 Time P o p u la ti o n I Fig. 3. Time series of infected population with parameter values and initial conditions from Table 1 and 2 during 1/4/2020 to 29/5/2020. 0 50 100 150 200 250 300 350 400 450 500 0 2 4 6 8 10 12 14 16 18 x 10 5 Time P o p u la ti o n I α 1 =0.0000000002 α 1 =0.00000000019 α 1 =0.00000000018 Fig. 4. Time series of the infected population with α2 = 1 × 10−10 for different values of α1 and other input values taken from Table 1 and 2. proposed model curve from 1st April to 29th May 2020. Therefore, the proposed Covid-19 model is best fitted to the current situation of India. For fixed α2 if we gradually decrease α1, the infected individuals is also reduces steadily, which is presented via Figure 4. Therefore, practically if we strictly follow the home quarantined restriction, then naturally α1 de- crease, and also the pick of the infected individual reduces. It is also observed from Figure 4 that for α1 = 2 × 10−10, the pick of the disease reached almost after 160 days from 1st April 2020, and the height number of infected cases around 1700000. For α1 = 1.9 × 10−10 and α1 = 1.8 × 10−10 the pick of infection reached almost after 175 and 220 days, respectively from 1st April 2020, and the corresponding height number of infected cases may be around 1200000 and 800000, respectively. Again if we fixed α1 at 2×10−10 and the values of α2 gradually increase, then the infected number of individuals is also gradually increasing, which Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 12 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... 0 50 100 150 200 250 300 350 0 0.5 1 1.5 2 2.5 3 3.5 x 10 6 Time P o p u la ti o n I α 2 =0.0000000001 α 2 =0.0000000002 α 2 =0.0000000003 Fig. 5. Time series of infected population with parameter values and initial conditions from Table 1 and 2 during 1/4/2020 to 29/5/2020. 0 20 40 60 80 100 120 140 160 180 200 0 1000 2000 3000 4000 5000 6000 7000 8000 Time P o p u la ti o n I α 1 =0.00000000001 α 1 =0.00000000005 α 1 =0.0000000001 Fig. 6. Time series of infected population with α1 = α2 using data from Table 1, 2 starting from 1st April to 29th May 2020. is depicted in Figure 5. This situation arises as we increase α2, then the government quarantined technique is slackly ap- plied to the population. In this case α2 = 3×10−10, the pick of infection reached almost 125 days after 1st April 2020, and the total number of highest infected individuals will be around 3000000. Also, we making α1 = α2 = 2 × 10−10, i.e., if Government take a policy to convert all home quarantined individuals into government quaran- tined. In that case, the value of R0 = 1.6368 < 3.0909 is less than the previous value of R0. Therefore, infected individuals are automatic de- creases, which are depicted in Figure 6. Figure 6 shows that if the Government takes said policy, then the maximum number of infected is Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 13 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... 30/4 20/5 10/6 30/6 20/7 7/8 0 0.2 0.4 0.6 0.8 1 Time u Fig. 7. The graph of u with respect to time t based on Table 1, 2 and G1 = 0.005 and G2 = 1000 starting from 30/4/2020 to 7/8/2020. 30/4 20/5 10/6 30/6 20/7 7/8 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 5 Time P o p u la ti o n I Fig. 8. The control diagram for the infected population (I) using data given in Table 1 along with G1 = 0.005 and G2 = 1000 from 30/4/2020 to 7/8/2020. restricted to about 8000. Unfortunately, to arrange this type of quarantined system is not possible for Government since India is a country with large populations. Therefore, the Government has to think above other possible ways to restrict Covid- 19 infection. Therefore, this paper provides a way to restrict the infection by optimal control policy. We try to recover the infected patients by using the mini- mum drug. The present section explores the idea to solve the control problem numerically and will interpret the findings graphically. The boundary value problem in this paper estranged boundary conditions ranges between t = 0 to t = tf . The optimality problem is solved intended for 100 days. Actually, given time tf = 100 represent the period at which the given treatment is stopped. TABLE III INITIAL DENSITIES FOR THE OPTIMAL CONTROL PROBLEM (17) S(0) E(0) H(0) G(0) I(0) R(0) 11.76×108 2×106 5×106 5×105 89987 16×106 The collocation method is the best technique to solve two-point BVPs numerically. The current optimization problem solves numerically using MATLAB for our control problem based on Table 1 and Table 3. Here, we choose weight constants G1 = 0.005 and G2 = 1000, respectively in the objective function given in (17). Figure 7 represents the optimal control graphs for treatment control u. It shows that treatment control is very much neces- sary when the disease prevails. Also, this control function minimizes the cost function J. Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 14 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... 0 10 20 30 40 50 60 70 80 90 100 0 1 2 3 4 5 6 7 8 x 10 6 Time P o p u la ti o n I Fig. 9. Diagram for the infected populace (I) without control using data from Tables 1, 2 and G1 = 0.005, G2 = 1000 from 30/4/2020 to 7/8/2020. 0 10 20 30 40 50 60 70 80 90 100 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 10 8 Time P o p u la ti o n R Fig. 10. The control diagram of the recovered population (R) based on Tables 1, 2 and G1 = 0.005, G2 = 1000 from 30/4/2020 to 7/8/2020. The graphs of the infected population (I) and recovered population (R) with treatment control and without treatment control with respect to time t are presented in Figure 8 to Figure 11, respec- tively. From these figures, we can predict that treat- ment control is exceptionally efficient in reducing COVID infection. Therefore, control acquiesces the best result to control the Covid-19 epidemic outbreak. VIII. DISCUSSIONS AND CONCLUSIONS This paper explores the idea of a six- compartmental Covid-19 infection model fitted for the India scenario. The present model has exhib- ited the effects of different precautions proposed by the administration to control India’s infectious disease. This study has also presented the impact of home and Government quarantined technique Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 15 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... 0 10 20 30 40 50 60 70 80 90 100 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 10 8 Time P o p u la ti o n R Fig. 11. Diagram of the recovered population (R) without of control based on Table 1, 2 and G1 = 0.005, G2 = 1000 during 30th April to 7th August 2020. on the Covid-19 epidemiological model via a non- linear differential equation system. The dynamical behavior of our proposed model is presented at DFE and endemic equilibrium points. We have calculated the BNR and verify that it behaves a crucial role in predicting the stability nature of all possible equilibrium points and the existence of the disease soon. Also, a sensitivity analysis of R0 is carried out to α1 and α2 through Figure 2. This analysis has shown that the parameters α1 and α2 are vital to restrict the spread of infection. The most crucial part related to public health importance is that this paper has built up a suitable optimal control problem to reduce the number of infected individuals. Though the infection may be controlled by reducing the parameters’ values α1 and α2, it is not a long term solution to restrict the spread of the disease. Therefore, we deem the treatment of infected individuals by medicine as a control to diminish the spread of Covid-19 infection. In this paper, we include a quadratic control to quantify this goal. To minimize the objective functional (9), the control function u(t) is considered. This study also numerically verified the theo- retical analysis by using MATLAB software to validate scientific findings through plot compar- ative figures of infected populations with different values of α1 and α2. We have observed from Figure 4 for α1 = 2 × 10−10 and α1 = 1 × 10−10 that the pick of infection will be attained in the mid of September 2020 and around 1700000 may be affected in that time. Again if the administration has taken the policy to cover all the populations under the Government quarantined process (which is impossible for a country like India), i.e., α1 = α2, then Figure 6 show that a maximum number of infected individuals are around 80000, the in- fection may be eradicated within October 2020. The proposed optimal control strategies are ben- eficial to reduce the number of infected popu- lations, which is presented through comparative Figure 8 to Figure 11. It may also be concluded from these figures that the only treatment of an infected individual by medicine may not be the possible way to die out the disease from India and the globe. Therefore, vaccination should be necessary as early as possible to protect the world from the Covid-19 endemic. Our mathematical model on the Covid-19 epi- demic diseases gives some consequences of pub- lic health policies. Many of our proposed model Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 16 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... parameters are assumed or estimated, but it de- pends on many factors; these parameters may be considered as fuzzy or stochastic rather than deterministic. Consequently, it may include fuzzy or stochastic differential equations in the proposed model for future work consideration. The progress of treating Covid-19 disease by different medicines in a cost-effective way is the main objective of health administrators, policy-makers, and scientists until a vaccine is discovered. The present paper gives a little effort to reach this objective to restrict Covid-19 infection. Acknowledgments: The authors would like to express their gratitude to the Editor Hristo V Kojouharov and Referees for their encouragement and constructive comments in revising the paper. REFERENCES [1] S. E. Adler, Why Coronaviruses Hit Older Adults Hard- est, AARP2020. [2] J. Cohen, D. Normile, New SARS-like virus in China triggers alarm, Science 367 (6475) (2020) 234–235. [3] P. Wang, X. Zheng, J. Li, B. Zhu, Prediction of epidemic trends in COVID-19 with logistic model and machine learning technics, Chaos Soliton Fract., 139 (2020), 110058. [4] C. Pai, A. Bhaskar, V. Rawoot, Investigating the dynam- ics of COVID-19 pandemic in India under lockdown, Chaos Soliton Fract., 138 (2020), 109988. [5] S. Djilali, B. Ghanbari, Coronavirus pandemic: A pre- dictive analysis of the peak outbreak epidemic in South Africa, Turkey, and Brazi, Chaos Soliton Fract., 138 (2020), 109971. [6] K. Ayinde, A. F. Lukman, R. I. Rauf, O. O. Alabi, C. E. Okon, O. E. Ayinde, Modeling Nigerian Covid-19 cases: A comparative analysis of models and estimators, Chaos Soliton Fract., 138 (2020), 109911. [7] M. Yousaf, S. Zahir, M. Riaz, S. M. Hussain, K. Shah, Statistical analysis of forecasting COVID-19 for upcoming month in Pakistan, Chaos Soliton Fract., 138 (2020), 109926. [8] I. Kırbaş, A. Sözen, A. D. Tuncer, F. Ş. Kazancıoğlu, Comparative analysis and forecasting of COVID-19 cases in various European countries with ARIMA, NARNN and LSTM approaches, Chaos Soliton Fract., 138 (2020), 110015. [9] W. E. Alnaser, M. Abdel-Aty, O. Al-Ubaydli, Mathe- matical prospective of coronavirus infections in bahrain, saudi arabia and egypt, Inf. Sci. Lett., 09 (2020), 51-64. [10] C. M. Păcurar, B. R. Necula, An analysis of COVID- 19 spread based on fractal interpolation and fractal dimension, Chaos Soliton Fract., 139 (2020), 110073. [11] P. Arora, H. Kumar, B. K. Panigrahi, Prediction and analysis of COVID-19 positive cases using deep learn- ing models: A descriptive case study of India, Chaos Soliton Fract., 139 (2020), 110017. [12] India covid-19 tracker. https://www.covid19india.org/, 2020. [Retrieved: 08/05/2020] [13] Q. Pan, T. Gao, M. He, Influence of isolation measures for patients with mild symptoms on the spread of COVID-19, Chaos Soliton Fract., 139 (2020), 110022. [14] B. K. Sahoo, B. K. Sapra, A data driven epidemic model to analyse the lockdown effect and predict the course of COVID-19 progress in India, Chaos Soliton Fract., 139 (2020), 110034. [15] H. Panwar, P. K. Gupta, M. K. Siddiqui, R. Morales- Menendez, V. Singh, Application of deep learning for fast detection of COVID-19 in X-Rays using nCOVnet, Chaos Soliton Fract., 138 (2020), 109944. [16] Y. Huang, Y. Wu, W. Zhang, Comprehensive identifica- tion and isolation policies have effectively suppressed the spread of COVID-19, Chaos Soliton Fract., 139 (2020), 110041. [17] P. C. L. Silva, P. V. C. Batista, H. S. Lima, M. A. Alves, F. G. Guimarães, R. C. P. Silva, An agent-based model of COVID-19 epidemic to simulate health and economic effects of social distancing interventions, Chaos Soliton Fract., 139 (2020), 110088. [18] B. K. Mishra, A. K. Keshri, Y. S. Rao, B. K. Mishra, B. Mahato, S. Ayesha, B. P. Rukhaiyyar, D. K. Saini, A. K. Singh, COVID-19 created chaos across the globe: Three novel quarantine epidemic models, Chaos Soliton Fract., 138 (2020), 109928. [19] V. M. Marquioni, M. A. M. de Aguiar, Quantifying the effects of quarantine using an IBM SEIR model on scalefree networks, Chaos Soliton Fract., 138 (2020), 109999. [20] H. B. Fredj, F. Chérif, Novel Corona virus disease infection in Tunisia: Mathematical model and the impact of the quarantine strategy, Chaos Soliton Fract., 138 (2020), 109969. [21] N. Ferguson, D. Laydon, G. N. Gilani, N. Imai, K. Ainslie, M. Baguelin, S. Bhatia, A. Boonyasiri, Z. C. Perez, G. Cuomo-Dannenburg, et al. Report 9: Impact of non-pharmaceutical interventions (npis) to reduce covid19 mortality and healthcare demand. 2020. [22] B. Tang, N. L. Bragazzi, Q. Li, S. Tang, Y. Xiao, J. Wu, An updated estimation of the risk of transmission of the novel coronavirus (2019-ncov). Infectious Disease Modelling, 5 (2020) 248-255. [23] R. Anguelov, J. Banasiak, C. Bright, J. Lubuma, R. Ouifki, The big unknown: The asymptomatic spread of COVID-19, Biomath. 9 (2020), 2005103. [24] A. Behnood, E. Mohammadi Golafshani, S. M. Hos- seini, Determinants of the infection rate of the COVID- 19 in the U.S. using ANFIS and virus optimization algorithm (VOA), Chaos Soliton Fract., 139 (2020), 110051. [25] R. G. da Silva, M. H. D. M. Ribeiro , V. C. Mariani, Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 17 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... L. D. S. Coelho, Forecasting Brazilian and American COVID-19 cases based on artificial intelligence coupled with climatic exogenous variables, Chaos Soliton Fract., 139 (2020), 110027. [26] P. C. L. Silva, P. V. C Batista, H. S. Lima, M. A. Alves, F. G. Guimarães, R. C. P. Silva, An agent-based model of COVID-19 epidemic to simulate health and economic effects of social distancing interventions, Chaos Soliton Fract., 139 (2020), 110088. [27] D. Okuonghae, A. Omame, Analysis of a mathematical model for COVID-19 population dynamics in Lagos, Nigeria, Chaos Soliton Fract., 139 (2020), 110032. [28] B. J. Quilty, S. Cli ord, et al. Effectiveness of air- port screening at detecting travellers infected with novel coronavirus (2019-ncov), Eurosurveillance, 25(5) (2020). [29] M. Shen, Z. Peng, Y. Xiao, L. Zhang, Modelling the epidemic trend of the 2019 novel coronavirus outbreak in china. bioRxiv, 2020. [30] B. Tang, X. Wang, Q. Li, N. L. Bragazzi, S. Tang, Y. Xiao, J. Wu, Estimation of the transmission risk of the 2019-ncov and its implication for public health interventions. Journal of Clinical Medicine, 9(2) (2020) 462. [31] V. Soukhovolsky, A. Kovalev, A. Pitt, B. Kessel, A new modelling of the COVID 19 pandemic, Chaos Soliton Fract., 139 (2020), 110039. [32] Sk S. Nadim, I. Ghosh, J. Chattopadhyay, Short-term predictions and prevention strategies for covid-2019: A model based study. arXiv preprint arXiv:2003.08150, 2020. [33] M. Cadoni, G. Gaeta, Size and timescale of epidemics in the SIR framework, Physica D, 411 (2020), 132626. [34] Y. Zhang, X. Yu, H. Sun, G. R. Tick, W. Wei, B. Jin, Applicability of time fractional derivative models for simulating the dynamics and mitigation scenarios of COVID-19, Chaos Soliton Fract., 138 (2020), 109959. [35] M. Higazy, Novel fractional order SIDARTHE mathe- matical model of COVID-19 pandemic, Chaos Soliton Fract., 138 (2020), 110007. [36] J. T. Wu, K. Leung, G. M Leung, Nowcasting and fore- casting the potential domestic and international spread of the 2019-ncov outbreak originating in wuhan, china: a modelling study. The Lancet, 395(10225) (2020) 689- 697. [37] J. M. Read, J. RE Bridgen, D. AT Cummings, A. Ho, C. P Jewell, Novel coronavirus 2019-ncov: early estimation of epidemiological parameters and epidemic predictions, MedRxiv, 2020. [38] A. Paul, S. Chatterjee, N. Bairagi, Prediction on Covid- 19 epidemic for different countries: Focusing on South Asia under various precautionary measures, medRxiv, 2020. [39] T. Sardar, Sk S. Nadimb, J. Chattopadhyayb, As- sessment of 21 days lockdown effect in some states and overall India: A predictive mathemati- cal study on COVID-19 outbreak, arXiv preprint arXiv:2004.03487V1, 2020. [40] D. Pal, D. Ghosh, P.K. Santra, G.S. Mahapatra, Mathe- matical analysis of a COVID-19 epidemic model by us- ing data driven epidemiological parameters of diseases spread in India, medRxiv, 2020. [41] P. Van den Driessche, J. Watmough, Reproduction num- bers and sub-threshold endemic equilibria for compart- mental models of disease transmission, Math Biosci 180 (2002) 29–48. [42] S. Çakan, Dynamic analysis of a mathematical model with health care capacity for COVID-19 pandemic, Chaos Soliton Fract., 139 (2020), 110033. [43] K. Y. Ng, M. M. Gui, COVID-19: Development of a robust mathematical model and simulation package with consideration for ageing population and time delay for control action and resusceptibility, Physica D, 411 (2020), 132599. [44] N. Al-Asuoad, S. Alaswad, L. Rong, M. Shillor, Mathe- matical model and simulations of MERS outbreak: Pre- dictions and implications for control measures, Biomath 5 (2016), 1612141 [45] H. Zhao, Z. Feng, Staggered release policies for COVID-19 control: Costs and benefits of relaxing re- strictions by age and risk, Chaos Soliton Fract., 138 (2020), 108405. [46] Z. Abbasi, I. Zamani, A. H. A. Mehra, M. Shafieirad, A. Ibeas, Optimal Control Design of Impulsive SQEIAR Epidemic Models with Application to COVID-19, Chaos Soliton Fract., 139 (2020), 110054. [47] V. Mbazumutima, C. Thron, L. Todjihounde, E. nu- merical solution for optimal control using treatment and vaccination for an SIS epidemic model, Biomath 8 (2019), 1912137. [48] S. Lenhart, J. T. Workman, Optimal control applied to biological models. Mathematical and computational biology. Boca Raton (Fla.), London: Chapman & Hall/CRC, 2007. [49] A. C. Chiang, Elements of dynamic optimization. New York, NY: McGraw-Hill international editions, McGraw-Hill. [u.a.], internat. ed. edition, 1992. [50] R. M. Anderson, R. M. May, Population biology of infectious diseases. Part I, Nature, 280 (1979) 361–367. [51] F. Brauer, C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology (Springer, Berlin, 2001). [52] A. Vivas-Barber, C. Castillo-Chavez, E. Barany, Dynam- ics of an “SAIQR” Influenza Model, Biomath 3 (2014), 1409251 [53] F. Nyabadza, S. D. Hove-Musekwa, From heroin epi- demics to methamphetamine epidemics: Modeling sub- stance abuse in a South African province, Math. Biosci. 225 (2010) 132–140. [54] P. van den Driessche, J. Watmough, Reproduction num- bers and sub-threshold endemic equilibria for compart- mental models of disease transmission, Math. Biosci. 180 (2002) 29–48. Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 18 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... [55] Carlos Castillo-Chavez, Zhilan Feng, Wenzhang Huang, On the computation of R0 and its role on Global Stability, Mathematical approaches for emerging and reemerging infectious diseases: an introduction, 1:229, 2002. [56] K. Blayneh, Y. Cao, H. D. Kwon, Optimal control of vectorborne disease: treatment and prevention, Discrete Continuous Dyn. Syst. Ser. B 11 (2009) 1–31. [57] C. Castillo-Chevez, Z. Feng, Global stability of an agestructure model for TB and its applications to opti- mal vaccination strategies, Math Biosci 151(1998) 135– 154. [58] D. T. Fleming, J. N. Wasserheit, From epidemiological synergy to public health policy and practice: the contri- bution of other sexually transmitted diseases to sexual transmission of HIV infection, Sex Transm. Infect. 75(4) (1999) 3–17. [59] H. R. Joshi, Optimal control of an HIV immunology model, Optim. Control. Appl. Methods 23 (2002) 199– 213. [60] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, E. F. Mishchenko, The mathematical theory of optimal processes. Gordon and Breach Science, London, 1986. [61] knoema, https://knoema.com/atlas/India/Birth-rate, Re- trieved : 2020-04-08. [62] M. Kot, Elements of Mathematical Ecology, Cambridge University Press, Cambridge, 2001. Appendix A. Proof of Theorem 4 The variational matrix of system (3) at DFE E0 is given by ME0 =   −d1 0 −α1Λd1 − α2Λ d1 0 0 0 −A α1Λ d1 α2Λ d1 0 0 0 β1 −B 0 0 0 0 β3 0 −C 0 0 0 β2 γ2 σ2 −D 0 0 0 γ1 σ1 ∈ −d1   Therefore, eigenvalues of the characteristic equa- tion of ME0 are −d1,−d1 and −D and and the solution of the cubic equation, P (λ) ≡ λ3 + A1λ2 + A2λ + A3 = 0 (A.1) where A1 = (A + B + C), A2 = (AB + AC)(1−R0)+BC + ( α1β1C B + α2β3B C ) Λ d1 , A3 =ABC(1 −R0). Now, it is easily noted that, A1 > 0, A3 > 0 if R0 < 1. After some simplifications, we get A1A2−A3 = (A+B)AB [ (1−R0)+ Λα2β3 d1AC ] +(A + C)AC [ (1 −R0) + Λα1β1 d1AB ] +BC(B + C) + 2ABC Here, we can notice that, if R0 < 1 then A1A2 − A3 > 0 if R0 < 1. Therefore, by the Routh– Hurwitz Routh–Hurwitz criterion [62] it follows that P (λ) = 0 has negative real roots if R0 < 1, i.e., the system (3) at DFE E0 when R0 < 1.This completes the proof. Appendix B. Proof of Theorem 6 The variational matrix of system (3) at E1(S ∗,E∗,H∗,G∗,I∗,R∗) is given by, ME1 =   b11 0 b13 b14 0 0 b21 b22 b23 b24 0 0 0 b32 b33 0 0 0 0 b42 0 b44 0 0 0 b52 b53 b54 b55 0 0 0 b63 b64 b65 b66   where, b11 = −d1R0, b13 = −α1S∗, b14 = −α2S∗, b21 = d1(R0−1), b22 = −A, b23 = α1S∗, b24 = α2S ∗, b32 = β1, b33 = −B, b42 = β3, b44 = −C, b52 = β2, b53 = γ2, b54 = σ2, b55 = −D, b63 = γ1, b64 = σ1, b65 =∈, b66 = −d1. Therefore, eigenvalues of the characteristic equation of ME1 are −D, −d1and the solution of the equation, Q (λ) ≡ λ4 +B1λ3 +B2λ2 +B3λ+B4 = 0 (B.1) where B1 = A + B + C + d1R0, B2 = Λ d1R0 [ α1β1C B + α2β3B C ] +(A + B + C)d1R0 + BC, B3 = (α1β1 + α2β3)d1(R0 − 1) + BCd1R0 + Λα1β1C B + Λα2β3B C B4 = ABCd1(R0 − 1). Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 19 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 D Pal, D Ghosh, P K Santra, G S Mahapatra, Mathematical modeling and analysis of Covid-19 ... Now, it is easily noted that Bi > 0 (i = 1, 2, 3) and B4 > 0 if R0 > 1. By the Routh–Hurwitz criterion [62], it follows that Q (λ) = 0 has negative real roots if Bi > 0 for i = 1, 2, 3, 4, D1 = B1 > 0, D2 = ∣∣∣∣B1 B31 B2 ∣∣∣∣ = B1B2 −B3 > 0, D3 = ∣∣∣∣∣∣ B1 B3 0 1 B2 B4 0 B1 B3 ∣∣∣∣∣∣ = B1B2B3 −B21B4 −B 2 3 > 0. Therefore the system (3) shows local asymptotic stability at E1 when R0 > 1, B1B2 −B3 > 0 and B1B2B3 − B21B4 − B 2 3 > 0. This completes the proof. Biomath 10 (2021), 2106147, http://dx.doi.org/10.11145/j.biomath.2021.06.147 Page 20 of 20 http://dx.doi.org/10.11145/j.biomath.2021.06.147 Introduction Multi-objective optimal control Derivation and Preliminaries of Covid-19 Model Fundamental Properties Positivity of the solutions Invariant region Existence of Equilibrium and Stability Analysis The basic reproduction number Existence of endemic equilibrium E1(S,E,H,G,I,R) Asymptotic behavior Proposed Covid-19 Model with Control Numerical Simulations Discussions and Conclusions References