1 Introduction

COVID-19 which is a threatful and terrible disease has been identified initially in Wuhan city of China in December 2019. The said infection transmitted in all over the world in coming few months. The spreading of this little and quickly transmissible virus in the recent time is due to corona virus [1, 2]. In 2020, the disease of COVID-19 is the world big threat that affected nearly all countries and continents around the globe. By the data given by the worldo-meter [3] and WHO [4, 5] shows that nearly 150 million cases of infection occurred while more than five million of population has been died. In the past history of the said virus, it is started in 1965 by the Tyrrell and Bynoe for identification of a virus due to B814 [6]. Such types of viruses had been identified in human organs of embryonic tracheal organ moves taken from the respiration vain of an aged person [7].

Most of scientists, scholars and politicians are trying to stabilize and control the transmission of the pandemic, because the aforesaid disease has killed millions of people around the world in last 2 years. One of the factors of transmitting the said pandemic so rapidly is the migration of infectious population from one place to another. Therefore, locally and globally, some precautionary measures have been implemented. Most of the countries have stopped their air traffic and avoid unnecessary traveling [8]. They also banned crowds and lockdown in the cities to minimize the loss of human lives. In this scenario, the researchers and policies makers are searching to discover a cure or vaccine for the mention disease to stabilize and control it in the coming days. In preparation of vaccine some countries have got succession and now vaccines are available.

To properly controlled this pandemic, it is important to know much about the transmission, symptoms and features of this disease. Implementation of a proper method against the disease out breaking which is the big type of challenge faced by the human population in past history. Therefore, scientists and researchers are trying continuously to model this disease mathematically. In the past, different mathematical models have been developed for infectious diseases, for instance (see the references [9,10,11,12,13,14,15]). Recently, a lot of research work has been published in the form of mathematical models. For reference, we give for instance some published work as [1, 2, 16,17,18,19,20,21,22].

Most of mathematical models which have been investigated in the past were either the system of differential, difference and integration equation having natural or discrete-order. But after fractional calculus has got attention in last few decades, fractional order differential equations (FDEs) applied in excessive numbers to model various real-world problems. FDEs have many applications in various fields of engineering and medical laboratories like physics, business, controlling phenomena, accounting and in biological problems. Therefore, the scientists and researchers increasingly have used FDEs to formulate the real-globe phenomena. Because of the extra degree of choices in fractional derivative which is not present in traditional order operator. Further traditional order derivatives of integer order are not generalized as compared to fractional order which is generalized. Hence fractional order derivative is non-locale in nature and preserves the memory properties which makes it better. Further fractional order derivative of a function produces accumulation of the function which include the corresponding integer order counter part as a special case. Further geometrically it gives spectrum of the function and hence produce the whole density of the function on whom it applies. This is consider the best one, in the conditions where the coming states of models not only related to the present state but may also depend on the past timing of each quantity. For some significance applications see [23,24,25,26]. Due to these properties FDEs not only formulate the problems containing the non-Gaussian nature but can also describe the dynamics for the non-Markovian conditions also. As the natural order derivative and its constituting equations do not give knowledge lying between any of the two consecutive different natural numbers. Therefore, FDEs have been introduced to overcome these limitations. Fractional differential operators have many applications in different areas of mathematical and physical sciences. Liouville, Euler, Reimann and Fourier established some definitions for fractional order derivative during. After that the area has given much more attention. Modern calculus has a lot of applications in the area of mathematical modeling where hereditary characteristic and memorization properties have been studied very easily. Integer order derivative is rarely used to study such behaviors. Non-integer order derivative is the generalization of the natural order derivative having extra degree of freedom as compared the natural order derivative (see [23,24,25,26,27,28,29,30]). Keeping these properties scholars and researchers have taken much interest to study FDEs from different aspects. In the definition of arbitrary order operators, theirs lie a definite integration which predicts physically the area under the function curve or spectrum to generalize it. Integer order differentiation is a specific class of the non-integer order derivative. Although, sufficient contributions have been made by the researchers to analyze the solution of various problems (see [31,32,33,34,35,36,37]). Remarkably, arbitrary order operators have been formulated in different mathematical forms. Fractional differential operators can be classified in two major classes. One is devoted to singular kernel type fractional order differential operators like Reimann-Liouville, Caputo, etc. While the other class is devoted to non-singular type operator where exponential or Mittag-Leffler function play the role of kernel. One of the famous operator of fractional derivative with Mittag-Leffler type kernel is known as \(\mathcal {A}{\mathscr{B}}\mathcal {C}\) introduced by Atangana, Baleanu and Caputo [38] in 2016. This operator replaced the singular kernel by non-singular one [39,40,41]. But this classification has own merits and de-merits. But researchers increasingly used these operators to investigate various real-world problems.

To treat FDEs under various operators for their numerical solution, optimization and numerical analysis, the traditional techniques have been extended for these purposes. For instance decomposition and homotopy perturbation techniques have been previously used to investigate various problems of FDEs (see [42,43,44]). For numerical solution mostly, RKM methods have been applied to various fractional order models. Here, in our work, the fractional Adams Bash-forth method is used for numerical simulation as applied in [45, 46]. This technique is an easy bi-step method which is more powerful than Taylor series, Euler method, and RKM techniques. Moreover, it is rapidly convergent and stable.

The investigation of epidemiological models of infected disease have gained great attention from research point of view. Several scholars have investigated the solution existence and uniqueness of many fractional order models [47,48,49,50]. For the description of the mathematical formulation of COVID-19, and to observe that how this disease impacts the susceptible, infected and quarantined people have been investigated. Some of the researchers have focused on the mathematical perspective of COVID-19. For knowing the dynamics structure and physical behavior of the outbreak of COVID-19, Mathematical models are playing important roles. In the problem presented in [51, 52] contains the susceptible people Sp(t), exposed population Ep(t), infectious density Ip(t), asymptotically infectious people Ap(t), humans recovery population Rp(t), reservoir M(t) and the their interactions have been modeled as

$$\begin{array}{@{}rcl@{}}\left\{\begin{array}{lllllll}& \mathcal{D}_{t}(S_{p}(t))=n_{p}-S_{p}m_{p}-b_{p}S_{p}(I_{p}+\kappa A_{p})-b_{w}MS_{p}, \\& \mathcal{D}_{t}(E_{p}(t))=(I_{p}+\kappa A_{p})b_{p}S_{p}+b_{w}S_{p}M-\omega_{p}E_{p}(1-\delta_{p})-\delta_{p}E_{p}{\omega}^{\prime}_{p}-E_{p}m_{p}, \\& \mathcal{D}_{t}(I_{p}(t))=\omega_{p}E_{p}(1-\delta_{p})-I_{p}(\gamma_{p}+m_{p}),\\& \mathcal{D}_{t}(A_{p}(t))=\delta_{p}{\omega}^{\prime}_{p}E_{p}-(\gamma^{\prime}_{p}+m_{p})A_{p},\\& \mathcal{D}_{t}(R_{p}(t))=\gamma_{p}I_{p}+\gamma^{\prime}_{p}A_{p}-m_{p}\gamma_{p},\\& \mathcal{D}_{t}(M(t))=\varepsilon I_{p}+\sigma A_{p}-\vartheta M,\\& S_{p}(0)=S_{0}, \ \ E_{p}(0)=E_{0},\ \ I_{p}(0)=I_{0}, \ \ A_{p}(0)=A_{0}, \ \ R_{p}(0)=A_{0},\ M(0)=M_{0}. \end{array}\right.\end{array}$$
(1)

The detail of parameters applied in the problem (1), with full explanation is provided in Table 1.

Table 1 Parameters description given in the model (1)
Full size table

Here authors have established some global, local stability by computing basic reproductive numbers. Also using simple integral transform method, they have presented some numerical results.

Motivated from the aforesaid literature and work published in the corresponding area, we consider Model (1) under the \(\mathcal {A}{\mathscr{B}}\mathcal {C}\) fractional order derivative as

$$\begin{array}{@{}rcl@{}} \left\{\begin{array}{llllllll}&{~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(S_{p}(t))=n_{p}-m_{p}S_{p}-b_{p}S_{p}(I_{p}+\kappa A_{p})-b_{w}S_{p}M, \\&{~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(E_{p}(t))=(I_{p}+\kappa A_{p})b_{p}S_{p}+b_{w}S_{p}M-\omega_{p}E_{p}(1-\delta_{p})-\delta_{p}{\omega}^{\prime}_{p}E_{p}-m_{p}E_{p}, \\&{~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(I_{p}(t))=\omega_{p}E_{p}(1-\delta_{p})-I_{p}(\gamma_{p}+m_{p}),\\&{~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(A_{p}(t))=\delta_{p}{\omega}^{\prime}_{p}E_{p}-(\gamma^{\prime}_{p}+m_{p})A_{p},\\&{~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(R_{p}(t))=\gamma_{p}I_{p}+\gamma^{\prime}_{p}A_{p}-m_{p}\gamma_{p},\\&{~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(M(t))=\varepsilon I_{p}+\sigma A_{p}-\vartheta M,\\& S_{p}(0)=S_{0}, \ \ E_{p}(0)=E_{0},\ \ I_{p}(0)=I_{0}, \ \ A_{p}(0)=A_{0}, \ \ R_{p}(0)=A_{0},\\& M(0)=M_{0},\ \ \ 0<r\leq1. \end{array}\right. \end{array}$$
(2)

We establish some appropriate results for existence theory of solution via fixed point approach. Further, we attempt on stability results for the suggested model. Some sensitivity results about the parameters of the model are also discussed. Further, numerical technique of Adams-Bashforth method is used to handle this model (2) for the approximate solution and numerical simulations. Further we testify the numerical interpretation by two sets of data one of Wuhan city reported and other one is reported in about Pakistan. Further we also compared our simulated data and real data in case of infected cases to see the validity of the numerical scheme.

Here we remark some limitations of using mathematical models to understand the mechanism of infections disease or other real-world problems. For instance, models that establish for addressing forecasts are usually designed to produce either short-term or long-term forecasts. Some times models designed for long-term forecasting often do not produce good short-term forecasts and vice versa. Also, the associated factors, assumptions and structure, required for the one purpose often make the model less suitable for the other. To construct an appropriate model is a crucial job, because it is often the only link between the model and the model user. Also in majority cases to verify the model by real data, we often have no access to the afore data or information. In short we say that models are abstractions of reality, because, real-world systems are complex and composed of many interrelated components. For a modeler it is impossible or tedious to include all the comments (see detail in [57]). On the other hand for simulations, different numerical schemes are using to deal mathematical models. The concerned schemes have some short comings. For instance often numerical scheme is stable that we are using but on the other hands it will suffer from convergency. In same fashion it is not necessary that a scheme we use is convergent then it must also be stable. Here we use Adams-Bashforth method to simulate our results. The advantage of the proposed method is that it uses only one additional function evaluation per step and produces preserve high-order accuracy. But the limitation of the said method is the necessity of using another method to start.

2 Method, feasibility and stability analysis

Here in this part, we have to find feasibility and stability analysis of the proposed model. We first here re-collect some required results, definitions from [39, 40].

Definition 2.1

\(\mathcal {A}{\mathscr{B}}\mathcal {C}\) fractional operator for a function Ψ(t) and \({\Psi }(t)\in {\mathscr{H}}^{1}(0,\tau )\) is formulated as:

$$\begin{array}{@{}rcl@{}} {~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathbf{D}_{0}^{r}{\Psi}(t)=\frac{\mathcal{A}\mathcal{B}\mathcal{C}(r)}{1-r}{{\int}_{0}^{t}}\frac{d}{dz}{\Psi}(z){\kappa}_{r}\bigg[\frac{-r}{1-r}\bigg(t-z\bigg)^{r}\bigg]dz. \end{array}$$
(3)

If we change \({\kappa }_{r}\bigg [\frac {-r}{1-r}\bigg (t-z\bigg )^{r}\bigg ]\) to \({\kappa }_{1}=\exp \bigg [\frac {-r}{1-r}\bigg (t-z\bigg )\bigg ],\) in (3), then we will obtain the Caputo-Fabrizo (CF) operator of fractional orders. Further, it is to be noted that

$${~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathbf{D}^{r}[constant]=0.$$

\(\mathcal {A}{\mathscr{B}}\mathcal {C}(r)\) is called normalized mapping as \(\mathcal {A}{\mathscr{B}}\mathcal {C}(0)=\mathcal {A}{\mathscr{B}}\mathcal {C}(1)=1\). Also \({\kappa }_{r}\) represents specific mapping known as Mittag-Leffler which is the general form of the exponential mapping [28,29,30].

Definition 2.2

Consider Ψ ∈ L[0,T], then the fractional order integration in the sense of \(\mathcal {A}{\mathscr{B}}\mathcal {C}\) is as follows:

$$\begin{array}{@{}rcl@{}} {~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathbf{I}_{0}^{r}{\Psi}(t)=\frac{1-r}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}{\Psi}(t)+\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r) {\Gamma}(r)}{{\int}_{0}^{t}} (t-z)^{r-1}{\Psi}(z)dz. \end{array}$$
(4)

Lemma 2.3

([54]) If \(Y(t)\rightarrow 0\) as t = 0, then the solution for 0 < r < 1 of the problem

$$\begin{array}{@{}rcl@{}} {~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathbf{D}_{0}^{r} {\Psi}(t)&=&Y(t),\ t\in [0, T],\\ {\Psi}(0)&=&{\Psi}_{0} \end{array}$$

is given by

$${\Psi}(t)={\Psi}_{0}+\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}Y(t)+\frac{r}{\Gamma(r)\mathcal{A}\mathcal{B}\mathcal{C}(r)}{{\int}_{0}^{t}} (t-z)^{r-1}Y(z)dz.$$

Note: For existence of solution, closed norm space is defined by:

$$\mathbf{Y}=\mathbf{Z}=C([0,T]\times R^{6},R),$$

where Z = C[0,T] under the norm:

$$\|\mathbf{Y}\|=\|{\Psi}\|=\sup_{t\in [0, T]}\left[|S_{p}(t)|+|E_{p}(t)|+|I_{p}(t)|+|A_{p}(t)|+|R_{p}(t)|+|M(t)|\right].$$

The Krasnosilkii’s theorem of fixed point theory is applied for the main result.

Theorem 2.4

[55] Consider A be any convexed subset of Y and consider that F,G are two different operators in the integral equations with

  1. 1.

    Gw + Fw ∈A for all w ∈A;

  2. 2.

    F is contracted operator;

  3. 3.

    G is compact and continuous operator.

Then equation Fw + Gw = w in operator form, has one or more than one solution.

Lemma 2.5

The solution of the proposed problem (2) is bounded in the region of feasibility given by

$${\Psi}=\bigg\{(S_{p}(t), E_{p}(t), I_{p}(t), A_{p}(t), R_{p}(t), M(t))\in \mathbf{R}^{6}_{+}:\ 0\leq N(t)\leq \frac{n_{p}}{m_{p}}\bigg\}.$$

Proof

Let consider

$$N_{p}=S_{p}(t)+E_{p}(t)+I_{p}(t)+A_{p}(t)+R_{p}(t)+M(t).$$

By adding all the equations of (2) we get as

$$\begin{array}{@{}rcl@{}}{ll} \frac{{~}^{\mathcal{A}\mathcal{B}\mathcal{C}}d^{r}{(N_{p})}}{dt^{r}}&\leq& n_{p}-m_{p}N_{p}. \end{array}$$
(5)

Solving Eq. (5), we get

$$\begin{array}{@{}rcl@{}} N_{p}&\leq&\frac{n_{p}}{m_{p}}-C\exp(-m_{p}t), \end{array}$$

or

$$N_{p}(t)\leq\frac{n_{p}}{m_{p}},$$

hence proved.

Next we find the disease free and the pandemic equilibrium points by setting all the equation of system (2) equal zero as

$$\begin{array}{@{}rcl@{}} {~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(S_{p}(t))=0, \\ {~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(E_{p}(t))=0, \\ {~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(I_{p}(t))=0,\\ {~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(A_{p}(t))=0,\\ {~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(R_{p}(t))=0,\\ {~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(M(t))=0,\\ \end{array}$$

or

$$E_{0}(\frac{n_{p}}{m_{p}}, 0, 0, 0, 0, 0).$$

Theorem 2.6

The basic reproductive number is computed as

$$R_{0}=\frac{\delta_{p}{\omega}^{\prime}_{p}(m_{p}+\gamma_{p})(\vartheta \kappa b_{p} n_{p}+n_{p}\delta b_{w})+(1-\delta_{p})\omega_{p}(\gamma^{\prime}_{p}+m_{p})(\vartheta b_{p} n_{p}+n_{p}\varepsilon b_{w})}{\vartheta m_{p}(m_{p}+\gamma_{p})(\gamma^{\prime}_{p}+m_{p})(\delta_{p}(\omega^{\prime}_{p}-\omega_{p})+m_{p}+w_{p})}$$

Proof

For this we take the four equations of model (2) as

$$\begin{array}{@{}rcl@{}} \frac{{~}^{\mathcal{A}\mathcal{B}\mathcal{C}}d^{r}{(N_{p})}}{dt^{r}}&=&\left(\begin{array}{c} (I_{p}+\kappa A_{p})b_{p}S_{p}+b_{w}S_{p}M-\omega_{p}E_{p}(1-\delta_{p})-\delta_{p}{\omega}^{\prime}_{p}E_{p}-m_{p}E_{p} \\ \omega_{p}E_{p}(1-\delta_{p})-I_{p}(\gamma_{p}+m_{p}) \\ \delta_{p}{\omega}^{\prime}_{p}E_{p}-(\gamma^{\prime}_{p}+m_{p})A_{p} \\ \varepsilon I_{p}+\sigma A_{p}-\vartheta M \end{array} \right). \end{array}$$

We define F and V as follows

$$\begin{array}{@{}rcl@{}} F=\left(\begin{array}{c} (I_{p}+\kappa A_{p})b_{p}S_{p}+b_{w}S_{p}M \\ 0 \\ 0 \\ 0 \end{array} \right) \end{array}$$

and

$$\begin{array}{@{}rcl@{}} V=\left(\begin{array}{c} \omega_{p}E_{p}(1-\delta_{p})+\delta_{p}{\omega}^{\prime}_{p}E_{p}+m_{p}E_{p} \\ \omega_{p}E_{p}(1-\delta_{p})-I_{p}(\gamma_{p}+m_{p}) \\ \delta_{p}{\omega}^{\prime}_{p}E_{p}-(\gamma^{\prime}_{p}+m_{p})A_{p} \\ \varepsilon I_{p}+\sigma A_{p}-\vartheta M \end{array} \right). \end{array}$$

Next, taking the Jacobian of F and V w.r.t, and putting the value of E0, we get

$$\begin{array}{@{}rcl@{}} \mathcal{F}=\left(\begin{array}{cccc} 0 & \frac{b_{p} n_{p}}{m_{p}} & \frac{\kappa b_{p} n_{p}}{m_{p}} & \frac{b_{w} n_{p}}{m_{p}} \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array} \right) \end{array}$$

and

$$\begin{array}{@{}rcl@{}} \mathcal{V}=\left(\begin{array}{cccc} \delta_{p}{\omega}^{\prime}_{p}+\omega_{p}(1-\delta_{p})+m_{p} & 0 & 0 & 0 \\ \omega_{p}(\delta_{p}-1) & m_{p}+\gamma_{p} & 0 & 0 \\ -\delta_{p}{\omega}^{\prime}_{p} & 0 & {\gamma}^{\prime}_{p}+m_{p} & 0 \\ 0 & -\varepsilon_{p} & -\delta & \vartheta \end{array} \right). \end{array}$$

Then the dominant eigen value of \(\mathcal {F}\mathcal {V}^{-1}=\rho (\mathcal {F}\mathcal {V}^{-1})\) is called basic reproduction number R0 and hence equal to

$$\begin{array}{@{}rcl@{}} R_{0}=\frac{\delta_{p}{\omega}^{\prime}_{p}(m_{p}+\gamma_{p})(\vartheta \kappa b_{p} n_{p}+n_{p}\delta b_{w})+(1-\delta_{p})\omega_{p}(\gamma^{\prime}_{p}+m_{p})(\vartheta b_{p} n_{p}+n_{p}\varepsilon b_{w})}{\vartheta m_{p}(m_{p}+\gamma_{p})(\gamma^{\prime}_{p}+m_{p})(\delta_{p}(\omega^{\prime}_{p}-\omega_{p})+m_{p}+w_{p})}. \end{array}$$
(6)

Hence proved

Theorem 2.7

E0 is locally asymptotically stable if R0 < 1.

Proof

The derivation of the theorem can be obtained by taking Jacobian of the system (2) and putting \(E_{0}=(\frac {n_{p}}{m_{p}},0,0,0,0,0)\), one has

$$\begin{array}{@{}rcl@{}} \mathcal{J}=\left(\begin{array}{cccccc} -m_{p} & 0 & \frac{-b_{p} n_{p}}{m_{p}} & \frac{-\kappa b_{p} n_{p}}{m_{p}} & 0 & \frac{-b_{w} n_{p}}{m_{p}} \\ 0 & -m_{p}-\delta_{p}\omega^{\prime}_{p}-(1-\delta_{p})\omega_{p} & \frac{b_{p} n_{p}}{m_{p}}& \frac{\kappa b_{p} n_{p}}{m_{p}} & 0 & \frac{b_{w} n_{p}}{m_{p}} \\ 0 & (1-\delta_{p})\omega_{p} & -m_{p}-\gamma_{p} & 0 & 0 & 0 \\ 0 & \delta_{p} \omega^{\prime}_{p}& 0 & -m_{p}-\gamma^{\prime}_{p} & 0 & 0 \\ 0 & 0 & \gamma_{p} & \gamma^{\prime}_{p} & -m_{p} & 0 \\ 0 & 0 & \varepsilon & \delta & 0 & -\vartheta \end{array} \right). \end{array}$$

In the above matrix, two of the eigen values on the main diagonal are negative, while the rest of the eigen values can be computed by characteristic equation as

$$\begin{array}{@{}rcl@{}} {\Lambda}^{4}+b_{1}{\Lambda}^{3}+b_{2}{\Lambda}^{2}+b_{3}{\Lambda}+b_{4}=0. \end{array}$$
(7)

Here

$$\begin{array}{@{}rcl@{}} b_{1}&=&\gamma^{\prime}_{p}+\delta+\delta_{p}\omega^{\prime}_{p}+(1-\delta_{p})\omega_{p}+3m_{p}+\gamma_{p},\\ b_{2}&=&\underbrace{(m_{p}+\gamma_{p})(\delta_{p}\omega^{\prime}_{p}+(1-\delta_{p})\omega_{p}+m_{p})-b_{p}(1-\delta_{p})\omega_{p}}\\ &+&\underbrace{(m_{p}+\gamma_{p})(\delta_{p}(\omega^{\prime}_{p}-\omega_{p})+m_{p}-\omega_{p})-\kappa b_{p}\delta_{p}\omega^{\prime}_{p}}\\ &+&\vartheta(\gamma^{\prime}_{p}+m_{p})+(m_{p}+\gamma_{p})(m_{p}+\gamma^{\prime}_{p})+\vartheta(\delta_{p}(\omega^{\prime}_{p}-\omega_{p})+m_{p}+\omega_{p})+\vartheta(m_{p}+\gamma_{p}),\\ b_{3}&=&\vartheta(\delta_{p}\omega^{\prime}_{p}+(1-\delta_{p})\omega_{p}+m_{p})[(\gamma^{\prime}_{p}+m_{p})(1-R_{1})+(m_{p}+\gamma_{p})(1-R_{2})]\\ &+&(m_{p}+\gamma_{p})\underbrace{(\delta(\gamma^{\prime}_{p}+m_{p})-(\kappa b_{p}\delta_{p}\omega^{\prime}_{p}))}+b_{p}\delta_{p}\omega_{p}(\gamma^{\prime}_{p}+m_{p})\\ &+&(\gamma^{\prime}_{p}+m_{p})((\gamma_{p}+m_{p})\delta_{p}(\omega^{\prime}_{p}-\omega_{p})+m_{p}+\omega_{p})-b_{p}\omega_{p},\\ b_{4}&=&\vartheta(\gamma_{p}+m_{p})(\gamma^{\prime}_{p}+m_{p})(\delta+\delta_{p}\omega^{\prime}_{p}+(1-\delta_{p})\omega_{p}+m_{p})(1-R_{0}), \end{array}$$

where R0 = R1 + R2 as follows

$$R_{1}=\frac{\delta_{p}{\omega}^{\prime}_{p}(\vartheta \kappa b_{p} n_{p}+n_{p}\delta b_{w})}{\vartheta m_{p}(\gamma^{\prime}_{p}+m_{p})(\delta_{p}(\omega^{\prime}_{p}-\omega_{p})+m_{p}+w_{p})}$$

and

$$R_{2}=\frac{(1-\delta_{p})\omega_{p}(\gamma_{p}+m_{p})(\vartheta b_{p} n_{p}+n_{p}\varepsilon b_{w})}{\vartheta m_{p}(\gamma_{p}+m_{p})(\delta_{p}(\omega^{\prime}_{p}-\omega_{p})+m_{p}+w_{p})}.$$

In the above characteristic equation, the terms which are underlines are less than R0, also b4 is positive if R0 < 1. Further if R1 < 1 and R2 < 1, then b3 will be positive. Hence all the coefficients are positive being the conditions for Routh-Hurwitz criteria [53]. Hence E0 is locally asymptotically stable.

Next we have to find the pandemic equilibrium point as

$$E^{*}=(S^{*}_{P}, E^{*}_{P}, I^{*}_{P}, A^{*}_{P}, R^{*}_{P},M^{*})$$

and

$$\begin{array}{@{}rcl@{}} S^{*}_{p}&=&\frac{n_{p}}{\Lambda+m_{p}},\\ E^{*}_{p}&=&\frac{\Lambda S^{*}_{p}}{\delta_{p}\omega^{\prime}_{p}-\delta_{p}\omega_{p}+m_{p}+\omega_{p}},\\ I^{*}_{p}&=&\frac{E^{*}_{P}(1-\delta_{p})\omega_{p}}{\gamma_{p}+m_{p}},\\ A^{*}_{P}&=&\frac{E^{*}_{p}\delta_{p}\omega^{\prime}_{p}}{\gamma^{\prime}_{p}+m_{p}},\\ R^{*}_{P}&=&\frac{A^{*}_{p}\gamma^{\prime}_{p}+I^{*}_{p}\gamma_{p}}{+m_{p}},\\ M^{*}&=&\frac{A^{*}_{p}\delta+I^{*}_{p}\varepsilon}{\vartheta}. \end{array}$$

Here

$${\Lambda}=\frac{b_{p}(\kappa n_{p} A^{*}_{p}+m_{p}I^{*}_{p})}{m_{p}(S^{*}_{p}+E^{*}_{p}+I^{*}_{p}+A^{*}_{p}+R^{*}_{p})}+\frac{b_{p}M^{*}}{m_{p}},$$

satisfying the given equation

$$\begin{array}{@{}rcl@{}} P({\Lambda}^{*})=a_{1}({\Lambda}^{*})2+a_{2}({\Lambda}^{*})=0. \end{array}$$

Here

$$a_{1}=\vartheta(m_{p}+\gamma_{p})(m_{p}+\gamma^{\prime}_{p})(\delta_{p}(\omega^{\prime}_{p}-\omega_{p})+m_{p}+\omega_{p})$$
$$a_{2}=\vartheta m_{p}(m_{p}+\gamma_{p})(m_{p}+\gamma^{\prime}_{p})(\delta_{p}(\omega^{\prime}_{p}-\omega_{p})+m_{p}+\omega_{p})(1-R_{0})$$

As a1 > 0,a2 ≥ 0 if R0 < 1, then \({\Lambda }^{*}=\frac {-a_{2}}{a_{1}}\leq 0\). Hence no pandemic equilibrium will lie if R0 ≤ 1. This implies that the endemic equilibrium exists and stable if R0 > 1.

3 Existence, uniqueness of solution and numerical simulations

It is natural to ask whether a dynamical system that we are investigating exists or not in reality. Fixed point theory answer this question. We examine our considered problem (2) for existence results about the solution. Regarding this, we write the right sides of our problem (2) as:

$$\begin{array}{@{}rcl@{}}\left\{\begin{array}{lllllll}& \mathrm{F}_{1}(t, S_{p}, E_{p}, I_{p}, A_{p}, R_{p}, M)=n_{p}-S_{p}m_{p}-b_{p}S_{p}(I_{p}+\kappa A_{p})-S_{p}Mb_{w}, \\& \mathrm{F}_{2}(t, S_{p},E_{p}, I_{p},A_{p}, R_{p}, M)=(I_{p}+\kappa A_{p})b_{p}S_{p}+b_{w}S_{p}M-(1-\delta_{p})\omega_{p}E_{p}-E_{p}\delta_{p}{\omega}^{\prime}_{p}-E_{p}m_{p}, \\& \mathrm{F}_{3}(t, S_{p},E_{p}, I_{p},A_{p}, R_{p}, M)=(1-\delta_{p})\omega_{p}E_{p}-(\gamma_{p}+m_{p})I_{p},\\& \mathrm{F}_{4}(t, S_{p},E_{p}, I_{p},A_{p}, R_{p}, M)=\delta_{p}{\omega}^{\prime}_{p}E_{p}-(\gamma^{\prime}_{p}+m_{p})A_{p},\\& \mathrm{F}_{5}(t, S_{p},E_{p}, I_{p},A_{p}, R_{p}, M)=\gamma_{p}I_{p}+\gamma^{\prime}_{p}A_{p}-m_{p}\gamma_{p},\\& \mathrm{F}_{6}(t, S_{p},E_{p}, I_{p},A_{p}, R_{p}, M)=\varepsilon I_{p}+\sigma A_{p}-\vartheta M,\\& S_{p}(0)=S_{0}, \ E_{p}(0)=E_{0},\ I_{p}(0)=I_{0}, \ A_{p}(0)=A_{0}, \ R_{p}(0)=A_{0},\ M(0)=M_{0}. \end{array}\right. \end{array}$$
(8)

To symbolize the system (2) by using (8) as follows

$$\begin{array}{@{}rcl@{}} {~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathbf{D}_{+0}^{r} \mathcal{Y}(t)&=&{\Omega}(t, \mathcal{Y}(t)),\ t\in [0, \tau], \ 0<r\leq1,\\ \mathcal{Y}(0)&=&\mathcal{Y}_{0}. \end{array}$$
(9)

By applying integral in sense of \(\mathcal {A}{\mathscr{B}}\mathcal {C}\) and by using lemma 2.3 we get

$$\mathcal{Y}(t)=\mathcal{Y}_{0}(t)+\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\bigg[{\Omega}(t, \mathcal{Y}(t))\bigg]+\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}{{\int}_{0}^{t}} (t-z)^{r-1} {\Omega}(z, \mathcal{Y}(z))dz,$$
(10)

where,

$$\mathcal{Y}(t)=\left\{\begin{array}{llllll}&S_{p}(t)\\& E_{p}(t)\\& I_{p}(t)\\& A_{p}(t)\\& R_{p}(t)\\& M(t) \end{array} \right., \ \ \ \mathcal{Y}_{0}(t)=\left\{\begin{array}{llllll}&S_{0}\\& E_{0}\\& I_{0}\\& A_{0}\\& R_{0}\\& M_{0} \end{array} \right., \ \ \ {\Omega}(t, \mathcal{Y}(t))=\left\{\begin{array}{llllll}&\mathrm{F}_{1}(S_{p},E_{p}, I_{p},A_{p}, R_{p}, M,t)\\& \mathrm{F}_{2}(t,S_{p},E_{p}, I_{p},A_{p}, R_{p}, M)\\& \mathrm{F}_{3}(t,S_{p},E_{p}, I_{p},A_{p}, R_{p}, M)\\& \mathrm{F}_{4}(t,S_{p},E_{p}, I_{p},A_{p}, R_{p}, M)\\& \mathrm{F}_{5}(t,S_{p},E_{p}, I_{p},A_{p}, R_{p}, M)\\& \mathrm{F}_{6}(t,S_{p},E_{p}, I_{p},A_{p}, R_{p}, M). \end{array} \right.$$
(11)

Using (9) and define operators F,G by using (10) as

$$\begin{array}{@{}rcl@{}} \mathbf{F}(\mathcal{Y})&=& \mathcal{Y}_{0}(t)+\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\bigg[{\Omega}(t, \mathcal{W}(t))\bigg],\\ \mathbf{G}(\mathcal{Y})&=& \frac{r}{\Gamma(r)\mathcal{A}\mathcal{B}\mathcal{C}(r)}{{\int}_{0}^{t}} (t-z)^{r-1}{\Omega}(Z, \mathcal{Y}(z))dz. \end{array}$$
(12)

Witting the growth condition and Lipschitz condition for solution’s existence and uniqueness as given below.

(B1):

Let we have a constants AΩ,EΩ, as:

$$|{\Omega}(t, \mathcal{Y}(t))|\leq A_{\Omega}|\mathcal{Y}|+E_{\Omega}.$$
(B2):

Let we have a constants LΩ > 0, as for all \(\mathcal {Y},\ {\bar {Y}}\in \mathbf {Y}\) as:

$$|{\Omega}(t, \mathcal{Y})-{\Omega}(t, {\bar{Y}})|\leq L_{\Omega}[|\mathcal{Y}|- {\bar{Y}}|].$$

Theorem 3.1

Under hypothesis (B1,B2), the problem (10) has at least one solution which implies that the proposed model (2) has at least one solution if \(\frac {(1-r)}{\mathcal {A}{\mathscr{B}}\mathcal {C}(r)}L_{\Omega }<1\).

Proof

The theorem can be proved by using the following two steps.

  1. Step I:

    Consider \({\bar {Y}}\in \mathbf {B}\) and \(\mathbf {B}=\{\mathcal {Y}\in \mathbf {Y}:\ \|\mathcal {Y}\|\leq \sigma , \sigma >0\}\) is convex and close set. Then by F in (12), we obtain

    $$\begin{array}{@{}rcl@{}} \|\mathbf{F}(\mathcal{Y})-\mathbf{F}({\bar{Y}})\| &=& \frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\max_{t\in [0, \tau]}\bigg|{\Omega}(t, \mathcal{Y}(t))-{\Omega}(t, {\bar{Y}}(t))\bigg| \\ &\leq& \frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}L_{\Omega}\|\mathcal{Y}-{\bar{Y}}\|. \end{array}$$
    (13)

    Hence, F is contracted.

  2. Step-II:

    To show that G is compact relative, it is enough to show that G is bounded and equi-continuous. Clearly, G is defined on their domain as Ω is defined on domain and for any \(\mathcal {Y}\in \mathbf {B}\), we follow

    $$\begin{array}{@{}rcl@{}} \|\mathbf{G}(\mathcal{Y})\|&=&\max_{t\in [0, \tau]}\|\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}{{\int}_{0}^{t}} (t-z)^{r-1}{\Omega}(z, \mathcal{Y}(z))dz\|\\ &\leq& \frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}{\int}_{0}^{\tau}(\tau-z)^{r-1}|{\Omega}(z, \mathcal{Y}(z))|dz\\ &\leq&\frac{\tau^{r}}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}[A_{\Omega}\sigma+E_{\Omega}]. \end{array}$$
    (14)

So, from (14) it is clear that G have bounds. Further, for equi-continuous let t1 > t2 ∈ [0,τ], we continue as

$$\begin{array}{@{}rcl@{}} |\mathbf{G}(\mathcal{Y}(t_{2})-\mathbf{G}(\mathcal{Y}(t_{1})|&=&\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}\bigg|{\int}_{0}^{t_{2}}(t_{2}-z)^{r-1}{\Omega}(z, \mathcal{Y}(z))dz-{\int}_{0}^{t_{1}}(t_{1}-z)^{r-1}{\Omega}(z, \mathcal{Y}(z))dz\bigg|\\ &\leq&\frac{[A_{\Omega}\sigma+E_{\Omega}]}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}[{t_{2}^{r}}-{t_{1}^{r}}]. \end{array}$$
(15)

Equation (15) implies that as \(t_{2}\rightarrow t_{1}\) then the right side will approaches to zero. As, G is continuous and hence

$$|\mathbf{G}(\mathcal{Y}(t_{2})-\mathbf{G}(\mathcal{Y}(t_{1})|\rightarrow 0,\ \text{as}\ t_{2}\rightarrow t_{1}.$$

Hence as G have bounds and are continuous so

$$\|\mathbf{G}(\mathcal{Y}(t_{2})-\mathbf{G}(\mathcal{Y}(t_{1})\|\rightarrow 0,\ \text{as}\ t_{2}\rightarrow t_{1}.$$

Thus, G have bounds and equi-continuous operator. Also, from theorem of Arzelá-Ascoli, the operator G is relative compact and hence continuous completely. Thus, from Theorem 3.1, the integration Equation (10) has Als at least one solution and therefore the proposed problem has at least one solution.

For unique solution we give the given theorem.

Theorem 3.2

By hypothesis (B2), the integral Equation (10) has unique solution which yields that the system under consideration (2) has unique solution if:

$$\bigg[\frac{(1-r)L_{\Omega}}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}+\frac{\tau^{r} L_{\Omega}}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}\bigg]<1.$$

Proof

Consider the mapping \(\mathbf {T}:\mathbf {Y}\rightarrow \mathbf {Y}\) defined by

$$\begin{array}{@{}rcl@{}} \mathbf{T}\mathcal{Y}(t)&=&\mathcal{Y}_{0}(t)+\bigg[{\Omega}(t, \mathcal{Y}(t))-{\Omega}_{0}(t)\bigg]\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\\&& +\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}{{\int}_{0}^{t}} (t-z)^{r-1} {\Omega}(z, \mathcal{Y}(z))dz,\ t\in [0, \tau]. \end{array}$$
(16)

Let \(\mathcal {Y}, {\bar {Y}} \in \mathbf {Y}\), then one can take

$$\begin{array}{@{}rcl@{}} \|\mathbf{T}\mathcal{Y}-\mathbf{T}{\bar{Y}}\|&\leq&\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\max_{t\in [0, \tau]}\bigg|{\Omega}(t, \mathcal{Y}(t))-{\Omega}(t, {\bar{Y}}(t))\bigg|\\ &+&\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}\max_{t\in [0, \tau]}\bigg|{{\int}_{0}^{t}} (t-z)^{r-1} {\Omega}(z, \mathcal{Y}(z))dz-{{\int}_{0}^{t}} (t-z)^{r-1} {\Omega}(z, {\bar{Y}}(z))dz\bigg|\\ &\leq&{\Upsilon}\|\mathcal{Y}-{\bar{Y}}\|, \end{array}$$
(17)

where

$$\begin{array}{@{}rcl@{}} {\Upsilon}=\bigg[\frac{(1-r)L_{\Omega}}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}+\frac{\tau^{r} L_{\Omega}}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}\bigg]. \end{array}$$
(18)

Hence, T is contracted from (17). So, the integration Equation (10) has one root. This implies that the problem (2) has one solution.

For approximate solution, we continue this part of manuscript to the proposed fractional order (2) model in sense of ABC operator. The iterative technique are then simulated on different fractional orders. For this, we use the arbitrary order AB iterative technique [56] to find the numerical scheme for the graphical representation of the problem (2). For model (8) we develop a numerical scheme as

$$\begin{array}{@{}rcl@{}}\left\{\begin{array}{lllllll}&{~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(S_{p}(t))=\mathrm{F}_{1}(S_{p}(t),t)=n_{p}-m_{p}S_{p}-b_{p}S_{p}(I_{p}+\kappa A_{p})-b_{w}S_{p}M, \\&{~}^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(E_{p}(t))= \mathrm{F}_{2}(E_{p}(t),t)=(I_{p}+\kappa A_{p})b_{p}S_{p}+b_{w}S_{p}M-\omega_{p}E_{p}(1-\delta_{p})-\delta_{p}{\omega}^{\prime}_{p}E_{p}-E_{p}m_{p}, \\& ^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(I_{p}(t))=\mathrm{F}_{3}(I_{p}(t),t)=\omega_{p}E_{p}(1-\delta_{p})-I_{p}(\gamma_{p}+m_{p}),\\& ^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(A_{p}(t))=\mathrm{F}_{4}(A_{p}(t),t)=\delta_{p}{\omega}^{\prime}_{p}E_{p}-(\gamma^{\prime}_{p}+m_{p})A_{p},\\& ^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(R_{p}(t))=\mathrm{F}_{5}(R_{p}(t), t)=\gamma_{p}I_{p}+\gamma^{\prime}_{p}A_{p}-m_{p}\gamma_{p},\\& ^{\mathcal{A}\mathcal{B}\mathcal{C}}\mathcal{D}_{t}^{r}(M(t))=\mathrm{F}_{6}(M(t), t)=\varepsilon I_{p}+\sigma A_{p}-\vartheta M,\\& S_{p}(0)=S_{0}, \ E_{p}(0)=E_{0}, \ I_{p}(0)=I_{0}, \ A_{p}(0)=A_{0}, \ R_{p}(0)=A_{0},\ M(0)=M_{0}. \end{array}\right. \end{array}$$
(19)

Integrating first equation of (19) in \(\mathcal {A}{\mathscr{B}}\mathcal {C}\) approach, we get

$$\begin{array}{@{}rcl@{}} S_{p}(t)-S_{p}(0)=\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\bigg[\mathrm{F}_{1}(S_{p}(t), t)\bigg]+\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}{{\int}_{0}^{t}} (t-z)^{r-1} \mathrm{F}_{1}(S_{p}(z),z)dz. \end{array}$$

Consider t = ti+ 1 for i = 0,1,2⋯, it follows that

$$\begin{array}{@{}rcl@{}} S_{p}(t_{i+1})-S_{p}(0)&=&\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\bigg[\mathrm{F}_{1}(S_{p}(t_{i}), t_{i})\bigg]\\&& +\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}{\int}_{0}^{t_{i+1}} (t_{i+1}-z)^{r-1} \mathrm{F}_{1}(S_{p}(z),z)dz,\\ &=&\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\bigg[\mathrm{F}_{1}(S_{p}(t_{i}), t_{i})\bigg]\\&& +\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)} \sum\limits_{q=0}^{i}{\int}_{q}^{t_{q+1}} (t_{i+1}-z)^{r-1} \mathrm{F}_{1}(S_{p}(z),z)dz. \end{array}$$

Next, approximating the mapping F1(Sp(t),t) on time interval [tq,tq+ 1], by the interpolating expression as follows:

$$F_{1}(S_{p}(t), t)\cong \frac{F_{1}(S_{p}(t_{q}), t_{q})}{\Delta}(t-t_{q-1})+ \frac{F_{1}(S_{p}(t_{q-1}), t_{q-1})}{\Delta}(t-t_{q})$$

or

$$\begin{array}{@{}rcl@{}} S_{p}(t_{i+1})&=&S_{p}(0)+\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\bigg[\mathrm{F}_{1}(S_{p}(t_{i}), t_{i})\bigg]\\&& +\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}{\sum}_{q=0}^{i}\bigg(\frac{F_{1}(S_{p}(t_{q}), t_{q})}{\Delta}{\int}_{q}^{t_{q+1}}(t-t_{q-1}) (t_{i+1}-t)^{r-1}dt\\&& -\frac{F_{1}(S_{p}(t_{q-1}), t_{q-1})}{\Delta}{\int}_{q}^{t_{q+1}}(t-t_{q}) (t_{i+1}-t)^{r-1}dt\bigg)\\ &=&S_{p}(0)+\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\bigg[\mathrm{F}_{1}(S_{p}(t_{i}), t_{i})\bigg]\\&& +\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r){\Gamma}(r)}{\sum}_{q=0}^{i}\bigg(\frac{F_{1}(S_{p}(t_{q}), t_{q})}{\Delta}I_{q-1,r}-\frac{F_{1}(S_{p}(t_{q-1}), t_{q-1})}{\Delta}I_{q,r}\bigg). \end{array}$$
(20)

Now the integrals Iq− 1,r and Iq,r can be calculated as follow:

$$\begin{array}{@{}rcl@{}} I_{q-1,r}&=&{\int}_{q}^{t_{q+1}}(t-t_{q-1}) (t_{i+1}-t)^{r-1}dt\\ &=&-\frac{1}{r}\bigg[(t_{q+1}-t_{q-1})(t_{i+1}-t_{q+1})^{r}-(t_{q}-t_{q-1})(t_{i+1}-t_{q})^{r}\bigg]\\&& -\frac{1}{r(r-1)}\bigg[(t_{i+1}-t_{q+1})^{r+1}-(t_{i+1}-t_{q})^{r+1}\bigg], \end{array}$$

and

$$\begin{array}{@{}rcl@{}} I_{q,r}&=&{\int}_{q}^{t_{q+1}}(t-t_{q}) (t_{i+1}-t)^{r-1}dt\\ &=&-\frac{1}{r}\bigg[(t_{q+1}-t_{q})(t_{i+1}-t_{q+1})^{r}\bigg]\\&& -\frac{1}{r(r-1)}\bigg[(t_{i+1}-t_{q+1})^{r+1}-(t_{i+1}-t_{q})^{r+1}\bigg], \end{array}$$

put tq = qΔ, we get

$$\begin{array}{@{}rcl@{}} I_{q-1,r}&=&-\frac{{\Delta}^{r+1}}{r}\bigg[(q+1-(q-1))(i+1-(q+1))^{r}-(q-(q-1))(i+1-q)^{r}\bigg]\\&& -\frac{{\Delta}^{r+1}}{r(r-1)}\bigg[(i+1-(q+1))^{r+1}-(i+1-q)^{r+1}\bigg],\\ &=&\frac{{\Delta}^{r+1}}{r(r-1)}\bigg[-2(r+1)(i-q)^{r}+(r+1)(i+1-q)^{r}-(i-q)^{r+1}+(i+1-q)^{r+1}\bigg],\\ &=&\frac{{\Delta}^{r+1}}{r(r-1)}\bigg[(i-q)^{r}(-2(r+1)-(i-q))+(i+1-q)^{r}(r+1+i+1-q)\bigg],\\ &=&\frac{{\Delta}^{r+1}}{r(r-1)}\bigg[(i+1-q)^{r}(i-q+2+r)-(i-q)^{r}(i-q+2+2r)\bigg], \end{array}$$
(21)

and

$$\begin{array}{@{}rcl@{}} I_{q,r}&=&-\frac{{\Delta}^{r+1}}{r}\bigg[(q+1-q)(i+1-(q+1))^{r}\bigg]-\frac{{\Delta}^{r+1}}{r(r-1)}\bigg[(i+1-(q+1))^{r+1}-(i+1-q)^{r+1}\bigg],\\ &=&\frac{{\Delta}^{r+1}}{r(r-1)}\bigg[-(r+1)(i-q)^{r}-(i-q)^{r+1}+(i+1-q)^{r+1}\bigg],\\ &=&\frac{{\Delta}^{r+1}}{r(r-1)}\bigg[(i-q)^{r}(-(r+1)-(i-q))+(i+1-q)^{r+1}\bigg],\\ &=&\frac{{\Delta}^{r+1}}{r(r-1)}\bigg[(i+1-q)^{r+1}-(i-q)^{r}(i-q+1+r)\bigg], \end{array}$$
(22)

substitute (21) and (22) in (20), we get

$$\begin{array}{@{}rcl@{}} S_{p}(t_{i+1})&=&S_{p}(0)+\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\bigg[\mathrm{F}_{1}(S_{p}(t_{i}), t_{i})\bigg]\\&& +\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r)} \sum\limits_{q=0}^{i}\bigg(\frac{F_{1}(S_{p}(t_{q}), t_{q})}{\Gamma{(r+2)}}{\Delta}^{r}\bigg[(i+1-q)^{r}(i-q+2+r)\\&& -(i-q)^{r}(i-q+2+2r)\bigg]\\&& -\frac{F_{1}(S_{p}(t_{q-1}), t_{q-1})}{\Gamma{(r+2)}}{\Delta}^{r}[(i+1-q)^{r+1}-(i-q)^{r}(i-q+1+r)]\bigg). \end{array}$$

Similarly, the proposed method can be used for the remaining five equations of (19) to form general algorithms as

$$\begin{array}{@{}rcl@{}} E_{p}(t_{i+1})&=&E_{p}(0)+\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\bigg[\mathrm{F}_{2}(E_{p}(t_{i}), t_{i})\bigg]\\&& +\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r)} \sum\limits_{q=0}^{i}\bigg(\frac{F_{2}(E_{p}(t_{q}), t_{q})}{\Gamma{(r+2)}}{\Delta}^{r}\bigg[(i+1-q)^{r}(i-q+2+r)\\&& -(i-q)^{r}(i-q+2+2r)\bigg]\\&& -\frac{F_{2}(E_{p}(t_{q-1}), t_{q-1})}{\Gamma{(r+2)}}{\Delta}^{r}[(i+1-q)^{r+1}-(i-q)^{r}(i-q+1+r)]\bigg). \end{array}$$
$$\begin{array}{@{}rcl@{}} I_{p}(t_{i+1})&=&I_{p}(0)+\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\bigg[\mathrm{F}_{3}(I_{p}(t_{i}), t_{i})\bigg]\\&& +\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r)} \sum\limits_{q=0}^{i}\bigg(\frac{F_{3}(I_{p}(t_{q}), t_{q})}{\Gamma{(r+2)}}{\Delta}^{r}\bigg[(i+1-q)^{r}(i-q+2+r)\\&& -(i-q)^{r}(i-q+2+2r)\bigg]\\&& -\frac{F_{3}(I_{p}(t_{q-1}), t_{q-1})}{\Gamma{(r+2)}}{\Delta}^{r}[(i+1-q)^{r+1}-(i-q)^{r}(i-q+1+r)]\bigg). \end{array}$$
$$\begin{array}{@{}rcl@{}} A_{p}(t_{i+1})&=&A_{p}(0)+\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\bigg[\mathrm{F}_{4}(A_{p}(t_{i}), t_{i})\bigg]\\&& +\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r)} \sum\limits_{q=0}^{i}\bigg(\frac{F_{4}(A_{p}(t_{q}), t_{q})}{\Gamma{(r+2)}}{\Delta}^{r}\bigg[(i+1-q)^{r}(i-q+2+r)\\&& -(i-q)^{r}(i-q+2+2r)\bigg]\\&& -\frac{F_{4}(A_{p}(t_{q-1}), t_{q-1})}{\Gamma{(r+2)}}{\Delta}^{r}[(i+1-q)^{r+1}-(i-q)^{r}(i-q+1+r)]\bigg). \end{array}$$
$$\begin{array}{@{}rcl@{}} R_{p}(t_{i+1})&=&R_{p}(0)+\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\bigg[\mathrm{F}_{5}(R_{p}(t_{i}), t_{i})\bigg]\\&& +\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r)} \sum\limits_{q=0}^{i}\bigg(\frac{F_{5}(R_{p}(t_{q}), t_{q})}{\Gamma{(r+2)}}{\Delta}^{r}\bigg[(i+1-q)^{r}(i-q+2+r)\\&& -(i-q)^{r}(i-q+2+2r)\bigg]\\&& -\frac{F_{5}(R_{p}(t_{q-1}), t_{q-1})}{\Gamma{(r+2)}}{\Delta}^{r}[(i+1-q)^{r+1}-(i-q)^{r}(i-q+1+r)]\bigg). \end{array}$$
$$\begin{array}{@{}rcl@{}} M(t_{i+1})&=&M(0)+\frac{(1-r)}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\bigg[\mathrm{F}_{6}(M(t_{i}), t_{i})\bigg]\\&& +\frac{r}{\mathcal{A}\mathcal{B}\mathcal{C}(r)}\sum\limits_{q=0}^{i}\bigg(\frac{F_{6}(M(t_{q}), t_{q})}{\Gamma{(r+2)}}{\Delta}^{r}\bigg[(i+1-q)^{r}(i-q+2+r)\\&& -(i-q)^{r}(i-q+2+2r)\bigg]\\&& -\frac{F_{6}(M(t_{q-1}), t_{q-1})}{\Gamma{(r+2)}}{\Delta}^{r}[(i+1-q)^{r+1}-(i-q)^{r}(i-q+1+r)]\bigg). \end{array}$$

We apply the above procedure to simulate the proposed model in the next section.

4 Experimental results and discussion

The numerical simulation is obtained from the results using data given in Table 2, from [51] with some updating assuming data. Further, the initial data are in million, we have taken in percentage as Sp(0) = 8.465518, Ep(0) = 0.2,Ip(0) = 0.0002, Ap(0) = 0.0002, Rp(0) = 0.0002, M(0) = 0.055, and the values of parameters are given in Table 2.

Table 2 Description of the parameters given in model (1) for Wuhan
Full size table

From Fig. 1, it is observed that increasing rate of transmission will decrease the number of susceptible individuals, and in return, it will increase number of infected population. In other words isolation and keeping social distance will greatly help in controlling the current outbreak for further spreading. The decrease of susceptible has been shown on different fractional order. The order of derivatives has also produced certain impact on the process, initially at smaller order the process of decay is faster than the higher order and vice versa. After 60 days the susceptible population decreases and then going towards convergency and stability.

Fig. 1
figure 1

Behavior of Susceptible population Sp(t) at various arbitrary order r of the proposed system (2) for h = 0.1,bp = 0.05

Full size image

From Fig. 2, we conclude that initially the number of exposed individuals is growing up to 50 days for about all orders of derivatives. After that the number decreases gradually, but at this time the decrease occurs differently at different order of derivative. It means that exposed population increases as the symptoms is recognized initially, when the outbreak of pandemic starts. After 40 or 45 days the exposed class begins to decay and then become constant or stable

Fig. 2
figure 2

Behavior of Exposed individuals Ep(t) at various arbitrary order r of the proposed system (2) for h = 0.1,bp = 0.05

Full size image

From Fig. 3, we see that at the given data the number of infected cases are less than that of susceptible and exposed ones. Here up to 50 days the rate of increase is very high and same for all about orders of derivatives, but after that as the transmission of people from place to place decreases the infection is decreasing. The concerned increase and decrease in population are different due to different fractional orders. The infection reached to the peak value on 50th day and then decays to some certain value of 10 infected cases per day, showing stability or convergency.

Fig. 3
figure 3

Behavior of total infected population Ip(t) at various arbitrary order r of the proposed system (2) for h = 0.1,bp = 0.05

Full size image

From Fig. 4, we conclude that initially the asymptotic infected population slightly increases up to 50 days, as the number of infected individuals at this stage is on peak. After that the number decreases gradually or very slowly at different orders of derivative. The number of this class after the 100th day shows stability and became constant. As decrease and increase in this class are very very low or small as compared to other classes, therefore it is known as asymptotic class.

Fig. 4
figure 4

Behavior of asymptotically infectious population Ap(t) at various arbitrary order r of the proposed system (2) for h = 0.1,bp = 0.05

Full size image

Figure 5 shows that the recovery from disease at the beginning is low as the numbers of infected and exposed are very low. But after some protective actions and precautionary measures, the number of recovered population increases. Here, the recovered population is different for different order of derivatives up to 300 days. After that the increases that occur in the recovered case are different at different fractional order.

Fig. 5
figure 5

Behavior of recovered population Rp(t) at various arbitrary order r of the proposed system (2) for h = 0.1,bp = 0.05

Full size image

Figure 6 demonstrates that the number of reservoir population decreases up to 35 days at different fractional order. After that the number of this class increases as compared to other classes with the passage of time. This means that as no precautionary actions are taken in the society more people will be infected but they will be unaware of their infection which will become cause or reservoir for infection in the future. It means that large number of population will be reserved (infection lies in their bodies) but with the passage of time it is also controlled.

Fig. 6
figure 6

Behavior of reservoir population M(t) at various arbitrary order r of the proposed system (2) for h = 0.1,bp = 0.05

Full size image

Furthermore, to check the sensitivity of the fractional order model by changing the values of contact rate bp and step size h in Fig. 7a–f.

Fig. 7
figure 7

Behavior of all populations at various arbitrary order r of the proposed system (2) for h = 0.05,bp = 0.5

Full size image

Another set has been given from Fig. 8a–f by taking h = 0.01,bp = 0.1.

Fig. 8
figure 8

Behavior of all populations at various arbitrary order r of the proposed system (2) for h = 0.01,bp = 0.1

Full size image

Here in Fig. 9, we have compared the reported real data [58] of infected cases in Wuhan city from 4th January 2020 to 8th March 2020 for 67 days as [6,12,19,25,31,38,44,60,80,131,131,259,467,688,776,1776,1460,1739,1984,2101,2590,2827,3233,3892,3697,3151,3387,2653,2984,2473,2022,1820,1998,1506,1278,2051,1772,1891,399,894,397,650,415,518,412,439,441,435,579,206,130,120,143,146,102,46,45,20,31,26,11,18,27,29,39,39]

Fig. 9
figure 9

Comparison between real and simulated data at given different fractional orders for our proposed model

Full size image

We see that the simulated data has close agreement with the plot of the real data. This phenomenon demonstrates the efficiency of our numerical results (Fig. 10).

Next we use the second data for Pakistan as in [59]. The initial population is Sp(0) = 220, Ep(0) = 120,Ip(0) = 1.30, Ap(0) = 0.3, Rp(0) = 1.02, M(0) = 1.055. The parameters values are given in Table 2 [59] (Table 3).

Table 3 Description of the parameters given in model (1) for Pakistan
Full size table
Fig. 10
figure 10

Behavior of all populations at various arbitrary order r of the proposed system (2) for h = 0.01

Full size image

Hence present comparison between real data and simulated data Fig. 11. The confirmed cases in Pakistan per day reported in [60] from the 1 March 2021 to 15th of September 2021 for 200 days as [4,4,5,5,5,5,5,6,15,17,18,19,19,31,51,182,245,331,439,485,629,758,856,962,1034,1171,1139,1454,1554,1836,19972262,2520,2646,2899,3058,3549,3735,3852,3902,4162,4150,4307,4362,4824,5143,5122,5660,6043,6742,7286,7703,8479,8925,9438,10103,10586,11058,11747,11996,12380,12900,13818,14498,14814,15716,16370,17574,18003,20267,21587,22037,23268,25609,26003,26230,27054,27904,29266,30503,31775,32578,34386,34642,36228,37657,38150,38900,39690,40358,40880,42687,44777,47607,50234,53300,56144,59394,63400,57170,60470,75053,78699,83182,79700,84762,85321,89583,93233,97690,100324,104648,105087,106142,107733,107270,107607,107460,107784,106023,106775,106361,108100,108466,103543,95388,95241,95219,94522,91408,90358,89250,87345,86770,84234,77418,77360,73536,60234,57668,53431,53333,52203,51057,50080,40242,29274,29626,27189,26191,25279,24983,24941,24912,24908,24935,24827,20597,19230,18253,17573,17548,17555,17588,17103,16229,16685,16014,16001,13706,13385,12464,11697,11542,10378,10446,9940,9356,8739,8555,8585,8500,8553,8623,8633,8564,8512,8660,8883,6020,6234,6477,6545,5291,5546,5979,5786,5582,5525]

Fig. 11
figure 11

Comparison between real and simulated data at given different fractional orders for our proposed model

Full size image

5 Concluding remarks

Assigning different experimental values taken from [51, 59] to the parameters of (2), we have performed the required simulations as compared to integer order simulation of the system (1). We noticed that by increasing rate of protection, cure and decrease rate of transmission, the minimization and stabling in the numbers of infected individuals can be achieved. By studying such dynamical system, one can know how to control the population from being infected and isolation of infected ones from transmission (immigration). This will be very easy for policy makers and health sector to implement precautionary measures. We can predict for future on the basis of basic reproductive number. From epidemiological point of view it will be very interesting for medical science researchers to know about the history (past), present and future of infection by investigating such type of fractional mathematical model for the pandemic. This model can be applied to the population where social gathering occurs locally or globally. Further, by using fixed point theory the solution of the considered fractional dynamical system has been proved for the existence and uniqueness, while the rate of decaying and growth has been shown through global ways. Hence, fractional calculus can be used for comprehensive explanation of various dynamical models. Further, we observe that increasing precautionary measures will increase the recovered population. The data of Wuhan and Pakistan have been used to demonstrate the model. Also we have compared our simulated data with some reported real data of Wuhan and Pakistan for infected population. We have observed that both simulated and real data plots closely agreed. This shows that the established results are true and applicable. These types of models usually provide interesting indications for future planning and understanding the transmission dynamics of the disease.