Dynamics of SEPAIVRD model for COVID-19 in Bangladesh

1. Introduction

The Covid-19 pandemic, caused by the novel coronavirus SARS-CoV-2, has brought significant changes to the world in the past few years. The virus first emerged in Wuhan, China, in late 2019 and has since spread rapidly worldwide, resulting in numerous deaths and overwhelming healthcare systems. Covid-19 is primarily transmitted through respiratory droplets, and its symptoms can range from mild flu-like symptoms to severe respiratory illness, leading to hospitalization and death. Governments worldwide have implemented various measures to curb the spread of Covid-19, including social distancing, face mask mandates, and vaccination campaigns. Despite these efforts, Covid-19 remains a significant global health challenge, with new variants emerging and the virus continuing to spread. Understanding the virus's transmission dynamics and developing effective prevention and treatment strategies are critical in managing the ongoing Covid-19 pandemic.

According to Barlow & Weinstein's [1] research, they employed asymptotic approximations to develop a closed-form analytical solution of the SIR model and utilized this solution to describe the characteristics of COVID-19. By comparing the SIR model approximations with real-time COVID-19 data from Japan, they established a strong correlation between the two. Tomchin & Fradkov [2] focused on predicting the behavior of COVID-19 in Russia using discrete-time SIR and SEIR models. Initially, they utilized a discrete-time SIR model that excluded the asymptomatic class. However, they later improved their model by incorporating asymptomatic individuals in the exposed class and presented a discrete-time SEIR model. Through the application of the inverse method, Marinov & Marinova [3] investigated the dynamics of the SIR model across 12 different countries by comparing and analyzing real-time data. While acknowledging the constraints of the epidemic model, they made short-term predictions regarding the progression of the disease. Neves & Guerrero [4] focused on forecasting the potential outcomes of COVID-19 endemics in Lombardy, Italy, and Sao Paulo, Brazil. To achieve this, they employed the SAIR model, which includes an asymptomatic class.

The work of Alenezi et al. [5] involved the development of SIR and SEIR models to approximate the COVID-19 outbreak in Kuwait. They discovered that the basic reproduction number varied from 3 to 4, indicating that the failure to implement effective measures to curb the spread of COVID-19 would result in exponential growth of cases. Din & Algehyne [6] presented a novel version of the SIR model and evaluated its performance by comparing it with real COVID-19 data from Pakistan. In contrast, Malavika et al. [7] utilized the SIR model with a logistic growth curve to make short-term predictions regarding COVID-19 behavior in India. According to their study, lockdown measures may not have a significant impact on reducing the infection rate. Kozyreff [8] expanded on the SIR model by including individuals infected with COVID-19 who require hospitalization. They evaluated their findings by comparing them to real-world data from various countries, and concluded that models with curved fitting achieved higher accuracy. Conversely, Campos et al. [9] developed a discrete-time SIR model to predict the age-based vaccination strategy for COVID-19 in Brazil. Their research indicated that vaccinating different age groups could significantly reduce the number of deaths caused by the virus. Kudryashov et al. [10] utilized an analytical approach to explore the potential application of the SIR model in predicting the behavior of COVID-19. They suggested that while the SIR model may be useful for short-term predictions and small populations, modifications are required to make accurate long-term forecasts. Nonetheless, the model's simplicity has contributed to its widespread use in research for several decades. Researchers such as Cooper et al. [11], Nesteruk [12,13], Turkyilmazoglu [14], Anand et al. [15], and Lounis and Bagal [16] have also used the SIR model in various studies related to COVID-19.

Ramos et al. [17] developed the θ-SEIHQRD model to better understand the COVID-19 outbreak in Italy. In their model, they extended the classical Susceptible-Infectious-Recovered (SIR) model to include additional compartments such as Exposed (E), Hospitals (H), Quarantine (Q), and Deaths (D). By incorporating these compartments, they were able to obtain a more accurate representation of the real situation of COVID-19 in Italy. Their findings showed good agreement with the real-world data, and they also provided insights on potential future observations. Shahrear et al. [18] developed a mathematical model called the SEQIRP model to forecast the spread of COVID-19 in Bangladesh. They solved the system of first-order nonlinear ordinary differential equations with the given initial conditions using the Runge-Kutta method of order four. Their numerical investigation showed that temperature plays a significant role in the transmission of COVID-19, and that other countries have also faced similar issues. Uddin et al. [19] suggested that the winter season may not be favorable for people in Bangladesh as low temperatures could contribute to the rapid spread of COVID-19. However, Ahmed et al. [20] recently discussed that temperature and humidity have very little impact on the fast transmission of COVID-19. They argued that high population density, poor air quality, and inadequate COVID-19 testing have a more significant effect on the increase of COVID-19 cases.

The COVID-19 situation in Malaysia and Pakistan was studied by Alrabaiah et al. [21] using the fractional-order SEIQRDC model, which included a compartment for the protected population class called C. They attributed the rapid spread of COVID-19 in Pakistan to factors such as high population density, poor economic conditions, and inadequate safety measures. This finding is consistent with the research outcomes of Ahmed et al. [20]. A modified SEIR model was presented by López & Rodó [22] to forecast the spread of COVID-19 in Italy and Spain. Their model included a confined susceptible class, quarantine class, and dead population class. The modified model displayed stable behavior and was able to forecast future outcomes. Meanwhile, Hussain et al. [23] utilized the SEIQR model to conduct an analytical study and predict the future situation of COVID-19 in Pakistan. A SEIAR model was proposed by Li et al. [24] to highlight the importance of human awareness and social distancing in mitigating the spread of COVID-19. On the other hand, Youssef et al. [25] examined the spread of COVID-19 in Saudi Arabia using the SEIQR model. They conducted an analysis of reproduction numbers, equilibrium points, and stability, and utilized the Lyapunov stability theorem to demonstrate global stability. Their model was validated using actual COVID-19 data, and they found it to be in good agreement. Yu et al. [26] proposed a delay-SEIRS model to analyze the COVID-19 pandemic. They argued that social distancing, quarantine, and lockdown measures may not significantly reduce the transmission of COVID-19, but could have a substantial negative impact on the global economy. In addition to their work, Adhikary & Pal [27], Akindeinde et al. [28], Basnarkov [29], Khedher et al. [30], Pekpe et al. [31], Zhu & Shen [32], Chen et al. [33], Contreras et al. [34], and León et al. [35] have also studied modified and extended SIR models to investigate COVID-19.

To accurately portray the COVID-19 situation in Bangladesh through our research, we must construct a fresh model and ensure that it adheres to certain guidelines:

Step 1: Our research requires us to validate the compartments and ensure that each class holds significant implications for presenting accurate scenarios. For instance, the presymptomatic class holds immense significance since a large number of individuals in Bangladesh fail to comply with COVID-19 guidelines until they test positive. This is a key reason why the presymptomatic class tends to infect densely populated areas more easily. Additionally, lockdown measures are not very effective in Bangladesh since individuals with low income need to venture outside in search of work to sustain themselves. Due to a disregard for health regulations, a large number of individuals in the vulnerable class frequently venture outdoors. It is presumed that many of these individuals have obtained herd immunity against COVID-19, and consequently, do not undergo testing since they exhibit no symptoms. A recent study conducted by icddr,b, and Bangladesh Health Watch revealed that approximately 72% of low-income residents in Dhaka and Chattogram have developed antibodies. Some of these individuals experienced mild symptoms, while the majority remained asymptomatic. Hence, this particular segment of the population is categorized as the asymptomatic class. It should also be noted that the vaccinated class has emerged recently, following the Bangladesh government's decision to provide free vaccines to the public from February 7, 2021, onwards. Despite being vaccinated, individuals within this group remain susceptible to contracting COVID-19.

Step 2: The time has come to deliberate on the model's parameters. To proceed with caution, we will adopt a three-pronged approach. Firstly, we will gather data from various sources. Secondly, we will derive some parameter values from real-world data. Thirdly, we will make certain assumptions. Once we have estimated the parameter values, we can conduct simulations and compare the results with actual data. Furthermore, we will verify the precision of the present code by cross-checking it against other researchers' findings.

Step 3: Our objective is to forecast the future outcomes of the COVID-19 situation in Bangladesh and also illustrate the potential future scenarios of other South Asian countries.

2. Mathematical Modelling

Let us define each compartment of the model:

S: The susceptible class comprises individuals who have either never been infected by COVID-19 or have fully recovered from the virus.

E: The exposed class consists of individuals who are at risk of being infected by COVID-19.

P: The presymptomatic class includes individuals who have been infected with COVID-19 but have yet to display any symptoms. These individuals will eventually develop symptoms and can transmit the virus to others during this asymptomatic period. They may wait for symptoms to appear before getting tested and, in the interim, unwittingly spread COVID-19 within their local communities. Once they receive a positive COVID-19 test result, they are moved to the infected class.

A: The asymptomatic class comprises individuals who have been infected with COVID-19 but do not display any symptoms. Despite being asymptomatic, they can still transmit the virus to others. At times, the asymptomatic class is also referred to as the mild-symptomatic class, where individuals exhibit symptoms similar to seasonal flu and recover naturally without requiring any treatment. As a result, they may not opt for COVID-19 testing and remain within the asymptomatic class. This class can offer insights into the development of herd immunity against COVID-19. It's worth noting that individuals in the asymptomatic class are not considered part of the infected class.

I: The infected class consists of individuals who have been infected with COVID-19 and received a positive test result. These individuals may receive treatment either in hospitals or at home.

V: The vaccinated class refers to individuals who have received one or two doses of COVID-19 vaccines and have either fully recovered from COVID-19 or have never been infected before. While they have some protection against COVID-19, there is still a possibility of getting infected between doses or even after completing the vaccination regimen, depending on the COVID-19 variants. Therefore, individuals in this class should strictly adhere to all health protocols and precautions as they still have a chance of contracting the virus.

R: The Recovered class comprises people who have fully recovered from COVID-19 and moved to the susceptible class. These individuals have successfully gone through the infected and quarantine stages, overcome the severity of COVID-19, and received post-treatment care, resulting in their recovery from the virus and the development of immunity. After being exposed to COVID-19, this class of people generates antibodies within a few days or weeks.

D: The Death class refers to the group of people who have passed away solely due to COVID-19 infection.

In this paper, we have considered a deterministic system of non-linear differential equations to formulate the model. The total human population, N_h (t), is divided into the following eight compartments: susceptible individuals (S), vaccinated individuals (V), individuals exposed with COVID-19 (E), pre-symptomatic individuals (P), asymptomatic individuals (A), COVID-19 infected individuals (I), recovered individuals (R).

$$\lambda =\beta \left(I+\tau P\right)$$

(1)

Assuming a constant recruitment rate (Λ) during the model formulation, the rate of susceptible individuals getting infected with COVID-19 is denoted by $\lambda$. Another assumption is that individuals who have recovered from COVID-19 may become susceptible again. Based on these assumptions, we can derive the following set of differential equations.

$${\displaystyle \frac{dS}{dt}}=\Lambda -\beta \left(I+\tau P\right)S+\sigma R-K_{1}S$$

(2a)

$${\displaystyle \frac{dV}{dt}}=\rho S-\left(1-\epsilon \right)\beta \left(I+\tau P\right)V-K_{2}V$$

(2b)

$${\displaystyle \frac{dE}{dt}}=\beta \left(I+\tau P\right)S+\left(1-\epsilon \right)\beta \left(I+\tau P\right)V-K_{3}E$$

(2c)

$${\displaystyle \frac{dP}{dt}}=\omega E-K_{4}P$$

(2d)

$${\displaystyle \frac{dA}{dt}}=\omega E-K_{5}A$$

(2e)

$${\displaystyle \frac{dI}{dt}}=\omega E+{\kappa }_{1}P-K_{6}I$$

(2f)

$${\displaystyle \frac{dR}{dt}}={\kappa }_{2}I-K_{7}R$$

(2g)

$${\displaystyle \frac{dD}{dt}}=\gamma I$$

(2e)

where

$$K_1=\rho +\mu ,K_2=\mu ,K_3=3\omega +\mu ,K_4={\kappa }_{1}I+\mu ,K_5=\mu ,K_6=\gamma I+{\kappa }_2+\mu ,K_7=\sigma +\mu .$$

Here $\beta$ is the rate of contact for COVID-19 infection. The modification parameters $\tau <1$ indicate low infectiousness of the individuals in the P class compared to the individuals in the I class. ρ is the vaccination rate provided to susceptible individuals. Susceptible populations move to the exposed class at a rate $\lambda$ upon contact with infected individuals. Due to immunity, among the susceptible individuals those who are vaccinated moved to the exposed class at a lower rate $(1-\epsilon )\lambda$. Individuals exposed to COVID-19 progress to the P, A and I classes at a rate $\omega$. Individuals in the pre-symptomatic class (P) move to the infected class (I) at a rate ${\kappa }_1$. ${\kappa }_2$ represents the recovery rate of infected individuals. $\gamma$ is the disease-induced death rate for infected individuals in I class. $\mu$ corresponds to the natural death rate. Recovered individuals move to the susceptible class at a rate $\sigma$.

Initial Conditions:

$$\begin{gathered} S\left(0\right)=S_0,E\left(0\right)=E_{0}P\left(0\right)=P_0,A\left(0\right)=A_0,I\left(0\right)=I_0, \\ V\left(0\right)=V_0,R\left(0\right)=R_0,D\left(0\right)=D_0 \end{gathered}$$

3. Mathematical analysis of the COVID-19 model

3.1. Positivity and Boundedness of the Solution of the Model

We have the following theorem.

Theorem 1: If $S\left(0\right)\ge 0,V\left(0\right)\ge 0,E\left(0\right)\ge 0,P\left(0\right)\ge 0,A\left(0\right)\ge 0,I\left(0\right)\ge 0\mathrm{a}\mathrm{n}\mathrm{d}R\left(0\right)\ge 0$ and all the parameters are positive then solutions $\left\{S,V,E_c,I_c,H\mathrm{a}\mathrm{n}\mathrm{d}R\right\}$ of the system (2) will remain positive and bounded for all time $t>0$.

Proof. From the first equation of model (2) we

$${\displaystyle \frac{d}{dt}}\left[S\left(t\right)\exp \left\{{\int }_0^t\lambda \left(r\right)dr+K_{1}t\right\}\right]\ge \mathsf{\Lambda }\left[\exp \left\{{\int }_0^t\lambda \left(r\right)dr+K_{1}t\right\}\right]$$

(3)

Solving the differential equation (3), the following expression can be obtained

$$S\left(t\right)\exp \left\{{\int }_0^t\lambda \left(r\right)dr+K_{1}t\right\}-S\left(0\right)\ge {\int }_0^t\mathsf{\Lambda }\left[\exp \left\{{\int }_0^x\lambda \left(r\right)dr+K_{1}t\right\}\right]dx$$

which can also be written as,

$$S\left(t\right)\ge S\left(0\right)\exp \left\{-{\int }_0^t\lambda \left(r\right)dr+K_{1}t\right\}+\exp \left\{-{\int }_0^t\lambda \left(r\right)dr+K_{1}t\right\}{\int }_0^t\mathsf{\Lambda }\left[\exp \left\{{\int }_0^x\lambda \left(r\right)dr+K_{1}t\right\}\right]dx>0$$

Similar expressions can be found for other variables and hence, we have

$V\ge 0,E\ge 0,P\ge 0,A\ge 0,I\ge 0,\mathrm{a}\mathrm{n}\mathrm{d}R\ge 0$ for all $t\ge 0$.

Now adding all the equations of model (2), we can obtain

$${\displaystyle \frac{dN}{dt}}=\mathsf{\Lambda }-\mu N-\gamma I$$

(4)

Again, we can write

$$\mathsf{\Lambda }-\left(\mu +\gamma \right)N_h\le {\displaystyle \frac{dN}{dt}}\le \mathsf{\Lambda }-\mu N$$

(5)

Hence,

$${\displaystyle \frac{\mathsf{\Lambda }}{\mu +{\delta }_1}}\le \underset{t\to \infty }{\lim } \inf N\le \underset{t\to \infty }{\lim } \sup N\le {\displaystyle \frac{\mathsf{\Lambda }}{\mu }}$$

which implies

$$\underset{t\to \infty }{\lim } \sup N\le {\displaystyle \frac{\mathsf{\Lambda }}{\mu }}$$

Again equations (4) and (5) implies

$${\displaystyle \frac{dN}{dt}}\le \mathsf{\Lambda }-\mu N$$

$$\Rightarrow N\left(t\right)\le N\left(0\right)e^{-\mu t}+{\displaystyle \frac{\mathsf{\Lambda }}{\mu }}\left(1-e^{-\mu t}\right)$$

(6)

When

$$t=0,N\left(t\right)\le {\displaystyle \frac{\mathsf{\Lambda }}{\mu }}ifN\left(0\right)\le {\displaystyle \frac{\mathsf{\Lambda }}{\mu }}$$

implying the region $\Omega =\left\{(S,V,E,P,A,I,R)ϵR_{+}^7\right\}:N\le {\displaystyle \frac{\mathsf{\Lambda }}{\mu }}$ is positively invariant [39]. Hence all the solutions of model (2) remains in $\Omega$ for all time $t>0$.

1.2. Local Asymptotic Stability of the DFE

From the COVID-19 model (2), the disease-free equilibrium (DFE), $\mathcal{D}$, is obtained as

$$\mathcal{D}=\left(S^{*},V^{*},E^{*},P^{*},A^{*},I^{*},R^{*}\right)=({\displaystyle \frac{\mathsf{\Lambda }}{\rho +\mu }},{\displaystyle \frac{\mathsf{\Lambda }\rho }{\mu \left(\rho +\mu \right)}},0,0,0,0,0)$$

(7)

The local asymptotic stability of the DFE can be established using the next-generation operator method [36,37] on the system (2). For system (2), we have

$$F=\left(\begin{array}{ccc}0& \tau \mathcal{F}& \begin{array}{cc}0& \mathcal{F}\end{array}\\ 0& 0& \begin{array}{cc}0& 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}\begin{array}{cc}0& 0\end{array}\\ \begin{array}{cc}0& 0\end{array}\end{array}\end{array}\right),V=\left(\begin{array}{ccc}{K}_3& 0& \begin{array}{cc}0& 0\end{array}\\ -\omega & K_4& \begin{array}{cc}0& 0\end{array}\\ \begin{array}{c}-\omega \\ -\omega \end{array}& \begin{array}{c}0\\ {\kappa }_1\end{array}& \begin{array}{c}\begin{array}{cc}{K}_5& 0\end{array}\\ \begin{array}{cc}0& K_6\end{array}\end{array}\end{array}\right)$$

where, F associates with new infection terms and V represents the remaining transfer terms and

$$\mathcal{F}=\beta \left[S+\left(1-\epsilon \right)V\right].$$

Following the approach described in Chavez et al. [38] and Hethcote [39], it can be shown that the control reproduction number, denoted by ${\mathcal{R}}_c$, is given by

$${\mathcal{R}}_c={\displaystyle \frac{\tau \mathcal{F}\omega }{K_3{K}_4}}+{\displaystyle \frac{\mathcal{F}\omega K_4}{{\mu K}_3{K}_5}}$$

(8)

Using Theorem 2 of Van den Driessche and Watmough [37] the following result can be established.

Lemma 1. The DFE of the COVID-19 model (2), is locally-asymptotically stable (LAS) if ${\mathcal{R}}_c<1$ and unstable if ${\mathcal{R}}_c>1$.

3.3. Endemic Equilibrium Point (EEP)

Let ${\mathcal{E}}_c=\left(S^{\mathrm{*}},V^{\mathrm{*}},E^{\mathrm{*}},P^{\mathrm{*}},A^{\mathrm{*}},I^{\mathrm{*}},R^{\mathrm{*}}\right)$ be any arbitrary equilibrium of the model (2) and let

$${\lambda }^{*}=\beta \left(I^{*}+\tau P^{*}\right)$$

(9)

be the force of infection at steady-state. Now, setting the left-hand side of each equation of the model (2), the following expressions can be obtained:

$$\begin{gathered} S^{*}={\displaystyle \frac{\mathsf{\Lambda }+\sigma R^{*}}{{\lambda }^{*}+K_1}},V^{*}={\displaystyle \frac{\rho S^{*}}{\left(1-\epsilon \right){\lambda }^{*}+K_2}},E^{*}={\displaystyle \frac{{\lambda }^{*}\left(a_1{\lambda }^{*}+a_2\right)}{a_3{{\lambda }^{*}}^2+a_4{\lambda }^{*}+a_5}}, \\ {P}^{*}=a_p{E}^{*},A^{*}=a_a{E}^{*},I^{*}=a_i{E}^{*},R^{*}=a_r{E}^{*} \end{gathered}$$

(10)

where

$$\begin{gathered} a_1=\mathsf{\Lambda }\left(1-\epsilon \right),a_2=\mathsf{\Lambda }\left\{\left(1-\epsilon \right)\rho +K_2\right\},a_3=\left(1-\epsilon \right)K_3,a_4=K_3\left[K_2+\left(1-\epsilon \right)K_1\right] \\ {a}_5=K_1{K}_2{K}_3,a_P={\displaystyle \frac{2\omega }{K_4-K_1}},a_a={\displaystyle \frac{\omega }{K_5}},a_i={\displaystyle \frac{\omega +{\kappa }_1{a}_P}{K_6}},a_r={\displaystyle \frac{{\kappa }_2}{K_7}}{a}_i \end{gathered}$$

Substituting (10) into (9) gives

$${\lambda }^{*}=\beta \left(a_i+\tau a_p\right){\displaystyle \frac{{\lambda }^{*}\left(a_1{\lambda }^{*}+a_2\right)}{a_3{{\lambda }^{*}}^2+a_4{\lambda }^{*}+a_5}}$$

(11)

After some algebraic calculation the following polynomial equation can be obtained in terms of ${\lambda }^{*}$:

$${\lambda }^{*}\left\{P_2{{\lambda }^{*}}^2+P_1{\lambda }^{*}+P_0\right\}=0$$

(12)

where

$$P_2=a_3$$

$$P_1=a_4-\beta \left(a_i+\tau a_p\right)a_1$$

$$P_0=a_5-\beta \left(a_i+\tau a_p\right)a_2$$

Out of the three roots, the root ${\lambda }^{*}=0$, of Eq. (12), corresponds to the DFE ($\mathcal{D})$. Equation (12) says that the non-zero equilibria of the model satisfy

$$f\left({\lambda }^{*}\right)=P_2{{\lambda }^{*}}^2+P_1{\lambda }^{*}+P_0=0$$

(13)

We have the following result.

Theorem 2. The COVID-19 only model (2) has

(a) A unique endemic equilibrium if $P_0<0$ (i.e. ${\mathcal{R}}_c>1$),

(b) A unique endemic equilibrium if $P_1<0 and P_0=0 or P_1^2-4P_2{P}_0=0$,

(c) Two endemic equilibria if $P_1<0 and P_0>0\left(i.e.{\mathcal{R}}_c<1\right) and P_1^2-4P_2{P}_0>0$

(d) No endemic equilibrium otherwise.

3.4. Backward Bifurcation Analysis of the COVID-19 Model

From condition-(c) of Theorem 2, we have the possibility of having two endemic equilibrium whenever ${\mathcal{R}}_c<1$ and hence there may exist backward bifurcation phenomenon. Now we will apply center manifold theory [37,40] to examine the existence of backward bifurcation. For this let us consider the changes of variable as $S=x_1,V=x_2,E=x_3,P=x_4,A=x_5,I=x_6,andR=x_7,$ so that in vector form the model (2) can be written as

$${\displaystyle \frac{dX}{dt}}={\left(f_1,f_2,f_3,f_4,f_5,f_6,f_7\right)}^T$$

where $X={(x_1,x_2,x_3,x_4,x_5,x_6,x_7)}^T$ and then we have

$${\displaystyle \frac{d{x}_1}{dt}}=f_1=\Lambda -\beta \left(x_6+\tau x_4\right)x_1+\sigma x_7-K_1{x}_1$$

(14a)

$${\displaystyle \frac{d{x}_2}{dt}}=f_2=\rho x_1-\left(1-\epsilon \right)\beta \left(x_6+\tau x_4\right)x_2-K_2{x}_2$$

(14b)

$${\displaystyle \frac{d{x}_3}{dt}}=f_3=\beta \left(x_6+\tau x_4\right)x_1+\left(1-\epsilon \right)\beta \left(x_6+\tau x_4\right)x_2-K_3{x}_3$$

(14c)

$${\displaystyle \frac{d{x}_4}{dt}}=f_4=\omega x_3-K_4{x}_4$$

(14d)

$${\displaystyle \frac{d{x}_5}{dt}}=f_5=\omega x_3-K_5{x}_5$$

(14e)

$${\displaystyle \frac{d{x}_6}{dt}}=f_6=\omega x_3+{\kappa }_1{x}_4-K_6{x}_6$$

(14f)

$${\displaystyle \frac{d{x}_7}{dt}}=f_7={\kappa }_2{x}_6-K_7{x}_7$$

(14g)

The Jacobian of the system (14) is given by:

where

$$J_1={\displaystyle \frac{\beta \mathsf{\Lambda }}{\rho +\mu }},J_2={\displaystyle \frac{\left(1-\epsilon \right)\rho \beta \mathsf{\Lambda }}{\mu \left(\rho +\mu \right)}},J_3={\displaystyle \frac{\beta \left\{\left(1-\epsilon \right)\rho +\mu \right\}\mathsf{\Lambda }}{\mu \left(\rho +\mu \right)}}$$

Now consider, ${\mathcal{R}}_c=1and\beta ={\beta }^{*}$ is a bifurcation parameter. Thus we get,

$$\beta ={\beta }^{*}={\displaystyle \frac{{\mu K}_3{K}_4{K}_5}{\tau \omega \mu K_5+\omega K_4^2\left\{S^{*}+\left(1-\epsilon \right)V^{*}\right\}}}$$

The Jacobian $\mathcal{J}\left({\mathcal{E}}_c\right)$ of (14) with $\beta ={\beta }^{*}$, denoted by ${\mathcal{J}}_{{\beta }^{*}}$, has a simple zero eigenvalue (with all other eigenvalues having negative real part). Hence, the center manifold theory can be used to analyze the dynamics of the model (2).

Eigenvector of ${\mathcal{J}}_{{\beta }^{\mathrm{*}}}={\left.\mathcal{J}\left({\mathcal{E}}_c\right)\right|}_{\beta ={\beta }^{\mathrm{*}}}$:

When ${\mathcal{R}}_c=1$, the Jacobian (${\mathcal{J}}_{{\beta }^{*}}$) of (14) has a right eigenvector given by $W={(w_1,w_2,w_3,w_4,w_5,w_6,w_7)}^T$ where,

$$w_1={\displaystyle \frac{J_1{w}_4+\tau J_1{w}_6-\sigma w_7}{K_1}},w_2={\displaystyle \frac{\rho w_1-J_2{w}_4+\tau {J_{2}w}_6}{K_2}},w_3=w_3,$$

$$w_4={\displaystyle \frac{\omega w_3}{K_4}},w_5={\displaystyle \frac{\omega w_3}{K_5}},w_6={\displaystyle \frac{\omega w_3+{\kappa }_1{w}_4}{K_6}},w_7={\displaystyle \frac{{\kappa }_2{w}_6}{K_7}}$$

Further, ${\mathcal{J}}_{{\beta }^{\mathrm{*}}}$ has a left eigenvector $v=[v_1,v_2,v_3,v_4,v_5,v_6{,v}_7]$, where,

$$v_1=v_2=0,v_3={\displaystyle \frac{\omega \left(v_4+v_6\right)}{K_3}},v_4=v_4,v_5=0,v_6=v_6,v_7=0$$

Computation of $a$ and $b$ [40,41]:

After some tedious manipulations it can be shown that

$$a=2v_3\beta \left(\tau w_6+w_4\right)\left\{w_1+w_2\left(1-\epsilon \right)\right\}$$

$$b={\displaystyle \frac{\mathsf{\Lambda }{v}_3\left(\tau w_6+w_4\right)\left\{\mu +\rho \left(1-\epsilon \right)\right\}}{\mu \left(\mu +\rho \right)}}$$

Thus, it follows from Theorem 4.1 of Castillo-Chavez & Song [41] that the COVID-19 model (2) undergoes backward bifurcation at ${\mathcal{R}}_c=1$ whenever $a>0\mathrm{a}\mathrm{n}\mathrm{d}b>0.$

3.4.1. Local asymptotic stability of Endemic Equilibrium Point (EEP) for the COVID-19 model

Theorem 3. The EEP of the COVID-19 model (2) given by ${\mathcal{E}}_c$, is locally asymptotically stable (LAS) whenever ${\mathcal{R}}_c>1$.

3.5. Global Stability of DFE When There Is No Re-Infection

Theorem 4. The DFE of the COVID-19 model (2) with no re-infection, given by $\mathcal{D}$, is globally asymptotically stable (GAS) whenever ${\mathcal{R}}_c\le 1$.

Proof. We consider the following linear Lyapunov function:

$$\mathcal{L}={\mathcal{l}}_{v_1}E+{\mathcal{l}}_{v_2}E+{\mathcal{l}}_{v_3}I,$$

where,

$${\mathcal{l}}_{v_1}={\displaystyle \frac{\tau \omega }{K_3{K}_4}}+{\displaystyle \frac{\omega }{K_3{K}_5}}+{\displaystyle \frac{\omega {\kappa }_1}{K_3{K_{4}K}_5}},{\mathcal{l}}_{v_2}={\displaystyle \frac{\tau }{K_4}}+{\displaystyle \frac{{\kappa }_1}{{K_{4}K}_5}},{\mathcal{l}}_{v_3}={\displaystyle \frac{1}{K_5}}$$

Differentiating the above Lyapunov function and substituting the expression of $\dot{E},\dot{P},\mathrm{a}\mathrm{n}\mathrm{d}\dot{I},$ we have the following

$$\dot{\mathcal{L}}\le \tau P\left({\mathcal{R}}_c-1\right)$$

Also $\dot{\mathcal{L}}=0$ if and only if ${\mathcal{R}}_c=1$. Hence $\dot{\mathcal{L}}\le 0$. Therefore, $\mathcal{L}$ is a Lyapunov function on $\mathcal{D}$ and thus it follows by the LaSalle’s invariance principle [42] that, the DFE of the model (2) is globally asymptotic stable whenever ${\mathcal{R}}_c\le 1$.

3.6. ${\mathcal{R}}_c$ Sensitivity Analysis

In order to manage the disease, it is necessary to ensure that ${\mathcal{R}}_c$ is less than 1. This can be achieved by regulating the model parameters that influence the value of ${\mathcal{R}}_c$. If the differentiation of ${\mathcal{R}}_c$ with respect to a specific parameter yields a positive sign, it means that increasing (decreasing) the parameter value will cause ${\mathcal{R}}_c$ to increase (decrease). Conversely, if the differentiation of ${\mathcal{R}}_c$ with respect to a parameter yields a negative sign, it implies that decreasing (increasing) the parameter value will lead to an increase (decrease) in ${\mathcal{R}}_c$. It can be easily verified that

$$\begin{array}{l}{\displaystyle \frac{\partial {\mathcal{R}}_c}{\partial \tau }}={\displaystyle \frac{\left\{\beta \mathsf{\Lambda }\left(\mu +\rho -\epsilon \rho \right)\omega \right\}}{\mu \left(\mu +\rho \right)\left(\mu +3\omega \right)\left(\mu +{\kappa }_1\right)}}>0\\ {\displaystyle \frac{\partial {\mathcal{R}}_c}{\partial \omega }}={\displaystyle \frac{\beta \mathsf{\Lambda }\left(\mu +\rho -\epsilon \rho \right)\left\{\mu \tau \gamma +2\mu {\kappa }_1+{\kappa }_1^2+\mu \left(\mu +\mu \tau +\tau {\kappa }_2\right)\right\}}{\mu \left(\mu +\rho \right){\left(\mu +3\omega \right)}^2\left(\mu +{\kappa }_1\right)\left(\mu +\gamma +{\kappa }_2\right)}}>0\\ {\displaystyle \frac{\partial {\mathcal{R}}_c}{\partial \rho }}=-{\displaystyle \frac{\beta \mathsf{\Lambda }\mathsf{\omega }\mathsf{\epsilon }\left\{\mu \tau \gamma +2\mu {\kappa }_1+{\kappa }_1^2+\mu \left(\mu +\mu \tau +\tau {\kappa }_2\right)\right\}}{\mu {\left(\mu +\rho \right)}^2{\left(\mu +3\omega \right)}^2\left(\mu +{\kappa }_1\right)\left(\mu +\gamma +{\kappa }_2\right)}}<0\\ {\displaystyle \frac{\partial {\mathcal{R}}_c}{\partial \gamma }}=-{\displaystyle \frac{\beta \mathsf{\Lambda }\left(\mu +\rho -\epsilon \rho \right)\left(\mu +{\kappa }_1\right)}{{\mu }^2\left(\mu +\rho \right)\left(\mu +3\omega \right){\left(\mu +\gamma +{\kappa }_2\right)}^2}}<0\\ {\displaystyle \frac{\partial {\mathcal{R}}_c}{\partial {\kappa }_2}}=-{\displaystyle \frac{\beta \mathsf{\Lambda }\left(\mu +\rho -\epsilon \rho \right)\left(\mu +{\kappa }_1\right)}{{\mu }^2\left(\mu +\rho \right)\left(\mu +3\omega \right){\left(\mu +\gamma +{\kappa }_2\right)}^2}}<0\\ {\displaystyle \frac{\partial {\mathcal{R}}_c}{\partial \beta }}={\displaystyle \frac{\omega \mathsf{\Lambda }\left(\mu +\rho -\epsilon \rho \right)\left\{\mu \tau \gamma +2\mu {\kappa }_1+{\kappa }_1^2+\mu \left(\mu +\mu \tau +\tau {\kappa }_2\right)\right\}}{{\mu }^2\left(\mu +\rho \right)\left(\mu +{\kappa }_1\right)\left(\mu +3\omega \right){\left(\mu +\gamma +{\kappa }_2\right)}^2}}>0\end{array}$$

(15)

Table 1. Estimated parameters for model (2).

Parameter	Baseline values	Units	References
$\mathsf{\Lambda }$	10,000	Day-1	Assumed
$\beta$	0.0005	Day-1	[43]
$\omega$	0.05	-	[44]
$\rho$	0.80	Day-1	[45]
$\epsilon$	0.85	-	Assumed
$\tau$	0.8	-	Assumed
$\sigma$	0.0001	Day-1	[43]
$\gamma$	0.0095	Day-1	[44]
${\kappa }_1,{\kappa }_2$	0.15, 0.2	Day-1	[43]
$\mu$	0.0001	Day-1	[44]

Table 1. Estimated parameters for model (2).

Parameter	Baseline values	Units	References
$\mathsf{\Lambda }$	10,000	Day-1	Assumed
$\beta$	0.0005	Day-1	[43]
$\omega$	0.05	-	[44]
$\rho$	0.80	Day-1	[45]
$\epsilon$	0.85	-	Assumed
$\tau$	0.8	-	Assumed
$\sigma$	0.0001	Day-1	[43]
$\gamma$	0.0095	Day-1	[44]
${\kappa }_1,{\kappa }_2$	0.15, 0.2	Day-1	[43]
$\mu$	0.0001	Day-1	[44]

4. Numerical Simulations

In this section, numerical simulations of the model (2) are performed in this section using the parameter values as given in table 1.

4.1. Numerical Method and Validations

The well-known Euler method was utilized in this research to solve system (2) along with the initial conditions. The numerical method's specifics are described in Burden et al. [46]. To validate the numerical code, multiple validation steps were undertaken. Initially, three COVID-19 models, namely SIR, SEIR, and SIQR, were examined, which were presented by Turkyilmazoglu [14], Weinstein et al. [47], and Odagaki [48], respectively. The outcomes are illustrated in Figure 1(a, b, c), and they demonstrate that our results are in excellent agreement with those of Turkyilmazoglu [14], Weinstein et al. [47], and Odagaki [48]. Our findings suggest that we can continue the research with our proposed model, as there is good agreement.

4.2. Effect of Vaccines

4.2.1. Effect of COVID-19 Vaccine Coverage

The results are presented in Figure 2(a, b), which demonstrates the effects of vaccination coverage on reducing the number of infected individuals and deaths. According to the results, there is a significant decrease in the number of infected individuals and deaths as the number of vaccinated individual’s increases. Figure 2(a) specifically shows that the peak of daily infected cases decreases by 11% when vaccination coverage increases from 10% to 40%. This means that a higher vaccination coverage could lead to fewer individuals being infected daily, thereby reducing the overall burden on the healthcare system and decreasing the risk of transmission to others. Overall, results suggest that increasing vaccination coverage can significantly reduce the number of COVID-19 infections and deaths, which is critical in controlling the pandemic's spread.

4.2.2. Effect of COVID-19 Vaccine Efficacy

Results are presented in Figure 3(a, b), which shows the effects of vaccine efficacy on reducing the number of daily infected cases and deaths. According to the results, there is a considerable reduction in the number of daily infected cases and deaths as vaccine efficacy increases. For example, Figure 3(a) shows that when vaccine efficacy increases from 60% to 80%, the peak of daily infected cases drops from 16750 to 15210, representing a 10% decrease in the number of daily infected cases. This indicates that vaccines with higher efficacy can significantly reduce the spread of the virus, leading to fewer infections and less severe disease outcomes. These results are crucial in determining the effectiveness of various vaccines and their potential impact on the control of the COVID-19 pandemic. The simulation results suggest that vaccines with higher efficacy could play a vital role in reducing the number of COVID-19 cases and deaths, especially in populations with lower vaccination coverage. In summary, results presented in Figure 3(a, b) indicate that the impact of vaccine efficacy is significant in reducing the number of daily infected cases and deaths, emphasizing the importance of using vaccines with higher efficacy in controlling the spread of the COVID-19 virus.

4.3. Effect of τ

The simulation results for different values of τ are displayed in Figure 4(a, b). These figures demonstrate that both the daily infected cases and cumulative deaths rise with an increase in τ. For example, in Figure 4(a), when τ equals 0.1, the peak of daily infected cases is 13,400. However, when τ is 0.4, the peak of daily infected cases increases to 14,430, resulting in an 8% rise in the peak of daily infected cases. Similar patterns can be observed in Figure 4(b).

4.4. Effect of ω

Numerical simulation is conducted to evaluate the impact of ω on disease transmission. The simulation results for model (2) with varying values of ω are presented in Figure 5(a, b). These figures illustrate that the complexity of COVID-19 increases as ω rises. For example, Figure 5(a) indicates a 27% increase in the peak of daily infected cases when ω increases from 0.01 to 0.05.

4.5. Effect of ${\kappa }_1$

Furthermore, numerical simulation is carried out to investigate the impact of ${\kappa }_1$. The simulation results are displayed in Figure 6(a, b), indicating that an increase in ${\kappa }_1$ leads to a rise in both daily infected cases and cumulative deaths. For instance, Figure 6(a) reveals a 16% decrease in the peak value of daily infected cases when ${\kappa }_1$ increases from 0.2 to 0.4. Moreover, the peak value of daily infected cases further increases by 5% when ${\kappa }_1$ changes from 0.4 to 0.7. Figure 6(b) demonstrates a similar trend.

4.6. Effect of ${\kappa }_2$

To evaluate the impact of ${\kappa }_2$, numerical simulation is performed using model (2) with varying values of ${\kappa }_2$. The simulation results are presented in Figure 7, indicating that an increase in ${\kappa }_2$ leads to a decrease in both daily infected cases and cumulative deaths. For example, in Figure 7(a), when ${\kappa }_2$ = 0.2, the peak of daily infected cases is 17540. However, when ${\kappa }_2$ = 0.4, the peak of daily infected cases reduces to 10570, resulting in a 40% decrease in the peak value of daily infected cases. The same trend is observed in Figure 7(b).

4.7. Effect of Transmission Rate

The figures 8(a, b) indicate that strict adherence to social distancing measures can play an important role in controlling the spread of COVID-19 infections and associated deaths. Figure 8(a) shows that when the contact rate is reduced from 0.0005 to 0.0002, the peak of daily infected cases decreases from 17,600 to 17,300, representing a 2% decrease in the peak value of daily infected cases. A similar pattern is observed in Figure 8(b). These results suggest that reducing the effective contact rate by adhering to social distancing measures can be an effective strategy for controlling the spread of COVID-19. The numerical simulations provide important insights into the impact of different variables on disease transmission and can inform public health policies and interventions aimed at reducing the spread of the virus.

5. Conclusion

Controlling the transmission of COVID-19, a highly infectious and deadly disease, is crucial. Mathematical models can provide insights into disease dynamics and potential methods for control. Thus, a new mathematical model has been developed to address these concerns. The first step in our analysis was to examine the model theoretically. We established the positivity and boundedness of the model's solutions, and calculated its basic reproduction number. Additionally, we used a theorem to demonstrate the local asymptotic stability of the DFE when the basic reproduction number is less than one. The study further demonstrated the local asymptotic stability of the EEP when the reproduction number is greater than one. The existence of the backward bifurcation phenomenon was examined using center manifold theory, revealing that it may occur under specific conditions. Additionally, the model's globally asymptotic stability was demonstrated for the DFE under certain conditions when the basic reproduction number is less than one.

In addition to the theoretical analysis, we conducted numerical simulations to further validate our findings and to explore the impact of various parameters on the dynamics of the model. The simulations confirmed that widespread and effective vaccination can play a significant role in controlling the spread of COVID-19. The results also demonstrated that reducing the contact rate can be an effective strategy for decreasing the number of infections and deaths related to COVID-19. Specifically, our simulations showed that when highly effective vaccines were administered to a large portion of the population, the number of daily infected cases and deaths decreased significantly. Conversely, when the effective contact rate was reduced, the number of daily infected cases also decreased. These findings highlight the importance of implementing measures to reduce contact rates and promote vaccination as part of an effective strategy for controlling the spread of COVID-19.