Motivated by viral persistence in HIV+ patients on long-term anti-retroviral treatment (ART), we present a stochastic model of HIV viral dynamics in the blood stream. We consider the hypothesis that the residual viremia in patients on ART can be explained principally by the activation of cells latently infected by HIV before the initiation of ART and that viral blips (clinically-observed short periods of detectable viral load) represent large deviations from the mean. We model the system as a continuous-time, multi-type branching process. Deriving equations for the probability generating function we use a novel numerical approach to extract the probability distributions for latent reservoir sizes and viral loads. We find that latent reservoir extinction-time distributions underscore the importance of considering reservoir dynamics beyond simply the half-life. We calculate blip amplitudes and frequencies by computing complete viral load probability distributions, and study the duration of viral blips via direct numerical simulation. We find that our model qualitatively reproduces short small-amplitude blips detected in clinical studies of treated HIV infection. Stochastic models of this type provide insight into treatment-outcome variability that cannot be found from deterministic models.
There are many locations in the body from which viruses could re-emerge during drug treatment; for a review, see [6]. Here, we will focus on the important possibility that viruses may emerge from a reservoir of latently infected cells. Usually when HIV infects target cells (such as CD4+ T cells and macrophages) the result is rapid virus production and cell death. However, a fraction of infected cells are known to enter a state of latent infection [12] where virus has integrated into the host cell DNA, but there is little, if any, viral gene expression. While cells are in this state, they are unaffected by ART and viral cytopathicity, and are effectively invisible to the host immune response [13]. However, upon re-activation, latent cells begin the normal processes of viral replication and production, and become immune targets [14]. A large fraction of the latent reservoir consists of resting memory CD4+ cells [15] and therefore, reactivation could occur as part of the normal immune response to a secondary pathogen [16]. However, we do not completely understand the reasons for activation of latently infected cells and it is likely that there is a pathogen-independent component as well. Indeed, the mechanisms for latency are generally poorly understood; there are differing opinions, but no consensus to date [12], [17]–[19]. The population of latently infected cells is established as early as 10 days after symptoms of seroconversion, within a few weeks of initial infection [14]. Estimates of reservoir size differ but consistently show that latently infected cells constitute only a small fraction of the total number of T-cells [12], [20]. Unfortunately, and in spite of its small relative size, the decline of this population is slow and it is estimated that it can persist for up to 70 years [21]. This long lifetime probably arises from the intrinsic stability of resting memory CD4+ cells which is an important part of immune memory [22]. Recent evidence also indicates that latently infected cells can undergo cell division [15], potentially increasing the lifetime of the reservoir. These factors, in combination with long lifetime of the reservoir, indicate that latently infected cells are an important factor that must be addressed in the search for therapies to eradicate HIV infection [23]. While on successful anti-retroviral treatment (ART) for HIV, an infected individual' s viral load remains non-zero [3], though it is very low and usually undetectable using standard assays that have a detection limit of 50 copies/mL in plasma. Occasionally, however, regular blood tests show viral blips: periods of detectable viral load, preceded and followed by undetectable loads. At one time there was a concern that blips might signal imminent drug failure [24], including the emergence of new, drug-resistant variants of virus [25]. However, there is a substantial body of evidence dating from the early 2000s, indicating that most blips are not associated with virological failure [26]–[28]. With that said, a recent large-scale study of 3530 Canadian patients refined these results by showing a two-fold increase in the risk of drug failure following viral blips that exceeded 500 copies/mL, but importantly, smaller blips were not associated with drug failure [29]. In this study, blips were detected at a frequency of about 0. 1/patient/year. This rate is compatible with data taken from the UK during 2006–2007 [30] and is lower than the rate estimated from earlier data [31]. The reduction in blip frequency over the last decade is likely a result of improved drug efficacy. The underlying cause of viral blips remains controversial. There is some evidence that immune activation, through secondary infection or vaccination, may be correlated with viral blips, [32], [33]. However, there have been observations of blips not associated with clinical or demographic variables. In an intensive 90-day study of 10 patients, Nettles et al. found that blips were fairly common, smaller in amplitude (mean 79 copies/mL) and short in duration (median less than 3 days), and that blip frequency was unrelated to illness, vaccination, or drug concentrations [34]. Finally, we must acknowledge that accurate detection of small-amplitude blips is bedevilled by assay variability and sensitivity [30], [34]. There has been extensive modeling work done to characterize viral load and pathogen-immune system interaction in HIV-infected individuals. However, standard viral dynamics models do not capture residual viremia in treated patients and various modeling approaches have been applied. Residual viremia can be captured by adding a latent cell reservoir to the standard model. Perelson et al. (1997) proposed the first model that included latent cell activation, in order to better understand decay characteristics of HIV-1-infected compartments during combination therapy [35]. This model was expanded to include a varying decay rate in the latent reservoir, and bystander proliferation in the latent reservoir, with the finding that a constant long-term activation rate for the latent reservoir, maintained through cell division, could explain residual viremia in treated patients [36]. We will use these elements in the development of our model of the latent reservoir. Careful modelling of viral blips has been fruitful in analyzing different mechanisms of blip generation. The focus of previous models has been on blips associated with immune system activation due secondary infection or vaccination. One successful approach has been to consider short periods of sustained viral replication driven by stochastic activation of CD4+ and CD8+ T cells [37]–[39]. A further series of models including T cell expansion due to vaccination and secondary infection showed episodes of detectable viremia of long duration (2–3 months), with amplitudes in the range of several hundred copies/mL [40], [41]. Viral blips can also result from production of virus following immune activation and clonal expansion of latently infected cells [33]; latent cell activation caused by sporadic immune activation has also been modeled as a source of viral blip generation. Indeed, antigen-induced latent cell activation has been modeled and shown to be a plausible source of viral blips [42]. Most recently, Rong and Perelson (2009) proposed a model with antigen-induced asymmetric activation and division of latently infected cells [43]. Blips produced by this model are of short duration, directly related to the length of stimulation, and of variable amplitude, consistent with observations. These models produce blips of larger (100 copies/mL) amplitude, with variable durations, and frequency depending directly on user-controlled periods of immune system activation in the model. The base mechanism in these models of the production of blips is immune system activation. As noted above, there have been observations of small-amplitude blips not associated with clinical or demographic variables [34]. Such blips can be imagined as random biological or statistical variation around a mean undetectable viral load. In order to capture this kind of stochastic effect, continuous (differential equation based) models are inadequate. Here, we propose a continuous-time branching process model of within-host viral dynamics for a patient undergoing successful treatment. We use this formulation to derive probability distribution functions for viral load as a function of time and examine the contribution of varying latent cell activation and proliferation to viral load. Using this methodology we first consider extinction times for the latent reservoir and examine the role of limited ongoing viral replication in replenishing the reservoir. We then examine the hypothesis that stochastic activation of latently infected cells can maintain low-level plasma viremia and generate small intermittent viral blips. Finally, via Gillespie simulation of the branching process, we calculate viral blip durations.
We consider a simple model of latent cell reactivation, presented schematically in Figure 1. Our model has three compartments: the number of latently infected cells, which can replicate at rate die at rate and activate at rate to become productively infected cells; the number of productively infected cells which die at rate; and the number of virions, produced by productively infected cells at rate, which can die at a rate. We allow for infection of new cells at rate, of which a fraction become latently infected. The efficacy of ART is given by. We will assume that this efficacy is very high, and that therefore the number of uninfected T-cells remains approximately constant and equal to. Clinical findings on viral blips show differing amplitudes [34], [44], [45]; though the small-amplitude blips were shown to be unassociated with clinical variables [34], it is possible that the larger-amplitude blips may be due to an immune response, increasing the activation rate for a period of time. We will initially restrict ourselves to constant activation rate, but we consider variable in a later section. The mean behaviour of the system shown in Figure 1 is given by the linear system of ordinary differential equations (1) where, and represent the mean numbers of latently infected cells, productively infected cells, and virions respectively. Our goal is to obtain probability distributions for the viral load and the size of the latent reservoir at time. We assume that the events in the model can be described by a multi-type continuous time branching process with the rates given in Figure 1. Importantly, the model does not scale up and so any computations we perform must be over the total number of and in the patient. The transition probabilities for each process in the model are stationary in time. We therefore know that our desired probability distributions depend only on the time since the initial state, and considerthe conditional probability that there are latently infected cells at time, productively infected cells at time, and virions at time given that there were initially, and of each species respectively. Then, given a joint initial distribution on these species, we can compute the joint probability distribution on each of theseNote that, as the latent reservoir must be of finite size, and as. By considering each process in Figure 1 in turn, we can derive the backwards Chapman-Kolmogorov differential equation for [46]: (2) with initial condition (is the Kronecker delta function). Multiplying through by and re-indexing yields an infinite-dimensional system of nonlinear ordinary differential equations for the conditional probability generating functionWe reduce the infinite dimensional system to a system of three equations by exploiting the branching property [46], with initial conditions and. To our knowledge we cannot solve this nonlinear system analytically. Therefore to calculate, we solve numerically using a standard differential equation integrator. Once and are calculated we can compute the full probability generating function, accounting for the initial distributionsOur goal is to recover the probability distributions of latently infected cells, productively infected cells, and virions at times. These can be recovered from the probability generating function by taking derivatives. For example, the probability that there are virions at time is given byBecause the distributions do not scale, we must perform computations over the total number of virions. Given a mean viral load of 25 copies/ml (henceforth abbreviated as 5 c/mL) within 5L of total blood volume, we must compute 125000 derivatives to get! Direct numerical differentiation would be difficult, so we exploit the Cauchy-Euler formula: where is a closed curve in complex space which contains, and is analytic on a simply connected domain containing. The probability generating function is a polynomial in, and and therefore satisfies the analyticity requirement. We want to evaluate integrals at so our contour must contain the origin and it is simplest to use the unit circle, , Thenwhere we have used the fact that, where indicates complex conjugate. By this method we can calculate our probabilities via straightforward and reliable numerical integration. The same approach can be used to compute joint probability distributions. We can also use this formulation to directly calculate cumulative probabilities. As, we can writeby interchanging the order of integration and summation. This final formula will be useful in calculating blip probabilities at a time, . To our knowledge this is a novel method for computing probability distributions from single- or multi-type continuous time branching processes. We thoroughly tested our method and its implementation; see Figure S1 for comparisons with Gillespie simulations. We also wish to calculate the distribution of times to extinction for the latent reservoir. We choose parameters so that the probability of extinction of the latent reservoir is 1 as. However, as clearing the latent reservoir is considered a major hurdle in clearing HIV, the distribution of times to this inevitable extinction is of interest. We obtain the cumulative probability of latent reservoir extinction directly from the probability generating function. Since Note that the marginal probability. We then find the probability distribution of extinction times by differentiating, If we assume that no newly infected cells become latently infected () or that treatment is completely effective (), we can obtain the extinction probability analytically. In this case, the latent cell dynamics decouples from the rest of the model and can be represented as a pure birth-and-death process with master equation (3) where is the probability that at time there are latently infected cells. This probability has the conditional probability generating functionwhere is the initial reservoir size. The cumulative distribution is and thus we obtain (4) Then given the initial distribution on the latent reservoir we have can compute the extinction probability. We use the analytic expression to compute latent reservoir extinction time distributions for or. Otherwise, we work numerically. The parameters used for simulation results presented below are given in Table 1. In our simulations the parameters, and are chosen based on estimates from [36] and based on estimates from [47]. The decay rate of the latent reservoir is chosen so that its half-life is, as measured in patients exhibiting viral blips [21]. For, the death rate of productively infected cells, and the virion clearance rate we set to estimates from [48] () and [49] (), respectively. The in vivo viral production rate is not well established and therefore we will consider a range for this parameter. The fraction of new viral infections that result in latency is also hard to estimate, but given the small size of the latent reservoir, it is likely to be rather small. For simplicity, we choose a baseline value of. We choose the initial mean latent reservoir size consistent with the experimental estimates [21], [50]. The activation rate and replication rate of latently infected cells remain unknown. We calculate values from the mean behaviour equations (1), taking. Since the dynamics of the productively infected cells and virus are more rapid than those of the latent reservoir, we can make a quasi-steady approximation to find. Then for an initial latent reservoir size of, . To calculate, we choose the mean decay rate so that the half-life of the latent reservoir is. Thus we can set. The resulting () values for each production rate and mean viral loads of or are given in Tables 2 and 3 respectively. As noted above, we must perform calculations over the entire blood volume, which we take to be. When presenting results below we report viral loads in copies per mL, as this is the standard measurement, but they are always obtained by re-scaling the axes for results over the entire blood volume. In order to correctly simulate viral blips and latent reservoir extinction in patients with established treated infection, we should carefully choose the initial joint distribution so that it is close to the (moving) equilibrium of the ongoing dynamics. Otherwise, transient effects will pollute our results. In the mean, the dynamics of the latent reservoir are very slow compared to those of the productively infected cells or virions. We therefore focus on getting the initial latent reservoir distribution correct since errors in the other two compartments will resolve themselves quickly. Indeed, for a constant latent reservoir size, and our parameters, the distributions on and converge to stationary distributions in less than a month (results not shown). In order to calculate a reasonable initial latent reservoir distribution we isolate its dynamics and consider the marginal probability distribution only, as in equation (3). We choose the marginal latent reservoir probability distribution at time such that the variance is maximized. We reason that transient dynamics on the latent reservoir are dominated by the spreading of the distribution about the decaying mean, and that at maximum variance the probabilities are sufficiently spread for our purposes. For birth-and-death processes maximum variance occurs at the half life. Therefore, in order to create the initial distribution on the latent reservoir, we solve (3) out to, starting with latently infected cells, where is the desired mean latent reservoir size. The resulting distributions for different parameter sets are shown in Figure 2. Notice that results based on a virus production rate have larger standard deviation. This is because lower production rates are associated with higher activation rates (cf. Tables 2 and 3). The higher activation rate speeds the dynamics of the latent reservoir, increasing the spread of its probability distribution function. Finally, we combine the computed initial latent cell distribution with single initial numbers of productively infected cells and virus, to obtain the whole initial joint probability distribution:
The reservoir of latently infected cells is considered a major obstacle to clearing HIV infection [19]. Within our model, when the reservoir goes extinct, viral load quickly goes to zero, since ongoing viral replication is too small to sustain the virus population. We are therefore interested in examining the reservoir lifetime after the onset of ART. To do this, we extend our approach to find the probability of reservoir extinction over time. We examine the reservoir lifetime using baseline parameters (), and then allowing for the possibility of latent reservoir replenishment (). Furthermore, since anti-retroviral treatments have improved substantially over the last 15 years, we also examine how the reservoir lifetime behaves as drug efficacy improves (). We now focus on the time evolution of viral load and the likelihood of small-amplitude viral blips. We interpret viral load above the threshold of detection of 50 c/mL as a viral blip. Note that unless otherwise specified, the following calculations and computations assume the fraction of newly infected cells that become latently infected is 0. In Figures 6 and 7, we plot full viral load distributions over time, assuming initial mean viral loads of 25 c/mL and 35 c/ml. As time advances, the mean viral load decreases as expected in all cases but the viral load distributions widen more significantly when is smaller (e. g. distributions in Figure 6A are widest, those in Figure 6C are narrowest). This is because the lower values of are associated with higher values of and, and the resulting dynamics on the latent reservoir cause the latent reservoir size probability distribution (not shown) to be wider. As a consequence the associated viral load distributions are wider, and this also causes higher blip amplitudes. This effect is more clearly understood by examining the insets in Figures 6 and 7, which represents a magnified view of the given probability distribution curves above the blip threshold (50 c/mL), using a log scale to more clearly distinguish the curves. We observe that viral blips occur with very small probability regardless of the production rate. The blip amplitudes vary between parameter sets but remain approximately within the range of blips unassociated with clinical variables shown in [34], i. e. 50–100 c/ml. Over three years, the range of reasonably likely detectable viral loads decays slowly, but small blips remain possible throughout that time (Figure 6). In this section we consider the following question: given a blip, defined as a detectable viral load measurement, how long should we expect the viral load to remain above the threshold of detection? This question is of clinical interest, since a repeat measurement following a measurable viral load should be performed after enough time that a second positive result might have clinical significance, such as suggesting drug failure. Different from the previous sections, all results in this section are computed via 10000 direct simulations of the branching process using the Gillespie algorithm, beginning with an initial “blip” condition. The initial conditions are chosen as follows: we set the latent reservoir size and viral load. Since dynamics on the viral load are so much faster than on the productively infected cells (cf. Table 1), we then use a quasi-steady approximation to set the initial number of productively infected cells.
Clinical results on latent reservoir decay (e. g. [21]) make predictions on latent reservoir lifetimes that are based on purely exponential decay. Our model results showed that, for the same mean decay rate, the time distribution - and the mean time to extinction - is sensitive to dynamics on the latent reservoir. Further, assuming some reservoir replenishment due to latency in newly infected cells, our model predicted only limited lifetime reduction associated with improving drug efficacy. Eradication of the latent reservoir is considered a major hurdle in eradicating HIV infection [19], [23], and these results demonstrate the importance of understanding the underlying dynamics on the latent reservoir: the reservoir half-life is only a small part of the equation. From a clinical point of view our results on the latent reservoir lifetime are quite depressing - in our model we essentially study perfectly drug-adherent patients, and even with perfect drugs, decades of drug treatment are needed to clear the latent cell reservoir. Our model predicts that drug treatments that increase the activation rate of latently infected cells should reduce the lifetime of the reservoir. This approach has been tried several times but so far without real success (reviews in [1], [53]). A candidate drug would need to work on the whole heterogeneous population of cells seeded during initial infection, and in particular on the longest-lived subpopulation. We would predict that early treatment with such a drug, along with aggressive ART, would be most likely to reduce the size of the latent reservoir. This is in line with current research in treating HIV infection earlier, to enhance survival on an individual level [54] and limit transmission on the population level [55]. Of concern with earlier treatment is the possibility of emergent drug resistance (DR). In future work we plan to expand our stochastic model to examine the likelihood of different mechanisms of acquired DR in patients on treatment, such as mutation during ongoing viral replication and activation of a cell latently infected with a DR strain [56]. We also examined the evolution of viral load over time, finding that as time progresses, the viral load distributions become more asymmetric, with a long tail towards higher viral loads. This can be explained by viewing our model as an extended subcritical birth-and-death process. Such processes produce asymmetric distributions (see Figures 6 and 7). The asymmetry is more pronounced for smaller production rates (vs), associated with larger activation rates, and for larger initial viral loads (vs). When examining blip probability we found that our model predicts that these probabilities decay exponentially over time. This decay is more dramatic for larger production rates, associated with smaller activation (‘birth’) rates, and smaller initial mean viral loads. In only one case is the decay so slow (with mean initial viral load) that our model predictions are broadly consistent with previous observations that blip probabilities don' t decay over time [45]. We also observed that blip probabilities show great sensitivity to model parameters, varying by orders of magnitude. Given sufficient high-quality data on blips (which does not currently exist), these would be the ideal results to compare with data for the purpose of parameter fitting, in order to gain some insight into latent reservoir dynamics - into the activation rate, for example. With improving drug efficacy () our model predicts a significant decrease in baseline viral load, down to, in accordance with recent viral load observations [10], [11]. From this low baseline, we found blip probabilities smaller than we can calculate. However, we found that by increasing the activation rate (roughly simulating immune system activation due to secondary infection) viral loads exceeding the threshold of detection were attainable. Therefore, our model supports the hypothesis that, for patients adhering to modern ART, viral blips signal an underlying secondary condition. We examined the duration of viral blips through direct (Gillespie) simulation of the model. We sought to answer the question “Given a patient measurement of X c/mL, how long can we expect the viral load to remain detectable? ” We found, unsurprisingly, that blips of larger initial amplitude have longer mean duration and larger standard deviation in duration. Perhaps more interestingly, we found very strong dependence on the production rate. Given an initial blip of amplitude 90 c/mL, doubling the production rate from to, and changing other parameters accordingly, more than triples the predicted mean blip duration (for parameters as in Figure 9). We also considered blip duration as a function of the latent reservoir size, anticipating longer durations for larger reservoir sizes, since associated with these is a higher quasi-steady mean viral load. Our expectations were confirmed by simulation results (see Figure 11). The sensitivity to production rate and associated parameters was also recovered. Repeat-blip measurements in patients are, predictably, rather rare, since blips are already quite unusual events. Across all the parameter sets we examined, we found that detectable viremia should be expected to vanish within 8–10 days at most. This result is in general agreement with previous reports [34], [45] and indicates that repeat low-level detectable viremia within 8–10 days could be due to a statistical fluctuation rather than drug resistance or other pathology. Over the last 15 years, enormous numbers of differential-equation models have been generated to study different aspects of various viral infections. We believe that stochastic models of the kind described here have an important role to play in certain situations where viral or cell populations are small enough that random effects still play a role. The obvious settings are during the first few days of any new infection (see also [57] and, very recently, [58]), during drug treatment of a chronic infection, and during the extinction phase of an acute infection. One issue with stochastic modeling of rare events (such as viral blips in our model) is that simulation-driven studies can require enormous numbers of simulations to reliably sample the rare events. The method we describe here is an alternative to simulation (or methods to capture rare-events) and provides a direct and relatively straightforward way to calculate probability distribution functions. We hope to adapt this method to other situations in viral dynamics in future work.
While on successful drug treatment, routine testing does not usually detect virus in the blood of an HIV patient. However, more sensitive techniques can detect extremely low levels of virus. Occasionally, routine blood tests show “viral blips”: short periods of elevated, detectable viral load. We explore the hypothesis that residual low-level viral load can be largely explained by re-activation of cells that were infected before the initiation of treatment, and that viral blips can be viewed as occasional statistical events. To do this, we propose a mathematical model of latently-infected cells, activated cells, and virus. The model captures random fluctuations of the system as well as the mean behaviour. We estimate the time it takes for all the latently-infected cells to be eradicated. Eradication of these cells is considered a major hurdle in eliminating infection. We predict a wide range of eradication times, highlighting the importance of studying latently-infected cells. We also estimate the frequency and duration of viral blips, and find qualitative agreement with clinical studies. By refining our models, we hope to find guidelines that can be used in practise to distinguish between clinically insignificant statistical blips, and instances of drug failure.
medicine
population dynamics
mathematics
population modeling
population biology
stochastic processes
infectious diseases
theoretical biology
hiv
biology
infectious disease modeling
probability theory
hiv diagnosis and management
ecosystem modeling
disease dynamics
viral diseases
computational biology
A Stochastic Model of Latently Infected Cell Reactivation and Viral Blip Generation in Treated HIV Patients