Natasha K Martin
email: natasha.martin@bristol.ac.uk

Peter Vickerman

Matthew Hickman

Mathematical modelling of Hepatitis C treatment for injecting drug users

Keywords: HCV, antiviral treatment, liver disease, IDUs, mathematical model

Hepatitis C virus (HCV) is a blood-borne infection that can lead to progressive liver failure, cirrhosis, hepatocellular carcinoma and death. In developing countries, the majority of HCV infections are transmitted via injecting drug users (IDUs). Despite effective antiviral treatment for HCV, very few active IDUs are treated. Reluctance to treat is partially due to the risk of reinfection. We develop a mathematical model of HCV transmission amongst active IDUs, and examine the potential effect of antiviral treatment. As most mathematical models of interventions utilise a treatment function proportional to the infected population, but many policy implementations set fixed yearly targets for specific numbers treated, we study the effects of using two different treatment terms: annually treating a proportion of infecteds or a fixed number of infecteds. We examine the behaviour of the two treatment models and find different bifurcation behaviours in each case. We calculate analytical solutions for the treatment level needed for disease clearance or control, and observe that achievable levels of treatment can result in control or eradication across a wide range of prevalence levels. Finally, we calculate the sensitivity of the critical treatment threshold to the model parameters, and find that for a given observed prevalence, the injecting duration and infection risk play the most important role in determining the treatment level needed. By contrast, the sensitivity analysis indicates presence (or absence) of immunity does not alter the treatment threshold. We conclude by discussing the public health implications of this work, and comment on the importance and feasibility of utilising treatment as prevention for HCV spread amongst IDUs.

Introduction

Hepatitis C virus (HCV) is a blood-borne disease with an estimated global prevalence of 2-3%, or 130-170 million people, and is one of the leading causes of chronic liver disease [START_REF] Shepard | Global epidemiology of hepatitis C virus infection[END_REF].

If left untreated, about 7-18% of those infected will progress to liver disease within 20 years, which can result in progressive liver failure, cirrhosis, hepatocellular carcinoma and death [START_REF] Seeff | The history of the "natural history" of hepatitis C (1968-2009)[END_REF].

In developed countries, the primary mode of transmission is amongst injecting drug users (IDUs) through needle and syringe sharing, with over 80% of new cases in the UK attributed to injecting drugs [START_REF] Acmd | The Primary Prevention of Hepatitis C Among Injecting Drug Users[END_REF]. HCV is easily transmitted amongst IDUs, with 15-90% of IDUs testing positive for HCV antibodies [START_REF] Page-Shafer | Testing Strategy To Identify Cases of Acute Hepatitis C Virus (HCV) Infection and To Project HCV Incidence Rates[END_REF][START_REF] Judd | Incidence of hepatitis C virus and HIV among new injecting drug users in London: prospective cohort study[END_REF][START_REF] Hahn | Hepatitis C Virus Seroconversion among Young Injection Drug Users: Relationships and Risks[END_REF]. Current preventative measures to reduce HCV transmission such as health education and advice, needle and syringe exchange, and opiate substitution therapy aim to prevent transmission by reducing unsafe injecting [START_REF] Acmd | The Primary Prevention of Hepatitis C Among Injecting Drug Users[END_REF]. However, public health surveillance indicates substantial decreases in prevalence have not been achieved [START_REF] Palmateer | Evidence for the Effectiveness of Sterile Injecting Equipment Provision in Preventing Hepatitis C and HIV Transmission among Injecting Drug Users: A Review of Reviews[END_REF].

HCV antiviral treatment (peginterferon-α and ribavirin) is effective, resulting in viral clearance in 45-80% of cases, depending on HCV genotype [START_REF] Nice | Peginterferon alfa and ribavirin for the treatment of mild chronic hepatitis C. NIH[END_REF]. Prior to 2002, guidelines in the US and UK recommended against treating active IDUs. However, current guidelines now do not exclude IDUs from treatment eligibility, given mounting evidence that IDUs exhibit a similar response to treatment, and are just as compliant with treatment as ex or non-IDUs [START_REF] Hellard | Hepatitis C Treatment for Injection Drug Users: A Review of the available evidence[END_REF][START_REF] Nice | Peginterferon alfa and ribavirin for the treatment of mild chronic hepatitis C. NIH[END_REF][START_REF] Shepherd | Interferon alfa (pegylated and non-pegylated) and ribavirin for the treatment of mild chronic hepatitis C: a systematic review and economic evaulation[END_REF]NIH, 2002). Nevertheless, despite these recommendations and the high numbers of IDUs infected, very few (<3-4%) active IDUs have ever been treated [START_REF] Grebely | Uptake of hepatitis C virus (HCV) treatment among injection drug users (IDUS) in Vancouver, Canada[END_REF][START_REF] Seal | Among injection drug users, interest is high, but access low to HCV antiviral therapy[END_REF]. Studies on treatment barriers have indicated a reluctance to treat active IDUs due to the possibility of subsequent reinfection [START_REF] Booth | Clinical guidelines on the management of Hepatitis C[END_REF][START_REF] Reimer | Guidelines for the Treatment of hepatitis C Virus Infection in Injection Drug Users: Status Quo in the European Union Countries[END_REF][START_REF] Foster | Injecting drug users with chronic hepatitis C: should they be offered antiviral therapy?[END_REF].

We examine the potential of antiviral treatment as a prevention strategy for HCV amongst IDUs. By using antiviral treatment to reduce prevalence amongst active IDUs, the treatment can act to reduce the risk of infection for other IDUs. But to what extent? This paper examines the potential impact of HCV treatment on prevalence and transmission, including the possibility of reinfection. We incorporate two treatment scenarios (treating a proportion of infected IDUs, and a fixed number of IDUs) and examine the resulting dynamics and treatment needed for eradication. Treating a constant proportion of the population is the function most commonly used in infectious disease modelling. However, annually treating a fixed number of IDUs would be more likely in the initial stages of a treatment delivery program, or in situations with budget constraints. Hence, we analyse both situations.

Background and assumptions for the model

Infection with HCV leads to a brief acute stage, which is relatively short (on the order of weeks to months) in comparison to the prolonged chronic stage (on the order of decades) [START_REF] Ecmdda | Hepatitis C and injecting drug use: impact, costs and policy options[END_REF]). In the first few weeks, viral levels may be undetectable, increasing but possibly remaining low during the remainder of the acute stage. A fraction of people (about 26%) spontaneously clear the acute infection [START_REF] Micallef | Spontaneous viral clearance following hepatitis C infection: a systematic review of longitudinal studies[END_REF]. The specifics of spontaneous clearance are not well known, although women and young adults exhibit higher spontaneous recovery rates. Due to the relatively short duration of the acute stage and the small fraction who spontaneously clear, we neglect the small contribution towards infections from the acute IDUs who spontaneously clear. Those who spontaneously clear either become susceptible again, or may become immune. The concept of sterilising immunity following exposure to HCV is uncertain. We assume a low proportion become immune, and explore the sensitivity of the model with respect to immunity in the sensitivity analysis. The remaining fraction which do not spontaneously clear the acute infection progress to the chronic infection stage.

There are six identified HCV genotypes (numbered 1-6), with different distributions among geographically distinct IDU populations. In the UK, for example, genotype 1 comprises about 50% and genotype 2 and 3 together comprise about 50% [START_REF] Nice | Peginterferon alfa and ribavirin for the treatment of mild chronic hepatitis C. NIH[END_REF]. In the US, the proportion of genotype 1 is slightly higher (about 70%). The differences in disease progression between the genotypes is not yet clear, but they do show differences in response rates to therapy, with genotypes 2 and 3 exhibiting higher cure rates than genotype 1. Treatment with peginterferon and ribavirin results in a sustained viral response 6 months after treatment in 40-50% of people with genotype 1, and 75-85% with genotype 2 or 3 [START_REF] Nice | Peginterferon alfa and ribavirin for the treatment of mild chronic hepatitis C. NIH[END_REF]. In this model, we do not explicitly model infections with different genotypes, and instead track total infections and use a weighted average cure rate. Additionally, we examine a worst-case scenario with a population comprised entirely of the harder to treat genotype 1, in case the differential treatment success rates alters the population genotype distribution. Further, countries such as the United States have a higher proportion of genotype 1 and so would tends towards these scenarios (Klevens et al. , n.d.).

Antiviral treatment leads to a substantial reduction in viral load in the first few weeks (even among some eventual nonresponders). Hence, we assume that IDUs currently on treatment are non-infectious.

Due to the lack of evidence to suggest otherwise, we assume that the chances of spontaneous clearance and immunity are equal for naive and re-infected IDUs. Furthermore, we assume that the probability of treatment success is the same between naive and re-infecteds, which is supported by experimental evidence [START_REF] Litwin | Successful treatment of chronic hepatitis C with pegylated interferon in combination with ribavirin in a methadone maintenance treatment program[END_REF]. Most importantly, we assume that people who fail treatment (and return to the chronically infected pool) can be retreated with the same chance of success. This assumption is based on recent data showing that novel drugs (specifically Teleprevir) may have high success rates (50%) amongst non-responders with genotype 1, and the anticipation that other future drugs will have similar effects [START_REF] Mchutchison | PROVE3 Final Results and 1-Year Durability of SVR with Telaprevir-Based Regimen in Hepatitis C Genotype 1-Infected Patients with Prior Non-response[END_REF].

Details and explanation of the model

We use a system of ordinary differential equations to describe the transmission of HCV amongst active IDUs. We utilise a four compartment model, tracking susceptible, chronically infected, treated, and immune IDUs. Susceptible IDUs become infected through sharing of needles with an infected IDU.

About one quarter spontaneously clear the infection, and become susceptible or immune. The remaining three-quarters progress to chronic infection. Chronic infecteds can be treated, with a certain chance of success, and either fail treatment and return to the infection compartment, or clear the disease and become susceptible again or immune.

In our model, X denotes susceptible IDUs (including those who have cleared the infection), C denotes both chronically infected and acutely infected IDUs which will proceed to chronic infection, T denotes IDUs in treatment, Z denoting immune IDUs, τ is time in years, and where N=total population=X+Z+C+T. The equations describing the HCV transmission are,

dX dτ = θ -π(1 -δ + δξ) C N X + ωασT -μX (1) dC dτ = π(1 -δ) C N X -f (C) + ω(1 -α)T -μC (2) dT dτ = f (C) -ωT -μT (3) dZ dτ = πδξ C N X + ωα(1 -σ)T -μZ, ( 4
)
with initial conditions X(0) = X 0 , C(0) = C 0 , T (0) = 0, and Z(0) = 0.

Equation (1) represents the susceptible population, where new IDUs enter at a fixed rate θ. The second term in Equation 1 models the infection of a susceptible IDU, which is proportional to the number of susceptibles, the fraction of the population chronically infected, and the infection rate, π.

The acute infection spontaneously clears in a proportion δ, a fraction of which become immune at a proportion ξ. The remaining infected fraction which do not spontaneously clear, 1δ, progress to chronic infection. The third term in Equation 1 represents IDUs who exit treatment at a rate ω, with successful treatment proportion α, and who are part of the proportion not immune, σ. Due to the short duration of the acute stage, the number of infections caused by people with acute HCV who spontaneously clear or become immune is small, and we neglect it for model simplicity.

In each of the Equations ( 1)-(4), IDUs leave (due to death or ceasing injection) proportional to the rate μ.

Equation (2) models chronically infected IDUs. The first term represents those who enter from the susceptible pool, which is proportional to the number of susceptibles, the fraction of the population chronically infected, the infection rate, π, and the fraction who do not spontaneously clear the acute infection 1δ. The fraction of nonresponders to treatment, 1α, return from treatment proportional to rate ω.

The second term in Equation (2), f (C), represents the movement of infected IDUs into treatment.

In this paper, we examine two forms of the treatment recruitment function, which we describe in

Treatment functions: proportional versus fixed

We examine two alternative forms of the treatment function, treating either a fixed proportion of chronic infecteds per year (which we term the proportional treatment case) or treating a fixed number of chronic infecteds per year (which we call the fixed treatment case). The analysis of these two cases allows us to examine how changing the strategy of treatment delivery (attaining a specific coverage versus hitting a target number) alters the impact. In most infectious disease modelling, a proportional treatment term is used. However, we believe an annual fixed number of treatments term is more realistic, particularly in settings with limited resources (which effectively cap the number of possible treatments per year due to cost), or when targets are set to increase the total number treated per year, regardless of the proportion on infections. Hence, it is likely that if a treatment programme were to be initiated for active IDUs, the yearly targets would be number of IDUs treated.

Hence, for the proportional term, we examine a situation where the treatments initiated are a fixed proportion of chronic infecteds,

f (C) = φC, ( 5
)
where the annual recruitment rate of chronic infecteds on to treatment, φ, is in units per 1000 IDUs annually. For the fixed treatment term, we utilise a treatment function where f (C) is a fixed number of treatments per year,

f (C) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ Φ if Φ < C, C if 0 ≤ C < Φ. (6)
Hence, infected IDUs are recruited onto treatment at a fixed rate, Φ treatments initiated per 1000

IDUs per year. If the infected prevalence is driven below Φ people, all the infected IDUs are treated.

Nondimensionalisation

Employing the following substitutions: 1)-( 4) with ( 5) become

x = X N , c = C N , t = T N , z = Z N , then Equations (
dx dτ = θ N -π(1 -δ + δξ)cx + ωασt - θx N (7) dc dτ = π(1 -δ)cx -Φc + ω(1 -α)t - θc N (8) dt dτ = g(c) -ωt - θt N (9) dz dτ = πδξcx + ωα(1 -σ)t - θz N , ( 10
)
with the equation for N ,

dN dτ = θ -μN (11) with initial conditions x(0) = X 0 /N 0 , c(0) = C 0 /N 0 , t(0) = 0, z(0) = 0, and N 0 = X 0 + C 0 . Here x, c, t,
and z are the proportions of the total population of the susceptibles, infecteds, treated, and immune, respectively.

For the proportional treatment case,

g(c) = φc. ( 12
)
For the fixed treatment case,

g(c) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ Φ N if Φ N < c, c if 0 ≤ c < Φ N .
(13)

Discussion of parameter estimates

We obtain the model parameters from the relevant biological literature on injecting drug use and HCV treatment, a list of which is provided in Table 2. The exit rate is determined by the sum of the cessation of injecting rate and the IDU death rate, which are about 7.75% and 0.75% per year, respectively [START_REF] Sweeting | Estimating the prevalence of ex-injecting drug use in the population[END_REF][START_REF] Hickman | Hepatitis C virus (HCV) prevalence, and injecting risk behaviour in multiple sites in England in 2004[END_REF][START_REF] Nordt | Incidence of heroin use in Zurich, Switzerland: a treatment case register analysis[END_REF][START_REF] Hickman | Assessing IDU prevalence and health consequences (HCV, overdose and drug-related mortality) in a primary care trust: implications for public health action[END_REF]. The new IDU rate is calculated from the exit rate to retain 1000 IDUs in the population.

Treatment duration and success depends on the specific genotype of HCV being treated. In general, the treatment success is high for IDUs with genotype 2 or 3 (75-85%), and lower for genotype 1 (40-50%) [START_REF] Nice | Peginterferon alfa and ribavirin for the treatment of mild chronic hepatitis C. NIH[END_REF]NIH, 2002). The recommended duration of treatment for genotype 2/3 is 24 weeks for both responders or nonresponders, and for genotype 1 the duration is 48 weeks for responders and 12 weeks for nonresponders [START_REF] Nice | Peginterferon alfa and ribavirin for the treatment of mild chronic hepatitis C. NIH[END_REF]NIH, 2002). In the United Kingdom, about half the infections are genotype 1, with the remaining half 2 or 3 [START_REF] Nice | Peginterferon alfa and ribavirin for the treatment of mild chronic hepatitis C. NIH[END_REF]. Hence, we take an average between the genotype 1 and 2/3 parameters for the treatment success parameter (α), as well as the treatment duration (ω). [START_REF] Micallef | Spontaneous viral clearance following hepatitis C infection: a systematic review of longitudinal studies[END_REF] performed a meta-analysis and found that 26% of infections lead to spontaneous clearance, which we use for the parameter δ. The proportion of infections resulting in immunity after spontaneous clearance (ξ) and treatment (1σ) is not well known. We use a conservative estimate [START_REF] Mehta | Protection against persistence of hepatitis C[END_REF] for the immune proportion after spontaneous clearance (25%), and assume the same proportion after treatment.

In our simulations we vary the infection rate, π, in order to obtain projections at different endemic prevalences.

Analytical Results

Steady states and local stability

Proportional treatment

At steady state, N = θ μ . Setting the left hand side of Equations ( 7)-( 10) to zero (with g(c) defined by ( 12)) and solving for the equilibrium values we find two steady states. One is the trivial disease free steady state,

F P 1 = (x * 1 , c * 1 , t * 1 , z * 1 ) = (1, 0, 0, 0). ( 14
)
The second equilibrium is the infected endemic steady state, defined by

F P 2 = (x * 2 , c * 2 , t * 2 , z * 2 )
, where

t * 2 = φc * 2 ω + μ (15) x * 2 = μ(ω + μ) + ωασφc * 2 π(1 -δ + δξ)(ω + μ)c * 2 + μ(ω + μ) (16) z * 2 = πδξωασφ(c * 2 ) 2 + πδξμ(ω + μ)c * 2 μπ(1 -δ + δξ)c * 2 + μ 2 + ωα(1 -α)φc * 2 μ(ω + μ) (17) c * 2 = -μ(ω -μ) ω(1-α)φ ω+μ -φ -μ -π(1 -δ)μ(ω + μ) ω(1-α)φ ω+μ -φ -μ π(1 -δ + δξ)(ω + μ) + π(1 -δ)ωασφ . ( 18
)
By rearranging, the equation for c * 2 (Equation ( 18)) becomes,

c * 2 = Δ 1 φβ -Γ + Δ 2 β -Γ φ (19)
where

Δ 1 = μ 2 (ω + μ) -π(1 -δ)μ(ω + μ) Δ 2 = μ(ω + μ) -μω(1 -α) β = ω(1 -α)π(1 -δ + δξ) -π(1 -δ + δξ)(ω + μ) + π(1 -δ)ωασ Γ = μπ(1 -δ + δξ)(ω + μ).
Hence, increasing φ decreases both terms on the right hand side of Equation ( 19), leading to a monotonic decrease in c * .

At the trivial disease-free steady state, FP1, the eigenvalues (λ) are:

λ 1 = -μ λ 2 = -μ λ 3,4 = - 1 2 (ω + 2μ + φ -π(1 -δ))± 1 2 2ωφ + ω 2 + φ 2 + 2πδφ + 2πω + π 2 δ 2 + π 2 -2πφ -2π 2 δ -2πδω -4ωαφ.
The largest eigenvalue is λ 3 , and when λ 3 < 0 all the eigenvalues are negative and the disease-free state (FP1) is stable, and no infections will exist. If λ 3 > 0, the disease-free state is unstable and the endemic infected state (FP2) is stable. The numerical results indicate that if λ 3 > 0, the endemic infected state (FP2) is stable and a transcritical bifurcation of the chronic infected prevalence occurs with respect to the treatment rate, φ. This bifurcation is shown in Figure 1a. The steady state population fractions for varying infection rates (π) are shown in Figure 1b.

Figure 1c projects the steady-state relative reduction in prevalence for varying levels of treatment and across a range of untreated equilibrium infected prevalences. Even small treatment rates (less than 6% annually) can result in substantial reductions in infected prevalence (if not eradication) at low-medium untreated equilibrium prevalences (less than 40%). The high untreated prevalence scenario (60%) exhibits a reduced treatment effect at low treatment rates, but higher treatment rates (10-20%

annually) could result in prevalence reductions of 50-80%.

Fixed treatment

At steady state, N = θ μ . Setting the left hand side of Equations ( 7)-( 10) to zero (with g(c) defined by ( 13)) we can solve for the equilibrium values for the different treatment regions depending on Equation (13).

Case A, when g(c) = Φ N : In this region there are two fixed points,

F P 1b = μ(ωθ + μθ + ωασΦ) θ(πc * + (ω + μ)(1 -δ + δξ) + ωμ + μ 2 ) , c * + , Φμ θ(ω + μ) , πδξθc * + (ω + μ) + ωαΦπc * + (1 -δ + δξ) + ωασΦπc * + (1 -δ) + ωαΦμ(1 -σ) θ(πc * + (ω + μ)(1 -δ + δξ) + ωμ + μ 2 ) F P 2b = μ(ωθ + μθ + ωασΦ) θ(πc * -(ω + μ)(1 -δ + δξ) + ωμ + μ 2 ) , c * -, Φμ θ(ω + μ) , πδξθc * -(ω + μ) + ωαΦπc * -(1 -δ + δξ) + ωασΦπc * -(1 -δ) + ωαΦμ(1 -σ) θ(πc * -(ω + μ)(1 -δ + δξ) + ωμ + μ 2 )
where c * ± are the roots of a quadratic,

c * ± = -B± √ B 2 -4AD 2A
where

A = πθ(ω + μ)(1 -δ + δξ) (20) B = Φπ[(ωα + μ)(1 -δ + δξ) -ωασ(1 -δ)] -θ(μ + ω)[π(1 -δ) -μ] (21) D = Φμ(μ + ωα) (22)
With our parameters, A, D ≥ 0 for all treatment levels (even if Φ = 0). Hence, by Decartes' Rules of Signs, if B < 0 there are either 2 or 0 positive roots (equilibrium values for c), and if B > 0 there are either 2 or 0 negative roots. If B 2 > 4AD there are 2 real roots. If B 2 < 4AD there are no real roots and 2 imaginary roots.

If Φ = where << 1, then with our parameters B < 0 and B 2 > 4AD, so there are 2 positive real roots. As Φ increases, B remains negative, but when B 2 < 4AD the positive roots become imaginary.

With Φ very large (in fact, larger than our population of infecteds), B becomes positive, and B 2 > 4AD, which leads to 2 negative real roots. With our model however, a level of Φ this high would result in the case where g(c) = c, and the dynamics would proceed as in the following section.

Therefore, realistically, there is no biological situation resulting in negative real roots.

In the situation with 2 positive real roots, the local stability of the equilibrium values can be calculated by substituting FP1b and FP2b into the Jacobian matrix of Equations ( 7)-( 10) and determining the sign of the eigenvalues. Due to the complexity of the eigenvalues, we calculated them analytically using the mathematical software program Maple. For all values of our parameters where two positive roots exist, the numerical results indicate that FP2b (containing c * -) is an unstable steady state, and FP1b (containing c * + ) is stable. Therefore, there exists a saddle node bifurcation with the treatment parameter, Φ as shown in Figure 2a. Hence, with treatment, the only endemic steady state present is F P 1b.

Case B, when g(c) = c: This is a special case of the proportional treatment solution, found in Equations ( 14)-( 18) with Φ = 1. As before, the disease is always driven to eradication in this case. An outline of the possible cases with varying treatment levels is shown in Table 1.

Figure 2b shows the steady state projected relative reduction in prevalence for varying levels of treatment and across a range of untreated infected prevalences. The prevalence reduction curves are vertical near the prevalence where treatment causes eventual eradication, which is as a result of the saddle-node bifurcation present at the treatment threshold. As before, even low levels of treatment could result in large reductions in prevalence, although areas with high untreated prevalence are harder to cause reductions with treatment.

Treatment threshold for eradication

Proportional treatment

By setting λ 3 = 0 it is possible to solve for the critical value of treatment, φ c , which makes the endemic infected state unstable and the disease-free state stable, thus driving the disease to eradication. We find

φ c = πω(1 -δ) + πμ(1 -δ) -μω -μ 2 μ + ωα . ( 23
)
Figure 1d shows the treatment threshold, given by Equation ( 23), for varying untreated equilibrium infected prevalences and for both a mixed genotype scenario and all genotype 1 scenario. The threshold treatment level needed for clearance increases monotonically, and with an increasingly steep slope, for higher untreated prevalences. This indicates that with a high untreated prevalence, small increases in the equilibrium infection prevalence require large increases in treatment level to control. Additionally, the treatment threshold needed is higher in the all genotype 1 scenario, reflecting the reduced treatment success rates with this genotype. Despite this, settings with lower untreated prevalence (below 50%) could require comparatively low and achievable levels of treatment coverage (less than 20 % and 25%

annually in the mixed genotype and genotype 1 scenarios, respectively) to eradicate the disease, indicating that, all things being equal, treatment could be a valuable strategy for control.

Fixed treatment

From the stability analysis, the threshold for disease elimination occurs when B < 0 and B 2 = 4AD (B, A, and D defined in Equations ( 20)-( 22)), where the 2 positive real roots disappear. Setting B 2 = 4AD and solving for Φ c , we find that Φ c is itself the root of a quadratic, .

Φ c = -b± √ b 2 -4ad 2a where a = [π(ωα + μ)(1 -δ + δξ) -ωασπ(1 -δ)] 2 b = 2[π(ωα + μ)(1 -δ + δξ) -ωασπ(1 -δ)][-πθ(μ + ω)(1 -δ) + μθ(μ + ω)] -4μπθ(μ + ω)(1 -δ + δξ)(μ + ωα) d = [-πθ(μ + ω)(1 -δ) + μθ(μ + ω)] 2
Furthermore, by substituting in our parameter values, the following situation holds:

• When Φ < Φ c-, then with our parameters B < 0 and B 2 > 4AD, so there are 2 positive real roots in c * .

• When Φ c-< Φ < Φ c+ , then B 2 < 4AD, so there are no real roots in c * .

• When Φ c+ < Φ, then B > 0 and B 2 > 4AD, so there are 2 real negative roots in c * . However, as mentioned previously, for treatment values this high, the system would instead fall under the g(c) = c case, so no populations would actually go negative with our system.

Hence, the treatment level to clearance is the negative root, defined by

Φ c = Φ c-= -b - √ b 2 -4ac 2a . ( 24
)
Figure 2c shows the threshold treatment number needed for eradication (given by Equation ( 24)) for a range of untreated infected prevalences and both mixed genotype and all genotype 1 scenarios. As before, treatment levels needed are higher for the all genotype 1 scenario due to the reduced treatment efficacy. For the mixed treatment scenario and 40% prevalence, eradication could be possible by annually treating 9 infections per 1000 IDUs. At a high untreated prevalence of 60%, annually treating 29 infections per 1000 IDUs would be needed for clearance for the mixed genotype, raising to 40 treatments per 1000 IDU for the all genotype 1 scenario.

Numerical results

Numerical methods

Numerical simulations of Equations ( 7)-( 11) were performed using the MATLAB ODE solver, ode45.

Because HCV is not in the breakout epidemic phase in the majority of real-world settings, the simulations were run until steady state (τ = 600 years), and then treatment initiated. The parameter values are given in Table 2. For specific predictions at untreated equilibrium infected prevalences of 20%, 40%, and 60%, the values of π used were π = 0.1468, π = 0.2033, and π = 0.3307, respectively.

The calculation of the sensitivity coefficients, S, tells us what effect a change in the parameter has on a chosen variable. Specifically, it describes the factor relative change in the target variable relative to a factor change in a parameter. Positive coefficients indicate increasing the parameter increases our target variable. Conversely, negative coefficients mean increasing the parameter decreases the target variable. The larger the magnitude, the greater effect of the parameter change.

For our purposes, we are most interested in how sensitive the threshold treatment level needed for eradication is to the model parameters. Hence, the sensitivity coefficient for the threshold treatment level in the fixed treatment case, Φ c ), with respect to each parameter, p, is calculated as

S Φc,p = p Φ c dΦ c dp . ( 25
)
The sensitivity coefficients were calculated analytically using the software Maple.

Numerical simulations: Reductions in prevalence

Figure 3a presents the numerical simulations of the change in prevalence through time with an untreated equilibrium prevalence of 40% and treatment with the proportional treatment term at three treatment rates (2%, 4%, and 6% of chronically infected IDUs treated annually). Even though the treatment rates are below the threshold eradication level, they result in measurable reductions in infected prevalence within the first 10-20 years. Indeed, the majority of the prevalence reduction occurs within the first few decades, stabilizing at the new endemic prevalence within approximately 60 years. For example, treating only 2% of chronic infections each year may reduce the prevalence by over 15% (to below 34%) within 20 years. Treating at a rate of 4% annually reduces prevalence by nearly one-third in 20 years, and increasing treatment to 6% nearly halves the chronic infection prevalence in the same time span.

Figure 3b shows the timecourse of infection for the same untreated equilibrium prevalence (40%) but with the fixed treatment term. In this figure, the annual number of infections treated are 8, 16, and 24 per 1000 IDUs, which equates, at time=0, to treating 2%, 4%, and 6% of the chronically infected population. Note, however, that because the treatment number remains fixed, this results in treating an increased proportion of chronic infecteds as the prevalence is driven down through time. Figure 3b exhibits the abrupt collapse of the infected population (saddle-node bifurcation discussed in Section 4.1.2) with treatment above the treatment threshold, Φ c . Furthermore, the model shows that initiating 16 treatments per 1000 IDU annually can result in clearance within 60 years for a baseline prevalence of 40%. By contrast, treating at a proportional rate of 4% (shown in Figure 3a) does not result in eradication, only reducing prevalence from 40% to 23% by year 60.

Sensitivity analysis

The sensitivity coefficients of the parameters with respect to the threshold treatment level needed for eradication, φ c ,are shown for the proportional treatment model in Figure 4a for three untreated infected prevalence levels without treatment (20%, 40%, 60%). At all prevalence levels the threshold treatment level is most sensitive to the infection rate, π. The higher the prevalence, the more sensitive the treatment level is to the exit rate, μ, and the fraction of infecteds progressing to chronic infection, 1δ.

All prevalences are equally sensitive to the treatment success rate, α. Interestingly, the treatment threshold is not sensitive to the immunity parameters or treatment duration.

The sensitivity coefficients for the fixed treatment model are shown in Figure 4b for untreated infection prevalences 20%, 40%, and 60%. As in the previous treatment model, the threshold treatment level, Φ c , is most sensitive to the infection rate, π, with higher prevalences more sensitive to the exit rate μ, and the fraction of infecteds that progress to the chronic infection 1δ. Also as before, the treatment threshold is not sensitive to the immunity parameters or treatment duration.

Discussion and conclusions

In this paper we analyse a mathematical model of HCV transmission amongst active IDUs, examining the potential for antiviral treatment to reduce HCV transmission. Despite guidelines stating that current IDUs should not be excluded from obtaining treatment, very few are treated, with the risk of reinfection used as justification for withholding treatment. However, our model indicates that antiviral treatment could act as a prevention measure for the wider IDU community by reducing prevalence and therefore infection risk. Despite the possibility of reinfection after treatment, our models show that low levels of treatment could lead to large reductions in HCV prevalence, or even eradication.

The treatment rates predicted for substantial reductions in prevalence (10-20%), although much higher than those currently implemented, are feasible and achievable. For example, over half of current IDUs are in regular contact with services and hence accessible for treatment [START_REF] Hickman | Assessing IDU prevalence and health consequences (HCV, overdose and drug-related mortality) in a primary care trust: implications for public health action[END_REF].

Preliminary studies on treatment uptake indicate that willingness to enter treatment is high among active injectors [START_REF] Grebely | Treatment uptake and outcomes among current and former injection drug users receiving directly observed therapy within a multidisciplinary group model for the treatment of hepatitis C virus infection[END_REF]. Further, a recent review of treatment of active IDUs notes that treatment compliance and success rates are equal between active IDU and ex-or non-IDU [START_REF] Hellard | Hepatitis C Treatment for Injection Drug Users: A Review of the available evidence[END_REF].

Due to the logistical difficulties in initiating a comprehensive, nationwide treatment programme, we have examined two possible treatment scenarios: annually treating a proportion of infecteds, and annually treating a fixed number of infections. As there will likely be challenges to finding, testing/diagnosis, and recruiting IDUs onto treatment, it is possible that targeting a certain percentage of IDUs for treatment is the most reasonable. However, in situations where budget or service provider constraints restrict the absolute number of treatments per year, a fixed number treatment model could be more appropriate. Additionally, this type of model is more likely to fit the scenario in the early stages of a treatment initiative where targets are set by absolute numbers of treatment rather than a proportion of infected. Furthermore, even in the later stages of a treatment roll-out, advances in easy and inexpensive diagnostic techniques can lead to the ability to identify successively larger proportions of infected IDUs, thus enabling clinics to maintain a fixed number of treatments as the number of infections declines. Additionally, as a greater percentage of IDUs receive treatment, reluctance against treatment may decline, and as a larger number of IDUs become aware of the possibility, they will be more likely to access and opt for treatment. Hence, a positive 'snowball' effect can aid in the ability to remain in the fixed number treatment model, effectively accessing a higher proportion of infections as prevalence declines.

The two treatment models (proportional and fixed treatment number) exhibit different behaviour.

The proportional model has a transcritical bifurcation of the infection prevalence with respect to the treatment parameter. Increasing levels of treatment drives the steady state infected prevalence progressively lower to zero. This scenario is similar to that of 'constant effort' harvesting found in fisheries models, although in our case we would like to minimise both the population and the harvest.

By contrast, the fixed treatment number exhibits a saddle-node bifurcation of prevalence with respect to treatment, indicating that steady state prevalence may not be an accurate indicator of how close the disease is to eradication. This situation is similar to the situation present in 'constant yield harvesting' fisheries models, where small increases in harvesting over a threshold can cause rapid collapse of the fish populations, resulting in what is described as a catastrophe [START_REF] May | Thresholds and breakpoints in ecosystems with a multiplicity of stable states[END_REF]. In both constant yield fish harvesting and the fixed treatment number HCV model, steady state prevalence is not necessarily a good indicator of how close the population/infection is to collapse/eradication. In these situations, mathematical models can be an important way to track the efficacy and potential of control methods.

A sensitivity analysis of the models indicate that the treatment level needed for eradication is most sensitive to infection rate π, exit rate μ, the fraction of acute infections which spontaneously clear, δ, and treatment success rate α. Decreasing the infection rate, increasing the exit rate, or increasing the spontaneously clearing proportion will decrease the treatment level needed for eradication.

Interestingly, the treatment threshold is not sensitive to the immunity parameters, indicating that despite the uncertainty surrounding the presence and extent of immunity, it does not substantially alter the treatment thresholds needed for eradication.

It is important to note the limitations of this model. In this paper, we have constructed a deliberately simplified model to capture the basic transmission dynamics. In doing so, we have made a number of simplifying assumptions. First, we assume that the population of injectors mixes homogeneously, which is commonly invoked in infectious disease modelling but likely oversimplifies the complexities of sharing networks. Second, we assume that residence times in the compartments are exponentially distributed. Third, we assume that non-responders to treatment can be retreated with the same chance of success as treatment naive patients. As there are a number of new HCV treatments in development which show good success rates amongst non-responders, this is not an unrealistic assumption. However, limitations on re-treatment may play an important role when the disease is near to eradication, and we examine an extended model in another manuscript [START_REF] Martin | Can antiviral therapy for hepatitis C reduce the prevalence of HCV among injecting drug user populations? A modelling analysis of its prevention utility[END_REF].

Fourth, the question of immunity from spontaneous clearance or treatment is controversial. Our sensitivity analysis indicates that at the treatment levels examined the steady state prevalence is not very sensitive to the immunity parameters, but if circumstances manage to achieve higher rates, the role of immunity may be important. As information becomes available about immunity, the model can be adapted appropriately. Fifth, the use of a continuum model means that our model is not valid for very small population sizes where stochastic effects are likely to play a significant role in system dynamics.

As such, our model is best applied to urban areas with several thousand injectors (and hence the per 1000 IDU treatment rates should be scaled up accordingly). Finally, although there is no indication that antiviral resistance has developed, should this develop in the future it would be important to extend the model to include this aspect.

Overall, this model shows the potential of HCV antiviral therapy to reduce prevalence, even with reinfection and potential treatment failure. The model projects that achievable levels of treatment can result in large reductions in prevalence across a wide range of baseline prevalences, or possibly even eradication. It is likely that for sustained and substantial control of the HCV infection amongst IDUs a combination of prevention measures will need to be employed. Furthermore, the model characteristics might also be applicable to other diseases such as syphilis and chlamydia and could be explored.

Notably, there is emerging interest in the use of antiretroviral treatment for HIV as a means of reducing HIV transmission [START_REF] Blower | A Tale of Two Futures: HiV and Antiretroviral Therapy in San Francisco[END_REF][START_REF] Granich | Universal voluntary HiV testing with immediate antiretroviral therapy as a strategy for elimination of HiV transmission: a mathematical model[END_REF]. Our models indicate that using treatment for prevention of HCV spread amongst current IDUs should be explored as a realistic public health measure.

Acknowledgements

Grant Support

This manuscript was supported by the NCCRCD/NIHR CRDHB, MRC New Investigator Award, and the Scottish Government Hepatitis C Action Plan. Increasing π to 0.3, the infected fraction rapidly increases to about 60% prevalence. Increasing π further to 0.6 only increases the prevalence slightly, to about 70%. Parameters used are as in Table 2.

Figure captions

(c) Threshold level of treatment (φ c , vertical axis) needed to eradicate HCV with varying untreated equilibrium infected prevalence fractions (horizontal axis). At high prevalences, small increases in prevalence result in the need for substantially larger treatment rates. Here, π is varied to produce different untreated equilibrium infected prevalences, and φ c calculated from Equation 23. Parameters used are as in Table 2. (d) Relative reduction in prevalence (vertical axis) for varying untreated prevalences (horizontal axis) and treatment rates (φ, curves as noted in legend). Parameters used are as in Table 2. initiation of treatment at varying annual rates (φ= 0, 0.1, 0.2, and 0.4) with an untreated infection prevalence of 40%. Numerical simulations are of Equations ( 7)-( 12), with parameters as in Table 2, with π = 0.2033. (black). Sensitivity coefficients are calculated from Equations ( 7)-( 12) with parameters as in Table 2, with π = 0.1468 for 20% prevalence, π = 0.2033 for 40% prevalence, π = 0.3307 for 60% prevalence.

Table 1: Equilibrium values and stability for the fixed treatment model, Equations ( 7)-( 11) and Equation ( 13) with the parameters found in Table 2.

  represents IDUs currently in treatment. Infected IDUs enter treatment at the rate f (C) as discussed in Equation 2. Due to the reduction of viral loads during treatment, we assume that IDUs on treatment are not infectious during this stage. IDUs exit treatment proportional to the rate ω. Equation (4) tracks the immune population of IDUs. The first term represents those who become immune after spontaneous clearance of the acute infection. Susceptible IDUs can become infected proportional to the number of susceptibles, the fraction of the population chronically infected, and the infection rate, π. The acute infection spontaneously clears in the proportion δ, a fraction of which become immune at a proportion ξ. The second term represents those who become immune after successfully completing treatment. The IDUs exit treatment at the rate ω, with a proportion successful α, and with a further proportion who become immune 1σ.

  It is clear that a, d > 0. With our parameters, b < 0 and b 2 > 4ad, so by Decartes' Rules of Signs, there are two positive real roots, Φ c+ =

Figure 1 :

 1 Figure 1: (a) Transcritical bifurcation of equilibrium infected prevalence with respect to annual treatment rate, φ, for the proportional treatment model, Equations (7)-(12). (b) Steady state population fractions with varying infection rates (π) and no treatment. If π < 0.12, no outbreak occurs.

Figure 2 :

 2 Figure 2: (a)The saddle node bifurcation present of the equilibrium infected prevalence and annual treatment number, Φ for the fixed treatment model, Equations (7)-(11) and Equation (13). (b) Reduction in relative infection prevalence (%) with varying treatment numbers per year, Φ. (c) The threshold treatment number (Φ c , in treatments initiated per 1000 IDUs annually) needed for eradication at varying untreated equilibrium prevalences.

Figure 3 :

 3 Figure 3: (a) Proportional treatment model numerical simulations of prevalence through time after

Figure 4 :

 4 Figure 4: Sensitivity coefficients of the treatment threshold needed for eradication with respect to the model parameters for the (a) proportional treatment model and (b) fixed treatment model. Sensitivity coefficients are shown for baseline prevalences of 20% (light gray), 40% (dark gray), and 60%

  Figure1

  Figure 1

  Figure 2

  Figure 3

Figure 4

 4 Figure 4

Table 2 :

 2 Parameter values used in the numerical simulations. a Average of the genotype 1 cure rate (α 1 = 0.45) and the genotype 2/3 cure rate (α 2/3 = 0.8). b Exit rate calculated from the average of the genotype 1 treatment length for responders and nonresponders: (α 1 × 48 + (1α 1 ) × 12) weeks and the genotype 2 treatment length, 24 weeks. c We assume the same proportion become immune after treatment as after spontaneous clearance. d Based on a cessation rate of 7.75% per year, and an IDU death rate of 0.75%.

Treatment

Steady-States Stability

Table 2