After infusion of bNAbs in a patient, the antibodies bind and neutralize the susceptible strains of HIV. The neutralized subpopulation of HIV-1 no longer infects T-cells, and the plasma RNA copy-number associated with this neutralized population decays. The dynamics of viremia in HIV-1 patients off ART following a bNAb therapy with 3BNC117 (Caskey et al., 2015), 10–1074 (Caskey et al., 2017), and their combination (Baron et al., 2018 ) are shown in Figure 3—figure supplements 1–3. With competition of the neutralized strains removed, the resistant subpopulation grows until the viral load typically recovers to a level close to the pretreatment state (i.e. the carrying capacity); see Figure 1A. The time it takes for the viral load to recover is the rebound time—a key quantity that characterizes treatment efficacy within a patient. Although the details of the viremia dynamics, especially at beginning and at the end of the therapy, may be complex (Lu et al., 2016; Reeves et al., 2020; Saha and Dixit, 2020; Meijers et al., 2021), the rebound time can be approximately modeled using a logistic growth after bNAb infusion (t > 0),

with the initial condition set for pre-treatment fraction of resistant subpopulation $x=N_r(0)/(N_r(0)+N_s(0))$, where $N_r(0)$ and $N_s(0)$ denote the size of resistant and susceptible subpopulations at time $t=0$, respectively. Here, $\gamma$ is the growth rate of the resistant population, $r$ is the neutralization rate impacting the susceptible subpopulation, and $N_k$ is the carrying capacity (Figure 1A, Methods). In our analysis, we set $\gamma =1/3{\text{days}}^{-1}$ or a doubling time of $\sim 2$ days, the known HIV-1 growth rate in patients (Perelson et al., 1996). We infer the neutralization rate $r$ as a global parameter for each trial, since it depends on the neutralization efficacy of a bNAb at the concentration used in the trial. We will infer the patient-specific pre-treatment fraction of resistant subpopulation $x$, using a population genetics based approach based on which we characterize the mutational target size and selection cost of escape in the absence of a bNAb (see below). The maximum-likelihood fits of $N(t)$ to the viremia measurements in Figures 2 and 3, Figure 3—figure supplements 1–3 specifies the initial resistant fraction x_p and the rebound time $T_p$ in each patient $p$, which in this simple model, is given by $T_p=-{\gamma }^{-1}\log x_p$ (Methods).

Figure 1

Download asset Open asset

Figure 2 with 1 supplement see all

Download asset Open asset

Figure 3 with 3 supplements see all

Download asset Open asset

The rebound times following passive infusion of 3BNC117 (Caskey et al., 2015 and 10–1077 Caskey et al., 2017) bNAbs range from 1 to 4 weeks, with a small fraction of patients exceeding the monitoring time window in the studies (late rebounds past 56 days). The distribution of rebound times summarizes the escape response of the virus to a therapy and directly relates to the distribution for the pre-treatment fraction of resistant variants $P(x)$ across patients $P(T)\sim x^{-1}P(x)$.

The fate of an HIV-1 population subject to bNAb therapy depends on the composition of the pre-treatment population with resistant and susceptible variants, and the establishment of resistant variants following the treatment. To capture these effects, we construct an individual-based stochastic model for viral rebound (Figure 1B). We specify a coarse-grained phenotypic model, where a viral strain of type $a$ is defined by a binary state vector ${\overrightarrow{\rho }}^a=[{\rho }_1^a,\mathrm{\dots },{\rho }_{\mathrm{\ell }}^a]$, with $\mathrm{\ell }$ entries for potentially escape-mediating epitope sites; the binary entry of the state vector at the epitope site $i$ represents the presence (${\rho }_i^a=1$) or absence (${\rho }_i^a=0$) of an escape mediating mutation against a specified bNAb at this site of variant $a$. We assume that a variant is resistant to a given antibody if at least one of the entries of its corresponding state vector is non-zero. For multivalent treatment, a virus must be resistant to all antibodies comprising the treatment to have a positive growth after infusion.

At each generation, a virus with phenotype $a$ can undergo one of three processes: birth, death and mutation to another type $b$, with rates ${\beta }_a$, ${\delta }_a$, and ${\mu }_{a\to b}$, respectively (Figure 1B). The net growth rate of the viral subpopulation with phenotype $a$ is the birth rate minus the death rate, ${\gamma }_a={\beta }_a-{\delta }_a$ (Figure 1C). The total rate of events (birth and death) per virion $\lambda ={\beta }_i+{\delta }_i$ modulates the amount of stochasticity in this birth-death process (Methods), which we assume to be constant across phenotypic variants. The continuous limit for this birth-death process results in a stochastic evolutionary dynamics for the sub-population size $N_a$,

where $\eta (t)$ is a Gaussian random variable with mean $⟨\eta (t)⟩=0$ and correlation $⟨\eta (t)\eta (t^{\prime })⟩=\delta (t-t^{\prime })$ (Methods). Here, f_a denotes the intrinsic fitness of variant $a$ and its net growth rate $\gamma$ is mediated by a competitive pressure $\varphi =\frac{1}{N_k}{\sum }_b{N}_b{f}_b$ with the rest of the population constrained by the carrying capacity $N_k$, such that ${\gamma }_a=f_a-\varphi$. In the presence of a bNAb, birth is effectively halted for susceptible variants and their death rate is set by the neutralization rate of the antibody, resulting in a net growth rate, ${\gamma }_{\text{sus.}}=-r$. At the carrying capacity, the competitive pressure is the mean population fitness $\varphi =\overline{f}$, making the net growth rate of the whole population zero (Figure 1C). When susceptible variants are neutralized by a bNAb, the competitive pressure $\varphi$ drops, and as a result, the resistant variants can rebound to carrying capacity, at growth rates near their intrinsic fitness.

To connect the birth-death model (Equation 2) with data, we should relate the simulation parameters of a birth-death process to molecular observables. We have already made a connection between the birth and death rates of a variant and its intrinsic fitness in Equation 2. In addition, the neutral diversity $\theta$ of a population at steady-state can be expressed as $\theta =2N_k\mu /\lambda$, where $\mu \approx {10}^{-5}/\text{day}$ is the per-nucleotide mutation rate, which we infer from intra-patient longitudinal HIV-1 sequence data (Zanini et al., 2015) (Methods). For consistency, we set the total rate of events $\lambda$ to be at least as large as the fastest process in the dynamics, which in this case is the growth rate of resistant viruses $\gamma \approx {(3\text{days})}^{-1}$, we choose $\lambda ={(0.5\text{days})}^{-1}$ (Methods). Therefore, the key parameters of the birth-death model, that is, $\beta ,\delta ,\text{and}{N}_k$ can be expressed in terms of the intrinsic fitness of the variants f_a and the neutral diversity $\theta$, which we will infer from data.

The neutral genetic diversity $\theta =2N_k\mu /\lambda$ (i.e. the number of segregating alleles) is an observable that relates to key population genetics parameters, that is, the per-nucleotide mutation rate $\mu$, the population carrying capacity $N_k$, and the total number of events per virus in the birth-death process $\lambda$, which determines the noise amplitude in the evolutionary dynamics (Methods). $N_k$ and $\lambda$ together determine the effective population size $N_e=N_k/\lambda$. We use synonymous changes as a proxy for diversity associated with the neutral variation in an HIV-1 population at a given time point within a patient. By developing a maximum-likelihood approach based on the multiplicities of different synonymous variants, we can accurately infer the neutral diversity of a population from the large survey of synonymous sites in the HIV-1 genome (Methods and Figure 2—figure supplement 1). Importantly, we infer the neutral diversity of transition ${\theta }_{\text{ts}}$ and transversion ${\theta }_{\text{tv}}$ mutations separately, and consistent with previous work (Feder et al., 2016; Zanini et al., 2017), find that transitions occur with a rate of about 8 times larger than transversions (Figure 2C, Figure 2—figure supplement 1B-E).

Our inference indicates that the neutral diversity grows over the course of an infection in untreated HIV-1 patients from Zanini et al., 2015 (Figure 2D). The patients enrolled in the three bNAb trials (Caskey et al., 2015; Caskey et al., 2017; Baron et al., 2018) show a broad range of neutral diversity prior to bNAb therapy (Figure 2D). In addition to the circulating viruses in a patient’s sera, the viral reservoir, which consists of replication-competent HIV-1 in latently infected cells or un-sampled tissue, can also contribute to a bNAb escape in a patient. Evidence that the latent reservoir can contribute to HIV escape from bnAbs is directly visible in trials as the failure of pre-trial sequencing to exclude patients who do not harbor escape variants. We model the effect of the reservoir as augmenting the neutral diversity by a constant multiplicative factor $r_{\text{resv.}}$, so that patients with more diverse sera, representing usually longer infections, are also expected to have correspondingly more diverse reservoir populations. By fitting the observed rebound data, we infer the reservoir factor $r_{\text{resv.}}\simeq 2.07$ (Methods, Figure 4—figure supplement 1). We use the augmented genetic diversity of HIV-1 prior to the bNAb therapy in each trial to generate the rebound time and the probability of HIV-1 escape in patients.

We define the mutational target size for escape from a bNAb as the number of trajectories that connect the susceptible codon to codons associated with escape variants, weighed by their probability of occurrences (Methods). The connecting paths with only single nucleotide transitions or transversions dominate the escape and can be represented as connected graphs shown in Figure 1D. To characterize the target size of escape for each bNAb, we determine the forward mutation rate $\mu \equiv {\mu }_{\text{sus.}\to \text{res.}}$ from the susceptible codons to the resistant (escape) codons, and the reverse mutation rate ${\mu }^{†}\equiv {\mu }_{\text{sus.}\leftarrow \text{res.}}$ back to the susceptible variant (Figure 1D, Methods). The mutational target sizes vary across bNAbs, with HIV-1 escape being most restricted from 10E8 ($\mu /{\mu }_{\text{ts}}=1.8$ where ${\mu }_{\text{ts}}$ is the single-nucleotide transition substitution rate) and most accessible in the presence of 10–1074 ($\mu /{\mu }_{\text{ts}}=4.9$); see Figure 4C and Appendix 1—table 1 for the list of mutational target size for escape against all bNAbs in this study.

Figure 4 with 4 supplements see all

Download asset Open asset

Since bNAbs target highly conserved regions of the virus, we expect HIV-1 escape mutations to be intrinsically deleterious for the virus (Ferguson et al., 2013; Meijers et al., 2021), and incur a fitness cost relative to pre-treatment baseline f₀. We assume that fitness cost associated with escape mutations are additive and background-independent so the fitness of a variant $a$ in the absence of bNAb follows, $f_a=f_0-{\sum }_i{\mathrm{\Delta }}_i{\sigma }_i^a$, where ${\mathrm{\Delta }}_i$ is the cost associated with the presence of a escape mutation at site $i$ (i.e., for ${\sigma }_i^a=1$); see Figure 1D.

Interestingly, we observe the escape variants against different bNAbs to be circulating in the HIV-1 populations from the cohort of ART- and bNAb-naive patients (Zanini et al., 2015, Figure 2B). We use this data (Zanini et al., 2015) and extract the multiplicity of susceptible and escape variants in HIV-1 populations at each sampled time point from a given patient. We use a single locus approximation under strong selection to represent the stationary distribution of the underlying frequency of escape alleles $x$ in each patient from (Zanini et al., 2015), $P(x;\sigma ,\theta ,{\theta }^{†})\sim x^{-1+\theta }{(1-x)}^{-1+{\theta }^{†}}\exp [-\sigma x]$, given the (scaled) fitness difference between the susceptible and the escape variants $\sigma =2N_e(f_{\text{sus}}-f_{\text{mut}})$; see Methods.

Based on the statistics of escape and susceptible variants in all patients, we define a likelihood function that determines a Bayesian posterior for selection $\sigma$ associated with escape at each site (Methods). We found that it is statistically more robust to infer the strength of selection relative to a reference diversity measure $\sigma /{\theta }_{\text{ts}}=(f_{\text{sus.}}-f_{\text{res.}})/{\mu }_{\text{ts}}$, for which we choose the transition rate (Methods). This approach generates unbiased selection estimates in simulations and is robust to effects of linkage and recombination (Methods and Model robustness Figure 4—figure supplement 1, Figure 4—figure supplement 2). The inferred values of the scaled fitness costs $\sigma /{\theta }_{\text{ts}}$ are shown for the escape-mediating sites of the trial bNAbs in Figure 2E, and are reported in Appendix 1—table 1.

In this study, we considered four clinical trials for passive therapy with bNAbs:

• Monoclonal therapy with 3BNC117 bNAb Caskey et al., 2015 and sequences reported in Schoofs et al., 2016. There were 16 patients enrolled, 13 of whom were off anti-retroviral therapy (ART).

• Monoclonal therapy with 10–1074 bNAb Caskey et al., 2017, with 19 patients enrolled, 16 of whom were off ART.

• Combination therapy with 10–1074+3BNC117 bNAb Baron et al., 2018, with 7 patients off ART.

• Monoclonal therapy with PGT121 with 17 patients off ART. Stephenson et al., 2021. Far fewer viral sequences were available for this trial, and therefore, we could not reliably infer the neural diversity of viruses within patients, as we did for other trials. We only used this trial to compare our predictions based on the DMS-inferred with the trial-inferred escape pathways in Appendix 4.

All sequenced patients across all trials were infected with distinct HIV-1 clade B viral strains. We limited our analyses to those patients not on ART at the time of treatment initiation. In these studies, the injected bNAb level falls off over time within patients and therefore, we only considered dynamics within an 8-week window since infusion. This assures that rebound is not confounded by a drop in bNAb below sensitive-strain neutralizing levels of ${\displaystyle {\text{IC}}_{{50}_S}<2\mu \text{g/m}}$.

We used single-genome sequence data of env collected from all patients in each trial to characterize the diversity of HIV-1 population within each patient shown in Figure 2 available from Zanini et al., 2015 and accessible through European Nucleotide Archive (Accession no: PRJEB9618). The patient sequence data for each trial is available through Caskey et al., 2015 GeneBank accession number: KX016803, Caskey et al., 2017 GenBank accession numbers KY323724.1 - KY324834.1, and Baron et al., 2018, GenBank accession numbers MH632763 - MH633255.

Single nucleotide polymorphism (SNP) data was obtained from Zanini et al., 2015 and aligned to the HXB2 reference using HIV-align tool (Gaschen et al., 2001). The dataset includes 11 patients observed for 5–8 years of infection, with HIV-1 sequence data sampled over 6–12 time points per patient (Figure 2).

Patient 4 and patient 7 were excluded from the original analysis done in Zanini et al., 2015 because of suspected superinfection and failure to amplify early samples, respectively. These patients were included in our analysis, since (i) super-infection poses no additional difficulties for our tree-free procedure, and (ii) only time points with measurable viral diversity entered into our selection likelihood, which automatically limits our analysis to samples with high-quality sequences.

All patients were infected with clade B of HIV, except for patient 6 (clade C), and patient 1 (clade 01_AE). We assessed the robustness of our inference to exclusion of these patients from our analysis in Effect of genomic linkage on the inference of selection. Overall, our inference was not strongly affected by this choice (Model robustness Figure 4—figure supplement 1), and therefore, we included these patients in our main analyses to enhance the statistics with larger data.

For our analysis we considered only data reported in single nucleotide polymorphism (SNP) counts. The number of raw SNP counts are the result of amplification and must be converted into estimates for the number of pre-amplification template fragments. Zanini et. al. reported that on average about 10² templates of amplicon were associated with the fragments of the envelope (env) protein Zanini et al., 2015. We converted raw SNP counts into templates by normalizing to 120 counts and rounding the resulting number to the nearest integer.

The concentration of viral RNA copies in blood serum is a delayed reflection of the total viral population size $N(t)$, containing a resistant and susceptible subpopulations, with respective sizes $N_r(t)$ and $N_s(t)$. After infusion of bNAbs in a patient, the susceptible sub-population decays due to neutralization by bNAbs and the resistant sub-population grows and approaches the carrying capacity $N_k$, with the dynamics,

Here, $\gamma$ is the growth rate of the resistant population, and $r$ is the neutralization rate impacting the susceptible subpopulation. By setting the initial condition for fraction of resistant subpopulation prior to treatment (at time $t=0$) $x=N_r(0)/(N_r(0)+N_s(0))$, we can characterize the evolution of the total viremia in a patient. This dynamics is governed by the combined processes of neutralization by the infused bNAb and the viral rebound (Figure 1A), which entails,

We use Equation 5 to define the evolution of blood concentration of viral RNA sequences which is observed indirectly via noisy viremia measurement data from Caskey et al., 2015; Caskey et al., 2017; Baron et al., 2018. To connect the data with the simple model of viral dynamics in Equation 5, we fit the initial frequency of resistant mutants $x$ for each patient separately, and fit a global estimate for the decay rate of susceptible variants $r$ shared across all patients in a trial, using a joint maximum-likelihood procedure. In addition, we fix the growth rate $\gamma$ to 1/3 days, corresponding to a doubling time of approximately 2 days García et al., 1999. Our analyses indicate that the initial viremia decline lagged treatment by about 1 day (Figure 3—figure supplements 1–3), consistent with previous findings (Ioannidis et al., 2000), and therefore, we included a 1 day lag between the fitted viremia response model and the treatment.

$t$ The number of viral RNA copies in a blood sample is subject to count fluctuations with respect to the true number of circulating virions in a given volume of the blood. We use a Poisson sampling model to define a likelihood for our model of viral population. The likelihood of observing viral counts in a sample collected from patient at time is given by a Poisson distribution,

with rate parameter set by the model value of the viral multiplicity $N_p(t)$ (Equation 5). We use the Poisson likelihood in Equation 6 to characterize an error model to fit the parameters of the viral dynamics in Equation 5. However, since the mean and variance of the Poisson distribution are related, combining data with different mean values $N_p(t)$ at different times and from different patients can cause inconsistencies in evaluations of errors in our fits. To overcome this problem, we use a variance stabilizing transformation (McCullagh and Nelder, 2019) and define a change in variable ${\widehat{n}}_p(t)=\sqrt{n_p(t)}$. This transformed variable has a constant variance, and in the limit of large-sample size, it is Gaussian distributed with a mean and variance given by, ${\displaystyle {\hat{n}}_p(t)\sim \mathcal{N}(\sqrt{\lambda },1/4)}$. The constant variance of the transformed variable enables us to combine data from all patients and time points, irrespective of the sample’s viral loads, and fit the model parameters $(r,x_p)$ using (non-linear) least-squares fitting of the function

Here, $N_p(t|r,x_p)$ is the model estimate of viremia in patient $p$ at time $t$ (Equation 5), given the pre-treatment fraction of resistant variants x_p, and the decay rate $r$.

Note that the viremia measurements have a minimum sensitivity threshold of 20 RNA copies per ml. We treat the data points below the threshold of detection as missing data and if $n_p(t)$ is below the threshold of detection we impute $n_p(t)=\text{min}(20,N_t)$.

The fitted viremia curves for patients enrolled in the three bNAb trials under consideration are shown in Figure 3—figure supplements 1–3, and the respective decay rates $r$ for each experiment are,

To encode for different viral variants, we specify a coarse-grained phenotypic model, where a viral strain of type $a$ is defined by a binary state vector ${\overrightarrow{\rho }}^a=[{\rho }_1^a,\mathrm{\dots },{\rho }_{\mathrm{\ell }}^a]$, with $\mathrm{\ell }$ entries for potentially escape-mediating epitope sites; the binary entry of the state vector at the epitope site $i$ represents the presence (${\rho }_i^a=1$) or absence (${\rho }_i^a=0$) of a escape mediating mutation against a specified bNAb at site $i$ of variant $a$. We assume that a variant is resistant to a given antibody if at least one of the entries of its corresponding state vector is non-zero.

We define an individual-based stochastic birth-death model (Wilkinson, 2019) to capture the competitive dynamics of different HIV variants within a population. This dynamic model will allow us to predict the distribution of rebound times under any combination of antibodies.

We assume that a viral strain of type $a$ can undergo one of three processes: birth, death and mutation to another type $b$ with rates ${\beta }_a$, ${\delta }_a$, and ${\mu }_{a\to b}$, respectively:

We specify an intrinsic fitness f_a for a given variant a, defined as the growth rate of the virus in the absence of neutralizing antibody or competition. Since bNAbs target highly vulnerable regions of the virus, we expect that HIV-1 escape mutations to be in trinsically deleterious for the virus and to confer a fitness cost relative to the susceptible viral variants prior to the infusion of bNAbs. Assuming that fitness cost of escape is additive across sites and background-independent, we can express the fitness of a variant as , where Δ is the cost associated with the presence of an escape mutation at site of variant (i.e., for = 1).

We assume that growth is self-limiting via a competition for host T-cells. This competition enforces a carrying capacity, which sets the steady-state population size $N_k$. Competition is mediated through a competitive pressure term $\varphi =\frac{{\sum }_a{N}_a{f}_a}{N_K}$ which attenuates the net growth rate ${\gamma }_a$ so that ${\gamma }_a=f_a-\varphi$. At the carrying capacity, the competitive pressure equals the mean population fitness $\varphi =\overline{f}$, making the net growth rate of the population zero.

The net growth rate of a variant $a$ is given by its birth rate minus the death rate: ${\gamma }_a={\beta }_a-{\delta }_a$. We assume that the total rate of events (i.e. the sum of birth and death events) is equal for all types, that is, $\lambda ={\beta }_a+{\delta }_a,\forall a$. Assuming that $\lambda$ is constant is to be agnostic about the mechanism of a fitness decrease, attributing fitness loss equally to (i) an increase in the death rate and (ii) a decrease in the birth rate.

Because the absolute magnitude of $\beta$ and $\delta$ asymptotically converge in the continuum limit for a surviving population, that is, ${lim}_{N\to \mathrm{\infty }}\beta /\delta =1$, it is impossible to distinguish between (i) and (ii) in the continuous limit. Choosing constant $\lambda$ simplifies both theoretical calculations and the simulation algorithm.

This leads to the following equations for the birth and the death rates:

In the presence of an antibody, birth is effectively halted for susceptible variants, resulting in birth and death rate for a susceptible variant $s$,

so that the susceptible phenotype decays at rate $r$.

We assume that mutations occur independently at each site,

where $1_s$ is the vector which has only one non-zero entry at site $i$, and ${\mu }_i$ and ${\mu }_i^{†}$ are the forward the backward mutation rates at site $i$, respectively. We characterize the state of a population by vector ${\displaystyle \mathit{n}=(n_1,\dots n_M)}$, where n_a is the number of type $a$ variants within the population. We approximate the evolution of the population state distribution using a Fokker-Planck equation (see Appendix 2).

Of special interest for our analysis is the steady-state density of frequencies, which we use to describe the initial state of the population (before treatment) and to infer selection intensity. In the steady state, the population is fluctuating around carrying capacity ${\sum }_a{n}_a\approx N_k$ and we can represent the population state via allele frequencies $x_a=n_a/N_k$. In the simple case of a bi-allelic problem, the equilibrium allele frequency distribution $P_{\text{eq}}(x)$ follows the Wright-equilibrium distribution (Crow and Kimura, 2010) with modified rates,

where $Z$ is the normalization factor, f₁ is the intrinsic fitness of the variant of interest, f₀ is the fitness of the competing variant, $\mu$ and ${\mu }^{†}$ are the forward and backward mutation rates, $N_k$ is the carrying capacity, and $\lambda$ is the total rate of events in the birth-death process, which sets a characteristic time scale over which the impact of selection and mutations can be measured. In this case, we can define an ‘effective population size’ that sets the effective size of a bottleneck and the natural time scale of evolution as $N_e=N_k/\lambda$, and specify a scaled selection factor $\sigma =N_{e}s=N_e(f_0-f_1)$, and scaled forward mutation and backward mutation rates (diversity) $\theta =2N_e\mu$, and ${\theta }^{†}=2N_e{\mu }^{†}$. The normalization factor is given by,

where ${}_{1}F_1(\cdot )$ denotes a Kummer confluent hypergeometric function and $ℬ(\theta ,{\theta }^{†})=\frac{\mathrm{\Gamma }[\theta ]\mathrm{\Gamma }[{\theta }^{†}]}{\mathrm{\Gamma }[\theta +{\theta }^{†}]}$ is the Euler beta function.

The logistic dynamics describing a patient’s viremia over time in Equation 5 is the deterministic approximation to the underlying birth-death process. However, the resistant population can also go extinct due stochastic effects, which in turn contribute to the probability of late rebound in a population. To capture this effect, we derive an approximate closed form expression for the probability of extinction.

Using the standard birth-death process generating function theory Allen, 2010 the probability $P(\text{extinct}|n_i)$ that a population consisting of n_i resistant variants of type $i$ go extinct can be expressed as,

To characterize the probability of extinction for a population of size $N_k$ with pre-treatment fraction of $i^{th}$ resistant variants x_i, we can convolve the extinction probability in Equation 14 with a Binomial probability density for sampling n_i resistant variants from $N_k$ trials. Given that the pre-treatment fraction of resistant variants is small $x_i\ll 1$ and $N_k$ is large, this Binomial distribution can be well approximated by a Poisson distribution, $\text{Poiss}(n_i;N_k{x}_i)$ with rate $N_k{x}_i$, resulting in an extinction probability,

Using the expressions for the growth in the absence of competition, ${\beta }_i-{\delta }_i={\gamma }_i=f_i$ (since $\varphi =0$), and assuming that fitness is small relative to the total rate of birth and death events $f_i\ll \lambda$, we can use the approximation ${\beta }_i=(\lambda +f_i)/2\approx \lambda /2$, to arrive at,

$x_{\text{ext}}$ where the characteristic escape threshold can be written in terms of concrete genetic observables,

Figure 3 shows that this threshold can well separate the fate of stochastic evolutionary trajectories, simulated with relevant parameters for intra-patient HIV evolution.

To treat the full viral dynamics including mutations, and transient competition effects, we can exactly simulate the viral dynamics defined by our individual based model. Below are the key steps in this simulation.

The nucleotide triplets which encode for amino acids at an escape site undergo substitutions which can change the amino acid type and create an escape variant. The changes in the state of an amino acid codon can be modeled as a Markov jump process and can be visualized as a weighted graph where the nodes represent codon states, and edges represent single nucleotide substitutions linking two codon states (Figure 1D). In our mutational model, these edges have weights associated with either the mutation rates for transitions ${\mu }_{\text{ts}}$ or transversions ${\mu }_{\text{tv}}$. We call this the codon substitution graph.

The codon states can be clustered into three distinct classes: (i) codons which are fatal $F$, (ii) wild-type (i.e. susceptible to neutralization by the bNAb) $W$, and escape mutants (i.e. resistant to the bNAb) $M$. We expect the escape mutants to be at a selective disadvantage compared to the resistant wild-type, and that the most common escape codons to be those which are adjacent to wild-type states.

The mutational target size is determined by the density of paths from the wild-type $W$ to the escape mutants $M$.

The functions $[c-d=\text{ts}]$ or $[c-d=\text{tv}]$ are 1 when the two codons are separated by a transition or transversion, and are zero otherwise. Note that since we only have an estimate for the ratio of the transition to transversion rates ${\mu }_{\text{tv}}/{\mu }_{\text{ts}}$, we can only determine the scaled mutational target sizes, $\widehat{\mu }=\mu /{\mu }_{\text{ts}}$ and ${\widehat{\mu }}^{†}={\mu }^{†}/{\mu }_{\text{ts}}$, which are sufficient for inference of selection in the next section. The full list of mutational target sizes inferred for the bNAbs in this study are shown in Appendix 1—table 1.

When discussing the mutational target size of escape from a given bNAb, we refer to the total mutation rate from the susceptible (wild type) to the escape variant as,

where the sum runs over all the mediating escape sites $i$, and can be interpreted as the average number of accessible escape variants. In the strong selection regime, we can write the frequency of escape mutants as,

Here, we develop an approximate likelihood approach to infer the selection ratio $\widehat{\sigma }=\frac{\sigma }{{\theta }_{\text{ts}}}$, using the high-throughput sequence data from bNAb-naive HIV-1 patients from Zanini et al., 2017. The quantity $\widehat{\sigma }$ is a dimensionless ratio which is independent of the coalescence timescale $N_e$, and therefore, represents a stable target for inference.

We assume that the probability to sample $m$ escape mutants and $s$ susceptible (wild-type) alleles at a given site in the genome of HIV in a population sampled from a patient at a given time point follows a binomial distribution $\text{Binom}(m,s|x)$, governed by the underlying frequency $x$ of the mutant allele. In addition, the frequency $x$ of the allele of interest itself is drawn from the equilibrium distribution $P_{\text{eq}}(x|\sigma ,\theta ,{\theta }^{†})$ (Equation 12), governed by the diversity ${\theta }_{\text{ts}}$ inferred from the neutral sites, the estimated mutational target sizes $\widehat{\mu }=\mu /{\mu }_{\text{ts}},{\widehat{\mu }}^{†}={\mu }^{†}/{\mu }_{\text{ts}}$, and the unknown selection ratio $\widehat{\sigma }=\sigma /{\theta }_{\text{ts}}$. As a result, we can characterize the probability $P(m,s|{\theta }_{\text{ts}},\widehat{\mu },{\widehat{\mu }}^{†},\widehat{\sigma })$ to sample $m$ escape mutants and $s$ susceptible-type alleles, given the scaled selection and diversity parameters as,

Here, $Z(\sigma ,\theta ,{\theta }^{†})$ is a confluent hypergeometric function of the model parameters that sets the normalization factor for the allele frequency distribution $P_{\text{eq}}(x)$ (Equation 13). It should be noted that the viral population is in fact out of equilibrium, due to constant changes in immune pressure evolution from the B-cell and T-cell populations (Nourmohammad et al., 2016). Although we are ignoring these significant complications, we later use the same equilibrium distribution in a consistent way to generate standing variation in simulations. For the model to make accurate predictions, it is not necessary that the equilibrium model be exactly correct, but only that it is rich enough to provide a consistent description for the distribution of mutant frequencies observed across viral populations.

We will use the probability density in Equation 29 to define a log-likelihood function in order to infer the scaled selection $\widehat{\sigma }=\sigma /{\theta }_{\text{ts}}$ from data. To do so, we first express the logarithm of this probability density as,

where the constant factors (const.) are independent of selection, and ${\displaystyle \mathbb{E}[\cdot {]}_{\text{Beta}(\cdot )}}$ denotes the expectation of the argument over a Beta distribution with parameters specified in the subscript. The expression in Equation 30 implies that we can evaluate the likelihood of selection strength by computing the difference between the logarithms of the expectation for $e^{-\sigma x}$ over allele frequencies drawn from two neutral distributions (Beta distributions), with parameters $(\theta ,{\theta }^{†})$ and $(\theta +m,{\theta }^{†}+s)$, respectively. This approach is more attractive as it would not require direct evaluation of the confluent hypergeometric functions for the normalization factors in Equation 29. Estimating these normalization factors is computationally intensive for large values of $\sigma$, since many terms in the underlying hypergeometric series should be taken into account to stably compute them. However, evaluating the expectations via sampling from these two neutral distributions has the disadvantage that it is subject to variations across simulations. We reduce the variance of our estimate of $\log P(m,s|\sigma ,\theta ,{\theta }^{†})$ in Equation 30 by using a mixture-importance sampling scheme (Owen and Zhou, 2000) with the details shown in Algorithm 4 (Appendix 5).

We use data collected across all time points and from all patients to infer reliable estimates for selection strengths. However, allele frequencies are correlated across time points within patients (Figure 2B), and thus, sequential measurements are not independent data points. Our estimates indicate a coalescence time of about $N_e\sim {10}^3$ days based on the estimates for the mutation rate ${\mu }_{\text{ts}}={10}^{-5}$ /nt/day, and the neutral diversity ${\theta }_{\text{ts}}=2N_e{\mu }_{\text{ts}}=0.01$. This coalescence time is much longer that the typical separation between sampled time points within a patient ($\sim {10}^2$ days), suggesting that sequential samples collected from each individual in this data are correlated. Therefore, we treat each patient as effectively a single observation, using the time-averaged likelihood for the (scaled) selection factor $\widehat{\sigma }$:

where $p$ and $t$ denote patient identity and sampled time points, respectively, and $T_p$ is the total number of time points sampled in patient $p$. We use the likelihood in Equation 31 to generate samples from the posterior distribution for selection strengths under a flat prior, with a standard Metropolis-Hastings algorithm Hastings, 1970. Since the prior is constant, this procedure amounts to simply accepting or rejecting samples based on the likelihood ratio of Equation 31. We used the centered-normal distribution with standard deviation of 50 ($\times {\mu }_{\text{ts}}$ in absolute units) as the proposal density for the jumps in the Markov chain.

Prior work has also inferred the fitness effect of mutations in HIV, but our approach differs in important aspects. For instance, maximum entropy models have been used to infer the preference of different amino acids in the Gag and the env proteins of HIV-1 from their prevalences across sequences sampled from different patients( Louie et al., 2018; Ferguson et al., 2013). The inferred preference values can explain the in-vitro growth rate (fitness) of the associated viral strains, especially for sites that are relatively conserved and are not the drivers of antigenic evolution in HIV-1. However, further modeling would be needed to quantitively map these inferred amino acid preferences onto population genetics measures that can be used to characterize the evolutionary dynamics of an HIV population. Other work has used longitudinal HIV-1 sequence data to infer selection and to characterize the role of selective sweeps due to immune pressure, genetic hitchhiking and recombination in the turnover of HIV-1 populations within patients (Zanini et al., 2017; Neher and Leitner, 2010; Illingworth et al., 2020; Haddox et al., 2018). In contrast, our work focuses on the expected composition of the population in a viremic patient prior to the start of treatment as opposed to the history of a viral population. Our approach uses a self-consistent formalism for inference of the population genetics parameters (e.g. population diversity and selection strength) and for the evolutionary simulations used to predict outcomes. Therefore, we can directly interpret the fitted parameters in terms of both the viral dynamics and the pre-treatment state of HIV-1 populations within patients.

We expect different distributions of patient outcomes depending on whether they have been recently infected and thus have relatively low viral diversity, or whether their infection is longstanding with a diverse viral population. To construct the distribution of initial population diversities ${\theta }_{\text{ts}}$ for simulating trial outcomes, we apply the ${\theta }_{\text{ts}}$ inference procedure (Equation 22) to pre-treatment sequence datasets available from the clinical trials under consideration. In Figure 2D this set of pre-trial ${\theta }_{\text{ts}}$ is compared to the longitudinal in-patient ${\theta }_{\text{ts}}$ from Zanini et al., 2017. We used random draws from the inferred ${\theta }_{\text{ts}}$ values for patients to generate ${\theta }_{\text{ts}}$ for simulations.

We found that there was considerably more viral escape and non-responders in our simulations than in the observed data as shown in Figure 4—figure supplement 3A. This is in addition to the fact that the patients were screened to have only susceptible variants to the antibodies used in trials (Caskey et al., 2015; Caskey et al., 2017; Baron et al., 2018). In theory, there should be zero non-responders, as such patients should have been excluded by screening. The over-prediction of both non-responders and late rebounds is a signature of undercounting the effective diversity of the viral populations.

The failure of both screening and our naive prediction in undercounting the diversity in the viral population can be explained by an effective viral reservoir. Viral variants which mediate rebound can come from compartments such as including bone marrow, lymph nodes, and organ tissues, and can be genetically distinct from those sample from the plasma T-cells during screening (Chaillon et al., 2020; Wong and Yukl, 2016; Chun et al., 2005). This reservoir of viral diversity can reappear in plasma after infusion of a bNAb and could in part contribute to treatment failure (Avettand-Fenoel et al., 2007; Shan and Siliciano, 2013; Sharkey et al., 2011; Tagarro et al., 2018).

We model the effect of reservoirs as a simple inflation of the diversity observed by a multiplicative factor $\xi$. We fit $\xi$ directly to trial observations, using a disparity-based approach by minimizing an empirical divergence estimator Ekström, 2008 between the observed and simulated data. To do so, we characterize the Hellinger distance Lindsay, 1994; Simpson, 1987 between the true distribution of rebound times $P(t)$ and the rebound times $Q(t|\xi )$ generated by simulations with a given reservoir factor $\xi$,

Algorithm 5 (Appendix 5) defines the procedure that we use to estimate the Hellinger distance D_H(Q(t)||P(t|ξ)). Specifically, we use n_q quantiles of the observed data x_(i) ~ Q to partition the space of observations into discrete outcomes

where $P_{(i)}$ is a constant by construction, and $P_{(i)}(\xi )$ is estimated by simulations, and $[\cdot ]$ is the Iverson bracket Graham et al., 1989; see Algorithm 5 (Appendix 5).

To simulate data for this analysis, we generate $S$ rebound times ($T_{1:S}$) by simulations, given the scaled diversity values $\xi {\theta }_{\text{ts}}$. We then find the optimal value ${\xi }^{*}$ by minimizing the disparity with the observed rebound times $t_{1:p}$ by brute-force search,

Here, ReboundDisparity is the function defined by Algorithm 5 (Appendix 5); see Ekström, 2008 for details. We find the optimal reservoir factor to be ${\xi }^{*}=2.1$, which we use in subsequent therapy prediction. The disparity over various values of $\xi$ for different trials is shown in Figure 4—figure supplement 3.

Given the reservoir-corrected estimate of the diversity $\xi \theta$ and the posterior samples for selection factors $\widehat{\sigma }$, we now summarize how we simulated the outcome of clinical trials.

For a given bNAb, we draw the selection factor at each of the escape mediating sites from the corresponding Bayesian posterior on $\widehat{\sigma }$; the posterior distributions are shown in Figure 4A, and summarized in Appendix 1—table 1. We also use the mutational target size (forward $\mu$ and backward ${\mu }^{†}$ rates) associated with each of the escape mediating sites of a given bNAb; see Appendix 1—table 1. The result can be summarized in a mutation / selection matrix ${\widehat{M}}_a$ for a given bNAb $a$, where each column corresponds to an escape mediating site $i$ against the bNAb,

The elements of the matrix ${\widehat{M}}_a$ are the scaled mutation and selection factors, i.e., ${\widehat{\mu }}_i={\mu }_i/{\mu }_{\text{ts}},{\widehat{\mu }}_i^{†}={\mu }_i/{\mu }_{\text{ts}},\text{and}{\widehat{\sigma }}_i={\sigma }_i/{\theta }_{\text{ts}}$, where the absolute value of mutation rate is set to ${\mu }_{\text{ts}}=1.11\times {10}^{-5}/\text{nt}/\text{day}$ from Zanini et al., 2017.

For each patient in our simulated trial, we then draw diversity ${\theta }_{\text{ts}}$ from the patient pool, and scale it by our fitted $\xi =2.1$, resulting in patient-specific selection and mutation factors,

These parameters are then used to initialize the state of an HIV population within a patient according to (Appendix 5), and to determine the absolute rates in Algorithm 3 (Appendix 5) for the population evolution according to Equation 38. The decay rate is set to the fitted trial average of $r=0.31{\text{days}}^{-1}$ (Equation 8). The carrying capacity $N_k$ is set according to Equation 21. This determines all parameters of the birth-death process simulating the intra-patient evolution of HIV, which are used in Algorithm 3 (Appendix 5).

We evolve a population through time until 56 days have elapsed since treatment, or until the escape fraction relative to the carrying capacity x_t is above 0.8; see Algorithm 3 (Appendix 5). After ${\displaystyle x_t>.8}$ the evolution is governed by the deterministic equations, and the stochastic simulation ends. The rebound time $T$, defined as the intersection of the exponential envelope and the carrying capacity, can then be calculated analytically as,

The resulting distribution for rebound times are shown as model predictions in Figure 3B–D.

Rebound times generated in this fashion were also used to estimate the probability of late rebound to characterize the efficacy of a given bNAb in curbing viral rebound. The probability of late rebound was estimated from 10⁴ simulated patients. The interdecile quantiles (0.1–0.9) of early rebound (lt₅₆ days) probability over 200 values of scaled selection coefficients $\widehat{\sigma }$ drawn from the posteriors in Figure 4A are shown in Figure 4C.

In our inference of selection (Equation 31, Figure 4A), we assume that the escape-mediating sites are at linkage equilibrium and that the distribution of allele frequencies can be approximated by a skewed Beta distribution (Equation 12), reflecting the equilibrium of allele frequencies. In reality, despite recombination, the HIV genome exhibits linkage effects, especially at nearby sites Zanini et al., 2015, and the viral populations experience changing selective pressures by the immune system Feder et al., 2021; Theys et al., 2018; Nourmohammad et al., 2019, and the transient population bottlenecks during therapy Feder et al., 2016.

To test the limits on the validity of our inference procedure, we applied it to in silico populations generated by full-genome forward-time simulations (Appendix 5, Algorithm 3) in the presence and absence of recombination. To do so, we considered an ensemble of ten patients with 100 genomes sampled at 10 time points, and used two diversity parameters ${\theta }_{\text{ts}}=0.01$ and ${\theta }_{\text{ts}}=0.1$, to cover the range reflected in patient data (Figure 2D, Figure 4—figure supplement 2).

One relevant scenario to consider is the impact of other selected sites in the genome on the distribution of alleles at the escape mediating sites against bNAbs. The sites under a strong constant selection are likely to be already fixed (or at a high a frequency) at their favorable state in the population. However, the strong selection on a large fraction of antigenic sites can be thought as time-varying, due to the changing pressure imposed by the immune system or therapy. To capture this effect, we simulated whole genome evolution in which linked sites were under strong selection ($0.1\times \text{growth rate}$), and where the sign of selection changed after exponentially distributed waiting times (i.e. as a Poisson process); this model of fluctuating selection has been used in the context of influenza evolution Strelkowa and Lässig, 2012, and for somatic evolution of B-cell repertoires in HIV patients Nourmohammad et al., 2019. The resulted evolutionary dynamics in this case can involve strong selective sweeps and clonal interference due to the continuous rise of beneficial mutations (in the linked sites) within a population.

To test the robustness of our selection inference, we evaluated the distribution of maximum likelihood estimates (MLEs) for the selection values $\widehat{\sigma }=\sigma /{\theta }_{t}s$ at the escape mediating sites, inferred from the ensemble of sequences obtain from simulated data with linkage. Figure 4—figure supplement 2 shows that even for fully linked genomes (zero recombination) our MLE estimate of selection has little bias relative to the true values used in the simulations. Adding recombination into the simulations only further attenuates the effect of linkage (Figure 4—figure supplement 2), making the estimates more accurate.

The reason that selective sweeps of linked beneficial mutations have only minor effects on our inference of selection for the escape mediating sites is two-fold: First, the primary mechanism by which selective sweeps change the strength of selection at linked sites is via a reduction in the effective population size Hill and Robertson, 2009; Comeron et al., 2008. However, variations in the effective population size are already accounted for in our inference procedure: The selection likelihood in Equation 31 is conditioned on the measured neutral-site diversity, $\theta =2N_e\mu$. The change in the effective population size impacts the selection coefficient $\sigma =2N_e\mathrm{\Delta }$ and the diversity $\theta =2N_e\mu$ in the same way, and therefore, the (scaled) selection parameter $\widehat{\sigma }=\sigma /\theta$ that we infer from data remains insensitive to changes in the effective population size.

The second reason for the robustness of our selection inference to linkage is due to the fact that a beneficial allele in a linked locus can appear on a genetic background with or without a susceptible variant, leading to the rise of either variants in the population. As a results, the impact of such hitchhiking remains a secondary issue in inference of selection at the escape sites, for which an ensemble of populations from different patients with distinct evolutionary histories of HIV-1are used.

To infer the selection effect of mutations from the longitudinal deep sequencing data of Zanini et al., 2015, we use time averaging of the likelihood (Equation 31) to avoid conflating our results due to temporal correlations between the circulating alleles within patients (Figure 2B). We can view this choice as being one choice among two extremes: (i) to treat each patient as effectively one independent data point so that all patients are given the same weight or (ii) to treat each time point as independent, giving patients with more time points a higher weight. These two choices correspond to different log-likelihood functions for the (scaled) selection factor $\widehat{\sigma }$:

where $T_p$ is the total number of time points from patient $p$, and $P_p(m_t,s_t|{\theta }_{\text{ts}}(t),\widehat{\sigma })$ is the probability to observe $m$ escape mutants, and $s$ susceptible variants at time $t$ in patient $p$, given the neutral diversity ${\theta }_{\text{ts}}(t)$ and the scaled selection factor $\widehat{\sigma }$.

We find that both of these approaches result in similar posteriors for selection $\widehat{\sigma }$ (Model robustness Figure 4—figure supplement 1A and C), although the $t$-averaged likelihood has a higher uncertainty due to fewer independent time points. Thus, our inference of selection is insensitive to the exact choice of the likelihood function given in Equation 40, yet our time-averaged approach remains the more conservative choice between the two.

In Equation 36, we introduced the reservoir factor ${\xi }^{*}=2.1$ to account for the diversity of HIV-1that is not sampled from a patient’s plasma prior to therapy, which resulted in a better fit of the rebound time distributions (Figure 3) compared to a reservoir-free model (Figure 4—figure supplement 3A). Here, we quantify the importance of the reservoir factor with a statistical test on the null hypothesis, ${\xi }_0=1$. Specifically, we perform a hypothesis test to test the necessity of using an inflated diversity $\xi {\theta }_{\text{ts}}$ relative to using the bare diversity observed in pre-trial sequence data ${\theta }_{\text{ts}}$. To do so, we construct a disparity-based test statistic Ekström, 2008, which is analogous to the likelihood ratio test statistic.

Recall that the optimal reservoir factor ${\xi }^{*}=\arg {\min }_{\xi }R(\xi |t_{1:p})$ was obtained by minimizing the disparity function $R(\xi |t_{1:p})$ across measurements of rebound times $t_{1:p}$ from all the $p$ patients in data (Equation 36). We can estimate the test statistic for the reduction in disparity between the null hypothesis, ${\xi }_0=1$ and the fitted reservoir factor ${\xi }^{*}$ as,

We can then determine the p-value by estimating the quantile of the observed test statistic ${\displaystyle {\mathrm{\Delta }}_R(t_{1:p})}$ relative to that inferred from the distribution of ${\displaystyle {\mathrm{\Delta }}_R(T_{1:p})}$ obtained from simulations under null hypothesis ${\displaystyle {\xi }_0=1}$ (Fisher, 1956). Specifically,

where ${\displaystyle \underset{T_{1:p}|{\xi }_0}{\mathbb{E}}(\cdot )}$ denotes the expectation over the rebound times $T_{1:p}$ obtained from 1000 realizations of simulated populations each with $p$ patients, and under the null hypothesis ${\displaystyle {\xi }_0=1}$ is the Iverson bracket that takes value 1 when its argument is true and 0, otherwise (Equation 34). The observed ${\mathrm{\Delta }}_R$ (Equation 41), the distribution of simulated values of ${\mathrm{\Delta }}_R(T_{1:p}|{\xi }_0)$ under the null hypothesis, and the resulting $\text{p-value}=0.004$ are shown in Figure 4—figure supplement 3B, C.

It should be noted that here we use the disparity measure because the corresponding likelihood function for the reservoir factor is inaccessible through forward simulations of populations. However, a general analogy exists between our approach and the more commonly used likelihood approach. Specifically, in an analogous likelihood-ratio test, the test statistic ${\mathrm{\Delta }}_L={\max }_{\xi }\log p(\xi |t_{1:p})-\log p({\xi }_0|t_{1:p})$ would be asymptotically ${\chi }^2$-distributed with one degree of freedom under the null hypothesis Fisher, 1954, and the quantile under the null-hypothesis (p-value) would be estimated by inverting the ${\chi }^2$ cumulative distribution function (i.e. a ${\chi }^2$ test).

The longitudinal deep sequencing data of Zanini et al., 2015 is collected from 11 patients, 9 of whom are infected with clade B strains of HIV-1, which is the dominant clade circulating in Europe Spira et al., 2003. All of the clinical trials we considered Caskey et al., 2015; Caskey et al., 2017; Baron et al., 2018 are from patients carrying clade B strains. For the results presented in the main text, we included all the 11 patients in our analysis. Here, we test wether our inference of selection is sensitive to the choice of including or excluding non-clade B patients in our analysis. We therefore repeated our inference procedure for selection by excluding the two non-clade B patients. Model robustness Figure 4—figure supplement 1A,B shows a strong agreement between the Bayesian posterior for selection factors in the two cases, with a slight increase in uncertainty for the case with only clade B patients. This increased uncertainty is related to the reduction in sample size by excluding the non-clade B patients from data. Nonetheless, the richness of the intra-patient diversity makes the inference robust to the exclusion of one or two patients.

How sensitive are the outcomes of our predictions for the rebound time distributions (Figure 1) to the exact values of inferred selection strengths we used for our simulations? We addressed this question by performing a disparity analysis similar to that for the diversity $\theta$ in Equation 41. Specifically, we assessed whether we might need to rescale our inferred selection strength $\sigma /{\theta }_{\text{ts}}$ by a multiplicative factor ${\xi }_s$ (Figure 4—figure supplement 4).

In contrast to diversity, the reduction in disparity for adjustment of selection with a factor ${\xi }_s$ is small (Figure 4—figure supplement 4A) and not statistically significant ($\text{p-value}=0.49$; Figure 4—figure supplement 4B), and could be attributed to count noise. Still, we cannot discount the possibility that selection was slightly overestimated, possibly due to the effect of compensatory mutations in linked genome, the interplay between the reservoir and the inference of selection, or other biological factors. Nonetheless, in absolute terms, the null hypothesis (i.e. ${\xi }_s=1$) cannot be rejected and we have no statistical justification for adding an adjustment factor for selection inferred from untreated patients.

Our predictions rely on identifying escape variants against each bNAb, either based on in vivo trial data and patient surveillance, or in vitro DMS assays. To test the sensitivity of our results to these methods, we compare the predictions of rebound time distributions when identifying escape sites from the DMS data versus the trial data. In Appendix 4, we perform this comparison for the 10–1074 and the PGT121 bNAbs; we do not include the 3BNC117 bNAb in this analysis since it targets the CD4 binding site of HIV-1and the DMS data is unreliable for identifying escape variants against it. As shown in Appendix 4—figure 1, both trial-inferred and DMS-inferred escape sites result in good predictions for rebound time distributions. However, its appears that using the DMS-inferred escape sites could lead to a more optimistic prediction for treatment success (i.e. a later rebound). This inconsistency may be due to the fact that DMS data is collected in vitro, and other biological factors could be influencing the in vivo escape patterns. Moreover, the differences in the genetic composition of the HIV-1strains circulating in patients enrolled in trials and the strains used for the DMS experiments could lead to different epistatic interactions that can enhance or reduce the chances of escape. For a systematic understanding of these differences and limitations, more experiments would be necessary. Nonetheless, DMS data can provides baseline to gauge the efficacy of different bNAbs for therapy and further in vivo investigation.