#### Abstract

This article is devoted to investigate a mathematical model consisting on susceptible, exposed, infected, quarantined, vaccinated, and recovered compartments of COVID-19. The concerned model describes the transmission mechanism of the disease dynamics with therapeutic measures of vaccination of susceptible people along with the cure of the infected population. In the said study, we use the fractal-fractional order derivative to understand the dynamics of all compartments of the proposed model in more detail. Therefore, the first model is formulated. Then, two equilibrium points disease-free (DF) and endemic are computed. Furthermore, the basic threshold number is also derived. Some sufficient conditions for global asymptotical stability are also established. By using the next-generation matrix method, local stability analysis is developed. We also attempt the sensitivity analysis of the parameters of the proposed model. Finally, for the numerical simulations, the Adams–Bashforth method is used. Using some available data, the results are displayed graphically using various fractal-fractional orders to understand the mechanism of the dynamics. In addition, we compare our numerical simulation with real data in the case of reported infected cases.

#### 1. Introduction

Various viral infectious diseases have greatly affected human life in the past. Some famous viral infectious diseases were known as swine flu, influenza, SARs, MERs, etc. Recently, coronavirus disease 2019 (COVID-19) has greatly affected human life all over the world. The said infection has been reported for the first time in Wuhan, China, in December 2019. After three months, the disease took the form of a pandemic which was announced by the WHO in April 2020 (see [1]). Recent estimates issued by the WHO showed that the full death toll associated directly or indirectly with the COVID-19 pandemic between 1 January 2020 and 31 December 2021 was approximately 14.9 million (see details in [2]). Approximately, more than six hundred million people have got infections (see [3, 4]). Since the outbreak surrounded the whole globe, each nation in the world has taken its own measures. Some countries have imposed a strict curfew and precautionary measures including wearing face masks, keeping social distance, and avoiding large gatherings. Europe has been seriously affected along with USA, Brazil, Iran, and India. Initially it was reported that the disease has been transmitted from some animals to humans because the virus was first identified in a man who was working in a fish market. After that, various researchers proved that bats are not a unique source and the mentioned disease may be transmitted from other animals, person to person, and dogs (see [5, 6]).

Epidemiology is the most important branch of medical science and it has been very well developed in recent times. In the said area, various infectious diseases are investigated for treatment, controlling, curing, etc. Moreover, virology is a sub-branch, particularly dealing with virus and their transmitted diseases. Here, we remark that infectious diseases are investigated via various tools and methods. Biomathematics and bio-informatics engineering have attracted researchers recently more than in the past. It is worth mentioning that mathematical models are powerful tools to describe various real-world problems. Therefore, modeling infectious diseases is a hot area of research at the present time. In the said area, the dynamics of various infectious diseases are explained by using various differential or integral equations (some studies we refer to as [7–9]). Mathematical models help us to investigate various diseases for prediction and controlling procedures, to save society from great loss. Inspired by the aforesaid importance of mathematical models, researchers have formulated various real-world infectious diseases, we refer some as (see [10–15] and [16]). COVID-19 has been transmitted nearly to every part of the Earth and affected every nation of the world. The health conditions and economical situation of well-developed counties have been disturbed very well. In recent times, the world needs effective medication for the said diseases to save the lives of people. In this regard, great trials have been performed on vaccine preparation so far (for such detail, we refer [17]). Some authors have identified that COVID-19 vaccination in low and middle-income nations is extremely cost-effective and even price-saving (see [18]). The infection can be minimized all over the world by providing a vaccine to all nations of the world without any obstacles. Because by doing so, we can secure the next generation from the deadly virus. Moreover, some measures, like treatment, and vaccination, nonpharmaceutical precautionary measures, for instance, quarantine of confirmed cases, isolation, face masks, hand washing, social distancing, and avoiding gathering, should be imposed to reduce the transmission. Currently, some versions of the virus have been identified which indicate a new threat to human life. However, researchers have proved that the currently available vaccine is useful in medication. Moreover, it has been proved that vaccinated people if infected can easily be recovered compared to nonvaccinated people. However, this is not a permanent solution. Therefore, preparing proper vaccines and their availability without any bound in every society at a low price will be the best solution. Although, in some people, the vaccine has shown some adverse reactions. However, with the passage of time, things will become relaxed.

Most real-world problems are nonlinear in nature. In the same line, most epidemiological processes or phenomena are modeled in the form of nonlinear equations under traditional order derivatives or difference equations. As ordinary differential operators are local in nature and cannot produce the global dynamical behavior of the phenomenon, researchers are increasingly using fractional order derivatives in mathematical models (for instance, see [19–22], and [23]). A significant amount of work has been performed in this regard. Some well-known studies are referred here as [24–28], and [26]. In current times, in almost every discipline including fluid dynamics, bio-engineering, control theory, epidemiology, and rheology, the concept of fractional calculus is increasingly used to study various processes and phenomenons (see [29–32]). On the other hand, the said area has been exploited in engineering sides also recently (see [33–35]). Here it is worth mentioning that fractal geometry can be used as a powerful tool to investigate complex phenomena and every irregular picture in nature properly. Furthermore, fractal curves and surfaces are treated by using fractal dimensions to investigate their roughness. Because natural graphs are analyzed by using fractal interpolation (see some details in [36]). In previous times, due to the complex nature of fractional calculus, researchers have studied fractal surfaces through classical integer order calculus. Recently, many researchers have used the fractal-fractional order concept to investigate various real-world problems relating to epidemiology and physics, and chemical sciences. Here, we refer to some good work as [37].

Due to the ability to explain complex nature and preserve memory concepts, the idea of fractal calculus is widely used in the mathematical modeling of various diseases (details can be seen in [38]). For dealing with problems involving the fractal-fractional order derivatives, the traditional analytical and numerical techniques have been updated (see [39, 40], and [41]). Moreover, various numerical methods used for fractional derivatives have been applied to treat fractal-fractional problems. For instance, Adams–Bashforth procedure has been updated for dealing classical and fractal-fractional order problems (see [42–49]). Also, some authors have investigated other infectious diseases models like [50–54], and [55]. Authors [56] studied the psychological effects of staying home due to the COVID-19 pandemic. For instance, the authors applied a novel fractal-fractional order operator to the mathematical model of optimal control for malaria, where the derivative is defined in the Caputo sense. Two control variables have been introduced, as well as the necessary optimality condition in order to minimize the low-risk and high-risk infectious humans. Furthermore, the optimal control has been studied and results were simulated using a numerical scheme (see [57, 58]). In the same fashion, the author [59] proposed a compartmental mathematical model to acquire a better understanding of COVID-19’s future dynamics. The problem has been described as a highly nonlinear coupled system of classical order ODEs, which has then been generalized using the fractal-fractional derivative with the Mittag-Leffler kernel. For recent work on fractional order models, see [60], and [61]. Moreover, the boundedness and non-negativity of the model have been established as well. Here it should be kept in mind that authors [62] studied the backward bifurcation of a vaccination model with nonlinear incidence. Authors [63] studied a COVID-19 model. Authors [64] studied some stability results in control process by using fractional derivative. For more applications of fractional derivative, we refer [65]. Using mathematical models and the concept of fractional calculus, we refer some recent work like [66–69], and [70]. Also authors investigated a fractional order model for HIV infection in [71].

Motivated by the aforesaid discussion, we formulate the proposed model under the fractal-fractional order derivative. We investigate the proposed mathematical model with the vaccinated class by considering the available therapeutic measures, vaccination of susceptible, and curing of infected/hospitalized people. Our considered model involves some important epidemiological and biological features of the said disease like inhibitory effect, death rates due to infection and nature, birth rate, and different vaccination rates. A flow chart of our study is given in Figure 1.

Considering the compartmental diagram, the model is formulated as

Symbols involved in model (1) are described in Table 1. Also, the nomenclature and their dimensions are described in Table 2.

Using various tools of mathematical analysis, we establish global and local stability. The said analysis is established by the Lyapunov method and next-generation matrix method. Moreover, by using the Pontryagin maximum principle, we develop some results about the optimal control procedure. Furthermore, sensitivity analysis is also investigated about the parameters involved in our proposed model. Some results for numerical stability are derived by using the Ulam–Hyers concept. We also perform numerical simulation for our considered model by using a numerical method based on the Adams–Bashforth method to investigate multiphase behaviors under various fractal-fractional order derivatives. Here, we remark that multistep methods are increasingly used because these procedures are more efficient than the earlier methods like the Euler method. Moreover, the derivation of the Adams–Bashforth method can be performed in a number of ways. By using numerical interpolation and numerical integration, one can easily derive the aforesaid method. Therefore, we use the aforesaid scheme to simulate our model by using some numerical values of the parameters and initial data given in Table 1. A graphical presentation of some real data and discussion are provided.

#### 2. Preliminaries

In this section, we recall the fractional-fractional operators as given below.

*Definition 1. *(see [61]). Let the function be differentiable in the opened interval , then the fractional-fractional derivative with and are fractal and fractional orders, respectively, in the Caputo sense with power law is given asHere, is continuous over the interval .

*Definition 2. *(see [61]). The fractal-fractional integral of a function with fractal order and fractional order is defined aswhere is continuous over the interval .

*Definition 3. *Consider the fractal-fractional nonlinear ODE, such thatFrom [45], the resultant Adams–Bashforth scheme for (4) can be written as

#### 3. Equilibrium Points and Stability

This section is enriched with various stability results and equilibrium points.

##### 3.1. Equilibrium Points

The disease-free equilibrium point is given bywherewhile the positive disease-endemic equilibrium is computed in terms of one class such thatwhere

##### 3.2. Basic Reproduction Number

The nonlinear and linear terms of the infected classes in matrices and , respectively, are given as

Now, the Jacobian matrix of and is given by

Calculating the inverse of matrix and the next-generation matrix , such that

Thus, the nonzero and largest eigenvalue is the basic reproduction number , which is

##### 3.3. Stability Analysis

Lemma 1 (see [60]). *Let be a matrix, then the eigenvalues of the matrix are negative in real part if , and .*

Theorem 1. *The COVID-19 model at the disease-free equilibrium point , is locally asymptotically stable if , otherwise unstable.*

*Proof 1. *The Jacobian matrix of system (1) at disease-free equilibrium point is given byFollowing the characteristic equation of the Jacobian matrix (14), we haveThus, the eigenvalues of the characteristic (16) are given byAs a result, all eigenvalues of the Jacobian matrix (15) are negative for , such that . Hence, model (1) is locally asymptotically stable around a disease-free equilibrium point.

Theorem 2. *The COVID-19 model at the disease-endemic equilibrium point is locally asymptotically stable if , otherwise unstable.*

*Proof 2. *The Jacobian matrix of system (1) at disease-endemic equilibrium point is given byThus, eigenvalues of the characteristic (16) are given byHowever, the reduced matrix takes the formThe matrix (20) possesses negative eigenvalues if the trace is negative and determinant is positive, such thatThe determinant is positive if . Hence, model (1) is locally asymptotically stable around the disease-endemic equilibrium point with negative eigenvalues of and .

#### 4. Globally Asymptotical Stability

Theorem 3. *For nonautonomous fractional order system, let be an equilibrium point, such that*

Let be a domain which contains . Let be a continuously differentiable function such that and , for ,where continuous positive definite functions , and on and is Lyapunov candidate function, then, , is globally asymptotically stable.

Lemma 2. *For the fractal-fractional order, we use the lemma defined for fractional order ODE from [42, 60]. At any instant of time , assume a continuous differentiable function as*

Theorem 4. *Equilibrium point is globally asymptotically stable.*

*Proof 3. *We consider the quadratic Lyapunov function to derive the Lyapunov candidate function for fractional order differential equation.We define the Lyapunov candidate function asThe linearity property is given asUsing Lemma 2, we havewhich further implies that

##### 4.1. Global Stability at Disease-Free Equilibrium Point

Using disease-free equilibrium point, in the last equation of the proposed model, one has

The last inequality (29) is negative if , such that if . Hence, model (1) is globally asymptotically stable around disease-free equilibrium point, with .

##### 4.2. Global Stability at Disease-Endemic Equilibrium Point

For global stability around disease-endemic equilibrium point and , we havewhere

The last inequality (32) is negative if is positive such that, if . Hence, model (1) is globally asymptotically stable around disease-endemic equilibrium point, with .

#### 5. Sensitivity Analysis

Now, according to the (33) relation, we have

Here, in Figures 2–7, we have described graphically the dynamics of sensitivity analysis of various parameters, respectively.

Furthermore, we present Table 3 for the sensitivity index of each parameter associated with based on system (34).

Considering Table 3, the change in the value of each parameter in the basic reproduction number causes an increase or decrease in the value of the basic reproduction number . Hence, it directly affects the spread in the population while keeping the values of the remaining parameters constant. Furthermore, in Table 3, if the sensitivity index for parameter goes negatively, then there is an inverse effect on . The parameters show the negative effect must be minimized for the sake of the spread of infection in the population. Moreover, the sensitivity index of is positive at its peak.

#### 6. Ulam–Hyers Stability

Stability theory is an important branch of the qualitative theory of differential equations. As we know, the computation of exact solutions to some problems is quite challenging to obtain. Therefore, various numerical techniques were developed to find a solution. In this regard, we check the stability of the given problem. We can find various types of stability in literature, including Lyapunov, exponential, and asymptotic. But the most important type of stability, which is first introduced by Ulam in 1940 is called Ulam stability. He posed a problem about the stability of functional equations. The proper introduction was given by Hyers in 1941. Therefore, this stability was named Ulam–Hyers stability (see [72, 73]). The said stability results have been extended and generalized by researchers for difference and functional equations in different directions. From a numerical and optimization point of view, Ulam–Hyers stability is essential because it provides a bridge between the exact and numerical solutions. The said stability is most useful and also easy to establish for the approximated solution *o* model (1). Therefore, in this section, by the use of nonlinear functional analysis, some adequate conditions are constructed for the mentioned stability of the proposed model (1).

We express model (35) as

The proposed problem (35) can be reformulated in the following form:

With the help of (36), we can write the considered system (33) as

On applying the integral, we get the solution of (38) aswhere

For Ulam–Hyers stability, let , which is a small perturbation independent of solution. In addition,

, for ; .

Lemma 3. *The solution of the perturbed problem,satisfies the following relation*

*Proof 4. *According to (37), the solution of (39) is given byUsing in (43), one can easily get the relation (41).

Theorem 5. *Inview of the inequality (42), and hypothesis given in (44), if the condition holds, then the solution of problem (33) is Ulam–Hyers stable, where*

*Proof 5. *Suppose is a unique solution and be any other solution of (37) in Banach space ; then,From (45), we can writewhere . Thus, from the result (43), we conclude that (33) is Ulam–Hyers stable.

#### 7. Numerical Scheme

Nonlinear problems are more difficult than linear ones to solve analytically. Therefore, researchers have implemented efficient and accurate numerical schemes for the treatment of nonlinear problems to explore the dynamics of real-world problems more precisely. To tackle this nonlinear problem, from definition (3), we use the Adams–Bashforth scheme for the justification of our work.

*Remark 1. *Here it should be kept in mind that the Adams–Bashforth method preserves the basic type of numerical stability associated with the usual one-step numerical methods including Euler, backward Euler, and trapezoidal. Furthermore, it is a -step explicit method, and whose truncation error is of size . For detailed convergence and stability of the Adams–Bashforth method, we refer for the readers to see [74].

#### 8. Numerical Results and Discussion

To simulate our model, we use the numerical values given in Table 4.

In this section, we explain the dynamics of the proposed model from graphical results such that the bar chart (2) shows the quantity of sensitivity index of parameters associated with , in which two parameters and have a large effect on . From Figures 3–7, we see the dynamics of by using numerical values of parameters , , , , and . Furthermore, Figures 8–13 show the stable behavior of susceptible, exposed, infected, quarantine, vaccinated, and recovered population, respectively, with a variation in fractal order and using value of fractional order . Moreover, Figures 14–19 show the behavior of the model with each compartment under different fractional order and using value of fractal order . While in Figures 20, 21, 22, 23, 24, and 25, we plot the solution for different fractal-fractional order. Here, we compared the real available reported cases of infection in Pakistan for 200 days from 15th March 2021 to 30th September 2021. The details about the COVID-19 situation in Pakistan can be found in [75, 76], and [77]. Also, some real data in comparison with fractional order simulation has been plotted in [78] recently. We see that the numerical simulation at two different fractional orders coincides very well with the plot of real data as shown in Figure 26. The graphical results also reveal that when the fractal and fractional orders approach 1, then model 1 is reduced to the classical order model.

From Figure 26, it is clear that the real data plot and the simulated data graph are closely related. Moreover, by fitting the real data, we can obtain the numerical values of the parameters of model 1 given in Table 4. Moreover, the model is numerically stable and takes less time and memory during simulation using the Adams–Bashforth scheme due to the nonlinearity and complexity of the problem.

#### 9. Conclusion

We have studied the compartmental model of COVID-19 which has been consisting of susceptible, infected, quarantined, vaccinated, and recovered human populations. The model has been specifically proposed for the implementation of a vaccinated class. The fractal-fractional calculus has been used to understand the dynamics of fractal order and fractional order . We have also investigated some results devoted to local and global stability for the proposed model based on the Jacobian matrix and the Lyapunov function method. Both local and global type stabilities have been demonstrated under the disease-free and endemic equilibrium points by showing that and , respectively. For the sensitivity analysis of each parameter, we have investigated global sensitivity analysis to justify the said feature. Furthermore, some results necessary for numerical stability based on the Ulam–Hyers concept have also been studied. Numerical simulations have been presented by means of the Adams–Bashforth scheme. Some discussion about the convergence and numerical stability of the proposed method has been given in Remark 1. Numerical results have been displayed in different fractional orders graphically to understand the dynamics of the model. We have given a comparison between real reported and simulated data in the case of the infected class. Both curves are closely agreed which shows the efficiency of the proposed numerical method. The given detail can be extended to more complex dynamical systems in the future.

#### Data Availability

All the data used in the paper are included in the manuscript.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest regarding this work.

#### Authors’ Contributions

All authors have contributed equally.

#### Acknowledgments

Kamal Shah, Thabet Abdeljawad, and Bahaaeldin Abdalla are thankful to the Prince Sultan University for paying the APC and support through TAS research lab.