Hepatitis C virus modelled as an indirectly transmitted infection highlights the centrality of injection drug equipment in disease dynamics

1. Introduction

While the ecology of infectious disease is a rich field with decades worth of empirical evidence and theory, there are aspects that remain relatively under-explored. One example is the importance of the free-living survival stage of certain pathogens, where diseases are transmitted indirectly between hosts through an environmental reservoir intermediate. These include infections transmitted indirectly between hosts via a surface or reservoir intermediate—often abiotic—where the pathogen lives freely and independently of a host [1–18], sometimes described as ‘sit and wait’ infections [19]. Other studies have focused on systems where pathogens are growing in the environment [9], or have explored indirectly transmitted infections in theoretical terms [12,15]. While frameworks already exist for studying indirect environmental transmission, most are engineered with constraints that render their application necessarily narrow [6], limiting their relevance for a wider number of indirectly transmitted infections.

One class of diseases where the indirect transmission paradigm has been scarcely applied are those spread through injection drug use in urban settings, such as the human immunodeficiency virus (HIV) and hepatitis C virus (HCV). HIV has been the object of many important mathematical models [20,21], some of which have implemented injection drug use effectively, even focusing on the specific dynamics of injection equipment [22–25]. HCV has also been studied using modelling methods, many focusing on treatment [26–28] and others on the particulars of transmission in injection drug-user communities [29–35]. Importantly, none of these existing dynamical models consider the peculiar ecology of HCV transmission, where transmission events occur through an environmental reservoir (injection equipment) that resembles a disease vector [36,37]. Unlike an insect vector, however, injection equipment is not an organism and is more realistically considered an abiotic reservoir for infection, similar to the role that the water supply serves in an outbreak of cholera or other waterborne diseases [38]. As HCV continues to pose a serious public health challenge in many communities, there is a need to understand how the dynamics of injection equipment influence HCV transmission. This is especially important for informing the utility of harm reduction programmes, such as needle exchange, which have been effective in decreasing transmission of HIV and HCV [39,40]. Lastly, but perhaps most importantly, the urgency for understanding these dynamics has increased dramatically in recent years with the growth of the modern opioid epidemic, much of it involving injection drug use [41,42]. The lack of models of HCV that specifically consider injection equipment, and increased social urgency related to the modern opioid epidemic implore more adaptable mathematical models of injection-drug use that could facilitate a better understanding of and predictions for the trajectory of modern HCV infections.

In this study, we model HCV as an indirectly (or environmentally) transmitted infection, where the drug paraphernalia serves as the environmental reservoir. As HCV epidemics are partly defined by injection drug users and injection drug equipment, we argue that this indirectly transmitted lens captures aspects that prior models have not. In §2, we introduce a theoretical iteration of an indirectly transmitted infection using a standard epidemiological model imbued with an environmental reservoir compartment. We describe analytic equations of such a system, and derive the reproductive number (R₀) using analytic methods. Then in §3, we introduce the full HCV mathematical model, demonstrating how it allows one to examine several otherwise-overlooked features of disease dynamics. We pontificate on these results in light of the ecology of infectious diseases, and in terms of public health policies, especially as they relate to the modern opioid epidemic.

2. An elementary adapted SIR indirectly transmitted iteration

2.1. Description

While the emphasis of our examination will reside in how we analyse a HCV epidemic, for explanatory purposes we will begin by describing how an environmental reservoir modifies very basic concepts in a classic, purposefully prosaic susceptible–infected–recovered (SIR) mathematical model. We will explain the basic structure of a model of indirect transmission, after which the HCV-specific iteration will be discussed.

While there are several existing frameworks that can be used to describe infections spreading through an environmental reservoir, we have conveniently labelled ours the waterborne, abiotic and indirectly transmitted (W.A.I.T.) infection model. Many diseases can be modelled using this kind of approach, but this study applies it to HCV in a community of injection drug users, which has not been previously modelled in this manner.

We use a standard SIR framework, where dynamics are defined by changes in a population of susceptible (S), infected (I) and recovered (R) hosts. Classically, flow of infection through the system is defined by contact between susceptible and infected individuals, often driven by a β factor, or transmission coefficient. Figure 1 is a compartmental model that depicts this interaction, and adds two additional compartments, labelled with a W (for W.A.I.T.), which influence the flow of hosts from the susceptible to infected compartments—indicated by the dashed lines in the figure.

• Download figure
• Open in new tab
• Download PowerPoint

2.2. The adapted SIR compartmental diagram

The S, I and R compartments represent the usual susceptible, infected and recovered populations of hosts, respectively. W_u and W_i represent uninfected and infected populations of environmental agents, respectively.

In traditional SIR models, the rate of new infection (arrow from the S compartment to the I) is generally proportional to the product of the susceptible and the infected populations, i.e. proportional to SI. In the W.A.I.T. framework, the environmental compartment plays a role analogous to the infected host compartment in driving the rate of infection. In this specific example, the W_i compartment contributes to the rate of infection as a fraction, W_i/(W_i + W_u), which appears as a factor in the rate terms.

The epidemic is then driven by a series of interactions: between uninfected (susceptible) hosts S and the infected (transmitting) environmental compartment W_i, and interactions between infected individuals I and the uninfected environmental compartment W_u. The epidemic is sustained through infected hosts I depositing pathogen into the environmental reservoir, creating new infections, which can then infect more susceptible hosts S (in a process resembling a feedback loop). These dynamics can be captured by the set of dynamical equations and visualized with the diagram in figure 1:

Equations (2.1)–(2.5) define an extension of the prosaic SIR model. π_S is the birthrate of new susceptible hosts and μ is the fractional death rate of hosts. In this context, β represents the strength of the interaction between the susceptible hosts S and the environmental reservoir. This will generally be proportional to the rate of contact between the two. Similarly, α characterizes the strength of interaction between infected hosts I and the environmental reservoir, and is also generally proportional to the contact rate between the two. α and β, while both generally proportional to the contact rate between environmental agents and living hosts, are distinguished by factors that indicate the probabilities of spreading the infection either from host to environment, as in the case of α, or from the environment to host, as in the case of β. Note that α and β could be replaced with the same parameter in settings where the infection is guaranteed to spread at any encounter with an infected host or environmental agent. ν represents the fractional recovery rate, π_W is the birthrate of new uninfected environmental agents and k is the fractional death or discard rate of environmental agents. Note that the discard rate k (which includes any force that removes environmental agents from the system) can be split into two discard rates, one for the infected compartment, and one for the uninfected compartment (we do, in fact make this distinction in the full HCV model). For simplicity, we will tend to set the discard rates of these two compartments equal, where there are two parameters, as we have no current mathematical grounds to distinguish them.

Also note that this model resembles vector-borne transmission models such as those used to study malaria [37]. In fact, prior studies have explored the utility of applying vector-borne transmission models to the spread of infection with needles as a proxy for vectors—although, only in the context of HIV [36]. In this paper, however, we wish to emphasize peculiarities of the spread of HCV and, further, to elaborate some features of the R₀ expression in the context of these shared dynamical settings—between hosts and agents—which we believe have not been rigorously addressed in the existing literature.

2.3. W.A.I.T. framework influences the basic reproductive number in a standard SIR model

Next, we briefly consider how the value of the basic reproductive ratio R₀ in this model compares to its SIR counterpart. While R₀ can have different theoretical formulations, we rely on definitions as provided by Jones [43] and Diekmann et al. [44]. In a density-dependent SIR model with constant birth of susceptible hosts π_S and death rate proportional to the host population −μS, the R₀ value is given by

or sometimes, more simply, $R_0^{SIR}=\beta /\nu$, depending on the form of the SIR equations used, e.g. frequency-dependent, constant population, etc. In this equation, β is the traditional transmission coefficient. It represents the coupling strength between infected and uninfected hosts, two non-environmental agents. Whereas, in the W.A.I.T. model, what is analogous to β is a pair of parameters α and β, which govern the interaction strengths between hosts and the environment. π_S, μ and ν have the same interpretation as in the W.A.I.T. model.

In the case of the W.A.I.T. iteration, the value of R₀ takes the form

There are some notable differences in the R₀ formulae of the SIR and W.A.I.T. models: the square root in the W.A.I.T. version arises as a consequence of implementing two infected agents (I and W_i) into the model, as opposed to just one in the SIR case. Next, one notices that the β factor in the SIR formula is augmented by the additional factor α in the W.A.I.T. formula, representing a kind of shared dependence between the couplings controlling the I-interaction (α) and the S-interaction (β) with the environment. Analogously, what was the responsibility of π_S in the SIR formula now presents itself as a shared dependence, π_S/π_W, the ratio of the birthrate of susceptible hosts to that of uninfected environmental agents. In this case, the two appear as a ratio under the square root, as opposed to a product in the αβ case, indicating that whereas α and β contribute to R₀ in the same way, π_S and π_W contribute in opposite ways: when π_S is increased, R₀ increases, but when π_W is increased, R₀ decreases.

It is possible to view $R_0^{\mathrm{WAIT}}$ as a geometric mean of two R₀ values. Namely, there is the reproductive ratio associated with the number of secondary host infections caused by a single infected environmental agent, and there is the reproductive ratio associated with the number of secondary environmental agent infections caused by a single infected host. We denote the former by $R_0^H$ and the latter by $R_0^W$ (H for hosts and W for the W.A.I.T. compartment). From equations (2.1)–(2.5), one can see that the rate of new host infection due to infected environmental agents W_i is given by βSW_i/(W_i + W_u). Near the disease-free equilibrium (DFE), S ≈ π_S/μ and W_i/(W_i + W_u) ≈ kW_i/π_W (near the DFE, W_i ≪ W_u), which implies that near the DFE, the rate of new host infection per infected environmental agent is ≈βπ_Sk/(μπ_W). The average amount of time an infected environmental agent remains infected is 1/k, i.e. the reciprocal of the exit rate of the infected state. Thus, the number of new host infections caused by an infected environmental agent in the time that the agent is infected, and while the system is near the DFE, is given by βπ_Sk/(μπ_W) × 1/k = βπ_S/(μπ_W). That is

Similarly, the rate of new infection of environmental agents, caused by infected hosts, is given by αIW_u/(W_i + W_u). Near the DFE, this rate, per infected host, is ≈α (since W_u/(W_i + W_u) ≈ 1), and the average time that an infected host remains infected is given by 1/(μ + ν), the reciprocal of the exit rate of the infected state. Thus, the number of new environmental agent infections caused by an infected host in the time that the host is infected (near the DFE) is given by

One can see that the value of R₀ given in equation (2.7) is the geometric mean of the two R₀ values calculated above:

From this perspective, one can observe how a characteristic feature of the epidemic is modified by indirect transmission.

3. The hepatitis C virus model

3.1. Description

Our HCV model represents an adaptation of the SIR W.A.I.T. model outlined in §2, but engineered around the particulars of HCV. Our model simulates a population of approximately 170 000 individuals—based on estimates of the size of the people who inject drugs (PWID) community in New York City [45]—where infected injection drug users may migrate into the population. In this model, injection paraphernalia serve as the environmental reservoir for HCV and the sharing of this equipment will constitute the means of transmitting new infections. While the entirety of injection paraphernalia might contain other components, many parameters in this model are based on the use of needle and syringe as the instrument of injection and sharing. Consequently, we use the term ‘needle’ in this paper as a synecdoche for the entire injection apparatus. It is also important to note that HCV can be transmitted sexually [46], but in this study we restrict our attention to transmission through infected needles. This main text focuses on the main structure and dynamical properties of the model. Further model details and discussion can be found in the electronic supplementary material, appendix.

3.2. HCV W.A.I.T. model: compartmental diagram

We model the dynamics of needle populations and injection drug users through a series of five ordinary differential equations. The compartments labelled S, I_E, I_L, N_u and N_i represent the populations of susceptible individuals, early-stage infected individuals (acute HCV infection), late-stage-infected individuals (chronic HCV infection), uninfected needles and infected needles, respectively (figure 2). Here, we refer to all needles in circulation within the entire PWID community. This model is defined by several features:

—	The susceptible compartment refers to individuals who are injecting drugs and who are sharing needles with other members in the PWID community.
—	The needle population is divided into two compartments: infected and uninfected, and we model the dynamics of each compartment separately. This is analogous to the Wi and Wu terms discussed in the preliminary model.
—	New infections (of both hosts and needles) will depend on the fraction of infected or uninfected needles in circulation.
—	Newly infected individuals enter the early stage compartment IE first before either spontaneously clearing the infection or moving into the late-stage compartment IL, from which we assume no spontaneous clearance occurs—individuals may leave IL either by treatment or death only, since cases of spontaneously clearing chronic HCV are rare.
—	There are various estimates for the ability of HCV to survive in needles [47,48]. We incorporate HCV free-living survival via the parameter ε, which quantifies the rate at which the virus decays on infected needles.

• Download figure
• Open in new tab
• Download PowerPoint

3.3. HCV W.A.I.T. model: analytic equations and parameters

The dynamics of the HCV transmission process are governed by equations (3.1)–(3.5). The population of individuals that are being treated and those who have recovered are not explicitly modelled in this iteration, as the dynamics of treatment and recovery are not central to the questions explored in this study. There are, however, several modelling studies of HCV that focus on treatment [26–28,49], and their effects are not ignored in the HCV W.A.I.T. model.

Entering treatment (and re-entering the susceptible population, as in case of drug relapse in the PWID population) are incorporated via removal terms −τI_L and −τI_E and the susceptible ‘birth’ term π_S.

where π_S is the birthrate of new members into the PWID community either via migration, first-time use or recovery from treatment—not from spontaneous self-clearance. ϕ represents the daily fractional rate that acutely infected individuals (or early infected I_E) spontaneously clear the infection—i.e. without treatment. α represents the per capita injection rate, scaled by the fraction of injection events by infected users that render a needle infectious. β represents the per capita injection rate, scaled by the fraction of injection events with an infected needle that leave a susceptible host infectious. μ is the combined fractional death and PWID cessation rate (individuals who leave the PWID community). ω is the daily fractional rate that early stage infected individuals progress to the late stage of infection. τ is the daily fractional rate that infected individuals go into treatment. π_N is the rate of introduction of uninfected needles into the PWID population. k_u is the daily fractional discard rate of uninfected needles. k_i is the daily fractional discard rate of infected needles. Lastly, $ϵ$ is the daily fractional rate that infected needles clear the infection due to de-activation (or ‘death’) of virus populations on the needle. Parameter values and sources can be seen in table 1.

3.4. HCV W.A.I.T. model parameters influence R₀

Having constructed and elaborated on the details of the HCV W.A.I.T. model, we now explore how parameters related to the environmental reservoir (in this case, those framing the population of infected needles) influence R₀. We directly measured the influence of parameters on R₀ by considering the partial rank correlation coefficient (PRCC), discussed below. The value of R₀ was calculated using established methods [43,44] and is outlined in the electronic supplementary material, appendix:

We emphasize that in a manner analogous to our example discussed in §2, we can regard our R₀ value as a geometric mean of two other R₀ values:

The left-most factor (under the square root) can be interpreted as the number of secondary infections of needles in the average time that a host is infected (near the DFE), and the right-most factor can be regarded as the number of secondary infections of hosts in the average time that a needle remains infected. Further discussion of this result can be found in the electronic supplementary material, appendix. As with traditional values of R₀, we find that our value is consistent with the statement that sign(R₀ − 1) = sign(λ), where λ is the maximal eigenvalue of the Jacobian of the infected subsystem—composed of the infected compartments of the ODE system: I_E, I_L and N_i—calculated at the DFE (all eigenvalues of the Jacobian were real-valued). This shows that the DFE is unstable when R₀ > 1.

We determine the sensitivities of our parameters on the value of R₀ by calculating the PRCC with respect to equation (3.6)—we base our calculation of PRCC on methods used in prior studies [60]. We find that parameters related to an interaction with the environmental reservoir (the population of needles) such as α and β, the couplings between hosts and needles, are at least as central to HCV dynamics as parameters traditionally associated with an epidemic, such as π_S, the birthrate of susceptibles, μ, the combined death and cessation rate of PWID, and τ, the rate of progressing to treatment (figure 3). This fortifies the notion that W.A.I.T.-specific properties dictate the spread of HCV, providing opportunities to explore more precise targeting by public health interventions.

• Download figure
• Open in new tab
• Download PowerPoint

3.5. HCV W.A.I.T. model and simulated interventions: needle exchange programmes

Having demonstrated the relevance of injection drug equipment in terms of how it influences the basic reproductive number, we can consider the utility of the model with respect to other properties, including how it offers insight into potential interventions.

In figure 4b, we demonstrate how changing k_u and k_i modifies the value of R₀. Notice that R₀ is reduced by increasing k_i across fixed values of k_u, and the opposite effect—increasing R₀—is observed when increasing k_u along fixed values of k_i. That is, removing infected needles at an increased rate may decrease infection risk in a population of PWID, while removing uninfected needles can increase the risk. One can also see that increasing k_u and k_i simultaneously, along the dashed line—where k_u = k_i—will increase R₀. This suggests that if a distinction between infected and uninfected needles cannot be established (as is often the case) then discarding needles indiscriminately can potentially exacerbate the spread of the infection.

• Download figure
• Open in new tab
• Download PowerPoint

This is a result of the fact that when k := k_u = k_i, the expression for R₀ is proportional to

We highlight this to show the explicit dependence on k in the R₀ expression. Notice that this factor is monotonically increasing in k, indicating that no matter the values of other parameters in the model, as long as they are all positive, R₀ will necessarily increase with k. Notice also that when $ϵ$ = 0, the k dependence cancels out entirely. This indicates that when $ϵ$ = 0, meaning that there is no flow of infected needles back to the uninfected compartment, then R₀ is not modified by the discard rate of needles. We point the reader to the electronic supplementary material, appendix, for a more thorough discussion of this point.

Next, we considered how certain interventions can modify the transfer rate of needles from infected to uninfected states, through modifying the $ϵ$ parameter in our study (figure 5). A high $ϵ$ value would indicate a scenario where needles move quickly from an infected state to an uninfected state. This would apply to settings where viral decay on a needle is high, or when infected needles are directly exchanged for uninfected ones (as in certain needle exchange programmes). The model is run with all uninfected populations initialized at their DFE values (S = 170 000 and N_u = 220 000), and we initialize I_E = N_u = 1, and I_L = 0. In the high $ϵ$ scenario, we observe generally slower dynamics and higher overall susceptible population sizes, along with lower infected populations (on long time scales).

• Download figure
• Open in new tab
• Download PowerPoint

4. Discussion

While diseases transmitted through injection drug use have been the object of prior modelling efforts, none have specifically investigated how injection equipment plays a role in the dynamics of HCV. Prior models of injection equipment have focused on HIV [24,25], and/or been so complicated that their structure is not easily translated to any other settings [23]. In this study, we model HCV as an indirectly transmitted infection, where the injection equipment is modelled as the environmental reservoir, just as a water source might be modelled in a waterborne infection [8,19]. We label our approach as the W.A.I.T. model, one that incorporates features of other approaches to studying environmentally transmitted pathogens [6,11], but grounding them in a flexible model that can be neatly applied to HCV. Our approach offers several specific insights. For example, we demonstrate that the composite R₀ that defines the entire dynamical system is the geometric mean of R₀ values used to describe each of two sub-components: disease flow through the hosts and flow through the injection equipment (equation (3.7)). This observation offers a practical suggestion for studying diseases like HCV: epidemiologists and modellers must understand, through empirical studies, properties of all major actors in the system (hosts and environmental injection drug equipment in the case of HCV).

The mathematical model of HCV presented in this paper (described as a W.A.I.T. model; see §§2 and 3) also offers nuanced findings about the dynamics of disease. Firstly, our model highlights the differing roles of uninfected and infected injection on disease dynamics. Specifically, the model speaks to the potential utility of harm reduction policies: indiscriminately removing injection equipment from a system—without an overall shift in needle populations from infected to uninfected—might increase the rate of infection. In order to attenuate an epidemic, intervention strategies should focus on steering the population of needles towards being more uninfected. Therefore, ideal intervention efforts should aim to decrease sharing events on an infected needle. This helps to explain why programmes like safe injection might be effective [61]: they do not change the number of infected needles in the system directly, but can alter the sharing rate, and consequently, the probability of sharing an infected needle.

Finally, understanding the dynamical properties of disease transmitted through injection drug use is now especially relevant as a result of the modern opioid epidemic. This epidemic is typified by recreational use of prescription and illicit opioids, with injection drug use being a major route through which drugs are consumed [42]. The relevance of viral diseases among opioid users gained national attention during a 2015 outbreak of HIV in rural Indiana that was driven by an injected opioid called oxymorphone [62,63]. This outbreak raised alarms in the public health community, and officials are increasingly aware of the potential for future outbreaks. However, it was not until relatively recently that the role of the opioid crisis in HCV transmission has been examined [64,65]. We propose, in closing, that modelling approaches (in general, and not relegated to the methods proposed in this study) are crucial for understanding, attenuating or preventing explosive outbreaks of HCV in an age when a new opioid epidemic has emerged.

Data accessibility

Code for the mathematical models presented in this paper are available on GitHub (https://github.com/ogplexus/WAIT).

Authors' contributions

All authors contributed to the development of the model, study design, analysis and interpretation of results. C.B.O. conceived the study, and created the original iteration of the model. M.D.M.D. and V.A.M. wrote the code, ran simulations, performed calculations and developed visualizations. All authors contributed to the writing and editing of the manuscript.

Competing interests

We declare we have no competing interests.

Funding

This work was supported by Brown University.

Footnotes

^†These authors contributed equally to this study.

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.4616084.