# A Fractional-Order Model for HIV Dynamics in a Two-Sex Population

Authors: Fatmawati; Endrik Mifta Shaiful; Mohammad Imam Utoyo
Journal: International Journal of Mathematics and Mathematical Sciences
 (2018)

---

## Abstract

Human Immunodeficiency Virus (HIV) is a virus that attacks or infects cells in the immune system that causes immune decline. Acquired Immunodeficiency Syndrome (AIDS) is the most severe stage of HIV infection. AIDS is the rapidly spreading and becoming epidemic diseases in the world of almost complete influence across the country. A mathematical model approach of HIV/AIDS dynamic is needed to predict the spread of the diseases in the future. In this paper, we presented a fractional-order model of the spread of HIV and AIDS diseases which incorporates two-sex population. The fractional derivative order of the model is in the interval(
0,1
]. We compute the basic reproduction number and prove the stability of the equilibriums of the model. The sensitivity analysis also is done to determine the important factor controlling the spread. Using the Adams-type predictor-corrector method, we then perform some numerical simulations for variation values of the order of the fractional derivative. Finally, the effects of various antiretroviral therapy (ART) treatments are studied and compared with numerical approach.


---

## Body

## 1. Introduction

Human Immunodeficiency Virus (HIV) is a virus that attacks or infects cells in the immune system that causes immune decline. Acquired Immunodeficiency Syndrome (AIDS) is the most severe stage of HIV infection, which can take from 2 to 15 years to develop depending on the individual. AIDS is defined by the development of certain cancers, infections, or other severe clinical manifestations. HIV can be transmitted via the exchange of a variety of body fluids from infected individuals, such as blood, breast milk, semen, and vaginal secretions. There were approximately 36.7 million people living with HIV at the end of 2016 with 1.0 million people died from HIV-related causes globally. In 2015, an estimated 44% of new infections occurred among key populations and their partners [1].Until now, there is no cure for HIV infection. However, effective antiretroviral therapy (ART) can inhibit HIV progression in immune defects. Since 1996, ART has begun to be used by people living with HIV in the world because it can prevent death early [2]. The benefits of ART for people with HIV/AIDS are to improve quality of life, prevent mother-to-child transmission of HIV, prolong survival, and restore the immune system [3]. In 2016, 19.5 million people with HIV have been receiving ART globally [1].Mathematical models are needed to understand the dynamics of epidemic infection [4–8]. At present many models have been proposed to describe the dynamics of HIV/AIDS infection [9–11]. For instance, authors in [9] formulated a mathematical model for the transmission dynamics of HIV/AIDS in a two-sex population considering counseling and antiretroviral therapy. Authors in [10] studied the impact of the optimal control on the treatment of HIV/AIDS incorporating use of condom, screening of unaware infective, and treatment of HIV individuals. In recent year, Yang et al. [11] formulated a two-group (female sex workers and senior male clients) compartmental model to study the impact of senior male clients on the transmission dynamics, the containment, and the elimination of the HIV.Memory effect plays an important role in the spread of disease. The presence of memory effects on past events will affect the spread of disease in the future so that the spread of disease in the future can be controlled. The distance of memory effect indicates the history of disease spread. Thus, memory effects on the spread of infectious diseases can be investigated using fractional derivatives. Fractional derivatives have been used in the literature to observe the effect of memory on a system dynamics by replacing the ordinary derivative order with the fractional derivative order [12–15]. Sardar et al. [13] formulated the dengue model with memory in the transmission process by using fractional differential operator and order of the fractional derivative as an index of memory. Huo et al. [12] analyzed the effects of vaccines on mathematical models of fractional order of HIV disease spread. They also performed a local stability analysis on the fractional-order framework of HIV disease spread and bifurcation behavior of the system. In 2017, Saeedian et al. [14] studied the evolution of the SIR epidemic model, considering memory effects. Using the fractional calculus technique, the authors in [14] show that the dynamics of such a system depend on the strength of memory effects, controlled by the order of fractional derivatives. In 2017, Pinto and Carvalho [15] derived a fractional-order model for the dynamics of the coinfection of HIV and TB in the presence of TB resistant strains.In this paper, we proposed a fractional order of two-sex mathematical model for dynamic HIV transmission, as a generalization of an integer model, proposed by Kimbir et al. [9]. We also modified the basic model in [9] by distinguishing populations infected with HIV and infected with AIDS. By using fractional ordinary differential equation, we hope that the fractional-order model can accommodate the real phenomenon of the spread of HIV/AIDS. The structure of the paper is organized as follows. In Section 2, we introduce the description of the fractional order of the HIV/AIDS model. In Section 3, the stability analysis of the equilibriums of the model is proven. In Section 4, we carry out the sensitivity analysis of the reproduction number to the parameters in the model. Section 5 shows some numerical results for distinct values of the fractional order α. Finally, in Section 6, we give a brief conclusion.


## 2. Description of the Model

In this section we propose a mathematical model of two-sex HIV/AIDS transmission based on [9]. The model is constructed under the following assumptions:(1)
The total population is divided into two groups, namely, the population of males and females.
(2)
The population of males is divided into three compartments: the susceptible males(
S
m
), the HIV infected males (
I
m
), and the AIDS infected males (
A
m
). The population of females is also divided into three compartments: the susceptible females (
S
f
), the HIV infected females (
I
f
), and the AIDS infected females (
A
f
). Moreover, the total population of males and females is given by N
=
S
m
+
I
m
+
A
m
+
S
f
+
I
f
+
A
f.
(3)
The HIV transmission is by heterosexual contact only.
(4)
Vertical transmission and age-structure are ignored.
(5)
The population is homogeneous mixing.
(6)
Male and female AIDS populations were isolated and then do not contribute to HIV infection.The transmission diagram of the model is shown in Figure1.Figure 1
HIV/AIDS transmission diagram.The basic model is derived as follows:(1)
d
S
m
d
t
=
p
Λ
-
β
f
I
f
S
m
-
μ
S
m
.
d
I
m
d
t
=
β
f
I
f
S
m
-
1
-
r
γ
m
I
m
-
μ
I
m
,
d
A
m
d
t
=
1
-
r
γ
m
I
m
-
μ
+
δ
A
m
,
d
S
f
d
t
=
1
-
p
Λ
-
β
m
I
m
S
f
-
μ
S
f
,
d
I
f
d
t
=
β
m
I
m
S
f
-
1
-
r
γ
f
I
f
-
μ
I
f
,
d
A
f
d
t
=
1
-
r
γ
f
I
f
-
μ
+
δ
A
f
,The description of the parameter for model (1) could be seen in Table 1. The biologically feasible region of model (1) is (2)
Ω
=
S
m
,
I
m
,
A
m
,
S
f
,
I
f
,
A
f
∈
R
+
6
:
0
≤
N
≤
Λ
μ
,and all of the parameters used in model (1) are nonnegative.Table 1
Parameters of model (1).
Description
Parameter
Recruitment rate
Λ
Proportion of the recruitment rate
p
Natural death rate
μ
Transmission rate by an infected male
β
m
Transmission rate by an infected female
β
f
Progression rate from male HIV infection to AIDS
γ
m
Progression rate from female HIV infection to AIDS
γ
f
AIDS disease induced death rate
δ
Proportion of the efficacy of ART treatment for HIV infection
rThe regionΩ is positively invariant. In this region, model (1) is well-posed. So, if it is given an initial condition in the region, then the solution is defined for all time t
≥
0 and remains in the region.Next, we consider a fractional-order model of system (1). The fractional model corresponding to system (1) is as follows: (3)
d
α
S
m
d
t
α
=
p
Λ
-
β
f
I
f
S
m
-
μ
S
m
,
d
α
I
m
d
t
α
=
β
f
I
f
S
m
-
1
-
r
γ
m
I
m
-
μ
I
m
,
d
α
A
m
d
t
α
=
1
-
r
γ
m
I
m
-
μ
+
δ
A
m
,
d
α
S
f
d
t
α
=
1
-
p
Λ
-
β
m
I
m
S
f
-
μ
S
f
,
d
α
I
f
d
t
α
=
β
m
I
m
S
f
-
1
-
r
γ
f
I
f
-
μ
I
f
,
d
α
A
f
d
t
α
=
1
-
r
γ
f
I
f
-
μ
+
δ
A
f
,where α
∈
(
0,1
] is the order of the fractional derivative. Fractional derivative of model (3) is in the sense of Caputo. The Caputo approach is mostly used in real applications. The main advantages of Caputo approach are the initial values for fractional differential equations with the Caputo derivatives taking on the same form as for integer order differential equations [16]. The Caputo fractional derivative is defined as follows.Definition 1 (see [16]).
The Caputo fractional differential operator of orderα
>
0, with n
-
1
<
α
<
n
,
n
∈
N, is defined by (4)
d
α
f
t
d
t
α
≔
I
n
-
α
d
n
f
t
d
t
n
=
1
Γ
n
-
α
∫
0
t
t
-
s
n
-
α
-
1
f
n
s
d
s
,where Γ
(
·
) is the gamma function.


## 3. Model Analysis

In this section, we study the stability of the equilibrium of the fractional-order model (3). We begin by computing the basic reproduction number (
R
0
) of model (3). The basic reproduction number is defined as the number of secondary cases of primary case during the infectious period due to the type of infection [17, 18].Now, we recall the properties of the stability of the fractional-order systems. The stability theorem on fractional-order system is as follows.Theorem 2 (see [19, 20]).
Consider the following autonomous nonlinear fractional-order system:(5)
d
α
x
d
t
α
=
f
x
,
x
0
=
x
0
,with 0
<
α
<
1 and x
∈
R
n. The equilibrium points of the above system are solutions to the equation f
(
x
)
=
0. An equilibrium is locally asymptotically stable if all eigenvalues (
λ
j
) of the Jacobian matrix J
=
∂
f
/
∂
x evaluated at the equilibrium satisfy arg
⁡
(
λ
j
)
>
α
π
/
2.Based on Theorem2, the equilibria are obtained by setting the right-hand sides of the equations in model (3) to zero. The disease-free equilibrium of model (3) to the coordinate (
S
m
,
I
m
,
A
m
,
S
f
,
I
f
,
A
f
) is given by E
0
=
(
Λ
p
/
μ
,
0,0
,
Λ
(
1
-
p
)
/
μ
,
0,0
).Then, the basic reproduction number(
R
0
) is computed by using the next-generation method [21, 22]. For the next-generation matrix method [22], we take the infected compartments (
I
m
,
I
f
). The Jacobian matrices F and V for the new infection in the compartment and the transfer of individuals between the compartment respectively, evaluated at E
0, are given by(6)
F
=
0
β
f
Λ
p
μ
β
m
Λ
1
-
p
μ
0
,
V
=
1
1
-
r
γ
m
+
μ
0
0
1
1
-
r
γ
f
+
μ
.The basic reproduction number of model (3) is the spectral radius of the matrix F
V
-
1 such that we have (7)
R
0
=
Λ
μ
β
m
β
f
p
1
-
p
1
-
r
γ
m
+
μ
1
-
r
γ
f
+
μ
.The following theorem provides the local stability of the disease-free equilibrium.Theorem 3.
The disease-free equilibriumE
0 is locally asymptotically stable if R
0
<
1.Proof.
The Jacobian matrix of model (3) around the disease-free equilibrium, E
0, is given by (8)
J
=
-
μ
0
0
0
-
β
f
p
Λ
μ
0
0
-
1
-
r
γ
m
-
μ
0
0
β
f
p
Λ
μ
0
0
1
-
r
γ
m
-
μ
-
δ
0
0
0
0
β
m
Λ
p
-
1
μ
0
-
μ
0
0
0
-
β
m
Λ
p
-
1
μ
0
0
-
1
-
r
γ
f
-
μ
0
0
0
0
0
1
-
r
γ
f
-
μ
-
δ
.The eigenvalues of matrix J are λ
1
=
λ
2
=
-
μ and λ
3
=
λ
4
=
-
(
μ
+
δ
) and the roots of quadratic equation x
2
+
a
1
x
+
a
2
=
0, where a
1
=
2
(
1
-
r
)
μ
(
γ
m
+
γ
f
) and a
2
=
[
(
1
-
r
)
2
γ
m
γ
f
+
μ
(
1
-
r
)
(
γ
m
+
γ
f
)
+
μ
2
]
(
1
-
T
0
), with R
0
=
T
0. Thus, we have arg
⁡
(
λ
1
)
=
arg
⁡
(
λ
2
)
=
arg
⁡
(
λ
3
)
=
arg
⁡
(
λ
4
)
=
π
>
α
π
/
2. Next, we check the roots of the quadratic equation x
2
+
a
1
x
+
a
2
=
0. Authors in [23] show that the Routh-Hurwitz criteria, a
1
,
a
2
>
0, are necessary and sufficient for arg
⁡
(
λ
i
)
>
α
π
/
2. It is clear that all of the eigenvalues are negative (arg
⁡
(
λ
j
)
>
α
π
/
2, for j
=
1,2
,
…
,
6) if T
0
<
1 or equivalently R
0
<
1. Hence, the disease-free equilibrium E
0 is locally asymptotically stable for α
∈
(
0,1
] if R
0
<
1.We continue with the computing of the endemic equilibrium(
E
1
) of model (3). The endemic equilibrium E
1 is given by (9)
E
1
=
S
m
∗
,
I
m
∗
,
A
m
∗
,
S
f
∗
,
I
f
∗
,
A
f
∗
,where (10)
S
m
∗
=
p
Λ
β
f
I
f
∗
+
μ
,
I
m
∗
=
1
-
r
γ
m
+
μ
1
-
r
γ
f
+
μ
μ
2
T
0
-
1
β
m
1
-
r
γ
m
+
μ
β
f
Λ
1
-
p
+
1
-
r
γ
f
+
μ
μ
,
A
m
∗
=
1
-
r
γ
m
I
m
∗
μ
+
δ
,
S
f
∗
=
1
-
p
Λ
β
m
I
m
∗
+
μ
,
I
f
∗
=
β
m
I
m
∗
1
-
p
Λ
1
-
r
γ
f
+
μ
β
m
I
m
∗
+
μ
,
A
f
∗
=
1
-
r
γ
f
I
f
∗
μ
+
δ
,with R
0
=
T
0. The endemic equilibrium E
1 exists if T
0
>
1 or equivalently R
0
>
1.The stability of the endemic equilibriumE
1 is difficult to prove analytically, because it involves a quartic equation which depend on the variables I
m and I
f. Numerical simulations show that the endemic equilibrium is locally asymptotically stable if R
0
>
1. This can be seen in Figures 2 and 3. Using three different initial conditions for the simulation, these orbits converge to the same point as time evolves.Figure 2
Phase portrait of model (3) in S
m-S
f plane for α
=
0.8.Figure 3
Phase portrait of model (3) in S
m-S
f plane for α
=
0.9.


## 4. Sensitivity Analysis

In this section we present the sensitivity analysis of the reproduction numberR
0 to the parameters in model (3). The aim of this analysis was to measure the parameters that have the most effects on the reproduction number. We derived analytically the sensitivity index of R
0 to each parameter following the approach in [24].Definition 4 (see [24]).
The normalized forward sensitivity index of a variable,R
0, that depends differentially on a parameter, l, is defined as (11)
Υ
l
R
0
≔
∂
R
0
∂
l
×
l
R
0
.Based on Definition4, the sensitivity indices of R
0 with respect to each parameter such as Λ
,
μ
,
β
m
,
β
f
,
γ
m
,
γ
f, and p can be computed in the same way as (11). For example, the sensitivity index of R
0 with respect to Λ is (12)
Υ
Λ
R
0
≔
∂
R
0
∂
Λ
×
Λ
R
0
=
1
.Thus, we compute the sensitivity indexes of the remaining parameters using the parameter values in the Table 2. The results are given in Table 3.Table 2
Parameter values for simulations.
Parameter
Value
Ref.
p
0.5
Assumed
Λ
20
[12]
β
m
2
×
1
0
-
3
Assumed
β
f
5
×
1
0
-
3
Assumed
μ
0.02
[12]
δ
0.125
[12]
γ
m
8
×
1
0
-
3
Assumed
γ
f
9
×
1
0
-
3
Assumed
r
0.5
AssumedTable 3
Sensitivity indices to parameter for model (3).
Parameter
Sensitivity index
Λ
1
β
m
0.5
β
f
0.5
r
0.175
μ
−1.825
γ
f
−0.092
γ
m
−0.083
p
0The sensitivity index can be analyzed as follows. The positive sensitivity index shows that an increase in the parameters will lead to increase in the basic reproduction number, while a negative sensitivity index means that an increase in the parameter will lead to a decrease in the basic reproduction number. For example, forΥ
β
m
R
0
=
0.5, increasing the value β
m by 10% increases the reproduction number R
0 by 5%. Thus, increasing natural death rate μ by 10% decreases R
0 by 18.25%.We also perform sensitivity simulation to verify our sensitivity analysis. The parameter values used in the simulations are given in Table2. In Figures 4 and 5, we can see that, for the parameter chosen with distinct values of β
f and β
m, respectively, R
0 increases monotonically with both β
m and β
f. This results indicate that increasing β
m and β
f will increase the basic reproduction number R
0.Figure 4
Sensitivity ofR
0 with respect to β
m for different values of β
f.Figure 5
Sensitivity ofR
0 with respect to β
f for different values of β
m.


## 5. Numerical Simulation

In this section, we conduct several numerical simulations of model (3). An Adams-type predictor-corrector method [25–27] is applied to solve the numerical solution of the fractional-order model (3). Parameters values used in these simulations could be seen in Table 2. The adopted initial conditions are (
S
m
(
0
)
,
I
m
(
0
)
,
A
m
(
0
)
,
S
f
(
0
)
,
I
f
(
0
)
,
A
f
(
0
)
)
=
(
150,10,6
,
300,7
,
2
). Here, we take 500 days for the time horizon. The simulations are carried out with varying values of the order of the fractional derivative α
∈
[
0.5,1.0
].In Figures6–11, we display the dynamic of the populations when R
0
<
1 and in Figures 12–17 when R
0
>
1. In each figure six different values of α are employed.Figure 6
Dynamic of nonendemicS
m.Figure 7
Dynamic of nonendemicS
f.Figure 8
Dynamic of nonendemicI
m.Figure 9
Dynamic of nonendemicI
f.Figure 10
Dynamic of nonendemicA
m.Figure 11
Dynamic of nonendemicA
f.Figure 12
Dynamic of endemicS
m.Figure 13
Dynamic of endemicS
f.Figure 14
Dynamic of endemicI
m.Figure 15
Dynamic of endemicI
f.Figure 16
Dynamic of endemicA
m.Figure 17
Dynamic of endemicA
f.Now, we setβ
m
=
2
×
1
0
-
5 and β
f
=
5
×
1
0
-
5 and the remaining of parameters as in Table 2 to simulate Figures 6–11. In this case, the value of R
0 is R
0
=
0.6521
<
1 which means that the infection will die out in the population. This condition is confirmed by simulation results in Figures 6 and 7, Figures 8 and 9, and Figures 10 and 11 for susceptible male and female populations, respectively, infected by HIV and infected by AIDS. These figures show that the solutions of model (3) are convergent to the disease-free equilibrium. Moreover, we observe that the solutions with higher order α have faster convergence speed compared to the smaller.Next, we plot in Figures12–17 the numerical simulations using the parameters as in Table 2. In this condition, the value of R
0 is R
0
=
65.2051
>
1 which means that the infection will persist in the population. In Figures 12 and 13, Figures 14 and 15, and Figures 16 and 17 the dynamics of the susceptible male and female populations, respectively, are seen, infected by HIV and infected by AIDS with different values of α. The figures show that the solutions of model (3) converge to the endemic equilibrium when R
0
>
1. Similar to the nonendemic condition, we see that as the order α increases, the convergence of solutions is faster for the endemic condition.In Figures18–23, we depict the dynamic of both male and female populations infected by AIDS for different values of the efficacy of ART treatment (
r
) for α
=
1, α
=
0.7, and α
=
0.5. For α
=
1, we observe that the AIDS infected both male and female populations decrease when the value of r is increase. On the contrary, the number of both AIDS infections increases when the efficacy of ART treatment, r, is smaller. The similar behavior is seen for α
=
0.7 and α
=
0.5. It is well known that the ART treatment could improve the quality of HIV infected both male and female patients.Figure 18
Dynamic ofA
m with α
=
1 for various r.Figure 19
Dynamic ofA
f with α
=
1 for various r.Figure 20
Dynamic ofA
m with α
=
0.7 for various r.Figure 21
Dynamic ofA
f with α
=
0.7 for various r.Figure 22
Dynamic ofA
m with α
=
0.5 for various r.Figure 23
Dynamic ofA
f with α
=
0.5 for various r.


## 6. Conclusion

In this paper, we have investigated a fractional order of two-sex mathematical model for dynamic HIV, as a generalization of an integer order model, proposed by Kimbir et al. [9]. The basic model in [9] is modified by distinguishing populations infected with HIV and AIDS. We have computed the basic reproduction number (
R
0
) and proved the stability of equilibriums of the fractional-order model of the HIV infection. Based on the mathematical analysis, the disease-free equilibrium is locally asymptotically stable when R
0
<
1 that means the infection will die out in the population. Numerically, the endemic equilibrium tends to be locally asymptotically stable when R
0
>
1 which means that the infection will persist in the population. We also studied analytically and numerically the sensitivity analysis to measure the parameters that have high impact on R
0. Finally, we have carried out the numerical simulations for different values of the order (
α
) of the fractional derivative. The simulations results show that the solutions with higher order α have faster convergence compared to the smaller α. We also found that as the efficacy of ART treatment (
r
) increases there is a corresponding decrease in the number of the AIDS infected both male and female populations for three values of α. These results indicate the effectiveness of the ART treatment to reduce the AIDS infected.



---

*Source: 6801475-2018-04-18.xml*