Preprint
Article

Fractional Order Mathematical Modelling of HFMD Transmission via Caputo Derivative

Altmetrics

Downloads

99

Views

50

Comments

0

A peer-reviewed article of this preprint also exists.

Submitted:

14 December 2023

Posted:

15 December 2023

You are already at the latest version

Alerts
Abstract
In this paper, we propose to explore a nonlinear fractional mathematical model of HFMD- Hand Foot Mouth Disease with vaccinated compartment. Firstly, we prove the non negativity and boundedness for the fractional order dynamical model. The system’s existence and uniqueness are studied using the formulation of the Caputo derivative operator. The fixed-point approach was used to get results indicating the presence of at least one solution. In numerical simulation, we obtained the approximate solution by using fractional Euler method.
Keywords: 
Subject: Computer Science and Mathematics  -   Applied Mathematics

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 S V E I R , 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 α 0 , and that after recovery, some people will take vaccination and therefore enter the vaccinated compartment at a rate of α 1 . As a result of the aforesaid assumptions, a basic nonlinear mathematical model was developed in [21]. That is,
d S d t = A δ 1 S I δ 2 S E μ S z S + α 0 R , d V d t = z S δ 3 V I δ 4 V E μ V + α 1 R , d E d t = δ 1 S I + δ 2 S E + δ 3 V I + δ 4 V E η E μ E , d I d t = η E γ I μ I μ 1 I , d R d t = γ I μ R α 0 R α 1 R .
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.
  C D t α [ S ( t ) ] = A δ 1 S I δ 2 S E μ S z S + α 0 R ,   C D t α [ V ( t ) ] = z S δ 3 V I δ 4 V E μ V + α 1 R ,   C D t α [ E ( t ) ] = δ 1 S I + δ 2 S E + δ 3 V I + δ 4 V E η E μ E ,   C D t α [ I ( t ) ] = η E γ I μ I μ 1 I ,   C D t α [ R ( t ) ] = γ I μ R α 0 R α 1 R .
Here, 0 < α ≤ 1 and S ( 0 ) = S 0 , V ( 0 ) = V 0 , E ( 0 ) = E 0 , I ( 0 ) = I 0 and R ( 0 ) = 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 ∈ C m , then   C D t α in (m-1,m) where m ∈ N is provided by,
  C D t α ( g ( t ) ) = 1 Γ m α 0 t ( t s ) ( m α 1 ) g ( m ) ( s ) d s .
Here,   C D t α ( g ( t ) ) g ( t ) as α tends to 1.
Definition 3.2. 
For α> 0, the corresponding integral of g : R + R is expressed as follows,
I t α ( g ( t ) ) = 1 Γ α 0 t ( t s ) ( α 1 ) g ( s ) d s ,
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 d N d t = d S d t + d V d t + d E d t + d I d t + d R d t . Also let consider μ > 0 .
  C D t α [ N ( t ) ] + μ N ( t ) = A δ 1 S I δ 2 S E μ S z S + α 0 R + z S δ 3 V I δ 4 V E μ V + α 1 R + δ 1 S I + δ 2 S E + δ 3 V I + δ 4 V E η E μ E + η E γ I μ I μ 1 I + γ I μ R α 0 R α 1 R + μ S + μ V + μ E + μ I + μ R = A μ 1 I A
By using standard comparison theorem [19] for (2), we get,
N ( t ) N ( 0 ) E α μ ( t ) α + A ( t ) α E α , α + 1 μ ( t ) α
where the Mittag Leffler function is represented by E α . As stated by [19], we have
N ( t ) A μ , t α .
Hence, all solutions of equation (2) are restricted to the domain Ω ,
Ω = { ( S , V , E , I , R ) R + 5 | N ( t ) A μ + ϵ 0 , f o r a n y ϵ 0 > 0 } .
Now, we show the solutions of (2) are positive. From first equation of (2), we have,
  C D t α [ S ( t ) ] = A δ 1 S I δ 2 S E μ S z S + α 0 R [ δ 1 I + δ 2 E + μ + z ] S [ ( δ 1 + δ 2 ) A μ + μ + z ] S d 1 S .
where d 1 = ( δ 1 + δ 2 ) A μ + μ + z .
As per the traditional comparison theorem [19] and the non negativity of Mitag Leffler function, we have, E α , 1 ( t ) > 0 for any α ( 0 , 1 ) ,
S S ( 0 ) E α , 1 ( d 1 t α ) S 0 .
From (2), we get
  C D t α [ V ( t ) ] = z S δ 3 V I δ 4 V E μ V + α 1 R [ δ 3 I + δ 4 E + μ ] V [ ( δ 3 + δ 4 ) A μ + μ ] V d 2 V .
where d 2 = ( δ 3 + δ 4 ) A μ + μ .
Therefore,
V V ( 0 ) E α , 1 ( d 2 t α ) V 0 .
For third equation of (2), we get
  C D t α [ E ( t ) ] = δ 1 S I + δ 2 S E + δ 3 V I + δ 4 V E η E μ E [ δ 2 S δ 4 V + η + μ ] E [ δ 2 A μ δ 4 A μ + η + μ ] E d 3 E .
Here d 3 = δ 2 A μ δ 4 A μ + η + μ .
E E ( 0 ) E α , 1 ( d 3 t α ) E 0 .
For fourth equation of (2), we get
  C D t α [ I ( t ) ] = η E γ I μ I μ 1 I [ γ + μ + μ 1 ] I d 4 I .
where d 4 = γ + μ + μ 1 .
I I ( 0 ) E α , 1 ( d 4 t α ) I 0 .
For fifth equation of (2), we get
  C D t α [ R ( t ) ] = γ I μ R α 0 R α 1 R [ μ + α 0 + α 1 ] I d 5 R .
where d 5 = μ + α 0 + α 1 .
Therefore,
R R ( 0 ) E α , 1 ( d 5 t α ) R 0 .
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).
S ( t ) = S ( 0 ) +   C I 0 α S ( t ) { A δ 1 S I δ 2 S E μ S z S + α 0 R } , V ( t ) = V ( 0 ) +   C I 0 α V ( t ) { z S δ 3 V I δ 4 V E μ V + α 1 R } , E ( t ) = E ( 0 ) +   C I 0 α E ( t ) { δ 1 S I + δ 2 S E + δ 3 V I + δ 4 V E η E μ E } , I ( t ) = I ( 0 ) +   C I 0 α I ( t ) { η E γ I μ I μ 1 I } , R ( t ) = R ( 0 ) +   C I 0 α R ( t ) { γ I μ R α 0 R α 1 R } .
By using Definition 3.2,
S ( t ) = S ( 0 ) + 1 Γ ( α ) 0 t ( t σ ) α 1 K 1 ( σ , S ) d σ , V ( t ) = V ( 0 ) + 1 Γ ( α ) 0 t ( t σ ) α 1 K 2 ( σ , V ) d σ , E ( t ) = E ( 0 ) + 1 Γ ( α ) 0 t ( t σ ) α 1 K 3 ( σ , E ) d σ , I ( t ) = I ( 0 ) + 1 Γ ( α ) 0 t ( t σ ) α 1 K 4 ( σ , I ) d σ , R ( t ) = R ( 0 ) + 1 Γ ( α ) 0 t ( t σ ) α 1 K 5 ( σ , R ) d σ .
Here,
K 1 ( σ , S ) = A δ 1 S I δ 2 S E μ S z S + α 0 R , K 2 ( σ , V ) = z S δ 3 V I δ 4 V E μ V + α 1 R , K 3 ( σ , E ) = δ 1 S I + δ 2 S E + δ 3 V I + δ 4 V E η E μ E , K 4 ( σ , I ) = η E γ I μ I μ 1 I , K 5 ( σ , R ) = γ I μ R α 0 R α 1 R .
Assume all population groups are positive bounded function that is ∃ some positive constants Π 1 , Π 2 , Π 3 , Π 4 , Π 5 , such that
S ( t ) Π 1 , V ( t ) Π 2 , E ( t ) Π 3 , I ( t ) Π 4 , R ( t ) Π 5 .
Theorem 5.1. 
If 0 M = m a x { Δ 1 , Δ 2 , Δ 3 , Δ 4 , Δ 5 } < 1 , then K i for i = 1 , 2 , 3 , 4 , 5 satisfies Lipchitz conditions.
Proof. 
Consider K 1 , for any S and S 1 ,
K 1 ( t , S ) K 1 ( t , S 1 ) = δ 1 I S δ 2 E S μ S z S + δ 1 I S 1 + δ 2 E S 1 + μ S 1 + z S 1 = δ 1 I ( S 1 S ) + δ 2 E ( S 1 S ) + μ ( S 1 S ) + z ( S 1 S ) δ 1 I ( t ) + δ 2 E ( t ) + μ + z S 1 S Δ 1 S S 1
where Δ 1 = δ 1 Π 4 + δ 2 Π 3 + μ + z . Thus K 1 satisfy Lipchitz condition. Similarly we can find Δ j for j = 2 , 3 , 4 , 5 . So that K i for i = 2 , 3 , 4 , 5 the lipchitz’s requirements are met under the condition, 0 M = m a x { Δ 1 , Δ 2 , Δ 3 , Δ 4 , Δ 5 } < 1 , the functions are contractions. Hence the proof.
Now, we rewrite (4) recursively as,
S n ( t ) = 1 Γ ( α ) 0 t ( t σ ) α 1 K 1 ( σ , S n 1 ) d σ , V n ( t ) = 1 Γ ( α ) 0 t ( t σ ) α 1 K 2 ( σ , V n 1 ) d σ , E n ( t ) = 1 Γ ( α ) 0 t ( t σ ) α 1 K 3 ( σ , E n 1 ) d σ , I n ( t ) = 1 Γ ( α ) 0 t ( t σ ) α 1 K 4 ( σ , I n 1 ) d σ , R n ( t ) = 1 Γ ( α ) 0 t ( t σ ) α 1 K 5 ( σ , R n 1 ) d σ .
The difference between two terms can be represented as
Λ n ( t ) = S n ( t ) S n 1 ( t ) = 1 Γ ( α ) 0 t ( t σ ) α 1 K 1 ( σ , S n 1 ) K 1 ( σ , S n 2 ) d σ , Ω n ( t ) = V n ( t ) V n 1 ( t ) = 1 Γ ( α ) 0 t ( t σ ) α 1 K 2 ( σ , V n 1 ) K 2 ( σ , V n 2 ) d σ , χ n ( t ) = E n ( t ) E n 1 ( t ) = 1 Γ ( α ) 0 t ( t σ ) α 1 K 3 ( σ , E n 1 ) K 3 ( σ , E n 2 ) d σ , ψ n ( t ) = I n ( t ) I n 1 ( t ) = 1 Γ ( α ) 0 t ( t σ ) α 1 K 4 ( σ , I n 1 ) K 4 ( σ , I n 2 ) d σ , ϑ n ( t ) = R n ( t ) R n 1 ( t ) = 1 Γ ( α ) 0 t ( t σ ) α 1 K 5 ( σ , R n 1 ) K 5 ( σ , R n 2 ) d σ .
where,
S n ( t ) = i = 0 n Λ i ( t ) , V n ( t ) = i = 0 n Ω i ( t ) , E n ( t ) = i = 0 n χ i ( t ) , I n ( t ) = i = 0 n ψ i ( t ) , R n ( t ) = i = 0 n ϑ i ( t ) .
Now consider,
Λ n ( t ) = S n ( t ) S n 1 ( t ) = 1 Γ ( α ) 0 t ( t σ ) α 1 K 1 ( σ , S n 1 ) K 1 ( σ , S n 2 ) d σ = Δ 1 Γ ( α ) 0 t S n 1 S n 2 d σ = Δ 1 Γ ( α ) 0 t Λ n 1 ( t ) d σ .
In the same fashion, we obtain
Ω n ( t ) = Δ 2 Γ ( α ) 0 t Ω n 1 ( t ) d σ , χ n ( t ) = Δ 3 Γ ( α ) 0 t χ n 1 ( t ) d σ , ψ n ( t ) = Δ 4 Γ ( α ) 0 t ψ n 1 ( t ) d σ , ϑ n ( t ) = Δ 5 Γ ( α ) 0 t ϑ n 1 ( t ) d σ .
Theorem 5.2.(i) The functions stated in (8) exists and are smooth, (ii) If t 0 > 1 Δ i Γ ( α ) 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 K i , where i = 1 , 2 , 3 , 4 , 5 satisfy Lipchitz’s conditions, then we will get the below relation.
Λ n ( t ) S ( 0 ) Δ 1 Γ ( α ) t n , Ω n ( t ) V ( 0 ) Δ 2 Γ ( α ) t n , χ n ( t ) E ( 0 ) Δ 3 Γ ( α ) t n , ψ n ( t ) I ( 0 ) Δ 4 Γ ( α ) t n , ϑ n ( t ) R ( 0 ) Δ 5 Γ ( α ) t n .
Therefore equation (9) shows the existence and smoothness of all populations which was defined in (8).
(ii) Here we show, S n , V n , E n , I n and R n converge to solutions of (2). Now let, u n * , w n * , x n * , y n * , z n * , are remainder terms after n iterations, ∋
S S ( 0 ) = S n u n * , V V ( 0 ) = V n w n * , E E ( 0 ) = E n x n * , I I ( 0 ) = I n y n * , R R ( 0 ) = R n z n * .
By using triangle inequality along with K 1 , we get
u n * ( t ) = 1 Γ ( α ) 0 t ( t σ ) α 1 K 1 ( σ , S ) K 1 ( σ , S n 1 ) d σ Δ 1 Γ ( α ) S S n 1 t
We obtain the below equation by iteratively using the previous process,
u n * ( t ) Δ 1 Γ ( α ) t n + 1 Π 1
Then at t 0 , we have
u n * ( t ) Δ 1 Γ ( α ) t 0 n + 1 Π 1 .
Taking limit as n ,
lim n u n * ( t ) lim n Δ 1 Γ ( α ) t 0 n + 1 Π 1 .
Using hypothesis, Δ 1 Γ ( α ) t 0 < 1 . Equation (11) becomes,
lim n u n * ( t ) = 0
As n , we get
w n * ( t ) 0 , x n * ( t ) 0 , y n * ( t ) 0 , z n * ( t ) 0 .
Thus, the system has at least one solution.
Theorem 5.3. 
A unique solution exists for (2) if 1 Δ i Γ ( α ) t > 0 for i = 1 to 5 .
Proof. 
Let S 1 ( t ) , V 1 ( t ) , E 1 ( t ) , I 1 ( t ) , R 1 ( t ) is another set of solution of (2) thus,
S ( t ) S 1 ( t ) = 1 Γ ( α ) 0 t ( t σ ) α 1 K 1 ( s , S ) K 1 ( s , S 1 ) d s Δ 1 Γ ( α ) t S ( t ) S 1 ( t )
By rearranging,
S ( t ) S 1 ( t ) Δ 1 Γ ( α ) t S ( t ) S 1 ( t ) 0
S ( t ) S 1 ( t ) 1 Δ 1 Γ ( α ) t 0
using hypothesis 1 Δ 1 Γ ( α ) t is strictly greater than zero, then (12) takes the structure S ( t ) S 1 ( t ) = 0 .
It means that S ( t ) = S 1 ( t ) .
By repeating the similar procedure to each solution for i = 2 , 3 , 4 , 5 , we get,
V ( t ) = V 1 ( t ) , E ( t ) = E 1 ( t ) , I ( t ) = I 1 ( t ) , R ( t ) = R 1 ( t ) .
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 = ( A θ 1 , z A θ 1 μ 0, 0, 0). Also for endemic equilibrium we obtain S * = η A θ 4 + α 0 γ η I * θ 4 ( η θ 1 + δ 1 η I * + δ 2 θ 3 η I * ) , V * = η z θ 4 ( η A θ 4 + α 0 γ η I * θ 4 ( η θ 1 + δ 1 η I * + δ 2 θ 3 η I * ) + η α 1 γ I * ) θ 4 ( δ 3 η I * + δ 4 θ 3 I * + η μ ) , E * = θ 3 I * η and R * = γ I * θ 4 . Here I * will be in the form of g ( I * ) = 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 = η S δ 1 + η V δ 3 + S δ 2 θ 3 + θ 3 V δ 4 θ 2 θ 3 . By substituting disease free equilibrium points we get,
R 0 = η ( A θ 1 ) δ 1 + η ( z A θ 1 μ ) δ 3 + ( A θ 1 ) δ 2 θ 3 + θ 3 ( z A θ 1 μ ) δ 4 θ 2 θ 3 .
where,
θ 1 = μ + z ,
θ 2 = η + μ ,
θ 3 = γ + μ + μ 1 .

6.3. Local stability

The Jacobian matrix for (2) is presented by,
J = δ 1 I δ 2 E θ 1 0 δ 2 S δ 1 S α 0 z δ 3 I δ 4 E μ δ 4 V δ 3 V α 1 δ 1 I + δ 2 E δ 3 I + δ 4 E δ 2 S + δ 4 V θ 2 δ 3 V + δ 1 S 0 0 0 η θ 3 0 0 0 0 γ θ 4 .
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( λ i )|> a π 2 .
The Jacobi matrix at E 0 is obtained as follows,
J 0 = θ 1 0 δ 2 A θ 1 δ 1 A θ 1 α 0 z μ δ 4 z A θ 1 μ δ 3 z A θ 1 μ α 1 0 0 δ 2 A θ 1 + δ 4 z A θ 1 μ θ 2 δ 3 z A θ 1 μ + δ 1 A θ 1 0 0 0 η θ 3 0 0 0 0 γ θ 4
Clearly, from the above matrix we can say, λ 1 = θ 1 , λ 2 = μ , λ 5 = θ 4 are the three eigen values. Here λ 3 and λ 4 can be computed from the below matrix.
δ 2 A θ 1 + δ 4 z A θ 1 μ θ 2 δ 3 z A θ 1 μ + δ 1 A θ 1 η θ 3
They are λ 3 = y 1 + y 2 + y 3 and λ 4 = y 1 y 2 + y 3 , where,
y 1 = A δ 2 μ θ 1 θ 2 μ θ 1 θ 3 μ + A δ 4 z , y 2 = A 2 δ 2 2 μ 2 + θ 1 2 θ 2 2 μ 2 + 2 θ 1 2 θ 2 θ 3 μ 2 + θ 1 2 θ 3 2 μ 2 + 2 A 2 δ 2 δ 4 μ z + A 2 δ 4 2 z 2 2 A δ 2 θ 1 θ 2 μ 2 2 A δ 2 θ 1 θ 3 μ 2 2 A δ 4 θ 1 θ 2 μ z 2 A δ 4 θ 1 θ 3 μ z , y 3 = 4 ( A δ 2 θ 1 θ 3 μ 2 + δ 1 A θ 1 η μ 2 + A δ 4 θ 1 θ 3 μ z + A δ 3 θ 1 η μ z θ 1 2 θ 2 θ 3 μ 2 ) .
y 1 < 0 when θ 1 θ 2 μ + θ 1 θ 3 μ > A δ 2 μ + A δ 4 z , y 2 > 0 when A 2 δ 2 2 μ 2 + θ 1 2 θ 2 2 μ 2 + 2 θ 1 2 θ 2 θ 3 μ 2 + θ 1 2 θ 3 2 μ 2 + 2 A 2 δ 2 δ 4 μ z + A 2 δ 4 2 z 2 > 2 A δ 2 θ 1 θ 2 μ 2 + 2 A δ 2 θ 1 θ 3 μ 2 + 2 A δ 4 θ 1 θ 2 μ z + 2 A δ 4 θ 1 θ 3 μ z , y 3 > 0 when 4 ( A δ 2 θ 1 θ 3 μ 2 + δ 1 A θ 1 η μ 2 + A δ 4 θ 1 θ 3 μ z + A δ 3 θ 1 η μ z ) > 4 θ 1 2 θ 2 θ 3 μ 2 .
Therefore λ 3 and λ 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 = ( A θ 1 , z A θ 1 μ 0, 0, 0) is stable whenever BRN is less than one.
Table 2. Necessary components for the stability of the equilibrium points.
Table 2. Necessary components for the stability of the equilibrium points.
Eigen value Sign Conditions Stability
λ 3 y 1 < 0 , y 2 > 0 and y 3 > 0 and hence y 1 > y 2 + y 3 Stable
λ 4 y 1 < 0 , y 2 > 0 and y 3 > 0 and hence y 1 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( λ i )|> a π 2 .
For (2), the Jacobi matrix at E 1 is obtained as,
J 1 = δ 1 I * δ 2 E * θ 1 0 δ 2 S * δ 1 S * α 0 z δ 3 I * δ 4 E * μ δ 4 V * δ 3 V * α 1 δ 1 I * + δ 2 E * δ 3 I * + δ 4 E * δ 2 S * + δ 4 V * θ 2 δ 3 V * + δ 1 S * 0 0 0 η θ 3 0 0 0 0 γ θ 4
Let us consider,
b 11 = δ 1 I * δ 2 E * θ 1 ,
b 13 = δ 2 S * ,
b 14 = δ 1 S * ,
b 15 = α 0 ,
b 21 = z ,
b 22 = δ 3 I * δ 4 E * μ ,
b 23 = δ 4 V * ,
b 24 = δ 3 V * ,
b 25 = α 1 ,
b 31 = δ 1 I * + δ 2 E * ,
b 32 = δ 3 I * + δ 4 E * ,
b 33 = δ 2 S * + δ 4 V * θ 2 ,
b 34 = δ 3 V * + δ 1 S * ,
b 43 = η ,
b 44 = θ 3 ,
b 54 = γ ,
b 55 = θ 4 .
The characteristic equation of the J 1 is,
λ 5 + x 1 λ 4 + x 2 λ 3 + x 3 λ 2 + x 4 λ + x 5 = 0 ,
where,
x 1 = b 11 b 22 b 33 b 44 b 55 , 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 b 34 b 43 b 23 b 32 b 13 b 31 , 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 + 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 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 b 13 b 21 b 32 , 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 + 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 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 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 , 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 + 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 b 15 b 21 b 32 b 43 b 54 .
Here we consider,
x 1 > 0 , if b 11 + b 22 + b 33 + b 44 + b 55 < 0 , x 2 > 0 , if b 11 b 33 + b 11 b 22 > b 34 b 43 + b 23 b 32 + b 13 b 31 , x 3 > 0 , 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 + 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 + b 11 b 22 b 33 + b 13 b 21 b 32 , x 4 > 0 , 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 + 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 + b 1 b 22 b 34 b 43 + b 14 b 21 b 32 b 43 , x 5 > 0 , 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 + b 14 b 22 b 31 b 43 b 55 + b 15 b 21 b 32 b 43 b 54 .
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
D t α [ S ( t ) ] = τ 1 ( t , S ( t ) ) , D t α [ V ( t ) ] = τ 2 ( t , V ( t ) ) , D t α [ E ( t ) ] = τ 3 ( t , E ( t ) ) , D t α [ I ( t ) ] = τ 4 ( t , I ( t ) ) , D t α [ R ( t ) ] = τ 5 ( t , R ( t ) ) .
where,
τ 1 ( t , S ( t ) ) = A δ 1 S I δ 2 S E μ S z S + α 0 R , τ 2 ( t , V ( t ) ) = z S δ 3 V I δ 4 V E μ V + α 1 R , τ 3 ( t , E ( t ) ) = δ 1 S I + δ 2 S E + δ 3 V I + δ 4 V E η E μ E , τ 4 ( t , I ( t ) ) = η E γ I μ I μ 1 I , τ 5 ( t , R ( t ) ) = γ I μ R α 0 R α 1 R .
From first equation of (2),
D t α [ S ( t ) ] = τ 1 ( t , S ( t ) ) , S ( 0 ) = S 0 , t > 0 .
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 ( t ) 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 = r h , r = 0 , 1 , 2 , 3 , . . . k , we partition the interval [ 0 , d ] into k subintervals [ t r , t r + 1 ] with equal width h = d k . Assume S ( t ) , D t α [ S ( t ) ] and D t 2 α [ S ( t ) ] are continuous on [ 0 , d ] . By using the generalised Taylor formula, extend S ( t ) to t = t 0 = 0 . There is a value c 1 for each value t, therefore
S ( t ) = S ( t 0 ) + D t α [ S ( t ) ] t 0 t α Γ ( α + 1 ) + D t 2 α [ S ( t ) ] c 1 t 2 α Γ ( 2 α + 1 )
Substitute D t α [ S ( t ) ] t 0 = τ 1 ( t 0 , S ( t 0 ) ) and h = t 1 in (16)
S ( t 1 ) = S ( t 0 ) + τ 1 ( t 0 , S ( t 0 ) ) h α Γ ( α + 1 ) + D t 2 α [ S ( t ) ] c 1 h 2 α Γ ( 2 α + 1 ) ,
If h is small, thus
S ( t 1 ) = S ( t 0 ) + τ 1 ( t 0 , S ( t 0 ) ) h α Γ ( α + 1 ) ,
A general formula, t r + 1 = t r + h can be written as
S ( t r + 1 ) = S ( t r ) + τ 1 ( t r , S ( t r ) ) h α Γ ( α + 1 ) .
Using the fractional integral on (15), we obtain
[ S ( t ) ] = S ( 0 ) + I α [ τ 1 ( t , S ( t ) ) ]
To obtain ( t 1 , S ( t 1 ) ) , we replace t = t 1 into (19), then
S ( t 1 ) = S ( 0 ) + ( I α [ τ 1 ( t , S ( t ) ) ] ) ( t 1 ) .
With the help modified trapezoidal rule, we approximate ( I α [ τ 1 ( t , S ( t ) ) ] ) ( t 1 ) with h = t 1 t 0 then (20) becomes
S ( t 1 ) = S ( 0 ) + α h α [ τ 1 ( t 0 , S ( t 0 ) ) ] Γ ( α + 2 ) + h α [ τ 1 ( t 1 , S ( t 1 ) ) ] Γ ( α + 2 )
The expression S ( t 1 ) appears in the formula on the RHS of (21). As a result, we estimate S ( t 1 ) . For this purpose, the fractional Euler approach can be used. Substitute (17) for (21).
S ( t 1 ) = S ( 0 ) + α h α [ τ 1 ( t 0 , S ( t 0 ) ) ] Γ ( α + 2 ) + h α [ τ 1 ( t 1 , S ( t 0 ) ) ] + h α Γ ( α + 1 ) τ 1 ( t 0 , S ( t 0 ) ) Γ ( α + 2 )
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,
S ( t r ) = S ( 0 ) + h α Γ ( α + 2 ) [ ( r 1 ) α + 1 ( r α 1 ) r α ] τ 1 ( t 0 , S ( t 0 ) ) + h α Γ ( α + 2 ) i = 1 r 1 [ ( r i + 1 ) α + 1 2 ( r 1 ) α + 1 + ( r i 1 ) α + 1 ] τ 1 ( t i