Pre-existing resistance in the latent reservoir can compromise VRC01 therapy during chronic HIV-1 infection

Passive immunization with broadly neutralizing antibodies (bNAbs) of HIV-1 appears a promising strategy for eliciting long-term HIV-1 remission. When administered concomitantly with the cessation of antiretroviral therapy (ART) to patients with established viremic control, bNAb therapy is expected to prolong remission. Surprisingly, in clinical trials on chronic HIV-1 patients, the bNAb VRC01 failed to prolong remission substantially. Identifying the cause of this failure is important for improving VRC01-based therapies and unraveling potential vulnerabilities of other bNAbs. In the trials, viremia resurged rapidly in most patients despite suppressive VRC01 concentrations in circulation, suggesting that VRC01 resistance was the likely cause of failure. ART swiftly halts viral replication, precluding the development of resistance during ART. If resistance were to emerge post ART, virological breakthrough would have taken longer than without VRC01 therapy. We hypothesized therefore that VRC01-resistant strains must have been formed before ART initiation, survived ART in latently infected cells, and been activated during VRC01 therapy, causing treatment failure. Current assays preclude testing this hypothesis experimentally. We developed a mathematical model based on the hypothesis and challenged it with available clinical data. The model integrated within-host HIV-1 evolution, stochastic latency reactivation, and viral dynamics with multiple-dose VRC01 pharmacokinetics. The model predicted that single but not higher VRC01-resistant mutants would pre-exist in the latent reservoir. We constructed a virtual patient population that parsimoniously recapitulated inter-patient variations. Model predictions with this population quantitatively captured data of VRC01 failure from clinical trials, presenting strong evidence supporting the hypothesis. We attributed VRC01 failure to single-mutant VRC01-resistant proviruses in the latent reservoir triggering viral recrudescence, particularly when VRC01 was at trough levels. Pre-existing resistant proviruses in the latent reservoir may similarly compromise other bNAbs. Our study provides a framework for designing bNAb-based therapeutic protocols that would avert such failure and maximize HIV-1 remission.

Introduction Antiretroviral therapy (ART) for HIV-1 infection rapidly suppresses viremia to undetectable levels and curtails disease progression but is unable to eradicate the virus [1]. Discontinuation of ART, even long after viremic control is established, typically leads to rapid viral rebound, often within days to weeks of discontinuation, and to progressive disease [2]. The rebound is caused by a reservoir of latently infected cells [3] that is formed soon after infection [4]. The reservoir is sustained by cell proliferation [5,6], which can continue even when ART has stopped new infections, allowing the reservoir to exist long-term [7]. The reservoir can get reactivated stochastically and reignite infection once ART is stopped [8,9]. ART must therefore be administered lifelong. In a remarkable breakthrough, the VISCONTI study showed that when ART is administered early in infection, some individuals can maintain viremic control for many years after the cessation of ART [10]. This study has raised hopes of a functional cure, or long-term remission, of HIV-1, where the viremic control once established by ART can be maintained without lifelong treatment [11]. The success of ART in inducing post-treatment control, however, is small: only~5-15% of the patients treated achieve lasting post-treatment control [12]. Enormous efforts are underway to improve this success rate [11].
One strategy that holds promise is to administer broadly neutralizing antibodies (bNAbs) of HIV-1 for a short period post-ART, i.e., during an analytical treatment interruption (ATI) [13,14]. bNAbs target diverse viral genomic variants [15,16], and are expected to suppress viremia arising from the reactivation of latently infected cells [17,18]. Simultaneously, they may engage the host immune system [19], potentiating it to maintain the viremic control longterm [20,21]. Two recent clinical trials tested this strategy using the bNAb VRC01 [17]. VRC01 targets the CD4 binding site on the HIV-1 envelope with high breadth and potency [15,22]. When it was administered to chronically infected patients who had achieved undetectable viremia with ART, the duration of viremic control was observed to increase only marginally, by a median of~2-4 weeks, beyond historical controls [17]. In the historical control group, which was not treated post-ART, viral rebound occurred in~2.6 weeks on average after ART interruption [23]. Why VRC01 was ineffective in maintaining remission longer is unclear. Unravelling the causes of this ineffectiveness is important to optimizing VRC01 usage, which is also in large clinical trials for preventing the transmission of infection [24,25], and to expose potential vulnerabilities of other bNAbs.
VRC01 exhibits potent antiviral activity in vivo [26], including at its trough levels in patients [27]. Yet, in the trials above [17], viral rebound was observed in most patients when VRC01 levels in circulation were significantly higher than its in vitro suppressive concentration (or IC 50 ), implicating the role of resistance. Mutations that confer resistance to bNAbs are well documented [28][29][30][31]. Indeed, VRC01-resistant strains were detected in the breakthrough viral populations in the trials above [17]. The rapid virological breakthrough during treatment, especially given the absence of circulating virions at the time of ART cessation and the potent activity of VRC01 against the wild-type, suggests that VRC01 resistance might have existed before VRC01 therapy. Here, we therefore hypothesized that pre-existing VRC01-resistant proviruses formed before ART and contained in the latent reservoir could underlie the failure of VRC01 therapy. Minority viral variants as well as latently infected cells are difficult to detect using current assays [32,33]. Testing our hypothesis experimentally would require detecting minority variants 'within' latently infected cells, highlighting the challenge involved. Current assays can rarely detect variants below a frequency of~1% [34], implying that, given the prevalent estimates of the latent reservoir size of 10 5 -10 8 cells [4,35], variants present in as many as 10 3 latently infected cells may go undetected and be responsible for therapy failure. Indeed, in many individuals who failed rapidly in the trials above [17], viral genome sequencing could not detect VRC01 resistance pre-treatment [17,36]. As an alternative approach, therefore, we resorted to mathematical modeling.
A number of mathematical models have been developed in recent years to describe latent cell infection and dynamics and its role in the outcomes of treatments [37][38][39][40][41][42][43][44][45][46][47][48]. Models have also been constructed, independently, for describing viral evolution and drug resistance, especially in the context of ART [49][50][51][52][53]. Describing the failure of VRC01 therapy required integrating these two independent formalisms, of viral evolution and latent reservoir reactivation, a task not accomplished so far because of the complexity involved. HIV-1 evolution involves mutation, recombination, fitness selection, and random genetic drift, which together define the timing and speed of the development of drug resistance during ART [51][52][53][54][55]. Latent cell reactivation is an intrinsically stochastic process [8,45,46], following which the virions released must establish lasting infection, which is not guaranteed [37,39,42], especially in the presence of bNAbs. Here, we developed a framework that integrates these processes by recognizing that the dynamics of viral evolution and latent cell reactivation could be decoupled in the context of post-ART bNAb therapy. Viral evolution primarily occurs pre-ART, where viral and productively infected cell populations are large and the contribution from latently infected cells to the dynamics can be ignored. Latent cell reactivation leading to treatment failure occurs during bNAb therapy, when active viral replication is small and so viral evolution can be ignored. These processes are linked by VRC01 resistant strains, which are predominantly formed before therapy and are harbored in latently infected cells and could get reactivated during therapy. Developing a model with this strategy, we were able, for the first time, to capture data from human clinical trials involving bNAb-based interventions quantitatively, offering an explanation of the inadequate effectiveness of VRC01 in maintaining remission, and providing a framework for rational treatment optimization.

Mathematical model
We considered the scenario where chronically infected individuals with viremic control established with long-term ART are administered VRC01 during an ATI, as in recent clinical trials [17]. We developed a model to predict the ensuing remission times based on the hypothesis that viral strains resistant to VRC01 harbored in latently infected cells were responsible for virological breakthrough (Fig 1). We provide an overview of the model here; details are in Methods.
We first considered a single infected individual ( Fig 1A). We estimated the diversity of the viral population in the individual at ART initiation using a detailed model of viral dynamics and evolution that considered target cells, free virions, and productively infected cells, and included mutation, cellular superinfection, recombination, and fitness selection based on the relative fitness of dominant VRC01 resistance mutations. From the diversity, we obtained the frequencies of productively infected cells containing different mutant proviruses resistant to VRC01. We let the latent reservoir harbor proviruses with the same frequencies, as has been done previously [55]. We recognized that of the latter cells, those harboring the most frequent mutant provirus were the most likely to re-establish infection. We employed a stochastic model to estimate the waiting time for the reactivation of such latent cells and tracked the ensuing dynamics during VRC01 therapy until the infection grew to detectable levels, at which point the therapy was deemed to have failed. During the simulations, we let the efficacy of VRC01 vary continuously with time based on its pharmacokinetic profile, which we estimated from independent fits to data. The parameter values employed for estimating the proviral frequencies and latent cell reactivation are listed in Tables 1 and 2, respectively. Next, we constructed a large virtual patient population based on inter-patient variations in the size of the latent reservoir and the fitness of mutant viral strains, to reflect variations in host and viral factors, respectively, that could influence the outcomes of therapy (Fig 1B). We applied our model above to each virtual patient and estimated the time of therapy failure. From these simulations, we estimated the distribution of remission times in the population and constructed Kaplan-Meier survival plots. We used data from one clinical trial to estimate the parameters defining the inter-patient variations in the virtual patient population and used them without adjustable parameters to predict the outcomes of another clinical trial, validating our model and the parameter estimates.

Frequencies of VRC01-resistant mutants pre-existing in the latent reservoir
To estimate the frequencies of VRC01 resistant proviruses that may exist in the latent reservoir and cause VRC01 failure, we considered viral evolution before the initiation of ART in a chronically infected individual. During ART, viral replication is quickly halted [5,6], leaving little scope for further viral diversification. We considered four mutations in the HIV-1 envelope region reported to be highly resistant to VRC01: N279K, N280D, R456W and G458D [28]. Resistance could come from strains carrying these mutations singly, in pairs, in triplets, or all together. We considered all these strains in addition to the wild-type, or VRC01 sensitive, strain in our model. The frequencies of the strains would depend on their relative fitness, which, following previous studies [56][57][58], is defined in our model by two components: infectivity and replicative ability. The relative infectivity of each of the strains involved has been estimated in independent experiments [28], which we employed (S1 Table). The replicative abilities can be estimated using competitive growth assays [59], which have been reported for some of the strains [28]. From the assays, the two fittest single mutant strains, N279K and N280D, appeared to have replicative fitness not significantly different from the wild-type [28]. We analysed data available from the assays for the other two single mutants, R456W and G458D, and the quadruple mutant [28] using a previously developed formalism [59] and estimated their replicative fitness (S1 Fig). For all the double (except N279K-R456W, which had replicative fitness similar to the wild type [28]) and triple mutants, as an approximation, we set the replicative fitness to values predicted assuming zero epistasis. (This assumption is not critical to our findings, which, as we show below, depend primarily on the single mutant frequencies.) The fitness values are in S1 Table. The fitness values were consistent with recent deep mutational scanning analyses, which suggest that most mutations in the HIV-1 envelope are under purifying selection [60]. With these fitness values and all other parameters representative of chronic HIV-1 infection (Table 1), we estimated the frequencies of productively infected cells harbouring the different mutant proviral genomes (Fig 2).
Our model predicted that the fittest mutant strain, N279K, would exist at a frequency of approximately 2.3×10 −4 , the highest among the different mutants (Fig 2A). Other single mutants were at frequencies~5-fold lower. The double mutants had frequencies of~10 −9 ( Fig  2B), triple mutants of~10 −14 (Fig 2C), and the quadruple mutant of~10 −19 (S1 Table). These frequencies are similar but not identical to those expected from the mutation-selection balance, where a strain with n mutations would exist at a frequency of~(μ/a) n , with a the constant per mutation fitness penalty and μ the mutation rate [49]. For N279K, for instance, where a~0.14, considering both relative replicative fitness and relative infectivity (S1 Table), and with μ~3×10 −5 [52], the mutation-selection balance would yield a frequency of 2.14×10 −4 , slightly lower than that predicted by our model.
Following previous models, where a constant fraction of infection events is assumed to lead to latency [39,42,47], the frequencies of the mutants are expected to be similar in productively and latently infected cells. The latent reservoir is not affected by ART directly, and even in the absence of active viral replication, the latent reservoir would decrease in size extremely slowly, taking years [7]. Thus, one could safely assume that the latent reservoir at the initiation of ART would exist nearly intact post-ART and at the start of VRC01 therapy, as also suggested by recent observations [61]. We do recognize that the latent reservoir can harbour genomes from early in the infection [61], so that the estimate diversity in the plasma would be an upper bound on the diversity in the latent reservoir. Current technologies have been unable to We used a model of within-host HIV-1 evolution to estimate the pre-ART frequencies of wild-type and VRC01-resistant mutants (left), letting them be identical in productively and latently infected cells. ART eliminates the former but not the latter cells (middle). We then used a stochastic model of latent cell reactivation and viral growth to estimate the time of virological failure following VRC01 therapy (right) for a given size of the latent reservoir and the fitness of the VRC01 resistant strain. (B) Dynamics at the patient population level and the outcomes of clinical trials. We created a virtual patient population by sampling the latent cell pool size and mutant viral fitness during VRC01 therapy from defined distributions (left). For each individual, we performed stochastic simulations as in (A) and estimated the time to virological failure (middle), from which we obtained the distribution of breakthrough times and Kaplan-Meier survival plots (right), which we compared with clinical data. estimate the diversity in the latent reservoir comprehensively [36]. We therefore employ the upper bound as a conservative estimate. Given the latent reservoir size, L 0 , of~10 5 -10 8 cells in chronically infected individuals [35,39], the expected number of cells infected with the N279K mutant proviruses would be~23-23000. The corresponding numbers would be~10-5000 for the other single mutants (Fig 2). For the double and triple mutants, however, the numbers would be <0.1 and <10 −6 , respectively. On average, thus, our calculations predicted that most latently infected cells would carry wild-type, or VRC01-sensitive proviruses. A small number, 10-10 4 cells, would carry single mutants resistant to VRC01. Cells carrying higher mutants were unlikely to exist. The single mutants, too, may not be detectable in most cases, given current assay limits (see Introduction).
With this description of the frequencies of VRC01 resistant proviruses in the latent cell reservoir, we examined the dynamics of VRC01 failure due to latent cell reactivation. & Estimates indicate a value much larger than 50 μg/mL [25,28], but a precise value is lacking. # 95% confidence limits are in brackets.

de novo generation of mutants during VRC01 therapy
From the above estimates, the lower end of the spectrum of latent cell numbers carrying single mutant proviruses,~10 cells, is small enough that it is possible that in some individuals, due to stochastic variations in the mutant frequencies and/or infected cell numbers, the latent reservoir contains no resistant proviruses. In such a scenario, de novo mutation after the reactivation of latent cells carrying VRC01 sensitive proviruses would have to give rise to resistance during VRC01 therapy. Even in the presence of latent cells carrying single mutant proviruses, the large majority of cells carrying VRC01 sensitive strains may result in de novo mutations following reactivation of latent cells carrying VRC01 sensitive proviruses being the predominant mechanism of VRC01 failure. We developed a model to test this possibility (S1 Text) and found that the estimated virological breakthrough times were far larger (over~100 days) than those observed clinically (~20 days) (S2 Fig). Although reactivation of such cells was more frequent, given their larger numbers, than the cells carrying mutant proviruses, such reactivation did not lead to lasting infection in our predictions because VRC01 successfully neutralized the viruses produced. Mutations, being intrinsically rare, did not lead to rapid enough de novo development of resistance. The VRC01 failure seen clinically was thus unlikely to be due to the reactivation of latently infected cells carrying VRC01 sensitive proviruses. The more likely mechanism therefore was the reactivation of latently infected cells carrying single mutant proviruses resistant to VRC01. We estimated remission times based on the latter mechanism next.

Growth of pre-existing mutants during VRC01 therapy
We focussed on latently infected cells carrying proviruses with the N279K mutation, which, with their 5-fold higher prevalence than other single mutants, were the most likely to be reactivated. Stochastic simulations (Methods) with a constant VRC01 efficacy against the mutant, indicated that the latent cell pool did not vary significantly over the durations considered ( Fig  3A), consistent with experiments [36]. Reactivation leading to the growth of productively infected cells carrying proviruses with the N279K mutation occurred over a duration of a few days to weeks ( Fig 3B). Detectable viremia, however, took longer given that the viral levels had to rise to 20 copies/ml ( Fig 3C). Defining the time for viremia to become detectable as the time of the failure of VRC01 therapy, or the breakthrough time, the simulations yielded a distribution of breakthrough times ranging from 20-100 d when L 0 was 10 6 cells ( Fig 3D), consistent with observations [17], suggesting that virological breakthrough was likely to be due to the reactivation of cells infected latently with single mutant VRC01-resistant strains. The breakthrough time depended on the time for latent cell reactivation as well as for the subsequent establishment of successful infection. Thus, higher L 0 , which decreased reactivation times, led to earlier breakthrough ( Fig 3E). Increasing the viral production rate (Fig 3F), lowering drug efficacy (( Fig 3G) or increasing viral fitness (Fig 3H), which improved the chances of establishment of successful infection after reactivation, all led to more rapid virological breakthrough. These factors are likely to vary across individuals. We examined next how their influence would manifest in clinical trials and whether our simulations could capture the VRC01 failure seen in the trials.

Multimodal distribution of breakthrough times in a virtual patient population
We focused here on the A5340 trial [17] where VRC01 therapy was initiated a week before the end of ART on patients with well-controlled viremia. 3 doses of VRC01 were administered, with a gap of 3 weeks between successive doses. Virological breakthrough was detected when the viral load crossed 20 copies/ml. To mimic the trials, it was necessary not only to account for inter-patient variations but also to consider time-varying efficacies of drugs within an individual, which can significantly affect the development of drug resistance [51,62]. We therefore first considered VRC01 pharmacokinetics and fit a model of biphasic decay to data of the VRC01 plasma concentration profile following a single intravenous dose [26] (Fig 4A inset). Using the resulting best-fit parameters, we predicted the multiple dose pharmacokinetics for the dosing protocol above and estimated the time-varying efficacy, ε m (Fig 4A), which we employed in our stochastic simulations.
Next, we created a virtual patient population to mimic inter-patient variations expected in clinical trials. A number of host factors, including HLA types, are known to affect the ability of individuals to control viremia [63]. Similarly, viral factors have also been argued to determine the level of viremia in chronic infection [64]. A convolution of host and viral factors is expected to determine clinical outcomes. The specific factors involved and how they vary across individuals, however, is not fully established [46,63]. Here, we employed a parsimonious approach where we let two parameters, one reflective of variations in host factors and the other viral factors, define the inter-patient variations. We thus constructed a virtual patient population with different initial latent pool sizes, L 0 , subsuming variations in all host-factors, and viral production rates, p, subsuming variations in all viral factors. For each individual, we sampled L 0 and p independently from pre-defined distributions (Methods) and ran a stochastic simulation to describe the ensuing dynamics.
The simulated dynamics showed virological breakthrough times varying from a few days to a few months post cessation of ART across individuals ( Fig 4B). Following breakthrough, the viral load rose sharply. It was suppressed partially following the administration of a VRC01 dose, but then rose again once the VRC01 level fell. The patterns were similar to the viral load resurgence patterns seen in patients [17]. From the breakthrough data, we estimated the distribution of virological breakthrough times in this virtual patient population ( Fig 4C). We found that the distribution was multimodal. Until a week or so following the cessation of ART, no breakthrough was expected based on the distribution because of high VRC01 levels in circulation. As the VRC01 concentration declined, breakthrough began. The distribution of breakthrough times peaked when the VRC01 concentration was at its trough level, just before the administration of the second VRC01 dose. Following dosing, the distribution dropped steeply, and rose again as the VRC01 concentration waned. It then attained a peak that was smaller than the first peak and began to fall subsequently. The smaller size of the peak and the subsequent fall was due to a convolution of the effect of VRC01 and the natural distribution of breakthrough times in the absence of treatment. Post ART cessation, it has been shown that the distribution of mean rebound times is unimodal and declines following its peak [41]. In other words, after the peak, the population of individuals that suffers breakthrough decreases with the breakthrough time. In our simulations, this decline is what explains the smaller second peak compared to the first and the drop in the distribution after the second peak although the VRC01 levels were low. Of course, with the third dose, the distribution dropped sharply again due to the rise in the VRC01 level and then rose as the VRC01 level waned. The distribution attained a third peak that was even smaller than the second peak and ended in a long tail representing the small fraction of individuals who experienced longer remission times than studied in our simulations.
Based on the distribution, we constructed a Kaplan-Meier survival plot, which at each time point marked the percentage of the virtual patient population that was still under remission, defined as viremia <200 copies/mL [17]. The plot, as expected, indicated no failure for a short period,~1-2 weeks, after ART, then dropped sharply, reaching 50% failures in about 4 weeks, and displayed a long tail with a small percentage,~10%, maintaining remission for longer than the duration of our simulations (100 days) (Fig 4D). We examined next whether these predictions could recapitulate clinical observations.

Model calibration to recapitulate the A5340 trial
Comparing our simulations with data required knowledge of the distributions of L 0 and p in patients. We chose p to mimic viral growth rates after rebound. The growth rate can vary from 0.4-1.5/day [39,42]. We therefore chose the mean value of p to yield a growth rate of 1/day. Further, we let p follow a normal distribution with a standard deviation that ensured positive growth rate across two standard deviations from the mean. The resulting distribution of p is  [26] (symbols) (Inset). Best-fit parameter estimates are in Table 2. Corresponding multiple dose concentration profiles (blue) and the VRC01 efficacy against a mutant strain with IC 50 = 800 μg/mL (green) are shown. Black arrows indicate VRC01 infusions. (B) Stochastic realizations of the dynamics in a virtual population of 10000 patients, manifesting as changes in plasma viremia. Each grey line represents a virtual patient. Some randomly chosen trajectories are colored to aid visualization of the dynamics. Note that ART was continued after the first infusion of VRC01 for 1 week. The detection limit of 20 copies/mL is marked as a red dashed line, crossing which marks clinical viral rebound. (C) The corresponding distribution of rebound times (orange). Rebound times of the participants in the A5340 trial with a 1-week uncertainty period, representing the gap between successive viral load measurements, are also marked. (D) Kaplan-Meier plot for the A5340 trial based on the percentage of patients with viremia �200 copies/mL. Each light purple line is a survival curve generated by randomly choosing 20 patients from the virtual population above. The dashed purple line is the mean of all these survival curves. Each light green line is an analogous survival curve if the failure were to occur by the recrudescence of latently infected cells infected by wild-type (VRC01 sensitive) proviruses. (IC 50 for the wild-type was 1 μg/mL [25].) The green dashed line represents the average of the latter survival curves. The blue dashed line marks the average survival curve due to resistant strains alone. It is indistinguishable from the purple line until~100 days, indicating that the predominant mode of failure is via resistance. The data from the trial [17] is shown as a red solid line.
https://doi.org/10.1371/journal.pcbi.1008434.g004 mentioned in Table 2. Although estimates of the variations in L 0 exist [35,39], in our simulations, they had to accurately mimic the variations in the pool carrying resistant proviruses, which are not known. We therefore adopted the following approach. We decided to employ data from the A5340 trial to estimate parameters characterizing the distributions of L 0 and then validate them using an independent trial, the NIH trial. Both the trials involved small sample sizes,~10-15 patients [17]. Non-linear mixed effects modeling, designed particularly to estimate parameter distributions using clinical data from small sample sizes, works with deterministic but not stochastic models [65,66]. Consequently, we adopted a heuristic approach to estimate the distributions of L 0 . We recognized that the product a c L 0 , with a c the latency reactivation rate, determines the waiting time for viral recrudescence; the larger the product, the earlier would be the reactivation of the latent reservoir. We fixed a c based on previous estimates [39,47]. (Note that previous studies report a wide range for a c [37,39,42].) Through small test simulations, we identified approximate values of L 0 that mimicked the mean waiting times seen in patients in the A5340 trial. We then performed more detailed simulations by varying the distribution of L 0 around the approximate parameters and identified those distributions that best described the Kaplan-Meier survival data from the A5340 trial (S3 Fig). The resulting parameters are also in Table 2. Similarly, we explored the implications of variations in the IC 50 of VRC01 against the mutant and found that values �800 μg/mL captured the data well (S4 Fig). (We recognize that IC 50 values of rebound virus in individual patients were lower [17]. These could be due to differences between in vitro and in vivo estimates [67], possible evolution post-breakthrough in response to autologous antibodies [36], and/or inter-patient variations (S5 Fig).) Simulations with the resulting distributions recapitulated data from the A5340 trial. We found that 11 of the 12 patients who experienced treatment failure despite high concentrations of VRC01 in circulation [17] had breakthrough times close to the peaks in the multimodal distribution of breakthrough times we predicted (Fig 4C). The number of patients with breakthrough times associated with the different peaks was also proportional to the area under the peaks. Note that the area under a peak is a measure of the probability, and hence the frequency, of failure corresponding to the times associated with the peak. 6 patients failed during the times associated with the first peak, 4 with the second peak, and one with the third peak. The areas under the peaks from our simulations yielded failure percentages of~35%,~44%, and 21%, respectively. The first two peaks, thus, appeared to have comparable failure percentages, whereas the third was substantially smaller. Given the small sample size in the clinical trial, the distribution of patients into the three peaks appeared to be consistent with the estimated failure percentages. One patient (A01) appeared to fail at the trough in the distribution after the first peak, and this could be due to stochastic effects or variations not captured in our virtual population.
Kaplan-Meier plots based on breakthrough times, i.e., times for the viral load to reach 200 copies/ml, from the virtual patient population also closely recapitulated the clinical data ( Fig  4D). Here, to account for the small sample size in the trials, we chose many samples of 20 individuals each, selected randomly from our virtual patient population, and constructed Kaplan-Meier curves for each sample. In addition to failure due to VRC01 resistant strains, for each virtual patient we also estimated the breakthrough time due to the wild-type, or VRC01 sensitive, strains (S1 Text), and set the breakthrough time of the individual to be the smaller of the two. The clinical data fell within the ranges defined by these curves (purple lines in Fig 4D). Further, the mean of these curves was in close agreement with the data. Accordingly, 50% of the treated population exhibited a breakthrough time of >4 weeks from the end of ART (or 5 weeks from the start of VRC01 therapy), consistent with the clinical data. Importantly, the simulations predicted that the failure was predominantly due to VRC01 resistant strains. Ignoring failure due to the wild-type strains made little difference to the survival plot (compare blue and purple lines in Fig 4D). The failure due to the wild-type strains became important only after VRC01 levels had declined significantly (green lines in Fig 4D). The simulations tended to marginally over-predict the clinical data for long remission times; whereas all the patients had shown rebound by 80 d, the simulations predicted that~10% patients would retain control still and take much longer (up to~120 d) to experience failure. We attributed this to the presence of individuals with strong immune responses and/or small latent cell populations, including post-treatment controllers [63], who may not be seen in the small population of 13 individuals in the A5340 trial. Notwithstanding these deviations, our simulations closely recapitulated the data from the A5340 trial.

Model validation with the NIH trial
To test and validate our model and the parameter estimates, we applied our simulations to describe a second, independent clinical trial, the NIH trial [17], where 10 individuals were subjected to VRC01 therapy during an ATI. Treatment commenced 3 days before ART cessation. Subsequent doses were administered on weeks 2 and 4 and then every month until 6 months. To describe the resulting breakthrough data, we created a virtual population exactly as above and subjected it to therapy following the clinical protocol. VRC01 pharmacokinetics was predicted based on the corresponding dosing times (Fig 5A). All the parameters were kept the same as those as in our simulations of the A5340 trial.
We found that viral rebound trajectories were similar to those in the A5340 trial, with sharp rises in viral load following breakthrough and wide inter-patient variations (Fig 5B). The distribution of breakthrough times again followed a multimodal distribution with the peak widths broadening progressively with each dose and culminating in a long tail (Fig 5C). Based on the areas under the peaks, the percentages of failure at times corresponding to the four peaks were~16%,~39%,~34%, and~9%, respectively. Indeed, in agreement, 7 of the 10 patients failed at times corresponding to the second and third peaks. 2 patients were associated with the first peak and one with the last. More precise timing of failure of the patients was not possible given the much larger intervals between successive viral load measurements in the NIH trial compared to the A5340 trial. Nonetheless, the agreement between our simulations and the observations of the distribution of failure times between peaks, given the small sample size, was remarkable. We also computed Kaplan-Meier survival plots, using samples of 20 virtual patients from our virtual population, and considering failure due to VRC01 resistant strains. (Failure due to the wild-type was unlikely, as observed above.) The Kaplan-Meier plots captured the clinical data quantitatively (Fig 5D), indicating that our simulations accurately mimicked the response of patients to VRC01 therapy in the NIH trial.
The difference between the A5340 trial and the NIH trial was in the dosing protocol alone. Indeed, once our model was calibrated with data from the A5340 trail, it captured the data from the NIH trial by simply changing the dosing protocol accordingly, and required no other adjustments, thereby providing a strong test and validation of our model. The model may be applied to predict the outcomes of other possible dosing protocols, which could involve changing dosages or half-lives (S6 Fig), providing a framework for rational therapy optimization.

Discussion
Passive immunization with HIV-1 bNAbs holds promise as a strategy to achieve long-term remission of HIV-1 infection [20,68]. Following its success in macaques [20,21], enormous efforts are underway to translate it to humans [69]. Trials with VRC01 administered to chronically infected patients following the cessation of successful ART saw rapid virological failure despite the presence of suppressive concentrations of VRC01 in circulation [17], signaling a potential vulnerability of such bNAb therapies. Here, using mathematical modeling and analysis of clinical data, we elucidated the likely cause of this rapid VRC01 failure. Model predictions attribute this failure to the reactivation during therapy of cells latently infected with VRC01 resistant proviruses before ART initiation. To arrive at this inference, we constructed a mathematical model that integrated within-host viral evolution, latency reactivation, and viral dynamics with VRC01 pharmacokinetics, and applied it to simulate the outcomes of therapy in a virtual patient population. Our simulations recapitulated data from clinical trials [17], giving us confidence in the inference. Accounting for pre-existing resistant strains in the latent reservoir would be important to the success of bNAb-based therapies.
That mutation-driven resistance can be an important cause of bNAb failure has been recognized earlier [28][29][30][31]. Indeed, efforts are ongoing to identify bNAbs, including in the VRC01 class, that are less vulnerable to such resistance [70]. For bNAb therapies that target maintenance of HIV remission post ART cessation during chronic infection, our study prescribes the requisite genetic barrier to resistance, i.e., the number of mutations HIV-1 must accumulate to develop significant resistance to the therapy. Using a detailed description of within-host HIV-1 evolution, our study predicts that with the estimated latently infected reservoir size of 10 5 -10 8 cells in a chronically infected individual [35,39], proviral strains carrying single but not more resistance mutations would pre-exist in the latent reservoir. Consequently, a genetic barrier of 2 would be necessary to prevent pre-existing resistance from compromising therapy. A bNAb with a genetic barrier of one, like VRC01, would be predestined to fail, as was observed in clinical trials [17]. When used in combination with another bNAb, however, the overall barrier would cross the threshold of 2, diminishing the chances of such failure. Resistance would then have to develop by de novo mutation during therapy, which according to our model would be unlikely as long as adequate bNAb concentrations exist in circulation and restrict the replication of the wild-type virus. Further, by explicitly incorporating bNAb pharmacokinetics, our model provides a framework with which optimal dosing protocols could be identified that would ensure adequate bNAb concentrations throughout.
Early ART initiation has been argued to restrict the viral reservoir size and improve the chances of post-treatment control [10,47]. Our study predicts an additional advantage of early ART initiation: The chances of failure due to resistance are reduced. This reduction happens in multiple ways. First, most HIV-1 infections involve a single founder virus, which gradually evolves into the diverse quasispecies seen in chronic infection [54,71,72]. Thus, the frequency of mutants at the time of ART initiation is expected to be smaller, the earlier the initiation [54]. Second, given the limited viral diversity early in infection, the latent reservoir is likely to have a much higher representation of bNAb-sensitive strains than estimated using our model of chronic infection. Finally, if the reservoir size is restricted, due to early ART initiation, far fewer mutant proviruses, or none at all, may exist in the latent reservoir. Together, thus, early initiation of ART is expected to reduce the chances of bNAb failure. Indeed, in trials where ART was initiated early, in the acute phase of infection, no genotypic resistance to VRC01 was  [73]. Similarly, a combination of two bNAbs, including one in the VRC01 class, administered early in infection saw no resistance to therapy in SHIV-infected macaques [20]. Surprisingly, however, the breakthrough times for VRC01 therapy following ART in acute infection were found to be similar to those seen in the trials in chronic infection we studied [17,73]. Additional mechanisms that are not predominant in the chronic phase thus appear to cause VRC01 therapy failure in individuals treated with ART early. One possibility could be that CD8 T cell exhaustion is weaker and/or more reversible in the acute phase than in the chronic phase because of the much longer duration of antigen exposure in the latter scenario [74,75]. The greater associated immune activation levels in the acute phase could imply more rapid reactivation of latently infected cells, which could offset the advantage from the lower frequencies of mutants. To test this possibility, models that incorporate CD8 T cell exhaustion [39,[76][77][78][79] would have to be integrated with our model of stochastic latency reactivation, a promising avenue for future study.
An important question in HIV cure research is how early should ART be initiated to maximize post-treatment control [63]. While starting early would restrict the reservoir size and improve the chances of post-treatment control, starting it too early would not allow enough time for the development of an immune response, compromising post-treatment control. If bNAb therapy were to be used post-ART to improve the chances of post-treatment control, the development of resistance to bNAbs would have to be factored in along with the latter trade-off between reservoir size and immune response strength to arrive at the optimal timing of ART initiation. Our study provides a framework that could be used to test whether preexisting resistance in the latent reservoir would compromise therapies with other bNAbs, especially those belonging to the VRC01-class [80][81][82], which are under trial for both preventive [24,83,84] and therapeutic [20] vaccination and are known to fail via diverse mutation-driven resistance pathways [29].
A modeling study several years ago compared the likelihood of the failure of therapy due to pre-existing resistance versus de novo generation of mutants in the context of ART and found that ART was more likely to fail due to pre-existing resistance [50]. A more recent modeling study examined the likelihood of the failure of different antiretroviral drug combinations and explained how patient adherence influences such failure [55]. Our findings are consistent with these studies and show that bNAbs too are vulnerable to pre-existing resistance, but the resistance now is restricted to the latent reservoir. Because the reservoir is small compared to the pool of infected cells pre-ART, the required genetic barrier for bNAb therapies is estimated to be lower than for ART. Thus, a combination of 2 bNAbs is predicted to overcome resistance, whereas first-line ART necessarily contains 3 drugs.
Previous studies have used either fully deterministic or fully stochastic frameworks to describe treatment failure [50,51,55], making their predictions approximate or computationally expensive, respectively. Here, we devised a strategy that retained both accuracy and computational tractability. We used a deterministic framework to estimate the frequencies of mutant proviral genomes in the latent reservoir pre-treatment and then a stochastic framework to estimate latency reactivation and treatment failure times. The deterministic framework was shown previously to agree well with stochastic population genetics-based simulations [52], giving us confidence in the strategy. In the stochastic simulations, we considered only latently infected cells containing the dominant resistant provirus, which accurately described the development of drug resistance without requiring expensive computations. Indeed, we were able to capture clinical trials of VRC01 failure with this hybrid framework, reiterating its applicability. Such hybrid deterministic-stochastic frameworks have been successfully applied in other settings, such as in predicting the pre-existing frequency of hepatitis C virus strains resistant to drugs [85], but not, to our knowledge, with HIV-1 infection.
The sample sizes involved in the clinical trials we studied were small,~10-15 patients each [17]. Nonlinear mixed-effects models have been used in recent studies to infer the effects of interventions by analyzing data from such trials [86][87][88]. For instance, bNAb therapy was argued to improve immune responses 8-fold over controls and not to synergize with the TLR7-agonist [86]. The strength of such an approach lies in its robust parameter estimation and the ability to infer effects at the population level using data from small sample sizes. The approach, however, does not work when the underlying model is stochastic, which is the case in our study. Also, when the effects are explicitly parameterized, as is often done [86,88], their quantitative estimates are restricted to the specific conditions studied. Thus, for instance, how a change in the dosage or dosing protocol would alter the effect becomes difficult to predict. Our approach, being fully mechanistic, is not similarly limited. Indeed, with the parameters that captured the A5340 trial, by simply changing the dosing protocol, our model captured data from the NIH trial without any adjustable parameters. Our model could therefore be used potentially to comparatively evaluate alternative treatment protocols and suggest optimal ones.
We recognize that in the trials we considered [17], the dominant mutant seen post treatment failure was not the same across individuals. The mutations, however, were all typically concentrated in the same genomic regions [17], suggesting that structural or conformational modifications that could drive VRC01 resistance could be produced by many mutations in the same genomic region. While we have focused on the most frequent mutant, as it is the most likely to have caused resistance, and as has been done in earlier studies [55], stochastic variations could result in latent cells carrying other mutants being reactivated and reestablishing infection. The mutants selected could also differ based on inter-host variations and the genetic backgrounds of the infecting strains. Our simulations must thus be viewed as reflecting therapy failure arising from the dominant mutant which could differ across hosts and which in our study is accounted for by the inter-host variation in viral fitness.
bNAbs are known to engage the immune system via multiple mechanisms [20,87,[89][90][91]. The result could be heightened activation, which for HIV-1 infection, could mean greater susceptibility of target CD4 T cells. At the same time, it could imply greater reactivation rates of latently infected cells. Indeed, the latency reactivation rates employed in our study were higher than those estimated for historical controls, the latter based on viral recrudescence post ART and in the absence of further intervention [39]. When bNAb therapy is administered post ART, cells latently infected with wild-type strains, which would be in a vast majority, would get frequently reactivated and produce virions but without causing sustained infection. The virions produced, being bNAb sensitive, would get neutralized and cleared by the bNAbs. Previous studies have argued that in the process bNAbs can stimulate CD8 T cells [20,87], improving the overall immune activation status, possibly explaining the higher latency reactivation rate we required to describe VRC01 failure than previously used to describe historical controls [38]. Our parameters would thus tend to underpredict the remission times in the historical controls. Previously, too, differences in the viral growth rates before and after ART have been assumed in order to capture viral recrudescence accurately and have been attributed to different strengths of the immune responses in the respective periods [39]. Mechanisms that could explain the differences in these latency reactivation rates are yet to be identified. A unified framework that describes remission both with and without bNAb therapy awaits future studies that would quantify how bNAbs influence immune responses.
Non-nucleoside reverse transcriptase inhibitors (NNRTIs), administered as part of ART, have been found to delay breakthrough post ART [23]. Patients administered NNRTIs were excluded from the A5340 trial and were switched to an integrase-inhibitor based regimen 2 weeks before VRC01 therapy in the NIH trial to eliminate the confounding effects of NNRTIs. We therefore did not consider the effect of NNRTIs in our study. A number of host [63] and viral [64] factors are thought to be involved in determining disease progression and treatment outcomes. Our model subsumed inter-patient variations in these factors into variations in two factors, the latent pool size and the viral production rate. Virtual patient populations that we created based on variations in these minimal factors captured clinical data from clinical trials, justifying the approximation, and suggesting that a small subset of factors may be adequate to capture outcomes of such bNAb therapies. Future studies may consider the variations in other factors explicitly to ascertain the validity of the approximation, especially for other kinds of bNAb-based interventions.
In summary, our study presents an explanation of the failure of VRC01 therapy to sustain HIV remission post ART, captures clinical data quantitatively, highlights the importance of accounting for pre-existing resistance in designing effective bNAb-based therapies, and facilitates rational optimization of such therapies.

Methods
We considered the scenario where a chronically infected individual maintains viral suppression with ART and is then subjected to VRC01 therapy concomitantly with cessation of ART. Because active viral replication is halted by successful ART, virological breakthrough during VRC01 therapy must arise from the reactivation of latently infected cells. We developed a model to estimate the timing of the failure of VRC01 therapy via the growth of resistant viral mutants pre-existing in the latent reservoir (Fig 1). Our approach was to estimate the population of latently infected cells carrying the resistant strains and then to follow their reactivation leading to successful infection. We then considered a virtual patient population subjected to the same therapy and applied our model to recapitulate data from clinical trials.

Mathematical model of virological breakthrough in a single infected individual
Within-host HIV-1 evolution and the frequencies of pre-existing resistant strains. To estimate the population of cells latently infected with VRC01-resistant proviral genomes, we reasoned that the frequencies would be the same as in productively infected cells before treatment initiation because the probability that a particular cell becomes latently infected is not known to depend on the nucleotides at the VRC01 resistance loci. We estimated the frequencies of VRC01-resistant strains in productively infected cells using an approach developed previously to quantify the pre-existence of resistance to antiretroviral drugs [52]. The approach combines virus dynamics with evolution, incorporating mutation, cellular superinfection, recombination, and fitness selection. We present details below.
Virus dynamics. We considered n = 4 positions in the env region where mutations with high level resistance to VRC01 have been identified, namely N279K, N280D, R456W and G458D [28]. The time evolution of the populations of cells and virions containing the different single, double, triple, and quadruple mutants were described using the following equations: Here, uninfected target CD4+ T cells, U, are produced from the thymus at the rate λ and are lost at the per capita death rate d T . They are infected by virions V jh containing genomes j and h with the second order rate constant k 0 β jh , where k 0 is the infectivity of wild-type virions and β jh is the relative infectivity of virions V jh . We assumed that β jh = (β j +β h )/2, where β j is the infectivity of genome j relative to the wild-type. The indices j and h denoting viral genomes range from 0 to S = 2 n −1. Thus, S = 15 here, with the different genomes including the wild-type (or sensitive strain), the 4 single mutants, 6 double mutants, 4 triple mutants, and the quadruple mutant, amounting to a total of 16 types. We denoted the genomes serially, starting with '0' for the wildtype and ending with S for the quadruple mutant. Each viral particle contains two genomes, not necessarily identical. The number of distinct viral particle types, V jh , where j2{0,1,..,S} and h2{j,j +1,..,S}, is thus (S+1)(S+2)/2, which here would be 136. (The range of values of h is to ensure that virions V 12 and V 21 , for instance, which are identical, are not counted separately.) Summing over all these viral types yields the total rate of loss of U due to infections in Eq 1.
Following infection with a virion V jh , reverse transcription, which includes mutation and recombination, yields genome i with a probability Q i (jh), where i2{0,1,..,S}. On average, thus, a fraction Q i (jh) of the infections with virions V jh yield productively infected cells carrying single proviruses i, which we denote T i . Summing over all V jh yields the total rate at which cells T i are produced. These cells die at the per capita rate δ. They can also be infected again, but with a lower infectivity k 1 because infected cells downregulate their CD4 receptors, rendering further infections difficult [92,93]. The net effect of these processes defines the dynamics of cells T i in Eq 2.
Doubly infected cells T ii are produced when cells T i are infected with virions V jh , following which reverse transcription again yields the provirus i. When a different provirus j(6 ¼i) is produced, the result is the doubly infected cell T ij carrying distinct proviruses i and j. Of course, cells T ij can also be produced by the infections of cells T j with another virion yielding the provirus i. Doubly infected cells too die with the rate constant δ. These processes are contained in Eqs 3 and 4. We neglected cells infected more than twice, following experiments that suggest a rare occurrence of such cellular superinfection [94].
Cells T i and T ii produce virions V ii at the per capita rate pξ i , where p is the production rate of wild-type virions and ξ i the production rate of virions V ii relative the wild-type. ξ i is thus the relative replicative fitness of genome i. Cells T hi (h6 ¼i) can also produce virions V ii . We assumed that viral RNA of the types h and i are present in cells T hi in proportion to ξ h and ξ i , respectively, and that they are randomly assorted into pairs and packaged into progeny virions. Thus, cells T hi produce virions V ii at the rate and V ih at the rate pξ h ξ i /(ξ h +ξ i ). Virions are cleared at the rate c. These processes determine the dynamics of viral populations in Eqs 5 and 6. Mutation and recombination. We next describe the formalism to compute the probability Q i (jh). Following previous studies [51,52,95,96], we let R k (jh) be the probability that genome k is produced by the recombination of genomes j and h, and P ik the probability that genome k mutates to genome i. Thus, P ik R k (jh) is the probability that genome i is produced from genomes j and h via the intermediate k. We recognize next that if genomes j and h differ in d positions, then recombination could produce a total of 2 d different genomes, depending on whether at each of the d positions, the nucleotide chosen is either from genome j or genome h. Summing over these different intermediates k yields the total probability of producing genome i from genomes j and h during reverse transcription: To compute R k (jh), we consider the desired path of the enzyme reverse transcriptase on the two genomes so that the enzyme is on the appropriate genome, j or h, at each of the d distinctive sites, so that the genome k is produced. We let the separation between the x th and x+1 st distinctive sites be l x . If the enzyme has to be on the same genome (j or h) at both these sites, then it must perform an even number of crossovers in the length l x . Else, it must perform an odd number of crossovers. If the enzyme has a probability ρ of crossover per site, then the probabilities of even and odd crossovers over a length l are P even (l) = (1+(1−2ρ) l )/2 and P odd (l) = (1−(1 −2ρ) l )/2, assuming that the crossover at any position is independent of the others and all crossovers happen with the same probability [95]. If we write P des (x+1|x) as the probability that the enzyme arrives on to the desired genome at the x+1 st distinctive site given that it was on the desired genome at the x th site, then depending on whether the associated crossovers must be even or odd, we write P des (x+1|x) = P even (l x ) or P des (x+1|x) = P odd (l x ). Note that P des (1) = 1/2 because the enzyme could be on either genome at the start of the reverse transcription process. Thus, multiplying the probabilities over all the distinctive sites yields R k (jh): Next, we estimated the probability of producing genome i from genome k. If the two genomes differ at u sites, then genome i is produced from genome k by mutating genome k at the distinctive sites and nowhere else. Thus, where μ is the per site mutation probability and n is the number of sites of interest.
Frequencies. Eqs 1-9 yield a model of viral dynamics that predicts the growth of the populations of different mutants. We solved the equations for their steady states and obtained the corresponding frequencies of all proviruses, ϕ i , contained in productively infected cells: The approximation in Eq 10 is justified by the small doubly infected cell population compared to the singly infected cell population. The parameter values employed are in Table 1. ϕ i yield the frequencies of various VRC01-resistant mutants before the start of treatment. Previous studies have shown that this deterministic formalism yields mutant frequencies in agreement with stochastic population genetics-based simulations of HIV evolution [52]. The frequencies are expected to hold also for the proviruses in latently infected cells. Further, we expect ART not to influence the latter frequencies; standard first-line ART drugs target HIV reverse transcriptase and protease, whereas VRC01 targets the HIV envelope. We assumed therefore that following ART, the frequencies of VRC01-resistant strains contained in the latent reservoir are given by Eq 10.
We focused next on the reactivation of latently infected cells carrying the resistant genomes.
Latency reactivation and viral rebound. We developed a stochastic framework to describe the reactivation of latently infected cells carrying mutant proviruses resistant to VRC01. Because reactivation is most likely of the cells carrying the fittest mutant, which are the most prevalent, we considered cells L m carrying the fittest mutant strain above. Previous studies on ART resistance too have considered the dominant mutant alone [55]. The reactivation of these cells and the ensuing growth of mutant virions is then described by the following events: T m À À À À À!
Here, cells L m proliferate at the per capita rate ρ l (Eq 11), die at the per capita rate d l (Eq 12), and get activated to productively infected cells, T m , at the per capita rate a c (Eq 13). Cells T m carry the mutant provirus and produce resistant virions, V mm , at the per capita rate pξ m , which in turn infect uninfected cells, U, with the rate constant kβ mm . Here, ξ m and β mm are the relative replicative ability and infectivity, respectively, of the strain m. Free virions are cleared at the per capita rate c. Viral production and clearance are typically rapid [97] compared to changes in cell populations, so that following an approximation widely used (e.g., see [50]), V mm can be assumed to be in pseudo-state with T m . Thus, V nm �pξ m T m /c. The rate, kβ mm VU, of the infection of U thus becomes kβ mm pξ m T m U/c. We recognized that U is not altered significantly due to new infections, especially following ART when the viremia is small. We let a fraction f of the new infections lead to latency and the remaining to productive infection. Further, we let a VRC01-resistant strain mutate back to the wild-type with a probability μ, the point mutation rate of HIV-1. The resulting description would capture scenarios where a single mutation is adequate to develop resistance, as is the case with VRC01 [28]. Thus, the rate kβ mm pξ m T m U(1−f)(1−μ)/c becomes the rate of the growth of T m in the absence of therapy. If VRC01 were to block infections with efficacy ε m , the rate would become X m T m (1−f)(1−μ)(1 −ε m ) with X m = kβ mm pξ m U/c. Thus, we let T m (effectively) double with the per capita rate X m (1 −f)(1−μ)(1−ε m ) (Eq 14), and yield new latently infected cells at the per capita rate X m f(1−μ)(1 −ε m ) (Eq 15). Cells T m die with the rate constant δ. These events together describe the growth of resistance to VRC01 and therapy failure. We solved the model equations using the Gillespie algorithm [98] with parameter values representative of HIV-1 infection in the presence of VRC01 therapy. The model was implemented using a program written in MATLAB (S2 Text). The parameter estimates are listed in Tables 1 and 2. The initial population of L m was set by the frequencies of mutants estimated. With each parameter setting, we performed 5000 realizations to examine the dynamics of treatment failure.

Comparison with clinical data
We applied our model to predict the outcomes of trials with VRC01 [17]. For this, we explicitly considered VRC01 pharmacokinetics. Further, we created a virtual patient population to account for the inter-patient variations seen in the trials. In each patient, we performed stochastic simulations and identified the VRC01 failure time following the model above. We compared the resulting distribution of failure times with those observed in the trials. We describe the methods we used here.
Data. We considered data of viral resurgence following VRC01 therapy during ATIs from two clinical trials, the A5340 trial and the NIH trial [17]. In the A5340 trial, 14 adult chronic HIV patients with a median ART duration of 4.7 years and viremia maintained below detection were administered 3 doses of VRC01 (40 mg/kg of body weight) intravenously, starting a week before the cessation of ART and with 3 week intervals. Plasma viremia was measured weekly to check for rebound (>20 copies/ml). Of the 14 participants, 1 participant stopped ART before VRC01 administration and was not part of the trial data analysis. We considered data from the remaining 13 patients. For each patient, we assumed the breakthrough time was anywhere within the week from the last undetectable to the first detectable viral load measurement.
The NIH trial had 10 participants with similar characteristics as the A5340 trial. The participants had been on ART much longer, however, with a median duration of 10 years, before entry into the trial. ART was stopped 3 days after the first VRC01 administration. VRC01 doses, at the same dosage as above, were administered subsequently on days 14 and 28 post ART cessation and once every month thereafter. Plasma viremia was measured every week for 1 month and then every 2 weeks until 6 months. The patients experienced virological breakthrough between 2 and 12 weeks from the discontinuation of ART. We considered this time of virological failure with uncertainties as defined above.
We note that both trials did not screen for pre-existing VRC01 resistance. While the primary endpoints of the trials were safety and tolerance, secondary endpoints were viral remission, based on which the trial data report Kaplan-Meier survival curves [17]. We considered the latter data too in our analysis.
VRC01 pharmacokinetics and pharmacodynamics. We let the efficacy of VRC01 against the drug-resistant mutant strain, ε m , be related to its plasma concentration, A, by the Hill function [51,57,99], where IC 50,m is the value of A at which VRC01 is 50% efficacious. For simplicity, we set the Hill coefficient to 1. The antibody concentration has been observed to decline in a biphasic manner upon dosing [26]. We therefore described the time course of the antibody concentration using the following expression when successive doses are administered at time points ω 1 , ω 2 , and so on: Here, A 1 and A 2 are the portions of a dose that decay in the first and second phases, respectively, with decay rates η 1 and η 2 , and W is the total number of doses. We estimated the pharmacokinetic parameters using fits of the above equation to data of VRC01 concentrations following a single dose [26].
ART efficacy. For the duration that ART is used simultaneously with VRC01, we replaced the term (1−ε m ) in Eqs 14 and 15 by 1−ε comb = (1−ε m )(1−ε ART ), where ε ART is the efficacy of ART and ε comb is the combined efficacy of VRC01 and ART. We let ε ART be constant while on ART and set it to zero thereafter.
Virtual patient population. To capture inter-patient variations in the response to bNAb therapy, we constructed a virtual population of clinical trial participants as follows. For simplicity, we considered variations in two factors, the initial latent pool size, L 0 , and the viral production rate, p, across individuals. We assumed that variations in all host factors could be subsumed in the variation in L 0 and that variations in viral factors could be subsumed in p. Introducing variations in other parameters did not affect our findings, as we found by varying the IC 50 of VRC01 across patients (S5 Fig), justifying the latter assumption. We let L 0 vary lognormally and p normally across individuals. Thus, the corresponding density functions were where μ l and σ l are the mean and standard deviation of lnL 0 , and μ p and σ p are the mean and standard deviation of p, respectively. The parameters in the distributions were chosen based on earlier studies or to fit the results of the A5340 trial. From the resulting distributions, we sampled 10000 pairs of values of L 0 and p, with each pair representing an individual. We thus created a virtual patient population, which we then subjected to VRC01 therapy according to the protocols in the respective trials. For each individual, we ran a stochastic simulation (Eqs 1-18) and examined the dynamics of virological breakthrough. From the resulting dynamics, we constructed the distribution of breakthrough times and used it to build Kaplan-Meier survival plots sampling virtual patients mimicking the clinical trial protocols.
Supporting information S1 Text. VRC01 therapy failure is unlikely to be due to de novo mutation.
(DOCX) S1 Table. Relative replicative fitness, relative infectivity, and pre-existing frequencies of various mutants. Replicative fitness was estimated by analyzing competitive growth assays where available (S1 Fig) and by assuming zero epistasis otherwise. For the N279K and N280D mutants, competitive growth assays indicated replicative fitness indistinguishable from the wild-type. The relative infectivity of each strain was calculated as the ratio of its % infectivity to that of the wild type, reported elsewhere [28]. The frequencies were obtained using our model (Eqs (1)- (10)).
(XLSX) S1 Fig. Estimation of fitness using data from competitive growth assays. We applied a previous formalism [59] to analyze competitive growth assays, where the selective advantage, s, of a mutant relative to the wild-type is given by , the corresponding probability density function calculated using the deterministic formalism (Eqs. S12-S17) is shown as solid lines. (H) Expected waiting time, τm, for the formation of a productively infected cell carrying a VRC01-resistant provirus, calculated using the deterministic formalism as a function of the initial latent pool size for the VRC01 efficacy indicated, which is the efficacy against the wild-type averaged over the dosing interval. Note that the modified bNAb VRC01LS has a >4-fold longer half-life than VRC01, and has been tested in healthy and uninfected adults for safety and pharmacokinetics [100]. The terminal half-life, η 1 , has been varied accordingly. The initial total antibody concentration (A 1 +A 2 ) has been varied while keeping the ratioA 1 /A 2 fixed, with the maximum (A 1 +A 2 ) (in μg/mL) corresponding to the 40 mg/kg dosage of VRC01 (fitted from data [26]), which appears to be the maximum dosage of VRC01 reported. The other parameters are the same as in Fig 4. (TIF)