Computational and Mathematical Methods in Medicine

Volume 2016, Article ID 6757928, 14 pages

http://dx.doi.org/10.1155/2016/6757928

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

Department of Mathematics and Statistics, University of Strathclyde, Glasgow G1 1XH, UK

Received 8 October 2015; Revised 7 December 2015; Accepted 22 December 2015

Academic Editor: Chuangyin Dang

Copyright © 2016 Yanfeng Liang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.

#### 1. Introduction

HIV (human immunodeficiency virus) is a deadly and infectious* Lentivirus* 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 time , 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 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 (, , , ) be a complete probability space with filtration satisfying the usual conditions (i.e., it is increasing and right continuous while contains all -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: : shooting gallery visiting rate for susceptible PWIDs and the PWIDs who are infected but do not know they are infected. : shooting gallery visiting rate for infected PWIDs who know that they are infected. : probability that the needle is flushed and the PWID is infected. : probability that the needle is flushed and the PWID remains uninfected. : probability that the PWID becomes infected without the needle being flushed. : probability that the PWID remains uninfected and the needle is not flushed. : probability that an infected PWID leaves uninfected a syringe that was initially uninfected. : 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 and represents the PWID population and represents the number of shooting galleries or syringes that each PWID visits at random. : 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 .

Define the following new composite parameters:In the expression for the factor represents the probability that an initially uninfected syringe is left infected and not cleaned by an infected PWID. The term represents the average rate at which an infected PWID visits syringes. Hence , where 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.

represents the rate at which a susceptible PWID visits syringes and 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 and denote the proportion of infected PWIDs and proportion of infected needles, respectively. Thus the absolute numbers of infected PWIDs and infected needles are and . The spread of the disease amongst syringes can be described by the following differential equation:Dividing by ,whereis 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:Dividing by ,

So in summary the equations describing the deterministic HIV model are

Greenhalgh and Hay define the basic reproduction number for the modified Kaplan model to bewhere 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 if or and . If there are two possible equilibria, one with no disease present and the other with disease present. Thus this value of clearly satisfies the usual properties of the deterministic threshold value in epidemic models. There is a unique endemic equilibrium If then the unique endemic equilibrium is locally stable. If and , then and as provided that or .

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. , 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 is approximately constant in the last equation of (7). Hence will approach its equilibrium value from this equationBy substituting (10) into the second equation in (7) we deduce thatA 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 , where is a Brownian motion and 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 is normally distributed with mean, and variance,

Notice that both of this mean and variance tend to zero as 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 and are approximately constant. We consider the changes in a small time interval and divide it into a series of equal width subintervals , where and is very large. Then the numbers of new infections caused by a single susceptible PWID visiting one infected syringe during each of the subintervals are identically distributed random variables, say, with common mean and common variance . We assume that . So by the Central Limit Theorem the total number of PWIDs who visit infected syringes and become infected in is approximately normally distributed with mean and variance . Moreover keeping fixed and doubling doubles thus the mean and variance of the number of susceptible PWIDs who become infected from visiting an infected syringe in are both proportional to . Hence it is appropriate to consider simple white noise where the mean number of infections in caused by a given susceptible PWID visiting a given infected syringe is (hence proportional to ) the same as in the deterministic model and the variance of this number is also proportional to .

As a result, we obtain the following SDE HIV model:

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 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 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 if it starts there.

Theorem 1. *For any given initial value , the SDE HIV model (14) has a unique global nonnegative solution for all with probability one; namely,*

*Proof. *For any given initial value , there is a unique global solution for . Let be sufficiently large so that lies within the interval . Then for each integer , define the stopping time where . It is easy to see that is increasing as . Let us also define . 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 and such thatHence, there is an integer such thatLet us define a function ,Now by Itô’s formula, we have that, for any and ,where is defined byFurthermore, it is easy to see thatwhere Here denotes the maximum of and . By substituting this into (20), we have that for any Then by using the Gronwall inequality we have that Let us set for , and so, by (18), we have that For every equals either or and thus Consequently we have thatLetting , we have a contradiction, where . 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 .

#### 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]:where all the parameters are defined as before.

For the stochastic model we define the stochastic basic reproduction numberThis is the deterministic basic reproduction number corrected for the effect of stochastic noise and plays a role in the stochastic model with many similarities to in the deterministic one.

Theorem 2. *If the stochastic reproduction numberthen, for any given initial value , the solution of (14) obeys In other words, will tend to zero exponentially a.s. Thus the fraction of population that is infected with HIV at time will approach zero.*

*Proof. *Let us define a function , where by Itô’s formula we have thatHere is defined aswhere Moreover Hence is a monotone decreasing function of for , and thus we must have that where is the stochastic reproduction number defined in Theorem 2. As a result, (31) becomes This implies thatHowever, since then, by the large number theorem of martingales (e.g., [28]), we have that Hence, we have arrived at our desired result, where In other words, tends to zero exponentially a.s.

In Theorem 2, 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 . 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 .

Theorem 3. *Ifthen, for any given initial value , the solution of (14) obeysIn other words, will tend to zero exponentially a.s. Thus the fraction of the population that are infected with HIV at time 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 atNote that, by substituting (43) back into the expression , we could easily obtain the same result as we would if we decided to work with the alternative functionwhere . Note also that by (41). Furthermore, by substituting the maximum turning point given in (43) into (33), we have that which is negative by condition (41). Therefore, arguing as before in Theorem 2, we have that which similarly implies thatIn other words, will also tend to zero exponentially a.s. for and thus we have completed the proof.

Note that , 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 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 for . We have not shown what will happen if andbut 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. *Ifthen, for any given initial value , the solution of (14) satisfieswherewhich is the unique root in of the functiondefined in (32). In other words, the solution will persist and oscillate around the level infinitely often with probability one.*

*Proof. *Let us recall the function defined in (33). Throughout this proof, we will be working with this function in order to simplify the computation.

By setting , we obtain one positive and one negative root, where the positive root iswhere is the maximum turning point of (33) defined in (43). For the purpose of consistency, we will now substitute (53) into the expression to get that where , and that is the equivalent maximum turning point of (32) defined in (44). Moreover, it is easy to see that As a result we have thatwhere and are defined as before. Let us now prove result (49) is true by contradiction. Assume that (49) is false and thus there must exist small enough such that where Hence, for every , there is such thatClearly we can choose so small such that . Therefore, from (56), (57), and (60), we have that for . Let us now recall that, for , and then arguing as before, by the large number theorem of martingales, there is with , such that, for every , Therefore by fixing any , then, for ,which implies that and thus we have that . This is clearly a contradiction to (60). Thus, our assumption at the beginning must be wrong and therefore we obtained our desired result that Similarly, we will prove (50) by assuming again that it is false and thus there must exist such that where Hence, for every , there is such thatThus, we have that for . Let us now fix any ; then similarly to before, we would get that, for ,and thus 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 that . Consider as defined by (51) as a function of forand then is strictly decreasing andwhich is the equilibrium state of the deterministic HIV model (11) andIn other words, the noise intensity lies between the deterministic equilibrium value for , namely, , and . Furthermore, if the noise intensity decreases to zero, then will increase to the deterministic equilibrium value. If is large, then will be close to but beneath the deterministic equilibrium value for *

*Proof. *Let us recall thatThen,Clearly, since 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 The numerator of the above expression is equal toAs a result, if , then it is obvious that . On the other hand, if , then . 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 of such that for all , where*

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, is bounded away from zero for . However, recalling from the proof of Theorem 1, we have already shown that the denominator, , is bounded away from zero (it is at least ). Thus it is therefore clear that this condition holds for our model.

Theorem 7. *If , then the SDE HIV model (14) has a unique stationary distribution.*

*Proof. *Let us fix any . From conditions (56)–(58) in the proof for Theorem 4 we can see thatLet us now define the stopping time as we did in Lemma 6. Recall that and then, by using (79), we have that, for all and for any , and then By letting , we have that for all Similarly, for any , we have thatand then By letting , we have thatClearly, 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 in . 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 and assume that all PWIDs visit shooting galleries at the same rate whether or not they are infected and thus . In addition, we take 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 Theorems 2 and 3.

*Example 8 (). *Let us choose realistic parameter values /year = 7.06849 10^{−4}/day [7], , , , (based on [5]), and [10]; then, from (1) and (4), we have that , , , and . Then by choosing , we have that where while . Therefore, by Theorem 2, we would expect the solution to reach zero with probability one.

The computer simulation produced in using the Euler-Maruyama method ([21, 28]) with the above parameter values is given in Figure 1, which clearly illustrates that hits zero in finite time a.s. The numerical simulations were repeated numerous times with different initial value of and similar results were obtained each time.