System of Time Fractional Models for COVID-19: Modeling, Analysis and Solutions

Abstract

The emergence of the COVID-19 outbreak has caused a pandemic situation in over 210 countries. Controlling the spread of this disease has proven difficult despite several resources employed. Millions of hospitalizations and deaths have been observed, with thousands of cases occurring daily with many measures in place. Due to the complex nature of COVID-19, we proposed a system of time-fractional equations to better understand the transmission of the disease. Non-locality in the model has made fractional differential equations appropriate for modeling. Solving these types of models is computationally demanding. Our proposed generalized compartmental COVID-19 model incorporates effective contact rate, transition rate, quarantine rate, disease-induced death rate, natural death rate, natural recovery rate, and recovery rate of quarantine infected for a holistic study of the coronavirus disease. A detailed analysis of the proposed model is carried out, including the existence and uniqueness of solutions, local and global stability analysis of the disease-free equilibrium (symmetry), and sensitivity analysis. Furthermore, numerical solutions of the proposed model are obtained with the generalized Adam–Bashforth–Moulton method developed for the fractional-order model. Our analysis and solutions profile show that each of these incorporated parameters is very important in controlling the spread of COVID-19. Based on the results with different fractional-order, we observe that there seems to be a third or even fourth wave of the spike in cases of COVID-19, which is currently occurring in many countries.

1. Introduction

In 2019, the COVID-19 outbreak was sudden and rapidly spreading beyond understanding at the time, with more than 213 countries declaring the disease as a pandemic. The coronavirus disease, formally known as COVID-19 or SARS-CoV-2, originated in the city of Wuhan, in the Hubei province of China. The source is still unknown. It is believed that an animal reservoir caused the appearance of the virus in the human population. In recent times, many countries have started an investigation into the exact source of this new coronavirus. COVID-19 has spawned a global pandemic that has substantially affected millions around the world. The coronavirus disease causes an upper-respiratory infection that is easily transmitted between hosts with no known vaccine. The Centers for Disease Control and Prevention have declared this outbreak a public health risk. This disease spreads through droplets caused by coughs, sneezes, or talking within a range of approximately six feet. It is known that social distancing and face masks are effective ways to prevent the spread of the virus between hosts. COVID-19 has affected people all over the world and continues to spread. Throughout history, major pandemics have occurred, such as the Spanish Flu (1981–1920) and the recent Ebola epidemic (2014–2016), that have had a significant impact on human health and global economics. As for recently (March 2021), the coronavirus has infected more than 127.3 million people worldwide, killing over 2.7 million people. Specifically, in the USA, over 31 million cases have already been recorded, with over 540,000 deaths [1].

The importance of modeling and tracking the severity of the virus is not only necessary to prepare for future outbreaks, but to inform governments and citizens of what they should be doing now to mitigate transmission and infection. By modeling and tracking the trend of the current outbreak, we can prepare for outbreaks in the upcoming seasons. Doing this will help reduce the severity of next coronavirus outbreaks. There have been many proposed mathematical models and analyses by a large number of infectious disease researchers on COVID-19 and similar diseases see [2,3,4,5,6,7,8,9,10,11]. Furthermore, some authors have proposed a SIR or related models for COVID-19. In [3,5,6], the Susceptible-Exposed-Infected-Recovered (SEIR) model was used to analyze data generated from Italy and China. Furthermore, in [8], Anastassopouloua et al. used the Susceptible-Infected-Recovered-Dead (SIRD) model to calculate the basic reproduction number, infection, and per day recovery rate for the data generated from China. Recent approaches in dealing with infectious diseases such as COVID-19 using data-driven methods such as machine learning and deep learning, hybrid artificial intelligence, and principal component analysis can be found in [12,13,14,15,16,17]. However, COVID-19 is rare, complex, and many things are yet unknown, which set limitations to what known models (especially the integer-order models) could capture.

Non-integer (fractional) order derivatives can be used to model most phenomena in science, engineering, and social science. This is due to the property of non-locality, which is inherent in many complex structures. Some example applications are seen in biological modeling, financial modeling, and random walk, dynamical system control theory, viscoelasticity, nanotechnology, anomalous transport, and anomalous diffusion [18,19,20,21,22,23,24,25,26,27,28]. It should be noted here that solving fractional differential equations is very difficult due to the memory effect. Beyond numerical methods, which are often computationally demanding, different analytical methods have been proposed to handle linear and nonlinear fractional differential equations, see [19,29,30,31,32,33,34,35,36].

Due to the complex nature of COVID-19, this paper explores the use of fractional calculus to develop mathematical models in order to understand and answer the following pertinent questions: (1) what parameters are essential (most impactful) for control measure purposes; (2) dynamics of the spread and long time impact of the disease in different population densities; and (3) what to expect at the resurgence of the outbreak and necessary control measures to reduce the impact. This work is focused on advancing research to understand better factors that could help curtailing, or possible control measures for eradication for a long time benefit. Different control measures are investigated using effective, efficient, and more general mathematical models to reveal the intrinsic nature of the COVID-19. Modeling, analysis, and numerical solutions are the focus of the current paper on the disease transmission and spread that includes understanding the delay to ’flattening the curve’ as experienced in many cities and countries.

Specifically, our proposed generalized compartmental COVID-19 model takes into account many key parameters to gain insight into the complexity of the disease. The parameters considered are effective contact rate, transition rate (from exposed quarantine and recovered to susceptible and infected quarantined individuals), quarantine rate, disease-induced death rate, natural death rate, natural recovery rate, and recovery rate of quarantine infected. We simulate and discuss the parameters on the dynamics of the solution profiles. In this paper, a detailed analysis of the proposed model is carried out, including the existence and uniqueness of solutions, local and global stability analysis, endemic equilibrium analysis, and sensitivity analysis. In addition, numerical solutions of the proposed fractional SEQIR model are obtained using the generalized Adam–Bashforth–Moulton method developed for the fractional-order model.

The rest of the paper is structured as follows: In Section 2, preliminaries are presented to include basic definitions, notations used in this present investigation, and valuable results from the literature. The generalized mathematical model is formulated in Section 3, which accounts for the loss of immunity of the infected recovered individuals. Analysis of the model to include the existence and uniqueness of solutions, local and global stability analysis, endemic equilibrium (symmetry) analysis, and sensitivity analysis is presented with detailed proofs in Section 4. The numerical simulation of the proposed model is presented in Section 5. Finally, Section 6 is devoted to summary and recommendations.

2. Preliminaries

In this section, we present some definitions, notations, and some known results needed in sequel. Caputo’s fractional derivative is adopted in this work.

Definition 1.

A real valued function u is said to be in the space $C_{\eta }$, $\eta \in \mathbb{R}$, $t>0$, if there exists a real number p with $p>\eta$ such that

$$u\left(t\right)=t^{p}g\left(t\right),$$

where $g\in C[0,\infty )$, and it is said to be in the space $C_{\eta }^m$ iff $u^{\left(m\right)}\in C_{\eta }$, $m\in \mathbb{N}$.

Definition 2.

The Laplace transform is defined by

$$\mathcal{N}\left(s\right)=\mathcal{L}\left(N\left(t\right)\right)\left(s\right)={\int }_0^{\infty }{e}^{-st}N\left(t\right)dt,$$

(1)

where $N\left(t\right)$ is an n-dimensional vector-valued function.

Definition 3.

The Riemann–Liouville’s $\left(RL\right)$ fractional integral operator of order $\alpha \ge 0$, of a function $u\in L^1(a,b)$ is given as

$$J^{\alpha }u\left(t\right)=\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^t{(t-\tau )}^{\alpha -1}u\left(\tau \right)d\tau ,t>0,\alpha >0,$$

(2)

where Γ is the Gamma function, and $J^{0}u\left(t\right)=f\left(t\right)$.

Definition 4.

The fractional derivative in the Caputo’s sense is defined as [37],

$${\mathcal{D}}^{\alpha }u\left(t\right)=J^{n-\alpha }{D}^{n}u\left(t\right)=\frac{1}{\mathsf{\Gamma }(n-\alpha )}{\int }_0^t{(t-\tau )}^{n-\alpha -1}{u}^{\left(n\right)}\left(\tau \right)d\tau ,$$

(3)

where $n-1<\alpha \le n$, $n\in \mathbb{N}$, $t>0$.

Lemma 1.

Let $t\in (a,b]$. Then

$$\left[J_a^{\alpha }{(t-a)}^{\beta }\right]\left(t\right)=\frac{\mathsf{\Gamma }(\beta +1)}{\mathsf{\Gamma }(\beta +\alpha +1)}{(t-a)}^{\beta +\alpha },\alpha ⩾0,\beta >0.$$

(4)

Definition 5.

The generalized Mittag–Leffler function (also called two-parameter Mittag–Leffler function) is defined as:

$$E_{\alpha ,\beta }\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\mathsf{\Gamma }(\alpha k+\beta )},\beta ,z\in \mathbb{C},\mathrm{Re}\left(\alpha \right)>0.$$

(5)

When $\beta =1$, a special case of this function is obtained, called a one-parameter Mittag–Leffler function, given as:

$$E_{\alpha }\left(z\right)=E_{\alpha ,1}\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\mathsf{\Gamma }(\alpha k+1)},z\in \mathbb{C},\mathrm{Re}\left(\alpha \right)>0.$$

(6)

Taking different values for $\alpha$ and $\beta$, we see that this function represents many important special cases, such as:

$$\begin{array}{ccc}\hfill E_{1,1}\left(z\right)& =& e^z,E_{1,2}\left(z\right)=\frac{e^z-1}{z},\hfill \\ \hfill E_{2,1}\left(z\right)& =& cosh\left(\sqrt{z}\right),E_{2,2}\left(z\right)=\frac{sinh\left(\sqrt{z}\right)}{\sqrt{z}},\hfill \\ \hfill E_{2,1}(-z^2)& =& cos\left(z\right),E_{0,1}\left(z\right)=\frac{1}{1-z},\left|z\right|<1.\hfill \end{array}$$

The generalized Mittag–Leffler function is an entire function [38,39], so bounded in any finite interval.

Lemma 2.

[40] Let $0<\alpha \le 1$, and $\lambda >0$. Then, the function $E_{\alpha ,\beta }(-\lambda z^{\alpha })$, $\alpha \le \beta$ is positive monotonically decreasing functions of $z>0$.

Definition 6.

Let $n-1<\alpha \le n$. Then, the Laplace transform of the Caputo’s fractional derivative is given as

$$\mathcal{L}\left\{{\mathcal{D}}^{\alpha }N\left(t\right)\right\}\left(s\right)=s^{\alpha }\mathcal{N}\left(s\right)+\sum _{k=0}^{n-1}{s}^{n-k-1}{N}^{\left(k\right)}\left(0_{+}\right),Re\left(s\right)>0,$$

(7)

where $Re\left(\alpha \right)>0$.

Lemma 3.

Considering the two-parameter Mittag–Leffler function, the Laplace transform formula is given as

$$\mathcal{L}\left\{t^{\beta -1}{E}_{\alpha ,\beta }\left({\lambda }^{*}{t}^{\alpha }\right)\right\}\left(s\right)=\frac{s^{\alpha -\beta }}{s^{\alpha }-{\lambda }^{*}},Re\left(z\right)>0\left|{\lambda }^{*}{s}^{-\alpha }\right|<1,$$

(8)

where $\alpha ,\beta ,{\lambda }^{*}\in \mathbb{C}$, $Re\left(\alpha \right)>0$, $Re\left(\beta \right)>0$.

Lemma 4.

[41] For $0<\alpha \le 1$, let $\phi \in C[a,b]$ and ${\mathcal{D}}^{\alpha }\in (a,b]$. Then

$$u\left(t\right)=u\left(a\right)+\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\mathcal{D}}^{\alpha }u\left(\eta \right){(t-a)}^{\alpha },0\le \eta \le t,\forall t\in (a,b].$$

(9)

Lemma 4 is called the generalized mean value theorem.

3. Model Formulation

In this section, we develop a compartmental COVID-19 model where individuals are classified at time t as: susceptible (S), exposed (E), exposed but quarantined ($Q_1$), infected (I), infected but quarantined ($Q_2$), and recovered (R). The infected compartment, $I,$ contains both asymptomatic and symptomatic carriers who have not been clinically confirmed to be COVID-19 positive. Transmission of COVID-19 occurs as a result of touching a contaminated surface or sharing very close space (less than 6 ft/2m) with an infectious host coughing or sneezing [42]. We assume that exposed and infected individuals are able to transmit the disease to the general public, but quarantined individuals, either exposed or infected, are unable to transmit the virus to general public because they stay in isolation. Susceptible individuals contract the disease at a rate $\beta .$ Once individuals are confirmed positive, they are placed in self-quarantine (or hospitalize) at rate $k_4.$ We incorporate the parameters $k_1,k2,k_3$ for quarantine rate of exposed individuals, progression rate from exposure to infected class, and progression rate from exposed quarantined to infected quarantined class, respectively. Individuals in I and $Q_2$ decrease by diseased-induced death rate $\delta .$ There are several reports that recovered individuals are susceptible to re-infection. On 24 April 2020, the World Health Organization (WHO) reported that there is currently no proof that individuals who have recovered from COVID-19 and have tested negative for the virus are immune to re-infection [1]. Because of this, we take into consideration loss of immunity at rate ${\rho }_2.$ See Figure 1 for an illustration of the model and

Table 1. Parameters of the disease model and their meanings.

Parameter	Likely Range (Sources)	Default Value
b (Recruitment rate of the population)	$\mu \times N\left(0\right)$	78,680
$\beta$ (Effective contact rate)	$2.1011\times {10}^{-8}-9.11\times {10}^{-8}$ [5,43]	$4.1011\times {10}^{-8}$ day${}^{-1}$
${\rho }_1$(Transition rate from exposed quarantined to susceptible individuals)	1/14 [43]	1/14
${\rho }_2$(Transition rate from recovered to susceptible individuals)	Assumed	5/100
$k_1$(Quarantine rate of exposed individuals)	1/10 [5]	1/10 day${}^{-1}$
$k_2$ (Transition rate from exposed to infected class $I\left(t\right)$)	1/14–1/3 [43,44,45]	1/7 day${}^{-1}$
$k_3$ (Transition from exposed quarantined to infected quarantined)	0.1259 [43]	0.1259
$k_4$ (Quarantining rate of individuals in the $I\left(t\right)$ class)	0.2– 1 [43,46]	0.3654 day${}^{-1}$
$\delta$(Disease-induced death rate)	$1.7826\times {10}^{-5}$ [5,43]	$1.7826\times {10}^{-5}$
$\mu$ (Natural death rate)	7.1/1000	7.1/1000
${\sigma }_1$ (Natural recovery rate)	1/30–1/3 day${}^{-1}$ [43,47]	1/20 day${}^{-1}$
${\sigma }_2$ (Recovery rate of quarantine infected)	0.11624 [43,48]	0.11624 day${}^{-1}$
$S\left(0\right)$ (Initial value of the susceptible)	[5]	11,081,000
$E\left(0\right)$ (Initial value of the expose)	[5]	399
$Q_1\left(0\right)$(Initial value of the expose, quarantined)	[5]	159
$I\left(0\right)$ (Initial value of the infected)	[5]	54
$Q_2\left(0\right)$ (Initial value of the infected, quarantined)	[5]	28
$R\left(0\right)$ (Initial value of the recovered)	[5]	12

for the model parameters.

A nonlinear system under the above transmission assumptions is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{d^{\alpha }S}{d{t}^{\alpha }}& =b-\beta (E+I)S+{\rho }_1{Q}_1+{\rho }_{2}R-\mu S\hfill \\ \hfill \frac{d^{\alpha }E}{d{t}^{\alpha }}& =\beta (E+I)S-(k_1+k_2+\mu )E\hfill \\ \hfill \frac{d^{\alpha }{Q}_1}{d{t}^{\alpha }}& =k_{1}E-({\rho }_1+k_3+\mu )Q_1\hfill \\ \hfill \frac{d^{\alpha }I}{d{t}^{\alpha }}& =k_{2}E-(k_4+{\sigma }_1+\delta +\mu )I\hfill \\ \hfill \frac{d^{\alpha }{Q}_2}{d{t}^{\alpha }}& =k_3{Q}_1+k_{4}I-({\sigma }_2+\delta +\mu )Q_2\hfill \\ \hfill \frac{d^{\alpha }R}{d{t}^{\alpha }}& ={\sigma }_{1}I+{\sigma }_2{Q}_2-({\rho }_2+\mu )R\hfill \end{array}\end{array}$$

(10)

and initial points to be

$$\begin{array}{cc}\hfill S\left(0\right)& =S_0\hfill \\ \hfill E\left(0\right)& =E_0\hfill \\ \hfill Q_1\left(0\right)& =Q_{1,0}\hfill \\ \hfill I\left(0\right)& =I_0\hfill \\ \hfill Q_2\left(0\right)& =Q_{2,0}\hfill \\ \hfill R\left(0\right)& =R_0\hfill \end{array}$$

(11)

with $0<\alpha \le 1$, $\frac{d^{\alpha }}{d{t}^{\alpha }}$ the Caputo fractional derivative of order $\alpha$. The parameters $b,\mu ,{\sigma }_1,{\sigma }_2,{\rho }_1$ are birth, per capital death rates of individuals, natural recovery rate, recovery rate of quarantine infected, and transition rate from recovered to susceptible individuals, respectively. All the parameters are defined and values provided in Table 1 below. Note that $N=S+E+Q_1+I+Q_2+R.$ This means that

$$\frac{d^{\alpha }N}{d{t}^{\alpha }}=b-\delta (I+Q_2)-\mu N$$

4. Analysis of the Time-Fractional COVID-19 Model

Our focus in this Section is on detailed analysis of the proposed COVID-19 model of time-fractional type. The existence, uniqueness, and non-negativity results of solutions to the system (10)–(11) are presented and proven. We further investigate conditions and results on existence, stability (local and global), and equilibrium of the proposed model.

Let us define a feasible region $\mathsf{\Delta }$ of system (10)–(11) as:

$$\nabla =\left\{x=(S,E,Q_1,I,Q_2,R)\in {\mathbb{R}}_{+}^6:S\ge 0,N\le \frac{b}{\mu }\right\}.$$

The region is an area which contains all biologically and mathematically relevant solutions of system (10)–(11).

We denote the following parameter for easy calculation.

$$\begin{array}{ccc}\hfill {\gamma }_1& =& k_1+k_2+\mu \hfill \\ \hfill {\gamma }_2& =& {\rho }_1+k_3+\mu \hfill \\ \hfill {\gamma }_3& =& {\sigma }_2+\delta +\mu \hfill \\ \hfill {\gamma }_4& =& k_4+{\sigma }_1+\delta +\mu \hfill \\ \hfill {\gamma }_5& =& k_2+k_4+{\sigma }_1+\delta +\mu .\hfill \end{array}$$

4.1. Results on the Existence and Uniqueness Of Solutions

The following Lemma is needed to prove the existence and uniqueness of solutions to the system considered in (10)–(11). Our results include the forward-invariant property of the domain considered.

Lemma 5.

For $0<\alpha \le 1$, let $u\in C[0,b]$ and ${\mathcal{D}}^{\alpha }u\in (0,b]$. Then

(a)

the function u is non-decreasing if ${\mathcal{D}}^{\alpha }u\left(t\right)\ge 0$, $\forall t\in (0,b)$.

(b)

the function u is non-increasing if ${\mathcal{D}}^{\alpha }u\left(t\right)\le 0$, $\forall t\in (0,b)$.

Proof. 

The proof is a direct consequence of Lemma 4. □

Theorem 1.

The initial valued problem for the COVID-19 model of fractional-order type in (10) and (11) has a unique solution in ∇.

Proof. 

Using Theorem (3.1) and Remark (3.2) in [49], the existence and uniqueness of the solution in $(0,\infty )$ is obtained. Furthermore, we need to show that the closed set ∇ is forward-invariant w.r.t (10). To prove this, we first consider the fractional differential:

$$\begin{array}{c}\hfill \frac{d^{\alpha }N}{d{t}^{\alpha }}=b-\mu N.\end{array}$$

Using the Laplace transform given in 6 in both sides, we obtain

$$s^{\alpha }\mathcal{N}\left(s\right)-s^{\alpha -1}N\left(0\right)=\frac{b}{s}-\mu \mathcal{N}\left(s\right).$$

(12)

Re-arranging (12), we get

$$\mathcal{N}\left(s\right)=b\frac{1}{s(s^{\alpha }+\mu )}+N\left(0\right)\frac{s^{\alpha -1}}{(s^{\alpha }+\mu )}.$$

(13)

Applying Laplace inverse transform, and Lemma 3 after some simplification, we obtain the solution as

$$N\left(t\right)=N\left(0\right)E_{\alpha }(-\mu t^{\alpha })+\frac{b}{\mu }\left[1-E_{\alpha }(-\mu t^{\alpha })\right].$$

(14)

Recall that

$$\begin{array}{c}\hfill \frac{d^{\alpha }N}{d{t}^{\alpha }}=b-\delta (I+Q_2)-\mu N\le b-\mu N.\end{array}$$

Using (14) combined with Lemma 2, it follows that

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim}supN\left(t\right)\le \frac{b}{\mu }.\end{array}$$

Thus, ∇ is positive and bounded. So, all solutions of the system (10) and (11) are confined in the set ∇. □

4.2. Existence, Stability, and Equilibrium Results for COVID-19 Model of Fractional Type

From above, the region ∇ is positively-invariant under models (10) and (11). Then the system is both epidemiologically and mathematically well-posed. Therefore, it is sufficient to study the dynamics of the model in ∇.

In the absence of COVID-19, we have $E=Q_1=I=Q_2=0$ and at equilibrium

$$\begin{array}{c}\hfill S^{*}=N^{*}\to \frac{b}{\mu }.\end{array}$$

There exists a disease-free equilibrium of system (10) given by

$$\begin{array}{c}\hfill DFE=\left(\frac{b}{\mu },0,0,0,0,0\right).\end{array}$$

(15)

4.2.1. Computation of ${\mathcal{R}}_0$

The basic reproduction number of infectious disease models, denoted by ${\mathcal{R}}_0$, is a critical epidemiological quantity of public health concerns. It measures the intensity of an outbreak of diseases and plays a significant role in evaluating control strategies. We will use a method called next generation described in [50,51,52] to compute ${\mathcal{R}}_0$ of our model. This method is defined as the largest eigenvalue of the matrix $F{V}^{-1}$, where the matrices F and $V^{-1}$ are determined as follows:

$$F_i=\left[\begin{array}{c}{F}_E\\ F_{Q_1}\\ F_I\\ F_{Q_2}\end{array}\right]=\left[\begin{array}{cc}& \beta (E+I)S\\ & 0\\ & 0\\ & 0\end{array}\right]$$

and

$$V_i=\left[\begin{array}{c}{V}_E\\ V_{Q_1}\\ V_I\\ V_{Q_2}\end{array}\right]=\left[\begin{array}{cc}& {\gamma }_{1}E\\ & -k_{1}E+{\gamma }_2{Q}_1\\ & -k_{2}E+{\gamma }_{4}I\\ & -k_3{Q}_1-k_{4}I+{\gamma }_3{Q}_2\end{array}\right].$$

Therefore

$$F=\left[\begin{array}{cccc}\beta S^{*}& 0& \beta S^{*}& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right].$$

and

$$V=\left[\begin{array}{cccc}{\gamma }_1& 0& 0& 0\\ -k_1& {\gamma }_2& 0& 0\\ -k_2& 0& {\gamma }_4& 0\\ 0& -k_3& -k_4& {\gamma }_3\end{array}\right],$$

where $S^{*}=\frac{b}{\mu }.$

The basic reproduction number (${\mathcal{R}}_0$) is the spectral radius of $\left|F{V}^{-1}-\lambda I\right|$ and is given by

$$\begin{array}{c}\hfill {\mathcal{R}}_0=\frac{b\beta {\gamma }_5}{\mu {\gamma }_1{\gamma }_4}\end{array}$$

(16)

4.2.2. Stability Results

Theorem 2

(Local Stability). The disease-free equilibrium $\left(\frac{b}{\mu },0,0,0,0,0\right)$ of system (10) is locally asymptotically stable if ${\mathcal{R}}_0<1$, and unstable if ${\mathcal{R}}_0>1$.

Proof. 

To determine the local stability of the disease-free equilibrium, we use the eigenvalues of the associated Jacobian matrix. The Jacobian matrix of system (3) at

$$DFE=\left(\frac{b}{\mu },0,0,0,0,0\right)$$

is

$$J=\left(\begin{array}{cccccc}-\mu & -\frac{b\beta }{\mu }& {\rho }_1& -\frac{b\beta }{\mu }& 0& {\rho }_2\\ 0& \frac{b\beta }{\mu }-{\gamma }_1& 0& \frac{b\beta }{\mu }& 0& 0\\ 0& k_1& -{\gamma }_2& 0& 0& 0\\ 0& k_2& 0& -{\gamma }_4& 0& 0\\ 0& 0& k_3& k_4& -{\gamma }_3& 0\\ 0& 0& 0& {\sigma }_1& {\sigma }_2& -({\rho }_2+\mu )\end{array}\right).$$

The characteristic equation of the Jacobian matrix is

$$\begin{array}{c}\hfill \left(\lambda +\mu \right)\left(\lambda +{\rho }_2+\mu \right)\left(\lambda +{\gamma }_2\right)\left(\lambda +{\gamma }_3\right)\left({\lambda }^2+a_1\lambda +a_0\right)=0,\end{array}$$

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill a_1& =-{\gamma }_1\left({\mathcal{R}}_0-1\right)+{\gamma }_4+\frac{b\beta k_2}{\mu {\gamma }_4}\hfill \\ \hfill a_0& =-{\gamma }_1{\gamma }_4\left({\mathcal{R}}_0-1\right).\hfill \end{array}\end{array}$$

For ${\mathcal{R}}_0<1,$ all the roots of the characteristic equation have negative real parts. □

The public health implication of Theorem 2 is that the spread of the virus can be reduced or controlled when ${\mathcal{R}}_0<1,$ if the initial sizes of the sub-populations of the model are in the basin of attraction of the DFE (∇). Global stability of the $DFE$ is necessary to ensure that the disease elimination is independent of the initial sizes of the sub-populations [53].

Theorem 3

(Global Stability). The disease-free equilibrium $\left(\frac{b}{\mu },0,0,0,0,0\right)$ of system 10 is globally asymptotically stable if ${\mathcal{R}}_0<1$.

Proof. 

Consider the Lyapunov function $V=(S,E,Q_1,I,Q_2,R):{\mathbb{R}}_{+}^6$ defined as

$$V=E+\frac{\beta b}{\mu {\gamma }_4}I.$$

Then

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{d^{\alpha }V}{d{t}^{\alpha }}=& \frac{d^{\alpha }E}{d{t}^{\alpha }}+\frac{\beta b}{\mu {\gamma }_4}\frac{d^{\alpha }I}{d{t}^{\alpha }}\hfill \\ \hfill =& \beta (E+I)S-{\gamma }_{1}E+\frac{\beta b}{\mu {\gamma }_4}\left[k_{2}E-{\gamma }_{4}I\right]\hfill \\ \hfill =& \beta ES-{\gamma }_{1}E+\frac{\beta b}{\mu {\gamma }_4}{k}_{2}E+\beta IS-\frac{\beta b}{\mu }I\hfill \\ \hfill \le & \beta E{S}^{*}-{\gamma }_{1}E+\frac{\beta b}{\mu {\gamma }_4}{k}_{2}E+\beta I{S}^{*}-\frac{\beta b}{\mu }I\hfill \\ \hfill =& \left(\frac{\beta b}{\mu }-{\gamma }_1\right)E+\frac{\beta b{k}_{2}E}{\mu {\gamma }_4}\hfill \\ \hfill =& {\gamma }_1\left(\frac{\beta b{\gamma }_4}{\mu {\gamma }_1{\gamma }_4}-1\right)E+\frac{\beta b{k}_{2}E}{\mu {\gamma }_4}\hfill \\ \hfill =& {\gamma }_1\left(\frac{\beta b{\gamma }_5}{\mu {\gamma }_1{\gamma }_4}-1\right)E-\frac{\beta b{k}_{2}E}{\mu {\gamma }_4}+\frac{\beta b{k}_{2}E}{\mu {\gamma }_4}\hfill \\ \hfill =& {\gamma }_1\left({\mathcal{R}}_0-1\right)E.\hfill \end{array}\end{array}$$

We have ${\displaystyle \frac{d^{\alpha }V}{d{t}^{\alpha }}=0}$ when $E=0,$ and ${\displaystyle \frac{d^{\alpha }V}{d{t}^{\alpha }}<0}$ when $E>0$ provided that ${\mathcal{R}}_0<1.$ It follows from the the generalized Lasalle invariance principle [54] that the disease-free equilibrium is globally asymptotically stable. □

4.2.3. Computation of an Endemic Equilibrium

Next, we summarize the predictions of endemic equilibrium of model 10 in the following Lemma.

Lemma 6 (Equilibrium).

Let

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill A_1& =\frac{\left({\rho }_1(k_2+\mu )+{\gamma }_1{\gamma }_4(k_3+\mu )\right)}{k_2{\gamma }_2}>0\hfill \\ \hfill A_2& ={\sigma }_1+\frac{{\sigma }_2{k}_3{k}_1{\gamma }_4}{k_2{\gamma }_2{\gamma }_3}+\frac{{\sigma }_2{k}_4}{{\gamma }_3}>0.\hfill \end{array}\end{array}$$

(a)

If ${\mathcal{R}}_0>1$ and

(i)

${\rho }_2=0,$ then the model (10) has a unique endemic equilibrium.

(ii)

${\displaystyle \frac{{\rho }_2}{{\rho }_2+\mu }<\frac{A_1}{A_2},}$ then the model (10) has a unique endemic equilibrium.

(b)

If ${\mathcal{R}}_0<1$ and ${\displaystyle \frac{{\rho }_2}{{\rho }_2+\mu }>\frac{A_1}{A_2},}$ then the model (10) has a unique endemic equilibrium.

Proof. 

The endemic equilibrium is the solution of

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill b-\beta (E^{**}+I^{**})S^{**}+{\rho }_1{Q}_1^{**}+{\rho }_2{R}^{**}-\mu S& =0\hfill \\ \hfill \beta (E^{**}+I^{**})S^{**}-{\gamma }_1{E}^{**}& =0\hfill \\ \hfill k_1{E}^{**}-{\gamma }_2{Q}_1^{**}& =0\hfill \\ \hfill k_2{E}^{**}-{\gamma }_4{I}^{**}& =0\hfill \\ \hfill k_3{Q}_1^{**}+k_4{I}^{**}-{\gamma }_3{Q}_2^{**}& =0\hfill \\ \hfill {\sigma }_1{I}^{**}+{\sigma }_2{Q}_2^{**}-({\rho }_2+\mu )R^{**}& =0.\hfill \end{array}\end{array}$$

(17)

It follows that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill E^{**}& =\frac{{\gamma }_4}{k_2}{I}^{**}\hfill \\ \hfill Q_1^{**}& =\frac{k_1}{{\gamma }_2}{E}^{**}=\frac{k_1{\gamma }_4}{k_2{\gamma }_2}{I}^{**}\hfill \\ \hfill Q_2^{**}& =\frac{k_3{Q}_1^{**}+k_4{I}^{**}}{{\gamma }_3}=\frac{k_3{k}_1{\gamma }_4+k_4{k}_2{\gamma }_2}{k_2{\gamma }_2{\gamma }_3}{I}^{**}\hfill \\ \hfill R^{**}& =\frac{{\sigma }_1{I}^{**}+{\sigma }_2{Q}_2^{**}}{({\rho }_2+\mu )}=\frac{{\sigma }_1{k}_2{\gamma }_2{\gamma }_3+{\sigma }_2\left(k_3{k}_1{\gamma }_4+k_4{k}_2{\gamma }_2\right)}{k_2{\gamma }_2{\gamma }_3({\rho }_2+\mu )}{I}^{**}\hfill \\ \hfill S^{**}& =\frac{{\gamma }_1{E}^{**}}{\beta (E^{**}+I^{**})}\hfill \end{array}\end{array}$$

(18)

Substitute the last line of Equation (18) into the first line of Equation (17) to obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & b-\beta (E^{**}+I^{**})\frac{{\gamma }_1{E}^{**}}{\beta (E^{**}+I^{**})}+{\rho }_1{Q}_1^{**}+{\rho }_2{R}^{**}-\frac{\mu {\gamma }_1{E}^{**}}{\beta (E^{**}+I^{**})}=0\hfill \\ \hfill & b-{\gamma }_1{E}^{**}+{\rho }_1{Q}_1^{**}+{\rho }_2{R}^{**}-\frac{\mu {\gamma }_1{E}^{**}}{\beta (E^{**}+I^{**})}=0\hfill \\ \hfill & b\beta (E^{**}+I^{**})-{\gamma }_1{E}^{**}\beta (E^{**}+I^{**})+{\rho }_1{Q}_1^{**}\beta (E^{**}+I^{**})+{\rho }_2{R}^{**}\beta (E^{**}+I^{**})\hfill \\ \hfill & -\mu {\gamma }_1{E}^{**}=0\hfill \end{array}\end{array}$$

(19)

Let

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill B_1& =\frac{{\gamma }_4}{k_2}\hfill \\ \hfill B_2& =\frac{k_1{\gamma }_4}{k_2{\gamma }_2}=\frac{k_1{B}_1}{{\gamma }_2}\hfill \\ \hfill B_3& =\frac{{\sigma }_1{k}_2{\gamma }_2{\gamma }_3+{\sigma }_2\left(k_3{k}_1{\gamma }_4+k_4{k}_2{\gamma }_2\right)}{k_2{\gamma }_2{\gamma }_3({\rho }_2+\mu )}\hfill \\ \hfill & =\frac{1}{{\rho }_2+\mu }\left({\sigma }_1+\frac{{\sigma }_2{k}_3{k}_1{B}_1}{{\gamma }_2{\gamma }_3}+\frac{{\sigma }_2{k}_4}{{\gamma }_3}\right).\hfill \end{array}\end{array}$$

(20)

Then

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & b\beta (B_1{I}^{**}+I^{**})-{\gamma }_1{B}_1{I}^{**}\beta (B_1{I}^{**}+I^{**})+{\rho }_1{B}_2{I}^{**}\beta (B_1{I}^{**}+I^{**})+\hfill \\ \hfill & {\rho }_2{B}_3{I}^{**}\beta (B_1{I}^{**}+I^{**})-\mu {\gamma }_1{B}_1{I}^{**}=0\hfill \end{array}\end{array}$$

(21)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \left[b\beta (B_1+1)-{\gamma }_1{B}_1\beta (B_1+1)I^{**}+{\rho }_1{B}_2\beta (B_1+1)I^{**}+{\rho }_2{B}_3\beta (B_1+1)I^{**}-\mu {\gamma }_1{B}_1\right]I^{**}=0\end{array}\end{array}$$

(22)

For $I^{**}\ne 0,$ we get

$$\begin{array}{c}\hfill \begin{array}{c}\hfill b\beta (B_1+1)-\mu {\gamma }_1{B}_1-\left[{\gamma }_1{B}_1\beta (B_1+1)-{\rho }_1{B}_2\beta (B_1+1)-{\rho }_2{B}_3\beta (B_1+1)\right]I^{**}=0,\end{array}\end{array}$$

which leads to

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill I^{**}& =\frac{b\beta (B_1+1)-\mu {\gamma }_1{B}_1}{{\gamma }_1{B}_1\beta (B_1+1)-{\rho }_1{B}_2\beta (B_1+1)-{\rho }_2{B}_3\beta (B_1+1)}\hfill \\ \hfill & =\frac{b\beta +B_1\left(b\beta -\mu {\gamma }_1\right)}{\beta (B_1+1)\left[{\gamma }_1{B}_1-{\rho }_1{B}_2-{\rho }_2{B}_3\right]}\hfill \\ \hfill & =\frac{\mu {\gamma }_1{B}_1\left({\mathcal{R}}_0-1\right)}{\beta (B_1+1)\left[{\gamma }_1{B}_1-{\rho }_1{B}_2-{\rho }_2{B}_3\right]}\hfill \\ \hfill & =\frac{\mu {\gamma }_1{B}_1\left({\mathcal{R}}_0-1\right)}{\beta (B_1+1)\left[{\gamma }_1{B}_1-\frac{{\rho }_1{k}_1{B}_1}{{\gamma }_2}-\frac{{\rho }_2}{{\rho }_2+\mu }\left({\sigma }_1+\frac{{\sigma }_2{k}_3{k}_1{B}_1}{{\gamma }_2{\gamma }_3}+\frac{{\sigma }_2{k}_4}{{\gamma }_3}\right)\right]}\hfill \\ \hfill & =\frac{\mu {\gamma }_1{B}_1\left({\mathcal{R}}_0-1\right)}{\beta (B_1+1)\left[\frac{{\rho }_1(k_2+\mu )+(k_3+\mu ){\gamma }_1}{{\gamma }_2}{B}_1-\frac{{\rho }_2}{{\rho }_2+\mu }\left({\sigma }_1+\frac{{\sigma }_2{k}_3{k}_1{B}_1}{{\gamma }_2{\gamma }_3}+\frac{{\sigma }_2{k}_4}{{\gamma }_3}\right)\right]}.\hfill \end{array}\end{array}$$

Further, let

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill A_1& =\frac{{\rho }_1(k_2+\mu )+{\gamma }_1(k_3+\mu )}{{\gamma }_2}{B}_1=\frac{\left({\rho }_1(k_2+\mu )+{\gamma }_1(k_3+\mu )\right){\gamma }_4}{k_2{\gamma }_2}\hfill \\ \hfill A_2& ={\sigma }_1+\frac{{\sigma }_2{k}_3{k}_1{B}_1}{{\gamma }_2{\gamma }_3}+\frac{{\sigma }_2{k}_4}{{\gamma }_3}={\sigma }_1+\frac{{\sigma }_2{k}_3{k}_1{\gamma }_4}{k_2{\gamma }_2{\gamma }_3}+\frac{{\sigma }_2{k}_4}{{\gamma }_3}.\hfill \end{array}\end{array}$$

Then

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill I^{**}& =\frac{\mu {\gamma }_1{B}_1\left({\mathcal{R}}_0-1\right)}{A_2\beta (B_1+1)\left(\frac{A_1}{A_2}-\frac{{\rho }_2}{{\rho }_2+\mu }\right)}.\hfill \end{array}\end{array}$$

(23)

It follows from (23) that $I^{**}>0$ if

• ${\mathcal{R}}_0>1$ and ${\rho }_2=0.$
• ${\mathcal{R}}_0>1$ and ${\displaystyle \frac{{\rho }_2}{{\rho }_2+\mu }<\frac{A_1}{A_2}.}$
• ${\mathcal{R}}_0<1$ and ${\displaystyle \frac{{\rho }_2}{{\rho }_2+\mu }>\frac{A_1}{A_2}.}$

□

4.2.4. Sensitivity Analysis

In this subsection, we perform sensitivity analysis to determine which model parameters have the most significant impact on the threshold ${\mathcal{R}}_0$ and the spread of the virus. When designing control strategies, sensitivity analysis can help identify parameters that should be controlled. We compute a sensitive index for the model parameters; the index predicts the relative change in ${\mathcal{R}}_0$ for the relative change in the parameters [55,56,57]. The parameters $b,\beta$ have a positive impact on the ${\mathcal{R}}_0,$ meaning that an increase in these parameters implies an increase in ${\mathcal{R}}_0$, while the $k_1,k_4,{\sigma }_1,\theta ,\mu$ have negative impact on (decrease) the ${\mathcal{R}}_0.$

Definition 7.

The normalized forward sensitivity index of a variable, L, that depends differentially on a parameter, y, is defined as:

$${\xi }_y=\frac{y}{L}\frac{\partial L}{\partial y}.$$

The sensitivity indices with respect to the model parameters are as follows.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\xi }_b& =\frac{b}{{\mathcal{R}}_0}\frac{\partial {\mathcal{R}}_0}{\partial b}=1\hfill \\ \hfill {\xi }_{\beta }& =\frac{\beta }{{\mathcal{R}}_0}\frac{\partial {\mathcal{R}}_0}{\partial \beta }=1\hfill \\ \hfill {\xi }_{k_1}& =\frac{k_1}{{\mathcal{R}}_0}\frac{\partial {\mathcal{R}}_0}{\partial k_1}=-\frac{k_1}{{\gamma }_1}\hfill \\ \hfill {\xi }_{k_2}& =\frac{k_2}{{\mathcal{R}}_0}\frac{\partial {\mathcal{R}}_0}{\partial k_2}=\frac{k_2\left(k_1-k_4-{\sigma }_1-\delta \right)}{{\gamma }_5{\gamma }_1}\hfill \\ \hfill {\xi }_{k_4}& =\frac{k_4}{{\mathcal{R}}_0}\frac{\partial {\mathcal{R}}_0}{\partial k_4}=-\frac{k_2{k}_4}{{\gamma }_4{\gamma }_5}\hfill \\ \hfill {\xi }_{{\sigma }_1}& =\frac{{\sigma }_1}{{\mathcal{R}}_0}\frac{\partial {\mathcal{R}}_0}{\partial {\sigma }_1}=-\frac{k_2{\sigma }_1}{{\gamma }_4{\gamma }_5}\hfill \\ \hfill {\xi }_{\delta }& =\frac{\delta }{{\mathcal{R}}_0}\frac{\partial {\mathcal{R}}_0}{\partial \delta }=-\frac{k_2\delta }{{\gamma }_4{\gamma }_5}\hfill \\ \hfill {\xi }_{\mu }& =\frac{\mu }{{\mathcal{R}}_0}\frac{\partial {\mathcal{R}}_0}{\partial \mu }=-\frac{k_2\left(\mu {\gamma }_1+\left({\gamma }_1+\mu \right){\gamma }_4\right)+{\gamma }_4^2\left({\gamma }_1+k_4+\mu \right)}{{\gamma }_1{\gamma }_4{\gamma }_5}.\hfill \end{array}\end{array}$$

(24)

It is important to note that a sensitive index can be a constant or depend on other parameters of the model. An ${\xi }_y=\pm 1$ implies that increasing or decreasing y by a given percentage always increases (decreases) ${\mathcal{R}}_0$ by that same percentage [58].

Based on Equation (24) and the sensitivity indices in

Table 2. Sensitivity indices of ${\mathcal{R}}_0$ based on the parameter values given in Table 1 and Equation (24).

Parameter	Sensitivity Index
b	1
$\beta$	1
$k_1$	−0.4001
$k_2$	−0.3188
$k_4$	−0.2185
$\delta$	$-1.0660\times {10}^{-5}$
$\mu$	−1.3878
${\sigma }_1$	−0.0299

, the most influential parameters are $\mu ,b,\beta ,k_1,k_2.$ We make the following recommendations.

• Quarantining infected people may be an essential prevention measure to reduce the spread of COVID-19 infection.
• Reducing the transmission parameter is essential in reducing the transmission of the infection. To reduce the transmission parameter, compliance with protocols such as social distancing, wearing a face mask, and washing hands are necessary.

5. Numerical Analysis and Simulation

For completeness, we present in this Section the basic idea behind the method of solution employed to obtain solutions to the proposed system of time fractional model for COVID-19 given in (10).

5.1. Adams–Bashforth–Moulton Method for Fractional-Order

Adams–Bashforth–Moulton Method has been applied to integer order differential equations for many years. Here, we give brief description of this method applied to fractional-order differential equations. We refer the readers to see [59,60,61] for detailed results on stability and convergence. Given the initial valued fractional differential system of the form

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{d^{\alpha }}{d{t}^{\alpha }}u\left(t\right)& =\phi (t,u(t\left)\right)\hfill \end{array}\end{array}$$

(25)

and the initial condition

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill u^{\left(j\right)}\left(0\right)& =u_0^j,j=0,1,2,\cdots ,m-1\hfill \end{array}\end{array}$$

(26)

with $\alpha >0$, $m=\lceil \alpha \rceil$, and ${\displaystyle u_0^j}$, $j=0,1,2,\cdots ,m-1$ are given constants which are real numbers. Define the Volterra integral equation as

$$u\left(t\right)=\sum _{j=0}^{\lceil \alpha \rceil -1}\frac{t^j}{j!}{u}_0^j+\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^t{(t-\tau )}^{\alpha -1}\phi (\tau ,u\left(\tau \right))d\tau$$

(27)

Any continuous solution to the initial value fractional differential system (26) is also a solution of the Volterra integral equation given in (27) and vice versa. With the assumption of uniform grid with $t_n=nh$, $n=0,1,2,\cdots ,N$, $N\in \mathbb{N}$, $h>0$ is the step size, we give Adams–Bashforth–Moulton method of fractional-order type. We further assume that we first computed the approximation $u_h\left(t_i\right)\approx u\left(t_k\right)$, $k=1,2,\cdots ,n$ then find the approximation $u_h\left(t_{n+1}\right)$ using (27). The derivation for the case $\alpha =1$ of the one-step Adams–Bashforth–Moulton scheme is similar to the approach used here with some modification making use of the Volterra integral equation in (27). This is a Predictor–Corrector method. The main idea for the corrector is to approximate the integral term in (27) with product trapezoidal quadrature formula with nodes $t_k$, $k=0,1,2,\cdots ,n+1$. To get the predictor for this case, product rectangle rule is employed. Define

$$\begin{array}{c}\hfill {\gamma }_{k,n+1}=\left\{\begin{array}{ll}{n}^{\alpha +1}-(n-\alpha ){(n+1)}^{\alpha }& \mathrm{if}k=0,\\ \\ 1\hfill & \mathrm{if}k=n+1,\\ \\ {(n+2-k)}^{\alpha +1}+{(n-k)}^{\alpha +1}-2{(n+1-k)}^{\alpha +1}\hfill & \mathrm{if}1\le k\le n.\end{array}\right.\end{array}$$

(28)

and

$$\begin{array}{c}\hfill {\epsilon }_{k,n+1}=\frac{h^{\alpha }}{\alpha }\left[{(n+1-k)}^{\alpha }-{(n-k)}^{\alpha }\right].\end{array}$$

(29)

Then, the Adams–Bashforth–Moulton Method for fractional-order type which solves system (26) is given as follows:

$$\begin{array}{cc}{u}_h\left(t_{n+1}\right)\hfill & =\sum _{j=0}^{\lceil \alpha \rceil -1}\frac{t_{n+1}^j}{k!}{u}_0^j\hfill \\ \hspace{1em}\hfill & +\frac{h^{\alpha }}{\mathsf{\Gamma }(\alpha +2)}\phi \left(t_{n+1},u_h^p\left(t_{n+1}\right)\right)+\frac{h^{\alpha }}{\mathsf{\Gamma }(\alpha +2)}\sum _{k=0}^n{\gamma }_{k,n+1}\phi \left(t_k,u_h\left(t_k\right)\right)\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill u_h^p\left(t_{n+1}\right)& =& \sum _{j=0}^{\lceil \alpha \rceil -1}\frac{t_{n+1}^j}{j!}{u}_0^j+\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}\sum _{k=0}^n{\epsilon }_{k,n+1}\phi \left(t_k,u_h\left(t_k\right)\right),\hfill \end{array}$$

(31)

R. Roberto in [62] developed MATLAB subroutine for the implementation of the above method in MATLAB code fde12. With slight modification, we have employed the MATLAB function fde12 for numerical simulation of our proposed time-fractional COVID-19 model given in (26).

5.2. Numerical Simulation and Analysis

Our numerical simulations of the proposed COVID-19 model of time-fractional type are presented here with various analysis. The values of parameters used are as given in Table 1 unless otherwise stated. Effects of the fractional-order in the dynamics of the solution profiles obtained for various variables are investigated with meaningful observations. This is an important aspect of our study since the study of the nature and dynamics of COVID-19 (spread, resurgence etc.) is still ongoing, and not much is yet known. These intrinsic properties of COVID-19 could be captured using fractional derivative in time compared with integer order. Furthermore, we simulate with different values of fractional-order with different values of key parameters. The effects of these key parameters are expedient as well be able to tell what measure is effective in handling the spread of the disease and for crucial decisions given different values of fractional-orders as well.

5.2.1. Impact of Time-Fractional-Order on the Solution Profiles for the COVID-19 Model

The proposed COVID-19 model is simulated for different fractional-order values $0<\alpha \le 1$, using the baseline values of the parameters provided in Table 1. Figure 2a represents the solution of susceptible individuals $S\left(t\right)$, Figure 2b represents the solution of exposed individuals $E\left(t\right)$, Figure 2c represents the solution of exposed but quarantined individuals $Q_1\left(t\right)$, Figure 2d represents the solution of infected individuals $I\left(t\right)$, Figure 2e represents the solution of infected but quarantined individuals $Q_2\left(t\right)$, and Figure 2f represents the solution of recovered individuals $R\left(t\right)$. As shown in Figure 2a–f, the intrinsic properties of the fractional-orders are unique; the solution curves display delays in the epidemic peak for different fractional-order and flatten faster. For smaller $\alpha$ values, the impact is more pronounced. From a public health perspective, this observation is vital, as hospitalization would not be overloaded due to the flattening of the curve. Although the number of infected individuals decreases dramatically for smaller fractional-orders, the number of susceptible individuals increases, compare Figure 2a,d for $\alpha =0.6$.

5.2.2. Impact of the Effective Contact Rate on the Solution Profiles for the Coivd-19 Model

We investigated the contributions of the effective contact rate on the propagation of the COVID-19. We set the fractional-order as $\alpha =1,0.9,0.8,0.7$, and varied the $\beta$. The dynamics of our proposed model for the different values of the transmission parameter are displayed in Figure 3a–f, Figure 4a–f, Figure 5a–f and Figure 6a–f, respectively. As the transmission parameter rises, the infected individuals attain a higher peak level, Figure 3d,e, Figure 4d,e, Figure 5d,e, and Figure 6d,e. These effects are consistent with the sensitivity analysis result that shows the influence of the effective contact rate.

5.2.3. Impact of Quarantining Exposed Individuals on the Solution Profiles for the Coivd-19 Model

We ran simulations of the model to access the population-level impact of quarantining on the outbreak. The simulations results are displayed in Figure 7a–f, for $\alpha =1.$ The results for $\alpha =0.9$, $\alpha =0.8$, $\alpha =0.7$ are given in Figure 8a–f, Figure 9a–f and Figure 10a–f, respectively. The results show that quarantining exposed individuals reduces the number of infected individuals. For example, the default scenario result, with $\alpha =1$, shows a projected 800,000 cases at the peak time (Figure 7a). Without quarantining, about 1,400,000 cases will be recorded at the peak time. This result is approximately a $75\%$ increment of cases. A $40\%$ increase in quarantine baseline value results in approximately a 25% (600,000) reduction in infected cases, Figure 7a. Similarly, with the default value $k_1=0.1$, the projected peak daily infected cases was 700,000 observed for $\alpha =0.9$, 625,000 cases reached for $\alpha =0.8$, and 600,000 cases attained for $\alpha =0.7.$ A $40\%$ increase in quarantining rate from the baseline value results in about 20–25% reduction of infected cases, Figure 8d, Figure 9d and Figure 10d.

5.2.4. The Impact of the Loss of Immunity

There are several reports of re-infection of the COVID-19 cases. We simulated and investigated the impacts of transition rate from recovered to susceptible individuals (loss of immunity) on the dynamics of the solution profiles. The results are displayed in Figure 11a–f, Figure 12a–f, Figure 13a–f, and Figure 14a–f. We observed that loss of immunity increases the number of infected cases. Without it, the number of individuals in the recovered compartment rises.

5.2.5. The Impact of Quarantine Infected Individuals on the Solution Profiles for the Coivd-19 Model

The impact of quarantining infected individuals was also investigated. The simulation results for different $k_4$ are given in Figure 15a–f for $\alpha =1$, Figure 16a–f for $\alpha =0.9$, Figure 17a–f, for $\alpha =0.8$, and Figure 18a–f for $\alpha =0.7$. The results indicate that a high quarantining rate decreases the number of individuals in the infected (I) and exposed (E) classes. As displayed in Figure 15b,d, Figure 16b,d, and Figure 17b,d, without quarantine intervention, the exposed and infected populations grow faster with very high peak levels. In Figure 15, d, with a $40\%$ increase in quarantining rate from the baseline value, about 50,000 cases are observed at the peak time. Without quarantine, approximately 3,800,000 people will be infected at the peak time; thus, an over 400% increase of infected cases.

6. Conclusions

The COVID-19 spread has caused many damages, and it is still posing a significant challenge to individuals and countries. Many studies have been conducted to understand the transmission dynamics and control of the outbreak. The complexity of the disease is very high, and the rate at which it spreads is alarming. Based on this, we have proposed a fractional SEQIR model that captured many intrinsic properties of the transmission of the virus (COVID-19). Deep and extensive simulations in this paper provided insights and exposed hidden properties in the transmission of the coronavirus.

The proposed generalized compartmental COVID-19 model incorporates many key parameters to gain insight into the complexity of the disease. Effective contact rate, transition rate (from exposed quarantine and recovered to susceptible and infected quarantined individuals), quarantine rate, disease-induced death rate, natural death rate, natural recovery rate, and recovery rate of quarantine infected are part of our study. We simulated and discussed the parameters on the dynamics of the solution profiles. Detailed analysis of the proposed model was done, and the results on the existence and uniqueness of solutions, local and global stability analysis of the disease-free equilibrium, and sensitivity analysis are rigorously proven. An expression for the reproduction number, the threshold that measures the contagiousness of the outbreak, was established. For a numerical solution, we employed the generalized Adam–Bashforth–Moulton method developed for the fractional-order model. The effects of the critical parameters are also discussed in detail.

Based on the results, quarantining exposed and infected individuals is sensitive and significantly impacts the transmission of the virus. Thus, the importance of quarantine in mitigating the pandemic cannot be overlooked, and an effective strategy to control the spread of the disease. The results also suggest that the effective contact parameter influences the reproduction number and the spread of the disease significantly. The results also indicate that the effective contact parameter influences the reproduction number and the spread of the disease significantly. We highly recommend using both non-pharmaceutical and pharmaceutical measures to substantially reduce the effective contact rate and minimize the transmission of the virus. New surges in waves have started to occur currently in many places globally, including the USA. Our results indicate such surges with some appropriate choice of parameters and specific fractional-order.