Computational and Mathematical Methods in Medicine
Volume 2016 (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

Yanfeng Liang, David Greenhalgh, and Xuerong Mao

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 ,

where

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:

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 be

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 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 equation

By substituting (10) into the second equation in (7) we deduce that

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 , 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 that

Hence, there is an integer such that

Let us define a function ,

Now by Itô’s formula, we have that, for any and ,

where is defined by

Furthermore, it is easy to see that

where

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 that

Letting , 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 number

This 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 number

then, 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 that

Here is defined as

where 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 that

However, 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. If

then, for any given initial value , the solution of (14) obeys

In 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 at

Note 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 function

where . 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 that

In 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 and

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

then, for any given initial value , the solution of (14) satisfies

where

which is the unique root in of the function

defined 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 is

where 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 that

where 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 that

Clearly 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 that

Thus, 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 for

and then is strictly decreasing and

which is the equilibrium state of the deterministic HIV model (11) and

In 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 that

Then,

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 to

As 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 that

Let 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 that

and then

By letting , we have that

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 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⁻⁴/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) becomes

with solution

Also a.s. Hence we expect that if

then 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 that

where 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

World Health Organization, “Global Health Observatory Data,” September 2015, http://www.who.int/gho/hiv/en/

World Health Organization, Media Centre, September 2015, http://www.who/int/mediacentre/factsheets/fs360/en/

World Health Organization, “HIV/AIDS,” September 2015, http://www.who.int/hiv/topics/idu/en/

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

Q. Lu, “Stability of SIRS system with random perturbations,” Physica A: Statistical Mechanics and Its Applications, vol. 388, no. 18, pp. 3677–3686, 2009.

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.

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.

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.

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.

X. Mao, Stochastic Differential Equations and Applications, Horwood, Chichester, UK, 2nd edition, 2008.

R. Z. Khasminskii, Stochastic Stability of Differential Equations, Sijthoff & Noordhoff, Alphen aan den Rijn, The Netherlands, 1980.

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.

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.

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.

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.