# A Mathematical Model for Assessing the Impact of Intravenous Drug Misuse on the Dynamics of HIV and HCV within Correctional Institutions

Authors: S. Mushayabasa; Claver P. Bhunu; Alexander G. R. Stewart
Journal: ISRN Biomathematics
 (2012)

---

## Abstract

Unsafe injecting practices, blood exchange, the use of nonsterile needles, and other cutting instruments for tattooing are common in correctional institutions, resulting in a number of blood transmitted infections. A mathematical model for assessing the dynamics of HCV and HIV coinfection within correctional institutions is proposed and comprehensively analyzed. The HCV-only and HIV-only submodels are first considered. Analytical expressions for the threshold parameter in each submodel and the cointeraction are derived. Global dynamics of this coinfection shows that whenever the threshold parameter for the respective submodels and the coinfection model is less than unity, then the epidemics die out, the reverse condition implies disease persistence within correctional institutions. Numerical simulations using a set of plausible parameter values are provided to support analytical findings.


---

## Body

## 1. Introduction

Prior studies suggests that prevalence of both human immunodeficiency virus (HIV) and hepatitis C virus (HCV) infections is known to be higher among incarcerated populations than the general population, since a high proportion of incarcerated individuals originate from high-risk environments and have high-risk behaviors, especially drug use [1–6]. HCV and HIV coinfections represent a public health problem of growing importance; because of similar modes of spread, many people are coinfected with HCV and HIV or HCV and HBV and in some cases with all three viruses at the same time [5]. In particular, HCV/HIV coinfections are common and are known as “twin epidemics” [6, 7]. Both HIV and HIV are blood born RNA (ribonucleic acid) viruses that replicate rapidly. Unsafe injecting practices, blood exchange and the use of non-sterile needles are the most efficient means of transmitting both viruses [6, 8].Correctional institutions host a disproportionately high prevalence of HCV infection and coinfections. The prevalence of HCV among prisoners approaches 57.5% and far exceeds that of HIV in prison [9, 10]. Transmission of these infections is believed to be a rare consequence of blood or body fluid exposures in the prison. However, if transmission does occur, the consequences are permanent and potentially fatal [5–7]. Coinfection with the two viruses (HCV/HIV) is associated with an accelerated course of hepatitis C disease [5–7]. Prison populations constitute a very high-risk group; they have high levels of HCV infection and HIV or HBV coinfection [11, 12]. HCV positive inmates are at exceptional risk for coinfection with HIV because of the association of injectable substance abuse [9–12].Mathematical models have become invaluable management tools for epidemiologists, both shedding light on the mechanisms underlying the observed dynamics as well as making quantitative predictions on the effectiveness of different control measures. The literature and development of mathematical epidemiology are well documented and can be found in [13–15]. Mathematical modeling provides an alternative means to define our problems, organize our thoughts, understand our data, communicate and test our understanding, and make predictions. The deterministic compartmental model provides means of obtaining insight into on the dynamics of HIV and HCV within correctional institutions. As with most models for disease transmission and control, our model is based on the simple SIR model [16]. The main parameter of the SIR model is the basic reproduction number, ℛ0. The basic reproductive number provides a quantitative framework to address the question of level of risk posed by a disease to a population. It allows the identification, first, of what needs to be known in order to assess risk, and, second, which of these factors are most important in determining the magnitude of risk. If this parameter (ℛ0) is below unity, then the disease dies out, whereas if this parameter is above unity, any small introduction of infected individuals in the population results in an oscillatory approach to an endemic equilibrium. Mathematically, there is a trivial equilibrium, known as the disease-free equilibrium, which is globally asymptotically stable whenever ℛ0<1 [17]. A brief survey reveals that a number of statistical and mathematical studies for assessing the transmission dynamics of infectious diseases within correctional institutions have been proposed (see [10, 18–24]) to mention a few. This paper seeks to use a mathematical model to asses the impact intravenous drug misuse on the dynamics of HIV and HCV within correctional institutions.The paper is structured as follows. The HIV-HCV transmission model is formulated in the next section. The HCV submodel is presented in Section3, followed by the HIV submodel. Section 5 presents analytical results for the HIV-HCV model formulated in Section 2. Simulation results and projection profiles of HIV-HCV model are presented in Section 6. Summary and concluding remarks round up the paper.


## 2. Model Formulation

Based on the individual’s epidemiological status, the total populationN has been subdivided into the following classes or subgroups: susceptible nonintravenous drug users (Sn), susceptible intravenous drug users (IDUs) (Sd), individuals singly infected with HCV (Ic) (both infective and chronic), individuals singly infected with HIV (Ih), and dually infected individuals (Ich).  We assume that AIDS patients dually or singly infected are not participants due to their health status. Prisoners are recruited at rate Λ. Although there is need for real demographic data for one to know the proportion of inmates recruited into each of the aforementioned epidemiological classes, we have assumed that a fraction π0,  π1, π2, π3, and π4 are recruited into Sn, Sd, Ic, Ih, and Ich, respectively. Natural mortality μ is assumed to be constant in all classes. Susceptible intravenous drug users acquire either HCV or HIV but not both at rates λc, λh, respectively, where λc=β(Ic+ηIch)/N, and λh=θ(Ih+ηIch)/N. Parameters β and θ denotes the probability of getting infected with either HCV or HIV, respectively. The parameter η>1 captures the assumed increased probability for individuals dually infected with HCV and HIV to infect their partners. The model takes the following form:
(1)Sn′=Λπ0-(α+μ+ω)Sn,Sd′=Λπ1+αSn-(λc+λh)Sd-(μ+ω)Sd,Ic′=Λπ2+λcSd-σλhIc-(μ+ω)Ic,Ih′=Λπ3+λhSd-σλcIh-(μ+ω+δ)Ih,Ich′=Λπ4+σ(λhIc+λcIh)-(μ+ω+ϕδ)Ich.
Upon completion of their sentence, inmates are released at rate ω, susceptible non-IDUs acquire IDUS behavior at rate α,  and σ(σ>1) captures the assumed increased likelihood for singly infected individuals to become dually infected due to suppressed immune system. Individuals infected with HIV alone and those dually infected with HIV and HCV progress to AIDS stage at rate δ, and ϕδ, respectively. The modification factor ϕ>1 captures the increased likelihood for dually infected individuals to progress to AIDS stage compared to individuals singly infected with HIV. The model flow diagram is depicted in Figure 1.Figure 1
Flow diagram of the HCV transmission.


### 2.1. Basic Properties of the Model

In this section, we study the basic properties of the solutions of model system (1), which are essential in the proofs of stability.Lemma 1.
The equations preserve positivity of solutions.Proof.
The vector field given by the right-hand side of (1) points inward on the boundary of ℝ+5∖{0}. For example, if Sn=0, then, Sn′=Λπ0≥0. In an analogous manner, the same result can be shown for the other model components (variables).Lemma 2.
All solutions of system (1) are bounded.Proof.
Using system (1) we have N′=Λ-(μ+ω)N-δ(Ih+ϕIch)≤Λ-(μ+ω)N. Assume that N(t)≤M for all t≥0 where M=(Λ/(μ+ω))+ϵ,ϵ>0. Suppose the assumption is not true then there exists a t1>0 such that
(2)N(t1)=Λ(μ+ω)+ϵ,N(t)<Λ(μ+ω)+ϵ,t<t1,N′(t1)≥0,N′(t1)≤Λ-(μ+ω)N(t1)=-(μ+ω)<0,
which is a contradiction meaning the assumption is true. This means N(t)≤M for all t≥0.Therefore all feasible solutions of system (1) enter the region
(3)Ω={(Sn,Sd,Ic,Ih,Ich)∈ℝ+5:N≤Λ(μ+ω)}.
Thus, Ω is positively invariant and it is sufficient to consider solutions of system (1) in Ω. Existence, uniqueness, and continuation results for system (1) hold in this region and all solutions of system (1) starting in Ω remain in Ω for all t≥0. All parameters and state variables for model system (1) are assumed to be nonnegative (for biological relevance) for all t≥0 since it monitors human population.


## 2.1. Basic Properties of the Model

In this section, we study the basic properties of the solutions of model system (1), which are essential in the proofs of stability.Lemma 1.
The equations preserve positivity of solutions.Proof.
The vector field given by the right-hand side of (1) points inward on the boundary of ℝ+5∖{0}. For example, if Sn=0, then, Sn′=Λπ0≥0. In an analogous manner, the same result can be shown for the other model components (variables).Lemma 2.
All solutions of system (1) are bounded.Proof.
Using system (1) we have N′=Λ-(μ+ω)N-δ(Ih+ϕIch)≤Λ-(μ+ω)N. Assume that N(t)≤M for all t≥0 where M=(Λ/(μ+ω))+ϵ,ϵ>0. Suppose the assumption is not true then there exists a t1>0 such that
(2)N(t1)=Λ(μ+ω)+ϵ,N(t)<Λ(μ+ω)+ϵ,t<t1,N′(t1)≥0,N′(t1)≤Λ-(μ+ω)N(t1)=-(μ+ω)<0,
which is a contradiction meaning the assumption is true. This means N(t)≤M for all t≥0.Therefore all feasible solutions of system (1) enter the region
(3)Ω={(Sn,Sd,Ic,Ih,Ich)∈ℝ+5:N≤Λ(μ+ω)}.
Thus, Ω is positively invariant and it is sufficient to consider solutions of system (1) in Ω. Existence, uniqueness, and continuation results for system (1) hold in this region and all solutions of system (1) starting in Ω remain in Ω for all t≥0. All parameters and state variables for model system (1) are assumed to be nonnegative (for biological relevance) for all t≥0 since it monitors human population.


## 3. HCV-Only Submodel

Before analyzing the full model (system (1)), it is instructive to gain insights on the dynamics of the HCV-only submodel, obtained by setting Ih=Ich=0, so that system (1) reduces to
(4)Sn′=Λπ0-(α+μ+ω)Sn,Sd′=Λπ1+αSn-(λc+μ+ω)Sd,Ic′=Λπ2+λcSd-(μ+ω)Ic.
For system (4), the first octant in the state space is positively invariant and attracting; that is, solutions that start where all the variables are nonnegative remain there. Thus, system (4) will be analyzed in a suitable region
(5)Ωc={(Sn,Sd,Ic)∈ℝ+3:N≤Λ(μ+ω),
which is positively invariant and attracting. Existence, uniqueness, and continuation results for system (4) hold in this region.


### 3.1. Disease-Free Equilibrium and Its Stability

System (4) has a disease-free equilibrium (DFE) given by
(6)ℰc0=(Sn0,Sd0,Ic0)=(Λπ0(α+μ+ω),Λ[απ0+π1(μ+α+ω)](μ+ω)(μ+α+ω),0).
Now, we compute the basic reproductive number, ℛ0, for system (4). Following Van Den Driessche and Watmough [17], the reproductive number for system (4) is given by
(7)ℛc=β[απ0+(α+μ+ω)π1](π0+π1)(μ+ω)(α+μ+ω),whereℛc measures the average number of new secondary cases generated by a single HCV infective individual during his/her entire infectious period when he/she is introduced into a susceptible population within a correctional setting in the absence of HCV intervention strategies. Using [17, Theorem  2], the following result is established.Theorem 3.
ℰ
c
0 is locally asymptotically stable (LAS) if ℛc<1, and unstable if ℛc>1.In order to examine the global stability ofℰc0 we use the method proposed by Kamgang and Sallet (2008). By closely following Kamgang and Sallet (2008) [25–27], we write system (1) in the form
(8)x1′=A1(x)·(x1-x1*)+A12·(x2),x2′=A2(x)x2,
on the positively invariant set Ωc⊂ℝ+n1+n2. Here x1=(Sn,Sd) and x2=(Ic). Here x1∈ℝ+2 denotes (its components) the number of uninfected individuals and x2∈ℝ+1 denotes the number of infected components, x1*=ℰc0. We have to prove that the following conditions are satisfied. (H1)
The system is defined on a positively invariant setΩc of the nonnegative orthant. The system is dissipative on Ωc.(H2)
The subsystemx1′=A1·(x1,0)·(x1-x1*) is globally asymptotically stable at the equilibrium x1* on the canonical projection of Ωc on ℝ+n1.(H3)
The matrixA2(x) is Metzler (A Metzler matrix is a matrix with off-diagonal nonnegative entries [28]) and irreducible for any given x∈Ωc.(H4)
There exists an upper-bound matrixA-2 for 𝕄={A2(x)/x∈Ωc} with the property that either A2∉𝕄 or if A2∉𝕄, (i.e., A2=maxΩc𝕄), then for any x-∈Ωc such that A-2=A2(x-), x-∈ℝ+n1×{0} (i.e., the points where the maximum is realized are contained in the disease-free submanifold).(H5)
ρ(A-2)≤0.If conditions (H1–H5) are satisfied, then ℰ0 is globally asymptotically stable in Ωc.We express the subsystemx1′=A1·(x1,0)·(x1-x1*) as follows:
(9)Sn′=Λπ0-(α+μ+ω)Sn,Sd′=Λπ1+αSn-(μ+ω)Sd.
System (9) is a linear system which is globally asymptotically stable at the equilibrium (Λπ0/(α+μ+ω),Λ[απ0+π1(μ+α+ω)]/(μ+ω)(μ+α+ω)), which corresponds to ℰc0, satisfying conditions (H1) and (H2). Hence, matrix A2(x) is given by
(10)A2=[-(μ+ω)(1-βSdN(μ+ω))],
Theorem 3 suggests that βSd/N<(μ+ω), whenever ℛc<1, hence matrix A2(x) is a Metzler matrix for any x∈Ωc and this satisfies conditions (H3) and (H4). The upper bound of A2(x) is given by
(11)A2=[-(μ+ω)(1-β[απ0+(α+μ+ω)π1](π0+π1)(μ+ω)(α+μ+ω))].
Condition (H5) requires that ρ(A2(x))≤0, that is,
(12)β[απ0+(α+μ+ω)π1](π0+π1)(μ+ω)(α+μ+ω)≤1,
which is the reproductive number for system (4). Thus, ℰc0 is globally asymptotically stable whenever ℛc≤1. We summarize the result in Theorem 4.Theorem 4.
The disease-free equilibrium(ℰc0) of model system (1) is globally asymptotically stable (GAS) if ℛc≤1 and unstable if ℛc>1.


### 3.2. Permanence of the Model

We first present the following definitions that are similar to those in [29, 30].Definition 5.
Model system (4) is said to be uniformly persistent if there is an ϵ>0 (independent of the initial data) such that every solution with positive initial conditions satisfies
(13)ϵ≤liminft→∞Sn(t),ϵ≤liminft→∞Sd(t),ϵ≤liminft→∞Ic(t).Definition 6.
Model system (4) is said to be permanent if there exists a compact region M⊂Ω∘c (the interior of Ωc) such that every solution of system (4) with positive initial conditions will eventually enter and remain in the region M.Since a test for permanence amounts to testing a suitable Lyapunov function, a suitable candidate for the model system (4) is
(14)ℒc(Sn,Sd,Ic)=Snk1Sdk2Ick3,
for k1,k2,k3>0 [31]. This function takes value zero on ∂Ωc (the boundary of Ωc) and is strictly positive in the interior of Ωc (Ω∘c). We show that model system (4) is uniformly permanent. The function defined in (14) vanishes for any values of (Sn,Sd,Ih) on the boundary and is strictly positive in Ω∘c, but we need to show that its derivative with respect to time for all points in Ω∘c close to the boundary is also positive. Differentiating (14) and writing ℒc for short, we obtain
(15)ℒc′ℒc=k1Sn′Sn+k2Sd′Sd+k3Ic′Ic∶=ℏ(Sn,Sd,Ic),(16)∫0T‍ℏ(Sn,Sd,Ic)dt=lnSnk1Sdk2Ick3≥0,
for some choice of parameters such that
(17)Snk1Sdk2Ick3>1.
Condition (17) is strong because there might be some parts of ∂Ω on which ℒc′=0. Therefore, we consider the time average of the derivative a condition which ensures that the boundary is a uniform repeller to orbits not starting on ∂Ω. For T large, the time average of the derivative is given by
(18)0<k1TlnSn(T)Sn(0)+k2TlnSd(T)Sd(0)+k3TlnIc(T)Ic(0).
It is actually this weaker condition that makes ℒc an average Lyapunov function [31]. During epidemics, Sn(T)>Sn(0),Sd(T)>Sd(0), consequently, the above inequality holds. Therefore, for orbits not starting on the boundary move away from it. Thus, model system (4) is permanent. Even though what orbits do away from the vicinity of the boundary of Ωc is immaterial, we need to show in our case that they tend to the endemic fixed point of the system. This is accomplished by showing that ℒc′ is negative semidefinite. From (15), on substituting the values of the derivatives therein, we obtain
(19)ℒc′ℒc=k1(μ+α+ω)(π0Λ(μ+α+ω)Sn-1)+k2(λc+μ+ω)(Λ[απ0+(α+μ+ω)π1](α+μ+ω)(λc+μ+ω)Sd-1)+k3(μ+ω)(Λπ2+λcSd(μ+ω)Ic-1).
The right- hand side of the above expression is zero only at ℰc*, the largest positive compact invariant set Ωc. Also,
(20)ℒc′≥-[k1(μ+α+ω)+k2(λc+μ+ω)+k3(μ+ω)]ℒc≤0,
and since
(21)k1(μ+α+ω)+k2(λc+μ+ω)+k3(μ+ω):=K>0,
then
(22)ℒc′ℒc≥-K<0,
and consequently, ℒc′<0. Hence, by LaSalle invariance principle [32] any trajectory starting in Ωc moves towards the maximal invariant set ℰc*, which is asymptotically stable in Ω∘c.


### 3.3. Endemic Equilibrium and Its Stability

From model system (4), we note that N=Λ/(μ+ω), and Sn=Λπ0/(α+μ+ω); hence system (4) can be rewritten as
(23)Sd′=Λ[απ0+(α+μ+ω)π1](α+μ+ω)-βIcSdΛ/(μ+ω)-(μ+ω)Sd,Ic′=Λπ2+βIcSdΛ/(μ+ω)-(μ+ω)Ic.
We now study the global stability of system (23) using the Poincare-Bendixson Theorem [33]. Denote the right-hand side of (23) by f and g for Sd and Ic, respectively, and choose a Dulac function as D(Sd,Ic)=1/SdIc. Then we have
(24)∂(Df)∂Sd+∂(Dg)∂Ic=-ΛSdIc(Λ[απ0+(α+μ+ω)π1](α+μ+ω)Sd+π2Ic)<0.
Thus, by Dulac’s criterion, there are no periodic orbits in Ωc. Since Ωc is positively invariant, and the endemic equilibrium exists whenever ℛc>1, then, it follows from the Poincare-Bendixson Theorem [33] that all solutions of the limiting system originating in Ωc remain in Ωc, for all t≥0. Further, the absence of periodic orbits in Ωc implies that the endemic equilibrium ℰc* of the HCV-only model is globally asymptotically stable whenever ℛc>1. We summarize the result in Theorem 7.Theorem 7.
The endemic equilibrium of the HCV-only model system (4) is GAS in Ωc whenever ℛc>1.


### 3.4. Sensitivity Analysis ofℛc

Sensitivity analysis of model parameters is very important to design and control strategies as well as a direction to future research. There are many methods [34] available for conducting sensitivity analysis such as differential analysis, response surface methodology, the Fourier amplitude sensitivity test (FAST) and other variance decomposition procedures, fast probability integration, and sampling-based procedures. In this section, sensitivity analysis based on Latin hypercube sampling (LHS) [34–39] has been performed using relevant parameters’ values in Table 1. Sensitivity analysis assesses the amount and type of change inherent in the model as captured by the terms that define the reproductive number. If the reproductive number is very sensitive to a particular parameter, then a perturbation of the conditions that connect the dynamics to such a parameter may prove useful in identifying policies or intervention strategies that reduce epidemic prevalence. In this section the Partial rank correlation coefficients (PRCCs) were calculated to estimate the correlation between values of ℛc and the model parameters across 1000 random draws from the empirical distribution of ℛc and its associated parameters.Table 1
Model parameters and their interpretations.
Parameter description
Symbol
Point estimate
Range
Source
Proportion of prisoners recruited intoSn
π
0
0.4
0–0.4
[24]
Proportion of prisoners recruited intoSd
π
1
0.3
0–0.3
[24]
Proportion of prisoners recruited intoIc
π
2
0.1
0–0.1
Assumed
Proportion of prisoners recruited intoIh
π
3
0.1
0–0.1
Assumed
Proportion of prisoners recruited intoIch
π
4
0.1
0–0.1
Assumed
Rate of progression to AIDS stage
δ
0.125 yr−1
0.01–0.125
[26]
Recruitment rate for prisoners
Λ
123 per 100 000
—
[24, 40]
Modification parameter
σ
1.2
≥1
[24]
Modification parameter
ϕ
1.2
≥1
[24]
Modification parameter
η
1.2
≥1
[24]
Natural mortality rate
μ
0.0142 yr−1
0.01–0.02
[24]
HCV transmissibility
β
0.1
0.0084–0.1
[41, 42]
HIV transmissibility
θ
0.05
0.0084–0.1
[41, 42]
Behaviour change
α
0.1 yr−1
0–0.2
[24]
Release rate
ω
0.2  yr−1
0.125–0.33
[24, 43]Figure2 illustrates the PRCCs using ℛc as an output variable. Results here suggest that HCV transmission is most sensitive to influence the magnitude of ℛc than any other parameter. An increase in the magnitude of HCV transmission will result in an increase in the magnitude of the reproductive number, while an increase in prison release rate will lead to a decrease in ℛc. We now examine the dependence of the six model parameters, namely, HCV transmissibility, prison release rate, incoming susceptible IDUs, incoming susceptible non-IDUs, IDU adoption rate, and natural mortality rate using the Latin hypercube sampling technique.Figure 2
Partial rank correlation coefficients showing the effects of parameter variation onℛc using ranges in the table. Parameters with positive PRCCs will increase ℛc when they are increased, whereas parameters with negative PRCCs will decrease ℛc when they are increased.A close analysis of results depicted in Figure3 shows that HCV transmission and prison release rate are highly correlated on ℛc, positively and negatively, respectively. Results displayed here are in agreement with numerical findings in Figure 2 and analytical results on (25).Scatter plots for the basic reproductive numberℛc and six parameter values (β, ω, π1, π0, α, and μ). These results were obtained from Latin hypercube sampling using a sample size of 1000.
(a)
(b)
(c)
(d)
(e)
(f)


### 3.5. Sensitivity Indices ofℛc Bases on Perturbation of Fixed Point Estimates

We now investigate the sensitivity ofℛc with respect to each of the parameters, following Arriola and Hyman [44]. The normalized forward sensitivity index with respect to each of the parameters is presented below:
(25)βℛc∂ℛc∂β=1,ωℛc∂ℛc∂ω=-ω(π1(μ+ω)2+(π0+π1)(α+2(μ+ω)))(μ+ω)(α+μ+ω)(π1(μ+ω)+(π0+π1)α)=-0.996725,π0ℛc∂ℛc∂π0=-π0π1(μ+ω)(π0+π1)(π1(μ+ω)+(π0+π1)α)=-0.169235,π1ℛc∂ℛc∂π1=π0π1(μ+ω)(π0+π1)(π1(μ+ω)+(π0+π1)α)=0.169235,αℛc∂ℛc∂α=απ0(μ+ω)(α+μ+ω)(π1(μ+ω)+(π0+π1)α)=0.131372,μℛc∂ℛc∂μ=-μ(π1(μ+ω)2+(π0+π1)(α+2(μ+ω)))(μ+ω)(α+μ+ω)(π1(μ+ω)+(π0+π1)α)=-0.141535.
Results in (25) suggest that HCV transmission rate β and prison release rate ω are highly correlated with the reproductive number ℛc. An increase in β will bring about an increase of the same proportion in ℛc and a decrease in β will result in a decrease in ℛc with about an equivalent magnitude. Thus, a 20% increase in β will result in 20% increase in ℛc. An increase in ω will lead to a decrease in ℛc not by a similar magnitude, but close to it. Although π0 (the proportion of incoming susceptible nonintravenous drug users) and π1 (the proportion of incoming susceptible intravenous drug users) have equal influence on the magnitude of ℛc, π0 will lead to a decrease in ℛc and the reverse occurs for π1. Results obtained here (on (25)), are in total agreement with the earlier findings in Figures 2 and 3.


## 3.1. Disease-Free Equilibrium and Its Stability

System (4) has a disease-free equilibrium (DFE) given by
(6)ℰc0=(Sn0,Sd0,Ic0)=(Λπ0(α+μ+ω),Λ[απ0+π1(μ+α+ω)](μ+ω)(μ+α+ω),0).
Now, we compute the basic reproductive number, ℛ0, for system (4). Following Van Den Driessche and Watmough [17], the reproductive number for system (4) is given by
(7)ℛc=β[απ0+(α+μ+ω)π1](π0+π1)(μ+ω)(α+μ+ω),whereℛc measures the average number of new secondary cases generated by a single HCV infective individual during his/her entire infectious period when he/she is introduced into a susceptible population within a correctional setting in the absence of HCV intervention strategies. Using [17, Theorem  2], the following result is established.Theorem 3.
ℰ
c
0 is locally asymptotically stable (LAS) if ℛc<1, and unstable if ℛc>1.In order to examine the global stability ofℰc0 we use the method proposed by Kamgang and Sallet (2008). By closely following Kamgang and Sallet (2008) [25–27], we write system (1) in the form
(8)x1′=A1(x)·(x1-x1*)+A12·(x2),x2′=A2(x)x2,
on the positively invariant set Ωc⊂ℝ+n1+n2. Here x1=(Sn,Sd) and x2=(Ic). Here x1∈ℝ+2 denotes (its components) the number of uninfected individuals and x2∈ℝ+1 denotes the number of infected components, x1*=ℰc0. We have to prove that the following conditions are satisfied. (H1)
The system is defined on a positively invariant setΩc of the nonnegative orthant. The system is dissipative on Ωc.(H2)
The subsystemx1′=A1·(x1,0)·(x1-x1*) is globally asymptotically stable at the equilibrium x1* on the canonical projection of Ωc on ℝ+n1.(H3)
The matrixA2(x) is Metzler (A Metzler matrix is a matrix with off-diagonal nonnegative entries [28]) and irreducible for any given x∈Ωc.(H4)
There exists an upper-bound matrixA-2 for 𝕄={A2(x)/x∈Ωc} with the property that either A2∉𝕄 or if A2∉𝕄, (i.e., A2=maxΩc𝕄), then for any x-∈Ωc such that A-2=A2(x-), x-∈ℝ+n1×{0} (i.e., the points where the maximum is realized are contained in the disease-free submanifold).(H5)
ρ(A-2)≤0.If conditions (H1–H5) are satisfied, then ℰ0 is globally asymptotically stable in Ωc.We express the subsystemx1′=A1·(x1,0)·(x1-x1*) as follows:
(9)Sn′=Λπ0-(α+μ+ω)Sn,Sd′=Λπ1+αSn-(μ+ω)Sd.
System (9) is a linear system which is globally asymptotically stable at the equilibrium (Λπ0/(α+μ+ω),Λ[απ0+π1(μ+α+ω)]/(μ+ω)(μ+α+ω)), which corresponds to ℰc0, satisfying conditions (H1) and (H2). Hence, matrix A2(x) is given by
(10)A2=[-(μ+ω)(1-βSdN(μ+ω))],
Theorem 3 suggests that βSd/N<(μ+ω), whenever ℛc<1, hence matrix A2(x) is a Metzler matrix for any x∈Ωc and this satisfies conditions (H3) and (H4). The upper bound of A2(x) is given by
(11)A2=[-(μ+ω)(1-β[απ0+(α+μ+ω)π1](π0+π1)(μ+ω)(α+μ+ω))].
Condition (H5) requires that ρ(A2(x))≤0, that is,
(12)β[απ0+(α+μ+ω)π1](π0+π1)(μ+ω)(α+μ+ω)≤1,
which is the reproductive number for system (4). Thus, ℰc0 is globally asymptotically stable whenever ℛc≤1. We summarize the result in Theorem 4.Theorem 4.
The disease-free equilibrium(ℰc0) of model system (1) is globally asymptotically stable (GAS) if ℛc≤1 and unstable if ℛc>1.


## 3.2. Permanence of the Model

We first present the following definitions that are similar to those in [29, 30].Definition 5.
Model system (4) is said to be uniformly persistent if there is an ϵ>0 (independent of the initial data) such that every solution with positive initial conditions satisfies
(13)ϵ≤liminft→∞Sn(t),ϵ≤liminft→∞Sd(t),ϵ≤liminft→∞Ic(t).Definition 6.
Model system (4) is said to be permanent if there exists a compact region M⊂Ω∘c (the interior of Ωc) such that every solution of system (4) with positive initial conditions will eventually enter and remain in the region M.Since a test for permanence amounts to testing a suitable Lyapunov function, a suitable candidate for the model system (4) is
(14)ℒc(Sn,Sd,Ic)=Snk1Sdk2Ick3,
for k1,k2,k3>0 [31]. This function takes value zero on ∂Ωc (the boundary of Ωc) and is strictly positive in the interior of Ωc (Ω∘c). We show that model system (4) is uniformly permanent. The function defined in (14) vanishes for any values of (Sn,Sd,Ih) on the boundary and is strictly positive in Ω∘c, but we need to show that its derivative with respect to time for all points in Ω∘c close to the boundary is also positive. Differentiating (14) and writing ℒc for short, we obtain
(15)ℒc′ℒc=k1Sn′Sn+k2Sd′Sd+k3Ic′Ic∶=ℏ(Sn,Sd,Ic),(16)∫0T‍ℏ(Sn,Sd,Ic)dt=lnSnk1Sdk2Ick3≥0,
for some choice of parameters such that
(17)Snk1Sdk2Ick3>1.
Condition (17) is strong because there might be some parts of ∂Ω on which ℒc′=0. Therefore, we consider the time average of the derivative a condition which ensures that the boundary is a uniform repeller to orbits not starting on ∂Ω. For T large, the time average of the derivative is given by
(18)0<k1TlnSn(T)Sn(0)+k2TlnSd(T)Sd(0)+k3TlnIc(T)Ic(0).
It is actually this weaker condition that makes ℒc an average Lyapunov function [31]. During epidemics, Sn(T)>Sn(0),Sd(T)>Sd(0), consequently, the above inequality holds. Therefore, for orbits not starting on the boundary move away from it. Thus, model system (4) is permanent. Even though what orbits do away from the vicinity of the boundary of Ωc is immaterial, we need to show in our case that they tend to the endemic fixed point of the system. This is accomplished by showing that ℒc′ is negative semidefinite. From (15), on substituting the values of the derivatives therein, we obtain
(19)ℒc′ℒc=k1(μ+α+ω)(π0Λ(μ+α+ω)Sn-1)+k2(λc+μ+ω)(Λ[απ0+(α+μ+ω)π1](α+μ+ω)(λc+μ+ω)Sd-1)+k3(μ+ω)(Λπ2+λcSd(μ+ω)Ic-1).
The right- hand side of the above expression is zero only at ℰc*, the largest positive compact invariant set Ωc. Also,
(20)ℒc′≥-[k1(μ+α+ω)+k2(λc+μ+ω)+k3(μ+ω)]ℒc≤0,
and since
(21)k1(μ+α+ω)+k2(λc+μ+ω)+k3(μ+ω):=K>0,
then
(22)ℒc′ℒc≥-K<0,
and consequently, ℒc′<0. Hence, by LaSalle invariance principle [32] any trajectory starting in Ωc moves towards the maximal invariant set ℰc*, which is asymptotically stable in Ω∘c.


## 3.3. Endemic Equilibrium and Its Stability

From model system (4), we note that N=Λ/(μ+ω), and Sn=Λπ0/(α+μ+ω); hence system (4) can be rewritten as
(23)Sd′=Λ[απ0+(α+μ+ω)π1](α+μ+ω)-βIcSdΛ/(μ+ω)-(μ+ω)Sd,Ic′=Λπ2+βIcSdΛ/(μ+ω)-(μ+ω)Ic.
We now study the global stability of system (23) using the Poincare-Bendixson Theorem [33]. Denote the right-hand side of (23) by f and g for Sd and Ic, respectively, and choose a Dulac function as D(Sd,Ic)=1/SdIc. Then we have
(24)∂(Df)∂Sd+∂(Dg)∂Ic=-ΛSdIc(Λ[απ0+(α+μ+ω)π1](α+μ+ω)Sd+π2Ic)<0.
Thus, by Dulac’s criterion, there are no periodic orbits in Ωc. Since Ωc is positively invariant, and the endemic equilibrium exists whenever ℛc>1, then, it follows from the Poincare-Bendixson Theorem [33] that all solutions of the limiting system originating in Ωc remain in Ωc, for all t≥0. Further, the absence of periodic orbits in Ωc implies that the endemic equilibrium ℰc* of the HCV-only model is globally asymptotically stable whenever ℛc>1. We summarize the result in Theorem 7.Theorem 7.
The endemic equilibrium of the HCV-only model system (4) is GAS in Ωc whenever ℛc>1.


## 3.4. Sensitivity Analysis ofℛc

Sensitivity analysis of model parameters is very important to design and control strategies as well as a direction to future research. There are many methods [34] available for conducting sensitivity analysis such as differential analysis, response surface methodology, the Fourier amplitude sensitivity test (FAST) and other variance decomposition procedures, fast probability integration, and sampling-based procedures. In this section, sensitivity analysis based on Latin hypercube sampling (LHS) [34–39] has been performed using relevant parameters’ values in Table 1. Sensitivity analysis assesses the amount and type of change inherent in the model as captured by the terms that define the reproductive number. If the reproductive number is very sensitive to a particular parameter, then a perturbation of the conditions that connect the dynamics to such a parameter may prove useful in identifying policies or intervention strategies that reduce epidemic prevalence. In this section the Partial rank correlation coefficients (PRCCs) were calculated to estimate the correlation between values of ℛc and the model parameters across 1000 random draws from the empirical distribution of ℛc and its associated parameters.Table 1
Model parameters and their interpretations.
Parameter description
Symbol
Point estimate
Range
Source
Proportion of prisoners recruited intoSn
π
0
0.4
0–0.4
[24]
Proportion of prisoners recruited intoSd
π
1
0.3
0–0.3
[24]
Proportion of prisoners recruited intoIc
π
2
0.1
0–0.1
Assumed
Proportion of prisoners recruited intoIh
π
3
0.1
0–0.1
Assumed
Proportion of prisoners recruited intoIch
π
4
0.1
0–0.1
Assumed
Rate of progression to AIDS stage
δ
0.125 yr−1
0.01–0.125
[26]
Recruitment rate for prisoners
Λ
123 per 100 000
—
[24, 40]
Modification parameter
σ
1.2
≥1
[24]
Modification parameter
ϕ
1.2
≥1
[24]
Modification parameter
η
1.2
≥1
[24]
Natural mortality rate
μ
0.0142 yr−1
0.01–0.02
[24]
HCV transmissibility
β
0.1
0.0084–0.1
[41, 42]
HIV transmissibility
θ
0.05
0.0084–0.1
[41, 42]
Behaviour change
α
0.1 yr−1
0–0.2
[24]
Release rate
ω
0.2  yr−1
0.125–0.33
[24, 43]Figure2 illustrates the PRCCs using ℛc as an output variable. Results here suggest that HCV transmission is most sensitive to influence the magnitude of ℛc than any other parameter. An increase in the magnitude of HCV transmission will result in an increase in the magnitude of the reproductive number, while an increase in prison release rate will lead to a decrease in ℛc. We now examine the dependence of the six model parameters, namely, HCV transmissibility, prison release rate, incoming susceptible IDUs, incoming susceptible non-IDUs, IDU adoption rate, and natural mortality rate using the Latin hypercube sampling technique.Figure 2
Partial rank correlation coefficients showing the effects of parameter variation onℛc using ranges in the table. Parameters with positive PRCCs will increase ℛc when they are increased, whereas parameters with negative PRCCs will decrease ℛc when they are increased.A close analysis of results depicted in Figure3 shows that HCV transmission and prison release rate are highly correlated on ℛc, positively and negatively, respectively. Results displayed here are in agreement with numerical findings in Figure 2 and analytical results on (25).Scatter plots for the basic reproductive numberℛc and six parameter values (β, ω, π1, π0, α, and μ). These results were obtained from Latin hypercube sampling using a sample size of 1000.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.5. Sensitivity Indices ofℛc Bases on Perturbation of Fixed Point Estimates

We now investigate the sensitivity ofℛc with respect to each of the parameters, following Arriola and Hyman [44]. The normalized forward sensitivity index with respect to each of the parameters is presented below:
(25)βℛc∂ℛc∂β=1,ωℛc∂ℛc∂ω=-ω(π1(μ+ω)2+(π0+π1)(α+2(μ+ω)))(μ+ω)(α+μ+ω)(π1(μ+ω)+(π0+π1)α)=-0.996725,π0ℛc∂ℛc∂π0=-π0π1(μ+ω)(π0+π1)(π1(μ+ω)+(π0+π1)α)=-0.169235,π1ℛc∂ℛc∂π1=π0π1(μ+ω)(π0+π1)(π1(μ+ω)+(π0+π1)α)=0.169235,αℛc∂ℛc∂α=απ0(μ+ω)(α+μ+ω)(π1(μ+ω)+(π0+π1)α)=0.131372,μℛc∂ℛc∂μ=-μ(π1(μ+ω)2+(π0+π1)(α+2(μ+ω)))(μ+ω)(α+μ+ω)(π1(μ+ω)+(π0+π1)α)=-0.141535.
Results in (25) suggest that HCV transmission rate β and prison release rate ω are highly correlated with the reproductive number ℛc. An increase in β will bring about an increase of the same proportion in ℛc and a decrease in β will result in a decrease in ℛc with about an equivalent magnitude. Thus, a 20% increase in β will result in 20% increase in ℛc. An increase in ω will lead to a decrease in ℛc not by a similar magnitude, but close to it. Although π0 (the proportion of incoming susceptible nonintravenous drug users) and π1 (the proportion of incoming susceptible intravenous drug users) have equal influence on the magnitude of ℛc, π0 will lead to a decrease in ℛc and the reverse occurs for π1. Results obtained here (on (25)), are in total agreement with the earlier findings in Figures 2 and 3.


## 4. HIV-Only Submodel

Consider the HIV-only submodel (obtained by settingIc=Ich=0) in system (1), so that we have
(26)Sn′=Λπ0-(α+μ+ω)Sn,Sd′=Λπ1+αSn-(λc+μ+ω)Sd,Ih′=Λπ3+λhSd-(μ+ω+δ)Ih.
For system (26), the first octant in the state space is positively invariant and attracting; that is, solutions that start where all the variables are nonnegative remain there. Thus, system (26) will be analyzed in a suitable region
(27)Ωh={(Sn,Sd,Ih)∈ℝ+3:N≤Λ(μ+ω),
which is positively invariant and attracting. Existence, uniqueness, and continuation results for system (4) hold in this region.


### 4.1. Disease-Free Equilibrium and Its Stability

System (26) has a disease-free equilibrium (DFE) given by
(28)ℰh0=(Sn0,Sd0,Ih0)=(Λπ0(α+μ+ω),Λ[απ0+π1(μ+α+ω)](μ+ω)(μ+α+ω),0).
Following Van Den Driessche and Watmough [17], the reproductive number for system (26) is given by
(29)ℛh=θ[απ0+(α+μ+ω)π1](π0+π1)(μ+ω+δ)(α+μ+ω),whereℛh, measures the average number of new secondary cases generated by a single HIV infective individual during his/her entire infectious period when he/she is introduced into a susceptible population within a correctional setting in the absence of HIV intervention strategies. Using [17, Theorem  2], the following result is established.Theorem 8.
The DFEℰh0 of system (26) is globally asymptotically stable provided ℛh≤1 and unstable if ℛh>1.Proof.
By closely following the approach on Section4.1 we set
(30)x1=(Sn,Sd),x2=(Ih),
so that the subsystem x˙1=A1(x1,0)·(x1-x1*),
(31)Sn′=Λπ0-(α+μ+ω)Sn,Sd′=Λπ1+αSn-(μ+ω)Sd,
and (31) is a linear system which is globally asymptotically stable at the equilibrium (Λπ0/(α+μ+ω),Λ[απ0+π1(μ+α+ω)]/(μ+ω)(μ+α+ω)) corresponding to the disease-free equilibrium ℰh0 which satisfies conditions (H1) and (H2). The matrix A2(x) is given by
(32)A2(x)=[-(μ+δ+ω)(1-θSdN(μ+δ+ω))].
The upper bound for x∈Ω is given by
(33)A2=[(1-θ[απ0+(α+μ+ω)π1](π0+π1)(μ+ω+δ)(α+μ+ω))-(μ+δ+ω)×(1-θ[απ0+(α+μ+ω)π1](π0+π1)(μ+ω+δ)(α+μ+ω))].
Thus, conditions (H3) and (H4) are satisfied. Then condition (H5) is equivalent to ℛh≤1,
(34)θ[απ0+(α+μ+ω)π1](π0+π1)(μ+ω+δ)(α+μ+ω)≤1.
Thus, ℰh0 is globally asymptotically stable whenever ℛh≤1.


### 4.2. Endemic Equilibrium and Its Stability

Following the approach on (23), system (26) can be rewritten as
(35)Sd′=Λ[απ0+(α+μ+ω)π1](α+μ+ω)-θIhSdN-(μ+ω)Sd,Ih′=Λπ2+θIhSdN-(μ+ω+δ)Ih.
In order to investigate the global stability of the endemic equilibrium, we adopt the approach by Korobeinikov (2006) [45] (and closely follow the approach in [46]). Assume that there exists ℰh* for all Sd, Ih>ϵ, for some ϵ>0. Let θIhSd/N=g(Sd,Ih) be a positive and monotonic function, and define the following continuous function in ℝ+2 (for more details, see Korobeinikov, 2006 [45]). A function
(36)V(Sd,Ih)=Sd-∫ϵSd‍g(Sd*,Ih*)g(τ,Ih*)dτ+Ih-∫ϵIh‍g(Sd*,Ih*)g(Sd*,τ)dτ.
If g(Sd,Ih) is monotonic with respect to its variables, then the endemic state ℰh* is the only extremum and the global minimum of this function. Indeed
(37)∂V∂Sd=1-g(Sd*,Ih*)g(Sd,Ih*),∂V∂Ih=1-g(Sd*,Ih*)g(Sd*,Ih),
grow monotonically, then the function g(Sd,Ih) has only one stationary point. Further, since
(38)∂2V∂Sd2=g(Sd*,Ih*)[g(Sd,Ih*)]2·∂g(Sd,Ih*)∂Sd,∂2V∂Ih2=g(Sd*,Ih*)[g(Sd*,Ih)]2·∂g(Sd*,Ih)∂Ih
are nonnegative, then the point ℰh* is a minimum. That is, V(Sd,Ih)≥V(Sd*,Ih*) and, hence, V is a Lyapunov function, and its time derivative is given by
(39)dVdt=Sd′-Sd′(g(Sd*,Ih*)g(Sh,Ih*))+Ih′-Ih′(g(Sd*,Ih*)g(Sd*,Ih)).
Recall that
(40)Λ[απ0+(α+μ+ω)π1](α+μ+ω)=g(Sd*,Ih*)+(μ+ω)Sd*,Λπ2=(μ+ω)Ih*-g(Sd*,Ih*).
We have
(41)dVdt=(μ+ω)Sd*(1-SdSd*)(1-g(Sd*,Ih*)g(Sd,Ih*))+(μ+ω+δ)Ih*(1-IhIh*)(1-g(Sd*,Ih*)g(Sd*,Ih))+g(Sd*,Ih*)(1-g(Sd*,Ih*)g(Sd,Ih*))(1-g(Sd*,Ih*)g(Sd,Ih*))+g(Sd,Ih)(1-g(Sd*,Ih*)g(Sd*,Ih))(1-g(Sd*,Ih*)g(Sd,Ih)).
Since ℰh*>0, the function g(Sd,Ih) is concave with respect to Ih, and ∂2g(Sd,Ih)/∂Ih2≤0, then dV/dt≤0 for all Sd, Ih>0. Also, the monotonicity of g(Sd,Ih) with respect to Sd and Ih ensures that
(42)(1-SdSd*)(1-g(Sd*,Ih*)g(Sd,Ih*))≤0,(1-IhIh*)(1-g(Sd*,Ih*)g(Sd*,Ih))≤0,(1-g(Sd*,Ih*)g(Sd,Ih*))(1-g(Sd*,Ih*)g(Sd,Ih*))≤0,(1-g(Sd*,Ih*)g(Sd*,Ih))(1-g(Sd*,Ih*)g(Sd,Ih))≤0
holds for all Sd, Ih>0. Thus, we establish the following result.Theorem 9.
The endemic equilibriumℰh* is globally asymptotically stable whenever conditions outlined in (42) are satisfied.


### 4.3. Sensitivity Analysis ofℛc

We now perform the sensitivity analysis ofℛh and the seven model parameters, namely, HIV transmission, prison release rate, incoming susceptible IDUs, incoming susceptible non-IDUs, IDU adoption rate, progression to AIDS stage, and natural mortality rate. Results in Figure 4 suggests that prison release rate is the more sensitive to ℛh than any of the aforementioned model parameters. It is also worth noting that an increase in prison release will lead to a decrease in the magnitude of ℛh. Although HIV transmission is not highly correlated to ℛh (compared to prison release rate), its contribution on increasing the magnitude of ℛh when it is increased cannot be neglected. Overall, the results depicted in Figures 2 and 4 suggest the need for HIV and HCV intervention strategies, in order to reduce the prevalence of the two diseases among prisoners.Figure 4
Partial rank correlation coefficients showing the effects of parameter variation onℛh using ranges in Table 1. Parameters with positive PRCCs will increase ℛh when they are increased, whereas parameters with negative PRCCs will decrease ℛh when they are increased.Results depicted in Figure5 support the earlier findings in Figure 4, that prison release rate and HIV transmission have a significant impact on influencing the magnitude of ℛh negatively and positively, respectively.Scatter plots for the basic reproductive numberℛh and seven mode parameter values (θ, δω, π1, π0, α, and μ). These results were obtained from Latin hypercube sampling using a sample size of 1000.
(a)
(b)
(c)
(d)
(e)
(f)


### 4.4. Sensitivity Indices ofℛh Bases on Perturbation of Fixed Point Estimates

Following Arriola and Hyman [44] we present the normalized forward sensitivity indices of ℛh with respect to θ, δμ, and ω(note that sensitivity indices for ℛh in relation to the model parameters π0,  π1, and α yield the same results as those on (25))
(43)θℛh∂ℛh∂θ=1,δℛh∂ℛh∂δ=-δμ+δ+ω=-0.522575,μℛh∂ℛh∂μ=-μ[(μ+ω)(π1+2α(π0+π1))+α(π0δ+α(π0+π1))](μ+α+ω)(δ+μ+ω)(π1(μ+ω)+α(π0+π1))=-0.0765568,ωℛh∂ℛh∂ω=-ω[(μ+ω)(π1+2α(π0+π1))+α(π0δ+α(π0+π1))](μ+α+ω)(δ+μ+ω)(π1(μ+ω)+α(π0+π1))=-0.53912.Analytical results on (43) support the numerical findings in Figures 4 and 5.


## 4.1. Disease-Free Equilibrium and Its Stability

System (26) has a disease-free equilibrium (DFE) given by
(28)ℰh0=(Sn0,Sd0,Ih0)=(Λπ0(α+μ+ω),Λ[απ0+π1(μ+α+ω)](μ+ω)(μ+α+ω),0).
Following Van Den Driessche and Watmough [17], the reproductive number for system (26) is given by
(29)ℛh=θ[απ0+(α+μ+ω)π1](π0+π1)(μ+ω+δ)(α+μ+ω),whereℛh, measures the average number of new secondary cases generated by a single HIV infective individual during his/her entire infectious period when he/she is introduced into a susceptible population within a correctional setting in the absence of HIV intervention strategies. Using [17, Theorem  2], the following result is established.Theorem 8.
The DFEℰh0 of system (26) is globally asymptotically stable provided ℛh≤1 and unstable if ℛh>1.Proof.
By closely following the approach on Section4.1 we set
(30)x1=(Sn,Sd),x2=(Ih),
so that the subsystem x˙1=A1(x1,0)·(x1-x1*),
(31)Sn′=Λπ0-(α+μ+ω)Sn,Sd′=Λπ1+αSn-(μ+ω)Sd,
and (31) is a linear system which is globally asymptotically stable at the equilibrium (Λπ0/(α+μ+ω),Λ[απ0+π1(μ+α+ω)]/(μ+ω)(μ+α+ω)) corresponding to the disease-free equilibrium ℰh0 which satisfies conditions (H1) and (H2). The matrix A2(x) is given by
(32)A2(x)=[-(μ+δ+ω)(1-θSdN(μ+δ+ω))].
The upper bound for x∈Ω is given by
(33)A2=[(1-θ[απ0+(α+μ+ω)π1](π0+π1)(μ+ω+δ)(α+μ+ω))-(μ+δ+ω)×(1-θ[απ0+(α+μ+ω)π1](π0+π1)(μ+ω+δ)(α+μ+ω))].
Thus, conditions (H3) and (H4) are satisfied. Then condition (H5) is equivalent to ℛh≤1,
(34)θ[απ0+(α+μ+ω)π1](π0+π1)(μ+ω+δ)(α+μ+ω)≤1.
Thus, ℰh0 is globally asymptotically stable whenever ℛh≤1.


## 4.2. Endemic Equilibrium and Its Stability

Following the approach on (23), system (26) can be rewritten as
(35)Sd′=Λ[απ0+(α+μ+ω)π1](α+μ+ω)-θIhSdN-(μ+ω)Sd,Ih′=Λπ2+θIhSdN-(μ+ω+δ)Ih.
In order to investigate the global stability of the endemic equilibrium, we adopt the approach by Korobeinikov (2006) [45] (and closely follow the approach in [46]). Assume that there exists ℰh* for all Sd, Ih>ϵ, for some ϵ>0. Let θIhSd/N=g(Sd,Ih) be a positive and monotonic function, and define the following continuous function in ℝ+2 (for more details, see Korobeinikov, 2006 [45]). A function
(36)V(Sd,Ih)=Sd-∫ϵSd‍g(Sd*,Ih*)g(τ,Ih*)dτ+Ih-∫ϵIh‍g(Sd*,Ih*)g(Sd*,τ)dτ.
If g(Sd,Ih) is monotonic with respect to its variables, then the endemic state ℰh* is the only extremum and the global minimum of this function. Indeed
(37)∂V∂Sd=1-g(Sd*,Ih*)g(Sd,Ih*),∂V∂Ih=1-g(Sd*,Ih*)g(Sd*,Ih),
grow monotonically, then the function g(Sd,Ih) has only one stationary point. Further, since
(38)∂2V∂Sd2=g(Sd*,Ih*)[g(Sd,Ih*)]2·∂g(Sd,Ih*)∂Sd,∂2V∂Ih2=g(Sd*,Ih*)[g(Sd*,Ih)]2·∂g(Sd*,Ih)∂Ih
are nonnegative, then the point ℰh* is a minimum. That is, V(Sd,Ih)≥V(Sd*,Ih*) and, hence, V is a Lyapunov function, and its time derivative is given by
(39)dVdt=Sd′-Sd′(g(Sd*,Ih*)g(Sh,Ih*))+Ih′-Ih′(g(Sd*,Ih*)g(Sd*,Ih)).
Recall that
(40)Λ[απ0+(α+μ+ω)π1](α+μ+ω)=g(Sd*,Ih*)+(μ+ω)Sd*,Λπ2=(μ+ω)Ih*-g(Sd*,Ih*).
We have
(41)dVdt=(μ+ω)Sd*(1-SdSd*)(1-g(Sd*,Ih*)g(Sd,Ih*))+(μ+ω+δ)Ih*(1-IhIh*)(1-g(Sd*,Ih*)g(Sd*,Ih))+g(Sd*,Ih*)(1-g(Sd*,Ih*)g(Sd,Ih*))(1-g(Sd*,Ih*)g(Sd,Ih*))+g(Sd,Ih)(1-g(Sd*,Ih*)g(Sd*,Ih))(1-g(Sd*,Ih*)g(Sd,Ih)).
Since ℰh*>0, the function g(Sd,Ih) is concave with respect to Ih, and ∂2g(Sd,Ih)/∂Ih2≤0, then dV/dt≤0 for all Sd, Ih>0. Also, the monotonicity of g(Sd,Ih) with respect to Sd and Ih ensures that
(42)(1-SdSd*)(1-g(Sd*,Ih*)g(Sd,Ih*))≤0,(1-IhIh*)(1-g(Sd*,Ih*)g(Sd*,Ih))≤0,(1-g(Sd*,Ih*)g(Sd,Ih*))(1-g(Sd*,Ih*)g(Sd,Ih*))≤0,(1-g(Sd*,Ih*)g(Sd*,Ih))(1-g(Sd*,Ih*)g(Sd,Ih))≤0
holds for all Sd, Ih>0. Thus, we establish the following result.Theorem 9.
The endemic equilibriumℰh* is globally asymptotically stable whenever conditions outlined in (42) are satisfied.


## 4.3. Sensitivity Analysis ofℛc

We now perform the sensitivity analysis ofℛh and the seven model parameters, namely, HIV transmission, prison release rate, incoming susceptible IDUs, incoming susceptible non-IDUs, IDU adoption rate, progression to AIDS stage, and natural mortality rate. Results in Figure 4 suggests that prison release rate is the more sensitive to ℛh than any of the aforementioned model parameters. It is also worth noting that an increase in prison release will lead to a decrease in the magnitude of ℛh. Although HIV transmission is not highly correlated to ℛh (compared to prison release rate), its contribution on increasing the magnitude of ℛh when it is increased cannot be neglected. Overall, the results depicted in Figures 2 and 4 suggest the need for HIV and HCV intervention strategies, in order to reduce the prevalence of the two diseases among prisoners.Figure 4
Partial rank correlation coefficients showing the effects of parameter variation onℛh using ranges in Table 1. Parameters with positive PRCCs will increase ℛh when they are increased, whereas parameters with negative PRCCs will decrease ℛh when they are increased.Results depicted in Figure5 support the earlier findings in Figure 4, that prison release rate and HIV transmission have a significant impact on influencing the magnitude of ℛh negatively and positively, respectively.Scatter plots for the basic reproductive numberℛh and seven mode parameter values (θ, δω, π1, π0, α, and μ). These results were obtained from Latin hypercube sampling using a sample size of 1000.
(a)
(b)
(c)
(d)
(e)
(f)


## 4.4. Sensitivity Indices ofℛh Bases on Perturbation of Fixed Point Estimates

Following Arriola and Hyman [44] we present the normalized forward sensitivity indices of ℛh with respect to θ, δμ, and ω(note that sensitivity indices for ℛh in relation to the model parameters π0,  π1, and α yield the same results as those on (25))
(43)θℛh∂ℛh∂θ=1,δℛh∂ℛh∂δ=-δμ+δ+ω=-0.522575,μℛh∂ℛh∂μ=-μ[(μ+ω)(π1+2α(π0+π1))+α(π0δ+α(π0+π1))](μ+α+ω)(δ+μ+ω)(π1(μ+ω)+α(π0+π1))=-0.0765568,ωℛh∂ℛh∂ω=-ω[(μ+ω)(π1+2α(π0+π1))+α(π0δ+α(π0+π1))](μ+α+ω)(δ+μ+ω)(π1(μ+ω)+α(π0+π1))=-0.53912.Analytical results on (43) support the numerical findings in Figures 4 and 5.


## 5. Analysis of the Full Model

Having analysed the dynamics of the two submodels, the full HCV-HIV model is now considered (1). It is evident that system (1) has a disease-free equilibrium (DFE) given by
(44)ℰch0=(Sn0,Sd0,Ic0,Ih0,Ich0)=(Λπ0(α+μ+ω),Λ[απ0+π1(μ+α+ω)](μ+ω)(μ+α+ω),0,0,0).
Using the next generation method, the reproductive number for system (1) is given by
(45)ℛch=max(β[απ0+(α+μ+ω)π1](π0+π1)(μ+ω)(α+μ+ω),θ[απ0+(α+μ+ω)π1](π0+π1)(μ+ω+δ)(α+μ+ω))=max(ℛc,ℛh).
By closely following Kamgang and Sallet (2008) [25], and the notations of Section 4.1 we set
(46)x1=(Sn,Sn),x2=(Ic,Ih,Ich).
We express the subsystem x˙1=A1(x1,0)·(x1-x1*),
(47)Sn′=Λπ0-(α+μ+ω)Sn,Sd′=Λπ1+αSn-(μ+ω)Sd,
and (47) is a linear system which is globally asymptotically stable at the equilibrium (Λπ0/(α+μ+ω),Λ[απ0+π1(μ+α+ω)]/(μ+ω)(μ+α+ω)) corresponding to the disease-free equilibrium ℰch0 which satisfies conditions (H1) and (H2). The matrix A2(x) is given by(48)A2(x)=[-((μ+ω)-βSdN)0ηβSdN0-((μ+δ+ω)-θSdN)ηθSdN00-(μ+ϕδ+ω)].The upper bound for x∈Ω is given by(49)A2=[-((μ+ω)-β[απ0+(α+μ+ω)π1](π0+π1)(α+μ+ω))0ηβ[απ0+(α+μ+ω)π1](π0+π1)(α+μ+ω)0-((μ+δ+ω)-θ[απ0+(α+μ+ω)π1](π0+π1)(α+μ+ω))ηθ[απ0+(α+μ+ω)π1](π0+π1)(α+μ+ω)00-(μ+ϕδ+ω)].Thus, conditions (H3) and (H4) are satisfied. Then condition (H5) is equivalent to ℛ0≤1. The maximum of two quantities in (50) is the reproductive number for system (1)
(50)β[απ0+(α+μ+ω)π1](π0+π1)(μ+ω)(α+μ+ω)≤1,θ[απ0+(α+μ+ω)π1](π0+π1)(μ+ω+δ)(α+μ+ω)≤1.
Thus, ℰ0 is globally asymptotically stable whenever ℛ0≤1. We summarize the result in Theorem 10.Theorem 10.
The disease-free equilibrium(ℰ0) of model system (1) is globally asymptotically stable (GAS) if ℛ0≤1 and unstable if ℛ0>1.Due to the complex nature of the endemic equilibrium for system (1), we shall use the following symmetric conditions, established from the results of Theorems 7 and 9:(a)
ℛc<1 (HCV dies out) and ℛh>1 (HIV persists),(b)
ℛh<1 (HIV persists) and ℛc>1 (HCV persists),(c)
ℛh>1 (HIV persists) and ℛc>1 (HCV persists). The symmetric conditions above may be summarized diagrammatically as shown in Figure 6.Figure 6
Symmetric conditions forℛ0>1, n∈ℵ,(n>1).With the aid of the above symmetric conditions and results deduced from Theorems7, 9, and 10, we deduce that for ℛ0>1 the endemic equilibrium may be globally asymptotically stable.


### 5.1. Impact of HCV on HIV Prevalence and Vice Versa

In many epidemiological models, the magnitude of the reproductive number is associated with the level of infection. The same is true for model (1). Expressions ℛc and ℛh can be rewritten as
(51)ℛc=β(π0+π1)(μ+ω)απ0+π1(α+μ+ω)(α+μ+ω)⇒απ0+π1(α+μ+ω)(α+μ+ω)=(π0+π1)(μ+ω)ℛcβ.
Making use of the results on (51) we have
(52)ℛh=θ(π0+π1)(μ+ω+δ)απ0+π1(α+μ+ω)(α+μ+ω)⇒ℛh(ℛc)=(μ+ω)θℛc(μ+ω+δ)β.
Taking the partial derivative of ℛh(ℛc) with respect to ℛc one gets
(53)∂ℛh(ℛc)∂ℛc=(μ+ω)θ(μ+ω+δ)β>0.
Results in (53) suggest that an increase in ℛc may result in an increase in ℛh, and this further suggests that whenever the two diseases coexist and are endemic, then they are high chances that they can influence the growth of each other. Following the above procedure in (53) it can be established that ∂ℛc(ℛh)/∂ℛh>0, which again suggests that whenever the two diseases coexist they may fuel one another.


### 5.2. Population Level Effects

In order to support analytical findings in this study, as well as demonstrating the impact of HIV and HCV coexistence on correctional institutions, numerical simulations of the model system (1) were performed using the MATLAB ODE solver, ode45, and parameter values in Table 1. It is worth noting that the scarcity and limited data on HIV and HCV within correctional institution limit our ability to calibrate and fit our model to the original data for future forecasts; nevertheless, we assume some of the parameters and initial conditions in the realistic range for illustrative purpose. These parsimonious assumptions reflect the lack of information currently available on the dynamics of the aforementioned epidemics with the correctional institutions. Reliable data on the two diseases (HIV and HCV) transmission within correctional institutions would enhance our understanding and aid in the possible interventions to be implemented.Figure7 illustrates the long term dynamics of single and dual cumulative infection cases. In Figure 7(a) we note that if Ic(0)>(Ih(0)>Ich(0)), then after a period of about 15 years cumulative dual cases will be more than any of the cumulative single infection cases, though cumulative HCV cases will remain higher than cumulative single HIV cases. In Figure 7(b) we observe that when Ih(0)>(Ic(0)>Ich(0)), then in a period of 5 years or less cumulative single HCV infections will be more than cumulative single HIV cases, while for a period of about 15 years cumulative dual cases will be more than either of the cumulative single infections. If Ih(0)=Ic(0)=Ich(0), then for the period 0–10 years cumulative HCV single infections will be dominant followed by cumulative dual infections; however, cumulative dual infections will outnumber cumulative single HCV infections 5 years later. Overall, results in Figure 7 suggest that in the presence of HCV and HIV within correctional institutions, the first periods will be dominated by HCV cases, but time dual cases will outnumber all single infections. This clearly shows that coexistence of HCV and HIV within correctional institutions is not a pleasant news to the public health domain.Simulations of model (1) showing the long term dynamics of Ic, Ih, and Ich for different initial conditions. The rest of the parameter values are in Table 1. In all the three figures Sn=500 and Sd=350, while (a) Ic=180, Ih=100, and Ich=50; (b) Ic=100, Ih=180, and Ich=50; (c) Ic=100, Ih=100, and Ich=100.
(a)
(b)
(c)Figure8(a) shows that increasing η will increase cumulative single HCV infections during the period 0–10 years which will also be associated with increased cumulative single HCV cases; thereafter the reverse will be true. The impact of η on cumulative single HIV and cumulative dual cases is shown in Figures 8(b) and 8(c), respectively; here we note that an increase in η is associated with an increase in cumulative infection cases.Simulations of model (1) showing the effects of varying the parameter η on cumulative single HCV cases, cumulative single HIV cases, and cumulative dual cases over a period of 40 years. The following assumed initial conditions were used. Sn=400, Sd=350, and Ic=Ih=Ich=100, and parameter values used are as in Table 1.
(a)
(b)
(c)The effects of increasingσ are shown in Figure 9. Here we observe that increasing σ may lead to a decrease in cumulative single cases, while leading to an increase in cumulative dual cases for the period 0–15 years, and thereafter high values of σ will be associated with low cumulative dual cases.Simulations of model (1) showing the effects of varying the parameter σ on cumulative single HCV cases, cumulative single HIV cases, and cumulative dual cases over a period of 40 years. The following assumed initial conditions were used. Sn=400, Sd=350, and Ic=Ih=Ich=100, and parameter values used are as in Table 1.
(a)
(b)
(c)


## 5.1. Impact of HCV on HIV Prevalence and Vice Versa

In many epidemiological models, the magnitude of the reproductive number is associated with the level of infection. The same is true for model (1). Expressions ℛc and ℛh can be rewritten as
(51)ℛc=β(π0+π1)(μ+ω)απ0+π1(α+μ+ω)(α+μ+ω)⇒απ0+π1(α+μ+ω)(α+μ+ω)=(π0+π1)(μ+ω)ℛcβ.
Making use of the results on (51) we have
(52)ℛh=θ(π0+π1)(μ+ω+δ)απ0+π1(α+μ+ω)(α+μ+ω)⇒ℛh(ℛc)=(μ+ω)θℛc(μ+ω+δ)β.
Taking the partial derivative of ℛh(ℛc) with respect to ℛc one gets
(53)∂ℛh(ℛc)∂ℛc=(μ+ω)θ(μ+ω+δ)β>0.
Results in (53) suggest that an increase in ℛc may result in an increase in ℛh, and this further suggests that whenever the two diseases coexist and are endemic, then they are high chances that they can influence the growth of each other. Following the above procedure in (53) it can be established that ∂ℛc(ℛh)/∂ℛh>0, which again suggests that whenever the two diseases coexist they may fuel one another.


## 5.2. Population Level Effects

In order to support analytical findings in this study, as well as demonstrating the impact of HIV and HCV coexistence on correctional institutions, numerical simulations of the model system (1) were performed using the MATLAB ODE solver, ode45, and parameter values in Table 1. It is worth noting that the scarcity and limited data on HIV and HCV within correctional institution limit our ability to calibrate and fit our model to the original data for future forecasts; nevertheless, we assume some of the parameters and initial conditions in the realistic range for illustrative purpose. These parsimonious assumptions reflect the lack of information currently available on the dynamics of the aforementioned epidemics with the correctional institutions. Reliable data on the two diseases (HIV and HCV) transmission within correctional institutions would enhance our understanding and aid in the possible interventions to be implemented.Figure7 illustrates the long term dynamics of single and dual cumulative infection cases. In Figure 7(a) we note that if Ic(0)>(Ih(0)>Ich(0)), then after a period of about 15 years cumulative dual cases will be more than any of the cumulative single infection cases, though cumulative HCV cases will remain higher than cumulative single HIV cases. In Figure 7(b) we observe that when Ih(0)>(Ic(0)>Ich(0)), then in a period of 5 years or less cumulative single HCV infections will be more than cumulative single HIV cases, while for a period of about 15 years cumulative dual cases will be more than either of the cumulative single infections. If Ih(0)=Ic(0)=Ich(0), then for the period 0–10 years cumulative HCV single infections will be dominant followed by cumulative dual infections; however, cumulative dual infections will outnumber cumulative single HCV infections 5 years later. Overall, results in Figure 7 suggest that in the presence of HCV and HIV within correctional institutions, the first periods will be dominated by HCV cases, but time dual cases will outnumber all single infections. This clearly shows that coexistence of HCV and HIV within correctional institutions is not a pleasant news to the public health domain.Simulations of model (1) showing the long term dynamics of Ic, Ih, and Ich for different initial conditions. The rest of the parameter values are in Table 1. In all the three figures Sn=500 and Sd=350, while (a) Ic=180, Ih=100, and Ich=50; (b) Ic=100, Ih=180, and Ich=50; (c) Ic=100, Ih=100, and Ich=100.
(a)
(b)
(c)Figure8(a) shows that increasing η will increase cumulative single HCV infections during the period 0–10 years which will also be associated with increased cumulative single HCV cases; thereafter the reverse will be true. The impact of η on cumulative single HIV and cumulative dual cases is shown in Figures 8(b) and 8(c), respectively; here we note that an increase in η is associated with an increase in cumulative infection cases.Simulations of model (1) showing the effects of varying the parameter η on cumulative single HCV cases, cumulative single HIV cases, and cumulative dual cases over a period of 40 years. The following assumed initial conditions were used. Sn=400, Sd=350, and Ic=Ih=Ich=100, and parameter values used are as in Table 1.
(a)
(b)
(c)The effects of increasingσ are shown in Figure 9. Here we observe that increasing σ may lead to a decrease in cumulative single cases, while leading to an increase in cumulative dual cases for the period 0–15 years, and thereafter high values of σ will be associated with low cumulative dual cases.Simulations of model (1) showing the effects of varying the parameter σ on cumulative single HCV cases, cumulative single HIV cases, and cumulative dual cases over a period of 40 years. The following assumed initial conditions were used. Sn=400, Sd=350, and Ic=Ih=Ich=100, and parameter values used are as in Table 1.
(a)
(b)
(c)


## 6. Discussion

HIV and hepatitis C virus (HCV) infections are among the most costly consequences of illicit drug use, having a high impact on individuals and on health care systems. Prison populations are considered to be at high risk for blood borne infections due to the high proportion of intravenous drug users, commercial sex workers and homeless people, high-risk sexual behaviors before and during incarceration, and tattooing among inmates [47]. A mathematical model for assessing the impact of HIV and HCV coexistence within correctional institutions has been developed and analyzed. The basic reproductive number for the model has been computed and qualitatively used to gain insights on the long term dynamics on cumulative single and dual cases when the two diseases coexist with correctional institution. Sensitivity analysis on the associated reproductive number has been carried out, and the influence of each of the parameters which define it (the reproductive number) has been clearly illustrated. Analysis of the two submodels (the HIV submodel and HCV submodel) suggest that each submodel has globally asymptotically disease-free equilibrium whenever the associated reproductive number is less than unity, and a globally asymptotically stable endemic when the associated reproductive number is greater than unity. Comprehensive numerical simulations performed using the MATLAB ODE solver, ode45, have been provided to support analytical results. At its best the study suggests that whenever HIV and HCV coexist within correctional institutions, then the long-term cumulative cases will be dominated by dual infections, followed by cumulative single HCV cases, and this clearly suggests the need for resources and attention in order to control the aforementioned epidemics since the consequences will also be felt by the population outside the correction institutions.



---

*Source: 919502-2012-12-17.xml*