Modeling the spread of hepatitis C virus amongst people who inject drugs

This article discusses two models, with two different needle assumptions for the transmission of hepatitis C virus (HCV) between people who inject drugs (PWIDs) who share needles and syringes. Our analysis demonstrates that the basic reproduction number R0 determines how the model behaves. R0=1 is a crucial threshold parameter which divides two qualitatively different scenarios. It has been shown that if R0≤1 only a disease‐free equilibrium point exists. On the other hand if R0>1 a unique endemic equilibrium point exists as well as the disease‐free one. The disease‐free equilibrium point is globally asymptotically stable if R0≤1 otherwise unstable. We look at an approximation to this model by using the fact that the timescale on which the individual PWIDs inject with needles or syringes is much faster than the timescale of epidemiological change. Also, we showed that if R0>1 the endemic equilibrium point is locally asymptotically stable for our approximation to this model. Additionally, we perform some simulations with realistic parameter values which confirm that if R0>1 then the solutions to this model tend to this endemic equilibrium value giving a steady‐state nonzero HCV prevalence. With realistic parameter values for both assumptions we find that R0=2.9987 and the total endemic equilibrium fraction of infectious PWIDs is 0.5348 and the total equilibrium fraction of antibody positive PWIDs is 0.65. The total fraction of infectious needles is 0.275.

both sexes and every age are infected. To understand the worldwide HCV problem governments of each nation globally need to accurately record HCV disease prevalence. Information about the global spread must be determined using data from all over the world, but for most countries this information is not present. 2 At present there are around 58 million individuals globally who have long-lasting persistent or chronic HCV. 1 HCV is a virus spread by blood, which affects the liver cells, 1 and is mostly transmitted by different pathways, including common use of injection paraphernalia by individual people who inject drugs (PWIDs), insufficient cleaning of medical syringes and needles and other materials in health environments. Moreover it is possible for HCV to be spread by the transfusion of blood and blood products which have not been screened. However this mode of transmission has been uncommon in recent times because blood screening has improved considerably. A woman who is infected while expecting a child can also pass HCV to her baby. HCV can also be spread via sex. 1 In developed countries the most common way of acquiring HCV infection is intravenous drug use. The disease can occur in both acute and chronic forms, ranging from mild illness to a severe lifelong disease. Hepatitis C is an important cause of liver cancer. At present no vaccine for HCV exists although research into HCV vaccines is in progress. Antiviral medicines can reduce the risk of death from cirrhosis or liver cancer. It is important to mathematically model the spread of HCV then we can predict the number of people infected and the long-term health consequences.
Vickerman et al. 3 performed mathematical modeling to explain the manner in which HCV spread among PWIDs who live in London. The model helped in understanding changes in HCV infection among PWIDs. Also, the model was employed to assess the impact of intervention measures that were aimed at reducing the sharing of needles and syringes among all PWIDs, including those who have shared needles and syringes for more than a year. In addition, the model was used to determine how the intervention measures have impacted PWIDs with either a low or a high rate of sharing needles and syringes. Hutchinson et al., 4 on the other hand, relied on a computer model to stochastically mimic and understand how HCV would spread amongst Glasgow PWIDs particularly by common use of syringes and needles. Hutchinson and her co-authors specifically used the model to determine both the frequency and incidence of HCV among PWIDs in Glasgow, particularly between 1960 and 2000.
In the next section we give an outline and description of our model. We then derive an expression for the basic reproduction number R 0 . This is followed by an equilibrium analysis of our model equations. For R 0 ≤ 1 there is only the disease free equilibrium (DFE) point, whereas for R 0 > 1 the disease-free equilibrium (DFE) point exists and there is a unique steady state. We then show that if R 0 ≤ 1 then the system is globally asymptotically stable and always tends to the DFE point. This is done by bounding the limsup of all the other infection related variables by the limsup of the fraction of acutely infected PWIDs ∞ h and then deducing that if R 0 ≤ 1 this leads to a contradiction unless ∞ h = 0. Next we show that for R 0 > 1 the DFE point is unstable. We then give an argument to show that if the disease is initially present in PWIDs or needles and R 0 > 1 then the disease will persist. We bound the liminf of all other infection related variables below by terms involving the liminf of one of the acutely infected classes, then we use this result to bound the time that the proportion of this type of acutely infectious PWID can be beneath an arbitrarily low level and use this to show disease persistence.
We note that in general PWIDs inject every day or every few days, whereas the time in each disease class (acutely infectious, chronically infected and immune) is usually months or years. We use this to obtain a simpler approximation to our model. For this simpler approximation we show that the endemic equilibrium point is locally asymptotically stable if R 0 > 1. Some simulations with realistic parameter values conclude the paper.

PWIDs population
We develop a differential equation model for the spread of HCV by the common use of needles and syringes. The model is a modification of that of Corson et al. 5 which itself loosely used the model of Vickerman et al. 3 as its basis in that PWIDs were divided into susceptibles, two classes of acutely infectious, chronic and immune except that we have only one class of immune and we introduce a new class of second and subsequent time susceptibles. Instead of all the second acute infectious class individuals progressing to the immune class at the end of their acutely infectious phase some of these individuals may enter the new second and subsequent time susceptible PWID class.
A major novel innovation in our model compared to Corson et al. 5 and the other previous models is that we introduce a different set of needle equations. In Corson et al. 5 when an infected individual PWID uses a syringe or needle after use the syringe or needle takes on the infectious state of that PWID. However this may not necessarily be the case as data from HIV 6 demonstrates that HIV can be obtained from needles at very low dilutions. So once a syringe is infected it remains infected indefinitely. Based on this we assume that needles and syringes can only move up the levels of infectiousness. Also a second important novel aspect of our model compared to these models is that we have also introduced the effect of treatment of chronically infected PWIDs.
We take to be the average rate per individual PWID at which PWIDs use needles and syringes and the proportion of individual PWIDs who clean their needle before use is . For cleaning to be effective there must be no HCV virus in the syringe before it is used. Successful cleaning needs a single PWID to use bleach or alcohol to disinfect the needles or syringes that they use to inject drugs so that no HCV virus is left before the syringe is used.
Suppose that the total PWIDs population has fixed size n, a large constant number. The PWIDs can exit the population at per capita rate when susceptible PWIDs replace them immediately. There are two types of acute HCV infection, one leading to temporary infection and the other leading to permanent infection. x denotes the number of PWIDs who have never been infected and x 1 denotes the number of PWIDs who are currently not infected but have previously been infected. y is the number of PWIDs in the chronic stage and there are z PWIDs immune to infection. h is the probability that if a susceptible individual injects with a syringe in the state of Acute Infectivity (denoted h 1 or h 2 ) without cleaning the needle the individual PWID becomes infected. y denotes the corresponding probability of infection if the syringe is the state of chronic infectivity. Data shows that h > y . 3 Including separate chances of HCV spread for acute and chronic infection, similarly to Hutchinson et al. 4 is similar to the supposition that different HIV stages have different viral loads. 7 Moreover Greenhalgh and Lewis 8 use high and low HIV transmission probabilities for respectively acute and chronic disease. In our model we use a similar varying HCV transmission probability.
When individual susceptible PWIDs are infected they move on to acute infection stages (h 1 and h 2 ). Of these new infections a fraction will move into the h 2 acute category. When they leave this compartment these individual PWIDs get rid of their infection. A proportion 1 − become susceptible to being reinfected with HCV and a fraction become HCV immune. 9 The other fraction 1 − of newly infectious individual PWIDs move into the acute h 1 class: these individual PWIDs move on to the chronic disease class where they stay until death or stopping shared injection.

Needles and syringes
As well as modeling the individual PWIDs, we include needles and syringes (according to their state of infection). These can be uninfected, in the state of acute infection, or in the state of chronic infection. m denotes the total population size of needles and syringes. Only needles and syringes that have been used by individual PWIDs who are chronically or acutely infected can spread HCV. The infectiousness of each needle depends on the highest level of infectiousness that it has come into contact with. Unused needles are uninfectious. For example, suppose that an acute stage h 1 infected syringe is used for injection by a chronic y disease state individual PWID. The needle remains in the acute h 1 infectious state. On the other hand a needle in the chronic y diseased state used by an individual PWID in the acute h 1 infectious stage will be left in the acute h 1 infectious state. This assumption agrees with that made by Lewis and Greenhalgh 10 and makes this model pessimistic compared to other possible assumptions and so we expect that it may be possible to use our model to find a lower limit on the fraction of HCV infected syringes and needles and PWIDs.
The needle turnover rate is the mean rate that individual PWIDs swap used needles or syringes for unused ones. This is assumed as per year. We also suppose that needles and syringes which are not used do not lose their infectiousness. Chronically infected individual PWIDs can be treated with antiviral therapy which is successful over 95% of the time 11 although successfully treated individual PWIDs can catch HCV again. Cousien et al. 11 used a dynamic agent-based model of transmission amongst PWIDs to study optimal interventions to best manage HCV amongst PWIDs in France, in order to eliminate HCV. They take the new DAAs into consideration. Jia et al. 12 consider a mathematical model for hepatitis C in China. Infected individuals are divided into acute and chronically infected. After the acute stage individuals can progress to the chronic, immune or treated stages. We suppose that the per capita rate at which chronically infected individual PWIDs are successfully treated is .

Governing equations
The differential equations which describe HCV spread in the PWID population will now be explained. Individual PWIDs move through the different HCV disease stages described previously and HCV infection is caused by sharing the three types of infectious needle described also in the previous Section 2.2. We derive a total of eleven equations: six equations for PWIDs and five for needles. For needles we have two alternative assumptions for how needles in state h 1 and h 2 interact, Assumptions 1 and 2. Assumption 1 means that we regard state h 1 infectious needles as being slightly more infectious than state h 2 infectious needles and after use the level of infectivity of a needle originally in state h 1 or state h 2 infectivity is the highest of its state immediately before use and the level of infectivity in the blood of the individual PWID who last used it (as in Lewis and Greenhalgh 10 ). Assumption 2 means that we regard state h 1 infectious needles and state h 2 infectious needles as being equally infectious and assume that a needle remains uninfectious until first used then adopts the infectious state of the last individual PWID to use it (as in Corson et al. 5  At the time t let x (t) be the proportion of PWIDs in the x-susceptible class, x 1 (t) the proportion in class x 1 , h 1 (t) the proportion in the acute h 1 class, h 2 (t) the proportion in the acute h 2 class, y (t) the proportion in the chronic y class, and z (t) the proportion in the immune class. h 1 (t) is the proportion of syinges in HCV infectious stage h 1 , h 2 (t) is the proportion of syinges in HCV infectious stage h 2 and y (t) is the proportion of syinges in HCV infectious stage y. The gallery ratio = n∕m is the PWIDs to syringe ratio. Each susceptible individual PWID injects at rate . Recall that is the chance that an individual PWID will effectively cleanse the injection equipment before injection. The probability that he or she becomes infected at injection is So the per capita rate of attack of a single susceptible individual PWID is So the set of equations (shown in Figure 1) that give transmission of HCV amongst PWIDs is Recall that under Assumption 1 we regard state 1 infectious needles as being slightly more infectious than state 2 infectious needles. Under Assumption 1 for needles F I G U R E 1 Flow diagrams of the equations. The top part is for numbers of PWIDs and the bottom part is for needles (Assumption 1).
Both are given in terms of absolute numbers not fractions. In these figures = h 1 + h 2 + y is the total fraction of infectious PWIDs, For both assumptions we have that Equations (1)-(6) explain how the PWID fraction at each HCV infectious stage evolves with time, while (7)-(11) explain how the infected needle fractions change over time. For either assumption the initial conditions are x (0) ≥ 0, x 1 (0) ≥ 0, h 1 (0) ≥ 0, h 2 (0) ≥ 0, y (0) ≥ 0 and z (0) ≥ 0; h 1 (0) ≥ 0, h 2 (0) ≥ 0 and y (0) ≥ 0 with x (0) + x 1 (0) + h 1 (0) + h 2 (0) + y (0) + z (0) = 1 and h 1 (0) + h 2 (0) + y (0) < 1. Note that as the population is constant x (t) + x 1 (t) + h 1 (t) + h 2 (t) + y (t) + z (t) = 1 for all t. The arguments work mathematically provided that the initial condition of the system lies in this region. However we note that not all of these conditions are biologically realistic, for example it is virtually impossible that a population would start off with x (0) = 0 as there would be no first time susceptible individuals. This completes our derivation of the fundamental model. Next we will focus on the basic reproduction number R 0 which is a basic determinant of the behavior of the system.

THE BASIC REPRODUCTION NUMBER R 0
The basic reproduction number, R 0 , is defined as the amount of secondary cases caused when a single infected individual enters an equilibrium population with no disease (Corson et al. 5 ). A secondary infection is when someone is infected after using a needle which was infected by the original infectious individual. In other words for our model the number of secondary cases in PWIDs caused via only one infected needle by a single newly infected PWID entering a disease-free population at equilibrium, equivalently the number of secondary needles infected via only one infected PWID from a single infected needle entering an equilibrium population with no disease. Note that because of the treatment each chronically infected individual is infected for time on average 1∕( + ). Arguing as in Corson et al. 5

an individual newly infected PWID generates on average
acute h 1 , acute h 2 and chronic y infectious needles respectively. If we consider a single newly infectious needle then, writinĝ= this needle must remain infected rather than being exchanged or cleaned before the next injection with probability (1− ) Hence the probability that the needle remains infectious for PWIDs 1,2,3, … i − 1 but does not remain infectious for PWID i is Hence, E, the total expected number of PWIDs infected by the needle is Here * is the probability that the next PWID is infected so * = h if the needle is initially in infectious state h 1 or h 2 and * = y if the needle is initially in infectious state y. Multiplying this equation by 1− 1+̂a nd subtracting from (13) Hence summing the geometric progression, writing * = 1− +̂, and solving for E, Thus using (12) and (14) Here the basic reproduction number is defined as the number of HCV cases in PWIDs caused via exactly one infected syringe from a single newly infected PWID entering the DFE point. [13][14][15][16][17] An alternative approach [17][18][19][20][21][22] counts individual cases of disease to be either vectors (here needles) or hosts (here individual PWIDs). If this method is applied the new basic reproduction number, R * 0 say, is the unique positive real root of a cubic equation. As R * 0 > 1 if and only if R 0 > 1 and R * 0 < 1 if and only if R 0 < 1 this R * 0 has the same threshold value as R 0 . Define Applying the next generation matrix method there are six sorts of disease entities, h 1 , h 2 , and y for individual PWIDs, and m h 1 , m h 2 , and m y for needles. The next generation matrix is taken to be the matrix M = {M ij ∶ i = 1, 2, 3, 4, 5, 6, j = 1, 2, 3, 4, 5, 6}. Here M ij is taken as the average number of infections in disease state i resulting from one infectious disease stage j individual PWID who has just caught disease coming into the disease free steady state. Hence M is given by (16). It is straightforward to show that three of the eigenvalues of the characteristic equation of this matrix are zero and the remaining three satisfy f ( ) = 0 where As f (0) < 0 and lim →∞ f ( ) = ∞ the equation f ( ) = 0 has either one or three strictly positive real roots. By Descartes' rule of signs it has either zero or one strictly positive real roots, hence has exactly one positive real root. Moreover as the next generation matrix (NGM) M is a positive irreducible matrix it has a real strictly positive eigenvalue corresponding to its spectral radius (Lemma 2.1 in Nold 23 ), which is R * 0 , the basic reproduction number by the NGM method. Moreover It is clear that the expression for R 0 is different than the corresponding expression in Corson et al. 5 because of the expressions involving and the factor +̂in the denominator. Moreover the derivation using the NGM was not discussed in Corson et al.. 5 We next analyze the behavior of our system. We are especially focusing on conditions which cause HCV to go extinct or remain persistent.

ANALYTICAL RESULTS
Now we look at how our differential equation model behaves. We look at when the disease will persist and when it goes extinct. We facilitate this by carrying out an equilibrium and stability analysis. There is always a disease-free equilibrium and there may additionally be a unique endemic equilibrium. The latter is a possibility only if R 0 > 1. We shall demonstrate that the solution with no disease present will always be approached if R 0 < 1 but for R 0 > 1 this equilibrium is unstable. We derive an approximation to our differential equation model which has the same basic reproduction number and equilibria. We show that for this approximation model the nonzero solution is locally asymptotically stable when R 0 > 1. The general structure of the argument follows Corson et al. 5 but there are significant differences due to the different needle equations and the introduction of treatment. Because there is the PWID treatment term and the needle equations are different in the two models the results of this paper do not follow directly from the results of Corson et al. 5 During the course of this analysis we assume that the probability of effectively cleansing the syringe or needle , is greater than or equal to zero but strictly less than one but cannot take the value one. If we allow = 1 then R 0 = 0 and no disease transmission will take place, = 0 is allowed since this corresponds to a scenario where PWIDs do not engage in cleaning practices prior to injecting. In addition, we assume that the remaining model parameters are strictly positive.

4.1
Equilibrium solutions Theorem 1. In either system (1)- (8) and (11) or (1)-(6) and (9)- (11), if R 0 is less than or equal to 1, there is a unique disease-free equilibrium with no disease present in either PWIDs or syringes. If R 0 > 1 as before we have the disease-free equilibrium (DFE), but additionally there is exactly one steady state with disease present.
to be the fraction of PWIDs in class i and * j to be the fraction of needles in infectious class j at the steady state with disease present. Setting the right side of the system equal to zero in either system (1)-(8) and (11) or (1)-(6) and (9)-(11) and writing * (3) and (4), we find that * .

For notational simplicity write
and Hence from the equilibrium versions of these equations * .
Adding the equilibrium versions of Equations (3) and (4) we deduce that k = 0 is always a possible solution of (19) corresponding to the DFE solution. Any other nonzero solution must satisfy Let g(k) denote the right-hand side of (20). g(k) is strictly monotone decreasing in k with g(0) = R 0 . Hence if R 0 ≤ 1, g lies beneath one for k > 0. So (20) has no strictly positive solutions for k.
It is straightforward to show that this strictly positive solution corresponds to an endemic equilibrium for either Assumptions 1 or 2. Hence it has been demonstrated that if R 0 > 1 there is the DFE and exactly one feasible steady state with disease present with all of * , and * y strictly positive and if R 0 ≤ 1 there is only the DFE. ▪

Global stability of the disease free equilibrium
The case 0 ≤ R 0 ≤ 1 will now be analyzed. Here our result will be that HCV is eventually eliminated in every PWID and every needle and syringe. . Hence we write these results in terms of just one limit supremum. R 0 < 1 then ensures that this limit supremum is zero. Theorem 2 then follows. Note that this proof does not use the Lyapunov method. 24 We were unable to use the Lyapunov method for our model because we were unable to find a suitable Lyapunov function. ▪ Proof. From Equation (5), we have that To solve by using an integrating factor, we have Let t ≥ t 1 ( ) and integrating over [t 1 ( ), t] gives Note there exists t 2 ( ) > t 1 ( ) such that exp(−( + )(t − t 1 ( ))) ≤ . For t ≥ t 2 ( ) Taking the limsup we deduce that As is arbitrary so is 1 and letting 1 tend to zero we deduce that ∞ y ≤ Proof.
Arguing as in Lemma 1 Corollary 1 is true. ▪ We shall next use the same method and Equations (7)- (8) and (11) for Assumption 1 and Equations (9) (7) we have that For Assumption 2 from Equation (9) we have that

So in both cases
So multiplying by exp( ( +̂)t) we have Proceeding as in Lemma 1

+̂.
Similarly for h 2 for Assumption 1 from Equation (8) we have that And from Equation (10) we have that Again proceeding as in the proof of Lemma 1 we deduce that

+̂.
For y for both Assumptions 1 and 2 we have that So once again we have that ∞ y ≤ ∞ y +̂.

Now for Assumptions 1 and 2, we have that
The next step in our argument is to derive an equation Lemma 2 permits us to express the results of Lemma 1, Corollary 1 and Equations (21) in terms of ∞ h to obtain: and ∞ y ≤ Adding Equations (3) and (4) together we have Using the upper limits for ∞ , and ∞ y given by inequalities (25)- (27) yields where 2 = (1− ) ( + ) . Using Equation (15) we obtain Integrating over [t 5 ( ), t] and applying a similar argument to Lemma 1 we find that , and the positive constant 5 can be made as small as we want. Letting 5 tend to zero and dividing by ∞ h we deduce that 1 + ∞ h R 0 < R 0 which contradicts 0 ≤ R 0 ≤ 1 and so ∞ h = 0 if 0 ≤ R 0 ≤ 1. Note that this proof is still valid even if , ∞ y are zero using (23)- (27) and the comments immediately before these inequalities. Therefore, we must have It is then straightforward to show that lim t→∞ x 1 (t) = 0 and lim t→∞ x (t) = 1. We have now shown that for 0 ≤ R 0 ≤ 1 all solutions must tend to the DFE which is globally stable.
Thus we have shown that if R 0 ≤ 1 the DFE is globally asymptotically stable. We next turn our attention to the situation where R 0 > 1, recall that in this case there is a unique endemic equilibrium. Our first result is the following lemma which shows that if R 0 exceeds one the DFE is unstable.
Proof. We look at the linearized system of Equations (1)- (11) which are evaluated at the DFE. When linearizing the population dynamics in the neighborhood of the DFE point (1,0,0,0,0,0,0,0,0) we find that it can be written as a matrix and, for both Assumptions 1 and 2, J is given by We want to look at the local asymptotic stability of the DFE, determined by the eigenvalues of this matrix. If it is unstable we desire to demonstrate that there is one or more eigenvalue with strictly positive real part. The characteristic polynomial of J has leading term ± 9 . If the equilibrium is unstable it is enough to demonstrate that the constant term a 9 < 0 (multiplying the polynomial by −1 if necessary to ensure the first term is 9 ). It is easy to demonstrate that It is clear that this term is negative when R 0 > 1 and the result follows. ▪ We shall now look at HCV persistence when R 0 > 1. h 1 and h 2 as the initial infectious stages drive the course of the infection, and hence, if they become small then so do the other terms. This statement is shown mathematically in the proof of Theorem 2. Our main persistence result is Theorem 4 which states that if R 0 > 1 then it is not possible that these terms can become as small as we like, hence they can be bounded away from zero. So if the disease is present at the start in either needles or PWIDs and R 0 > 1 then infection will persist indefinitely.
This theorem needs several steps like Theorem 2. For a function f write f ∞ = liminf t→∞ f (t). Our first step in the persistence proof is to show that the persistence of HCV in both PWIDs and needles will follow if we have h 1 ,∞ > 0 instead of needing to bound all separate disease classes below. Hence we need some preliminary lemmas and inequalities (Lemmas 3-5 and inequalities 28).  . In a similar way, to prove Theorem 4 we need to establish a connection between h 1 ,∞ and h 2 ,∞ .

Lemma 5.
( Using the same method as Lemma 3 and Equations (7), (9), (17), and (18) for both Assumptions 1 and 2 it is easy to demonstrate that . (28) Note that this argument is different from the argument in Corson et al. 5 because of the way we have to combine the 's. From Lemmas 3-5 and inequalities (28) we see that it is enough to demonstrate that h 1 ,∞ > 0 to bound trajectories from the origin. , for t > .
Proof. Note that this argument is different to Corson et al. 5 because it is more complicated and we have to deal with h 1 , h and instead of h 1 , h 2 and y . Also the bounds for the 's are different. We take to be a small constant. We have either . Then using Lemma 5 we find that Arguing similarly to the above there is T 2 with for t ≥ T 2 then } is the next instant following t = t 0 where h 1 rises then by the definition of t 0 , we have h 1 (t 0 + ) < 1 2 * h 1 for some small and positive.
Hence, t 1 > t 0 and by continuity h 1 (t 0 ) = 1 just after t 1 . We shall next demonstrate that as h 1 becomes small then so do the rest of the variables. ▪ Lemma 7. Suppose that * > 0 is small and positive, then there is U * 1 > 0 such that 0 < y < ( 1 2 + * ) * y if t 0 + U * 1 < t 1 and t ∈ [t 0 + U * 1 , t 1 ].U * 1 is a function only of the model constants, * and .

Proof. See Appendix A. ▪
In Lemma 7 we have shown that h 1 becomes small then this causes y to become small. We shall now prove similar results for h 2 , z , h 1 , h and .
is a function only of the model constants, * and .
is a function only of the model constants, * and .
for Assumption 1 and for Assumption 2.U * 4 is a function only of the model constants, * and .
for both Assumptions 1 and 2. U * 5 is a function only of the model constants, * and .
Proof. See Appendix A. ▪

Lemma 12.
There is U * 6 > 0 with for both Assumptions 1 and 2. U * 6 is a function only of the model constants, * and .

Proof. See Appendix A. ▪
We have demonstrated that if h 1 approaches zero then the proportions of all other infectious individual PWID and needle stages do so as well. Now we demonstrate h 1 cannot decrease below a fixed limit. We do this by showing that there is a constant which is an upper bound for t 1 − t 0 and so, h 1 cannot be beneath 1 2 * h 1 sufficiently long to approach zero. Note that We want to show that t 1 − t 0 < T. Here T < ∞ is a function only of * , and the model parameters. In option 2 our result is true. We next discuss option 1 where all the terms are small by time t 1 . As the DFE is unstable when R 0 > 1 we shall demonstrate that it is not possible for h 1 to be small for an unbounded length of time. Although the idea is similar there are significant differences to the proof of Corson et al., 5 in particular a different co-ordinate system is used (see the proof of Lemma 14). (10) . Let F 1 ( , ) be a polynomial of nth degree in and . Denote the roots (possibly complex) of F 1 ( , ) = 0 by j for j = 1, … , n. Then in a neighborhood of = 0 each j is defined and continuous.

at a time t 0 , then h 1 (t) comes back to this value by time no later than t
where T is a function which depends on the system constants, , * and nothing else. The same argument holds whenever h 1 drops beneath 1 . A modification of this argument also works whenever h 1 drops beneath 1

Local stability of the endemic equilibrium
To show that the endemic equilibrium is locally stable we need to look at the Routh-Hurwitz criteria for the characteristic polynomial whose leading term is 9 . This is impractical due to level of complexity and difficulty. In its place, a similar approximation model will be used that has five dimensions. The first concept in reducing the complexity of the system is that individual PWID demographic processes are much slower than the timescale on which individual PWIDs inject and the former are measured in years whilst the latter are measured in days. This is similar to the way that Kaplan It is possible to simplify the system further by eliminating the x equation. It is known that , y → * y , and z → * z then x → * x . So if everything else is known it is possible to determine the limiting behavior of x . Therefore our approximate model is as follows: and It is important to note that for Equations (29)-(33) the equilibrium values and R 0 are the same as in our original model with needles included. Recall that * = 1− +̂. If we denote R * * to be the basic reproduction number obtained by using the NGM it satisfies characteristic equation f ( ) = 0 where and a similar argument as for the original model shows that R 0 and R * * 0 have the same threshold value. Theorem 5. The model (29)-(33) has a locally stable endemic equilibrium for R 0 > 1.
where k satisfies the equation As in Corson et al., 5 one eigenvalue is = −( + ) and the remaining eigenvalues satisfy the characteristic equation where Substituting Equation (34) into A 1 we can write Substituting (35) in (39) we can express A 2 in terms of a 3 as Expressing + − Ka 1 in terms of K and a 3 , we have Similarly writing and substituting for + − Ka 1 and a 2 in terms of K and a 3 we see that To show stability of the endemic equilibrium we need to show that the Routh-Hurwitz conditions are satisfied. These are A 1 , A 2 , A 3 > 0, and A 1 A 2 > A 3 . From (41), (42) and (44) it is clear that A 1 , A 2 , and A 3 > 0. To show A 1 A 2 > A 3 note that: Now we shall show that (2 + ) 2 > 3 . Note that (2 + )(2 + + (1 − )) > ( + + (1 − )).
Hence if we consider (2 + ) multiplied by the first term in the second pair of square brackets in (43), it exceeds multiplied by the second term in the second pair of square brackets in (46). Also as (2 + ) multiplied by the second term in the second pair of square brackets in (43), exceeds multiplied by the first term in the second pair of square brackets in (46). Similarly as the terms multiplying ( h − y ) in (2 + ) 2 exceed those multiplying ( h − y ) in 3 . So The result that A 1 A 2 > A 3 follows from (47) using (45) and (48). ▪

Parameter combinations resulting in R 0 > 1
These analytical findings will now be shown numerically. We wish to verify that our extended model has a globally stable DFE if R 0 ≤ 1 and explore whether either extended system (1)-(8) and (11) or (1)-(6) and (9)-(11) tends to the unique endemic equilibrium if R 0 > 1, and there is initially HCV in either syringes or PWIDs. We now use Wolfram Mathematica version 11.1 to give simulations of HCV prevalence in Glasgow, Scotland, over time. The results were checked in detail by examining the output from extensive simulations. We do not have new data to calculate the epidemiological constants of the model. So we decided to use the existing scientific literature: Corson et al., 5 Hutchinson et al., 4 Micallef et al., 26 Vickerman et al., 3 and Vickerman et al. 27 Using the parameters in Table 1 we calculated R 0 = 2.9987 > 1. The parameters apart from were used by Corson et al. 5 in their simulations. The parameters , and were estimated from Glasgow survey data. is estimated to be 0.26. 26 was conservatively estimated to be 0.25. We use parameter estimates for h and y taken from the literature (De Carli et al., 28 Gerlach et al., 29 Kamal et al., 30 and Vilano et al. 31 ). is estimated to be 2 per annum following Vickerman et al. 3,27 and supposing that the acute HCV stage lasts half a year.
Martin et al. 32 consider a model in which up to 6% of continuing PWIDs were treated annually. This is in line with observed treatment rates in the UK up to 2015 where at most 3% of PWIDs were treated annually. 33 Scott et al. 34 discuss a mathematical model to investigate HCV elimination in line with WHO targets in Iceland, using a complex mathematical model. The model estimated that an 80% reduction in domestic HCV incidence was achievable by 2030, 2025, or 2020 if at least 5.5%, 7.5%, and 18.8% of PWIDs were treated per year. However treatments for HCV have improved dramatically in recent years moving from relatively expensive and ineffective interferon based treatments to much more effective and cheaper directly acting antiviral (DAA) treatments. So recently there has been a huge increase in the number of PWIDs being treated. 35,36 Pitcher et al. 37 review the epidemic modeling literature on HCV transmission and prevention amongst PWIDs. They describe the basic mathematical modeling approaches. They survey the use of dynamic models incorporating a similar deterministic susceptible-infected-susceptible SIS compartmental structure although some use a stochastic model. They present levels of combinations of various control strategies (DAA treatments with 0-90 per 1000 PWIDs treated per annum, opiate substitution therapy and needle exchange programs) to achieve WHO elimination targets (80% reduction by 2030). These are new simulations but using a previously published model. 38 Traeger et al. 36 state that in Australia DAAs have moved from a cumulative total of less than 1% of RNA tested PWIDs receiving treat-T A B L E 1 Parameter values. Sources for these parameter estimates are given in the text. The parameters apart from were used in Corson et al. 5 with different PWID needle interaction assumptions and no treatment of PWIDs

Time (Years) HCV Prevalence
F I G U R E 2 HCV amongst Glasgow PWIDs who share injections (given by the unbroken black line) and infectious syringes (given by the broken red line). By taking the parameters as in Table 1 we calculate R 0 = 2.9987 > 1 ment to a cumulative total of around 45% of RNA tested PWIDs (note that this is RNA tested PWIDs not all PWIDs). Moreover World Health Organization forward targets are even more ambitious with 80% of PWIDs to receive treatment by 2030. 39 Hence we take the number of continuing PWIDs treated annually to be 10% as it lies in the middle of these estimates.
We estimate the spread of HCV in our model when R 0 = 2.9987 over a period of 25 years. First we took 1% of the PWIDs to have acute HCV (h 1 ) disease and the other PWIDs or needles and syringes were free of disease. Hence, x (0) = 0.99, x 1 (0) = 0, h 1 (0) = 0.01, h 2 (0) = 0, y (0) = 0, z (0) = 0 and h 1 (0) = h 2 (0) = y (0) = 0 for both assumptions. We follow Corson et al. 5 and Health Protection Scotland and take x 1 + h 1 + h 2 + y + z to be the prevalence. So we count PWIDs who have already experienced infection in the definition as these will be positive in an antibody test. The total fraction of infected needles is also given. The results in Figure 2

5.2
Parameter combinations resulting in R 0 ≤ 1 In this section we focus on combining parameters and finding the range of values that allow for HCV elimination. We begin by assuming that = 48.10 per annum, = 0.763 and = 637.152 per annum. We then simulate our model behavior using these values and evaluate the time taken for our system to reach a disease free state. Figure 3 shows the resulting model behavior when we combine these parameters. Initially we assume that 70% of PWIDs are infected with HCV: 10% are in the acute h 1 disease stage, 10% belong to the state h 2 disease class and 50% belong to the state y PWID disease class. That is, x (0) = 0.3, h 1 (0) = 0.1, h 2 (0) = 0.1, y (0) = 0.5 and x 1 (0) = z (0) = 0. In addition, we assume that 30% of needles and syringes are infected with HCV with the disease spread evenly among the three infectious needle classes. Hence, h 1 (0) = h 2 (0) = y (0) = 0.1. HCV dies out in both PWIDs and needles since R 0 = 0.0492. Figure 3 shows that HCV elimination in our population is possible in approximately 25 years.
We performed other simulations with other parameter estimates and other starting points for both R 0 less than or equal to one and R 0 greater than one. In each case the results of the simulations verified our analytical results. If R 0 was less than or equal to one HCV always died out and for R 0 greater than one HCV approached the steady state with disease present, which was the unique steady state with disease present, provided that the infection was there at the start. This suggests that if R 0 > 1 the steady state with disease present is globally asymptotically stable in the sense that if infection is present at the start the system will tend to the steady state with disease present.

CONCLUSIONS AND DISCUSSION
This article has taken the model for the spread of HCV amongst PWIDs studied by Corson et al. 5 and introduced two sets of alternative PWID needle interaction assumptions. Instead of assuming that the needle always adopts the infectious state of the last user the models assume that needles cannot lose infectivity and get successively more infectious over time. Additionally in recent years treatment with DAAs have offered other options for treatment of HCV in addition to needle cleaning and needle exchange, so this also has been included in the model. Both the different PWID needle interaction assumptions result in more disease transmission, and as the PWID needle interaction assumptions are different and more complicated, new mathematical ideas are needed. This is particularly true of the existence and uniqueness of the equilibria and the local stability of the endemic equilibrium in the approximate model. An expression for the basic reproduction number R 0 was derived. Another novel aspect of our work compared with previous work is that we used the NGM approach to find R 0 . We performed an equilibrium analysis by analyzing the solutions of the Equation (20) where the right-hand side is monotone decreasing in k. The dynamics of the model are given by the basic reproduction number R 0 . R 0 = 1 is a crucial boundary above which HCV will become endemic. It has been shown that if R 0 ≤ 1 then the model progresses toward the DFE. This was done by bounding the limsup of all the nonsusceptible class variables in terms of the limsup of h , the fraction of acutely infectious PWIDs, as these initially infectious PWIDs drive the disease propagation. By examining the Jacobian of the linearized equation about the DFE point we show instability of the DFE when R 0 > 1. We then outline an argument to show the persistence of HCV when R 0 > 1 and disease is initially present. We obtain lower bounds for the liminf of each of the nonsusceptible disease class proportions in terms of the liminf of h 1 . So it is enough to show that h 1 ,∞ > 0. We show that h 1 becoming small, in the sense h 1 ,∞ < 1 2 * h 1 for some > 0, implies that all of the other nonsusceptible disease class proportions become small in a time which is a function only of the model constants, * and , where * is an arbitrarily small quantity. We then use the instability of the DFE point when R 0 > 1 to show that h 1 must rise again to the level 1 2 * h 1 if is small enough and the time for this to happen is a function only of the model constants, * and . From this we can deduce a contradiction to h 1 ,∞ = 0 and hence the result follows.
We then obtained an approximation to the model by using the fact that the timeframe in which individual PWIDs inject is much quicker than the time for the PWIDs to move between disease states and used the Routh-Hurwitz criteria to demonstrate local asymptotic stability of the endemic steady state point in the approximation to our model if R 0 > 1.
The results of the simulations confirmed the results of our analyses and enabled the model to be used to estimate how to intervene to ensure that R 0 ≤ 1 and consequently to drive out HCV among PWIDs and syringes. With the parameters as in Table 1 we get R 0 = 2.9987. By keeping all but one of these parameters fixed and varying the remaining parameter we can find the critical value of the latter which will just eliminate the disease. Applying this method to the critical control parameters which are the per capita needle common injection rate , the syringe disinfection probability and the individual syringe exchange rate gives, = 48.10 per annum, = 0.763 and = 637.152 per annum. With all of the other parameters fixed as in the table and the per individual chronic PWID treatment rate varying it is just possible to eliminate HCV but the value of required to do this is so high as not to be realistically achievable. We have performed further simulations using this model to explore the effect of different control measures such as needle cleaning, needle exchange and treatment on both disease prevalence and disease elimination. We have also investigated sensitivity of the disease prevalence to model parameters and numerically compared these results with the results of the simpler model of Corson et al. 5 but do not have space to present these. We have also modified the model to introduce heterogeneity amongst drug users, both in injection rate and time since onset of injection, but these results will be discussed elsewhere. Now note that exp(−( + )(t − t 6 ( ))) → 0 as t → ∞. So there exists t 7 ( ) > t 6 ( ) such that for t ≥ t 7 ( ), As is arbitrary the result of Lemma 3 follows. So for any > 0 there is t 7 ( ) > 0 so if t ≥ t 7 ( ), Moreover we can find t 8 > t 7 such that for t( ) ≥ t 8 ( ) then ) .