An epidemiological model for analysing pandemic trends of novel coronavirus transmission with optimal control

ABSTRACT

Symptomatic and asymptomatic individuals play a significant role in the transmission dynamics of novel Coronaviruses. By considering the dynamical behaviour of symptomatic and asymptomatic individuals, this study examines the temporal dynamics and optimal control of Coronavirus disease propagation using an epidemiological model. Biologically and mathematically, the well-posed epidemic problem is examined, as well as the threshold quantity with parameter sensitivity. Model parameters are quantified and their relative impact on the disease is evaluated. Additionally, the steady states are investigated to determine the model's stability and bifurcation. Using the dynamics and parameters sensitivity, we then introduce optimal control strategies for the elimination of the disease. Using real disease data, numerical simulations and model validation are performed to support theoretical findings and show the effects of control strategies.

KEYWORDS:

• Epidemiological model
• SARS-CoV-2 virus
• backward bifurcation
• centre manifold theory
• steady states
• sensitivity
• optimal control
• numerical simulations

MATHEMATICS SUBJECT CLASSIFICATIONS:

• 34D08
• 34D20
• 34H05
• 37N25
• 37N35

1. Introduction

Infectious or communicable diseases are disorders produced by fungi and bacteria, viruses, or parasites. Various infections can be transferred from individual to individual, or from some contaminated water, food being exposed to organisms in the environment. Symptoms and signs are depending on the organisms causing the infection, but fatigue, fever, etc are common symptoms. Hepatitis A and B, Tuberculosis, flu, SARS, and influenza are some common examples. Coronaviruses are zoonotic diseases transmitted between animals and humans. The novel Coronavirus or SARS-CoV-2, initially, started in a seafood market in Wuhan, China. The pandemic initially began with a hospital admitting patients diagnosed with pneumonia. The sickness of patients was associated with the wet animal and seafood market [1]. The infection of SARS-CoV-2 virus transmits via multiple ways, e.g. due to human interaction via exposure to spoil surfaces directly and by respiratory droplets inhalation as well as from an infected person, however, the pandemic rises due to human interaction [2]. The various infection phases of the novel disease are asymptomatic and symptomatic. Especially, the asymptomatically infected individuals are crucial as they do not exhibit symptoms but transmit the infection. Therefore, a small amount of asymptomatically infected individuals will lead to a major disaster. Various pharmaceutical and non-pharmaceutical control strategies have been exercised for the control of the novel disease of the Coronavirus. For example, imposing mandatory and strict lockdown, maintaining social distancing, avoiding crowds, and masks use especially in public places, and imposing the ban on any type of gathering (social and religious) [3,4]. Similarly, the detected cases are placed in isolation for quick interventions [4]. Coronavirus has affected many countries around the world, including China, India, Canada, France, Germany, the United States, Pakistan, Oman, the United Arab Emirates, Africa, Southeast Asia, Europe, and the Western Pacific. As a result, many people suffered and died [5–7].

Mathematical modelling and optimal control of infectious disease is an emerging area that plays a significant role in exploring disease dynamics and predicting control measures [8–13]. The temporal dynamics of infectious diseases can be well understood with the aid of mathematical models. Usually, the models represent the interacting components arising in epidemiology and population biology [14–16]. These models can help to understand the dynamics and propose control strategies for disease eradication which are useful for health control authorities. With these guidelines, different steps are easily taken to control infection. Several models have been created to depict the time-based behaviour of contagious diseases. For example, a mathematical model has been presented to discuss the stability analysis with optimal vaccination using the concept of the classical SIR model in Ref. [17]. A model representing the dynamics of hepatitis B has been studied by Khan et al. [18]. Other mathematical models can be found in the literature for investigating the epidemiology of infectious diseases (see for detail, Refs [19–21]). Since, the epidemiology of infectious diseases is characterized by dynamic and evolving scenarios, so optimal control theory to design adaptive controls that can respond to change the conditions of disease transmissibility, the emergence of new variants and population behaviour. By continuously updating and optimizing control measures based on real-time data and feedback, optimal control theory enables a more proactive and agile approach to disease management [22,23].

The spreading and control of SARS-CoV-2 was difficult but in fact still a burning issue, and therefore the attention of various researchers has been attracted (see for more detail, Refs [24–33]). Since it is very much clear that the newly reported disease of SARS-CoV-2 virus has a complex nature. Especially, the different infection phases are very significant because the asymptomatically infected individuals have no symptoms and transmit the disease to others. It could be also noted that some of the individuals having no history of the asymptomatic stage while direct leads to the symptomatic stage. In this case, the probability-based transmission co-efficient is more reasonable because a susceptible population leads to asymptomatic and symptomatic at the same time. Moreover, vaccines for SARS-CoV-2 are now available and provide protection up to a certain degree of realism. Thus, we focus to develop an epidemiological model that contains the various infected groups with probability-based transmission and vaccination of susceptible in the form of an extraordinarily complex dynamical system to investigate the time dynamics and optimal control of the novel disease.

By assuming both symptomatic and asymptomatic individuals as sources of disease transmission, in this paper, we develop an epidemiological model to study the pandemic trend and optimal control of the novel disease. In more detail, the different populations are divided into four compartments, s, c, a, and r, which denote susceptible, symptomatically and asymptomatically infected, and recovered populations. When susceptible individuals interact successfully with infected individuals, they will enter both symptomatic and asymptomatic compartments. Due to the fact that susceptible individuals move to both infected compartments simultaneously, we use the probability-based transmission co-efficient. We then discuss the biological and mathematical feasibility of the considered epidemic problem. A detailed sensitivity analysis will be performed for every parameter in the model as well as the role that it plays in disease transmission and control, as well as the model equilibria will be examined. In addition, we analyse the bifurcation analysis of the model by using linear stability theory. The proposed model exhibits backward bifurcation based on the centre manifold theory. To eliminate the contagious infection of novel Coronavirus transmission, we designed an optimal control mechanism using optimal control theory, based on the temporal dynamics and parameter sensitivity of the model. Using real data of the novel disease, we validate the model and parameterize the parameters involved in its epidemic process. As a result, we use the estimated value for long-run numerical simulations to support our theoretical findings and show how the control measure affects disease transmission. This article is structured as follows. In Section 2, we discuss the detailed model formulation with positivity and boundedness. Section 3 shows how we derive the threshold quantity and calculate the equilibria in order to test the sensitivity of the model parameters and show how backward bifurcations occur. In Section 4, we proposed an optimal control mechanism for the presentation of Coronavirus disease. Section 5 presents the model calibration to support our theoretical results and the validity of the model. In the final Section 6, we provide a brief conclusion.

2. Formulation of the model

We present the model formulation to describe the time dynamics of the novel Coronavirus (SARS-CoV-2) disease. Since the symptomatically and asymptomatically infected individuals are very significant in the transmission of the novel disease of Coronavirus, therefore due to the characteristics of the novel disease, we distribute the entire community in the vulnerable or susceptible (those who are at risk), symptomatic, asymptomatic, and removed/recovered individuals, who are respectively characterized by $s(t),c(t),a(t)$, and $r(t)$. The susceptible individuals $s(t)$ at time t are those who are yet not infected but have a chance to catch the infection at time $t+\delta t$, where $\delta t$ is a small increment in t. The asymptomatic individuals $a(t)$ have no symptoms while becoming the source of the disease transmission. Similarly, the individuals showing symptoms of SARS-CoV-2 are known as symptomatic individuals, and those who become healthy after getting infected are the recovered population. We assume the newborn and immigrants from other compartments are susceptible because there is still no evidence of vertical transmission, although the immigrants may or may not be infected. Further, we assume a homogeneously mixed population, then the susceptible individuals enter the asymptomatic and symptomatic compartment after successful interaction with symptomatically and asymptomatically infected individuals. Two types of recovery for infected individuals are natural and based on treatment. The asymptomatic individuals enter the symptomatic compartment after producing symptoms of the SARS-CoV-2 virus. Moreover, the individuals leave the asymptomatic and symptomatic population group if they are fully recovered or die. Based on the effectiveness of the novel Coronavirus vaccination, those susceptible individuals which are successfully vaccinated lead to the recovered compartment. The complete transition of the novel disease is presented in the flowchart as presented in Figure 1.

Figure 1. The flowchart characterizes the transfer mechanism of the novel disease.

Therefore, we illustrate the time dynamics of disease as follows: (1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) with the initial compartmental population sizes (2) $s(0),r(0)>0,a(0),c(0)\ge 0.$(2) Here, the inflow of newborns is denoted by Π, and β is the disease transmission coefficient. It is assumed that γ represents the reduced rate of spreading, whereas the rate ν is indicative of the vaccination. Moreover, the natural as well as the SARS-CoV-2 related death rates are symbolized by ${\mu }_0$, ${\mu }_1$ and ${\mu }_2$, respectively. Similarly, ${\gamma }_2$ and ${\gamma }_1$ denote the recovery of asymptomatic (natural) and symptomatic individuals (due to treatment). Also, it is assumed that p is the probability of those who lead to the asymptomatic population group while q denotes the probability of those who recover from symptomatic directly.

We assume that $q_1=\nu +{\mu }_0$, $q_2={\mu }_0+{\mu }_1+{\gamma }_1$ and $q_3={\mu }_0+{\mu }_2+{\gamma }_2$ to shorten our calculations. We also assume that all parameters and constants used in the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) )–(2(2) $s(0),r(0)>0,a(0),c(0)\ge 0.$(2) ) have non-negative values because of the population dynamical problem. In the upcoming section, we will prove the well-posedness of the proposed epidemic problem.

2.1. Positivity and bounded-ness

To prove the well-posedness of the model, it is essential to investigate all the state variables having non-negative values.

Theorem 2.1

For all non-negative initial conditions, the solutions, $(s,a,c,r)\in {\mathcal{R}}_{+}^4$, of the proposed model are non-negative and uniformly bounded that is $n(t)\le \frac{\mathrm{\Pi }}{{\mu }_0}$, whenever the time grows without bounds. Further, the proposed epidemic problem is well-defined and positively invariant in the feasible biological region given by (3) $\mathrm{\Omega }=\{(s,a,c,r)\in {\mathcal{R}}_{+}^4:0\le s+a+c+r\le \frac{\mathrm{\Pi }}{{\mu }_0}\}.$(3)

Proof.

Let $(s,a,c,r)\in {\mathcal{R}}_{+}^4$ be solution of the epidemiological system (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ), then by integrating the first equation, we obtain (4) $\begin{array}{rl}s(t)& =s(0)\exp \{-{\int }_0^t\{\beta a(x)+\gamma \beta c(x)+q_1\}\mathrm{d}x\}\\ & +\mathrm{\Pi }\exp \{-{\int }_0^t\{\beta a(x)+\gamma \beta c(x)+q_1\}\mathrm{d}x\}\\ & \times {\int }_0^t\exp \left\{{\int }_0^x\{\beta a(y)+\gamma \beta c(y)+q_1\}\mathrm{d}y\right\}\mathrm{d}x.\end{array}$(4) According to the second equation in the epidemiological model (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ), it can be written as $\frac{\mathrm{d}a(t)}{\mathrm{d}t}\ge -\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),$which leads to (5) $a(t)\ge a(0)\exp \{-{\int }_0^t\{{\mu }_0+{\mu }_1+{\gamma }_1\}\mathrm{d}x\}\ge 0.$(5) In a similar fashion, the last two equations of the model look like $\frac{\mathrm{d}c(t)}{\mathrm{d}t}\ge -\{{\mu }_0+{\mu }_2+{\gamma }_2\}a(t),\frac{\mathrm{d}r(t)}{\mathrm{d}t}\ge -{\mu }_{0}r(t).$Integrating, we then obtain (6) $c(t)\ge c(0)\exp \{-{\int }_0^t\{{\mu }_0+{\mu }_2+{\gamma }_2\}\mathrm{d}x\}\ge 0,r(t)\ge r(0)\exp \{-{\mu }_{0}t\}\ge 0.$(6) It could be noted from Equations (4(4) $\begin{array}{rl}s(t)& =s(0)\exp \{-{\int }_0^t\{\beta a(x)+\gamma \beta c(x)+q_1\}\mathrm{d}x\}\\ & +\mathrm{\Pi }\exp \{-{\int }_0^t\{\beta a(x)+\gamma \beta c(x)+q_1\}\mathrm{d}x\}\\ & \times {\int }_0^t\exp \left\{{\int }_0^x\{\beta a(y)+\gamma \beta c(y)+q_1\}\mathrm{d}y\right\}\mathrm{d}x.\end{array}$(4) ), (5(5) $a(t)\ge a(0)\exp \{-{\int }_0^t\{{\mu }_0+{\mu }_1+{\gamma }_1\}\mathrm{d}x\}\ge 0.$(5) ) and (6(6) $c(t)\ge c(0)\exp \{-{\int }_0^t\{{\mu }_0+{\mu }_2+{\gamma }_2\}\mathrm{d}x\}\ge 0,r(t)\ge r(0)\exp \{-{\mu }_{0}t\}\ge 0.$(6) ) that the solutions $(s,a,c,r)$ of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) satisfying the initial conditions (2(2) $s(0),r(0)>0,a(0),c(0)\ge 0.$(2) ) remain non-negative for every $t\ge 0$.

To proceed further, let $n(t)$ denote the sum of the compartmental population, then $\frac{\mathrm{d}n(t)}{\mathrm{d}t}+{\mu }_{0}n(t)\le \mathrm{\Pi },$which yields that (7) $0<n(t)\le \frac{\mathrm{\Pi }}{{\mu }_0}\{1-e^{-{\mu }_{0}t}\}+n(0)e^{-{\mu }_{0}t}.$(7) From Equation (7(7) $0<n(t)\le \frac{\mathrm{\Pi }}{{\mu }_0}\{1-e^{-{\mu }_{0}t}\}+n(0)e^{-{\mu }_{0}t}.$(7) ), it is observed that $n(t)$ becomes less than or equal to $\frac{\mathrm{\Pi }}{{\mu }_0}$ whenever time grows. Thus it is follows that the total population in the region $\{{\mathcal{R}}_{+}^4\mathrm{\setminus }\{0\}\}$ is bounded by $\frac{\mathrm{\Pi }}{{\mu }_0}$ with growing time ($t\to \mathrm{\infty }$), therefore, the solution trajectories of the model satisfying the initial conditions are bounded.

Since the model state variables are non-negative and $n(t)\le \frac{\mathrm{\Pi }}{{\mu }_0}$ implies that the epidemic problem (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) )–(2(2) $s(0),r(0)>0,a(0),c(0)\ge 0.$(2) ) is well-posed. Further, we assume that $n(0)\le \frac{\mathrm{\Pi }}{{\mu }_0}$, then from Equation (7(7) $0<n(t)\le \frac{\mathrm{\Pi }}{{\mu }_0}\{1-e^{-{\mu }_{0}t}\}+n(0)e^{-{\mu }_{0}t}.$(7) ), we conclude that $n(t)\le \frac{\mathrm{\Pi }}{{\mu }_0}$, thus every solution of the proposed epidemic problem with initial conditions in ${\mathcal{R}}_{+}^4$ remains in Ω and is well-defined dynamical system and positively invariant.

3. Steady states analysis and reproductive number

Here, we discuss the steady states of the model and find the reproductive number. First, we calculate the SARS-CoV-2 virus-free equilibrium of the model and find the threshold parameter to perform the sensitivity analysis. Further, we use the threshold parameter to calculate the endemic state of the model and its detailed temporal dynamics.

3.1. SARS-CoV-2 virus-free equilibrium and reproductive number

The model reported by the system (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) is examined for equilibrium. We calculate the infection-free and disease-endemic states of the proposed model by assuming that ${\mathcal{E}}_0$ is the disease-free state as given by ${\mathcal{E}}_0=(s_0,0,0,r_0)$, where $s_0=\frac{\mathrm{\Pi }}{q_1}$, $r_0=\frac{\nu \mathrm{\Pi }}{{\mu }_0{q}_1}$. We calculate the reproductive quantity by following Driessche and Watmough [34] and Diekmann et al. [35], and accordingly used by Khajanchi et al. [36]. For this, let us assume that $X=(a,c{)}^T$, then $\frac{\mathrm{d}X}{\mathrm{d}t}{|}_{{\mathcal{E}}_0}=F-V,$where $F=\left[\begin{array}{cc}p\beta s_0& p\beta \gamma s_0\\ (1-p)p\beta s_0& (1-p)\beta \gamma s_0\end{array}\right],V=\left[\begin{array}{cc}{q}_2& 0\\ -{\gamma }_{1}q& q_3\end{array}\right].$Calculating the dominant eigenvalue ρ of the matrix $F{V}^{-1}$ implies (8) ${\mathcal{R}}_0=\frac{p\mathrm{\Pi }\beta }{q_1{q}_2}+\frac{\mathrm{\Pi }\gamma {\gamma }_1\beta pq}{q_1{q}_2{q}_3}+\frac{(1-p)\mathrm{\Pi }\gamma \beta }{q_1{q}_3}.$(8) For the shake of simplicity, let us assume that ${\mathcal{R}}_1=\frac{p\mathrm{\Pi }\beta }{q_1{q}_2}$, ${\mathcal{R}}_2=\frac{\mathrm{\Pi }\gamma {\gamma }_1\beta pq}{q_1{q}_2{q}_3}$ and ${\mathcal{R}}_3=\frac{(1-p)\mathrm{\Pi }\gamma \beta }{q_1{q}_3}$, then the threshold quantity becomes ${\mathcal{R}}_0={\mathcal{R}}_1+{\mathcal{R}}_2+{\mathcal{R}}_3$.

3.2. Sensitivity analysis of the reproductive number

We discuss the sensitivity analysis of the proposed model and its associated impact on the threshold parameter. The reproduction number plays an essential role to study the robustness of model prediction as the initial transmission of any epidemic is directly associated with this quantity. Because a small change in the value of the most sensitive parameter can result in a great effect on any epidemiological process [30]. The normalized sensitivity analysis helps to figure out the most sensitive parameter with respect to the threshold quantity. We carry out the normalized sensitivity to determine the most sensitive parameters about the persistence, reduction, and prevalence of SARS-CoV-2 virus transmission.

Definition 3.1

Let Φ be any epidemic parameter and ${\mathcal{R}}_0$ is the threshold quantity, the normalized forward sensitivity index of Φ is defined as (9) $S_{\mathrm{\Phi }}^{{\mathcal{R}}_0}=\frac{\mathrm{\Phi }}{{\mathcal{R}}_0}\frac{\mathrm{\partial }{\mathcal{R}}_0}{\mathrm{\partial }\mathrm{\Phi }}.$(9)

To calculate the association of the model parameters with the basic reproductive quantity and its relative effect, we calculate the following partial derivatives: $\begin{array}{rl}\frac{\mathrm{\partial }{\mathcal{R}}_0}{\mathrm{\partial }\beta }& =\frac{\mathrm{\Pi }p}{q_1{q}_2}+\frac{\mathrm{\Pi }{\gamma }_1\gamma qp}{q_3{q}_1{q}_2}+\frac{(1-p){\gamma }_0\mathrm{\Pi }}{q_1{q}_3},\frac{\mathrm{\partial }{\mathcal{R}}_0}{\mathrm{\partial }\gamma }=\frac{{\gamma }_1\mathrm{\Pi }\beta qp}{q_1{q}_2{q}_3}+\frac{(1-p)\beta \mathrm{\Pi }}{q_1{q}_3},\\ \frac{\mathrm{\partial }{\mathcal{R}}_0}{\mathrm{\partial }\nu }& =-\frac{p\beta \mathrm{\Pi }}{q_1^2{q}_2}-\frac{\gamma \mathrm{\Pi }{\gamma }_1\beta qp}{q_1^2{q}_2{q}_3},\frac{\mathrm{\partial }{\mathcal{R}}_0}{\mathrm{\partial }{\gamma }_2}=-\frac{qp\beta \mathrm{\Pi }\gamma {\gamma }_1}{q_1{q}_2{Q}_3^2}-\frac{(1-p)\mathrm{\Pi }\beta \gamma }{q_1{q}_3^2{q}_3}.\end{array}$It is evident from the above relationships that certain parameters exert a positive influence on the basic reproductive number, while others have a negative influence. The + and − signs of the indices respectively describe the direct and inverse relation of the model parameters with the reproductive quantity. If the parameters that directly affect the threshold quantity are decreased or increased, the value of the threshold quantity will also decrease or increase accordingly. But in contrast, the basic reproductive number will decrease or increase with the increase or decrease in the value of parameters which are in inverse relation. We calculate the sensitivity indices and their relative effect with the aid of the normalized sensitivity index formula given in the Equation (9(9) $S_{\mathrm{\Phi }}^{{\mathcal{R}}_0}=\frac{\mathrm{\Phi }}{{\mathcal{R}}_0}\frac{\mathrm{\partial }{\mathcal{R}}_0}{\mathrm{\partial }\mathrm{\Phi }}.$(9) ) to compute the parameters which are more sensitive to the novel disease of Coronavirus transmission and its prevention. In this case, we use the parameter values as: $\mathrm{\Pi }=158.43$, $\nu =0.9$, $\beta =0.079$, $\gamma =0.0457$, ${\gamma }_1=0.89$, ${\gamma }_2=0.0081$, p = 0.26, q = 0.79, ${\mu }_1=0.8$, ${\mu }_0=0.0000350$ and ${\mu }_2=0.63$. We accordingly calculate the indices with sensitivity as given in Table 1.

Table 1. Sensitivity indices and its relative impact on the basic reproductive number (${\mathcal{R}}_0$).

Parameters	Indices	% Increase or Decrease	Impact on ${\mathcal{R}}_0$
ν	$-0.999$	10	9.99 %
β	$+1.00$	10	10.0 %
γ	$+0.283$	10	2.83 %
${\gamma }_2$	$-0.003$	10	0.03 %
${\gamma }_1$	$-0.360$	10	3.60 %
p	$+0.036$	10	0.36 %
q	$+0.036$	10	0.36 %

The parameters which have in direct proportion with the basic reproductive number are $A_1=\{\beta ,\gamma ,q,p\}$, while the set of those parameters which have inverse relation is listed as $A_2=\{{\gamma }_2,{\gamma }_1,\nu \}$. We observe that an increase or decrease in the value of β, for example, by 10 % would increase or decrease the value of ${\mathcal{R}}_0$ by 10.0 %. In a similar fashion, if we perturb the value of γ, q, and p by 10 %, increase or decrease the value of ${\mathcal{R}}_0$ by 2.8 %, 0.36 % and 2.02 %, respectively. On the other hand, the sensitivity indices of the parameters ν, ${\gamma }_2$, and ${\gamma }_1$ reveal that if their values are increased by 10 % would collectively decrease ${\mathcal{R}}_0$ by 13.94  

The parameters that significantly affect reproduction number are β, γ, ${\gamma }_1$, and ν. In the case of the novel Coronavirus disease, the disease-spreading coefficient is the most influential epidemic parameter. Furthermore, the normalized sensitivity index of the vaccination rate indicates vaccination is crucial to reducing Coronavirus disease propagation. In addition, the parameters ${\gamma }_1$ and ${\gamma }_2$ are also essential for decreasing the reproductive number by increasing the recovery of the infected population. Consequently, we recommend control measures to slow the progress of β and γ and to strengthen the capacity of ν, ${\gamma }_1$, and ${\gamma }_2$.

3.3. Endemic state of the model

We use this the threshold quantity (${\mathcal{R}}_0$) to calculate the endemic state of the epidemiological model (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) by assuming ${\mathcal{E}}^{\ast }=(s^{\ast },a^{\ast },c^{\ast },r^{\ast })$, where (10) $\begin{array}{rl}{s}^{\ast }& =\frac{q_2{q}_3}{\beta \{q_{3}p+{\gamma }_{1}q\gamma p+q_2\gamma (1-p)\}},a^{\ast }=\frac{p{q}_1{q}_2{q}_3\{{\mathcal{R}}_0-1\}}{q_2\beta \{q_{3}p+{\gamma }_{1}q\gamma p+q_2\gamma (1-p)\}},\\ c^{\ast }& =\frac{q_1\{pq{\gamma }_1+q_2(1-p)\}\{{\mathcal{R}}_0-1\}}{\beta \{q_{3}p+{\gamma }_{1}q\gamma p+(1-p)q_2{\gamma }_1\}},r^{\ast }=\frac{\{{\gamma }_2{c}^{\ast }+a^{\ast }(1-q)+\nu s^{\ast }\}}{{\mu }_0}.\end{array}$(10) It can be observed clearly from the Equation (10(10) $\begin{array}{rl}{s}^{\ast }& =\frac{q_2{q}_3}{\beta \{q_{3}p+{\gamma }_{1}q\gamma p+q_2\gamma (1-p)\}},a^{\ast }=\frac{p{q}_1{q}_2{q}_3\{{\mathcal{R}}_0-1\}}{q_2\beta \{q_{3}p+{\gamma }_{1}q\gamma p+q_2\gamma (1-p)\}},\\ c^{\ast }& =\frac{q_1\{pq{\gamma }_1+q_2(1-p)\}\{{\mathcal{R}}_0-1\}}{\beta \{q_{3}p+{\gamma }_{1}q\gamma p+(1-p)q_2{\gamma }_1\}},r^{\ast }=\frac{\{{\gamma }_2{c}^{\ast }+a^{\ast }(1-q)+\nu s^{\ast }\}}{{\mu }_0}.\end{array}$(10) ) that the endemic states of the model exist whenever ${\mathcal{R}}_0>1$. Thus we state a lemma as follows.

Lemma 3.2

The proposed epidemiological model has disease endemic states whenever ${\mathcal{R}}_0>1$. Otherwise, the proposed epidemic problem has only a disease-free equilibrium.

3.4. Stability analysis and existence of bifurcation

We investigate the model stabilities and existence of bifurcation to show the dynamics of the epidemic problem that is under consideration. For this, we state the results in the following.

Theorem 3.3

The infection-free state (${\mathcal{E}}_0$) of the epidemiological model (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) is stable locally asymptotically, whenever ${\mathcal{R}}_0<1$.

Proof.

To determine the local dynamics of the epidemic problem at ${\mathcal{E}}_0$, we compute the Jacobian matrix $J{|}_{{\mathcal{E}}_0}$ of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) around infection-free equilibrium as (11) $J{|}_{{\mathcal{E}}_0}=\left(\begin{array}{cccc}-q_1& -\beta \frac{\mathrm{\Pi }}{q_1}& -\gamma \beta \frac{\mathrm{\Pi }}{q_1}& 0\\ 0& \frac{p\beta \mathrm{\Pi }-q_1{q}_2}{q_1}& \frac{p\beta \gamma \mathrm{\Pi }}{q_1}& 0\\ 0& \frac{(1-p)\beta \mathrm{\Pi }+q_{1}q{\gamma }_1}{q_1}& \frac{(1-p)\gamma \mathrm{\Pi }\beta -q_1{q}_3}{q_1}& 0\\ \nu & {\gamma }_1(1-q)& {\gamma }_2& -{\mu }_0\end{array}\right).$(11) Clearly, $J{|}_{{\mathcal{E}}_0}$ have two eigenvalues such that ${\lambda }_1=-{\mu }_0$ and ${\lambda }_2=-q_1$. For the rest of the eigenvalues we take the following matrix in reduced form $A=\left(\begin{array}{cc}\frac{p\beta \mathrm{\Pi }-q_1{q}_2}{q_1}& \frac{p\gamma \beta \mathrm{\Pi }}{q_1}\\ \frac{(1-p)\beta \mathrm{\Pi }+q_1{q}_2{\gamma }_1}{q_1}& \frac{(1-p)\gamma \mathrm{\Pi }\beta -q_1{q}_3}{q_1}\end{array}\right).$It could be noted from the matrix A, that $\mathrm{trace}(A)=-\frac{1}{q}\{q_1{q}_2(1-{\mathcal{R}}_1)+q_1{q}_3(1-{\mathcal{R}}_3)\},$and the determinant of the matrix (A) gives $\det (A)=q_2{q}_3\{(1-{\mathcal{R}}_1)(1-{\mathcal{R}}_3)\}+q_2{q}_3\{{\mathcal{R}}_1{\mathcal{R}}_3+{\mathcal{R}}_2\}.$Clearly, the $trace(A)<0$ and the $\det (A)>0$ whenever ${\mathcal{R}}_0<1$, which implies that Routh-Hurwitz criteria satisfied and the matrix A having negative eigenvalues or having negative real parts of the eigenvalues, thus all the eigenvalues of $J{|}_{{\mathcal{E}}_0}$ at ${\mathcal{E}}_0$ are negative.

Now to discuss the bifurcation analysis of the proposed model, we assume the bifurcation parameter is ${\beta }^{\ast }$, and at ${\mathcal{R}}_0=1$, the value of the bifurcation parameter, ${\beta }^{\ast }$ gives (12) ${\beta }^{\ast }=\frac{q_2{q}_3}{S^0}\{p{q}_3+\gamma {\gamma }_{1}pq+(1-p)\gamma \beta q_2\},$(12) then regarding the existence of bifurcation, we establish the following result.

Theorem 3.4

The proposed SARS-CoV-2 virus epidemiological model undergoes a backward bifurcation at ${\mathcal{R}}_0=1$.

Proof.

The Jacobian matrix of the considered model evaluated at $\beta ={\beta }^{\ast }$ and at the disease-free equilibrium which looks like $J{|}_{\beta ={\beta }^{\ast }}=\left(\begin{array}{cccc}-q_1& -\frac{{\beta }^{\ast }\mathrm{\Pi }}{q_1}& -\frac{\gamma {\beta }^{\ast }\mathrm{\Pi }}{q_1}& 0\\ 0& \frac{p{\beta }^{\ast }\mathrm{\Pi }-q_1{q}_2}{q_1}& \frac{p{\beta }^{\ast }\gamma \mathrm{\Pi }}{q_1}& 0\\ 0& \frac{(1-p){\beta }^{\ast }\mathrm{\Pi }+q_{1}q{\gamma }_1}{q_1}& \frac{(1-p)\gamma \mathrm{\Pi }{\beta }^{\ast }-q_1{q}_3}{q_1}& 0\\ \nu & {\gamma }_1(1-q)& {\gamma }_2& -{\mu }_0\end{array}\right).$Calculated the eigenvalues of the matrix, $J{|}_{\beta ={\beta }^{\ast }}$, we obtain $\begin{array}{rl}& {\lambda }_1=0,{\lambda }_2=-{\mu }_0,{\lambda }_3=-q_1,\\ & {\lambda }_4=-\frac{1}{\{\gamma q_2(1-p)+p{q}_3+{\gamma }_1\gamma pq\}}\{\gamma q_2^2(1-p)+p{q}_3^2+{\gamma }_1\gamma pq{q}_2+{\gamma }_1\gamma pq{q}_3\}.\end{array}$We observed that three of the eigenvalues are negative or having negative real parts and one zero eigenvalue. Thus, the proposed epidemiological model has hyperbolic solution which greentree the existence of bifurcation. To find the direction of bifurcation, we assume that $l=(l_1,l_2,l_3,l_4)$ and $r=(r_1,r_2,r_3,r_4)$ are the left and right eigenvector associated with the zero eigenvalue, then $\begin{array}{rl}& l_2=l_2,l_1=-\frac{q_2{l}_2}{p{q}_1},l_3=\frac{l_2}{p{q}_3}\{q_2(1-p)+pq{\gamma }_1\},\\ & l_4=\frac{l_2}{p{q}_1{q}_3{\mu }_0}\{p{q}_1{q}_3{\gamma }_1(1-q)+q_1{q}_2{\gamma }_2(1-p)+pq{q}_1{\gamma }_1{\gamma }_2-v{q}_2{q}_3)\},\end{array}$and $r_1=r_4=0,r_3=r_3,r_2=\frac{r_3}{q_2\gamma }\{q_3+q\gamma {\gamma }_1\}.$To proceed further, let us assume that ${\psi }_k$ represent the right-hand side of the k-th equation of the proposed model and $z_k$ be the state variables whose derivative is described by the k-th equation for k = 1, 2, 3, 4. Following centre manifold theory and the results stated by Theorem 4.1 in Ref. [37], the local bifurcation at around the bifurcation parameter $\beta ={\beta }^{\ast }$ is then investigated by calculating the two constants given as ${\alpha }_1=\sum _{k,i,j=1}^n{r}_k{l}_i{l}_j\frac{{\mathrm{\partial }}^2{\psi }_k(00)}{\mathrm{\partial }{y}_i\mathrm{\partial }{y}_j},{\alpha }_2=\sum _{k,i,j=1}^n{r}_k{l}_i\frac{{\mathrm{\partial }}^2{\psi }_k(0,0)}{\mathrm{\partial }{y}_i\mathrm{\partial }\varphi },$where $\varphi =\beta -{\beta }^{\ast }$ and the first zero in ${\psi }_k(0,0)$ correspond to the SARS-CoV-2 virus-free equilibrium of the model. To calculate the above constants for the proposed model, let as assume that $y_1=s-s_0$, $y_2=a$, $y_3=c$ and $y_4=r-r_0$. Thus, after some algebraic manipulation, we arrive at $\begin{array}{rl}{\alpha }_1& =\frac{2{\beta }^{\ast }{r}_3{l}_2^2}{p^2\gamma q_1{q}_2}\{(1-p)q_2\gamma +p(q_3+q\gamma {\gamma }_1){\}}^2,\\ {\alpha }_2& =\frac{4r_3{l}_2{s}_0}{p\gamma q_2{q}_3}\{(1-p)\gamma q_2+p(q_3+q\gamma {\gamma }_1){\}}^2.\end{array}$It should be noted that both the co-efficient, ${\alpha }_1>0$ and ${\alpha }_2>0$, thus the system exhibits the phenomenon of backward bifurcation.

Theorem 3.5

If ${\mathcal{R}}_0>1$, then the proposed model (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) at the disease endemic state (${\mathcal{E}}^{\ast }$) is stable.

Proof.

In a similar fashion, the linearized matrix of the system (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) at the endemic state looks like $J{|}_{{\mathcal{E}}^{\ast }}=\left(\begin{array}{cccc}-\beta a^{\ast }-\gamma \beta c^{\ast }-q_1& -\beta s^{\ast }& -\gamma \beta s^{\ast }& 0\\ p\{\beta a^{\ast }+\gamma \beta c^{\ast }\}& -q_2& 0& 0\\ (1-p)\{\beta a^{\ast }+\gamma \beta c^{\ast }\}& q{\gamma }_1& -q_3& 0\\ \nu & {\gamma }_1(1-q)& {\gamma }_2& -{\mu }_0\end{array}\right).$Clearly, $J{|}_{{\mathcal{E}}^{\ast }}$ has one negative eigenvalue i.e. ${\lambda }_1=-{\mu }_0$. We take the reduced matrix to find the other eigenvalues as follows $J{|}_{{\mathcal{E}}^{\ast }}=\left(\begin{array}{ccc}-\beta a^{\ast }-\gamma \beta c^{\ast }-q_1& -\beta s^{\ast }& -\gamma \beta s^{\ast }\\ p\{\beta a^{\ast }+\gamma \beta c^{\ast }\}& -q_2& 0\\ (1-p)\{\beta a^{\ast }+\gamma \beta c^{\ast }\}& q{\gamma }_1& -q_3\end{array}\right).$Calculating the auxiliary equation $\det (J{|}_{E^{\ast }}-\lambda I)=0$, we obtain (13) ${\lambda }^3+A{\lambda }^2+B\lambda +C=0,$(13) where $\begin{array}{rl}A& =q_1+q_2+q_3+a^{\ast }\beta +c^{\ast }\gamma \beta ,\\ B& =q_1{q}_2+q_1{q}_3+q_2{q}_3+a^{\ast }{q}_2\beta +a^{\ast }{q}_3\beta +a^{\ast }p{s}^{\ast }{\beta }^2+c^{\ast }{q}_2\beta \gamma +c^{\ast }{q}_3\beta \gamma \\ & +a^{\ast }{s}^{\ast }{\beta }^2\gamma -a^{\ast }p{s}^{\ast }{\beta }^2\gamma +c^{\ast }p{s}^{\ast }{\beta }^2\gamma +c^{\ast }{s}^{\ast }{\beta }^2{\gamma }^2(1-p),\\ C& =q_1{q}_2{q}_3+a^{\ast }{q}_2{q}_3\beta +a^{\ast }p{q}_3{s}^{\ast }{\beta }^2+c^{\ast }{q}_2{q}_3\beta \gamma +a^{\ast }{q}_2{s}^{\ast }{\beta }^2\gamma (1-p)\\ & +c^{\ast }p{q}_3{s}^{\ast }{\beta }^2\gamma +c^{\ast }{q}_2{s}^{\ast }{\beta }^2{\gamma }^2(1-p)+a^{\ast }pq{s}^{\ast }{\beta }^2\gamma {\gamma }_1+c^{\ast }pq{s}^{\ast }{\beta }^2{\gamma }^2{\gamma }_1.\end{array}$Moreover $\begin{array}{rl}AB-C& =\{q_1+q_2+q_3+a^{\ast }\beta +c^{\ast }\gamma \beta \}\{q_1{q}_2+q_1{q}_3+q_2{q}_3+a^{\ast }{q}_2\beta \\ & +a^{\ast }{q}_3\beta +a^{\ast }p{s}^{\ast }{\beta }^2+c^{\ast }{q}_2\beta \gamma +c^{\ast }{q}_3\beta \gamma +a^{\ast }{s}^{\ast }{\beta }^2\gamma (1-p)\\ & +c^{\ast }p{s}^{\ast }{\beta }^2\gamma +c^{\ast }{s}^{\ast }{\beta }^2{\gamma }^2(1-p)\}-q_1{q}_2{q}_3-a^{\ast }{q}_2{q}_3\beta -a^{\ast }p{q}_3{s}^{\ast }{\beta }^2\\ & -c{q}_2{q}_3\beta \gamma -a^{\ast }{q}_2{s}^{\ast }{\beta }^2\gamma (1-p)-c^{\ast }p{q}_3{s}^{\ast }{\beta }^2\gamma -c^{\ast }{q}_2{s}^{\ast }{\beta }^2{\gamma }^2(1-p)\\ & -a^{\ast }pq{s}^{\ast }{\beta }^2\gamma {\gamma }_1-c^{\ast }pq{s}^{\ast }{\beta }^2{\gamma }^2{\gamma }_1.\end{array}$We noted that the above coefficients A, B, and C of the auxiliary Equation (13(13) ${\lambda }^3+A{\lambda }^2+B\lambda +C=0,$(13) ) are positive, and if AB−C is positive then it ensures that the characteristic equation has all negative roots and so the proposed model around disease endemic state is stable.

Theorem 3.6

The disease-free state (${\mathcal{E}}_0$) of the epidemiological system (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) is stable globally asymptotically, if and only if ${\mathcal{R}}_0<1$.

Proof.

We investigate the global analysis of the considered problem as stated by Equation (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) around ${\mathcal{E}}_0$ by using the concept of Lyapunov function theory. On the basis of this, we define a real-valued function ${\mathcal{H}}_1:{\mathcal{R}}_{+}^4\to \mathcal{R}$ is given by (14) ${\mathcal{H}}_1(s,a,c)=\{s-s_0\}+a+c.$(14) Thus, the temporal differentiation of Equation (14(14) ${\mathcal{H}}_1(s,a,c)=\{s-s_0\}+a+c.$(14) ) and making the substitution of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) gives that (15) $\frac{\mathrm{d}{\mathcal{H}}_1(s,a,c)}{\mathrm{d}t}=\mathrm{\Pi }-q_{1}s(t)-q_{2}a(t)+q{\gamma }_{1}a(t)-q_{3}c(t).$(15) Further, the substitution of the components of the disease-free equilibrium and then a little re-arrangement leads to the following assertion (16) $\frac{\mathrm{d}{\mathcal{H}}_1(s,a,c)}{\mathrm{d}t}=-q_1(s-s_0)-\{{\mu }_0+{\mu }_1\}a(t)-{\gamma }_1(1-q)a(t)-q_{2}c(t).$(16) It can be noted from Equation (16(16) $\frac{\mathrm{d}{\mathcal{H}}_1(s,a,c)}{\mathrm{d}t}=-q_1(s-s_0)-\{{\mu }_0+{\mu }_1\}a(t)-{\gamma }_1(1-q)a(t)-q_{2}c(t).$(16) ), that $\frac{\mathrm{d}{\mathcal{H}}_1}{\mathrm{d}t}$ is negative and consequently equal to zero at (${\mathcal{E}}_0$), which satisfied the properties of Lyapunov function. Hence ${\mathcal{H}}_1$ is a Lyapunov function and therefore the invariance principle [38] then implies that the proposed epidemiological model (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) is globally stable at the infection-free state.

Theorem 3.7

The endemic disease equilibrium (${\mathcal{E}}^{\ast }$) of the system (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) is globally asymptotically stable, if and only if ${\mathcal{R}}_0>1$.

Proof.

We discuss the global analysis of the model around the endemic equilibrium. For this, we define a function given by (17) ${\mathcal{H}}_2(s,a,c)=\frac{1}{2}{\{(s-s^{\ast })+(a-a^{\ast })+(c-c^{\ast })\}}^2.$(17) Taking the temporal differentiation of ${\mathcal{H}}_2$ and making use of Equation (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) gives $\begin{array}{rl}\frac{\mathrm{d}{\mathcal{H}}_2(s,a,c)}{\mathrm{d}t}& =-\{(s-s^{\ast })+(a-a^{\ast })+(c-c^{\ast })\}\{q_1(s-s_0)\\ & +({\mu }_0+{\mu }_1)a(t)+{\gamma }_1(1-q)a(t)+q_{2}c(t)\}.\end{array}$From the above equation, we observe that $\frac{\mathrm{d}{\mathcal{H}}_2}{\mathrm{d}t}<0$, and consequently equal to zero around the endemic state ${\mathcal{E}}^{\ast }$, which justifies the axioms of Lyapunov function, so ${\mathcal{H}}_2$ is a Lyapunov function, and hence according to the LaSalle's invariance principle [38], the proposed problem (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) is globally stable at the disease endemic state (${\mathcal{E}}^{\ast }$).

4. Optimal control analysis

Optimal control theory plays a pivotal role in the epidemiology of infectious diseases by providing a rigorous framework for designing effective control strategies and guiding decision-making processes. Optimal control theory enables researchers and policymakers to identify control measures that maximize the desired public health outcomes. By formulating the epidemiological problem as an optimization task, one can seek control strategies that minimize disease burden, reduce transmission rates, limit the number of cases, or achieve other specific objectives. This approach ensures that limited resources and interventions are utilized in the most efficient and effective manner, ultimately leading to improved public health outcomes. Moreover, optimal control theory takes into account the inherent constraints and limitations of resource availability. With the finite availability of control measures in the form of vaccines, treatment facilities, healthcare personnel, and other resources, it becomes essential to allocate them optimally to achieve the greatest impact. By incorporating resource constraints into the optimization framework, optimal control theory aids in determining the optimal distribution and utilization of resources to curb the spread of infectious diseases.

In addition, the epidemiology of infectious diseases is characterized by dynamic and evolving scenarios. Optimal control theory allows for the design of adaptive control strategies that can respond to changing conditions, such as variations in disease transmissibility, the emergence of new variants, or alterations in population behaviour. By continuously updating and optimizing control measures based on real-time data and feedback, optimal control theory enables a more proactive and agile approach to disease management. Similarly, optimal control theory provides a framework to explore trade-offs and balance competing objectives in disease control. It allows for the quantification and comparison of the costs and benefits associated with different control measures, considering both short-term and long-term perspectives. This analysis assists policymakers in making informed decisions by identifying cost-effective interventions that yield the greatest health benefits within the available resources. Optimal control theory facilitates the assessment of the potential impact of different control strategies. By mathematically modelling the epidemiological dynamics and incorporating various intervention measures, researchers can simulate different scenarios and predict their outcomes. This predictive capability enables the evaluation of the effectiveness and efficiency of different interventions before their implementation, helping to inform policy decisions and prioritize interventions based on their expected impact.

We formulate an optimal control problem to mitigate the spread of SARS-CoV-2 virus transmission by assuming a set of three control measures $u_1(t)$, $u_2(t)$ and $u_3(t)$ that minimize the infected population over a given time frame. Physically, $u_1(t)$ represents the non-pharmaceutical control intervention strategy, implementing personal protection (using mask, avoiding gathering and keeping social distances etc.) that minimizes the spreading of novel Coronavirus transmission and slow down the intensity of β and γ. To optimize the interaction of infected and susceptible, we assume that $u_1(t)\in [0,1]$ and introduce the control factor, $1-u_1(t)$. This implies that whenever $u_1(t)\equiv 1$, the control is fully responsive, while having no response at $u_1(t)\equiv 0$. By full response, we mean that the interaction between susceptible and infected individuals vanishes. Moreover, the current scenario justifies that there is sufficient treatment available for those individuals which are clinically ill, therefore, we incorporate treatment as a control measure for clinically infected individuals by introducing the term $-u_2(t)c(t)$ to minimize the infected population by keeping in view the normalized sensitivity indices of ${\gamma }_1$ and ${\gamma }_2$, where the controlling attempt $u_2(t)$ represents the treatment intensity lies between zero and unity. In addition, the sensitivity index of vaccination parameter reveals that vaccination plays a vital role in the reduction of novel Coronavirus, and vaccine for the prevention of SARS-CoV-2 virus are available which provide safety up to a certain degree of realism, so to incorporate vaccination is a control measure, we incorporate $-u_3(t)s(t)$ to maximize the number of recovered population. Thus, the clear goal of our control strategies is to enhance the response of susceptible and recovered by increasing the intensities of the control function $u_1(t)$ and $u_3(t)$ to overcome the novel Coronavirus spreading, while reducing the infected population with the aid of controlling effect of $u_2(t)$. More precisely, the goal of our control mechanism is to maximize the recovered and to minimize the infected population, thus we define the following control set (18) $U=\{(u_1(t),u_2(t),u_3(t)):0\le u_i(t)\le 1,i=1,2,3,0\le t\le T\}.$(18) To describe the control problem, we formulate the objective functional to be minimized is as (19) $J={\int }_0^T\{{\xi }_{1}a(t)+{\xi }_{2}c(t)+\frac{1}{2}({\xi }_3{u}_1^2(t)+{\xi }_4{u}_2^2(t)+{\xi }_5{u}_3^2(t))\}\mathrm{d}t,$(19) subject to the control system, an extended version of the proposed model, which takes the following form (20) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\{\beta a(t)s(t)+\beta \gamma s(t)c(t)\}\{1-u_1(t)\}-\{u_3(t)+{\mu }_0\}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =p\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}\{1-u_1(t)\}-\{{\mu }_1+{\gamma }_1+{\mu }_0\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =(1-p)\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}\{1-u_1(t)\}+q{\gamma }_{1}a(t)\\ & -\{{\mu }_2+{\mu }_0+{\gamma }_2+u_2(t)\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& =\{{\gamma }_2+u_2(t)\}c(t)+u_3(t)s(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t),\end{array}$(20) with (21) $s(0)>0,a(0),c(0)\ge 0,r(0)>0.$(21) In the above, the quantities ${\xi }_i,i=1,2,3,4,5$ respectively represent the weight constants of asymptomatic, symptomatic, and the weight constants associated with the personal protection, treatment and vaccination against SARS-CoV-2 virus disease. The quadratic term ${\xi }_1{u}_1^2$ represent the cost corresponding to the control of susceptible and infected individuals interaction by personal protection, while ${\xi }_1{u}_2^2$ and ${\xi }_3{u}_3^2$ are the costs associated with the treatment and vaccination control measures. Basically, the cost correspond to the control measure $u_1(t)$ comes from the implementation of personal protection, educating the health official and public to properly wear mask, use sanitizer and keeping social distances etc., while the cost associated with the other two controls directly come from the treatment and vaccination of the individuals. The costs is assumed to be in proportion to the square of the control measure they are associated with. Thus, the objective of our control problem is to find the optimal value $(u_1^{\ast },u_2^{\ast },u_3^{\ast })$ for the control measure such that $J(u_1^{\ast },u_2^{\ast },u_3^{\ast })=min\{J(u_1,u_2,u_3):(u_1,u_2,u_3)\in U\}.$First, we prove the existence analysis for the proposed optimal control problem and then investigate the optimality system. For this, we assume that $\phi =\left(\begin{array}{c}s\\ a\\ c\\ r\end{array}\right),L=\left(\begin{array}{cccc}-q_1& 0& 0& 0\\ 0& -q_2& 0& 0\\ 0& q{\gamma }_1& -q_3& 0\\ u_3(t)& {\gamma }_1(1-q)& {\gamma }_2+u_2(t)& -{\mu }_0\end{array}\right),$and $N(\phi )=\left(\begin{array}{c}\mathrm{\Pi }-\{1-u_1(t)\}\{\beta s(t)a(t)+\beta \gamma c(t)s(t)\}\\ p\{1-u_1(t)\}\{\beta s(t)a(t)+\beta s(t)c(t)\}\\ (1-p)\{\beta s(t)a(t)+\beta s(t)c(t)\}\{1-u_1(t)\}\\ 0\end{array}\right),$then the control system (20(20) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\{\beta a(t)s(t)+\beta \gamma s(t)c(t)\}\{1-u_1(t)\}-\{u_3(t)+{\mu }_0\}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =p\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}\{1-u_1(t)\}-\{{\mu }_1+{\gamma }_1+{\mu }_0\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =(1-p)\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}\{1-u_1(t)\}+q{\gamma }_{1}a(t)\\ & -\{{\mu }_2+{\mu }_0+{\gamma }_2+u_2(t)\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& =\{{\gamma }_2+u_2(t)\}c(t)+u_3(t)s(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t),\end{array}$(20) ) can be re-written as $\frac{\mathrm{d}\phi }{\mathrm{d}t}=L\phi +N(\phi ).$To proceed further, we also define that $G(\phi )=L\phi +N(\phi )$, then for any ${\phi }_1$ and ${\phi }_2$, we have $\Vert G({\phi }_1)-G({\phi }_2)\Vert \le M\Vert {\phi }_1-{\phi }_2\Vert ,$where $M={\mu }_0+{\mu }_2$, is the Lipschitz constants implies that the function $G(\phi )$ is Lipschitz continuous, which ensure the existence analysis of the control system (20(20) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\{\beta a(t)s(t)+\beta \gamma s(t)c(t)\}\{1-u_1(t)\}-\{u_3(t)+{\mu }_0\}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =p\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}\{1-u_1(t)\}-\{{\mu }_1+{\gamma }_1+{\mu }_0\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =(1-p)\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}\{1-u_1(t)\}+q{\gamma }_{1}a(t)\\ & -\{{\mu }_2+{\mu }_0+{\gamma }_2+u_2(t)\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& =\{{\gamma }_2+u_2(t)\}c(t)+u_3(t)s(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t),\end{array}$(20) ). In addition, we give the existence analysis for optimal value of the control measure in the result given below.

Theorem 4.1

For the objective functional J and the control $(u_1,u_2,u_3)\in U$ subject to the control system, there exists optimal control value $(u_1^{\ast },u_2^{\ast },u_3^{\ast })$, such that $J(u_1^{\ast },u_2^{\ast },u_3^{\ast })=minJ(u_1,u_2,u_3)$.

Proof.

The integrand of the objective functional J is convex on the closed bounded control set U. Since, the model is linear in the controls function and is bounded by a linear system in the model state variables [39], which shows that the optimal control exist.

We now find the optimality system with the aid of Pontryagin's maximum principle by following Pontryagin et al. [40], thus we define the Hamiltonian for the proposed control system as (22) $\begin{array}{rl}\mathcal{H}& ={\xi }_{1}a(t)+{\xi }_{2}c(t)+\frac{1}{2}\{{\xi }_3{u}_1^2(t)+{\xi }_4{u}_2^2(t)+{\xi }_5{u}_3^2(t)\}\\ & +{\lambda }_1[\mathrm{\Pi }-\{1-u_1(t)\}\{\beta s(t)a(t)+\beta \gamma c(t)s(t)\}-\{{\mu }_0+u_3(t)\}s(t)]\\ & +{\lambda }_2[p\{1-u_1(t)\}\{\beta a(t)s(t)+\gamma \beta s(t)c(t)\}-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t)]\\ & +{\lambda }_3[(1-p)\{1-u_1(t)\}\{\beta s(t)a(t)+\beta \gamma c(t)s(t)\}+q{\gamma }_{1}a(t)\\ & -\{{\mu }_2+{\mu }_0+{\gamma }_2+u_2(t)\}c(t)]+{\lambda }_4[\{{\gamma }_2+u_2(t)\}c(t)+u_3(t)s(t)\\ & -{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)].\end{array}$(22) In the next theorem, we use the Hamiltonian H to present the characterization of adjoint variables and optimal control measures.

Theorem 4.2

For the given control measures $(u_1(t),u_2(t),u_3(t))$ and associated states s, a, c, r of the control system (20(20) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\{\beta a(t)s(t)+\beta \gamma s(t)c(t)\}\{1-u_1(t)\}-\{u_3(t)+{\mu }_0\}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =p\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}\{1-u_1(t)\}-\{{\mu }_1+{\gamma }_1+{\mu }_0\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =(1-p)\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}\{1-u_1(t)\}+q{\gamma }_{1}a(t)\\ & -\{{\mu }_2+{\mu }_0+{\gamma }_2+u_2(t)\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& =\{{\gamma }_2+u_2(t)\}c(t)+u_3(t)s(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t),\end{array}$(20) ), there exists ${\lambda }_i$, i = 1, 2, 3, 4, satisfying (23) $\begin{array}{rl}\frac{\mathrm{d}{\lambda }_1(t)}{\mathrm{d}t}& =\{(\beta a(t)+\beta \gamma c(t))(1-u_1(t))\}\{{\lambda }_1(t)-p{\lambda }_2(t)-(1-p){\lambda }_3(t)\}\\ & +\{{\lambda }_1(t)-{\lambda }_4(t)\}{u}_3(t)+{\mu }_0{\lambda }_1(t),\\ & \frac{\mathrm{d}{\lambda }_2(t)}{\mathrm{d}t}=-{\xi }_1+\beta s(t)\{1-u_1(t)\}\{{\lambda }_1(t)-p{\lambda }_2(t)-(1-p){\lambda }_3(t)\}\\ & +q_2{\lambda }_2(t)-q\gamma {\lambda }_3(t)+{\gamma }_1(1-q){\lambda }_4(t),\\ & \frac{\mathrm{d}{\lambda }_3(t)}{\mathrm{d}t}=-{\xi }_2+\{\beta \gamma s(t)(1-u_1(t))\}\{{\lambda }_1(t)-p{\lambda }_2(t)-(1-p){\lambda }_3(t)\}\\ & +\{q_2+u_2(t)\}{\lambda }_3(t)-\{{\gamma }_2+u_2(t)\}{\lambda }_4(t),\\ & \frac{\mathrm{d}{\lambda }_4(t)}{\mathrm{d}t}=-{\mu }_0{\lambda }_4(t).\end{array}$(23) Moreover, the terminal conditions and optimal values of the control measures are (24) ${\lambda }_i(T)=0,\mathrm{for}i=1,2,3,4$(24) and (25) $\begin{array}{rl}& u_1^{\star }(t)=max\{0,min\{1,\frac{1}{{\xi }_3}\{\beta sa+\gamma \beta sc\}\{(1-p){\lambda }_3+p{\lambda }_2-{\lambda }_1\}\}\},\\ & u_2^{\star }(t)=max\{0,min\{1,\frac{1}{{\xi }_4}\{{\lambda }_3-{\lambda }_4\}c\}\},\\ & u_3^{\star }(t)=max\{0,min\{1,\frac{1}{{\xi }_5}\{{\lambda }_1-{\lambda }_4\}s\}\}.\end{array}$(25)

Proof.

The adjoint system (23(23) $\begin{array}{rl}\frac{\mathrm{d}{\lambda }_1(t)}{\mathrm{d}t}& =\{(\beta a(t)+\beta \gamma c(t))(1-u_1(t))\}\{{\lambda }_1(t)-p{\lambda }_2(t)-(1-p){\lambda }_3(t)\}\\ & +\{{\lambda }_1(t)-{\lambda }_4(t)\}{u}_3(t)+{\mu }_0{\lambda }_1(t),\\ & \frac{\mathrm{d}{\lambda }_2(t)}{\mathrm{d}t}=-{\xi }_1+\beta s(t)\{1-u_1(t)\}\{{\lambda }_1(t)-p{\lambda }_2(t)-(1-p){\lambda }_3(t)\}\\ & +q_2{\lambda }_2(t)-q\gamma {\lambda }_3(t)+{\gamma }_1(1-q){\lambda }_4(t),\\ & \frac{\mathrm{d}{\lambda }_3(t)}{\mathrm{d}t}=-{\xi }_2+\{\beta \gamma s(t)(1-u_1(t))\}\{{\lambda }_1(t)-p{\lambda }_2(t)-(1-p){\lambda }_3(t)\}\\ & +\{q_2+u_2(t)\}{\lambda }_3(t)-\{{\gamma }_2+u_2(t)\}{\lambda }_4(t),\\ & \frac{\mathrm{d}{\lambda }_4(t)}{\mathrm{d}t}=-{\mu }_0{\lambda }_4(t).\end{array}$(23) ) directly comes from the Pontryagin's Principle [40], such that $\frac{\mathrm{d}{\lambda }_1(t)}{\mathrm{d}t}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }s},\frac{\mathrm{d}{\lambda }_2(t)}{\mathrm{d}t}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }a},\frac{\mathrm{d}{\lambda }_3(t)}{\mathrm{d}t}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }c},\frac{\mathrm{d}{\lambda }_4(t)}{\mathrm{d}t}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }r},$with zero final time conditions.

In addition to get the optimal value of the controls measures, we solve the following equations $\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{u}_1}=0,\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{u}_2}=0,\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{u}_3}=0.$To get the compact form of the control measures, we then use the bound on the controls, we get the desired characterization of the control measures.

5. Numerical simulations

The numerical simulations of the proposed SARS-CoV-2 epidemiological model (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) will be performed to present the graphical visualization of the analytical findings. We first present the graphical representation of the sensitivity for reproductive number and then the parameterizations of the model from real data of the novel disease of the Coronavirus. To visualize the results and to verify the theoretical results, we use the estimated values of the parameters retrieved from real data of SARS-CoV-2 virus reported cases.

5.1. Parameterization

For parameterization purposes, reported data of confirmed SARS-CoV-2 infected cases in Oman are used between January 2021 to May 2021. We estimate the value of parameters with the aid of Ordinary Least Square (OLS). The procedure of OLS is applied to minimize the cost representing the daily registered cases as well as the computed given in Equation (26(26) $min\left\{\sum _{i=1}^n\frac{c_i-{\hat{c}}_i}{c_i^2}\right\}.$(26) ), with the corresponding relative error in the goodness of fit (26) $min\left\{\sum _{i=1}^n\frac{c_i-{\hat{c}}_i}{c_i^2}\right\}.$(26) In Equation (26(26) $min\left\{\sum _{i=1}^n\frac{c_i-{\hat{c}}_i}{c_i^2}\right\}.$(26) ), $c_i$ is used for the cumulative reported cases, while ${\hat{c}}_i$ represents the cumulative simulated cases. In order to estimate the model parameters, we first calculate the values of Π and µ. In 2021, Oman's average life expectancy was 78.16 years, while the total population was 4.52 million. Assuming that the limiting population is $\mathrm{\Pi }/{\mu }_0=4520,000$ in the absence of disease, the birth rate is $\mathrm{\Pi }=4520,000\times {\mu }_0$, where ${\mu }_0=1/(78.16\times 365)$ per day, so we can obtain the death rate and birth rate as ${\mu }_0=0.0000350$ and $\mathrm{\Pi }=158.43$ per day, respectively. Then, we solve the proposed model using the fourth-order Runge–Kutta method in order to estimate the model's other epidemic parameters. Using the optimization algorithm, we update the parameter values to obtain a better agreement with the real data after evaluating the error. The process is repeated until the desired parameter value is reached, which meets the convergence criteria. Here, Figure 2 shows the comparison between real data and model solution, while Table 2 shows the associated values. The parameters are all measured in days.

Figure 2. The curves represent the model fitting against novel Coronavirus reported cases in Oman.

Table 2. The fitted and assumed value of the epidemic parameters using (OLS) approach.

Parameter	Value	Source	Parameter	Value	Source
ν	0.9000	Fitted	${\mu }_2$	0.0633	Estimated
β	0.0790	Fitted	γ	0.0457	Fitted
${\gamma }_2$	0.0081	Fitted	${\gamma }_1$	0.8900	Fitted
p	0.2600	Fitted	q	0.7900	Estimated
${\mu }_1$	0.8000	Estimated

5.2. Sensitivity analysis

The main purpose of this section is to provide graphical representations of the sensitivity analysis of the basic reproductive number that was conducted using Equation (9(9) $S_{\mathrm{\Phi }}^{{\mathcal{R}}_0}=\frac{\mathrm{\Phi }}{{\mathcal{R}}_0}\frac{\mathrm{\partial }{\mathcal{R}}_0}{\mathrm{\partial }\mathrm{\Phi }}.$(9) ) as demonstrated in Table 1. As shown in Figure 3(a), we plot the most influential parameters along with their sensitivity indices on x-axis. It is shown in Figure 3(b) that β and γ have positive incidences and are the most influential parameters associated with the basic reproductive number. This contour plot describes the relationship between the reproductive number and β and γ. The graphical representation of the sensitivity analysis reveals that the substantial diminishing of the parameters with positive indices while improving the parameters with non-positive indices can lead to a decrease in the basic reproductive number to control the outbreak of the SARS-CoV-2 virus epidemic.

Figure 3. The plots represent the normalized sensitivity indices and contour plot of the epidemic parameters related with the reproductive number. In this case the parametric values are: $\mathrm{\Pi }=158.43$, $\nu =0.9$, $\beta =0.079$, $\gamma =0.0457$, ${\gamma }_1=0.89$, ${\gamma }_2=0.0081$, p = 0.26, q = 0.79, ${\mu }_1=0.8$, ${\mu }_0=0.0000350$, ${\mu }_2=0.63$. (a) Normalized Sensitivity Indices, (b) Contour plot of ${\mathcal{R}}_0$ as a function of β and γ.

5.3. Validation of the model

To present the time dynamics of novel disease with the aid of the proposed epidemiological model, we use a numerical procedure of the fourth-order Runge–Kutta (RK) method [41,42] and consequently obtain the dynamics of the compartmental population of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) as given in Figure 4 which respectively demonstrate the dynamics of susceptible, asymptomatic, symptomatic and the recovered compartments. We use the estimated parametric values and obtain the value of the reproductive number is 2.98. In this case, since ${\mathcal{R}}_0>1$, the dynamics of various compartmental population of the model attain their associated endemic position which proves that the model stability at the disease endemic state as shown in Figure 4. This reveals that the infected population persists (see Figure 4(b,c)).

Figure 4. The graph visualizing the temporal dynamics of the compartmental model (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ), which shows that infected population will always exist whenever $R_0>1$. (a) Susceptible Population, (b) Asymptomatic Population, (c) Symptomatic Population and (d) Recovered Population.

On the other hand, to show the model dynamics at the infection-free state, we assume the value of model parameters as $\mathrm{\Pi }=0.25$, v = 0.01, $\beta =0.0044$, ${\gamma }_1=0.059$, $\gamma =0.045$, ${\gamma }_2=0.0081$, p = 0.36, q = 0.59, ${\mu }_1=0.09$, ${\mu }_0=0.0079$ and ${\mu }_2=0.063$. We calculate the reproductive number is ${\mathcal{R}}_0=0.1661987904$ which is less than one. Therefore, we perturb the initial sizes of the various compartments s, a, c, r of the proposed model from ${\mathcal{E}}_0$, so that the solutions curves reach disease-free equilibrium irrespective of its initial sizes in the long run as shown in Figure 5. This investigates that the model is stable around disease-free equilibrium ${\mathcal{E}}_0$. The biological interpretation reveals that in the case of disease-free equilibrium, each solution curve of s, i.e. susceptible population will tend to its equilibrium position analyses that susceptible individuals will always exist as shown in Figure 5(a). The dynamics of asymptomatic and symptomatic compartments are depicted in Figure 5(b,c), which describes that the solution curves will tend to zero, whenever ${\mathcal{R}}_0<1$ implies that the infected population vanishes. Moreover, the dynamics of recovered individuals are shown in Figure 5(d) which shows that there will be always a recovered population. Thus, we conclude that the eradication of the novel disease continuously depends on the value of the threshold quantity. In addition, the system exhibits regular and irregular periodic oscillation, if we vary the value of bifurcation parameter from 2 to 6 as shown in Figure 6, which implies that in this case there exist an unstable equilibrium.

Figure 5. The graph visualizing the temporal dynamics of the compartmental population of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\beta a(t)s(t)-\beta \gamma s(t)c(t)-\{{\mu }_0+\nu \}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}p-\{{\mu }_0+{\mu }_1+{\gamma }_1\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}(1-p)+q{\gamma }_{1}a(t)-\{{\mu }_2+{\mu }_0+{\gamma }_2\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& ={\gamma }_{2}c(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t)+\nu s(t),\end{array}$(1) ) around disease-free equilibrium, where the parameters values are as: $\mathrm{\Pi }=0.25$, $\nu =0.01$, $\beta =0.0044$, $\gamma =0.045$, ${\gamma }_1=0.059$, ${\gamma }_2=0.0081$, p = 0.36, q = 0.59, ${\mu }_1=0.09$, ${\mu }_0=0.0079$, ${\mu }_2=0.063$ and ${\mathcal{R}}_0=0.1661987904<1$. (a) Susceptible Population, (b) Asymptomatic Population, (c) Symptomatic Population and (d) Recovered Population.

Figure 6. The graphs represent the bifurcation analysis of compartmental population against the bifurcation parameter $(\beta )$ for the parameters value: $\mathrm{\Pi }=0.1$, $\gamma =0.52$, ${\mu }_0=0.02$, v = 0.0025, p = 0.2, ${\mu }_1=0.9$, ${\gamma }_1=1$, q = 0.3, ${\mu }_2=0.002$ and ${\gamma }_2=0.9$. (a) Susceptible Population, (b) Asymptomatic Population and (c) Symptomatic Population.

5.4. Optimal control implementation

To visualize the control effect graphically, we use a similar procedure of numerical analysis with the help of fourth-order Runge–Kutta method. First, we solve the state system (20(20) $\begin{array}{rl}\frac{\mathrm{d}s(t)}{\mathrm{d}t}& =\mathrm{\Pi }-\{\beta a(t)s(t)+\beta \gamma s(t)c(t)\}\{1-u_1(t)\}-\{u_3(t)+{\mu }_0\}s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}& =p\{\beta a(t)s(t)+\gamma \beta c(t)s(t)\}\{1-u_1(t)\}-\{{\mu }_1+{\gamma }_1+{\mu }_0\}a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}& =(1-p)\{\beta s(t)a(t)+\beta \gamma s(t)c(t)\}\{1-u_1(t)\}+q{\gamma }_{1}a(t)\\ & -\{{\mu }_2+{\mu }_0+{\gamma }_2+u_2(t)\}c(t),\\ \frac{\mathrm{d}r(t)}{\mathrm{d}t}& =\{{\gamma }_2+u_2(t)\}c(t)+u_3(t)s(t)-{\mu }_{0}r(t)+{\gamma }_1(1-q)a(t),\end{array}$(20) ) with initial compartmental population sizes by using the forward RK procedure, while for the solution of adjoint system, we use backward-RK method and the transversality conditions. Moreover, we used the parametric values as reported by Table 1 and while the value of weight constants are assumed as: ${\xi }_1=100$, ${\xi }_2=0.6$, ${\xi }_3=10$ and ${\xi }_4=0.9$ and ${\xi }_5=0.44$. In addition, the time units are taken to be 0 to 100. By using the Matlab software, executing the proposed algorithm of Runge–Kutta method, we obtain the graphs as demonstrated in Figure 7, which describe the temporal dynamics of the susceptible, asymptomatic, symptomatic and recovered population with the implementation (blue dashed) and without the implementation (red dashed) of the control measure. The visualization clearly illustrate our objective in implementation: to reduce the number of infected population and to maximize the non-infected individuals (recovered & susceptible) as shown in Figure 7( a–d), respectively. We noted that collectively, the implementation of the optimal control measures, $u_i(t)$, i = 1, 2, 3 slow down the intensity of the transmission co-efficient β and γ, and strengthen the recovery of infected individuals which leads to the elimination of the current pandemic of novel disease. The difference between the two cases: with the implementation and without the implementation of the control measures are clearly visible and one can eliminate the contagious infection of Coronavirus, if we can apply the proposed control strategies in a true sense. Particularly, the disease transfer coefficient and rate of vaccination are very important in the reduction of the value of ${\mathcal{R}}_0$. So using non-pharmaceutical and pharmaceutical interventions (face mask, avoiding gathering, vaccination, treatment and social distancing) are the best possible control measure to prevent the pandemic of the novel Coronavirus disease.

Figure 7. The graphs visualizing the time dynamics of the compartmental population of the proposed model with the optimal control measures implementation and with out the implementation of the control measure. (a) Susceptible Population, (b) Asymptomatic Population, (c) Symptomatic Population and (d) Recovered Population.

6. Conclusion

This study presented an epidemiological analysis and analysed the pandemic trend of SARS-CoV-2 virus transmission with optimal control strategies. We discussed the feasibility of the epidemic problem. It has been demonstrated that the considered epidemiological system is well-posed and biologically meaningful. The basic reproductive number of the model was also determined by examining the equilibria. Specifically, we discussed the positivity and boundedness of the model and calculated the infection-free state to obtain the $({\mathcal{R}}_0)$ infection-free reproduction number. Using normalized sensitivity analysis, we quantified the most sensitive parameters and their impact on disease prevention using normalized sensitivity analysis. Furthermore, we calculated the endemic equilibrium and performed a detailed stability analysis. A backward bifurcation phenomenon has been observed in the proposed model. To verify the analytical results, the proposed model was parameterized based on real data of the novel Coronavirus and simulated the problem. The disease will die out when ${\mathcal{R}}_0<1$, but persists if ${\mathcal{R}}_0>1$. We calculated the basic reproduction number as 1.406 using the estimated parameter values, which indicates that SARS-CoV-2 is endemic with the current control measures and immunization programme. Likevise, the sensitivity analysis shows that the transmission coefficient and vaccination got the highest sensitivity indices among other epidemic parameters, so vaccination and non-pharmaceutical interventions should be combined as an optimal control mechanism. As a result, using optimal control theory, we developed a control mechanism based on normalized sensitivity and stability and proposed a variety of non-pharmaceutical and pharmaceutical intervention strategies in the form of vaccinations, treatments, and personal protection. Combining these three control interventions can help eradicate novel Coronavirus contagious infections. To show the significance of the model and the work done, all theoretical results and the effects of controls are visualized.

While the model proposed yielded some interesting results, we will incorporate the stochastic nature by using stochastic differential equations in the future to capture more valuable information. Various fractional operators will also be used to discuss the dynamics of SARS-CoV-2 disease using the proposed model's associated fractional version.

Acknowledgments

The authors would like to thank the editor and reviewers for their excellent comments, which helped improve the quality of the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The work has been supported by UAE University, fund # 12S107.