Advancing COVID-19 Understanding: Simulating Omicron Variant Spread Using Fractional-Order Models and Haar Wavelet Collocation

Abstract

This study presents a novel approach for simulating the spread of the Omicron variant of the SARS-CoV-2 virus using fractional-order COVID-19 models and the Haar wavelet collocation method. The proposed model considers various factors that affect virus transmission, while the Haar wavelet collocation method provides an efficient and accurate solution for the fractional derivatives used in the model. This study analyzes the impact of the Omicron variant and provides valuable insights into its transmission dynamics, which can inform public health policies and strategies that are aimed at controlling its spread. Additionally, this study’s findings represent a significant step forward in understanding the COVID-19 pandemic and its evolving variants. The results of the simulation showcase the effectiveness of the proposed method and demonstrate its potential to advance the field of COVID-19 research. The COVID epidemic model is reformulated by using fractional derivatives in the Caputo sense. The existence and uniqueness of the proposed model are illustrated in the model, taking into account some results of fixed point theory. The stability analysis for the system is established by incorporating the Hyers–Ulam method. For numerical treatment and simulations, we apply the Haar wavelet collocation method. The parameter estimation for the recorded COVID-19 cases in Pakistan from 23 June 2022 to 23 August 2022 is presented.

1. Introduction

The global pandemic caused by SARS-CoV-2, a spike protein virus, has resulted in widespread infections of coronavirus. The severe acute respiratory syndrome (SARS) epidemic in China in 2003 and the Middle East respiratory syndrome (MERS) epidemic in Saudi Arabia in 2012 [1] were both coronavirus epidemics. Coronaviruses are a large family of viruses. The World Health Organization (WHO) first noted a SARS coronavirus 2 (COVID-19) outbreak in Wuhan, China, in December 2019. In March 2020, the WHO declared the outbreak to be a global pandemic. Over 6.6 million COVID-19 deaths [2] have been reported as of 8 January 2023, with over 659 million confirmed cases worldwide. Despite the fact that the disease is primarily transmitted through respiratory droplets produced by breathing, coughing, sneezing, and talking [3,4,5], more research indicates that it may also be transmitted through the air [6,7,8,9]. Additionally, COVID-19 can be acquired by coming into contact with contaminated objects. The virus can cause symptoms such as coughing, muscle pain, vertigo, high temperature, loss of smell, throat irritation, weakness, and nasal congestion and has an incubation period of 2–14 days.

The COVID-19 Omicron variant, which was first identified in South Africa in November 2021, has since raised serious concerns among public health authorities throughout the world, including in Pakistan. The COVID-19 virus, SARS-CoV-2, is thought to be more contagious and potentially more resistant to current vaccinations than earlier SARS-CoV-2 variants. The possibility of the virus spreading further in communities as a result has raised concerns, especially in regions with high population densities and low vaccination rates. The Omicron variant has been found in Pakistan’s major cities such as Karachi, Lahore, and Islamabad as well as other regions of the nation. In response to the threat posed by the variant, the government has increased testing and contact tracing efforts, increased vaccination efforts, and implemented stricter public health measures such as social seclusion and mask use. Despite these initiatives, the situation is still worrying because COVID-19 cases are still on the rise across much of the nation. This emphasizes the requirement for ongoing watchfulness and a thorough, multifaceted strategy to stop the virus’s spread. Along with public health initiatives, this also entails work intended to improve vaccination rates and access, as well as ongoing studies into the biology and epidemiology of the Omicron variant to better understand its potential effects on society. In the end, it will take a concerted effort from all facets of society, including the government, healthcare professionals, and the general public, to stop the spread of COVID-19 and its variants in Pakistan. By working together, it may be possible to reduce the impact of the Omicron variant and protect the health and well-being of communities in Pakistan.

Modeling the spread of the COVID-19 pandemic has been essential in guiding public health policies and interventions as it has been a significant global health crisis. Multiscale modeling techniques have become an effective method for understanding the intricate dynamics of virus pandemics and have been heavily utilized in the case of COVID-19 [10,11]. Epidemiological models have been extensively used to study the dynamics of COVID-19 transmission at the macroscopic level. These models span a spectrum ranging from straightforward compartmental models, such as the susceptible-infected-recovered (SIR) model, to more intricate models that take into account spatial and temporal heterogeneity, demographic factors, and interventions, such as vaccination and social isolation [12,13]. Additionally, the behavior of individual agents and their interactions during the COVID-19 outbreak have been simulated using agent-based models. Network models have been used to study the dissemination of COVID-19 within social networks and communities at the mesoscopic level. These models take into account the variation in contact patterns as well as the effects of social isolation policies on the dynamics of transmission. Large-scale datasets have also been analyzed to find patterns and trends in the spread of COVID-19 [14,15] using data-driven approaches, such as machine learning and artificial intelligence techniques. Models have been created to study the molecular interactions between the virus and host cells at the microscopic level. These models include the structural modeling of virus proteins and their interactions with host cell receptors (see, for example, [16,17]), as well as molecular dynamics simulations. These methods have been applied to comprehend the molecular mechanisms of virus replication and infection and to find potential therapeutic targets. Along with these conventional modeling techniques, there is growing interest in the use of multiscale models, which combine various modeling scales to capture the intricate relationships between virus spread, individual behavior, and public health policies. These models seek to offer a more thorough understanding of the COVID-19 dynamics and assist in the creation of efficient interventions and policies. Multiscale modeling has generally been shown to be an effective method for analyzing the COVID-19 epidemic and informing public health policies and interventions. Multiscale modeling techniques offer a promising way to deepen our understanding of this intricate global health crisis, even though there is still much to learn about the dynamics of COVID-19.

Fractional order differential equations (FODs) have been used to gain a deeper understanding of diseases, and mathematical models have been formulated and studied for various diseases [18,19,20,21,22,23]. The process of differentiation and integration in fractional calculus is generalized to non-integer orders, making it a valuable tool for research in various fields [24]. Caputo–Fabrizio (CF) operators [25], which use non-singular kernels, are one type of fractional derivative that has been developed to overcome the limitations of the ordinary operator. However, the CF operator has a locality problem, leading to the proposal of Mittag-Leffler kernels as a novel type of fractional derivative by Atangana and Baleanu [26]. Overall, the use of fractional calculus and mathematical models can provide insights into the dynamics of infectious diseases and inform public health policy. Omame et al. [27] explore the potential impact of COVID-19 on the dynamics of dengue and HIV transmission in their paper titled “Assessing the impact of SARS-CoV-2 infection on the dynamics of dengue and HIV via fractional derivatives”. Using fractional calculus, the study creates mathematical models and suggests that the pandemic could significantly affect the transmission and control of these diseases. The authors emphasize the importance of considering multiple infectious diseases during pandemics. In their paper “Backward bifurcation and optimal control in a co-infection model for SARS-CoV-2 and ZIKV” [28], the authors present a mathematical model for the co-infection of SARS-CoV-2 and ZIKV. The study examines the role of backward bifurcation in the dynamics of the co-infection and evaluates the potential benefits of optimal control strategies. The study highlights the importance of considering multiple infectious diseases and optimal control strategies in public health policy.

Furthermore, fractional operators with non-singular and singular kernels have been proposed in several works [29,30,31], and research related to these topics and their applications can be found in a number of recent publications [32,33,34,35,36]. As of late, a number of studies have appeared in mathematical modeling that investigate addressing social issues, such as criminal issues, using FCs. A time lag exists between the individual’s offense and the judgment, which is why Bansal et al. [37] introduced the time-delay coefficient to extend the proposed fractional-order crime transmission model to the delayed model. With regard to analyzing crime congestion, in their study, Pritam et al. [38] examined a mathematical model of crime transmission using a fractional-order derivative that includes memory effects, allowing for the previous input’s impact to be taken into account when forecasting the growth rate of crime. Using the iterative fractional-order Adams–Bashforth approach, ref. [39] found the approximate solution and numerically simulated it for various control strategies in different fractional orders.

In recent years, there has been a growing interest in the use of Haar wavelet numerical methods for solving problems related to the COVID-19 pandemic. These methods have been used to model the spread of the virus over time, taking into account the various factors that influence its transmission. The results of these studies have provided valuable insights into the dynamics of the pandemic and its evolution and have informed public health policies and strategies aimed at controlling its impact. Overall, Haar wavelet numerical methods have proven to be a valuable tool for solving problems with fractional derivatives and have been widely adopted in various fields due to their efficiency and accuracy. With ongoing research and development, they are expected to play an increasingly important role in solving complex problems in the future.

Similarly, wavelet analyses have been extensively applied in numerical analyses, statistical applications, image digital processing, quantum field theory, and many other fields. A wide range of applications has been made for Haar wavelets, including in communication and physics research as well as more mathematically-based research on differential equations and nonlinear problems [40]. There is an emphasis on Haar wavelets among all wavelet families. In mathematics, they are the simplest wavelet family because they consist of pairs of piecewise constant functions. In addition, Haar wavelets can also be integrated analytically in random times. Recently, researchers applied this technique for solving different fractional-order mathematical models [41,42]. In addition to being a fast method, the method is also more stable.

2. Preliminaries

The definition of the fractional derivative has developed greatly in recent years [43,44], ranging from the non-singular kernel derivatives and Riemann–Liouville (RL) fractional derivative without a singular kernel to the two-parameter derivative with non-singular and non-local kernels [45]. The two most commonly used definitions are described below:

Definition 1.

The fractional derivative of y of order δ is defined by Riemann–Liouville as follows [42]:

$${\mathbb{D}}_{*}^{\delta }\mathcal{F}(t)=\left\{\begin{array}{cc}\frac{1}{\mathsf{\Gamma }(s-\delta )}{\left(\frac{d}{dt}\right)}^s{\int }_0^t\frac{\mathcal{F}(\nu )}{{(t-\nu )}^{\delta -s+1}}d\nu ,\hfill & 0\le s-1<\delta <s,s\in \mathbb{N},\hfill \\ {(d/dt)}^s\mathcal{F}(t),\hfill & \delta =s,s\in \mathbb{N}.\hfill \end{array}\right.$$

(1)

Definition 2.

The Caputo fractional derivative of the function $\mathcal{F}$ of order δ is defined as follows [42]:

$${\mathbb{D}}_{*}^{\delta }\mathcal{F}(t)=\left\{\begin{array}{cc}\frac{1}{\mathsf{\Gamma }(s-\delta )}{\int }_0^t\frac{{(d/d\nu )}^s\mathcal{F}(\nu )}{{(t-\nu )}^{\delta -s+1}}d\nu ,\hfill & 0\le s-1<\delta <s,s\in \mathbb{N},\hfill \\ {(d/dt)}^s\mathcal{F}(t),\hfill & \delta =s,s\in \mathbb{N}.\hfill \end{array}\right.$$

(2)

Also used in this study is the RL form of the fractional integral operator ${\mathbb{D}}_{*}^{-\delta }$ of order $\delta$, which is defined as follows:

$${\mathbb{D}}_{*}^{-\delta }\mathcal{F}(t)=\frac{1}{\mathsf{\Gamma }(\delta )}{\int }_0^t\mathcal{F}(\nu ){(t-\nu )}^{\delta -1}d\nu$$

(3)

Haar Wavelets

If $\psi (t)$ and ${\tilde{\psi }}_0(t)$ represent the mother Haar wavelet function (on the real line) and Haar scaling function, respectively, then they are given by [46,47]:

$$\begin{array}{cc}\hfill \psi (t)& =\left\{\begin{array}{cc}1,\hfill & \mathrm{if}t\in \left[0,\frac{1}{2}\right),\hfill \\ -1,\hfill & \mathrm{if}t\in \left[\frac{1}{2},1\right),\hfill \\ 0,\hfill & \mathrm{elsewhere},\hfill \end{array}\right.\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\tilde{\psi }}_0(t)& =1,\mathrm{if}t\in [0,1).\hfill \end{array}$$

(5)

As a result, if the various Haar wavelets that are produced on the interval $[0,1)$ using multiresolution analysis are ${\tilde{\psi }}_m(t)$, then:

$${\tilde{\psi }}_m(t)=2^{j/2}\psi \left(2^{j}t-p\right),m=1,2,\dots ;$$

(6)

where $m=2^j+p:p=0,1,\dots ,2^j-1;j=0,1,\dots .$ Furthermore, we can translate the Haar functions on $u-1\le t<u$ as

$${\tilde{\psi }}_{u,m}(t)={\tilde{\psi }}_m(t+1-u),m=0,1,2,\dots ,u=1,2,\dots ,\varrho ,\varrho \in \mathbb{N}.$$

(7)

The resulting sequence ${\left\{{\tilde{\psi }}_m(t)\right\}}_{m=0}^{\infty }$ forms a complete orthonormal system [47] in ${\mathcal{L}}^2[0,1).$ Similarly, the sequence ${\left\{{\tilde{\psi }}_{u,m}(t)\right\}}_{m=0}^{\infty },u=1,2,\dots ,\varrho$, forms a complete orthonormal system in ${\mathcal{L}}^2[0,\varrho )$. Therefore, any function $\mathcal{F}(t)\in {\mathcal{L}}^2[0,\varrho )$ can be expanded in terms of Haar orthonormal basis functions as

$$\mathcal{F}(t)=\sum _{u=1}^{\varrho }\sum _{m=0}^{\infty }{\mathcal{G}}_{u,m}{\tilde{\psi }}_{u,m}(t).$$

(8)

Additionally, after truncating the series $\mathcal{F}(t)$, we obtain the equivalent approximation $y_p(t)$ of $\mathcal{F}(t)$ as

$$\mathcal{F}(t)\approx y_p(t)=\sum _{u=1}^{\varrho }\sum _{m=0}^{p-1}{\mathcal{G}}_{u,m}{\tilde{\psi }}_{u,m}(t)=B_{\varrho p\times 1}^T{\tilde{\psi }}_{\varrho p\times 1}(t),$$

(9)

where the coefficients ${\mathcal{G}}_{u,m}$ can be expressed by the inner product

$$\begin{array}{c}〈\mathcal{F}(t),{\tilde{\psi }}_{u,m}(t)〉={\int }_{u-1}^u\mathcal{F}(t){\tilde{\psi }}_{u,m}(t)dt,m=1,2,\dots ,(p-1),u=1,2,\dots ,\varrho ,\hfill \\ B_{\varrho p\times 1}={\left[{\mathcal{G}}_{1,0},\dots ,{\mathcal{G}}_{1,p-1},{\mathcal{G}}_{2,0},\dots ,{\mathcal{G}}_{2,p-1},\dots ,{\mathcal{G}}_{\varrho ,0},\dots ,{\mathcal{G}}_{\varrho ,p-1}\right]}^T,\hfill \\ {\tilde{\psi }}_{\varrho p\times 1}={\left[{\tilde{\psi }}_{1,0},\dots ,{\tilde{\psi }}_{1,p-1},{\tilde{\psi }}_{2,0},\dots ,{\tilde{\psi }}_{2,p-1},\dots ,{\tilde{\psi }}_{\varrho ,0},\dots ,{\tilde{\psi }}_{\varrho ,p-1}\right]}^T,\hfill \end{array}$$

(10)

and superscript T indicates the transpose of a matrix.

3. Mathematical Model

Understanding and forecasting the spread of infectious illnesses are greatly aided by mathematical modeling. Researchers can model how illnesses behave and spread within communities using mathematical equations and algorithms. As a result, they can forecast how the disease will develop in the future and can try different intervention techniques. In addition to helping to identify risk factors and guide public health policy, mathematical models can offer important insights into the fundamental mechanisms of disease transmission. It is crucial to keep in mind that the correctness of these models depends on the caliber of the data input and model assumptions, and they should always be utilized in conjunction with other information sources. Regarding the work of [6,7], the model used consists of the following ODEs:

$$\left\{\begin{array}{c}\frac{d\mathbf{S}(t)}{dt}=B-\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))-({\epsilon }_1+\rho +\eta )\mathbf{S}(t)+\mathcal{K}\mathbf{V}(t),\hfill \\ \frac{d\mathbf{V}(t)}{dt}=\pi \mathbf{S}(t)-\psi \mathbf{V}(t)-(\mathcal{K}+\rho )\mathbf{V}(t),\hfill \\ \frac{d\mathbf{E}(t)}{dt}=\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))+\psi \mathbf{V}(t)-({\epsilon }_2+\rho +\phi )\mathbf{E}(t),\hfill \\ \frac{d\mathbf{I}(t)}{dt}=\phi \mathbf{E}(t)-(\lambda +ϵ+\rho +{\epsilon }_3)\mathbf{I}(t),\hfill \\ \frac{d\mathbf{Q}(t)}{dt}={\epsilon }_1\mathbf{S}(t)+{\epsilon }_2\mathbf{E}(t)+{\epsilon }_3\mathbf{I}(t)-(\rho +\sigma )\mathbf{Q}(t),\hfill \\ \frac{d\mathbf{R}(t)}{dt}=\eta \mathbf{S}(t)+\sigma \mathbf{Q}(t)+\lambda \mathbf{I}(t)-\rho \mathbf{R}(t),\hfill \\ \mathbf{S}(t)\ge 0,\mathbf{V}(t)\ge 0,\mathbf{E}(t)\ge 0,\mathbf{I}(t)\ge 0,\mathbf{Q}(t)\ge 0,\mathbf{R}(t)\ge 0.\hfill \end{array}\right.$$

(11)

A representation of COVID-19 transmission is provided by the system of ordinary differential equations in mathematical model refe1; the model is divided into six compartments, each representing a different group of people in the population that is affected by the disease. These compartments are as follows: susceptible people $\mathbf{S}(t)$, vaccinated people $\mathbf{V}(t)$, exposed people $\mathbf{E}(t)$, infectious people $\mathbf{I}(t)$, quarantined people $\mathbf{Q}(t)$, and recovered people $\mathbf{R}(t)$. Susceptible people have not been infected but can contract the virus, whereas vaccinated people have received a vaccine and are immune to the virus. Exposed people have been infected but have not yet developed symptoms, whereas infectious people are infected and can spread the virus. Individuals quarantined have been identified as infected and are being isolated to prevent further transmission. Finally, recovered individuals have recovered from the disease and are now immune.

Numerous parameters that affect the dynamics of the disease are included in the model for predicting the spread of COVID-19. The following parameters make up this list: birth rate (B), transmission rate ($\theta$), infectivity coefficient ($\tau$), quarantine rates for exposed and susceptible individuals (${\epsilon }_2$) and (${\epsilon }_1$) and for natural deaths ($\rho$), vaccination efficacy rate ($\eta$), vaccination rate ($\mathcal{K}$), vaccination loss rate ($\psi$), rate of susceptibility to vaccination ($\pi$), progression rate ($\phi$), recovery rate ($\lambda$), disease-induced death rate ($ϵ$), quarantine rate for infectious individuals (${\epsilon }_3$), and quarantine loss rate ($\sigma$). The model’s differential equations predict how these parameters and variables will change over time to affect how many people are in each compartment. It is possible to forecast the spread of COVID-19 and assess the effects of different interventions, such as vaccination and quarantine, by simulating the model.

3.1. Formulation of Fractional Model

The ability to model with fractional-order derivatives, which can more precisely capture the dynamics of some diseases that exhibit non-integer order behavior, is one benefit of using the Caputo derivative over the classical derivative in the context of disease modeling. Because the fractional-order derivatives can capture both memory effects and the power-law decay characteristic of many disease models, this can result in more accurate predictions and better control strategies. The Caputo derivative is also a more economical option for disease modeling because it frequently requires fewer computational resources and data. Additionally, the majority of natural phenomena, including epidemiological dynamics, exhibit the time memory effect. Model (11) is expressed in integral form as:

$$\left\{\begin{array}{c}\frac{d\mathbf{S}(t)}{dt}={\int }_{t_0}^t\varsigma (t-\varpi )\left[B-\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))-({\epsilon }_1+\rho +\eta )\mathbf{S}(t)+\mathcal{K}\mathbf{V}(t)\right]d\varpi ,\hfill \\ \frac{d\mathbf{V}(t)}{dt}={\int }_{t_0}^t\varsigma (t-\varpi )\left[\pi \mathbf{S}(t)-\psi \mathbf{V}(t)-(\mathcal{K}+\rho )\mathbf{V}(t)\right]d\varpi ,\hfill \\ \frac{d\mathbf{E}(t)}{dt}={\int }_{t_0}^t\varsigma (t-\varpi )\left[\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))+\psi \mathbf{V}(t)-({\epsilon }_2+\rho +\phi )\mathbf{E}(t)\right]d\varpi ,\hfill \\ \frac{d\mathbf{I}(t)}{dt}={\int }_{t_0}^t\varsigma (t-\varpi )\left[\phi \mathbf{E}(t)-(\lambda +ϵ+\rho +{\epsilon }_3)\mathbf{I}(t)\right]d\varpi ,\hfill \\ \frac{d\mathbf{Q}(t)}{dt}={\int }_{t_0}^t\varsigma (t-\varpi )\left[{\epsilon }_1\mathbf{S}(t)+{\epsilon }_2\mathbf{E}(t)+{\epsilon }_3\mathbf{I}(t)-(\rho +\sigma )\mathbf{Q}(t)\right]d\varpi ,\hfill \\ \frac{d\mathbf{R}(t)}{dt}={\int }_{t_0}^t\varsigma (t-⫋)\left[\eta \mathbf{S}(t)+\sigma \mathbf{Q}(t)+\lambda \mathbf{I}(t)-\rho \mathbf{R}(t)\right]d\varpi .\hfill \end{array}\right.$$

(12)

Incorporating the Caputo derivative, we get

$$\left\{\begin{array}{c}{}^C{D}_t^{\delta -1}\left[\frac{d\mathbf{S}(t)}{dt}\right]{=}^C{D}_t^{\delta -1}{I}^{-(\delta -1)}\left[B-\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))-({\epsilon }_1+\rho +\eta )\mathbf{S}(t)+\mathcal{K}\mathbf{V}(t)\right],\hfill \\ {}^C{D}_t^{\delta -1}\left[\frac{d\mathbf{V}(t)}{dt}\right]{=}^C{D}_t^{\delta -1}{I}^{-(\delta -1)}\left[\pi \mathbf{S}(t)-\psi \mathbf{V}(t)-(\mathcal{K}+\rho )\mathbf{V}(t)\right],\hfill \\ {}^C{D}_t^{\delta -1}\left[\frac{d\mathbf{E}(t)}{dt}\right]{=}^C{D}_t^{\delta -1}{I}^{-(\delta -1)}\left[\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))+\psi \mathbf{V}(t)-({\epsilon }_2+\rho +\phi )\mathbf{E}(t)\right],\hfill \\ {}^C{D}_t^{\delta -1}\left[\frac{d\mathbf{I}(t)}{dt}\right]{=}^C{D}_t^{\delta -1}{I}^{-(\delta -1)}\left[\phi \mathbf{E}(t)-(\lambda +ϵ+\rho +{\epsilon }_3)\mathbf{I}(t)\right],\hfill \\ {}^C{D}_t^{\delta -1}\left[\frac{d\mathbf{Q}(t)}{dt}\right]{=}^C{D}_t^{\delta -1}{I}^{-(\delta -1)}\left[{\epsilon }_1\mathbf{S}(t)+{\epsilon }_2\mathbf{E}(t)+{\epsilon }_3\mathbf{I}(t)-(\rho +\sigma )\mathbf{Q}(t)\right],\hfill \\ {}^C{D}_t^{\delta -1}\left[\frac{d\mathbf{R}(t)}{dt}\right]{=}^C{D}_t^{\delta -1}{I}^{-(\delta -1)}\left[\eta \mathbf{S}(t)+\sigma \mathbf{Q}(t)+\lambda \mathbf{I}(t)-\rho \mathbf{R}(t)\right].\hfill \end{array}\right.$$

(13)

After calculations, we reach

$$\left\{\begin{array}{c}{}^C{D}_t^{\delta }\mathbf{S}(t)=B-\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))-({\epsilon }_1+\rho +\eta )\mathbf{S}(t)+\mathcal{K}\mathbf{V}(t),\hfill \\ {}^C{D}_t^{\delta }\mathbf{V}(t)=\pi \mathbf{S}(t)-\psi \mathbf{V}(t)-(\mathcal{K}+\rho )\mathbf{V}(t),\hfill \\ {}^C{D}_t^{\delta }\mathbf{E}(t)=\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))+\psi \mathbf{V}(t)-({\epsilon }_2+\rho +\phi )\mathbf{E}(t),\hfill \\ {}^C{D}_t^{\delta }\mathbf{I}(t)=\phi \mathbf{E}(t)-(\lambda +ϵ+\rho +{\epsilon }_3)\mathbf{I}(t),\hfill \\ {}^C{D}_t^{\delta }\mathbf{Q}(t)={\epsilon }_1\mathbf{S}(t)+{\epsilon }_2\mathbf{E}(t)+{\epsilon }_3\mathbf{I}(t)-(\rho +\sigma )\mathbf{Q}(t),\hfill \\ {}^C{D}_t^{\delta }\mathbf{R}(t)=\eta \mathbf{S}(t)+\sigma \mathbf{Q}(t)+\lambda \mathbf{I}(t)-\rho \mathbf{R}(t).\hfill \end{array}\right.$$

(14)

3.2. Basic Reproductive Number $R_0$

The basic reproductive number $R_0$, which expresses the typical number of secondary infections caused by a single infected person in a population that is fully susceptible, is an essential component of epidemiological modeling. The advantages of $R_0$ in epidemiological modeling are extensive: it offers a precise measurement of a disease’s transmissibility, where a disease is more contagious if its $R_0$ is high, whereas a low $R_0$ suggests that it is less contagious; predicting the potential spread of an outbreak, which helps epidemiologists forecast the size and length of an outbreak and create the most efficient control strategies by calculating $R_0$; it can be employed to assess the efficacy of interventions, and epidemiologists can assess the success of an intervention in reducing transmission by comparing the $R_0$ before and after the implementation of a vaccine or quarantine measures; it is useful for locating the crucial control points, and the amount of infected people in a population, for example, is one of the crucial factors in the transmission of an infection that $R_0$ can be used to pinpoint; and the herd immunity threshold can be predicted using the number, and the percentage of the population that must be immune to a disease for there to be herd immunity is known as the herd immunity threshold. The herd immunity threshold and subsequently the overall efficacy of vaccination programs are predicted with the aid of $R_0$. Overall, $R_0$ is a useful tool for understanding the spread of infectious diseases and for designing effective public health interventions.

The DFE of Model (14) is denoted by $E^0=({\mathbf{S}}_0,{\mathbf{V}}_0,0,0,{\mathbf{Q}}_0,{\mathbf{R}}_0)$, where

$${\mathbf{S}}_0=\frac{B}{\eta +\rho +{\epsilon }_1},{\mathbf{V}}_0=\frac{\pi {\mathbf{S}}_0}{\mathcal{K}+\rho },{\mathbf{Q}}_0=\frac{{\epsilon }_1{\mathbf{S}}_0}{\rho +\tau },{\mathbf{R}}_0=\frac{\eta {\mathbf{S}}_0+\tau {\mathbf{Q}}_0}{\rho }.$$

Our proposed model is split into two matrices [48].

$$\tilde{U}=\left[\begin{array}{cc}0& \theta {\mathbf{S}}_0\\ 0& 0\end{array}\right],\tilde{V}=\left[\begin{array}{cc}\phi +\rho +{\epsilon }_2& 0\\ -\phi & \rho +ϵ+\lambda +{\epsilon }_3\end{array}\right],$$

$${\tilde{V}}^{-1}=\frac{1}{(\phi +\rho +{\epsilon }_2)(\rho +ϵ+\lambda +{\epsilon }_3)}\left[\begin{array}{cc}\rho +ϵ+\lambda +{\epsilon }_3& 0\\ \phi & \phi +\rho +{\epsilon }_2\end{array}\right],$$

$$\tilde{U}{\tilde{V}}^{-1}=\left[\begin{array}{cc}\frac{\theta {\mathbf{S}}_0\phi }{(\phi +\rho +{\epsilon }_2)(\rho +ϵ+\lambda +{\epsilon }_3)}& \frac{\theta {\mathbf{S}}_0}{\rho +ϵ+\lambda +{\epsilon }_3}\\ 0& 0\end{array}\right].$$

Hence,

$$R_0=\frac{\phi \theta B}{(\eta +\rho +{\epsilon }_1)(\phi +\rho +{\epsilon }_2)(\rho +ϵ+\lambda +{\epsilon }_3)}.$$

4. Existence and Uniqueness

The solution for the system (14) using the Caputo operator will be described below, along with its existence and uniqueness. Assume that the continuous real-valued function $\mathcal{A}(Y)$, which has the sup-norm property, is a Banach space on $J=[0,b]$ and that $Y=[0,\kappa ]$ and $P=\mathcal{A}(Y)\times \mathcal{A}(Y)\times \mathcal{A}(Y)\times \mathcal{A}(Y)\times \mathcal{A}(Y)$ with norm $∥\left(\mathbf{S},\mathbf{V},\mathbf{E},\mathbf{Q},\mathbf{I},\mathbf{R}\right)∥=\parallel \mathbf{S}\parallel +\parallel \mathbf{V}\parallel +\parallel \mathbf{E}\parallel +∥\mathbf{Q}∥+∥\mathbf{I}∥+\parallel \mathbf{R}\parallel$, where $\parallel \mathbf{S}\parallel ={sup}_{t\in Y}\left|\mathbf{S}(t)\right|,\parallel \mathbf{V}(t)\parallel ={sup}_{t\in Y}\left|\mathbf{V}(t)\right|,\parallel \mathbf{E}(t)\parallel ={sup}_{t\in Y}\left|\mathbf{E}(t)\right|,∥\mathbf{Q}∥={sup}_{t\in Y}\left|\mathbf{Q}(t)\right|,∥\mathbf{I}∥={sup}_{t\in Y}\left|\mathbf{I}(t)\right|,∥\mathbf{R}∥={sup}_{t\in Y}\left|\mathbf{R}(t)\right|$. The following equation is obtain by using the Caputo fractional integral operator on both sides of (14):

$$\left\{\begin{array}{c}\mathbf{S}(t)-\mathbf{S}(0){=}^C{\mathbb{D}}_{0,t}^{\delta }\mathbf{S}(t)\left\{B-\theta \mathbf{S}\mathbf{I}(t)(1+\tau \mathbf{I})-({\epsilon }_1+\rho +\eta )\mathbf{S}(t)+\mathcal{K}\mathbf{V}(t)\right\},\hfill \\ \mathbf{V}(t)-\mathbf{V}(0){=}^C{\mathbb{D}}_{0,t}^{\delta }\mathbf{V}(t)\left\{\pi \mathbf{S}(t)-\psi \mathbf{V}(t)-(\mathcal{K}+\rho )\mathbf{V}(t)\right\},\hfill \\ \mathbf{E}(t)-\mathbf{E}(0){=}^C{\mathbb{D}}_{0,t}^{\delta }\mathbf{E}(t)\left\{\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))+\psi \mathbf{V}(t)-({\epsilon }_2+\rho +\phi )\mathbf{E}(t)\right\},\hfill \\ \mathbf{I}(t)-\mathbf{I}(0){=}^C{\mathbb{D}}_{0,t}^{\delta }\mathbf{I}(t)\left\{\phi \mathbf{E}(t)-(\lambda +ϵ+\rho +{\epsilon }_3)\mathbf{I}(t)\right\},\hfill \\ \mathbf{Q}(t)-\mathbf{Q}(0){=}^C{\mathbb{D}}_{0,t}^{\delta }\mathbf{Q}(t)\left\{{\epsilon }_1\mathbf{S}(t)+{\epsilon }_2\mathbf{E}(t)+{\epsilon }_3\mathbf{I}(t)-(\rho +\sigma )\mathbf{Q}(t)\right\},\hfill \\ \mathbf{R}(t)-\mathbf{R}(0){=}^C{\mathbb{D}}_{0,t}^{\delta }\mathbf{R}(t)\left\{\eta \mathbf{S}(t)+\sigma \mathbf{Q}(t)+\lambda \mathbf{I}(t)-\rho \mathbf{R}(t)\right\}.\hfill \end{array}\right.$$

(15)

After calculation,

$$\begin{array}{c}\mathbf{S}(t)-\mathbf{S}(0)=\mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }{\mathcal{B}}_1(\delta ,\varpi ,\mathbf{S}(\varpi ))d\varpi ,\hfill \\ \\ \mathbf{V}(t)-\mathbf{V}(0)=\mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }{\mathcal{B}}_2(\delta ,\varpi ,\mathbf{V}(t)(\varpi ))d\varpi ,\hfill \\ \\ \mathbf{E}(t)-\mathbf{E}(0)=\mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }{\mathcal{B}}_2(\delta ,\varpi ,\mathbf{E}(t)(\varpi ))d\varpi ,\hfill \\ \\ \mathbf{I}(t)-\mathbf{I}(0)=\mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }{\mathcal{B}}_4\left(\delta ,\varpi ,\mathbf{I}(t)(\varpi )\right)d\varpi ,\hfill \\ \\ \mathbf{Q}(t)-\mathbf{Q}(0)=\mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }{\mathcal{B}}_3\left(\delta ,\varpi ,\mathbf{Q}(t)(\varpi )\right)d\varpi ,\hfill \\ \\ \mathbf{R}(t)-\mathbf{R}(0)=\mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }{\mathcal{B}}_5(\delta ,\varpi ,\mathbf{R}(t)(\varpi ))d\varpi ,\hfill \end{array}$$

(16)

where

$$\left\{\begin{array}{c}\mathbf{S}(t)-\mathbf{S}(0){=}^C{\mathbb{D}}_{0,t}^{\delta }\mathbf{S}(t)\left\{B-\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))-({\epsilon }_1+\rho +\eta )\mathbf{S}(t)+\mathcal{K}\mathbf{V}(t)\right\},\hfill \\ \mathbf{V}(t)-\mathbf{V}(0){=}^C{\mathbb{D}}_{0,t}^{\delta }\mathbf{V}(t)\left\{\pi \mathbf{S}(t)-\psi \mathbf{V}(t)-(\mathcal{K}+\rho )\mathbf{V}(t)\right\},\hfill \\ \mathbf{E}(t)-\mathbf{E}(0){=}^C{\mathbb{D}}_{0,t}^{\delta }\mathbf{E}(t)\left\{\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))+\psi \mathbf{V}(t)-({\epsilon }_2+\rho +\phi )\mathbf{E}(t)\right\},\hfill \\ \mathbf{I}(t)-\mathbf{I}(0){=}^C{\mathbb{D}}_{0,t}^{\delta }\mathbf{I}(t)\left\{\phi \mathbf{E}(t)-(\lambda +ϵ+\rho +{\epsilon }_3)\mathbf{I}(t)\right\},\hfill \\ \mathbf{Q}(t)-\mathbf{Q}(0){=}^C{\mathbb{D}}_{0,t}^{\delta }\mathbf{Q}(t)\left\{{\epsilon }_1\mathbf{S}(t)+{\epsilon }_2\mathbf{E}(t)+{\epsilon }_3\mathbf{I}(t)-(\rho +\sigma )\mathbf{Q}(t)\right\},\hfill \\ \mathbf{R}(t)-\mathbf{R}(0){=}^C{\mathbb{D}}_{0,t}^{\delta }\mathbf{R}(t)\left\{\eta \mathbf{S}(t)+\sigma \mathbf{Q}(t)+\lambda \mathbf{I}(t)-\rho \mathbf{R}(t)\right\}.\hfill \end{array}\right.$$

(17)

$$\begin{array}{c}{\mathcal{B}}_1(\delta ,t,\mathbf{S}(t))=B-\theta \mathbf{S}\mathbf{I}(t)(1+\tau \mathbf{I}(t))-({\epsilon }_1+\rho +\eta )\mathbf{S}(t)+\mathcal{K}\mathbf{V}(t),\hfill \\ \\ {\mathcal{B}}_2(\delta ,t,\mathbf{V}(t))=\pi \mathbf{S}(t)-\psi \mathbf{V}(t)-(\mathcal{K}+\rho )\mathbf{V}(t),\hfill \\ \\ {\mathcal{B}}_2(\delta ,t,\mathbf{E}(t))=\theta \mathbf{S}\mathbf{I}(t)(1+\tau \mathbf{I}(t))+\psi \mathbf{V}(t)-({\epsilon }_2+\rho +\phi )\mathbf{E}(t),\hfill \\ \\ {\mathcal{B}}_3\left(\delta ,t,\mathbf{I}(t)\right)=\phi \mathbf{E}(t)-(\lambda +ϵ+\rho +{\epsilon }_3)\mathbf{I}(t),\hfill \\ \\ {\mathcal{B}}_4\left(\delta ,t,\mathbf{Q}(t)\right)={\epsilon }_1\mathbf{S}(t)+{\epsilon }_2\mathbf{E}(t)+{\epsilon }_3\mathbf{I}(t)-(\rho +\sigma )\mathbf{Q}(t),\hfill \\ \\ {\mathcal{B}}_5(\delta ,t,\mathbf{R}(t))=\eta \mathbf{S}(t)+\sigma \mathbf{Q}(t)+\lambda \mathbf{I}(t)-\rho \mathbf{R}(t).\hfill \end{array}$$

(18)

If $\mathbf{S}(t),\mathbf{V}(t),\mathbf{E}(t),\mathbf{Q}(t),\mathbf{I}(t)$, and $\mathbf{R}(t)$ have an upper bound, then the symbols ${\mathcal{B}}_1,{\mathcal{B}}_2,{\mathcal{B}}_3$, ${\mathcal{B}}_4,{\mathcal{B}}_5$, and ${\mathcal{B}}_6$ are necessary for the Lipschitz condition. It should be noted that $\mathbf{S}(t)$ and ${\mathbf{S}}^{*}(t)$ are paired functions, and we reach

$$\begin{array}{c}∥{\mathcal{B}}_1(\delta ,t,\mathbf{S}(t))-{\mathcal{B}}_1\left(\delta ,t,{\mathbf{S}}^{*}(t)\right)∥=∥-\left(\theta \mathbf{I}(t)(1+\nsim \mathbf{I}(t))+\eta +\rho +{\epsilon }_1\right)\left(\mathbf{S}(t)-{\mathbf{S}}^{*}(t)\right)∥.\hfill \end{array}$$

(19)

Taking into account ${\mathsf{\Lambda }}_1=∥-\left(\theta \mathbf{I}(t)(1+\tau \mathbf{I}(t))+\eta +\rho +{\epsilon }_1\right)∥$, one reaches

$$∥{\mathcal{B}}_1(\delta ,t,\mathbf{S}(t))-{\mathcal{B}}_1\left(\delta ,t,{\mathbf{S}}^{*}(t)\right)∥\le {\mathsf{\Lambda }}_1∥\mathbf{S}(t)-{\mathbf{S}}^{*}(t)∥.$$

(20)

Similarly,

$$\begin{array}{c}∥{\mathcal{B}}_2(\delta ,t,\mathbf{V}(t))-{\mathcal{B}}_2\left(\delta ,t,{\mathbf{V}}^{*}(t)\right)∥\le {\mathsf{\Lambda }}_2∥\mathbf{V}(t)-{\mathbf{V}}^{*}(t)∥,\hfill \\ \\ ∥{\mathcal{B}}_3(\delta ,t,\mathbf{E}(t))-{\mathcal{B}}_3\left(\delta ,t,{\mathbf{E}}^{*}(t)\right)∥\le {\mathsf{\Lambda }}_3∥\mathbf{E}(t)-{\mathbf{E}}^{*}(t)∥,\hfill \\ \\ ∥{\mathcal{B}}_4\left(\delta ,t,\mathbf{Q}(t)\right)-{\mathcal{B}}_4\left(\delta ,t,{\mathbf{Q}}^{*}(t)\right)∥\le {\mathsf{\Lambda }}_4∥\mathbf{Q}(t)-{\mathbf{Q}}^{*}(t)∥,\hfill \\ \\ ∥{\mathcal{B}}_5\left(\delta ,t,\mathbf{I}(t)\right)-{\mathcal{B}}_5\left(\delta ,t,{\mathbf{I}}^{*}(t)\right)∥\le {\mathsf{\Lambda }}_5∥\mathbf{I}(t)-{\mathbf{I}}^{*}(t)∥,\hfill \\ \\ ∥{\mathcal{B}}_6(\delta ,t,\mathbf{R}(t))-{\mathcal{B}}_6\left(\delta ,t,{\mathbf{R}}^{*}(t)\right)∥\le {\mathsf{\Lambda }}_6∥\mathbf{R}(t)-{\mathbf{R}}^{*}(t)∥.\hfill \end{array}$$

(21)

where

$$\begin{array}{cc}\hfill {\mathsf{\Lambda }}_3=& ∥-\left(\phi +\rho +{\epsilon }_2\right)∥\hfill \\ & \\ \hfill {\mathsf{\Lambda }}_4=& ∥-\left(\rho +ϵ+\lambda +{\epsilon }_3\right)∥\hfill \\ & \\ \hfill {\mathsf{\Lambda }}_5=& ∥-\left(\rho +\tau \right)∥\hfill \\ & \\ \hfill {\mathsf{\Lambda }}_6=& ∥-\left(\rho \right)∥.\hfill \end{array}$$

This indicates that for each of the five functions, the Lipschitz condition is true. Recursively applying the expressions in (16), we obtain

$$\begin{array}{c}{\mathbf{S}}_n(t)=\mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }{\mathcal{B}}_1\left(\delta ,\varpi ,{\mathbf{S}}_{n-1}(\varpi )\right)d\varpi ,\hfill \\ \\ {\mathbf{V}(t)}_n(t)=\mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }{\mathcal{B}}_2\left(\delta ,\varpi ,{\mathbf{V}(t)}_{n-1}(\varpi )\right)d\varpi ,\hfill \\ \\ {\mathbf{E}(t)}_n(t)=\mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }{\mathcal{B}}_3\left(\delta ,\varpi ,{\mathbf{E}(t)}_{n-1}(\varpi )\right)d\varpi ,\hfill \\ \\ {\mathbf{Q}(t)}_n(t)=\mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }{\mathcal{B}}_4\left(\delta ,\varpi ,{\mathbf{Q}(t)}_{n-1}(\varpi )\right)d\varpi ,\hfill \\ \\ {\mathbf{I}(t)}_n(t)=\mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }{\mathcal{B}}_5\left(\delta ,\varpi ,{\mathbf{I}(t)}_{n-1}(\varpi )\right)d\varpi ,\hfill \\ \\ {\mathbf{R}(t)}_n(t)=\mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }{\mathcal{B}}_6\left(\delta ,\varpi ,{\mathbf{R}(t)}_{n-1}(\varpi )\right)d\varpi ,\hfill \end{array}$$

(22)

We can obtain the successive terms difference by considering ${\mathbf{S}}_0(t)={\mathbf{S}}_0,{\mathbf{V}}_0(t)={\mathbf{V}}_0,{\mathbf{E}}_0(t)={\mathbf{E}}_0,{\mathbf{Q}}_0(t)={\mathbf{Q}}_0,{\mathbf{I}}_0(t)={\mathbf{I}}_0$ and ${\mathbf{R}}_0(t)={\mathbf{R}}_0$ in conjunction with other relevant information.

$$\begin{array}{cc}\hfill {\mathsf{\Xi }}_{\mathbf{S},n}(t)=& {\mathbf{S}}_n(t)-{\mathbf{S}}_{n-1}(t)\hfill \\ \hfill =& \mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }\left({\mathcal{B}}_1\left(\delta ,\varpi ,{\mathbf{S}}_{n-1}(\varpi )\right)\right.\left.-{\mathcal{B}}_1\left(\delta ,\varpi ,{\mathbf{S}}_{n-2}(\varpi )\right)\right)d\varpi ,\hfill \\ \hfill {\mathsf{\Xi }}_{\mathbf{V},n}(t)=& {\mathbf{V}}_n(t)-{\mathbf{V}}_{n-1}(t)\hfill \\ \hfill =& \mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }\left({\mathcal{B}}_2\left(\delta ,\varpi ,{\mathbf{V}}_{n-1}(\varpi )\right)\right.\left.-{\mathcal{B}}_2\left(\delta ,\varpi ,{\mathbf{V}}_{n-2}(\varpi )\right)\right)d\varpi ,\hfill \\ \hfill {\mathsf{\Xi }}_{\mathbf{E}(t),n}(t)=& {\mathbf{E}(t)}_n(t)-{\mathbf{E}(t)}_{n-1}(t)\hfill \\ \hfill =& \mathcal{H}(\delta ){\int }_0^l{(t-\varpi )}^{-\delta }\left({\mathcal{B}}_3\left(\delta ,\varpi ,{\mathbf{E}(t)}_{n-1}(\varpi )\right)\right.\left.-{\mathcal{B}}_3\left(\delta ,\varpi ,{\mathbf{E}(t)}_{n-2}(\varpi )\right)\right)d\varpi ,\hfill \\ \hfill {\mathsf{\Xi }}_{\mathbf{I}(t),n}(t)=& {\mathbf{I}(t)}_{2n}(t)-{\mathbf{I}(t)}_{n-1}(t)\hfill \\ \hfill =& \mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }\left({\mathcal{B}}_4\left(\delta ,\varpi ,{\mathbf{I}(t)}_{n-1}(\varpi )\right)\right.\left.-{\mathcal{B}}_4\left(\delta ,\varpi ,{\mathbf{I}(t)}_{n-2}(\varpi )\right)\right)d\varpi ,\hfill \\ \hfill {\mathsf{\Xi }}_{\mathbf{Q}(t),n}(t)=& {\mathbf{Q}(t)}_{1n}(t)-{\mathbf{Q}(t)}_{n-1}(t)\hfill \\ \hfill =& \mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }\left({\mathcal{B}}_5\left(\delta ,\varpi ,{\mathbf{Q}(t)}_{n-1}(\varpi )\right)\right.\left.-{\mathcal{B}}_5\left(\delta ,\varpi ,{\mathbf{Q}(t)}_{n-2}(\varpi )\right)\right)d\varpi ,\hfill \\ \hfill {\mathsf{\Xi }}_{\mathbf{R}(t),n}(t)=& {\mathbf{R}(t)}_n(t)-{\mathbf{R}(t)}_{n-1}(t)\hfill \\ \hfill =& \mathcal{H}(\delta ){\int }_0^t{(t-\varpi )}^{-\delta }\left({\mathcal{B}}_6\left(\delta ,\varpi ,{\mathbf{R}(t)}_{n-1}(\varpi )\right)\right.\left.-{\mathcal{B}}_6\left(\delta ,\varpi ,{\mathbf{R}(t)}_{n-2}(\varpi )\right)\right)d\varpi .\hfill \end{array}$$

(23)

It is vital to observe that

$$\begin{array}{cc}\hfill {\mathbf{S}}_n(t)& =\sum _{m=0}^n{\mathsf{\Xi }}_{\mathbf{S},m}(t),{\mathbf{V}}_n(t)=\sum _{m=0}^n{\mathsf{\Xi }}_{\mathbf{V},m}(t),{\mathbf{V}}_n(t)=\sum _{m=0}^n{\mathsf{\Xi }}_{\mathbf{E},m}(t),\hfill \\ \hfill {\mathbf{Q}}_n(t)& =\sum _{m=0}^n{\mathsf{\Xi }}_{\mathbf{Q},m}(t),{\mathbf{I}}_n(t)=\sum _{m=0}^n{\mathsf{\Xi }}_{\mathbf{I},m}(t),{\mathbf{R}}_n(t)=\sum _{m=0}^n{\mathsf{\Xi }}_{\mathbf{R},m}(t).\hfill \end{array}$$

Additionally, by using Equations (20) and (21) and considering that

$$\begin{array}{cc}& {\mathsf{\Xi }}_{\mathbf{S},n-1}(t)={\mathbf{S}}_{n-1}(t)-{\mathbf{S}}_{n-2}(t),{\mathsf{\Xi }}_{\mathbf{V},n-1}(t)={\mathbf{V}}_{n-1}(t)-{\mathbf{V}}_{n-2}(t),{\mathsf{\Xi }}_{\mathbf{E},n-1}(t)={\mathbf{E}}_{n-1}(t)-{\mathbf{E}}_{n-2}(t),\hfill \\ & {\mathsf{\Xi }}_{\mathbf{Q},n-1}(t)={\mathbf{Q}}_{n-1}(t)-{\mathbf{Q}}_{n-2}(t),{\mathsf{\Xi }}_{\mathbf{I},n-1}(t)={\mathbf{I}}_{n-1}(t)-{\mathbf{I}}_{n-2}(t),{\mathsf{\Xi }}_{\mathbf{R},n-1}(t)={\mathbf{R}(t)}_{n-1}(t)-{\mathbf{R}}_{n-2}(t),\hfill \end{array}$$

we reach

$$\begin{array}{c}∥{\mathsf{\Xi }}_{\mathbf{S},n}(t)∥\le \mathcal{H}(\delta ){\mathsf{\Lambda }}_1{\int }_0^t{(t-\varpi )}^{-\delta }∥{\mathsf{\Xi }}_{\mathbf{S},n-1}(\varpi )∥d\varpi ,\hfill \\ \\ ∥{\mathsf{\Xi }}_{\mathbf{V},n}(t)∥\le \mathcal{H}(\delta ){\mathsf{\Lambda }}_2{\int }_0^t{(t-\varpi )}^{-\delta }∥{\mathsf{\Xi }}_{\mathbf{V},n-1}(\varpi )∥d\varpi ,\hfill \\ \\ ∥{\mathsf{\Xi }}_{\mathbf{E},n}(t)∥\le \mathcal{H}(\delta ){\mathsf{\Lambda }}_3{\int }_0^t{(t-\varpi )}^{-\delta }∥{\mathsf{\Xi }}_{\mathbf{E},n-1}(\varpi )∥d\varpi ,\hfill \\ \\ ∥{\mathsf{\Xi }}_{\mathbf{I},n}(t)∥\le \mathcal{H}(\delta ){\mathsf{\Lambda }}_4{\int }_0^t{(t-\varpi )}^{-\delta }∥{\mathsf{\Xi }}_{\mathbf{I},n-1}(\varpi )∥d\varpi ,\hfill \\ \\ ∥{\mathsf{\Xi }}_{\mathbf{Q},n}(t)∥\le \mathcal{H}(\delta ){\mathsf{\Lambda }}_5{\int }_0^t{(t-\varpi )}^{-\delta }∥{\mathsf{\Xi }}_{\mathbf{Q},n-1}(\varpi )∥d\varpi ,\hfill \\ \\ ∥{\mathsf{\Xi }}_{\mathbf{R},n}(t)∥\le \mathcal{H}(\delta ){\mathsf{\Lambda }}_6{\int }_0^t{(t-\varpi )}^{-\delta }∥{\mathsf{\Xi }}_{\mathbf{R},n-1}(\varpi )∥d\varpi .\hfill \end{array}$$

(24)

Theorem 1.

If the following condition holds,

$$\frac{\mathcal{H}(\delta )}{\delta }{\kappa }^{\delta }{\mathsf{\Lambda }}_m<1,m=1,2,\dots ,6.$$

(25)

Then, (14) has a unique solution for $t\in [0,\kappa ]$.

Proof. 

It has been demonstrated that the functions $\mathbf{S}(t),\mathbf{V}(t),\mathbf{E}(t),\mathbf{Q}(t),\mathbf{I}(t)$, and $\mathbf{R}(t)$ are bounded. Furthermore, it can be observed from Equations (20) and (21) that the symbols ${\mathcal{B}}_1,{\mathcal{B}}_2,{\mathcal{B}}_3,{\mathcal{B}}_4,{\mathcal{B}}_5$, and ${\mathcal{B}}_6$ are applicable to the Lipschitz condition. Thus, by using Equation (24) along with a recursive hypothesis, we can derive:

$$\begin{array}{c}∥{\mathsf{\Xi }}_{\mathbf{S},n}(t)∥\le ∥{\mathbf{S}}_0(t)∥{\left(\frac{\mathcal{H}(\delta )}{\delta }{\kappa }^{\delta }{\mathsf{\Lambda }}_1\right)}^n,\hfill \\ \\ ∥{\mathsf{\Xi }}_{\mathbf{V},n}(t)∥\le ∥{\mathbf{V}}_0(t)∥{\left(\frac{\mathcal{H}(\delta )}{\delta }{\kappa }^{\delta }{\mathsf{\Lambda }}_2\right)}^n,\hfill \\ \\ ∥{\mathsf{\Xi }}_{\mathbf{E},n}(t)∥\le ∥{\mathbf{E}}_0(t)∥{\left(\frac{\mathcal{H}(\delta )}{\delta }{\kappa }^{\delta }{\mathsf{\Lambda }}_3\right)}^n,\hfill \\ \\ ∥{\mathsf{\Xi }}_{\mathbf{I},n}(t)∥\le ∥{\mathbf{I}}_0(t)∥{\left(\frac{\mathcal{H}(\delta )}{\delta }{\kappa }^{\delta }{\mathsf{\Lambda }}_4\right)}^n,\hfill \\ \\ ∥{\mathsf{\Xi }}_{\mathbf{Q},n}(t)∥\le ∥{\mathbf{Q}}_0(t)∥{\left(\frac{\mathcal{H}(\delta )}{\delta }{\kappa }^{\delta }{\mathsf{\Lambda }}_5\right)}^n,\hfill \\ \\ ∥{\mathsf{\Xi }}_{\mathbf{R},n}(t)∥\le ∥{\mathbf{R}}_0(t)∥{\left(\frac{\mathcal{H}(\delta )}{\delta }{\kappa }^{\delta }{\mathsf{\Lambda }}_6\right)}^n.\hfill \end{array}$$

(26)

As a result, it is evident that the sequences fulfill and exist

$$∥{\mathsf{\Xi }}_{\mathbf{S},n}(t)∥\to 0,∥{\mathsf{\Xi }}_{\mathbf{V},n}(t)∥\to 0,∥{\mathsf{\Xi }}_{\mathbf{E},n}(t)∥\to 0,∥{\mathsf{\Xi }}_{\mathbf{I},n}(t)∥\to 0,∥{\mathsf{\Xi }}_{\mathbf{Q},n}(t)∥\to 0,∥{\mathsf{\Xi }}_{\mathbf{R},n}(t)∥\to 0\mathrm{as}n\to \infty .$$

Additionally, using Equation (26) and the triangle inequality, for any s, we have

$$\begin{array}{c}∥{\mathbf{S}}_{n+s}(t)-{\mathbf{S}}_n(t)∥\le \sum _{j=n+1}^{n+s}{X}_1^j=\frac{X_1^{n+1}-X_1^{n+s+1}}{1-X_1},\hfill \\ \\ ∥{\mathbf{V}}_{n+s}(t)-{\mathbf{V}}_n(t)∥\le \sum _{j=n+1}^{n+s}{X}_2^j=\frac{X_2^{n+1}-X_2^{n+s+1}}{1-X_2},\hfill \\ \\ ∥{\mathbf{E}}_{n+s}(t)-{\mathbf{E}}_n(t)∥\le \sum _{j=n+1}^{n+s}{X}_2^j=\frac{X_3^{n+1}-X_3^{n+s+1}}{1-X_3},\hfill \\ \\ ∥{\mathbf{I}}_{n+s}(t)-{\mathbf{I}}_n(t)∥\le \sum _{j=n+1}^{n+s}{X}_4^j=\frac{X_4^{n+1}-X_4^{n+s+1}}{1-X_4},\hfill \\ \\ ∥{\mathbf{Q}}_{n+s}(t)-{\mathbf{Q}}_n(t)∥\le \sum _{j=n+1}^{n+s}{X}_5^j=\frac{X_5^{n+1}-X_5^{n+s+1}}{1-X_5},\hfill \\ \\ ∥{\mathbf{R}}_{n+s}(t)-{\mathbf{R}}_n(t)∥\le \sum _{m=n+1}^{n+s}{X}_6^j=\frac{X_6^{n+1}-X_6^{n+s+1}}{1-X_6},\hfill \end{array}$$

(27)

with $X_m=\frac{\mathcal{H}(\delta )}{\delta }{\kappa }^{\delta }{\mathsf{\Lambda }}_m<1$ by hypothesis. Therefore, with ${\mathbf{S}}_n,{\mathbf{V}}_n,{\mathbf{E}}_n,{\mathbf{I}}_n,{\mathbf{Q}}_n$, it is possible to think of ${\mathbf{R}(t)}_n$ as a Cauchy sequence in the $\mathcal{A}(Y)$ Banach space. This has demonstrated that they are uniformly convergent [49]. □

5. Parameter Estimation

Parameter estimation in epidemic models involves identifying the model parameters that best fit the observed data. One widely used approach for parameter estimation is least-squared curve-fitting, which entails minimizing the difference between the observed data and the model predictions by adjusting the model parameters. This method assumes that the residuals of the model are normally distributed and have a constant variance. Least-squared curve-fitting tools estimate the parameters of the epidemic model by minimizing the sum of squared residuals. This method is computationally efficient and provides reliable estimates when the assumptions of normality and constant variance are met. However, it is critical to verify that the model assumptions are satisfied and that the goodness-of-fit is assessed using appropriate statistical tests. In this study, we employed least-squared curve-fitting methods to analyze the COVID-19 cases reported in Pakistan between 23 June 2022 and 23 August 2022 [50]. The estimated parameters of the system were based on Pakistan’s confirmed cases and fatalities. The ordinary least square solution (OLS) was used to minimize the error terms in the daily reports, and the goodness-of-fit was evaluated using the relative error.

$$min\left(\frac{{\sum }_{\iota =1}^n{\left({\mathbf{I}}_{\iota }-{\widehat{\mathbf{I}}}_l\right)}^2}{{\sum }_{\iota =1}^n{\mathbf{I}}_{\iota }^2}\right).$$

(28)

The simulated cumulative number of infected individuals, denoted by ${\mathbf{I}}_{\iota }$, is computed by adding the total number of individuals transitioning from the infected compartment to the recovered compartment each day. Meanwhile, the reported total number of infected individuals is represented by ${\widehat{\mathbf{I}}}_{\iota }$.

Table 1. The table contains descriptions and estimated values for the parameters.

Symbol of Parameters	Values of Parameters	References
B	$50.4057$	[51]
$\theta$	$0.0019$	Estimated
$\tau$	$0.0205$	Estimated
$ϵ$	$0.1571$	Estimated
${\epsilon }_1$	$0.01772$	Estimated
${\epsilon }_2$	$0.0805$	Estimated
${\epsilon }_3$	$0.0876$	Estimated
$\rho$	$0.0205$	[51]
$\eta$	$0.1460$	Estimated
$\alpha$	$0.3506$	Estimated
$\mathcal{K}$	$0.1530$	Estimated
$\sigma$	$0.0059$	Estimated
$\pi$	$0.0029$	Estimated
$\psi$	$0.0532$	Estimated
$\lambda$	$0.0105$	Estimated

presents the estimated parameter values, while a comparison between the model predictions and the reported cases is illustrated in Figure 1.

6. Sensitivity Analysis

A sensitivity analysis identifies the parameters that are most effective in curbing COVID-19 spread. Even though forward sensitivity analysis becomes tedious for complex biological models, it is an essential component of phenomena modeling. Ecologists and epidemiologists have taken an active interest in $R_0$ sensitivity analysis.

Definition 3.

The normalized forward sensitivity index of the $R_0$ that depends differentiability on a parameter ϰ is defined as

$${\mathsf{{\rm Y}}}_{\varkappa }=\frac{\varkappa }{R_0}\frac{\partial R_0}{\partial \varkappa }.$$

A vital tool for assessing how uncertain input parameters affect a system’s output is sensitivity analysis. Sensitivity indices can be computed using a variety of techniques, including direct differentiation, Latin hypercube sampling, and system linearization. In this study, the analytical expressions for the indices were provided using the direct differentiation method, allowing for a more thorough understanding of the behavior of the system. This study investigated the effects of various factors and gained a critical understanding of the comparative variability of the fundamental reproduction number, $R_0$, and other parameters by applying this method to a COVID-19 model. These sensitivity indices provide policymakers and law enforcement agencies with useful data to create strategies that effectively combat COVID-19, and they are therefore essential in battling the pandemic. The graphical results are displayed in Figure 2 and Figure 3.

7. Numerical Scheme and Graphical Results

The Haar wavelet collocation method is a powerful numerical technique used for solving various types of differential equations. One of the main advantages of this method is its high accuracy, as it is capable of producing highly precise solutions for a wide range of problems. Additionally, the Haar wavelet collocation method is highly efficient, having a relatively low computational cost compared with other numerical methods. This is due to the compact support of the Haar wavelet, which allows for the use of fewer basis functions to represent the solution, resulting in faster calculations. Overall, the Haar wavelet collocation method is a reliable and efficient tool for numerical analysis, making it a popular choice among researchers and practitioners in various fields of science and engineering.

Consider the square integrable function space $L_2[0,1)$, where $\dot{\mathbf{S}}(t),\dot{\mathbf{V}}(t),\dot{\mathbf{E}}(t),\dot{\mathbf{I}}(t),\dot{\mathbf{Q}}(t)$, and $\dot{\mathbf{R}}(t)$ can be expressed as a Haar series, given by:

$$\begin{array}{c}\hfill \dot{\mathbf{S}}(t)=\sum _{m=1}^N{\alpha }_m{\tilde{\psi }}_m(t)\\ \hfill \dot{\mathbf{V}(t)}=\sum _{m=1}^N{\gamma }_m{\tilde{\psi }}_m(t)\\ \hfill \dot{\mathbf{E}(t)}=\sum _{m=1}^N{\theta }_m{\tilde{\psi }}_m(t)\\ \hfill \dot{\mathbf{I}(t)}=\sum _{m=1}^N{\lambda }_m{\tilde{\psi }}_m(t),\\ \hfill \dot{\mathbf{Q}(t)}=\sum _{m=1}^N{\varpi }_m{\tilde{\psi }}_m(t)\\ \hfill \dot{\mathbf{R}(t)}=\sum _{m=1}^N{\sigma }_m{\tilde{\psi }}_m(t).\end{array}$$

Here, $\mathbf{S}(0)={\mathbf{S}}_0$ represents the initial susceptible compartment, $\mathbf{V}(0)={\mathbf{V}}_0$ represents the initial vaccinated compartment, $\mathbf{E}(0)={\mathbf{E}}_0$ represents the initial exposed compartment, $\mathbf{I}(0)={\mathbf{I}}_0$ represents the initial infected compartment, $\mathbf{Q}(0)={\mathbf{Q}}_0$ represents the initial infected compartment, and $\mathbf{R}(0)={\mathbf{R}}_0$ represents the initial recovered compartment. Integration of these equations leads to the following relation:

$$\begin{array}{cc}\hfill & \mathbf{S}(t)={\mathbf{S}}_0+\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t),\mathbf{V}(t)={\mathbf{V}}_0+\sum _{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t),\mathbf{E}(t)={\mathbf{E}}_0+\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t),\hfill \\ \hfill & \mathbf{I}(t)={\mathbf{I}}_0+\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t),\mathbf{Q}(t)={\mathbf{Q}}_0+\sum _{m=1}^N{\varpi }_m{\vartheta }_{m,1}(t),\mathbf{R}(t)={\mathbf{R}}_0+\sum _{m=1}^N{\sigma }_m{\vartheta }_{m,1}(t).\hfill \end{array}$$

(29)

By using the Caputo derivative, we have

$$\left\{\begin{array}{c}\frac{1}{\mathsf{\Gamma }(n-\delta )}{\int }_0^t{\mathbf{S}}^{(n)}(\varpi ){(t-\varpi )}^{n-\delta -1}d\varpi =B-\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))-({\epsilon }_1+\rho +\eta )\mathbf{S}(t)+\mathcal{K}\mathbf{V}(t),\hfill \\ \frac{1}{\mathsf{\Gamma }(n-\delta )}{\int }_0^t{\mathbf{V}}^{(n)}(\varpi ){(t-\varpi )}^{n-\delta -1}d\varpi =\pi \mathbf{S}(t)-\psi \mathbf{V}(t)-(\mathcal{K}+\rho )\mathbf{V}(t),\hfill \\ \frac{1}{\mathsf{\Gamma }(n-\delta )}{\int }_0^t{\mathbf{E}}^{(n)}(\varpi ){(t-\varpi )}^{n-\delta -1}d\varpi =\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))+\psi \mathbf{V}(t)-({\epsilon }_2+\rho +⫌)\mathbf{E}(t),\hfill \\ \frac{1}{\mathsf{\Gamma }(n-\delta )}{\int }_0^t{\mathbf{I}}^{(n)}(\varpi ){(t-\varpi )}^{n-\delta -1}d\varpi =\phi \mathbf{E}(t)-(\lambda +ϵ+\rho +{\epsilon }_3)\mathbf{I}(t),\hfill \\ \frac{1}{\mathsf{\Gamma }(n-\delta )}{\int }_0^t{\mathbf{Q}}^{(n)}(\varpi ){(t-\varpi )}^{n-\delta -1}d\varpi ={\epsilon }_1\mathbf{S}(t)+{\epsilon }_2\mathbf{E}(t)+{\epsilon }_3\mathbf{I}(t)+-(\rho +\sigma )\mathbf{Q}(t),\hfill \\ \frac{1}{\mathsf{\Gamma }(n-\delta )}{\int }_0^t{\mathbf{R}}^{(n)}(\varpi ){(t-\varpi )}^{n-\delta -1}d\varpi =\eta \mathbf{S}(t)+\sigma \mathbf{Q}(t)+\lambda \mathbf{I}(t)-\rho \mathbf{R}(t).\hfill \end{array}\right.$$

(30)

As we have assumed that $0<\delta <1$, therefore $n=1$, and we have

$$\left\{\begin{array}{c}\frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t\dot{\mathbf{S}}(t)(\varpi ){(t-\varpi )}^{-\delta }d\varpi =B-\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))-({\epsilon }_1+\rho +\eta )\mathbf{S}(t)+\mathcal{K}\mathbf{V}(t),\hfill \\ \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t\dot{\mathbf{V}(t)}(\varpi ){(t-\varpi )}^{-\delta }d\varpi =\pi \mathbf{S}(t)-\psi \mathbf{V}(t)-(\rho +\mathcal{K})\mathbf{V}(t),\hfill \\ \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t\dot{\mathbf{E}(t)}(\varpi ){(t-\varpi )}^{-\delta }d\varpi =\theta \mathbf{S}(t)\mathbf{I}(t)(1+\tau \mathbf{I}(t))+\psi \mathbf{V}(t)-({\epsilon }_2+\rho +⫌)\mathbf{E}(t),\hfill \\ \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t\dot{\mathbf{I}(t)}(\varpi ){(t-\varpi )}^{-\delta }d\varpi =\phi \mathbf{E}(t)-(\lambda +ϵ+\rho +{\epsilon }_3)\mathbf{I}(t),\hfill \\ \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t\dot{\mathbf{Q}(t)}(\varpi ){(t-\varpi )}^{-\delta }d\varpi ={\epsilon }_1\mathbf{S}(t)+{\epsilon }_2\mathbf{E}(t)+{\epsilon }_3\mathbf{I}(t)-(\rho +\sigma )\mathbf{Q}(t),\hfill \\ \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t\dot{\mathbf{R}(t)}(\varpi ){(t-\varpi )}^{-\delta }d\varpi =\eta \mathbf{S}(t)+\sigma \mathbf{Q}(t)+\lambda \mathbf{I}(t)-\rho \mathbf{R}(t).\hfill \end{array}\right.$$

(31)

Haar approximations are used, and we have

$$\begin{array}{cc}\hfill \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t& \sum _{m=1}^N{\alpha }_m{\tilde{\psi }}_m(\varpi ){(t-\varpi )}^{-\delta }d\varpi =B-\theta \left({\mathbf{I}}_0+\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)\left({\mathbf{S}}_0+\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)\hfill \\ & \left(1+\tau \left({\mathbf{I}}_0+\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)\right)-({\epsilon }_1+\rho +\eta )\left({\mathbf{S}}_0+\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)+\mathcal{K}\left({\mathbf{S}}_0+\sum _{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t& \sum _{m=1}^N{\gamma }_m{\tilde{\psi }}_m(\varpi ){(t-\varpi )}^{-\delta }d\varpi =\pi \left({\mathbf{S}}_0+\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)\hfill \\ & -\psi \left({\mathbf{V}}_0+\sum _{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t)\right)-(\rho +\mathcal{K})\left({\mathbf{V}}_0+\sum _{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t& \sum _{m=1}^N{\theta }_m{\tilde{\psi }}_m(\varpi ){(t-\varpi )}^{-\delta }d\varpi =\theta \left({\mathbf{I}}_0+\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)\left({\mathbf{S}}_0+\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)+\psi \left({\mathbf{V}}_0+\sum _{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t)\right)\hfill \\ & \left(1+\tau \left({\mathbf{I}}_0+\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)\right)-({\epsilon }_2+\rho +\phi )\left({\mathbf{E}}_0+\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t& \sum _{m=1}^N{\lambda }_m{\tilde{\psi }}_m(\varpi ){(t-\varpi )}^{-\delta }d\varpi =\hfill \\ & \phi \left({\mathbf{E}}_0+\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)-(\lambda +ϵ+\rho +d_1)\left({\mathbf{I}}_0+\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t& \sum _{m=1}^N{\varpi }_m{\tilde{\psi }}_m(\varpi ){(t-⫋)}^{-\delta }d\varpi ={\epsilon }_3\left({\mathbf{I}(t)}_0+\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)+{\epsilon }_2\left({\mathbf{E}}_0+\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ & +{\epsilon }_1\left({\mathbf{S}}_0+\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)-(\rho +\sigma )\left({\mathbf{Q}}_0+\sum _{m=1}^N{\varpi }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t& \sum _{m=1}^N{\sigma }_m{\tilde{\psi }}_m(\varpi ){(t-\varpi )}^{-\delta }d\varpi =\eta \left({\mathbf{S}}_0+\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)\hfill \\ & +\sigma \left({\mathbf{Q}}_0+\sum _{m=1}^N{\sigma }_m{\vartheta }_{m,1}(t)\right)+\lambda \left({\mathbf{I}}_0+\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)+\rho \left({\mathbf{R}}_0+\sum _{m=1}^N{\zeta }_m{\vartheta }_{m,1}(t)\right).\hfill \end{array}$$

(32)

Upon simplification, we have

$$\left\{\begin{array}{c}\frac{1}{\mathsf{\Gamma }(1-\delta )}{\sum }_{m=1}^N{\alpha }_m{\tilde{\psi }}_m(\varpi ){(t-\varpi )}^{-\delta }d\varpi -B+\theta \left(1+\tau {\mathbf{I}}_0\right)\times \hfill \\ \\ \left({\mathbf{I}}_0{\mathbf{S}}_0+{\mathbf{I}}_0{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)+{\mathbf{S}}_0{\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \\ +\theta \left[{\mathbf{I}}_0{\mathbf{S}}_0\tau {\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\mathbf{I}}_0\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\mathbf{S}}_0\tau {\left({\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)}^2\right.\hfill \\ \\ \left.+\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\left({\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)}^2\right]+({\epsilon }_1+\rho +\eta ){\mathbf{S}}_0+({\epsilon }_1+\rho +\eta ){\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\hfill \end{array}\right\}=0,$$

(33)

$$\begin{array}{cc}\hfill \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t\sum _{m=1}^N{\gamma }_m{\tilde{\psi }}_m(\varpi )& {(t-\varpi )}^{-\delta }d\varpi -\pi \left({\mathbf{S}}_0+\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill & +\psi \left({\mathbf{V}}_0+\sum _{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t)\right)+(\rho +\mathcal{K})\left({\mathbf{V}}_0+\sum _{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t)\right)=0,\hfill \end{array}$$

(34)

$$\left\{\begin{array}{c}\frac{1}{\mathsf{\Gamma }(1-\delta )}{\sum }_{m=1}^N{\theta }_m{\tilde{\psi }}_m(\varpi ){(t-\varpi )}^{-\delta }d\varpi +\theta \left(1+\tau {\mathbf{I}}_0\right)\times \hfill \\ \\ \left({\mathbf{I}}_0{\mathbf{S}}_0+{\mathbf{I}}_0{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)+{\mathbf{S}}_0{\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \\ +\theta \left[{\mathbf{I}}_0{\mathbf{S}}_0\tau {\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\mathbf{I}}_0\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\mathbf{S}}_0\tau {\left({\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)}^2\right.\hfill \\ \\ \left.+\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\left({\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)}^2\right]+({\epsilon }_2+\rho +\phi ){\mathbf{E}}_0+({\epsilon }_2+\rho +\phi ){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\hfill \end{array}\right\}=0,$$

(35)

$$\begin{array}{cc}\hfill \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t& \sum _{m=1}^N{\lambda }_m{\tilde{\psi }}_m(\varpi ){(t-\varpi )}^{-\delta }d\varpi -\phi {\mathbf{E}}_0-\phi \left(\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill & +(\lambda +ϵ+\rho +{\epsilon }_3){\mathbf{I}}_0+(\lambda +ϵ+\rho +{\epsilon }_3)\left(\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)=0,\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t& \sum _{m=1}^N{\varpi }_m{\tilde{\psi }}_m(\varpi ){(t-⫋)}^{-\delta }d\varpi -{\epsilon }_3{\mathbf{I}(t)}_0-{\epsilon }_3\left(\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)+{\epsilon }_2{\mathbf{E}}_0+{\epsilon }_2\left(\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill & +{\epsilon }_1{\mathbf{S}}_0+{\epsilon }_1\left(\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)+(\rho +\sigma ){\mathbf{Q}}_0+(\rho +\sigma )\left(\sum _{m=1}^N{\varpi }_m{\vartheta }_{m,1}(t)\right)=0,\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill \frac{1}{\mathsf{\Gamma }(1-\delta )}{\int }_0^t& \sum _{m=1}^N{\zeta }_m{\tilde{\psi }}_m(\varpi ){(t-\varpi )}^{-\delta }d\varpi -\eta {\mathbf{S}}_0-{\mathbf{S}}_0\left(\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)-\sigma {\mathbf{Q}}_0\hfill \\ \hfill & -\sigma \left(\sum _{m=1}^N{\sigma }_m{\vartheta }_{m,1}(t)\right)-\lambda {\mathbf{I}}_0-\lambda \left(\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)+\rho {\mathbf{R}}_0+\rho \left(\sum _{m=1}^N{\zeta }_m{\vartheta }_{m,1}(t)\right)=0.\hfill \end{array}$$

(38)

Using the method of Haar integration [52], the integral in the aforementioned system is approximately calculated as

$$\begin{array}{c}\hfill {\int }_{\varkappa }^{\kappa }f(t)dt\approx \frac{\kappa -\varkappa }{N}\sum _{p=1}^{N}f\left(t_p\right)=\sum _{p=1}^{N}f\left(\varkappa +\frac{(\kappa -\varkappa )(p-0.5)}{N}\right)\end{array}$$

(39)

$$\left\{\begin{array}{c}\frac{t}{N\mathsf{\Gamma }(1-\delta )}{\sum }_{s=1}^N{\sum }_{m=1}^N{\alpha }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t-{\varpi }_s\right)}^{-\delta }-B+\theta \left(1+\tau {\mathbf{I}}_0\right)\times \hfill \\ \\ \left({\mathbf{I}}_0{\mathbf{S}}_0+{\mathbf{I}}_0{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)+{\mathbf{S}}_0{\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \\ +\theta \left[{\mathbf{I}}_0{\mathbf{S}}_0\tau {\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\mathbf{I}}_0\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\mathbf{S}}_0\tau {\left({\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)}^2\right.\hfill \\ \\ \left.+\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\left({\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)}^2\right]+({\epsilon }_1+\rho +\eta ){\mathbf{S}}_0+({\epsilon }_1+\rho +\eta ){\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\hfill \end{array}\right\}=0$$

(40)

$$\begin{array}{cc}\hfill \frac{1}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N{\gamma }_m{\tilde{\psi }}_m(\varpi )& {(t-\varpi )}^{-\delta }d\varpi -\pi \left({\mathbf{S}}_0+\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill & +\psi \left({\mathbf{V}}_0+\sum _{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t)\right)+(\rho +\mathcal{K})\left({\mathbf{V}}_0+\sum _{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t)\right)=0,\hfill \end{array}$$

(41)

$$\left\{\begin{array}{c}\frac{t}{N\mathsf{\Gamma }(1-\delta )}{\sum }_{s=1}^N{\sum }_{m=1}^N{\theta }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t-{\varpi }_s\right)}^{-\delta }-\theta \left(1+\tau {\mathbf{I}}_0\right)\times \hfill \\ \\ \left({\mathbf{I}}_0{\mathbf{S}}_0+{\mathbf{I}}_0{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)+{\mathbf{S}}_0{\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \\ +\theta \left[{\mathbf{I}}_0{\mathbf{S}}_0\tau {\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\mathbf{I}}_0\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\mathbf{S}}_0\tau {\left({\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)}^2\right.\hfill \\ \\ \left.+\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\left({\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)}^2\right]+({\epsilon }_2+\rho +\phi ){\mathbf{E}}_0+({\epsilon }_2+\rho +\phi ){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\hfill \end{array}\right\}=0,$$

(42)

$$\begin{array}{cc}\hfill \frac{t}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N& {\lambda }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t-{\varpi }_s\right)}^{-\delta }-\phi {\mathbf{E}}_0-\phi \left(\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill & +(\lambda +ϵ+\rho +{\epsilon }_3){\mathbf{I}}_0+(\lambda +ϵ+\rho +{\epsilon }_3)\left(\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)=0,\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill \frac{t}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N& {\varpi }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t-{\varpi }_s\right)}^{-\delta }-{\epsilon }_3{\mathbf{I}}_0-{\epsilon }_3\left(\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)+{\epsilon }_2{\mathbf{E}(t)}_0+{\epsilon }_2\left(\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill & +{\epsilon }_1{\mathbf{S}}_0+{\epsilon }_1\left(\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)+(\rho +\sigma ){\mathbf{Q}}_0+(\rho +\sigma )\left(\sum _{m=1}^N{\varpi }_m{\vartheta }_{m,1}(t)\right)=0,\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill \frac{t}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N& {\zeta }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t-{\varpi }_s\right)}^{-\delta }-\eta {\mathbf{S}}_0-{\mathbf{S}}_0\left(\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)-\sigma {\mathbf{Q}}_0\hfill \\ \hfill & -\sigma \left(\sum _{m=1}^N{\sigma }_m{\vartheta }_{m,1}(t)\right)-\lambda {\mathbf{I}}_0-\lambda \left(\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)+\rho {\mathbf{R}}_0+\rho \left(\sum _{m=1}^N{\zeta }_m{\vartheta }_{m,1}(t)\right)=0.\hfill \end{array}$$

(45)

Let

$${\mathsf{\Phi }}_{1,j}=\left\{\begin{array}{c}\frac{t}{N\mathsf{\Gamma }(1-\delta )}{\sum }_{s=1}^N{\sum }_{m=1}^N{\alpha }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t-{\varpi }_s\right)}^{-\delta }-B+\theta \left(1+\tau {\mathbf{I}}_0\right)\times \hfill \\ \\ \left({\mathbf{I}}_0{\mathbf{S}}_0+{\mathbf{I}}_0{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)+{\mathbf{S}}_0{\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \\ +\theta \left[{\mathbf{I}}_0{\mathbf{S}}_0\tau {\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\mathbf{I}}_0\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\mathbf{S}}_0\tau {\left({\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)}^2\right.\hfill \\ \\ \left.+\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\left({\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)}^2\right]+({\epsilon }_1+\rho +\eta ){\mathbf{S}}_0+({\epsilon }_1+\rho +\eta ){\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t).\hfill \end{array}\right.$$

(46)

Let

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{2,j}=\frac{1}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N{\gamma }_m{\tilde{\psi }}_m(\varpi )& {(t-\varpi )}^{-\delta }d\varpi -\pi \left({\mathbf{S}}_0+\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill & +\psi \left({\mathbf{V}}_0+\sum _{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t)\right)+(\rho +\mathcal{K})\left({\mathbf{V}}_0+\sum _{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t)\right).\hfill \end{array}$$

(47)

Let

$${\mathsf{\Phi }}_{3,j}=\left\{\begin{array}{c}\frac{t}{N\mathsf{\Gamma }(1-\delta )}{\sum }_{s=1}^N{\sum }_{m=1}^N{\theta }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t-{\varpi }_s\right)}^{-\delta }-\theta \left(1+\tau {\mathbf{I}}_0\right)\times \hfill \\ \\ \left({\mathbf{I}}_0{\mathbf{S}}_0+{\mathbf{I}}_0{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)+{\mathbf{S}}_0{\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \\ +\theta \left[{\mathbf{I}}_0{\mathbf{S}}_0\tau {\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\mathbf{I}}_0\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)+{\mathbf{S}}_0\tau {\left({\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)}^2\right.\hfill \\ \\ \left.+\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t){\left({\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)}^2\right]+({\epsilon }_2+\rho +\phi ){\mathbf{E}}_0+({\epsilon }_2+\rho +\phi ){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t).\hfill \end{array}\right.$$

(48)

Let

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{4,j}=\frac{t}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N& {\lambda }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t-{\varpi }_s\right)}^{-\delta }-\phi {\mathbf{E}}_0-\phi \left(\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill & +(\lambda +ϵ+\rho +{\epsilon }_3){\mathbf{I}}_0+(\lambda +ϵ+\rho +{\epsilon }_3)\left(\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right).\hfill \end{array}$$

(49)

Let

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{5,j}=\frac{t}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N& {\varpi }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t-{\varpi }_s\right)}^{-\delta }-{\epsilon }_3{\mathbf{I}}_0-{\epsilon }_3\left(\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)+{\epsilon }_2{\mathbf{E}(t)}_0+{\epsilon }_2\left(\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t)\right)\hfill \\ \hfill & +{\epsilon }_1{\mathbf{S}}_0+{\epsilon }_1\left(\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)+(\rho +\sigma ){\mathbf{Q}}_0+(\rho +\sigma )\left(\sum _{m=1}^N{\varpi }_m{\vartheta }_{m,1}(t)\right).\hfill \end{array}$$

(50)

Let

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{6,j}=\frac{t}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N& {\zeta }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t-{\varpi }_s\right)}^{-\delta }-\eta {\mathbf{S}}_0-{\mathbf{S}}_0\left(\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t)\right)-\sigma {\mathbf{Q}}_0\hfill \\ \hfill & -\sigma \left(\sum _{m=1}^N{\sigma }_m{\vartheta }_{m,1}(t)\right)-\lambda {\mathbf{I}}_0-\lambda \left(\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t)\right)+\rho {\mathbf{R}}_0+\rho \left(\sum _{m=1}^N{\zeta }_m{\vartheta }_{m,1}(t)\right).\hfill \end{array}$$

(51)

The nodal points are placed to create the system of nonlinear algebraic equations shown below:

$${\mathsf{\Phi }}_{1,j}=\left\{\begin{array}{c}\frac{t_m}{N\mathsf{\Gamma }(1-\delta )}{\sum }_{s=1}^N{\sum }_{m=1}^N{\alpha }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t_m-{\varpi }_s\right)}^{-\delta }-B+\theta \left(1+\tau {\mathbf{I}}_0\right)\times \hfill \\ \\ \left({\mathbf{I}}_0{\mathbf{S}}_0+{\mathbf{I}}_0{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m)+{\mathbf{S}}_0{\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)+{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m){\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right)\hfill \\ \\ +\theta \left[{\mathbf{I}}_0{\mathbf{S}}_0\tau {\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)+{\mathbf{I}}_0\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m){\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right.\hfill \\ \\ \left.+{\mathbf{S}}_0\tau {\left({\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right)}^2+\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m){\left({\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right)}^2\right]\hfill \\ \\ +({\epsilon }_1+\rho +\eta ){\mathbf{S}}_0+({\epsilon }_1+\rho +\eta ){\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m)\hfill \end{array},\right.$$

(52)

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{2,j}=\frac{1}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N{\gamma }_m{\tilde{\psi }}_m(\varpi )& {(t_m-\varpi )}^{-\delta }d\varpi -\pi \left({\mathbf{S}}_0+\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m)\right)\hfill \\ \hfill & +\psi \left({\mathbf{V}}_0+\sum _{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t_m)\right)+(\rho +\mathcal{K})\left({\mathbf{V}}_0+\sum _{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t_m)\right).\hfill \end{array}$$

(53)

$${\mathsf{\Phi }}_{3,j}=\left\{\begin{array}{c}\frac{t_m}{N\mathsf{\Gamma }(1-\delta )}{\sum }_{s=1}^N{\sum }_{m=1}^N{\theta }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t_m-{\varpi }_s\right)}^{-\delta }-B+\theta \left(1+\tau {\mathbf{I}}_0\right)\times \hfill \\ \\ \left({\mathbf{I}}_0{\mathbf{S}}_0+{\mathbf{I}}_0{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m)+{\mathbf{S}}_0{\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)+{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m){\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right)\hfill \\ \\ +\theta \left[{\mathbf{I}}_0{\mathbf{S}}_0\tau {\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)+{\mathbf{I}}_0\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m){\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right.\hfill \\ \\ \left.+{\mathbf{S}}_0\tau {\left({\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right)}^2+\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m){\left({\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right)}^2\right]\hfill \\ \\ +\psi \left({\mathbf{V}}_0+{\sum }_{m=1}^N{\gamma }_m{\vartheta }_{m,1}(t_m)\right)+({\epsilon }_2+\rho +\phi ){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t_m)\hfill \end{array},\right.$$

(54)

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{4,j}=\frac{t_m}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N& {\lambda }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t_m-{\varpi }_s\right)}^{-\delta }-\phi {\mathbf{E}}_0-\phi \left(\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t_m)\right)\hfill \\ \hfill & +(\lambda +ϵ+\rho +{\epsilon }_3){\mathbf{I}}_0+(\lambda +ϵ+\rho +{\epsilon }_3)\left(\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right).\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{5,j}=\frac{t_m}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N& {\varpi }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t_m-{\varpi }_s\right)}^{-\delta }-{\epsilon }_3{\mathbf{I}}_0-{\epsilon }_3\left(\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right)\hfill \\ \hfill & +{\epsilon }_2{\mathbf{E}}_0+{\epsilon }_2\left(\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t_m)\right)+{\epsilon }_1{\mathbf{S}}_0+{\epsilon }_1\left(\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m)\right)\hfill \\ \hfill & +(\rho +\sigma ){\mathbf{Q}}_0+(\rho +\sigma )\left(\sum _{m=1}^N{\varpi }_m{\vartheta }_{m,1}(t_m)\right).\hfill \end{array}$$

(56)

Let

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{6,j}=\frac{t_m}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N& {\zeta }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t_m-{\varpi }_s\right)}^{-\delta }-\eta {\mathbf{S}}_0-{\mathbf{S}}_0\left(\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m)\right)-\sigma {\mathbf{Q}}_0\hfill \\ \hfill & -\sigma \left(\sum _{m=1}^N{\sigma }_m{\vartheta }_{m,1}(t_m)\right)-\lambda {\mathbf{I}}_0-\lambda \left(\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right)+\rho {\mathbf{R}}_0\hfill \\ \hfill & +\rho \left(\sum _{m=1}^N{\zeta }_m{\vartheta }_{m,1}(t_m)\right).\hfill \end{array}$$

(57)

Using Broyden’s approach, this system is resolved. The Jacobian is given by

$$\mathbf{J}={\left[J_{jp}\right]}_{6N\times 6N}$$

(58)

where

$$\left\{\begin{array}{cc}\hfill \frac{\partial {\mathsf{\Phi }}_{1,j}}{\partial {\alpha }_m}& =\left\{\begin{array}{c}\frac{t_m}{N\mathsf{\Gamma }(1-\delta )}{\sum }_{s=1}^N{\sum }_{m=1}^N{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t_m-{\varpi }_s\right)}^{-\delta }+\theta \left(1+\tau {\mathbf{I}}_0\right)\times \hfill \\ \left(+{\mathbf{I}}_0{\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)+{\sum }_{m=1}^N{\vartheta }_{m,1}(t_m){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t_m)\right)\hfill \\ +\theta \left[+{\mathbf{I}}_0\tau {\sum }_{m=1}^N{\vartheta }_{m,1}(t_m){\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)+\tau {\sum }_{m=1}^N{\vartheta }_{m,1}(t_m){\left({\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right)}^2\right]\hfill \\ +({\epsilon }_1+\rho +\eta ){\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)\hfill \end{array}\right.\hfill \\ & \frac{\partial {\mathsf{\Phi }}_{1,j}}{\partial {\gamma }_m}=\sum _{m=1}^N{\vartheta }_{m,1}(t_m),\frac{\partial {\mathsf{\Phi }}_{1,j}}{\partial {\theta }_m}=0,\frac{\partial {\mathsf{\Phi }}_{1,j}}{\partial {\varpi }_m}=0,\frac{\partial {\mathsf{\Phi }}_{1,j}}{\partial {\sigma }_m}=0,\hfill \\ \hfill \frac{\partial {\mathsf{\Phi }}_{1,j}}{\partial {\lambda }_m}& =\left\{\begin{array}{c}\theta \left(1+\tau {\mathbf{I}}_0\right)\times \left({\mathbf{S}}_0{\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)+{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m){\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)\right)\hfill \\ +\theta \left[{\mathbf{I}}_0{\mathbf{S}}_0\tau {\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)+{\mathbf{I}}_0\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m){\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)\right.\hfill \\ \left.+2{\lambda }_m{\mathbf{S}}_0\tau {\left({\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)\right)}^2+\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m){\left({\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)\right)}^2\right]\hfill \end{array}\right.\hfill \end{array},\right.$$

(59)

$$\begin{array}{cc}\hfill & \frac{\partial {\mathsf{\Phi }}_{2,j}}{\partial {\alpha }_m}=-\pi \left(\sum _{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m)\right),\frac{\partial {\mathsf{\Phi }}_{2,j}}{\partial {\theta }_m}=0,\frac{\partial {\mathsf{\Phi }}_{2,j}}{\partial {\lambda }_m}=0,\frac{\partial {\mathsf{\Phi }}_{1,j}}{\partial {\varpi }_m}=0,\frac{\partial {\mathsf{\Phi }}_{1,j}}{\partial {\sigma }_m}=0,\hfill \\ \hfill & \frac{\partial {\mathsf{\Phi }}_{2,j}}{\partial {\gamma }_m}=\frac{1}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N{\gamma }_m{\tilde{\psi }}_m(\varpi ){(t_m-\varpi )}^{-\delta }d\varpi +\psi \left({\mathbf{V}}_0+\sum _{m=1}^N{\vartheta }_{m,1}(t_m)\right)+(\rho +\mathcal{K})\left({\mathbf{V}}_0+\sum _{m=1}^N{\vartheta }_{m,1}(t_m)\right).\hfill \end{array}$$

(60)

$$\left\{\begin{array}{cc}\hfill \frac{\partial {\mathsf{\Phi }}_{3,j}}{\partial {\alpha }_m}& =\left\{\begin{array}{c}\theta \left(1+\tau {\mathbf{I}}_0\right)\times \left(+{\mathbf{I}}_0{\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)+{\sum }_{m=1}^N{\vartheta }_{m,1}(t_m){\sum }_{m=1}^N{\theta }_m{\vartheta }_{m,1}(t_m)\right)\hfill \\ +\theta \left[+{\mathbf{I}}_0\tau {\sum }_{m=1}^N{\vartheta }_{m,1}(t_m){\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right.\left.+\tau {\sum }_{m=1}^N{\vartheta }_{m,1}(t_m){\left({\sum }_{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right)}^2\right]\hfill \end{array}\right.\hfill \\ & \frac{\partial {\mathsf{\Phi }}_{3,j}}{\partial {\gamma }_m}=\sum _{m=1}^N{\vartheta }_{m,1}(t_m),\frac{\partial {\mathsf{\Phi }}_{3,j}}{\partial {\theta }_m}=\frac{t_m}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t_m-{\varpi }_s\right)}^{-\delta }+({\epsilon }_2+\rho +\phi )\sum _{m=1}^N{\vartheta }_{m,1}(t_m),\hfill \\ & \frac{\partial {\mathsf{\Phi }}_{3,j}}{\partial {\varpi }_m}=0,\frac{\partial {\mathsf{\Phi }}_{3,j}}{\partial {\sigma }_m}=0,\hfill \\ \hfill \frac{\partial {\mathsf{\Phi }}_{3,j}}{\partial {\lambda }_m}& =\left\{\begin{array}{c}\theta \left(1+\tau {\mathbf{I}}_0\right)\times \left({\mathbf{S}}_0{\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)+{\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m){\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)\right)\hfill \\ +\theta \left[{\mathbf{I}}_0{\mathbf{S}}_0\tau {\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)+{\mathbf{I}}_0\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m){\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)\right.\hfill \\ \left.+2{\lambda }_m{\mathbf{S}}_0\tau {\left({\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)\right)}^2+\tau {\sum }_{m=1}^N{\alpha }_m{\vartheta }_{m,1}(t_m){\left({\sum }_{m=1}^N{\vartheta }_{m,1}(t_m)\right)}^2\right]\hfill \end{array}\right.\hfill \end{array},\right.$$

(61)

$$\begin{array}{cc}\hfill & \frac{\partial {\mathsf{\Phi }}_{4,j}}{\partial {\alpha }_m}=0,\frac{\partial {\mathsf{\Phi }}_{4,j}}{\partial {\gamma }_m}=0,\frac{\partial {\mathsf{\Phi }}_{4,j}}{\partial {\varpi }_m}=0,\frac{\partial {\mathsf{\Phi }}_{4,j}}{\partial {\sigma }_m}=0,\frac{\partial {\mathsf{\Phi }}_{4,j}}{\partial {\theta }_m}=-\phi \left(\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t_m)\right),\hfill \\ \hfill & \frac{\partial {\mathsf{\Phi }}_{4,j}}{\partial {\lambda }_m}=\frac{t_m}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t_m-{\varpi }_s\right)}^{-\delta }(\lambda +ϵ+\rho +{\epsilon }_3)\left(\sum _{m=1}^N{\vartheta }_{m,1}(t_m)\right).\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill & \frac{\partial {\mathsf{\Phi }}_{5,j}}{\partial {\alpha }_m}={\epsilon }_1\left(\sum _{m=1}^N{\vartheta }_{m,1}(t_m)\right),\frac{\partial {\mathsf{\Phi }}_{5,j}}{\partial {\gamma }_m}=0,\frac{\partial {\mathsf{\Phi }}_{5,j}}{\partial {\varpi }_m}=\frac{t_m}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N{\varpi }_m{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t_m-{\varpi }_s\right)}^{-\delta },\hfill \\ \hfill & +(\rho +\sigma )\left(\sum _{m=1}^N{\varpi }_m{\vartheta }_{m,1}(t_m)\right)\frac{\partial {\mathsf{\Phi }}_{5,j}}{\partial {\sigma }_m}=0,\frac{\partial {\mathsf{\Phi }}_{5,j}}{\partial {\theta }_m}={\epsilon }_2\left(\sum _{m=1}^N{\theta }_m{\vartheta }_{m,1}(t_m)\right),\hfill \\ \hfill & \frac{\partial {\mathsf{\Phi }}_{5,j}}{\partial {\lambda }_m}=-{\epsilon }_3\left(\sum _{m=1}^N{\vartheta }_{m,1}(t_m)\right).\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill & \frac{\partial {\mathsf{\Phi }}_{6,j}}{\partial {\alpha }_m}={\mathbf{S}}_0\left(\sum _{m=1}^N{\vartheta }_{m,1}(t_m)\right),\frac{\partial {\mathsf{\Phi }}_{6,j}}{\partial {\gamma }_m}=0,\frac{\partial {\mathsf{\Phi }}_{6,j}}{\partial {\varpi }_m}=\frac{t_m}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t_m-{\varpi }_s\right)}^{-\delta },\hfill \\ \hfill & +(\rho +\sigma )\left(\sum _{m=1}^N{\vartheta }_{m,1}(t_m)\right)\frac{\partial {\mathsf{\Phi }}_{6,j}}{\partial {\theta }_m}=0,\frac{\partial {\mathsf{\Phi }}_{6,j}}{\partial {\lambda }_m}=-\lambda \left(\sum _{m=1}^N{\lambda }_m{\vartheta }_{m,1}(t_m)\right).\hfill \\ \hfill & \frac{\partial {\mathsf{\Phi }}_{6,j}}{\partial {\sigma }_m}=\frac{t_m}{N\mathsf{\Gamma }(1-\delta )}\sum _{s=1}^N\sum _{m=1}^N{\tilde{\psi }}_m\left({\varpi }_s\right){\left(t_m-{\varpi }_s\right)}^{-\delta }-\sigma \left(\sum _{m=1}^N{\sigma }_m{\vartheta }_{m,1}(t_m)\right)+\rho \left(\sum _{m=1}^N{\vartheta }_{m,1}(t_m)\right).\hfill \end{array}$$

(64)

This system’s solution yields the values of ${\alpha }_m$’s, ${\gamma }_m$’s, ${\lambda }_m$’s, ${\varpi }_m$’s, and ${\sigma }_m$’s unknown coefficients. By entering ${\alpha }_m$’s, ${\gamma }_m$’s, ${\lambda }_m$’s, ${\varpi }_m$’s, and ${\sigma }_m$’s unknown coefficients into Equation (29), it is possible to calculate the necessary solutions $\mathbf{S}(t),\mathbf{V}(t),\mathbf{E}(t),\mathbf{I}(t),\mathbf{Q}(t)$, and $\mathbf{R}(t)$ at nodal locations. The experimental rate of convergence, denoted by the formula r${}_{\varrho }(N)$ [53], can be calculated as follows:

$${\mathbf{r}}_{\varrho }(N)=\frac{1}{log2}log\left[\frac{\mathrm{Maximum}\mathrm{absolute}\mathrm{error}\mathrm{at}\frac{N}{2}}{\mathrm{Maximum}\mathrm{absolute}\mathrm{error}\mathrm{at}N}\right].$$

(65)

Graphical Results

This section presents graphical results for the fractional-order model (14). The fractional model was numerically solved based on the method outlined in [47,54] and using the information from Table 1. The resulting figures, namely Figure 4, Figure 5, Figure 6 and Figure 7, illustrate the dynamics of the susceptible $\mathbf{S}(t)$, vaccinated $\mathbf{V}(t)$, exposed $\mathbf{E}(t)$, infected $\mathbf{I}(t)$, quarantined $\mathbf{Q}(t)$, and recovered $\mathbf{R}(t)$ individuals. Figure 5a shows that susceptible individuals are characterized by fractional-order derivatives ranging between $0.55$ and $0.95$, and as time progresses, the number of susceptible individuals decreases due to exposure to the virus. This behavior is expected and observed in other epidemiological models. The population of vaccinated individuals is shown in Figure 5b, and it grows both steadily and quickly as the fractional-order derivative approaches its classical counterpart. Similarly, the population of exposed individuals is shown in Figure 5c, and it greatly decreases in the first 15 days while both steadily and quickly increasing as the fractional-order derivative approaches its classical counterpart. This increase is brought on by more susceptible people becoming infected during the first few weeks of the outbreak and joining the exposed class. An increased risk of transmission during the early stages of the epidemic may be indicated by the rise in exposed individuals. Figure 5c displays the number of infected people, which rises as the fractional order gets closer to one. The fractional order, which becomes more sensitive as it gets closer to one, is what is causing this increase. Within a few weeks of exposure, the majority of people in the quarantined and infectious stages of the infection leave the exposed class. Figure 5d illustrates the dynamics of the quarantined individuals, who exhibit a similar behavior as the exposed population. The population of exposed individuals increases as the fractional-order derivative approaches the integer order. Figure 5e shows how the fractional order affects the number of recovered individuals, which grows steadily as the fractional-order derivative approaches the classical value. This is due to the infected individuals recovering, which aids in disease containment. Raising the fractional order can cause the population of the recovered class to grow considerably more quickly. These results suggest that the exposed population is similarly affected by the fractional order regardless of the strain of infection.

8. Conclusions

The $SVEIQR$ COVID-19 epidemic model under fractional derivatives in the Caputo sense has been studied in this article. Some findings from fixed point theory are taken into consideration when determining the existence and uniqueness of the proposed model. This study simulates the spread of the COVID-19 Omicron variant in Pakistan using fractional-order models and the Haar wavelet collocation method. The study makes use of actual data from Pakistan, which improves the model’s precision and applicability. To determine the factors that most influence the spread of the virus, a sensitivity analysis was also carried out. This study sheds light on the dynamics of the pandemic and emphasizes the significance of taking fractional order into account when modeling infectious diseases. In terms of future directions, this study can be extended to incorporate additional elements such as vaccination rates, travel restrictions, and socio-economic conditions. To evaluate the efficacy of various intervention strategies and compare the effects of various COVID-19 variants, the model can also be applied to other regions or nations. Additionally, this study can be expanded to assess the effects of various compliance levels with public health measures, such as mask use and social distancing. To determine the most effective and precise method, this study can be expanded to compare the performance of various numerical solution techniques for fractional-order models.