# On Fractional Order Dengue Epidemic Model

Authors: Hamed Al-Sulami; Moustafa El-Shahed; Juan J. Nieto; Wafa Shammakh
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper deals with the fractional order dengue epidemic  model. The stability of disease-free and positive fixed points is studied. Adams-Bashforth-Moulton algorithm  has been used to solve and simulate the system of differential equations.


---

## Body

## 1. Introduction

Dengue is a major public health problem in tropical and subtropical countries. It is a vector-borne disease transmitted byAedes aegypti andAedes albopictus mosquitoes. Four different serotypes can cause dengue fever. A human infected by one serotype, when recovers, gains total immunity to that serotype and only partial and transient immunity with respect to the other three.Dengue can vary from mild to severe. The more severe forms of dengue include shock syndrome and dengue hemorrhagic fever (DHF). Patients who develop these more serious forms of dengue fever usually need to be hospitalized. The full life cycle of dengue fever virus involves the role of the mosquito as a transmitter (or vector) and humans as the main victim and source of infection. Preventing or reducing dengue virus transmission depends entirely on the control of mosquito vectors or interruption of human vector contact [1, 2].In this paper we study the fractional order dengue epidemic model. The stability of equilibrium points is studied. Numerical solutions of this model are given. We like to argue that fractional order equations are more suitable than integer order ones in modeling biological, economic, and social systems (generally complex adaptive systems) where memory effects are important. Adams-Bashforth-Moulton algorithm has been used to solve and simulate the system of differential equations.


## 2. Model Derivation

Esteva and Vargas [3] developed a dengue fever transmission model by assuming that, once a person recovers from the disease, he or she will not be reinfected by the disease. The model also assumes that the host population N
h is constant, that is, the death rate and the birth rate equal μ
H. The host-vector model for the dengue transmission of Esteva and Vargas [3] is as follows:
(1)
d
S
h
d
t
=
A
-
β
h
b
N
h
S
h
I
v
-
μ
h
S
h
,
d
I
h
d
t
=
β
h
b
N
h
S
h
I
v
-
(
μ
h
+
γ
)
I
h
,
d
R
h
d
t
=
γ
I
h
-
μ
h
R
h
,
d
S
v
d
t
=
B
-
β
v
b
N
h
S
v
I
h
-
μ
v
S
v
,
d
I
v
d
t
=
β
v
b
N
h
S
v
I
h
-
μ
v
I
v
,
whereA is the recruitment rate of the host population,
B is the recruitment rate of the vector population,
S
h is the number of susceptible in the host population,
I
h is the number of infective in the host population,
R
h is the number of immunes in the host population,
N
v is the vector population,
S
v is the number of susceptible in the vector population,
I
v is the number of infective in the vector population,
μ
v is the death rate in the vector population,
β
h is the transmission probability from vector to host,
β
v is the transmission probability from host to vector,
γ is the recovery rate in the host population,
b is the biting rate of the vector.The notion of fractional calculus was anticipated by Leibniz, one of the founders of standard calculus, in a letter written in 1695. Recently great considerations have been made to the models of FDEs in different aria of researches. The most essential property of these models is their nonlocal property which does not exist in the integer order differential operators. We mean by this property that the next state of a model depends not only upon its current state but also upon all of its historical states. There are many definitions of fractional derivatives [4, 5]. Perhaps the best-known is the Riemann-Liouvile definition. The Riemann-Liouville derivative of order α is defined as
(2)
RL
D
0
+
α
f
(
t
)
=
1
Γ
(
n
-
α
)
(
d
d
t
)
n
∫
0
t
f
(
s
)
(
t
-
s
)
α
-
n
+
1
d
s
,
000000000000000000000000000000
n
=
[
α
]
+
1
,
where Γ() is the gamma function and n is an integer. An alternative definition was introduced by Caputo as follows, which is a sort of regularization of the Riemann-Liouville derivative:
(3)
D
t
α
f
(
t
)
=
1
Γ
(
n
-
α
)
∫
0
t
f
(
n
)
(
s
)
(
t
-
s
)
α
-
n
+
1
d
s
.Pooseh et al. [6] introduced the notion of fractional derivative in the sense of Riemann-Liouville to reformulate the dynamics of the classical model (1) in terms of fractional derivatives. They applied a recent approximate technique to obtain numerical solutions to the fractional model. The system in this paper will be in the sense of Caputo fractional derivative by the following set of fractional order differential equations:
(4)
D
t
α
S
h
=
A
-
β
h
b
N
h
S
h
I
v
-
μ
h
S
h
,
D
t
α
I
h
=
β
h
b
N
h
S
h
I
v
-
(
μ
h
+
γ
)
I
h
,
D
t
α
R
h
=
γ
I
h
-
μ
h
R
h
,
D
t
α
S
v
=
B
-
β
v
b
N
h
S
v
I
h
-
μ
v
S
v
,
D
t
α
I
v
=
β
v
b
N
h
S
v
I
h
-
μ
v
I
v
.
Because model (4) monitors the dynamics of human populations, all the parameters are assumed to be nonnegative. Furthermore, it can be shown that all state variables of the model are nonnegative for all time t
≥
0 (see, for instance, [7–9]).Lemma 1.
The closed setΩ
=
{
(
S
h
,
I
h
,
R
h
,
S
v
,
I
v
)
∈
R
+
5
:
S
h
+
I
h
+
R
h
=
A
/
μ
h
,
S
v
+
I
v
=
B
/
μ
v
} is positively invariant with respect to model (4).Proof.
The fractional derivative of the total population, obtained by adding all the equations of model (4), is given by
(5)
D
t
α
N
h
(
t
)
=
A
-
μ
h
N
h
(
t
)
.
The solution to (5) is given by N
h
(
t
)
=
N
h
(
0
)
E
α
,
1
(
-
μ
h
t
α
)
+
A
t
α
E
α
,
α
+
1
(
-
μ
h
t
α
), where E
α
,
β is the Mittag-Leffler function. Considering the fact that the Mittag-Leffler function has an asymptotic behavior [4, 10],
(6)
E
α
,
β
(
z
)
~
-
∑
k
=
1
ω
z
-
k
Γ
(
β
-
α
k
)
+
O
(
|
z
|
-
1
-
ω
)
,
000000
(
|
z
|
⟶
∞
,
α
π
2
<
|
arg
⁡
(
z
)
|
≤
π
)
.
One can observe thatN
h
(
t
)
→
A
/
μ
h as t
→
∞. The proof of vector population case is completely similar to that of host population and is therefore omitted. One can observe that N
v
(
t
)
→
B
/
μ
v. Therefore, all solutions of the model with initial conditions in Ω remain in Ω for all t > 0. Thus, region Ω is positively invariant with respect to model (4).In the following, we will study the dynamics of system (4).


## 3. Equilibrium Points and Stability

To evaluate the equilibrium points let(7)
D
t
α
S
h
=
0
,
D
t
α
I
h
=
0
,
D
t
α
R
h
=
0
,
D
t
α
S
v
=
0
,
D
t
α
I
v
=
0
.ThenE
0
=
(
A
/
μ
h
,
0,0
,
(
B
/
μ
v
)
0
). By (4), a positive equilibrium E
1
=
(
S
h
1
,
I
h
1
,
R
h
1
,
S
v
1
,
I
v
1
) satisfies
(8)
S
h
1
=
N
h
μ
v
(
A
b
β
v
+
N
h
(
γ
+
μ
h
)
μ
v
)
b
β
v
(
b
B
β
h
+
N
h
μ
h
μ
v
)
,
I
h
1
=
A
b
2
B
β
h
β
v
-
N
h
2
μ
h
(
γ
+
μ
h
)
μ
v
2
b
β
v
(
γ
+
μ
h
)
(
b
B
β
h
+
N
h
μ
h
μ
v
)
,
R
h
1
=
A
b
2
B
γ
β
h
β
v
-
γ
N
h
2
μ
h
(
γ
+
μ
h
)
μ
v
2
b
β
v
μ
h
(
γ
+
μ
h
)
(
b
B
β
h
+
N
h
μ
h
μ
v
)
,
S
v
1
=
N
h
(
γ
+
μ
h
)
(
b
B
β
h
+
N
h
μ
h
μ
v
)
b
β
h
(
A
b
β
v
+
N
h
(
γ
+
μ
h
)
μ
v
)
,
I
v
1
=
A
b
2
B
β
h
β
v
-
N
h
2
μ
h
(
γ
+
μ
h
)
μ
v
2
b
β
h
μ
v
(
A
b
β
v
+
N
h
(
γ
+
μ
h
)
μ
v
)
.
The Jacobian matrix J
(
E
0
) for the system given in (4) evaluated at the disease-free equilibrium is as follows:
(9)
J
(
E
0
)
=
(
-
μ
h
0
0
0
-
A
b
β
h
N
h
μ
h
0
-
γ
-
μ
h
0
0
A
b
β
h
N
h
μ
h
0
γ
-
μ
h
0
0
0
-
b
B
β
v
N
h
μ
v
0
-
μ
v
0
0
b
B
β
v
N
h
μ
v
0
0
-
μ
v
)
.Theorem 2.
The disease-free equilibriumE
0 is locally asymptotically stable if R
0
<
1 and is unstable if R
0
>
1.Proof.
The disease-free equilibrium is locally asymptotically stable if all the eigenvalues,λ
i, i
=
1,2
,
3,4
,
5 of the Jacobian matrix J
(
E
0
) satisfy the following condition [11–14]:
(10)
|
arg
(
λ
i
)
|
>
α
π
2
.
The eigenvalues of the Jacobian matrixJ
(
E
0
) are λ
1
=
-
μ
h, λ
2
=
-
μ
h, and λ
3
=
-
μ
v; the other two roots are determined by the quadratic equation
(11)
λ
2
+
λ
(
γ
+
μ
h
+
μ
v
)
+
μ
v
(
γ
+
μ
h
)
(
1
-
R
0
)
=
0
,
where R
0
=
(
A
B
b
2
β
h
β
v
)
/
(
μ
h
μ
v
2
N
h
2
(
γ
+
μ
h
)
). Hence E
0 is locally asymptotically stable if R
0
<
1 and is unstable if R
0
>
1.
The quantityR
0
*
=
R
0 is called the basic reproductive number of the disease, since it represents the average number of secondary cases that one case can produce if introduced into a susceptible population.We now discuss the asymptotic stability of the endemic (positive) equilibrium of the system given by (4). The Jacobian matrix J
(
E
1
) evaluated at the endemic equilibrium is given as
(12)
J
(
E
1
)
(
-
b
β
h
I
v
1
N
h
-
μ
h
0
0
0
-
b
S
h
1
β
h
N
h
b
β
h
I
v
1
N
h
-
γ
-
μ
h
0
0
b
S
h
1
β
h
N
h
0
γ
-
μ
h
0
0
0
-
b
S
h
1
β
v
N
h
0
-
b
β
v
I
h
1
N
h
-
μ
v
0
0
b
S
v
1
β
v
N
h
0
b
β
v
I
h
1
N
h
-
μ
v
)
.The characteristic equation ofJ
(
E
1
) is
(13)
(
λ
+
μ
h
)
(
λ
+
μ
v
)
(
λ
3
+
a
1
λ
2
+
a
2
λ
+
a
3
)
=
0
,
where
(14)
a
1
=
γ
+
μ
v
+
M
2
(
K
θ
+
M
1
μ
v
)
K
θ
M
1
N
h
μ
v
,
a
2
=
(
M
2
2
+
K
θ
M
1
N
h
2
μ
h
μ
v
(
θ
+
μ
v
)
0
+
M
2
N
h
(
K
θ
2
+
(
K
θ
+
M
1
(
θ
+
μ
h
)
)
μ
v
)
)
×
(
K
θ
M
1
N
h
2
μ
v
)
-
1
,
a
3
=
M
2
(
M
2
+
N
h
(
K
θ
+
M
1
μ
h
)
μ
v
)
K
M
1
N
h
2
μ
v
,
M
1
=
A
b
2
B
β
h
β
v
-
N
h
2
μ
h
(
γ
+
μ
h
)
μ
v
2
,
M
2
=
A
b
β
v
+
N
h
(
γ
+
μ
h
)
μ
v
,
K
=
b
B
β
h
+
N
h
μ
h
μ
v
,
θ
=
γ
+
μ
h
.Ifp
(
x
)
=
x
3
+
a
1
x
2
+
a
2
x
+
a
3. Let D
(
p
) denote the discriminant of a polynomial p
(
x
); then
(15)
D
(
p
)
=
-
|
1
a
1
a
2
a
3
0
0
1
a
1
a
2
a
3
3
2
a
1
a
2
0
0
0
3
a
1
a
2
0
0
0
3
2
a
1
a
2
|
=
18
a
1
a
2
a
3
+
(
a
1
a
2
)
2
-
4
a
3
a
1
3
-
4
a
2
3
-
27
a
3
2
.Following [14–18], we have Proposition 3.Proposition 3.
One assumes thatE
1 exists in R
+
3.(i)
If the discriminant ofp
(
x
), D
(
p
), is positive and Routh-Hurwitz are satisfied, that is, D
(
p
)
>
0, a
1
>
0, a
3
>
0, and a
1
a
2
>
a
3, then E
1 is locally asymptotically stable.
(ii)
IfD
(
p
)
<
0, a
1
>
0, a
2
>
0, a
1
a
2
=
a
3, and α
∈
[
0,1
), then E
1 is locally asymptotically stable.
(iii)
IfD
(
p
)
<
0, a
1
<
0, a
2
<
0, and α
>
2
/
3, then E
1 is unstable.
(iv)
The necessary condition for the equilibrium pointE
1, to be locally asymptotically stable, is a
3
>
0.


## 4. Numerical Methods and Simulations

Since most of the fractional order differential equations do not have exact analytic solutions, so approximation and numerical techniques must be used. Several analytical and numerical methods have been proposed to solve the fractional order differential equations. For numerical solutions of the system (4) one can use the generalized Adams-Bashforth-Moulton method. To give the approximate solution by means of this algorithm, consider the following nonlinear fractional differential equation [19]:
(16)
000000000000.000000
D
t
α
y
(
t
)
=
f
(
t
,
y
(
t
)
)
,
0
≤
t
≤
T
,
y
(
k
)
(
0
)
=
y
0
k
,
k
=
0,1
,
2
,
…
,
m
-
1
,
where
m
=
[
α
]
.
This equation is equivalent to Volterra integral equation:
(17)
y
(
t
)
=
∑
k
=
0
m
-
1
y
0
(
k
)
t
k
k
!
+
1
Γ
(
α
)
∫
0
t
(
t
-
s
)
α
-
1
f
(
s
,
y
(
s
)
)
d
s
.Diethelm et al. used the predictor-correctors scheme [15, 16, 20] based on the Adams-Bashforth-Moulton algorithm to integrate (17). By applying this scheme to the fractional order dengue epidemic model and setting h
=
T
/
N, t
n
=
n
h, and n
=
0,1
,
2
,
,
…
,
N
∈
Z
+, (17) can be discretized as follows [19]:
(18)
S
n
+
1
=
S
0
+
h
α
Γ
(
α
+
2
)
(
A
-
β
h
b
h
N
h
S
n
+
1
p
Y
n
+
1
p
-
μ
h
S
n
+
1
p
)
+
h
α
Γ
(
α
+
2
)
∑
j
=
0
n
a
j
,
n
+
1
(
A
-
β
h
b
h
N
h
S
j
Y
j
-
μ
h
S
j
)
,
I
n
+
1
=
I
0
+
h
α
Γ
(
α
+
2
)
(
β
h
b
h
N
h
S
n
+
1
p
Y
n
+
1
p
-
(
μ
h
+
γ
)
I
n
+
1
p
)
+
h
α
Γ
(
α
+
2
)
∑
j
=
0
n
a
j
,
n
+
1
(
β
h
b
h
N
h
S
j
Y
j
-
(
μ
h
+
γ
)
I
j
)
,
R
n
+
1
=
R
0
+
h
α
Γ
(
α
+
2
)
(
γ
I
n
+
1
p
-
μ
h
R
n
+
1
p
)
+
h
α
Γ
(
α
+
2
)
∑
j
=
0
n
a
j
,
n
+
1
(
γ
I
j
-
μ
h
R
j
)
,
X
n
+
1
=
X
0
+
h
α
Γ
(
α
+
2
)
(
B
-
β
v
b
h
N
h
X
n
+
1
p
I
n
+
1
p
-
μ
h
X
n
+
1
p
)
+
h
α
Γ
(
α
+
2
)
∑
j
=
0
n
a
j
,
n
+
1
(
B
-
β
v
b
h
N
h
X
j
I
j
-
μ
v
X
j
)
,
Y
n
+
1
=
Y
0
+
h
α
Γ
(
α
+
2
)
(
β
v
b
h
N
h
X
n
+
1
p
I
n
+
1
p
-
μ
h
Y
n
+
1
p
)
+
h
α
Γ
(
α
+
2
)
∑
j
=
0
n
a
j
,
n
+
1
(
β
v
b
h
N
h
X
j
I
j
-
μ
v
Y
j
)
,
where
(19)
S
n
+
1
p
=
S
0
+
1
Γ
(
α
)
∑
j
=
0
n
b
j
,
n
+
1
(
A
-
β
h
b
h
N
h
S
j
Y
j
-
μ
h
S
j
)
,
I
n
+
1
p
=
I
0
+
1
Γ
(
α
)
∑
j
=
0
n
b
j
,
n
+
1
(
β
h
b
h
N
h
S
j
Y
j
-
(
μ
h
+
γ
)
I
j
)
,
R
n
+
1
p
=
R
0
+
1
Γ
(
α
)
∑
j
=
0
n
b
j
,
n
+
1
(
γ
I
j
-
μ
h
R
j
)
,
X
n
+
1
p
=
X
0
+
1
Γ
(
α
)
∑
j
=
0
n
b
j
,
n
+
1
(
B
-
β
v
b
h
N
h
X
j
I
j
-
μ
v
X
j
)
,
Y
n
+
1
p
=
Y
0
+
1
Γ
(
α
)
∑
j
=
0
n
b
j
,
n
+
1
(
β
v
b
h
N
h
X
j
I
j
-
μ
v
Y
j
)
,
a
j
,
n
+
1
=
{
n
α
+
1
-
(
n
-
α
)
(
n
+
1
)
,
j
=
0
,
(
n
-
j
+
2
)
α
+
1
+
(
n
-
j
)
α
+
1
-
2
(
n
-
j
+
1
)
α
+
1
,
0
.
1
≤
j
≤
n
,
1
,
j
=
n
+
1
,
b
j
,
n
+
1
=
h
α
α
(
(
n
-
j
+
1
)
α
-
(
n
-
j
)
α
)
,
0
≤
j
≤
n
.


## 5. Discussion

In this paper, we have considered a fractional calculus model for dengue disease. Following [21], Figure 1 shows that S
h drops significantly in a relatively small period of time. Both I
h and I
v increase significantly during the period of 30 days and then eventually oscillate around the endemic state (0.09529, 0.0.00029, and 0.00058). This seems unrealistic in the nature. With constant population of mosquitoes, this fluctuation (in a short period of time) cannot be shown to happen in the nature [21]. As mentioned by [6], Figures 2 and 3 show that even a simple fractional model may give surprisingly good results. However, the transformation of a classical model into a fractional one makes it very sensitive to the order of differentiation α: a small change in α may result in a big change in the final result. From the numerical results in Figures 2 and 3, it is clear that the approximate solutions depend continuously on the fractional derivative α.Figure 1
S
h
(
t
), I
h
(
t
), and I
v
(
t
) for α
=
1 and  μ
h
=
0.0000457; μ
v
=
0.25; b
=
0.5; β
h
=
0.75; β
v
=
1; γ
=
0.1428; N
h
=
10000; B
=
5000.Figure 2
S
h
(
t
), I
h
(
t
),and I
v
(
t
) for α
=
0.9 and μ
h
=
0.0000457; μ
v
=
0.25; b
=
0.5; β
h
=
0.75; β
v
=
1; γ
=
0.1428; N
h
=
10000; B
=
5000.Figure 3
S
h
(
t
), I
h
(
t
), and I
v
(
t
) for α
=
0.8 and μ
h
=
0.0000457; μ
v
=
0.25; b
=
0.5; β
h
=
0.75; β
v
=
1; γ
=
0.1428; N
h
=
10000; B
=
5000.The approximate solutionsS
h
(
t
), I
h
(
t
), and I
v
(
t
) are displayed in Figures 4 and 5 with different values of α. In each figure three different values of α are considered. When α
=
1, system (4) is the classical integer-order system (1). In Figure 4, the variation of S
h
(
t
) versus time t is shown for different values of α
=
1
,
0.9
,
0
.
8 by fixing other parameters. It is revealed that S
h does not drop significantly in a relatively small period of time for small values. Figure 5 depicts I
h
(
t
) versus time t. As mentioned by [22, 23], one should note that although the equilibrium points are the same for both integer order and fractional order models the solution of the fractional order model tends to the fixed point over a longer period of time. One also needs to mention that when dealing with real life problems, the order of the system can be determined by using the collected data.Figure 4
S
h
(
t
) for α
=
1
,
0.9
,
0.8 and μ
h
=
0.0000457; μ
v
=
0.25; b
=
0.5; β
h
=
0.75; β
v
=
1; γ
=
0.1428; N
h
=
10000; B
=
5000.Figure 5
I
h
(
t
) for α
=
1
,
0.9
,
0.8, μ
h
=
0.0000457; μ
v
=
0.25; b
=
0.5; β
h
=
0.75; β
v
=
1; γ
=
0.1428; N
h
=
10000; B
=
5000.



---

*Source: 456537-2014-08-18.xml*