Quantitative analysis of a fractional order of the $SEI_{c}\, I_{\eta} VR$ epidemic model with vaccination strategy

1.   Introduction

Mathematical modeling plays an important role in understanding a lot of phenomena and situations in the real world. It can be applied to describe the relationships among the components of a problem by using mathematical tools. In addition, it is a useful technique that is applied to introduce predictions for real problems. Moreover, it deals with a different field of sciences, for instance, see ^{[1,2,3,4,5,6,7]}. One of the mathematical tools applied to analyze and discuss the problems is differential equations. There is more than one operator of differential equation that helps study the problem and gives accurate representations for the problems.

Mathematical modeling is commonly used in studying the outbreak of infectious diseases. Recently, many articles were published to investigate the transmission of various diseases among the people such as influenza ^[8,9], Zika virus ^[10], Covid-19 ^[11,12], and measles ^[13]. In the field of mathematical modeling of infectious diseases, the researchers are interested in presenting the qualitative analysis of their proposed model. They focus on finding the conditions that reduce the pandemic by illustrating the stability theory. Therefore, they employ the ordinary differential equations (ODEs), partial differential equations (PDEs), discrete differential equations, and fractional order differential equations (FDEs) in their study. The earliest modeling of infectious disease considered the population of the outbreak has only three groups: susceptible group $S$, infected group $I$, and recovered group $R$ ^[14]. After that, the scientists extended a classic model $S\, I\, R$ into $S\, E\, I\, R$ to incorporate the exposed individual ^[15,16,17]. A lot of modifications for the classical epidemic models were introduced by including other categories of population components such as: vaccination class, treatment control, or new factors that affect on the transmission of the infectious disease ^[18,19,20].

Vaccine is an appropriate way and it is important to reduce the spread of infectious diseases. Therefore, the researchers considered the impact of vaccination in their models, and many of the published articles describe the transmission of infectious diseases through the population by assuming the individuals who are vaccinated as a component of the population ^{[21,22,23,24,25]}.

One of the significant examples of infectious diseases is Coronavirus. Coronavirus disease (Covid-19) was discovered in 2019 and then spread worldwide, leading to a continuing pandemic outbreak worldwide ^[26,27]. Recently, many published articles studied the outbreak of Covid-19 by considering more than one class of infected people, since a person can be infected by Covid-19 without displaying any symptoms. Also, some of these articles show the role of a vaccine in fighting the diffusion of Covid-19 ^[28,29,30]. Therefore, it is important to develop a mathematical model to represent the spread of infectious diseases by including the different types of infected people and vaccinated people at the same time. Accordingly, this work aims to propose mathematical modeling in which the population involves two types of infected groups and vaccinated individuals, and it also considers the fractional differential equations as an operator for this system. The fractional order of differential equations was used in many biological systems because the memory feature, which is the ability to include all previous conditions to explain the state condition, and the non-locality effect distinguishes the fractional derivative, which means it is valuable for the pandemic models depend on the history ^[31,32,33]. Furthermore, the fractional order of differential equations is more similar to real-life problems, and it expresses the whole time domain for physical stages, not as in the ODEs systems that are related to the local properties of a certain position ^[32]. Hence, fractional calculus investigates the arbitrary-order derivative, that is, it is a generalization of ODEs ^[34].

In ^[34], Sun et al. constructed an SEQIR epidemic model with saturated incidence and vaccination. They investigated the stability analysis of equilibrium points by presenting different theorems. Soulaimani and Kaddar ^[35] introduced a study of an optimal control of fractional order SEIR model with general incidence and vaccination. The SVEIR epidemic model was considered and verified by Nabti and Ghanbari ^[36], and they proved properties of their model such as existence, boundedness, and global stability of equilibria. In ^[37], the authors proposed a new fractional infectious disease model under the non-singular Mittage-Leffler derivative. They assumed the infectious group has been divided into two classes: acute and chronically infectious people. In ^[38], Ali et al. investigated the transmission dynamics of a fractional-order mathematical model of COVID-19 including susceptible, exposed, asymptomatic infected, symptomatic infected, and recovered.

However, these previous models neglected the effect of vaccination on some infectious individual cases. Hence, this work distributes a new modification of infection dynamics that contains two types of infected individuals and people who are vaccinated. It is organized as follows: In Section 2, we introduce the description of our model. Some definitions of fractional derivatives, existence, positivity, and boundedness of the proposed model are shown in Section 3. In Section 4, we prove the global stability of the equilibria. The numerical simulations are given in Section 5, to validate the analysis results proved in previous sections.

2.   Model formulation

This paper is interested in modifying the SEIR epidemic model. The proposed model consists of the following components: the susceptible individuals are denoted by $\mathrm{S}(t)$, the exposed individuals are given by $\mathrm{E}(t)$, $\mathrm{I_{c}}(t)$ and $\mathrm{I_{\eta}}(t)$ are defined as the asymptomatic infected and symptomatic infected respectively, $\mathrm{V}(t)$ is described the individuals that are vaccinated, and $\mathrm{R}(t)$ are recovered population.

2.1.   System description

Our proposed model is described as

with initial conditions

Here, the total population is $\mathbb{N}(t) = \mathrm{S}(t)+\mathrm{E}(t)+\mathrm{I_{c}}(t)+\mathrm{I_{\eta}}(t)+\mathrm{V}(t)+\mathrm{R}(t)$, where $D^{\vartheta}$ is fractional derivative in the Caputo sense and $\vartheta$ is a parameter that describes the order of the fractional time-derivative with $0 < \vartheta\leq1$. Each of the proposed model's parameters is assumed to be positive and constant, and they are shown in Table 1. Figure 1 describes all elements of system (2.1).

3.   Mathematical preliminaries and basic properties

In our work, we will apply the Caputo-derivative because it has important features such as nonlocal and nonsingular exponential kernel ^[30]. Consequently, in this section, we introduce some of the related definitions. Also, we will investigate the solution's following characteristics: existence of the solution, uniqueness, positivity, and boundedness.

Definition 3.1. ^[39] The Caputo derivative of order $\vartheta\, > \, 0$ for any function $\chi\, \in\, \mathbb{C}^{n}([t_{0}, \infty), \mathbb{R})$ is given as

where $t\, \geq \, t_{0}$, ${\Gamma}$ is Gamma function, and $n\, \in\, \mathbb{Z}^+$, such that $\vartheta\, \in\, (n-1, n)$.

Lemma 3.1. ^[40] Consider the following fractional system

where $\vartheta \in (0, 1], \Upsilon:[t_{0}, \infty)\times {\ddot{\Xi}}\longrightarrow\mathbb{R}^{6}$. If $\Upsilon$ satisfied the Lipschitz condition, then the system (2.1) has a unique solution on $[t_{0}, \infty)$.

3.1.   Existence, positivity, and boundedness of solutions

Assume that $\ddot{\Xi} = \{\mathrm{S}, \mathrm{E}, \mathrm{I_{c}}, \mathrm{I_{\eta}}, \mathrm{V}, \mathrm{R}\in \mathbb{R}^{6}_{+}:\mathrm{S}\geq 0, \mathrm{E}\geq 0, \mathrm{I_{c}}\geq 0, \mathrm{I_{\eta}} \geq 0, \mathrm{V} \geq 0, \mathrm{R} \geq 0, \text{max}(|\mathrm{S}|, |\mathrm{E}|, |\mathrm{I_{c}}|, |\mathrm{I_{\eta}}|, |\mathrm{V}|, |\mathrm{R}|)\leq \phi\}$. Then, we will indicate that for each initial value in $\ddot{\Xi}$, the solution of system (2.1) is unique.

Theorem 3.2. The solution of system (2.1) is non-negative and ultimately uniformly bounded for all $t > 0$.

Proof. First, we will illustrate that the solution of system (2.1) is non-negative. From system (2.1), we get

Hence, we can see that $\mathrm{S}(t), \mathrm{E}(t), \mathrm{I_{c}}(t), \mathrm{I_{\eta}}(t), \mathrm{V}(t), \mathrm{R}(t)\, \geq 0$ for any $t\geq 0$. Now, we figure out the solution of system (2.1) is uniformly bounded. Let $\mathbb{N}(t) = \mathrm{S}(t)+\mathrm{E}(t)+\mathrm{I_{c}}(t)+\mathrm{I_{\eta}}(t)+\mathrm{V}(t)+\mathrm{R}(t)$. Then,

Applying the Laplace transform on Eq (3.1), we get

Then, by using the inverse Laplace transform, we find that

where $\mathbb{E}$ is the Mittag-Leffler function. Thus, according to Lemma 5 and Corollary 6 in ^[41], we have

Therefore, we find $\limsup_{t\longrightarrow\infty} \mathbb{N}(t)\, \leq \, {\Delta\over{\rho}}$, implies that $\mathrm{S}(t), \mathrm{E}(t), \mathrm{I_{c}}(t), \mathrm{I_{\eta}}(t), \mathrm{V}(t), \text{and}\, \mathrm{R}(t)$ are bounded. Thus,

is a positive invariant set with respect to system (2.1). □

Theorem 3.3. For any given initial conditions in $\ddot{\Xi}$, such that $\mathrm{S}(0), \, \mathrm{E}(0), \, \mathrm{I}_{c}(0), \, \mathrm{I}_{\eta}(0), \, \mathrm{V}(0), \, \mathrm{R}(0)\, > \, 0$, the model (2.1) has a unique solution on $[0, \infty)$ and it remains non-negative and bounded for all $t\, > \, 0$.

Proof. Let $\Psi(\chi) = (\Psi_{1}(\chi), \Psi_{2}(\chi), \Psi_{3}(\chi), \Psi_{4}(\chi), \Psi_{5}(\chi), \Psi_{6}(\chi))$ be a mapping, where

and $\chi = (\mathrm{S}, \mathrm{E}, \mathrm{I_{c}}, \mathrm{I_{\eta}}, \mathrm{V}, \mathrm{R})$. Suppose that $\varphi$ and $\bar{\varphi}$ are two arbitrary solutions of system (2.1), such that $\varphi = (\mathrm{S}, \mathrm{E}, \mathrm{I_{c}}, \mathrm{I_{\eta}}, \mathrm{V}, \mathrm{R}), \text{and}\, \bar{\varphi} = (\bar{\mathrm{S}}, \bar{\mathrm{E}}, \bar{\mathrm{I_{c}}}, \bar{\mathrm{I_{\eta}}}, \bar{\mathrm{V}}, \bar{\mathrm{R}})$. Then, we obtain

where,

As a result, $\Psi(\chi)$ satisfies Lipschitz condition, implying that the solution of system (2.1) exists and unique. □

4.   Analysis of the model

In this section, we calculate the basic reproduction number $\mathcal{R}_{0}$. Then, we explore the global stability analysis of the disease-free equilibrium point and endemic-equilibrium point by constructing the appropriate Lyaounov functions.

4.1.   The basic reproduction number

It is observed that our system (2.1) always admits a disease-free equilibrium point $U_{0} = (\mathrm{S}_{0}, 0, 0, 0, \mathrm{V}_{0}, \mathrm{R}_{0})$, and to study the existence of the endemic equilibrium point $U^{*}$ of the system (2.1), we establish the reproduction number $\mathcal{R}_{0}$ by employing the well-known method of the next generation matrix ^[42], which is a number of newly infected individuals generated from an infected individual at the beginning of the infectious process. Mathematically, it is a spectral radius of the next-generation matrix $\mathfrak{J}^{-1}\, \Pi$, where $\mathfrak{J}$ represents the positive matrix of new infection cases, which is a derivative of the non-linear terms at $U_{0}$, and $\Pi$ denotes the matrix of the transition of the infections, which is a derivative of the linear terms at $U_{0}$. For a system (2.1), we obtain

Thus, we get

Then, the basic reproduction number $\mathcal{R}_{0}$ is the spectral radius of the matrix $\mathfrak{J}^{-1}\Pi$, which is defined as

Since, $\mathcal{R}_{0, 1}$ represents the average number of secondary cases caused by contact with the infected people without symptoms while $\mathcal{R}_{0, 2}$ determines the average number of secondary cases related to infected people with symptoms.

4.2.   Steady states

Lemma 4.1. The model (2.1) has a positive basic reproduction number $\mathcal{R}_{0}$ such that:

(ⅰ) if $\mathcal{R}_{0}\, \leq\, 1$, then there exists only one fixed point.

(ⅱ) if $\mathcal{R}_{0}\, > \, 1$, then there exist two equilibrium points.

Proof. To calculate the equilibria of system (2.1), we set $(\mathrm{S}, \mathrm{E}, \mathrm{I_{c}}, \mathrm{I_{\eta}}, \mathrm{V}, \mathrm{R})$ as an equilibrium point that satisfies the following equations

Thus, from Eq (4.3) and Eq (4.4), we get

Substituting Eq (4.7) in Eq (4.1), we have

From substituting Eqs (4.7) and (4.8) in Eq (4.5), we obtain

Substituting Eqs (4.7) and (4.9) in Eq (4.6), we find

Now, from substituting Eqs (4.7)–(4.9) into Eq (4.2), we have

where,

Hence, from Eq (4.11) we see that,

(ⅰ) if $\mathrm{E}\, = \, 0$, then from Eqs (4.7)–(4.10), we get the disease-free equilibrium point $U_{0} = (S_{0}, 0, 0, 0, V_{0}, R_{0}) = \left({\Delta\over{\kappa+\rho}}, 0, 0, 0, {\Delta\, \kappa\over{(\varrho+\rho)(\kappa+\rho)}}, {\varrho\, \Delta\, \kappa\over{\rho(\varrho+\rho)(\kappa+\rho)}}\right).$

(ⅱ) if $\mathrm{E}\neq 0$, then we have $\mathcal{K}_{1}\mathrm{E^{2}}+\mathcal{K}_{2}\mathrm{E}+\mathcal{K}_{3} = 0$. Since $\mathcal{K}_{2}^{2}-4\, \mathcal{K}_{1}\, \mathcal{K}_{3}\, > 0$ and $\mathcal{K}_{3}\, < \, 0$ if and only if $\mathcal{R}_{0}\, > \, 1$, that means there exist a positive real root $\mathrm{E}^*$ when $\mathcal{R}_{0}\, > 1$. By substituting $\mathrm{E}^*$ into Eqs (4.7)–(4.10), we have

and

It is clear that the endemic equilibrium point $U^* = (\mathrm{S}^*, \mathrm{E}^*, \mathrm{I}_{c}^*, \mathrm{I}^*_{\eta}, \mathrm{V}^*, \mathrm{R}^*)$ exists if $\mathcal{R}_{0}\, > \, 1$. □

4.3.   Stability of equilibria

Clearly, from system (2.1), the equation of recovered group $\mathrm{R}$ can be ignored without any loss of generality since it does not appear in the rest of the equations. Furthermore, recovered equation $\mathrm{R}$ and the other equations have no transmission rates between them. Thus, in studying the stability of the equilibria, we will reduce the model (2.1) to the following system

Hence, we will clarify the global stability of the disease-free equilibrium point and endemic equilibrium point by establishing the Lyapunov function.

Let the function $\mathcal{T}(u(t)):\mathbb{R}^{+}\longrightarrow \mathbb{R}^{+}$ be defined as $\mathcal{T}(u(t)) = u(t)-u^{c}-u^{c}\, \ln({u(t)\over{u^{c}}})$, and note that $\mathcal{T}(u(t))$ is non-negative for any $u(t)\, > \, 0$. In addition, we set

Lemma 4.2. ^[43] Assume $\chi(t)\, \in\, \mathbb{R}^{+}$ is a continuous and derivable function. Then, for any $t\, \geq t_{0}$,

Theorem 4.3. The disease-free equilibrium point $U_{0}$ is globally asymptotically stable (GAS) if $\mathrm{R}_{0}\leq 1$.

Proof. We construct a Lyapunov function $\Theta\longrightarrow \mathbb{R}^{5}_{+}$ as follows

where, $\omega_{1} = {(\rho+d_{1})\, (\rho+d_{2})\over{(\mathrm{S}_{0}+\varepsilon\, \mathrm{V}_{0})}}, \, \, \omega_{2} = \beta_{1}\, (\rho+d_{2}), \, \text{and}\, \omega_{3} = \beta_{2}\, (\rho+d_{1}).$ It is observe that $\mathcal{Z}$ is positive definite. The time derivative of $\mathcal{Z}$ along solutions of system (4.12) is given as

By using $\Delta = (\kappa+\rho)\, \mathrm{S}_{0}$ and $\kappa\, \mathrm{S}_{0} = (\varrho+\rho)\, \mathrm{V}_{0}$, we obtain

Then,

Therefore,

Hence,

Thus, since the geometric mean is less than or equal arithmetical mean, we obtain

Moreover, when $\mathcal{R}_{0}\, \leq \, 1$, $D^{\vartheta}\mathcal{Z}\, \leq \, 0$ for all $\mathrm{S}, \mathrm{E}, \mathrm{I}_{c}, \mathrm{I}_{\eta}, \mathrm{V}\, > 0$. Note that $D^{\vartheta}\mathcal{Z}\, = \, 0$ if and only if $\mathrm{S} = \mathrm{S}_{0}, \, \mathrm{E} = 0,$ and $\mathrm{V} = \mathrm{V}_{0}$. Hence, the subset $\mathcal{H}'$ is the largest invariant set of $\mathcal{H} = \left\lbrace (\mathrm{S}, \mathrm{E}, \mathrm{I}_{c}, \mathrm{I}_{\eta}, \mathrm{V})\, \in\, \Theta:{D^{\vartheta}\mathcal{Z}} = 0\right\rbrace$. By applying La Salle's invariant principle ^[44], we observe that all solutions converge to $\mathcal{H}'$. In $\mathcal{H}'$, the elements are equal to $\mathrm{S} = \mathrm{S}_{0}, \, \mathrm{V} = \mathrm{V}_{0}$, and $\mathrm{E} = 0$. Also, from the system (4.12), we obtain $\mathrm{I}_{c} = 0$ and $\mathrm{I}_{\eta} = 0$ when $\mathrm{E} = 0$. Thus, $\mathcal{H}' = \lbrace(\mathrm{S}, \mathrm{E}, \mathrm{I}_{c}, \mathrm{I}_{\eta}, \mathrm{V})\, \in\, \mathcal{H}: \mathrm{S} = \mathrm{S}_{0}, \mathrm{V} = \mathrm{V}_{0}, \mathrm{E} = \mathrm{I}_{c} = \mathrm{I}_{\eta} = 0\rbrace\, = \, \lbrace U_{0} \rbrace$. Therefore, $U_{0}$ is GAS when $\mathcal{R}_{0}\, < 1$. □

Theorem 4.4. The endemic equilibrium point $U^*$ is GAS when $\mathcal{R}_{0}\, > \, 1$.

Proof. We define a Lyapunov function as follows:

Note that $\bar{\mathcal{Z}}$ is positive definite. By calculating the time derivative of $\bar{\mathcal{Z}}$ along the solutions of system (4.12), we obtain

From the equilibria, we have the following relationships

After substituting the above relationships and then eliminating some terms, we obtain

Hence,

Therefore,

Then,

We obtain,

Thus,

Therefore,

Finally, we have

Now, from the geometrical and arithmetical means relationship, we have

As a result, $D^\vartheta \bar{\mathcal{Z}}\, \leq 0$ if $\mathcal{R}_{0}\, > 1$ for all $\mathrm{S}, \mathrm{E}, \mathrm{I}_{c}, \mathrm{I}_{\eta}, \mathrm{V}\, > \, 0$. Moreover, $D^\vartheta \bar{\mathcal{Z}}\, = 0$ if and only if $\mathrm{S} = \mathrm{S}^*, \, \mathrm{E} = \mathrm{E}^*, \, \mathrm{I}_{c} = \mathrm{I}_{c}^*, \, \mathrm{I}_{\eta} = \mathrm{I}_{\eta}^*$ and $\mathrm{V} = \mathrm{V}^*$. Suppose that $\bar{\mathcal{H}}'$ is the largest subset of $\bar{\mathcal{H}} = \lbrace(\mathrm{S}, \mathrm{E}, \mathrm{I}_{c}, \mathrm{I}_{\eta}, \mathrm{V}):D^{\vartheta}\bar{\mathcal{Z}} = 0\rbrace$. Then, $\bar{\mathcal{H}}'\, = \, \lbrace U^*\rbrace$. Therefore, from La Salle's invariant principle ^[44], the endemic equilibrium point $U^*$ is GAS when $\mathcal{R}_{0}\, > 1$. □

5.   Numerical simulation and sensitivity analysis

5.1.   Sensitivity analysis

A valuable insight into how changes in various parameters impact the overall dynamics of the system can be provided by the sensitivity analysis of the model with respect to the basic reproduction number $\mathcal{R}_{0}$. Now, to study the sensitivity, we calculate the partial derivatives of $\mathcal{R}_{0}$ with respect to each parameter. Thus, we have

Then, we see that

and

where $\dot{\xi} = (d_{1}+\rho)\, (d_{2}+\rho)\, (\kappa+\rho)\, (\gamma+\rho)+(\varrho+\rho)\, (d_{2}+\rho)\, (\kappa+\rho)\, (\gamma+\rho)+(d_{1}+\rho)\, (\varrho+\rho)\, (\kappa+\rho)\, (\gamma+\rho)+(d_{1}+\rho)\, (d_{2}+\rho)\, (\varrho+\rho)\, (\gamma+\rho)+(d_{1}+\rho)\, (d_{2}+\rho)\, (\kappa+\rho)\, (\varrho+\rho)$. Therefore, an increase in the value of $\beta_{1}, \, \beta_{2}, \, \Delta, \, \varepsilon$, and $\gamma$ would result in a rise in $\mathcal{R}_{0}$ while no change occurs if $\rho$ varies. Conversely, the rest parameters would indicate the opposite effect.

5.2.   Numerical simulation

In this subsection, we will illustrate that our theoretical investigations are confirmed with numerical results. Therefore, we set three groups of initial conditions:

● IC1:$\mathrm{S}(0) = 25, \, \mathrm{E}(0) = 3, \, \mathrm{I}_{c}(0) = 6, \, \mathrm{I}_{\eta}(0) = 5, \, \mathrm{V}(0) = 5, \, \mathrm{R}(0) = 10.$

● IC2:$\mathrm{S}(0) = 50, \, \mathrm{E}(0) = 6, \, \mathrm{I}_{c}(0) = 12, \, \mathrm{I}_{\eta}(0) = 1, \, \mathrm{V}(0) = 5, \, \mathrm{R}(0) = 5.$

● IC3:$\mathrm{S}(0) = 60, \, \mathrm{E}(0) = 15, \, \mathrm{I}_{c}(0) = 1, \, \mathrm{I}_{\eta}(0) = 1, \, \mathrm{V}(0) = 10, \, \mathrm{R}(0) = 8.$

Now, we assume that values of the parameters are $\Delta = 10, \, \kappa = 0.4, \, \gamma = 0.8, \, \, \rho = 0.2, \, \varrho = 0.2, \, d_{1} = d_{2} = 0.2, \, \theta = 0.5, \, \varepsilon = 0.3$ while $\beta_{1}$ and $\beta_{2}, \,$ varied. Additionally, we plot the solutions of the system (2.1) at different values of $\vartheta$ ($0.95, \, 0.85, \, 0.75,$ and $1$) by using predictor-corrector PECE method of Adams-Bashforth-Moulton ^[45]. Then, we will study the following cases:

Ⅰ: The effect of $\beta_{1}$ and $\beta_{2}$ on the stability

By setting $\beta_{1} = 0.02, \, \beta_{2} = 0.02$ and using ${\bf{IC1}}$, we obtain $\mathcal{R}_{0} = 0.6933\, < \, 1$, and the result is given in Figure 2. The solution trajectories converge to $U_{0}\, = (\mathrm{S}_{0}, \mathrm{E}_{0}, \mathrm{I}_{c, 0}, \mathrm{I_{\eta, 0}}, \mathrm{V}_{0}, \mathrm{R}_{0})\, = \, (16.67, 0, 0, 0, 16.67, 16.67)$. This ensures that $U_{0}$ is GAS based on the result of Theorem 4.3. This means the spread of the infectious disease is decreasing.

By setting $\beta_{1} = 0.2, \, \beta_{2} = 0.02$ and using ${\bf{IC1}}$, we obtain $\mathcal{R}_{0} = 3.8133\, > \, 1$, and the result is given in Figure 3. The solution trajectories converge to $U^*\, = (\mathrm{S}^*, \mathrm{E}^*, \mathrm{I}_{c}^*, \mathrm{I}_{\eta}^*, \mathrm{V}^*, \mathrm{R}^*) = (4.96, 8.04, 6.43, 6.43, 2.40, 21.71)$. Also, when $\beta_{1} = 0.1, \, \beta_{2} = 0.1$, we get $\mathcal{R}_{0} = 3.4667 > 1$, as shown in Figure 4, and it is clear that the solutions tend to $U^*\, = (\mathrm{S}^*, \mathrm{E}^*, \mathrm{I}_{c}^*, \mathrm{I}_{\eta}^*, \mathrm{V}^*, \mathrm{R}^*) = (5.41, 7.79, 6.23, 6.23, 2.79, 21.51)$. This prove that, from Figures 3 and 4, the equilibrium point $U^*$ is GAS based on the result of Theorem 4.4. Biologically, this indicates there is an outbreak of infectious disease.

In Figure 5, we plot the solutions at different initial conditions and $\vartheta = 0.95$, and we also assume that $\beta_{1} = 0.02, \, \beta_{2} = 0.02$. Then we have $\mathcal{R}_{0} = 0.6933$ and $\beta_{1} = 0.2, \, \beta_{2} = 0.02$, which implies $\mathcal{R}_{0} = 3.8133$. Thus, we observe that the results support the theoretical results of Theorems 4.3 and 4.4. Finally, we notice that the proposed fractional system converts to the following ODE system at $\vartheta = 1$

Then, we solve the ODE system (5.12) by using the Runge-Kutta method (rk4). Figures 2–4 show that the solution of system (2.1) is consistent with the solution of system (5.12) when $\vartheta$ is equal to one. In addition, from Figures 2–4, we conclude that the increasing of fractional order reduces the outbreak of the infectious disease.

Ⅱ: The effect of $\kappa$ on equilibrium point

In this part, we select different values of $\kappa$ to examine the effect of the vaccination rate $\kappa$ on the transmission of the infectious disease. We fix the values $\Delta = 10, \, \varepsilon = 0.3, \, \gamma = 0.8, \, \, \rho = 0.2\, \varrho = 0.2, \, d_{1} = d_{2} = 0.2, \, \theta = 0.5,$ and $\beta_{1} = \, \beta_{2} = 0.1$ with IC2. Table 2 illustrates that increasing $\kappa$ leads to decreasing $\mathcal{R}_{0}$ and the number of infected individuals. This result is shown in Figure 6.

6.   Conclusions

In this work, we presented a $SEI_{c}I_{\eta}VR$ epidemic model with two types of infected individuals and a vaccination strategy. The model has six compartments: susceptible $\mathrm{S}(t)$, exposed $\mathrm{E}(t)$, asymptomatic infected $\mathrm{I_{c}}(t)$, symptomatic infected $\mathrm{I_{\eta}}(t)$, vaccinated $\mathrm{V}(t)$, and recovered $\mathrm{R}(t)$. We proved the existence, positivity, and boundedness of all solutions in this model. By applying the next-generation method, we calculated the basic reproduction number $\mathcal{R}_{0}$ that controls the existence and stability of the equilibria. By applying the Lyapunov method together with LaSalle's invariance principle, we derived that if the basic reproduction number $\mathcal{R}_{0}$ is less than one, then the free-disease equilibrium point is GAS and if $\mathcal{R}_{0}$ is greater than one, then the endemic equilibrium point is GAS. In addition, the effect of the vaccination rate $\kappa$ is illustrated numerically and we conclude that, as $\kappa$ increases, $\mathcal{R}_{0}$ is decreased. This means the vaccine can be useful in reducing the spread of infectious diseases. Finally, some of the numerical simulations were introduced to support our theoretical results by using MATLAB.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The author declare no conflict of interest.