The effect of delayed death in HIV/AIDS models

ABSTRACT

HIV-infected patients who receive treatment survive for some years after they have acquired the disease. The received treatment causes sustained reduction of viral reproduction by improving the immune function, leading to prolonged progression period to AIDS development. This prolonged progression period has created variability in survival times that affects estimates produced using mathematical models that do not include delay in disease related mortality. This paper investigates the effect of including delay in AIDS death occurrence in HIV/AIDS transmission models. A simple mathematical model with two stages of HIV progression is developed and extended to include time delay in the occurrence of AIDS deaths. Numerical simulations indicate that time delay changes the mortality curves considerably but has less effect on the proportion of infectives. The study highlights the importance of incorporating delay in models of HIV/AIDS for the production of accurate HIV/AIDS estimates.

KEYWORDS:

• Mathematical model
• delay
• HIV/AIDS
• survival
• numerical simulation

1. Introduction

Developments and improvements in care and treatment services to HIV/AIDS patients, genetic heterogeneity, and socio-economic aspects of life contribute to increased differences in survival times between individuals. These differences in survival times lead to a situation where no common survival period is available to all HIV-infected individuals.

The survival times of a HIV/AIDS-infected individual depend on factors such as gender, treatment, economic status, nutrition, biological make up, and the region where the person lives. In sub-Saharan Africa, a region where the large majority of HIV is transmitted heterosexually and the epidemic is more mature, there is a substantial difference in age of infection between men and women, with women becoming infected at an earlier age [13, 20, 21]. Furthermore, an approximate of 9.4 years survival time for women and 8.6 years for men [5] has been estimated by the UNAIDS. Treatment of infected persons has also shown to delay individuals developing AIDS through reduction of the number of copies of viruses in the body [9]. This has led to an increase of survival rates of HIV positive individuals.

Mathematical models used for prediction of the spread of HIV/AIDS are constrained by the differences in survival times that exist in HIV-infected individuals. The Weibull distribution has commonly been used as an incubation period distribution of HIV survival times [14, 18]. Other forms of distribution functions have been used to explain the survival probabilities of HIV-infected individuals under which this probability decreases as one progresses from HIV infection to death [3, 5]. These distribution functions fit well to data within a short period from infection. However, due to the introduction of an effective antiretroviral therapy (ART), the incubation period defined with approximately parameterized Weibull distribution could form a different shape if survival is prolonged. Modelling work has also assumed an exponential distribution of the infectious period [1, 2, 17]. Exponential distribution functions assume that the chance of dying due to AIDS within a given time interval is constant, regardless of time since infection.

To accurately model survival, delay equations have been used by some authors in population studies and epidemiology [6, 10] and with application to HIV spread models [4, 7, 8, 12]. However, none of them have included the delay in HIV death occurrences or have provided a general and simple framework for modelling arbitrary survival times of newly infected individuals. In this paper, a mathematical model that includes the above features is formulated and simulated to investigate the impact of delayed death in HIV estimates. In particular, we investigate the effect of including time delay in the occurrence of AIDS death of HIV-infected individuals by comparing to a model without delay in HIV/AIDS death.

2. Methods

A continuous-time mathematical model of HIV/AIDS transmission that includes the effect of stages of HIV progression [11] is extended to incorporate delay in AIDS death occurrences. The model consists of one susceptible class, $Z\left(t\right)$, which forms the larger part of the population, and two infected classes: those in the primary stage, $Y_1\left(t\right)$, and those who are in the asymptomatic or late stage, $Y_2\left(t\right)$. The model assumes that AIDS-related mortality occur exponentially and that the rate of dying in the late stage of the disease is constant [11]. This property is shown by the term $\gamma Y_2\left(t\right)$ in the third equation of the model as shown below: (1) $\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }-r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}-r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Z\left(t\right),\\ \frac{\mathrm{d}{Y}_1\left(t\right)}{\mathrm{d}t}& =r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}+r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Y_1\left(t\right)-\rho Y_1\left(t\right),\\ \frac{\mathrm{d}{Y}_2\left(t\right)}{\mathrm{d}t}& =\rho Y_1\left(t\right)-\mu Y_2\left(t\right)-\gamma Y_2\left(t\right).\end{array}$(1) with initial conditions $Z\left(0\right),Y_1\left(0\right)$, and $Y_2\left(0\right)$. The total population $N\left(t\right)$ is given by $Z\left(t\right)+Y_1\left(t\right)+Y_2\left(t\right).$

As described in more detail in [11], the population dynamics of the model in Equation (1(1) $\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }-r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}-r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Z\left(t\right),\\ \frac{\mathrm{d}{Y}_1\left(t\right)}{\mathrm{d}t}& =r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}+r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Y_1\left(t\right)-\rho Y_1\left(t\right),\\ \frac{\mathrm{d}{Y}_2\left(t\right)}{\mathrm{d}t}& =\rho Y_1\left(t\right)-\mu Y_2\left(t\right)-\gamma Y_2\left(t\right).\end{array}$(1) ) are governed by the recruitment rate into the susceptible population, Λ, of new individuals into the susceptible population, exit through natural mortality at a rate μ and through the rate at which infected individuals in $Y_2\left(t\right)$ die of the disease, γ. The parameters $r_1$ and $r_2$ are rates of transmission when the interaction between the susceptible and infected individuals, $Y_1\left(t\right)$ and $Y_2\left(t\right)$ occurs. $r_1>r_2$ always [16], and ρ is the progression rate from $Y_1\left(t\right)$ to $Y_2\left(t\right)$. In the absence of the disease, the model assumes that the total population grows or decreases exponentially.

2.1. Modelling delay in HIV/AIDS deaths

To extend model (1(1) $\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }-r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}-r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Z\left(t\right),\\ \frac{\mathrm{d}{Y}_1\left(t\right)}{\mathrm{d}t}& =r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}+r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Y_1\left(t\right)-\rho Y_1\left(t\right),\\ \frac{\mathrm{d}{Y}_2\left(t\right)}{\mathrm{d}t}& =\rho Y_1\left(t\right)-\mu Y_2\left(t\right)-\gamma Y_2\left(t\right).\end{array}$(1) ) to include delay in AIDS-related death, we assume that mortality of an infected individual is dependent on the survival period, τ. We assume that individuals in the asymptomatic stage of HIV infection, $Y_2\left(t\right)$ survive for a period τ years before they die. In this case, we modify the death term $\gamma Y_2\left(t\right)$ in the third equation of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }-r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}-r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Z\left(t\right),\\ \frac{\mathrm{d}{Y}_1\left(t\right)}{\mathrm{d}t}& =r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}+r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Y_1\left(t\right)-\rho Y_1\left(t\right),\\ \frac{\mathrm{d}{Y}_2\left(t\right)}{\mathrm{d}t}& =\rho Y_1\left(t\right)-\mu Y_2\left(t\right)-\gamma Y_2\left(t\right).\end{array}$(1) ) to become $\gamma Y_2(t-\tau )$.

The last equation in model (1(1) $\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }-r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}-r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Z\left(t\right),\\ \frac{\mathrm{d}{Y}_1\left(t\right)}{\mathrm{d}t}& =r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}+r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Y_1\left(t\right)-\rho Y_1\left(t\right),\\ \frac{\mathrm{d}{Y}_2\left(t\right)}{\mathrm{d}t}& =\rho Y_1\left(t\right)-\mu Y_2\left(t\right)-\gamma Y_2\left(t\right).\end{array}$(1) ) becomes (2) $\frac{\mathrm{d}{Y}_2\left(t\right)}{\mathrm{d}t}=\rho Y_1\left(t\right)-\mu Y_2\left(t\right)-\gamma Y_2(t-\tau ),$(2) which describes the dynamics of the group of individuals in the asymptomatic stage. The dynamics of the disease in the population with a time delay in the AIDS-related mortality is given by the following system of nonlinear differential equations: (3) $\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }-r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}-r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Z\left(t\right),\\ \frac{\mathrm{d}{Y}_1\left(t\right)}{\mathrm{d}t}& =r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}+r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Y_1\left(t\right)-\rho Y_1\left(t\right),\\ \frac{\mathrm{d}{Y}_2\left(t\right)}{\mathrm{d}t}& =\rho Y_1\left(t\right)-\mu Y_2\left(t\right)-\gamma Y_2(t-\tau ).\end{array}$(3) All parameters are positive. The dynamics of the population change according to (4) $\frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}=\mathrm{\Lambda }-\mu N\left(t\right)-\gamma Y_2(t-\tau ).$(4)

All parameters are positive (Table 1). Since the total population is variable, we define $z\left(t\right)=Z\left(t\right)/N\left(t\right)$, $y_1\left(t\right)=Y_1\left(t\right)/N\left(t\right)$, and $y_2\left(t\right)=Y_2\left(t\right)/N\left(t\right)$. Converting the system in (3(3) $\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }-r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}-r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Z\left(t\right),\\ \frac{\mathrm{d}{Y}_1\left(t\right)}{\mathrm{d}t}& =r_1\frac{Z\left(t\right)Y_1\left(t\right)}{N\left(t\right)}+r_2\frac{Z\left(t\right)Y_2\left(t\right)}{N\left(t\right)}-\mu Y_1\left(t\right)-\rho Y_1\left(t\right),\\ \frac{\mathrm{d}{Y}_2\left(t\right)}{\mathrm{d}t}& =\rho Y_1\left(t\right)-\mu Y_2\left(t\right)-\gamma Y_2(t-\tau ).\end{array}$(3) ) to proportions gives (5) $\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }(1-z(t\left)\right)\zeta \left(t\right)-r_{1}z\left(t\right)y_1\left(t\right)-r_{2}z\left(t\right)y_2\left(t\right)+\gamma z\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_1\left(t\right)}{\mathrm{d}t}& =r_{1}z\left(t\right)y_1\left(t\right)+r_{2}z\left(t\right)y_2\left(t\right)-(\rho -\mathrm{\Lambda })y_1\left(t\right)\zeta \left(t\right)+\gamma y_1\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_2\left(t\right)}{\mathrm{d}t}& =\rho y_1\left(t\right)-\mathrm{\Lambda }{y}_2\left(t\right)-\gamma (1-y_2(t\left)\right)y_2(t-\tau )\alpha \left(t\right),\end{array}$(5) where (6) $y_2(t-\tau )=\frac{Y_2(t-\tau )}{N(t-\tau )}$(6) (7) $\alpha \left(t\right)=\frac{N(t-\tau )}{N\left(t\right)},$(7) and (8) $\zeta \left(t\right)=\frac{1}{N\left(t\right)}.$(8)

Table 1. Description and values of parameters of the model.

Parameter	Description	Units	Baseline	Source
Λ	Recruitment rate of individuals into the susceptible population	people year^−1	200	Estimate
$r_1$	Transmission rate in the primary stage	year^−1	4.0	Estimate
$r_2$	Transmission rate in the asymptomatic stage	year^−1	$r_1/12$	[16, 19, 22]
μ	Natural mortality rate	year^−1	$0.02$	[24]
ρ	Progression rate from primary to asymptomatic stage	year^−1	$6.0$	[15, 23]
γ	Mortality rate due to disease	year^−1	$0.1$, variable	Estimate
τ	Survival period of infected individuals	year	10, variable	Estimate

Equation (4(4) $\frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}=\mathrm{\Lambda }-\mu N\left(t\right)-\gamma Y_2(t-\tau ).$(4) ) also becomes (9) $\frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}=\mathrm{\Lambda }-\{\mu -\gamma y_2(t-\tau \left)\alpha \right(t\left)\right\}N\left(t\right)$(9) integrating to (10) $N\left(t\right)=N_0\exp \left({\int }_0^t(\mathrm{\Lambda }-\{\mu -\gamma y_2(s-\tau )\alpha \left(s\right)\left\}\right)\mathrm{d}s\right).$(10)

Model (5(5) $\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }(1-z(t\left)\right)\zeta \left(t\right)-r_{1}z\left(t\right)y_1\left(t\right)-r_{2}z\left(t\right)y_2\left(t\right)+\gamma z\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_1\left(t\right)}{\mathrm{d}t}& =r_{1}z\left(t\right)y_1\left(t\right)+r_{2}z\left(t\right)y_2\left(t\right)-(\rho -\mathrm{\Lambda })y_1\left(t\right)\zeta \left(t\right)+\gamma y_1\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_2\left(t\right)}{\mathrm{d}t}& =\rho y_1\left(t\right)-\mathrm{\Lambda }{y}_2\left(t\right)-\gamma (1-y_2(t\left)\right)y_2(t-\tau )\alpha \left(t\right),\end{array}$(5) ) is epidemiologically and mathematically well posed in the region ${\mathrm{\Omega }}_2$ defined by (11) ${\mathrm{\Omega }}_2=\left\{\left(\begin{array}{c}z\\ y_1\\ y_2\end{array}\right)\in {\mathbb{R}}^3\right.\left|\begin{array}{c}z\ge 0,\\ y_1\ge 0,\\ y_2\ge 0,\\ z+y_1+y_2\le 1\end{array}\right\}.$(11)

Theorem 2.1.

If the initial conditions of system (5(5) $\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }(1-z(t\left)\right)\zeta \left(t\right)-r_{1}z\left(t\right)y_1\left(t\right)-r_{2}z\left(t\right)y_2\left(t\right)+\gamma z\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_1\left(t\right)}{\mathrm{d}t}& =r_{1}z\left(t\right)y_1\left(t\right)+r_{2}z\left(t\right)y_2\left(t\right)-(\rho -\mathrm{\Lambda })y_1\left(t\right)\zeta \left(t\right)+\gamma y_1\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_2\left(t\right)}{\mathrm{d}t}& =\rho y_1\left(t\right)-\mathrm{\Lambda }{y}_2\left(t\right)-\gamma (1-y_2(t\left)\right)y_2(t-\tau )\alpha \left(t\right),\end{array}$(5) ) lie in region ${\mathrm{\Omega }}_2,$ then (5(5) $\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }(1-z(t\left)\right)\zeta \left(t\right)-r_{1}z\left(t\right)y_1\left(t\right)-r_{2}z\left(t\right)y_2\left(t\right)+\gamma z\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_1\left(t\right)}{\mathrm{d}t}& =r_{1}z\left(t\right)y_1\left(t\right)+r_{2}z\left(t\right)y_2\left(t\right)-(\rho -\mathrm{\Lambda })y_1\left(t\right)\zeta \left(t\right)+\gamma y_1\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_2\left(t\right)}{\mathrm{d}t}& =\rho y_1\left(t\right)-\mathrm{\Lambda }{y}_2\left(t\right)-\gamma (1-y_2(t\left)\right)y_2(t-\tau )\alpha \left(t\right),\end{array}$(5) ) has a unique solution $z\left(t\right),$ $y_1\left(t\right),$ and $y_2\left(t\right)$ which remains in ${\mathrm{\Omega }}_2$ for all time $t\ge 0$.

Proof.

The right-hand side of the system (5(5) $\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }(1-z(t\left)\right)\zeta \left(t\right)-r_{1}z\left(t\right)y_1\left(t\right)-r_{2}z\left(t\right)y_2\left(t\right)+\gamma z\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_1\left(t\right)}{\mathrm{d}t}& =r_{1}z\left(t\right)y_1\left(t\right)+r_{2}z\left(t\right)y_2\left(t\right)-(\rho -\mathrm{\Lambda })y_1\left(t\right)\zeta \left(t\right)+\gamma y_1\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_2\left(t\right)}{\mathrm{d}t}& =\rho y_1\left(t\right)-\mathrm{\Lambda }{y}_2\left(t\right)-\gamma (1-y_2(t\left)\right)y_2(t-\tau )\alpha \left(t\right),\end{array}$(5) ) is continuous with continuous partial derivatives in ${\mathrm{\Omega }}_2$; therefore (5(5) $\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }(1-z(t\left)\right)\zeta \left(t\right)-r_{1}z\left(t\right)y_1\left(t\right)-r_{2}z\left(t\right)y_2\left(t\right)+\gamma z\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_1\left(t\right)}{\mathrm{d}t}& =r_{1}z\left(t\right)y_1\left(t\right)+r_{2}z\left(t\right)y_2\left(t\right)-(\rho -\mathrm{\Lambda })y_1\left(t\right)\zeta \left(t\right)+\gamma y_1\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_2\left(t\right)}{\mathrm{d}t}& =\rho y_1\left(t\right)-\mathrm{\Lambda }{y}_2\left(t\right)-\gamma (1-y_2(t\left)\right)y_2(t-\tau )\alpha \left(t\right),\end{array}$(5) ) has a unique solution that exist for all time. It remains to show that ${\mathrm{\Omega }}_2$ is forward-invariant. We see from system (5(5) $\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }(1-z(t\left)\right)\zeta \left(t\right)-r_{1}z\left(t\right)y_1\left(t\right)-r_{2}z\left(t\right)y_2\left(t\right)+\gamma z\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_1\left(t\right)}{\mathrm{d}t}& =r_{1}z\left(t\right)y_1\left(t\right)+r_{2}z\left(t\right)y_2\left(t\right)-(\rho -\mathrm{\Lambda })y_1\left(t\right)\zeta \left(t\right)+\gamma y_1\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_2\left(t\right)}{\mathrm{d}t}& =\rho y_1\left(t\right)-\mathrm{\Lambda }{y}_2\left(t\right)-\gamma (1-y_2(t\left)\right)y_2(t-\tau )\alpha \left(t\right),\end{array}$(5) ) that if z=0, then $z^{\mathrm{\prime }}=\mathrm{\Lambda }\zeta \left(t\right)\ge 0$; if $y_1=0$, then $y_1^{\mathrm{\prime }}=r_{2}z\left(t\right)y_2\left(t\right)\ge 0$; if $y_2=0$, then $y_2^{\mathrm{\prime }}=\rho y_1\left(t\right)-\gamma y_2(t-\tau )\alpha \left(t\right)\ge 0$. Therefore, all solutions of the system of equations (Equation (5(5) $\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }(1-z(t\left)\right)\zeta \left(t\right)-r_{1}z\left(t\right)y_1\left(t\right)-r_{2}z\left(t\right)y_2\left(t\right)+\gamma z\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_1\left(t\right)}{\mathrm{d}t}& =r_{1}z\left(t\right)y_1\left(t\right)+r_{2}z\left(t\right)y_2\left(t\right)-(\rho -\mathrm{\Lambda })y_1\left(t\right)\zeta \left(t\right)+\gamma y_1\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_2\left(t\right)}{\mathrm{d}t}& =\rho y_1\left(t\right)-\mathrm{\Lambda }{y}_2\left(t\right)-\gamma (1-y_2(t\left)\right)y_2(t-\tau )\alpha \left(t\right),\end{array}$(5) )) are contained in the region ${\mathrm{\Omega }}_2$.

Model (5(5) $\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }(1-z(t\left)\right)\zeta \left(t\right)-r_{1}z\left(t\right)y_1\left(t\right)-r_{2}z\left(t\right)y_2\left(t\right)+\gamma z\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_1\left(t\right)}{\mathrm{d}t}& =r_{1}z\left(t\right)y_1\left(t\right)+r_{2}z\left(t\right)y_2\left(t\right)-(\rho -\mathrm{\Lambda })y_1\left(t\right)\zeta \left(t\right)+\gamma y_1\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_2\left(t\right)}{\mathrm{d}t}& =\rho y_1\left(t\right)-\mathrm{\Lambda }{y}_2\left(t\right)-\gamma (1-y_2(t\left)\right)y_2(t-\tau )\alpha \left(t\right),\end{array}$(5) ) is quite complex for analysis. We calculate the equilibrium points and carry out analysis of certain behaviour numerically.

2.2. Equilibrium points

Definition 2.1.

If as $t\to \mathrm{\infty }$ the system in equation (5(5) $\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }(1-z(t\left)\right)\zeta \left(t\right)-r_{1}z\left(t\right)y_1\left(t\right)-r_{2}z\left(t\right)y_2\left(t\right)+\gamma z\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_1\left(t\right)}{\mathrm{d}t}& =r_{1}z\left(t\right)y_1\left(t\right)+r_{2}z\left(t\right)y_2\left(t\right)-(\rho -\mathrm{\Lambda })y_1\left(t\right)\zeta \left(t\right)+\gamma y_1\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_2\left(t\right)}{\mathrm{d}t}& =\rho y_1\left(t\right)-\mathrm{\Lambda }{y}_2\left(t\right)-\gamma (1-y_2(t\left)\right)y_2(t-\tau )\alpha \left(t\right),\end{array}$(5) ) attains an equilibrium state, then we can define $y_1\left(t\right)\to y_1^{\ast },$ $y_2\left(t\right)\to y_2^{\ast },$ $y_2(t-\tau )\to y_2^{\ast }$ and $\alpha \left(t\right)\to {\alpha }^{\ast }$. Since $N\left(t\right)\to N^{\ast }$ and $N(t-\tau )\to N^{\ast }$ at equilibrium, then, by Equation (7(6) $y_2(t-\tau )=\frac{Y_2(t-\tau )}{N(t-\tau )}$(6) ), ${\alpha }^{\ast }\to 1.$

Setting the system in Equation (5(5) $\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }(1-z(t\left)\right)\zeta \left(t\right)-r_{1}z\left(t\right)y_1\left(t\right)-r_{2}z\left(t\right)y_2\left(t\right)+\gamma z\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_1\left(t\right)}{\mathrm{d}t}& =r_{1}z\left(t\right)y_1\left(t\right)+r_{2}z\left(t\right)y_2\left(t\right)-(\rho -\mathrm{\Lambda })y_1\left(t\right)\zeta \left(t\right)+\gamma y_1\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_2\left(t\right)}{\mathrm{d}t}& =\rho y_1\left(t\right)-\mathrm{\Lambda }{y}_2\left(t\right)-\gamma (1-y_2(t\left)\right)y_2(t-\tau )\alpha \left(t\right),\end{array}$(5) ) to zero at the equilibrium and applying Definition 2.1, the disease-free equilibrium $(1,0,0)$ is found.

From the third equation of (5(5) $\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }(1-z(t\left)\right)\zeta \left(t\right)-r_{1}z\left(t\right)y_1\left(t\right)-r_{2}z\left(t\right)y_2\left(t\right)+\gamma z\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_1\left(t\right)}{\mathrm{d}t}& =r_{1}z\left(t\right)y_1\left(t\right)+r_{2}z\left(t\right)y_2\left(t\right)-(\rho -\mathrm{\Lambda })y_1\left(t\right)\zeta \left(t\right)+\gamma y_1\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_2\left(t\right)}{\mathrm{d}t}& =\rho y_1\left(t\right)-\mathrm{\Lambda }{y}_2\left(t\right)-\gamma (1-y_2(t\left)\right)y_2(t-\tau )\alpha \left(t\right),\end{array}$(5) ), and by Definition 2.1, $y_1^{\ast }$ is given (12) $y_1^{\ast }=\frac{(\mathrm{\Lambda }-\gamma )y_2^{\ast }-\gamma y_2^{\ast 2}}{\rho }.$(12)

Substituting Equation  (12(12) $y_1^{\ast }=\frac{(\mathrm{\Lambda }-\gamma )y_2^{\ast }-\gamma y_2^{\ast 2}}{\rho }.$(12) ) in the first equation of (5(5) $\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }(1-z(t\left)\right)\zeta \left(t\right)-r_{1}z\left(t\right)y_1\left(t\right)-r_{2}z\left(t\right)y_2\left(t\right)+\gamma z\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_1\left(t\right)}{\mathrm{d}t}& =r_{1}z\left(t\right)y_1\left(t\right)+r_{2}z\left(t\right)y_2\left(t\right)-(\rho -\mathrm{\Lambda })y_1\left(t\right)\zeta \left(t\right)+\gamma y_1\left(t\right)y_2(t-\tau )\alpha \left(t\right),\\ \frac{\mathrm{d}{y}_2\left(t\right)}{\mathrm{d}t}& =\rho y_1\left(t\right)-\mathrm{\Lambda }{y}_2\left(t\right)-\gamma (1-y_2(t\left)\right)y_2(t-\tau )\alpha \left(t\right),\end{array}$(5) ) gives (13) $z^{\ast }=\frac{\mathrm{\Lambda }\rho {\zeta }^{\ast }}{(\mathrm{\Lambda }\rho {\zeta }^{\ast }+r_1(\mathrm{\Lambda }-\gamma )y_2^{\ast }-\gamma y_2^{\ast 2})+r_2\rho y_2^{\ast }-\gamma \rho y_2^{\ast }}$(13) and $y_2^{\ast }=1-z^{\ast }-y_1^{\ast }$ which approximates the endemic equilibrium point.

3. Model parameterization

We parameterize the model with data from the literature. Because of the variations in data from studies, a single value was chosen from a range of values and used for the numerical simulation of the model.

Perelson and Nelson [15] and Witten and Perelson [23] have shown that an infected individual stays 10–15 years in the asymptomatic stage of HIV progression before developing AIDS. In the AIDS stage, an infected individual stays for about 1–2 years [22]. In this regard,the survival period of a HIV-infected person can extent to more than 10 years.

Understanding how infectious a person can be when infected and estimating the rate at which this person is able to transmit HIV to others has been difficult in the scientific community since one cannot perform experiments. Rapatski et al. [19] have shown using mathematical models (finding the best fit to data) that transmission rate of the virus by stage differs when an infected individual progresses from one stage to another. The findings from the model showed that infected persons in the primary stage are 12 times more likely to infect susceptible than those in the asymptomatic stage. On the other hand, in communities such as sub-Saharan Africa where heterosexual transmission is the main mode of HIV transmission, a Rakai study by Wawer et al. [22] and Quinn et al. [16] presented an analysis that provided the first empirical data on the substantial variation in transmission by stage of HIV infection after seroconversion. These studies also showed that the rate of HIV transmission within the first two-and-a-half months was almost 12 times higher than that observed in chronic couples. This presents observed evidence that the HIV transmission rate for those individuals in the primary stage of HIV progression ($r_1$) is higher than that of individuals in the asymptomatic stage ($r_2$).

The rate of progression from $Y_1\left(t\right)$ to $Y_2\left(t\right)$, ρ, is assumed to be constant and is derived from time spent in the primary stage by a recently infected individual ($1/\rho$). Perelson and Nelson [15] and Witten and Perelson [23] described the waiting times in the first stage of HIV progression to be between 2 and 10 weeks. Asymptomatic individuals stay in $Y_2\left(t\right)$ for 10–15 years [15, 23].

In sub-Saharan Africa, majority of the young adults are expected to live an average of 50 years [24], which makes the mortality rate, $\mu =0.02$ years⁻¹. We also assume in the simulations that 200 individuals are recruited into the population each year. We summarize these parameter values in Table 1.

4. Results

The model that incorporated delay in the occurrence of AIDS mortality was simulated to numerically study the differences in survival times. Results from this model were compared to numerical results of the model that did not consider delay to gain insights into the differences between them.

Figure 1 presents results for different survival periods. Results show oscillatory behaviour in all numerical results generated using the model with delay in AIDS mortality. However, this oscillatory behaviour decreases with increasing time as the equilibrium point is approached. In addition, the survival period has an impact on the duration to reach the equilibrium points. The longer the infected individual survives the asymptomatic stage, the longer it takes for the equilibrium point to be attained. The amplitude of the oscillations increases with increasing survival period. This behaviour was evident in all sub-plots.

Figure 1. Effect of the survival periods, τ, for the model with delay in death occurrence. (a) Susceptible people (b) Proportion of infected people in the primary stage, $y_1$ (c) Proportion of infected people in the asymptomatic stage, $y_2$ (d) AIDS mortality rates and (e) Total population. Parameters: $r_1=4.0$, $r_2=r_1/12$, $\mu =0.02$, $\mathrm{\Lambda }=200$, $\rho =6.0$ year⁻¹, and $\gamma =0.1$ year⁻¹.

Given that the transmission rate was the same for the different scenarios of the survival periods, the rise of the curves in $Y_1$, $Y_2$, and the total population does not differ for the different values of τ. Similarly, the proportion of the susceptible population decreased at the same rate for all the different values of the survival times.

Results show that, as the survival period is increased, predictions show roughly similar results within a short period of about t=20 years measured from the start of the epidemic for the proportion of infected individuals who are in the primary stage, asymptomatic stage, and the total population. During this period, deaths due to HIV/AIDS do not occur and people still survive although they are infected. Differences arise in the proportion of infected individuals who are in the primary stage, asymptomatic stage, and the total population after t=20 years. As τ is increased, the proportion of infected individuals in the primary stage decreases, a dynamic behaviour that is different from that of those in the asymptomatic stage.

As expected, after the initial rise, in the proportion of infected population in the asymptomatic stage, more individuals accumulate when the survival period in that particular group is high. Similarly, the mortality curve indicates that, although there is a delay in the occurrence of deaths, more people die when τ is high. Similar dynamics were also observed for different values of γ in the model without delay (Figure 2).

Figure 2. Effect of the death rate, γ, for the model without delay in the occurrence of deaths. (a) Susceptible people (b) Proportion of infected people in the primary stage, $y_1$ (c) Proportion of infected people in the asymptomatic stage, $y_2$ (d) AIDS mortality rates and (e) Total population. Parameters: $r_1=4.0$, $r_2=r_1/12$, $\mu =0.02$, $\mathrm{\Lambda }=200$, and $\rho =6.0$ year⁻¹.

The effect of a delay in AIDS mortality is assessed by comparing the results of the model with and without a delay (Figure 3). The model with delay shows oscillatory behaviour than the model without delay. This is due to the time lag introduced in the $Y_2$, a sub-population that accumulates infected individuals. The initial rise of infected individuals in the primary and asymptomatic stage is not the same for the two models. The rise is sharp for the model with delay compared to that without delay. Unlike infected groups, mortality is initially high for the model without delay. At later times, disease-related mortality is high, though with an oscillatory behaviour. Although the model with delay reaches the equilibrium point later when compared to the model without delay, numerical results in Figure 3 show that the equilibrium points for both models is the same, confirming that Definition 2.1 is correct.

Figure 3. A comparison between the model with and without delay. (a) Susceptible people (b) Proportion of infected people in the primary stage, $y_1$ (c) Proportion of infected people in the asymptomatic stage, $y_2$ (d) AIDS mortality rates and (e) Total population. Parameters: $r_1=4.0$, $r_2=r_1/12$, $\mu =0.02$, $\mathrm{\Lambda }=200$, $\rho =6.0$year⁻¹, $\gamma =0.1$, and $\tau =10$.

5. Discussion and concluding remarks

This work has attempted to describe the dynamics of HIV in a population and the role played by introducing a time delay in the occurrence of mortality. A mathematical model that divides the infected group into two groups – primary and asymptomatic stages – were formulated. This model was extended to include time delay in HIV/AIDS-related mortality. In the extended model, we assumed that AIDS deaths have a profound effect on both the spread of the disease and predictions. We altered the mechanism of AIDS deaths occurrence and consider an explicit constant time delay for all infected individuals from their initial infection to death.

We have proved that both models are mathematically and epidemiologically well-posed. Furthermore, we have analyzed the existence of disease-free equilibrium points and proved the existence of endemic equilibrium points.

The results from this work indicate that models of death strongly determine disease-related mortality. AIDS mortality for the model with delay was found to be far higher than that in the model without delay. These results suggest that, if one is to do projection of AIDS mortality using mathematical models, then, including a delay in the occurrence of mortality would make a significant impact on the results. The proportion of new infections and the prevalence of the disease are mainly determined by transmission rates in the early stages of the epidemic rather than the extension of the survival period.

The models with and without delay in survival periods did not show large differences in prevalence and in the proportion of infected individuals in the primary stage. However, small oscillatory behaviour occurred in the solutions. The hypothesis that increasing HIV-infected individuals' survival period increases the proportion of infected individuals in the population and facilitates continued spread of HIV in the population was also investigated. Findings have shown that HIV survival has a low impact on the spread of HIV. Increasing survival time of HIV-infected individuals may increase prevalence and new infections but decrease new infections and mortality in a long term. This happens due to the accumulation of infected individuals in the asymptomatic group. Similar findings are also reported by Mukandavire et al. [12], who observed that prolonged incubation periods of HIV/AIDS enlarges the epidemic in the model delay.

Results have been based on estimated parameter values which may however not apply in some cases. A more realistic model would allow for the possibility that parameter values are not constant over time, and are different in different regions. The assumption that the total population is not constant has led to increased nonlinearity of the terms in the proportion forms of the models. This study is also limited by having few parameters in the models. For example, instead of studying the effect of changing survival periods of infected individuals in a general case, it could be interesting to include parameters that describe the effects due to treatment programmes directly. The time delay introduced in the model assumes that people die following a step function distribution which is an extreme case.

In general, this work provides an understanding of how difficult it can be to produce quality HIV/AIDS estimates using mathematical models. However, a model that incorporates delay in mortality would be an added advantage as the coverage of care and treatment programs is now on the rise.

Acknowledgments

This work was developed from my MSc dissertation submitted to the University of Stellenbosch, South Africa. The author acknowledges the helpful discussions she had with Fritz Hahne and grateful to the anonymous reviewers for their reviews that improved the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

Part of this work was partially supported by a bursary from the African Institute for Mathematical Sciences.