Biomedical Research

DOI: 10.35841/biomedicalresearch.30-19-142

Abstract

Background: Markov models are widely used in studying progression of chronic diseases using the time homogeneity approach. Relatively few studies use the time non-homogeneous approach to model the progression of HIV/AIDS based particularly in the use of viral load monitoring mainly due to fact that that viral load is expensive to measure frequently hence not readily available in the format required. Non-homogeneous models reveal the interval in which viral rebound occurs. Methods: This study compares the time-homogeneous Markov approach to the time nonhomogeneous approach in predicting HIV/AIDS progression. A piece-wise constant approach to non-homogeneous modelling is used. Results: This study reveals that a piece-wise constant Markov approach to non-homogeneous modelling which allows variation in transition rates within intervals [0,0.5) years and [0.5,∞) years is better for predicting HIV/AIDS progression compared to the time-homogeneous approach. Failure to combination antiretroviral therapy (cART) increases significantly the rate of viral rebound and deaths 6 months post treatment initiation. Conclusion: Hence, there is need to closely monitor the uptake of antiretroviral drugs in the first six months of treatment to increase the chances of survival for HIV/AIDS patients.

Keywords

Non-homogeneous Markov model, Time-homogeneous Markov model, HIV/AIDS progression, Viral load.

Introduction

Combination antiretroviral therapy (cART) reduces HIV viral load by blocking replication of the virus [1-3]. Hence, constant monitoring of viral load counts is beneficial over routine monitoring of CD4⁺ cell count since immunodeficiency can be prevented by early detection of virologic failure thereby enabling proper management of HIV/AIDS patients by interested stakeholders [4].

The suppression of HIV viral load from very high levels, greater than 500 000 copies/mL, to suppressed levels due to the effects of cART can be considered as a multi-state process based on the viral load including the endpoint, death state. Patients are then monitored only at visit times resulting in the exact time the transition occurred to be unknown [5]. Homogeneous and non-homogeneous Markov processes are an important field of research when dealing with interval censored data in which the exact time of transition is not known. A Markov process is a stochastic process in which the probability of transition to a future state given the current state is independent of the past history. Thus, a Markov model possesses the memoryless property involving random movements between states that occur at regular or irregular intervals [6].

The constant hazard assumption of the time-homogeneous models is unrealistic when modelling evolution in chronic diseases [7] because they put severe limitations on the previous behaviour of the disease [8]. However, most applications to HIV/AIDS disease assumed a homogeneous process where transition probabilities only depend on the elapsed time between observations [6,9-11].

Application of non-homogeneous Markov models to study the evolution of HIV allows computations of estimates of the time in which HIV patients spend in each state but also the randomness of different states in which the states can evolve known as a stochastic process [12]. Saint-Pierre et al. used a piece-wise approach to non-homogeneous Markov modelling on asthma patients. They argued that piece-wise constant transition intensities help in preserving the tractability of constant intensities [8]. However, non-homogeneity is not treated in survival analysis as extensively as homogeneous models and in particular for HIV/AIDS disease progression when based on viral load monitoring [13].

This paper models the virology of HIV using a piece-wise constant approach to form a time non-homogeneous Markov process. The virology of HIV is split into mutually exclusive states defined by viral load including the absorbing state, death. The use of Markov models with piece-wise constant transition intensities helps in preserving the tractability of constant intensities [8].

In the next section we describe the HIV/AIDS data that is used in this study and also how the model is formulated basing on the data. In section 3, we have results and analysis where three different models are fitted and their appropriateness is assessed using the Akaike Information Criteria (AIC), likelihood ratio tests as well as the plots of percentage prevalence in each state. The last section discusses and concludes the findings.

Materials and Methods

Data description

The data set used in this study has been described in previous studies [14-17]. At treatment initiation the variables age, Viral Load Baseline (VLBL) and CD4 baseline (CD4BL) are described in Table 1.

Table 1: Descriptive Statistics for variables age, baseline viral load and baseline CD4+ cell count.

For each and every visit, blood samples were obtained for each patient and stored frozen until assayed. Plasma HIV RNA was measured using an amplicator HIV-1 monitor assay kit which has a limit of sensitivity ranging from 50 copies/μL to 500 000 copies/μL.

Table 2 shows the possible transitions between the defined states from t=0 to t=1.5 years of treatment uptake.

Table 2: Number of HIV/AIDS patients in each viral load state from t=0 to t=0.5 years.

At enrolment, most patients were in state 3 followed by state 4. State 1 had the least number of patients but during the first 0.25 years, it had the highest number of patients followed by state 2. From t=0.25 to t=0.5 years, the number of patients in state 3 increased from 6 to 13, an indication of viral rebound since the number of patients with higher viral loads (state 3 and state 4) reduced collectively only by 3. There is need to investigate the causes of viral rebound for HIV/AIDS patients since they are on treatment.

Model formulation

The non-homogeneous continuous time Markov model is formulated by splitting HIV/AIDS progression into 5 viral load defined states followed by the end point, death, which is state 6. The HIV/AIDS progression states are:

These 6 states define the HIV/ AIDS disease process is denoted by X (t) as shown above. The machine that is used to detect the viral load used thresholds of [50; 500 000) such that viral load below 50 copies /mL is classified as undetectable. The other covariates included in the model are.

Formulation of the model is based on the approach by Saint- Pierre and others [8] where the study period [τ_(l-1),τ_l) for, l=1,….,r+1, and considering a vector of artificial timevarying covariates

such that:

for

(1)

And the model with transition intensities;

(2)

q_ij0 is the baseline transition intensity corresponding to the interval [τ₀,τ₁], β*_ij is a vector of regression coefficients associated with the artificial time-varying covariates. For this model the transition intensities are step-functions of time defined for each interval as follows:

Incorporating the effects of covariates, the model becomes

(3)

Where q_ij0 is the baseline transition intensities for the interval 0 ≤ t<0.5 with covariates set to their means, r is the log-linear effect of the r interval on the baseline transition intensity and Z*t is the artificial time-dependent covariate, β_ij is the log-linear effect of the relating the instantaneous rate of transitions from state i to state j to the covariates

Statistical analysis

All the analysis is done using the Multi-State Modelling [MSM] package for multistate modelling in R software developed by Jackson [18]. The process of identifying the appropriate model started off by fitting time homogeneous Markov model for the data. The time homogeneous model converges and has -2Log-likelihood of 2799.465. From this time homogeneous model, percentage in each state are estimated and plotted as shown in Figure 1.

Figure 1: Prevalence plot for the time homogeneous Markov model.

A comparison of the expected prevalence and observed percentages in each state from Figure 1 show that the predicted number of HIV infected individuals who die (State 6) is underestimated by the model from about 1 year onwards. The number of individuals alive and in state 1 is overestimated by the model from 1 year onwards. Such discrepancies, according to Jackson, indicate a possibility that the transition rates vary with time since the beginning of the process. In this particular case, it could be the time on treatment therapy [19,20]. This calls for the need to fit a Markov process that is non-homogeneous which can be handled by fitting a piecewise constant function.

Individuals in this study were followed up after every 6 months, thus we start off with a 9-segment non-homogeneous Markov model with 0.5 year intervals. Although the 9-segment model did not converge to maximum likelihood, the prevalence plots for each state in the model helped in revealing intervals with constant transition intensities.

We further investigated the best piece-wise constant function for the data by considering different change points. Most of the models did not converge to a maximum likelihood except for the models;

1. With 2-segments and having a change point at 0.5 years.

2. With 2-segments and having a change point at 1 year.

3. With 3-segments and having change points at 0.5 and 1 year.

In the next section, we base our analysis on these three models and also the time homogeneous model. We further investigate the effects of covariates; gender, non-adherence and CD4 baseline on the chosen model.

Results

For the 2-segment model, the transition intensities are defined for the following intervals.

Where q_ij,0 is the baseline transition and q_ij,1 is the transition intensity matrix for the interval t ≥ 0.5years and β*_ij1 is the regression coefficient associated with the artificial timedependent covariate Z₁ *(t). The point estimates of parameters and their corresponding confidence intervals are shown in Table 3.

Table 3: Baseline transition intensities for the 2-segment non-homogeneous model and the time varying log-linear effects.

The results from Table 3 also narrow confidence intervals for transitions from state 1 and from state 2 except for the transition from state 2 to 6 (death). This shows that estimates of transitions from these states are statistically significant. However, most of the point estimates have wide confidence intervals and therefore are not statistically significant. This could be due to the fact that within 6 months of treatment uptake, most patients had made transitions to either state 1 or state 2 as shown by the prevalence plots in Figure 1.

Next we compute estimates of parameters for the 3-segment non-homogeneous Markov model with half-year and oneyear change points, with transition intensities defined for each interval as follows:

is the baseline transition intensities for the interval, 0 ≤ t <0.5, β*_ijr tis the log-linear effect of the r^th interval on the baseline transition intensity and Z*(t) is the artificial time-dependent covariate. The estimated parameters are shown in Table 4.

Table 4: Effects of half-year and one-year changes in time on the baseline transition

Results from Table 4 show an improvement on the estimated parameter values. This is shown by the confidence intervals that are narrower compared to the confidence intervals shown in Table 3 for the 2-segment model. However, transitions to death from states 2, 3, 4 and 5 have got very wide confidence interval. Only mortality rates from state 1 (undetectable viral load) decrease significantly with time as indicated by the narrow confidence intervals and exclusion of zero in the confidence intervals.

Next we present the results from a 2-segment piece-wise Markov model together with the effects of covariates. The model is given as follows:

q_ij0 is the baseline transition intensities for the interval 0 ≤ t <0.5 with covariates set to their means, β*_ijk is the log-linear effect of the k^th interval on the baseline transition intensity and Z*(t) is the artificial time-dependent covariate, β*_ij is the log-linear effect of the relating the instantaneous rate of transitions from state i to state j to the covariates.

The results are shown in Table 5.

The results from Table 5 show that during the first 0.5 years of treatment uptake (baseline) the rates of viral load suppression are higher than the rates of viral rebound regardless of the state. The highest transition rates are noted from a viral load above 500 000 copies/μL to a viral load between 100 000 and 499 999 copies/μL. During the first six months most transitions to death occurred from a viral load state between 100 000 and 500 000 copies/μL. This is mainly attributed to initiating treatment with a CD4 baseline below 200 cells/ m^3.

Having a CD4 baseline below 200 cells/m^3 at treatment initiation contribute significantly to the increase in viral suppression to undetectable level as shown by a narrow confidence interval and exclusion of zero in the confidence interval. For these patients, there is also a significant increase in viral suppression from state 4 to state 3. Though not significant, there is reduction in viral suppression from state 3 to state 2.

For non-adherent patients, the results reveal a significant reduction in viral suppression to undetectable levels (state 2 to state 1). Thus, adherence to treatment increases significantly the transition rates to undetectable levels. Non-adherence increases transitions to viral rebound from undetectable viral load (state 1 to state 2). Although some of the estimated parameters are not significant, the results show that nonadherence to treatment contribute more to the increased rates of viral rebound than viral suppression.

Males in the study experience a reduction in the rates of viral suppression to undetectable levels, but once the undetectable viral load is reached they have reduced rates of viral rebound. This means that for their female counterparts, rates of viral suppression to undetectable levels are higher than that of males. Deaths of males are more likely to occur from viral load states between 10 000 and 500 000 copies/μL compared to their female counterparts.

The results also show reduction in the rates of transition to better states (viral suppression) for patients who initiated treatment with a viral load level above 10 000 copies/μL. These patients experienced increased transitions to death regardless of the viral load state. Thus, starting treatment with a viral load above 10 000 copies/μL increases the mortality rates and also accelerates diseases progression.

From 6 months (0.5 years) of treatment uptake onwards, most of the transitions between live states are significant except for the transition from state 3 to state 4 and state 5 to state 4. The results show a significant reduction in viral suppression to undetectable levels from 6 months onwards but once the undetectable viral load is reached there is reduction in viral rebound and also reduction in death occurrences. However, if the viral load is still detectable after 6 months, it increases the occurrence of deaths.

Assessment of the fitted models

In this subsection, we assess the fitted models using several techniques: prevalence plots, contour plots and estimates of log-likelihoods and Akaikie Information Criteria (AIC). Prevalence plots give a rough indication of the goodness of the fitted models and contours plot the likelihood surface with respect to 2 parameters default to plus or minus 2 standard errors obtained from the Hessian matrix at the maximum likelihood estimate.

In Figure 2 we compare the percentage prevalence plots for each of the states for the a two-segment non-homogeneous Markov model with change point at 0.5 years, 3-segment non-homogeneous Markov model and the 2-segment model with change point at 1 year respectively.

Figure 2: Percentage Prevalence plots for: (a) the 2-segment model with change point at 0.5 years, (b) the 3-segment model with change points at 0.5 and 1 year and (c) the 2-segment model with change point at 1 year.

Results from the percentage prevalence plots show that the fitted model gives a perfect fit of the observed data for patients in state 2,3,4 and 5. However, the expected percentage prevalence plot for the 2-segment non-homogeneous model with change point at 1 year overestimates the observed plot for patients in state 1. From 3 years onwards, expected percentage prevalence for mortality (state 6) underestimate the observed mortality from all fitted models although the 2-segment model with change point at 1 year underestimates mortality from 2 years onwards. Compared to the prevalence plot for the time homogeneous model shown in Figure 1, these plots show an improved fit to the observed data.

The prevalence of State 1 is down sloping from half a year on. This is due to non-adherence to treatment as shown by the results in Table 4.

Before computing estimates of AICs we plot the contours for the fitted models. These contours help to diagnose any irregularities in the likelihood surface. Thus, they present a graphical visualisation of the fitted surface. For biological data, these contours are expected to be elliptical. Figure 3 shows contours for the fitted models.

Figure 3: Contour plots for: (a) the time homogeneous Markov model, (b) the 2-segment model with change point at 1 year, (c) the 2-segment model with change point at 0.5 years and (d) the 3-segment model with change points at 0.5 and 1 year.

Figure 3 shows that contours for the time homogeneous Markov (top left) are not elliptical but the contour for the three non-homogeneous models are elliptical. So the nonhomogeneous model give a better fit the data compared to the homogeneous model. However, the contours for the time homogeneous model (top left) and the contours for the 2-segment with change point at 1 year (top right) are not symmetric; hence the models cannot explain this data fully. The contours for the non-homogeneous models; 2-segment with change point at 0.5 years (bottom left) and 3-segment with change points at 0.5 and 1 year (bottom right) are both elliptical and evenly distributed contours plots with symmetrical surfaces that peak at the centre (indicated by the white region). Hence these models give an adequate explanation of the data.

Table 6 shows estimates of the -2xlog-likelihood (-2 × LL), log-likelihoods (LL) are shown in brackets, the degrees of freedom for each of the fitted models and the Akaike information criteria (AIC). The AICs are calculated using the

formula;

AIC = −2log(likelihood) + 2× df

For example, for the homogeneous model AIC=2799.465+2x13=2825.465 as shown in the Table 6.

The results show that the time homogeneous Markov model has the highest AIC and log-likelihood compared to the non-homogeneous Markov models. From the fitted non-homogeneous Markov model, the 3-segment model, with change points at 0.5 and 1 year, has the highest loglikelihood compared to all the other fitted models followed by the 2-segment model with change points at 0.5 years. However, a further assessment basing on the AICs reveal that a 2-segment model with change points at 0.5 years fits the data better than the 3-segment model since it has the lowest AIC. Effects of the covariates were included in the 2-segment non-homogeneous model and the model yielded the lowest AIC.

Discussion

In this study, a comparison of the time-homogeneous Markov model and the non-homogeneous Markov model to estimate the progression of HIV/AIDS on viral load monitoring is done. For this model, the expected percentages in state 1 overestimated the observed percentages and the predicted number of individuals who died was underestimated by the model from 1 year of treatment uptake onwards. This indicated the need to fit a non-homogeneous model in which transition intensities are piece-wise constant. Nonhomogeneous models with different change points were fitted. However, most of these models did not converge to a maximum likelihood except for the 2-segment model with change point at 0.5 years, 2-segment model with change point at 1 year and 3-segment model with change points at 0.5 and 1 year. From these models, non-homogeneous models and the time homogeneous model, assessment of the best model was done using the percentage prevalence plots, contour plots, perspective plots, CIs, log-likelihoods and the AICs. The model with the lowest AIC, the 2-segment model with change point at 0.5 years, was considered as the model that best explains HIV progression in patient on treatment under viral load monitoring for this data set.

This means that HIV progression is best described by a non-homogeneous Markov model with a change point at 0.5 years. This is mainly influenced by the results from viral load monitoring which are that, most of the patients on antiretroviral therapy reach the undetectable viral load within the first six months of treatment uptake (14,16). Thus, when monitoring HIV/AIDS patients, one should ensure that viral load suppression is reached as this reduces mortality rates.

Results from the fitted models show wider confidence intervals for the estimated transition rates to death. However, transitions within the 5 live states defined by viral load had narrow confidence intervals indicating significance in the estimated parameters, particularly the 2-segment model with change point at 0.5 years and covariates, in predicting transitions between live states.

The results also revealed a slower response to treatment when patients start treatment when viral load levels are above 10 000 copies/μL than when treatment is initiated when the viral load levels are below 10 000 copies/μL in the first 6 months [17]. Patients who initiated therapy with a viral load level above 10 000 copies/μL experienced accelerated rates of transition to death [14]. In 2018, Shoko et al. observed the same finding although they used a time homogeneous Markov modelling approach [17].

The results from this study revealed reduction with time in viral suppression to undetectable levels for males and patients who were non-adherent to treatment [14]. For males, once the undetectable viral load is reached, there is reduction in viral rebound particularly from 0.5 years onwards.

This study shows a significant reduction in transitions to death from 0.5 years onwards for patients who have achieved an undetectable viral load despite the challenges of nonadherence. Thus, time non-homogeneity is the best approach to modelling HIV/AIDS progression for patients receiving cART.

Ethics and Consent to Participate

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), in accordance with 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 prior to their involvement; and data obtained was stripped of personal identifiers to ensure the confidentiality of the participants.

Funding

The dataset used was generated from funding awards to POB by the South African Medical Research Council (RCDI) through funding received from the South African National Treasury, the South African National Research Foundation (GUN109312), and the University of Venda. The contents are solely the responsibility of the authors and do not necessarily represent the official views of the South African Medical Research Council , the National Research Foundation, or the University of Venda.

Availability of Data and Materials

The data is available upon a submitted request.

Author Contributions

CS drafted the manuscript. POB provided the data used for the analysis. CS, DC and POB contributed to the analysis and interpretation of the data. All authors participated in critical revision of the manuscript drafts and approved the final version.

Competing Interests

The authors declare that they have no competing interests.

Consent for Publication

Not applicable.

References

1. Brennan AT, Maskew M, Sanne I, Fox MP. The Interplay between CD4 Cell Count, Viral Load Suppression and Duration of Antiretroviral Therapy on Mortality in a Resource-Limited Setting. Trop Med Int Health 2013; 18: 619-631.
2. Cohen MS, Chen YQ, McCauley M, Gamble T, Hosseinipour MC, Kumarasamy N, Hakim JG, Kumwenda J, Grinsztejn B, Pilotto JHS, Godbole SV, Mehendale S, Chariyalertsak S, Santos BR, Mayer KH, Hoffman IF, Eshleman SH, Piwowar-Manning E, Wang L, Makhema J, Mills LA, Bruyn GD, Sanne I, Eron J, Gallant J, Havlir D, Swindells S, Ribaudo H, Elharrar, Burns D, Taha TE, Nielsen-Saines K, Celentano D, Essex M, Fleming TR. Prevention of HIV-1 infection with early antiretroviral therapy. N Engl J Med 2011; 365: 493-505.
3. Cole SR, Hernan MA, Anastos K, Jamieson BD, Robins JM. Determining the effect of highly active antiretroviral therapy on changes in human immunodeficiency virus type 1 RNA viral load using a marginal structural left-censored mean model. Am J Epidemiol. 2007; 166: 219-227.
4. Lecher S, Williams J, Fonjungo PN, Kim AA, Ellenberger D, Zhang G, Toure CA, Agolory S, Appiah-Pippim G, Beard S, Borget MY, Carmona S, Chipungu G, Diallo K, Downer M, Edgil D, Haberman H, Hurlston M, Jadzak S, Kiyaga C, MacLeod W, Makumb B, Muttai H, Mwangi C, Mwangi JW, Mwasekaga M, Naluguza M, NgAngA LW, Nguyen S, Sawadogo S, Sleeman K, Stevens W, Kuritsky J, Hader S, Nkengasong J. Progress with Scale-Up of HIV Viral Load Monitoring ? Seven Sub-Saharan African Countries, January 2015?June 2016. R Morb Mortal Wkly Rep 2016; 65: 1332-1335.
5. Chenand B, Zhou XH. Non-homogeneous Markov process models with informative observations with an application to Alzheimer?s disease. Biom J 2011; 53: 444-463.
6. Lee S, Ko J, Tan X, patel I, Balakrishnan R, Chang J. Markov Chain Modelling Analysis of HIV/AIDS Progression: A Race Based Forecast in the United States. Indian J Pharm Sci 2014; 76: 107-115
7. Lange JM, Hubbard RA, Inoue LYT, Minin V. A joint model for multistate disease processes and random informative observation times, with applications to electronic medical records Data. UW Biostatistics Working Paper Series, 2014: 401.
8. Saint-Pierre P, Combescure C, Daures JP, Godard P. The analysis of asthma control under a Markov assumption with use of covariates. Stat Med. 2003; 22: 3755-3770.
9. Shoko C, Chikobvu D. Time-homogeneous Markov process for HIV/AIDS progression under a combination treatment therapy: cohort study, South Africa. Theor Biol Med Model 2018; 15: 3.
10. Grover G, Gadpayle KA, Swain PK, Deka B. A Multistate Markov Model Based on CD4 cell count for HIV/AIDS Patients on Antiretroviral Therapy (ART). Int J Stat Med Res 2013; 2: 144-151.
11. Mathieu E, Foucher Y, Dellamanica P, Daures JP. Parametric and non homogeneous Semi-Markov process for HIV control. Methodol Comput Appl Probab 2007; 9: 389-397.
12. Ocarn-Riola R. Non-homogeneous Markov Processes for Biomedical data Analysis. Biom J 2005; 47: 369-376.
13. Dessie ZG. Modeling of HIV/AIDS dynamics evolution using non-homogeneous semi-Markov process. SpringerPlus 2014; 3: 537.
14. Shoko C and Chikobvu D. Determinants of viral load rebound on HIV/AIDS patients receiving antiretroviral therapy: Results from South Africa. Theor Biol Med Model 2018; 15: 10.
15. Shoko C, Chikobvu D. A superiority of viral load over CD4 cell count when predicting mortality in HIV patients on therapy. BMC Infectious Disease. 2019; 19: 169.
16. Chikobvu D and Shoko C. A Markov Model to estimate Mortality due to HIV/AIDS using CD4 cell counts based states and viral load: A Principal Component Analysis approach. Biomedical Research. 2018; 29 (15): 3090-3098
17. Shoko C, Chikobvu D, Pascal O. Bessong. A Markov Model to estimate Mortality due to HIV/AIDS using viral load levels based states and CD4 cell counts: A Principal Component Analysis approach. Infect Dis Ther 2018; 4: 457-471.
18. Jackson C, Multi-state modelling with R: the msm package Version 1.6.4. MRC Biostatistics Unit. 2016
19. Silveira MP, De Lourdes DM, Leite JC, Pinheiro CA, Da Silveira VL. Predictors of undetectable plasma viral load in HIV-positive adults receiving antiretroviral therapy in Southern Africa, Braz J Infect Dis 2002; 6: 164-171.
20. Azia IN, Mukumbang FC, Van Wyk, B. Barriers to adherence to antiretroviral treatment in a regional hospital in Vredenburg, Western Cape, South Africa. South Afr J HIV Med 2016; 17: 476.