Table of Contents Author Guidelines Submit a Manuscript
Computational and Mathematical Methods in Medicine
Volume 2016, Article ID 6757928, 14 pages
http://dx.doi.org/10.1155/2016/6757928
Research Article

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 [611]. 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 [1315]. 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, 2528] 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.

Figure 1: Computer simulations of the path for the SDE HIV model (14) and its corresponding deterministic HIV model (11) with step size with parameter values given in Example 8 with initial value .

Example 9 (). By using the same parameter values as in Example 8 but choosing to be 0.07 and thus , we have that where while . As a result, by Theorem 3, we could conclude that, for any initial value , the solution obeys Clearly Figure 2 supports this result by showing that the solution 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 for the SDE HIV model (14) and its corresponding deterministic HIV model (11) with step size using parameter values given in Example 9 with initial value .
7.2. Simulation on Persistence

We will now move on to the numerical simulations for results given in Theorem 4 and Proposition 5.

Example 10 (). Let us use the same parameter values as in Example 8 but changing to /year and thus /day [4]. Let us define and thus and . Therefore by Theorem 4, for any given initial value , the solution for the SDE HIV model (14) should obey Figure 3 clearly supports our analytical results given in Theorem 4 by showing the solution path of 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 for the SDE HIV model (14) and its corresponding deterministic HIV model (11) with step size using parameter values given in Example 10 with initial value .

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 Example 10 and reducing to , we have that , and By Theorem 4 and Proposition 5, we would expect the solution to persist and oscillate around the level . Furthermore by Proposition 5, as , we would expect to tend towards the deterministic equilibrium value for the corresponding deterministic model given by (11); namely, .
From Figure 4, we can clearly see that the solution path 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 for the SDE HIV model (14) and its corresponding deterministic HIV model (11) with step size using parameter values given in Example 11 with initial value .

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 Example 10 and choose to be , , , , and . We then let the simulations run for 1 million iterations but disregarding the first 800,000 iterations in order to allow to reach its recurrent level.
From Figure 5, 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, , , , , and corresponding to , , , , and , 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 for the SDE HIV model (14) with step size using parameter values given in Example 12 with initial value and , and .
Figure 6: QQ plot for the solution path for the SDE HIV model (14) corresponding to the histograms shown in Figure 5 for , and .

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 for any given initial value . Furthermore, we have constructed the basic reproduction number for the stochastic model, namely, , and the conditions required for extinction and persistence for our solution In general, if , the solution will almost surely go extinct as shown in Theorems 2 and 3. There is a gap in our results if and 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 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 that 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) becomeswith solution

Also a.s. Hence we expect that ifthen the disease dies out whereas if , the disease takes off. Thus this is a natural biological interpretation of the stochastic basic reproduction number .

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 thatwhere represents the basic reproduction number in the deterministic model. So in the deterministic model 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 number , then it is known that the disease will die out, whereas if , then the disease will persist. The results in this paper show that, in the stochastic model, if , then the disease dies out (almost surely), whereas if , 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 and take off if .

These findings provide new insights into the spread of HIV amongst PWIDs. Because the stochastic basic reproduction number is less than the deterministic one, it is possible for the noise to drive the disease to extinction, that is, if , 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 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.

Disclosure

The research data associated with this paper which was funded by EPSRC is available at http://dx.doi.org/10.15129/bf5ba8b5-d484-43c7-8006-14fb76819be2.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

Yanfeng Liang is grateful to the EPSRC for a studentship to support this work (EPSRC DTP Grant Reference no. EP/K503174/1) and the authors are grateful for the R coding provided by Drs. Alison Gray and Jiafeng Pan, Department of Mathematics and Statistics, University of Strathclyde, Glasgow, UK, which allowed the authors to modify the code to create the simulations shown in this paper. David Greenhalgh is grateful to the Leverhulme Trust for support from a Leverhulme Research Fellowship (RF-2015-88). Xuerong Mao is grateful to the Leverhulme Trust (RF-2015-385) and the Royal Society of London (IE131408) for their financial support.

References

  1. World Health Organization, “Global Health Observatory Data,” September 2015, http://www.who.int/gho/hiv/en/.
  2. World Health Organization, Media Centre, September 2015, http://www.who/int/mediacentre/factsheets/fs360/en/.
  3. World Health Organization, “HIV/AIDS,” September 2015, http://www.who.int/hiv/topics/idu/en/.
  4. E. H. Kaplan, “Needles that kill: modeling human immunodeficiency virus transmission via shared drug injection equipment in shooting galleries,” Clinical Infectious Diseases, vol. 11, no. 2, pp. 289–298, 1989. View at Publisher · View at Google Scholar
  5. D. Greenhalgh and G. Hay, “Mathematical modelling of the spread of HIV/AIDS amongst injecting drug users,” IMA Journal of Mathemathics Applied in Medicine and Biology, vol. 14, no. 1, pp. 11–38, 1997. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  6. S. M. Blower, D. Hartel, H. Dowlatabadi, R. M. Anderson, and R. M. May, “Drugs, sex and HIV: a mathematical model for New York City,” Philosophical Transactions of the Royal Society B: Biological Sciences, vol. 331, no. 1260, pp. 171–187, 1991. View at Publisher · View at Google Scholar · View at Scopus
  7. J. P. Caulkins and E. H. Kaplan, “AIDS impact on the number of intravenous drug users,” Interfaces, vol. 21, no. 3, pp. 50–63, 1991. View at Publisher · View at Google Scholar
  8. F. Lewis and D. Greenhalgh, “Three stage AIDS incubation period: a worst case scenario using addict-needle interaction assumptions,” Mathematical Biosciences, vol. 169, no. 1, pp. 53–87, 2001. View at Publisher · View at Google Scholar · View at MathSciNet
  9. D. Greenhalgh and F. Lewis, “The general mixing of addicts and needles in a variable-infectivity needle-sharing environment,” Journal of Mathematical Biology, vol. 44, no. 6, pp. 561–598, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet · View at Scopus
  10. E. H. Kaplan and E. O'Keefe, “Let the needles do the talking! evaluating the new haven needle exchange,” Interfaces, vol. 23, no. 1, pp. 7–26, 1993. View at Publisher · View at Google Scholar
  11. F. Lewis and D. Greenhalgh, “Three stage AIDS incubation period: a worst case scenario using addict-needle interaction assumptions,” Mathematical Biosciences, vol. 169, no. 1, pp. 53–87, 2001. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  12. T. Rhodes, M. Singer, P. Bourgois, S. R. Friedman, and S. A. Strathdee, “The social structure production of HIV risk among injecting drug users,” Journal of Social Science and Medicine, vol. 61, no. 5, pp. 1026–1044, 2005. View at Publisher · View at Google Scholar
  13. Y. Ding, M. Xu, and L. Hu, “Risk analysis for AIDS control based on a stochastic model with treatment rate,” Human and Ecological Risk Assessment, vol. 15, no. 4, pp. 765–777, 2009. View at Publisher · View at Google Scholar · View at Scopus
  14. G. R. Gupta, J. O. Parkhurst, J. A. Ogden, P. Aggleton, and A. Mahal, “Structural approaches to HIV prevention,” The Lancet, vol. 372, no. 9640, pp. 764–775, 2008. View at Publisher · View at Google Scholar · View at Scopus
  15. S. A. Strathdee, T. B. Hallett, N. Bobrova et al., “HIV and risk environment for injecting drug users: the past, present, and future,” The Lancet, vol. 376, no. 9737, pp. 268–284, 2010. View at Publisher · View at Google Scholar · View at Scopus
  16. N. Dalal, D. Greenhalgh, and X. Mao, “A stochastic model for internal HIV dynamics,” Journal of Mathematical Analysis and Applications, vol. 341, no. 2, pp. 1084–1101, 2008. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  17. H. C. Tuckwell and E. Le Corfec, “A stochastic model for early HIV-1 population dynamics,” Journal of Theoretical Biology, vol. 195, no. 4, pp. 451–463, 1998. View at Publisher · View at Google Scholar · View at Scopus
  18. N. Dalal, D. Greenhalgh, and X. Mao, “A stochastic model of AIDS and condom use,” Journal of Mathematical Analysis and Applications, vol. 325, no. 1, pp. 36–53, 2007. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  19. D. Peterson, K. Willard, M. Altmann, L. Gatewood, and G. Davidson, “Monte Carlo simulation of HIV infection in an intravenous drug user community,” Journal of Acquired Immune Deficiency Syndromes, vol. 3, no. 11, pp. 1086–1095, 1990. View at Google Scholar · View at Scopus
  20. D. Greenhalgh and F. Lewis, “Stochastic models for the spread of HIV amongst in-travenous drug users,” Stochastic Models, vol. 17, no. 4, pp. 491–512, 2001. View at Publisher · View at Google Scholar · View at MathSciNet
  21. A. Gray, D. Greenhalgh, L. Hu, X. Mao, and J. Pan, “A stochastic differential equation SIS epidemic model,” SIAM Journal on Applied Mathematics, vol. 71, no. 3, pp. 876–902, 2011. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  22. Z. Huang, Q. Yang, and J. Cao, “Complex dynamics in a stochastic internal HIV model,” Chaos, Solitons & Fractals, vol. 44, no. 11, pp. 954–963, 2011. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  23. S. W. Corson, D. Greenhalgh, and S. Hutchinson, “Mathematically modelling the spread of hepatitis C in injecting drug users,” Mathematical Medicine and Biology, vol. 29, no. 3, pp. 205–230, 2012. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  24. S. W. Corson, D. Greenhalgh, and S. J. Hutchinson, “A time since onset of injection model for hepatitis C spread amongst injecting drug users,” Journal of Mathematical Biology, vol. 66, no. 4-5, pp. 935–978, 2013. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  25. Q. Lu, “Stability of SIRS system with random perturbations,” Physica A: Statistical Mechanics and Its Applications, vol. 388, no. 18, pp. 3677–3686, 2009. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  26. X. Mao, G. Marion, and E. Renshaw, “Environmental Brownian noise suppresses explosions in population dynamics,” Stochastic Processes and Their Applications, vol. 97, no. 1, pp. 95–110, 2002. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  27. Q. Yang, D. Jiang, N. Shi, and C. Ji, “The ergodicity and extinction of stochastically perturbed SIR and SEIR epidemic models with saturated incidence,” Journal of Mathematical Analysis and Applications, vol. 388, no. 1, pp. 248–271, 2012. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  28. Q. Yang and X. Mao, “Extinction and recurrence of multi-group SEIR epidemic models with stochastic perturbations,” Nonlinear Analysis: Real World Applications, vol. 14, no. 3, pp. 1434–1456, 2013. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  29. E. Tornatore, S. M. Buccellato, and P. Vetro, “Stability of a stochastic SIR system,” Physica A: Statistical Mechanics and its Applications, vol. 354, pp. 111–126, 2005. View at Publisher · View at Google Scholar · View at Scopus
  30. X. Mao, Stochastic Differential Equations and Applications, Horwood, Chichester, UK, 2nd edition, 2008. View at Publisher · View at Google Scholar · View at MathSciNet
  31. R. Z. Khasminskii, Stochastic Stability of Differential Equations, Sijthoff & Noordhoff, Alphen aan den Rijn, The Netherlands, 1980.
  32. D. Greenhalgh and F. Lewis, “Modelling the spread of AIDS among intravenous drug users including variable infectivity, HIV testing and needle exchange,” Journal of Medical Informatics and Technologies, vol. 2, no. 1, pp. IP3–IP11, 2001. View at Google Scholar
  33. F. I. Lewis, Assessing the impact of variable infectivity on the transmission of HIV amongst intravenous drug users [Ph.D. thesis], University of Strathclyde, Glasgow, UK, 2000.
  34. F. Lewis and D. Greenhalgh, “HIV testing as an effective control strategy against the spread of AIDS among intravenous drug users,” Archives of Control Sciences, vol. 9, no. 1-2, pp. 69–96, 1999. View at Google Scholar
  35. S. W. Corson, Modelling the spread of hepatitis C virus infection amongst injecting drug users in Glasgow [Ph.D. thesis], University of Strathclyde, Glasgow, UK, 2012.