In the following transmission via sexual contact is assumed. Thus, we consider a sexually active population of varying size N. AIDS patients are assumed to be too sick to be sexually active. Thus, the total population can be divided into a susceptible class of size S and an infectious class before the onset of AIDS. As the infectious period is very long (⩾10 years), it is further divided into several stages. The widely-used CDC-classification suggests 2–6 stages of infection before AIDS with the classes defined by either viral load or CD$4^{+}$T cell count, and clinical symptoms (see [5]). In our model, we use two stages according to clinical stages, that is the asymptomatic and the symptomatic phases as defined in [18]. The numbers of individuals in the infectious classes are given by $I_1$ and $I_2$, respectively. All individuals entering the population are assumed to be susceptible, because vertical transmission can be almost excluded by means of antiretroviral treatment [5]. The number of individuals entering the population per unit of time is bN with birth-rate constant b, where ‘birth’ means that an individual becomes sexually active. Similarly, d is the natural death rate constant for all subpopulations. For simplicity reasons, the total population size is assumed to grow exponentially with rate constant $b-d$ in the absence of the disease. It is further assumed that infected individuals do not change their behaviour, no matter if they are treated or not. Thus, the average number of contacts of an individual per unit of time is a constant c for all subpopulations. The disease progression rate constant from stage j to stage $j+1$ or the AIDS-stage is given by $k_j$ for $j=1\text{,}2$, respectively. To model treatment effects as described in Section 1, it is assumed that an infected individual in the second stage can move back to the first stage by means of successful drug treatment. The corresponding transfer-rate constant is denoted by α. By β and $a\beta$ we denote the probability of disease transmission during one contact with an infected individual in stages 1 and 2, respectively.

See Table 1 for an overview of the parameters and variables used in the formulation of the model equations. Here also the parameter values used in the simulations below are given. The values for the disease progression rates $k_j$ are taken from [18]. The reasons for using these quite old estimates are obvious: as today antiretroviral treatments are widely used, new estimates on these transfer rates always include treatment effects. The main goal of our work was to study the effects of treatment on the population growth and the spread of the epidemic. This is why we have to use these very old estimates for the $k_j$ that hopefully represent the spread of the epidemic in an untreated population. In case $b<d$, an epidemic of HIV infection can only speed up extinction of the population, whereas in case $b⩾d$, the disease of HIV infection can slow down the population growth or even reverse growth to decay. As the latter case is much more interesting, the simulations shown in this paper use parameter values b and d chosen such that the population grows in the absence of the disease. According to [5], the probability of disease transmission per contact is between 11% and 50%. It is highest right after infection and in the last stage before AIDS, thus we take $\beta =0.5$ and $a\beta =0.6\beta$. In the simulations we analyse two treatment levels – medium level and high level treatment with $\alpha =0.1$ and $\alpha =0.2$, respectively. It has to be pointed out that there is no data available for the size of parameter α, thus our values are only reasonable estimates for different treatment levels. Parameter c reflects the level of risky behaviour in the population.

Let $\dot{X}$ denote the derivative of X with respect to time. The model equations are then given by:

The treatment reproduction number R is used to define a threshold condition for the existence and stability of equilibrium points. It can be interpreted as the average number of infections a typical infective causes in a completely susceptible population. In [32], a formula for R was derived for a very general model with n infectious stages. The complexity of this formula increases significantly with increasing n. In the case of two infectious compartments, it yields the following expression for R:

To discuss the local stability of the endemic steady state $E^{*}$, we consider the linearized system of (1) at $E^{*}$. The Jacobian at $E^{*}$ is given by:

By using the Routh–Hurwitz criterion [35], it follows that all eigenvalues of the characteristic equation (9) have negative real parts if and only if the following condition is satisfied:

For the parameter values as given in Table 1 and $M_1\text{,}{M}_2\text{,}{M}_3$ as defined in (6)–(8), we get:

In [32], it was shown in detail that $E^{*}$ is asymptotically stable for n infectious stages whenever it exists.

Let constant τ (in months) be the time delay from the start of treatment in the symptomatic stage $I_2$ until treatment effects become visible. As with every virus infection, one has to distinguish between infected individuals having the virus in their bloodstream and infectious individuals, who also can infect others. When a patient is infected with HIV, he typically can transmit the disease not at once, but after a short timespan (in the range of a few days, see [5]) called the latency period. In the following the latency period will be modelled by a time-delay constant λ (in months).

We incorporate these two delays into our ordinary differential-equation model (1) and use the same parameters as in Table 1. We then get the following delay differential-equation model:

The characteristic equation of (16) is given by:

The characteristic equation (17) of the linearised system at the endemic steady state $E^{*}$ can be written as a third-degree polynomial equation of the following form:

Here we have (again with $M_1$, $M_2$, $M_3$ as given in (6)–(8)):

The location of the roots of Eq. (17) has been studied by many authors (see [38–46]). The following result, giving necessary and sufficient conditions for the absolute stability of (16), was proved by Hale et al. in [37,47].

To study the delay case, we consider the ordinary differential-equation model (1) first (with $\lambda =\tau =0$ for model (15)), and assume that condition (13) for the asymptotic stability of $E^{*}$ is satisfied. Thus, all the roots of (1) have negative real parts, and assumption (1) in Lemma 1 is fulfilled. Assumption (2) in Lemma 1 means that iω, $\omega >0$ is not a root of Eq. (17). Then by using Lemma 1, the delay system is asymptotically stable if and only if the ODE system is asymptotically stable and the characteristic equation (17) has no purely imaginary roots.

Let $\theta (\tau \text{,}\lambda )=\eta (\tau \text{,}\lambda )+\mathrm{i}\omega (\tau \text{,}\lambda )$, with $\omega >0$, be the eigenvalues of the characteristic equation (17), where $\eta (\tau \text{,}\lambda )$ and $\omega (\tau \text{,}\lambda )$ depend on the delays τ and λ.

If, on the other hand, the characteristic equation (17) has a pair of purely imaginary roots iω, $\omega >0$, then the delay system (15) is not absolutely stable, but could be conditionally stable. For example in the case of one delay, say for τ, and we suppose that ${\omega }_0$ is given, for some value ${\tau }_0$. Hence, when $\tau <{\tau }_0$, the real parts of all roots of the characteristic equation (17) still remain negative and the system is conditionally stable. When $\tau ={\tau }_0$, the characteristic equation (17) has a purely imaginary root iω, $\omega >0$, and the system (15) loses its stability. By using Rouché's theorem (Dieudonné [48], Theorem 9.17.4) and continuity in τ, and if the transversality condition holds ($\frac{\mathrm{d}\eta (\tau )}{\mathrm{d}\tau }{|}_{{\tau }_0}>0$) at $\tau ={\tau }_0$, then the characteristic equation (17) will have at least one root with positive real part when $\tau >{\tau }_0$ and becomes unstable.

To be able to find a condition for all eigenvalues to have negative real parts, we will determine if (17) has purely imaginary roots.

Since the positive equilibrium $E^{*}$ of the ODE model is stable, we have $\eta (0)<0$ when $\tau =\lambda =0$. By continuity, if $\tau >0$, $\lambda >0$ are sufficiently small, we still have $\eta (\tau \text{,}\lambda )<0$ and $E^{*}$ still remains stable. If $\eta ({\tau }_0\text{,}{\lambda }_0)=0$ for certain values ${\tau }_0>0$ and ${\lambda }_0>0$, the steady state $E^{*}$ loses its stability and becomes unstable when $\eta (\tau \text{,}\lambda )>0$. On the other hand, if such an $\omega ({\lambda }_0\text{,}{\tau }_0)$ does not exist for all delays, that is the characteristic equation (17) does not have purely imaginary roots, the system will always be stable.

If $\mathrm{i}\omega (\tau \text{,}\lambda )=\mathrm{i}\omega$, $\omega >0$ is a root of the characteristic equation (17), then we have:

The resulting differential-equation model including delayed onset of treatment effects but no latency period is formulated as follows:

Notice that in the case of a population with high-level treatment ($\alpha =0.2$) and all parameter values as given in Table 1, we have: $\nu \cong -7.9222\mathrm{e}-004<0$ and $\iota \cong 1.5178\mathrm{e}-009>0$. Thus, the equation has a positive root ${\omega }_0$. Hence the characteristic equation (17) with $\lambda =0$ has a purely imaginary root $\mathrm{i}\omega \cong 0.1983\mathrm{i}$.

Again, let $\theta (\tau )=\eta (\tau )+\mathrm{i}\omega (\tau )$ be the eigenvalues of Eq. (17) – with $\lambda =0$ – such that $\eta ({\tau }_0)=0\text{,}\omega ({\tau }_0)={\omega }_0\cong 0.1983$. From (21) and (22), we have

Then the real part of $\theta (\tau )$ becomes positive when $\tau >{\tau }_0\cong 8.82035$, and the steady state becomes unstable (see Fig. 4). From that, it follows that ${\tau }_0$ is a bifurcation value, and a Hopf bifurcation (see [49]) occurs when τ passes through the critical value ${\tau }_0$ in Fig. 3. We can see in Fig. 2 that the steady state $E^{*}$ still is stable for delay $\tau <{\tau }_0$, allthough the delay causes transient oscillations. Thus, we have the following theorem.

Biologically, this means that there is a critical value for the treatment-induced delay ${\tau }_0$, that determines the stability of the endemic equilibrium point $E^{*}$. With the parameter values given in Table 1, the endemic equilibrium point is asymptotically stable when antiretroviral drugs in average show positive effects in the patients within less than ${\tau }_0\cong 8.82035$ months. In this case (see Fig. 2) the population will reach the infected steady state $E^{*}$ as given in Eqs. (3)–(5). As soon as it takes more than ${\tau }_0\cong 8.82035$ months for the patients to feel better, the infected steady state $E^{*}$ loses its stability (see Fig. 4).

The underlying assumption in the present paper is that people who are infected with HIV do not change their behaviour despite the knowledge of the risks. In [50], S. Blower states that incidence rates of HIV will fall as more HIV-positive individuals gain access to treatment (HAART), but that this public health benefit will only occur if the levels of risky behaviour do not increase. Computer simulations based on a simple ordinary differential-equation model for an HIV epidemic (see [33]), on the contrary, show that treatment without reduction of risky behaviour may even increase the proportion of infected individuals. This effects can also be seen in the delay differential-equation models presented in the present paper. The simulations for the ODE model in [33] show that the combination of reduction of risky behaviour together with antiretroviral drug treatment is a very promising strategy in fighting the epidemic of HIV infection. This has to be taken into account, when anti-HIV regimens are distributed in the poorest countries of the world, where behaviour change is hard to achieve because of low educational levels or cultural reasons.

In our delay differential-equation models we do have two threshold parameters: The treatment reproduction number R serves as a threshold for the existence of an endemic equilibrium point, and the size of the treatment-induced delay τ has enormous effects on the dynamics of the spread of the disease of HIV infection in the population. This treatment-induced delay τ can be seen as the minimum time span that patients have to be under medical treatment until positive treatment effects occur. If and only if the treatment reproduction number R is above one, there is a positive steady state. If in this case, the treatment induced delay τ is below a critical value, the total sexually active population will reach the infected equilibrium point (see Fig. 2). Note that at this positive steady state there are still people that have not been infected with HIV. The proportion of these susceptibles is given by $\frac{1}{R}$. As soon as τ increases and reaches the critical value, a Hopf bifurcation occurs (see Fig. 3), and, finally, the endemic equilibrium point loses its stability above the critical value (see Fig. 4).

The present work wants to be understood as a starting point for further research on analysing treatment effects including time delay on the spread of an HIV epidemic. The next step will be to analyse the general model with n infected stages. We also have to point out that we assumed exponential growth for the total population in the absence of the disease, what is a very restrictive assumption. The simulations with logistic growth for the ODE model given in [32] show that the disease of HIV infection also could lower the carrying capacity or give rise to a second positive steady state below the carrying capacity. It remains to analyse the delay differential-equation model with other than exponential population growth and the role of treatment and behaviour change under these circumstances.