Mathematical modelling of COVID-19 outbreak using caputo fractional derivative: stability analysis

Abstract

The novel coronavirus SARS-Cov-2 is a pandemic condition and poses a massive menace to health. The governments of different countries and their various prohibitory steps to restrict the virus's expanse have changed individuals' communication processes. Due to physical and financial factors, the population's density is more likely to interact and spread the virus. We establish a mathematical model to present the spread of the COVID-19 in worldwide. In this article, we propose a novel mathematical model (‘$\mathbb{S}\mathbb{L}\mathbb{I}{\mathbb{I}}_q{\mathbb{I}}_h\mathbb{R}\mathbb{P}$’) to assess the impact of using hospitalization, quarantine measures, and pathogen quantity in controlling the COVID-19 pandemic. We analyse the boundedness of the model's solution by employing the Laplace transform approach to solve the fractional Gronwall's inequality. To ensure the uniqueness and existence of the solution, we rely on the Picard-Lindelof theorem. The model's basic reproduction number, a crucial indicator of epidemic potential, is determined based on the greatest eigenvalue of the next-generation matrix. We then employ stability theory of fractional differential equations to qualitatively examine the model. Our findings reveal that both locally and globally, the endemic equilibrium and disease-free solutions demonstrate symptomatic stability. These results shed light on the effectiveness of the proposed interventions in managing and containing the COVID-19 outbreak.

KEYWORDS:

• Caputo fractional derivatives
• Karsnosels'kil's fixed point theorem
• Arzela Ascoli theorem
• Lyapunov function technique
• trace-determinant approach

MATHEMATICS SUBJECT CLASSIFICATIONS:

• 26A33
• 37C25
• 34D08

1. Introduction

The Wuhan Health Commission in the province of China alerted the World Health Organization (WHO) of a group of 27 cases of pneumonia with an unknown etiological cause, which was regularly reported to the Wuhan city market around the end of December 2019. Chinese officials, On January 7, 2020, a new strain of the Corona virus family was identified as the source of the epidemic. The virus's medical observation was given the designation COVID-19. WHO declared a pandemic on the 12th of March, [1–4]. Symptoms include sore throat, fever, dry cough, fatigue, and shortness of breath. In general, infections such as parasite infection by Cryptosporidia may occur without any symptoms.

The aim of the research was to construct a new ‘$\mathbb{S}\mathbb{L}\mathbb{I}{\mathbb{I}}_q{\mathbb{I}}_h\mathbb{R}\mathbb{P}$’ compartment model for predicting corona virus epidemiological dynamics. A model that presents the spreading dynamic of COVID-19 should be developed. In this case, Chen and Zhao [5] established a ‘SUQC’ mathematical model to describe the transmission dynamics of coronavirus. Similarly, López et al. [6] developed a new model of COVID-19, including healthy People, close contacts, and infected people.

Khajanchi et al. [7] developed the COVID-19 model with quarantine for the transmission dynamics of COVID-19 and calibrated the model for the four provinces of India with daily and accumulated cases. Further, they have conducted a detailed theoretical analysis of the basic reproductive numbers and predicted accumulated cases. Moreover, the study suggested that quarantine, unreported and reported individuals as well as intervention policies like lock down, and face masks can play a significant role in controlling the spread of corona virus. Also, S. Choi and Moran. Ki [8] developed a deterministic SEIHR mathematical model, in which they evaluated the effectiveness of control measures. Also, Lin et al. [9] modelled the outbreak in China with governmental action and individual reaction (city lock down, holiday extension, quarantine, and hospitalization). Also Wang and Yang [10] developed a model to investigate the corona virus disease outbreak.

Non-locality is the key areas of interest in the fractional calculus (FC) applications. There are fascinating phenomena called memory effects, this means that their conditions do not depend only on time and position even on the previous condition [11]. Fractional differential equation (FDE) has numerous important in different fields of study including physics, biophysics, aerodynamics, electron-analytical chemistry, capacitor theory, biology, and electrical circuits [12–20]. The Caputo derivative is one of the widely used definitions of fractional derivatives, and it differs from other fractional derivative definitions, such as the Riemann-Liouville derivative and the Grünwald-Letnikov derivative, in how it handles initial conditions and initial values.The Caputo derivative is most commonly used when we are dealing with fractional differential equations that involve initial conditions, as it leads to more physically meaningful solutions by considering the initial state of the system. The Riemann-Liouville derivative is used when initial values are involved in fractional differential equations. However, it can lead to solutions that are not physically meaningful in some cases, especially when dealing with real-world systems with initial conditions.The Grünwald-Letnikov derivative is commonly used in numerical methods for approximating fractional derivatives. It is straightforward to implement and provides accurate results when used with small values of h. However, it may not be suitable for all types of fractional differential equations, especially those involving initial conditions.

There have been several recent studies on epidemiological modelling using the FD equation [21–26]. The most used fractional derivative by mathematicians is Rieman-Liouville fractional derivative, but this procedure is not worthy for most real world-problems. Therefore, an important definition was introduced by Caputo [27]. Since the disease transmission rate is uncertain, we develop a mathematical model to present the spread of the COVID-19 in India and worldwide. This model is helpful for further modelling studies [28–30], and we can compare our model with previous mathematical models with actual data values (MoHFW India) (Corona virus disease (COVID-19) 2020) for the authentication of the model. The paper is organised as follows: Section 2 is related to the mathematical model formulation and some basic definitions and theorems. In Sections 3 and 4, we study the boundedness and uniqueness of the solution of the model, respectively. Section 5 present the disease free and endemic equilibrium points. Section 6 describe the basic reproduction number of the model.Local and Global stability are presented in Sections 7 and 8, respectively. Some numerical simulations are executed to illustrate the analytical results in Section 9. Finally, conclusions are presented in Section 10.

2. Model formulation

A fractional ‘$\mathbb{S}\mathbb{L}\mathbb{I}{\mathbb{I}}_q{\mathbb{I}}_h\mathbb{R}\mathbb{P}$’ mathematical model of the COVD-19 was developed in this part. We consider the whole population $\mathbb{N}$ into seven classes: susceptible individuals $\mathbb{S}$, latent individuals $\mathbb{L}$, infected individuals $\mathbb{I}$, quarantined individuals ${\mathbb{I}}_q$, hospitalized individuals ${\mathcal{I}}_h$, and the recovered individuals as $\mathbb{R}$. $\mathbb{P}$ represents the quantity of pathogens present in the environment.

The COVID-19 transmission dynamics are represented in Figure 1 flow chart.

Figure 1. Flow diagram of the COVID-19 Model (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ).

The virus can spread in two ways, according to scientists. Specifically, direct human-to-human and environmental-to-human transmission. Infected persons spread COVID-19 at a rate of π, according to the epidemic statistics, and are assigned to the latent individuals $\mathbb{L}$. When infected individuals sneeze or cough, the virus spreads to the surroundings they are in. Because pathogen $\mathbb{P}$ may survive in the environment for several days, $\mathbb{S}$ persons who are in close touch with it are likely to get infected. The rate of transmission from $\mathbb{S}$ to $\mathbb{L}$ due to contact with $\mathbb{P}$ is represented by π. Infected people with serious illnesses are admitted to hospitals at a rate of φ. In addition, individuals that were quarantined or hospitalised recovered from COVID-19 naturally at rates of ϵ and λ, respectively, and were placed in the recovered class $\mathbb{R}$.

The Caputo derivative is used to develop a mathematical model is given by (1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) with the following initial conditions $\mathbb{S}\left(0\right)>0,\mathbb{L}\left(0\right)\ge 0,\mathbb{I}\left(0\right)\ge 0,{\mathbb{I}}_q\left(0\right)\ge 0,{\mathbb{I}}_h\left(0\right)\ge 0,\mathbb{R}\left(0\right)\ge 0,\mathbb{P}\left(0\right)>0.$

Detailed description of relevant parameters and state variables of the newly proposed COVID-19 model are well explained in Table 1.

The paper is arranged as: Section 1 gives the introduction. Section 2 includes model formulation and basic definitions. Boundedness and uniqueness for the solution are considered in Sections 3 and 4, respectively. Section 5 is devoted to the disease free equilibrium point (DFE) point. Basic reproduction is found in Section 6. The local and global stabilities are discussed in Sections 7 and 8, respectively. The numerical discussion is considered in Section 9. The last section is the conclusion of our work.

Table 1. Biological interpretations of parameters.

Parameter	Biological interpretations
σ	The human population recruited rate.
υ	Elaborated the natural deaths rate.
$\upsilon {}_i$	Represents the mortality rate in the $\mathbb{I}$ due to the coronavirus.
$\upsilon {}_q$	Describe the death rate due to coronavirus in the ${\mathbb{I}}_q$ class.
$\upsilon {}_h$	Denotes the death rate due to COVID-19 in the ${\mathbb{I}}_h$ class.
β	The rate of transmission from $\mathbb{S}$ to $\mathbb{L}$ is due to contact with pathogens in the environment.
τ	The rate of progression from latent individuals ($\mathbb{L}$) to infected ($\mathbb{I}$).
ϕ and φ	${\mathbb{I}}_q$ and ${\mathbb{I}}_h$ rates for infectious class were described, respectively.
ϵ	Shows the rates of recovery of the class ${\mathbb{I}}_q$.
λ	Describe the rates of recovery of the class ${\mathbb{I}}_h$.
δ	Shows the rates of recovery of the class $\mathbb{I}$.
ϖ	Denotes the recovery rate of the class $\mathbb{L}$.
θ	Mild cases transmission rate.
α	Rate of virus spread to environment by infected class.
π	Progression rate from $\mathbb{L}$ back to $\mathbb{S}$ due to robust immune system.
$\upsilon {}_p$	Natural death rate of pathogens in the environment.

Definition 2.1

[31]

Let $\Psi \left(t\right)$ be a smooth function, then the Caputo fractional derivative of $\Psi \left(t\right)$ is ${}_0^c{\mathcal{D}}_t^{\gamma }\left[\Psi \right(t\left)\right]=\frac{1}{\Gamma (\mathcal{K}-\gamma )}{\int }_0^t{(t-\zeta )}^{\mathcal{K}-\gamma -1}{\Psi }^{\left(\mathcal{K}\right)}\left(\zeta \right)\mathrm{d}\zeta ,$ where, $\mathcal{K}-1\le \gamma \le \mathcal{K}$.

Definition 2.2

[32]

Let $\Psi :{\mathbb{R}}^{+}\to \mathbb{R}$ be a smooth function, then the fractional integral of order $\gamma >0$ for $\Psi \left(t\right)$ is defined by ${}_0{\mathcal{I}}_t^{\gamma }\Psi \left(t\right)=\frac{1}{\Gamma \left(\gamma \right)}{\int }_0^t{(t-\zeta )}^{\gamma -1}\Psi \left(\zeta \right)\mathrm{d}\zeta ,$ and ${\mathcal{I}}_t^{-\gamma }{\mathcal{D}}_t^{\gamma }\left(\Psi \right(t\left)\right)={\mathcal{D}}_t^{\gamma }\Psi \left(t\right),$ where $0<\gamma \le 1$.

The following is an expression of a general n-dimension autonomous system: (2) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\alpha }\left[\Psi t\right]=\mathcal{G}\left(\Psi \right(t\left)\right)\\ \Psi \left(t_0\right)={\Psi }_0.\end{array}$(2)

Theorem 2.3

[31]

For an equilibrium point ${\mathcal{E}}^{\ast }$ for the system in Caputo sense (2(2) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\alpha }\left[\Psi t\right]=\mathcal{G}\left(\Psi \right(t\left)\right)\\ \Psi \left(t_0\right)={\Psi }_0.\end{array}$(2) ) and $\mathrm{\Omega }\in {\mathcal{R}}^n$be the domain such that ${\mathcal{F}}^{\ast }\in \Xi$ and let $\mathcal{Z}(t,\mathcal{E}):[0,\mathrm{\infty })\times \Xi \to \mathcal{R}$, be a continuously differentiable function. If (3) ${\Psi }_1\left(\mathcal{E}\right)\le \mathcal{Z}(t,\mathcal{E})\le {\Psi }_2\left(\mathcal{E}\right).$(3) Then, (4) ${}_0^C{\mathcal{D}}_t^{\alpha }\mathcal{Z}(t,\mathcal{E})\le -{\Psi }_3\left(\mathcal{E}\right).$(4) ∀$z\in \Xi$ and $\alpha \in (0,1)$, where $\Psi {}_1\left(\mathcal{E}\right)$, $\Psi {}_2\left(\mathcal{E}\right)$ and $\Psi {}_3\left(\mathcal{E}\right)$ are positive definite continuous functions in Ξ, then the point ${\mathcal{F}}^{\ast }$ of (2(2) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\alpha }\left[\Psi t\right]=\mathcal{G}\left(\Psi \right(t\left)\right)\\ \Psi \left(t_0\right)={\Psi }_0.\end{array}$(2) ) is stable asymptotically.

3. Boundedness of the system

The size of each population class varies over time in our suggested model (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ), and the overall population size $\mathbb{N}$ is provided by $\mathbb{N}=\mathbb{S}+\mathbb{L}+\mathbb{I}+{\mathbb{I}}_h+{\mathbb{I}}_q+\mathbb{R}+\mathbb{P}.$ Using the Caputo fractional derivative and the Caputo derivative linearity property, we get (5) $\begin{array}{rl}{}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{N}\right(t\left)\right]& {=}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{S}\right(t\left)\right]+{}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{L}\right(t\left)\right]+{}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{I}\right(t\left)\right]+{}_0^c{\mathcal{D}}_t^{\alpha }\left[{\mathbb{I}}_q\right(t\left)\right]\\ & {+}_0^c{\mathcal{D}}_t^{\alpha }\left[{\mathbb{I}}_h\right(t\left)\right]+{}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{R}\right(t\left)\right]+{\mathcal{D}}_t^{\alpha }\left[\mathbb{P}\right],\\ {}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{N}\right(t\left)\right]& =\tau -\upsilon \mathbb{S}\left(t\right)-\upsilon \mathbb{L}\left(t\right)-\upsilon \mathbb{I}\left(t\right)-{\upsilon }_i\mathbb{I}\left(t\right)-\upsilon {\mathbb{I}}_q\left(t\right)-{\upsilon }_q{\mathbb{I}}_q\left(t\right)\\ & -\upsilon {\mathbb{I}}_h\left(t\right)+{\upsilon }_h{\mathbb{I}}_h\left(t\right)-\upsilon \mathbb{R}\left(t\right)-{\upsilon }_p\mathbb{P}\left(t\right),\\ {}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{N}\right(t\left)\right]& =\tau -\upsilon \left(\mathbb{S}\right(t)+\mathbb{L}(t)+\mathbb{I}(t)+{\mathbb{I}}_q(t)+{\mathbb{I}}_h(t)+\mathbb{R}(t)+\mathbb{P}(t\left)\right)\\ & -{\upsilon }_i\mathbb{I}\left(t\right)-{\upsilon }_q{\mathbb{I}}_q\left(t\right)-{\upsilon }_h{\mathbb{I}}_h\left(t\right),\\ {}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{N}\right(t\left)\right]& =\tau -\upsilon \mathbb{N}\left(t\right)-{\upsilon }_i\mathbb{I}\left(t\right)-{\upsilon }_q{\mathbb{I}}_q\left(t\right)-{\upsilon }_h{\mathbb{I}}_h\left(t\right).\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{N}\right]+\upsilon \mathbb{N}& \le \sigma .\end{array}$(5) We solve the Gronwall's inequality (5(5) $\begin{array}{rl}{}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{N}\right(t\left)\right]& {=}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{S}\right(t\left)\right]+{}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{L}\right(t\left)\right]+{}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{I}\right(t\left)\right]+{}_0^c{\mathcal{D}}_t^{\alpha }\left[{\mathbb{I}}_q\right(t\left)\right]\\ & {+}_0^c{\mathcal{D}}_t^{\alpha }\left[{\mathbb{I}}_h\right(t\left)\right]+{}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{R}\right(t\left)\right]+{\mathcal{D}}_t^{\alpha }\left[\mathbb{P}\right],\\ {}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{N}\right(t\left)\right]& =\tau -\upsilon \mathbb{S}\left(t\right)-\upsilon \mathbb{L}\left(t\right)-\upsilon \mathbb{I}\left(t\right)-{\upsilon }_i\mathbb{I}\left(t\right)-\upsilon {\mathbb{I}}_q\left(t\right)-{\upsilon }_q{\mathbb{I}}_q\left(t\right)\\ & -\upsilon {\mathbb{I}}_h\left(t\right)+{\upsilon }_h{\mathbb{I}}_h\left(t\right)-\upsilon \mathbb{R}\left(t\right)-{\upsilon }_p\mathbb{P}\left(t\right),\\ {}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{N}\right(t\left)\right]& =\tau -\upsilon \left(\mathbb{S}\right(t)+\mathbb{L}(t)+\mathbb{I}(t)+{\mathbb{I}}_q(t)+{\mathbb{I}}_h(t)+\mathbb{R}(t)+\mathbb{P}(t\left)\right)\\ & -{\upsilon }_i\mathbb{I}\left(t\right)-{\upsilon }_q{\mathbb{I}}_q\left(t\right)-{\upsilon }_h{\mathbb{I}}_h\left(t\right),\\ {}_0^c{\mathcal{D}}_t^{\alpha }\left[\mathbb{N}\right(t\left)\right]& =\tau -\upsilon \mathbb{N}\left(t\right)-{\upsilon }_i\mathbb{I}\left(t\right)-{\upsilon }_q{\mathbb{I}}_q\left(t\right)-{\upsilon }_h{\mathbb{I}}_h\left(t\right).\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{N}\right]+\upsilon \mathbb{N}& \le \sigma .\end{array}$(5) ) using the Laplace transform with initial conditions ${}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{N}\right]\ge 0$, (6) $\begin{array}{rl}& \mathcal{L}\left\{{}_0^c{\mathcal{D}}_t^{\gamma }\right[\mathbb{N}\left]\right\}+\upsilon \mathcal{L}\left[\mathbb{N}\right]\le \mathcal{L}\left\{\sigma \right\}.\\ & S^{\gamma }\mathcal{L}\left[\mathbb{N}\right]-\sum _{k=0}^{n-1}{S}^{\gamma -k-1}{\mathbb{N}}^k\left(t_0\right)+\upsilon \mathcal{L}\left[\mathbb{N}\right]\le \frac{\sigma }{S}.\\ & \mathcal{L}\left\{\mathbb{N}\right\}\le \frac{\sigma }{S(\upsilon +S^{\gamma })}+\sum _{k=0}^{n-1}\frac{S^{\gamma -k-1}{\mathbb{N}}^k\left(t_0\right)}{(\upsilon +S^{\gamma })}.\end{array}$(6) We obtain, when we split (6(6) $\begin{array}{rl}& \mathcal{L}\left\{{}_0^c{\mathcal{D}}_t^{\gamma }\right[\mathbb{N}\left]\right\}+\upsilon \mathcal{L}\left[\mathbb{N}\right]\le \mathcal{L}\left\{\sigma \right\}.\\ & S^{\gamma }\mathcal{L}\left[\mathbb{N}\right]-\sum _{k=0}^{n-1}{S}^{\gamma -k-1}{\mathbb{N}}^k\left(t_0\right)+\upsilon \mathcal{L}\left[\mathbb{N}\right]\le \frac{\sigma }{S}.\\ & \mathcal{L}\left\{\mathbb{N}\right\}\le \frac{\sigma }{S(\upsilon +S^{\gamma })}+\sum _{k=0}^{n-1}\frac{S^{\gamma -k-1}{\mathbb{N}}^k\left(t_0\right)}{(\upsilon +S^{\gamma })}.\end{array}$(6) ) into partial fractions. (7) $\begin{array}{rl}\mathcal{L}\left\{\mathbb{N}\right\}& \le \sigma (\frac{1}{S}-\frac{S^{\gamma -1}}{(S^{\gamma }+\upsilon )})+\sum _{k=0}^{n-1}\frac{S^{\gamma -k-1}{\mathbb{N}}^k\left(0\right)}{\upsilon +S^{\gamma }}.\\ \mathcal{L}\left\{\mathbb{N}\right\}& \le \sigma (\frac{1}{S}-\frac{1}{S(1+\frac{\upsilon }{S^{\gamma }})})+\sum _{k=0}^{n-1}\frac{{\mathbb{N}}^k\left(0\right)}{S^{k+1}(1+\frac{\upsilon }{S^{\gamma }})}.\end{array}$(7) By expanding through the Talyor series, we get $\frac{1}{S(1+\frac{\upsilon }{S^{\gamma }})}=\sum _{n=0}^{\mathrm{\infty }}{\left(\frac{-\upsilon }{S^{\gamma }}\right)}^n.$ Therefore, (7(7) $\begin{array}{rl}\mathcal{L}\left\{\mathbb{N}\right\}& \le \sigma (\frac{1}{S}-\frac{S^{\gamma -1}}{(S^{\gamma }+\upsilon )})+\sum _{k=0}^{n-1}\frac{S^{\gamma -k-1}{\mathbb{N}}^k\left(0\right)}{\upsilon +S^{\gamma }}.\\ \mathcal{L}\left\{\mathbb{N}\right\}& \le \sigma (\frac{1}{S}-\frac{1}{S(1+\frac{\upsilon }{S^{\gamma }})})+\sum _{k=0}^{n-1}\frac{{\mathbb{N}}^k\left(0\right)}{S^{k+1}(1+\frac{\upsilon }{S^{\gamma }})}.\end{array}$(7) ) becomes (8) $\begin{array}{rl}\mathcal{L}\left\{\mathbb{N}\right\}& \le \sigma (\frac{1}{S}-\frac{1}{S}\sum _{n=0}^{\mathrm{\infty }}{\left(\frac{-\upsilon }{S^{\gamma }}\right)}^n)+\sum _{k=0}^{n-1}\frac{1}{S^{k+1}}{\mathbb{N}}^k\left(0\right)\sum _{n=0}^{\mathrm{\infty }}{\left(\frac{-\upsilon }{S^{\gamma }}\right)}^n.\\ \mathcal{L}\left\{\mathbb{N}\right\}& \le \sigma (\frac{1}{S}-\sum _{n=0}^{\mathrm{\infty }}\frac{(-\upsilon {)}^n}{S^{\mathrm{n\gamma }+1}})+\sum _{k=0}^{n-1}\sum _{n=0}^{\mathrm{\infty }}\frac{(-\upsilon {)}^n{\mathbb{N}}^k(0)}{S^{\mathrm{n\gamma }+k+1}}.\end{array}$(8) The inverse Laplace transform is applied to (8(8) $\begin{array}{rl}\mathcal{L}\left\{\mathbb{N}\right\}& \le \sigma (\frac{1}{S}-\frac{1}{S}\sum _{n=0}^{\mathrm{\infty }}{\left(\frac{-\upsilon }{S^{\gamma }}\right)}^n)+\sum _{k=0}^{n-1}\frac{1}{S^{k+1}}{\mathbb{N}}^k\left(0\right)\sum _{n=0}^{\mathrm{\infty }}{\left(\frac{-\upsilon }{S^{\gamma }}\right)}^n.\\ \mathcal{L}\left\{\mathbb{N}\right\}& \le \sigma (\frac{1}{S}-\sum _{n=0}^{\mathrm{\infty }}\frac{(-\upsilon {)}^n}{S^{\mathrm{n\gamma }+1}})+\sum _{k=0}^{n-1}\sum _{n=0}^{\mathrm{\infty }}\frac{(-\upsilon {)}^n{\mathbb{N}}^k(0)}{S^{\mathrm{n\gamma }+k+1}}.\end{array}$(8) ) (9) $\mathbb{N}\le \sigma {\mathcal{L}}^{-1}\left\{\frac{1}{S}\right\}-\sigma \sum _{n=0}^{\mathrm{\infty }}(-\upsilon {)}^n{\mathcal{L}}^{-1}\left\{\frac{1}{S^{\mathrm{n\gamma }+1}}\right\}+\sum _{k=0}^{n-1}\sum _{n=0}^{\mathrm{\infty }}(-\upsilon {)}^n{\mathbb{N}}^k\left(0\right){\mathcal{L}}^{-1}\left\{\frac{1}{S^{\mathrm{n\gamma }+k+1}}\right\}.$(9) Since $\mathcal{L}\left\{t^k\right\}=\frac{k!}{S^{k+1}}=\frac{\Gamma (k+1)}{S^{k+1}}$, or ${\mathcal{L}}^{-1}\left\{\frac{1}{S^{k+1}}\right\}=\frac{t^k}{\Gamma (k+1)}$. Therefore, (10) $\mathcal{N}\le \sigma -\tau \sum _{n=0}^{\mathrm{\infty }}\frac{(-\upsilon t^{\gamma }{)}^n}{\mathrm{\Gamma }(\mathrm{n\alpha }+1)}+\sum _{k=0}^{n-1}\sum _{n=0}^{\mathrm{\infty }}\frac{(-\upsilon {)}^n{\mathbb{N}}^k(0)t^{\mathrm{\gamma n}+k}}{\mathrm{\Gamma }(\mathrm{n\gamma }+k+1)}.$(10) We are now utilising the Mittag–Leffler function definition (11) $\mathbb{N}\le \sigma (1-{\mathbb{E}}_{\gamma ,1}(-\upsilon t^{\gamma }\left)\right)+{\mathbb{E}}_{\gamma ,k+1}(-\upsilon t^{\gamma })t^k{\mathbb{N}}^k\left(0\right),$(11) where, ${\mathbb{E}}_{\alpha ,1}(-\mu t^{\alpha })$ and ${\mathbb{E}}_{\alpha ,k+1}(-\mu t^{\alpha })$ are the series of entire function which converges for any argument, hence, the solution of Model (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) is bounded. $\left\{\right(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})\in {\mathcal{R}}_{+}^7,0\le \mathbb{N}\le \sigma (1-E_{\gamma ,1}(-\upsilon t^{\gamma }\left)\right)+E_{\gamma ,k+1}(-\upsilon t^{\gamma }\left){\mathbb{N}}^k\right(0\left)t^k\right\}.$

4. Uniqueness of the solution

Using fixed point theory, we discuss the uniqueness of the model's solution. Let us rewrite the proposed model (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) in the following way. (12) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]={\mathrm{{\rm Y}}}_1(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P}),\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\mathrm{{\rm Y}}}_2(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P}),\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]={\mathrm{{\rm Y}}}_3(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P}),\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]={\mathrm{{\rm Y}}}_4(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P}),\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]={\mathrm{{\rm Y}}}_5(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P}),\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]={\mathrm{{\rm Y}}}_6(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P}),\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]={\mathrm{{\rm Y}}}_7(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P}),\end{array}$(12) where (13) $\{\begin{array}{l}{\mathrm{{\rm Y}}}_1(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})=\sigma -\upsilon \mathbb{S}-{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{M})-\pi \mathbb{S}\mathbb{P},\\ {\mathrm{{\rm Y}}}_2(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L}(t\left)\right)-(\upsilon +\varpi +\tau )\mathbb{L}\left(t\right)+\pi \mathcal{S}\mathbb{P},\\ {\mathrm{{\rm Y}}}_3(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {\mathrm{{\rm Y}}}_4(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {\mathrm{{\rm Y}}}_5(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {\mathrm{{\rm Y}}}_6(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})=\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+\delta \mathbb{I}+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {\mathrm{{\rm Y}}}_7(t,\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P}.\end{array}$(13) We can write the model (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) as (14) $\{\begin{array}{ll}{}_0^c{\mathcal{D}}_t^{\gamma }\Psi \left(t\right)=\mathcal{Y}(t,\Psi );& \Psi \left(0\right)={\Psi }_0,\\ \mathcal{Y}(t,\Psi )=\mathrm{B\Psi }+A\left(\Psi \right)+c,& \end{array}$(14) where $\begin{array}{rl}\Psi & =\Psi \left(t\right)={(\mathbb{S},\mathbb{L},\mathbb{I},\mathbb{I}{}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})}^T,\\ A& =\left(\begin{array}{cc}-{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\theta \mathbb{L}+\mathbb{I})& -\pi \mathbb{S}\mathbb{P}\\ {\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\theta \mathbb{M}+\mathbb{I})& \pi \mathbb{S}\mathbb{P}\\ 0& 0\\ 0& 0\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right),\\ B& =\left(\begin{array}{ccccccc}-\upsilon & 0& 0& 0& 0& 0& 0\\ 0& -a_1& 0& 0& 0& 0& 0\\ 0& \tau & -a_2& 0& 0& 0& 0\\ 0& 0& \varphi & -a_3& 0& 0& 0\\ 0& 0& \phi & 0& -a_4& 0& 0\\ 0& \varpi & \delta & ϵ& \lambda & -\upsilon & 0\\ 0& 0& \alpha & 0& 0& 0& -{\upsilon }_p\end{array}\right),\\ c& =(\sigma ,0,0,0,0,0{)}^T,a_1=(\upsilon +\varpi +\tau ),a_2=(\phi +\delta +\upsilon +{\upsilon }_i+\varphi +\alpha ),\\ a_3& =(\upsilon +ϵ+{\upsilon }_q)\mathrm{and}{a}_4=(\upsilon +{\upsilon }_h+\lambda ).\end{array}$

Theorem 4.1

Lipschitz continuity is satisfied by System (14(14) $\{\begin{array}{ll}{}_0^c{\mathcal{D}}_t^{\gamma }\Psi \left(t\right)=\mathcal{Y}(t,\Psi );& \Psi \left(0\right)={\Psi }_0,\\ \mathcal{Y}(t,\Psi )=\mathrm{B\Psi }+A\left(\Psi \right)+c,& \end{array}$(14) ).

Proof.

$\begin{array}{rl}\left|\mathcal{Y}\right(t,\Psi )-\mathcal{Y}(t,{\Psi }^{\ast }\left)\right|& =|A(\Psi -{\Psi }^{\ast })+g\left(\Psi \right)-g\left({\Psi }^{\ast }\right)|\\ & \le (\Vert A\Vert +1)\Vert \Psi -{\Psi }^{\ast }\Vert .\end{array}$ This implies (15) $\Vert \mathcal{Y}(t,\Psi )-\mathcal{Y}(t,{\Psi }^{\ast })\Vert \le L\Vert \Psi -{\Psi }^{\ast }\Vert ,L=\Vert A\Vert +1<\mathrm{\infty }.$(15) Hence, the function $\mathcal{Y}$ is continuous and bounded.

The following theorem was constructed using the Picard–Lindelof theorem [33].

Theorem 4.2

Let $0<\gamma <1,\mathrm{\Lambda }=[0,l^{\ast }]\subseteq R$ and $J=|\Psi \left(t\right)-\Psi \left(0\right)|\le \mathcal{K}$.

Let $g:\mathrm{\Lambda }\times J\to R$ be a continuous and bounded function, that is ∃ $M>0$ such that $|g(t,\mathcal{Y})|\le M$, since $\mathcal{Y}$ satisfies Lipschitz continuity. If $\mathcal{K}L<M$, then there exist a unique $\Psi \in C[0,l^{\ast }]$ that satisfied the problem (14(14) $\{\begin{array}{ll}{}_0^c{\mathcal{D}}_t^{\gamma }\Psi \left(t\right)=\mathcal{Y}(t,\Psi );& \Psi \left(0\right)={\Psi }_0,\\ \mathcal{Y}(t,\Psi )=\mathrm{B\Psi }+A\left(\Psi \right)+c,& \end{array}$(14) ), where $l^{\ast }=min\{l,(\frac{\mathcal{K}\Gamma (\gamma +1)}{M}{)}^{\frac{1}{\gamma }}\}$.

Proof.

Let $T=\{\Psi \in C[o,l^{\ast }]:\Vert \Psi (t)-\Psi (0)\Vert \le \mathcal{K}\}$, clearly T is complete metric space because $T\subseteq R$ and its closed set.

The following is an equivalent integral equation for the system (14(14) $\{\begin{array}{ll}{}_0^c{\mathcal{D}}_t^{\gamma }\Psi \left(t\right)=\mathcal{Y}(t,\Psi );& \Psi \left(0\right)={\Psi }_0,\\ \mathcal{Y}(t,\Psi )=\mathrm{B\Psi }+A\left(\Psi \right)+c,& \end{array}$(14) ): (16) $\begin{array}{rl}{}_0^c{\mathcal{D}}_t^{-\gamma }\left[{}_0^c{\mathcal{D}}_t^{\gamma }\Psi \right(t\left)\right]& ={}_0^c{\mathcal{D}}_t^{-\gamma }\mathcal{Y}(t,\Psi );\\ \Psi & =\Psi \left(0\right)+I_t^{\gamma }\mathcal{Y}(t,\Psi )\\ & ={\Psi }_0+\frac{1}{\Gamma \left(\gamma \right)}{\int }_0^t{(t-\zeta )}^{\gamma -1}\mathcal{Y}(\zeta ,\Psi (\zeta \left)\right)\mathrm{d}\zeta .\end{array}$(16) We define ${\mathcal{T}}^{\ast }$ in T as (17) ${\mathcal{T}}^{\ast }\left[\Psi \right]={\Psi }_0+\frac{1}{\Gamma \left(\gamma \right)}{\int }_0^t{(t-\zeta )}^{\gamma -1}\mathcal{Y}(\zeta ,\Psi )\mathrm{d}\zeta .$(17) We need to verify that (17(17) ${\mathcal{T}}^{\ast }\left[\Psi \right]={\Psi }_0+\frac{1}{\Gamma \left(\gamma \right)}{\int }_0^t{(t-\zeta )}^{\gamma -1}\mathcal{Y}(\zeta ,\Psi )\mathrm{d}\zeta .$(17) ) satisfies the contraction mapping principle hypothesis.

Firstly, we show that ${\mathcal{T}}^{\ast }:T\to T$ (18) $\begin{array}{rl}|{\mathcal{T}}^{\ast }\left[\Psi \right(t\left)\right]-\Psi (0\left)\right|& =\left|\frac{1}{\Gamma \left(\gamma \right)}{\int }_0^t{(t-\zeta )}^{\gamma -1}\mathcal{Y}(\zeta ,\Psi )\mathrm{d}\zeta \right|\\ & \le \frac{1}{\Gamma \left(\gamma \right)}{\int }_0^t{(t-\zeta )}^{\gamma -1}\left|\mathcal{Y}(\zeta ,\Psi )\right|\mathrm{d}\zeta \\ & \le \frac{1}{\Gamma \left(\gamma \right)}{\int }_0^t{(t-\zeta )}^{\gamma -1}M\mathrm{d}\zeta \\ & =\frac{M}{\Gamma (\gamma +1)}{t}^{\gamma }\\ & =\frac{M}{\Gamma (\gamma +1)}(l^{\ast }{)}^{\gamma }\\ & =\frac{M}{\Gamma (\gamma +1)}\frac{\mathcal{K}\Gamma (\gamma +1)}{M}.\\ |{\mathcal{T}}^{\ast }\left[\Psi \right(t\left)\right]-\Psi (0\left)\right|& \le \mathcal{K}.\end{array}$(18) Equivalently, (19) $\Psi \left(0\right)-\mathcal{K}\le {\mathcal{T}}^{\ast }\left[\Psi \right]\left(t\right)\le \Psi \left(0\right)+\mathcal{K},\mathrm{\forall }t\in [0,l^{\ast }].$(19) Clearly, the operator ${\mathcal{T}}^{\ast }$ maps T onto itself.

We need to demonstrate that T is a contraction. $\begin{array}{rl}\left|{\mathcal{T}}^{\ast }\right[\Psi ]-{\mathcal{T}}^{\ast }[{\Psi }^{\ast }\left]\right|& =\left|\frac{1}{\Gamma \left(\gamma \right)}{\int }_0^t{(t-\zeta )}^{\gamma -1}[\mathcal{Y}(\zeta ,\Psi (\zeta \left)\right)-\mathcal{Y}(\zeta ,{\Psi }^{\ast }(\zeta \left)\right)]\mathrm{d}\zeta \right|\\ & \le \frac{1}{\Gamma \left(\gamma \right)}{\int }_0^t{(t-\zeta )}^{\gamma -1}\left|[\mathcal{Y}(\zeta ,\Psi (\zeta \left)\right)-\mathcal{Y}(\zeta ,{\Psi }^{\ast }(\zeta \left)\right)]\right|\mathrm{d}\zeta .\\ & \le \frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\zeta )}^{\gamma -1}L\Vert \Psi -{\Psi }^{\ast }\Vert \mathrm{d}\zeta \\ & =\frac{L}{\Gamma \left(\gamma \right)}\Vert \Psi -{\Psi }^{\ast }\Vert {\int }_0^t{(t-\zeta )}^{\gamma -1}\mathrm{d}\zeta \\ & =\frac{L}{\Gamma (\gamma +1)}\Vert \Psi -{\Psi }^{\ast }\Vert t^{\gamma }\\ & =\frac{L}{\Gamma (\gamma +1)}\Vert \Psi -{\Psi }^{\ast }\Vert (l^{\ast }{)}^{\gamma }\\ & =\frac{L}{\Gamma (\gamma +1)}\Vert \Psi -{\Psi }^{\ast }\Vert \frac{\mathcal{K}\Gamma (\gamma +1)}{M}\\ & =\frac{\mathcal{K}L}{M}\Vert \Psi -{\Psi }^{\ast }\Vert .\end{array}$ Thus (20) $\Vert {\mathcal{T}}^{\ast }\left[\Psi \right]\left(t\right)-{\mathcal{T}}^{\ast }\left[{\Psi }^{\ast }\right]\left(t\right)\Vert \le \frac{L\mathcal{K}}{M}\Vert \Psi -{\Psi }^{\ast }\Vert .$(20) According to hypothesis $\frac{L\mathcal{K}}{M}<1$, T is contraction and has a single fixed point. Hence, system (14(14) $\{\begin{array}{ll}{}_0^c{\mathcal{D}}_t^{\gamma }\Psi \left(t\right)=\mathcal{Y}(t,\Psi );& \Psi \left(0\right)={\Psi }_0,\\ \mathcal{Y}(t,\Psi )=\mathrm{B\Psi }+A\left(\Psi \right)+c,& \end{array}$(14) ) has unique solution.

5. Disease free and endemic equilibrium points

We set the right-hand sides of all equations in (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) to zero to determine $E^0$ of the system (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ). (21) $\{\begin{array}{l}\sigma -\mu \mathbb{S}-\frac{\pi \mathbb{S}}{\mathbb{N}}(\theta \mathbb{L}+\mathbb{I})-\text{β}\mathbb{S}\mathbb{P}=0,\\ {\displaystyle \frac{\pi \mathbb{S}}{\mathbb{N}}}(\theta \mathbb{L}+\mathbb{I})-(\upsilon +\varpi +\tau )\mathbb{L}+\text{β}\mathbb{S}\mathbb{P}=0,\\ \sigma \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I}=0,\\ \varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q)\mathbb{I}{}_q=0,\\ \phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h)\mathbb{I}{}_h=0,\\ \delta \mathbb{I}+\varpi \mathbb{L}+ϵ{\mathbb{I}}_q+\lambda \mathbb{I}{}_h-\upsilon \mathbb{R}=0,\\ \alpha \mathbb{I}-{\upsilon }_p\mathbb{P}=0.\end{array}$(21) Let $\mathbb{L}=\mathbb{I}={\mathbb{I}}_q={\mathbb{I}}_h=0$, which means that $\mathbb{P}=0$ and $\mathbb{R}=0$ too. Hence we have (22) $\sigma -\upsilon \mathbb{S}=0\Rightarrow \mathbb{S}=\frac{\sigma }{\upsilon }.$(22) Therefore, $E^0$ of the considered model (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) is given by $E^0=(\frac{\sigma }{\upsilon },0,0,0,0,0,0)$, which is always feasible.

To find the model endemic equilibrium point (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ), where at least one disease class is not zero. The steps below must be followed. Let $E^{\ast }=({\mathbb{S}}^{\ast },{\mathbb{L}}^{\ast },{\mathbb{I}}^{\ast },{\mathbb{I}}_q^{\ast },{\mathbb{I}}_h^{\ast },{\mathbb{R}}^{\ast },{\mathbb{P}}^{\ast })$ be the endemic equilibrium point of the model (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) and is given by (23) $\{\begin{array}{l}{\mathbb{S}}^{\ast }={\displaystyle \frac{1}{\upsilon }}(\sigma -{\displaystyle \frac{(\upsilon +\varpi +\tau )(\upsilon +\delta +\phi +{\upsilon }_i+\varphi +\alpha ){\mathbb{I}}^{\ast }}{\tau }}),\\ {\mathbb{L}}^{\ast }={\displaystyle \frac{(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi ){\mathbb{I}}^{\ast }}{\tau }},\\ {\mathbb{I}}^{\ast }={\displaystyle \frac{\mathrm{\tau \sigma }}{\left(\upsilon +\varpi +\tau \right)(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )}}-{\displaystyle \frac{\mathrm{\tau \upsilon }{\mathbb{N}}^{\ast }}{\text{β}(\theta (\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )+\tau )}},\\ {\mathbb{I}}_q^{\ast }={\displaystyle \frac{\varphi {\mathbb{I}}^{\ast }}{(\upsilon +ϵ+{\upsilon }_q)}},\\ {\mathbb{I}}_h^{\ast }={\displaystyle \frac{\phi {\mathbb{I}}^{\ast }}{(\upsilon {}_h+\lambda +\upsilon )}},\\ {\mathbb{R}}^{\ast }=\left[{\displaystyle \frac{\delta }{\upsilon }}+{\displaystyle \frac{\varpi (\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )+\mathrm{\lambda \phi }(\upsilon +ϵ+{\upsilon }_q)+\mathrm{ϵ\varphi }(\upsilon +{\upsilon }_h+\lambda )}{\mathrm{\upsilon \tau }(\upsilon {}_h+ϵ+\upsilon )({\upsilon }_q+ϵ+\upsilon )}}\right]{\mathbb{I}}^{\ast },\\ {\mathbb{P}}^{\ast }={\displaystyle \frac{\alpha {\mathbb{I}}^{\ast }}{\upsilon }}.\end{array}$(23) where, $\begin{array}{rl}{\mathbb{N}}^{\ast }& =\frac{\sigma }{\upsilon }-{\mathbb{I}}^{\ast }\frac{(\upsilon +\tau +\varpi )(\upsilon +\phi +\delta +{\upsilon }_i+\varphi +\alpha )}{\mathrm{\tau \upsilon }}+\frac{(\upsilon +\varphi +\alpha +{\upsilon }_i+\phi +\delta )}{\tau }\\ & +\frac{\varphi }{({\upsilon }_q+\upsilon +ϵ)}+\frac{\phi }{(\upsilon +{\upsilon }_h+\lambda )}\\ & +[\frac{\delta }{\upsilon }+\frac{\varpi ({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )+\mathrm{\lambda \phi }({\upsilon }_q+\upsilon +ϵ)+\mathrm{ϵ\varphi }(\upsilon +{\upsilon }_h+ϵ)}{\mathrm{\upsilon \tau }(\upsilon +{\upsilon }_h+\lambda )({\upsilon }_q+\upsilon +ϵ)}].\end{array}$

6. Basic reproduction number ${\mathcal{R}}_0$

In this section, we use the next generation matrix approach to compute the basic reproduction number (${\mathcal{R}}_0$).

Let the matrices for transition and new infections terms of the model (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) be V and F, respectively. Let $p=(\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{P}{)}^T$ then; $F\left(p\right)=\left(\begin{array}{c}{\displaystyle \frac{\pi \mathbb{S}}{\mathbb{N}}}(\theta \mathbb{L}+\mathbb{I})\\ 0\\ 0\\ 0\end{array}\right),V\left(p\right)=\left(\begin{array}{c}-(\tau +\upsilon +\varpi )\mathbb{L}-\pi \mathbb{S}\mathbb{P}\\ \tau \mathbb{L}-({\upsilon }_i+\mu +\delta +\phi +\varphi +\alpha )\mathbb{I}\\ \varphi \mathbb{I}-({\upsilon }_q+\upsilon +\lambda ){\mathbb{I}}_q\\ \phi \mathbb{I}-({\upsilon }_h+\upsilon +\lambda ){\mathbb{I}}_h\end{array}\right).$ The Jacobian matrix of $V\left(p\right)$ and $F\left(p\right)$ calculated at the disease free equilibrium ($E^0)$ are $\begin{array}{rl}\mathrm{DF}\left(E^0\right)& =\left(\begin{array}{cccc}{\displaystyle \frac{\sigma \text{β}\theta }{\upsilon }}& {\displaystyle \frac{\sigma \text{β}}{\upsilon }}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right),\\ \mathrm{DV}\left(E^0\right)& =(\begin{array}{cc}-(\tau +\varpi +\upsilon )& 0\\ \tau & -(\phi +{\upsilon }_i+\delta +\upsilon +\varphi +\alpha )\\ 0& \varphi \\ 0& \phi \end{array}\\ & \begin{array}{cc}0& 0\\ 0& 0\\ -({\upsilon }_q+\upsilon +ϵ)& 0\\ 0& -({\upsilon }_h+\upsilon +\lambda )\end{array}).\end{array}$ Let, $\mathcal{V}=\mathrm{DV}\left(E^0\right)$ and $\mathcal{F}=\mathrm{DF}\left(E^0\right)$.

Now, $\begin{array}{rl}{\mathcal{V}}^{-1}& =(\begin{array}{c}-{\displaystyle \frac{1}{(\tau +\upsilon +\varpi )}}\\ -{\displaystyle \frac{\tau }{(\sigma +\upsilon +\varpi )({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )}}\\ -{\displaystyle \frac{\mathrm{\tau \varphi }}{(\tau +\upsilon +\varpi )({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )({\upsilon }_q+\upsilon +ϵ)}}\\ -{\displaystyle \frac{\mathrm{\tau \phi }}{(\tau +\upsilon +\varpi )({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )({\upsilon }_h+\upsilon +\lambda )}}\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ -{\displaystyle \frac{1}{({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )}}& 0& 0\\ -{\displaystyle \frac{\varphi }{(\tau +\upsilon +\varpi )({\upsilon }_q+\upsilon +ϵ)}}& -{\displaystyle \frac{1}{({\upsilon }_q+\upsilon +ϵ)}}& \\ -{\displaystyle \frac{\phi }{({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )({\upsilon }_h+\upsilon +\lambda )}}& 0& -{\displaystyle \frac{1}{({\upsilon }_h+\upsilon +ϵ)}}\end{array}).\end{array}$ Therefore, the next generation matrix of model (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) is $\begin{array}{rl}\mathcal{F}{\mathcal{V}}^{-1}& =(\begin{array}{c}-{\displaystyle \frac{\sigma \text{β}\theta ({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )+\mathrm{\sigma \tau \pi }}{\upsilon (\tau +\upsilon +\varpi )({\upsilon }_i+\phi +\delta +\varphi +\phi +\alpha )}}\\ 0\\ 0\\ 0\end{array}\\ & \begin{array}{ccc}-{\displaystyle \frac{\sigma \text{β}}{\upsilon ({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )}}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}),\end{array}$ and the characteristic equation of $\left(\mathcal{F}{\mathcal{V}}^{-1}\right)$ is given by ${\lambda }^3(-a-\lambda )=0,$ where, $a=\frac{\sigma \text{β}\theta (\upsilon {}_i+\upsilon +\delta +\varphi +\alpha +\phi )+\mathrm{\sigma \tau \pi }}{\upsilon (\tau +\upsilon +\varpi )(\upsilon {}_i+\upsilon +\delta +\varphi +\phi +\alpha )}$. Eigenvalues of the above characteristic equation are ${\lambda }_1=-\frac{\mathrm{\sigma \pi \theta }({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )+\mathrm{\sigma \tau \pi }}{\upsilon (\tau +\upsilon +\varpi )({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )},{\lambda }_2=0,{\lambda }_3=0,{\lambda }_4=0.$ Therefore, the largest eigenvalue of $\left(\mathcal{F}{\mathcal{V}}^{-1}\right)$ is $\begin{array}{rl}\rho \left(\mathcal{F}{\mathcal{V}}^{-1}\right)& =|-\frac{\mathrm{\sigma \theta \pi }({\upsilon }_i+\upsilon +\varphi +\delta +\alpha +\phi )+\mathrm{\sigma \tau \pi }}{\upsilon (\tau +\upsilon +\varpi )({\upsilon }_i+\upsilon +\delta +\alpha +\varphi +\phi )}|\\ & =\frac{\sigma \text{β}\theta ({\upsilon }_i+\upsilon +\delta +\phi +\alpha +\varphi )+\mathrm{\sigma \tau \pi }}{\upsilon (\tau +\upsilon +\varpi )({\upsilon }_i+\upsilon +\delta +\alpha +\phi +\varphi )}.\end{array}$ Thus, ${\mathcal{R}}_0$ of the considered model is ${\mathcal{R}}_0=\frac{\text{β}\sigma (\theta ({\upsilon }_i+\upsilon +\delta +\varphi +\alpha +\phi )+\tau )}{\upsilon (\tau +\upsilon +\varpi )({\upsilon }_i+\upsilon +\delta +\phi +\varphi +\alpha )}.$

Figure 2 represents the dynamics of $R_0$ with respect to the transmission rates β and θ.

7. Local stability

From system (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ), we construct the following matrix $J=\left(\begin{array}{ccccccc}{\rho }_1& -{\displaystyle \frac{\text{β}\theta \mathbb{S}}{\mathbb{N}}}& -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}& 0& 0& 0& -\pi \mathbb{S}\\ {\rho }_2& {\rho }_3& {\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}& 0& 0& 0& \pi \mathbb{S}\\ 0& \tau & {\rho }_4& 0& 0& 0& 0\\ 0& 0& \varphi & -({\upsilon }_q+\upsilon +ϵ)& 0& 0& 0\\ 0& 0& \phi & 0& -({\upsilon }_h+\upsilon +\lambda )& 0& 0\\ 0& \varpi & \delta & ϵ& \lambda & -\upsilon & 0\\ 0& 0& \alpha & 0& 0& 0& -{\upsilon }_p\end{array}\right),$ where $\begin{array}{rl}{\rho }_1& =-\frac{\text{β}}{\mathbb{N}}(\mathbb{I}+\theta \mathbb{L})-\upsilon -\pi \mathbb{P},{\rho }_2=\frac{\text{β}}{\mathbb{N}}(\mathbb{I}+\theta \mathbb{L})+\pi \mathbb{P},\\ {\rho }_3& =\frac{\text{β}\theta \mathbb{S}}{\mathbb{N}}-(\upsilon +\tau +\varpi ),\\ {\rho }_4& =-(\upsilon +{\upsilon }_i+\delta +\phi +\varphi +\alpha ).\end{array}$

Figure 2. Dynamics of ${\mathcal{R}}_0$ with respect to transmission rates β and θ. (a) Impact of β on ${\mathcal{R}}_0$. (b) Impact of θ on ${\mathcal{R}}_0$ and (c) Impact of θ and β on ${\mathcal{R}}_0$.

Theorem 7.1

When ${\mathcal{R}}_0<1$, the disease free equilibrium $E^0$ is locally asymptotically stable.

Proof.

The Jacobian matrix of the considered model at $E^0$ is given by $\begin{array}{rl}{J}_{E^0}& =(\begin{array}{cccc}-\upsilon & -\text{β}\theta & -\text{β}& 0\\ 0& \text{β}\theta -(\tau +\upsilon +\varpi )& \text{β}& 0\\ 0& \tau & -{\rho }_3& 0\\ 0& 0& \varphi & -({\upsilon }_q+\upsilon +ϵ)\\ 0& 0& \phi & 0\\ 0& \varpi & \delta & ϵ\\ 0& 0& \alpha & 0\end{array}\\ & \begin{array}{ccc}0& 0& -{\displaystyle \frac{\mathrm{\pi \sigma }}{\upsilon }}\\ 0& 0& {\displaystyle \frac{\mathrm{\pi \sigma }}{\upsilon }}\\ 0& 0& 0\\ 0& 0& 0\\ -({\upsilon }_h+\upsilon +\lambda )& 0& 0\\ n& -\upsilon & 0\\ 0& 0& -{\upsilon }_p\end{array})\end{array}$ We get the eigenvalues $\begin{array}{rl}{\lambda }_1& =-\upsilon ,{\lambda }_2=-({\upsilon }_q+\upsilon +\lambda ),\\ {\lambda }_3& =-({\upsilon }_h+\upsilon +\lambda ),{\lambda }_4=-\upsilon ,\\ {\lambda }_5& =({\mathcal{R}}_0-1)\frac{3\sqrt{2}{s}_5}{{(s_1+s_4+\sqrt{4{\left(s_3\right)}^3+{\left(s_2\right)}^2})}^{1/3}}\\ & -\frac{{(s_4+s_1+\sqrt{4{\left(s_3\right)}^3+{\left(s_2\right)}^2})}^{1/3}}{6(2{)}^{1/3}}-\frac{(m_2-\upsilon {}_p-m_3)}{3},\\ {\lambda }_6& =(\sqrt{3}-1)\frac{{(s_4+s_1+\sqrt{4{\left(s_3\right)}^3+{\left(s_2\right)}^2})}^{1/3}}{3(2{)}^{1/3}}\\ & -\frac{(1+\sqrt{3})s_5}{3(2{)}^{2/3}{(s_1+\sqrt{{(s_2+4{\left(s_3\right)}^3)}^2}+s_4)}^{1/3}}-\frac{(m_2-{\upsilon }_p-m_3)}{3},\\ {\lambda }_7& =-(1+\sqrt{3})\frac{{(s_1+\sqrt{4{\left(s_3\right)}^3+{\left(s_2\right)}^2}+s_4)}^{1/3}}{3(2{)}^{1/3}}\\ & -\frac{(1-\sqrt{3})s_5}{3(2{)}^{2/3}{(s_1+\sqrt{4{\left(s_3\right)}^3+{\left(s_2\right)}^2}+s_4)}^{1/3}}-\frac{(m_2-m_3-{\upsilon }_p)}{3},\end{array}$ with $\begin{array}{rl}{s}_1& =2m_2^3+3m_2^2{m}_3+3m_2^2{\upsilon }_p,\\ s_2& =(2m_2^3+3m_2^2{m}_3+3m_2^2{\mu }_p+9m_2\pi m_1-3m_2^2{m}_3+12m_2{m}_3{\mu }_p-3m_2^2{\mu }_p\\ & {+9\pi m_1{m}_3+18\pi m_1{\upsilon }_p+27m_4{b}_1\alpha -2m_3^3+3m_3^2{\mu }_p-2{\upsilon }_p^3)}^2,\\ s_3& =3(-m_2{m}_3-m_2{\mu }_p-\pi m_1+m_3{\upsilon }_p)-(-m_2+m_3-{\upsilon }_p),\\ s_4& =9m_2\pi m_1-3m_2{m}_3^2+12m_2{m}_3-p-3m_2{\mu }_p^2-9\pi m_1{m}_3+18\pi m_1{\upsilon }_p\\ & +27c{b}_1\alpha -2m_3^3+3m_3^2{\upsilon }_p+3m_3{\upsilon }_p^2-2{\upsilon }_p^3,\\ s_5& =3\left(\right(-m_2{m}_3-m_2{\mu }_p-\text{β}\tau +m_3\upsilon {}_p)-(-m_2+m_3+{\upsilon }_p{)}^2),\end{array}$ where, $m_1=\tau ,m_2=\text{β}\theta -(\upsilon +\tau +\varpi ),m_3=(\upsilon +\upsilon {}_i+\delta +\phi +\varphi +\alpha ),m_4=\frac{\mathrm{\pi \sigma }}{\upsilon }$. The imaginary portion of all eigenvalues is obviously zero. Therefore, $|\arg ({\lambda }_i\left)\right|=\text{β}>\frac{\gamma \text{β}}{2},0<\gamma <1.$ Hence, $E^0$ is LAS.

Theorem 7.2

If ${\mathcal{R}}_0>1$, endemic equilibrium Point $E^{\ast }$ is locally asymptotically stable.

Proof.

$E^{\ast }$ can be describe as (24) $\{\begin{array}{l}{\mathbb{I}}^{\ast }={\displaystyle \frac{\mathrm{\tau \sigma }}{(\tau +\upsilon +\varpi )({\upsilon }_i+\upsilon +\alpha +\phi +\delta +\varphi )}}[1-{\displaystyle \frac{{\mathbb{N}}^{\ast }}{{\mathcal{R}}_0}}],\\ {\mathbb{S}}^{\ast }={\displaystyle \frac{\tau {\mathbb{N}}^{\ast }}{\upsilon {\mathcal{R}}_0}},\\ {\mathbb{L}}^{\ast }={\displaystyle \frac{\sigma ({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )}{(\tau +\upsilon +\varpi )}}[1-{\displaystyle \frac{{\mathbb{N}}^{\ast }}{{\mathcal{R}}_0}}],\\ {\mathbb{I}}_q^{\ast }={\displaystyle \frac{\mathrm{\tau \varphi \sigma }}{({\upsilon }_q+\upsilon +ϵ)(\tau +\upsilon +\varpi )({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )}}[1-{\displaystyle \frac{{\mathbb{N}}^{\ast }}{{\mathcal{R}}_0}}],\\ {\mathbb{I}}_h^{\ast }={\displaystyle \frac{\mathrm{\tau \phi \sigma }}{(\tau +\upsilon +\varpi )({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )({\upsilon }_h+\upsilon +\lambda )}}[1-{\displaystyle \frac{{\mathbb{N}}^{\ast }}{{\mathcal{R}}_0}}],\\ {\mathbb{R}}^{\ast }=[{\displaystyle \frac{\delta }{\upsilon }}+{\displaystyle \frac{\varpi ({\upsilon }_i+\upsilon +\delta +\phi +\varphi +\alpha )+\mathrm{\phi ϵ}({\upsilon }_q+\upsilon +\lambda )+\mathrm{ϵ\varphi }({\upsilon }_h+\lambda +\upsilon )}{\mathrm{\upsilon \tau }(\upsilon {}_h+\lambda +\upsilon )({\upsilon }_q+\upsilon +ϵ)}}]\\ {\displaystyle \frac{\mathrm{\tau \sigma }}{(\tau +\varpi +\upsilon )({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )}}[1-{\displaystyle \frac{{\mathbb{N}}^{\ast }}{{\mathcal{R}}_0}}],\\ {\mathbb{P}}^{\ast }={\displaystyle \frac{\mathrm{\tau \alpha \sigma }}{\upsilon (\tau +\upsilon +\varpi )({\upsilon }_i+\upsilon +\delta +\alpha +\phi +\varphi )}}[1-{\displaystyle \frac{{\mathbb{N}}^{\ast }}{{\mathcal{R}}_0}}].\end{array}.$(24) After the Linearizing around $E^{\ast }$ we obtain $J_{E^{\ast }}=\left(\begin{array}{ccccccc}{f}_1& -{\displaystyle \frac{\theta \text{β}\sigma }{\upsilon {\mathcal{R}}_0}}& -{\displaystyle \frac{\pi \text{β}\sigma }{\upsilon {\mathcal{R}}_0}}& 0& 0& 0& -{\displaystyle \frac{\mathrm{\sigma \pi }{\mathbb{N}}^{\ast }}{\upsilon {\mathcal{R}}_0}}\\ f_2& f_3& {\displaystyle \frac{\sigma \text{β}}{\upsilon {\mathcal{R}}_0}}& 0& 0& 0& {\displaystyle \frac{\mathrm{\sigma \pi }{\mathbb{N}}^{\ast }}{\upsilon {\mathcal{R}}_0}}\\ 0& \tau & -f_3& 0& 0& 0& 0\\ 0& 0& \varphi & -({\upsilon }_q+\upsilon +ϵ)& 0& 0& 0\\ 0& 0& \phi & 0& -({\upsilon }_h+\upsilon +\lambda )& 0& 0\\ 0& \varpi & \delta & \lambda & \lambda & -\upsilon & 0\\ 0& 0& \alpha & 0& 0& 0& -{\upsilon }_p\end{array}\right),$ where $\begin{array}{rl}{f}_1& =-\frac{\upsilon {\mathcal{R}}_0}{{\mathbb{N}}^{\ast }}+\frac{\text{β}\mathrm{\tau \alpha }}{(\theta ({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )+\tau )\text{β}}({\mathbb{N}}^{\ast }-{\mathcal{R}}_0),\\ f_2& =\frac{\upsilon {\mathcal{R}}_0}{{\mathbb{N}}^{\ast }}+\frac{\mathrm{\alpha \sigma \pi }\text{β}}{(\theta ({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )+\tau )}({\mathcal{R}}_0-{\mathbb{N}}^{\ast })-\upsilon ,\\ f_3& =\frac{\mathrm{\theta \pi }\text{β}\sigma }{\upsilon {\mathcal{R}}_0}-(\tau +\upsilon +\varpi ),\\ f_4& =({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha ).\end{array}$ Characteristic equation of the above Jacobian matrix (25) $\begin{array}{rl}P\left(\lambda \right)& =(ϵ+({\upsilon }_q+\mu +ϵ))(ϵ+\upsilon )(ϵ+({\upsilon }_h+\upsilon +\lambda ))\\ & \times (c_0{\lambda }^4+c_1{\lambda }^3+c_2{\lambda }^2+c_3\lambda +c_4)=0,\end{array}$(25) where $\begin{array}{rl}{c}_0& =1,\\ c_1& =-\frac{\mathrm{\theta \pi }\text{β}\sigma }{\upsilon {\mathcal{R}}_0}+(\tau +\upsilon +\varpi )+{\upsilon }_p+f_4-f_1,\\ c_2& =f_1{f}_3+{\upsilon }_p{f}_3-f_4(d_3-f_1-\upsilon {}_p)+{\upsilon }_p{f}_1-\frac{\text{β}\mathrm{\tau \sigma }}{\upsilon {\mathcal{R}}_0}-\frac{\theta \text{β}\sigma }{\upsilon R_0}(f_2-\upsilon ),\\ c_3& =f_1{f}_3(f_4-{\upsilon }_p)+f_1(\frac{\mathrm{\tau \pi \sigma }}{\upsilon {\mathcal{R}}_0}+{\upsilon }_p{f}_4)(\frac{\tau \text{β}\sigma }{\upsilon {\mathcal{R}}_0}+{\upsilon }_p{f}_4)\\ & -\frac{\mathrm{\theta \pi \sigma }}{\upsilon {\mathcal{R}}_0}{f}_4{f}_2+{\upsilon }_p{f}_3{f}_4+\frac{\tau \text{β}\sigma {\upsilon }_p}{\upsilon {\mathcal{R}}_0}-\frac{\sigma \text{β}{\mathbb{N}}^{\ast }\mathrm{\tau \alpha }}{\upsilon {\mathcal{R}}_0},\\ c_4& =\frac{\mathrm{\pi \theta \sigma }{f}_2}{\upsilon {\mathcal{R}}_0}({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha ){\upsilon }_p-\frac{f_1{\upsilon }_p\text{β}\mathrm{\tau \sigma }}{\upsilon {\mathcal{R}}_0}+\frac{f_1\mathrm{\sigma \pi \tau \sigma }{\mathbb{N}}^{\ast }}{\upsilon {\mathcal{R}}_0}\\ & -f_1{f}_3{f}_4{\upsilon }_p+\frac{f_2\mathrm{\sigma \pi \alpha \tau }{\mathbb{N}}^{\ast }}{\upsilon {\mathcal{R}}_0}.\end{array}$ The stability criterion is defined by Rouh-Hurtwiz as follows: the equation's equilibrium point is stable if these inequalities are satisfied $c_0>0,{\mathrm{△}}_1=c_1>0,{\mathrm{△}}_2=c_1{c}_3-c_1{c}_2>0{\mathrm{△}}_3=c_3>0,{\mathrm{△}}_4=c_4>0.$ To begin, use $c_0=1>0$. Now, $c_1>0=\mathrm{△}{}_1$, we can do this if $(\tau +\upsilon +\varpi )+{\upsilon }_p+f_4>\frac{\mathrm{\theta \pi }\text{β}\sigma }{\upsilon {\mathcal{R}}_0}+f_1,$ then $\mathrm{△}{}_1=c_1>0$. For $\mathrm{△}{}_2=c_0{c}_3-c_1{c}_2>0$, it is also clear, if $\begin{array}{rl}& f_1{f}_3(f_4-{\upsilon }_p)+f_1(\frac{\mathrm{\tau \pi \sigma }}{\upsilon {\mathcal{R}}_0}+{\mu }_p{f}_4)(\frac{\tau \text{β}\sigma }{\mu {\mathcal{R}}_0}+{\upsilon }_p{f}_4)-\frac{\theta \text{β}\sigma }{\upsilon {\mathcal{R}}_0}{f}_4{f}_2\\ & +{\upsilon }_p{f}_3{f}_4+\frac{\tau \text{β}\sigma {\upsilon }_p}{\upsilon {\mathcal{R}}_0}-\frac{\mathrm{\sigma \pi }{\mathbb{N}}^{\ast }\mathrm{\alpha \tau }}{\upsilon {\mathcal{R}}_0}>(-\frac{\theta \text{β}\mathrm{\pi \tau }}{\upsilon {\mathcal{R}}_0}+(\alpha +\mu +\varpi )+{\upsilon }_p+f_4-f_1)\\ & \times (d_1{f}_3+{\mu }_p{f}_3-f_4(f_3-f_1-{\upsilon }_p)+{\upsilon }_p{f}_1-\frac{\text{β}\mathrm{\tau \sigma }}{\upsilon {\mathcal{R}}_0}-\frac{\theta \text{β}\sigma }{\upsilon {\mathcal{R}}_0}(f_2-\upsilon )),\end{array}$ then $\mathrm{△}{}_2>0$.

Now, for $\mathrm{△}{}_3>0$, if $\begin{array}{rl}& f_1{f}_3(f_4-{\upsilon }_p)+f_1(\frac{\tau \text{β}\sigma }{\upsilon {\mathcal{R}}_0}+{\upsilon }_p{f}_4)(\frac{\tau \text{β}\sigma }{\upsilon {\mathcal{R}}_0}+{\upsilon }_p{f}_4)\\ & +{\mu }_p{f}_3{f}_4+\frac{\tau \text{β}\sigma {\upsilon }_p}{\upsilon R_0}>\frac{\text{β}\mathrm{\theta \sigma }}{\upsilon {\mathcal{R}}_0}{f}_4{f}_2+\frac{\mathrm{\sigma \pi }{\mathbb{N}}^{\ast }\mathrm{\alpha \tau }}{\upsilon {\mathcal{R}}_0},\end{array}$ then $\mathrm{△}{}_3>0$.

Lastly, $\mathrm{△}{}_4>0$, $\begin{array}{rl}& \frac{\theta \text{β}\sigma f_2}{\upsilon {\mathcal{R}}_0}({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha ){\upsilon }_p+\frac{f_1\mathrm{\alpha \pi \tau \sigma }{\mathbb{N}}^{\ast }}{\upsilon {\mathcal{R}}_0}\\ & +\frac{f_2\mathrm{\sigma \pi \alpha \tau }{\mathbb{N}}^{\ast }}{\upsilon {\mathbb{R}}_0}>\frac{f_1{\upsilon }_p\sigma \text{β}\tau }{\upsilon {\mathcal{R}}_0}+f_1{f}_3{f}_4{\upsilon }_p,\end{array}$ then $\mathrm{△}{}_4>0$.

The Routh–Hurwitz stability criteria are satisfied if ${\mathcal{R}}_0>1$ because all of the conditions are dependent on ${\mathcal{R}}_0$. Hence,the endemic equilibrium is asymptomatically stable when ${\mathcal{R}}_0>1$.

8. Global stability

Theorem 8.1

When ${\mathcal{R}}_0<1$, the disease free equilibrium point $E^0$ is globally asymptomatically stable.

Proof.

To derive the Lyapunov candidate function for fractional order as in [16], consider the family of quadratic Lyapunov function $L(f_1,f_2,\dots ,f_n)=\sum _{i=1}^n\frac{c_i}{2}{\left(f_i\right(t)-f^0)}^2,$ where $c_i=t-1-\ln t$. We defined the Lyapunov function as (26) $\begin{array}{rl}L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})& =\frac{1}{2}{(\mathbb{S}-{\mathbb{S}}^0)}^2+\frac{1}{2}{(\mathbb{L}-{\mathbb{L}}^0)}^2+\frac{1}{2}{(\mathbb{I}-{\mathbb{I}}^0)}^2\\ & +\frac{1}{2}{({\mathbb{I}}_q-{\mathbb{I}}_q^0)}^2+\frac{1}{2}{({\mathbb{I}}_h-{\mathbb{I}}_h^0)}^2+\frac{1}{2}{(\mathbb{R}-{\mathbb{R}}^0)}^2\\ & +\frac{1}{2}{(\mathbb{P}-{\mathbb{P}}^0)}^2.\end{array}$(26) Using the linearity property of the Caputo operator: (27) $\begin{array}{rl}{}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})& =\frac{1}{2}[{}_0^c{\mathcal{D}}_t^{\gamma }{(\mathbb{S}-{\mathbb{S}}^0)}^2+{}_0^c{\mathcal{D}}_t^{\gamma }{(\mathbb{L}-{\mathbb{L}}^0)}^2\\ & {+}_0^c{\mathcal{D}}_t^{\gamma }{(\mathbb{I}-{\mathbb{I}}^0)}^2{+}_0^c{\mathcal{D}}_t^{\gamma }{({\mathbb{I}}_q-{\mathbb{I}}_q^0)}^2\\ & {+}_0^c{\mathcal{D}}_t^{\gamma }{({\mathbb{I}}_h-{\mathbb{I}}_h^0)}^2{+}_0^c{\mathcal{D}}_t^{\gamma }{(\mathbb{R}-{\mathbb{R}}^0)}^2\\ & {+}_0^c{\mathcal{D}}_t^{\gamma }{(\mathbb{P}-{\mathbb{P}}^0)}^2].\end{array}$(27) Using Lemma 2 [15], we get (28) $\begin{array}{rl}{}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})& \le {}_0^c{\mathcal{D}}_t^{\gamma }(\mathbb{S}-{\mathbb{S}}^0)+{}_0^c{D}_t^{\gamma }(\mathbb{L}-{\mathbb{L}}^0)\\ & {+}_0^c{\mathcal{D}}_t^{\gamma }(\mathbb{I}-{\mathbb{I}}^0){+}_0^c{D}_t^{\gamma }({\mathbb{I}}_q-{\mathbb{I}}_q^0)+{}_0^c{\mathcal{D}}_t^{\gamma }({\mathbb{I}}_h-{\mathbb{I}}_h^0)\\ & {+}_0^c{\mathcal{D}}_t^{\gamma }(\mathbb{R}-{\mathbb{R}}^0){+}_0^c{D}_t^{\gamma }(\mathbb{P}-{\mathbb{P}}^0),\\ & =\sigma +\upsilon (\mathbb{S}-{\mathbb{S}}^0+\mathbb{L}-{\mathbb{L}}^0+\mathbb{I}-{\mathbb{I}}^0+{\mathbb{I}}_q-{\mathbb{I}}_q^0\\ & +{\mathbb{I}}_h-{\mathbb{I}}_h^0+\mathbb{R}-{\mathbb{R}}^0+\mathbb{P}-{\mathbb{P}}^0)-{\upsilon }_i(\mathbb{I}-{\mathbb{I}}^0)\\ & -{\upsilon }_q({\mathbb{I}}_q-{\mathbb{I}}_q^0)-{\upsilon }_h({\mathbb{I}}_h-{\mathbb{I}}_h^0)-{\upsilon }_v(\mathbb{P}-{\mathbb{P}}^0).\end{array}$(28) Thus (29) $\begin{array}{rl}{}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})& \le \sigma -\upsilon ({\mathbb{N}}^0-\mathbb{N})-{\upsilon }_i(\mathbb{I}-{\mathbb{I}}^0)\\ & -{\upsilon }_q({\mathbb{I}}_q-{\mathbb{I}}_q^0)-{\upsilon }_h({\mathbb{I}}_h-{\mathbb{I}}_h^0)-{\upsilon }_v(\mathbb{P}-{\mathbb{P}}^0).\end{array}$(29) Substituting the point $E^0$ into (29(29) $\begin{array}{rl}{}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})& \le \sigma -\upsilon ({\mathbb{N}}^0-\mathbb{N})-{\upsilon }_i(\mathbb{I}-{\mathbb{I}}^0)\\ & -{\upsilon }_q({\mathbb{I}}_q-{\mathbb{I}}_q^0)-{\upsilon }_h({\mathbb{I}}_h-{\mathbb{I}}_h^0)-{\upsilon }_v(\mathbb{P}-{\mathbb{P}}^0).\end{array}$(29) ), we get ${}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})\le -\upsilon \mathbb{N}-{\upsilon }_i\mathbb{I}-{\upsilon }_h{\mathbb{I}}_h-{\upsilon }_q{\mathbb{I}}_q-{\upsilon }_v\mathbb{P}.$ Therefore (30) ${}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})\le -Y\left(x\right(t\left)\right),$(30) where $\upsilon \mathbb{N}+{\upsilon }_i\mathbb{I}+{\upsilon }_h{\mathbb{I}}_h+{\upsilon }_q{\mathbb{I}}_q=Y\left(x\right(t\left)\right).$ Thus, ${}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})\le 0$. The equilibrium point $E_0$ is thus globally asymptotically stable, according to Theorem 2.3.

Theorem 8.2

When ${\mathcal{R}}_0>1$, the point $E^{\ast }$ is globally asymptomatically stable.

Proof.

To derive the Lyapunov candidate function for fractional order as in [16], consider the family of quadratic Lyapunov function $L(f_1,f_2,\dots ,f_n)=\sum _{i=1}^n\frac{c_i}{2}{\left(f_i\right(t)-f^0)}^2,$ where $c_i=t-1-\ln t$. We defined the Lyapunov function as (31) $\begin{array}{rl}L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})& =\frac{1}{2}{(\mathbb{S}-{\mathbb{S}}^{\ast })}^2+\frac{1}{2}{(\mathbb{L}-{\mathbb{L}}^{\ast })}^2+\frac{1}{2}{(\mathbb{I}-{\mathbb{I}}^{\ast })}^2\\ & +\frac{1}{2}{({\mathbb{I}}_q-{\mathbb{I}}_q^{\ast })}^2+\frac{1}{2}{({\mathbb{I}}_h-{\mathbb{I}}_h^{\ast })}^2+\frac{1}{2}{(\mathbb{R}-{\mathbb{R}}^{\ast })}^2\\ & +\frac{1}{2}{(\mathbb{P}-{\mathbb{P}}^{\ast })}^2.\end{array}$(31) Using the linearity property Caputo operator, we get (32) $\begin{array}{rl}{}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})& =\frac{1}{2}[{}_0^c{\mathcal{D}}_t^{\gamma }{(\mathbb{S}-{\mathbb{S}}^{\ast })}^2+{}_0^c{\mathcal{D}}_t^{\gamma }{(\mathbb{L}-{\mathbb{L}}^{\ast })}^2\\ & {+}_0^c{\mathcal{D}}_t^{\gamma }{(\mathbb{I}-{\mathbb{I}}^{\ast })}^2{+}_0^c{\mathcal{D}}_t^{\gamma }{({\mathbb{I}}_q-{\mathbb{I}}_q^{\ast })}^2\\ & {+}_0^c{\mathcal{D}}_t^{\gamma }{({\mathbb{I}}_h-{\mathbb{I}}_h^{\ast })}^2{+}_0^c{\mathcal{D}}_t^{\gamma }{(\mathbb{R}-{\mathbb{R}}^{\ast })}^2\\ & {+}_0^c{\mathcal{D}}_t^{\gamma }{(\mathbb{P}-{\mathbb{P}}^{\ast })}^2].\end{array}$(32) Using Lemma 2 [15], we get (33) $\begin{array}{rl}{}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})& \le {}_0^c{\mathcal{D}}_t^{\gamma }(\mathbb{S}-{\mathbb{S}}^{\ast })+{}_0^c{\mathcal{D}}_t^{\gamma }(\mathbb{L}-{\mathbb{L}}^{\ast })\\ & {+}_0^c{\mathcal{D}}_t^{\gamma }(\mathbb{I}-{\mathbb{I}}^{\ast }){+}_0^c{\mathcal{D}}_t^{\gamma }({\mathbb{I}}_q-{\mathbb{I}}_q^{\ast })+{}_0^c{\mathcal{D}}_t^{\gamma }({\mathbb{I}}_h-{\mathbb{I}}_h^{\ast })\\ & {+}_0^c{\mathcal{D}}_t^{\gamma }(\mathbb{R}-{\mathbb{R}}^{\ast }){+}_0^c{\mathcal{D}}_t^{\gamma }(\mathbb{P}-{\mathbb{P}}^{\ast }),\\ & =\sigma +\upsilon (\mathbb{S}-{\mathbb{S}}^{\ast }+\mathbb{L}-{\mathbb{L}}^{\ast }+\mathbb{I}-{\mathbb{I}}^{\ast }+{\mathbb{I}}_q-{\mathbb{I}}_q^{\ast }\\ & +{\mathbb{I}}_h-{\mathbb{I}}_h^{\ast }+\mathbb{R}-{\mathbb{R}}^{\ast }+\mathbb{P}-{\mathbb{P}}^{\ast })-{\upsilon }_i(\mathbb{I}-{\mathbb{I}}^{\ast })\\ & -{\upsilon }_q({\mathbb{I}}_q-{\mathbb{I}}_q^{\ast })-{\upsilon }_h({\mathbb{I}}_h-{\mathbb{I}}_h^{\ast })-{\upsilon }_v(\mathbb{P}-{\mathbb{P}}^0).\end{array}$(33) (34) $\begin{array}{rl}{}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})& \le \sigma +\upsilon ({\mathbb{N}}^{\ast }-\mathbb{N})-{\upsilon }_i(\mathbb{I}-{\mathbb{I}}^{\ast })\\ & -{\upsilon }_q({\mathbb{I}}_q-{\mathbb{I}}_q^{\ast })-{\upsilon }_h({\mathbb{I}}_h-{\mathbb{I}}_h^{\ast })-{\upsilon }_v(\mathbb{P}-{\mathbb{P}}^0).\end{array}$(34) At the endemic equilibrium, (34(34) $\begin{array}{rl}{}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})& \le \sigma +\upsilon ({\mathbb{N}}^{\ast }-\mathbb{N})-{\upsilon }_i(\mathbb{I}-{\mathbb{I}}^{\ast })\\ & -{\upsilon }_q({\mathbb{I}}_q-{\mathbb{I}}_q^{\ast })-{\upsilon }_h({\mathbb{I}}_h-{\mathbb{I}}_h^{\ast })-{\upsilon }_v(\mathbb{P}-{\mathbb{P}}^0).\end{array}$(34) ) becomes $\begin{array}{rl}{}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})& \le \sigma -\upsilon \mathbb{N}-{\upsilon }_i\mathbb{I}-{\upsilon }_h{\mathbb{I}}_h-{\upsilon }_q{\mathbb{I}}_q-{\upsilon }_v\mathbb{P}\\ & -\frac{\mathrm{\varphi \tau \sigma }{\upsilon }_q}{(\upsilon +{\upsilon }_q+ϵ){\delta }_1}({\mathcal{R}}_0-\frac{\mathbb{N}(t{)}^{\ast }}{{\mathcal{R}}_0})\\ & -\frac{\mathrm{\tau \sigma }{\upsilon }_i}{{\delta }_1}({\mathcal{R}}_0-\frac{{\mathbb{N}}^{\ast }}{{\mathcal{R}}_0})\\ & -\frac{\mathrm{\phi \tau \sigma }}{{\delta }_1({\upsilon }_h+\upsilon +\lambda )}({\mathcal{R}}_0-\frac{{\mathbb{N}}^{\ast }}{{\mathcal{R}}_0})\\ & -\upsilon {\mathbb{N}}^{\ast },\end{array}$ where, $\delta {}_1=\frac{1}{(\upsilon +\tau +\varpi )(\upsilon +\upsilon {}_i+\delta +\phi +\varphi +\alpha )}$, $\begin{array}{rl}& {}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{M},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})\\ & \le -\sigma -\upsilon \mathbb{N}-{\upsilon }_i\mathbb{I}-{\upsilon }_h{\mathbb{I}}_h-{\upsilon }_q{\mathbb{I}}_q-{\upsilon }_v\mathbb{P}\\ & -({\mathcal{R}}_0-\frac{{\mathbb{N}}^{\ast }}{{\mathcal{R}}_0})\left[\frac{\mathrm{\tau \varphi \sigma }{\upsilon }_q+\mathrm{\sigma \tau }{\upsilon }_i({\upsilon }_h+\upsilon +\lambda )+\mathrm{\tau \phi \tau }}{({\upsilon }_q+\upsilon +ϵ){\delta }_1}\right]\\ & -\upsilon {\mathbb{N}}^{\ast }.\\ & {}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})\\ & \le -\tau -\upsilon \mathbb{N}-{\upsilon }_i\mathbb{I}-{\upsilon }_h{\mathbb{I}}_h-{\upsilon }_q{\mathbb{I}}_q-{\upsilon }_v\mathbb{P}\\ & -{\mathcal{R}}_0\left[\frac{\mathrm{\tau \varphi \sigma }{\upsilon }_q+\mathrm{\sigma \tau }{\upsilon }_i({\upsilon }_h+\lambda +\upsilon )+\mathrm{\tau \phi \sigma }}{({\upsilon }_q+\upsilon +ϵ){\delta }_1}\right]\\ & +{\mathbb{N}}^{\ast }[\frac{\mathrm{\tau \varphi \sigma }{\upsilon }_q+\mathrm{\sigma \tau }{\upsilon }_i({\upsilon }_h+\lambda +\upsilon )+\mathrm{\tau \phi \sigma }}{{\mathcal{R}}_0({\upsilon }_q+\upsilon +ϵ){\delta }_1}-\upsilon ].\end{array}$ The positive equilibrium is stable if ${\mathcal{R}}_0>1$, as we previously proved.

Now, we look the relations $\begin{array}{rl}& \text{β}\sigma (\theta (\upsilon +{\upsilon }_i+\delta +\varphi +\phi +\alpha )+\tau )\\ & >\upsilon (\tau +\upsilon +\varpi )({\upsilon }_i+\upsilon +\delta +\phi +\varphi +\alpha )\\ & .\frac{\sigma \text{β}(\theta ({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )+\tau )}{\upsilon (\tau +\upsilon +\varpi )({\upsilon }_i+\upsilon +\delta +\phi +\varphi +\alpha )}>1.\end{array}$ Back substituting the above mentioned equations into (34(34) $\begin{array}{rl}{}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})& \le \sigma +\upsilon ({\mathbb{N}}^{\ast }-\mathbb{N})-{\upsilon }_i(\mathbb{I}-{\mathbb{I}}^{\ast })\\ & -{\upsilon }_q({\mathbb{I}}_q-{\mathbb{I}}_q^{\ast })-{\upsilon }_h({\mathbb{I}}_h-{\mathbb{I}}_h^{\ast })-{\upsilon }_v(\mathbb{P}-{\mathbb{P}}^0).\end{array}$(34) ), we get $\begin{array}{rl}& {}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})\\ & \le -\sigma -\upsilon \mathbb{N}-{\upsilon }_q{\mathbb{I}}_q-{\upsilon }_i\mathbb{I}-{\upsilon }_h{\mathbb{I}}_h-{\upsilon }_v\mathbb{P}\\ & -{\mathcal{R}}_0\left[\frac{\mathrm{\tau \varphi \sigma }{\upsilon }_q+\mathrm{\sigma \tau }{\upsilon }_i({\upsilon }_h+\lambda +\upsilon )+\mathrm{\phi \sigma \tau }}{({\upsilon }_q+\upsilon +ϵ){\delta }_1}\right]\\ & -{\mathbb{N}}^{\ast }[1-{\mathcal{R}}_0]\end{array}$ Therefore (35) ${}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})\le -\mathcal{W}\left(x\right(t\left)\right),$(35) where $\begin{array}{rl}\mathcal{W}\left(x\right(t\left)\right)& =-\upsilon \mathbb{N}-{\upsilon }_q{\mathbb{I}}_q-{\upsilon }_i\mathbb{I}-{\upsilon }_h{\mathbb{I}}_h-{\upsilon }_v\mathbb{P}\\ & -{\mathcal{R}}_0\left[\frac{\mathrm{\tau \varphi \sigma }{\upsilon }_q+\mathrm{\sigma \tau }{\upsilon }_i(\upsilon +{\upsilon }_h+\lambda )+\mathrm{\phi \tau \sigma }}{({\upsilon }_q+\upsilon +\lambda )(\sigma +\upsilon +\varpi )({\upsilon }_i+\upsilon +\delta +\varphi +\phi +\alpha )}\right]\\ & -[1-{\mathcal{R}}_0]\end{array}$ Thus, ${}_0^c{\mathcal{D}}_t^{\gamma }L(\mathbb{S},\mathbb{L},\mathbb{I},{\mathbb{I}}_q,{\mathbb{I}}_h,\mathbb{R},\mathbb{P})\le 0$ for ${\mathcal{R}}_0>1$ According to LaSalle's invariance principle, the endemic equilibrium $E^{\ast }$ is globally asymptotically stable.

9. Graphical representation

In this part, we looked at how each variable in the system (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) behaved when applied to the Caputo fractional operators. To build graphical profiles for each variable, values of some of the parameters were obtained from Table 2 and variations in the fractional order γ were employed.

Table 2. Parameters and their corresponding values.

Variables and parameters	Values	References
${\mathbb{S}}_0$	$1.0\times {10}^{-7}$	[34]
${\mathbb{L}}_0$	$4.0\times {10}^3$	[34]
${\mathbb{I}}_0$	$1.0\times {10}^3$	[34]
${\mathbb{I}}_q\left(0\right)$	$1.0\times {10}^3$	Assumed
${\mathbb{I}}_h\left(0\right)$	$1.0\times {10}^3$	Assumed
${\mathbb{R}}_0$	$1.0\times {10}^2$	[34]
${\mathbb{P}}_0$	$1.0\times {10}^2$	Assumed
σ	0.5	[34]
υ	0.1791	[34]
$\upsilon {}_i$	0.5002	[34]
$\upsilon {}_q$	0.0329	Fitted
$\upsilon {}_h$	0.0373	Fitted
β	0.0082	Fitted
π	0.006	[34]
τ	$\frac{1}{12}$	[34]
ϕ	0.75	Fitted
ϵ	0.097	Fitted
λ	0.325	Fitted
δ	$\frac{1}{30}$	[34]
α	0.399	[34]
θ	0.42	[34]
ϖ	0.006	Fitted
$\upsilon {}_p$	0.00453	Fitted

Figures 3–9 depict the dynamical behaviour of each compartment of the system (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) for various values of the fractional parameter γ. Figure 4 shows that the asymptomatic population decrease when the quantity of pathogen decreases. Figure 3 describes the susceptible population's diminishing behaviour as γ increases. Similar results are obtained for hospitalized and quarantined classes shown in Figures 7 and 6, respectively. Using the Caputo operator, Figure 5 depicts the pandemic's long-term behaviour. At the conclusion of 50 days, the overall number of mild individuals seem to be decreasing. The significance of the order parameter γ for the pathogen class in the model is seen in Figure 9. The decrease in the values of γ has resulted in a considerable reduction in the number of infected people, as shown in Figure 7, and an increase in the numbers of infected people, as depicted in Figure 7.

Figure 3. Dynamics of Susceptible individuals for various values of fractional order $\gamma =0.40,0.50,0.60,0.70,0.80$ and $0.90$.

Figure 4. Dynamics of latent individuals for various values of fractional order $\gamma =0.40,0.50,0.60,0.70,0.80$ and $0.90$.

Figure 5. Dynamics of infected individuals for various values of fractional order $\gamma =0.40,0.50,0.60,0.70,0.80$ and $0.90$.

Figure 6. Dynamics of quarantined individuals for various values of fractional order $\gamma =0.40,0.50,0.60,0.70,0.80$ and $0.90$.

Figure 7. Dynamics of hospitalized individuals for various values of fractional order $\gamma =0.40,0.50,0.60,0.70,0.80$ and $0.90$.

Figure 8. Dynamics of recovered individuals for various values of fractional order $\gamma =0.60,0.70,0.80$ and $0.90$.

Figure 9. Dynamics of pathogen individuals for various values of fractional order $\gamma =0.40,0.50,0.60,0.70,0.80$ and $0.90$.

We have shown illustrations for seven different values of γ as can be found in plots 10–16. The purpose of these graphical results has to explore the effects of control measures like hospitalization, quarantine, contact rate, and quantity of pathogens. These results show that taking fractional orders has a significant effect on the dynamics of the disease. In Figure 10, We analyse the dynamic behaviour of the model variable for the fractional order $\gamma =0.90$. It can be seen in Figure 11 that when the number of infected individuals decreases, the number of recovered individuals increases. This observation is also evident in Figures 12 – 15, where a decrease in the number of infected people is associated with a significant decrease in the number of mild, hospitalized, and quarantined individuals. The results presented in Figure 16 suggest that as the order γ decreases, the graph of infected individuals also decreases significantly for the same value of the contact rate β. Furthermore, peaks in the pandemic can be observed to decrease with an increase in the hospitalization and quarantine rates. Interestingly, a similar effect is observed for the same values of β and γ. This analysis demonstrates that an increase in quarantine and hospitalization rates could significantly reduce the epidemic peak and decrease the total number of infected cases, as depicted in Figure 16.

10. Conclusion

The ongoing COVID-19 pandemic, lasting for more than two years now, remains one of the most significant challenges humanity has faced in combating coronavirus pathogens. The causative agent, SARS-CoV-2, possesses distinctive biological and transmission properties compared to its predecessors, SARS-CoV and MERS-CoV. Despite the implementation of numerous control strategies worldwide, the number of infected cases continues to rise, necessitating a deeper understanding of COVID-19's dynamics and the development of effective intervention measures. To gain insights into the transmission dynamics of COVID-19, a mathematical model called ‘$\mathbb{S}\mathbb{L}\mathbb{I}{\mathbb{I}}_q{\mathbb{I}}_h\mathbb{R}\mathbb{P}$’ is formulated and analysed in this study. This mathematical model aims to provide a comprehensive representation of the interaction between susceptible ($\mathbb{S}$), latent ($\mathbb{L}$), infected ($\mathbb{I}$), quarantined infected (${\mathbb{I}}_q$), hospitalized infected (${\mathbb{I}}_h$), recovered ($\mathbb{R}$), and deceased ($\mathbb{P}$) individuals within the population. To ensure the validity and reliability of the model, its uniqueness and existence are established through the application of the Karsnoselskil's fixed-point theorem. Additionally, the basic reproduction number, denoted as ${\mathcal{R}}_0$, is calculated using the next-generation matrices approach. This parameter provides crucial information about the disease's potential for spread and helps in determining the necessary control measures. The analysis of the model reveals both disease-free and endemic equilibrium points. It is proven that these equilibrium points are locally and globally asymptotically stable, indicating the stability of the system under certain conditions. Such insights are vital for understanding the long-term behavior of the COVID-19 transmission dynamics and devising appropriate strategies for disease control and containment. To visually represent the model's dynamics and gain a better understanding of the population's changes over time, graphical representations are employed. Different scenarios are considered by varying the value of the parameter γ, which is a measure of the recovery rate. These graphical representations enable researchers and policymakers to observe how the populations of susceptible, infected, quarantined, hospitalized, recovered, and deceased individuals evolve under different conditions of recovery rate. In conclusion, this study presents a comprehensive mathematical model for understanding the transmission dynamics of COVID-19. The analysis of this model, including the calculation of ${\mathcal{R}}_0$ and stability analysis of equilibrium points, contributes to a deeper understanding of the pandemic's behavior. The graphical representations further enhance the comprehension of population dynamics and can aid in devising effective strategies to control the spread of COVID-19 and mitigate its impact on communities worldwide.

Figure 10. Numerical solution of system (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) with fractional order $\gamma =0.90$.

Figure 11. Numerical solution of system (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) with fractional order $\gamma =0.80$.

Figure 12. Numerical solution of system (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) with fractional order $\gamma =0.70$.

Figure 13. Numerical solution of the system with fractional order $\gamma =0.60$.

Figure 14. Numerical solution of system (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) )with fractional order $\gamma =0.50$.

Figure 15. Numerical solution of system (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) ) with fractional order $\gamma =0.40$.

Figure 16. Numerical solution of system (1(1) $\{\begin{array}{l}{}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{S}\right]=\sigma -{\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-\upsilon \mathbb{S}-\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{L}\right]={\displaystyle \frac{\text{β}\mathbb{S}}{\mathbb{N}}}(\mathbb{I}+\theta \mathbb{L})-(\upsilon +\varpi +\tau )\mathbb{L}+\pi \mathbb{S}\mathbb{P},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{I}\right]=\tau \mathbb{L}-(\upsilon +\delta +{\upsilon }_i+\phi +\alpha +\varphi )\mathbb{I},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_q\right]=\varphi \mathbb{I}-(\upsilon +ϵ+{\upsilon }_q){\mathbb{I}}_q,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[{\mathbb{I}}_h\right]=\phi \mathbb{I}-(\upsilon +\lambda +{\upsilon }_h){\mathbb{I}}_h,\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{R}\right]=\delta \mathbb{I}+\varpi \mathbb{L}+\lambda {\mathbb{I}}_h+ϵ{\mathbb{I}}_q-\upsilon \mathbb{R},\\ {}_0^c{\mathcal{D}}_t^{\gamma }\left[\mathbb{P}\right]=\alpha \mathbb{I}-{\upsilon }_p\mathbb{P},\end{array}$(1) )with fractional order $\gamma =0.30$.

Authors' contributions

Authors approved the final manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Data sharing is not applicable to this paper as no datasets were generated or analysed during the current study.

Additional information

Funding

This research was funded by Deanship of Scientific Research at King Khalid University, Abha, Saudi Arabia, Project under Grant Number RGP.2/27/44.