# Epidemic Model of Leptospirosis Containing Fractional Order

Authors: Muhammad Altaf Khan; S. F. Saddiq; Saeed Islam; Ilyas Khan; Dennis Ling Chuan Ching
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We study an epidemic model of leptospirosis in fractional order numerically. The multistep generalized differential transform method is applied to find the accurate approximate solution of the epidemic model of leptospirosis disease in fractional order. A unique positive solution for the epidemic model in fractional order is presented. For the integer case derivative, the approximate solution of MGDTM is compared with the Runge-Kutta order four scheme. The numerical results are presented for the justification purpose.


---

## Body

## 1. Introduction

Leptospirosis disease is an important infectious disease. This kind of infection occurs in urban areas of industrialized and developed countries and also in the rural areas. The people of the city who walk in dirty water are mostly infected. Workers planting rice, sewer cleaners, cleaning canals workers, and agriculture labor get the disease easily. The disease flourishes due to delay in diagnosis and unavailability of clinical infrastructure. The cause of the disease is bacteria. It is potentially fatal infection of brain, kidney, liver, heart, and lung. The people who can get infection are those who have contact with infected animals, soil, or water in which the bacteria is present. The outdoor people, who work with animals, face the risk of leptospirosis infection, similarly workers in farms, sewer, mine, slaughter houses, dairy farmers, and animal caretakers and those who work with fishes and military personnel. Those people who work outdoors like swimming, rafting, and kayaking also face the risk of infection [1].The mathematical modeling dealing with the complex biological structures is a great challenge for the researchers. For simple biological models, the study of integer order differential equations is enough for the description of their dynamics while for those systems contains complexities are characterized due to its variability structured properties, like nonlinearity of the systems, multiscale behavior, and the mathematical relation between the parameters [2]. To study such nonlinear biological models, with complex mathematical structures, the fractional derivatives provide the tool to handle the dynamical behavior of such complex systems. The fundamental property related to these models is the nonexistence of the differential operator of integer order. These characteristics show information about the present but also show its historical states.The fractional derivative and fractional integrals were first used by Magin [3] in his model which is related to the relation between stress-strain and biomaterials. A new fractional analytical approach via a modified Riemann Liouville derivative has been used by [4]. For solution of the 2D region with obstacles, the author [5] used the fractional differential equations. For the giving up smoking model, the fractional differential equation was used by [6]. The author used the differential transforms method and Pade approximants for a fractional population growth model [7]. For more information about fractional calculus, see [2, 8–10]. The author had made many efforts for modeling of leptospirosis epidemic disease since 2011; see [11–16]. The author proposed a mathematical model that describes the epidemic leptospirosis disease [11]. The model proposed by [11] is given by
(1)
d
S
h
d
t
=
b
1
-
μ
h
S
h
-
β
2
S
h
I
v
-
β
1
S
h
I
h
+
λ
h
R
h
,
d
I
h
d
t
=
β
2
S
h
I
v
+
β
1
S
h
I
h
-
μ
h
+
δ
h
+
γ
h
I
h
,
d
R
h
d
t
=
γ
h
I
h
-
μ
h
R
h
-
λ
h
R
h
,
d
S
v
d
t
=
b
2
-
γ
v
S
v
-
β
3
S
v
I
h
,
d
I
v
d
t
=
β
3
S
v
I
h
-
(
γ
v
+
δ
v
)
I
v
,
and the related initial conditions:
(2)
S
h
(
0
)
=
e
1
,
I
h
(
0
)
=
e
2
,
R
h
(
0
)
=
e
3
,
S
v
(
0
)
=
e
4
,
I
v
(
0
)
=
e
5
.
Here S
h
(
t
), I
h
(
t
), R
h
(
t
), S
v
(
t
), and I
v
(
t
), respectively, represent the population of susceptible human, infected human, recovered human, susceptible vector, and infected vector at time t. The rate at which the population of human increases is shown by b
1. The natural mortality rate for the human population is μ
h; β
1, β
2, and β
3 represent the transmission coefficients. The parameter λ
h shows the individuals who become susceptible again. The death from the disease that occurs to humans is shown by δ
h. The rate of recovery from infection for the human is denoted by γ
h. The growth rate of the vector population is represented by b
2; γ
v is the natural death rate for vector and disease related death rate for the vector is δ
v.The model presented by [11] is given by system (1) which represents a system of nonlinear ODE; their fractional order differential equation is given as follows:
(3)
D
t
α
S
h
t
=
b
1
-
μ
h
S
h
-
β
2
S
h
I
v
-
β
1
S
h
I
h
+
λ
h
R
h
,
D
t
α
I
h
t
=
β
2
S
h
I
v
+
β
1
S
h
I
h
-
μ
h
+
δ
h
+
γ
h
I
h
,
D
t
α
R
h
t
=
γ
h
I
h
-
μ
h
R
h
-
λ
h
R
h
,
D
t
α
S
v
t
=
b
2
-
γ
v
S
v
-
β
3
S
v
I
h
,
D
t
α
I
v
(
t
)
=
β
3
S
v
I
h
-
(
γ
v
+
δ
v
)
I
v
,
where D
t
α represents the fractional derivative in Caputo derivative sense and the order of the fractional derivative is shown by the parameter α and 0
<
α
<
1, with the related initial conditions (2). The parameter which describes the general response expression in representing the order of fractional derivative gives different results for different values. Obviously, the integer-order system can be viewed as a special case of the fractional-order system by putting the time-fractional order of the derivative equal to one. To put it simple, for the higher in order, the behavior of the fractional order system is the same in the case of integer order. According to research, this is the first work, which is available on the epidemic model of leptospirosis in fractional order. Due to this reason, we express the approximate solutions to the problems (3) by the technique of (MSGDTM). The solution to such nonlinear fractional equation can also be available by using an efficient approach [17].The MSGDTM is a modified form of the GDTM [6, 8, 18–20], in which it is treated as an algorithm in a sequence of small steps, to obtain the accurate approximate solution to the desire models. By using GDTM, the obtained solution is valid for short interval of time while the solution obtained from multistep generalized differential transform method (MSGDTM) is valid for a long time. To obtain the solution, using MSGDTM is more accurate and valid for long interval of time and agrees well with the classical Runge-Kutta numerical solution method, with the unity order derivative [6, 21].The paper has been organized/structured as follows. Section1 is introduction; Section 2 contains basic definitions of fractional calculus and their notations; Section 3 gives nonnegative solution of the model; Section 4 contains application of MSGDTM; Section 5 gives numerical results; and Section 6 is the conclusion.


## 2. Basic Definitions and Notations

This section contains some fundamental definitions and notations of fractional calculus that will help us in the proceeding sections [3, 22–24].Definition 1.
A functiong
(
x
) with the positive values of x is known in the space C
α  (
α
∈
R
) if it is written in the form g
(
x
)
=
x
q
g
1
(
x
) and for some q
>
α, where g
1
(
x
) is continuous in [
0
,
∞
), and it is known to be in the space C
α
n if g
(
n
)
∈
C
α, n
∈
N.Definition 2.
The Riemann Liouville integral operator of orderα
>
0 with a
≥
0 is defined as
(4)
J
a
α
g
x
=
1
Γ
α
∫
a
x
‍
x
-
t
α
-
1
g
t
d
t
,
x
>
α
,
(
J
a
0
g
)
(
x
)
=
g
(
x
)
.
For the properties of the operator [22], we need only the following. For g
∈
C
α, α
>
0, β
>
0, c
∈
R, and γ
>
-
1, one has
(5)
J
a
α
J
a
β
g
x
=
J
a
β
J
a
α
g
x
=
J
a
α
+
β
g
x
,
J
α
α
x
γ
=
x
γ
+
α
Γ
(
α
)
B
(
x
-
a
)
/
x
(
α
,
γ
+
1
)
,
where B
τ
(
α
,
γ
+
1
) represents the incomplete beta function defined as
(6)
B
τ
(
α
,
γ
+
1
)
=
∫
0
τ
t
α
1
1
-
t
γ
d
t
,
J
a
α
e
c
x
=
e
a
c
x
-
a
α
∑
k
=
0
∞
c
x
-
a
k
Γ
(
α
+
k
+
1
)
.
The Riemann Liouville derivative has some disadvantages when applying to real world model with fractional differential equations. Therefore, here we use a modified fractional differential operator D
a
α which is used in his work by Caputo on the theory of viscoelasticity.Definition 3.
The Caputo fractional derivative ofg
(
x
) of order α
>
0 with a
≥
0 is defined as
(7)
D
a
α
g
x
=
J
a
m
-
α
g
m
(
x
)
=
1
Γ
(
m
-
α
)
∫
a
x
‍
g
(
m
)
(
t
)
x
-
t
(
α
+
1
-
m
)
d
t
,
for m
-
1
<
α
≤
m, m
∈
N, x
≥
a, f
(
x
)
∈
C
-
1
m. Many authors investigated the Caputo fractional order derivative; for m
-
1
<
α
≤
m, g
(
x
)
∈
C
α
m, and α
≥
-
1, one gets
(8)
J
a
α
D
a
α
g
x
=
J
(
m
)
D
(
m
)
g
(
x
)
=
g
(
x
)
-
∑
k
=
0
m
-
1
g
(
k
)
(
a
)
x
-
a
k
k
!
.


## 3. Nonnegativity of the Model

SupposeR
+
5
=
{
X
♮
∈
R
5
:
X
♮
≥
0
} and X
♮
(
t
)
=
(
S
h
t
,
I
h
t
,
R
h
t
,
S
v
t
,
I
v
t
)
T. To prove the theorem, first, we state the following lemma [25].Lemma 4.
Supposeg
(
x
)
∈
C
[
a
,
b
] and D
α
g
(
x
)
∈
C
[
a
,
b
] for α
∈
(
0,1
]. Then one has
(9)
g
(
x
)
=
g
(
a
)
+
1
Γ
(
α
)
D
α
g
(
η
)
x
-
a
η
,
with 0
≤
η
≤
x, for all x
∈
[
a
,
b
]. This is also called generalized mean value theorem.Remark 5.
Supposeg
(
x
)
∈
C
[
a
,
b
] and D
α
g
(
x
)
∈
C
[
a
,
b
] for 0
<
α
≤
1. It is obvious from Lemma 4 that if D
α
g
(
x
)
≥
0  ∀
x
∈
(
0
,
b
), then the g is nonincreasing function.Theorem 6.
For the initial value problems (3), a unique solution exists and remains in R
+
5.Proof.
The existence and uniqueness of the solution of (3) in (
0
,
∞
) can be obtained from [26]. We just need to prove that the domain R
+
5 is positive invariant. Since
(10)
D
t
α
S
h
t
S
h
=
0
=
b
1
+
λ
h
R
h
t
≥
0
,
D
t
α
I
h
t
I
h
=
0
=
β
2
S
h
t
I
v
t
≥
0
,
D
t
α
R
h
t
R
h
=
0
=
γ
h
I
h
t
≥
0
,
D
t
α
S
v
t
S
v
=
0
=
b
2
≥
0
,
D
t
α
I
v
(
t
)
I
v
=
0
=
β
3
S
v
(
t
)
I
h
(
t
)
≥
0
,
on each hyperplane bounding the nonnegative orthant, the vector field points into R
+
5.


## 4. MSGDTM: Its Application

The general differences transform method (GDTM) is used to provide approximate solutions for nonlinear problems with respect to the convergent series with easily computable components. It has shown that the solution obtained is estimated to be invalid for larget for some systems [21, 27]. We therefore use multistep generalized differential transform method (MSGDTM) to approximate the solution of (3), which provides accurate solution for a period that is longer than the standard generalized difference transform method [21, 27]. Several authors used the MSGDTM; see, for example, [28, 29]. We take the differential transform of (3) with respect to time t as
(11)
S
h
k
+
1
=
Γ
α
k
+
1
Γ
α
k
+
1
+
1
×
b
1
-
∑
k
2
=
0
k
∑
k
1
=
0
k
2
B
2
k
1
S
h
k
2
-
k
1
I
v
k
-
k
2
llllllllllllll
-
μ
h
S
h
(
k
)
+
λ
h
R
h
(
k
)
llllllllllllll
-
∑
k
2
=
0
k
∑
k
1
=
0
k
2
B
1
k
1
S
h
k
2
-
k
1
I
h
k
-
k
2
,
I
h
k
+
1
=
Γ
α
k
+
1
Γ
α
k
+
1
+
1
×
∑
k
2
=
0
k
∑
k
1
=
0
k
2
B
2
k
1
S
h
k
2
-
k
1
I
v
k
-
k
2
llllllllllllll
-
(
μ
h
+
δ
h
+
γ
h
)
I
h
(
k
)
llllllllllllll
+
∑
k
2
=
0
k
‍
∑
k
1
=
0
k
2
‍
B
1
k
1
S
h
k
2
-
k
1
I
h
k
-
k
2
,
R
h
k
+
1
=
Γ
α
k
+
1
Γ
α
k
+
1
+
1
γ
h
I
h
k
-
μ
h
+
λ
h
R
h
k
,
S
v
k
+
1
=
Γ
α
k
+
1
Γ
α
k
+
1
+
1
×
-
∑
k
2
=
0
k
∑
k
1
=
0
k
2
B
3
k
1
S
v
k
2
-
k
1
I
h
k
-
k
2
b
2
-
γ
v
S
v
k
lllllllllll
-
∑
k
2
=
0
k
∑
k
1
=
0
k
2
B
3
k
1
S
v
k
2
-
k
1
I
h
k
-
k
2
,
I
v
k
+
1
=
Γ
(
α
k
+
1
)
Γ
(
α
(
k
+
1
)
+
1
)
×
∑
k
2
=
0
k
∑
k
1
=
0
k
2
B
3
k
1
S
v
k
2
-
k
1
I
h
k
-
k
2
llllllllll
-
γ
v
+
δ
v
I
v
k

∑
k
2
=
0
k
∑
k
1
=
0
k
2
B
3
k
1
S
v
k
2
-
k
1
I
h
k
-
k
2
,
where S
h
(
k
), I
h
(
k
), R
h
(
k
), S
v
(
k
), I
v
(
k
), B
1
(
k
), B
2
(
k
), and B
3
(
k
) represent the differential transformations of S
h
(
t
), I
h
(
t
), R
h
(
t
), S
v
(
t
), I
v
(
t
), β
1
(
t
), β
2
(
t
), and β
3
(
t
). The initial conditions in terms of differential transform form are represented as S
h
(
0
)
=
e
1, I
h
(
0
)
=
e
2, R
h
(
0
)
=
e
3, S
v
(
0
)
=
e
4, and I
v
(
0
)
=
e
5. From differential inverse transforms point of view, the differential transform series solution for the systems (3) can be obtained as
(12)
S
h
t
=
∑
k
=
0
K
S
h
k
t
α
k
,
I
h
t
=
∑
k
=
0
K
I
h
k
t
α
k
,
R
h
t
=
∑
k
=
0
K
R
h
k
t
α
k
,
S
v
t
=
∑
k
=
0
K
S
v
k
t
α
k
,
I
v
t
=
∑
k
=
0
K
I
v
k
t
α
k
.
The series solution for the systems (3), according to MSGDTM is suggested as
(13)
S
h
t
=
∑
k
=
0
K
S
1
h
(
k
)
t
α
k
,
t
∈
[
0
,
t
1
]
∑
k
=
0
K
S
2
h
(
k
)
t
-
t
1
α
k
,
t
∈
[
t
1
,
t
2
]
⋮
∑
k
=
0
K
S
M
h
(
k
)
t
M
-
t
M
-
1
α
k
,
t
∈
t
M
-
1
,
t
M
,
I
h
t
=
∑
k
=
0
K
I
1
h
(
k
)
t
α
k
,
t
∈
[
0
,
t
1
]
∑
k
=
0
K
I
2
h
(
k
)
t
-
t
1
α
k
,
t
∈
[
t
1
,
t
2
]
⋮
∑
k
=
0
K
I
M
h
(
k
)
t
M
-
t
M
-
1
α
k
,
t
∈
t
M
-
1
,
t
M
,
R
h
(
t
)
=
∑
k
=
0
K
R
1
h
(
k
)
t
α
k
,
t
∈
[
0
,
t
1
]
∑
k
=
0
K
R
2
h
(
k
)
t
-
t
1
α
k
,
t
∈
[
t
1
,
t
2
]
⋮
∑
k
=
0
K
R
M
h
(
k
)
t
M
-
t
M
-
1
α
k
,
t
∈
t
M
-
1
,
t
M
,
S
v
(
t
)
=
∑
k
=
0
K
S
1
v
(
k
)
t
α
k
,
t
∈
[
0
,
t
1
]
∑
k
=
0
K
S
2
v
(
k
)
t
-
t
1
α
k
,
t
∈
[
t
1
,
t
2
]
⋮
∑
k
=
0
K
S
M
v
(
k
)
t
M
-
t
M
-
1
α
k
,
t
∈
t
M
-
1
,
t
M
,
I
v
(
t
)
=
∑
k
=
0
K
I
1
v
(
k
)
t
α
k
,
t
∈
[
0
,
t
1
]
∑
k
=
0
K
I
2
v
(
k
)
t
-
t
1
α
k
,
t
∈
[
t
1
,
t
2
]
⋮
∑
k
=
0
K
I
M
v
(
k
)
t
M
-
t
M
-
1
α
k
,
t
∈
t
M
-
1
,
t
M
,
where S
i
h
(
k
), I
i
h
(
k
), R
i
h
(
k
), S
i
v
(
k
), and I
i
v
(
k
) for i
=
1,2
,
3
,
…
,
M satisfy the following recurrence relations:
(14)
S
i
h
k
+
1
=
Γ
α
k
+
1
Γ
α
k
+
1
+
1
×
b
1
-
∑
k
2
=
0
k
∑
k
1
=
0
k
2
B
2
k
1
S
i
h
k
2
-
k
1
I
i
v
k
-
k
2
llllllllllllll
-
μ
h
S
i
h
(
k
)
+
λ
h
R
i
h
(
k
)
llllllllllllll
-
∑
k
2
=
0
k
∑
k
1
=
0
k
2
B
1
k
1
S
i
h
k
2
-
k
1
I
i
h
k
-
k
2
,
I
i
h
k
+
1
=
Γ
α
k
+
1
Γ
α
k
+
1
+
1
×
∑
k
2
=
0
k
∑
k
1
=
0
k
2
B
2
k
1
S
i
h
k
2
-
k
1
I
i
v
k
-
k
2
llllll
llllllll
-
(
μ
h
+
δ
h
+
γ
h
)
I
i
h
(
k
)
llllllllllllll
+
∑
k
2
=
0
k
∑
k
1
=
0
k
2
B
1
k
1
S
i
h
k
2
-
k
1
I
i
h
k
-
k
2
,
R
i
h
k
+
1
=
Γ
α
k
+
1
Γ
α
k
+
1
+
1
γ
h
I
i
h
k
-
μ
h
+
λ
h
R
i
h
k
,
S
i
v
k
+
1
=
Γ
α
k
+
1
Γ
α
k
+
1
+
1
×
∑
k
2
=
0
k
∑
k
1
=
0
k
2
B
3
k
1
S
i
v
k
2
-
k
1
I
i
h
k
-
k
2
b
2
-
γ
v
S
i
v
k
llllllllll
-
∑
k
2
=
0
k
∑
k
1
=
0
k
2
B
3
k
1
S
i
v
k
2
-
k
1
I
i
h
k
-
k
2
,
I
i
v
k
+
1
=
Γ
(
k
+
1
)
Γ
α
α
k
+
1
+
1
×
∑
k
2
=
0
k
‍
∑
k
1
=
0
k
2
‍
B
3
k
1
S
i
v
k
2
-
k
1
I
i
h
k
-
k
2
lllllllllll
-
γ
v
+
δ
v
I
i
v
k

∑
k
2
=
0
k
‍
∑
k
1
=
0
k
2
‍
B
3
,
such that S
i
h
(
t
i
-
1
)
=
S
i
-
1
h
(
t
i
-
1
), I
i
h
(
t
i
-
1
)
=
I
i
-
1
h
(
t
i
-
1
), R
i
h
(
t
i
-
1
)
=
R
i
-
1
h
(
t
i
-
1
), S
i
v
(
t
i
-
1
)
=
S
i
-
1
v
(
t
i
-
1
), and I
i
v
(
t
i
-
1
)
=
I
i
-
1
v
(
t
i
-
1
). Initially starting from S
0
h
=
e
1, I
0
h
=
e
2, R
0
h
=
e
3, S
0
v
=
e
4, and I
0
v
=
e
5, with the use of (14). The solution can be obtained by the MSGDTM given in (13). Table 1 represents the relation between the original function and the differential transform.Table 1
The relation between the original function and the differential transform.
The original function
The transformed function
g
(
t
)
=
x
(
t
)
+
y
(
t
)
G
(
k
)
=
X
(
k
)
+
Y
(
k
)
g
(
t
)
=
α
x
(
t
)
G
(
k
)
=
α
X
(
k
)
g
(
t
)
=
x
(
t
)
y
(
t
)
G
(
k
)
=
∑
l
=
0
k
X
(
l
)
Y
(
k
-
l
)
g
(
t
)
=
d
x
(
t
)
d
(
t
)
G
(
k
)
=
(
k
+
1
)
X
(
k
+
1
)
g
(
t
)
=
d
m
x
(
t
)
d
t
m
G
(
k
)
=
(
k
+
1
)
(
k
+
2
)
⋯
(
k
+
m
)
X
(
k
+
m
)


## 5. Numerical Results

In this section, we present the numerical solution of the proposed model (3) numerically by the classical Runge-Kutta order four method for the integer case derivative and analytically we use the MSGDTM. We apply the algorithm on the interval [
0,30
] for the approximate solution of nonlinear fractional differential equations (3). For K
=
10 and M
=
3000, the results are obtained. The results are obtained by using computer algebra package Mathematica. The initial conditions e
1
=
40, e
2
=
20, e
3
=
10, e
4
=
50, and e
5
=
30 are used. The values for the parameters we assumed are presented in Table 2. In Figures 1, 2, 3, 4, and 5, the approximate solution is obtained by MSGDTM and classical Runge-Kutta method order four scheme for α
=
1. The approximate solution of MSGDTM has good agreement with the Runge-Kutta order four method. Figures 6, 7, 8, 9, and 10 represent the approximate solution obtained by MSGDTM and classical Runge-Kutta order four scheme considering different values of α. We conclude from our graphical results that the MSGDTM results behave as closely as Runge-Kutta method.Table 2
Parameter values used in the numerical simulations of the optimal control.
Notation
Parameter description
Value
Reference
b
1
Recruitment rate for human
1.6
[11]
λ
h
Proportionality constant
0.066
[30]
μ
h
Natural death rate of human
4.6
×
1
0
-
5
[11]
γ
v
Natural death rate of vector
1.8
×
1
0
-
3
[11]
δ
h
Death rate due to disease at human class
1.0
×
1
0
-
5
[31]
γ
h
Recovery rate of the infection
2.7
×
1
0
-
3
[31]
b
2
Birth rate of vector
1.9
Assumed
β
2
Transmission betweenS
h and I
v
0.0089
Assumed
β
3
Transmission betweenS
v and I
h
0.0079
Assumed
β
1
Transmission coefficient betweenS
h and I
h
0.00013
Assumed
γ
v
Natural death rate of vector
0.0027
[11]Figure 1
The dotted line and the solid line connectively show the Runge-Kutta method and MSGDTM.Figure 2
Comparison of RK method (dotted line) and MSGDTM (solid line).Figure 3
Comparison of RK method (dotted line) and MSGDTM (solid line).Figure 4
Comparison of RK method (dotted line) and MSGDTM (solid line).Figure 5
Comparison of RK method (dotted line) and MSGDTM (solid line).Figure 6
Susceptible human: solid lineα
=
1.0, +
+ line (
α
=
0.95
), and *
* line (
α
=
0.90
).Figure 7
Infected human: solid lineα
=
1.0, +
+ line (
α
=
0.95
), and *
* line (
α
=
0.90
).Figure 8
Recovered human: solid lineα
=
1.0, +
+ line (
α
=
0.95
), and *
* line (
α
=
0.90
).Figure 9
Susceptible vector: solid lineα
=
1.0, +
+ line (
α
=
0.95
), and *
* line (
α
=
0.90
).Figure 10
Infected vector: solid line(
α
=
1.0, +
+ line (
α
=
0.95
), and *
* line (
α
=
0.90
).


## 6. Conclusion

In this paper, we have presented an approximate solution to an epidemic model of leptospirosis in fractional order. First, we obtained the analytical solution to the fractional order of epidemic model. Then, we solved the fractional order system and the results were compared with classical Runge-Kutta method. The results obtained from MGDTM have a good agreement with Runge-Kutta order four method. The MGDTM results are valid for a largert.



---

*Source: 317201-2014-12-11.xml*