Fractional Order Mathematical Modelling of HFMD Transmission via Caputo Derivative

1. Introduction

Mathematical modeling involves using equations or algorithms to represent real-world systems across various disciplines such as physics, engineering, biology, economics, and social sciences. These models serve to comprehend, evaluate, and predict system behavior, playing a vital role in scientific research, engineering design, and decision-making. Mathematical modeling allows the study and control of systems in ways often impractical or costly in real-world experiments. The process typically includes problem formulation, model construction, parameter estimation, model validation, and model analysis. Proficiency in both system understanding and mathematics is crucial for creating and solving these equations. Validation and verification of models are essential before using them for predictions or decision-making. This interdisciplinary field combines mathematical concepts with domain-specific knowledge, providing insightful analyses of complex phenomena. Mathematical modeling is indispensable for advancing scientific understanding, addressing real-world issues, and making well-informed decisions. Researchers, engineers, and decision-makers leverage mathematical modeling to gain knowledge, predict outcomes, and optimize systems, contributing significantly to various fields [1,2,3,4].

The viral illness known as hand, foot, and mouth disease (HFMD) predominantly affects children under the age of ten and occasionally even adults. There are several enterovirus types that cause it, but Coxsackievirus A16 and Enterovirus 71 are the most frequent ones. The manifestation of diminutive and distressing blisters or ulcers on the hands, feet, and oral cavity, coupled with other indicators like fever and rash, typically points to the likelihood of Hand, Foot, and Mouth Disease (HFMD). It is highly contagious and spreads through close contact with infected individuals, contaminated surfaces, or respiratory droplets from coughs and sneezes. It is most commonly found in crowded places such as day care centres, schools, and playgrounds, and outbreaks often occur in communities, especially during the warmer months. The symptoms of HFMD usually start with fever, sore throat, and a general feeling of malaise. Within a day or two, small blisters will develop in the hands, feet and mouth. These blisters may be painful and can sometimes become ulcers. Other symptoms may include loss of appetite, headache, and in some cases, a skin rash. Although HFMD is often a temporary disease that cure by itself within seven to ten days, it can cause discomfort and sometimes complications, especially in severe cases. Complications can include viral meningitis, encephalitis (inflammation of the brain), and in rare cases, more serious neurological complications. Treatment for HFMD is usually supportive, focusing on relieving symptoms such as fever and discomfort. To manage discomfort and lower fever, over-the-counter pain medicines such as acetaminophen or ibuprofen may be used, and maintaining excellent hygiene habits, such as frequent handwashing and avoiding close contact with sick persons, is essential for preventing the virus’s spread. The important prevention is to maintain good oral hygiene and minimising contact with ill peoples. Vaccines for HFMD are not widely available, but some countries have developed vaccines against Enterovirus 71, which is one of the common causes of severe HFMD cases. Some recent research are presented in [5,6,7].

Fractional differential equations (FDEs) are differential equations that incorporate non-integer order derivatives. This permit the use of fractional orders, which can be any real or complex number, in contrast to ordinary differential equations (ODEs), where the order of the derivatives is always a positive integer [8]. When a system’s behaviour depends on the fractional derivatives of its variables, FDEs are an effective mathematical tool for describing a wide range of physical, biological, and engineering phenomena. In 18th and 19th centuries, mathematicians especially Leibniz and Riemann established the concept of fractional calculus, that is associated with non-integer order derivatives and integrals. However, because of its applications in numerous disciplines, including physics, engineering, finance, and biology, it has recently attracted new interest. Among other complex phenomena, FDEs have been used to model and examine viscoelastic materials, diffusion in fractal media, electrochemical processes, and biological systems. Most of the models stated above are integer order models, however they are unable to depict real life difficulties because there is no non-locality impact in local differentiation. Mathematicians thus developed the idea of differentiation with non-local operators to address this issue. The memory effect of fractional calculus has been demonstrated to produce more accurate results when predicting physical systems, including mathematical models [9]. In order to solve systems of differential equations of both integer and non-integer orders resulting from real life issues, such as applications related to integrodifferential equations, [10,11], mathematical epidemiology [12,13], economic and financial [14,15], new discoveries have resulted from advances in constructing new fractional order operators, especially Caputo, Atangana-Baleanu and Caputo-Fabrizio. Most researchers are optimistic on the fractional-order differential equation, especially in biological modelling [16]. A unique interpretation of the fractional operator without a single kernel was proposed by Caputo and Fabrizio [17]. In this paper, we choose Hand foot mouth disease as our social problem and formulated a five compartmental nonlinear mathematical model. Further, we extend the nonlinear ODE model into fractional order model by using Caputo derivative.

An overview of our article is provided below. In Section 2, the fractional order model was formulated, followed by some essential preliminaries in Section 3. In Section 4, we addressed the system’s positivity and boundedness. Section 5 discusses the suggested model’s existence and uniqueness, followed by the qualitative analysis in Section 6. The 7th section explains how to use the modified Euler’s approach to get the general solution to the suggested system. In Section 8, a brief summary was presented.

2. Model Formulation

We consider a $SVEIR$, five compartment HFMD model. The basic structure of the model comprises distinct population groups, including those that are susceptible, vaccinated, exposed, infected and recovered. In this case, we presume that not all vaccinated persons are immune to this virus, and that some may be exposed to it. We also believe that persons who have recovered from HFMD do not have a lifetime immunity to the condition. As a result, we believe that recovered people are susceptible at a rate of ${\alpha }_0$, and that after recovery, some people will take vaccination and therefore enter the vaccinated compartment at a rate of ${\alpha }_1$. As a result of the aforesaid assumptions, a basic nonlinear mathematical model was developed in [21]. That is,

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =A-{\delta }_{1}SI-{\delta }_{2}SE-\mu S-zS+{\alpha }_{0}R,\hfill \\ \hfill \frac{dV}{dt}& =zS-{\delta }_{3}VI-{\delta }_{4}VE-\mu V+{\alpha }_{1}R,\hfill \\ \hfill \frac{dE}{dt}& ={\delta }_{1}SI+{\delta }_{2}SE+{\delta }_{3}VI+{\delta }_{4}VE-\eta E-\mu E,\hfill \\ \hfill \frac{dI}{dt}& =\eta E-\gamma I-\mu I-{\mu }_{1}I,\hfill \\ \hfill \frac{dR}{dt}& =\gamma I-\mu R-{\alpha }_{0}R-{\alpha }_{1}R.\hfill \end{array}$$

(1)

It is vital to investigate viral disease mathematical models in order to have a deeper comprehension of their assessment, presence, stability, and control. Because typical mathematical models do not provide a great level of precision in describing these diseases, FDE which have many uses in practical domains, were designed to manage such situations. So, here ODE model for HFMD in [21] was extended to the fractional order model by using the Caputo derivative.

$$\begin{array}{cc}\hfill { }^C{\mathbb{D}}_t^{\alpha }\left[S\left(t\right)\right]& =A-{\delta }_{1}SI-{\delta }_{2}SE-\mu S-zS+{\alpha }_{0}R,\hfill \\ \hfill { }^C{\mathbb{D}}_t^{\alpha }\left[V\left(t\right)\right]& =zS-{\delta }_{3}VI-{\delta }_{4}VE-\mu V+{\alpha }_{1}R,\hfill \\ \hfill { }^C{\mathbb{D}}_t^{\alpha }\left[E\left(t\right)\right]& ={\delta }_{1}SI+{\delta }_{2}SE+{\delta }_{3}VI+{\delta }_{4}VE-\eta E-\mu E,\hfill \\ \hfill { }^C{\mathbb{D}}_t^{\alpha }\left[I\left(t\right)\right]& =\eta E-\gamma I-\mu I-{\mu }_{1}I,\hfill \\ \hfill { }^C{\mathbb{D}}_t^{\alpha }\left[R\left(t\right)\right]& =\gamma I-\mu R-{\alpha }_{0}R-{\alpha }_{1}R.\hfill \end{array}$$

(2)

Here, 0 <$\alpha$≤ 1 and $S\left(0\right)$ = $S_0$, $V\left(0\right)$ = $V_0$, $E\left(0\right)$ = $E_0$, $I\left(0\right)$ = $I_0$ and $R\left(0\right)$ = $R_0$ are the initial conditions. Table 1 contains a description of the model (2) parameters.

3. Preliminaries

In this, we present the essential definitions of fractional calculus within the context of Caputo significance.

Definition 3.1. 

Consider a function y ∈${\mathbb{C}}^m$, then ${ }^C{\mathbb{D}}_t^{\alpha }$ in (m-1,m) where m ∈$\mathbb{N}$ is provided by,

$${ }^C{\mathbb{D}}_t^{\alpha }\left(g\left(t\right)\right)=\frac{1}{\Gamma \left(m-\alpha \right)}{\int }_0^t{(t-s)}^{(m-\alpha -1)}{g}^{\left(m\right)}\left(s\right)ds.$$

Here, ${ }^C{\mathbb{D}}_t^{\alpha }\left(g\left(t\right)\right)$→$g^{\prime }\left(t\right)$ as α tends to 1.

Definition 3.2. 

For α> 0, the corresponding integral of $g:{\mathbb{R}}^{+}\to \mathbb{R}$ is expressed as follows,

$$I_t^{\alpha }\left(g\left(t\right)\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-s)}^{(\alpha -1)}g\left(s\right)ds,$$

where, 0 <α< 1, t> 0.

4. Positivity and Boundedness

Here we demonstrate the system’s boundedness and non-negative solutions of (2).

Theorem 4.1. 

Every solution of (2) are non negative and bounded.

Proof. 

By using the same approach from [18], We define $\frac{dN}{dt}$ = $\frac{dS}{dt}+\frac{dV}{dt}+\frac{dE}{dt}+\frac{dI}{dt}+\frac{dR}{dt}$. Also let consider $\mu >0.$

$$\begin{array}{cc}\hfill { }^C{\mathbb{D}}_t^{\alpha }\left[N\left(t\right)\right]+\mu N\left(t\right)& =A-{\delta }_{1}SI-{\delta }_{2}SE-\mu S-zS+{\alpha }_{0}R+zS-{\delta }_{3}VI-{\delta }_{4}VE-\mu V+{\alpha }_{1}R\hfill \\ & +{\delta }_{1}SI+{\delta }_{2}SE+{\delta }_{3}VI+{\delta }_{4}VE-\eta E-\mu E+\eta E-\gamma I-\mu I-{\mu }_{1}I\hfill \\ & +\gamma I-\mu R-{\alpha }_{0}R-{\alpha }_{1}R+\mu S+\mu V+\mu E+\mu I+\mu R\hfill \\ \hfill & =A-{\mu }_{1}I\hfill \\ \hfill & \le A\hfill \end{array}$$

By using standard comparison theorem [19] for (2), we get,

$$\begin{array}{cc}\hfill N\left(t\right)& \le N\left(0\right)E_{\alpha }\left(-\mu {\left(t\right)}^{\alpha }\right)+A{\left(t\right)}^{\alpha }{E}_{\alpha ,\alpha +1}\left(-\mu {\left(t\right)}^{\alpha }\right)\hfill \end{array}$$

where the Mittag Leffler function is represented by $E_{\alpha }$. As stated by [19], we have

$$\begin{array}{cc}\hfill N\left(t\right)& \le \frac{A}{\mu },t\to \alpha .\hfill \end{array}$$

Hence, all solutions of equation (2) are restricted to the domain $\Omega$,

$$\Omega =\{(S,V,E,I,R)\in {\mathbb{R}}_{+}^5|N\left(t\right)\le \frac{A}{\mu }+{ϵ}_0,forany{ϵ}_0>0\}.$$

(3)

Now, we show the solutions of (2) are positive. From first equation of (2), we have,

$$\begin{array}{cc}\hfill { }^C{\mathbb{D}}_t^{\alpha }\left[S\left(t\right)\right]& =A-{\delta }_{1}SI-{\delta }_{2}SE-\mu S-zS+{\alpha }_{0}R\hfill \\ \hfill & \ge -[{\delta }_{1}I+{\delta }_{2}E+\mu +z]S\hfill \\ \hfill & \ge -[({\delta }_1+{\delta }_2)\frac{A}{\mu }+\mu +z]S\hfill \\ \hfill & \ge -d_{1}S.\hfill \end{array}$$

where $d_1=({\delta }_1+{\delta }_2)\frac{A}{\mu }+\mu +z$.

As per the traditional comparison theorem [19] and the non negativity of Mitag Leffler function, we have, $E_{\alpha ,1}\left(t\right)>0$ for any $\alpha \in (0,1)$,

$$\begin{array}{c}\hfill S\ge S\left(0\right)E_{\alpha ,1}(-d_1{t}^{\alpha })\Rightarrow S\ge 0.\end{array}$$

From (2), we get

$$\begin{array}{cc}\hfill { }^C{\mathbb{D}}_t^{\alpha }\left[V\left(t\right)\right]& =zS-{\delta }_{3}VI-{\delta }_{4}VE-\mu V+{\alpha }_{1}R\hfill \\ \hfill & \ge -[{\delta }_{3}I+{\delta }_{4}E+\mu ]V\hfill \\ \hfill & \ge -[({\delta }_3+{\delta }_4)\frac{A}{\mu }+\mu ]V\hfill \\ \hfill & \ge -d_{2}V.\hfill \end{array}$$

where $d_2=({\delta }_3+{\delta }_4)\frac{A}{\mu }+\mu$.

Therefore,

$$\begin{array}{c}\hfill V\ge V\left(0\right)E_{\alpha ,1}(-d_2{t}^{\alpha })\Rightarrow V\ge 0.\end{array}$$

For third equation of (2), we get

$$\begin{array}{cc}\hfill { }^C{\mathbb{D}}_t^{\alpha }\left[E\left(t\right)\right]& ={\delta }_{1}SI+{\delta }_{2}SE+{\delta }_{3}VI+{\delta }_{4}VE-\eta E-\mu E\hfill \\ \hfill & \ge -[-{\delta }_{2}S-{\delta }_{4}V+\eta +\mu ]E\hfill \\ \hfill & \ge -[-\frac{{\delta }_{2}A}{\mu }-\frac{{\delta }_{4}A}{\mu }+\eta +\mu ]E\hfill \\ \hfill & \ge -d_{3}E.\hfill \end{array}$$

Here $d_3=-\frac{{\delta }_{2}A}{\mu }-\frac{{\delta }_{4}A}{\mu }+\eta +\mu$.

$$\begin{array}{c}\hfill \therefore E\ge E\left(0\right)E_{\alpha ,1}(-d_3{t}^{\alpha })\Rightarrow E\ge 0.\end{array}$$

For fourth equation of (2), we get

$$\begin{array}{cc}\hfill { }^C{\mathbb{D}}_t^{\alpha }\left[I\left(t\right)\right]& =\eta E-\gamma I-\mu I-{\mu }_{1}I\hfill \\ \hfill & \ge -[\gamma +\mu +{\mu }_1]I\hfill \\ \hfill & \ge -d_{4}I.\hfill \end{array}$$

where $d_4=\gamma +\mu +{\mu }_1$.

$$\begin{array}{c}\hfill \therefore I\ge I\left(0\right)E_{\alpha ,1}(-d_4{t}^{\alpha })\Rightarrow I\ge 0.\end{array}$$

For fifth equation of (2), we get

$$\begin{array}{cc}\hfill { }^C{\mathbb{D}}_t^{\alpha }\left[R\left(t\right)\right]& =\gamma I-\mu R-{\alpha }_{0}R-{\alpha }_{1}R\hfill \\ \hfill & \ge -[\mu +{\alpha }_0+{\alpha }_1]I\hfill \\ \hfill & \ge -d_{5}R.\hfill \end{array}$$

where $d_5=\mu +{\alpha }_0+{\alpha }_1$.

Therefore,

$$\begin{array}{c}\hfill R\ge R\left(0\right)E_{\alpha ,1}(-d_5{t}^{\alpha })\Rightarrow R\ge 0.\end{array}$$

Thus, the solutions of (2) are positive. □

5. Existence and Uniqueness

In this, we established the formulated fractional order system’s existence and uniqueness. The following equations were found after utilising the equation (2).

$$\begin{array}{cc}\hfill S\left(t\right)& =S\left(0\right)+{ }^C{I}_0^{\alpha }S\left(t\right)\{A-{\delta }_{1}SI-{\delta }_{2}SE-\mu S-zS+{\alpha }_{0}R\},\hfill \\ \hfill V\left(t\right)& =V\left(0\right)+{ }^C{I}_0^{\alpha }V\left(t\right)\{zS-{\delta }_{3}VI-{\delta }_{4}VE-\mu V+{\alpha }_{1}R\},\hfill \\ \hfill E\left(t\right)& =E\left(0\right)+{ }^C{I}_0^{\alpha }E\left(t\right)\{{\delta }_{1}SI+{\delta }_{2}SE+{\delta }_{3}VI+{\delta }_{4}VE-\eta E-\mu E\},\hfill \\ \hfill I\left(t\right)& =I\left(0\right)+{ }^C{I}_0^{\alpha }I\left(t\right)\{\eta E-\gamma I-\mu I-{\mu }_{1}I\},\hfill \\ \hfill R\left(t\right)& =R\left(0\right)+{ }^C{I}_0^{\alpha }R\left(t\right)\{\gamma I-\mu R-{\alpha }_{0}R-{\alpha }_{1}R\}.\hfill \end{array}$$

By using Definition 3.2,

$$\begin{array}{cc}\hfill S\left(t\right)& =S\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}{\mathbb{K}}_1(\sigma ,S)d\sigma ,\hfill \\ \hfill V\left(t\right)& =V\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}{\mathbb{K}}_2(\sigma ,V)d\sigma ,\hfill \\ \hfill E\left(t\right)& =E\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}{\mathbb{K}}_3(\sigma ,E)d\sigma ,\hfill \\ \hfill I\left(t\right)& =I\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}{\mathbb{K}}_4(\sigma ,I)d\sigma ,\hfill \\ \hfill R\left(t\right)& =R\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}{\mathbb{K}}_5(\sigma ,R)d\sigma .\hfill \end{array}$$

(4)

Here,

$$\begin{array}{cc}\hfill {\mathbb{K}}_1(\sigma ,S)& =A-{\delta }_{1}SI-{\delta }_{2}SE-\mu S-zS+{\alpha }_{0}R,\hfill \\ \hfill {\mathbb{K}}_2(\sigma ,V)& =zS-{\delta }_{3}VI-{\delta }_{4}VE-\mu V+{\alpha }_{1}R,\hfill \\ \hfill {\mathbb{K}}_3(\sigma ,E)& ={\delta }_{1}SI+{\delta }_{2}SE+{\delta }_{3}VI+{\delta }_{4}VE-\eta E-\mu E,\hfill \\ \hfill {\mathbb{K}}_4(\sigma ,I)& =\eta E-\gamma I-\mu I-{\mu }_{1}I,\hfill \\ \hfill {\mathbb{K}}_5(\sigma ,R)& =\gamma I-\mu R-{\alpha }_{0}R-{\alpha }_{1}R.\hfill \end{array}$$

(5)

Assume all population groups are positive bounded function that is ∃ some positive constants ${\Pi }_1,{\Pi }_2,{\Pi }_3,{\Pi }_4,{\Pi }_5$, such that

$$\begin{gathered} \parallel S\left(t\right)\parallel \le {\Pi }_1, \\ \parallel V\left(t\right)\parallel \le {\Pi }_2, \\ \parallel E\left(t\right)\parallel \le {\Pi }_3, \\ \parallel I\left(t\right)\parallel \le {\Pi }_4, \\ \parallel R\left(t\right)\parallel \le {\Pi }_5. \end{gathered}$$

Theorem 5.1. 

If $0\le \mathbb{M}=max\{{\Delta }_1,{\Delta }_2,{\Delta }_3,{\Delta }_4,{\Delta }_5\}<1$, then ${\mathbb{K}}_i$ for $i=1,2,3,4,5$ satisfies Lipchitz conditions.

Proof. 

Consider ${\mathbb{K}}_1$, for any S and $S_1$,

$$\begin{array}{cc}\hfill \parallel {\mathbb{K}}_1(t,S)-{\mathbb{K}}_1(t,S_1)\parallel & =\parallel -{\delta }_{1}IS-{\delta }_{2}ES-\mu S-zS+{\delta }_{1}I{S}_1+{\delta }_{2}E{S}_1+\mu S_1+z{S}_1\parallel \hfill \\ \hfill & =\parallel {\delta }_{1}I(S_1-S)+{\delta }_{2}E(S_1-S)+\mu (S_1-S)+z(S_1-S)\parallel \hfill \\ \hfill & \le \left[{\delta }_1\parallel I\left(t\right)\parallel +{\delta }_2\parallel E\left(t\right)\parallel +\mu +z\right]\parallel S_1-S\parallel \hfill \\ \hfill & \le {\Delta }_1\parallel S-S_1\parallel \hfill \end{array}$$

where ${\Delta }_1={\delta }_1{\Pi }_4+{\delta }_2{\Pi }_3+\mu +z.$ Thus ${\mathbb{K}}_1$ satisfy Lipchitz condition. Similarly we can find ${\Delta }_j$ for $j=2,3,4,5.$ So that ${\mathbb{K}}_i$ for $i=2,3,4,5$ the lipchitz’s requirements are met under the condition, $0\le \mathbb{M}=max\{{\Delta }_1,{\Delta }_2,{\Delta }_3,{\Delta }_4,{\Delta }_5\}<1$, the functions are contractions. Hence the proof.

□

Now, we rewrite (4) recursively as,

$$\begin{array}{cc}\hfill S_n\left(t\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}{\mathbb{K}}_1(\sigma ,S_{n-1})d\sigma ,\hfill \\ \hfill V_n\left(t\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}{\mathbb{K}}_2(\sigma ,V_{n-1})d\sigma ,\hfill \\ \hfill E_n\left(t\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}{\mathbb{K}}_3(\sigma ,E_{n-1})d\sigma ,\hfill \\ \hfill I_n\left(t\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}{\mathbb{K}}_4(\sigma ,I_{n-1})d\sigma ,\hfill \\ \hfill R_n\left(t\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}{\mathbb{K}}_5(\sigma ,R_{n-1})d\sigma .\hfill \end{array}$$

(6)

The difference between two terms can be represented as

$$\begin{array}{cc}\hfill {\Lambda }_n\left(t\right)& =S_n\left(t\right)-S_{n-1}\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}\left[{\mathbb{K}}_1(\sigma ,S_{n-1})-{\mathbb{K}}_1(\sigma ,S_{n-2})\right]d\sigma ,\hfill \\ \hfill {\Omega }_n\left(t\right)& =V_n\left(t\right)-V_{n-1}\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}\left[{\mathbb{K}}_2(\sigma ,V_{n-1})-{\mathbb{K}}_2(\sigma ,V_{n-2})\right]d\sigma ,\hfill \\ \hfill {\chi }_n\left(t\right)& =E_n\left(t\right)-E_{n-1}\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}\left[{\mathbb{K}}_3(\sigma ,E_{n-1})-{\mathbb{K}}_3(\sigma ,E_{n-2})\right]d\sigma ,\hfill \\ \hfill {\psi }_n\left(t\right)& =I_n\left(t\right)-I_{n-1}\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}\left[{\mathbb{K}}_4(\sigma ,I_{n-1})-{\mathbb{K}}_4(\sigma ,I_{n-2})\right]d\sigma ,\hfill \\ \hfill {\vartheta }_n\left(t\right)& =R_n\left(t\right)-R_{n-1}\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}\left[{\mathbb{K}}_5(\sigma ,R_{n-1})-{\mathbb{K}}_5(\sigma ,R_{n-2})\right]d\sigma .\hfill \end{array}$$

(7)

where,

$$\begin{array}{cc}\hfill S_n\left(t\right)& =\sum _{i=0}^n{\Lambda }_i\left(t\right),\hfill \\ \hfill V_n\left(t\right)& =\sum _{i=0}^n{\Omega }_i\left(t\right),\hfill \\ \hfill E_n\left(t\right)& =\sum _{i=0}^n{\chi }_i\left(t\right),\hfill \\ \hfill I_n\left(t\right)& =\sum _{i=0}^n{\psi }_i\left(t\right),\hfill \\ \hfill R_n\left(t\right)& =\sum _{i=0}^n{\vartheta }_i\left(t\right).\hfill \end{array}$$

(8)

Now consider,

$$\begin{array}{cc}\hfill \parallel {\Lambda }_n\left(t\right)\parallel & =\parallel S_n\left(t\right)-S_{n-1}\left(t\right)\parallel \hfill \\ \hfill & =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}\left[{\mathbb{K}}_1(\sigma ,S_{n-1})-{\mathbb{K}}_1(\sigma ,S_{n-2})\right]d\sigma \hfill \\ \hfill & =\frac{{\Delta }_1}{\Gamma \left(\alpha \right)}{\int }_0^t\parallel S_{n-1}-S_{n-2}\parallel d\sigma \hfill \\ \hfill & =\frac{{\Delta }_1}{\Gamma \left(\alpha \right)}{\int }_0^t\parallel {\Lambda }_{n-1}\left(t\right)\parallel d\sigma .\hfill \end{array}$$

In the same fashion, we obtain

$$\begin{array}{cc}\hfill \parallel {\Omega }_n\left(t\right)\parallel & =\frac{{\Delta }_2}{\Gamma \left(\alpha \right)}{\int }_0^t\parallel {\Omega }_{n-1}\left(t\right)\parallel d\sigma ,\hfill \\ \hfill \parallel {\chi }_n\left(t\right)\parallel & =\frac{{\Delta }_3}{\Gamma \left(\alpha \right)}{\int }_0^t\parallel {\chi }_{n-1}\left(t\right)\parallel d\sigma ,\hfill \\ \hfill \parallel {\psi }_n\left(t\right)\parallel & =\frac{{\Delta }_4}{\Gamma \left(\alpha \right)}{\int }_0^t\parallel {\psi }_{n-1}\left(t\right)\parallel d\sigma ,\hfill \\ \hfill \parallel {\vartheta }_n\left(t\right)\parallel & =\frac{{\Delta }_5}{\Gamma \left(\alpha \right)}{\int }_0^t\parallel {\vartheta }_{n-1}\left(t\right)\parallel d\sigma .\hfill \end{array}$$

Theorem 5.2.(i) The functions stated in (8) exists and are smooth, (ii) If ∃$t_0>1$∋$\frac{{\Delta }_i}{\Gamma \left(\alpha \right)}{t}_0$ strictly less than one, $i=1,2,3,4,5$, then atleast one solution of the system exist.

Proof. (i) Since all populations are bounded and for each ${\mathbb{K}}_i$, where $i=1,2,3,4,5$ satisfy Lipchitz’s conditions, then we will get the below relation.

$$\begin{array}{cc}\hfill \parallel {\Lambda }_n\left(t\right)\parallel & \le \parallel S\left(0\right)\parallel {∥\frac{{\Delta }_1}{\Gamma \left(\alpha \right)}t∥}^n,\hfill \\ \hfill \parallel {\Omega }_n\left(t\right)\parallel & \le \parallel V\left(0\right)\parallel {∥\frac{{\Delta }_2}{\Gamma \left(\alpha \right)}t∥}^n,\hfill \\ \hfill \parallel {\chi }_n\left(t\right)\parallel & \le \parallel E\left(0\right)\parallel {∥\frac{{\Delta }_3}{\Gamma \left(\alpha \right)}t∥}^n,\hfill \\ \hfill \parallel {\psi }_n\left(t\right)\parallel & \le \parallel I\left(0\right)\parallel {∥\frac{{\Delta }_4}{\Gamma \left(\alpha \right)}t∥}^n,\hfill \\ \hfill \parallel {\vartheta }_n\left(t\right)\parallel & \le \parallel R\left(0\right)\parallel {∥\frac{{\Delta }_5}{\Gamma \left(\alpha \right)}t∥}^n.\hfill \end{array}$$

(9)

Therefore equation (9) shows the existence and smoothness of all populations which was defined in (8).

(ii) Here we show, $S_n^{\prime },V_n^{\prime },E_n^{\prime },I_n^{\prime }$ and $R_n^{\prime }$ converge to solutions of (2). Now let, ${u_n^{*}}^{\prime },$${w_n^{*}}^{\prime },$${x_n^{*}}^{\prime },$${y_n^{*}}^{\prime },$${z_n^{*}}^{\prime },$ are remainder terms after n iterations, ∋

$$\begin{array}{cc}\hfill S^{\prime }-S\left(0\right)& =S_n^{\prime }-{u_n^{*}}^{\prime },\hfill \\ \hfill V^{\prime }-V\left(0\right)& =V_n^{\prime }-{w_n^{*}}^{\prime },\hfill \\ \hfill E^{\prime }-E\left(0\right)& =E_n^{\prime }-{x_n^{*}}^{\prime },\hfill \\ \hfill I^{\prime }-I\left(0\right)& =I_n^{\prime }-{y_n^{*}}^{\prime },\hfill \\ \hfill R^{\prime }-R\left(0\right)& =R_n^{\prime }-{z_n^{*}}^{\prime }.\hfill \end{array}$$

(10)

By using triangle inequality along with ${\mathbb{K}}_1$, we get

$$\begin{array}{cc}\hfill \parallel u_n^{*}\left(t\right)\parallel & =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}\left[{\mathbb{K}}_1(\sigma ,S)-{\mathbb{K}}_1(\sigma ,S_{n-1})\right]d\sigma \hfill \\ \hfill & \le \frac{{\Delta }_1}{\Gamma \left(\alpha \right)}\parallel S-S_{n-1}\parallel t\hfill \end{array}$$

We obtain the below equation by iteratively using the previous process,

$$\begin{array}{cc}\hfill \parallel u_n^{*}\left(t\right)\parallel & \le {∥\frac{{\Delta }_1}{\Gamma \left(\alpha \right)}t∥}^{n+1}{\Pi }_1\hfill \end{array}$$

Then at $t_0$, we have

$$\begin{array}{cc}\hfill \parallel u_n^{*}\left(t\right)\parallel & \le {∥\frac{{\Delta }_1}{\Gamma \left(\alpha \right)}{t}_0∥}^{n+1}{\Pi }_1.\hfill \end{array}$$

Taking limit as $n\to \infty$,

$$\begin{array}{cc}\hfill \underset{n\to \infty }{lim}\parallel u_n^{*}\left(t\right)\parallel & \le \underset{n\to \infty }{lim}{∥\frac{{\Delta }_1}{\Gamma \left(\alpha \right)}{t}_0∥}^{n+1}{\Pi }_1.\hfill \end{array}$$

(11)

Using hypothesis, $\frac{{\Delta }_1}{\Gamma \left(\alpha \right)}{t}_0<1$. Equation (11) becomes,

$$\begin{array}{c}\hfill \underset{n\to \infty }{lim}\parallel u_n^{*}\left(t\right)\parallel =0\end{array}$$

As $n\to \infty$, we get

$$\begin{array}{c}\hfill \parallel w_n^{*}\left(t\right)\parallel \to 0,\\ \hfill \parallel x_n^{*}\left(t\right)\parallel \to 0,\\ \hfill \parallel y_n^{*}\left(t\right)\parallel \to 0,\\ \hfill \parallel z_n^{*}\left(t\right)\parallel \to 0.\end{array}$$

Thus, the system has at least one solution.

□

Theorem 5.3. 

A unique solution exists for (2) if $\left(1-\frac{{\Delta }_i}{\Gamma \left(\alpha \right)}t\right)>0$ for $i=1\mathrm{to}5$.

Proof. 

Let $S_1\left(t\right),V_1\left(t\right),E_1\left(t\right),I_1\left(t\right),R_1\left(t\right)$ is another set of solution of (2) thus,

$$\begin{array}{cc}\hfill \parallel S\left(t\right)-S_1\left(t\right)\parallel & =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\sigma )}^{\alpha -1}\left[{\mathbb{K}}_1(s,S)-{\mathbb{K}}_1(s,S_1)\right]ds\hfill \\ \hfill & \le \frac{{\Delta }_1}{\Gamma \left(\alpha \right)}t\parallel S\left(t\right)-S_1\left(t\right)\parallel \hfill \end{array}$$

By rearranging,

$$\begin{array}{cc}\hfill \parallel S\left(t\right)-S_1\left(t\right)\parallel -\frac{{\Delta }_1}{\Gamma \left(\alpha \right)}t\parallel S\left(t\right)-S_1\left(t\right)\parallel & \le 0\hfill \end{array}$$

$$\begin{array}{cc}\hfill \parallel S\left(t\right)-S_1\left(t\right)\parallel \left[1-\frac{{\Delta }_1}{\Gamma \left(\alpha \right)}t\right]& \le 0\hfill \end{array}$$

(12)

using hypothesis $\left(1-\frac{{\Delta }_1}{\Gamma \left(\alpha \right)}t\right)$ is strictly greater than zero, then (12) takes the structure $\parallel S\left(t\right)-S_1\left(t\right)\parallel =0$.

It means that $S\left(t\right)=S_1\left(t\right).$

By repeating the similar procedure to each solution for $i=2,3,4,5$, we get,

$$\begin{array}{cc}\hfill V\left(t\right)& =V_1\left(t\right),\hfill \\ \hfill E\left(t\right)& =E_1\left(t\right),\hfill \\ \hfill I\left(t\right)& =I_1\left(t\right),\hfill \\ \hfill R\left(t\right)& =R_1\left(t\right).\hfill \end{array}$$

Hence the proof. □

6. Model Analysis

This section examines the existence of equilibrium points and the local stability of the system.

6.1. Equilibrium Points

In our case, we found two equilibrium points: $E_0$= $(S^0,V^0,E^0,I^0,R^0)$ and $E_1$ = $(S^{*},V^{*},E^{*},I^{*},R^{*})$. We get the disease free equilibrium points for system (2) as $E_0$ = ($\frac{A}{{\theta }_1}$, $\frac{zA}{{\theta }_1\mu }$ 0, 0, 0). Also for endemic equilibrium we obtain $S^{*}$ = $\frac{\eta A{\theta }_4+{\alpha }_0\gamma \eta I^{*}}{{\theta }_4(\eta {\theta }_1+{\delta }_1\eta I^{*}+{\delta }_2{\theta }_3\eta I^{*})}$, $V^{*}$ = $\frac{\eta z{\theta }_4(\frac{\eta A{\theta }_4+{\alpha }_0\gamma \eta I^{*}}{{\theta }_4(\eta {\theta }_1+{\delta }_1\eta I^{*}+{\delta }_2{\theta }_3\eta I^{*})}+\eta {\alpha }_1\gamma I^{*})}{{\theta }_4({\delta }_3\eta I^{*}+{\delta }_4{\theta }_3{I}^{*}+\eta \mu )}$, $E^{*}$ = $\frac{{\theta }_3{I}^{*}}{\eta }$ and $R^{*}$ = $\frac{\gamma I^{*}}{{\theta }_4}$. Here $I^{*}$ will be in the form of $g\left(I^{*}\right)=l_1{I^{*}}^3+l_2{I^{*}}^2+l_3{I}^{*}+l_4=0$. By using Descartes rule of signs [21], we can able to the positive roots of $I^{*}$.

6.2. Basic Reproduction Number (BRN)

The spectral radius of the next generation matrix determines the $R_0$[4] for the system (2). Therefore, we have $R_0=\frac{\eta S{\delta }_1+\eta V{\delta }_3+S{\delta }_2{\theta }_3+{\theta }_{3}V{\delta }_4}{{\theta }_2{\theta }_3}$. By substituting disease free equilibrium points we get,

$$\begin{array}{c}\hfill R_0=\frac{\eta \left(\frac{A}{{\theta }_1}\right){\delta }_1+\eta \left(\frac{zA}{{\theta }_1\mu }\right){\delta }_3+\left(\frac{A}{{\theta }_1}\right){\delta }_2{\theta }_3+{\theta }_3\left(\frac{zA}{{\theta }_1\mu }\right){\delta }_4}{{\theta }_2{\theta }_3}.\end{array}$$

where,

${\theta }_1=\mu +z,$

${\theta }_2=\eta +\mu$,

${\theta }_3=\gamma +\mu +{\mu }_1$.

6.3. Local stability

The Jacobian matrix for (2) is presented by,

$$J=\left(\begin{array}{ccccc}-{\delta }_{1}I-{\delta }_{2}E-{\theta }_1& 0& -{\delta }_{2}S& -{\delta }_{1}S& {\alpha }_0\\ z& -{\delta }_{3}I-{\delta }_{4}E-\mu & -{\delta }_{4}V& -{\delta }_{3}V& {\alpha }_1\\ {\delta }_{1}I+{\delta }_{2}E& {\delta }_{3}I+{\delta }_{4}E& {\delta }_{2}S+{\delta }_{4}V-{\theta }_2& {\delta }_{3}V+{\delta }_{1}S& 0\\ 0& 0& \eta & -{\theta }_3& 0\\ 0& 0& 0& \gamma & -{\theta }_4\end{array}\right).$$

Theorem 6.1. 

The disease free equilibrium point $E_0$ of the fractional order model (2) is stable if $R_0<1$.

Proof. 

By the approach of Matignon’s condition in [20], $E_0$ of the fractional order (2) is locally asymptotically stable ⇔ all the eigen values of $J_0$ should fulfil |arg(${\lambda }_i$)|>$\frac{a\pi }{2}$.

□

The Jacobi matrix at $E_0$ is obtained as follows,

$$J_0=\left(\begin{array}{ccccc}-{\theta }_1& 0& -\frac{{\delta }_{2}A}{{\theta }_1}& -\frac{{\delta }_{1}A}{{\theta }_1}& {\alpha }_0\\ z& -\mu & -\frac{{\delta }_{4}zA}{{\theta }_1\mu }& -\frac{{\delta }_{3}zA}{{\theta }_1\mu }& {\alpha }_1\\ 0& 0& \frac{{\delta }_{2}A}{{\theta }_1}+\frac{{\delta }_{4}zA}{{\theta }_1\mu }-{\theta }_2& \frac{{\delta }_{3}zA}{{\theta }_1\mu }+\frac{{\delta }_{1}A}{{\theta }_1}& 0\\ 0& 0& \eta & -{\theta }_3& 0\\ 0& 0& 0& \gamma & -{\theta }_4\end{array}\right)$$

Clearly, from the above matrix we can say, ${\lambda }_1=-{\theta }_1,{\lambda }_2=-\mu ,{\lambda }_5=-{\theta }_4$ are the three eigen values. Here ${\lambda }_3$ and ${\lambda }_4$ can be computed from the below matrix.

$$\left(\begin{array}{cc}\frac{{\delta }_{2}A}{{\theta }_1}+\frac{{\delta }_{4}zA}{{\theta }_1\mu }-{\theta }_2& \frac{{\delta }_{3}zA}{{\theta }_1\mu }+\frac{{\delta }_{1}A}{{\theta }_1}\\ \eta & -{\theta }_3\end{array}\right)$$

They are ${\lambda }_3=y_1+\sqrt{y_2+y_3}$ and ${\lambda }_4=y_1-\sqrt{y_2+y_3}$, where,

$$\begin{array}{cc}\hfill y_1& =A{\delta }_2\mu -{\theta }_1{\theta }_2\mu -{\theta }_1{\theta }_3\mu +A{\delta }_{4}z,\hfill \\ \hfill y_2& =A^2{\delta }_2^2{\mu }^2+{\theta }_1^2{\theta }_2^2{\mu }^2+2{\theta }_1^2{\theta }_2{\theta }_3{\mu }^2+{\theta }_1^2{\theta }_3^2{\mu }^2+2A^2{\delta }_2{\delta }_4\mu z+A^2{\delta }_4^2{z}^2\hfill \\ \hfill & -2A{\delta }_2{\theta }_1{\theta }_2{\mu }^2-2A{\delta }_2{\theta }_1{\theta }_3{\mu }^2-2A{\delta }_4{\theta }_1{\theta }_2\mu z-2A{\delta }_4{\theta }_1{\theta }_3\mu z,\hfill \\ \hfill y_3& =4(A{\delta }_2{\theta }_1{\theta }_3{\mu }^2+{\delta }_{1}A{\theta }_1\eta {\mu }^2+A{\delta }_4{\theta }_1{\theta }_3\mu z+A{\delta }_3{\theta }_1\eta \mu z-{\theta }_1^2{\theta }_2{\theta }_3{\mu }^2).\hfill \end{array}$$

$$\begin{array}{cc}\hfill y_1& <0\mathrm{when}{\theta }_1{\theta }_2\mu +{\theta }_1{\theta }_3\mu >A{\delta }_2\mu +A{\delta }_{4}z,\hfill \\ \hfill y_2& >0\mathrm{when}{A}^2{\delta }_2^2{\mu }^2+{\theta }_1^2{\theta }_2^2{\mu }^2+2{\theta }_1^2{\theta }_2{\theta }_3{\mu }^2+{\theta }_1^2{\theta }_3^2{\mu }^2+2A^2{\delta }_2{\delta }_4\mu z\hfill \\ \hfill & +A^2{\delta }_4^2{z}^2>2A{\delta }_2{\theta }_1{\theta }_2{\mu }^2+2A{\delta }_2{\theta }_1{\theta }_3{\mu }^2+2A{\delta }_4{\theta }_1{\theta }_2\mu z+2A{\delta }_4{\theta }_1{\theta }_3\mu z,\hfill \\ \hfill y_3& >0\mathrm{when}4(A{\delta }_2{\theta }_1{\theta }_3{\mu }^2+{\delta }_{1}A{\theta }_1\eta {\mu }^2+A{\delta }_4{\theta }_1{\theta }_3\mu z+A{\delta }_3{\theta }_1\eta \mu z)>4{\theta }_1^2{\theta }_2{\theta }_3{\mu }^2.\hfill \end{array}$$

Therefore ${\lambda }_3$ and ${\lambda }_4$ have negative real numbers with the above conditions. From this clearly we can conclude that all eigenvalues are negative and does fulfil the Matignon’s condition [20]. Hence, $E_0$= ($\frac{A}{{\theta }_1}$, $\frac{zA}{{\theta }_1\mu }$ 0, 0, 0) is stable whenever BRN is less than one.

Table 2. Necessary components for the stability of the equilibrium points.

Eigen value	Sign	Conditions	Stability
${\lambda }_3$	−	$y_1<0$, $y_2>0$ and $y_3>0$ and hence $y_1>\sqrt{y_2+y_3}$	Stable
${\lambda }_4$	−	$y_1<0$, $y_2>0$ and $y_3>0$ and hence $y_1-\sqrt{y_2+y_3}<0$	Stable

Table 2. Necessary components for the stability of the equilibrium points.

Eigen value	Sign	Conditions	Stability
${\lambda }_3$	−	$y_1<0$, $y_2>0$ and $y_3>0$ and hence $y_1>\sqrt{y_2+y_3}$	Stable
${\lambda }_4$	−	$y_1<0$, $y_2>0$ and $y_3>0$ and hence $y_1-\sqrt{y_2+y_3}<0$	Stable

Theorem 6.2. 

When $R_0>1$, the endemic equilibrium $E_1$ of the model (2) is asymptotically stable point.

Proof. 

By using [20], the endemic equilibrium $E_1$ of (2) is locally asymptotically stable ⇔ all the eigenvalues of below mentioned matrix should fulfil |arg(${\lambda }_i$)|>$\frac{a\pi }{2}$.

For (2), the Jacobi matrix at $E_1$ is obtained as,

$$J_1=\left(\begin{array}{ccccc}-{\delta }_1{I}^{*}-{\delta }_2{E}^{*}-{\theta }_1& 0& -{\delta }_2{S}^{*}& -{\delta }_1{S}^{*}& {\alpha }_0\\ z& -{\delta }_3{I}^{*}-{\delta }_4{E}^{*}-\mu & -{\delta }_4{V}^{*}& -{\delta }_3{V}^{*}& {\alpha }_1\\ {\delta }_1{I}^{*}+{\delta }_2{E}^{*}& {\delta }_3{I}^{*}+{\delta }_4{E}^{*}& {\delta }_2{S}^{*}+{\delta }_4{V}^{*}-{\theta }_2& {\delta }_3{V}^{*}+{\delta }_1{S}^{*}& 0\\ 0& 0& \eta & -{\theta }_3& 0\\ 0& 0& 0& \gamma & -{\theta }_4\end{array}\right)$$

Let us consider,

$b_{11}=-{\delta }_1{I}^{*}-{\delta }_2{E}^{*}-{\theta }_1$,

$b_{13}=-{\delta }_2{S}^{*}$,

$b_{14}=-{\delta }_1{S}^{*}$,

$b_{15}={\alpha }_0$,

$b_{21}=z$,

$b_{22}=-{\delta }_3{I}^{*}-{\delta }_4{E}^{*}-\mu$,

$b_{23}=-{\delta }_4{V}^{*}$,

$b_{24}=-{\delta }_3{V}^{*}$,

$b_{25}={\alpha }_1$,

$b_{31}={\delta }_1{I}^{*}+{\delta }_2{E}^{*}$,

$b_{32}={\delta }_3{I}^{*}+{\delta }_4{E}^{*}$,

$b_{33}={\delta }_2{S}^{*}+{\delta }_4{V}^{*}-{\theta }_2$,

$b_{34}={\delta }_3{V}^{*}+{\delta }_1{S}^{*}$,

$b_{43}=\eta$,

$b_{44}=-{\theta }_3$,

$b_{54}=\gamma$,

$b_{55}=-{\theta }_4$.

The characteristic equation of the $J_1$ is,

$$\Rightarrow {\lambda }^5+x_1{\lambda }^4+x_2{\lambda }^3+x_3{\lambda }^2+x_4\lambda +x_5=0,$$

where,

$$\begin{array}{cc}\hfill x_1& =-b_{11}-b_{22}-b_{33}-b_{44}-b_{55},\hfill \\ \hfill x_2& =b_{44}{b}_{55}+b_{33}{b}_{55}+b_{22}{b}_{55}+b_{11}{b}_{55}+b_{33}{b}_{44}+b_{22}{b}_{44}+b_{11}{b}_{44}+b_{22}{b}_{33}+b_{11}{b}_{33}+b_{11}{b}_{22}\hfill \\ \hfill & -b_{34}{b}_{43}-b_{23}{b}_{32}-b_{13}{b}_{31},\hfill \\ \hfill x_3& =b_{34}{b}_{43}{b}_{55}+b_{23}{b}_{32}{b}_{55}+b_{13}{b}_{31}{b}_{55}+b_{23}{b}_{32}{b}_{44}+b_{13}{b}_{31}{b}_{44}+b_{22}{b}_{34}{b}_{43}+b_{11}{b}_{34}{b}_{43}\hfill \\ \hfill & +b_{11}{b}_{23}{b}_{32}+b_{13}{b}_{22}{b}_{31}-b_{33}{b}_{44}{b}_{55}-b_{22}{b}_{44}{b}_{55}-b_{11}{b}_{44}{b}_{55}-b_{22}{b}_{33}{b}_{55}-b_{11}{b}_{33}{b}_{55}\hfill \\ \hfill & -b_{11}{b}_{22}{b}_{55}-b_{22}{b}_{33}{b}_{44}-b_{11}{b}_{33}{b}_{44}-b_{11}{b}_{22}{b}_{44}-b_{24}{b}_{32}{b}_{43}-b_{14}{b}_{31}{b}_{43}-b_{11}{b}_{22}{b}_{33}\hfill \\ & -b_{13}{b}_{21}{b}_{32},\hfill \\ \hfill x_4& =b_{22}{b}_{33}{b}_{44}{b}_{55}+b_{11}{b}_{33}{b}_{44}{b}_{55}+b_{11}{b}_{22}{b}_{44}{b}_{55}+b_{24}{b}_{32}{b}_{43}{b}_{55}+b_{14}{b}_{31}{b}_{43}{b}_{55}+b_{11}{b}_{22}{b}_{33}{b}_{55}\hfill \\ \hfill & +b_{13}{b}_{21}{b}_{32}{b}_{55}+b_{11}{b}_{22}{b}_{33}{b}_{44}+b_{11}{b}_{23}{b}_{32}{b}_{44}+b_{13}{b}_{21}{b}_{32}{b}_{44}+b_{11}{b}_{24}{b}_{32}{b}_{43}+b_{14}{b}_{22}{b}_{31}{b}_{43}\hfill \\ \hfill & -b_{23}{b}_{32}{b}_{44}{b}_{55}-b_{13}{b}_{31}{b}_{44}{b}_{55}-b_{22}{b}_{34}{b}_{43}{b}_{55}-b_{11}{b}_{34}{b}_{43}{b}_{55}-b_{11}{b}_{23}{b}_{32}{b}_{55}-b_{13}{b}_{22}{b}_{31}{b}_{55}\hfill \\ \hfill & -b_{25}{b}_{32}{b}_{43}{b}_{54}-b_{15}{b}_{31}{b}_{43}{b}_{54}-b_{13}{b}_{22}{b}_{31}{b}_{44}-b_1{b}_{22}{b}_{34}{b}_{43}-b_{14}{b}_{21}{b}_{32}{b}_{43},\hfill \\ \hfill x_5& =b_{11}{b}_{23}{b}_{32}{b}_{44}{b}_{55}+b_{13}{b}_{22}{b}_{31}{b}_{44}{b}_{55}+b_{11}{b}_{22}{b}_{34}{b}_{43}{b}_{55}+b_{14}{b}_{21}{b}_{32}{b}_{43}{b}_{55}+b_{11}{b}_{25}{b}_{32}{b}_{43}{b}_{54}\hfill \\ \hfill & +b_{15}{b}_{22}{b}_{31}{b}_{43}{b}_{54}-b_{11}{b}_{22}{b}_{33}{b}_{44}{b}_{55}-b_{13}{b}_{21}{b}_{32}{b}_{44}{b}_{55}-b_{11}{b}_{24}{b}_{32}{b}_{43}{b}_{55}-b_{14}{b}_{22}{b}_{31}{b}_{43}{b}_{55}\hfill \\ & -b_{15}{b}_{21}{b}_{32}{b}_{43}{b}_{54}.\hfill \end{array}$$

Here we consider,

$$\begin{array}{cc}\hfill •x_1& >0,\mathrm{if}{b}_{11}+b_{22}+b_{33}+b_{44}+b_{55}<0,\hfill \\ \hfill •x_2& >0,\mathrm{if}{b}_{11}{b}_{33}+b_{11}{b}_{22}>b_{34}{b}_{43}+b_{23}{b}_{32}+b_{13}{b}_{31},\hfill \\ \hfill •x_3& >0,\mathrm{if}{b}_{11}{b}_{34}{b}_{43}+b_{11}{b}_{23}{b}_{32}+b_{13}{b}_{22}{b}_{31}>b_{33}{b}_{44}{b}_{55}+b_{22}{b}_{44}{b}_{55}+b_{11}{b}_{44}{b}_{55}+b_{22}{b}_{33}{b}_{55}\hfill \\ \hfill & +b_{11}{b}_{33}{b}_{55}+b_{11}{b}_{22}{b}_{55}+b_{22}{b}_{33}{b}_{44}+b_{11}{b}_{33}{b}_{44}+b_{11}{b}_{22}{b}_{44}+b_{24}{b}_{32}{b}_{43}+b_{14}{b}_{31}{b}_{43}\hfill \\ \hfill & +b_{11}{b}_{22}{b}_{33}+b_{13}{b}_{21}{b}_{32},\hfill \\ \hfill •x_4& >0,\mathrm{if}{b}_{11}{b}_{24}{b}_{32}{b}_{43}+b_{14}{b}_{22}{b}_{31}{b}_{43}>b_{23}{b}_{32}{b}_{44}{b}_{55}+b_{13}{b}_{31}{b}_{44}{b}_{55}+b_{22}{b}_{34}{b}_{43}{b}_{55}\hfill \\ \hfill & +b_{11}{b}_{34}{b}_{43}{b}_{55}+b_{11}{b}_{23}{b}_{32}{b}_{55}+b_{13}{b}_{22}{b}_{31}{b}_{55}+b_{25}{b}_{32}{b}_{43}{b}_{54}+b_{15}{b}_{31}{b}_{43}{b}_{54}+b_{13}{b}_{22}{b}_{31}{b}_{44}\hfill \\ & +b_1{b}_{22}{b}_{34}{b}_{43}+b_{14}{b}_{21}{b}_{32}{b}_{43},\hfill \\ \hfill •x_5& >0,\mathrm{if}{b}_{11}{b}_{23}{b}_{32}{b}_{44}{b}_{55}>b_{11}{b}_{22}{b}_{33}{b}_{44}{b}_{55}+b_{13}{b}_{21}{b}_{32}{b}_{44}{b}_{55}+b_{11}{b}_{24}{b}_{32}{b}_{43}{b}_{55}\hfill \\ \hfill & +b_{14}{b}_{22}{b}_{31}{b}_{43}{b}_{55}+b_{15}{b}_{21}{b}_{32}{b}_{43}{b}_{54}.\hfill \end{array}$$

From the above conditions, we can say that, $E_1=(S^{*},E^{*},I^{*},Q^{*},R^{*},D^{*})$ of equation (2) is stable whenever $R_0>1$.

□

7. Numerical Simulation

In this, we derive general fractional order Euler method approach for (2). By reformulating, we get

$$\begin{array}{cc}\hfill {\mathbb{D}}_t^{\alpha }\left[S\left(t\right)\right]& ={\tau }_1(t,S\left(t\right)),\hfill \\ \hfill {\mathbb{D}}_t^{\alpha }\left[V\left(t\right)\right]& ={\tau }_2(t,V\left(t\right)),\hfill \\ \hfill {\mathbb{D}}_t^{\alpha }\left[E\left(t\right)\right]& ={\tau }_3(t,E\left(t\right)),\hfill \\ \hfill {\mathbb{D}}_t^{\alpha }\left[I\left(t\right)\right]& ={\tau }_4(t,I\left(t\right)),\hfill \\ \hfill {\mathbb{D}}_t^{\alpha }\left[R\left(t\right)\right]& ={\tau }_5(t,R\left(t\right)).\hfill \end{array}$$

(13)

where,

$$\begin{array}{cc}\hfill {\tau }_1(t,S\left(t\right))& =A-{\delta }_{1}SI-{\delta }_{2}SE-\mu S-zS+{\alpha }_{0}R,\hfill \\ \hfill {\tau }_2(t,V\left(t\right))& =zS-{\delta }_{3}VI-{\delta }_{4}VE-\mu V+{\alpha }_{1}R,\hfill \\ \hfill {\tau }_3(t,E\left(t\right))& ={\delta }_{1}SI+{\delta }_{2}SE+{\delta }_{3}VI+{\delta }_{4}VE-\eta E-\mu E,\hfill \\ \hfill {\tau }_4(t,I\left(t\right))& =\eta E-\gamma I-\mu I-{\mu }_{1}I,\hfill \\ \hfill {\tau }_5(t,R\left(t\right))& =\gamma I-\mu R-{\alpha }_{0}R-{\alpha }_{1}R.\hfill \end{array}$$

(14)

From first equation of (2),

$$\begin{array}{cc}\hfill {\mathbb{D}}_t^{\alpha }\left[S\left(t\right)\right]& ={\tau }_1(t,S\left(t\right)),\hfill \\ \hfill S\left(0\right)& =S_0,t>0.\hfill \end{array}$$

(15)

Let $[0,d]$ be the collection of points for which we wish to discover the solution of (15). In fact, we are unable to assess $S\left(t\right)$ that will be corresponding to (15). Instead, a collection of $(t_r,t_r+1)$ is formed, and the points are used in our iterative approach. From $t_r=rh$, $r=0,1,2,3,...k$, we partition the interval $[0,d]$ into k subintervals $[t_r,t_r+1]$ with equal width $h=\frac{d}{k}$. Assume $S\left(t\right)$, ${\mathbb{D}}_t^{\alpha }\left[S\left(t\right)\right]$ and ${\mathbb{D}}_t^{2\alpha }\left[S\left(t\right)\right]$ are continuous on $[0,d]$. By using the generalised Taylor formula, extend $S\left(t\right)$ to $t=t_0=0.$ There is a value $c_1$ for each value t, therefore

$$\begin{array}{cc}\hfill S\left(t\right)& =S\left(t_0\right)+{\mathbb{D}}_t^{\alpha }\left[S\left(t\right)\right]t_0\frac{t^{\alpha }}{\Gamma (\alpha +1)}+{\mathbb{D}}_t^{2\alpha }\left[S\left(t\right)\right]c_1\frac{t^{2\alpha }}{\Gamma (2\alpha +1)}\hfill \end{array}$$

(16)

Substitute ${\mathbb{D}}_t^{\alpha }\left[S\left(t\right)\right]t_0={\tau }_1(t_0,S\left(t_0\right))$ and $h=t_1$ in (16)

$$\begin{array}{cc}\hfill S\left(t_1\right)& =S\left(t_0\right)+{\tau }_1(t_0,S\left(t_0\right))\frac{h^{\alpha }}{\Gamma (\alpha +1)}+{\mathbb{D}}_t^{2\alpha }\left[S\left(t\right)\right]c_1\frac{h^{2\alpha }}{\Gamma (2\alpha +1)},\hfill \end{array}$$

If h is small, thus

$$\begin{array}{cc}\hfill S\left(t_1\right)& =S\left(t_0\right)+{\tau }_1(t_0,S\left(t_0\right))\frac{h^{\alpha }}{\Gamma (\alpha +1)},\hfill \end{array}$$

(17)

A general formula, $t_{r+1}=t_r+h$ can be written as

$$\begin{array}{cc}\hfill S\left(t_{r+1}\right)& =S\left(t_r\right)+{\tau }_1(t_r,S\left(t_r\right))\frac{h^{\alpha }}{\Gamma (\alpha +1)}.\hfill \end{array}$$

(18)

Using the fractional integral on (15), we obtain

$$\begin{array}{cc}\hfill \left[S\right(t\left)\right]& =S\left(0\right)+I^{\alpha }\left[{\tau }_1(t,S\left(t\right))\right]\hfill \end{array}$$

(19)

To obtain $(t_1,S\left(t_1\right)),$ we replace $t=t_1$ into (19), then

$$\begin{array}{cc}\hfill S\left(t_1\right)& =S\left(0\right)+\left(I^{\alpha }\left[{\tau }_1(t,S\left(t\right))\right]\right)\left(t_1\right).\hfill \end{array}$$

(20)

With the help modified trapezoidal rule, we approximate $\left(I^{\alpha }\left[{\tau }_1(t,S\left(t\right))\right]\right)\left(t_1\right)$ with $h=t_1-t_0$ then (20) becomes

$$\begin{array}{cc}\hfill S\left(t_1\right)& =S\left(0\right)+\frac{\alpha h^{\alpha }\left[{\tau }_1(t_0,S\left(t_0\right))\right]}{\Gamma (\alpha +2)}+\frac{h^{\alpha }\left[{\tau }_1(t_1,S\left(t_1\right))\right]}{\Gamma (\alpha +2)}\hfill \end{array}$$

(21)

The expression $S\left(t_1\right)$ appears in the formula on the RHS of (21). As a result, we estimate $S\left(t_1\right)$. For this purpose, the fractional Euler approach can be used. Substitute (17) for (21).

$$\begin{array}{cc}\hfill S\left(t_1\right)& =S\left(0\right)+\frac{\alpha h^{\alpha }\left[{\tau }_1(t_0,S\left(t_0\right))\right]}{\Gamma (\alpha +2)}+\frac{h^{\alpha }\left[{\tau }_1(t_1,S\left(t_0\right))\right]+\frac{h^{\alpha }}{\Gamma (\alpha +1)}{\tau }_1(t_0,S\left(t_0\right))}{\Gamma (\alpha +2)}\hfill \end{array}$$

The technique is iterated until a sequence of points that closely approximate the solution is obtained. Therefore, the general formula for our approach can be written as,

$$\begin{array}{cc}\hfill S\left(t_r\right)& =S\left(0\right)+\frac{h^{\alpha }}{\Gamma (\alpha +2)}[{(r-1)}^{\alpha +1}-(r-\alpha -1)r^{\alpha }]{\tau }_1(t_0,S\left(t_0\right))\hfill \\ & +\frac{h^{\alpha }}{\Gamma (\alpha +2)}\sum _{i=1}^{r-1}[{(r-i+1)}^{\alpha +1}-2{(r-1)}^{\alpha +1}+{(r-i-1)}^{\alpha +1}]{\tau }_1(t_i,S\left(t_i\right))\hfill \\ & +\frac{h^{\alpha }}{\Gamma (\alpha +2)}{\tau }_1[t_r,S\left(t_{r-1}\right)+\frac{h^{\alpha }}{\Gamma (\alpha +1)}{\tau }_1(t_{r-1},S\left(t_{r-1}\right))].\hfill \end{array}$$

(22)

Using the same method, we construct numerical schemes for the model’s additional compartments.

$$\begin{array}{cc}\hfill V\left(t_r\right)& =V\left(0\right)+\frac{h^{\alpha }}{\Gamma (\alpha +2)}[{(r-1)}^{\alpha +1}-(r-\alpha -1)r^{\alpha }]{\tau }_2(t_0,V\left(t_0\right))\hfill \\ & +\frac{h^{\alpha }}{\Gamma (\alpha +2)}\sum _{i=1}^{r-1}[{(r-i+1)}^{\alpha +1}-2{(r-1)}^{\alpha +1}+{(r-i-1)}^{\alpha +1}]{\tau }_2(t_i,V\left(t_i\right))\hfill \\ & +\frac{h^{\alpha }}{\Gamma (\alpha +2)}{\tau }_2[t_r,V\left(t_{r-1}\right)+\frac{h^{\alpha }}{\Gamma (\alpha +1)}{\tau }_2(t_{r-1},V\left(t_{r-1}\right))],\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill E\left(t_r\right)& =E\left(0\right)+\frac{h^{\alpha }}{\Gamma (\alpha +2)}[{(r-1)}^{\alpha +1}-(r-\alpha -1)r^{\alpha }]{\tau }_3(t_0,E\left(t_0\right))\hfill \\ & +\frac{h^{\alpha }}{\Gamma (\alpha +2)}\sum _{i=1}^{r-1}[{(r-i+1)}^{\alpha +1}-2{(r-1)}^{\alpha +1}+{(r-i-1)}^{\alpha +1}]{\tau }_3(t_i,E\left(t_i\right))\hfill \\ & +\frac{h^{\alpha }}{\Gamma (\alpha +2)}{\tau }_3[t_r,E\left(t_{r-1}\right)+\frac{h^{\alpha }}{\Gamma (\alpha +1)}{\tau }_3(t_{r-1},E\left(t_{r-1}\right))],\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill I\left(t_r\right)& =I\left(0\right)+\frac{h^{\alpha }}{\Gamma (\alpha +2)}[{(r-1)}^{\alpha +1}-(r-\alpha -1)r^{\alpha }]{\tau }_4(t_0,I\left(t_0\right))\hfill \\ & +\frac{h^{\alpha }}{\Gamma (\alpha +2)}\sum _{i=1}^{r-1}[{(r-i+1)}^{\alpha +1}-2{(r-1)}^{\alpha +1}+{(r-i-1)}^{\alpha +1}]{\tau }_4(t_i,I\left(t_i\right))\hfill \\ & +\frac{h^{\alpha }}{\Gamma (\alpha +2)}{\tau }_4[t_r,I\left(t_{r-1}\right)+\frac{h^{\alpha }}{\Gamma (\alpha +1)}{\tau }_4(t_{r-1},I\left(t_{r-1}\right))],\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill R\left(t_r\right)& =R\left(0\right)+\frac{h^{\alpha }}{\Gamma (\alpha +2)}[{(r-1)}^{\alpha +1}-(r-\alpha -1)r^{\alpha }]{\tau }_5(t_0,R\left(t_0\right))\hfill \\ & +\frac{h^{\alpha }}{\Gamma (\alpha +2)}\sum _{i=1}^{r-1}[{(r-i+1)}^{\alpha +1}-2{(r-1)}^{\alpha +1}+{(r-i-1)}^{\alpha +1}]{\tau }_5(t_i,R\left(t_i\right))\hfill \\ & +\frac{h^{\alpha }}{\Gamma (\alpha +2)}{\tau }_5[t_r,R\left(t_{r-1}\right)+\frac{h^{\alpha }}{\Gamma (\alpha +1)}{\tau }_5(t_{r-1},R\left(t_{r-1}\right))].\hfill \end{array}$$

(26)

The numerical scheme described above is employed to obtain the numerical results. Furthermore, assumptions are introduced for the initial states of variables and input parameters to streamline the calculations. Figure 1, Figure 2, Figure 3 and Figure 4 show the results of the most effective cases. Figure 1a shows the local stability of $SVEIR$ individuals for the corresponding disease free equilibrium. Let us consider, $A=1000$, ${\delta }_1=0.0009$, ${\delta }_2=0.0008$, ${\delta }_3=0.009$, ${\delta }_4=0.0008$, $\mu =0.9$, ${\mu }_1=0.042$, $z=0.099$, ${\alpha }_0=0.05$, ${\alpha }_1=0.005$, $\gamma =0.0052$ and $\eta =0.7$. For this set of parameter values we obtain BRN $R_0<1$. So from this we clearly sum up that our $E_0$ is locally stable. For different set of parameter we get $R_0>1$, which was illustrated in the Figure 1b. So from this we sum up that our $E_1$ is locally stable. Figure 2a,b shows the phase diagram for vaccinated, infected and recovered population of both DFE and EE. Next, we checked the influence of progression rate for infected individuals, which was depicted in the Figure 3. For this simulation we consider progression rate as 0.10, 0.25, 0.50 respectively. Then we also checked the influence of recovery rate due to treatment for recovered individuals, which was shown in Figure 4. For this simulation we consider recovery rate as 0.10, 0.25, 0.50 respectively. Our analysis predicts that how the input factors influence the output of the system of HFMD transmission.

8. Conclusion

In this study, we focus on investigating the applicability of the CFO mathematical model in understanding the spread dynamics of Hand Foot Mouth Disease (HFMD). The focus of our investigation is to determine whether the system is non negative and bounded, which allowed us to analyse the stability of the model. Through nonlinear analysis, we demonstrated the existence and uniqueness of the model, which enhances its reliability in predicting disease spread. Furthermore, we calculated the numerical solution of the suggested system using the fractional Euler’s approach. This approach is widely used to solve fractional differential equations, and it provides accurate results that are vital in disease modelling. By utilizing this method, we were able to obtain precise solutions for the model, which allowed us to make better predictions about HFMD transmission. Overall, our study demonstrates that the Caputo fractional order mathematical model is a reliable tool for understanding the spread of HFMD. The establishment of positive invariance and boundedness of the system, coupled with the use of nonlinear analysis and numerical solution, enhances the accuracy of the model, which can be used to inform policies aimed at controlling and preventing the spread of HFMD.