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

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 (measured in RNA copies per milliliter) and the bnAb dosage (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 . Here, we infer analogous IC50 values for mutant strains and use them as a measure of bnAb resistance.

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; 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 (wt, mt1, or mt2) has an intrahost replication (birth) rate that depends on the antibody dosage and on an effective viral load:In this model, the basic replication rate, , is reduced by three factors characterizing the intrahost environment. First, we use a simple Hill model of neutralization; functional antibody binding prevents new cell infections by constraining replication to unbound virions. The fraction of unbound virions, , depends on the antibody dosage and the strain-specific antibody resistance , defined as the dissociation constant of binding that confers neutralization. Similar biophysical fitness models linking functional binding and growth have been established in other microbial and viral systems (25–30). Second, the replication rate of the virus is reduced by a saturation factor characterizing the autologous immune pressure from B cells that target free virions in the blood or cytotoxic T cells that target infected CD4+ 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) (). Autologous constraints act on the effective viral load , defined as the total number of infective virions, which are not neutralized by bnAbs. Neutralized, bnAb-bound virions do not count for these saturation effects because they do not contribute to cell infection. In the minimal model, these processes define a strain-independent parameter of autologous constraint; this simplification will be justified a posteriori from our inference procedure. We describe clearance of virions, either by active immune processes or by decay, by a single clearance (death) rate , which will turn out to be host independent.

Replication and clearance determine the net growth rate or absolute fitness of a given strain:The basic fitness in the low-density regime (i.e., in the absence of saturation effects and autologous immune pressure) is denoted by the shorthand . In summary, the fitness defined in Eq. 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, , and can be computed in the minimal model (). 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. and determines specific selection coefficientsImportantly, because of the saturation effect of autologous immune response, these selection coefficients depend on the viral load (). In ecology, this dependence is known as density-dependent selection (32, 33). For comparison, we also consider an alternative model with linear niche constraint, , which defines selection coefficients independent of the viral load (32). In the following, we will show that the viral escape data are best described by the saturation model Eqs. and with universal (host-independent) fitness parameters , , , and host-specific autologous constraint .

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 formwhich includes the fitness model defined by Eqs. and . The terms describe the mutational turnover between strains given by , , and , taking into account the multiplicity of amino acid changes and assuming a uniform point mutation rate, (19). These dynamics determine the time-dependent viral load, , as well as the evolution of strain frequencies, , which follow a Wright–Fisher population dynamics under time-dependent selection (34). We note that the specific resistance mutants, the corresponding mutation rates, and hence, the initial strain frequencies depend on the antibody used for treatment.

At a given antibody dosage, the strain dynamics of Eq. 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. and , the equilibrium viral load, or carrying capacity, at a fixed dosage increases slower than linearly with the dosage-dependent maximum fitness , reflecting the load-dependent increase of autologous immune pressure. In contrast, the linear constraint model sets a carrying capacity that is proportional to the fitness ( and ).

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 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 , which we identify with the carrying capacity of HIV in a drug-free environment. In our model, the variation of initial loads translates into host-specific constraint parameters , as given by Eq. in (). Following the infusion at time , the antibody concentration decays exponentially,with host-specific initial values in the range g/mL and a characteristic decay time d (). In contrast, the relative load shows a common initial response to bnAbs, , from which we infer a universal clearance rate of the wt strain, (Fig. 1). Consistently, antibody concentrations in the initial time interval are well above the average wt IC50 concentration g/mL (13), which implies that most virions are bound to antibodies and cannot replicate. The rebound of the viral load is driven by the bnAb-resistant strain mt1 in 9 of 11 hosts (see the data points in Fig. 3). In these hosts, shifting the time axis to a common initial load, we obtain a data collapse describing the rebound to % of the initial load (Fig. 1). This indicates a universal growth pattern that extends throughout the approach to the carrying capacity, as given by Eqs. and . The rebound is marked by an initially exponential increase of the load, , from which we estimate a universal basic growth rate (we will show below that the growth of mt1 is independent of since ). The analysis of viral load trajectories from the 3BNC117 dataset, which targets the more complex CD4 binding site, produces consistent results (). Model-based trajectories with host-independent reproduction rate and clearance rate consistently describe the viral load dynamics of nine individuals after the bnAb infusion. Together, this gives evidence that a universal growth pattern across different hosts may be a general feature of HIV-1 escape from bnAbs.

Fig. 1.

Viral load trajectories have universal growth parameters. (A) Observed time series of the viral load in 11 individuals (RNA copies per milliliter; red dots) are shown together with the load trajectory of the maximum-likelihood fitness model. (B) Collapse plot of the initial load decline. Individual hosts are indicated by color. Measured relative load (solid lines) and universal (host-independent) exponential fit with inferred clearance rate (dashed line). (C) Collapse plot of the load rebound. Measured relative load plotted against the time from a common initial value (solid lines), universal fit curve (dotted line), and exponential fit to the initial rebound with inferred mutant growth rate (dashed line). (D) Inference of the initial mutant frequency. The measured relative load is plotted against the time from the start of treatment at (solid lines). Extrapolation of the exponential rebound back to (dashed lines) provides estimates of the initial frequencies (intercept with the vertical axis).

Fig. 3.

Prediction of escape evolution. (A) Observed time series of strain frequencies (dots; bars indicate sampling errors) are shown together with predicted frequency trajectories for 11 validation protocols (lines). The host indicator of each trajectory is displayed directly above each axis. Fitness parameters used for predictions are obtained from complementary training sets. The first data point for each strain (open circles) is used as the initial condition, and the subsequent points (filled dots) are to be compared with predictions. When mt2 is first observed, the normalization of the predicted trajectories is updated. (B) Model predictions of strain frequency changes, , are plotted against the corresponding observed changes, . Frequency ratios with () below the sampling threshold are evaluated with pseudocounts. Frequency increase is correctly predicted in 22 of 28 instances (first quadrant), and frequency decline is correctly predicted in 21 of 28 instances (third quadrant).

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 , we can infer the mutant frequency in the unperturbed viral population at the start of the treatment with bnAb 10-1074 (Fig. 1). We find host-specific mutant frequencies of order to (), which are higher than frequency estimates obtained from deep sequencing data in two individuals () and estimates from mutation selection balance (). 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 (Fig. 1).

Simulations of the viral escape dynamics Eqs. – highlight the complementary roles of fitness parameters for viral load trajectories (). The host-independent basic fitness parameters and set the slopes of decline and rebound, respectively. For given values of and , the initial mutant frequency and the niche constraint independently determine the minimum load and the stationary load after rebound.

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 and .

The maximum-likelihood fitness model has strain-dependent basic replication rates , resistance parameters , and a strain-independent clearance rate (). In Fig. 2, the resulting basic fitness values in the absence of antibodies are plotted against for all strains; error bars indicate 95% CIs obtained from the posterior distribution. The universal maximum-likelihood fitness model reproduces the viral load trajectories across all 11 hosts (blue lines in Fig. 1). In particular, the maximum-likelihood fitness parameters and and initial frequencies agree well with the corresponding values obtained by fitting decline and rebound of the viral load ().

Fig. 2.

Resistance–cost trade-off and dosage-dependent fitness ranking of viral strains. (A) The inferred drug-free growth rate, , is plotted against the antibody resistance of the viral strains = wt, mt1, mt2 (maximum-likelihood values; error bars indicate 95% CIs). (B) Michaelis–Menten growth profiles of the strains = wt, mt1, mt2 interpolate between the basic growth rate and the clearance rate with an IC50 concentration (solid lines). The maximum-growth rate (dashed line) is the envelope of the growth profiles of individual strains. (C) Dosage-dependent selection between resistance mutants. Estimates of the selection coefficient obtained from relative frequency changes (dots; bars indicate sampling errors, and large bars indicate inequalities involving frequencies below the sampling threshold) are compared with the predicted Michaelis–Menten form (yellow line).

Apart from estimating the fitness parameters in Eqs. and , Bayesian inference can also serve to rank the fitness model components against alternative functional forms. First, the data support the Michaelis–Menten function linking growth to antibody density; an alternative model with a fitness cost linear in has a significantly lower likelihood (). Second, the data favor the saturation model, , against the alternative model with linear niche constraint. The saturation model is also supported by the observed nonlinearity of the carrying capacity; the observed load ratio is larger than the fitness ratio that would set the load ratio given a linear niche constraint ( and ). Importantly, good fits of the saturation model are obtained with a strain-independent parameter , 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, the drug-free growth rate and the antibody resistance vary between strains in a correlated way; higher resistance (i.e., weaker binding to bnAbs) implies slower growth in the absence of antibodies. This resistance–cost trade-off is intuitive; more drastic changes to a viral protein (here, gp120) can more effectively reduce bnAb binding but also, have larger impact on protein stability and/or reproductive functions, such as binding to host cells. Here, we infer a trade-off of considerable amplitude; the cost of escape, defined as , is about 50% for both escape mutants.

Fig. 2 shows the dosage-dependent growth profile of the wt and escape mutant strains, as obtained from our maximum-likelihood fitness model. For a given strain, the growth rate takes a sigmoid (Michaelis–Menten) form that interpolates between the asymptotic values in the low-binding regime () and in the strong-binding regime (), with a cross-over at the half-binding point (IC50 concentration, ). Specifically, mt1 is in the weak-binding regime throughout (i.e., in all hosts), and no dosage dependence of its fitness is inferred from this dataset.

The growth profiles of Fig. 2 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, , is the envelope of the growth profiles of individual strains (gray line in Fig. 2); this growth rate determines a dosage-dependent carrying capacity by Eq. . 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, , which changes sign at a specific dosage g/mL. In Fig. 2, we compare the model prediction (yellow line) with data points obtained from changes in the frequency ratio between subsequent sampling time points, plotted against the average bnAb concentration in between these points ( and ). 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 .

In Fig. 3, we show the predicted frequency trajectories given by Eq. 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 and a limitation of fitness universality across individuals. Second, at intermediate times and bnAb dosages, mt2 appears and qualitatively follows the predicted trajectories. However, we do not attempt to predict the appearance of mt2 in the observable frequency range because its low-frequency dynamics appear to be more susceptible to the host-specific environment, and more initial-state sequencing would be necessary to determine its initial frequency . Instead, when mt2 appears, this observation is integrated into the initial data for predictions of the subsequent time points (the corresponding renormalization of frequencies is seen as jumps of the model trajectories). Third, a rebound of the wt occurs at late times and low bnAb dosages. This rebound is observed in almost all protocols where strain frequencies have been tracked over a sufficiently long time interval ( d).

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, , with observed changes in the same time interval, (Fig. 3) (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 . 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). 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, , inferred in this paper determines the bnAb dosage-dependent carrying capacity resulting from the evolution of sensitive strains and escape mutants. Conceptually, this result underscores that the ecology and evolution of the antigen population are critical to describe clinical patterns of time-dependent viral load and resistance against antibodies. The maximum-growth curve also sets a necessary criterion for clinically suitable bnAbs; at therapeutically sustainable dosages , the growth of common resistance mutants must be suppressed, . Importantly, resistance–cost trade-offs and maximum-growth curves can be used to compare the performance of bnAbs targeting different viral epitopes, including antibody combinations.

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 , time series data of strain frequencies , and time series data of host-specific antibody dosages . These data are used to infer the fitness parameters (including the replication rates, clearance rate, and antibody resistances), initial strain frequency data , and host-dependent niche constraint parameters . To this end, we construct an error model that gives the log likelihood and minimizes the error with respect to , , and . Details are given in .

Model Comparison.

In the main text, we compare two ecological fitness models: the saturation model defined by Eqs. and and a model with a linear constraint of strength , where fitness takes the form . First, the log likelihood obtained from our Bayesian inference favors the saturation model (details are given in ). Second, we compare the stationary viral load. In the saturation model, we obtainwhere we assume that the population contains a dominant fittest strain with parameters , , and . This load depends in a sublinear way on the fitness . In contrast, the stationary load of the linear model,is proportional to the replication rate . This leads to different predictions for the load ratio between mutant and wt strains (); comparison with the data again favors the saturation model.

Dosage-Dependent Selection.

Selection coefficients can be inferred from observed frequency trajectories and . The log frequency ratio, , changes with the instantaneous selection coefficient,where the selection coefficient . Then, the observed log frequency change in a time interval () measures the average selection coefficient,which is plotted in Fig. 2 against the average bnAb dosage in the same time window. Details are given in .

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 that contains data of individuals and a test set that contains data of individual . Then, the fitness parameters are inferred from the training set using the Bayesian inference scheme of the fitness model. To predict the viral escape for individual , it is necessary to first obtain the host-specific parameters. The initial load (which sets the niche constraint ) and initial dosage are given by measurements, while the initial frequency is given by deep sequencing where available (in patients 1HB3 and 1HD1) or inferred from backward extrapolation of the load rebound of using a fitness parameter inferred from the training set (Fig. 1 and ); we set . The extrapolation procedure is a substitute for direct deep sequencing of the genetic variation in the initial strain population, which is not available for the whole dataset. The trajectories can be computed using Eqs. – using the fitness parameters inferred from the training set. This results in the predicted trajectories that are compared with the data. For validation, we compare the predicted frequency change between subsequent sampling points with the corresponding observed changes . Details are given in .