Differential characteristics of primary infection and re-infection can cause backward bifurcation in HCV transmission dynamics

Mathematical Biosciences 263 (2015) 51–69

Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Differential characteristics of primary infection and re-infection can cause backward bifurcation in HCV transmission dynamics F. Nazari a , A.B. Gumel a,b,∗, E.H. Elbasha c a

Simon A. Levin Mathematical, Computational and Modeling Sciences Center, Arizona State University, Tempe, AZ 85287-1904, USA School of Mathematical and Natural Sciences, Arizona State University, Phoenix, AZ 85069-7100, USA c Merk Research Laboratories, UG1C-60, P.O. Box 1000, North Wales, PA 19454-1099, USA b

a r t i c l e

i n f o

Article history: Received 25 July 2014 Revised 4 February 2015 Accepted 5 February 2015 Available online 14 February 2015 Keywords: HCV Backward bifurcation Reproduction number Stability

a b s t r a c t Backward bifurcation, a phenomenon typically characterized by the co-existence of multiple stable equilibria when the associated reproduction number of the model is less than unity, has been observed in numerous disease transmission models. This study establishes, for the first time, the presence of this phenomenon in the transmission dynamics of hepatitis C virus (HCV) within an IDU population. It is shown that the phenomenon does not exist under four scenarios, namely (i) in the absence of re-infection, (ii) in the absence of differential characteristics of HCV infection (with respect to infectivity, progression, treatment and recovery) between re-infected individuals and primary-infected individuals, (iii) when re-infected and treated individuals do not transmit HCV infection and (iv) when the average infectivity-adjusted duration of re-infection is less than that of primary infection. This study identifies, using sensitivity analysis, five parameters of the model that have the most influence on the disease transmission dynamics, namely: effective contact rate, progression rate from acute to chronic infection, recovery rate from acute infection, natural death rate and the relative infectiousness of chronically-infected individuals. Numerical simulations of the model show that the reinfection of recovered individuals has marginal effect on the HCV burden (as measured in terms of the cumulative incidence and the prevalence of the disease) in the IDU community. Furthermore, treatment of infected IDUs, even for small rate (such as 4%), offers significant impact on curtailing HCV spread in the community. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Hepatitis C, a blood-borne viral infectious disease caused by the single-strained RNA Hepatitis C virus (HCV) [2,10], continues to pose major public health challenges globally. The World Health Organization (WHO) recently estimated the HCV prevalence to be between 2 and 3% (i.e., 130–170 million people currently live with HCV infection globally) [31,38,49,50]. While countries in Africa and Asia have the highest reported prevalence, industrialized countries in North America, Northern and Western Europe and Australia have lower prevalence (Germany (0.6%), Canada (0.8%), France (1.1%), and Australia (1.1%) have relatively low rates of HCV seroprevalence, with the USA (1.8%), Japan (1.5–2.3%), and Italy (2.2%) having slightly higher seroprevalence rates) [38,49]. The primary mode of HCV transmission is through blood contact [1,50]. In particular, there are three main age-specific transmission ∗ Corresponding author at: Arizona State University, SAL-MCMSC, 1551 S. Forest, Avenue, Tempe, Arizona, USA. Tel.: +1 480 7272 690. E-mail address: [email protected] (A.B. Gumel).

http://dx.doi.org/10.1016/j.mbs.2015.02.002 0025-5564/© 2015 Elsevier Inc. All rights reserved.

patterns [1,42]. The first is the 30–49 year age (middle age) group in developed countries (such as the United Kingdom and USA) [31,42]. For this age group, injecting drug use (IDU) is the major cause of HCV infection (via needle and syringe-sharing; with over 80% of new cases attributed to injecting drugs use) [31,42]. The second transmission pattern is for the elderly (as is the case in Japan [26,35,42,48]). The third transmission pattern entails all age groups. For these two later categories, the causes (or risk factors) of HCV infection include unsafe therapeutic injections (performed by both healthcare professionals and non-professionals) and blood transfusion from unscreened donors. In addition to the aforementioned HCV risk factors, other factors, such as exposure to blood by the healthcare workers (mostly through contact with contaminated needles), mother-to-child transmission, sex with an infected partner, sex with multiple partners and other healthcare-related procedures, further contribute to HCV transmission [38]. The common symptoms of HCV infection include jaundice, dark urine, fatigue, nausea, vomiting, and abdominal pain [13,25]. While majority of patients with acute HCV will progress to chronic infection [2,10,13], about 25% of cases clear the virus and build

52

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

natural immunity against re-infection [18]. The mean incubation period for acute HCV infection is 7 weeks [24]. Unfortunately, up to 90% of HCV-infected individuals (acute or chronic) may not to be aware of their infection status (i.e., they are asymptomaticallyinfected). Consequently, if undetected and untreated, about 7–18% of these (asymptomatically-infected) individuals will progress to develop liver disease, such as liver fibrosis, cirrhosis, hepatocellular carcinoma, within 20 years (and about 5%–7% of these patients may ultimately die of HCV) [2,10,13,15,25,31,49,50]. HCV-infected individuals can be treated using a combination therapy with pegylated interferon and ribavirin (having a response rate of 40%–80%) [13,17] (furthermore, several new anti-HCV drugs have recently been approved and/or are undergoing various stages of clinical trials [29]). Although there is currently no safe and effective vaccine for use against HCV infection in humans, efforts are underway to develop one [13,15]. Another intervention strategy for controlling HCV transmission among IDUs is increasing the access to unused syringes and needles, aimed at reducing the frequency of sharing/unsafe injection needles (it should, however, be mentioned that although this approach may have a positive effect on reducing HCV transmission, there is no evidence for substantial reductions in HCV prevalence) [32,41,43]. Several mathematical and statistical models have been developed and used to gain insight into the transmission dynamics of HCV in an IDU population [8,9,27,31,40,45,46,52]. Corson et al. [9] developed a deterministic compartmental mathematical model for the spread of HCV in an IDU population that has been separated into two groups (naive and experienced) based on the time since the onset of injection (and includes measures that allow for the prevention of HCV infection). Sutton et al. [40] used statistical modelling to estimate the force of infection of HCV and hepatitis B virus in England and Wales (using saliva sample of IDUs) for 1998–2003. Vickerman et al. [45] used a deterministic model to simulate the transmission of HCV in IDUs in London, England and assessed the impact of intervention measures that reduced syringe and needle sharing in some of the targeted IDU populations. Elbasha [13] introduced the effect of the re-infection of recovered primary infected individuals (and associated differential infection characteristics between re-infected individuals and primary infected individuals) on HCV transmission dynamics via the use of a deterministic model. Furthermore, Elbasha [13] provided a rigorous qualitative analysis of a special case of the model in [13] (where re-infected individuals behave in the same manner as primary infected individuals, with respect to disease infectivity, recovery, progression and treatment). The main purpose of this study is to provide qualitative insight into the effect of treatment on the transmission dynamics of HCV in an IDU population. To achieve this objective, the treatment model developed in [13] will be considered (and fully analyzed, unlike the special case considered in [13]). Furthermore, the effect of uncertainties on the associated parameters of the model (on the overall simulation results obtained) will be assessed using Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCC) [3,30,36]. The paper is organized as follows. A treatment model for HCV transmission dynamics in an IDU population is formulated in Section 2. The analogous model in the absence of treatment is analyzed in Section 3. The full treatment model is analyzed in Section 4. Uncertainty and sensitivity analyses of the treatment model are also carried out. Numerical simulation results are also presented. 2. Model formulation The HCV transmission model to be considered in this section is developed in [13], and is formulated by splitting the total IDU population at time t, denoted by N(t), into mutually-exclusive compartments of susceptible, S(t), acutely-infected, I (t), chronically-infected,

P (t), treated chronically infected, T (t), recovered with partial immunity, R(t), acutely re-infected, V (t), untreated chronically re-infected, W (t) and treated chronically re-infected, Q (t), individuals, so that

N(t) = S(t) + I (t) + P (t) + R(t) + V (t) + W (t) + T (t) + Q (t). The population of susceptible individuals (S) is increased by the recruitment of new IDUs into the IDU population (at a rate ). It is further increased by the loss of infection-acquired immunity of recovered individuals (at a per capita rate γ ). It is decreased by infection, following effective contacts with infected individuals, at a rate λ, given by

λ=

β(I + π P + υ V + ωπ W + π χT T + π χQ Q ) N

,

(2.1)

In (2.1), β is the effective contact rate, π , υ , ω, χT and χQ are modification parameters accounting for the relative infectivity of chronically infected, acutely re-infected, chronically re-infected, treated infected and treated re-infected individuals, in comparison to acutely infected individuals, respectively. This population is further decreased by natural death (at a rate μ; this rate is assumed, for mathematical convenience, to be the same for all of the epidemiological compartments). Thus,

dS =  + γ R − λS − μS. dt The population of acutely-infected individuals (I) is increased by the infection of susceptible individuals (at the rate λ). It is decreased by recovery (at a rate σ ), progression to chronic stage (at a rate ε ) or natural death. Thus,

dI = λS − (σ + ε + μ)I. dt The population of chronically-infected individuals (P) is generated by progression of acutely-infected individuals to the chronicallyinfection (at the rate ε ). It is also increased by treatment failure of chronically-infected individuals (at a rate ρ ). This population is decreased by treatment (at a rate τ ), recovery (at a rate δ ) or natural death. Hence,

dP = ε I + ρ T − (δ + τ + μ)P. dt The population of recovered individuals (R) is generated by the recovery of acutely-infected individuals (at the rate σ ), chronically-infected individuals (at the rate δ ), acutely-re-infected individuals (at a rate ασ , where α > 1 is the modification parameter accounting for the assumption that acutely-re-infected individuals recover at a faster rate in comparison to acutely-infected individuals), and chronicallyre-infected individuals (at a rate ηδ where, η > 1 is the modification parameter accounting for the assumption that chronically-re-infected individuals recover at a faster rate in comparison to chronicallyinfected individuals) [13]. It is also increased by the successful treatment of chronically-infected and chronically-re-infected individuals (at a rate θ ). This population is decreased by infection (at a reduced rate ψλ, where 0 < ψ < 1 accounts for the assumption that recovered individuals acquire HCV infection at a rate lower than susceptible individuals). It is further decreased by the loss of infection-acquired immunity (at the rate γ ) and natural death. Thus,

dR = σ I + δ P + ασ V + ηδ W + θ T + θ Q − ψλR − (γ + μ)R. dt The population of acutely re-infected individuals (V) is increased by the re-infection of recovered individuals (at the rate ψλ). It is decreased by progression to the chronic re-infection stage (at a rate κε, where 0 < κ < 1 is the modification parameter accounting for the assumption that acutely-re-infected individuals progress to the chronically-re-infection stage (W) at a slower rate in comparison to

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

53

acutely-infected individuals), recovery (at the rate ασ ) or natural death. Hence,

dV = ψλR − (ασ + κε + μ)V. dt The population of chronically-re-infected individuals (W) is increased by the progression of acutely-reinfected individuals (at the rate κε ) and by the failure of treatment in chronically-re-infected individuals (at a rate ζ ). It diminishes by recovery (at the rate ηδ ), treatment of chronically-re-infected individuals (at a rate φ ) or natural death. Thus,

dW = κε V + ζ Q − (ηδ + φ + μ)W. dt The population of chronically-infected treated individuals (T) is increased by the treatment of the chronically-infected individuals (at the rate τ ). It is decreased by treatment failure or successful treatment of the chronically-infected individuals (at the rates ρ and θ , respectively) or by natural death. Thus,

dT = τ P − (ρ + θ + μ)T. dt Similarly, the population of chronically-re-infected individuals (Q) is generated by the successful treatment of chronically-re-infected individuals (at the rate φ ). It diminishes by the treatment failure or successful treatment of the chronically-re-infected individuals (at the rates ζ or θ , respectively). It is further decreased by the natural death. Thus,

Based on the above assumptions and derivations, the treatment model for HCV transmission dynamics within an IDU population is given by the following deterministic system of non-linear differential equations [13]:

=  + γ R − λS − μS, =





δ + τ + μ P,

(2.2) 



=

σ I + δ P + ασ V + ηδ W + θ T + θ Q − ψλR − γ + μ R,

=

ψλR − (ασ + κε + μ) V,

=

κε V + ζ Q − ηδ + φ + μ W,

=

τ P − ρ + θ + μ T,

=

where,

λ=

λS − (σ + ε + μ) I,

= εI + ρ T −













φ W − ζ + θ + μ Q, 

β I + π P + υ V + π ωW + π χT T + π χQ Q N

 ,

(2.3)

is the infection rate. A flow diagram of the model is given in Fig. 1, and the associated variables and parameters are tabulated in Table 1. It is worth stating that Elbasha [13] studied a special case of the HCV treatment model (2.2), where re-infection plays the same role as primary infection, with respect to infectivity, progression, treatment and recovery (i.e., the model (2.2) with κ = α = η = ω = υ = 1, χQ = χT , ζ = ρ and φ = τ ), given by (the reduced model)

dS =  + γ R − λS − μS, dt dI = λS − (μ + σ + ε) I, dt

Table 1 Description of variables and parameters of the treatment model (2.2) [13]. Symbol

Description

Variables S(t) I (t) P (t) R(t) V (t) W (t) T (t) Q (t)

Population of susceptible individuals Population of acutely-infected individuals Population of chronically-infected individuals Population of recovered individuals with partial immunity Population of acutely-reinfected individuals Population of chronically-reinfected individuals Population of chronically-infected treated individuals Population of chronically-reinfected treated individuals

Parameters

dQ = φ W − (ζ + θ + μ)Q. dt

dS dt dI dt dP dt dR dt dV dt dW dt dT dt dQ dt

Fig. 1. Schematic diagram of the treatment model (2.2) [13].

 μ β σ δ ε ψ α η κ

Recruitment rate Natural death rate Contact rate Recovery rate from acute infection Recovery rate from chronic infection Rate of progression from acute infection to chronic infection Relative susceptibility of recovered individuals Relative rate of recovery from acute re-infection Relative rate of recovery from chronic re-infection Relative rate of progression from acute re-infection to chronic re-infection Treatment rate of chronically-infected individuals Treatment rate of chronically-re-infected individuals Treatment cure rate Treatment failure rate of chronically-infected individuals Treatment failure rate of chronically-re-infected individuals Rate of waning immunity Relative infectivity of chronically-infected individuals Relative infectivity of acutely-re-infected individuals Relative infectivity of chronically-re-infected individuals Relative infectivity of treated infected individuals Relative infectivity of treated re-infected individuals

τ φ θ ρ ζ γ π υ ω χT χQ

dP dt dR dt dV dt dW dt dT dt dQ dt

= εI + ρ T −





μ + δ + τ P,

(2.4) 



=

σ I + δ P + σ V + δ W + θ T + θ Q − ψλR − μ + γ R,

=

ψλR − (μ + σ + ε) V,

= εV + ρ Q −





μ + δ + τ W, 



=

τ P − μ + ρ + θ T,

=

φ W − μ + ρ + θ Q,





where, now,

λ=

β(I + π P + V + π W + π χT T + π χT Q ) N

.

The above simplifications, made in [13], allow for the change of variables, I¯ = I + V, P¯ = W + P, T¯ = T + Q (and, for consistency, R¯ = R

54

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

and S¯ = S), so that the reduced model (2.4) can be re-written as [13]:

dS¯ dt dI¯ dt dP¯ dt dR¯ dt dT¯ dt

¯ S¯ − μS¯ + γ R, ¯ = −λ ¯ S¯ + ψ λ ¯ R¯ − (μ + σ + ε) I¯, =λ = ε I¯ + ρ T¯ −

Theorem 2.2. The closed set



¯ μ + δ + τ P,

(2.5) 



=

¯ σ I¯ + δ P¯ + θ T¯ − ψ λ¯ R¯ − μ + γ R,

=

¯ τ P¯ − μ + ρ + θ T,

with,

λ¯ =





β I¯ + π P¯ + π χT T

¯ and N¯ = S¯ + I¯ + P¯ + R¯ + T.





β π ε(θ + μ + ρ + τ χT ) 1+ , ε+μ+σ (δ + μ)(θ + μ + ρ) + (θ + μ)τ

is less than unity. It is further shown, for the special case of the treatment model (2.5) in the absence of re-infection (ψ = 0), that the unique endemic equilibrium (EEP) of the reduced model (2.5) is globally asymptotically stable whenever it exists [13]. 2.1. Basic properties Theorem 2.1. Let the initial data for the model (2.2) be S(0) > 0, I (0) > 0, P (0) > 0, R(0) > 0, V (0) > 0, W (0) > 0, T (0) > 0andQ (0) > 0. Then, the solutions

(S(t), I (t), P(t), R(t), V (t), W (t), T (t), Q (t)) of the model (2.2), with positive initial data, will remain positive for all time t > 0. Proof. Let

t1 = sup{t > 0 : S(t) > 0, I (t) > 0, P (t) > 0, R(t) > 0, V (t) > 0, W (t) > 0, T (t) > 0, Q (t) > 0} > 0. It follows from the first equation of the model (2.2) that

dS =  − λS − μS + γ R ≥  − λS − μS, dt which can be written as,

 

  

 t t d λ(τ )dτ ≥  exp μt + λ(τ )dτ . S(t)exp μt + dt 0 0

S(t1 )exp ≥

t1 0



μt1 + 



t1 0

 exp μy +

so that,





λ(τ )dτ − S(0)



S(t1 ) ≥ S(0)exp −μt1 −



y 0



  + exp −μt1 −

t1 0





μ

dN =  − μN, dt

It was shown in [13] that the disease-free equilibrium (DFE) of the model (2.5) is globally-asymptotically stable whenever the associated reproduction number, given by,

Thus,



   S, I, P, R, V, W, T, Q ∈ R8+ : N ≤

Proof. Adding the equations of the model (2.2) gives



N¯

¯c = R

DT =

is positively-invariant and attracts all positive solutions of the model (2.2).





Similarly, it can be shown that I (t) > 0, P (t) > 0, R(t) > 0, V (t) > 0, W (t) > 0, T (t) > 0, and Q (t) > 0 for all time t > 0. Hence, all solutions of the model (2.2) remain positive for all non-negative initial conditions, as required.

 exp μy +

λ(τ )dτ

t1 0 t1

0



λ(τ )dτ λ(τ )dτ y

0



dy,  

λ(τ )dτ



dy > 0.

(2.6)

from which it is clear that dN is negative if N(t) >  μ . It follows from dt the solution of equation (2.6)

N(t) =





  − μt + N(0) − e , μ μ

 that if N(0) <  μ , then N (t) ≤ μ for all t > 0. That is, all orbits of the model (2.2) with initial conditions in DT remain in DT for all t > 0. Thus, the region DT is positively-invariant. Furthermore, if N(0) >  μ , then either the solution enters DT in finite time or N (t) ap-

proaches to  μ as t → ∞. Hence, the region DT attracts all solutions in R8+ .

Since the region DT is positively-invariant, the unique solution of the treatment-free model (2.2) exists and depends continuously on the initial data of the model (hence, it is sufficient to study its asymptotic dynamics in the region DT [23]). As stated in Section 1, a combination therapy with pegylated interferon and ribavirin are known to be quite effective against HCV infection [31]. Although the prevailing opinion among public health and medical practitioners in the US and UK, prior to 2002, was against treating active IDUs [31], IDUs are now receiving anti-HCV treatment (owing to the increasing evidence showing that IDUs exhibit a similar (non-IDUs) response to anti-HCV treatment [32]). Nevertheless, despite this fact, and the high number of IDUs infected, very few active IDUs have been (and are being) treated [31,34,47]. Therefore, it is constructive to study the transmission dynamics of HCV among IDUs in the absence of treatment. This is done below. 3. Analysis of model without treatment The special case of the treatment model (2.2) in the absence of treatment (referred to as treatment-free model), obtained by T (t) = Q (t) = χT = χQ = τ = φ = ζ = ρ = θ = 0 in (2.2), is given by the following deterministic system of non-linear differential equations [13]:

dS dt dI dt dP dt dR dt dV dt dW dt

=  + γ R − λwt S − μS, =

λwt S − (σ + ε + μ)I,

= ε I − (δ + μ)P, =

σ I + δ P + ασ V + ηδ W − ψλwt R − (γ + μ)R,

=

ψλwt R − (ασ + κε + μ)V,

=

κε V − (ηδ + μ)W,

(3.1)

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

where,

λwt =

3.2.2. Existence of EEP In this section, the possible existence of an endemic equilibrium (that is, an equilibrium of the treatment-free model (3.1) when the infected components are non-zero will be explored). Let,

β(I + π P + υ V + ωπ W )

and Nwt Nwt (t) = S(t) + I (t) + P (t) + R(t) + V (t) + W (t).

It is worth noting that the treatment-free model (3.1) reduces to an SIR model in the absence of re-infection of recovered individuals (ψ = 0) and loss of infection-acquired immunity (γ = 0). Furthermore, it can be shown that it reduces to an SIS model if recovered individuals acquire HCV infection at the same rate as whollysusceptible individuals (ψ = 1). The treatment-free model (3.1) will now be analyzed to gain insight into its qualitative features. 3.1. Basic properties The following results can be proved using the approaches in Section 2.1. Theorem 3.1. Let the initial data for the treatment-free model (3.1) be S(0) > 0, I (0) > 0, P (0) > 0, R(0) > 0, V (0) > 0 and W (0) > 0. Then, the solutions

(S(t), I (t), P(t), R(t), V (t), W (t)) of the treatment-free model (3.1), with positive initial data, will remain positive for all time t > 0. Theorem 3.2. The closed set Dwt =

 

  S, I, P, R, V, W ∈ R6+ : Nwt ≤



3.2. Existence and stability of equilibria 3.2.1. Local asymptotic stability of DFE The DFE of the treatment-free model (3.1), obtained by setting the right-hand side of the equations in the model to zero, is given by







 , 0, 0, 0, 0, 0 . μ

(3.2)

Using the next generation operator method [11,44], the matrices Fwt (of the new infection terms) and Hwt (of the transition terms) associated with the model (3.1) are given, respectively, by

⎛

β π β υβ π ωβ

⎜ ⎜0

Fwt = ⎜ ⎜0 ⎝

0

0

0

0

0

0

0

0

0

0

⎞ ⎟ ⎟ ⎟, ⎟ ⎠

⎛

G1

0

⎜ ⎜−ε G2

0 0

Hwt = ⎜ ⎜0 ⎝

0

G3

0

0

−κε

0

⎞

⎟ 0⎟ ⎟, 0⎟ ⎠ G4

where, G1 = μ + σ + ε , G2 = μ + δ , G3 = μ + ασ + κε , and G4 = μ + ηδ . It follows, from Theorem 2 in [44], that the basic reproduction number of the treatment-free model (3.1), defined by −1 ) (where ρwt is the spectral radius of the next R0 = ρwt (Fwt Hwt −1 ), is given by generation matrix Fwt Hwt



(3.4)

be an arbitrary endemic equilibrium of the treatment-free model (3.1), where S∗∗ , I∗∗ , P ∗∗ , R∗∗ , V ∗∗ , and W ∗∗ are obtained from setting the right-hand-sides of the equations in the model (3.1) to zero, given by ∗∗ γ R∗∗ +  ∗∗ λ∗∗ wt S ,I = , σ + ε + μ λ∗∗ + μ wt ε I∗∗ σ I∗∗ + δ P∗∗ + ασ V ∗∗ + ηδ W ∗∗ P ∗∗ = , R∗∗ = , δ+μ μ + γ + ψλ∗∗ wt ∗∗ ψλ∗∗ κε V ∗∗ wt R V ∗∗ = , W ∗∗ = . ασ + κε + μ ηδ + μ

S∗∗ =

(3.5)

Furthermore, let

λ∗∗ wt =

βμ(I∗∗ + π P∗∗ + υ V ∗∗ + π ωW ∗∗ ) , 

(3.6)

∗ =  ) be the (where Nwt (t) is now replaced by its limiting value, Nwt μ force of infection of the treatment-free model (3.1) at an endemic steady-state. Substituting the equations in (3.5) into (3.6), and simplifying, gives the following quadratic equation (in terms of λ∗∗ wt )

(3.7)

where,

is positively-invariant and attracts all positive solutions of the model (3.1).



E1wt = (S∗∗ , I∗∗ , P ∗∗ , R∗∗ , V ∗∗ , W ∗∗ ),

2 ∗∗ c2 (λ∗∗ wt ) + c1 λwt + c0 = 0,

μ

E0wt = S∗ , I∗ , P ∗ , R∗ , V ∗ , W ∗ =

55



β επ 1+ . R0 = ε+μ+σ δ+μ

(3.3)

The result below follows from Theorem 2 of [44]. Theorem 3.3. The DFE, E0wt , of the treatment-free model (3.1), given by (3.2), is locally-asymptotically stable (LAS) if R0 < 1, and unstable if R0 > 1. The epidemiological implication of Theorem 3.3 is that HCV can be effectively controlled in the community (when R0 < 1) if the initial sizes of the sub-populations of the model (3.1) are in the basin of attraction of the DFE. The threshold quantity, R0 , represents the average number of secondary infections that one HCV-infected individual can generate if introduced into a completely-susceptible IDU population [23,44].

ψ(μ + δ)(μ + δη + κε)(μ + σ + ε), c1 = γ (δη + μ)(ασ + μ + κε)(δ + μ + ε) + μψ(δ + μ)(δη + μ + κε)(ε + μ + σ )(1 − R0 ) (3.8) + (μ + δ)(μ + σ + ε)(μ + ασ + κε)(μ + ηδ) − βψ [σ (μ + δ) + δε ] [υ(μ + ηδ) + ωπ κε ] , c0 = (γ + μ)(δ + μ)(δη + μ)(ε + μ + σ )(μ + ασ + κε)(1 − R0 ),

c2 =

The components of the endemic equilibrium are then obtained by solving for λ∗∗ wt from the quadratic (3.7), and substituting the positive values of λ∗∗ wt into the expressions in (3.5). Furthermore, it follows from (3.8) that the coefficient c2 , of the quadratic (3.7), is always positive, and c0 is positive (negative) if R0 is less (greater) than unity. Hence, it follows from (3.8) that the quadratic (3.7) has a unique positive equilibrium (an endemic equilibrium) whenever R0 > 1. These results are summarized below. Theorem 3.4. The treatment-free model (3.1) has: (i) a unique endemic equilibrium if c0 < 0 ⇔ R0 > 1; (ii) a unique endemic equilibrium if c1 < 0 and c0 = 0 or c12 − 4c0 c2 = 0; (iii) two endemic equilibria if c0 > 0, c1 < 0 and c12 − 4c0 c2 > 0; (iv) no endemic equilibrium otherwise. Item (iii) of Theorem 3.4 suggests the possibility of backward bifurcation [7,14,16,19–21,23,33,37,51] in the treatment-free model (3.1). The phenomenon of backward bifurcation is characterized by the co-existence of a stable DFE and a stable EEP when the associated reproduction number of the model (R0 ) is less than unity. The epidemiological consequence of backward bifurcation is that disease control (when R0 < 1) is dependent on the initial sizes of the sub-populations of the model (see, for example, [37]). Hence, the presence of backward bifurcation in the transmission dynamics of a disease makes its effective control (or elimination) difficult. Consequently, it is instructive to explore the possibility of backward bifurcation in the treatment-free model (3.1). Before doing so, it is worth checking for the existence of an EEP of the model (3.1) when R0 ≤ 1 (which is a signature for a backward bifurcation in disease transmission models [7,14,16,20,23,33,37,51]). This is done below.

56

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

σ (μ + δ) + δε ≤ (μ + δ)(μ + σ + ε), it can be shown from (3.8) that (noting that ψ ≤ 1, υ ≤ 1 and ω ≤ 1)

γ (δη + μ)(ασ + μ + κε)(δ + μ + ε)

c1 ≥

+ μψ(δ + μ)(δη + μ + κε)(ε + μ + σ )(1 − R0 )

(3.10)

+ (μ + δ)(μ + σ + ε){(μ + ασ + κε)(μ + ηδ) − βψ [υ(μ + ηδ) + ωπ κε ]}, γ (δη + μ)(ασ + μ + κε)(δ + μ + ε) + μψ(δ + μ)(δη + μ + κε)(ε + μ + σ )(1 − R0 ) + (μ + δ)(μ + σ + ε){(μ + ασ + κε)(μ + ηδ) − β [(μ + ηδ) + π κε ]}.

≥

It should, first of all, be noted from (3.3) that the inequality R0 < 1 implies that

β< Fig. 2. Backward bifurcation diagram for the treatment-free model (3.1), showing the prevalence as a function of the basic reproduction number (R0 ). Parameter values used are:  = 40 640, β ∗ = 3.324462329, μ = 0.091, σ = 0.255, δ = 0.0025, ε = 3.5, ψ = 0.70, α = 3.3, η = 4.050, κ = 0.0827, γ = 0.0225, π = 0.019, υ = 0.89999 and ω = 0.89999 (so that, a = 0.0001180441877 > 0; R0 = 1).

3.3. Backward bifurcation Theorem 3.5. The treatment-free model (3.1) undergoes a backward bifurcation at R0 = 1 whenever the Inequality (A.5), given in Appendix A, holds. The proof of Theorem 3.5, based on using center manifold theory [4–6,12,44], is given in Appendix A. A schematic description of the backward bifurcation phenomenon of the model (3.1) is depicted in Fig. 2. It is worth stating that, to the authors’ knowledge, this is the first time the phenomenon of backward bifurcation has been established in the transmission dynamics of HCV in a population. 3.3.1. Non-existence of backward bifurcation First of all, it should be mentioned that setting the re-infection parameter ψ , to zero reduces the treatment-free model (3.1) to an SIRS model, which is known not to undergo backward bifurcation [20]. That is, as in the case of other disease such as TB [20], the reinfection of recovered individuals causes backward bifurcation in HCV transmission dynamics. Case (i). Effect of relative rate of progression of acute re-infected individuals (κ ). Consider the treatment-free model (3.1) for the case when acutelyre-infected individuals (V) progress to chronic re-infection stage (W) at the same rate as acutely-infected (I) individuals (i.e., κ = 1). For this case (of the model (3.1) with κ = 1), the coefficients c0 , c1 and c2 , of the quadratic (3.7), reduce to:

ψ(μ + δ)(μ + δη + ε)(μ + σ + ε), (3.9) c1 = γ (δη + μ)(ασ + μ + ε)(δ + μ + ε) + μψ(δ + μ)(δη + μ + ε)(ε + μ + σ )(1 − R0 ) + (μ + δ)(μ + σ + ε)(μ + ασ + ε)(μ + ηδ) − βψ [σ (μ + δ) + δε ] [υ(μ + ηδ) + ωπ ε ] , c0 = (γ + μ)(δ + μ)(δη + μ)(ε + μ + σ )(μ + ασ + ε)(1 − R0 ),

c2 =

from which the result below follows. Theorem 3.6. Consider the treatment-free model (3.1) with κ = 1. The resulting model has a unique positive equilibrium if R0 > 1, and no positive endemic equilibrium otherwise. Proof. Consider the model (2.2) with κ = 1. It is clear from (3.9) that c2 > 0 and c0 ≥ 0 if R0 ≤ 1. Furthermore, using the fact that

(μ + σ + ε)(μ + δ) . πε + μ + δ

(3.11)

Thus, it follows from (3.10) and (3.11) that c1 > γ (δη + μ)(ασ + μ + κε)(δ + μ + ε) + μψ(δ + μ)(δη + ε + μ)(ε + μ + σ )(1−R0 ) + (μ + δ)(μ + σ + ε) 

(μ + σ + ε)(μ + δ)[(μ + ηδ) + π κε] (μ + ασ + κε)(μ + ηδ)− , (3.12) π ε + (μ + δ) = γ (δη + μ)(ασ + μ + κε)(δ + μ + ε) + μψ(δ + μ)(δη + ε + μ)(ε + μ + σ )(1−R0 ) + (μ + δ)(μ + σ + ε) 

(μ+ασ +κε)(μ+ηδ)[π ε +(μ+δ)]−(μ+σ +ε)(μ+δ)[(μ+ηδ)+π κε ] , π ε +(μ+δ) = γ (δη + μ)(ασ + μ + κε)(δ + μ + ε) + μψ(δ + μ)(δη + ε + μ)(ε + μ + σ )(1−R0 ) + (μ + δ)(μ + σ + ε)M,

where,

M=

M1 + M2

πε + μ + δ

,

with,

M1 = π ε [(μ + ασ + κε)(μ + ηδ) − (μ + σ + ε)(μ + δ)], M2 = (μ + δ)(μ + ηδ)(ασ + κε − σ − ε). The sign of coefficient c1 can be deduced from (3.12) as follows. Since η > 1 and α > 1, it follows that (μ + ασ + ε) > (μ + σ + ε), (μ + ηδ) > (μ + δ) and (ασ + ε) > (σ + ε). Consequently, if R0 < 1 and κ = 1, then c1 > 0 (and the quadratic (3.7) will have no endemic equilibrium for this special case). The above analysis (along with Theorem 3.4) reveals that the relative rate of progression of acute re-infected individuals to the chronic re-infection stage (κ ) plays a critical role in the existence of the phenomenon of backward bifurcation in the treatment-free model (3.1). In fact, setting κ = 1 in (3.8), and using the assumptions η ≥ 1, α ≥ 1 and υ ≤ 1, ω ≤ 1 and ψ ≤ 1, it can be shown that the backward bifurcation coefficient (a), given by (A.3) in Appendix A, becomes negative (and, consequently, in line with Theorem 4.1 of [6], the treatmentfree model (3.1) does not undergo backward bifurcation in this case). To do so, we just need to show that the term S2 , in the expression for the backward bifurcation coefficient a (given by (A.3) in Appendix A), is positive (in which case, a < 0). This is shown below.

G1 G2 ψ(π εκω + G4 υ) (π ε + G2 ) G3 G4 (π ε + G2 ) − G1 G2 ψ(π εκω + G4 υ) = , (G2 + π ε) G3 G4 (π ε + G2 ) − G1 G2 (π εκ + G4 ) > , (π ε + G2 ) π ε(G3 G4 − G1 G2 ) + G2 G4 (G3 − G1 ) = > 0, (π ε + G2 )

S2 = G 3 G 4 −

(3.13)

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

since, assuming κ = 1, G3 = μ + ασ + ε > G1 = μ + σ + ε and G4 = μ + ηδ > G2 = μ + δ (so that, G3 G4 − G1 G2 > 0 and G3 − G1 > 0). Thus, it follows, based on Theorem 3.6 and Item (iv) of Theorem 4.1 in [6], that the treatment-free model (3.1) does not undergo backward bifurcation in this case (with κ = 1). Hence, this study shows, for the first time, that the relative rate of progression from acute reinfection to chronic re-infection stage (κ ) induces the phenomenon of backward bifurcation in HCV transmission dynamics. Case (ii). Effect of infectiousness of acute and chronic re-infected individuals (υ , ω) Consider, now, the special case of the treatment-free model (3.1) where acute and chronic re-infected individuals do not transmit HCV infection (i.e., υ = ω = 0). For this special case, it follows from (3.13) that S2 = G3 G4 > 0. Consequently, the backward bifurcation parameter, a (given by (A.3) in Appendix A), is negative (so that backward bifurcation does not occur in this case, in line with Theorem 4.1 of [6]). It is worth emphasizing that, for this scenario (υ = ω = 0), backward bifurcation does not occur in the treatment-free model even in the presence of re-infection (ψ = 0). Hence, this study shows, for the first time, that disease transmission by acute and chronic re-infection individuals also induces the phenomenon of backward bifurcation in HCV transmission dynamics. Case (iii). Effect of average infectivity-adjusted duration of infection/re-infection The effect of average duration of infectivity for primary and reinfected individuals, on the backward bifurcation property of the model (3.1), will now be analyzed. The aim is to explore whether or not the model (3.1) has positive equilibria when R0 ≤ 1 and the average infectivity-adjusted duration of re-infection is less than or equal to that of primary infection. It should be recalled that, for the model (3.1), each type of infection consists of two stages: acute and chronic. Furthermore, the average duration of primary acute and chronic infection is, respectively, given by

1 ε 1 , . μ+σ +ε μ+σ +ε μ+δ Adjusting for infectivity (1 for acute and π for chronic infection), and adding up the two durations, it follows that the average infectivityadjusted duration of primary infection (denoted by DP) is given by:

DP =

(μ + δ + π ε) . (μ + σ + ε)(μ + δ)

(3.14)

Similarly, the average duration of acute and chronic re-infection is given, respectively, by

κε , , μ + ασ + κε μ + ασ + κε μ + ηδ 1

1

and adjusting for infectivity (υ for acute and π ω for chronic infection) gives the average infectivity-adjusted duration of re-infection (denoted by DR):

DR =

υ(μ + ηδ) + κεπ ω . (μ + ασ + κε)(μ + ηδ)

(3.15)

Theorem 3.7. Consider the model (3.1) with DR ≤ DP. If R0 ≤ 1, then the model has no positive endemic equilibrium. Proof. Consider the model (3.1) with DR ≤ DP and R0 ≤ 1. It is sufficient to show that the coefficient, c1 , of the quadratic (3.7), is positive (in which case, the DFE, E0wt , is the only equilibrium of the model (3.1)). Solving for β from setting R0 = 1 gives

β=

(μ + σ + ε)(μ + δ)R0 . μ + δ + πε

Using (3.16) in the expression for c1 in (3.8) gives

(3.16)

57

(δη + μ)(ασ + μ + κε)γ (δ + μ + ε) + μ(1 − R0 )ψ(δ + μ)(μ + σ + ε)(δη + μ + κε) + (δη + μ)(ασ + μ + κε)(μ + δ)(μ + σ + ε) (3.17) R0 ψ(δ + μ)(μ + σ + ε)[δ(σ + ε) + μσ ](δηυ + μυ + ωεκπ ) − . δ + μ + πε

c1 =

Since the first two terms in (3.17) are positive (for R0 ≤ 1), we focus on the last two terms and rewrite them as

(δη + μ)(ασ + μ + κε)(μ + δ)(μ + σ + ε) × 

R0 ψ [δ(σ + ε) + μσ ][(μ + δη)υ + ωεκπ ] 1− . (δ + μ + π ε)(δη + μ)(ασ + μ + κε)

(3.18)

Since R0 and ψ are less than or equal to unity and δ(σ + ε) + μσ < (μ + δ)(μ + σ + ε), the second term in (3.18) can be further simplified to (noting that DR ≤ DP) R0 ψ [δ(σ + ε) + μσ ][(μ + δη)υ + ωεκπ ]

(δ + μ + π ε)(δη + μ)(ασ + μ + κε) R0 ψ [δ(σ + ε) + μσ ]DR < 1. = (δ + μ)(μ + σ + ε)DP

Therefore, it follows from (3.19) that (3.18) is positive. Hence, c1 > 0, as required. Theorem 3.7 suggests that the model (3.1) does not undergo backward bifurcation at R0 = 1 when DR ≤ DP. This can be proven by showing that, for R0 ≤ 1 and DR ≤ DP, the backward bifurcation coefficient (a), given by (A.3) in Appendix A, becomes negative (and, consequently, in line with Theorem 4.1 of [6], the treatment-free model (3.1) does not undergo backward bifurcation in this case). To do so, it is sufficient to show that the term S2 , in the expression for the backward bifurcation coefficient a (given by (A.3) in Appendix A), is positive (in which case, a < 0). This is shown below. First of all, using the definitions of G1 , . . . , G4 , given in Section 3.2.1, the expressions for DP and DR can be rewritten, respectively, as:

DP =

G2 + π ε G1 G2

and DR =

π εκω + υ G4 G3 G4

.

(3.19)

Hence,

G1 G2 ψ(π εκω + G4 υ) S2 = G 3 G 4 − (π ε + G2 )   G1 G2 ψ(π εκω + G4 υ) = G3 G4 1 − (π ε + G2 )G3 G4  ψ DR . = G3 G4 1 − DP

(3.20)

Thus, since ψ < 1, it follows that S2 > 0 whenever DR < DP. That is, the model (3.1) does not undergo backward bifurcation whenever the average infectivity-adjusted duration of re-infected individuals is less than that of primary infected individuals. In summary, the treatment-free model (3.1) has the following dynamical features: (i) The DFE of the model is LAS when R0 < 1, and unstable if R0 > 1; (ii) The model undergoes backward bifurcation at R0 = 1 under certain conditions. This phenomenon does not exist if: (a) there is no re-infection (i.e., ψ = 0); (b) acutely-re-infected individuals progress to chronic-reinfection stage at the same rate as acute-primary infected individuals progress to chronic-primary infection stage (i.e., κ = 1); (c) acutely and chronically re-infected individuals do not transmit infection (in the presence of the aforementioned differential characteristics of infection);

58

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

(d) the average infectivity-adjusted duration of re-infected individuals is less than that of primary infected individuals (DR < DP). (iii) The model has a unique endemic equilibrium point whenever R0 > 1.

Solving the equations of the treatment model (2.2) at the endemic steady-state gives:

S∗∗ = R∗∗ =

4. Analysis of the treatment model

V ∗∗ =

Consider, now, the treatment model (2.2).

γ R∗∗ +  λ∗∗ S∗∗ ε I∗∗ + ρ T ∗∗ , I∗∗ = , P ∗∗ = , ∗∗ λ +μ K1 K2 σ I∗∗ + δ P∗∗ + ασ V ∗∗ + ηδ W ∗∗ + θ T ∗∗ + θ Q ∗∗ , (μ + γ + ψλ∗∗ ) ψλ∗∗ R∗∗ κε V ∗∗ + ζ Q ∗∗ ∗∗ τP

K3 ∗∗

,

W

φW

K4 ∗∗

,

4.1. Existence and asymptotic stability of equilibria

T ∗∗ =

4.1.1. Local asymptotic stability of DFE The DFE of the treatment model (2.2) is given by

Substituting the expressions in (4.5) into (4.4) gives:



E0T = S∗ , I∗ , P ∗ , R∗ , V ∗ , W ∗ , T ∗ , Q

 ∗





 = , 0, 0, 0, 0, 0, 0, 0 . (4.1) μ

Using the next generation operator method [44], as in Section 3, the matrices FT and HT , associated with the treatment model (2.2), are given, respectively, by

⎛

β

⎜0 ⎜ ⎜ ⎜0 ⎜ FT = ⎜ ⎜0 ⎜ ⎜0 ⎝ 0 ⎛ K1 ⎜ ⎜−ε ⎜ ⎜0 ⎜ HT = ⎜ ⎜0 ⎜ ⎜ ⎝0 0

π ωβ

υβ

π χT β

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0 ⎞

0

0

0

0

K2

0

0

−ρ

0

K3

0

0

0

−κε

K4

0

−τ

0

0

K5

0

0

−φ

0

0

⎟ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎠

RT = ρ FT HT−1



γ K3 (K2 K5 − ρτ )(K4 K6 − ζ φ)(μ + γ ) γ (K2 K5 − ρτ )[K3 ψ(K4 K6 − ζ φ) − ψ(ασ (K4 K6 − ζ φ) + K6 ηδκε + θ φκε)], = μγ K3 (K4 K6 − ζ φ)[K1 (K2 K5 − ρτ ) − β(K2 K5 − ρτ + π ε K5 + χ π ετ ), = γ K1 K3 (K2 K5 − ρτ )(K4 K6 − ζ φ), (4.7) = γ K3 (K2 K5 − ρτ )(K4 K6 − ζ φ) + βμψ(K2 K5 − ρτ )[υ(K4 K6 − ζ φ) + κεωπ K6 + χQ π φκε ], = γ K3 (K4 K6 − ζ φ)[σ (K2 K5 − ρτ ) + δε K5 + θ ετ ].

B2 =

B6

a2 (λ∗∗ )2 + a1 λ∗∗ + a0 = 0,

(4.8)

where,

γ 2 K1 K3 A1 (K4 K6 − ζ φ)(K2 K5 − ρτ )2 , a1 = μγ 2 K1 K3 A1 (K4 K6 − ζ φ)(K2 K5 − ρτ )2 (1 − RT ) + γ 2 K32 K1 (μ + γ )(K4 K6 − ζ φ)2 (K2 K5 − ρτ )2 (4.9) − A2 [γ K3 (K4 K6 − ζ φ)(σ (K2 K5 − ρτ ) + δε K5 + θ ετ )], a0 = γ 2 K1 K32 (μ + γ )(1 − RT )(K4 K6 − ζ φ)2 (K2 K5 − ρτ )2 ,

a2 =



(4.2)

It should be mentioned that the expression for RT is the same ¯ c in [13]. The result below follows from Theorem 2 as that of R of [44]. Theorem 4.1. The DFE, E0T , of the treatment model (2.2), given by (4.1), is LAS if RT < 1, and unstable if RT > 1. 4.1.2. Existence of EEP As in Section 3, let

(4.3)

be an arbitrary endemic equilibrium of the treatment model (2.2). Furthermore, let

λ∗∗ =

(4.6)

It follows that the non-zero (endemic) equilibria of the treatment model (2.2) satisfy the following polynomial (in terms of λ∗∗ ),

β K2 K5 − ρτ + π ε K5 + χT π ετ = . K1 (K2 K5 − ρτ )

E1T = (S∗∗ , I∗∗ , P ∗∗ , R∗∗ , V ∗∗ , W ∗∗ , T ∗∗ , Q ∗∗ ),

.

B1 =

B4

⎟ 0 ⎟ ⎟ 0 ⎟ ⎟ ⎟, −ζ ⎟ ⎟ ⎟ 0 ⎠ K6



K6

where,

B5

where, K1 = μ + σ + ε , K2 = μ + δ + τ , K3 = μ + ασ + κε , K4 = μ + ηδ + φ , K5 = μ + ρ + θ , and K6 = μ + ζ + θ . It follows, from Theorem 2 in [44], that the control reproduction number of the model (2.2), defined by RT = ρ(FT HT−1 ), is given by (where K2 K5 − ρτ = (μ + δ + τ )(μ + θ ) + μρ + δρ > 0, so that RT > 0)



,

(B1 + B2 λ∗∗ )(B3 + B4 λ∗∗ ) = B5 B6 λ∗∗ ,

B3

π χQ β ⎞

πβ

K5

Q ∗∗ =

=

(4.5)

βμ(I∗∗ + π P∗∗ + υ V ∗∗ + π ωW ∗∗ + π χT T ∗∗ + π χQ Q ∗∗ ) ,  (4.4)

(where the total population, N(t), is now replaced by its limiting value, N∗ =  μ ) be the force of infection at steady-state.

with,

ψ(K4 K6 − ζ φ)(μ + κε) + K6 ψηδκε + ψθ φκε , (4.10) A2 = (K2 K5 − ρτ ){γ K3 (K4 K6 − ζ φ) + βμ[υψ(K4 K6 − ζ φ) + κεψωπ K6 + χ π φκεψ ]}.

A1 =

The endemic equilibria of the treatment model (2.2) can then be obtained by solving for λ∗∗ from (4.4), and substituting the positive values of λ∗∗ into the steady-state expressions in (4.5). Furthermore, it follows from (4.9) that the coefficient a2 , of the quadratic (4.8), is always positive (it should be recalled, from Section 4.1.1, that K4 K6 − ζ φ = (μ + ηδ + φ)(μ + θ ) + ζ (μ + ηδ) > 0 and K2 K5 − ρτ > 0) and a0 is positive (negative) if RT is less (greater) than unity. The quadratic has a unique endemic equilibrium whenever RT > 1. These results are summarized below. Theorem 4.2. The treatment model (2.2) has: (i) a unique endemic equilibrium if a0 < 0 ⇔ RT > 1; (ii) a unique endemic equilibrium if a1 < 0 and a0 = 0 or a21 − 4a0 a2 = 0; (iii) two endemic equilibria if a0 > 0, a1 < 0 and a21 − 4a0 a2 > 0; (iv) no endemic equilibrium otherwise. Here, too, Item (iii) of Theorem 4.2 suggests the possibility of the backward bifurcation in the model (2.2). The backward bifurcation phenomenon in the treatment model (2.2) is explored below.

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

Fig. 3. Backward bifurcation diagram for the treatment model (2.2), showing the prevalence as a function of the control reproduction number (RT ). Parameter values used are:  = 40 640, β ∗ = 3.324462329, μ = 0.091, σ = 0.255, δ = 0.0025, ε = 3.5, ψ = 0.70, α = 3.3, η = 4.050, κ = 0.0827, τ = 0.70, φ = 0.36, θ = 0.7603, ζ = 0.902, ρ = 0.638, γ = 0.0225, π = 0.019, υ = 0.89999, ω = 0.89999χT = 0.35, and χQ = 0.745 (so that, a = 0.0001180441877 > 0; RT = 1).

Theorem 4.3. The treatment model (2.2) undergoes backward bifurcation at RT = 1 whenever the Inequality (B.5), given in Appendix B, holds. The proof of Theorem 4.3 is given in Appendix B (and a schematic description of the backward bifurcation phenomenon of the treatment model (2.2) is depicted in Fig. 3). It should be recalled that this (backward bifurcation) phenomenon does not occur in the reduced model (2.5), considered in [13]. In other words, the treatment model (2.2) has at least one dynamical feature (backward bifurcation) that is not present in the reduced model (2.5). Thus, like in the case of the treatment-free model (3.1), it is instructive to explore the possible causes of the backward bifurcation phenomenon in the treatment model (2.2). This is done below. 4.1.3. Non-existence of backward bifurcation Two main cases will be considered, as follows. Case 1. Absence of re-infection of recovered individuals (ψ = 0) Consider the treatment model (2.2) in the absence of re-infection (i.e., ψ = 0). Setting the re-infection parameter (ψ ) to zero in the expression of the backward bifurcation coefficient, a, given by (B.3) in Appendix B (it should be recalled that K4 K6 − ζ φ > 0 and K2 K5 − ρτ > 0) shows that

a=

−2μv2 w22 K1 [K (K K − ζ φ)(μ + γ ) K3 (K4 K6 − ζ φ)(K2 K5 − ρτ )(μ + γ ) 3 4 6

×(K2 K5 − ρτ + ε K5 + τ ε) + K3 (K4 K6 − ζ φ)F] < 0.

(4.11)

Thus, it follows from (4.11) and Appendix B (where the eigenvector v2 > 0 and F > 0) that, in the absence of re-infection (ψ = 0), the bifurcation coefficient, a (given by (B.3)), is negative. Hence, based on the Item (iv) of Theorem 4.1 in [6], it can be concluded that the treatment model (2.2) does not undergo backward bifurcation in the absence of re-infection. Hence, the re-infection of recovered individuals causes backward bifurcation in the treatment model (2.2). It is worth mentioning that the reason that the treatment model (2.2) exhibits backward bifurcation (in the presence of re-infection), while the model (2.5) with re-infection (considered in [13]) does not, is (clearly) the differential characteristics of primary infection and reinfection. In the model (2.5), no such differential characteristics exist (and, in such a case, unlike in the model (2.2) where such differential exists, re-infection does not cause backward bifurcation in the model

59

Fig. 4. Simulations of the treatment model (2.2), showing the total number of infected individuals as a function of time, using various initial conditions. Parameter values used are as given in Table 2, with β = 1.68 (so that, RT = 0.8323).

(2.5)). To further confirm the absence of backward bifurcation in the treatment model (2.2) without re-infection, the following result is proved for the DFE of the model (2.2) for this special case. Theorem 4.4. The DFE (E0T ) of the treatment model (2.2), with ψ = 0, is GAS in DT whenever RT < 1. The proof, based on using Comparison Theorem [28], is given in Appendix C. Fig. 4 depicts the solutions profile of the model (2.2) generated for the case when ψ = 0 and RT < 1, using various initial conditions, showing the convergence to the DFE (in line with Theorem 4.4). Case 2. Presence of re-infection (ψ = 0). Consider the treatment model (2.2) in the presence of re-infection of recovered individuals (i.e., ψ = 0). The following cases are considered. (i) No differential characteristics between primary infection and re-infection. For this case (with κ = α = η = ω = υ = 1, χQ = χT , ζ = ρ and φ = τ ), the treatment model (2.2) reduces to the model (2.5). Elbasha [13] proved the GAS property of the DFE of the model for this special case (ruling out the possibility of backward bifurcation in the model (2.5)). The result in [13] is further verified by applying the centre manifold theory on the model (2.5), as detailed in Appendix D (from which it is clear that the reduced model (2.5) does not undergo backward bifurcation). Thus, the analysis in Appendix D shows that the absence of the differential characteristics of infection between re-infected and primary-infected individuals removes the backward bifurcation property of the model (2.2), even in the presence of the re-infection of recovered individuals. It is worth stating that the aforementioned conditions can be relaxed to κ = η = 1, χQ = χT , ζ = ρ and φ = τ (i.e., for υ ≤ 1, ω ≤ 1 and α ≥ 1 ), as proved at the end of Appendix B. (ii) Effect of infectiousness of acute, chronic and treated reinfected individuals (υ = 0, ω = 0, χQ = 0). Consider the case of the treatment model (2.2) where acute, chronic and treated re-infected individuals do not transmit infection (i.e., υ = ω = χQ = 0). For this case,

ˆ = λ|υ =ω=χ =0 = λ Q

β(I + π P + χT π T ) N

,

60

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

active IDU population. In particular, upon primary infection, a susceptible individual enters the acute infection stage (I) and, assuming survival, proceeds to the chronic compartment P. This individual can move back-and-forth between P and T classes, until he/she leaves the active population. Of the newly-infected individuals entering into compartment I, a fraction μ+εσ +ε , survives and proceeds to stage P. For an individual in stage P, the probability of proceeding to stage (T), and then returning to stage P, is given by

GPT =

ρ τ . μ+τ +δ μ+ρ +θ

Thus, given that an individual has reached P initially, the expected number of visits to class P, denoted by vP , is given by:

vp = 1 + GPT + G2PT + · · · , 1 , 1 − GPT (μ + τ + δ)(μ + ρ + θ ) = . (μ + τ + δ)(μ + ρ + θ ) − ρτ =

Fig. 5. Simulations of the treatment model (2.2), showing the total number of infected individuals as a function of time, using various initial conditions. Parameter values used are as given in Table 2 (so that, RT = 1.4574).

where, N = S + I + P + R + V + W + T + Q. Setting υ = ω = χQ = 0 in the expression of the backward bifurcation coefficient, a (given by (B.3) in Appendix B), and simplifying, shows that

a=

−2μv2 w22 K1 [K (K K − ζ φ)(μ + γ ) K3 (K4 K6 ζ φ)(K2 K5 − ρτ )(μ + γ ) 3 4 6 × (K2 K5 − ρτ + ε K5 + τ ε) + K3 (K4 K6 − ζ φ)F] < 0. (4.12)

Additionally, since a fraction, μ+ττ +δ , of individuals in class P proceed to class T, the expected number of visits to class T, denoted by vT , given that an individual has survived and reached class P, is given by:

vT =

τ τ (μ + ρ + θ ) v = . μ + τ + δ P (μ + τ + δ)(μ + ρ + θ ) − ρτ

(4.13)

Thus, the expected total time that a newly-infected individual spends in compartments I, P and T is given, respectively, by

1

TItotal = TI =

,

(μ+σ +ε) ε ε (μ+ρ +θ ) = T v = . , (μ+σ +ε) 2 P (μ+σ +ε) (μ+τ +δ)(μ+ρ +θ ) − ρτ ε ε τ = T v = . . (μ+σ +ε) 3 T (μ+σ +ε) (μ+τ +δ)(μ+ρ +θ ) − ρτ

Thus, it follows, from (4.12) and Appendix B, that disease transmission by re-infected individuals causes backward bifurcation in HCV transmission dynamics. Therefore, it was shown that if individuals in the re-infected classes (V, W, Q ) are not able to transmit HCV infection, then the treatment model (2.2) does not undergo backward bifurcation even in the presence of re-infection (ψ = 0). This fact is further illustrated by proving the global asymptotic stability of the DFE of the treatment model (2.2) for this special case, as below.

Tptotal

Theorem 4.5. The DFE, E0T , of the treatment model (2.2), with υ = ω =

DPtreatment (μ + τ + δ)(μ + ρ + θ ) − ρτ + π ε(μ + ρ + θ ) + χT π τ ε . = (μ + σ + ε)[(μ + τ + δ)(μ + ρ + θ ) − ρτ ]

χQ = 0, is GAS in DT whenever RT < 1.

The proof, based on using Comparison Theorem [28], is given in Appendix E. As in Section 3, the result given in Theorem 4.5 shows that HCV can be effectively-controlled (or eliminated) in the IDU population if the infectivity of re-infected individuals is negligible (or, for instance, cured IDUs change their behaviour, and cease being IDUs). Moreover, numerical simulations of the treatment model (2.2), for the case when RT > 1 (Fig. 5), suggest that the associated unique endemic equilibrium (E1T ) is stable when it exists. (iii) Effect of infectivity-adjusted duration of infection/reinfection Next we verify how the infectivity-adjusted duration of infection/re-infection can affect the qualitative dynamics of the treatment model (2.2). It should be recalled that each type of infection consists of three stages. For example, primary infection consists of acute (I), chronic (P) and treated (T) stages. The average duration that a primary infected individual spends in these stages is given, respectively, by TI = μ+1σ +ε , TP = μ+1τ +δ and TT = μ+1ρ +θ . This individual

passes through the infective stages until he/she leaves the active IDU population. Associated with an infected individual is a Markov chain that describes the individual’s passage through the infective stages (see also [22]). This process continues until the individual leaves the

TTtotal

Adjusting for infectivity (1 for acute, π for chronic and π χT for treated individuals), and adding up the three durations (TItotal , Tptotal , TTtotal ), it follows that the average infectivity-adjusted duration of primary infection (for the treatment model (2.2)) is given by:

Similarly, the average duration in the acute (V), chronic (W) and treated (Q) re-infected classes is, respectively, given by:

TVtotal = total TW

TQtotal

1

,

(μ + ασ + κε) (μ + ζ + θ ) κε = , (μ + ασ + κε) (μ + φ + ηδ)(μ + ζ + θ ) − φζ φ κε = . (μ + ασ + κε) (μ + φ + ηδ)(μ + ζ + θ ) − φζ

Adjusting for infectivity (υ for acute, π ω for chronic and π χQ for total , TQtotal ), the treated), and adding up the three durations (TVtotal , TW average infectivity-adjusted duration in the re-infected classes (for the treatment model (2.2)) is given by: DRtreatment υ [(μ + φ + ηδ)(μ + ζ + θ) − φζ ] + π ωκε(μ + ζ + θ) + χQ π φκε = . (μ + ασ + κε)[(μ + φ + ηδ)(μ + ζ + θ) − φζ ]

We claim the following result. Theorem 4.6. The treatment model (2.2) with DRtreatment ≤ DPtreatment does not undergo backward bifurcation at RT = 1.

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

61

Proof. Consider the treatment model with RT = 1 and DRtreatment ≤ DPtreatment . Further, consider the terms C1 + C2 F − β ∗ ψ FC3 in the expression for the backward bifurcation coefficient, a, given in Eq. (B.3) of Appendix B. Using the definitions of C2 and C3 , given in Eq. (B.4), it follows that

C2 F − β ∗ ψ FC3 = K3 (K4 K6 −ζ φ)F−β ∗ ψ F[(K4 K6 −ζ φ)υ +K6 εκπ ω+φεκπ χQ ],   β ∗ ψ [(K4 K6 −ζ φ)υ +K6 εκπ ω+φεκπ χQ ] = K3 (K4 K6 −ζ φ)F 1− , K3 (K4 K6 −ζ φ)   β ∗ [(K4 K6 −ζ φ)υ +K6 εκπ ω+φεκπ χQ ] , ≥ K3 (K4 K6 −ζ φ)F 1− K3 (K4 K6 −ζ φ) since, by assumption, ψ ≤ 1. On the other hand, using definitions of K1 , . . . and K6 , given in Section 4.1.1, the expressions for DPtreatment and DRtreatment can be rewritten as:

(μ+τ +δ)(μ+ρ +θ ) − ρτ +π ε(μ+ρ +θ )+χT π τ ε , (μ+σ +ε)[(μ+τ +δ)(μ+ρ +θ ) − ρτ ] 1 K2 K5 − ρτ + π ε K5 + χT π τ ε = ∗, (4.14) = K1 (K2 K5 − ρτ ) β

DPtreatment =

where β ∗ is given by (B.2) in Appendix B, and DRtreatment υ [(μ + φ + ηδ)(μ + ζ + θ) − φζ ] + π ωκε(μ + ζ + θ) + χQ π φκε = , (μ + ασ + κε)[(μ + φ + ηδ)(μ + ζ + θ) − φζ ] =

υ(K4 K6 − ζ φ) + π ωκεK6 + χQ π κεφ . K3 (K4 K6 − ζ φ)

Substituting (4.14) and (4.15) into the last line of C2 F − β ∗ ψ FC3 gives: C2 F − β ∗ ψ FC3

  β ∗ [(K4 K6 −ζ φ)υ +K6 εκπ ω+φεκπ χQ ] , ≥ K3 (K4 K6 −ζ φ)F 1− K3 (K4 K6 −ζ φ)   DRtreatment . = K3 (K4 K6 − ζ φ)F 1 − DPtreatment

(4.15)

Therefore, since C1 > 0 (by definition given in (B.4)), it follows that C1 + C2 F − β ∗ ψ FC3 ≥ 0 if and only if DRtreatment ≤ DPtreatment . Hence, the backward bifurcation coefficient a, given by (B.3) in Appendix B, is negative whenever DRtreatment ≤ DPtreatment (or, equivalently, the treatment model (2.2) does not undergo backward bifurcation at RT = 1 if the average infectivity-adjusted duration of re-infection is less than that of primary infection).

The qualitative results for the treatment model (2.2) are summarized below (a) The model does not undergo backward bifurcation in the absence of re-infection of recovered individuals (ψ = 0). (b) In the presence of re-infection of recovered individuals (ψ = 0), the phenomenon of backward bifurcation does not exist under any of the following scenarios: (i) the absence of differential characteristics of infection between re-infected and primary-infected individuals (i.e., κ = α = η = ω = υ = 1, χQ = χT , ζ = ρ , and φ = τ ); (ii) if re-infected individuals do not transmit HCV infection (i.e., υ = ω = χQ = 0); (iii) if the average infectivity-adjusted duration of re-infected individuals is less than that of primary infected individuals (DRtreatment < DPtreatment ). (c) The disease-free equilibrium of the model is globallyasymptotically stable when RT < 1 for the following scenarios: (i) no differential characteristics between primary infection and re-infection [13]; (ii) acute, chronic and treated individuals do not transmit HCV infection;

Fig. 6. Backward bifurcation diagram for the treatment model (2.2), showing the prevalence as a function of the control reproduction number (RT ) when κ = 1. Parameter values used are:  = 40 640, μ = 0.00091, σ = 0.355, δ = 0.00025, ε = 3.5, ψ = 0.90, α = 1.3 , η = 1, τ = 0.836, φ = 0.26, θ = 0.7603, ζ := 0.756, ρ = 0.338, γ 1 = 0.0, π = 0.019, υ = 0.8999, ω = 0.8999, χT = 0.25, and χQ = 0.945 (so that, RT = 1).

(iii) no re-infection of recovered individuals. The analyses in Sections 3 and 4 show that both the treatmentfree and the treatment model, (3.1) and (2.2), respectively, exhibit the same qualitative dynamics with respect to the asymptotic stability of the respective disease-free equilibrium of the model as well as the backward bifurcation property of each of the models. However, unlike in the treatment-free model (where setting the modification parameter κ = 1 removes its backward bifurcation property), the treatment model (2.2) imposes additional conditions for the backward bifurcation not to exist when κ (the relative progression rate from acute re-infection to chronic re-infection stage) is set to unity. In this case, we , additionally, need to set α = η = ω = υ = 1, χQ = χT , ζ = ρ , and φ = τ ; which is equivalent to removing the differential in the characteristics of between primary infection and re-infection. Furthermore, Fig. 6 shows that, unlike in the treatment-free model (3.1), setting κ = 1 is not enough to remove the backward bifurcation property of the treatment model (2.2). 4.2. Assessment of treatment impact Following Elbasha [13], the reproduction threshold (RT ) is differentiated partially with respect to the treatment rate of chronicallyinfected individuals (τ ), giving

π ε(μ + ρ + θ )

∂ RT β =− , (4.16) ∂τ (μ + σ + ε) [(μ + δ)(μ + θ + ρ) + τ (μ + θ )]2 where,

= (μ + θ ) − χT (μ + δ). Thus, RT is a decreasing (increasing) function of τ whenever > 0 ∂R (<). Furthermore, ∂τ T = 0 if = 0. This leads to the following result (same result was also derived for the reduced model considered in [13]). Theorem 4.7. Consider the treatment model (2.2) in the absence of backward bifurcation. The treatment of chronically-infected individuals offers: (i) a positive population-level impact whenever > 0 ; (ii) no-population level impact if = 0;

62

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69 Table 2 Baseline values and ranges of the parameters of the treatment model (2.2). Parameters

Baseline values [13]

Ranges [13]

 μ β σ δ ε ψ α η κ τ φ θ ρ ζ γ π υ ω χT = χQ = χ

39, 600 year−1 0.09 year−1 2.68 year−1 0.5 year−1 0.002 year−1 1.5 year−1 0.5 3.3 3.3 1/3.3 0.04 year−1 0.04 year−1 0.67 year−1 0.82 year−1 0.82 year−1 0.025 year−1 0.01 1/6.5 1/6.5 0.5

[35640, 43560] [0.081, 0.099] [2.444, 2.948] [0.45, 0.55] [0.0018, 0.0022] [1.35, 1.65] [0.45, 0.55] [2.97, 3.36] [2.97, 3.36] [0.2727, 0.3333] [0.036, 0.044] [0.036, 0.044] [0.603, 0.737] [0.738, 0.902] [0.738, 0.902] [0.0225, 0.0275] [0.009, 0.011] [0.1386, 0.1694] [0.1386, 0.1694] [0.45, 0.55]

Fig. 7. Boxplots of the control reproduction number (RT ), as a function of the number of LHS runs (NR ) carried out, for the treatment model (2.2). Parameter values and ranges used are as given in Table 2.

(iii) a detrimental impact population-level (increase disease burden) if < 0. As noted by Elbasha [13], the threshold quantity, , is expected to always be positive since the cure rate (θ ) is expected to exceed the natural recovery for chronically-infected individuals (δ ), and that the relative infectiousness of treated individuals is small (χT < 1). Similarly, differentiating the reproduction threshold (RT ) partially with respect to the treatment failure rate for chronically-infected individuals (ρ ) gives:

π ετ

∂ RT β = , ∂ρ (μ + σ + ε) [(μ + δ)(μ + θ + ρ) + τ (μ + θ )]2

(4.17)

so that RT is an increasing function of ρ whenever > 0, as expected. 4.3. Uncertainty and sensitivity analysis The treatment model (2.2) contains 21 parameters, and the effect of the uncertainties in the estimates of the parameter values used in the numerical simulations of the model (2.2) [3,36] are assessed using Latin Hypercube Sampling (LHS). The LHS method involves defining baseline values and ranges for each of the parameters of the model (as in Table 2), where each parameter is assumed to obey a uniform distribution [13], and carrying out multiple runs (NR = 1000) of the sampled data for the response output (the control reproduction threshold, RT , in this case) [3,36]. Furthermore, the sensitivity of the parameters of the treatment model (2.2) is measured by finding PRCC between each parameter and control reproduction number (RT ). A boxplot of the control reproduction number (RT ), as a function of the number of LHS runs carried out, is depicted in Fig. 7, showing a range of RT from 1.44 to 1.48 (which is consistent with the range in [13]). The boxplot corresponding to an increased treatment rate (70%) is also depicted in Fig. 8, showing a marked decrease of RT range (of RT ∈ [1.33, 1.37]), as expected. Furthermore, Table 3 and Fig. 9 give the PRCC values of the parameters of the model, from which it follows that the parameters that most affect RT (hence, drive the HCV transmission dynamics) are the effective contact rate (β ), the rate of progression to chronic infection (ε ), the rate of recovery from acute infection (σ ), natural death (μ) and the relative infectivity of chronically-infected individuals (π ). Further numerical simulations of the treatment model (2.2) (Figs. 10 and 11) show that a 10% increase in the baseline values of these top-PRCC ranked parameters (β , ε , σ , μ and π ) increases the cumulative incidence and prevalence of HCV in the community (and

Fig. 8. Boxplots of the control reproduction number (RT ), as a function of the number of LHS runs (NR ) carried out, for the treatment model (2.2) with increased treatment rate (τ = 70%). Parameter values and ranges used are as given in Table 2.

Table 3 PRCC values of the parameters of the treatment model (2.2), with RT as the output (response function). Parameter values and ranges used are as given in Table 2. Parameters

PRCC(RT )

β ε σ μ π θ τ δ ρ  η γ κ υ χ ζ α ψ ω φ

0.989521805 −0.972702393 −0.852219599 −0.703177587 0.637106865 −0.11532655 −0.097537395 0.054609189 0.053632664 0.035866564 0.033461692 −0.031457596 −0.017800737 −0.014450844 0.011501931 −0.010004813 −0.006900158 0.004427847 0.004084257 0.002943497

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

63

Fig. 9. Distribution of PRCC values for the parameters of the treatment model (2.2). Parameter values and ranges used are as given in Table 2. Fig. 11. Simulations of the treatment model (2.2), showing the prevalence of HCV as a function of time, for various values of the top-five PRCC-ranked parameters in Table 3 (β , ε , σ , μ and π ): green curve (10% decrease in the baseline values of the top-five PRCC-ranked parameters); blue curve (baseline values); red curve (10% increase in the baseline values of the top-five PRCC-ranked parameters). Parameter values used are as given in Table 2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 10. Simulations of the treatment model (2.2), showing the cumulative number of new infected individuals as a function of time, for various values of the top-five PRCCranked parameters in Table 3 (β , ε , σ , μ and π ): green curve (10% decrease in the baseline values of the top-five PRCC-ranked parameters); blue curve (baseline values); red curve (10% increase in the baseline values of the top-five PRCC-ranked parameters). Parameter values used are as given in Table 2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

that a 10% decrease lead to a decrease in the cumulative incidence and prevalence of the disease during the first few years (about 10 years), and a marginal increase shortly thereafter).

Fig. 12. Simulations of the treatment model (2.2), showing the cumulative number of new infected individuals as a function of time, for various values of the re-infection parameter (ψ ): green curve (ψ = 0.0), blue curve (ψ = 0.5) and red curve (ψ = 1.0). Parameter values used are as given in Table 2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

4.4. Numerical simulations The treatment model (2.2) is simulated using the parameters values given in Table 2 (unless otherwise stated). The following initial data (relevant to HCV dynamics in an IDU population [13]) is also used:

(S(0), I (0), P(0), R(0), V (0), W (0), T (0), Q (0)) = (439000, 550, 350, 50, 20, 10, 10, 10). The cumulative number of new cases, as a function of time, for various values of the re-infection rate (ψ ) is shown in Fig. 12. It is evident from Fig. 12 that re-infection has little or no effect on the cumulative incidence of HCV for the first few years. The effect is, however, noticeable after about 7 years. In other words, this study shows that re-infection has marginal effect on HCV burden (as measured in terms of cumulative number of new cases) in the community. Similar results are obtained for the prevalence of HCV (Fig. 13).

The effect of treatment of chronically-infected individuals on the cumulative incidence and prevalence of the disease is depicted in Figs. 14 and 15, respectively. It follows from these figures that, as expected, the treatment of chronically-infected individuals significantly reduces the cumulative incidence and prevalence of HCV in the community. Furthermore, Fig. 16 shows that the treatment of active (chronically-infected and -re-infected) IDUs, even if only a tiny percentage of the IDUs is treated, has a positive population-level impact (i.e., minimizes HCV burden in the IDU population). Furthermore, this figure shows that the treatment of a sizeable proportion of chronically-infected individuals (e.g., τ = 70%) significantly reduces the prevalence of HCV in the population. It is also evident from Fig. 16 that while a 4% treat rate leads to about 5% reduction in prevalence of HCV, the 70% treat rate leads to about 20% reduction in disease prevalence.

64

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

Fig. 13. Simulations of the treatment model (2.2), showing the prevalence of total infected individuals as a function of time, for various values of the re-infection parameter (ψ ): green curve (ψ = 0.0), blue curve (ψ = 0.5) and red curve (ψ = 1.0). Parameter values used are as given in Table 2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 15. Simulations of the treatment model (2.2), showing the prevalence of HCV as a function of time, for various values of the treatment rate of chronically-infected individuals (τ ): green curve (τ = 0.04), blue curve (τ = 0.4) and red curve (τ = 0.7). Parameter values used are as given in Table 2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 14. Simulations of the treatment model (2.2), showing the cumulative number of new infected individuals as a function of time, for various values of the treatment rate of chronically-infected individuals (τ ): green curve (τ = 0.04), blue curve (τ = 0.4) and red curve (τ = 0.7). Parameter values used are as given in Table 2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 16. Simulations of the model (2.2), showing the prevalence of HCV as a function of time, in presence and absence of anti-HCV treatment: blue curve shows the prevalence of HCV without treatment (τ = φ = ζ = ρ = θ = χT = χQ = 0), green curve exhibits the prevalence where only 4% of chronically-infected IDUs are treated (τ = 0.04) and red curve shows the prevalence of HCV for the case where 70% of chronically-infected IDUs are treated (τ = 0.7). Other parameter values used are as given in Table 2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

4.5. Discussions and conclusions The study considers the HCV treatment model, presented by Elbasha [13]. The aim was to extend the qualitative analyses in [13] (where only a special case of the treatment model was analyzed) and explore the possibility of new dynamical features (of the model) not observed or established in [13]. It was shown that the treatment model undergoes the phenomenon of backward bifurcation when the associated control reproduction number (RT ) is less than unity. One of the notable contributions of this study is that two main scenarios (based on whether or not a difference exists in characteristics of primary and re-infected individuals) where such phenomenon does not exist are identified. They are discussed below. The backward bifurcation property, of the HCV transmission model considered in this study, does not exist in the absence of differential characteristics of primary infected and re-infected recovered

individuals (with respect to infectivity, progression, treatment and recovery). In other words, the model will not undergo backward bifurcation if newly-infected and re-infected recovered individuals behave the same way (with respect to the rates of the infectivity, recovery, disease progression and treatment). This result supports, and extends, the results reported in [13], where a special case (5-dimensional) of the 8-dimensional model (2.2) was analyzed. Hence, this study shows that the differential characteristics of infection between primary and re-infected individuals can induce the phenomenon of backward bifurcation in the transmission dynamics of a disease (such as HCV). It is further shown that in the presence of such differential characteristics, the backward bifurcation phenomenon does not exist in three cases, namely, (i) when re-infection of recovered individuals

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

does not occur (ψ = 0); (ii) when acute, chronic and treated reinfected individuals do not transmit infection (i.e., υ = ω = χQ = 0) and (iii) when the average infectivity-adjusted duration of re-infected individuals is less than that of primary infected individuals (i.e., DRtreatment < DPtreatment ). For the first of these two cases, comparison theorem is used to prove the global asymptotic stability of the disease-free equilibrium of the model when the associated reproduction number is less than unity. It should, however, be stressed that the condition for the existence of backward bifurcation, established in this study, involves a multitude of parameters. The values of many of these parameters cannot, at the present time, be determined with any reasonable degree of confidence. In other words, in the absence of robust empirical evidence, we cannot determine the likelihood that the assumptions needed to produce backward bifurcation in HCV dynamics are satisfied. The treatment model is shown to have a unique endemic equilibrium when the associated control reproduction number (RT ) exceeds unity. In such a case, numerical simulations show that the disease will persist in the IDU population. Further simulations show that the re-infection of recovered individuals has marginal effect on the disease burden (as measured in terms of HCV cumulative incidence and prevalence in the community). Moreover, it is shown (as expected) that the treatment of chronically-infected individuals significantly reduces the cumulative incidence and prevalence of HCV in the IDU community. Important developments in screening, diagnosis, and treatment of chronic HCV patients have emerged or continue to emerge. The area of treatment of chronic HCV patients is witnessing unprecedented developments. For example, data presented at the 2014 Conference on Retroviruses and Opportunistic Infections (CROI) show new interferon-free drugs achieving cure rates of more than 90%, with simple and short durations of therapy. Such high efficacy, coupled with treatment scale up, brings the possibility of HCV elimination within reach. This makes the analysis of the phenomenon of backward bifurcation especially important, so as to find the minimum treatment level needed to eliminate HCV transmission in a given population. Unlike models with forward bifurcation, where it is sufficient to have the control reproduction number be below unity to achieve elimination, this model requires the control reproduction number to be below a threshold, that is further smaller than unity.

dx5 = f5 = dt dx6 = f6 = dt

65

ψλx4 − (μ + ασ + κε) x5 , 



κε x5 − μ + ηδ x6 ,

where,

β(x2 + π x3 + υ x5 + π ωx6 )

λwt =

6

i=1

xi

,

and f = [f1 , · · · , f6 ]T represents the vector field of the model (3.1). Evaluating the Jacobian of the system (A.1) at the DFE (E0wt ) gives:

⎛

−β

−βπ

γ

−βυ

0

βπ

0

βυ

−G2

0

0

0

0

β − G1 ε σ

δ

−μ − γ

ασ

ηδ

0

0

0

0

−G3

0

0

0

0

0

κε

−G4

⎜ ⎜ ⎜ ⎜ ⎜ wt J(E0 ) = ⎜ ⎜ ⎜ ⎜ ⎝

−μ 0

−βωπ

⎞

βωπ ⎟ ⎟

⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎟ ⎠

Consider the case of the model (A.1) with R0 = 1. Suppose, also, that β is chosen as the bifurcation parameter. Solving for β from R0 = 1 gives (where G1 and G2 are as defined in Subsection 3.2.1)

G1 G2

β∗ =

π ε + G2

.

(A.2)

The transformed system (A.1), with β = β ∗ , has a simple eigenvalue with zero real part (and all other eigenvalues have negative real parts). Hence, the centre manifold theory [6] can be used to analyzed the dynamics of (A.1) near β ∗ . To apply the theory, the following computations are necessary. Eigenvectors of J(E0wt )|β =β ∗ Let J(E0wt )|β =β ∗ = Jβ ∗ . In order to apply the method described in [6], the following computations are necessary. The matrix Jβ ∗ has a left eigenvector (associated with the zero eigenvalue) given by,

v = [v1 , v2 , v3 , v4 , v5 , v6 ], where,

v1 = 0, v2 = v2 > 0, v3 =

β ∗ (π κεω + υ G4

β ∗π G2

β ∗π ω

Acknowledgments

v5 =

One of the authors (FN) acknowledges, with thanks, the support of the University of Manitoba and the Manitoba Graduate Scholarship. ABG acknowledges the support of NSERC of Canada. The authors are grateful to the anonymous reviewers for their constructive comments.

Furthermore, the matrix Jβ ∗ has a right eigenvector (associated with the zero eigenvalue) given by,

G3 G4

v2 .

where,

w3 =

Proof. It is convenient to let

G4

w = [w1 , w2 , w3 , w4 , w5 , w6 ]T ,

w1 = −

Appendix A. Proof of Theorem 3.4

v2 , v6 =

v2 , v4 = 0,

β ∗ (μ + γ )(π ε + G2 ) − γ (σ G2 + δε) w2 , w2 = w2 > 0, μG2 (μ + γ )

ε G2

w2 , w4 =

σ G2 + δε w2 , w5 = 0, w6 = 0. G2 (μ + γ )

S = x 1 , I = x 2 , P = x 3 , R = x 4 , V = x 5 , W = x6 , Computation of bifurcation coefficients, a and b

so that the treatment-free model (3.1) can be re-written as:

It can be shown, by computing the non-zero partial derivatives of the right-hand side functions in (A.1), that the associated backward bifurcation coefficients, a and b, are given, respectively, by (see Theorem 4.1 in [6]):

dx1 = f1 =  − λwt x1 − μx1 + γ x4 , dt dx2 = f2 = dt

λwt x1 − (μ + σ + ε) x2 , 



dx3 = f3 = dt

ε x2 − μ + δ x 3 ,

dx4 = f4 = dt

σ x2 + δ x3 + ασ x5 + ηδ x6 − ψλwt x4 − μ + γ x4 ,

(A.1) 



a=

6  k,i,j=1

=

vk wi wj

∂ 2 fk wt ∗ (E , β ) ∂ xi ∂ xj 0

−2μv2 w22 β ∗

G3 G4 G22 (μ + γ )

[S1 + (εδ + G2 σ )S2 ] ,

(A.3)

66

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

where,

Following the approach in Appendix A, it can be shown that the associated Jacobian at β ∗ has the left and right eigenvectors given, respectively, by

S1 = G3 G4 (μ + γ )(G2 + ε) > 0, S2 = G3 G4 −

G1 G2 ψ(π εκω + G4 υ) , G2 + π ε

and,

b=

6 

(A.4)

w = [w1 , w2 , w3 , w4 , w5 , w6 , w7 , w8 ]T , 



vk wi

k,i=1

where,

G2 + π ε ∂ 2 fk wt ∗ v2 w2 > 0. (E , β ) = ∂ xi ∂β 0 G2

β ∗ (π K5 + π χT τ ) v , v4 = 0, (K2 K5 − ρτ ) 2 β ∗ (K6 π κεω + υ(K4 K6 − ζ φ) + π φκεχQ ) v5 = v2 , K3 (K4 K6 − ζ φ) β ∗ (π ωK6 + π χQ φ) β ∗ (π χT K2 + π ρ) v6 = v2 , v7 = v , (K4 K6 − ζ φ) (K2 K5 − ρτ ) 2 β ∗ (π χQ K4 + π ζ ω) v2 . v8 = K4 K6 − ζ φ

v1 = 0,

It follows from (A.3) that the bifurcation coefficient, a, is positive whenever

(εδ + σ G2 )

v = [v1 , v2 , v3 , v4 , v5 , v6 , v7 , v8 ], and

G1 G2 ψ(π εκω + G4 υ) > S1 + (εδ + σ G2 )G3 G4 , G2 + π ε

(A.5)

Thus, it follows, from Theorem 4.1 of [6], that the treatment-free model (3.1) (or, equivalently, (A.1)) undergoes backward bifurcation at R0 = 1 whenever Inequality (A.5) holds.

v2 = v2 > 0,

v3 =

and, Appendix B. Proof of Theorem 4.3

w1 =−

Proof. Consider the model (2.2). It is convenient to let

S = x1 , I = x2 , P = x3 , R = x4 , V = x5 , W = x6 , T = x7 , Q = x8 ,

ε K5 w , (K2 K5 − ρτ ) 2 σ (K2 K5 − ρτ ) + δεK5 + τ εθ w4 = w2 , (K2 K5 − ρτ )(μ + γ ) τε w , w8 = 0. w5 = 0, w6 = 0, w7 = (K2 K5 − ρτ ) 2 w2 = w2 > 0,

so that the treatment model (2.2) can be re-written as:

dx1 dt dx2 dt dx3 dt dx4 dt dx5 dt dx6 dt dx7 dt dx8 dt

= f1 =  − λx1 − μx1 + γ x4 , = f2 = λx1 − (μ + σ + ε) x2 , = f 3 = ε x2 + ρ x7 −





μ + δ + τ x3 ,

(B.1)

= f4 = σ x2 + δ x3 + ασ x5 + ηδ x6 + θ x7   + θ x8 − ψλx4 − μ + γ x4 ,

= f 7 = τ x3 − = f 8 = φ x6 −

=





μ + ηδ + φ x6 ,

μ + ρ + θ x7 ,

C2 = K3 (K4 K6 − ζ φ),

β(x2 + π x3 + υ x5 + π ωx6 + π χT x7 + π χQ x8 ) i=1

xi

F =

,

⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ =⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝

−β

−βπ

−βυ −βωπ

βωπ

0

0

0

0

⎞

βχQ π ⎟ ⎟

−μ − γ

ασ

ηδ

0

0

0

−K3

0

0

0

0

0

0

0

κε

−K4

0

ζ

0

0

τ

0

0

0

−K5

0

0

0

0

0

0

φ

0

−K6

0

βυ

−βχQ π

0

0

0

−βχT π

βχT π ρ θ

0

β − K1 βπ ε −K2 σ δ

γ

θ

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

Consider the case of the model (B.1) with RT = 1. Suppose, also, that β is chosen as the bifurcation parameter. Solving for β from RT = 1 gives

β∗ =

σ (K2 K5 − ρτ ) + δ K5 ε + τ εθ ,

and (noting that K2 K5 − ρτ > 0 and K4 K6 − ζ φ > 0),

b=

8  k,i=1

−μ

K1 (K2 K5 − ρτ ) . K2 K5 − ρτ + π ε K5 + π χT τ ε

(B.4)

C3 = (K4 K6 − ζ φ)υ + K6 εκπ ω + φεκπ χQ ,

and f = [f1 , · · · , f8 ]T represents the vector field of the model (2.2). Evaluating the Jacobian of the system (B.1) at the DFE (E0T ) gives: J(E0T ) ⎛

  −2μv2 w22 K1 C + C2 F − β ∗ ψ FC3 , K3 (K4 K6 − ζ φ)(K2 K5 − ρτ )(μ + γ ) 1

C1 = K3 (K4 K6 − ζ φ)(μ + γ )(K2 K5 − ρτ + ε K5 + τ ε),

μ + ζ + θ x8 ,

8

∂ 2 fk T ∗ (E , β ), ∂ xi ∂ xj 0

where,





vk wi wj

(B.3)





8  k,i,j=1

= f5 = ψλx4 − (μ + ασ + κε) x5 , = f6 = κε x5 + ζ x8 −

w3 =

Furthermore, it can be shown (by computing the associated nonzero partial derivatives of (B.1) at the DFE (E0T )) that the associated backward bifurcation coefficients, a and b, are given, respectively, by:

a=

where,

λ=

β ∗ (μ+γ )(K2 K5 −ρτ +π εK5 +π χT τ ε)−γ (σ (K2 K5 −ρτ )+δεK5 +τ εθ ) w2 , μ(μ+γ )(K2 K5 −ρτ )

(B.2)

=

vk wi

∂ 2 fk T ∗ (E , β ) ∂ xi ∂β 0

  K2 K5 − ρτ + π ε K5 + π χT τ ε v2 w2 > 0. (K2 K5 − ρτ )

It follows from (B.3), with (B.4), that the bifurcation coefficient, a, is positive whenever

J 2 > J1 ,

(B.5)

where,

J1 = K3 (K4 K6 − ζ φ)F + K3 (μ + γ )(K4 K6 − ζ φ) × (K2 K5 − ρτ + ε K5 + τ ε), J2 =

β ∗ ψ [(K4 K6 − ζ φ)υ + K6 εκπ ω + φεκπ χQ ]F.

Thus, it follows, from Theorem 4.1 of [6], that the treatment model (2.2) (or, equivalently, (B.1)) undergoes backward bifurcation at RT = 1 whenever Inequality (B.5) holds.

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

B.1. Non-existence of backward bifurcation

where,

Consider backward bifurcation coefficient a in (B.3). Assume κ = η = 1, φ = τ , ρ = ζ and χT = χQ . Then, it follows from the definitions in Section 4.1.1 that K3 ≥ K1 , K2 K5 − ρτ = K4 K6 − ζ φ . Consequently, the expression C1 + C2 F − β ∗ ψ FC3 given in (B.3), can be re-written as:

C1 + C2 F − β ∗ ψ FC3 = K3 (K4 K6 − ζ φ)[(μ + γ )(K2 K5 − ρτ + ε K5 + τ ε) + F]

x = [I (t), P (t), W (t), T (t), Q (t)]T , ⎛

β

⎜0 ⎜ ⎜ F1 = ⎜ ⎜0 ⎜ ⎝0 ⎛

− β ∗ ψ [(K4 K6 − ζ φ)υ + K6 επ ω + φεπ χT ]F. Substituting β ∗ by its definition in (B.2), and defining N = (μ + γ )(K2 K5 − ρτ + εK5 + τ ε), gives

K3 (K4 K6 − ζ φ)(K2 K5 − ρτ + π ε K5 + τ επ χT )(N + F ) K2 K5 − ρτ + π ε K5 + τ επ χT

K1 (K4 K6 − ζ φ)(K2 K5 − ρτ + K5 επ + τ επ χT )F ≥ 0, K2 K5 − ρτ + π ε K5 + τ επ χT

Proof. It should, first of all, be mentioned that the system (C.1) satisfies the Type K condition [39] (hence, Comparison Theorem can be used [28]). Furthermore, consider the following reduced model (C.1) of the treatment model (2.2) for the special case with ψ = 0, given by:

=

0

0

0

0

0

0

0

0

0

0

0

0

0

0



μ + δ + τ P,

(C.1) 

0

0

0

−ρ

0

K4

0

−τ

0

K5

0 ⎟ ⎟ ⎟ −ζ ⎟ ⎟, ⎟ 0 ⎠

0

−φ

0

K6

⎛

β

⎜   ⎜0 μS(t) ⎜ ⎜0 J = 1−  ⎜ ⎜ ⎝0

πβ

ωπ β

χT π β

χQ π β

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

⎞

⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎠

S(t) ≤ N(t) ≤

 in DT . μ

Hence, it follows from (C.2) that

d x ≤ (F1 − H1 )x. dt

(C.3)

Since RT |ψ =0 = ρ(F1 H1−1 ) = RT ≤ 1 (or, equivalently, the eigenvalues of the matrix F1 − H1 all have negative real parts), it follows that the linearized differential inequality system (C.3) is stable whenever RT < 1. Thus, it follows, by Comparison Theorem [28], that

lim (I (t), P (t), W (t), T (t), Q (t)) → (0, 0, 0, 0, 0).

t→∞

Substituting I (t) = P (t) = W (t) = T (t) = Q (t) = 0 into the (C.1) and using the fact that V (t) → 0 as t → ∞ show that S(t) → S∗ as t → ∞ (for RT < 1). Therefore,

Proof. Consider the reduced model (2.5). For this model [13], the DFE is given by







=

τ P − μ + ρ + θ T,

=

φ W − μ + ζ + θ Q.

E¯0 = (S¯ ∗ , I¯∗ , P¯ ∗ , R¯ ∗ , T¯ ∗ ) =







 , 0, 0, 0, 0 . μ

¯c Furthermore, the associated reproduction number (denoted by R [13]) is given by

It follows from (C.1) that V (t) → 0 as t → ∞. Hence, the equation for dV can be temporal removed from (C.1). The infected components of dt the model (C.1) (with the equation for as:

Hence, the DFE (E0T ) of the reduced model (C.1) is GAS in DT whenever RT < 1. Appendix D. Proof of non-existence of backward bifurcation in the reduced model (2.5)

κε V + ζ Q − μ + ηδ + φ W,

d x = (F1 − H1 )x − Jx, dt

0 K2



σ I + δ P + ασ V + ηδ W + θ T + θ Q − μ + γ R,



0

⎞

= E0T |ψ =0 .

λS − (μ + σ + ε) I,



0

⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎠

t→∞

= −(μ + ασ + κε)V, =

0

lim (S(t), I (t), P (t), V (t), W (t), T (t), Q (t)) → (S∗ , 0, 0, 0, 0, 0, 0)

=  − λS − μS + γ R,

= εI + ρ T −

0

⎞

Thus, J is a non-negative matrix since

Appendix C. Proof of Theorem 4.4



χQ π β

0

since K3 ≥ K1 and N ≥ 0. Thus, under the aforementioned assumptions (κ = η = 1, φ = τ , ρ = ζ and χT = χQ ) the associated backward bifurcation coefficient, a, is negative. Hence, the treatment model (2.2) does not undergo backward bifurcation in this case. In other words, these analyses show that the requirements υ = 1, ω = 1 and α = 1 in Section 4.1.3 can be relaxed (i.e., backward bifurcation does not exist in the model (2.2) even if υ = 1, ω = 1 and α = 1).

=

χT π β

(with K1 , K2 , K4 , K5 and K6 as defined in Section 4.1) and,

K1 (K2 K5 − ρτ )ψ [(K4 K6 − ζ φ)υ + K6 επ ω + φεπ χT ]F − K2 K5 − ρτ + π ε K5 + τ επ χT

dS dt dI dt dP dt dR dt dV dt dW dt dT dt dQ dt

ωπ β

0

K3 (K4 K6 − ζ φ)(K2 K5 − ρτ + π ε K5 + τ επ χT )(N + F ) = K2 K5 − ρτ + π ε K5 + τ επ χT

−

πβ

K1 ⎜−ε ⎜ ⎜ H1 = ⎜ ⎜0 ⎜ ⎝0

C1 + C2 F − β ∗ ψ FC3

≥

67

dV dt

removed) can be re-written

(C.2)

¯c = R

β(K¯ 2 K¯ 3 − ρτ + π ε K¯ 3 + π χT τ ε) , K¯ 1 (K¯ 2 K¯ 3 − ρτ )

where K¯ 1 = μ + σ + ε , K¯ 2 = μ + δ + τ , K¯ 3 = μ + ρ + θ . It can be shown, as in Appendix A, that the associated backward bifurcation coefficients, a and b, are given, respectively, by (see

68

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

Theorem 4.1 in [6]):

a¯ =

Appendix E. Proof of Theorem 4.5

¯ 22 −2μv¯ 2 w

[K¯ 1 (K¯ 2 K¯ 3 − ρτ )(μ + γ ) (K¯ 2 K¯ 3 − ρτ )2 (μ + γ ) × (K¯ 2 K¯ 3 − ρτ + ε K¯ 3 + τ ε) + K¯ 1 (K¯ 2 K¯ 3 − ρτ )F¯ ¯ − β¯ ∗ ψ(K¯ 2 K¯ 3 − ρτ + π ε K¯ 3 + π χT τ ε)F],

(D.1)

Proof. Consider the treatment model (2.2) for the special case υ = ω = χQ = 0. It follows, from the next generation operator method [44] ˆ T (the associated control reproduction number of the treatment that R ˆ ) is given by: model (2.2) with λ = λ





ˆ T = ρ Fˆ H ˆ −1 = R

where,

K¯ 1 (K¯ 2 K¯ 3 − ρτ ) , ¯ ¯ K2 K3 − ρτ + π ε K¯ 3 + π χT τ ε F¯ = σ (K¯ 2 K¯ 3 − ρτ ) + δ K¯ 3 ε + τ εθ ,   K¯ 2 K¯ 3 − ρτ + π ε K¯ 3 + π χT τ ε ¯ 2 > 0, v¯ 2 w b¯ = (K¯ 2 K¯ 3 − ρτ )

(E.1)

β¯ ∗ =

(D.2)

−μ ⎜ ⎜ 0 ⎜ ⎜ J¯|β¯ ∗ (E¯0 ) = ⎜ 0 ⎜ ⎜ 0 ⎝

−β¯ ∗

β¯ ∗ − K¯ 1 ε σ

0

0

−β¯ ∗ π β¯ ∗ π

γ 0

−K¯ 2

0

δ τ

−μ − γ 0

ˆ (of the transiwhere, the matrices Fˆ (of new infection terms) and H tion terms) evaluated at the DFE (E0T ) are given, respectively, by:

⎛

β

(D.3)

¯ as the respective left and right eigenvectors correspondwith v¯ and w ing to zero eigenvalue of the Jacobian of the system (2.5), evaluated at the associated disease-free equilibrium (E¯0 ), given by:

⎛

⎞ −β¯ ∗ χT π ⎟ β¯ ∗ χT π ⎟ ⎟ ⎟ , ρ ⎟ ⎟ ⎟ θ ⎠ −K¯ 3

with v¯ = [v¯ 1 , v¯ 2 , v¯ 3 , v¯ 4 , v¯ 5 ], where,

β ∗ (π K¯ 3 + π χT τ ) v¯ 1 = 0, v¯ 2 = v¯ 2 > 0, v¯ 3 = v¯ 2 , v¯ 4 = 0, (K¯ 2 K¯ 3 − ρτ ) β¯ ∗ (π χT K¯ 2 + π ρ) v¯ 5 = v¯ 2 , (K¯ 2 K¯ 5 − ρτ )

πβ

χT π β

Fˆ = ⎝ 0

0

0

0

0

0

⎜

⎞

⎛

⎟ ⎠,

and

0

ˆ = ⎝−ε H

K2

⎟ 0 ⎠,

−τ

K5

⎜

0

0

⎞

K1

with, K1 , K2 and K5 as defined in Section 4.1. The proof is based on the Comparison Theorem [28]. First of all, as in Appendix C, it should be mentioned that the equations of the inˆ , satisfies fected components of the treatment model (2.2), with λ = λ the Type K condition [39] (hence, Comparison Theorem can be used [28]). Furthermore, the infected components of the treatment model ˆ , can be re-written as: (2.2), with λ = λ

d ˆ )x − Jˆx, x = (Fˆ − H dt

(E.2)

where,

x = [I (t), P (t), T (t)]T , and,

 Jˆ = 1 −

⎛

 β μS(t) ⎜ ⎝0  0

¯ = [w ¯ 1, w ¯ 2, w ¯ 3, w ¯ 4, w ¯ 5 ]T , with, and w

¯1 = − w

β [(K2 K5 − ρτ ) + π ε K5 + χT π ετ ] = RT , K1 (K2 K5 − ρτ )

πβ

χT π β

0

0

0

⎞ ⎟ ⎠.

0

β¯ ∗ (μ + γ )(K¯ 2 K¯ 3 − ρτ + π ε K¯ 3 + π χT τ ε) − γ [σ (K¯ 2 K¯ 3 − ρτ ) + δε K¯ 3 + τ εθ ] ¯ 2, w μ(μ + γ )(K¯ 2 K¯ 3 − ρτ ) Thus, Jˆ is a non-negative matrix since

ε K¯ 3 ¯ 2 > 0, w ¯3 = ¯2 = w ¯ 2, w w (K¯ 2 K¯ 3 − ρτ ) σ (K¯ 2 K¯ 3 − ρτ ) + δε K¯ 3 + τ εθ ¯4 = ¯ 2, w w (K¯ 2 K¯ 3 − ρτ )(μ + γ ) ¯5 = w

τε

(K¯ 2 K¯ 3 − ρτ )

¯ 2. w

Substituting (D.2) into (D.1) gives:

a¯ =

¯ 22 −2μv¯ 2 w

(K¯ 2 K¯ 3 − ρτ )2 (μ + γ )

[K¯ 1 (K¯ 2 K¯ 3 − ρτ )(μ + γ )

¯ − ψ K¯ 1 (K¯ 2 K¯ 3 − ρτ )F], ¯ 22 −2μv¯ 2 w

(K¯ 2 K¯ 3 − ρτ )2 (μ + γ )

 in DT . μ

Hence, it follows from (E.2) that

d ˆ )x. x ≤ (Fˆ − H dt

(E.3)

ˆ T = ρ(Fˆ H ˆ −1 ) = RT ≤ 1 (or, equivalently, the eigenvalues of Since R ˆ ˆ ) all have negative real-parts), it follows that the the matrix (F − H linearized differential inequality system (E.3) is stable whenever RT < 1. Thus, it follows, by Comparison Theorem [39], that

lim (I (t), P (t), T (t)) → (0, 0, 0).

t→∞

× (K¯ 2 K¯ 3 − ρτ + ε K¯ 3 + τ ε) + K¯ 1 (K¯ 2 K¯ 3 − ρτ )F¯

=

S(t) ≤ N(t) ≤

[K¯ 1 (K¯ 2 K¯ 3 − ρτ )(μ + γ )

¯ × (K¯ 2 K¯ 3 − ρτ + ε K¯ 3 + τ ε) + K¯ 1 (K¯ 2 K¯ 3 − ρτ )(1 − ψ)F]. Since ψ < 1, the bifurcation coefficient a¯ is always negative. Furthermore, the bifurcation coefficient b¯ is always positive. Thus, it follows from Item (iv) of Theorem 4.1 in [6], that the reduced model (2.5) does not undergo backward bifurcation in this case.

ˆ = 0. It should be noted that setting I (t) = P (t) = T (t) = 0 implies λ ˆ = 0 into the equation for dv in (2.2) gives V (t) → Hence, substituting λ dt 0, as t → ∞. Substituting (I (t), P (t), V (t), T (t)) = (0, 0, 0, 0) into the model (2.2) gives the following system of linear equations

dW = dt dQ = dt





ζ Q − μ + ηδ + φ W, 

(E.4)



φ W − μ + ζ + θ Q,

from which it is clear that the system (E.4) has only one equilibrium, (W ∗ , Q ∗ ) = (0, 0), whichis stable (hence, globally-asymptotically

F. Nazari et al. / Mathematical Biosciences 263 (2015) 51–69

stable, since the system is linear and eigenvalues of the associated Jacobin evaluated at (0, 0) have negative real part). Thus, (W (t), Q (t)) → (0, 0), as t → ∞. Finally, substituting

(I (t), P(t), V (t), W (t), T (t), Q (t)) = (0, 0, 0, 0, 0, 0), into the equations for that R(t) → 0 and S →

dR dt S∗

and

dS dt

in the treatment model (2.2) shows

=  μ , as t → ∞ when RT < 1. Therefore,

lim (S(t), I (t), P (t), R(t), V (t), W (t), T (t), Q (t))

t→∞

→ (S∗ , 0, 0, 0, 0, 0, 0, 0) = E0T |υ =ω=χQ =0 . ˆ , is GAS Hence, the DFE (E0T ), of the treatment model (2.2), with λ = λ in DT whenever RT < 1. References [1] M.J. Alter, Epidemiology of hepatitis C virus infection, World J. Gastroenterol. 13(17) (2007) 2436–2441. [2] D. Amarapurkar, Natural history of hepatitis C virus infection, J. Gastroenterol. Hepatol. 15 (2000) E105–E110. [3] S.M. Blower, H. Dowlatabadi, Sensitivity and uncertainty analysis of complex models of disease transmission, Int. Stat. Rev. 62(2) (1994) 229–243. [4] B. Boldin, Introducing a population into a steady community: the critical case, the centre manifold and the direction of bifurcation, SIAM J. Appl. Math. 66(4) (2006) 1424–1453. [5] J. Carr, Applications of Centre Manifold Theory, Springer-Verlag, New York, 1981. [6] C. Castillo-Chavez, B. Song, Dynamical models of tuberculosis and their applications, Math. Biosci. Eng. 2 (2004) 361–404. [7] C. Castillo-Chavez, K. Cooke, W. Huang, S.A. Levin, Results on the dynamics for models for the sexual transmission of the human immunodeficiency virus, Appl. Math. Lett. 2 (1989) 327–331. [8] S. Corson, D. Greenhalgh, Mathematical modelling the spread of hepatitis C in injecting drug users, Math. Med. Biol. 29(3) (2012) 205–230. [9] S. Corson, D. Greenhalgh, S.J. Hutchinson, A time since onset of injection model for hepatitis C spread amongst injecting drug users, J. Math. Biol. 66(4–5) (2013) 935–978. [10] A.M.D. Bisceglie, Natural history of hepatitis C: its impact on clinical management, Hepatology 31 (2000) 1014–1018. [11] O. Diekmann, J. Heesterbeek, J. Metz, On the definition and computation of the basic reproduction ratio R0 in models for infectious disease in heterogeneous population, Math. Biosci. 180 (2002) 29–48. [12] J. Dushoff, W. Huang, C. Castillo-Chavez, Backward bifurcations and catastrophe in simple models of fatal diseases, J. Math. Biol. 36 (1998) 227–248. [13] E.H. Elbasha, Model for hepatitis C virus transmissions, Math. Biosci. Eng. 10(4) (2013) 1045–1065. [14] E.H. Elbasha, A.B. Gumel, Theoretical assessment of public health impact of imperfect prophylactic HIV-1 vaccines with therapeutics benefits, Bull. Math. Biol. 68 (2006) 577–614. [15] P. Ferenci, K.R. Reddy, Impact of HCV protease-inhibitor-based triple therapy for chronic HCV genotype 1 infection. Antivir Ther. 16 (2011) 1187–1201. [16] S.M. Garba, A.B. Gumel, M.R.A. Bakar, Backward bifurcation in dengue transmission dynamics, Math. Biosci. 215 (2008) 11–25. [17] M.G. Ghany, D.B. Strader, D.L. Thomas, L.B. Seek, Diagnosis, management, and treatment of hepatitis C: an update, Hepatology 49 (2009) 1335–1374. [18] J. Grebely, M. Prins, M. Hellard, A.L. Cox, W.O. Osburn, G. Lauer, K. Page, A.R. Lloyd, G.J. Dore, Hepatitis c virus clearance, re-infection, and persistence, with insights from studies of injecting drug users: towards a vaccine, Lancet Infect. Dis. 12 (2012) 408–414. [19] D.D. Greenhalgh, O. Diekmann, M.C.M. de Jong, Subcritical endemic steady states in mathematical models for animal infections with incomplete immunity, Math. Biosci. 165(1) (2000) 1–25. [20] A.B. Gumel, Causes of backward bifurcations in some epidemiological models, J. Math. Anal. Appl. 395(1) (2012) 355–365. [21] A.B. Gumel, B. Song, Existence of multiple-stable equilibria for a multi- drugresistant model of mycobacterium tuberculosis, Math. Biosci. Eng. 5 (2008) 437–455. [22] A.B. Gumel, C.C. McCluskey, P. van den Driessche, Mathematical study of a stagedprogression HIV model with imperfect vaccine, Bull. Math. Biol. 68 (2006) 2105–2128. [23] H.M. Hethcote, The mathematics of infectious diseases, SIAM Rev. 42 (2000) 599–963.

69

[24] J.H. Hoofnagle, Hepatitis C: the clinical spectrum of disease, Hepatology 26(3) (1997) 15S–20S. [25] S. Kamal, Acute hepatitis C: a systematic review, Am. J. Gastroenterol. 103 (2008) 1283–1297. [26] K. Kiyosawa, E. Tanaka, T. Sodeyama, K. Yoshizawa, K. Yabu, K. Furuta, H. Imai, Y. Nakano, S. Usuda, K. Uemura, Transmission of hepatitis C in an isolated area in Japan: community-acquired infection. The South Kiso Hepatitis Study Group, J. Gastroenterol. 106(6) (1994) 1596–1602. [27] M. Kretzschmar, L. Wiessing, Modeling the transmission of hepatitis C in injecting drug users, European Monitoring Centre for Drug and Drug Addiction Monograph. [28] V. Lakshmikantham, S. Leela, A.A. Martynyuk, Stability Analysis of Nonlinear Systems. Marcel Dekker Inc., New York and Basel, 1989. [29] C.L. Leen, Hepatitis C and HIV co-infection, Int. J. STD AIDS. 15(5) (2004) 289– 295. [30] S. Marino, I.B. Hogue, C.J. Ray, D.E. Kirschner, A methodology for performing global uncertainty and sensitivity analysis in systems biology, J. Theo. Biol. 254(1) (2008) 178–196. [31] N.K. Martin, P. Vickerman, M. Hickman, Mathematical modelling of hepatitis C treatment for injecting drug users, J. Theo. Biol. 274(1) (2011) 58–66. [32] N.K. Martin, P. Vickerman, G.R. Foster, S.J. Hutchinson, D.J. Goldberg, M. Hickman, Can antiviral therapy for hepatitis C reduce the prevalence of HCV among injecting drug user populations? a modelling analysis of its prevention utility, J. Hepatol. 54 (2011) 1137–1144. [33] R.G. McLeod, A.B. Gumel, D.A. Slonowsky, Sensitivity and uncertainty analysis for a SARS model with time-varying inputs and outputs, Math. Biosci. Eng. 3(3) (2006) 527–544. [34] S.H. Mehta, B.L. Genberg, J. Astemborski, R. Kavasery, G.D. Kirk, D. Vlahov, S.A. Strathdee, D.L. Thomas, Limited uptake of hepatitis C treatment among injection drug users, J. Commun. Health 33 (2008) 126–133. [35] S. Noguchi, M. Sata, H. Suzuki, M. Mizokami, K. Tanikawa, Routes of transmission of hepatitis C virus in an endemic rural area of Japan; molecular epidemiologic study of hepatitis C virus infection, Scand. J. Infect. Dis. 29 (1997) 23– 28. [36] M.A. Sanchez, S.M. Blower, Uncertainty and sensitivity analysis of the basic reproductive rate, Am. J. Epidemiol. 145(12) (1997) 1127–1137. [37] O. Sharomi, C.N. Podder, A.B. Gumel, E. Elbasha, J. Watmough, Role of incidence function in vaccine-induced backward bifurcation in some HIV models, Math. Biosci. 210(2) (2007) 436–463. [38] C.W. Shepard, L. Finelli, M.J. Alter, Global epidemiology of hepatitis C virus infection, Lancet. Infect. Dis. 5(9) (2005) 558–567. [39] H.L. Smith, Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems, in: Mathematical Surveys and Monographs, vol. 41, American Mathematical Society, Providence, 1995. [40] A.J. Sutton, N.J. Gray, W.J. Edmunds, V.D. Hope, O.N. Gill, M. Hickman, Modelling the force of infection for hepatitis B and hepatitis C in injecting drug users in England and Wales, BMC Infect. Dis. 6 (2006) 93–102. [41] M.J. Sweeting, V.D. Hope, M. Hickman, J.V. Parry, F. Ncube, M.E. Ramsay, D. De Angelis, Hepatitis c infection among injecting drug users in England and Wales (1992–2006): there and back again?, Am. J. Epidemiol. 170(3) (2009) 352–360. [42] B.J. Thomson, R.G. Finch, Hepatitis C virus infection. Clin. Microbiol. Infect. Dis. 11(2) (2005) 83–165. [43] C.V.D. Berg, B. Colette, S.G.V. Brussel, R. Coutinho, M. Prins, Full participation in harm reduction program is associated with decreased risk for human immunodeficiency virus and hepatitis C virus, J. Addic. 102(9) (2007) 1454–1462. [44] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (2002) 29–48. [45] P. Vickerman, M. Hickman, A. Judd, Modeling the impact of hepatitis C transmission of reducing syringe sharing: London case study, Int. J. Epidemiol. 36(2) (2007) 396–405. [46] P. Vickerman, L. Platt, S. Hawkes, Modeling the transmission of HIV and HCV among injecting drug users in Rawalpindi, a low HCV prevalence setting in Pakistan, Sex. Transm. Infect. 85 (2009) ii23–ii30. [47] M.L. Volk, R. Tocco, S. Saini, A.S. Lok, Public health impact of antiviral therapy for hepatitis C in the united states, J. Hepatol. 50 (2009) 1750–1755. [48] A. Wasley, M.J. Alter, Epidemiology of hepatitis C: geographic differences and temporal trends, Semin. Liver. Dis. 20 (2000) 1–16. [49] World Health Organization (WHO), \Facts Sheet: Hepatitis C," http://www. who.int/mediacentre/factsheets/fs164/en/index.html (accessed 29-10-12). [50] J. Yuan, Z. Yang, Global dynamics of an SEI model with acute and chronic stages, J. Comput. App. Math. 213(2) (2008) 465–476. [51] C.K. Zaleta, J. Velasco-Hernandez, A simple vaccination model with multiple endemic state, Math. Biosci. 164 (2000) 183–201. [52] S. Zhang, Y. Zhou, Dynamics and application of an epidemiological model for hepatitis C, J. Math. Comp. Biol. 56(1–2) (2012) 36–42.