# Predicting in vivo escape dynamics of HIV-1 from a broadly neutralizing antibody

^{a}Institut für Biologische Physik, University of Cologne, 50937 Cologne, Germany;^{b}Laboratory of Experimental Immunology, Institute of Virology, Faculty of Medicine, University of Cologne, 50931 Cologne, Germany;^{c}Laboratory of Experimental Immunology, Institute of Virology, Faculty of Medicine and University Hospital Cologne, University of Cologne, 50931 Cologne, Germany;^{d}Partner Site Bonn–Cologne, German Center for Infection Research, 50931 Cologne, Germany;^{e}Center for Molecular Medicine, University of Cologne, 50931 Cologne, Germany

See allHide authors and affiliations

Edited by Stephen P. Goff, Columbia University Irving Medical Center, New York, NY, and approved June 22, 2021 (received for review March 11, 2021)

## Significance

Broadly neutralizing antibodies (bnAbs) show promise for antibody-mediated prevention or treatment of HIV-1 infection. Recent clinical trials, however, indicate that the virus can evolve resistance mutations that escape neutralization by bnAbs. Here, we establish a fitness model for the escape dynamics of HIV in humans. We show that this model can be applied universally across different human hosts, providing a proof of principle that the in vivo response of HIV to bnAb therapies is predictable. Our analysis identifies a fitness trade-off in viral evolution that can inform antibody design and therapy protocols.

## Abstract

Broadly neutralizing antibodies are promising candidates for treatment and prevention of HIV-1 infections. Such antibodies can temporarily suppress viral load in infected individuals; however, the virus often rebounds by escape mutants that have evolved resistance. In this paper, we map a fitness model of HIV-1 interacting with broadly neutralizing antibodies using in vivo data from a recent clinical trial. We identify two fitness factors, antibody dosage and viral load, that determine viral reproduction rates reproducibly across different hosts. The model successfully predicts the escape dynamics of HIV-1 in the course of an antibody treatment, including a characteristic frequency turnover between sensitive and resistant strains. This turnover is governed by a dosage-dependent fitness ranking, resulting from an evolutionary trade-off between antibody resistance and its collateral cost in drug-free growth. Our analysis suggests resistance–cost trade-off curves as a measure of antibody performance in the presence of resistance evolution.

HIV-1 infection is characterized by a high turnover rate in combination with high mutation rates, which contribute to the virus’s extraordinary capacity to evade the immune response of the host (1). In return, the immune system explores a broad set of responses, creating a coevolution of the two systems (2, 3). In recent years, the discovery of broadly neutralizing antibodies (bnAbs) created the prospect of antibody-mediated prevention (4, 5) or treatment (6⇓–8) of HIV-1 infection. Such antibodies have been isolated from a minority of HIV-1–infected individuals, termed elite neutralizers, and can potently neutralize a large spectrum of existing HIV-1 strains. They target relatively conserved sites on the envelope protein of the virus (9), blocking virus entry into cells. Additionally, they can engage the host immune system to target infected cells (10). Clinical studies show that in individuals infected with sensitive viral strains, antibody treatment leads to an initial decline in viral load. However, the decline is followed by a rebound on a timescale of weeks, indicating viral escape from antibody neutralization (11⇓⇓–14). Specific insight into the escape dynamics comes from two recent studies recording viral populations in cohorts of HIV-1–infected individuals following infusion of specific bnAbs (12, 13). These studies provide time-resolved in vivo data, including viral load, antibody concentration, and single-virion genome sequences. They reveal complex escape dynamics involving the turnover of several resistant mutant strains.

Here, we establish a biophysically grounded, predictive fitness model for the dominant HIV-1 variant and common escape mutants in response to a bnAb. In the first part of the paper, we use the data in refs. 12 and 13 to infer viral within-host replication rates of sensitive and resistant strains. These rates depend on two key factors: the antibody concentration, which determines the likelihood of antibody–antigen binding, and the virion density, which is subject to saturation effects setting the carrying capacity of the HIV-1 population in a given patient. Both of these densities vary by orders of magnitude in the course of a treatment protocol. Remarkably, we can quantify their effects on viral growth with a fitness model based on a few parameters, which are independent of the host-specific genomic background. This model captures escape patterns, including time-dependent viral load and antibody resistance levels, reproducibly across individuals.

Our fitness model contains two key parameters characterizing a given HIV-1 strain: its growth rate in the absence of antibodies and its resistance to bnAbs, defined as the reduction in the antibody–antigen binding that confers neutralization of viral growth. We infer an in vivo trade-off between these parameters; mutants with increased antibody resistance have a reduced reproductive rate in the absence of antibodies. This trade-off reflects a well-known general phenomenon; resistance mutations can involve collateral costs (15). Mutations can reduce the stability of protein folding and thereby, decrease fitness (16, 17). Moreover, since bnAbs target conserved sites of the virus, mutations in these regions can impair viral function. The cost of escape mutations can be quantified in several ways. Using virus replication assays, Lynch et al. (14) showed that escape from VRC01-class bnAbs targeting the CD4 binding site results in reduced viral replication. Computational methods based on deep mutational scanning (18), deep sequencing of longitudinal samples (19), or multiple sequence alignment (20) can also reveal a fitness cost at target sites. Similar trade-offs between evolution of resistance and function have been observed in microbial systems (21) and for cancer (22). However, to our knowledge, a resistance–growth rate trade-off arising in the in vivo escape from bnAbs has not been inferred so far.

The resistance–growth spectrum inferred here has an important consequence for viral dynamics; it determines a bnAb dosage-dependent fitness ranking of sensitive and resistant strains, which in turn, drives strain turnover during the escape process. We infer a maximum-growth curve that determines the dosage-dependent equilibrium viral load resulting from the sensitive strain and the observed escape mutants. Furthermore, the resistance–growth spectrum can be used to compare different bnAbs and their suitability for treatment.

In the second part of the paper, we show that the fitness model can successfully predict viral escape dynamics during bnAb treatment, given host-specific initial data of viral load and standing variation of resistance mutations prior to bnAb exposure. Predictable features of these dynamics include the strain turnover during the escape process, manifesting in the dosage-dependent ranking of strains in our fitness model. We discuss the consequences of our findings for the evolutionary optimization of antibodies and of time-dependent treatment protocols.

## Results

### Time-Resolved In Vivo Data of bnAb Escape.

Here, we use data from two studies of viral dynamics under bnAb treatment in humans. First, in the study of Caskey et al. (12), 17 individuals infected with HIV-1 received a single intravenous dose of the CD4 binding site–directed bnAb 3BNC117 at various dosages. Second, in the study of Caskey et al. (13), 19 individuals infected with HIV-1 received a single infusion of the V3 loop–directed antibody 10-1074. The total viral load N (measured in RNA copies per milliliter) and the bnAb dosage A (measured in micrograms per milliliter) were tracked over several weeks after the infusion in both studies. Additionally, in the 10-1074 study, single-virion genome sequencing was performed on plasma samples obtained at specific time points, providing frequency estimates of sensitive and mutant strains. The neutralizing power of the antibody against the sensitive strain was measured in terms of the half maximal inhibitory concentration (IC50), which we denote as

We use the data from the 10-1074 study for the 11 individuals who responded to the bnAb infusion, were not on antiretroviral therapy, and had single-genome sequencing performed at three or more time points. Similar pruning has been performed on the 3BNC117 data; *SI Appendix* has details and the full datasets. Sequence analysis and phenotype analyses with pseudoviruses show that the viral escape from neutralization with bnAb 10-1074 can be associated with amino acid changes at one of the gp120 epitope residues 334, 332, or 325 away from the wild-type (wt) allele that is common to all sensitive strains (23, 24). The escape mutations at residue 334 or 332 eliminate a glycosylation site, at which the antibody makes a critical contact to a glycan. Given this functional equivalence, we group them together as mutant 1 (mt1). The escape mutation at residue 325 (mt2) alters a different contact site of the antibody. We will show that the fitness of the sensitive (wt) and resistant (mt1, mt2) strains is largely independent of the genetic background, which differs between viral populations in different host.

### An Ecological Fitness Model for Viral Escape.

HIV-1 replicates in a complex intrahost environment under constraints set by the external bnAb and by the host’s intrinsic immune response. Here, we describe the growth of a viral strain by a continuous birth–death process. The birth term includes the entire replication cycle of virions, including cell entry, replication within host cells, and cell exit. The death term describes clearance of virions from circulation. In a minimal model, a given viral strain i (wt, mt1, or mt2) has an intrahost replication (birth) rate that depends on the antibody dosage and on an effective viral load:^{+} T cells; the level of immune activation depends strongly on the viral load. Third, another host factor, the local depletion of uninfected CD4^{+} T cells, affects replication in a similar way (31) (*SI Appendix*). Autologous constraints act on the effective viral load

Replication and clearance determine the net growth rate or absolute fitness of a given strain:**2** combines the intrinsic growth rate of the virus with ecological pressures from its intrahost environment. These ecological pressures are shaped by the antibody dosage and the autologous immune pressure.

### Ecology Shapes Viral Escape Dynamics.

How do these ecological pressures determine clinical patterns of viral load? First, the autologous constraint sets the equilibrium viral load, or carrying capacity, prior to the start of treatment. This point is defined by the equality of birth and death rates of the sensitive strain, *Methods*). Second, ecology impacts the load through viral evolution; the fitness model determines bnAb dosage-dependent growth rate differences (selection coefficients) that govern the competition between coexisting viral strains. Selection, in turn, shapes strain composition and viral load at future points of the treatment protocol. The minimal model defined by Eqs. **1** and **2** determines specific selection coefficients*SI Appendix*, Fig. S1). In ecology, this dependence is known as density-dependent selection (32, 33). For comparison, we also consider an alternative model with linear niche constraint, **1** and **2** with universal (host-independent) fitness parameters

### Model-Based Ecoevolutionary Dynamics.

The observed escape dynamics of HIV-1 after the infusion of a bnAb show a remarkably complex strain turnover, which involves a sensitive strain and multiple resistant strains. These dynamics depend on the fitness of strains and on the mutation rates that govern the generation of resistant variants from an initially predominant sensitive strain. In the following, we describe time-dependent strain population sizes by a multistrain mutation selection model of the form**1** and **2**. The terms

At a given antibody dosage, the strain dynamics of Eq. **4** lead to a mutation selection equilibrium that is dominated by the fittest strain (in the present system, we will show that this can be wt, mt1, or mt2, depending on dosage). In the saturation model Eqs. **1** and **2**, the equilibrium viral load, or carrying capacity, at a fixed dosage A increases slower than linearly with the dosage-dependent maximum fitness *Methods* and *SI Appendix*).

### Viral Load Trajectories.

A convenient starting point for data analysis is the time-dependent viral load data under bnAb treatment in 11 hosts (Fig. 1*A*). A host indicator (1HB1, etc.) is shown directly above each trajectory in the figure. Prior to the start of treatment, there is a host-specific load **6** in *Methods* (*SI Appendix*, Table S2). Following the infusion at time *SI Appendix*, Fig. S3 and Table S2). In contrast, the relative load shows a common initial response to bnAbs, *B*). Consistently, antibody concentrations in the initial time interval are well above the average wt IC50 concentration *A*). In these hosts, shifting the time axis to a common initial load, we obtain a data collapse describing the rebound to *C*). This indicates a universal growth pattern that extends throughout the approach to the carrying capacity, as given by Eqs. **4** and **6**. The rebound is marked by an initially exponential increase of the load, *SI Appendix*, Fig. S4). Model-based trajectories with host-independent reproduction rate

On the other hand, the initial viral population, which has evolved in the autologous immune environment, turns out to be strongly host specific. By extrapolating the measured (unshifted) rebound curves back to *D*). We find host-specific mutant frequencies of order *SI Appendix*, Table S2), which are higher than frequency estimates obtained from deep sequencing data in two individuals (*SI Appendix*, Tables S5 and S9) and estimates from mutation selection balance (*SI Appendix*). This difference may point to activation of resistant strains in the latent repertoire contributing to seeding the viral rebound (35); such seeding effects are included in our initial condition *D*).

Simulations of the viral escape dynamics Eqs. **1**–**5** highlight the complementary roles of fitness parameters for viral load trajectories (*SI Appendix*, Fig. S2). The host-independent basic fitness parameters

### Bayesian Inference of the Fitness Model.

To infer the full fitness model for the 10-1074 data, we use a Bayesian procedure based jointly on the time series data of viral load and strain frequencies. Optimal fitness parameters as well as host-specific initial mutant frequencies are inferred using a Markov Chain Monte Carlo algorithm that constructs a posterior distribution best fitting the in vivo escape data; details are given in *Methods* and *SI Appendix*.

The maximum-likelihood fitness model has strain-dependent basic replication rates *SI Appendix*, Table S1). In Fig. 2*A*, the resulting basic fitness values *A*). In particular, the maximum-likelihood fitness parameters *SI Appendix*, Fig. S5 and Table S1).

Apart from estimating the fitness parameters in Eqs. **1** and **2**, Bayesian inference can also serve to rank the fitness model components against alternative functional forms. First, the data support the Michaelis–Menten function *SI Appendix*). Second, the data favor the saturation model, *Methods* and *SI Appendix*). Importantly, good fits of the saturation model are obtained with a strain-independent parameter C, which models the host-specific autologous constraint. Because the niche constraint is strain independent, the autologous immune system exerts similar pressure on the wt and the escape mutant strains of the virus.

### Resistance–Cost Trade-Off and Dosage-Dependent Selection.

As shown in Fig. 2*A*, the drug-free growth rate

Fig. 2*B* shows the dosage-dependent growth profile

The growth profiles of Fig. 2*B* exhibit an important consequence of the resistance–cost trade-off; the fitness ranking of viral strains depends on the antibody dosage. The fittest strains are the wt at low dosage, mt2 at intermediate dosage, and mt1 at high dosage. The resulting dosage-dependent maximum-growth rate, *B*); this growth rate determines a dosage-dependent carrying capacity by Eq. **6**. The maximum-growth curve of the virus in the presence of the bnAb is the key immune feature inferred in this paper; it characterizes the potency of the antibody while accounting for the evolutionary response of the virus.

Another aspect of the fitness ranking is a dosage-dependent selection coefficient between the two escape mutants, *C*, we compare the model prediction *Methods* and *SI Appendix*). The observed frequency changes are seen to be in agreement with the predicted functional form (large error bars reflect instances where at least one of the mutants is observed at only one of the two time points, resulting in an inequality for the corresponding data point for

### Predicting Evolutionary Escape Trajectories.

The resistance–cost trade-off and the dosage-dependent fitness ranking of strains lead to an evolutionary prediction of our fitness model; in a bnAb treatment with time-dependent dosage, there is a reproducible turnover between prevalent strains. Here, we test this prediction by training the model on subsets of 10 hosts and using frequency trajectories of the 11th host for prediction and validation; this protocol is detailed in *Methods*.

In Fig. 3*A*, we show the predicted frequency trajectories given by Eq. **4** for 11 validation protocols together with the observed frequency data, which have been excluded from the model training step. The model successfully predicts salient qualitative features of the evolutionary strain turnover. First, the rebound of the viral load, which takes place while bnAb dosage is still high, is predominantly carried by mt1. Exceptions to this pattern are observed in hosts 1HD6K and 1HD10K, where the wt is not fully suppressed in the initial treatment phase, suggesting increased values of

To test the quantitative predictability of frequency trajectories, we predict validation trajectories from one sampling time point to the next and compare model predictions of strain frequency changes, *B*) (hats indicate predicted quantities). The model correctly predicts 22 of 28 instances of frequency increase and 21 of 28 instances of decline, which amount to an overall prediction accuracy of 77%. The main limiting factor of predictability appears to be the incomplete knowledge of low-frequency escape mutations, which can be mitigated by deeper sampling of the initial populations.

## Discussion

In this paper, we establish a fitness model for HIV-1 that successfully predicts the in vivo ecoevolutionary response of the viral system to bnAbs, given simple treatment protocols extending over limited periods. In a broader context, the predictability of evolution under realistic ecological conditions remains a largely open question. In the present system, we show that predictability emerges from a partial universality; viral fitness parameters governing the short-term escape from bnAbs depend only weakly on host environment and viral genomic background, while initial loads and mutant frequencies reflect autologous immune pressures that are variable across hosts. These ecological conditions enter our model as initial data but are beyond the scope of predictions.

Recent mechanistic and genomic fitness models for HIV (20, 36, 37) provide a more detailed picture of the intrahost reproductive dynamics and the viral genome sites relevant for resistance evolution. In contrast, our model gives a coarse-grained description of the viral dynamics geared to inference from limited in vivo time series data. Despite its simplicity, the model captures viral escape from bnAbs across sets of hosts with varying immunological responses and genetic variation of the viral population, as observed for two bnAbs included in our study.

Our inference of viral fitness factors suggests biological features of antigen–antibody interactions. First, the resistance mutations analyzed in this study operate largely independently of their genetic background, which differs between viral populations in different hosts. Second, the host-to-host variation in autologous immune response can, over the limited time intervals of our treatment protocols, be absorbed into a host-dependent constraint parameter C. This parameter is uniform across strains, indicating that wt and escape mutants are in the same niche of the host’s immune environment. In other words, we do not find evidence that the bnAb escape mutants studied here confer escape from autologous immune control. Third, we find evidence of density-dependent selection on escape mutations, suggesting a specific form of autologous immune suppression; a viral strain of increased basic fitness induces a supralinear increase in activation of the immune system, generating, in turn, a sublinear increase of the viral load as observed in our data.

A striking feature of bnAb treatment protocols is the complex strain turnover observed already over short timescales. Our fitness model explains this turnover in terms of time-dependent bnAb dosages together with a dosage-dependent fitness ranking between strains. This interplay is likely to be generic; bnAb treatment protocols generate antigenic fitness seascapes driving time-dependent strain prevalence. Given bnAb exposure over longer timescales, the strain dynamics is expected to become even more complex. In particular, compensatory mutations can stabilize mutant strains by reducing the fitness cost of antibody escape (14, 17).

The key fitness characteristic driving viral escape evolution is a trade-off between resistance against bnAb neutralization and growth in the absence of antibody challenge (Fig. 2*A*). This finding has implications for the optimization of bnAbs, which is a topic of high current interest (38⇓–40). The standard procedure to measure the power of bnAbs is based on neutralization assays against a panel of reference strains, which determine the potency (i.e., the mean IC50 against susceptible strains) and the breadth (i.e., the number of reference strains against which the IC50 is lower than a threshold value) of neutralization. However, these measures do not take into account the specific genetic changes that carry resistance evolution against that antibody. In contrast, the maximum-growth curve,

The predictability of the viral escape dynamics is also an important tool for designing optimal bnAb dosage protocols. Any such protocol has to balance medical limitations and physiological collateral costs of treatment with their effect of curbing the viral load. In particular, the effectivity of multistep protocols crucially depends on timing and strength of individual applications (41, 42). The optimization of protocols has to take into account the full ecoevolutionary response of the viral population, including the turnover in strain prevalence. Here, we have established a proof of principle that this response is computable for realistic treatment settings. Extending the method to other bnAbs and to more complex, resistance-limiting therapies with multiple bnAbs is an important avenue for future work.

## Methods

### Inference of Ecoevolutionary Dynamics.

In the Bayesian inference of the fitness model, we use the following input data: time series data of the total population size N, time series data of strain frequencies Y, and time series data of host-specific antibody dosages A. These data are used to infer the fitness parameters F (including the replication rates, clearance rate, and antibody resistances), initial strain frequency data X, and host-dependent niche constraint parameters C. To this end, we construct an error model that gives the log likelihood *SI Appendix*.

### Model Comparison.

In the main text, we compare two ecological fitness models: the saturation model defined by Eqs. **1** and **2** and a model with a linear constraint of strength *SI Appendix*). Second, we compare the stationary viral load. In the saturation model, we obtain*SI Appendix*); comparison with the data again favors the saturation model.

### Dosage-Dependent Selection.

Selection coefficients can be inferred from observed frequency trajectories *C* against the average bnAb dosage A in the same time window. Details are given in *SI Appendix*.

### Predicting Viral Escape Evolution.

In order to predict the escape evolution for a given individual α, we use the following scheme. First, we divide the data into a training set *C* and *SI Appendix*, Table S2); we set **1**–**3** using the fitness parameters inferred from the training set. This results in the predicted trajectories *SI Appendix*.

## Data Availability

All study data are included in the article and/or *SI Appendix*.

## Acknowledgments

We acknowledge discussions with Armita Nourmohammed and Colin Lamont. This work was supported by Deutsche Forschungsgemeinschaft Grant SFB 1310.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: mlaessig{at}uni-koeln.de.

Author contributions: M.M., F.K., and M.L. designed research; M.M., K.V., H.G., F.K., and M.L. performed research; M.M., K.V., H.G., F.K., and M.L. analyzed data; and M.M., K.V., H.G., F.K., and M.L. wrote the paper.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2104651118/-/DCSupplemental.

- Copyright © 2021 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- R. M. Lynch et al.

- ↵
- K. J. Bar et al.

- ↵
- ↵
- F. Klein et al.

- ↵
- Cl. Lu et al.

- ↵
- ↵
- ↵
- ↵
- R. M. Lynch et al.

- ↵
- ↵
- C. S. Wylie,
- E. I. Shakhnovich

- ↵
- ↵
- ↵
- ↵
- R. H. Y. Louie,
- K. J. Kaczorowski,
- J. P. Barton,
- A. K. Chakraborty,
- M. R. McKay

- ↵
- D. I. Andersson,
- D. Hughes

- ↵
- ↵
- H. Mouquet et al.

- ↵
- A. Dingens,
- D. Arenz,
- H. Weight,
- J. Overbaugh,
- J. D. Bloom

- ↵
- ↵
- ↵
- V. Mustonen,
- J. Kinney,
- C. G. Callan Jr,
- M. Lässig

- ↵
- ↵
- A. Rotem et al.

- ↵
- ↵
- ↵
- J. Mallet

- ↵
- ↵
- J. Crow,
- M. Kimura

- ↵
- Y. Z. Cohen et al.

- ↵
- ↵
- ↵
- J. S. Shaffer,
- P. L. Moore,
- M. Kardar,
- A. K. Chakraborty

- ↵
- K. G. Sprenger,
- J. E. Louveau,
- A. K. Chakraborty

- ↵
- V. Sachdeva,
- K. Husain,
- S. Wang,
- A. Murugan

- ↵
- ↵
- M. Lässig,
- V. Mustonen

## Citation Manager Formats

## Article Classifications

- Biological Sciences
- Evolution