 Research Article
 Open Access
 Open Peer Review
 Published:
Multistate models for the analysis of timetotreatment modification among HIV patients under highly active antiretroviral therapy in Southwest Ethiopia
BMC Infectious Diseases volume 17, Article number: 453 (2017)
Abstract
Background
Highly active antiretroviral therapy (HAART) has shown a dramatic change in controlling the burden of HIV/AIDS. However, the new challenge of HAART is to allow longterm sustainability. Toxicities, comorbidity, pregnancy, and treatment failure, among others, would result in frequent initial HAART regimen change. The aim of this study was to evaluate the durability of first line antiretroviral therapy and to assess the causes of initial highly active antiretroviral therapeutic regimen changes among patients on HAART.
Methods
A Hospital based retrospective study was conducted from January 2007 to August 2013 at Jimma University Hospital, Southwest Ethiopia. Data on the prescribed ARV along with start date, switching date, and reason for change was collected. The primary outcome was defined as the timetotreatment change. We adopted a multistate survival modeling approach assuming each treatment regimen as state. We estimate the transition probability of patients to move from one regimen to another.
Result
A total of 1284 ART naive patients were included in the study. Almost half of the patients (41.2%) changed their treatment during follow up for various reasons; 442 (34.4%) changed once and 86 (6.69%) changed more than once. Toxicity was the most common reason for treatment changes accounting for 48.94% of the changes, followed by comorbidity (New TB) 14.31%. The HAART combinations that were robust to treatment changes were tenofovir (TDF) + lamivudine (3TC)+ efavirenz (EFV), tenofovir + lamivudine (3TC) + nevirapine (NVP) and zidovudine (AZT) + lamivudine (3TC) + nevirapine (NVP) with 3.6%, 4.5% and 11% treatment changes, respectively.
Conclusion
Moving away from drugs with poor safety profiles, such as stavudine(d4T), could reduce modification rates and this would improve regimen tolerability, while preserving future treatment options.
Background
The implementation of HAART at a large scale has shown a dramatic change in controlling the burden of HIV/AIDS. Various studies from developed as well as developing countries have reported an improvement in CD4 cell counts following ART initiation [1–6] and decreases in mortality [7].
Currently, Ethiopia and most resourcelimited countries have adopted nonnucleoside reverse transcriptase inhibitors (NNRTIs) based therapy. They use either NVP or EFV plus two nucleoside reverse transcriptase inhibitors (NRTI) as a first line therapy for adults and adolescents. This combination has been shown to be efficacious, are generally less expensive, and have generic formulations [7]. However, the new challenge of HAART is to allow longterm durability. Many patients will be forced to modify or switch their treatment regimens for various reasons, including poor drug tolerance, drug toxicities, drugtodrug interactions, pregnancy and treatment failure [8–10].
Studies from developed and developing countries have shown that a substantial number of patients (up to 69%) may modify their regimen overtime, where 25%  44% of them modify their initial treatment within the first years of treatment [8–14]. Drug related toxicity was the most common reason for treatment modification [8, 9, 11, 13, 15–18] and this can be an important barrier to adherence and potentially lead to treatment failure [19]. The majority of these studies found that patients that receive d4T as a part of their treatment were at increased risk of treatment modification due to toxicity [13, 15–18], which raises questions about the continued role of d4T in firstline treatment. The WHO had revised its guidelines on the use of antiretroviral drugs several times. Recently, WHO recommends to move away from d4T giving preference to the use of TDF and AZT in standard firstline therapy when possible [20, 21]. Due to cost and management of toxicity, however, the transition from d4T to TDF has been slow in resourcelimited settings [20, 21].
The high rate of HAART switching emphasizes the complexity of managing these therapies. Given the limited number of secondline treatment options available in resourcelimited settings, maximizing regimen durability by minimizing the rate of treatment modification and rates of treatment failure amongst those on firstline regimens is vital to extend firstline treatment options and prevent premature initiation of secondline therapy. In order to achieve this goal, key reasons for changes in ART regimens should be studied and durable regimens should be identified for recommendations. Moreover, evaluating the influence of initial ART regimens on the likelihood of treatment modification has a vital role in determining what treatment to initiate and what treatment to preserve. Although relevant data on patients’ longterm experience on ART from resource limited settings are less commonly available, some investigators have described the reasons for modification of HAART and compared durability of individual ARV’s using routine clinical programme data [15–18, 22–27]. However, the majority of these studies have had short followup times and consider only first time regimen switching or first time single drug substitution with no distinction made between NNRTI and NRTI substitutions.
Therefore, this study aims to compare the durability of firstline ART regimens and investigate reasons for treatment modification in patients under HAART in Jimma university specialized hospital. For this purpose, we adopted a multistate survival modeling approach assuming each treatment regimen as state. We estimate the transition probability of patients to move from one regimen to another in general as well as due to a specific event that triggers the move. The proposed model allows modelling of both the occurrence of different event types (such as, single drug substitution or regimen switch) and the occurrence of subsequent events, the latter potentially of different types.
Methods
Data
Description of the cohort
The data used for this study were obtained from Jimma University Specialized Hospital HIV/AIDS clinic, located 352 Km Southwest of Addis Ababa, Ethiopia. The Hospital gives VCT, PMTCT and free ART service for people living in Jimma Town and Southwest Ethiopia. Patients begin ART after they have been checked for medical eligibility and are counseled for adherence for ART. Patients presenting with WHO stage 4 disease and/or a CD4 count lower than 200 were eligible to start ART. Those who started ART, have a regular followup for drug adverse effects, management of opportunistic infection, TB screen and counseling related to family planning. In addition, CD4 count is measured at each visit. Viral load measurement is not available. Adverse event monitoring is conducted by clinicians during medical visits in accordance with national guidelines.
Decisions on which treatment regimen to start or substitute are made by the clinician in consultation with the patient. During the study period, the standard firstline regimens to be initiated were d4T + 3TC + NVP (1), d4T + 3TC + EFV (2), AZT + 3TC + NVP (3), AZT + 3TC + EFV (4), TDF + 3TC + EFV (5), and TDF + 3TC + NVP (6). If a patient suffered from side effects/toxicities related to the NRTI’s (d4T, AZT or TDF), and was not in need of secondline therapy for virologic failure, the recommendation was to substitute d4T with either AZT or TDF, to substitute AZT either with d4T or TDF, and to substitute TDF with either d4T or AZT. Similarly, patients initiated on the NNRTI EFV could substitute with NVP, while those on NVP could substitute with EFV. In addition, upon the recommendation of WHO, in Oct 2012 the hospital started to phase out d4T backbone by replacing either AZT or TDF for patients who were on d4T based regimen.All data, including demographic, clinical conditions, laboratory test results and medications are recorded and entered in to the database by a data entry clerk at the clinic. In addition, data on prescribed ARV along with start and stop dates of the drug and reasons for discontinuation are documented. Use of Jimma University Hospital HIV/AIDS clinic data and analysis of deidentified data was approved by the Human Research Ethics Committee of Jimma University.
Study population
All ART naive patients, aged 18 years or older and who initiated a standard, publicsector, firstline ART regimen at the clinic between between January 1, 2007 and December 31, 2011 were eligible for this analysis. The data was closed for analysis at the end of August, 2013.
Outcome
The primary outcome was defined as the timetotreatment change (treatment modification or regimen switching). For the purpose of this study, treatment change is defined as changing at least one ARV in the regimen without initiating a secondline therapy. ARV dosage adjustments were not considered as treatment change. Time zero was defined as the day of ART initiation and each recurrent treatment change time was measured from the beginning of the patient’s ART initiation in months. Persontime of the study subject ended at the earliest of initiation on secondline therapy, lost to follow up, death, transfer or closure of the data set for analysis (August 25, 2013).
Multistate survival model
Model formulation
Possible transition between treatment combinations are presented in Fig. 1 which illustrates the treatment history of patients under ART. From here onward we use the term state to denote a specific treatment combination. The model has 6 transient states (which represents 6 firstline treatment combinations): d4T + 3TC + NVP (1), d4T + 3TC + EFV (2), AZT + 3TC + NVP (3), AZT + 3TC + EFV (4), TDF + 3TC + EFV (5), and TDF + 3TC + NVP (6). The model assumes that every patient can switch to all the regimen at one point or other. However, a patient can only switch to one regime at a time. For example, a patient who started treatment with d4T + 3TC + NVP State 1) is at risk of making one of the following transitions at a particular time; 1→2,1→3,1→4,1→5 and 1→6. If the patients made the transition 1→2 or 1→3, the subject has undergone a single drugsubstitution (treatment modification). Transition 1→2 implies that the patient has substituted their NNRTI’s NVP by EFV without changing their NRTI treatment (d4T). However, transition 1→3 implies that the patient has substituted their NRTI’s d4T by AZT without changing their NNRTI treatment. Transitions 1→4 or 1→5, imply regimen switching, substituting both NNRTI and NRTI at the same time. After making one of these possible transitions patients will be at risk of making further transition.
Let (X _{ t })_{ t>0} be a multistate process with a state space {1,2,…,6}. The stochastic process (X _{ t })_{ t>0} is defined as X _{ t }=ℓ, if the process is in state ℓ at time t (in months). As mentioned above, for the case study presented in this paper, there are 6 possible first lines regimens which implies that the initial state of the patient X _{0}∈{1,2,…,6}.
Our main interest is to model the transition from ℓ ^{th} regimen (state ℓ) to j ^{th} regimen (state j) at time t. A transition will be simply denoted by ℓ j. The distribution of this multistate process is characterized by the transition intensities, or hazard rate, a _{ ℓ j }(t), which expresses the instantaneous risk of a transition from state ℓ into state j at time t, that is
Here, \({F}_{t^{}}\phantom {\dot {i}\!}\) represents process history prior to time t. In our application, time t represents time since ART initiation. The cumulative transition hazards is defined as \(A_{\ell j}(t) = \int _{0}^{t} a_{\ell j}(u)du, u \leq t, \text {where}\ A_{\ell j}(t) = 0\) if a direct transition between state ℓ and j is impossible. These intensities can be gathered in to a 6×6 matrix A(t) with diagonal elements \(A_{\ell \ell }(t)=\sum _{\substack {j=1,j \neq \ell }}^{6} A_{\ell j}(t), \ell, j = \{1,2, \ldots, 6\}\). Note that individuals who have no transition should remain on their initial regimen (starting state) after ART initiation.
A central issue related to ART management is the ability to estimate the probability of the future treatment combination of the patient (i.e. the patient state) given all the information available until the present moment. For example, given a patient who substituted his NRTI’s d4T by AZT without changing his NNRTI after 6 months and (i.e., the current state of the patient is either in state 3 or state 4, depending on the initial NNRTI component) who has had no further events at one year post ART, one may be interested in estimating the probability of staying on this combination for additional 6 months as well as in comparing this probability to a patient who did not substitute their NRTI (d4T). We propose to use the transition probabilities for longterm prediction of a patient’s state. Let s be the time at which the prediction is made measured from the time origin of the patient (start of treatment) and let us denote the event history of the patient up to time s by X _{ u },0≤u≤s. Then, the transition probability from state ℓ to state j in the time interval [s,t], given information available until time s, is defined as
In order to estimate P _{ ℓ j }(s,t) we proposed to use a Markov model [28]. The model assumes that the future course of the patient depends on the patient’s state in the current time but not on the patient’s history before the current state. This means that, conditional on the present state, the past has no influence on the risk. This implies that
and
Analogous to A(t),these probabilities can be gathered in to a 6×6 matrix M(s,t) with (P _{ ℓ j }(s,t)) as its (ℓ,j)^{th} entry. A single element (P _{ ℓ j }(s,t)) combines both direct and indirect transition from state ℓ to state j. The matrix is fully presented in Additional file 1: Section S.2.
Inference
In this section we present nonparametric approaches for time continuous Markov models with finite state space and under independent right censoring. We consider n individual multistate processes \(\left (X^{(i)}_{t}\right)_{t \geq 0}, X^{(i)}_{t} \in \{1,...,6\},\ \text {and}\ i=1,2,...,n\). We assume that the n process are all, conditional on the initial state \(X^{(i)}_{0}\), independent multistate processes. Observation of the individual multistate data is subject to a right censoring time C _{ i }. Our notation and ideas are based on a counting process formulation [29, 30].
Let N _{ ℓ j;i }(t) be the counting process denoting individual i ^{′} s number of observed direct (without visiting another state in between) ℓ→j transition in [0,t],ℓ,j∈1,2,…,6,ℓ≠j. Here, time t refers to the time since the patient entered the initial state (i.e., the time since ART initiation). Let Y _{ ℓ;i }(t) be an indicator variable which represent the at risk process where we have
Let the aggregated process \(N_{\ell j}(t) = \sum _{i=1}^{n} N_{\ell j;i} (t), \ell \neq j\) and \(Y_{\ell }(t) = \sum _{i=1}^{n} Y_{\ell ;i}(t)\) denote, respectively, the number of observed direct ℓ→j transitions during the time interval [0,t] and the number of individuals to be observed at risk in state ℓ just prior to time t. We define the the increment of N _{ ℓ j }(t) as △N _{ ℓ j }(t)=N _{ ℓ j }(t)−N _{ ℓ j }(t ^{−}) for the increment of N _{ ℓ j }(t) which gives the number of ℓ→j transitions observed exactly at time t.
Nonparametric estimation of baseline hazards
From the definition of the transition intensities in Eq. (3) a _{ ℓ j } d t=P(X _{(t+d t)−}=j∣X _{ t−}=ℓ),ℓ≠j. Hence, if we observe no ℓ→j transition at t (i.e △N _{ ℓ j }(t)=0) we estimate the increment a _{ ℓ j }(t)d t of the cumulative hazard as 0. If we do observe ℓ→j transition at t (i.e △N _{ ℓ j }(t)>0), we estimate this conditional transition probability by
summing up over these increments yields the NelsonAalen estimators [29]
where summation is over all observed event times in [0,t] and its variance is given by
Nonparametric estimation of Transition probabilities
The transition probabilities are a complex function of the transition hazards, because the state occupied at some time t may potentially result from a complex nested series of competing risks experiments and there may also be more than one possible sequence of competing risks experiments leading to being in a certain state at a certain time [31]. Under the Markov model the transition probabilities defined in (4) are the solution of a set of differential equations [29]
where M(s,t) is transition probability matrix with (ℓ,j) element P _{ ℓ j }(s,t)=P(X _{ t }=j∣X _{ s }=ℓ) and A(t) is a matrix with off diagonal elements A _{ ℓ j }(t)=a _{ ℓ ℓ }(t) and diagonal elements \(A_{\ell \ell }(t)=\sum _{\substack {j=1,j \neq \ell }}^{6} a_{\ell j}(t)\). In coordinates, (8) is \(d/dt P_{\ell j}(s,t) =\sum _{k} P_{\ell k}(s,t)A_{kj}(t)\) and for any fixed initial state ℓ, the vector of transition probabilities from ℓ, (P _{ ℓ1}(s,t),P _{ ℓ2}(s,t),...,P _{ ℓ6}(s,t)) satisfies this equation. Even though this equation can not be solved in general due to the nonconstancy over time of the matrix A(t), under Markov assumption, the transition probabilities satisfy
and based on a partition s=t _{0}<t _{1}<t _{2}<…<t _{ k−1}<t _{ k }=t of the time interval [s,t], the matrix of transition probabilities M(s,t) can be approximated by [31]
where I is the (6×6) identity matrix, the (ℓ,j)^{th} element of △A(t _{ k }) is equal to A _{ ℓ j }(t _{ k })−A _{ ℓ j }(t _{ k−1}), and \(A_{\ell \ell }(t)=\sum _{\substack {j=1,j \neq l }}^{6} A_{lj}(t)\). Computing the approximation for ever finer partition [s,t] approaches a limit, namely a matrixvalued product integral [31, 32], which equals the matrix of transition probabilities,
where u ranges from s to t and d A(u) is defined as d A _{ ℓ j }(u)=a _{ ℓ j }(u)d u,ℓ,j∈1,2,...,6 [29]. Therefore, for Markov models, given A(t), the product integration is the mapping that switches from cumulative transition hazards to the matrix of transition probabilities while all cumulative transition hazards are involved.
An estimator of M(s,t) can be obtained by replacing A(u) with the NelsonAalen estimators, \(\hat {\mathbf {A}}(u)\), and by defining \(d\hat {\mathbf {A}}(u)\) as the matrix with entries \(\triangle \hat {A}(u)= \hat {A}_{\ell j}(u)  \hat {A}_{\ell j}(u^{})\) (i.e., the increment of the NelsonAalen estimators at time u). This results in the AalenJohansen type estimator [29],
in which u indicates all event times in (s,t]. Note that the transition probability matrix defined in (11) is calculated by means of a product integral, while its estimator in (12) is based on a finite product, which only changes at event times.
The transition probabilities can be used for two types of prediction: forward and fixed horizon [30, 33]. In the former case, at a given fixed time s the probabilities of possible future events are evaluated for varying time horizons t. In the latter case, the prediction is made from several starting points to one future fixed time point. In both cases, Aalentype or Greenwood type estimators of the variancecovariance matrix of M(s,t) can be calculated directly or through a recursion formula which can for instance be used to construct pointwise confidence intervals around the estimated transition probability curves [30]. In our application we use forward prediction type.
Robustness of Firstline HAART towards Treatment Modification/change
The primary aim of this study is to quantify the robustness of first line treatments to treatment modification. Given the individuals initial state ℓ at time s, the waiting time in state ℓ. i.e., the duration of stay at state ℓ, can be used as a summary measure of the model. The waiting time in state ℓ is generated with hazard
We define the total cumulative transition hazard out of state ℓ as \( A_{\ell.}(t)=\int _{0}^{t} a_{\ell.}(u)du=\sum _{j=1,j \neq \ell }^{6} A_{\ell j}(t)\). Using these quantities one can evaluate the probability of no events during a period. The survival function of the waiting time in the initial state ℓ, i.e., probability to stay in state ℓ until time t, given that the individual had already been in the respective state at time s,s≤t is given by,
These quantities are essentially common survival probabilities with all cause hazard a _{ ℓ.}(u), taking time s as time origin [31]. However, this can also be seen as a solution of the product integral in (11). Since, the ℓ ^{th} row of the AalenJohansen type estimator of \(\hat {\textbf {M}}(s,t)\) contains the estimates \(\hat {P}_{\ell j}(s, t)\) for ℓ≠j and the diagonal element is such that the sum over the ℓ ^{th} row equals 1, the AalenJohansen estimator of P(X _{ t }=ℓx _{ s }=ℓ) is just \(\hat {P}_{\ell \ell }(s,t)\).
The multistate model formulated above allows us to evaluate whether treatment modification reflect a substitution of NNRTI, substitution of NRTI or substitution of both NNRTI and NRTI by initial treatment combinations. We propose to use the following measures of HAART robustness to treatment modification:

1.
Probability of NNRTI substitution
P _{12} for state 1 P _{34} for state 3 P _{56} for state 5
P _{21} for state 2 P _{43} for state 4 P _{65} for state 6

2.
Probability of NRTI substitution
P _{13}+P _{16} for state 1 P _{35}+P _{36} for state 3 P _{54}+P _{52} for state 5
P _{24}+P _{25} for state 2 P _{42}+P _{45} for state 4 P _{63}+P _{61} for state 6

3.
Probability of regimen switching
P _{14}+P _{15} for state 1 P _{32}+P _{35} for state 3 P _{53}+P _{51} for state 5
P _{23}+P _{26} for state 2 P _{41}+P _{46} for state 4 P _{64}+P _{62} for state 6
Result
Of the 1453 eligible patients, 169 patients were excluded because of limited follow up (i.e those with at most 1 month followup data) and missing information (patients with missing information about prescribed ARV or start and stop dates of the drug). A total of 1284 subjects were included for the analysis presented in this paper. Patients persontime were cut at the earliest of changing to second line treatment, death, lost to followup, transfer or end of study (Aug 25, 2013). The median followup time was 37.40 months (IQR: 22.3256.15 months) and the average followup time was 38.25 months per person. At ART initiation, patients had a median CD4 cell count of 137c e l l s/m m ^{3} (IQR: 78201 c e l l s/m m ^{3}), were predominately female (68.81%) and had a median age of 30 years (IQR: 2635 years) (Table 1). The most common regimens initiated were d4T + 3TC + EFV consisting of 526(40.96%) patients while the treatment TDF + 3TC + EFV was administered to 401(31.23%) patients at initiation.
For the majority of the patients (58.8%) the first line treatment was not modified, 442 patients (34.4%) had their treatment changed only once while 86 patients (6.69%) had their treatment changed more than once. In total, 615 (32.4%) treatment changes occurred in this cohort over the period of followup. Of those 615 changes, 426 (69.27%) substituted NRTI only, 144 (23.41%) substituted NNRTI only and 45 (7.32%) substituted both the NRTI and NNRTI at the same time. The number of events (transition made from each state) and the number of patients in total that were at risk for treatment modification are presented in Table 2. Five hundred forty one patients were on d4T + 3TC + NVP combination, but only 125 (23%) remained on this regimen without any modification. 89% of the 483 patients who were on AZT + 3TC + NVP and 95% of the 89 patients who were on TDF + 3TC + NVP did not experience treatment modification. Among the regimens containing EFV, 36% of 157 patients on d4T + 3TC + EFV, 86% of 161 patients on AZT + 3TC + EFV and 96% of 471 patients on TDF + 3TC + EFV did not experience treatment modification. The frequency of treatment change was the lowest amongst those patients initiated on TDFbased regimens (3.75%) compared to those initiated on AZT (11.96%) and d4Tbased regimens (73.92%). It is interesting that regimens containing d4T were more prone to treatment modification than those containing AZT and TDF. Apart from d4T, patients who received NVP (42.5%) were more susceptible to treatment modification than patients who received EFV (17.87%).
As seen from Fig. 2 (Panel b), however, when we look at the time spent in the current treatment combination of the patients who modified their treatment, patients initiated on d4T had a tendency to stay longer (40.40 months; IQR: 14.6055.73) as compared to AZT patients (3.7 months; IQR: 1.9016.02) and TDF patients (12.60 months; IQR: 6.8420.40). Similarly, patients initiated on NVP had a tendency to stay longer (38.37 months; IQR: 5.4256.05) as compared to EFV patients (21.80 months; IQR: 8.7039.80) (Fig. 2 (Panel c)). The duration of stay in each treatment combination before the first change to another treatment combination is presented in Additional file 1: Figure S3.1.
The multistate model described in the previous section was estimated using the mstate packagee developed by [34]. Details about different R function used for the estimation is given in Additional file 1: Section S1. In particular, using the estimated transition hazards as described in Eq. (12), we calculate the transition probabilities P _{ ℓ j }(s,t) from all starting states to all possible states, between the starting time s=0 and all event times successively. Note that several probabilities estimates cannot be obtained due to limited information in some states. As shown in Table 2, treatment change was observed only in 4 out of the 89 patients initiated on TDF + 3TC + NVP; hence we have chosen to consider as inadmissible the occurrence treatment modification from this treatment combination (State 6). The model has 6states as before but with a different transition matrix. The transition matrix in Additional file 1: Section S1 shows the multistate structure which reflects this framework. We show in Fig. 3 the estimated transition probabilities from all starting states to all possible states, between the starting time s=0 and all event times successively. Treatment combinations containing d4T have the lowest probability of no treatment modification while the combination of TDF and EFV are the most robust to treatment modification.
As mentioned above, we are mainly interested in prediction of the four measures of HAART robustness to treatment modification: (1), probability of no treatment modification, (2) NRTI substitution, (3) NNRTI substitution and (4) regimen changes. The estimated probabilities are shown in Fig. 4. As was to be expected on the basis of the previous discussion, the prospects for a patient who received the regimens containing d4T are indeed worse than those patient who received the regimens containing AZT or TDF, the former having a far larger probability of NRTI substitution, regimen changes and the lowest treatment modificationfree survival probabilities.
As mentioned above, treatment modification occurs more frequently for AZT and TDF early after treatment initiation while it occurs later on in followup among patients on d4T. Thus, it is interesting to compare treatments based on the situation after some months to account for early ART complication. For this, the transition probabilities at 10 months post ART were estimated and the results are presented in Fig. 5. A comparison of Figs. 4 and 5 clearly shows that the fact that a patient on a treatment combination containing AZT or TDF has not had any early ART complication leading to treatment change or modification in the first 10 months postART has decreased his/her probability of future treatment change or modification considerably; notably, his/her probability of longterm NRTI substitutionfree survival has increased significantly. On the other hand, the longterm treatment modificationfree survival of a patient on a treatment combination containing d4T was unchanged by the fact that he/she has not experienced treatment modification in the first 10 months postART.
Reason for treatment modification
Table 3 shows the reason for treatment change for the total 615 observed treatment changes in the cohort, stratified by treatment combinations. We were able to obtain the reason for the majority of treatment changes (88.62%). Toxicity and comorbidity were was the main reasons for treatment modification accounting for 48.94% and 14.31% of the observed treatment changes, respectively. About 50% of the patients on all the regimens except TDF + 3TC + EFV reported toxicity or side effects. In addition, phasing out of d4T from the NRTI backbone accounts for 20.16% of the observed treatment changes.
In order to quantify the effect of toxicity on treatment change, we modify the definition of timetotreatment change to timetotreatment change due to toxicity. Here, treatment change related to other reasons during the followup period were censored at the time of their occurrence. As seen from Table 4, a large proportion of patients (26.25%) on NVP and d4T combination had NRTI substitution with d4T replaced by AZT due to toxicity. Similarly, 15.92% and 14.01% of patients on EFV combination with d4T had d4T replaced by AZT and TDF, respectively due to toxicity. Treatment combination of TDF and EFV was the most robust to treatment modification among the six first line regimens. As previously noted, this combination seems the least toxic. In a similar fashion, we calculate the toxicity driven transition probabilities P _{ ℓ j }(s,t) from all starting states to all possible states, between the starting time s=0 and all event times successively and extract the four measures of transition probabilities. Here again, Fig. 6 shows treatment combination containing d4T had highest probability for treatment modification due to toxicity. Where as the combination of TDF and EFV was the most robust to treatment modification.
Discussion
This study has provided unique and important data on durability of firstline ART and on reasons responsible for antiretroviral treatment modification in the setting of a tertiary care Hospital in a resource limited country. This work adds to the previous observational studies [15–18, 22–27] conducted in resourcelimited settings, where no distinction was made between NRTI substitution and NNRTI substitution, treatment modification due to all causes and toxicity driven treatment modification, and the incidence of subsequent treatment modification was not studied. Further, we show how a simple multistate survival model can be used to estimate the probability of the future treatment combination of the patient given all the information available up to the present moment.
In our cohort, a large proportion of patients (41.2%) changed their treatment during follow up time, where 34.4% patients changed once and 6.69% changed more than once, which is far higher to what has been previously reported [22, 27]. The short follow up time and consideration of only the first treatment modification as event of interest in their study may be the reason for such discrepancy. Among the regimens containing d4T, 77% of 541 patients on d4T + 3TC + NVP and 64% of 157 patients on d4T + 3TC + EFV experienced all cause treatment modification. In all cause analysis, regimens containing d4T had highest probability for treatment modification, NRTI substitution, and regimen switching as compared to those regimens containing AZT and TDF, consistent with previous findings [15–18, 22–27]. Whereas the combination of TDF and EFV was the most robust to treatment modification. Apart from d4T, patients on EFV were less susceptible to treatment modification than patients on NVP, similar to what has been reported previously [27].
We also found that treatment modification occurring more frequently for AZT and TDF early after treatment initiation while treatment modification occurs later on in followup amongst patients on d4T, consistent with previous findings [27]. The superiority of AZT and TDF over d4T, however, should not be shadowed by this finding. A further comparison of treatment combinations accounting for early ART complication shows that, if a patient on a treatment combination containing AZT or TDF has not had any early ART complication leading to modification in the first 10 months postART his probability of future treatment modification decreased considerably. On the other hand, the longterm treatment modificationfree survival of a patient on a treatment combination containing d4T was unchanged by the fact that he/she has not experienced treatment modification in the first 10 months postART. Further, no significant difference in the timing of treatment modification was observed among NVP and EFV.
The unique feature of this study is we manage to determine the type of treatment modification along with the reason for modification. Only in less than 7% of the treatment changes were we unable to determine the reason for treatment modification. Toxicityrelated treatment modification has been identified as the most common reason for treatment modification accounting for 48.94% of the changes, followed by comorbidity (New TB) 14.31%, similar to what has been reported previously [15–18, 22–27]. About 50% of the patients on all the regimens except TDF plus EFV reported toxicity or side effects. The largest number of treatment modification due to toxicity was in patients on d4T: approximately 27% of those who originally started with NVP and d4T combination had NRTI substitution with d4T replaced by AZT and 15.92% and 14.01% patients who originally started with EFV combination with d4T had d4T replaced by AZT and TDF,respectively. Treatment combination of TDF and EFV was the most robust to toxicity related treatment modification among the six first line regimens. As previously noted, this combination seems the least toxic. This is a significant finding because TDF is a WHO recommended preferred treatment, with AZT as alternative [20]. Phasing out of d4T also accounted for 20.16% of treatment changes observed in our study.
This study provided unique and important data on durability of firstline ART and on reasons responsible for antiretroviral treatment modification in a resource limited setting. Furthermore, the study shows the use of multistate models to study the evolution of patient’s state (treatment regimen) over time and to predict the probability of changing treatment. The proposed model allow us to model both the occurrence of different event types (such as, single drug substitution or regimen switch) and the occurrence of subsequent events, the latter potentially of different types in a unified way.
Our findings must be interpreted in light of some limitations. Our model assumes that the future course of a patient only depends on where you are at the current time, but not on how you got there. Deviations from this could have led to bias.
Conclusion
Our study shows the burden of toxicity/side effect related to d4T use is a matter of major concern, as it accounts for the majority of modifications. Safer and more tolerable regimens like a combination of TDF and EFV should be made more accessible to treatment programs in resourcelimited settings. Moving away from drugs with poor safety profiles, such as d4T, could reduce modification rates and this would improve regimen tolerability, while preserving future treatment options.
Abbreviations
 ART:

Antiretroviral therapy
 HAART:

Highly active antiretroviral therapy
 NRTI:

Nucleotide reverse transcriptase inhibitors
 NNRTI:

Non nucleotide reverse transcriptase inhibitors
 WHO:

World Health Organization
 VCT:

Voluntary counselling and testing
 PMTCT:

Prevention of mothertochild transmission
References
 1
Wolbers M, Battegay M, Hirschel B, Furrer H, Cavassini M, Hasse B, et al.CD4 ̂+ Tcell count increase in HIV1infected patients with suppressed viral load within 1 year after start of antiretroviral therapy. Antivir Ther. 2007; 12(6):889.
 2
Moore RD, Keruly JC. CD4+ cell count 6 years after commencement of highly active antiretroviral therapy in persons with sustained virologic suppression. Clin Infect Dis. 2007; 44(3):441–6.
 3
Mocroft A, Phillips AN, Gatell J, Ledergerber B, Fisher M, Clumeck N et al. Normalisation of CD4 counts in patients with HIV1 infection and maximum virological suppression who are taking combination antiretroviral therapy: an observational cohort study. The Lancet. 2007; 370(9585):407–13.
 4
Lawn SD, Myer L, Bekker LG, Wood R. CD4 cell count recovery among HIVinfected patients with very advanced immunodeficiency commencing antiretroviral treatment in subSaharan Africa. BMC Infect Dis. 2006; 6(1):1.
 5
Erhabor O, Ejele O, Nwauche C. The effects of highly active antiretroviral therapy (HAART) of stavudine lamivudine and nevirapine on the CD4 lymphocyte count of HIVinfected Africans. The Nigerians experience. Nigerian journal of clinical practice. 2006; 9(2):128–33.
 6
Seid A, Getie M, Birlie B, Getachew Y. Joint modeling of longitudinal CD4 cell counts and timetodefault from HAART treatment: a comparison of separate and joint models. Electronic J Appl Stat Anal. 2014; 7(2):292–314.
 7
Calmy A, Pinoges L, Szumilin E, Zachariah R, Ford N, Ferradini L, et al.Generic fixeddose combination antiretroviral treatment in resourcepoor settings: multicentric observational cohort. Aids. 2006; 20(8):1163–9.
 8
Van Roon E, Verzijl J, Juttmann J, Lenderink A, Blans M, Egberts A. Incidence of discontinuation of highly active antiretroviral combination therapy (HAART) and its determinants. JAIDS J Acquir Immune Defic Syndr. 1999; 20(3):290–4.
 9
Monforte Ad, Lepri AC, Rezza G, Pezzotti P, Antinori A, Phillips AN, et al.Insights into the reasons for discontinuation of the first highly active antiretroviral therapy (HAART) regimen in a cohort of antiretroviral naive patients. Aids. 2000; 14(5):499–507.
 10
Fellay J, Ledergerber B, Bernasconi E, Furrer H, Battegay M, Hirschel B, et al. Prevalence of adverse events associated with potent antiretroviral treatment: Swiss HIV Cohort Study. The Lancet. 2001; 358(9290):1322–7.
 11
Mocroft A, Youle M, Moore A, Sabin CA, Madge S, Lepri AC, et al. Reasons for modification and discontinuation of antiretrovirals: results from a single treatment centre. Aids. 2001; 15(2):185–94.
 12
Mocroft A, Phillips A, Soriano V, Rockstroh J, Blaxhult A, Katlama C, et al.Reasons for stopping antiretrovirals used in an initial highly active antiretroviral regimen: increased incidence of stopping due to toxicity or patient/physician choice in patients with hepatitis C coinfection. AIDS Res Hum Retrovir. 2005; 21(9):743–52.
 13
Cardoso SW, Grinsztejn B, Velasque L, Veloso VG, Luz PM, Friedman RK, et al.Incidence of modifying or discontinuing first HAART regimen and its determinants in a cohort of HIVinfected patients from Rio de Janeiro, Brazil. AIDS Res Hum Retrovir. 2010; 26(8):865–74.
 14
Teklay G, Legesse B, Legesse M. Adverse effects and regimen switch among patients on antiretroviral treatment in a resource limited setting in Ethiopia. J Pharmacovigilance. 2014;2013.
 15
Woldemedhin B, Wabe NT, et al. The reason for regimen change among HIV/AIDS patients initiated on first line highly active antiretroviral therapy in Southern Ethiopia. N Am J Med Sci. 2012; 4(1):19.
 16
Assefa D, Hussein N. Reasons for Regimen Change among HIV/AIDS Patients Initiated on First Line Highly Active Antiretroviral Therapy in Fitche Hospital, Oromia, Ethiopia. Adv Pharmacol Pharm. 2014; 2(5):77–83.
 17
Jima YT, Angamo MT, Wabe NT. Causes for antiretroviral regimen change among HIV/AIDS patients in Addis Ababa, Ethiopia. Tanzania J Health Res. 2013;15(1).
 18
Wube M, Tesfaye A, Hawaze S. Antiretroviral therapy regimen change among HIV/AIDS patients in Nekemt Hospital: a primary care Hospital in Oromia Regional State, Ethiopia. J Appl Pharm Sci. 2013; 3(8):36.
 19
O’Brien ME, Clark RA, Besch CL, Myers L, Kissinger P. Patterns and correlates of discontinuation of the initial HAART regimen in an urban outpatient cohort. JAIDS J Acquir Immune Defic Syndr. 2003; 34(4):407–14.
 20
World Health Organization. Antiretroviral therapy for HIV infection in adults and adolescents: recommendations for a public health approach2010 revision: Geneva: World Health Organization; 2010.
 21
World Health Organization and others. Consolidated guidelines on the use of antiretroviral drugs for treating and preventing HIV infection: recommendations for a public health approach: World Health Organization; 2016.
 22
Louwagie G, Zuma K, Okello V, Takuva S. Durability of first line antiretroviral therapy: reasons and predictive factors for modifications in a Swaziland cohort. J Antivirals Antiretrovirals. 2012;2012.
 23
Velen K, Lewis JJ, Charalambous S, Grant AD, Churchyard GJ, Hoffmann CJ. Comparison of tenofovir, zidovudine, or stavudine as part of firstline antiretroviral therapy in a resourcelimitedsetting: a cohort study. PLoS ONE. 2013; 8(5):e64459.
 24
Chi BH, Mwango A, Giganti M, Mulenga LB, TambatambaChapula B, Reid SE, et al. Early clinical and programmatic outcomes with tenofovirbased antiretroviral therapy in Zambia. J Acquir Immune Defic Syndr. 2010; 54(1):63.
 25
Bygrave H, Ford N, van Cutsem G, Hilderbrand K, Jouquet G, Goemaere E, et al. Implementing a tenofovirbased firstline regimen in rural Lesotho: clinical outcomes and toxicities after two years. JAIDS J Acquir Immune Defic Syndr. 2011; 56(3):e75–8.
 26
Njuguna C, Orrell C, Kaplan R, Bekker LG, Wood R, Lawn SD. Rates of switching antiretroviral drugs in a primary care service in South Africa before and after introduction of tenofovir. PLoS One. 2013; 8(5):e63596.
 27
Brennan AT, Maskew M, Ive P, Shearer K, Long L, Sanne I, et al. Increases in regimen durability associated with the introduction of tenofovir at a large publicsector clinic in Johannesburg, South Africa. J Intl AIDS Soc. 2013;16(1).
 28
Hougaard P. Analysis of multivariate survival data: Springer Science & Business Media; 2012.
 29
Andersen PK, Borgan O, Gill RD, Keiding N. Statistical models based on counting processes: Springer Science & Business Media; 2012.
 30
De Wreede LC, Fiocco M, Putter H. The mstate package for estimation and prediction in nonand semiparametric multistate and competing risks models. Comput Methods Prog Biomed. 2010; 99(3):261–74.
 31
Beyersmann J, Latouche A, Buchholz A, Schumacher M. Simulating competing risks data in survival analysis. Stat Med. 2009; 28(6):956–71.
 32
Gill RD, Johansen S. A survey of productintegration with a view toward application in survival analysis, The annals of statistics: JSTOR; 1990. pp. 1501–55.
 33
Putter H, van der Hage J, de Bock GH, Elgalta R, van de Velde CJ. Estimation and prediction in a multistate model for breast cancer. Biom J. 2006; 48(3):366–80.
 34
de Wreede LC, Fiocco M, Putter H, et al. mstate: an R package for the analysis of competing risks and multistate models. J Stat Softw. 2011; 38(7):1–30.
Acknowledgements
The authors are grateful to Jimma University Specialized Hospital ART clinic for the permission to use the data and the staff members of the clinic for their support in extracting the information from patients medical card. Financial support to the first authors for his research visit from the Institutional University Cooperation of the Council of Flemish Universities (VLIRIUC) is gratefully acknowledged. The authors are grateful to Kibralem Sisay for his support in the programming of the R function we used for data preparation.
Funding
This work was financially supported by Jimma University Interuniversity cooperation (IUCJU)
Availability of data and materials
The data sets analyzed in this study available from the corresponding author on reasonable request. The R code used to analyze the data provided as a supplement of the article.
Authors’ contributions
BB contributed to the study concept and design, performed the analysis on the data set as well as wrote the first draft of the paper. TA contributed to the analysis and interpretation of the data, in addition to drafting and critical revision of the manuscript. RB and ZS contributed to the study concept and design, the statistical methodology and finalization of the writing. AS contributed to the analysis and interpretation of the Application section and critical revision of the manuscript. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Human subject research approval for this study was received from Jimma University Research Ethics Committee and the medical director of the Hospital. As the study was retrospective, informed consent was not obtained from the study participants, but data were anonymous and kept confidential.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Author information
Additional file
Additional file 1
Supplementary Appendix: Evaluation of the Durability of Firstline Highly Active Antiretroviral Therapy in Southwest Ethiopia Using Multistate Survival Model. (PDF 218 KB)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 HIV/AIDS
 Highly active antiretroviral therapy
 Treatment modification
 Survial analysis
 Multistate models