The effect of the needle exchange program on the spread of some sexually transmitted diseases

1.   Introduction

Besides their principal transmission form, some sexually transmitted diseases (STD) can also be spread via several other ways, including donor tissue, blood, breastfeeding or during childbirth. In this paper we shall focus on STDs which can be transmitted among intravenous drug users (IDU) through sharing syringes or needles during drug use. Typical examples for pathogens spread this way are human immunodeficiency virus (HIV), hepatitis B, C, D virus, as well as the bacterium Treponema pallidum causing syphilis ^[1,2,3,4]. At the same time, some IDUs also engage in high-risk sexual behaviour facilitating the transmission of these diseases between different groups ^[5,6]. Hence, the application of infected needles not only elevates the risk of transmission among IDUs, but also creates a greater risk for non-drug users.

Recently, the European Centre for Disease Prevention and Control reported a sharp increase of HIV notifications in Greece and Romania due to the reduction of the sterile needle exchange program among drug users ^[7]. Several other European countries – Belgium, Cyprus, Hungary, Slovakia, Turkey – were rated as endangered based on the number of syringes given out per IDU per year.

Several studies have already been published on STDs with the possibility of transmission via blood ^[9,10,11,8]. In ^[9] and ^[8], the authors considered different compartments for people with different rates of addiction to intravenous drugs. Here we enhance the models developed by Kaplan et al. ^[12,13,14]. In these works, the authors established models to examine the dynamics of the population of susceptible and contaminated needles as well. Based on these models, further research was carried out by Greenhalgh et al. ^[15,16,17]. Several papers deal with a special subetaoup of the population, e.g. people visiting shooting galleries ^[18] or a rare transmission way, e.g. parental exposure and blood transfusion ^[19]. Ji et al. ^[11] studied optimal needle allocation to prevent HIV transmission.

Motivated by the above reports and predictions, our main purpose is to investigate the dynamics of a model including the effect of needle exchange and study the risk of an increased transmission among non-drug users, induced by the reduction of the needle exchange program. Our new approach to this problem consists in considering not just a needle population and the prevalence of infections among drug users, but beyond the latter one, a population of sexually active humans consisting of drug users and non-drug users.

The structure of the paper is the following. In the next section, we establish a system of differential equations to model the spread of a disease through sexual contact and needle sharing. In Section 3, we calculate the disease-free equilibrium and the basic reproduction number $\mathcal{R}_0$. In Section 4, we show that the disease-free equilibrium is globally asymptotically stable if $\mathcal{R}_0 < 1$. In Section 5, we show that the infected compartments are uniformly persistent if $\mathcal{R}_0 > 1$.

In Section 6, we provide some numerical simulations which suggest that the reduction of the needle exchange program can highly contribute to the spread of diseases both among drug users and non-drug users. We fit a time-dependent variant of our model to recent data from Hungary to give a prognosis for the expected number of infections in the next decade.

2.   The model

We establish a compartmental model, in which, besides humans, also needles are considered as a population. In our model, the human population is divided into four isolated groups: $S$ denotes susceptible individuals who do not use drugs, $S_D$ stands for susceptible IDUs. We denote by $I$ those who are infected and do not use drugs, while $I_D$ stands for infected IDUs. The needle population is separated to two disjoint classes, $S_N$ denotes susceptible, clean needles, while $I_N$ stands for infected needles. As in a developed country more than 99% of the cases happen via sexual transmission or contaminated needles, in our paper we concern only these two major factors, and investigate a model emphasizing their role (see e.g. ^[20]).

We assume that both human and needle populations are closed in the sense that birth and death rates (denoted by $\lambda$) are balanced ($S + S_D + I + I_D = \mathrm{const.}, \ S_N + I_N = \mathrm{const.}$) and that the disease spread by people-people or people-needles contacts.

The transitions between the compartments are as follows. All newborn individuals enter the $S$ compartment. Susceptibles become addicted to intravenous drugs and hence move from the compartment $S$ to $S_D$ with rate $a$. The other way around, $g$ denotes the rate of giving up drug use. Getting addicted to and giving up drugs is modelled similarly for the infected compartments $I$ and $I_D$.

The model includes two ways of contracting HIV. Sexual intercourse with an infected partner results in getting infected with a rate $\beta$. The other way of infection considered in our model is the use of contaminated needles among IDUs. We consider a needle population where new needles are distributed with a constant rate $\nu$, and we assume that old needles are discarded with the same rate $\nu$. Transmission rate from infected drug users to clean needles is denoted by $\omega$, while $\alpha$ denotes the transmission rate from infected needles to susceptible drug users. Used needles are sterilized with rate $\theta$. The transmission diagram of the model is shown in Figure 1.

Using the above notations and assumptions, we obtain the following system of ordinary differential equations:

3.   Disease-free equilibrium and basic reproduction number

In this section, we study the disease-free equilibrium which can be obtained by solving the following algebraic system of equations:

The unique disease-free equilibrium can be calculated as

The most important parameter to characterize an epidemic is the basic reproduction number $\mathcal{R}_0$ defined as the expected number of secondary infections caused by a single infectious individual over the course of his/her infectious period, introduced in a completely susceptible population.

We apply the method of next generation matrices introduced by Diekmann, Heesterbeek and Roberts ^[22]. According to this method, the basic reproduction number is given by the dominant eigenvalue of the next generation matrix $-FV^{-1}$, where $F$ is the transmission matrix and $V$ is the transition matrix of the infection subsystem at the disease-free equilibrium. In our case, this subsystem is given by

hence, the transmission and transition matrices are

and

Therefore, the next generation matrix is calculated as

The characteristic equation of the next generation matrix has the form

Remark 3.1. It is easy to see that there are three distinct real eigenvalues, two positive and one negative. Observe the cubic Vieta's formulas:

Let us suppose that $x_1 = \mathcal{R}_0$ and $x_2, x_3$ are complex conjugate roots. The product $\mathcal{R}_0x_2x_3$ is negative according to the third equation, which leads to a contradiction, since multiplying $\mathcal{R}_0$ by a complex conjugate pair is always positive. Moreover, the product of three real roots will be negative only if the number of the negative terms is odd.

Using Cardano's formula for cubic equations we obtain

The formula for $\mathcal{R}_0$ is far too complicated to be easily applicable in the analysis. In order to cope with this difficulty, we introduce the parameter $\mathcal{P}_0$ as follows:

Although this expression is much simpler, one can see that the parameter $\mathcal{P}_0$ is not equivalent to $\mathcal{R}_0$ as a threshold parameter: the smaller positive eigenvalue may pass through $1$ as well. It is easy to see that $\mathcal{R}_0 < 1$ implies $\mathcal{P}_0 < 0$ and $\mathcal{P}_0 > 0$ implies $\mathcal{R}_0 > 1$ (see Figure 2).

4.   Global asymptotic stability of the disease-free equilibrium

In this section, we will show that the basic reproduction number calculated in the previous section serves as a threshold parameter regarding the dynamics of our system. We need to show some simple assertions before we can state the main theorem of the section. Let $f\colon[0, \infty) \to \mathbb{R}$. We introduce the notations

for the limit superior, resp. the limit inferior of $f$ as $t \to \infty$.

Lemma 4.1. For the limit superior of $S, S_D, S_N$ the following inequalities hold:

Proof. Adding the first two equations of (2.1), we obtain

Since $1-S_D(t)-I_D(t) = S(t)+I(t)$ we get

According to the fluctuation lemma (see e.g. ^[28]), we have

Since $S$ and $I$ are nonnegative, $S^\infty \leq \left(S+I\right)^\infty$ holds.

Similarly, let us consider $\left(S_D(t)+I_D(t)\right)'$:

from which we have

Since $S(t)+I(t)+S_D(t)+I_D(t) = 1$ we get

Using the fluctuation lemma we obtain

Since $S_D$ and $I_D$ are nonnegative, $S_D^\infty \leq (S_D+I_D)^\infty$ holds.

$S_N^\infty\leq 1$ holds trivially.

In the proof of the global asymptotic stability of the disease-free equilbrium $E_0$ we will need [27,Theorem 2.1]. In order to apply this result, we first recall this theorem.

Theorem A. ([27,Theorem 2.1]). Let us consider a system

with $g = (g_1, \dots, g_m)^T$. The vectors $x = (x_1, \dots, x_n)\in\mathbb{R}^n$ and $y = (y_1, \dots, y_m)^T\in\mathbb{R}^m$ stand for the populations in the disease compartments and the non-disease compartments, respectively. The functions $\mathcal{F}$ and $\mathcal{V}$ are defined as $\mathcal{F} = (\mathcal{F}_1, \dots, \mathcal{F}_n)^T$ and $\mathcal{V} = (\mathcal{V}_1, \dots, \mathcal{V}_n)^T$ such that $\mathcal{F}_i$ denotes the rate of new infections in the $i$th disease compartment and $\mathcal{V}_i, \ i = 1, \dots, n$ denote transition terms. We assume that the disease-free system $y' = g(0, y)$ has a unique equilibrium $y_0 > 0$ which is locally asymptotically stable within the disease-free space. Let us define the two $n\times n$ matrices $F$ and $V$ as

and let

Let $\omega^T\geq 0$ be the left eigenvector of the matrix $V^{-1}F$ corresponding to the eigenvalue $\mathcal{R}_0$. If $f(x, y)\leq 0$ in $\Gamma\subset\mathbb{R}_+^{n+m}$, $F\geq 0$, $V^{-1}\geq 0$ and $\mathcal{R}_0\leq 1$, then $Q = \omega^TV^{-1}x$ is a Lyapunov function for (4.1) on $\Gamma$.

Theorem 4.2. The disease-free equilibrium of system (2.1) is globally asymptotically stable if $\mathcal{R}_0 < 1$ and it is unstable if $\mathcal{R}_0 > 1$.

Proof. The matrices $F$ and $V$ are given in (3.1) and (3.2), and one can easily calculate

Using the notations of Theorem 4.2, and introducing the notations $x = (S, S_D, S_N)$ and $y = (I, I_D, I_N)$, for model (2.1) we have

Let us now proceed as seen in Theorem 4.2: let $\omega^T$ be the left eigenvector of the matrix $V^{-1}F$ corresponding to the basic reproduction number $\mathcal{R}_0$ and define $Q = \omega^TV^{-1}x$. Then, proceeding as in the proof of Theorem 4.2 in ^[27], differentiating $Q$ along solutions of (2.1), we obtain

From the above, we can see that the conditions $F > 0$ and $V^{-1} > 0$ hold, while the condition $f(x, y)\leq 0$ does not always hold: the first two coordinates are positive if and only if the conditions

hold. However, from Lemma 4.1 we know that for any $\varepsilon > 0$, there exists a $T$ large enough such that

for all $t > T$. Hence, for $t > T$, the (possibly) positive terms in $Q'$ can be made arbitrarily small, i.e. the derivative is negative for $t$ large enough. From LaSalle's invariance principle (see e.g. [26]) we know that the limit set of each solution is a subset of the set $\big\{\frac{\mathrm{d}}{\mathrm{d}t} Q(t) = 0 \big\}$. Hence, on the limit set, the equation (2.1) reduces to

It is straightforward that $\lim_{t\to\infty}S_N(t) = 1$, while applying the Dulac–Bendixson criterion with the Dulac function $\frac{1}{S_D}$ shows that all solutions of the subsystem formed by the first two equations of the latter limit system tend to the equilibrium $\big(\frac{g+\lambda}{a+g+\lambda}, \frac{a}{a+g+\lambda}\big)$, hence, we have shown the first statement of the theorem.

Now, let $\mathcal{R}_0 > 1$. The Jacobian of the system (2.1) at the disease-free equilibrium is

The characteristic equation is

where

As the leading coefficient of $f(x)$ is positive, in order to show that there exists at least one eigenvalue with positive real part, according to the Routh–Hurwitz theorem, we need to show that at least one of the further coefficients is negative. We will show that this holds for at least one of the last two coefficients.

In the case $\mathcal{P}_0 > 0$, the constant term of (4.2) is negative.

In the case $\mathcal{P}_0 < 0$, suppose that the assertion does not hold, i.e. both the constant term and the coefficient of the first-order term, i.e.

are positive.

Taking into account that $\mathcal{R}_0 > 1$, $\chi\big(\frac{\beta}{\lambda}\big) > 0$ implies $\beta > \lambda$ (see Figure 3). From this it follows that apart from the third one, all terms of (4.3) are negative. Thus, if (4.3) is positive, then the inequalities

and

hold. Let us now suppose that $\mathcal{P}_0 < 0$, i.e.

holds (from which the positivity of the constant term would follow). From this we obtain

which yields

contradicting (4.4).

5.   Uniform persistence of the infected compartments

In this section, we focus on the persistence of the disease in the case $\mathcal{R}_0 > 1$. Before the main results of the section can be stated, we will need some definitions from ^[29].

Definition 5.1 ([29]). Let $X$ be an arbitrary nonempty set and $\rho \colon X \to \mathbb{R}_+$. A semiflow $\Phi: \mathbb{R}_+ \times X \to X$ is called uniformly weakly $\rho$-persistent, if there exists some $\varepsilon > 0$ such that

$\Phi$ is called uniformly (strongly) $\rho$-persistent, if there exists some $\varepsilon > 0$ such that

A set $M \subseteq X$ is called weakly $\rho$-repelling if there is no $x \in X$ such that $\rho(x) > 0$ and $\Phi(t, x) \to M$ as $t \to \infty$.

Lemma 5.2. The susceptible compartments $S$, $S_D$ and $S_N$ are always uniformly persistent. More precisely, the following inequalities hold for the limit inferior of $S, S_D$ and $S_N$ respectively:

Proof. Let us suppose that $I(t) \leq 1$, $I_D(t) \leq 1$, and $I_N(t) \leq 1$. First we consider $S'(t)$:

Then,

and hence

Similarly, let us consider $S_D(t)$:

Then,

Hence,

Finally, considering $S_N(t)$, from

we obtain

Hence,

System (2.1) generates a continuous flow $\Phi$ on the feasible state space

We cite some definitions and important theorems from ^[23,29] to state the main theorem of the section.

Definition 5.3. A matrix $A \in \mathbb{R}^{n \times n}$ is called a Metzler matrix if all of its off-diagonal entries are nonnegative.

Definition 5.4. Let $\lambda_1, \dotsc, \lambda_r$ denote the eigenvalues of a matrix $A \in \mathbb{R}^{n \times n}$. The spectral abscissa $\eta(A)$ is defined as

Definition 5.5. A point $p$ is called an $\omega$-limit point of a solution $x(t, x_0)$ if there exists a nonnegative sequence $t_1, \dotsc, t_k$ of time moments such that $t_k\to \infty$ as $k \to \infty$, for which $x(t_k, x_0) \to p$, $k\to\infty$ holds. The set of all such points of $x(t_k, x_0)$ is called the $\omega$-limit (or forward) set of $x(t, x_0)$.

Definition 5.6. Assume that a semiflow $\Phi\colon \mathbb{R}_+ \times X \to X$ is continuous, and $\rho\colon X \to \mathbb{R}_+$ is continuous and not identically zero. Let us introduce the extinction space

which is closed and forward invariant.

Theorem 5.7. ([29,Theorem 8.17]). Let $\Omega \subset \bigcup_{i = 1}^{k} M_i$ where each $M_i \subset X_0$ is isolated (in $X$), compact, invariant, and weakly $\rho$-repelling, $M_i \cap M_j = \emptyset$ if $i\neq j$. If $\{M_1, \dotsc, M_k\}$ is acyclic, then $\Phi$ is uniformly weakly $\rho$-persistent.

Theorem 5.8 (Perron–Frobenius theorem for Metzler matrices, see ^[23]). Let $A \in \mathbb{R}^{n \times n}$ be a Metzler matrix. Then the spectral abscissa $\eta(A)$ of $A$ is an eigenvalue of $A$ and there exists a nonnegative eigenvector $x \geq 0$, $x \neq 0$ such that $Ax = \eta(A)x$.

Now we can proceed to state the main result of this section.

Theorem 5.9. The infected compartments $I, I_D$ and $I_N$ are uniformly persistent if $\mathcal{R}_0 > 1$.

Proof. For simplicity we introduce the notation $x = (S, I, S_D, I_D, S_N, I_N)\in X$ for the state of the system. Let us choose $\rho(x) = I+k_1I_D+k_2I_N$ where $k_1, k_2$ are positive constants to be defined later. The three-dimensional invariant extinction space of the infected compartments is

Using the usual notation $\omega(x)$ for the $\omega$-limit set of a point in $X$, and analogously to [29, Chapter 8], we examine the set $\Omega := \bigcup_{x \in X_I}^{} \omega(x)$. It is clear that all solutions started from the extinction space $X_I$ converge to the disease-free equilibrium, as follows $\Omega = \left\{E_0\right\}$.

To prove uniform weak $\rho$-persistence we apply Theorem 5.7. According to the terminology of ^[29], we let $M_1 = \left\{E_0\right\}$, then we have $\Omega \subset M_1$, where $M_1$ is isolated, compact, invariant and acyclic. We show that $M_1$ is weakly $\rho$-repelling, then by [29, Theorem 8.17], this implies the weak persistence.

Let us suppose that $M_1$ is not weakly $\rho$-repelling, i.e. there exists a solution of (2.1) which converges to the disease-free equilibrium and $\rho(x) = I+k_1I_D+k_2I_N > 0$. Then for any $\varepsilon > 0$, for $t$ sufficiently large, $S(t) > \frac{g+\lambda}{a+g+\lambda}-\varepsilon$, $S_D(t)\geq\frac{a}{a+g+\lambda}-\varepsilon$ and $S_N(t)\geq 1-\varepsilon$ hold and for $t$ sufficiently large, we can give the estimate

with a right choice of positive constants $k_1, k_2$ and $K$. So if we can find positive constants that satisfy the above equations, then, for sufficiently large $t$, we can estimate the solution from below by the solution of the following equation:

which contradicts $I(t)+I_D(t)+I_N(t) \to 0$.

Finding appropriate constants $k_1, k_2$ and $K$ means finding a positive eigenvalue $K$ with positive corresponding eigenvector $\left(1, k_1, k_2\right)$ of the matrix

Because of continuity reasons, it is enough to find a positive eigenvalue with positive corresponding eigenvector of the matrix

The matrix $M$ is a Metzler matrix, so we can use Theorem 5.8. It follows that the spectral abscissa $\eta(M)$ of the matrix $M$ is nothing else but an eigenvalue of the matrix $M$ with a corresponding nonnegative eigenvector. Let us show that $\eta(M)$ is positive if $\mathcal{R}_0 > 1$ (in this case the corresponding eigenvector is clearly positive). The characteristic polynomial of the matrix $M$ is

This is exactly $-f(x)$ (see (4.2)). We proved that $f(x)$ has at least one negative coefficient thus it follows that $-f(x)$ has at least one positive. Then, according to the Routh–Hurwitz criterion, the existence of a positive root is proved.

Hence, we have shown uniform weak persistence of $I+I_D+I_N$, from which, using [29, Theorem 4.5], uniform persistence follows. This implies that at least one of $I$, $I_D$ and $I_N$ is uniformly persistent. Applying the fluctuation lemma, using that the noninfectious compartments have shown to be uniformly persistent, we obtain the uniform persistence of the remaining two infected compartments.

Remark 5.10. The existence of at least one endemic equilibrium follows from the uniform persistence proved above. Based on this and numerical simulations, we conjecture that system (2.1) has a unique endemic equilibrium which is globally asymptotically stable on $X\setminus X_I$ if $\mathcal{R}_0 > 1$.

6.   Numerical simulations

6.1.   The effect of reducing the needle distribution rate

In this subsection we present numerical experiments to study the effect of changes in the needle distribution rate on the number of HIV infections among non-drug users and drug users. It is important to note that several parameters are hard to estimate, just like the number of drug users and HIV infections in a given population. In this numerical study, we mainly base our estimations on ^[21,25] and consider a population of secondary school students exposed to the risk of drugs and HIV infection (hence, in this case, the birth and death rate $\lambda$ stands for the rate of joining, resp. leaving the subpopulation under risk of getting addicted to intravenous drugs or getting in contact with IDUs). Table 1 shows the parameter values used in our simulations. Both the constant human and needle populations are supposed to be equal to 1.

Figure 4 shows the number of infected non-drug users, resp. IDUs for different rates of sterile needle distribution. The first figure shows the situation after the reduction of the needle exchange program, when 28 needles/IDU were distributed yearly. In this case, the reproduction number takes the value $\mathcal{R}_0 = 1.294$. In the second figure, one can see the situation before the reduction of the program, when 81 needles/IDU were provided yearly. The basic reproduction number significantly decreases, though it is still above 1, taking the value $\mathcal{R}_0 = 1.11$. The third figure shows the case of 320 needles distributed/IDU. In this case, the reproduction number decreases below 1, with $\mathcal{R}_0 = 0.9035$. Hence, these simulations suggest that a significant number of new infections could have been spared with keeping the earlier level of sterile needle distribution, moreover, a further increase of the number of new needles could have even driven the reproduction number below 1.

6.2.   Fitting of the model with time-dependent needle distribution parameter—a case study for Hungary

According to the world-wide situation, the number of HIV infected individuals and intravenous drug users is increasing. The budget of the needle exchange program has been reduced recently in Hungary. The number of reported cases and distributed needles can be seen in Figure 5. Needle distribution is extrapolated from 2018 to 2035 at the current low level as well as at the highest level for further study.

Our model given in system (2.1) cannot be directly applied to this case, since the needle distribution is changing from year to year. Hence, instead of the constant needle exchange rate $\nu$, we introduce the rate function $\nu (t)$ interpolated to the given data. This non-autonomous system is fitted to the HIV cases from 2000 to 2017. Because of the high number of parameters, many of them are taken from the literature (see above) or fixed by some assumptions. In addition, by our experience, some kind of underreporting has to be taken into account, as well.

In the study we made the following assumptions. HIV infected and intravenous drug user populations are quite closed, we assume that the exposed population is about $100 000$ (estimated: population between 11 and 65 years, living in cities, having exposed life style). The death rate in Hungary is $\lambda\approx 0.01$, we keep it, due to the higher rate in this population.

Concerning the epidemiology and further parameters: we have the following assumptions: $\lambda = 0.01, \ a = 0.06, \ g = 0.0225, \ \alpha = 0.15, \ \beta = 0.06$. Due to the uncertainties, parameters concerning the needle infections/disinfections $\theta, \omega$ are used in the fitting. The model including the underreporting ratio $q$ is simple: $\operatorname{Model}(q, \theta, \omega)(t) = q (I(t)+I_D(t))$, where $I(t), \ I_D(t)$ are solutions of (2.1).

Note that the fitting is very sensitive to the initial values. They are taken from and estimated by the available data at 2000 ^[24,25,30], and corrected by the underreporting ratio $q$: $I(0) = 0.006/q, \ I_D(0) = 0.0004/q, \ S_D(0) = 0.0026, S(0) = 1-I(0)-I_D(0) - S_D(0)$.

Taking into account that the number of collected used needles is about $60\%$ of new ones ^[30] and that the number of HIV infected and intravenous drug users are increasing, we assume that altogether $10^6$ needles are in use, and only $15\%$ of needles were infected in 2000, hence: $S_N(0) = 0.85; \ I_N(0) = 0.15$. For numerical studies, we use Wolfram Mathematica 11.3 commands ParametricNDSolve and NonlinearModelFit. The command ParametricNDSolve numerically solves differential equations with undefined parameters by creating a ParametricFunction object, and enables to work with the solutions as if they were given symbolically. Fitting was performed by the sophisticated and well-tested command NonlinearModelFit, which can be applied to implicitly defined models such as numerical solutions of differential equations provided by ParametricNDSolve, and it can measure the goodness of fitting. Default options and $\mathtt{ConfidenceLevel} \rightarrow 0.95$ were used.

The result of the fitting is presented in Table 2. Although by changing our assumptions the result would certainly be different, the main problems are well visible now. Underreporting is high, due to the nature of HIV, the ratio of hidden cases is significant, and it fits the general opinion among epidemiologists. The sterilization rate of needles ($\theta$) is low as expected, and the high rate of infection of needles ($\omega$) is in correspondence with this.

Finally, Figures 6, 7 and 8 present graphically the fitting and forecast. In Figure 6 continuous red line shows the model $I(t)+I_D(t)$ fitted to the corrected data (blue dotted line) using the real needle distribution (see Figure 5). The dashed red line shows the trend of the model (using the same parameters) if the needle distribution were at the highest level, approximately 650000 needles/year. Figure 7 shows the accumulated number of infected non-drug users ($I$) and infected drug users ($I_D$). Figure 8 shows the ratio of the two values of the compartments $I$, $I_D$, resp. $I+I_D$ obtained in the real and in the imagined model. The figures illustrate very well that the decrease of needle exchange dramatically increase the ratio of drug user infected persons $I_d$, hence it is responsible for the high value of $I$ as well. The difference in 2030 is more than significant between these two estimates, even if these results are quite hypothetic under the assumptions of our models.

7.   Discussion

In this paper, we have established a new model to study the spread of a sexually transmitted disease through sexual transmission and the use of contaminated needles among intravenous drug users. Beyond determining the dynamics of the model, our aim was to study in what extent the recent reduction of needle exchange programs in several countries indirectly affects non-drug users via sexual contacts with infected intravenous drug users. We first determined the basic reproduction number $\mathcal{R}_0$ of the system and showed that for $\mathcal{R}_0 < 1$, the unique disease-free equilibrium is globally asymptotically stable. In the case $\mathcal{R}_0 > 1$, the disease was shown to persist in the population. Using numerical simulations, we showed that the reduction of the needle exchange program can highly contribute to the spread of diseases not only among drug users, but also among those who do not use drugs, while an increased program can contribute to the eradication of these diseases. Based on real world data from Hungary, we gave a forecast for the number of HIV infections for the next decade. We also compared the number of HIV infections with the current tendency of the needle exchange program and a situation where the original high level of distribution is maintained. This simulation suggests that a significant number of new infections could be avoided with the help of the needle exchange program. Of course, the model has certain limitations. For example, for a more realistic result, one should consider several subetaoups corresponding to different tendencies of drug use, needle sharing and risky sexual behaviour. The homogeneous mixing assumption may lead to an overestimation of the impact of needle exchange programs on HIV prevalence for non-drug users. HIV being a sexually transmitted disease, a more realistic model should differentiate sexes and subetaoups with different transmission probabilities. These generalizations of our model are the topic of a future work.

Acknowledgements

This research was supported by the Ministry of Human Capacities, Hungary grant 20391-3/2018/FEKUSTRAT. A. Dénes was supported by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences, by the project no. 125119, implemented with the support provided from the National Research, Development and Innovation Fund of Hungary, financed under the SNN_17 funding scheme and by the project no. 128363, implemented with the support provided from the National Research, Development and Innovation Fund of Hungary, financed under the PD_18 funding scheme. J. Karsai was supported by the project no. 125628, implemented with the support provided from the National Research, Development and Innovation Fund of Hungary, financed under the KH_17 funding scheme and by the EU-funded Hungarian grant EFOP-3.6.2-16-2017-00015. L. Székely was supported by the EU-funded Hungarian grant EFOP-3.6.1-16-2016-00008. The authors are grateful to Dr. Gergely Röst for his useful suggestions.

Conflict of interest

The authors declare there is no conflict of interest.