# Dynamics of COVID-19 Using SEIQR Epidemic Model

Authors: N. Avinash; G. Britto Antony Xavier; Ammar Alsinai; Hanan Ahmed; V. Rexma Sherine; P. Chellamani
Journal: Journal of Mathematics
 (2022)

---

## Abstract

The major goal of this study is to create an optimal technique for managing COVID-19 spread by transforming the SEIQR model into a dynamic (multistage) programming problem with continuous and discrete time-varying transmission rates as optimizing variables. We have developed an optimal control problem for a discrete-time, deterministic susceptible class (S), exposed class (E), infected class (I), quarantined class (Q), and recovered class (R) epidemic with a finite time horizon. The problem involves finding the minimum objective function of a controlled process subject to the constraints of limited resources. For our model, we present a new technique based on dynamic programming problem solutions that can be used to minimize infection rate and maximize recovery rate. We developed suitable conditions for obtaining monotonic solutions and proposed a dynamic programming model to obtain optimal transmission rate sequences. We explored the positivity and unique solvability nature of these implicit and explicit time-discrete models. According to our findings, isolating the affected humans can limit the danger of COVID-19 spreading in the future.


---

## Body

## 1. Introduction

Mathematical models are useful for determining how an infection behaves when it enters a population and determining whether it will be eradicated or continue under different settings. COVID-19 is currently causing tremendous concern among researchers, governments, and the general public due to its rapid spread and a high number of deaths [1]. The transmission of this disease is caused by the tiny particles or droplets called aerosols that carry the virus into the atmosphere caused by a contaminated person while sneezing, coughing, or exhaling. Many researchers and scientists are continuously working to reduce the transmission of this vicious disease throughout the world. Infectious diseases are the disciplines that focus on the study of the dynamics of infectious diseases as well as the relationship between these diseases and the various factors involved in their appearance and evolution, in order to implement a fight against this spread. Despite its youth, mathematical modeling is a valuable tool for understanding disease transmission mechanisms, is playing an increasingly important role in epidemiology, and has already contributed to significant successes. The most influential work in the field of mathematics epidemiology was first introduced by Kermack and McKendrick as the SIR model in the year 1927 [2]. Cao et al.[3] discussed a modified model of the SIR (susceptible, infected, recovered) epidemic introduced in order to detect the confirmed number of infected cases and consecutive burdens on isolation wards and ICUs. Also, Nesteruk [4] developed the variables used in the proposed model by introducing a SIR epidemic model and explaining how to dominate the spread of the disease. To restore the pandemic with the involvement of social distancing and lockdown, Gerberry and Milner presented a data-driven susceptible, exposed, infected, quarantined, and recovered (SEIQR) model in [5]. From the publication of Zeb et al. [6], epidemiological’s SEIQR model with isolation class in 2020 and their mathematical epidemiology has expanded in numerous directions, involving biology and computer science by Zima et al.[7], Zhou et al. [8], and Kermack and McKendrick [9, 10]. Some recent studies have focused on this area of research by He et al. [11], Rahimi et al. [12], Hussain et al. [13], Youssef et al. [14], Prabakaran et al. [15], and Youssef et al. [16]. In this current paper, we implement the discrete type of SEIQR model and discuss the solvability of both continuous and discrete type SEIQR model. We examined the behaviour of the time-continuous model. We have developed two time-discrete models: time-implicit and time-explicit. We looked at the theory and methods for solving the time-implicit model. Then, to control and anticipate the dynamics of COVID dissemination, we establish appropriate transmission rate limits. To do this, we devised a dynamic programming problem to optimize transmission rate sequences under arbitrary beginning conditions. We propose safety guidelines and essential precautionary measures based on the optimized rate sequences to control COVID spread. The article gives the technique for optimizing the transmission rate sequences. The epidemic models and their time-discrete variations have been studied by Allen [17] and Ghosh et al. [18]. Several approaches towards fractional-order mathematical models of COVID-19 were studied by the authors Alqhtani et al. [19], Valliammal and Ravichandran [20], Nisar et al. [21], Vijayakumar et al. [22], and Alderremy et al. [23]. However, the aforementioned studies and references mostly contain explicit approaches with respect to time-discrete epidemic models.


### 1.1. Review Literature

In 2020, COVID-19 is a worldwide emergency. The first cases occurred in December 2019, and as of 6 : 34 pm CEST, 28 July 2022, there have been 571,198,904 confirmed cases of COVID-19, including 6,387,863 deaths, reported to WHO. As of 25 July 2022, a total of 12,248,795,623 vaccine doses have been administered. The rapid spread of COVID-19 has already caused great public attention and many heated discussions, and the Chinese mass media have been reporting relevant information about the virus and the outbreak.Ming et al. [24] show that effective public health measures are required to be implemented in time to avoid the breakdown of the health system, and the media can certainly play a crucial role in conveying updated policies and regulations from authorities to the citizens. The finding that SARS-2-S exploits ACE2 for entry, which was also reported by Kermack and McKendrick [9] while the present manuscript was in revision, suggests that the virus might target a similar spectrum of cells as SARS-CoV. However, upon its outbreak, various research, including but not limited to Okhuese [25], began to predict the scale that the virus would hit the world; the ratio of the death to recovery rate has seemingly been a positive proportion. Allen [17] studied about time-discrete SI, SIR, and SIS epidemic models, and its properties. Kermack and McKendrick [10] analyzed an outbreak such as the one in Hubei is captured by SIR dynamics where the population is divided into three compartments that differentiate the state of individuals with respect to the contagion process: infected (I), susceptible (S) to infection, and removed (R) (i.e., not taking part in the transmission process). Mathematical modeling has been influential in providing a deeper understanding on the transmission mechanisms and burden of the ongoing COVID-19 pandemic, contributing to the development of public health policy and understanding. Most mathematical models of the COVID-19 pandemic can broadly be divided into either population-based, SIR (Kermack-McKendrick)-type models, driven by (potentially stochastic) differential equations proposed by Nesteruk [4] in which individuals typically interact on a network structure and exchange infection stochastically. This point emerges also clearly from a number of recent model-based contributions that have extended the basic SIR model to account for key insights from economic theory, namely by allowing for peoples’ (rational) adjustment of work, consumption, and leisure activities in the face of infection risk. More generally, the idea is to model explicitly the exposure to the virus (of those people who are susceptible), as in the susceptible-exposed-infectious-recovered (SEIR) model which has been analyzed extensively by He et al. [11] in the context of the COVID-19 pandemic. The Jacobian method used for the SEIR model yields a biologically reasonable R0, but for more complex compartmental models, especially those with more infected compartments, the method is hard to apply as it relies on the algebraic Routh Hurwitz conditions for stability of the Jacobian matrix. An alternative method proposed by the authors Van den Driessche and Watmough [26] gives a way of determining R0 for an ODE compartmental model by using the next generation matrix. Batista [27] applies logistic growth regression model to predict the final size of the Covid-19 epidemic. Basically, NSFD is an iterative method in which we get closer to solution through iteration was given by Mickens [28]. The authors Vijayakumar et al. [22] discussed about approximate controllability results for fractional Sobolev type Volterra-Fredholm integro-differential systems of order 1<r<2. Finding the variants that predict severe disease, we developed a collaboration of four international computational centers (Iran, Italy, Malaysia, and Greece). In [29], the authors Bairagi et al. have introduced a mathematical model for controlling the outbreak of COVID-19 by augmenting isolation and social distancing features of individuals and also solved the utility maximization problem by using a nonco-operative game. In 2021, the multidisciplinary approach was necessary to address the multidimensional aspects of COVID-19 infection by established collaborations discussed by the authors Rahimi et al. [12]. The study by the authors Prabakaran et al. [15] looked into the evolving geographic diversity of the SARS-CoV-2. We then consider how positive factors like social distancing measures and detrimental factors such as delays in testing onset affect optimal testing strategies and outbreak controllability. Throughout, Youssef et al. [14] focus their analyses on empirically supported parameter values including realistic testing rates. While many existing COVID-19 SIR-like compartmental models explore the effects of testing with forms of isolation like quarantine or hospitalization, the majority of these studies assume simple linear equations for the rates at which tests are administered and individuals are isolated. The authors Hussain et al. [13] discussed about the complex systems and network science approaches, along with technological advances and data availability, are becoming instrumental for the design of effective containment strategies. In a nonsense region, Hilfer’s neutral fractional derivative provided controllability results using Monch’s method, Banach’s contraction principle, fractional calculus, and semigroup property was studied by the authors Nisar et al. [21]. Some recent updates regarding the modeling of the coronavirus, the authors Alderremy et al. [23] constructed a mathematical model based on the fuzzy fractional derivative and obtained the results. The authors Valliammal and Ravichandran [20] are discussed in detail the fractional integro-differential equation with different conditions in various spaces. In [30], the authors Awal et al. proposed a framework that uses Bayesian optimization to optimize the hyperparameters of the classifier and adaptive synthetic (ADASYN) algorithm to balance the COVID and non-COVID classes of the dataset. In 2022, the authors Youssef et al. considered a modified model to analyze the disease dynamics of the coronavirus infection by taking the real cases from Saudi Arabia [16]. Alqhtani et al. [19] analyzed about spatiotemporal dynamical patterns arising from subdiffusion reaction-diffusion systems of predator-prey interaction are modeled in the sense of the Caputo fractional operator. Ghosh et al. [18] have studied about discrete-time epidemic model for the analysis of transmission of COVID-19 based upon data of epidemiological parameters.In this article, we have considered the epidemic model published in [6, 30]. Then, we have extended the idea of the article [1] to the considered model. As a result, we recap and extend certain conclusions on the features of the time-continuous classical SEIQR model, and we suggest an implicit time-discrete version of this classical SEIQR model, proving that it retains many of the qualities of the time-continuous version. As a result, the goal of this research is to propose a nonautonomous SEIQR model, investigate the properties of its time-continuous formulation, and design an implicit numerical solution approach that preserves the time-continuous variant’s primary properties. The goal of this article is to propose, analyze, and optimize COVID-19 using the SEIQR epidemic model. According to our investigations, COVID-19 outbreaks might be caused by human-to-human interaction. As a result, isolation of the infected humans can reduce the COVID-19 spread in the future. Literature review and comparison of various of these models are presented in Table 1.Table 1
Overview of models studied.
ModelsCompartmentsModel typeStudy areaSEIR/SLIRSusceptible (S), exposed/latent (E/L), infectious (I), removed (R)DeterministicEurope and North America, New York, Mexico, Zhejiang, Guangdong, Japan, India [31–36]SEIQRSusceptible (S), exposed (E), hospitalized infected (I), quarantine (Q), recovered or removed (R)DeterministicIndia [37]SIR-XInfected (I), susceptible (S), removed (R), quarantined (X)DeterministicChina [38]SIRDSusceptible (S), infected (I), recovered (R), dead (D)DeterministicChina, Italy, and France [39]SEIHARDSusceptible (S), exposed (E), symptomatic infectious (I), hospitalized (H), asymptomatic infectious (A), recovered (R), deaths (D)DeterministicWashington, New York [40]SIRUSusceptible (S), asymptomatic infectious (I), reported symptomatic infectious (R), unreported symptomatic infectious (U)DeterministicChina, Hubei, Wuhan [41]SEIPAHRFSusceptible (S), exposed (E), symptomatic (I), super-spreaders class (P), asymptomatic infectious (A), hospitalized (H), recovery (R), fatality (F)DeterministicWuhan [42]SEIRUSusceptible (S), asymptomatic noninfectious (E), asymptomatic infectious (I), reported symptomatic infectious (R), unreported symptomatic infectious (U)DeterministicChina [43]SEIHRSusceptible (S), exposed (E), symptomatic infectious (I), hospitalized (H), recovered or death (R)DeterministicSouth Korea [44]SEIRPSusceptible (S), exposed (E), infectious (I), removed (R), pathogens (P)DeterministicPakistan [45]More precisely, our main contributions can be summarized as follows:(i)
First, we suggest a time-continuous SEIQR model modification with time-varying transmission and recovery rates.(ii)
Second, we draw the conclusion that the formulation of our time-continuous problem is well-posed. This comprises continuous reliance on initial conditions and time-varying rates, global existence in time, and global uniqueness in time, all of which are based on an inductive application of Banach’s fixed point theorem.(iii)
In the case of the time-discrete implicit model, we provide unique solvability, monotonicity properties, and an upper error bound between the solution of the implicit time-discrete problem formulation and the solution of the time-continuous problem formulation.(iv)
In order to maximize transmission rate sequences under arbitrary beginning conditions, we have developed a dynamic programming problem. Based on the optimal rate sequences, we suggest safety guidelines and important safety precautions to control COVID spread.The paper is arranged as follows: Section1 is dedicated to the introduction. In section 2, we present the time-continuous and time-discrete SEIQR model. In section 3, we give the monotonicity properties and long-time behaviour. An error analysis is given in section 4. The conclusion of our research work is implemented in the last section 5.


## 1.1. Review Literature

In 2020, COVID-19 is a worldwide emergency. The first cases occurred in December 2019, and as of 6 : 34 pm CEST, 28 July 2022, there have been 571,198,904 confirmed cases of COVID-19, including 6,387,863 deaths, reported to WHO. As of 25 July 2022, a total of 12,248,795,623 vaccine doses have been administered. The rapid spread of COVID-19 has already caused great public attention and many heated discussions, and the Chinese mass media have been reporting relevant information about the virus and the outbreak.Ming et al. [24] show that effective public health measures are required to be implemented in time to avoid the breakdown of the health system, and the media can certainly play a crucial role in conveying updated policies and regulations from authorities to the citizens. The finding that SARS-2-S exploits ACE2 for entry, which was also reported by Kermack and McKendrick [9] while the present manuscript was in revision, suggests that the virus might target a similar spectrum of cells as SARS-CoV. However, upon its outbreak, various research, including but not limited to Okhuese [25], began to predict the scale that the virus would hit the world; the ratio of the death to recovery rate has seemingly been a positive proportion. Allen [17] studied about time-discrete SI, SIR, and SIS epidemic models, and its properties. Kermack and McKendrick [10] analyzed an outbreak such as the one in Hubei is captured by SIR dynamics where the population is divided into three compartments that differentiate the state of individuals with respect to the contagion process: infected (I), susceptible (S) to infection, and removed (R) (i.e., not taking part in the transmission process). Mathematical modeling has been influential in providing a deeper understanding on the transmission mechanisms and burden of the ongoing COVID-19 pandemic, contributing to the development of public health policy and understanding. Most mathematical models of the COVID-19 pandemic can broadly be divided into either population-based, SIR (Kermack-McKendrick)-type models, driven by (potentially stochastic) differential equations proposed by Nesteruk [4] in which individuals typically interact on a network structure and exchange infection stochastically. This point emerges also clearly from a number of recent model-based contributions that have extended the basic SIR model to account for key insights from economic theory, namely by allowing for peoples’ (rational) adjustment of work, consumption, and leisure activities in the face of infection risk. More generally, the idea is to model explicitly the exposure to the virus (of those people who are susceptible), as in the susceptible-exposed-infectious-recovered (SEIR) model which has been analyzed extensively by He et al. [11] in the context of the COVID-19 pandemic. The Jacobian method used for the SEIR model yields a biologically reasonable R0, but for more complex compartmental models, especially those with more infected compartments, the method is hard to apply as it relies on the algebraic Routh Hurwitz conditions for stability of the Jacobian matrix. An alternative method proposed by the authors Van den Driessche and Watmough [26] gives a way of determining R0 for an ODE compartmental model by using the next generation matrix. Batista [27] applies logistic growth regression model to predict the final size of the Covid-19 epidemic. Basically, NSFD is an iterative method in which we get closer to solution through iteration was given by Mickens [28]. The authors Vijayakumar et al. [22] discussed about approximate controllability results for fractional Sobolev type Volterra-Fredholm integro-differential systems of order 1<r<2. Finding the variants that predict severe disease, we developed a collaboration of four international computational centers (Iran, Italy, Malaysia, and Greece). In [29], the authors Bairagi et al. have introduced a mathematical model for controlling the outbreak of COVID-19 by augmenting isolation and social distancing features of individuals and also solved the utility maximization problem by using a nonco-operative game. In 2021, the multidisciplinary approach was necessary to address the multidimensional aspects of COVID-19 infection by established collaborations discussed by the authors Rahimi et al. [12]. The study by the authors Prabakaran et al. [15] looked into the evolving geographic diversity of the SARS-CoV-2. We then consider how positive factors like social distancing measures and detrimental factors such as delays in testing onset affect optimal testing strategies and outbreak controllability. Throughout, Youssef et al. [14] focus their analyses on empirically supported parameter values including realistic testing rates. While many existing COVID-19 SIR-like compartmental models explore the effects of testing with forms of isolation like quarantine or hospitalization, the majority of these studies assume simple linear equations for the rates at which tests are administered and individuals are isolated. The authors Hussain et al. [13] discussed about the complex systems and network science approaches, along with technological advances and data availability, are becoming instrumental for the design of effective containment strategies. In a nonsense region, Hilfer’s neutral fractional derivative provided controllability results using Monch’s method, Banach’s contraction principle, fractional calculus, and semigroup property was studied by the authors Nisar et al. [21]. Some recent updates regarding the modeling of the coronavirus, the authors Alderremy et al. [23] constructed a mathematical model based on the fuzzy fractional derivative and obtained the results. The authors Valliammal and Ravichandran [20] are discussed in detail the fractional integro-differential equation with different conditions in various spaces. In [30], the authors Awal et al. proposed a framework that uses Bayesian optimization to optimize the hyperparameters of the classifier and adaptive synthetic (ADASYN) algorithm to balance the COVID and non-COVID classes of the dataset. In 2022, the authors Youssef et al. considered a modified model to analyze the disease dynamics of the coronavirus infection by taking the real cases from Saudi Arabia [16]. Alqhtani et al. [19] analyzed about spatiotemporal dynamical patterns arising from subdiffusion reaction-diffusion systems of predator-prey interaction are modeled in the sense of the Caputo fractional operator. Ghosh et al. [18] have studied about discrete-time epidemic model for the analysis of transmission of COVID-19 based upon data of epidemiological parameters.In this article, we have considered the epidemic model published in [6, 30]. Then, we have extended the idea of the article [1] to the considered model. As a result, we recap and extend certain conclusions on the features of the time-continuous classical SEIQR model, and we suggest an implicit time-discrete version of this classical SEIQR model, proving that it retains many of the qualities of the time-continuous version. As a result, the goal of this research is to propose a nonautonomous SEIQR model, investigate the properties of its time-continuous formulation, and design an implicit numerical solution approach that preserves the time-continuous variant’s primary properties. The goal of this article is to propose, analyze, and optimize COVID-19 using the SEIQR epidemic model. According to our investigations, COVID-19 outbreaks might be caused by human-to-human interaction. As a result, isolation of the infected humans can reduce the COVID-19 spread in the future. Literature review and comparison of various of these models are presented in Table 1.Table 1
Overview of models studied.
ModelsCompartmentsModel typeStudy areaSEIR/SLIRSusceptible (S), exposed/latent (E/L), infectious (I), removed (R)DeterministicEurope and North America, New York, Mexico, Zhejiang, Guangdong, Japan, India [31–36]SEIQRSusceptible (S), exposed (E), hospitalized infected (I), quarantine (Q), recovered or removed (R)DeterministicIndia [37]SIR-XInfected (I), susceptible (S), removed (R), quarantined (X)DeterministicChina [38]SIRDSusceptible (S), infected (I), recovered (R), dead (D)DeterministicChina, Italy, and France [39]SEIHARDSusceptible (S), exposed (E), symptomatic infectious (I), hospitalized (H), asymptomatic infectious (A), recovered (R), deaths (D)DeterministicWashington, New York [40]SIRUSusceptible (S), asymptomatic infectious (I), reported symptomatic infectious (R), unreported symptomatic infectious (U)DeterministicChina, Hubei, Wuhan [41]SEIPAHRFSusceptible (S), exposed (E), symptomatic (I), super-spreaders class (P), asymptomatic infectious (A), hospitalized (H), recovery (R), fatality (F)DeterministicWuhan [42]SEIRUSusceptible (S), asymptomatic noninfectious (E), asymptomatic infectious (I), reported symptomatic infectious (R), unreported symptomatic infectious (U)DeterministicChina [43]SEIHRSusceptible (S), exposed (E), symptomatic infectious (I), hospitalized (H), recovered or death (R)DeterministicSouth Korea [44]SEIRPSusceptible (S), exposed (E), infectious (I), removed (R), pathogens (P)DeterministicPakistan [45]More precisely, our main contributions can be summarized as follows:(i)
First, we suggest a time-continuous SEIQR model modification with time-varying transmission and recovery rates.(ii)
Second, we draw the conclusion that the formulation of our time-continuous problem is well-posed. This comprises continuous reliance on initial conditions and time-varying rates, global existence in time, and global uniqueness in time, all of which are based on an inductive application of Banach’s fixed point theorem.(iii)
In the case of the time-discrete implicit model, we provide unique solvability, monotonicity properties, and an upper error bound between the solution of the implicit time-discrete problem formulation and the solution of the time-continuous problem formulation.(iv)
In order to maximize transmission rate sequences under arbitrary beginning conditions, we have developed a dynamic programming problem. Based on the optimal rate sequences, we suggest safety guidelines and important safety precautions to control COVID spread.The paper is arranged as follows: Section1 is dedicated to the introduction. In section 2, we present the time-continuous and time-discrete SEIQR model. In section 3, we give the monotonicity properties and long-time behaviour. An error analysis is given in section 4. The conclusion of our research work is implemented in the last section 5.


## 2. Time-Continuous SEIQR Model

The time-continuous SEIQR model is formulated, and its behaviours are described using the Lipchitz condition and Grownwall and Bellman’s inequality in this section.


### 2.1. Mathematical Background Material

Here, we revisit the Lipschitz continuity of a function on Euclidean spaces, the local Lipchitz condition, Banach’s fixed point theorem, and the method of variation of the parameter, which will be used in the subsequent sections.Definition 1.
(see [46]). Let q1 and q2 be two positive integers and D⊂ℝq1. A function H:D⟶ℝq2 is said to be Lipchitz continuous on D if there exists a nonnegative constant L≥0 such that Hx−Hyℝq2≤L⋅x−yℝq1 holds for x,y∈D.
(i) LetU⊂ℝq1 be an open set and H:U⟶ℝq2. Then, H is called as locally Lipchitz continuous if for every element y0∈U there exists a neighborhood V of y0 such that the restrictions of H to V are Lipchitz continuous on V.In a more general framework, we consider a nonlinear initial value problem (IVP) z′t = Ht,zt; z0 = z0, where zt = y1t,y2t,…,yntT is solution vector, and Ht,zt = h1t,zt,h2t,zt,…,hnt,ztT is vectorial function with initial point z0∈ℝn. The following theorem, which is a direct consequence of Gronwall’s lemma, can be used to prove global existence in time.Theorem 1 (see [17]).
LetG:ℝn⟶ℝn be locally Lipchitz conditions. If there exist nonnegative real constants B and k such that Gt,zt≤kzt∞+B holds for all zt∈ℝn, then the solution of the initial value problem z′t=Gt,zt,z0=z0, exists for all time t∈ℝ, and moreover, zt satisfies ztℝn≤z0ℝn⋅expk⋅t+B/kexpk⋅t−1,∀t∈ℝ.
We give the following: Banach’s fixed point theorem, which will be used to preserve the global uniqueness in time [47, 47].Theorem 2 (see [48]).
LetX,μ be a complete metric space. Let T:X⟶X be a strict contraction, that is, there exists a constant K∈0,1 such that μTx,Ty≤K⋅μx,y holds for all x,y∈X. Then, the mapping T has a unique fixed point.
In the following theorem, we present the Grownwall and Bellman inequality, which will be used in the subsequence theorems related to the continuous functions.Theorem 3 (see [49]).
LetI:=a,b, u,f:I⟶0,∞ be two continuous and nonnegative functions and g:I⟶0,∞ be a continuous, positive, and nondecreasing function. If the inequality(1)ut≤gt+∫atfs⋅usds,holds for all t∈I, we have(2)ut≤gt⋅exp∫atfs⋅usds.Theorem 4.
(Method of variation of parameter) For a first-order nonhomogeneous linear differential equation,y′t+ptyt=ft has the general solution yt=vtePt+AePt, where v′t=e−Ptft, Pt=−∫ptdt and A is an arbitrary constant.


### 2.2. Continuous Problem Formulation

At first, let us assume the following assumptions [50, 51] for the upcoming calculations.(i)
Let the population size varies over time beN is varying over time (i.e., population size = μNt for all t∈0,∞).(ii)
We divide the population into five homogeneous subgroups, namely susceptible people (S), exposed (E), infectious (I), quarantined (Q), and recovered (R). We can clearly assign every individual to exactly one subgroup. Hence, we obtain S,E,I,Q,R model satisfying the condition μN=St+Et+It+Qt+Rt for all t∈0,∞.(iii)
Each time-varying transmission ratex:0,∞⟶0,∞ is Lipchitz continuous and continuously differentiable, and there exist constants xmin and xmax such that 0<xmin≤xt≤xmax for all t≥0 and x∈π,β,γ,σ,θ,μ.The choice of time-dependent transmission rates is possible because the countermeasures such as lockdowns, social distancing, or other political actions like curfews and different medical treatments reduce possible contact between susceptible and infectious people.Our equations of the time-continuous SEIQR model read as follows:(3)dStdt=μN−μSt−βNStEt+It,dEtdt=βNStEt+It−πEt−μ+γEtdItdt=πEt−σIt−μIt,dQtdt=γEt+σIt−θQt−μQt,dRtdt=θQt−μRt,with initial conditions S0=S0>0,E0=E0>0,I0=I0>0,Q0=Q0>0 and R0=R0≥0. The detailed parameters and description are given in Table 2Table 2
Parameters and description.
ParametersDescriptionS(t)At timet, the number of susceptible peopleE(t)At timet, the number of exposed peopleI(t)At timet, the number of infected peopleQ(t)At timet, the number of quarantined peopleR(t)At timet, the number of recovered peopleβThe rate at which susceptible populations migrate to exposed and infected populationsπThe rate at which an exposed population moves to an infected populationγTransmission rate at which exposed people take outside as isolatedσTransmission rate at which infected people were added to isolated individualθTransmission rate at which isolated persons recoveredμNatural death rate and disease-related death rate


### 2.3. Nonnegativity and Boundedness of Solutions

Now, we prove the boundedness of the solution to (3). For this purpose, we modify ideas given in [51, 52] deriving the following lemmas, so consider the bounded, time-varying transmission rates given above.Lemma 1.
Each solution of system (3) is bounded below by zero.Proof.
Consider, the first relation of (3),(4)dStdt=μN−μSt−βNStEt+It.
By takingdSt/dt=S′t, equation (4) can be expressed as a first-order nonhomogeneous linear differential equation in St as(5)S′t+μ+βNEt+βNItSt=μN.
Applying Theorem4, and by applying the same procedure to the first-order nonhomogenous linear equation in Et,(6)dEtdt−βNSt−π−μ−γEt=βNStIt.
We can easily show thatEt≥0 for all t∈0,∞. Proceeding like this, we can show that It,Qt, and Rt≥0 for all t∈0,∞.
SinceμNt=St+Et+It+Qt+Rt, and total population is finite, St,Et,It,Qt, and Rt are bounded above, and hence, the proof is complete.Theorem 5.
For all solution functions of (3), we have 0≤X≤μN, where X∈St,Et,It,Qt,Rt.Proof.
The proof follows fromμNt=St+Et+It+Qt+Rt and Lemma 1.


### 2.4. Global Existence in Time

We arrive at a theorem regarding global existence in the time of (3) based on Theorem 1. For abbreviation, we use the supremum norm ft∞:=supt∈a,bft for an arbitrary continuous function f:a,b⟶ℝ. A similar definition holds for vector-valued bounded functions. In our case, using the boundedness of S,E,I,Q,R on 0,∞, obtain the following global existence theorem.Theorem 6.
The system of nonlinear first-order ODE (3) has at least one solution which exists for all t≥0.Proof.
By denotingzt=St,Et,It,Qt,Rt, we can set
G:0,∞×ℝ5⟶ℝ5 by(7)t,zt⟶A−μSt−βNStEt+ItβNStEt+It−πEt−μ+γEtπEt−σIt−μItγEt+σIt−θQt−μQtθQt−μRtT.
Clearly,G is Lipchitz continuous, due to the continuity of each components.
Assuming the supremum norm on our Euclidean space, and with the help of triangle inequality, we arrive(8)Gt,zt∞=supt∈0,∞A−μSt−βStEt+It,βNStEt+It−βEt−μ+γEt,πEt−σIt−μIt,γEt+σIt−θQt−μQt,θQt−μRt,Gt,zt∞≤supt∈0,∞A−μSt,βNStEt+It,πEt,γEt+σIt,θQt−μRt,≤supt∈0,∞μmaxN+μmaxSt,βmaxNStEt+It,πmaxEt,γmaxEt+σmaxIt,θmaxQt−μmaxRt,≤supt∈0,∞μmax+βmax+πmax+γmax+σmax+θmaxSt,Et,It,Qt,Rt,≤supt∈0,∞μmax+βmax+πmax+γmax+σmax+θmaxzt∞,Gt,zt∞≤kzt∞,where k=μmax+βmax+πmax+γmax+σmax+θmax.
From the boundedness of our solution functions and the boundedness of our time-varying transmission rates, all requirements of Theorem1 are fulfilled, and our proof is complete.


### 2.5. Global Uniqueness in Time

We present the global uniqueness theorem for (3) by utilizing the inductive application of Banach’s fixed point theorem.Theorem 7.
The nonlinear ODE system (3) has a unique solution that exists for all t≥0.Proof.
Consider the system of equations given in (3):(1)
Consider the time interval0,τ is applicable to Banach’s fixed point theorem(2)
Forx1,x2,y1,y2∈ℝ, by triangle inequality, we have x1⋅y1−x2⋅y2≤x1‖y1−y2+y2‖x1−x2(3)
We assume thatS,E,I,Q,R,S¯,E¯,I¯,Q¯,R¯:0,∞⟶0,∞ are two solutions of (3)
Beginning the proof by letting(9)supt∈0,τSt−S¯t=supS1ePt−S¯1ePtwherePt=−∫Ptdt,supt∈0,τSt−S¯t=supS1−S¯1‖ePt,supt∈0,τSt−S¯t≤S1−S¯1τzt−zt¯∞.
Then, the second equation in (3) becomes(10)dEtdt−βNSt−π−μ−γEt=βNStIt.
Since it is a first-order nonhomogenous linear equation inEt, we take p1t=−βNSt−π−μ−γ and ft=βNStIt. As P1t=−∫p1tdt and v1t=∫e−P1tβNStItdt. Hence, Et=v1tdteP1t+E1ePt and(11)supt∈0,τEt−E¯t≤E1−E¯1τzt−zt¯∞.
Similarly, we can easily show the following inequalities:(12)supt∈0,τQt−Q¯t≤Q1−Q¯1τzt−zt¯∞,supt∈0,τIt−I¯t≤I1−I¯1τzt−zt¯∞,supt∈0,τRt−R¯t≤R1−R¯1τzt−zt¯∞.
Summing implies(13)zt−zt¯∞≤∑S,E,I,Q,RS1−S¯1τzt−zt¯∞.
By choosingτ<1/21/∑S1−S¯1. As a result, we conclude that the solution is unique in 0,τ.


### 2.6. Time-Discrete Implicit SEIQR Model

Assume that our time interval0,τ can be divided by strictly increasing sequence tjj=1M for M∈ℕ with t1=0 and tM=T. We write ftj=fj, j∈1,2,3,…,M for an arbitrary time-dependent function f.


### 2.7. Discussion and Formulations

Here, we are transforming the continuous system (3) to the fully explicit discrete scheme (14) as given as follows:(14)ΔSjΔtj=μN−μj+1Sjt−βj+1NSjtEjt+Ijt,ΔEjΔtj=βj+1NSjtEj+Ij−πj+1Ej−μj+1+γj+1Ej,ΔIjΔtj=πj+1Ej−σj+1Ij−μj+1Ij,ΔQjΔtj=γj+1Ej+σj+1Ij−θj+1Qj−μj+1Qj,ΔRjΔtj=θj+1Qj−μj+1Rj,and a fully implicit scheme (15) as(15)ΔSjΔtj=μj+1N−μj+1Sj+1−βj+1NSj+1Ej+1+Ij+1t,ΔEjΔtj=βj+1NSj+1Ej+1+Ij+1−πj+1+μj+1+γj+1Ej+1,ΔIjΔtj=πj+1Ej+1−σj+1Ij+1−μj+1Ij+1,ΔQjΔtj=γj+1Ej+1+σj+1Ij+1−θj+1Qj+1−μj+1Qj+1,ΔRjΔtj=θj+1Qj+1−μj+1Rj+1,where ΔXj/Δtj = Xj+1−Xj/tj+1−tj for j∈1,2,3,…,M−1.Observe thatμj+1N=Sj+1+Ej+1+Ij+1+Qj+1+Rj+1 for all j∈1,2,3,…,M−1.Since the fully explicit scheme (14) simply reduces to a linear system, our main interest is in a fully implicit discrete scheme because it preserves the nonlinear structure of the continuous problem.


### 2.8. Implicit Time-Discrete Problem Formulation

In this section and subsequent sections, we derive the recurrence type of solutions for the implicit scheme (15). For that, we assume that 0<xmin≤xj≤xmax, x∈μ,β,π,γ,σ,θ, for j∈1,2,3,…,M and that 0<tj+1−tj≤1 for j∈1,2,3,…,M and that S0>0,E0>0,I0>0,Q0>0,R0≥0. Now, (15) can be expressed as(16)Sj+1=Sj+μj+1Ntj+1−tj1+μj+1+βj+1NEj+1+βj+1NIj+1tj+1−tj,Ej+1=Ej+βj+1NSj+1Ij+1tj+1−tj1+πj+1+μj+1+γj+1−βNSj+1tj+1−tj,Ij+1=Ij+πj+1Ej+1tj+1−tj1+σj+1+μj+1tj+1−tj,Qj+1=Qj+γj+1Ej+1+σj+1Ij+1tj+1−tj1+θj+1+μj+1tj+1−tj,Rj+1=Rj+θj+1Qj+1tj+1−tj1+μj+1tj+1−tj.


### 2.9. Unique Solvability

In this section, we provide the method for finding the solution of (16). Letting tj+1−tj=Δtj in Sj+1, we obtain(17)Sj+1=Sj+μj+1NΔtj1+μj+1+βj+1NEj+1+βj+1NIj+1Δtj,whenj=0, S1=S0+μ1NΔt0/1+μ1+β1NE1+β1NI1Δt0.whenj=1, S2=S1+μ2NΔt1/1+μ2+β2NE2+β2NI2Δt1.Now, substitutingS1 value in S2, we obtain(18)S2=Nμ2Δt11+μ1+Nβ1E1+Nβ1I1Δt0+Nμ1Δt0+S0∏m=121+μm+NβmEm+NβmImΔtm−1.Similarly, we can find the value of S3 as(19)S3=∑z=12Nμz+1Δtz∏m=1z1+μm+NβmEm+NβmImΔtm−1∏m=131+μm+NβmEm+NβmImΔtm−1+Nμ1Δt0+S0∏m=131+μm+NβmEm+NβmImΔtm−1.Finally, we obtain the general solution as(20)Sk=S0∏m=1k1+μm+NβmEm+NβmImΔtm−1+∑z=1kNμzΔtz−1∏m=0z−11+μm+NβmEm+NβmImΔtm−1∏m=1k1+μm+NβmEm+NβmImΔtm−1,where Ix=Ex=μx=βx=tx=0 for x≤0.In a similar way, we can easily find the other parameters as(21)Ek=E0∏m=1k1+πm+μm+γm−NβmSmΔtm−1+∑z=1kNβzSzIzΔtz−1∏m=0z−11+πm+μm+γm−NβmSmΔtm−1∏m=1k1+πm+μm+γm−NβmSmΔtm−1,Ik=I0+∑z=1kπzEzΔtz−1∏m=0z−11+σm+μmΔtm−1∏m=1k1+σm+μmΔtm−1,Qk=Q0+∑z=1kγzEz+σzIzΔtz−1∏m=0z−11+θm+μmΔtm−1∏m=1k1+θm+μmΔtm−1,Rk=R0+∑z=1kθzQzΔtz−1∏m=0z−11+μmΔtm−1∏m=1k1+μmΔtm−1,where Sx=Ex=Ix=Qx=Rx=0 and μx=βx=πx=γx=σx=θx=tx=0 for x≤0.Theorem 8.
Assume0<μmin≤μj≤μmax<1, 0<βmin≤βj≤βmax<1, 0<πmin≤πj≤πmax<1, 0<γmin≤γj≤γmax<1, 0<σmin≤σj≤σmax<1, and 0<θmin≤θj≤θmax<1 and 0<tj+1−tj≤1 holds for all j∈1,2,…,M−1, and S1>0, E1>0, I1>0, Q1>0 and R1>0. The implicit solution scheme (16) uniquely solvable for all j∈1,2,…,M−1, and we have(22)Ej+1=−B±B2−4AC2A,j∈1,2,…,M−1,where A, B, and C are given in upcoming equations in (27)–(29).Proof.
SubstitutingSj+1 in Ej+1 and taking Sj+1=S¯, Ej+1=E¯, Ij+1=I¯, Qj+1=Q¯ and Rj+1=R¯, we get(23)E¯=Ej+EjXT+I¯EjNβT+βNTNμT+Sj1+XT+ϕT+ϕT2X−NβNμT+SjT+I¯NβT+ϕT2Nβ,where π+μ+γ=ϕ, μ+NβE¯=X, and Δtj=T.
Again, substituting theI¯ value above, we arrive at(24)E¯=NumeratorN˜DenominatorD˜,(25)N˜=Ej+EjμTδ+EjNβT+βNTNμT+SjIj+E¯NβEjTδ+EjNβT+βNTNμT+SjπT,D˜=1+μT+ϕT+ϕT2μ−NβNμT+SjTδ+IjNβT+ϕT2Nβ+E¯NβT+ϕT2Nβδ+πTNβT+ϕT2Nβ,where Δtj=T,1+σ+μT=δ.
We get the quadratic equation in the form when we solve equation (24)(26)AE¯2+BE¯+C=0⋅i.eAEj+12+BEj+1+C=0,where(27)A=NβΔtj+π+μ+γΔtj2Nβ1+σ+μΔtj+πΔtjNβΔtj+π+μ+γΔtj2Nβ,(28)B=1+μΔtj+π+μ+γΔtj+π+μ+γμΔtj2−NβNμΔtj+SjΔtj1+σ+μΔtj+IjNβΔtj+π+μ+γΔtj2Nβ−NβEjΔtj1+σ+μΔtj−EjNβΔtj+βNΔtjNμΔtj+SjπΔtj,and(29)C=Ej+EjμΔtj1+σ+μΔtj+EjNβΔtjNβΔtjNμΔtj+SjIj.
The proof is completed by taking the roots of equation (26).
Similarly, by substitutingIj+1 in Qj+1 in the system, the equation Qj+1 is in the function of Ej+1. Then, replacing Ej+1 by equation (24), we will get an explicit formula for computing Q¯.
SinceSj+1+Ej+1+Ij+1+Qj+1+Rj+1=μj+1N, we can easily get an explicit formula for Rj+1 for j∈1,2,…,M−1.


## 2.1. Mathematical Background Material

Here, we revisit the Lipschitz continuity of a function on Euclidean spaces, the local Lipchitz condition, Banach’s fixed point theorem, and the method of variation of the parameter, which will be used in the subsequent sections.Definition 1.
(see [46]). Let q1 and q2 be two positive integers and D⊂ℝq1. A function H:D⟶ℝq2 is said to be Lipchitz continuous on D if there exists a nonnegative constant L≥0 such that Hx−Hyℝq2≤L⋅x−yℝq1 holds for x,y∈D.
(i) LetU⊂ℝq1 be an open set and H:U⟶ℝq2. Then, H is called as locally Lipchitz continuous if for every element y0∈U there exists a neighborhood V of y0 such that the restrictions of H to V are Lipchitz continuous on V.In a more general framework, we consider a nonlinear initial value problem (IVP) z′t = Ht,zt; z0 = z0, where zt = y1t,y2t,…,yntT is solution vector, and Ht,zt = h1t,zt,h2t,zt,…,hnt,ztT is vectorial function with initial point z0∈ℝn. The following theorem, which is a direct consequence of Gronwall’s lemma, can be used to prove global existence in time.Theorem 1 (see [17]).
LetG:ℝn⟶ℝn be locally Lipchitz conditions. If there exist nonnegative real constants B and k such that Gt,zt≤kzt∞+B holds for all zt∈ℝn, then the solution of the initial value problem z′t=Gt,zt,z0=z0, exists for all time t∈ℝ, and moreover, zt satisfies ztℝn≤z0ℝn⋅expk⋅t+B/kexpk⋅t−1,∀t∈ℝ.
We give the following: Banach’s fixed point theorem, which will be used to preserve the global uniqueness in time [47, 47].Theorem 2 (see [48]).
LetX,μ be a complete metric space. Let T:X⟶X be a strict contraction, that is, there exists a constant K∈0,1 such that μTx,Ty≤K⋅μx,y holds for all x,y∈X. Then, the mapping T has a unique fixed point.
In the following theorem, we present the Grownwall and Bellman inequality, which will be used in the subsequence theorems related to the continuous functions.Theorem 3 (see [49]).
LetI:=a,b, u,f:I⟶0,∞ be two continuous and nonnegative functions and g:I⟶0,∞ be a continuous, positive, and nondecreasing function. If the inequality(1)ut≤gt+∫atfs⋅usds,holds for all t∈I, we have(2)ut≤gt⋅exp∫atfs⋅usds.Theorem 4.
(Method of variation of parameter) For a first-order nonhomogeneous linear differential equation,y′t+ptyt=ft has the general solution yt=vtePt+AePt, where v′t=e−Ptft, Pt=−∫ptdt and A is an arbitrary constant.


## 2.2. Continuous Problem Formulation

At first, let us assume the following assumptions [50, 51] for the upcoming calculations.(i)
Let the population size varies over time beN is varying over time (i.e., population size = μNt for all t∈0,∞).(ii)
We divide the population into five homogeneous subgroups, namely susceptible people (S), exposed (E), infectious (I), quarantined (Q), and recovered (R). We can clearly assign every individual to exactly one subgroup. Hence, we obtain S,E,I,Q,R model satisfying the condition μN=St+Et+It+Qt+Rt for all t∈0,∞.(iii)
Each time-varying transmission ratex:0,∞⟶0,∞ is Lipchitz continuous and continuously differentiable, and there exist constants xmin and xmax such that 0<xmin≤xt≤xmax for all t≥0 and x∈π,β,γ,σ,θ,μ.The choice of time-dependent transmission rates is possible because the countermeasures such as lockdowns, social distancing, or other political actions like curfews and different medical treatments reduce possible contact between susceptible and infectious people.Our equations of the time-continuous SEIQR model read as follows:(3)dStdt=μN−μSt−βNStEt+It,dEtdt=βNStEt+It−πEt−μ+γEtdItdt=πEt−σIt−μIt,dQtdt=γEt+σIt−θQt−μQt,dRtdt=θQt−μRt,with initial conditions S0=S0>0,E0=E0>0,I0=I0>0,Q0=Q0>0 and R0=R0≥0. The detailed parameters and description are given in Table 2Table 2
Parameters and description.
ParametersDescriptionS(t)At timet, the number of susceptible peopleE(t)At timet, the number of exposed peopleI(t)At timet, the number of infected peopleQ(t)At timet, the number of quarantined peopleR(t)At timet, the number of recovered peopleβThe rate at which susceptible populations migrate to exposed and infected populationsπThe rate at which an exposed population moves to an infected populationγTransmission rate at which exposed people take outside as isolatedσTransmission rate at which infected people were added to isolated individualθTransmission rate at which isolated persons recoveredμNatural death rate and disease-related death rate


## 2.3. Nonnegativity and Boundedness of Solutions

Now, we prove the boundedness of the solution to (3). For this purpose, we modify ideas given in [51, 52] deriving the following lemmas, so consider the bounded, time-varying transmission rates given above.Lemma 1.
Each solution of system (3) is bounded below by zero.Proof.
Consider, the first relation of (3),(4)dStdt=μN−μSt−βNStEt+It.
By takingdSt/dt=S′t, equation (4) can be expressed as a first-order nonhomogeneous linear differential equation in St as(5)S′t+μ+βNEt+βNItSt=μN.
Applying Theorem4, and by applying the same procedure to the first-order nonhomogenous linear equation in Et,(6)dEtdt−βNSt−π−μ−γEt=βNStIt.
We can easily show thatEt≥0 for all t∈0,∞. Proceeding like this, we can show that It,Qt, and Rt≥0 for all t∈0,∞.
SinceμNt=St+Et+It+Qt+Rt, and total population is finite, St,Et,It,Qt, and Rt are bounded above, and hence, the proof is complete.Theorem 5.
For all solution functions of (3), we have 0≤X≤μN, where X∈St,Et,It,Qt,Rt.Proof.
The proof follows fromμNt=St+Et+It+Qt+Rt and Lemma 1.


## 2.4. Global Existence in Time

We arrive at a theorem regarding global existence in the time of (3) based on Theorem 1. For abbreviation, we use the supremum norm ft∞:=supt∈a,bft for an arbitrary continuous function f:a,b⟶ℝ. A similar definition holds for vector-valued bounded functions. In our case, using the boundedness of S,E,I,Q,R on 0,∞, obtain the following global existence theorem.Theorem 6.
The system of nonlinear first-order ODE (3) has at least one solution which exists for all t≥0.Proof.
By denotingzt=St,Et,It,Qt,Rt, we can set
G:0,∞×ℝ5⟶ℝ5 by(7)t,zt⟶A−μSt−βNStEt+ItβNStEt+It−πEt−μ+γEtπEt−σIt−μItγEt+σIt−θQt−μQtθQt−μRtT.
Clearly,G is Lipchitz continuous, due to the continuity of each components.
Assuming the supremum norm on our Euclidean space, and with the help of triangle inequality, we arrive(8)Gt,zt∞=supt∈0,∞A−μSt−βStEt+It,βNStEt+It−βEt−μ+γEt,πEt−σIt−μIt,γEt+σIt−θQt−μQt,θQt−μRt,Gt,zt∞≤supt∈0,∞A−μSt,βNStEt+It,πEt,γEt+σIt,θQt−μRt,≤supt∈0,∞μmaxN+μmaxSt,βmaxNStEt+It,πmaxEt,γmaxEt+σmaxIt,θmaxQt−μmaxRt,≤supt∈0,∞μmax+βmax+πmax+γmax+σmax+θmaxSt,Et,It,Qt,Rt,≤supt∈0,∞μmax+βmax+πmax+γmax+σmax+θmaxzt∞,Gt,zt∞≤kzt∞,where k=μmax+βmax+πmax+γmax+σmax+θmax.
From the boundedness of our solution functions and the boundedness of our time-varying transmission rates, all requirements of Theorem1 are fulfilled, and our proof is complete.


## 2.5. Global Uniqueness in Time

We present the global uniqueness theorem for (3) by utilizing the inductive application of Banach’s fixed point theorem.Theorem 7.
The nonlinear ODE system (3) has a unique solution that exists for all t≥0.Proof.
Consider the system of equations given in (3):(1)
Consider the time interval0,τ is applicable to Banach’s fixed point theorem(2)
Forx1,x2,y1,y2∈ℝ, by triangle inequality, we have x1⋅y1−x2⋅y2≤x1‖y1−y2+y2‖x1−x2(3)
We assume thatS,E,I,Q,R,S¯,E¯,I¯,Q¯,R¯:0,∞⟶0,∞ are two solutions of (3)
Beginning the proof by letting(9)supt∈0,τSt−S¯t=supS1ePt−S¯1ePtwherePt=−∫Ptdt,supt∈0,τSt−S¯t=supS1−S¯1‖ePt,supt∈0,τSt−S¯t≤S1−S¯1τzt−zt¯∞.
Then, the second equation in (3) becomes(10)dEtdt−βNSt−π−μ−γEt=βNStIt.
Since it is a first-order nonhomogenous linear equation inEt, we take p1t=−βNSt−π−μ−γ and ft=βNStIt. As P1t=−∫p1tdt and v1t=∫e−P1tβNStItdt. Hence, Et=v1tdteP1t+E1ePt and(11)supt∈0,τEt−E¯t≤E1−E¯1τzt−zt¯∞.
Similarly, we can easily show the following inequalities:(12)supt∈0,τQt−Q¯t≤Q1−Q¯1τzt−zt¯∞,supt∈0,τIt−I¯t≤I1−I¯1τzt−zt¯∞,supt∈0,τRt−R¯t≤R1−R¯1τzt−zt¯∞.
Summing implies(13)zt−zt¯∞≤∑S,E,I,Q,RS1−S¯1τzt−zt¯∞.
By choosingτ<1/21/∑S1−S¯1. As a result, we conclude that the solution is unique in 0,τ.


## 2.6. Time-Discrete Implicit SEIQR Model

Assume that our time interval0,τ can be divided by strictly increasing sequence tjj=1M for M∈ℕ with t1=0 and tM=T. We write ftj=fj, j∈1,2,3,…,M for an arbitrary time-dependent function f.


## 2.7. Discussion and Formulations

Here, we are transforming the continuous system (3) to the fully explicit discrete scheme (14) as given as follows:(14)ΔSjΔtj=μN−μj+1Sjt−βj+1NSjtEjt+Ijt,ΔEjΔtj=βj+1NSjtEj+Ij−πj+1Ej−μj+1+γj+1Ej,ΔIjΔtj=πj+1Ej−σj+1Ij−μj+1Ij,ΔQjΔtj=γj+1Ej+σj+1Ij−θj+1Qj−μj+1Qj,ΔRjΔtj=θj+1Qj−μj+1Rj,and a fully implicit scheme (15) as(15)ΔSjΔtj=μj+1N−μj+1Sj+1−βj+1NSj+1Ej+1+Ij+1t,ΔEjΔtj=βj+1NSj+1Ej+1+Ij+1−πj+1+μj+1+γj+1Ej+1,ΔIjΔtj=πj+1Ej+1−σj+1Ij+1−μj+1Ij+1,ΔQjΔtj=γj+1Ej+1+σj+1Ij+1−θj+1Qj+1−μj+1Qj+1,ΔRjΔtj=θj+1Qj+1−μj+1Rj+1,where ΔXj/Δtj = Xj+1−Xj/tj+1−tj for j∈1,2,3,…,M−1.Observe thatμj+1N=Sj+1+Ej+1+Ij+1+Qj+1+Rj+1 for all j∈1,2,3,…,M−1.Since the fully explicit scheme (14) simply reduces to a linear system, our main interest is in a fully implicit discrete scheme because it preserves the nonlinear structure of the continuous problem.


## 2.8. Implicit Time-Discrete Problem Formulation

In this section and subsequent sections, we derive the recurrence type of solutions for the implicit scheme (15). For that, we assume that 0<xmin≤xj≤xmax, x∈μ,β,π,γ,σ,θ, for j∈1,2,3,…,M and that 0<tj+1−tj≤1 for j∈1,2,3,…,M and that S0>0,E0>0,I0>0,Q0>0,R0≥0. Now, (15) can be expressed as(16)Sj+1=Sj+μj+1Ntj+1−tj1+μj+1+βj+1NEj+1+βj+1NIj+1tj+1−tj,Ej+1=Ej+βj+1NSj+1Ij+1tj+1−tj1+πj+1+μj+1+γj+1−βNSj+1tj+1−tj,Ij+1=Ij+πj+1Ej+1tj+1−tj1+σj+1+μj+1tj+1−tj,Qj+1=Qj+γj+1Ej+1+σj+1Ij+1tj+1−tj1+θj+1+μj+1tj+1−tj,Rj+1=Rj+θj+1Qj+1tj+1−tj1+μj+1tj+1−tj.


## 2.9. Unique Solvability

In this section, we provide the method for finding the solution of (16). Letting tj+1−tj=Δtj in Sj+1, we obtain(17)Sj+1=Sj+μj+1NΔtj1+μj+1+βj+1NEj+1+βj+1NIj+1Δtj,whenj=0, S1=S0+μ1NΔt0/1+μ1+β1NE1+β1NI1Δt0.whenj=1, S2=S1+μ2NΔt1/1+μ2+β2NE2+β2NI2Δt1.Now, substitutingS1 value in S2, we obtain(18)S2=Nμ2Δt11+μ1+Nβ1E1+Nβ1I1Δt0+Nμ1Δt0+S0∏m=121+μm+NβmEm+NβmImΔtm−1.Similarly, we can find the value of S3 as(19)S3=∑z=12Nμz+1Δtz∏m=1z1+μm+NβmEm+NβmImΔtm−1∏m=131+μm+NβmEm+NβmImΔtm−1+Nμ1Δt0+S0∏m=131+μm+NβmEm+NβmImΔtm−1.Finally, we obtain the general solution as(20)Sk=S0∏m=1k1+μm+NβmEm+NβmImΔtm−1+∑z=1kNμzΔtz−1∏m=0z−11+μm+NβmEm+NβmImΔtm−1∏m=1k1+μm+NβmEm+NβmImΔtm−1,where Ix=Ex=μx=βx=tx=0 for x≤0.In a similar way, we can easily find the other parameters as(21)Ek=E0∏m=1k1+πm+μm+γm−NβmSmΔtm−1+∑z=1kNβzSzIzΔtz−1∏m=0z−11+πm+μm+γm−NβmSmΔtm−1∏m=1k1+πm+μm+γm−NβmSmΔtm−1,Ik=I0+∑z=1kπzEzΔtz−1∏m=0z−11+σm+μmΔtm−1∏m=1k1+σm+μmΔtm−1,Qk=Q0+∑z=1kγzEz+σzIzΔtz−1∏m=0z−11+θm+μmΔtm−1∏m=1k1+θm+μmΔtm−1,Rk=R0+∑z=1kθzQzΔtz−1∏m=0z−11+μmΔtm−1∏m=1k1+μmΔtm−1,where Sx=Ex=Ix=Qx=Rx=0 and μx=βx=πx=γx=σx=θx=tx=0 for x≤0.Theorem 8.
Assume0<μmin≤μj≤μmax<1, 0<βmin≤βj≤βmax<1, 0<πmin≤πj≤πmax<1, 0<γmin≤γj≤γmax<1, 0<σmin≤σj≤σmax<1, and 0<θmin≤θj≤θmax<1 and 0<tj+1−tj≤1 holds for all j∈1,2,…,M−1, and S1>0, E1>0, I1>0, Q1>0 and R1>0. The implicit solution scheme (16) uniquely solvable for all j∈1,2,…,M−1, and we have(22)Ej+1=−B±B2−4AC2A,j∈1,2,…,M−1,where A, B, and C are given in upcoming equations in (27)–(29).Proof.
SubstitutingSj+1 in Ej+1 and taking Sj+1=S¯, Ej+1=E¯, Ij+1=I¯, Qj+1=Q¯ and Rj+1=R¯, we get(23)E¯=Ej+EjXT+I¯EjNβT+βNTNμT+Sj1+XT+ϕT+ϕT2X−NβNμT+SjT+I¯NβT+ϕT2Nβ,where π+μ+γ=ϕ, μ+NβE¯=X, and Δtj=T.
Again, substituting theI¯ value above, we arrive at(24)E¯=NumeratorN˜DenominatorD˜,(25)N˜=Ej+EjμTδ+EjNβT+βNTNμT+SjIj+E¯NβEjTδ+EjNβT+βNTNμT+SjπT,D˜=1+μT+ϕT+ϕT2μ−NβNμT+SjTδ+IjNβT+ϕT2Nβ+E¯NβT+ϕT2Nβδ+πTNβT+ϕT2Nβ,where Δtj=T,1+σ+μT=δ.
We get the quadratic equation in the form when we solve equation (24)(26)AE¯2+BE¯+C=0⋅i.eAEj+12+BEj+1+C=0,where(27)A=NβΔtj+π+μ+γΔtj2Nβ1+σ+μΔtj+πΔtjNβΔtj+π+μ+γΔtj2Nβ,(28)B=1+μΔtj+π+μ+γΔtj+π+μ+γμΔtj2−NβNμΔtj+SjΔtj1+σ+μΔtj+IjNβΔtj+π+μ+γΔtj2Nβ−NβEjΔtj1+σ+μΔtj−EjNβΔtj+βNΔtjNμΔtj+SjπΔtj,and(29)C=Ej+EjμΔtj1+σ+μΔtj+EjNβΔtjNβΔtjNμΔtj+SjIj.
The proof is completed by taking the roots of equation (26).
Similarly, by substitutingIj+1 in Qj+1 in the system, the equation Qj+1 is in the function of Ej+1. Then, replacing Ej+1 by equation (24), we will get an explicit formula for computing Q¯.
SinceSj+1+Ej+1+Ij+1+Qj+1+Rj+1=μj+1N, we can easily get an explicit formula for Rj+1 for j∈1,2,…,M−1.


## 3. Monotonicity Properties and Long-Time Behaviour

In this section, we develop a suitable atmosphere in which our implicit scheme obeys the monotonic properties as in the continuous case. For this, we give the following lemmas and finally provide a nonlinear programming problem to optimize the transmission sequences.Lemma 2.
IfSjβjμj−1Ej+Ij+N−1≥1, then Sj+1≤Sj, i.e., Snn=1∞ becomes a decreasing sequence.Proof.
Takingy=μj+1NΔtj and x=1+μj+1+βj+1NEj+1+Ij+1 in Sj+1 of (16). We arrive the relation(30)Sj+1=y+Sj1+x.
From (30), we obtain 1+xSj+1=y+Sj and which implies Sj+1−Sj=y−xSj+1.
Clearly,y−xSj+1≤0 if Sj+1βj+1μj+1−1Ej+1+Ij+1+N−1≥1. Thus,(31)Sj+1−Sj≤0,forj∈1,2,…,m,and the proof is complete.Lemma 3.
IfβjNSj1+IjEj−1≤πj+μj+γj, then Enn=1∞ is a decreasing sequence.Proof.
Takingx = βj+1NSj+1Ij+1Δtj and y = πj+1+μj+1+γj+1−βj+1NSj+1Δtj in Ej+1 of (16), we derive(32)Ej+1−Ej=x−yEj+1.
Here,x−yEj+1≤0 follows from βj+1NSj+11+Ij+1Ej+1−1≤πj+1+μj+1+γj+1, and hence, we obtain Ej+1≤Ej for all j, and the proof is complete.Lemma 4.
Ifπjσj+μj−1≤IjEj−1, then Inn=1∞ be a decreasing sequence.Proof.
TakingIj+1=Ij+x/1+y, where x=πj+1Ej+1Δtj and y=σj+1+μj+1Δtj in (16), then we find(33)Ij+1−Ij=x−yIj+1.
Sinceπj+1σj+1+μj+1−1≤Ij+1Ej+1−1, we get x−yIj+1≤0 and Ij+1≤Ij for all j.Lemma 5.
IfγjEj+σjIjθj+μj−1≥Qj, then Qjn=1∞ becomes an increasing sequence.Proof.
TakingQj+1=Qj+x/1+y, where x=γj+1Ej+1+σj+1Ij+1Δtj and y=θj+1+μj+1Δtj in (16), then it is easy to arrive(34)Qj+1Qj=x−yQj+1.From our assumption γj+1Ej+1+σj+1Ij+1θj+1+μj+1−1≥Qj+1, we derive(35)x−yQj+1≥0,and the proof is complete.Lemma 6.
IfθjQjμj−1≥Rj for all j, then Rj+1≥Rj for all j.Expressing Rj+1=Rj+x/1+y in (16), we can easily find(36)Rj+1−Rj=x−yRj+1,from the given condition, θj+1Qj+1μj+1−1≥Rj+1 and straightforward calculations, we obtain(37)x−yRj+1≥0,and the proof is complete.Remark 1.
SinceRj is monotonic and bounded above by total population, then it will converge and limj⟶∞Rj=R∗ exists, and Ej and Sj are decreasing sequence, and we easily observe that limj⟶∞Ij=0.


### 3.1. Formulation and Discussion

Due to the ongoing nature of the COVID-19 pandemic, it was impossible to fully comprehend the short- or long-term implications of this global disruption. In our study, when there is no quarantine, a single infected individual can spread the infection to about two other people; however, when quarantine is imposed, there is a chance of preventing further transmission of infection. However, some of the exposed individuals may avoid quarantine due to fear of stigma and death. In other words, this does not achieve zero infection in the population, implying that additional interventions are required to eradicate the virus. If there are no adequate interventions in place, the virus will remain in the population for a long time, but it will eventually drop over time. But, still, there will be a small number of sick people who have the ability to start another outbreak even after measures like quarantine and public health education/awareness raise the number of exposed and infected people dramatically but not to zero. According to this, COVID-19 will not be completely eradicated even with prompt development of measures.In order to study the effects of isolation, quarantine, and the percentage of exposed people who will be quarantined, we did numerical simulations. Many authors have developed numerous mathematical models to limit the spread of viruses. Here, we have developed the optimization technique to control the spread of the virus. This approach enables us to control a viruses future spread and predict how it will spread in the future.After summarizing all of the prior principles, the problem is transformed into a dynamic programming problem model with constraints imposed by the previous lemmas by keepingSk,Ek,Ik,Qk, and Rk are constants and βk,πk,μk,γk,σk, and θk are variables at each level of the optimization. Note that the following dynamic programming problem preserves the monotonic properties.The dynamic programming problem is given byMinβk+πk+μk−γk+σk+θk and subject to the constraints:(38)SkβkEk+Ik+μkN−1−μk≥0,πk+μk+γkEk−βkNSkEk+Ik≥0,Ikσk+μk−πkEk≥0,γkEk+σkIk−Qkθk+μk≥0,θkQk−Rkμk≥0,Sk+Ek+Ik+Qk+Rk−μkN=0,Sk,Ek,Ik,Qk,Rk≥0,βk,πk,μk,γk,σk,θk≥0.SinceS0,E0,I0,Q0, and R0 are known initial conditions, we get a dynamic programming problem for level-0 time if we keep k=0 in the above model. We shall obtain an optimal (feasible) solution β0∗,π0∗,μ0∗,γ0∗,σ0∗,θ0∗ by employing the optimization technique in operation research. Then, these values are assigned as β1=β0∗, π1=π0∗, μ1=μ0∗, γ1=γ0∗, σ1=σ0∗, θ1=θ0∗. We will receive the values for finding S1,E1,I1,Q1, and R1. We will receive a dynamic programming problem for level-1 time if we keep k=1 in this model. We will achieve the ideal solution as β1∗,π1∗,μ1∗,γ1∗,σ1∗,θ1∗ using the optimization technique. These values are assigned as β2=β1∗,π2=π1∗,μ2=μ1∗,γ2=γ1∗,σ2=σ1∗,θ2=θ1∗. We will receive the values for S2,E2,I2,Q2, and R2. If we keep going in this direction, we will end up with transmission rate sequences βkk=1∞,πkk=1∞,μkk=1∞,γkk=1∞,σkk=1∞, and θkk=1∞ which provide a sufficiently viable stable solution for the situation, correspond to S0,E0,I0,Q0, and R0. After a certain stage, each transmission rate becomes constant. We must follow the appropriate standard operating procedures and safety precautions in order to acquire these sequences in practice. The isolation class, it appears, plays a significant role in achieving this possible solution.


## 3.1. Formulation and Discussion

Due to the ongoing nature of the COVID-19 pandemic, it was impossible to fully comprehend the short- or long-term implications of this global disruption. In our study, when there is no quarantine, a single infected individual can spread the infection to about two other people; however, when quarantine is imposed, there is a chance of preventing further transmission of infection. However, some of the exposed individuals may avoid quarantine due to fear of stigma and death. In other words, this does not achieve zero infection in the population, implying that additional interventions are required to eradicate the virus. If there are no adequate interventions in place, the virus will remain in the population for a long time, but it will eventually drop over time. But, still, there will be a small number of sick people who have the ability to start another outbreak even after measures like quarantine and public health education/awareness raise the number of exposed and infected people dramatically but not to zero. According to this, COVID-19 will not be completely eradicated even with prompt development of measures.In order to study the effects of isolation, quarantine, and the percentage of exposed people who will be quarantined, we did numerical simulations. Many authors have developed numerous mathematical models to limit the spread of viruses. Here, we have developed the optimization technique to control the spread of the virus. This approach enables us to control a viruses future spread and predict how it will spread in the future.After summarizing all of the prior principles, the problem is transformed into a dynamic programming problem model with constraints imposed by the previous lemmas by keepingSk,Ek,Ik,Qk, and Rk are constants and βk,πk,μk,γk,σk, and θk are variables at each level of the optimization. Note that the following dynamic programming problem preserves the monotonic properties.The dynamic programming problem is given byMinβk+πk+μk−γk+σk+θk and subject to the constraints:(38)SkβkEk+Ik+μkN−1−μk≥0,πk+μk+γkEk−βkNSkEk+Ik≥0,Ikσk+μk−πkEk≥0,γkEk+σkIk−Qkθk+μk≥0,θkQk−Rkμk≥0,Sk+Ek+Ik+Qk+Rk−μkN=0,Sk,Ek,Ik,Qk,Rk≥0,βk,πk,μk,γk,σk,θk≥0.SinceS0,E0,I0,Q0, and R0 are known initial conditions, we get a dynamic programming problem for level-0 time if we keep k=0 in the above model. We shall obtain an optimal (feasible) solution β0∗,π0∗,μ0∗,γ0∗,σ0∗,θ0∗ by employing the optimization technique in operation research. Then, these values are assigned as β1=β0∗, π1=π0∗, μ1=μ0∗, γ1=γ0∗, σ1=σ0∗, θ1=θ0∗. We will receive the values for finding S1,E1,I1,Q1, and R1. We will receive a dynamic programming problem for level-1 time if we keep k=1 in this model. We will achieve the ideal solution as β1∗,π1∗,μ1∗,γ1∗,σ1∗,θ1∗ using the optimization technique. These values are assigned as β2=β1∗,π2=π1∗,μ2=μ1∗,γ2=γ1∗,σ2=σ1∗,θ2=θ1∗. We will receive the values for S2,E2,I2,Q2, and R2. If we keep going in this direction, we will end up with transmission rate sequences βkk=1∞,πkk=1∞,μkk=1∞,γkk=1∞,σkk=1∞, and θkk=1∞ which provide a sufficiently viable stable solution for the situation, correspond to S0,E0,I0,Q0, and R0. After a certain stage, each transmission rate becomes constant. We must follow the appropriate standard operating procedures and safety precautions in order to acquire these sequences in practice. The isolation class, it appears, plays a significant role in achieving this possible solution.


## 4. Error Analysis

Now, we will set an upper limit for error propagation. We need to construct certain assumptions for our convergence analysis before proving the required statements. The following is a list of them:(i)
Let0,T be the time interval under consideration with t1=0<t2<⋯<tM−1<tM=T(ii)
Allow the time-continuous and time-discrete models beginning circumstances to coincide(iii)
LetS,E,I,Q,R:0,T⟶0,∞ be twice continuously differentiable solution functions(iv)
Allow the time-varying transmission ratesμ,β,π,γ,σ,θ:0⋅T⟶0,∞ be continuously differentiable just once(v)
Allow the time-varying transmission and recovery rates to be bounded by 0 and 1(vi)
ChooseΔp<min1/12μmax+βmax+4πmax+γmax+σmax+θmax≤1 for all p∈N and Δ:=maxp∈NΔpWe get the following theorem under these conditions, in which we adopt notions from the error analysis of an explicit-implicit solution algorithm.Theorem 9.
The difference between the solutions of the time-continuous system formulation (3) and the time-discrete system (16) fulfills if the aforementioned assumptions are met, then(39)zp+1−ztp+1≤Cloc⋅Δ11−6μmax+βmax+2πmax+θmax+γmax+σmaxΔp−1.Proof.
Since this is technical proof, we will start with a brief description of our technique. The first step is the estimation of local errors between time-continuous and time-discrete solutions. After that, we look at error propagation over time. Finally, we look into the accumulation of these errors. At the same time, time-discrete solutions are expressed asSp at time tp, whereas time-continuous solutions are written as Stp.
1) For the purpose of examining local errors, we assume that(40)tp,SpT=tp,StpT,tp,EpT=tp,EtpT,tp,IpT=tp,ItpT,tp,QpT=tp,QtpT,tp,RpT=tp,RtpT,hold for arbitrary p∈1,2,…,M−1 on the time interval tp,tp+1. Here, we consider one time step and denote corresponding time-discrete solutions by Sp+1˜, Ep+1˜, Ip+1˜, Qp+1˜, and Rp+1˜.
a) For the time-discrete solutionSp at time tp, we take(41)Sp+1˜=Sp+μp+1NΔtp+11+μp+1+βp+1NEp+1˜+βp+1NIp+1˜Δtp+1,and in similar way for the time-continuous solutions Stp, we have(42)Sp+1˜=Stp−μp+1+βp+1NEp+1˜+βp+1NIp+1˜Δtp+1Stp1+μp+1+βp+1NEp+1˜+βp+1NIp+1˜Δtp+1+μp+1NΔtp+11+μp+1+βp+1NEp+1˜+βp+1NIp+1˜Δtp+1,which yields(43)Stp+1−Sp+1˜=Stp+1−Stp+μp+1+βp+1NEp+1˜+βp+1NIp+1˜Δtp+1Stp−μp+1NΔtp+11+μp+1+βp+1NEp+1˜+βp+1NIp+1˜Δtp+1,which implies(44)∫tptp+1S′τdτ+μp+1+βp+1NEp+1˜+βp+1NIp+1˜Δtp+1Stp−μp+1NΔtp+11+μp+1+βp+1NEp+1˜+βp+1NIp+1˜Δtp+1.
Now, add and subtractΔp+1⋅S′tp. And then, applying the triangle inequality in the above relation, we get(45)Stp+1−Sp+1˜≤∫tptp+1S′τdτ−Δp+1⋅S′tp+Δp+1⋅S′tp+μp+1+βp+1NEp+1˜+βp+1NIp+1˜Δtp+1Stp−μp+1NΔtp+11+μp+1+βp+1NEp+1˜+βp+1NIp+1˜Δtp+1IS,1≔=∫tptp+1S′τ−S′tpdτ=∫tptp+1τ−tpS′τ−S′tpτ−tpdτ.
By applying the mean value theorem, there existsξS,1∈tp,tp+1 such that(46)S″ξS,1=S′τ−S′tpτ−tp≤S″t∞.
This yields(47)IS,1≤S″t∞⋅∫tptp+1τ−tpdτ=Δp+122⋅S″t∞.Also,(48)IIS,1≔=Δp+1.S′tp+μp+1+βp+1NEp+1˜+βp+1NIp+1˜Δtp+1Stp−μp+1NΔtp+11+μp+1+βp+1NEp+1˜+βp+1NIp+1˜Δtp+1,IIS,1=Δtp+1⋅S′tp+μp+1+βp+1NEp+1˜+βp+1NIp+1˜Stp−μp+1N1+μp+1+βp+1NEp+1˜+βp+1NIp+1˜.
Substituting equation (42) in the above and then solving the above relation, we get IIS,1≤Δtp+1⋅IIIS,1+IVS,1, where(49)IIIS,1≔=Nμp−μp+1−Stpμp−μp+1+NβpEtp+Itp−βp+1Ep+1˜+Ip+1˜,IVS,1≔=μpN−Stpμp+βpNEtp+Itpμp+1+βp+1NEp+1˜+Ip+1˜Δtp+1.
Now, applying the triangle inequality inIIIS,1, we obtain ≤Nμp+1−μp+Stp⋅VS,1, where(50)VS,1≔=μp−μp+1+NβpEtp+Itp−βp+1Ep+1˜+Ip+1˜,letNμp+1−μp=Ntp+1−tpμp+1−μptp+1−tp,Nμp+1−μp=NΔtp+1μp+1−μptp+1−tp.
By the mean value theorem, there existsξμ,1∈tp,tp+1 such that(51)μ′ξμ,1=μp+1−μptp+1−tp≤μ′t∞.
Hence,(52)N⋅μp+1−μp≤N⋅Δtp+1⋅μ′t∞.
Using the triangle inequality again inVS,1 yields(53)VS,1≤μp+1−μp+NβpEtp−βp+1Ep+1˜+NβpItp−βp+1Ip+1˜.
From equation (52), we can easily find that(54)μp+1−μp≤Δtp+1⋅μ′t∞.
Now, substitutingEp+1˜ value (53), we arrive(55)NβpEtp−βp+1Ep+1˜=NβpEtp−βp+1Etp+βp+1NSp+1˜Ip+1˜Δtp+11+πp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1,≤N⋅IS,2+N⋅IIS,2,where(56)IS,2≔=βpEtp−βp+1Etp1+πp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1,IIS,2≔=βp+12NSp+1˜Ip+1˜Δtp+11+πp+1+μp+1+γp+1−Nβp+1Sp+1˜Δtp+1.
Now,IS,2≤Etp⋅(57)βp+βpπp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1−βp+11+πp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1≤N1−βmaxβp+1−βp+βpπp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1,IS,2≤N⋅Δtp+11−βmaxβ′t∞+βmaxπmax+μmax+γmax+βmax.
Similarly, forIIS,2, it is easy to arrive(58)IIS,2=Nβp+12Sp+1˜Ip+1˜Δtp+11+πp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1≤Nβp+12Sp+1˜Ip+1˜Δtp+11−βmax,IIS,2≤N⋅βmax2⋅Δtp+11−βmax.Applying equation (57) and (58) in (55), we get(59)NβpEtp−βp+1Ep+1˜≤N2Δtp+11−βmaxβ′t∞+βmaxπmax+μmax+γmax+βmax+βmax2.
Proceeding through the similar steps, we can find(60)NβpItp−βp+1Ip+1˜≤N2Δtp+11−βmaxβ′t∞+βmaxσmax+μmax+NΔtp+11−βmaxπmax.Applying equations (54), (59), and (60) in VS,1, we obtain(61)VS,1≤Δtp+1μ′t∞+N2Δtp+11−βmaxβ′t∞+βmaxπmax+μmax+γmax+βmax+βmax2+N2Δtp+11−βmaxβ′t∞+βmaxσmax+μmax+NΔtp+11−βmaxπmax.
Using equations (52) and (61), IIIS,1 will be(62)IIIS,1≤2NΔtp+1⋅μ′t∞+N3Δtp+11−βmaxβ′t∞+N3Δtp+11−βmaxβ′t∞+βmaxσmax+μmax+βmaxπmax+μmax+γmax+βmax+βmax2+NΔtp+11−βmaxπmax.
Now, the termIVS,1 satisfies the inequality(63)IVS,1≤IS,3+Stp⋅IIS,3,where(64)IS,3≔=μpNμp+1+βp+1NEp+1˜+Ip+1˜Δtp+1,IIS,3≔=μp+βpNEtp+Itpμp+1+βp+1NEp+1˜+Ip+1˜Δtp+1,which yields(65)IS,3=NΔtp+1μpμp+1+βp+1NEp+1˜+Ip+1˜,≤NΔtp+1μmaxμmax+βmax,and(66)IIS,3≤NΔtp+1μmax+βmaxμmax+βmax,IIS,3=NΔtp+1μmax+βmax2.
From (65) and (66), we get(67)IVS,1≤NΔtp+1μmaxμmax+βmax+N2Δtp+1μmax+βmax2.
Similarly, from (62) and (66),(68)IIS,1≤2NΔtp+12⋅μ′t∞+N3Δtp+121−βmaxβ′t∞+βmaxπmax+μmax+γmax+βmax+βmax2+N3Δtp+121−βmaxβ′t∞+βmaxσmax+μmax+NΔtp+121−βmaxπmax+NΔtp+12μmaxμmax+βmax+N2Δtp+12μmax+βmax2.Hence,(69)Stp+1−Sp+1˜≤Δp+122⋅S″t∞+2NΔtp+12⋅μ′t∞+N3Δtp+121−βmaxβ′t∞+βmaxπmax+μmax+γmax+βmax+βmax2+N3Δtp+121−βmaxβ′t∞+βmaxσmax+μmax+NΔtp+121−βmaxπmax+NΔtp+12μmaxμmax+βmax+N2Δtp+12μmax+βmax2.b) For the time-discrete solution for Ep at time tp, we take(70)Ep+1˜=Ep+βp+1NSp+1˜Ip+1˜Δtp+11+πp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1,as well as for time-continuous solutions Etp, we find(71)Ep+1˜=Etp−πp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1Etp1+πp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1+βp+1Sp+1˜Ip+1˜NΔtp+11+πp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1.
Now, applying the procedure similar to the casea, we can easily find(72)Etp+1−Ep+1˜≤∫tptp+1E′τdτ−Δtp+1E′tp+Δtp+1E′tp+πp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1Etp1+πp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1−βp+1Sp+1˜Ip+1˜NΔtp+11+πp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1,which becomes(73)Etp+1−Ep+1˜≤Δtp+122E″t∞+Δtp+1⋅IE,1,where(74)IE,1≔=E′tp+πp+1+μp+1+γp+1−βp+1NSp+1˜Etp−βp+1Sp+1˜Ip+1˜N1+πp+1+μp+1+γp+1−βp+1NSp+1˜Δtp+1,which can be written as(75)E′tp1+ϕp+1−βp+1NSp+1˜Δtp+11+ϕp+1−βp+1NSp+1˜Δtp+1+ϕp+1−βp+1NSp+1˜Etp−βp+1Sp+1˜Ip+1˜N1+ϕp+1−βp+1NSp+1˜Δtp+1,and πp+1+μp+1+γp+1=ϕp+1.Thus,(76)IE,1≤11−βmax⋅βpNStpItp−βp+1NSp+1˜Ip+1˜+EtpβpNStp−βp+1NSp+1˜+πp+1−πp+μp+1−μp+γp+1−γp+βpNStpItp+EtpβpNStp−ϕp+1ϕp+1−βp+1NSp+1˜Δp+1,where(77)IIE,2≔=βpNStpItp−βp+1NSp+1˜Ip+1˜,IIIE,2≔=EtpβpNStp−βp+1NSp+1˜+πp+1−πp+μp+1−μp+γp+1−γp,IVE,2≔=βpNStpItp+EtpβpNStp−ϕp+1ϕp+1−βp+1NSp+1˜Δp+1,IIE,2=NβpStpItp−βp+1Ip+1˜Stp+μp+1NΔtp+11+μp+1+βp+1NEp+1˜+βp+1NEp+1˜Δtp+1,≤NβpStpItp−βp+1Ip+1˜Stp1+μp+1+βp+1NEp+1˜+βp+1NEp+1˜Δtp+1+Nβp+1Ip+1˜μp+1NΔtp+11+μp+1+βp+1NEp+1˜+βp+1NEp+1˜Δtp+1,where(78)VE,2≔=βpStpItp−βp+1Ip+1˜Stp1+μp+1+βp+1NEp+1˜+βp+1NEp+1˜Δtp+1,VIE,2≔=βp+1Ip+1˜μp+1NΔtp+11+μp+1+βp+1NEp+1˜+βp+1NEp+1˜Δtp+1.
Similarly, substituting theIp+1˜ value in VE,2 and then solving the VE,2, we get(79)VE,2a=11−βmaxβpItp−βp+1Itp+11−βmaxβpItpσp+1+μp+1Δtp+1+μp+1+βp+1NEp+1˜+βp+1NIp+1˜+σp+1+μp+1μp+1+βp+1NEp+1˜+βp+1NIp+1˜Δ2tp+1,VE,2a≤NΔtp+11−βmaxβ′t∞+NΔtp+11−βmaxβmaxσmax+2μmax+βmax+μmax+2βmaxσmax+μmaxΔtp+1.
We can also easily find a(80)VE,2b≤πmaxΔtp+11−βmax.
Finally,(81)VE,2≤N2Δtp+11−βmaxβ′t∞+NΔtp+11−βmaxπmax+N2Δtp+11−βmaxβmaxσmax+2μmax+βmax+μmax+2βmaxσmax+μmaxΔtp+1.
We can also easily arrive(82)VIE,2≤NΔtp+11−βmaxβmaxμmax,and(83)IIE,2≤N3Δtp+11−βmaxβ′t∞+N3Δtp+11−βmaxβmaxσmax+2μmax+βmax+μmax+2βmaxσmax+μmaxΔtp+1+N2Δtp+11−βmaxβmaxμmax+πmax.
Using the same procedure given above, we get(84)IIIE,2≤N3Δtp+11−βmaxβ′t∞+βmaxμmax+2βmax+N2Δtp+11−βmaxμmax+π′t∞+γ′t∞Δtp+1,and(85)IVE,2≤N3Δtp+1βmaxπmax+μmax+γmax+βmax+NΔtp+1βmax+πmax+μmax+γmax2.
Substituting equations (83)–(85) in IE,1, we arrive(86)IE,1≤N3Δtp+11−βmax2β′t∞+N3Δtp+1βmax1−βmaxσmax+2μmax+βmax+μmax+2βmax+πmax+μmax+γmax+βmax+μmax+2βmaxσmax+μmaxΔtp+1+N2Δtp+11−βmaxβmaxμmax+πmax+μmax+NΔtp+11−βmaxβmax+πmax+μmax+γmax2+Δtp+1π′t∞+γ′t∞.
Substituting equation (86) in (73), we obtain the case.
c) For the time-discrete solutionIp at time tp, let(87)Ip+1˜=Ip+πp+1Ep+1˜Δtp+11+σp+1+μp+1Δtp+1,and for the time-continuous solutions Itp, let(88)Ip+1˜=Itp−σp+1+μp+1Δtp+1Itp1+σp+1+μp+1Δtp+1+πp+1Ep+1˜Δtp+11+σp+1+μp+1Δtp+1.
Using the triangle inequality and mean value theorem, we can easily find(89)Itp+1−Ip+1˜≤Δtp+122I″t∞+Δtp+1I′tp+σp+1+μp+1Δtp+1Itp−πp+1Ep+1˜Δtp+11+σp+1+μp+1Δtp+1,where(90)II,1≔=Δtp+1I′tp+σp+1+μp+1Δtp+1Itp−πp+1Ep+1˜Δtp+11+σp+1+μp+1Δtp+1.
ApplyingI′tp value in (89), we obtain(91)II,1=Δtp+1πpEtp−σp+μpItp+σp+1+μp+1Δtp+1Itp−πp+1Ep+1˜Δtp+11+σp+1+μp+1.
Solving the above equations, we get(92)II,1≤Δtp+1πpEtp−πp+1Ep+1˜+Itpσp+1−σp+μp+1−μp+πpEtp−Itpσp+μpσp+1+μp+1Δtp+1,which yields(93)II,1≤NΔtp+121−βmaxπ′t∞+πmaxπmax+μmax+γmax+2βmax+NΔtp+12σ′t∞+μ′t∞+πmax+σmax+μmaxσmax+μmax.
The case (c) is verified by substituting equation (93) in (89).
d) Let the time-discrete solution forQp at time tp is(94)Qp+1˜=Qp+γp+1Ep+1˜+σp+1Ip+1˜Δtp+11+θp+1+μp+1Δtp+1,and the time-continuous solutions Qtp, we take(95)Qp+1˜=Qtp−θp+1+μp+1Δtp+1Qtp1+θp+1+μp+1Δtp+1+γp+1Ep+1˜+σp+1Ip+1˜Δtp+11+θp+1+μp+1Δtp+1,which implies(96)Qtp+1−Qp+1˜≤Δtp+122Q″t∞+|Q′tpΔp+1+θp+1+μp+1QtpΔtp+1−γp+1Ep+1˜+σp+1Ip+1˜Δtp+11+θp+1+μp+1Δtp+1|,where(97)IQ,1≔=Q′tpΔp+1+θp+1+μp+1QtpΔtp+1−γp+1Ep+1˜+σp+1Ip+1˜Δtp+11+θp+1+μp+1Δtp+1≤Δtp+1γpEtp−γp+1Ep+1˜+σpItp−σp+1Ip+1˜+Qtp‖θp+1−θp+Qtp‖μp+1−μp+γpEtp+σpItp−Qtpθp+μpθp+1+μp+1Δtp+1.
Simplifying the above equation, we arrive(98)II,1≤NΔtp+121−βmaxγ′t∞+σ′t∞+NΔtp+12Q′t∞+μ′t∞+NΔtp+121−βmaxγmaxπmax+μmax+γmax+2βmax+σmaxσmax+μmax+NΔtp+12γmaxθmax+μmax+σmaxθmax+μmax+θmax+μmax2+Δtp+121−βmaxσmaxπmax.
Substituting equation (98) in (96), we get the quarantine case.
e) Let the time-discrete solution forRp at time tp is(99)Rp+1˜=Rp+θp+1Qp+1˜Δtp+11+μp+1Δtp+1,and the time-continuous solutions Rtp, let(100)Rp+1˜=Rtp−μp+1Δtp+1Rtp1+μp+1Δtp+1+θp+1Qp+1˜Δtp+11+μp+1Δtp+1,which gives(101)Rtp+1−Rp+1˜≤Δtp+122R″t∞+IR,1,where(102)IR,1≔=|Δtp+1R′tp+μp+1Δtp+1Rtp−θp+1Qp+1˜Δtp+11+μp+1Δtp+1|.The IR,1 expression can be simplified as follows:(103)IR,1≤NΔtp+121−βmaxθ′t∞+θmax+μmax+Δtp+121−βmaxθmaxγmax+σmax+Nμ′t∞+θmaxμmax+μmax2,
which yields recovered instances by substituting equation (103) into (101).
(69) can be rewritten in the following form:(104)Stp+1−Sp+1˜≤Δ2tp+1Cloc,S.Similarly, we can get the other cases as(105)Etp+1−Ep+1˜≤Δ2tp+1Cloc,E,Itp+1−Ip+1˜≤Δ2tp+1Cloc,I,Qtp+1−Qp+1˜≤Δ2tp+1Cloc,Q,Rtp+1−Rp+1˜≤Δ2tp+1Cloc,R.Definition 2.
DefineCloc:=maxCloc,S;Cloc,E;Cloc,I;Cloc,Q;Cloc,R which holds(106)ztp+1−zp+1˜≤Δp+12Cloc,for local errors on time interval tp,tp+1.
2) In reality,tp,SpT, tp,EpT, tp,IpT, tp,QpT, and tp,RpT do not exactly lie on the graph of the time-continuous solution. Thus, we should examine how procedural errors such as Sp−Stp, Ep−Etp, Ip−Itp, Qp−Qtp, and Rp−Rtp propagate to the p+1th time step. That is,(107)zp+1−ztp+1=zp−ztp+Δtp+1Gtp+1,zp+1−Gtp+1,ztp+1,which implies(108)zp+1−ztp+1∞≤zp−ztp∞+Δtp+1Gtp+1,zp+1−Gtp+1,ztp+1∞.Here,(109)Gtp+1,zp+1−Gtp+1,ztp+1∞=A−μp+1Stp+1−Sp+1−βp+1NStp+1Etp+1+tp+1−Sp+1Ep+1+Ip+1βp+1NStp+1Ep+1+Ip+1−Sp+1Ep+1+Ip+1−πp+1+μp+1+γp+1Etp+1−Ip+1πp+1Etp+1−Ep+1−σp+1+μp+1Itp+1−Ip+1γp+1Etp+1−Ep+1+σp+1Itp+1−Ip+1−θp+1+μp+1Qtp+1−Qp+1θp+1Qtp+1−Qp+1−μp+1Rtp+1−Rp+1∞,which is less than equal to(110)A+μp+1Sp+1−Stp+1∞+βp+1N2‖Sp+1−Stp+1‖∞+Ep+1−Etp+1∞+Ip+1−Itp+1∞βp+1N2Sp+1−Stp+1∞+Ep+1−Etp+1∞+‖Ip+1−Itp+1‖∞+πp+1+μp+1+γp+1Ep+1−Etp+1∞πp+1Ep+1−Etp+1∞+σp+1+μp+1Ip+1−Itp+1∞γp+1Ep+1−Etp+1∞+σp+1Ip+1−Itp+1∞−θp+1+μp+1Qp+1−Qtp+1∞θp+1Qp+1−Qtp+1∞−μp+1Rp+1−Rtp+1∞∞≤6μmax+βmax+2πmax+γmax+θmax⋅zp+1−ztp+1∞.Thus,(111)zp+1−ztp+1∞≤zp−ztp∞+6μmax+βmax+2πmax+γmax+θmax⋅zp+1−ztp+1∞.Hence, we conclude zp+1−ztp+1∞(112)≤zp−ztp∞1−6μmax+βmax+2πmax+γmax+θmaxΔ,with(113)Δ≕=maxp∈1,2,3,…,M−1Δtp+1<112μmax+βmax+4πmax+γmax+θmax.
3) The upper bound between the time-discrete solution and the time-continuous solution is verified below.
Forp=0 in (112), it becomes z1−zt1∞=0.
Forp=1, we get z2−zt2∞≤z2−z2˜∞+z2˜−zt2∞. By applying equations (106) and (112), we get z2−zt2∞(114)≤z1−zt1∞1−6μmax+βmax+2πmax+γmax+θmaxΔ+Cloc⋅Δ2,which becomes(115)z2−zt2∞≤Cloc⋅Δ2.For p=2, we get(116)z3−zt3∞≤z3−z3˜∞+z3˜−zt3∞≤z2−zt2∞1−6μmax+βmax+2πmax+γmax+θmaxΔ+Cloc⋅Δ2≤z2−zt2∞1−6μmax+βmax+2πmax+γmax+θmaxΔ+Cloc⋅Δ2≤∑j=01Cloc⋅Δ21−6μmax+βmax+2πmax+γmax+θmaxΔj.For p=3, we obtain(117)z4−zt4∞≤z4−z4˜∞+z4˜−zt4∞≤z3−zt3∞1−6μmax+βmax+2πmax+γmax+θmaxΔ+Cloc⋅Δ2,which yields z4−zt4∞(118)≤∑j=02Cloc⋅Δ21−6μmax+βmax+2πmax+γmax+θmaxΔj.Proceeding like this, we can write zp+1−ztp+1∞(119)≤∑j=0p−1Cloc⋅Δ21−6μmax+βmax+2πmax+γmax+θmaxΔj,for p∈1,2,3,…,M−2. From (119), we get zp+2−ztp+2∞(120)≤∑j=0pCloc⋅Δ21−6μmax+βmax+2πmax+γmax+θmaxΔj,by induction method.
Now, applying the geometric series in (119), we obtain zp+1−ztp+1∞.(121)≤Cloc⋅Δ2∑j=0p−111−6μmax+βmax+2πmax+γmax+θmaxΔj,=Cloc⋅Δ21/1−6μmax+βmax+2πmax+γmax+θmaxΔp−11/1−6μmax+βmax+2Zπmax+γmax+θmaxΔ−1.
AssumingΔ<1/12μmax+βmax+4πmax+γmax+θmax, we conclude Δ/1−6μmax+βmax+2πmax+γmax+θmaxΔ−1≤1.
Hence, we get (39).


## 5. Conclusion

The present paper is devoted to the analysis and optimization of the SEIQR epidemic model containing an isolation class. We derived both continuous and discrete schemes of the SEIQR model as well as the global existence of solutions and nonnegativity bounded properties for both schemes. Along with this, we have illustrated the solving technique for the discrete scheme and proposed a new optimization technique with the help of a dynamic programming problem. In addition, we have analyzed the error between continuous and discrete schemes in the last section. Finally, we conclude that the isolation class plays a vital role in controlling the COVID-19 pandemic situation. More interestingly, the results also reveal that COVID-19 can exhibit oscillatory behaviour in the future. On the other hand, social distancing methods, quarantine efficiency, and isolation can be used to keep it under control. Future research could look into the effects of current coronavirus mutations like Delta and Omicron on the COVID-19 pandemic’s dynamics. We also suggested an alternative dynamic model of the SEIQR class, which can theoretically be generalized to generate continuous and discrete-time models such as SEIR, SEIRS, SIRS, SEI, SEIS, SI, SIS, SEIQRS, SIDARTHE, and others in future work.



---

*Source: 2138165-2022-09-26.xml*