# A Stochastic Differential Equation Model for the Spread of HIV amongst People Who Inject Drugs

Authors: Yanfeng Liang; David Greenhalgh; Xuerong Mao
Journal: Computational and Mathematical Methods in Medicine
 (2016)

---

## Abstract

We introduce stochasticity into the deterministic differential equation model for the spread of HIV amongst people who inject drugs (PWIDs) studied by Greenhalgh and Hay (1997). This was based on the original model constructed by Kaplan (1989) which analyses the behaviour of HIV/AIDS amongst a population of PWIDs. We derive a stochastic differential equation (SDE) for the fraction of PWIDs who are infected with HIV at time. The stochasticity is introduced using the well-known standard technique of parameter perturbation. We first prove that the resulting SDE for the fraction of infected PWIDs has a unique solution in (0, 1) provided that some infected PWIDs are initially present and next construct the conditions required for extinction and persistence. Furthermore, we show that there exists a stationary distribution for the persistence case. Simulations using realistic parameter values are then constructed to illustrate and support our theoretical results. Our results provide new insight into the spread of HIV amongst PWIDs. The results show that the introduction of stochastic noise into a model for the spread of HIV amongst PWIDs can cause the disease to die out in scenarios where deterministic models predict disease persistence.


---

## Body

## 1. Introduction

HIV (human immunodeficiency virus) is a deadly and infectiousLentivirus which attacks and weakens the immune system by especially attacking the CD4 cells. As a result, HIV causes AIDS (acquired immune seficiency syndrome). Since the first discovery of HIV in 1981, it has already infected almost 78 million people with about 39 million lives having been taken [1]. Despite the massive improvement in technology and medical equipment, we are still unable to fully find a cure for the HIV virus. In 2014, according to the reports by the World Health Organization, there were still approximately 36.9 million people living with HIV, with around 2 million new cases globally [2]. In order to control the epidemic, it is crucial to understand the dynamical behaviour of HIV and how it spreads within our community. There are various routes by which HIV can be transmitted, for example, transmission via unprotected sexual intercourse, vertically from infected mothers to their unborn children and people who inject drugs (PWIDs) sharing contaminated needles. Amongst all the possible routes of HIV transmission, PWIDs have become a significant risk group with around 3 million of them living with HIV [3]. For every 10 new cases of HIV infection, on average, one of them is caused by injecting drug use. In regions of Central Asia and Eastern Europe, injecting drug use accounts for 80 percent of HIV infections [3]. As a result, in this paper, we will focus on looking at this particular risk group.Over the past years, mathematical models have been used successfully to analyse and predict the dynamical behaviour in biological systems. The first mathematical model for the spread of HIV and AIDS amongst PWIDs in shooting galleries was created by Kaplan [4], where a shooting gallery is a place for PWIDs to purchase and inject drugs. Kaplan incorporated many factors into his model such as the injection equipment sharing rate and the effect of cleaning injection equipment in order to better understand how HIV is transmitted within this type of community. Based on the original model created in [4], Greenhalgh and Hay [5] modified the model by changing some of the assumptions made by Kaplan to make them more realistic. These assumptions include having different visiting rates to the shooting galleries for PWIDs who have been diagnosed positive for the HIV virus and thus may have been advised to stop sharing injections and for those who either are not HIV positive or are but do not know it. The modified Kaplan model in [5] also allows for the possibility that an infected PWID may not always leave a needle infected before cleaning, as well as introducing different transmission probabilities for flushed and unflushed needles. The term “flushing” refers to the process where an infectious piece of injecting equipment is used by an uninfected PWID and thus after injecting the syringe is left uninfected. The HIV model that we will be looking at in this paper is based on the modified Kaplan model given in [5]. There have been many papers that have already looked at the connection between the spread of HIV and PWIDs [6–11]. However the models used are all deterministic models.The real world is not deterministic, and there are many factors that can influence the behaviour of a disease and thus it is not always possible to predict with certainty what would happen. Consequently, a stochastic model would be more appropriate. Furthermore, a stochastic model also possesses many useful unique properties such as being able to calculate the probability that an endemic will not occur and the expected duration of an endemic. By running a stochastic model many times, we can also build up a distribution of the possible outcomes which allows us to identify the number of infectives at a particular timet, whereas for a deterministic model we will only get a single outcome. HIV infection is a behavioural disease and thus there are many environmental factors that can influence the spread of HIV. Rhodes et al. [12] have mentioned in detail how factors such as injecting environments, social network, and neighbourhood deprivation and poverty can affect the spread of HIV amongst PWIDs. There are also other papers which emphasised how the dynamical behaviour of HIV is highly correlated with other factors [13–15]. Consequently, it is crucial for us to understand how HIV would spread under those environmental influences, especially amongst PWIDs. In this case, a stochastic model would be useful. There is also natural biological variation within people in their response to HIV. Using a stochastic model with environmental perturbation in the disease transmission parameter as we will do is one way to include this.The stochastic aspects of the HIV model have been studied by many authors. For example, in [16], Dalal et al. considered a stochastic model for internal HIV dynamics. They incorporated environmental stochasticity into their model by using the standard technique of parameter perturbation. They proved that the solution (representing the concentrations of uninfected cells, infected cells, and virus particles) is nonnegative and have looked at the stability aspect of their model by establishing the conditions required in order for the numbers of infected cells and virus particles to tend asymptotically to zero exponentially almost surely. Ding et al. [13] looked at a stochastic model for AIDS transmission and control taking into consideration the treatment rate of HIV patients. They have also examined the effect that knowledge, attitude, and behaviour of patients have on the spread of AIDS. Tuckwell and Le Corfec [17] used a stochastic model to analyse the behaviour of HIV-1 but focus only on the early stage after infection. Dalal et al. [18] have also used a stochastic model to look at another aspect of HIV. Once again, by using parameter perturbation, they introduced environmental randomness into their HIV model which allows them to examine the effect that condom use has on the spread of AIDS among a homogeneous homosexual population which is split into distinct risk groups according to the tendency of individuals to use condoms. Peterson et al. [19] constructed a population-based simulation of a community of PWIDs using the Monte Carlo technique. Greenhalgh and Lewis [20] modelled the spread of disease using a set of behavioural assumptions due to Kaplan and O’Keefe [10]. They use a branching process approximation to show that if the basic reproduction number R
0 is less than or equal to unity then the disease will always go extinct. They calculate an expression for the probability of extinction. They discuss an extended model which incorporates a three-stage incubation period and again examine a branching process approximation. They then compare them to investigate whether the deterministic model provides a good approximation to the simulated stochastic model. Although there have been many papers that looked at the stochastic aspect of the spread of HIV, as far as we know there are not many studies that focus on the stochastic aspect of the spread of HIV amongst PWIDs despite this particular risk group being responsible for many new HIV cases around the world. Thus it is crucial for us to examine the effect of environmental noise on this type of community.Inspired by the model constructed in [5], in this paper we will introduce environmental stochasticity into the model by parameter perturbation which is a standard technique in stochastic population modelling [16, 18, 21, 22]. To the best of our knowledge, this is the first paper which examines the effect that environmental stochasticity has on the dynamical behaviour of the modified Kaplan model [5]. The techniques used in this paper are inspired by the work done in [21]. The paper is organised as follows. In the next section, we will describe the formulation of the stochastic HIV model amongst PWIDs. In Section 3, we shall prove the existence of a unique nonnegative solution. In Sections 4 and 5, we will investigate two of the main important properties of any biological system, namely, the conditions required for extinction and persistence, respectively. Then in Section 6, we shall show that there exists a stationary distribution for our system. Finally, we will perform some numerical simulations with realistic parameter values to verify the results.


## 2. The Stochastic HIV Model

Throughout this paper, we let (Ω, F, {
F
t
}
t
≥
0, P) be a complete probability space with filtration {
F
t
}
t
≥
0 satisfying the usual conditions (i.e., it is increasing and right continuous while F
0 contains all P-null sets). Let us consider the following deterministic HIV model, which has been constructed by Greenhalgh and Hay [5] based on the model of Kaplan [4]. Define the following parameters:λ
1: shooting gallery visiting rate for susceptible PWIDs and the PWIDs who are infected but do not know they are infected.
λ
2: shooting gallery visiting rate for infected PWIDs who know that they are infected.
P
1: probability that the needle is flushed and the PWID is infected.
P
2: probability that the needle is flushed and the PWID remains uninfected.
P
3: probability that the PWID becomes infected without the needle being flushed.
P
4: probability that the PWID remains uninfected and the needle is not flushed.
ϕ
1: probability that an infected PWID leaves uninfected a syringe that was initially uninfected.
θ
1: probability that an infected PWID leaves uninfected a syringe that was initially infected.
ξ: fraction of all PWIDs (susceptible or not) who bleach their injection equipment after use.
γ: gallery ratio, where γ
=
n
/
m and n represents the PWID population and m represents the number of shooting galleries or syringes that each PWID visits at random.
p: probability that infected PWIDs know that they are infected.
μ: per capita rate at which infected PWIDs cease to share injection equipment (including those who cease sharing because of developing AIDS).Note that P
1
+
P
2
+
P
3
+
P
4
=
1.Define the following new composite parameters:(1)
σ
=
λ
1
1
-
p
+
λ
2
p
γ
1
-
ξ
1
-
ϕ
1
,
η
=
λ
1
1
-
p
+
λ
2
p
γ
ξ
+
θ
1
1
-
ξ
,
ρ
=
λ
1
γ
1
-
1
-
ξ
1
-
P
1
-
P
2
,
υ
=
λ
1
P
1
+
P
3
.In the expression for σ the factor (
1
-
ξ
)
(
1
-
ϕ
1
) represents the probability that an initially uninfected syringe is left infected and not cleaned by an infected PWID. The term λ
1
(
1
-
p
)
+
λ
2
p represents the average rate at which an infected PWID visits syringes. Hence σ
=
γ
σ
¯, where γ
=
n
/
m is the gallery ratio and σ
¯ is the rate at which an infected PWID visits syringes multiplied by the probability that he or she leaves an uninfected syringe infected after use. Similarly η
=
γ
η
¯, where η
¯ is the rate at which an infected PWID visits syringes multiplied by the probability that he or she leaves an infected syringe uninfected after use.λ
1 represents the rate at which a susceptible PWID visits syringes and 1
-
(
1
-
ξ
)
(
1
-
P
1
-
P
2
) represents the probability that an initially infected syringe is left uninfected after use by that PWID. Hence ρ
=
γ
ρ
¯, where ρ
¯ is the rate at which a susceptible PWID visits syringes multiplied by the probability that he or she leaves an infected syringe uninfected after use. υ represents the rate at which a susceptible PWID visits syringes multiplied by the probability that he or she becomes infected given that the syringe which they visit is infected. Thus υ
β represents the rate at which a susceptible PWID visits syringes and becomes infected. υ can thus be regarded as the “potential” infection rate of a susceptible PWID.Letπ
(
t
) and β
(
t
) denote the proportion of infected PWIDs and proportion of infected needles, respectively. Thus the absolute numbers of infected PWIDs and infected needles are n
π
(
t
) and m
β
(
t
). The spread of the disease amongst syringes can be described by the following differential equation:(2)
d
m
β
t
d
t
=
n
π
t
σ
¯
1
-
β
t
-
n
π
t
η
¯
β
t
-
n
1
-
π
t
ρ
¯
β
t
.Dividing by m,(3)
d
β
t
d
t
=
π
t
σ
1
-
β
t
-
π
t
η
β
t
-
1
-
π
t
ρ
β
t
=
π
t
σ
-
τ
β
t
-
1
-
π
t
ρ
β
t
,where(4)
τ
=
σ
+
η
=
γ
σ
¯
+
η
¯
,is the gallery ratio multiplied by the rate at which an infected PWID visits syringes multiplied by the sum of the probability that he or she leaves an uninfected syringe infected after use plus the probability that he or she leaves an infected syringe uninfected after use.The spread of the disease amongst PWIDs can be described by the differential equation:(5)
d
n
π
t
d
t
=
n
1
-
π
t
υ
β
t
-
μ
n
π
t
.Dividing by n,(6)
d
π
t
d
t
=
1
-
π
t
υ
β
t
-
μ
π
t
.So in summary the equations describing the deterministic HIV model are(7)
d
β
t
d
t
=
π
t
σ
-
τ
β
t
-
1
-
π
t
ρ
β
t
,
d
π
t
d
t
=
1
-
π
t
υ
β
t
-
μ
π
t
.Greenhalgh and Hay define the basic reproduction number for the modified Kaplan model to be(8)
R
0
D
=
υ
σ
ρ
μ
,where in Section  4.5 of their paper [5] they have shown in detail that it corresponds to the usual biological definition, that is, the expected number of secondary infected PWIDs (infected PWIDs who became infected from sharing a syringe with the original infected PWID) caused during his or her entire infectious period by a single newly infected PWID entering a disease-free population at equilibrium. They also point out that it is also the expected number of secondary infected needles caused by a single newly infected needle entering the disease-free population at equilibrium.Greenhalgh and Hay then show that the disease dies out ifR
0
D
<
1 or R
0
D
=
1 and τ
>
ρ. If R
0
D
>
1 there are two possible equilibria, one with no disease present and the other with disease present. Thus this value of R
0
D clearly satisfies the usual properties of the deterministic threshold value in epidemic models. There is a unique endemic equilibrium (9)
β
∗
=
σ
τ
1
-
ρ
μ
σ
υ
,
π
∗
=
σ
υ
-
ρ
μ
μ
τ
+
σ
υ
-
ρ
μ
.If R
0
D
>
1 then the unique endemic equilibrium is locally stable. If τ
>
ρ and R
0
D
>
1, then β
(
t
)
→
β
∗ and π
(
t
)
→
π
∗ as t
→
∞ provided that π
(
0
)
>
0 or β
(
0
)
>
0.Note that the average rate at which PWIDs leave the sharing, injecting population is around 0.25/year [5] so, on average, they each share for four years. P
1
+
P
3, the probability of HIV transmission to a susceptible PWID on making a single injection with an infected syringe, is quite small (e.g., one estimate is 0.01 [5]). On the other hand, PWIDs are injected every few days which is a much shorter timescale than the demographic epidemiological changes which take several years. Hence changes in the fraction of syringes infected will typically happen a lot faster than changes in the fraction of PWIDs infected.Thus over an intermediate timescale it may be reasonable to assume thatπ
(
t
) is approximately constant in the last equation of (7). Hence β
(
t
) will approach its equilibrium value from this equation(10)
π
t
σ
π
t
+
ρ
-
ρ
π
t
.By substituting (10) into the second equation in (7) we deduce that(11)
d
π
t
d
t
=
1
-
π
t
π
t
υ
σ
π
t
τ
+
ρ
-
π
t
ρ
-
μ
π
t
.A similar technique of reducing the dimensions of the model by assuming that the needle equations are at equilibrium is used in models of variable infectivity of spread of HIV amongst PWIDs discussed by Greenhalgh and Lewis [8, 11] and Corson et al. [23, 24].In this paper, we introduce environmental stochasticity into system (11) by replacing the parameter υ by υ
+
υ
¯
(
d
B
t
/
d
t
), where B
(
t
) is a Brownian motion and υ
¯
>
0 is the intensity of the noise which is associated with the potential rate of infection υ. It is therefore clear that the total number of new PWIDs infected during the small time interval [
t
,
t
+
d
t
) is normally distributed with mean,(12)
n
1
-
π
t
π
t
υ
σ
π
t
τ
+
ρ
-
π
t
ρ
d
t
, and variance,(13)
n
2
1
-
π
t
2
π
t
2
υ
¯
2
σ
2
π
t
τ
+
ρ
-
π
t
ρ
2
d
t
.Notice that both of this mean and variance tend to zero asd
t goes to zero which is a biologically desirable property. This is a standard technique of introducing random noise in stochastic modelling [16, 18, 25–28] and corresponds to some stochastic environmental factor acting on each individual in the population.To justify why simple white noise is appropriate for our model suppose that we consider a timescale on whichβ
(
t
) and π
(
t
) are approximately constant. We consider the changes in a small time interval [
t
,
t
+
T
0
) and divide it into a series of n
0 equal width subintervals t
,
t
+
T
,
t
+
T
,
t
+
2
T
,
…
,
[
t
+
(
n
0
-
1
)
T
,
t
+
n
0
T
), where n
0
T
=
T
0 and n
0 is very large. Then the numbers of new infections caused by a single susceptible PWID visiting one infected syringe during each of the subintervals t
,
t
+
T
,
t
+
T
,
t
+
2
T
,
…
,
[
t
+
(
n
0
-
1
)
T
,
t
+
n
0
T
) are identically distributed random variables, say, with common mean μ
0 and common variance σ
0
2. We assume that σ
0
2
<
∞. So by the Central Limit Theorem the total number of PWIDs who visit infected syringes and become infected in [
t
,
t
+
n
0
T
)
=
[
t
,
t
+
T
0
) is approximately normally distributed with mean n
0
μ
0 and variance n
0
σ
0
2. Moreover keeping T fixed and doubling T
0 doubles n
0 thus the mean and variance of the number of susceptible PWIDs who become infected from visiting an infected syringe in [
t
,
t
+
T
0
) are both proportional to T
0. Hence it is appropriate to consider simple white noise where the mean number of infections in [
t
,
t
+
d
t
) caused by a given susceptible PWID visiting a given infected syringe is υ
d
t (hence proportional to d
t) the same as in the deterministic model and the variance of this number is also proportional to d
t.As a result, we obtain the following SDE HIV model:(14)
d
π
t
=
1
-
π
t
π
t
υ
σ
π
t
τ
+
ρ
-
π
t
ρ
-
μ
π
t
d
t
+
1
-
π
t
π
t
υ
¯
σ
π
t
τ
+
ρ
-
π
t
ρ
d
B
t
.The reason why we chose to perturb the parameterυ corresponding to the total rate at which PWIDs visit syringes and potentially become infected is because as it multiplies the term π
(
1
-
π
) in (11) it is a key parameter in the transmission of HIV amongst PWIDs and we thought that this would be the most interesting and important parameter when analysing the effect that environmental noise would have on the spread of HIV.There are some environmental factors which can cause a perturbation inυ, for example, natural biological variation between people and between HIV viruses. These factors affect the probability P
1
+
P
3 of HIV transmission to a susceptible PWID. It is possible that environmental noise causes variation in other parameters too, but it would be quite complicated to include these as well. Analysis of the model with environmental stochasticity in υ provides theoretical insight into the behaviour of the model. A similar approach of introducing environmental stochasticity into only the disease transmission parameter was discussed in stochastic studies of epidemic models by Ding et al. [13], Gray et al. [21], Lu [25], Tornatore et al. [29], and others.For the rest of the paper, we shall focus on analysing the SDE HIV model (14). Throughout this paper, unless stated otherwise, we shall assume that the unit of time is one day.


## 3. Existence of Unique Nonnegative Solution

Before we begin to investigate the dynamical behaviour of the SDE HIV model (14), it is important for us to show whether this SDE has a unique global nonnegative solution. It is well-known that, in order for an SDE to have a unique global solution for any given initial value, the coefficients of the equation are generally required to satisfy the linear growth condition and the local Lipschitz conditions [30]. It is clear that our coefficients in (14) satisfy the linear growth condition and they are locally Lipschitz continuous. As a result there is a unique, nonexplosive solution to (14). The following theorem shows that the solution remains in (
0,1
) if it starts there.Theorem 1.
For any given initial valueπ
(
0
)
=
π
0
∈
(
0,1
), the SDE HIV model (14) has a unique global nonnegative solution π
(
t
)
∈
(
0,1
) for all t
≥
0 with probability one; namely,(15)
P
π
t
∈
0,1
,
∀
t
≥
0
=
1
.Proof.
For any given initial valueπ
0
∈
(
0,1
), there is a unique global solution π
(
t
) for t
≥
0. Let k
0
≥
0 be sufficiently large so that π
0 lies within the interval (
1
/
k
0
,
1
-
(
1
/
k
0
)
). Then for each integer k
≥
k
0, define the stopping time (16)
τ
k
=
inf
⁡
t
≥
0
:
π
t
∉
1
k
,
1
-
1
k
,where inf
⁡
∅
=
∞. It is easy to see that τ
k is increasing as k
→
∞. Let us also define τ
∞
=
l
i
m
k
→
∞
τ
k. To complete the proof, we need to show that τ
∞
=
∞ a.s. We will carry this proof out by contradiction. Let us therefore assume that the statement is false and thus there exists a pair of constants T
>
0 and ε
∈
(
0,1
) such that(17)
P
τ
∞
≤
T
>
ε
.Hence, there is an integer k
1
≥
k
0 such that(18)
P
τ
k
≤
T
>
ε
∀
k
≥
k
1
.Let us define a function V
:
(
0,1
)
→
R,(19)
V
x
=
1
x
+
1
1
-
x
.Now by Itô’s formula, we have that, for any t
∈
[
0
,
T
] and k
≥
k
1,(20)
E
V
π
t
∧
τ
k
=
V
π
0
+
E
∫
0
t
∧
τ
k
L
V
π
s
d
s
,where L
V
:
(
0,1
)
→
R is defined by(21)
L
V
x
=
-
1
-
x
υ
σ
x
x
τ
+
ρ
-
x
ρ
+
x
υ
σ
x
τ
+
ρ
-
x
ρ
1
-
x
+
μ
x
-
μ
x
1
-
x
2
+
1
-
x
2
υ
¯
2
σ
2
x
τ
+
ρ
-
x
ρ
2
x
+
x
2
υ
¯
2
σ
2
1
-
x
x
τ
+
ρ
-
x
ρ
2
.Furthermore, it is easy to see that(22)
L
V
x
≤
υ
σ
min
⁡
τ
,
ρ
1
-
x
+
μ
x
+
υ
¯
2
σ
2
min
⁡
τ
,
ρ
2
1
x
+
1
1
-
x
≤
C
V
x
,where (23)
C
=
υ
σ
min
⁡
τ
,
ρ
∨
μ
+
υ
¯
2
σ
2
min
⁡
τ
,
ρ
2
.Here a
∨
b denotes the maximum of a and b. By substituting this into (20), we have that for any t
∈
[
0
,
T
]
(24)
E
V
π
t
∧
τ
k
≤
V
π
0
+
C
∫
0
t
E
V
π
s
∧
τ
k
d
s
.Then by using the Gronwall inequality we have that (25)
E
V
π
t
∧
τ
k
≤
V
π
0
e
C
t
≤
V
π
0
e
C
T
.Let us set Ω
k
=
{
τ
k
≤
T
} for k
≥
k
1, and so, by (18), we have that P
(
Ω
k
)
≥
ε
. For every ω
∈
Ω
k
,
π
(
τ
k
,
ω
) equals either 1
/
k or 1
-
(
1
/
k
) and thus V
(
π
(
τ
k
,
ω
)
)
≥
k
. Consequently we have that(26)
V
π
0
ε
C
T
≥
E
1
Ω
k
ω
V
π
τ
k
,
ω
≥
k
P
Ω
k
≥
ε
k
.Letting k
→
∞, we have a contradiction, where ∞
>
V
π
0
ε
C
T
=
∞. Therefore, our assumption at the beginning must be false and thus we obtained our desired result that τ
∞
=
∞ a.s.In this section we have managed to show that there exists a unique nonnegative global solution for the SDE HIV model (14) which remains in (
0,1
).


## 4. Extinction

When studying the dynamical behaviour of a population system, it is important for us to consider the conditions required in order for the HIV amongst PWIDs to die out, in other words when the disease will become extinct. We will split this proof into two parts, each considering two different scenarios of the noise intensity, namely,υ
¯
. Before we begin the proof, let us recall the basic reproduction number for the deterministic model of Greenhalgh and Hay [5]:(27)
R
0
D
=
υ
σ
ρ
μ
,where all the parameters are defined as before.For the stochastic model we define the stochastic basic reproduction number(28)
R
0
S
=
R
0
D
-
υ
¯
2
σ
2
2
ρ
2
μ
.This is the deterministic basic reproduction number R
0
D corrected for the effect of stochastic noise and plays a role in the stochastic model with many similarities to R
0
D in the deterministic one.Theorem 2.
If the stochastic reproduction number(29)
R
0
S
=
R
0
D
-
υ
¯
2
σ
2
2
ρ
2
μ
<
1
,
υ
¯
2
≤
υ
ρ
σ
,then, for any given initial value π
(
0
)
=
π
0
∈
(
0,1
), the solution of (14) obeys(30)
l
i
m
t
→
∞
sup
⁡
1
t
log
⁡
π
t
≤
υ
σ
ρ
-
υ
¯
2
σ
2
2
ρ
2
-
μ
=
μ
R
0
S
-
1
<
0
a
.
s
. In other words, π
(
t
) will tend to zero exponentially a.s. Thus the fraction of population that is infected with HIV at time t will approach zero.Proof.
Let us define a functionV
(
x
)
=
log
⁡
x, where by Itô’s formula we have that(31)
log
⁡
π
t
=
log
⁡
π
0
+
∫
0
t
f
π
s
d
s
+
∫
0
t
1
-
π
s
υ
¯
σ
π
s
τ
+
ρ
-
π
s
ρ
d
B
s
.Here f
:
(
0,1
)
→
R is defined as(32)
f
x
=
1
-
x
υ
σ
x
τ
+
ρ
-
x
ρ
-
μ
-
1
-
x
2
υ
¯
2
σ
2
2
x
τ
+
ρ
-
x
ρ
2
(33)
=
υ
σ
φ
+
ρ
-
μ
-
υ
¯
2
σ
2
2
φ
+
ρ
2
,where φ
=
x
τ
/
(
1
-
x
)
. Moreover (34)
f
′
φ
=
-
υ
σ
φ
+
ρ
2
+
υ
¯
2
σ
2
φ
+
ρ
3
<
-
υ
σ
φ
+
ρ
2
+
υ
ρ
σ
φ
+
ρ
3
<
0
.Hence f
(
φ
) is a monotone decreasing function of φ for φ
>
0, and thus we must have that(35)
f
x
≤
f
x
φ
=
0
=
μ
R
0
S
-
1
<
0
, where R
0
S is the stochastic reproduction number defined in Theorem 2. As a result, (31) becomes(36)
log
⁡
π
t
≤
log
⁡
π
0
+
t
μ
R
0
S
-
1
+
∫
0
t
1
-
π
s
υ
¯
σ
π
s
τ
+
ρ
-
π
s
ρ
d
B
s
. This implies that(37)
l
i
m
t
→
∞
sup
⁡
1
t
log
⁡
π
t
≤
μ
R
0
S
-
1
+
l
i
m
t
→
∞
s
u
p
1
t
∫
0
t
1
-
π
s
υ
¯
σ
π
s
τ
+
ρ
-
π
s
ρ
d
B
s
.However, since (38)
0
≤
1
-
π
s
π
s
τ
+
ρ
-
π
s
ρ
≤
1
min
⁡
τ
,
ρ
,then, by the large number theorem of martingales (e.g., [28]), we have that(39)
l
i
m
t
→
∞
s
u
p
1
t
∫
0
t
1
-
π
s
υ
¯
σ
π
s
τ
+
ρ
-
π
s
ρ
d
B
s
=
0
a
.
s
. Hence, we have arrived at our desired result, where(40)
l
i
m
t
→
∞
s
u
p
1
t
log
⁡
π
t
≤
μ
R
0
S
-
1
<
0
a
.
s
. In other words, π
(
t
) tends to zero exponentially a.s.In Theorem2, we have focused on discussing the extinction conditions for our SDE HIV model (14); we have considered a partial case, where the noise intensity satisfies the condition υ
¯
2
≤
υ
ρ
/
σ. In order to get a better picture of the dynamical behaviour of our SDE HIV model (14), it is important for us to investigate what happens to the population system when υ
¯
2
>
υ
ρ
/
σ.Theorem 3.
If(41)
R
0
S
=
R
0
D
-
υ
¯
2
σ
2
2
ρ
2
μ
<
1
,
υ
¯
2
>
υ
ρ
σ
∨
υ
2
2
μ
,then, for any given initial value π
(
0
)
=
π
0
∈
(
0,1
), the solution of (14) obeys(42)
l
i
m
t
→
∞
sup
⁡
1
t
log
⁡
π
t
≤
υ
2
2
υ
¯
2
-
μ
<
0
a
.
s
.In other words, π
(
t
) will tend to zero exponentially a.s. Thus the fraction of the population that are infected with HIV at time t will become zero.Proof.
In order to simplify the computation, throughout this proof, we will be working with (33). It is easy to see that this function has a maximum turning point at(43)
φ
=
φ
^
=
υ
¯
2
σ
υ
-
ρ
.
Note that, by substituting (43) back into the expression φ
=
x
τ
/
(
1
-
x
), we could easily obtain the same result as we would if we decided to work with the alternative function(44)
x
=
x
^
=
υ
¯
2
σ
-
ρ
υ
υ
τ
-
ρ
+
υ
¯
2
σ
,where x
^
∈
(
0,1
). Note also that φ
^
>
0 by (41). Furthermore, by substituting the maximum turning point φ
^ given in (43) into (33), we have that f
(
x
)
|
φ
=
φ
^
=
υ
2
/
2
υ
¯
2
-
μ which is negative by condition (41). Therefore, arguing as before in Theorem 2, we have that(45)
log
⁡
π
t
≤
log
⁡
π
0
+
t
υ
2
2
υ
¯
2
-
μ
+
∫
0
t
1
-
π
s
υ
¯
σ
π
s
τ
+
ρ
-
π
s
ρ
d
B
s
, which similarly implies that(46)
l
i
m
t
→
∞
sup
⁡
1
t
log
⁡
π
t
≤
υ
2
2
υ
¯
2
-
μ
<
0
a
.
s
.
In other words,π
(
t
) will also tend to zero exponentially a.s. for υ
¯
2
>
υ
ρ
/
σ
∨
υ
2
/
2
μ and thus we have completed the proof.Note thatR
0
S
<
R
0
D, which implies that the condition for extinction is weaker in the stochastic case compared to the deterministic case. In addition, as υ
¯ increases, the stochastic reproduction number R
0
S will become smaller and thus it will be more likely for the HIV virus to die out for large noise intensity. As a result, this highlights the fact that environmental factors play an important role in the dynamical behaviour of HIV amongst PWIDs.There is a gap in our results forR
0
S
<
1. We have not shown what will happen if R
0
S
<
1 and(47)
υ
ρ
σ
<
υ
¯
2
<
υ
2
2
μ
,but we conjecture that in this case the disease will die out a.s. This was confirmed by simulation.


## 5. Persistence

Another very important aspect of the behaviour of a dynamical system is the conditions for persistence. In this section we will discuss the persistence conditions required for our SDE HIV model (14).Theorem 4.
If(48)
R
0
S
=
R
0
D
-
υ
¯
2
σ
2
2
ρ
2
μ
>
1
,then, for any given initial value π
(
0
)
=
π
0
∈
(
0,1
), the solution of (14) satisfies(49)
l
i
m
t
→
∞
sup
⁡
π
t
≥
η
a
.
s
,
(50)
l
i
m
t
→
∞
inf
⁡
π
t
≤
η
a
.
s
,where(51)
η
=
υ
σ
-
2
μ
ρ
+
υ
2
σ
2
-
2
μ
υ
¯
2
σ
2
2
μ
τ
+
υ
σ
-
2
μ
ρ
+
υ
2
σ
2
-
2
μ
υ
¯
2
σ
2
>
0
,which is the unique root in (
0,1
) of the function(52)
f
x
=
1
-
x
υ
σ
x
τ
+
ρ
-
x
ρ
-
μ
-
1
-
x
2
υ
¯
2
σ
2
2
x
τ
+
ρ
-
x
ρ
2
=
0
,defined in (32). In other words, the solution π
(
t
) will persist and oscillate around the level η infinitely often with probability one.Proof.
Let us recall the functionf
:
(
0,1
)
→
R defined in (33). Throughout this proof, we will be working with this function in order to simplify the computation.
By settingf
(
x
)
=
0, we obtain one positive and one negative root, where the positive root is(53)
φ
∗
=
1
2
μ
υ
σ
-
2
μ
ρ
2
+
4
μ
υ
σ
ρ
-
μ
ρ
2
-
0.5
υ
¯
2
σ
2
+
υ
σ
-
2
μ
ρ
>
φ
^
,where φ
^ is the maximum turning point of (33) defined in (43). For the purpose of consistency, we will now substitute (53) into the expression φ
∗
=
x
∗
τ
/
(
1
-
x
∗
) to get that(54)
x
∗
=
η
=
υ
σ
-
2
μ
ρ
+
υ
2
σ
2
-
2
μ
υ
¯
2
σ
2
2
μ
τ
+
υ
σ
-
2
μ
ρ
+
υ
2
σ
2
-
2
μ
υ
¯
2
σ
2
>
x
^
, where x
∗
∈
(
0,1
), and that x
^ is the equivalent maximum turning point of (32) defined in (44). Moreover, it is easy to see that (55)
f
0
=
υ
σ
ρ
-
μ
-
υ
¯
2
σ
2
2
ρ
2
>
0
,
f
1
=
-
μ
<
0
.As a result we have that(56)
f
x
>
0
i
s
s
t
r
i
c
t
l
y
i
n
c
r
e
a
s
i
n
g
o
n
x
∈
0,0
∨
x
^
,
(57)
f
x
>
0
i
s
s
t
r
i
c
t
l
y
d
e
c
r
e
a
s
i
n
g
o
n
x
∈
0
∨
x
^
,
x
∗
,
(58)
f
x
<
0
i
s
s
t
r
i
c
t
l
y
d
e
c
r
e
a
s
i
n
g
o
n
x
∈
x
∗
,
1
,where x
^ and x
∗ are defined as before. Let us now prove result (49) is true by contradiction. Assume that (49) is false and thus there must exist ε
∈
(
0,1
) small enough such that (59)
P
Ω
1
>
ε
,where Ω
1
=
{
ω
∈
Ω
:
l
i
m
t
→
∞
sup
⁡
π
(
t
)
≤
η
-
2
ε
}
. Hence, for every ω
∈
Ω
1, there is T
=
T
(
ω
)
>
0 such that(60)
π
t
,
ω
≤
η
-
ε
f
o
r
t
≥
T
ω
.
Clearly we can chooseε so small such that f
(
0
)
>
f
(
η
-
ε
). Therefore, from (56), (57), and (60), we have that f
(
π
(
t
,
ω
)
)
>
f
(
η
-
ε
) for t
≥
T
(
ω
). Let us now recall that, for t
≥
0, (61)
log
⁡
π
t
=
log
⁡
π
0
+
∫
0
t
f
π
s
d
s
+
∫
0
t
1
-
π
s
υ
¯
σ
π
s
τ
+
ρ
-
π
s
ρ
d
B
s
,and then arguing as before, by the large number theorem of martingales, there is Ω
2
⊂
Ω with P
(
Ω
2
)
=
1, such that, for every ω
∈
Ω
2, (62)
l
i
m
t
→
∞
1
t
∫
0
t
1
-
π
s
υ
¯
σ
π
s
τ
+
ρ
-
π
s
ρ
d
B
s
=
0
. Therefore by fixing any ω
∈
Ω
1
∩
Ω
2, then, for t
≥
T
(
ω
),(63)
log
⁡
π
t
,
ω
≥
log
⁡
π
0
+
∫
0
T
ω
f
π
s
,
ω
d
s
+
f
η
-
ε
t
-
T
ω
+
∫
0
t
1
-
π
s
υ
¯
σ
π
s
τ
+
ρ
-
π
s
ρ
d
B
s
,
ω
,which implies that(64)
l
i
m
t
→
∞
inf
⁡
1
t
log
⁡
π
t
,
ω
≥
f
η
-
ε
>
0
, and thus we have that l
i
m
t
→
∞
π
(
t
,
ω
)
=
∞. This is clearly a contradiction to (60). Thus, our assumption at the beginning must be wrong and therefore we obtained our desired result that (65)
l
i
m
t
→
∞
sup
⁡
π
t
≥
η
a
.
s
.Similarly, we will prove (50) by assuming again that it is false and thus there must exist δ
∈
(
0,1
) such that (66)
P
Ω
3
>
δ
,where Ω
3
=
{
ω
∈
Ω
:
l
i
m
t
→
∞
inf
⁡
π
(
t
)
≥
η
+
2
δ
}
. Hence, for every ω
∈
Ω
3, there is τ
=
τ
(
ω
)
>
0 such that(67)
π
t
,
ω
≥
η
+
δ
f
o
r
t
≥
τ
ω
.Thus, we have that f
(
π
(
t
,
ω
)
)
≤
f
(
η
+
δ
) for t
≥
τ
(
ω
). Let us now fix any ω
∈
Ω
2
∩
Ω
3; then similarly to before, we would get that, for t
≥
τ
(
ω
),(68)
log
⁡
π
t
,
ω
≤
log
⁡
π
0
+
∫
0
τ
ω
f
π
s
,
ω
d
s
+
f
η
+
δ
t
-
τ
ω
+
∫
0
t
1
-
π
s
υ
¯
σ
π
s
τ
+
ρ
-
π
s
ρ
d
B
s
,
ω
⟹
l
i
m
t
→
∞
sup
⁡
1
t
log
⁡
π
t
,
ω
≤
f
η
+
δ
<
0
,and thus (69)
⟹
lim
t
→
∞
⁡
π
t
,
ω
=
0
.
This is clearly a contradiction to (67) and thus we have completed our proof.In order to allow us to better understand the effect of the noise intensityυ
¯ on the dynamical behaviour of our SDE HIV model (14) and its connection to the corresponding deterministic model (11), we have the following proposition.Proposition 5.
Suppose thatR
0
S
>
1. Consider η as defined by (51) as a function of υ
¯ for(70)
0
<
υ
¯
<
2
ρ
υ
σ
-
μ
ρ
σ
=
υ
^
,and then η is strictly decreasing and(71)
lim
υ
¯
→
0
⁡
η
=
υ
σ
-
μ
ρ
υ
σ
-
μ
ρ
+
μ
τ
,which is the equilibrium state of the deterministic HIV model (11) and(72)
lim
υ
¯
→
υ
^
⁡
η
=
0
,
i
f
1
≤
R
0
D
≤
2
,
υ
σ
-
2
μ
ρ
υ
σ
-
2
μ
ρ
+
μ
τ
,
i
f
R
0
D
>
2
.In other words, the noise intensity υ
¯ lies between the deterministic equilibrium value for π
(
t
), namely, (
υ
σ
-
μ
ρ
)
/
(
υ
σ
-
μ
ρ
+
μ
τ
), and m
a
x
(
0
,
(
υ
σ
-
2
μ
ρ
)
/
(
υ
σ
-
2
μ
ρ
+
μ
τ
)
). Furthermore, if the noise intensity decreases to zero, then η will increase to the deterministic equilibrium value. If R
0
D is large, then η will be close to but beneath the deterministic equilibrium value for π
(
t
)
.Proof.
Let us recall that(73)
η
=
υ
σ
-
2
μ
ρ
+
υ
2
σ
2
-
2
μ
υ
¯
2
σ
2
2
μ
τ
+
υ
σ
-
2
μ
ρ
+
υ
2
σ
2
-
2
μ
υ
¯
2
σ
2
.
Then,(74)
d
η
d
υ
¯
=
-
4
μ
2
υ
¯
σ
2
τ
υ
2
σ
2
-
2
μ
υ
¯
2
σ
2
1
/
2
2
μ
τ
+
υ
σ
-
2
μ
ρ
+
υ
2
σ
2
-
2
μ
υ
¯
2
σ
2
1
/
2
2
.Clearly, d
η
/
d
υ
¯
<
0 since σ
>
0 and thus η is strictly deceasing as υ
¯ increases. By letting υ
¯ tend to zero in the function for η defined above, we have the desired result given in (71). Moreover, as υ
¯
→
υ
^, we have that(75)
lim
υ
¯
→
υ
^
⁡
η
=
υ
σ
-
2
μ
ρ
+
υ
2
σ
2
-
2
μ
υ
^
2
σ
2
2
μ
τ
+
υ
σ
-
2
μ
ρ
+
υ
2
σ
2
-
2
μ
υ
^
2
σ
2
. The numerator of the above expression is equal to(76)
υ
σ
-
2
μ
ρ
+
υ
σ
-
2
μ
ρ
.
As a result, if1
≤
R
0
D
≤
2, then it is obvious that l
i
m
υ
¯
→
υ
^
η
=
0. On the other hand, if R
0
D
>
2, then l
i
m
υ
¯
→
υ
^
η
=
(
υ
σ
-
2
μ
ρ
)
/
(
υ
σ
-
2
μ
ρ
+
μ
τ
). We have completed the proof.


## 6. Stationary Distribution

In this section, we will use the well-known Khasminskii theorem [31] to prove that there exists a stationary distribution for our stochastic HIV model (14). Before we begin, let us recall the conditions for the existence of a stationary distribution mentioned in [31].Lemma 6.
The SDE HIV model (14) has a unique stationary distribution if there is a strictly proper subinterval (
a
,
b
) of (
0,1
) such that E
(
τ
)
<
∞ for all π
0
∈
(
0
,
a
)
∪
(
b
,
1
), where(77)
τ
=
inf
⁡
t
≥
0
:
π
t
∈
a
,
b
,
sup
π
0
∈
a
¯
,
b
¯
⁡
E
τ
<
∞
f
o
r
e
v
e
r
y
i
n
t
e
r
v
a
l
a
¯
,
b
¯
⊂
0,1
.Note that, in the original Khasminskii theorem, there is an additional condition which states that the square of the diffusion coefficient of the SDE HIV model (14), namely, (78)
1
-
π
t
π
t
υ
¯
σ
π
t
τ
+
ρ
-
π
t
ρ
2
,is bounded away from zero for π
(
t
)
∈
(
a
,
b
). However, recalling from the proof of Theorem 1, we have already shown that the denominator, (
π
(
t
)
τ
+
ρ
-
π
(
t
)
ρ
), is bounded away from zero (it is at least min
⁡
(
τ
,
ρ
)). Thus it is therefore clear that this condition holds for our model.Theorem 7.
IfR
0
S
>
1, then the SDE HIV model (14) has a unique stationary distribution.Proof.
Let us fix any0
<
a
<
η
<
b
<
1. From conditions (56)–(58) in the proof for Theorem 4 we can see that(79)
f
x
≥
f
0
∧
f
a
>
0
i
f
0
<
x
≤
a
,
f
x
≤
f
b
<
0
i
f
b
≤
x
<
1
.Let us now define the stopping time τ as we did in Lemma 6. Recall that (80)
log
⁡
π
t
=
log
⁡
π
0
+
∫
0
t
f
π
s
d
s
+
∫
0
t
1
-
π
s
υ
¯
σ
π
s
τ
+
ρ
-
π
s
ρ
d
B
s
,and then, by using (79), we have that, for all t
≥
0 and for any π
0
∈
(
0
,
a
),(81)
log
⁡
a
≥
E
log
⁡
π
t
∧
τ
≥
log
⁡
π
0
+
f
0
∧
f
a
E
t
∧
τ
, and then(82)
log
⁡
a
π
0
≥
f
0
∧
f
a
E
t
∧
τ
. By letting t
→
∞, we have that for all π
0
∈
(
0
,
a
)
(83)
E
τ
≤
log
⁡
a
/
π
0
f
0
∧
f
a
.Similarly, for any π
0
∈
(
b
,
1
), we have that(84)
log
⁡
b
≤
E
log
⁡
π
t
∧
τ
≤
log
⁡
π
0
-
f
b
E
t
∧
τ
,
∀
t
≥
0
,and then(85)
log
⁡
b
π
0
≤
-
f
b
E
t
∧
τ
. By letting t
→
∞, we have that(86)
E
τ
≤
log
⁡
π
0
/
b
f
b
≤
log
⁡
1
/
b
f
b
∀
π
0
∈
b
,
1
.Clearly, the conditions required for existence of a unique stationary distribution mentioned in Lemma 6 are satisfied by (83) and (86) and thus we have completed our proof and our SDE HIV model (14) has a unique stationary distribution.


## 7. Simulations

In this section we will support our analytical results using numerical simulations produced inR. Throughout this section, various simulations are produced using realistic parameter values but our main objective is to verify the analytic results. Before we begin, let us make the same assumptions as in [5]. Without loss of generality, let us take p
=
0 and assume that all PWIDs visit shooting galleries at the same rate whether or not they are infected and thus λ
1
=
λ
2. In addition, we take ϕ
1
=
θ
1
=
0 as these probabilities are very small.


### 7.1. Simulations on Extinction

In this section, we will focus on looking at the numerical simulations produced which support the analytical results given in Theorems2 and 3.Example 8 (R
0
S
<
1
,
υ
¯
2
≤
υ
ρ
/
σ).
Let us choose realistic parameter valuesμ
=
0.258/year = 7.06849 × 10−4/day [7], λ
1
=
λ
2
=
0.143, α
=
(
P
1
+
P
3
)
=
0.01, θ
=
(
P
1
+
P
2
)
=
0.25, γ
=
1 (based on [5]), and ξ
=
0.6 [10]; then, from (1) and (4), we have that σ
=
0.0572, τ
=
0.143, ρ
=
0.1001, and υ
=
0.00143. Then by choosing υ
¯
=
0.046, we have that (87)
υ
¯
2
=
0.002116
<
υ
ρ
σ
=
0.0025025
,where R
0
S
=
0.66729
<
1 while R
0
D
=
1.156. Therefore, by Theorem 2, we would expect the solution π
(
t
) to reach zero with probability one.The computer simulation produced inR using the Euler-Maruyama method ([21, 28]) with the above parameter values is given in Figure 1, which clearly illustrates that π
(
t
) hits zero in finite time a.s. The numerical simulations were repeated numerous times with different initial value of π
0
∈
(
0,1
) and similar results were obtained each time.Figure 1
Computer simulations of the pathπ
(
t
) for the SDE HIV model (14) and its corresponding deterministic HIV model (11) with step size Δ
=
0.01 with parameter values given in Example 8 with initial value π
(
0
)
=
0.5.Example 9 (R
0
S
<
1
,
υ
¯
2
>
υ
ρ
/
σ
∨
υ
2
/
2
μ).
By using the same parameter values as in Example8 but choosing υ
¯ to be 0.07 and thus υ
¯
2
=
0.0049, we have that (88)
υ
¯
2
>
υ
ρ
σ
∨
υ
2
2
μ
,where R
0
S
=
0.02425249
<
1 while R
0
D
=
1.156035. As a result, by Theorem 3, we could conclude that, for any initial value π
(
0
)
=
π
0
∈
(
0,1
), the solution π
(
t
) obeys (89)
l
i
m
t
→
∞
sup
⁡
1
t
log
⁡
π
t
≤
-
0.000498186
<
0
a
.
s
.
Clearly Figure2 supports this result by showing that the solution π
(
t
) reaches zero at finite time. Again, the numerical simulations were repeated several times with different initial values and the same results were concluded.Figure 2
Computer simulations of the pathπ
(
t
) for the SDE HIV model (14) and its corresponding deterministic HIV model (11) with step size Δ
=
0.01 using parameter values given in Example 9 with initial value π
(
0
)
=
0.5.


### 7.2.  Simulation on Persistence

We will now move on to the numerical simulations for results given in Theorem4 and Proposition 5.Example 10 (R
0
S
>
1).
Let us use the same parameter values as in Example8 but changing μ to 0.125/year and thus 3.42466
×
1
0
-
4/day [4]. Let us define υ
¯
=
0.05 and thus R
0
D
=
2.386057 and R
0
S
=
1.1942
>
1. Therefore by Theorem 4, for any given initial value π
0
=
π
(
0
)
∈
(
0,1
), the solution π
(
t
) for the SDE HIV model (14) should obey (90)
l
i
m
t
→
∞
inf
⁡
π
t
≤
η
=
0.3206092
≤
l
i
m
t
→
∞
sup
⁡
π
t
a
.
s
.
Figure3 clearly supports our analytical results given in Theorem 4 by showing the solution path of π
(
t
) oscillates around the level η in finite time. Again the numerical simulations were repeated and the same conclusion can be drawn each time.Figure 3
Computer simulations of the pathπ
(
t
) for the SDE HIV model (14) and its corresponding deterministic HIV model (11) with step size Δ
=
0.01 using parameter values given in Example 10 with initial value π
(
0
)
=
0.3.In order to further illustrate the effect of the noise intensityυ
¯ has on the solution, in the next example we will keep all the parameter values the same as in Example 10 but reducing the noise intensity.Example 11.
By keeping the parameter values the same as in Example10 and reducing υ
¯ to 0.02, we have that R
0
D
=
2.38605
,
R
0
S
=
2.195363
>
1, and η
=
0.4770654
. By Theorem 4 and Proposition 5, we would expect the solution π
(
t
) to persist and oscillate around the level η. Furthermore by Proposition 5, as υ
¯
→
0, we would expect η to tend towards the deterministic equilibrium value for the corresponding deterministic model given by (11); namely, (
υ
σ
-
μ
ρ
)
/
(
υ
σ
-
μ
ρ
+
μ
τ
)
=
0.4924476.
From Figure4, we can clearly see that the solution path π
(
t
) does indeed oscillate about the level η. Moreover, by comparing Figures 3 and 4, we can also see that as we reduce the noise intensity from 0.05 to 0.02, the level η does indeed tend towards the deterministic equilibrium value as expected.Figure 4
Computer simulations of the pathπ
(
t
) for the SDE HIV model (14) and its corresponding deterministic HIV model (11) with step size Δ
=
0.01 using parameter values given in Example 11 with initial value π
(
0
)
=
0.3.In the next example we will use histograms to see how the solution of the SDE HIV model oscillates around the levelη as we vary the noise intensity υ
¯.Example 12.
Let us use the same parameter values as in Example10 and choose υ
¯ to be 0.05, 0.04, 0.03, 0.005, and 0.001. We then let the simulations run for 1 million iterations but disregarding the first 800,000 iterations in order to allow π
(
t
) to reach its recurrent level.
From Figure5, we can see from the histograms that, for larger υ
¯, the distribution of the solution is more skewed, while, for smaller υ
¯, the distribution is more normally distributed about the level η. This is further confirmed by the sample skewness coefficients, namely, 1.000705, 0.4264459, -
0.3415049, 0.2185976, and 0.003324568 corresponding to υ
¯
=
0.05, 0.04, 0.03, 0.005, and 0.001, respectively.
Furthermore, we have used the quantile-quantile plot (QQ plot) to further illustrate that, for the smaller values ofυ
¯, these data are not far from being normally distributed. The result is shown in Figure 6.Figure 5
Histogram for the solution pathπ
(
t
) for the SDE HIV model (14) with step size Δ
=
0.01 using parameter values given in Example 12 with initial value π
(
0
)
=
0.3 and υ
¯
=
0.05,0.04,0.03,0.005, and 0.001.Figure 6
QQ plot for the solution pathπ
(
t
) for the SDE HIV model (14) corresponding to the histograms shown in Figure 5 for υ
¯
=
0.03,0.005, and 0.001.


## 7.1. Simulations on Extinction

In this section, we will focus on looking at the numerical simulations produced which support the analytical results given in Theorems2 and 3.Example 8 (R
0
S
<
1
,
υ
¯
2
≤
υ
ρ
/
σ).
Let us choose realistic parameter valuesμ
=
0.258/year = 7.06849 × 10−4/day [7], λ
1
=
λ
2
=
0.143, α
=
(
P
1
+
P
3
)
=
0.01, θ
=
(
P
1
+
P
2
)
=
0.25, γ
=
1 (based on [5]), and ξ
=
0.6 [10]; then, from (1) and (4), we have that σ
=
0.0572, τ
=
0.143, ρ
=
0.1001, and υ
=
0.00143. Then by choosing υ
¯
=
0.046, we have that (87)
υ
¯
2
=
0.002116
<
υ
ρ
σ
=
0.0025025
,where R
0
S
=
0.66729
<
1 while R
0
D
=
1.156. Therefore, by Theorem 2, we would expect the solution π
(
t
) to reach zero with probability one.The computer simulation produced inR using the Euler-Maruyama method ([21, 28]) with the above parameter values is given in Figure 1, which clearly illustrates that π
(
t
) hits zero in finite time a.s. The numerical simulations were repeated numerous times with different initial value of π
0
∈
(
0,1
) and similar results were obtained each time.Figure 1
Computer simulations of the pathπ
(
t
) for the SDE HIV model (14) and its corresponding deterministic HIV model (11) with step size Δ
=
0.01 with parameter values given in Example 8 with initial value π
(
0
)
=
0.5.Example 9 (R
0
S
<
1
,
υ
¯
2
>
υ
ρ
/
σ
∨
υ
2
/
2
μ).
By using the same parameter values as in Example8 but choosing υ
¯ to be 0.07 and thus υ
¯
2
=
0.0049, we have that (88)
υ
¯
2
>
υ
ρ
σ
∨
υ
2
2
μ
,where R
0
S
=
0.02425249
<
1 while R
0
D
=
1.156035. As a result, by Theorem 3, we could conclude that, for any initial value π
(
0
)
=
π
0
∈
(
0,1
), the solution π
(
t
) obeys (89)
l
i
m
t
→
∞
sup
⁡
1
t
log
⁡
π
t
≤
-
0.000498186
<
0
a
.
s
.
Clearly Figure2 supports this result by showing that the solution π
(
t
) reaches zero at finite time. Again, the numerical simulations were repeated several times with different initial values and the same results were concluded.Figure 2
Computer simulations of the pathπ
(
t
) for the SDE HIV model (14) and its corresponding deterministic HIV model (11) with step size Δ
=
0.01 using parameter values given in Example 9 with initial value π
(
0
)
=
0.5.


## 7.2.  Simulation on Persistence

We will now move on to the numerical simulations for results given in Theorem4 and Proposition 5.Example 10 (R
0
S
>
1).
Let us use the same parameter values as in Example8 but changing μ to 0.125/year and thus 3.42466
×
1
0
-
4/day [4]. Let us define υ
¯
=
0.05 and thus R
0
D
=
2.386057 and R
0
S
=
1.1942
>
1. Therefore by Theorem 4, for any given initial value π
0
=
π
(
0
)
∈
(
0,1
), the solution π
(
t
) for the SDE HIV model (14) should obey (90)
l
i
m
t
→
∞
inf
⁡
π
t
≤
η
=
0.3206092
≤
l
i
m
t
→
∞
sup
⁡
π
t
a
.
s
.
Figure3 clearly supports our analytical results given in Theorem 4 by showing the solution path of π
(
t
) oscillates around the level η in finite time. Again the numerical simulations were repeated and the same conclusion can be drawn each time.Figure 3
Computer simulations of the pathπ
(
t
) for the SDE HIV model (14) and its corresponding deterministic HIV model (11) with step size Δ
=
0.01 using parameter values given in Example 10 with initial value π
(
0
)
=
0.3.In order to further illustrate the effect of the noise intensityυ
¯ has on the solution, in the next example we will keep all the parameter values the same as in Example 10 but reducing the noise intensity.Example 11.
By keeping the parameter values the same as in Example10 and reducing υ
¯ to 0.02, we have that R
0
D
=
2.38605
,
R
0
S
=
2.195363
>
1, and η
=
0.4770654
. By Theorem 4 and Proposition 5, we would expect the solution π
(
t
) to persist and oscillate around the level η. Furthermore by Proposition 5, as υ
¯
→
0, we would expect η to tend towards the deterministic equilibrium value for the corresponding deterministic model given by (11); namely, (
υ
σ
-
μ
ρ
)
/
(
υ
σ
-
μ
ρ
+
μ
τ
)
=
0.4924476.
From Figure4, we can clearly see that the solution path π
(
t
) does indeed oscillate about the level η. Moreover, by comparing Figures 3 and 4, we can also see that as we reduce the noise intensity from 0.05 to 0.02, the level η does indeed tend towards the deterministic equilibrium value as expected.Figure 4
Computer simulations of the pathπ
(
t
) for the SDE HIV model (14) and its corresponding deterministic HIV model (11) with step size Δ
=
0.01 using parameter values given in Example 11 with initial value π
(
0
)
=
0.3.In the next example we will use histograms to see how the solution of the SDE HIV model oscillates around the levelη as we vary the noise intensity υ
¯.Example 12.
Let us use the same parameter values as in Example10 and choose υ
¯ to be 0.05, 0.04, 0.03, 0.005, and 0.001. We then let the simulations run for 1 million iterations but disregarding the first 800,000 iterations in order to allow π
(
t
) to reach its recurrent level.
From Figure5, we can see from the histograms that, for larger υ
¯, the distribution of the solution is more skewed, while, for smaller υ
¯, the distribution is more normally distributed about the level η. This is further confirmed by the sample skewness coefficients, namely, 1.000705, 0.4264459, -
0.3415049, 0.2185976, and 0.003324568 corresponding to υ
¯
=
0.05, 0.04, 0.03, 0.005, and 0.001, respectively.
Furthermore, we have used the quantile-quantile plot (QQ plot) to further illustrate that, for the smaller values ofυ
¯, these data are not far from being normally distributed. The result is shown in Figure 6.Figure 5
Histogram for the solution pathπ
(
t
) for the SDE HIV model (14) with step size Δ
=
0.01 using parameter values given in Example 12 with initial value π
(
0
)
=
0.3 and υ
¯
=
0.05,0.04,0.03,0.005, and 0.001.Figure 6
QQ plot for the solution pathπ
(
t
) for the SDE HIV model (14) corresponding to the histograms shown in Figure 5 for υ
¯
=
0.03,0.005, and 0.001.


## 8. Conclusion

In this paper we have introduced environmental stochasticity into the extended Kaplan model for the spread of HIV amongst PWIDs constructed by Greenhalgh and Hay [5]. Inspired by the work done on introducing stochasticity by parameter perturbation into the SIS epidemic model in [21], we explored the properties for the resulting stochastic HIV model by first proving that there exists a unique nonnegative solution π
(
t
) for any given initial value π
0
∈
(
0,1
). Furthermore, we have constructed the basic reproduction number for the stochastic model, namely, R
0
S, and the conditions required for extinction and persistence for our solution π
(
t
)
. In general, if R
0
S
<
1, the solution will almost surely go extinct as shown in Theorems 2 and 3. There is a gap in our results if R
0
S and υ
ρ
/
σ
<
υ
¯
2
<
υ
2
/
2
μ but here we conjecture that disease will always die out. This conjecture was supported by simulation. On the other hand, the solution will almost surely persist and oscillate around the level η if R
0
S
>
1 as shown in Theorem 4. Most importantly, we have shown that, by altering the noise intensity υ
¯, it will affect the dynamical behaviour of our system.By using the well-known Khasminskii theorem, we have shown that the SDE HIV model has a unique stationary distribution. Lastly, numerical simulations using realistic parameter values are constructed to support our analytical results.Note thatR
0
S has a natural interpretation as follows: if we consider introducing a single newly infected individual into the disease-free equilibrium (DFE) and consider the number of secondary cases that he or she produces, then near the DFE equation (14) becomes(91)
d
π
t
=
π
υ
σ
ρ
-
μ
d
t
+
υ
¯
σ
ρ
π
d
B
,with solution(92)
π
t
=
π
0
exp
⁡
υ
σ
ρ
-
μ
-
1
2
υ
¯
2
σ
2
ρ
2
t
+
υ
¯
σ
ρ
B
t
.Alsol
i
m
t
→
∞
B
t
/
t
=
0 a.s. Hence we expect that if(93)
R
0
S
=
υ
σ
μ
ρ
-
υ
¯
2
σ
2
2
ρ
2
μ
<
1
,then the disease dies out whereas if R
0
S
>
1, the disease takes off. Thus this is a natural biological interpretation of the stochastic basic reproduction number R
0
S.Deterministic models have in the past proved very useful in describing the spread of HIV amongst PWIDs but they have their faults. The real world is stochastic and in general stochastic models are more realistic than deterministic ones. Recall that(94)
R
0
S
=
R
0
D
-
υ
¯
2
σ
2
2
ρ
2
μ
,where R
0
D represents the basic reproduction number in the deterministic model. So in the deterministic model R
0
D is the expected number of secondary cases caused by a single newly infected PWID entering a population consisting entirely of susceptible PWIDs and uninfected needles. The second term in (94) is an adjustment factor for the stochastic model.In the deterministic model we have a straightforward scenario where if the basic reproduction numberR
0
D
≤
1, then it is known that the disease will die out, whereas if R
0
D
>
1, then the disease will persist. The results in this paper show that, in the stochastic model, if R
0
S
<
1, then the disease dies out (almost surely), whereas if R
0
S
>
1, then the disease ultimately persists and oscillates about a nonzero level. These theoretical results are confirmed by numerical simulations. Moreover the argument above shows that if a single newly infected PWID enters the DFE, then we expect the disease to die out if R
0
S
<
1 and take off if R
0
S
>
1.These findings provide new insights into the spread of HIV amongst PWIDs. Because the stochastic basic reproduction numberR
0
S is less than the deterministic one, it is possible for the noise to drive the disease to extinction, that is, if R
0
D
>
1, so that in the deterministic model the disease will persist; then if the stochastic noise is large enough, in the stochastic model the disease will die out. This has important implications for control strategies. Deterministic models have often been used to predict control strategies, for example, the fraction of PWIDs who must clean their needles after use, the effects of HIV testing, or the amount that PWIDs need to decrease their syringe sharing rates in order to reduce R
0
D beneath one and eliminate disease. Examples of this applied to HIV amongst PWIDs include Greenhalgh and Lewis [32], Lewis [33] as well as Lewis and Greenhalgh [34]. Examples applied to hepatitis C virus (HCV) control include Corson [35] and Corson et al. [23].The analytical and numerical results of this paper provide new insight into this. If there is significant stochastic noise in the system, then these estimates will be overestimated; that is, a smaller fraction of PWIDs cleaning their needles or a smaller reduction in PWID syringe sharing rates will still be sufficient for elimination of disease transmission.



---

*Source: 6757928-2016-03-09.xml*