INTRODUCTION
The availability of new antiretroviral drugs and highly active antiretroviral therapy (HAART) has led to a significant reduction in HIV morbidity and mortality^{1}. The pharmacological activity of HAART has inhibitory effects on human immunodeficiency virus (HIV) duplication and has shown a significant reduction in acquired immunodeficiency syndrome (AIDS) epidemics as well as deaths^{2}. Inhibitory of HIV duplication is generally associated with a steady increase in the CD4 cell count and results in improved clinical outcomes^{3}. This increases the chances of survival for HIV infected patients.
Sieleunou et al.^{3} used a retrospective cohort study of HIV patients older than 15 years in a rural centre in the FarNorth Province of Cameroon to explore determinants of survival of HIV patients on antiretroviral therapy (ART). They used KaplanMeier analysis to estimate survival and Cox proportional hazard models to explain survival. Their findings show that CD4 cell count, gender and clinical stage at enrolment are the main predictors of mortality. A similar study was conducted by Zhang et al.^{4} in Liangshan Prefecture, Southwest of China using the same methods. Zhang et al.^{4} observed that male patients on ART were at a higher risk of death from AIDS than their female counterparts and that a CD4 baseline cell count below 350 cells/mm^{3} results in a higher risk of death for the patients than those with a CD4 cell count of at least 350.
Seyoum et al.^{5} carried out a retrospective cohort study with collected data from clinical records of adult HIV patients following (20062010) antiretroviral therapy (ART) in Southwestern Ethiopia. Their findings reveal the main factors associated with mortality as baseline age (below 35 years) and low adherence to ART.
This study is carried out to explore the survival of HIV patients on combination antiretroviral therapy (cART) using retrospective data from a Wellness clinic in the northern part of South Africa. A multistate HIV, AIDS, Death (HAD) model is developed and used to compare survival and mortality rates of patients receiving cART. Factors associated with survival or mortality of patients are analyzed using the continuoustime Markov modelling approach. Some of the factors explored include age at baseline, CD4 cell count at baseline, gender, peripheral neuropathy (PN), nonadherence (NA), and an orthogonal viral load residual covariate (RV).
In the next section, the methods used for the analysis of data are explained. Descriptive statistics for the data that is used for analysis are given. This is followed by a presentation of results from the analysis. The last section discusses and concludes the findings.
MATERIALS AND METHODS
Ethical considerations: The procedures used in this study were as approved by the Research Ethics Committee of the University of Venda, South Africa (Protocol number SMNS/13/ MBY/01/0625), following the 1964 Helsinki declaration and its subsequent amendments. Additionally, permission to access health facilities was obtained from the Limpopo Provincial Department of Health, South Africa and the collaborating health facilities. Informed consent was obtained from study participants before their involvement and data obtained were stripped of personal identifiers to ensure the anonymity and confidentiality of the participants.
Data: A retrospective cohort study was carried out at a HIV Wellness Clinic Wellness clinic in the Limpopo Province of South Africa on 320 patients living with HIV/AIDS who had been attending ART followup care from 2005 to 2009. At treatment commencement (time t = 0), there were 224 females and 96 males. About 50 and 65% of the female and male deaths, respectively occurred during the first 6 months of treatment uptake. The interquartile range of patient ages is (33, 48) years with a mean and median age of 40.62 and 41 years, respectively. The ages were negatively skewed (skew = 0.08) confirming that there were more younger patients than older patients in the cohort.
Timehomogeneous Markov jump model: A Markov jump model on a finite or countable set, S, is a family of random variables (Xc(t))_{t}_{≥0} ((right continuous), on a probability space (Ω, Fx_{c}(s), P). Fx_{c}(s), denote all the information pertaining to the history of X_{c} up to s<t and c represents the number of states. According to the continuoustime homogeneous Markov jump process assumptions, individuals transition independently among states. This means we can assume that the transition intensities are constant over time, that is, the transition intensities are independent of t^{6,7}. Thus, for the timehomogeneous Markov jump model we have:
where, Q is a c×c transition rate matrix and becomes:
For some constant matrix Q. This implies that sojourn time with a particular state, i, has an exponential distribution with rate parameter:
where, q_{ij} is the (i, j)th entry of Q. Thus, transition probabilities only depend on the interval between times t_{1} and t_{2} and not on t_{1} itself.
For a continuoustime homogeneous model, the transition probabilities satisfy the ChapmanKolmogorov’s equation^{8}:
According to Longini and Hudgens^{9}:
leads to the Kolmogorov’s forward equation:
The Kolmogorov’s forward differential equation is derived as:
Rearranging to derive the forward differential equation gives:
Taking the limits as Δt→0 gives the desired result:
This is the Kolmogorov forward equation for the process. In the biology literature, this system of equations is termed the chemical master equation^{10}.
The Kolmogorov backward equation is derived from ChapmanKolmogorov’s equation by substituting s = Δt as follows:
Now, since:
and:
We have:
If we then take P_{ij} (t) term to the lefthand side, divide by Δt and then taking limits.
As Δt→0, we obtain the differential equation:
Model formulation: Consider the following HIV, AIDS, Death (HAD) model with the given transition rates. The state space is S = {H, A, D}. These states are based on CD4 cell counts as follows:

Fig. 1: 
HIV, AIDS and Death (HAD) model for HIVinfected patients on combination antiretroviral therapy 
The CD4 state H represents CD4 cell counts above 200 cells mm^{–}^{3}, state A represents the AIDSdefining state, which is characterised by CD4 cell counts below 200 cells mm^{–}^{3} and state D is the DEATH state that can be reached from either state H or state A. At each stage, an individual is expected to be in either state A, state H or state D. The states H, A and D are mutually exclusive. States H, A and D are defined for patients receiving antiretroviral therapy, such that the transitions between states are bidirectional due to adherence or nonadherence to treatment as shown in Fig. 1.
The transition rate from the AIDS state to the DEATH state is denoted by ν. Life may be in the HIV state or the AIDS state on several separate occasions before making the oneway transition to the death state. Alternatively, life may pass from the HIV state to the DEATH state without ever having been in the AIDS state.
Using the notations of the HAD model, the expression for q_{AH}, q_{HH}, q_{HA}, q_{HD}, q_{AD}, q_{AA} and q_{DD }can be given as follows:
The generator matrix for the HAD model is:
where, the order of the states has been taken to be H, A then D.
Kolmogorov’s forward equation: For the HAD model, the differential equation for P_{HH} (t) can be defined by the general forward equation as a template. This gives:
Now substituting it for transition rates, we have:
The Kolmogorov’s forward equation for the transition probability P_{HA} (t) is:
Kolmogorov’s backward equation: For the HAD model, the Kolmogorov’s backward differential equation for P_{HH} (t) can be obtained using the general backward equation as a template. This gives:
Now substituting in the transition rates, we have:
The backward equation for the transition probability P
_{HA} (t) is given by:
Maximum likelihood estimators: From the HAD model defined:
Let:
T_{hi} 
= 
Waiting time of the ith life in the HIV state 
T_{ai} 
= 
Waiting time of the ith life in the AIDS state 
S_{i} 
= 
Number of transitions HIV→AIDS by the ith life 
R_{i} 
= 
Number of transitions AIDS→HIV by the ith life 
D_{i} 
= 
Number of transitions HIV→Death by the ith life 
U_{i} 
= 
Number of transitions AIDS→Death by the ith life 
We also need to define totals:
Using the lower case symbols for the observed samples, it can be shown that the likelihood for the parameters, μ, ν, σ, ρ for the HAD model is given by:
Where:
The likelihood function L (μ, ν, σ, ρ) for the ith life reflects:
• 
Probability of the life remaining in the HIV state for total time T_{Hi} and in the AIDS state for time T_{Ai}, giving the factors e^{(μ+σ)THi} and e^{(ν+ρ)TAi} respectively 
• 
Probability of life making the relevant number of transitions between states giving the factors μ^{di}, ν^{ui}, σ^{si}, and ρ^{ri} 
The likelihood factorises into functions of each parameter of the form e^{μTA}μ^{d}:
So the loglikelihood is:
Differentiating with respect to each of the four parameters gives:
Setting each of the derivatives to 0 and solving the resulting equations, we see that:
When there is more than one parameter to be estimated, the secondorder condition to check for maxima is that the Hessian matrix is negative definite, or equivalently, the eigenvalues of the Hessian matrix are all negative. The Hessian matrix is the matrix of the second derivatives. So, in this case, we consider the matrix:
Since this is a negative definite matrix, the maximum likelihood estimates of μ, ν, σ, ρ are:
Therefore, the maximum likelihood estimators are as given:
Coding of covariates: The effects of covariates on estimated transition intensities are analysed. This helps in determining variables that have a strong influence on the survival of HIV patients receiving antiretroviral therapy. These variables include age, peripheral neuropathy (PN), nonadherence (NA), gender, an orthogonal viral load residual covariate (RV). As proposed by Shoko et al.^{11} the effect of the orthogonal viral load variable is included in the continuoustime Markov model. The variables are coded as follows:
The residuals are obtained from a linear regression model of viral load on CD4 cell count. The residual covariate is independent of the CD4 count covariate correcting for collinearity effects. For more details on the collinearity argument, see Shoko et al.^{11}. The effects of the covariates on transition intensities, q_{ij}, for a patient h is given by the model:
For this model, the baseline transition intensities, , refer to a patient with age category 0 (over 45 years old), no LA, no PN, Gender = 0 (female), no NA and RV = 0 means lower viral load than expected. The transition intensities, q_{ij}, are presented in rates per year. q_{ij} are the elements of a 3×3 transition intensity matrix Q from a continuous timehomogeneous Markov process.
RESULTS
This section gives the results from the analysis of the data based on the HAD model given earlier in Fig. 1 of the methods section. Three states are defined namely: state HHIV state marked by CD4 cell count greater than 200 cell mm^{–}^{3}, state Athe AIDSdefining state, marked by CD4 cell count below 200 cells/mm^{3} and state DDeath. We start off by computing the transition counts. The transition counts are shown in Table 1.
Results from Table 1 above show fewer deaths recorded from CD4 cell count greater or equal to 200 (3.1%) than from the AIDSdefining state (15.6%). Once a patient has achieved a CD4 cell count greater or equal 200 (state H), there is a possibility of reverting to the AIDS state.
Next, the maximum likelihood estimators for the HAD model are computed using the continuoustime homogeneous Markov model.
Continuoustime homogeneous Markov model for the HAD model: In this subsection, we give estimates of the transition intensity matrix given by Q_{ij}, probability of each state being next, given by:
For example, in
Table 2:
The results are shown in Table 2.
Results from Table 2 show higher risks of death from the AIDS state (state A) than from state H where the CD4 cell count is above 200. The results also showed that for patients on antiretroviral therapy and in the AIDSdefining state, there are higher rates of recovery to a state of CD4 cell count above 200 cell mm^{3} compared to rates of transitions to the death state. Thus, confirmation that antiretroviral therapy improves the survival of HIV infected patients.
Next, we compute the total time spent in each state and the mean sojourn time for the HAD model as shown in Table 3.
The total times spent in each state reveal higher survival chances once the HIVinfected patient is alive and has transitioned out of the AIDSdefining state. Results from Table 3 show that patients with CD4 cell count above 200 and on cART can spend approximately 17 years before absorption into the death state. HIVinfected patients can also spend an average of about 4 years in state H before the transition to other states.
For the continuoustime homogeneous model with parameters estimated above, we plot the percentage prevalence in each of the states indicated in the HAD model for HIVinfected patients on combination antiretroviral therapy.
The prevalence plots from Fig. 2 show that the fitted model overestimates observed percentage prevalence in state 1 (state H). The fitted model gives a better fit of the observed patients in the AIDSdefining state 2 (state A). Observed death prevalence is underestimated by the fitted model. This suggests that transition intensities may not be constant throughout the whole study period. Hence, the need to fit a continuoustime nonhomogeneous Markov model. This is done below.
Table 1: 
Transition counts for the HAD model 

Table 2: 
Maximum likelihood estimates of the transition intensities, probability of each state being next and transition probabilities for the HAD model (95% CI are given in brackets) 

2* loglikelihood = 1977.895 
Table 3: 
Estimates of the total time spent in each state and sojourn times (95% CI are given in brackets) 

Table 4: 
Maximum likelihood estimates for the baseline transition intensities and hazard ratios for each time interval (95% CI are given in brackets) 

2* loglikelihood = 1935.603 

Fig. 2(ac): 
Prevalence for the continuoustime homogeneous Markov plots for the HAD model (a) State 1, (b) State 2 and (c) State 3 
Next, we fit a continuoustime inhomogeneous Markov model to assess the interval in which immune deterioration is experienced after patients have achieved an improvement in the immune system.
Results from Table 4 showed the maximum likelihood estimates for the baseline transition intensities and hazard ratios for the effects of time on transition intensities. The results show that in the second 0.5 years of treatment uptake, there is a reduction in the rates of immune suppression. The second 0.5 years of treatment uptake is also characterised by a reduction in deaths particularly from the AIDSdefining state. From 1 year of treatment uptake onwards, there is a significant increase in the rates of immune recovery from the AIDSdefining state. However, there is a slight increase in the rate of occurrences of deaths after 1 year of treatment compared to the 0.5 years to 1year interval. There is a further reduction in deaths for patients with CD4 cell count above 200, 1year posttreatment uptake compared to the 0.5 years to 1year interval. This means that chances of survival for HIVinfected patients increase once their CD4 cell count is above 200 and the effect of time for this to happen is crucial.
We now plot percentage prevalence for each of the states in the continuoustime nonhomogeneous model. The plots are shown in Fig. 3.
Results in Fig. 3 show a great improvement in fitting a nonhomogeneous model compared to the homogeneous model. The fitted model in Fig. 3 now gives an almost perfect fit of the observed percentage prevalence. In the next section, we analyse factors that contribute to the survival of HIVinfected patients on combination antiretroviral therapy.
Covariates that contribute to the survival of HIVinfected patients on cART: In this subsection, we analyse the effects of routinely collected viral load (an orthogonal viral load covariate (RV)), age, gender, peripheral neuropathy (PN) and nonadherence to treatment (NA) on the survival of HIVinfected patients on cART.
Results from Table 5 show that risks of immune suppression from CD4 cell count above 200 to the AIDSdefining state is attributed to nonadherence to treatment and having higher viral load than expected. Patients below the age of 40 years had significantly higher rates of immune recovery than their older counterparts. Results also show that male patients, patients below the age of 40 years, patients with a higher viral load than expected and nonadherent patients are at higher risk of death from the CD4 cell count level above 200.

Fig. 3(ac): 
Prevalence (%) for the continuoustime nonhomogeneous Markov model for the HAD model (a) State 1, (b) State 2 and (c) State 3 
Figure 4 shows that the inclusion of the effects of covariates in the Markov model results in a good fit of the observed percentage prevalence for each of the states in the HAD model. This has also resulted in better prediction of mortality (state 3). This is shown by death prevalence that accumulated to slightly below 20% by the end of the 4 years of treatment uptake. This gives the best prediction of mortality since 60 deaths out of 320 patients occurred during the study. This is approximately equal to 18.8%.

Fig. 4(ac): 
Prevalence (%) for the continuoustimehomogeneous Markov model with covariates (a) State 1, (b) State 2 and (c) State 3 
Assessment of the fitted models: In this section, we assess the fitted models by performing a likelihood ratio test as well as comparing the 2 *loglikelihood(2LL) ratios for the fitted models. The results are given in Table 6.
Results from Table 6 show that the likelihood ratio test for the timehomogeneous model versus the timehomogeneous model with covariates is in favour of a timehomogeneous model with covariates. A comparison of the nonhomogeneous model with the homogeneous model with covariates is in favour of the homogeneous model with covariates.
Table 5: 
Maximum likelihood estimates for the baselines transition intensities with hazard ratios for each covariate (95% CI are given in brackets) 

2 * loglikelihood = 1391.64 
Table 6: 
Likelihood ratio tests for the fitted models 

Timehomogeneous model: 2LL = 1977.895, Timenonhomogeneous model: 2LL = 1935.603, timehomogeneous model with covariates: 2LL = 1391.64 
The timehomogeneous model with covariates also has the lowest2*loglikelihood compared to all the other fitted models. This shows that the survival of patients on cART is best explained by a time homogenous model with covariates including the routinely collected viral load covariate.
DISCUSSION
In this study, the survival of HIVinfected patients after ART initiation is analysed. Continuoustime Markov models are fitted based on the states from an HIV, AIDS and DEATH (HAD) states model. These states are based on CD4 cell count. Factors that affect the survival of HIVinfected patients on ART were analysed. These, among others, include age, gender, routinely collected viral load, nonadherence, peripheral neuropathy and time on cART.
Results from Table 5 on estimated hazard ratios showed that patients below the age of 40 years achieved a normal CD4 cell count 1.15 times faster than their older counterparts. This is corroborated by findings from a previous study which concluded lower mean CD4 increases for older patients than younger patients^{12}. However, these patients have higher risks (Hazard ratio = 34.6, CI: (0.00001, 1756)) of mortality compared to their older counterparts.
Time on treatment also had a significant effect on the survival of patients. Thus, after one year of treatment uptake, patients have accelerated transitions from the AIDS state to the HIV state defined by CD4 cell count above 200 cell/mm^{3}.
Males have generally lower chances of survival than their female counterparts. Results showed that males also have slower rates of immune recovery than females. Sieleunou et al.^{3} also concluded that gender is a predictor of mortality for patients on antiretroviral therapy.
Nonadherence to treatment reduces the chances of survival for patients on cART. Nonadherence accelerates (Hazard ratio = 2.2) transitions from CD4 cell count state above 200 cell/mm^{3} to the AIDS state defined by a CD4 cell count below 200 cells/mm^{3}. Patients who were nonadherent to treatment are 3.8 times more likely to transit from the CD4 state above 200 cell/mm^{3} to death compared to patients who were adherent to treatment. This is corroborated by the findings from Seyoum et al.^{5} who observed that lowadherence to ART is associated with increased mortality.
Patients with higher viral loads than expected had higher risks (Hazard ratios = 11.147, CI: (2.6110, 47.588)) of HIV progression to the AIDS state and high risks (Hazard ratio = 1.091, CI: (0.0914, 13.026)) of mortality from a CD4 cell count above 200 cell/mm^{3}. Thus, due to the effects of treatment and time spent on cART, patients are expected to have suppressed viral loads (possibly undetectable viral load) leading to an improved immune system. However, if the viral load remains higher than expected, this increases risks of immune deterioration even after achieving normal CD4 cell counts and consequently, mortality risks are increased.
A continuoustime homogeneous model without covariates, a continuoustime nonhomogeneous model and a continuoustime homogeneous model with covariates were fitted for the data. These models were assessed to select the model that best predicts the survival of HIVinfected patients on cART. Form both the continuoustime homogeneous model and the nonhomogeneous model, cumulative deaths by the end of four years of treatment uptake had increased to close to 40%. For the model with covariates, cumulative deaths prevalence had increased to below 20%. However, from the data only 60 out of 320 deaths were recorded throughout the whole study period which is approximately equal to 18.8%. Therefore, we conclude that a continuoustime homogeneous Markov model with the effects of covariates, including the orthogonal viral load covariate, gives the best prediction of survival of HIVinfected patients on cART than the other two models. This corroborates with the findings from studies by Shoko et al.^{11} who proposed the inclusion of both viral load monitoring and CD4 cell count in one model to account for the aspect of mortality which one variable can fail to account for without the inclusion of the other.
CONCLUSION
Recovery from AIDS state by HIV infected patients on cART is likely to occur after one year of cART treatment. However, if the viral load remains higher than expected, this increases risks of immune deterioration even after having achieved normal CD4 cell counts and consequently, mortality risks are increased.
SIGNIFICANCE STATEMENT
This study concludes that constant monitoring of an HIV/AIDS patient’s viral load can be beneficial for reducing the risks of immune deterioration and mortality. This study will help the researchers to uncover the critical areas regarding HIV/AIDS progression and causes of mortality that many researchers were not able to explore. Thus, a new theory on HIV, AIDS, Death (HAD) Markov model may be arrived at.