Modeling the Effects of Vorinostat Reveals both Transient and Delayed HIV Transcriptional Activation and Minimal Killing of Latently Infected Cells

Summary

Combination antiretroviral therapy (cART) for HIV infection must be taken for life due to the existence of long lived latently infected cells. Recent efforts have focused on developing latency reversing agents to eliminate latently infected cells by activating HIV production. In this work, we assess the impact of a latency reversing agent, vorinostat, by fitting dynamic models to data from a clinical trial. Results show that vorinostat treatment induces HIV transcription transiently and that the sustained induction of HIV transcription may depend on the temporal impact of vorinostat on host gene expression. Our results also suggest that vorinostat treatment is not sufficient to induce killing of latently infected cells in a majority of HIV-infected individuals on cART.

Introduction

Treatment of HIV-infected individuals with combination antiretroviral therapy (cART) effectively suppresses HIV to levels below the limit of detection of conventional assays and substantially reduces morbidity and mortality of HIV infected patients [1]. However, it does not eradicate the virus and treatment is lifelong [2]. Therefore, developing novel therapeutics to cure HIV infection remains an important research priority [3,4]. A major barrier to cure is the presence of a population of long lived latently infected cells [4] that can persist indefinitely in patients treated with highly potent cART [5]. Recent efforts have focused on strategies that activate HIV production in latently infected cells. The idea, termed ‘shock and kill’ [6], is to first shock latently infected cells thereby activating HIV gene expression, such that the cells are then killed by viral cytopathic effects or immune-mediated cell death. Histone acetylation is one of several factors that regulate HIV transcription and is therefore important for establishing and maintaining latency [7]. Drugs such as histone deacetylase inhibitors (HDACi) enhance acetylation of both histones and proteins and thereby induce changes in gene transcription, including transcription of HIV [8].

Vorinostat, a histone deacetylase inhibitor licensed for the treatment of cutaneous T-cell lymphoma [9], has been shown to activate HIV transcription in resting memory CD4+ T-cells in vivo [10,11]. In a recent clinical trial, 20 HIV-1 infected individuals on suppressive cART were treated orally with 400 mg a day of vorinostat for 14 days and then followed for an additional 70 days. Overall, vorinostat induced a rapid and sustained increase of cell-associated unspliced (CA-US) HIV RNA [10]. However, the response pattern was highly variable among the participants. For example, in half of the participants, after an initial significant increase in CA-US HIV RNA, the level of CA-US HIV RNA decreased rapidly within 1–3 days before increasing again, and in 14 of the 20 participants, the level of CA-US HIV RNA continued to increase after vorinostat was stopped. These puzzling observations raise important questions about the temporal impact of vorinostat treatment on HIV transcription and the design of treatment strategies to eradicate the latent reservoir.

We constructed mathematical models to better understand the temporal changes in CA-US HIV RNA in individuals treated with vorinostat. Mathematical models have been widely applied to study viral dynamics in vivo [12–14]. They played an instrumental role in quantifying important parameters, such as the half-lives of virions and infected cells in vivo [14]. Recently, several models have been developed to understand the maintenance of the latent reservoir under cART treatment [15–17], the viral rebound time distribution after latency reversing agent (LRA) treatment [18], and the optimal time to start a LRA [19]. However, the dynamic response of HIV transcription in latently infected cells following treatment with a LRA has not been investigated. This question has important implications for future clinical trial design and optimizing treatment strategies to eliminate latently infected cells. Previous models have generally assumed that the HIV provirus in latently infected cells becomes fully activated following treatment with a LRA and that subsequent events will be identical to latently infected cells activated by normal immunological signals or through the T-cell receptor [18,19]. However, evidence suggests that current LRA treatments primarily activate HIV transcription and its impact on translation may be mild or minimal [20,21]. Here, we construct models that treat cells activated by a LRA and naturally activated cells separately. By fitting models to the clinical data, we show the complex dynamic response of latently infected cells to vorinostat can be explained. Furthermore, we use the models to quantify the extent to which vorinostat activates HIV transcription and induces cell death in vivo.

Results

Direct Activation Model

We first construct a mathematical model based on a previously published latency model by Rong et al. [15]. As the study participants were treated with suppressive cART for a medium of 5 years, we have chosen parameters in the Rong et al. model such that before vorinostat treatment the viral load is approximately 5 HIV RNA copies/ml, target cells levels are 750 T cells/μL similar to mean CD4 count in the 20 clinical trial subjects [10] and latently infected cell levels are 2/mL, which is approximately 2.7/million CD4 cells (roughly consistent with previous studies [22,23]). The major innovation in the direct activation model is that we assume latently infected cells become transcriptionally induced and express CA-US HIV RNA directly upon vorinostat treatment (Fig 1A; see Methods for full description of the model and Table 1 for parameter values). Thus, two equations are added to the Rong et al. model [15], one for the number of cells that have HIV transcription induced by vorinostat, L_A, and one for the average amount of CA-US HIV RNA per transcriptionally activated cell, R. Also, cells that are transcriptionally induced by vorinostat are assumed to be in a different state than cells that are naturally activated.

We fitted this model to the clinical data collected during the entire 84-day study period (see Methods for the fitting procedure and S1 Table for best-fit parameter values). In general, the direct activation model does not explain the data well especially during the first 1–3 days’ treatment (S1 Fig) and the period after treatment stops (S2 Fig). First, in 10 of the 20 participants, the level of CA-US HIV RNA first increased upon initiation of vorinostat, and then decreased rapidly after the first 1–3 days, whereas the best-fit model curves in these patients have CA-US HIV RNA increasing continuously during vorinostat treatment. Second, in 14 out of the 20 participants, the level of CA-US HIV RNA increased at variable time points after cessation of vorinostat at day 14, whereas the model predicts that the level always decreases over time after cessation of vorinostat.

To better understand the initial peaking pattern observed in half of the study participants soon after the initiation of vorinostat, we fitted the model to data obtained during the first 7 days of treatment only (Fig 2). The peaking pattern was well described by the model in five participants (VOR001, VOR004, VOR010, VOR016, VOR018). In these individuals, the estimated loss rate of transcriptionally activated cells, d_LA, ranged between 1.8 and 10 day^-1 (see S2 Table), and is much greater than the death rate of productively infected cells, i.e. 1.0 day^-1 [27]. Rates of decrease of CA-US HIV RNA levels higher than 1.0 day^-1 were also apparent in 5 other participants (VOR008, VOR019 and VOR021-023). This decrease can potentially result from either death of cells transcriptionally activated by vorinostat or from shutdown of HIV transcription and loss of CA-US HIV RNA. We reason that this decline is not due solely to cell death, as it is unlikely that cells activated by vorinostat die at a faster rate than productively infected cells. Then why would HIV transcription shut down during vorinostat treatment?

Previous in vitro studies have shown the activation of HIV transcription is a transient stochastic process, and that the duration of this transient process is dependent on the strength of Tat transcriptional feedback [34–36], as well as the availability and regulation of many host factors that are necessary for transcriptional activation, such as the NAD-dependent deacetylase sirtuin-1, NF-κB, Yin Yang 1 and the positive transcription elongation factor, P-TEFb [37–41]. In latently infected cells, mostly memory T cells, these transcription factors are likely to be at low levels [42–44], whereas many host enzymes such as Murr1 (a gene product that restricts HIV-1 replication), human schlafen 11 and the lipid raft associated protein tetherin, actively inhibit HIV transcription initiation [45], mRNA translation [46], and viral release [47]. Therefore, before vorinostat treatment, the host factors/enzymes required for full HIV gene activation are most likely limiting in latently infected cells. After vorinostat treatment initiation, host genes undergo rapid differential regulation at 2, 8 and 24 hours [10]; however, the immediate impact of vorinostat treatment may not be sufficient to induce HIV gene transcription sustainably. This unfavorable cellular environment and rapid changes in gene expression may lead to very short transcriptional pulses of Tat activity and CA-US HIV RNA production. Without further production of CA-US HIV RNA, the rapid decrease observed in the data may be a result of the loss of US HIV RNA by degradation and by splicing.

Delayed Activation Model

Vorinostat treatment not only induces rapid changes in host gene expression but also induces changes after treatment cessation [10]. It is therefore plausible that the late increase in HIV transcription after vorinostat treatment is due to a longer-term impact on host gene transcription. To test this hypothesis, we extended the direct activation model to include a ‘transiently activated’ state (L_T) and a ‘waiting’ state (L_W), and denote this as the ‘delayed activation’ model. We assume that upon vorinostat treatment latently infected cells first get activated transiently, i.e., enter the transient activation state, L_T, where CA-US HIV RNAs are produced for a short period of time. We assume the cells then enter a waiting state, L_W, in which there is no CA-US HIV RNA production before becoming sustainably activated cells, L_A (Fig 1B). The waiting state reflects the time needed for the transcriptional programs to produce sufficiently high levels of host factors necessary for transcriptional activation such that the cellular environment becomes favorable for sustained HIV transcription. We assume that the cells in the transiently activated state and in the waiting state die at the same rate as in the latent state (L) as transient activation is not likely to be strong enough to produce the shock needed for kill. The ordinary differential equations (ODEs) describing this model are given in the Methods.

Next, we tested whether the above hypotheses explain both the short-term and the long-term dynamics of CA-US HIV RNAs by fitting the delayed activation model to the full data set (see S3 Table for best-fit parameter values), and found that the delayed activation model describes the data much better than the direct activation model in a majority of participants (compare S3 and S4 Figs with S1 and S2 Figs, respectively). It successfully describes the initial pattern of CA-US HIV RNA change following initiation of vorinostat in most individuals as well as the dynamics of CA-US HIV RNA in 6 of the 14 individuals where the level of CA-US HIV RNA increased after cessation of vorinostat. For the other 8 patients, the delayed activation model does not predict the magnitude of the late increase in CA-US HIV RNA level at some time points (S4 Fig). We speculate that the discrepancy may arise from the assumption of an exponentially distributed residence time for latently infected cells in the waiting state before becoming sustainably activated (an assumption implicitly assumed in the ODE system). This assumption is valid when a single event is needed for the transition to sustained activation. However, it is likely that multiple events must occur before the transition to sustained activation, such as upregulation of several host factors, HIV RNA splicing and expression of tat and other regulatory proteins including rev.

Multistage Delayed Activation Model

We, thus, further modified our model to assume that cells in the waiting state have to go through several stages before becoming sustainably activated as in previous work describing the multiple events needed to drive an initially infected cell into viral production [48]. We denote this model the ‘multistage delayed activation’ model and the equations describing this model are given in the Methods. See Table 2 for a summary of the assumptions made with regard to the impact of vorinostat on latently infected cells in the three different models.

In this model, the L_W state is divided into n identical sub-states, i.e. L_W,1, L_W,2, … L_W,n. The transition rate from one sub-state to the next is set to nk_w such that the average residence time in the overall waiting state is 1/k_w. This model is equivalent to one in which we assume the transition out of the waiting state is stochastic with the delay described by a gamma probability distribution [49]. We let n change from 1 to 10, and fitted these 10 model variants to the clinical data from all 20 participants (S5 Fig). The fitting results show that this multistage delayed activation model describes the patterns of increases of CA-US HIV RNA after cessation of vorinostat as well as the initial peak following initiation of vorinostat (Fig 3 and S6 Fig; see S4 Table for best-fit parameter values).

We further performed model selection using the corrected Akaike information criterion (AICc) (see Methods). The direct activation model significantly underperformed compared to the delayed activation model and the multistage delayed activation model in 19 out of the 20 participants. The multistage delayed activation model was significantly better than the delayed activation model in 12 participants (Table 3). These results support the hypothesis that the immediate impact of vorinostat treatment is to activate HIV transcription for a short period of time (1–3 days), possibly due to the limited availability of many host cellular factors and that sustained activation may take longer to attain and a number of events (possibly in host cell transcriptional regulation) must occur before sustained HIV transcription becomes possible. As the number of events required varied among the participants, cells in different individuals may be in different states of latency. Also, host gene expression patterns, which can differ among individuals, may play a role in determining the length of the delay. Analyzing the effect of changing the number of waiting stages on the model fit to the data using AICc shows that for 17 out of the 20 participants, using a model with more than 7 stages would be a good choice in general (S5 Fig).

Quantitative Impact of Vorinostat Treatment on Latently Infected Cells

We next examined the best-fit parameter values of the multistage delayed activation model to assess the impact of vorinostat on latently infected cells in vivo. First, we find that the estimated values of α, the rate of CA-US HIV RNA production induced by vorinostat in latently infected cells, varies over a wide range (over 1.5 logs) among the 20 patients, suggesting the response to vorinostat is very heterogeneous across participants. In a majority of patients, the estimated values of α are smaller than the production rate of CA-US HIV RNA in productively infected cells, α_I (Fig 4A; see Methods for calculation of α_I). We then examined the estimated loss rate of transcriptionally activated latently infected cells, d_LA, and found that in 12 participants the estimated loss rates are extremely low, close to the death rate of latently infected cells, d_L (Fig 4B). Although the estimated loss rates are higher than d_L in other patients, we were not able to distinguish whether the loss is through shutdown of HIV transcription or through cell death. Nonetheless, the low estimates of the loss rate in most participants suggest that vorinostat treatment does not induce killing of transcriptionally activated latent cells in vivo in a majority of individuals, and thus, according to this model, the reductions in reservoir size were minimal or absent in most participants.

We further tested the robustness of the parameter estimates to variations in our assumptions. First, we varied the values of two fixed parameter values that describe the intracellular dynamics of CA-US HIV RNAs, i.e. the rate of US HIV RNA export from the cell in the form of virions, ρ, which we initially assumed to be 0, and the combined rate of RNA splicing and degradation, μ. We find the parameter estimates are robust to changes in the value of ρ (S7 Fig), and that the estimated production rate of CA-US HIV RNAs, α, decreases approximately linearly with decreases in μ (S8 Fig). Thus, if the rate of CA-US HIV RNA loss is lower in cells activated by vorinostat than we have estimated, the estimated production rate of CA-US HIV RNA would also be lower. Last, we tested the robustness of our results to the assumption in the Rong et al. (15) model about how the latently infected cell population is maintained by employing a different model based on the work of Kim and Perelson [50] in which the latent population is maintained by homeostatic proliferation rather than asymmetric division (see Methods). We found the model fits to the CA-US HIV RNA data and the estimates of α and d_L are largely unaffected (S9 Fig).

Discussion

We have constructed mathematical models to describe the dynamics of CA-US HIV RNA in HIV-infected individuals on ART who received multiple doses of the HDAC inhibitor vorinostat. By fitting these models to a clinical dataset, we have assessed the dynamic response of latently infected cells to vorinostat and estimated the quantitative impact of vorinostat on the latently infected cell population.

Model analyses show that the multistage delayed activation model, can describe both the short-term and the long-term patterns of change in CA-US HIV RNA induced by vorinostat in most individuals. This model assumes that in response to vorinostat treatment, HIV transcription in latently infected cells is induced transiently. Afterwards, the cells rather than returning to their original latent state go through several waiting stages where CA-US HIV RNAs are not produced but host gene expression patterns may change before becoming sustainably induced. The sustained induction of HIV transcription may even occur after vorinostat treatment is stopped. The induction of HIV gene expression depends on the availability of the HIV Tat protein as well as many host factors [34–36,42–44]. In latently infected cells, the number of Tat proteins [34] and the host factors necessary for inducing HIV transcription, such as P-TEFb, are likely to be at low levels [42–44], and at the same time, the presence of inhibitory molecules, such as Murr1, human schlafen 11 and tetherin, prevent transcriptional activation [45–47] before and at the early stage of response to vorinostat treatment. A recent proteomics and transcriptomics study showed that after 24 h of vorinostat treatment of primary CD4+ T cells the expression of a large number of host genes and proteins as well as genes and proteins previously reported to be involved in HIV transcription was modulated, with some effects appearing to be stimulatory and others inhibitory for HIV reactivation [51]. Therefore, it is likely that the immediate impact of vorinostat on histone acetylation and host gene transcription lead only to a transient induction of HIV RNA transcription and sustained HIV transcription may depend on the longer-term impact of vorinostat on host gene transcription [10]. This delay in sustained transcriptional induction may explain the later increase in the level of CA-US HIV RNA after cessation of vorinostat seen in this study, and the observed refractory periods in response to multiple doses of vorinostat in another study [52]. Note that the effect of vorinostat on host genes may also include the generation of read through transcripts containing HIV RNA [53], but a recent report suggests such transcripts are a minor fraction of total gag RNA [54].

Analyzing the model, we found that the number of stages latently infected cells goes through the waiting state and the total waiting period before sustained induced transcription varied among individuals. This suggests that latently infected cells in different individuals may be in different states, possibly due to variations in Tat protein copy number, host gene expression or alternatively different degrees of chromatin silencing or configuration potentially dependent on the sites of HIV integration. This, in turn, would cause responses to vorinostat to be heterogeneous. Interestingly, the maximal fold increase of CA-US HIV RNA was strongly correlated with the basal level of CA-US HIV RNA before vorinostat treatment [10]. Thus, it is plausible that the basal level of CA-US HIV RNA serves an indicator of the status of latency in a patient and the ease of induction of transcription using LRAs, suggesting that future treatment strategies may be able to be tailored to individual patients.

We further assessed the impact of vorinostat on the rate of loss of cells in the sustained activated state. This estimated loss rate, which serves as an upper bound on the death rate of activated cells (as cells could lose their activated state), is extremely low in most individuals, suggesting that vorinostat treatment does not induce killing of transcriptionally activated cells in most participants. This is in agreement with several previous in vitro and ex vivo studies showing vorinostat activates HIV transcription in only a subset of cells and that this level of HIV transcription and protein expression does not lead to cell death [21,52,53,55,56]. Interestingly, a recent in vitro study showed that vorinostat treatment only has significant impact on HIV transcriptional activation, with the impact on translation being minimal, suggesting that HIV proteins may not be produced sufficiently to lead to virion production or to induce viral cytopathic or cytotoxic T cell mediated cell death [20]. Thus, new treatment strategies aiming at both transcriptional and translational activation of HIV may be needed to induce efficient killing of latently infected cells.

In the model, we have assumed that in each participant, latently infected cells are a homogeneous population and respond to vorinostat by going through two activation steps. However, because of limited data sampling there could have been additional transient activation steps that we were unable to detect. In addition, the latent state of individual infected cells in vivo may differ [57] and it is possible that within an individual, some cells go through a different number of activation steps or have different waiting periods before becoming sustainably activated. This may give rise to the minor discrepancies between the data and the model seen in some participants (VOR010, VOR018, VOR019, VOR021, VOR023 in Fig 3). Although a model that accounts for different responses of latently infected cell subpopulations or has additional activation steps might explain the data, such a model would have more unknown parameters than our current model with only a marginal improvement in model fit. Nonetheless, the possibility that there exist different cell populations in individual patients in terms of their response to vorinostat or additional activation steps cannot be excluded. Further experiments examining the dynamics of host factors/enzymes that are responsible for transcriptional activation and inhibition under LRA treatment could validate our model, improve our understanding of the impact of vorinostat, and ultimately aid the design of treatment strategies to eradicate the latent reservoir.

To conclude, our results suggest that vorinostat induces both immediate transient induction and delayed sustained induction of HIV transcription. Similar dynamic patterns of CA-US HIV RNA were also observed in clinical trials of the LRAs panobinostat and romidepsin [58,59], suggesting LRAs may induce both transient and delayed transcription activation in latently infected cells in general. Therefore, designing clinical trials with frequent longitudinal sampling during both treatment and the follow-up period would help quantify the impact of LRAs. To our knowledge, our work represents the first mathematical model to assess the impact of a LRA on the dynamics of CA-US HIV RNA in vivo. Our model can be easily adapted to study other LRAs as well as combinations of these agents once data are available. In addition, our model or variants of it could be used to assess the efficacy of different candidate treatments, such as those using anti-HIV monoclonal antibodies combined with LRAs, and ultimately suggest optimal drug combinations to eliminate latently infected cells in HIV-infected individuals.

Methods

Clinical Data

As described by Elliot et al. [10], 20 chronically HIV-infected adults receiving at least three antiretroviral agents, having plasma HIV RNA < 50 copies per mL for at least three years (excluding single viral ‘blips’), a CD4+ T-cell count > 500 cells/μL and documented subtype B HIV-1 infection were recruited into a vorinostat trial (ClinicalTrials.gov, NCT01365065). Participants received vorinostat 400 mg orally once daily for 14 days. Levels of CA-US HIV RNA were measured in peripheral blood mononuclear cells at 0, 2, 8 and 24 hours, and on days 7, 14, 21, 28 and 84 (as well as on day 3 for participants VOR019-023). For each blood sample, four replicate q-PCR runs were performed to measure the levels of CA-US HIV RNA. The total number of data points in each patient ranges from 35 to 44.

The Direct Activation Model

The ordinary differential equations (ODEs) describing the model are:

In this model, target cells, T, are produced at a constant rate, s, and die at per capita rate, d_T. In the absence of cART, they are infected at per capita rate, βV_I, where β is a rate constant and V_I is the concentration of infectious virus. The effect of reverse transcriptase inhibitors (RTI) is to multiply the rate of infection by the factor 1-ε_RT, where ε_RT is the effectiveness of the RTI with 0 ≤ ε_RT ≤ 1. Under protease inhibitor treatment, a fraction, ε_PI, of produced viruses are non-infectious (V_NI). Upon infection with viruses, a fraction, f, of infected cells becomes latently infected, and the remaining fraction, 1-f, becomes productively infected. Productively infected cells, I, die at per capita rate, δI^1.44, where the power 1.44 models the effects of both viral cytopathicity and cell-mediated immune responses [60,61].

We assume that latently infected cells, L, when activated by antigen at rate a, undergo an asymmetric division in which a daughter cell remains in latency with probability, p_L, and becomes productively infected and produces virus with probability, 1-p_L [15]. Latently infected cells die at per capita rate, d_L. We assume that vorinostat causes latently infected cells to move into a transcriptionally activated state, L_A, at rate ν. The transcriptionally activated cells are lost at rate, d_LA. We also assume there is a pharmacological delay, t₀, in the effect of vorinostat upon treatment initiation such that the rate of activation remains 0 when t<t₀, and it becomes ν when t₀≤t≤14 days. After treatment is terminated, we assume the effectiveness declines exponentially as νe^-w(t-14) for t>14 days, where w is the rate at which vorinostat is cleared from the system. We also model the population average of the amount of CA-US HIV RNA (R) within transcriptionally activated cells. We assume that CA-US HIV RNA is produced at a constant rate, α, once vorinostat becomes effective, i.e. when t≥ t₀. CA-US HIV RNA is encapsidated and exported as virions at per capita rate, ρ, and lost by degradation and splicing at per capita rate, μ. Here, we assume that the transcriptionally activated cells do not produce mature viral particles, i.e. ρ = 0, because recent work shows that the production of virus from cells treated with vorinostat is minimal [20,53]. As other LRAs may induce viral production, we leave ρ in the model. Viruses are produced from productively infected cells at rate p_v and from transcriptionally activated cells at rate ρR (note that ρ is set to 0 for vorinostat). Viruses are cleared at per capita rate c. The amount of CA-US HIV RNAs, US, is calculated as the sum of the basal level of CA-US HIV RNA, US₀, and the average number of US HIV RNAs per activated cell, R, multiplied by the number of activated cells, L_A.

The Delayed Activation Model

In the delayed activation model, the equations describing the transient activation state, L_T, the waiting state, L_W, the transcriptionally activated state, L_A, and the total number of CA-US HIV RNAs, US, are

where k_T and k_w are transition rate constants. The ODEs describing other variables are kept the same as in the direct activation model.

The Multistage Delayed Activation Model

The ODEs describing the multistage waiting states are:

where the waiting state is divided into n identical sub-states, i.e. L_W,1, L_W,2, … L_W,n. The transition rate from one sub-state to the next is set to nk_w, so that the average residence time in the waiting state, 1/k_w, is the same as in the delayed activation model. The ODEs describing other variables are the same as in the delayed activation model.

A Model Assuming Constant Proliferation of Latently Infected Cells

To test the assumption about how the population of latently infected cells is maintained, we modified the multi-stage delayed activation such that the latently infected cells proliferate at a constant rate as in Kim and Perelson [50] instead of asymmetric division in the main text. The ordinary differential equations (ODEs) describing the productively infected population and latently infected population are shown below and other terms in the model are kept the same as in the main text.

In this model, the latently infected cells becomes activated at per capita rate a and proliferate at per capita rate r. We set a = 0.0088 day^-1 according to [62] and r = 0.0183 day^-1, such that the half-life of the latent reservoir is 44 months as estimated before [63].

Parameter Values

The fixed parameter values are based on prior work to yield a baseline state with a plasma viral load of 5 copies/mL and a latently infected cell population of ~2 cells/million cells (Table 1). The production and death rates of uninfected CD4⁺ T cells were chosen to yield a baseline CD4 count similar to the average level in patients in Elliot et al. [10]. Changes in the values of the fixed parameters that govern the dynamics of target cells and productively infected cells do not impact the estimation of the fitted parameters (S10 Fig). This is because cART is so effective that the contribution of new infections to the latent reservoir is negligible over the time period we study here. We also tested the robustness of the parameter estimates against changes in the fixed parameters that govern intracellular HIV transactivation in latently infected cells and to the structure of the model with regard to how latently infected cells proliferate (S7–S9 Figs).

The CA-US HIV RNA production rate in a productively infected cell, α_I, is calculated based on the derivations in Ref. [28]. The basal transcription rate in eukaryotic cells is estimated to be approximately 40 nucleotides/s [64]. The HIV genome has around 9500 nucleotides. Thus, the basal production rate of CA-US HIV RNA production can be estimated as 40/9500≈0.0047 s^-1≈406 day^-1. In the presence of Tat, the transcription is upregulated to 100 fold [65]. Then, the production rate of CA-US HIV RNA in a productively infected cell, α₀, can be calculated as approximately 40,000 transcripts day^-1. This maybe a minimal estimate if productively infected cells live approximately one day while producing virus, as Chen et al. [32] have estimated that a productively infected cell produces between 40,000 and 55, 000 virions over its lifespan.

Model Fitting and Model Selection

To fit the direct activation model to the data, we varied 5 unknown model parameters: the production rate of CA-US HIV RNA in a transcriptionally activated cell (α), the death rate of transcriptionally activated cells (d_LA), the rate of transcriptional activation (ν), the initial delay of vorinostat effectiveness (t₀) and the baseline US HIV RNA level (US₀). To fit the delayed activation model and the multistage delayed activation model, we allow 2 additional parameters to be estimated (together with the 5 parameters in the direct-activation model): the transition rate, k_T, from the transiently activated state to the waiting state, and the transition rate, k_w, from one sub-state of the waiting state to the next.

To estimate the parameter values in each model, we first calculate the residual sum of squares (RSS) between model predicted log CA-US RNA level and log transformed data, and then minimize the RSS using the Nelder-Mead algorithm [66]. 1,000 individual fits were performed for each model starting from parameter values randomly sampled within biologically plausible ranges. The parameter values with the smallest RSS among the 1,000 fits are taken as the best-fit parameter values for each model.

We perform model selection using the corrected Akaike information criterion (AICc) to account for the low number of data points (ranges from 35 to 44) for each patient [67]. The AICc score is calculated as

where n is the number of data points and K is the number of fitted parameters. When comparing models, the model with the lowest score is the best model, although small difference in AICc scores, e.g. ≤ 2, is not significant [67].

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10

Attachment 11

Attachment 12

Attachment 13

Attachment 14