Design and Analysis of a New COVID-19 Model with Comparative Study of Control Strategies

Abstract

The COVID-19 pandemic has become a worldwide concern and has caused great frustration in the human community. Governments all over the world are struggling to combat the disease. In an effort to understand and address the situation, we conduct a thorough study of a COVID-19 model that provides insights into the dynamics of the disease. For this, we propose a new $L_S{H}_{S}EAIHR$ COVID-19 model, where susceptible populations are divided into two sub-classes: low-risk susceptible populations, $L_S$, and high-risk susceptible populations, $H_S$. The aim of the subdivision of susceptible populations is to construct a model that is more reliable and realistic for disease control. We first prove the existence of a unique solution to the purposed model with the help of fundamental theorems of functional analysis and show that the solution lies in an invariant region. We compute the basic reproduction number and describe constraints that ensure the local and global asymptotic stability at equilibrium points. A sensitivity analysis is also carried out to identify the model’s most influential parameters. Next, as a disease transmission control technique, a class of isolation is added to the intended $L_S{H}_{S}EAIHR$ model. We suggest simple fixed controls through the adjustment of quarantine rates as a first control technique. To reduce the spread of COVID-19 as well as to minimize the cost functional, we constitute an optimal control problem and develop necessary conditions using Pontryagin’s maximum principle. Finally, numerical simulations with and without controls are presented to demonstrate the efficiency and efficacy of the optimal control approach. The optimal control approach is also compared with an approach where the state model is solved numerically with different time-independent controls. The numerical results, which exhibit dynamical behavior of the COVID-19 system under the influence of various parameters, suggest that the implemented strategies, particularly the quarantine of infectious individuals, are effective in significantly reducing the number of infected individuals and achieving herd immunity.

1. Introduction

COVID-19, a viral respiratory infection that was first reported in December 2019, is still a life threatening infection and has affected almost every part of the world [1,2,3]. Despite the vaccination being distributed, it will be necessary to continue taking preventive measures such as practicing social distancing and wearing masks until the prevalence of COVID-19 in a given area decreases to a safe level. This study deals with the mathematical analysis of the effectiveness of these preventive measures and strategies in minimizing the spread of infection. The virus causes a variety of symptoms. The most prevalent that may appear in affected individuals are a sore throat, loss of smell and taste, nausea, vomiting, diarrhea, and breathing difficulties [4,5]. COVID-19 was originally discovered in the city of Wuhan, China, where its origin was identified to be the Huanan seafood market, a live animal and seafood market. On 11 March 2020, the World Health Organization declared it a global pandemic [6]. The time to the onset of signs and symptoms, known as the incubation period, is about 2 to 14 days. Although during the incubation period, affected individuals may not suffer any illness, they can still transmit infection to other individuals. Persons with co-morbidities, including hypertension, diabetes, immune compromisation, respiratory infections, and age greater than 60 years, are at a high risk of developing a COVID-19 infection [7].

Mathematical modeling plays a crucial role in solving real-world problems across various fields, such as engineering, economics, physics, and biology. It is an essential tool in developing hybrid models for partially known intracellular signaling pathways and chemical processes [8,9,10]. Hybrid modeling approaches that combine physical models with machine learning techniques can be particularly effective in capturing complex, nonlinear relationships between different variables in the process. Through a hybrid modeling approach, researchers and engineers can gain a deeper understanding of the complex systems and make more informed decisions about how to optimize and control them [11,12]. Mathematical models have also been proven to be an efficient and powerful tool in the design, analysis and control of plans for infectious diseases [13,14,15], epidemics [16], and pandemics, such as Ebola [17], SARS, and MERS [18,19]. Mathematical models with integer- and fractional-order derivatives [20,21,22,23] have long been used to aid people in understanding how disease patterns evolve in populations, such as in dengue, malaria, TB, and other illnesses [24,25,26,27]. Such models attempt to account for a variety of important elements in disease transmission, including the presence of a disease vector, the occurrence of relapses and reinfections, clinical and sub-clinical cases, the examination of the impact of low-cost interventions, and so on.

To reduce the infection rate in the human population, a strong mechanism that provides the most effective and preventive strategy among the available strategies is critical. The best way to prevent infection is to use the best control strategies to treat the illness. Educational campaigns, lock-downs, therapies, vaccines, quarantines, treatments, and isolation are examples of such strategies [4,6,28]. Even though individuals worldwide have received vaccinations, it is necessary to continue taking preventive measures, such as maintaining a safe distance from others and wearing masks, until the number of people who contract COVID-19 reduces to a safe level. The primary goal of creating a COVID-19 model is to understand disease behaviour over time and to identify possible mandatory control strategies. As a result, there is an urgent need to examine disease prevention measures and devise a reasonable control strategy to limit disease spread. This research includes an examination of the mathematical and computational analysis of the newly proposed COVID-19 model, aiming to explore preventive measures and suggest an ideal control approach to restrict the transmission of the infection.

One of the proposed methods to investigate control strategies in the simulation of system dynamics, which may be used in various compartmental models [29,30,31], is to use isolation as a variable in mathematical models of COVID-19 dynamics. Because social exposure is a potential cause of COVID-19 outbreaks, isolating affected persons can lower the likelihood of COVID-19 spreading in the future. In this study, we propose a new $L_S{H}_{S}EAIHR$ COVID-19 model where susceptible populations are divided in to two sub-classes: low-risk susceptible populations, $L_S$, and high-risk susceptible populations, $H_S$. Based on our current knowledge of COVID-19, certain groups of people are at a higher risk of contracting the virus. These individuals are typically referred to as the high-risk susceptible population and include healthcare workers, providers (including all front-line workers), relatives of infected individuals, and individuals involved in burial processes. Conversely, the rest of the population is considered to be at a low risk of acquiring COVID-19 [32,33]. It is worth noting that an increase in the high-risk population can lead to a higher likelihood of COVID-19 transmission. Thus, the division of susceptible populations into low- and high-risk populations may help us study the disease transmission pattern more appropriately, and accordingly, a suitable disease control strategy maybe adopted. Regardless of the risk level, everyone should take appropriate precautions to prevent the spread of COVID-19. These measures include washing hands frequently, practicing social distancing, wearing masks in public, and becoming vaccinated when possible. Additionally, individuals who experience symptoms of COVID-19 should seek medical attention immediately.

In this article, we analyze the newly constructed model for the existence of a unique positive and bounded solution. To examine the validity of the model, we perform a stability analysis at the equilibrium points and establish the relevant important results. For disease control analysis, we reconstruct the proposed model to include an isolation compartment as a control strategy. First, we perform different simulations to observe the effect of different quarantine levels of infectious individuals on disease spread. Later, we formulate an optimal control problem to determine the best time-dependent quarantine rates that may be suggested to restrict the spread of disease. For this, we derive optimality conditions and design an algorithm to solve these conditions. For numerical solutions, we implement the RK4 method.

The rest of the paper is organized as follows. Section 2 discusses the formulation of a new COVID-19 model. Section 3 deals with the basic properties of the solution, e.g., the existence of a unique positive bounded solution. In Section 4, we determine the equilibrium points as well as the reproduction number. The local and global asymptotic stabilities of the model at equilibrium points are elaborated in Section 5. In Section 6, a comprehensive sensitivity analysis of the reproduction number to parameters is carried out. We investigate the effect of different time-independent isolation levels on disease spread in Section 7. An optimal control analysis for disease control is also part of this section. Finally, Section 8 summarizes the conclusion of the manuscript.

2. Model Formulation

A mathematical model can aid in the knowledge of the pathways for the transmission of a disease as well as illness prevention through the application of preventative measures [4,6]. We formulate a mathematical model of COVID-19 transmission to explore its dynamics. To mathematically describe the model, the whole population is segregated into seven classes: susceptible at low risk $L_S$, susceptible at high risk $H_S$, exposed $E\left(t\right)$, infected $I\left(t\right)$, asymptomatic $A\left(t\right)$, hospitalized $H\left(t\right)$, and recovered $R\left(t\right)$. Individuals with a high risk of infection are known as high-risk susceptible individuals, $H_S\left(t\right)$. This category includes clinicians, family members, and paramedical employees who deal with infected patients. Outside of this group, the population is considered susceptible at a low risk level, $L_S\left(t\right)$. Thus, at any time t, the entire population $N\left(t\right)$ can be represented through a mathematical equation as:

$$N\left(t\right)=H_S\left(t\right)+L_S\left(t\right)+E\left(t\right)+I\left(t\right)+A\left(t\right)+H\left(t\right)+R\left(t\right),$$

(1)

where the compartmental variables $H_S\left(t\right),L_S\left(t\right),E\left(t\right),I\left(t\right),A\left(t\right),H\left(t\right),$ and $R\left(t\right)$ are considered continuously differentiable functions of time $t\in [0,\infty )$. Figure 1 shows the pattern of disease flow through these compartments.

A mathematical model describing the disease flow pattern shown in Figure 1 is given by the following system of non-linear ordinary differential equations:

$$\begin{array}{cc}{\displaystyle \frac{d{L}_S}{dt}}=\hfill & \Pi (1-\omega )-{\displaystyle \frac{L_S({\alpha }_{1}A+{\alpha }_{2}I+{\alpha }_{3}H)}{N}}-\psi L_S,\hfill \end{array}$$

(2a)

$$\begin{array}{cc}{\displaystyle \frac{d{H}_S}{dt}}=\hfill & \Pi \omega -{\displaystyle \frac{H_S({\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}H)}{N}}-\psi H_S,\hfill \end{array}$$

(2b)

$$\begin{array}{cc}{\displaystyle \frac{dE}{dt}}=\hfill & {\displaystyle \frac{L_S({\alpha }_{1}A+{\alpha }_{2}I+{\alpha }_{3}H)}{N}}+{\displaystyle \frac{H_S({\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}H)}{N}}-(\psi +ϵ)E,\hfill \end{array}$$

(2c)

$$\begin{array}{cc}{\displaystyle \frac{dA}{dt}}=\hfill & ϵPE-({\sigma }_1+\psi )A,\hfill \end{array}$$

(2d)

$$\begin{array}{cc}{\displaystyle \frac{dI}{dt}}=\hfill & ϵ(1-P)E-({\alpha }_4+{\sigma }_2+{\delta }_I+\psi )I,\hfill \end{array}$$

(2e)

$$\begin{array}{cc}{\displaystyle \frac{dH}{dt}}=\hfill & {\alpha }_{4}I-({\sigma }_3+{\delta }_H+\psi )H,\hfill \end{array}$$

(2f)

$$\begin{array}{cc}{\displaystyle \frac{dR}{dt}}=\hfill & {\sigma }_{2}I+{\sigma }_{3}H+{\sigma }_{1}A-\psi R,\hfill \end{array}$$

(2g)

subject to the non-negative initial conditions:

$$\begin{array}{c}\begin{array}{cc}\hfill H_S\left(0\right)=& {\left(H_S\right)}_0,L_S\left(0\right)={\left(L_S\right)}_0,E\left(0\right)=E_0,A\left(0\right)=A_0,I\left(0\right)=I_0,\hfill \\ \hfill H\left(0\right)=& H_0,R\left(0\right)=R_0.\hfill \end{array}\hfill \end{array}$$

(2h)

A detailed description of the parameters of the proposed model (2) along with their numerical values are shown in

Table 1. Parametric description and the corresponding values.

Parameter	Description	Value	Source
$\Pi$	Recruitment rate	$0.5$	[7]
${\alpha }_1$	Contact rate of $L_S$ with A	$0.16$	[4,6]
${\gamma }_1$	Contact rate of $H_S$ with A	$0.2$	Assumed
${\alpha }_2$	Contact rate of $L_S$ with I	$0.25$	[30]
${\gamma }_2$	Contact rate of of $H_S$ with I	$0.3$	Estimated
${\alpha }_3$	Contact rate of $L_S$ with H	$0.1$	[30]
${\gamma }_3$	Contact rate of $H_S$ with H	$0.2$	Estimated
$ϵ$	Infectious rate	$0.1$	[31]
${\alpha }_4$	Transfer rate from I to H	$0.13266$	[30]
${\beta }_5$	Transfer rate from I to Q	$0.17$	[4,6]
${\beta }_6$	Transfer rate from A to Q	$0.18$	[4,6]
${\beta }_7$	Transfer rate from Q to H	$0.125$	[4,6]
${\sigma }_1$	Rate of recovery of A	$0.01$	Assumed
${\sigma }_2$	Rate of recovery of I	$0.03521$	[30]
${\sigma }_3$	Rate of recovery of H	$0.04255$	[30].
${\sigma }_4$	Rate of recovery of Q	$0.02$	Assumed
${\delta }_I$	Disease-induced mortality of I	$0.0079$	[31]
${\delta }_H$	Disease-induced mortality of H	$0.0068$	[30]
${\delta }_Q$	Disease-induced mortality of Q	$0.0068$	[31]
$\psi$	Natural death rate	$0.05$	Estimated

.

3. Physical Properties of the Model

In this section, we establish that the model (2) is epidemiologically meaningful in a feasible and bounded region. This can be established by proving that the model has a unique solution and that the solution set of the model is non-negative and bounded.

3.1. Existence and Uniqueness of the Solution:

We put the autonomous model (2) into the form:

$$\begin{array}{cc}{\displaystyle \frac{dy}{dt}}\hfill & =g\left(y\right(t\left)\right),\hfill \end{array}$$

(3a)

$$\begin{array}{cc}y\left(0\right)\hfill & =y_0,\hfill \end{array}$$

(3b)

where $y\left(t\right)\in C^1[0,\mathcal{T}]$ and $y\left(t\right):R_{+}\to R_{+}^7$ is defined by

$$y\left(t\right)={\left(L_S\left(t\right),H_S\left(t\right),E\left(t\right),I\left(t\right),A\left(t\right),H\left(t\right),R\left(t\right)\right)}^T,$$

with

$$y_0={\left(L_S\left(0\right),H_S\left(0\right),E\left(0\right),I\left(0\right),A\left(0\right),H\left(0\right),R\left(0\right)\right)}^T.$$

and

$$g\left(y\left(t\right)\right)={(g_1\left(y\left(t\right)\right),g_2\left(y\left(t\right)\right),g_3\left(y\left(t\right)\right),g_4\left(y\left(t\right)\right),g_5\left(y\left(t\right)\right),g_6\left(y\left(t\right)\right),g_7\left(y\left(t\right)\right))}^T,$$

is a vector-valued function from $R_{+}^7$ to $R_{+}^7$ with components

$$\begin{array}{cc}{g}_1\left(y\left(t\right)\right)\hfill & =\Pi (1-\omega )-{\displaystyle \frac{L_S({\alpha }_{1}A+{\alpha }_{2}I+{\alpha }_{3}H)}{N}}-\psi L_S,\hfill \\ g_2\left(y\left(t\right)\right)\hfill & =\Pi \omega -{\displaystyle \frac{H_S({\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}H)}{N}}-\psi H_S,\hfill \\ g_3\left(y\left(t\right)\right)\hfill & ={\displaystyle \frac{L_S({\alpha }_{1}A+{\alpha }_{2}I+{\alpha }_{3}H)}{N}}+{\displaystyle \frac{H_S({\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}H)}{N}}-(\psi +ϵ)E,\hfill \\ g_4\left(y\left(t\right)\right)\hfill & =ϵPE-({\sigma }_1+\psi )A,\hfill \\ g_5\left(y\left(t\right)\right)\hfill & =ϵ(1-P)E-({\alpha }_4+{\sigma }_2+{\delta }_I+\psi )I,\hfill \\ g_6\left(y\left(t\right)\right)\hfill & ={\alpha }_{4}I-({\sigma }_3+{\delta }_H+\psi )H,\hfill \\ g_7\left(y\left(t\right)\right)\hfill & ={\sigma }_{2}I+{\sigma }_{3}H+{\sigma }_{1}A-\psi R.\hfill \end{array}$$

Now, we state some fundamental theorems to prove the existence of a unique solution of the given model.

Theorem 1

([34]). Suppose $D=\left\{(t,z)\right|t\in R_{+},z\in R^n\}$, and let $h(t,z)$ be continuous on D and satisfy the Lipschitz condition there; then, the initial value problem

$$\begin{array}{c}{\displaystyle \frac{dz}{dt}}=h(t,z),z\left(t_0\right)=z_0.\hfill \end{array}$$

has a solution.

Theorem 2.

If the function g(y(t)) of (3) is continuously differentiable over $[0,\mathcal{T}]$, then it is Lipschitz-continuous.

Proof. 

Let $\mathbb{D}$ be a compact and convex subset of

$$\mathcal{V}=\left\{y\left(t\right)\right|t_0\le t\le \mathcal{T},y\left(t\right)\in R_{+}^7\}.$$

Let $y_1,y_2\in \mathbb{D}$; then, by mean value theorem, ∃ $\zeta \in (y_1,y_2)$ such that

$$\begin{array}{cc}{\displaystyle \frac{g\left(y_1\right)-g\left(y_2\right)}{y_1-y_2}}=\hfill & g^{{}^{\prime }}\left(\zeta \right),\hfill \end{array}$$

or

$$\begin{array}{cc}g\left(y_1\right)-g\left(y_2\right)=\hfill & g^{{}^{\prime }}\left(\zeta \right).(y_1-y_2).\hfill \end{array}$$

$$\begin{array}{cc}\mid g\left(y_1\right)-g\left(y_2\right)\mid \hfill & =\mid g^{{}^{\prime }}\left(\zeta \right).(y_1-y_2)\mid ,\hfill \\ \hfill & \le \parallel g^{{}^{\prime }}\left(\zeta \right)\parallel \parallel y_1-y_2\parallel .\hfill \end{array}$$

Since $g\in C^1[0,\mathcal{T}]$, over a convex compact set $\mathbb{D}$, there ∃ a constant $\tau >0$ such that

$$\begin{array}{c}\parallel g^{{}^{\prime }}\left(\zeta \right)\parallel \le \tau .\hfill \end{array}$$

Hence,

$$\begin{array}{cc}\mid g\left(y_1\right)-g\left(y_2\right)\mid \hfill & \le \tau \parallel \left(y_1-y_2\right)\parallel \hfill \\ \mid g\left(y_1\right)-g\left(y_2\right)\mid \hfill & \le \tau \parallel y_1-y_2\parallel \hfill \\ \parallel g\left(y_1\right)-g\left(y_2\right)\parallel \hfill & \le \tau \parallel y_1-y_2\parallel .\hfill \end{array}$$

Therefore, $g\left(y\right)$ is Lipschitz.    □

Since $g\left(y\right)$ given in (3) satisfies the Lipschitz conditions, Theorem 1 ensures the existence of a solution to System (3).

Theorem 3.

Suppose that the function $g\left(y\right)$ represents the state variables and satisfies the Lipschitz condition

$$\begin{array}{c}\parallel g\left(y_2\right)-g\left(y_1\right)\parallel \le \kappa \parallel y_2-y_1\parallel ;\hfill \end{array}$$

then, Problem (3) has a solution, and this solution is unique for

$$\begin{array}{c}k=\tau \mathcal{T}<1.\hfill \end{array}$$

Proof. 

By the fundamental theorem of calculus, the solution of an IVP (3) will be in the form:

$$\begin{array}{c}y\left(t\right)=y\left(0\right)+{\int }_0^{\mathcal{T}}g\left(y\right)dt.\hfill \end{array}$$

(4)

We will prove that the function $y\left(t\right)$ is a solution of (3) if and only if it satisfies the integral Equation (4). Let $y\left(t\right)$ be the solution of (3). Then, it satisfies Equation (4), i.e.,:

$$\begin{array}{c}y\left(t\right)=y\left(0\right)+{\int }_0^{\mathcal{T}}g\left(y\right)dt.\hfill \end{array}$$

(5)

For a converse implication, we let $y_n\left(t\right)$ be a sequence of solutions that converge to the solution of (3) with a successive iteration form, which can be defined as follows:

$$\begin{array}{c}{y}_n\left(t\right)=y_0+{\int }_0^{\mathcal{T}}g\left(y_{n-1}\left(t\right)\right)dt,i=1,2,.,.,n,\mathrm{with}{y}_0\left(t\right)=y\left(0\right).\hfill \end{array}$$

(6)

First of all, we show that the sequence (6) is contractive if $k=\tau T<1.$

$$\begin{array}{cc}\mid y_n\left(t\right)-y_{n-1}\left(t\right)\mid =\hfill & \mid {\int }_0^{\mathcal{T}}[g\left(y_{n-1}\left(t\right)\right)-g\left(y_{n-2}\left(t\right)\right)]dt\mid ,\hfill \\ \hfill & \le {\int }_0^{\mathcal{T}}\mid g\left(y_{n-1}\left(t\right)\right)-g\left(y_{n-2}\left(t\right)\right)\mid dt.\hfill \end{array}$$

Using the Lipchitz property of the function $g\left(y\right)$,

$$\begin{array}{cc}\mid y_n\left(t\right)-y_{n-1}\left(t\right)\mid \hfill & \le {\int }_0^{\mathcal{T}}\tau \mid y_{n-1}\left(t\right)-y_{n-2}\left(t\right)\mid dx,\hfill \\ \hfill & \le \tau \parallel y_{n-1}-y_{n-2}\parallel {\int }_0^{\mathcal{T}}dt,\hfill \\ \hfill & \le \tau \mathcal{T}\parallel y_{n-1}\left(t\right)-y_{n-2}\left(t\right)\parallel ,\hfill \\ \hfill \parallel y_n-y_{n-1}\parallel & \le \tau \mathcal{T}\parallel y_{n-1}-y_{n-2}\parallel ,\hfill \\ \hfill \parallel y_n-y_{n-1}\parallel & \le k\parallel y_{n-1}-y_{n-2}\parallel .\hfill \end{array}$$

This implies

$$d(y_n,y_{n-1})\le kd(y_{n-1},y_{n-2}).$$

(7)

Thus, the sequence (6) is contractive, and hence, it is a Cauchy sequence. Now, for $m,n\in N$ and $m>n$,

$$\begin{array}{cc}\parallel y_m-y_n\parallel =\hfill & \parallel y_m-y_{m-1}+y_{m-1}-y_{m-2}+y_{m-2}...-y_{n+1}+y_{n+1}-y_n+y_n-y_n\parallel ,\hfill \\ \hfill & \le \parallel y_m-y_{m-1}\parallel +\parallel y_m-1-y_{m-2}\parallel +...+\parallel y_{n+1}-y_n\parallel ,\hfill \\ & \le k^{m-1}\parallel y_1-y_0\parallel +k^{m-2}\parallel y_1-y_0\parallel +...+k^n\parallel y_1-y_0\parallel ,\hfill \\ \hfill & \le [k^{m-1}+k^{m-2}+...+k^n]\parallel y_1-y_0\parallel ,\hfill \end{array}$$

where

$$k=\tau \mathcal{T}<1.$$

Hence, the right-hand side is a geometric series which is always convergent for $\left|k\right|<1$.

$$\begin{array}{cc}\parallel y_m-y_n\parallel \le \hfill & k^n{\displaystyle \frac{1-k^{m-n}}{1-k}}\parallel y_1-y_0\parallel \le k^n{\displaystyle \frac{1}{1-k}}\parallel y_1-y_0\parallel .\hfill \end{array}$$

Since $0<k<1$, ${lim}_{n\to \infty }{k}^n=0$. Therefore, we infer that the sequence $\left(y_n\right)$ is Cauchy, and hence, it is convergent. Let ${lim}_{n\to \infty }{y}_n=y$; then, Equation (6) gives

$$\begin{array}{c}\underset{n\to \infty }{lim}{y}_n\left(t\right)=y\left(t\right)=y\left(0\right)+{\int }_0^{\mathcal{T}}g\left(y\left(t\right)\right)dt.\hfill \end{array}$$

(8)

Equation (8) is the required solution.

Uniqueness

To prove the uniqueness of the solution, we suppose, on the contrary, that the sequence $\left(y_n\right)$ converges to two different limits, $y^{{}^{\prime }}$ and $y^{{}^{″}}$. Then, there exist $n_{1}and{n}_2\in N$ such that

$$\begin{array}{c}\parallel y_n-y^{{}^{\prime }}\parallel <{ϵ}_1{n}_1\ge n,\hfill \\ \\ \parallel y_n-y^{{}^{″}}\parallel <{ϵ}_2{n}_2\ge n.\hfill \end{array}$$

Let $n^{*}=max\{n_1,n_2\}.$

$$\begin{array}{c}\parallel y^{{}^{\prime }}-y^{{}^{″}}\parallel =\parallel y^{{}^{\prime }}-y_n+y_n-y^{{}^{″}}\parallel \le \parallel y^{{}^{\prime }}-y_n\parallel +\parallel y_n-y^{{}^{″}}\parallel <{ϵ}_1+{ϵ}_2=ϵn^{*}\ge n,\hfill \end{array}$$

which implies

$$\begin{array}{c}\parallel y^{{}^{\prime }}-y^{{}^{″}}\parallel =0\Rightarrow y^{{}^{\prime }}=y^{{}^{″}}.\hfill \end{array}$$

Hence, we have proved that Solution (8) of Equation (3) exists and is unique.    □

3.2. Boundedness and Positivity

In order to examine the basic characteristics of the COVID-19 model, we shall demonstrate that the state variables $y=(L_S,H_S,E,I,A,H,R)$ are bounded and non-negative ∀ $t\ge 0$ in a feasible region $\Omega$ defined as

$$\Omega =\{L_S,H_S,E,I,A,H,R\in {\mathbb{R}}_{+}^7,0\le N\left(t\right)\le {\displaystyle \frac{\Pi }{\psi }}\}.$$

Theorem 4.

The variable $y\left(t\right)$ representing the states of the model (2) are bounded.

Proof. 

From Equation (1), we can write

$${\displaystyle \frac{dN}{dt}}={\displaystyle \frac{d{L}_S}{dt}}+{\displaystyle \frac{d{H}_S}{dt}}+{\displaystyle \frac{dE}{dt}}+{\displaystyle \frac{dA}{dt}}+{\displaystyle \frac{dI}{dt}}+{\displaystyle \frac{dH}{dt}}+{\displaystyle \frac{dR}{dt}}.$$

(9)

Using the equations of the model (2) in (9) and then simplifying them, we reach the following equation, which gives the rate of change for the entire population, i.e.,

$${\displaystyle \frac{dN}{dt}}=\Pi -{\delta }_{I}I-{\delta }_{H}H-\psi N.$$

(10)

Since ${\delta }_{I}I+{\delta }_{H}H\ge 0$, then from (10), we can write the inequality:

$${\displaystyle \frac{dN}{dt}}\le \Pi -\psi N.$$

(11)

Solving the inequality (11) by using grown wall’s inequality [6], we obtain

$$N\left(t\right)\le \frac{\Pi }{\psi }+\left(N\left(0\right)-\frac{\Pi }{\psi }\right)exp\left(-\psi t\right).$$

Consequently, for every initial condition

$$N\left(0\right)\le \frac{\Pi }{\psi },$$

we can write

$$\underset{t\to \infty }{lim}N\left(t\right)\le \frac{\Pi }{\psi }.$$

As a result, we deduced that the feasible region $\Omega$ is bounded for the state variables $(H_S,L_S,E,I,A,H,R)$.    □

Theorem 5.

Considering the initial conditions (2h), the solution $y\left(t\right)=(H_S,L_S,E,I,A,H,R)$ of model (2) is non-negative ∀$t\ge 0$.

Proof. 

Let us consider the Equation (2a) of model (2).

$${\displaystyle \frac{d{L}_S}{dt}}+\left(\mathrm{Y}\left(t\right)+\psi \right)L_S=\Pi (1-\omega ).$$

(12)

After being multiplied by the integrating factor $exp\left(\psi t+{\int }_0^t\mathrm{Y}\left(x\right)dx\right)$, Equation (12) can be written as

$${\displaystyle \frac{d}{dt}}\left[exp\left(\psi t+{\int }_0^t\mathrm{Y}\left(x\right)dx\right)L_S\right]=(1-\omega )\Pi exp\left(\psi t+{\int }_0^t\mathrm{Y}\left(x\right)dx\right),$$

which gives us

$$L_S\left(t_1\right)exp\left(\psi t_1+{\int }_0^{t_1}\mathrm{Y}\left(x\right)dx\right)=L_S\left(0\right)+\Pi (1-\omega ){\int }_0^{t_1}exp\left\{\psi t+{\int }_0^t\mathrm{Y}\left(x\right)dx\right\}dt.$$

Finding solutions for $L_S\left(t\right)$, it is clear that $L_S\left(t_1\right)>0forallt>0$. Using a similar method, we can demonstrate that $y\left(t\right)\ge 0$, $\forall t\ge 0$.    □

4. Equilibrium Points and Reproduction Number

4.1. Equilibrium Points

The equilibrium points of the model (2) are obtained by setting ${\displaystyle \frac{d{L}_s}{dt}}=0,{\displaystyle \frac{d{H}_s}{dt}}=0,{\displaystyle \frac{dE}{dt}}=0,{\displaystyle \frac{dI}{dt}}=0,{\displaystyle \frac{dA}{dt}}=0,{\displaystyle \frac{dH}{dt}}=0,$ and ${\displaystyle \frac{dR}{dt}}=0$ and then solving the resulting equations for the cases when the disease is present in the community and when the disease is absent in the community. The disease-free equilibrium point for model (2) is:

$$P_0=(L_S^{*},E^{*},H_S^{*},I^{*},A^{*},H^{*},R^{*})=\left({\displaystyle \frac{(1-\omega )\Pi }{\psi }},{\displaystyle \frac{\omega \Pi }{\psi }},0,0,0,0,0\right).$$

(13)

The endemic equilibrium is given as:

$$P_1=({\overline{L}}_S,{\overline{H}}_S,\overline{E},\overline{A},\overline{I},\overline{H},\overline{R}),$$

where

$$\begin{array}{cc}{\overline{L}}_S\hfill & ={\displaystyle \frac{(1-\omega )\Pi }{s_1+\psi }},{\overline{H}}_S={\displaystyle \frac{\omega \Pi }{s_2+\psi }},\overline{A}={\displaystyle \frac{Pϵ\overline{E}}{k_2}},\overline{I}={\displaystyle \frac{(1-P)ϵ\overline{E}}{k_3}},\overline{H}=\frac{{\alpha }_4(1-P)ϵ\overline{E}}{k_3{k}_4},\hfill \\ \overline{R}\hfill & ={\displaystyle \frac{1}{\psi }}[{\sigma }_1\overline{A}+{\sigma }_2\overline{I}+{\sigma }_3\overline{H}],\hfill \end{array}$$

and

$$s_1{s}_2{\overline{E}}^2+\left[k{s}_1{s}_2-\psi (s_1+s_2)\right]\overline{E}-k\psi (s_1-s_2)-{\psi }^2+k{s}_1\psi =0,$$

(14)

where

$$\begin{array}{cc}k\hfill & ={\displaystyle \frac{\Pi }{k_1}},s_1={\displaystyle \frac{k_3{k}_4{\alpha }_{1}Pϵ+k_2{k}_4{\alpha }_2(1-P)ϵ+k_2{k}_3{\alpha }_3{\alpha }_4(1-P)ϵ}{k_2{k}_3{k}_{4}N}},\hfill \\ s_2\hfill & ={\displaystyle \frac{k_3{k}_4{\gamma }_{1}Pϵ+k_2{k}_4{\gamma }_2(1-P)ϵ+k_2{\gamma }_3{\alpha }_4(1-P)ϵ}{k_2{k}_3{k}_{4}N}}.\hfill \\ k_1\hfill & =\psi +ϵ,k_2=\psi +{\sigma }_1,k_3={\alpha }_4+\psi +{\sigma }_2+{\delta }_I,k_4=\psi +{\sigma }_3+{\delta }_H.\hfill \end{array}$$

Equation (14) can be resolved for $\overline{E}$, and hence, the endemic point $P_1$ can be determined.

4.2. Reproduction Number

The reproduction number ${\mathcal{R}}_0$ is a criterion that estimates disease transmission in a population. It is the average proportion of individuals infected by an infectious agent in a susceptible population. The next-generation matrix approach established by Diekmann and Heesterbeek in 1990 [35] is used to compute ${\mathcal{R}}_0$. If ${\mathcal{R}}_0<1$, the disease will be eradicated. The disease will spread across the population if ${\mathcal{R}}_0>1$.

The recruitment of new infection $\mathbb{F}$ and transmission rates inside and outside in infected compartments $\mathbb{V}$ are given as:

$$\mathbb{F}=\left(\begin{array}{c}{\displaystyle \frac{L_S({\alpha }_{1}A+{\alpha }_{2}I+{\alpha }_{3}H)}{N}}+{\displaystyle \frac{H_S({\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}H)}{N}}\\ 0\\ 0\\ 0\end{array}\right),$$

$$\mathbb{V}=\left(\begin{array}{c}(\psi +ϵ)E\\ -ϵPE+({\sigma }_1+\psi )A\\ -ϵ(1-P)E+({\alpha }_4+{\sigma }_2+{\delta }_I+\psi )I\\ -{\alpha }_{4}I+({\sigma }_3+{\delta }_H+\psi )H\end{array}\right).$$

The Jacobians of $\mathbb{F}$ and $\mathbb{V}$ are computed at $P_0$ as

$$F={\left(\begin{array}{c}{\displaystyle \frac{\partial {\mathbb{F}}_j}{\partial {\mathbb{X}}_i}}\end{array}\right)}_{P_0},i,j=1,2,3,4.$$

$$V={\left(\begin{array}{c}{\displaystyle \frac{\partial {\mathbb{V}}_j}{\partial {\mathbb{X}}_i}}\end{array}\right)}_{P_0},i,j=1,2,3,4,$$

where $({\mathbb{X}}_1,{\mathbb{X}}_2,{\mathbb{X}}_3,{\mathbb{X}}_4)=(E,A,I,H)$. Thus, the resulting F and V are matrices that are given as

$$\begin{array}{cc}F\hfill & =\left(\begin{array}{cccc}0& (1-\omega ){\alpha }_1+\omega {\gamma }_1& (1-\omega ){\alpha }_2+\omega {\gamma }_2& (1-\omega ){\alpha }_3+\omega {\gamma }_3\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right).\hfill \\ V\hfill & =\left(\begin{array}{cccccc}{k}_1& 0& & 0& & 0\\ -Pϵ& k_2& & 0& & 0\\ -(1-P)ϵ& 0& & k_3& & 0\\ 0& 0& & -{\alpha }_4& & k_4\end{array}\right).\hfill \end{array}$$

Since ${\mathcal{R}}_0=$ the spectral radius of $F{V}^{-1}$,

Thus,

$$\begin{array}{cc}{\mathcal{R}}_0=\hfill & {\displaystyle \frac{k_3{k}_4ϵP\left((1-\omega ){\alpha }_1+\omega {\gamma }_1\right)+k_2{k}_4ϵ(1-P)\left((1-\omega ){\alpha }_2\right)+\omega {\gamma }_2}{k_1{k}_3{k}_4}}\hfill \\ & +{\displaystyle \frac{k_2(1-P){\alpha }_4ϵ\left((1-\omega ){\alpha }_3+\omega {\gamma }_3\right)}{k_1{k}_3{k}_4}}.\hfill \end{array}$$

5. Stability Analysis

In this part, we exploit ${\mathcal{R}}_0$ to discuss the local stability and global stability of the model (2) at DFE and EE. To investigate and discuss global stabilities, we implement theories by Lyapunov and Castillo-Chavez [36,37,38,39].

5.1. Local Stability

This section investigates the local stability of the COVID-19 model (2) at the DFE point $P_0$.

Theorem 6.

If ${\mathcal{R}}_0<1$, the model (2) is locally asymptotically stable at $P_0$, and it is unstable if ${\mathcal{R}}_0>1$.

Proof. 

The model (2) at the disease-free equilibrium $P_0$ corresponds to the following Jacobian matrix:

$$J\left(P_0\right)=\left(\begin{array}{ccccccccc}-\psi & & & 0& 0& -(1-\omega ){\alpha }_1& -(1-\omega ){\alpha }_2& -(1-\omega ){\alpha }_3& 0\\ 0& & & -\psi & 0& -\omega {\gamma }_1& -\omega {\gamma }_2& -\omega {\gamma }_3& 0\\ 0& & & 0& -k_1& q_1& q_2& q_3& 0\\ 0& & & 0& ϵP& -k_2& 0& 0& 0\\ 0& & & 0& ϵ(1-P)& 0& -k_3& 0& 0\\ 0& & & 0& 0& 0& {\alpha }_4& -k_4& 0\\ 0& & & 0& 0& {\sigma }_1& {\sigma }_2& {\sigma }_3& -\psi \end{array}\right).$$

We determine the following eigenvalues of the Jacobian matrix $J\left(P_0\right)$,

$$\begin{array}{cc}{\xi }_1\hfill & =-\psi ,\hfill \\ {\xi }_2\hfill & =-\psi ,\hfill \\ {\xi }_3\hfill & =-\psi ,\hfill \\ {\xi }_4\hfill & =-{\displaystyle \frac{(1-{\mathcal{R}}_0)k_2{k}_3{k}_4}{{\xi }_7{\xi }_6}},\hfill \\ {\xi }_5\hfill & =-k_1,\hfill \\ {\xi }_6\hfill & =-k_2[(1-{\mathcal{R}}_0)+{\displaystyle \frac{k_4ϵ(1-P)q_2+ϵ(1-P)q_3}{k_1{k}_3{k}_4}}],\hfill \\ {\xi }_7\hfill & ={\displaystyle \frac{[(1-{\mathcal{R}}_0)+{\displaystyle \frac{(1-P)ϵq_3}{k_1{k}_3{k}_4}}]k_2{k}_3}{{\xi }_6}},\hfill \end{array}$$

where $q_1={\alpha }_1(1-\omega )+{\gamma }_1\omega ,q_2={\alpha }_2(1-\omega )+{\gamma }_2\omega ,$ and $q_3={\alpha }_3(1-\omega )+{\gamma }_3\omega$.

All of the eigenvalues ${\xi }_i$ for $i=1,2,\dots ,7,$ are less than zero for ${\mathcal{R}}_0<1$. Therefore, at the disease-free equilibrium point $P_0$, the system of equations (2) is locally asymptotically stable. If ${\mathcal{R}}_0>1$, the eigenvalue ${\xi }_4$ with a positive sign demonstrates the model’s local instability at $P_0$.    □

5.2. Global Stability at DFE

At a DFE point, the Castillo-Chavez approach [40] is employed to ensure global stability. We replicate our model in the following equations by following the procedure established by Castillo-Chavez et al.

$$\begin{array}{c}\begin{array}{cc}{\displaystyle \frac{dX}{dt}}\hfill & =K(X,\varphi ),\hfill \\ {\displaystyle \frac{d\varphi }{dt}}\hfill & =G(X,\varphi ),G(X,0)=0,\hfill \end{array}\hfill \end{array}$$

(15)

where $X=(H_S,L_S)\in {\mathbb{R}}_{+}^2$ denotes the number of individuals who have not been infected, and $\varphi =(E,I,A,H)\in {\mathbb{R}}_{+}^4$ denotes the fraction of people who have been exposed, have been infected, are asymptomatic, or have been hospitalized. The last equation of model (2) is not taken into account because the others are independent of it. Here, $P_0=(X_0,0)$ is the disease-free equilibrium (DFE) point.

To verify the GAS of the DFE point in the Castillo-Chavez technique, the following two requirements must be met. That is

$$\begin{array}{cc}\left(M1\right)\mathrm{For}{\displaystyle \frac{dX}{dt}}\hfill & =K(X,0),X_0\mathrm{is}GAS,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\left(M2\right)G(X,\varphi )\hfill & =B\varphi -\overline{G}(X,\varphi ),\mathrm{where}\overline{G}(X,\varphi )\ge 0\mathrm{for}\mathrm{all}(X,\varphi )\in \Omega ,\hfill \end{array}$$

(17)

where $B=D_{\varphi }G(X_0,0)$ is an M-matrix, and $\Omega$ is the model’s feasible region. As a conclusion to Castillo-Chavez et al., we may state the following theorem.

Theorem 7.

If ${\mathcal{R}}_0<1$ and conditions (M1) and (M2) are fulfilled, then the DFE point $P_0$ of the proposed model is GAS.

Proof. 

Suppose $X=(H_S,L_S)$ stands in for uninfected individuals, $\varphi =(E,I,A,H)$ represents those who have been exposed to the disease, are infected, are asymptomatic, and are being hospitalized, and $P_0=(X_0,0)$ is the disease-free equilibrium. Then,

$$\begin{array}{c}K(X,\varphi )=\left[\begin{array}{c}\Pi (1-\omega )-\left({\alpha }_{1}A-{\alpha }_{2}I-{\alpha }_{3}H\right){\displaystyle \frac{L_S}{N}}-\psi L_S\\ \Pi \omega -\left({\gamma }_{1}A-{\gamma }_{2}I-{\gamma }_{3}H\right){\displaystyle \frac{H_S}{N}}-\psi H_S\end{array}\right].\hfill \end{array}$$

(18)

At $P_0=(X_0,0)$ and $K(X,0)$, we obtain

$$K(X,0)=\left[\begin{array}{c}(1-\omega )\Pi -\psi {L_S}_0\\ \omega \Pi -\psi {H_S}_0\end{array}\right]=0.$$

(19)

As $t\to \infty$, and $X\to X_0$. Therefore, $X=X_0=({L_S}_0,{H_S}_0)$ is globally asymptotically stable.

Now,

$$\begin{array}{cc}{\displaystyle \frac{d\varphi }{dt}}\hfill & =\left[\begin{array}{cccc}-(\psi +ϵ)& {\displaystyle \frac{{\alpha }_1{L_S}_0}{N}}+{\displaystyle \frac{{\gamma }_1{H_S}_0}{N}}& {\displaystyle \frac{{\alpha }_2{L_S}_0}{N}}+{\displaystyle \frac{{\gamma }_2{H_S}_0}{N}}& {\displaystyle \frac{{\alpha }_3{L_S}_0}{N}}+{\displaystyle \frac{{\gamma }_3{H_S}_0}{N}}\\ ϵP& -({\sigma }_1+\psi )& 0& 0\\ -(P-1)ϵ& 0& (-{\alpha }_4-{\sigma }_2-\psi -{\delta }_I)& 0\\ 0& 0& {\alpha }_4& -({\alpha }_4+{\sigma }_3+\psi +{\delta }_H)\end{array}\right]\left[\begin{array}{c}E\\ A\\ I\\ H\end{array}\right]\hfill \\ & -\left[\begin{array}{c}\left({\displaystyle \frac{{\alpha }_{1}A+{\alpha }_{2}I+{\alpha }_{3}H}{N}}\right)({L_S}_0-L_S)+\left({\displaystyle \frac{{\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}H}{N}}\right)({H_S}_0-H_S)\\ 0\\ 0\\ 0\end{array}\right],\hfill \end{array}$$

(20)

where

$$B=\left[\begin{array}{cccc}-(ϵ+\psi )& {\displaystyle \frac{{\alpha }_1{L_S}_0+{\gamma }_1{H_S}_0}{N}}& {\displaystyle \frac{{\alpha }_2{L_S}_0+{\gamma }_2{H_S}_0}{N}}& {\displaystyle \frac{{\alpha }_3{L_S}_0+{\gamma }_3{H_S}_0}{N}}\\ Pϵ& -({\sigma }_1+\psi )& 0& 0\\ (1-P)ϵ& 0& -({\alpha }_4+{\sigma }_2+\psi +{\delta }_I)& 0\\ 0& 0& {\alpha }_4& -({\alpha }_4+{\sigma }_3+\psi +{\delta }_H)\end{array}\right],\varphi =\left[\begin{array}{c}E\\ A\\ I\\ H\end{array}\right],$$

and

$$\overline{G}(X,\varphi )=\left[\begin{array}{c}\left({\displaystyle \frac{({\alpha }_{1}A+{\alpha }_{2}I+{\alpha }_{3}H)(L_{S_0}-L_S)}{N}}\right)+\left({\displaystyle \frac{({\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}H)({H_S}_0-H_S)}{N}}\right)\\ 0\\ 0\\ 0\end{array}\right].$$

B is obviously an M-matrix. Since $L_S\le L_{S_0}$ and $H_S\le H_{S_0}$, we have $\overline{G}(X,\varphi )\ge 0$. As a result, the DFE point $P_0$ is GAS.    □

5.3. Global Behavior at EE

Theorem 8.

The endemic equilibrium (EE) point denoted by $P_1$ is globally asymptotically stable if ${\mathcal{R}}_0>1$, and it is unstable if ${\mathcal{R}}_0<1$.

Proof. 

Assume that ${\mathcal{R}}_0>1$ so that the endemic equilibrium point exists. We now develop a Volterra-type Lyapunov functional $\Phi$, given as

$$\begin{array}{cc}\Phi \left(y\right)=\hfill & \left(L_S-\overline{L_S}-\overline{L_S}log\left({\displaystyle \frac{L_S}{\overline{L_S}}}\right)\right)-\left(-H_S+\overline{H_S}+\overline{H_S}log\left({\displaystyle \frac{H_S}{\overline{H_S}}}\right)\right)\hfill \\ & -\left(-E+\overline{E}+\overline{E}log\left({\displaystyle \frac{E}{\overline{E}}}\right)\right)-\left(-A+\overline{A}+\overline{A}log\left({\displaystyle \frac{A}{\overline{A}}}\right)\right)\hfill \\ \hfill & -\left(-I+\overline{I}+\overline{I}log\left({\displaystyle \frac{I}{\overline{I}}}\right)\right)-\left(-H+\overline{H}+\overline{H}log\left({\displaystyle \frac{H}{\overline{H}}}\right)\right)\hfill \\ \hfill & -\left(-R+\overline{R}+\overline{R}log\left({\displaystyle \frac{R}{\overline{R}}}\right)\right).\hfill \end{array}$$

(21)

Taking differentiation, inserting equations from model (2), and then rearranging them, we obtain

$$\begin{array}{cc}{\displaystyle \frac{d\Phi }{dt}}=\hfill & [(1-\omega )\Pi +(c1+\psi ){\displaystyle \frac{{\overline{L_S}}^2}{L_S}}+\omega \Pi +(c2+\psi ){\displaystyle \frac{{\overline{H_S}}^2}{H_S}}+c1L_S+c2H_S+k_1{\displaystyle \frac{{\overline{E}}^2}{E}}+PϵE\hfill \\ & +k2{\displaystyle \frac{{\overline{A}}^2}{A}}+(1-P)ϵE+k_3{\displaystyle \frac{{\overline{I}}^2}{I}}+{\alpha }_{4}I+k_4{\displaystyle \frac{{\overline{H}}^2}{H}}+{\sigma }_{1}A+\eta 2I+{\sigma }_{3}H+\psi {\displaystyle \frac{{\overline{R}}^2}{R}}]\hfill \\ & -[(c1+\psi ){\displaystyle \frac{{(L_S-\overline{L_S})}^2}{L_S}}+(c1+\psi )\overline{L_S}+(1-\omega )\Pi {\displaystyle \frac{\overline{L_S}}{L_S}}+(c2+\psi ){\displaystyle \frac{{(H_S-\overline{H_S})}^2}{H_S}}\hfill \\ & +(c2+\psi )\overline{H_S}+\omega \Pi {\displaystyle \frac{\overline{H_S}}{H_S}}+k_1{\displaystyle \frac{{(E-\overline{E})}^2}{E}}+(c1L_S+c2H_S){\displaystyle \frac{\overline{E}}{E}}+k_1\overline{E}+{\displaystyle \frac{{(A-\overline{A})}^2}{A}}{k}_2\hfill \\ & +PϵE{\displaystyle \frac{\overline{A}}{A}}+k_2\overline{A}+{\displaystyle \frac{{(I-\overline{I})}^2}{I}}{k}_3+(1-P)ϵE{\displaystyle \frac{\overline{I}}{I}}-k_3\overline{I}+{\displaystyle \frac{{(H-\overline{H})}^2}{H}}{k}_4+{\alpha }_{4}I{\displaystyle \frac{\overline{H}}{H}}+k_4\overline{H}\hfill \\ & +{\displaystyle \frac{{(R-\overline{R})}^2}{H}}\psi +\psi \overline{R}+{\displaystyle \frac{\overline{H}}{H}}({\sigma }_{1}A+\eta 2I+{\sigma }_{3}H)],\hfill \end{array}$$

where $c_1={\displaystyle \frac{{\alpha }_{1}A-{\alpha }_{2}I-{\alpha }_{3}H}{N}}{L}_S,$ and $c_2={\displaystyle \frac{{\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}H}{N}}{H}_S.$

Hence, this can be written as

$${\displaystyle \frac{d\Phi }{dt}}={\lambda }_1-{\lambda }_2,$$

where

$$\begin{array}{c}{\lambda }_1=(1-\omega )\Pi +(c1+\psi ){\displaystyle \frac{{\overline{L_S}}^2}{L_S}}+\omega \Pi +(c2+\psi ){\displaystyle \frac{{\overline{H_S}}^2}{H_S}}+c1L_S+c2H_S+k_1{\displaystyle \frac{{\overline{E}}^2}{E}}+PϵE\hfill \\ \hfill +k2{\displaystyle \frac{{\overline{A}}^2}{A}}+(1-P)ϵE+k_3{\displaystyle \frac{{\overline{I}}^2}{I}}+{\alpha }_{4}I+k_4{\displaystyle \frac{{\overline{H}}^2}{H}}+{\sigma }_{1}A+\eta 2I+{\sigma }_{3}H+\psi {\displaystyle \frac{{\overline{R}}^2}{R}},\end{array}$$

and

$$\begin{array}{cc}{\lambda }_2=\hfill & (c1+\psi ){\displaystyle \frac{{(L_S-\overline{L_S})}^2}{L_S}}+(c1+\psi )\overline{L_S}+(1-\omega )\Pi {\displaystyle \frac{\overline{L_S}}{L_S}}+(c2+\psi ){\displaystyle \frac{{(H_S-\overline{H_S})}^2}{H_S}}\hfill \\ & +(c2+\psi )\overline{H_S}+\omega \Pi {\displaystyle \frac{\overline{H_S}}{H_S}}+k_1{\displaystyle \frac{{(E-\overline{E})}^2}{E}}+(c1L_S+c2H_S){\displaystyle \frac{\overline{E}}{E}}+k_1\overline{E}+{\displaystyle \frac{{(A-\overline{A})}^2}{A}}{k}_2\hfill \\ & +PϵE{\displaystyle \frac{\overline{A}}{A}}+k_2\overline{A}+{\displaystyle \frac{{(I-\overline{I})}^2}{I}}{k}_3+(1-P)ϵE{\displaystyle \frac{\overline{I}}{I}}-k_3\overline{I}+{\displaystyle \frac{{(H-\overline{H})}^2}{H}}{k}_4+{\alpha }_{4}I{\displaystyle \frac{\overline{H}}{H}}+k_4\overline{H}\hfill \\ & +{\displaystyle \frac{{(R-\overline{R})}^2}{H}}\psi +\psi \overline{R}+{\displaystyle \frac{\overline{H}}{H}}({\sigma }_{1}A+\eta 2I+{\sigma }_{3}H).\hfill \end{array}$$

It is known that all the parameters used in the model are positives; therefore, ${\displaystyle \frac{d\Phi }{dt}}<0$ when ${\lambda }_1<{\lambda }_2$&${\displaystyle \frac{d\Phi }{dt}}=0$ if ${\lambda }_1={\lambda }_2$, which implies that $H_S=\overline{H_S}$, $L_S=\overline{L_S}$, $E=\overline{E}$, $I=\overline{I}$, $A=\overline{A}$, $H=\overline{H}$, and $R=\overline{R}$. Thus, the EE point $P_1$ is GAS according to LaSalle’s invariance principle.    □

6. Sensitivity Analysis

We investigated the sensitivity of the COVID-19 model to determine which parameter has a greater impact on ${\mathcal{R}}_0$. We calculated the sensitivity index with regards to each parameter of the reproduction number ${\mathcal{R}}_0$. The index has a greater value when the associated parameter is more sensitive to ${\mathcal{R}}_0$, meaning that ${\mathcal{R}}_0$ increases in the same proportion. The importance of the sensitivity index’s positive (or negative) sign is that ${\mathcal{R}}_0$ increases (or decreases) when the parameter increases. The sensitivity index is calculated utilizing the normalized forward sensitivity technique [41], which is defined as

$$S_{\delta }={\displaystyle \frac{\delta }{{\mathcal{R}}_0}}{\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial \delta }},$$

where $\delta$ is the parameter in ${\mathcal{R}}_0$ whose sensitivity must be determined. The calculated sensitivity indices for each of the parameters used in ${\mathcal{R}}_0$ are shown in

Table 2. Sensitivity indices for the reproduction number ${\mathcal{R}}_0$.

Parameter	Sensitivity Value/Index	Parameter	Sensitivity Value/Index
${\alpha }_1$	$0.5430$	$ϵ$	$0.1379$
${\gamma }_1$	$0.1579$	${\sigma }_1$	$-0.4124$
${\alpha }_2$	$0.1415$	${\sigma }_2$	$-0.04124$
${\gamma }_2$	$0.0424$	${\sigma }_3$	$-0.0485$
${\alpha }_3$	$0.0756$	${\delta }_I$	$-0.0104$
${\gamma }_3$	$0.0378$	${\delta }_H$	$-0.0078$
${\alpha }_4$	$-0.0163$	$\psi$	$-1.0418$

. The most sensitive parameter with a positive index is ${\alpha }_1$. Other parameters that have a greater effect on ${\mathcal{R}}_0$ than the rest of the parameters include ${\beta }_1$, ${\beta }_2$, ${\gamma }_1$, ${\sigma }_1$, and $ϵ$.

7. Disease Control Strategies

In this section, we study different strategies that may be implemented to control the disease. The COVID-19 outbreak may spread by direct exposure between infected and uninfected persons. Therefore, separating infected symptomatic and infected asymptomatic persons from susceptible populations is the simplest way to prevent COVID-19 from spreading. For this purpose, we updated the model (2) to add an isolation class, $Q\left(t\right)$, with some control rates.

We assume that infected (symptomatic) people are isolated at a rate of ${\beta }_5$. A rate of ${\beta }_6$ is also assumed for the isolation of asymptomatic infected individuals. Isolated individuals are hospitalised at a rate of ${\beta }_7$. Assume that ${\delta }_Q$, $\psi$ and ${\sigma }_4$ stand, respectively, for the death rate due to disease, the natural death rate, and the recovery rate of isolated individuals. The model (2) is modified based on the preceding assumptions to yield the following updated system.

$$\begin{array}{cc}{\displaystyle \frac{d{L}_S}{dt}}=\hfill & \Pi (1-\omega )-{\displaystyle \frac{({\alpha }_{1}A+{\alpha }_{2}I+{\alpha }_{3}H)L_S}{N}}-\psi L_S,\hfill \end{array}$$

(22a)

$$\begin{array}{cc}{\displaystyle \frac{d{H}_S}{dt}}=\hfill & \Pi \omega -{\displaystyle \frac{({\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}H)H_S}{N}}-\psi H_S,\hfill \end{array}$$

(22b)

$$\begin{array}{cc}{\displaystyle \frac{dE}{dt}}=\hfill & {\displaystyle \frac{({\alpha }_{1}A+{\alpha }_{2}I+{\alpha }_{3}H)L_S}{N}}+{\displaystyle \frac{({\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}H)H_S}{N}}-(ϵ+\psi )E,\hfill \end{array}$$

(22c)

$$\begin{array}{cc}{\displaystyle \frac{dA}{dt}}=\hfill & ϵPE-(\psi +{\sigma }_1+{\beta }_6)A,\hfill \end{array}$$

(22d)

$$\begin{array}{cc}{\displaystyle \frac{dI}{dt}}=\hfill & ϵ(1-P)E-({\alpha }_4+{\beta }_5+{\sigma }_2+\psi +{\delta }_I)I,\hfill \end{array}$$

(22e)

$$\begin{array}{cc}{\displaystyle \frac{dQ}{dt}}=\hfill & {\beta }_{6}A+{\beta }_{5}I-(\psi +{\beta }_7+{\sigma }_4+{\delta }_Q)Q,\hfill \end{array}$$

(22f)

$$\begin{array}{cc}{\displaystyle \frac{dH}{dt}}=\hfill & {\beta }_{7}Q+{\alpha }_{4}I-(\psi +{\sigma }_3+{\delta }_H)H,\hfill \end{array}$$

(22g)

$$\begin{array}{cc}{\displaystyle \frac{dR}{dt}}=\hfill & {\sigma }_{2}I+{\sigma }_{1}A+{\sigma }_{3}H+{\sigma }_{4}Q-\psi R,\hfill \end{array}$$

(22h)

along corresponding initial conditions given as:

$$\begin{array}{c}\begin{array}{cc}\hfill & {L_S}_0=L_S\left(0\right)\ge 0,{H_S}_0=H_S\left(0\right)\ge 0,E_0=E\left(0\right)\ge 0,A_0=A\left(0\right)\ge 0,\hfill \\ & I_0=I\left(0\right)\ge 0,Q_0=Q\left(0\right)\ge 0,H_0=H\left(0\right)\ge 0,R_0=R\left(0\right)\ge 0.\hfill \end{array}\hfill \end{array}$$

(22i)

The flow diagram of the modified model is shown in Figure 2. The model (22) will serve as a set of restrictions for the optimal control problem.

7.1. Effects of Different Isolation Levels

In this section, we adopt a qualitative approach to determine the effects of isolation on the disease dynamics governed by the proposed model (2). The RK-4 method was implemented to visualize these effects. We studied two different cases.

In the first case, we examined the effect of the isolation of asymptomatic infected individuals on disease control by considering different isolation levels, i.e., ${\beta }_6=0,0.1,0.3,0.5,$ and $0.7$. We can observe from Figure 3 and Figure 4 that the isolation of asymptomatic infected individuals has a great impact on controlling the pandemic COVID-19 disease. As we increase the isolation rate ${\beta }_6$ of asymptomatic individuals, solution curves for infected compartments move towards a disease-free state.

In the second case, we employed different quarantine rates of isolating infected (symptomatic) individuals to study the effect of this decision on disease control. We simulated the model (2) by considering ${\beta }_5=0,0.1,0.3,0.5,$ and $0.7$ and keeping all other parameters fixed. Figure 5 and Figure 6 give us a clear picture about the effects of the isolation of infected (symptomatic) individuals.

In both cases, we observe a decrease in the infected individuals with an increase in the isolation rates. However, a significant decrease in infected individuals is observed if asymptomatic infected individuals are isolated with higher quarantine levels. This study also reveals that ${\beta }_6$ may be considered a time-dependent control variable to control the spread of disease optimally.

7.2. Optimal Control Problem

Optimal control theory is an area of study concerned with identifying the best controls of a system over time in order to achieve a set of requirements. In the proposed model, we employ Pontryagin’s maximum principle [42] to derive the necessary conditions for identifying the possible optimal control. Pontryagin and Boltyansikii’s idea of optimal control has been employed in a variety of epidemic models [4,6,39] to forecast future outcomes. The goal of a control problem is to optimize an objective functional that consists of state and control variables.

To formulate an optimal control problem, the cost functional that takes into account the number of infected asymptomatic people A, the infected (symptomatic) people I, and the cost of executing control strategies is given as

$$\begin{array}{c}J(A,I,u)={\int }_0^{\mathcal{T}}\left[A+I+{\displaystyle \frac{1}{2}}{w}_1{u}_1^2+{\displaystyle \frac{1}{2}}{w}_2{u}_2^2+{\displaystyle \frac{1}{2}}{w}_3{u}_3^2\right]dt,\hfill \end{array}$$

(23)

where $w_1,w_2,w_3$ are the weights associated with each of the time-dependent controls, such that $u_1(={\beta }_5),u_2(={\beta }_6),$ and $u_3(={\beta }_7)$, and $\mathcal{T}$ is the final time. The goal is to determine the optimal controls $u^{*}\left(t\right)=(u_1^{*}\left(t\right),u_2^{*}\left(t\right),u_3^{*}\left(t\right))$ that minimize the objective functional (23). Consequently, the problem of optimal control is defined as follows:

$$\mathrm{Determine}\mathrm{minimizer}{u}^{*}\in U\mathrm{of}J(A,I,u)\mathrm{with}\mathrm{constraints}\left(22\right).$$

(24)

Here, U is a control space defined as

$$U=\{u:0\le u\le 1,0\le t\le \mathcal{T}\},$$

and the quarantine rates ${\beta }_5,{\beta }_6,$ and ${\beta }_7$ are, respectively, replaced with time-dependent controls $u_1\left(t\right),u_2\left(t\right),$ and $u_3\left(t\right)$ in an updated model (22).

To find the optimal control $u^{*}\left(t\right)$ for the optimal control problem (24), we first developed the necessary optimality conditions.

7.3. Necessary Conditions

In order to determine the conditions for the optimal control problem, Pontryagin’s maximum principle (PMP) was used. The conditions are derived from the Hamiltonian $\mathcal{H}$, which is defined as follows:

$$\begin{array}{c}\mathcal{H}(t,Z,u,\Lambda )=A+I+{\displaystyle \frac{1}{2}}{w}_1{u}_1^2+{\displaystyle \frac{1}{2}}{w}_2{u}_2^2+{\displaystyle \frac{1}{2}}{w}_3{u}_3^2+\sum _{n=1}^8{\Lambda }_n{g}_n(t,Z,u),\hfill \end{array}$$

where Z denotes the state variables, ${\Lambda }_n,n=1,\dots ,8,$ are the corresponding adjoint variables, and $g_n(t,Z,u),n=1,\dots ,8,$ are the right-hand sides of System (22).

The associated Hamiltonian $\mathcal{H}$ is given as follows:

$$\begin{array}{cc}\mathcal{H}(t,Z,u,\Lambda )=\hfill & A+I+{\displaystyle \frac{1}{2}}{w}_1{u}_1^2+{\displaystyle \frac{1}{2}}{w}_2{u}_2^2+{\displaystyle \frac{1}{2}}{w}_3{u}_3^2\hfill \\ & +{\Lambda }_1((1-\omega )\Pi -{\alpha }_1{\displaystyle \frac{A}{N}}{L}_S-{\alpha }_2{\displaystyle \frac{I}{N}}{L}_S-{\alpha }_3{\displaystyle \frac{H}{N}}{L}_S-\psi L_S)\hfill \\ & +{\Lambda }_2(\omega \Pi -{\gamma }_1{\displaystyle \frac{A}{N}}{H}_S-{\gamma }_2{\displaystyle \frac{I}{N}}{H}_S-{\gamma }_3{\displaystyle \frac{H}{N}}{H}_S-\psi H_S)\hfill \\ & +{\Lambda }_3({\alpha }_1{L}_S+{\gamma }_1{H}_S){\displaystyle \frac{A}{N}}-({\alpha }_2{L}_S+{\gamma }_2{H}_S){\displaystyle \frac{I}{N}}-({\alpha }_3{L}_S+{\gamma }_3{H}_S){\displaystyle \frac{H}{N}}-(\psi +ϵ)E)\hfill \\ & +{\Lambda }_4(PϵE-(\psi +u_2+{\sigma }_1)A)+{\Lambda }_5((1-P)ϵE-({\alpha }_4+u_1+\psi +{\delta }_I+{\sigma }_2)I)\hfill \\ \hfill & +{\Lambda }_6\left(u_{1}I+u_{2}A-(u_3+\psi +{\delta }_Q+{\sigma }_4)Q\right)+{\Lambda }_7({\alpha }_{4}I+u_{3}Q-(\psi +{\delta }_H+{\sigma }_3)H)\hfill \\ & +{\Lambda }_8({\sigma }_{1}A+{\sigma }_{2}I+{\sigma }_{3}H+{\sigma }_{4}Q-\psi R).\hfill \end{array}$$

(25)

The first optimality condition

$$\begin{array}{c}{\displaystyle \frac{\partial \mathcal{H}}{\partial u}}=0,\hfill \end{array}$$

of the PMP yields the following expressions for optimal control variables:

$$\begin{array}{c}{u}_1=\frac{({\Lambda }_5-{\Lambda }_6)I}{w_1},\hfill \\ u_2=\frac{({\Lambda }_4-{\Lambda }_6)A}{w_2},\hfill \\ u_3=\frac{({\Lambda }_6-{\Lambda }_7)Q}{w_3}.\hfill \end{array}$$

Under bounds, the controls are updated as follows:

$$\begin{array}{cc}{u}_1\hfill & =max\left[0,min\left(1,\frac{({\Lambda }_5-{\Lambda }_6)I}{w_1}\right)\right],\hfill \end{array}$$

(26a)

$$\begin{array}{cc}{u}_2\hfill & =max\left[0,min\left(1,\frac{({\Lambda }_4-{\Lambda }_6)A}{w_2}\right)\right],\hfill \end{array}$$

(26b)

$$\begin{array}{cc}{u}_3\hfill & =max\left[0,min\left(1,\frac{({\Lambda }_6-{\Lambda }_7)Q}{w_3}\right)\right].\hfill \end{array}$$

(26c)

The optimality condition

$$\begin{array}{c}{\displaystyle \frac{\partial \mathcal{H}}{\partial Z_n}}=-{\displaystyle \frac{d{\Lambda }_n}{dt}},n=1,\dots ,8,\hfill \end{array}$$

of PMP gives the following linear system of adjoint equations.

$$\begin{array}{cc}{\displaystyle \frac{d{\Lambda }_1}{dt}}=\hfill & \left({\displaystyle \frac{({\alpha }_{1}A+{\alpha }_{2}I+{\alpha }_{3}H)}{N}}+\psi \right){\Lambda }_1-\left({\displaystyle \frac{({\alpha }_{1}A+{\alpha }_{2}I+{\alpha }_{3}H)}{N}}\right){\Lambda }_3,\hfill \end{array}$$

(27a)

$$\begin{array}{cc}{\displaystyle \frac{d{\Lambda }_2}{dt}}=\hfill & \left({\displaystyle \frac{({\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}H)}{N}}+\psi \right){\Lambda }_2-\left({\displaystyle \frac{({\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}H)}{N}}\right){\Lambda }_3,\hfill \end{array}$$

(27b)

$$\begin{array}{cc}{\displaystyle \frac{d{\Lambda }_3}{dt}}=\hfill & \left(\psi ϵ\right){\Lambda }_3-ϵ(1-P){\Lambda }_5-ϵP{\Lambda }_4,\hfill \\ {\displaystyle \frac{d{\Lambda }_4}{dt}}=\hfill & \left({\sigma }_1+u_2+\psi \right){\Lambda }_4-{\sigma }_1{\Lambda }_8-u_2{\Lambda }_6+\left({\displaystyle \frac{{\alpha }_1{L}_S}{N}}\right){\Lambda }_1+\left({\displaystyle \frac{{\gamma }_1{H}_S}{N}}\right){\Lambda }_2\hfill \end{array}$$

(27c)

$$\begin{array}{cc}& -\left({\displaystyle \frac{{\alpha }_1{L}_S+{\gamma }_1{H}_S}{N}}\right){\Lambda }_3-1,\hfill \\ {\displaystyle \frac{d{\Lambda }_5}{dt}}=\hfill & \left({\alpha }_4+{\sigma }_2+u_1+\psi +{\delta }_I\right){\Lambda }_5-u_1{\Lambda }_6-{\alpha }_4{\Lambda }_7-{\sigma }_2{\Lambda }_8+\left({\displaystyle \frac{{\alpha }_2{L}_S}{N}}\right){\Lambda }_1\hfill \end{array}$$

(27d)

$$\begin{array}{cc}& +\left({\displaystyle \frac{{\gamma }_2{H}_S}{N}}\right){\Lambda }_2-\left({\displaystyle \frac{{\alpha }_2{L}_S+{\gamma }_2{H}_S}{N}}\right){\Lambda }_3-1,\hfill \end{array}$$

(27e)

$$\begin{array}{cc}{\displaystyle \frac{d{\Lambda }_6}{dt}}=\hfill & \left(u_3+\psi +{\delta }_Q+{\sigma }_4\right){\Lambda }_6-u_3{\Lambda }_7-{\sigma }_4{\Lambda }_8,\hfill \\ {\displaystyle \frac{d{\Lambda }_7}{dt}}=\hfill & \left({\sigma }_3+\psi +{\delta }_H\right){\Lambda }_7-{\sigma }_3{\Lambda }_8+\left({\displaystyle \frac{{\alpha }_3{L}_S}{N}}\right){\Lambda }_1+\left({\displaystyle \frac{{\gamma }_3{H}_S}{N}}\right){\Lambda }_2\hfill \end{array}$$

(27f)

$$\begin{array}{cc}& -\left({\displaystyle \frac{{\alpha }_3{L}_S+{\gamma }_3{H}_S}{N}}\right){\Lambda }_3,\hfill \end{array}$$

(27g)

$$\begin{array}{cc}{\displaystyle \frac{d{\Lambda }_8}{dt}}=\hfill & \psi {\Lambda }_8,\hfill \end{array}$$

(27h)

with conditions

$$\begin{array}{c}{\Lambda }_n\left(\mathcal{T}\right)=0,n=1,2,3,\dots ,8.\hfill \end{array}$$

(27i)

The derivatives of $\mathcal{H}$ with respect to adjoint variables ${\Lambda }_i,i=1,2,\dots ,8$ give us the state system (22).

7.4. Solution Algorithm

To solve optimality conditions for optimal solution, we follow steps of the Algorithm 1 [4] and use RK-4 method for numerical simulations.

Algorithm 1 Algorithm to find minimizer of the control problem (24)

Take $u_s\in U$ for $s=0$. Approximate the system (22) as well as the system (27) using control $u_s$. Determine $u^{*}$ using Equations (26). Reset control $u_s$ by computing $u_s=(u_s+u^{*})/2$. Stop the iterations when $\parallel {\theta }_s-{\theta }_{s-1}\parallel <\nu \parallel {\theta }_s\parallel$ for $s>0$,  otherwise $s+1⟵s$ and move to step 2.

Here $\theta$ is a representation for state variables, adjoint variables, control variables, and $\nu >0$ is the adjusted tolerance.

7.5. Optimal Solutions and Discussion

This section deals with the presentation and discussion of the optimal solutions to problem (24) with constant and time-dependent controls. MATLAB software was used to implement the RK-4 method to attain the optimal solution curves for the control problem (24). Simpson’s $\frac{1}{3}$ formula was employed to approximate the integral defined for the cost functional (23). We present four different cases with combinations of considered controls for the sake of analysis and to study the impact of controls on the eradication of diseases from society. A detailed discussion is presented in the following cases.

Case 1: In the first case, we considered two controls, $u_1$ and $u_2$, and studied their impact on reducing infection in compartments having infection. We took $u_3$ as a constant and ignored the term involving $u_3$ from the cost functional (23) by taking $w_3=0$.

The curves of optimal control variables and the corresponding objective functional are shown in Figure 7. We observe that, with each iteration, the curve for the cost functional reduces to attain its minimum. The state variables before and after optimization are shown in Figure 8. In this case, we find that the number of asymptomatic and infected individuals have decreased significantly, whereas the number of hospitalized patients increase initially and then decrease with time.

Case 2: Now, we take $u_1$ and $u_3$ as time-dependent controls and set $u_2$ as constant. We also take $w_2=0$ in the cost functional (23) to ignore the term involving $u_2$.

For this case, the optimal curves for control variables along with the corresponding objective functional are shown in Figure 9. It is observed that the optimal curve for the control $u_3$ does not change with time and stays almost zero for all time. Thus, the optimal hospitalization rate of quarantined individuals for the optimal solution is negligibly small. However, the optimal curve for $u_1$ changes with time and drops from its maximum value to its minimum value. The cost functional again reaches its smallest value under the effect of controls $u_1$ and $u_3$, as shown in Figure 9. Figure 10 represents solution curves for the state variables before and after optimization. We notice that the asymptomatic infected people decrease less as compared to the same individuals in case 1. Hospitalized individuals decrease after optimization in this case.

Case 3: In this case, we take $u_2$ and $u_3$ to be time-dependent controls and take $u_1$ as a constant. For this strategy, we take $w_1=0$ in the cost functional (23). Figure 11 shows curves for optimal control, which are the minimizers of the optimal control problem (24). The functional (23) reduces to its smallest value under the influence of these controls; see Figure 11. Figure 12 shows the graphs for curves of model variables before and after optimization. Compared to case 1, we see a relatively small reduction in the number of exposed and infected people, although this is better as compared to case 2.

Case 4: In the last case, we examine the impact of all three control variables $u_1,u_2,\mathrm{and}{u}_3$ together. Here, the combination of the isolation of asymptomatic and infected individuals and the hospitalization of isolated individuals is taken into account. The change in the intensity of control variables along with the cost functional over time are shown in Figure 13. The optimized hospitalization rate $u_3$ remains almost zero throughout the simulation, but the other rates fluctuate between the upper bound of 1 and the lower bound of zero. As shown in Figure 13, these control variables provide us the optimized isolation and hospitalization rates, reducing the cost functional to a minimum. Figure 14 shows the variations in state variables before and after optimization with respect to time. It is evident that the number of exposed, asymptomatic, infected, and hospitalized persons have decreased significantly under the implementation of all three control variables. Initially, the number of isolated persons increases, but then, it gradually decreases.

We can see that when all three control variables are considered together, the number of infectious persons decreases significantly compared to the scenarios when one of the control variables is considered time-independent. As a result, it is clear that Case 4 is the most efficient strategy for eliminating COVID-19 from the population.

8. Conclusions

In this study, a nonlinear COVID-19 model was examined for disease transmission prevention utilizing different control approaches. We developed a COVID-19 $L_S{H}_{S}EAIHR$ model in which the susceptible class was divided into two groups: a low-risk susceptible population, $L_S$, and a high-risk susceptible population, $H_S$. The existence and uniqueness of a solution is proved using the Picard iterative approach. We investigated the model’s fundamental characteristics and discovered that its solutions are bounded and non-negative. The reproduction number $R_0$ is evaluated, and the model’s local and global stability is examined at both of the equilibrium points. The DFE point is locally and globally asymptotically stable when $R_0<1$ and unstable when $R_0>1$. A sensitivity analysis of the fundamental reproduction number $R_0$ was also performed, revealing that the contact rate of hospitalized to susceptible individuals is the most sensitive parameter. We chose the isolation of infected individuals as the control strategy for limiting disease transmission. We not only examined the effect of different levels of time-independent controls on disease control but also developed an objective functional with the goal of decreasing asymptomatic and infected people optimally with specific time-dependent controls. Pontryagin’s maximum principle was employed to identify a solution to an optimal control problem. Moreover, the optimal control problem was numerically resolved with a comparison of the system without controls. Cases 1–4 were assessed with different control combinations, and we identified that activating all three controls at the same time reduced the number of infectious people significantly. We also observed a significant decrease in the asymptomatic infected and infected (symptomatic) curves while using a cost-effective approach with time-dependent controls.

To expand the scope of the research, vaccination can be utilized as an optimal control. In our forthcoming work, we aim to provide a precise depiction of COVID-19 disease by employing a fractional model with an ABC derivative operator and several intervention approaches. Additionally, we will utilize a fractional-order optimal control problem to determine the optimal methods for vaccination and quarantine.