Nicolas Chéron #

Adrian W R Serohijos #

Jeong-Mo Choi

Eugene I Shakhnovich
email: shakhnovich@chemistry.harvard.edu

Evolutionary dynamics of viral escape under antibody stress: A biophysical model

Keywords: Viral evolution, optimal mutation rate, neutralizing antibody, folding stability

Content of supplementary material: (1)

Table of experimentally observed mutation rates and burst size of some RNA viruses.

(2) Influence of biophysical parameters G and Kd on Pinf.. (3) Viral population dynamics. (4) Influence of the genome size. (5) Influence of the burst size (M) on the phase diagram.

Introduction

Viruses face the stress of various environments [START_REF] Meyers | Fighting change with change: adaptive variation in an uncertain world[END_REF] . It has been hypothesized that an adaptive strategy among viral populations to fight this pressure is to modulate their mutation rates [START_REF] Combe | Variation in RNA virus mutation rates across host cells[END_REF] and previous works have studied the existence of an optimal mutation rate for viruses [START_REF] Kamp | Viral evolution under the pressure of an adaptive immune system: Optimal mutation rates for viral escape[END_REF] . An upper bound to the rate of mutation exists due to lethal mutagenesis whereby mutational load leads to a population decline [START_REF] Bull | Theory of lethal mutagenesis for viruses[END_REF][START_REF] Chen | Lethal mutagenesis in viruses and bacteria[END_REF][START_REF] Zeldovich | Protein stability imposes limits on organism complexity and speed of molecular evolution[END_REF] . Lethal mutagenesis is the presumed mode of action of some anti-retroviral drugs, such as ribavirin, that increase the mutation rate to lethal levels [START_REF] Crotty | The broad-spectrum antiviral ribonucleoside ribavirin is an RNA virus mutagen[END_REF][START_REF] Crotty | Implications of high RNA virus mutation rates: lethal mutagenesis and the antiviral drug ribavirin[END_REF] . A lower bound for mutation rate has been hypothesized for several reasons. First, there is the requirement for an elevated mutation rate to increase the supply of beneficial mutations. These mutations could help the virus to adapt to the changing environmental stress that can come from the immune response of the host cell. This scenario is supported by the observed hypermutability of bacteria isolated from natural populations [START_REF] Oliver | High frequency of hypermutable Pseudomonas aeruginosa in cystic fibrosis lung infection[END_REF][START_REF] Bjorkholm | Mutation frequency and biological cost of antibiotic resistance in Helicobacter pylori[END_REF] . Another potential explanation of the lower bound on the mutation rate is the high cost of maintaining replication fidelity [START_REF] Combe | Variation in RNA virus mutation rates across host cells[END_REF][START_REF] Furio | The cost of replication fidelity in an RNA virus[END_REF][START_REF] Sniegowski | Evolution: constantly avoiding mutation[END_REF] , although such cost remains to be experimentally demonstrated. The existence of an optimal mutation rate due to the balance for the requirement for beneficial mutations and the avoidance of the lethal mutagenesis have been postulated in the past [START_REF] Kamp | Viral evolution under the pressure of an adaptive immune system: Optimal mutation rates for viral escape[END_REF] . There have been several theoretical models on how such an optimal mutation rate could arise [START_REF] Kamp | Viral evolution under the pressure of an adaptive immune system: Optimal mutation rates for viral escape[END_REF][START_REF] Sniegowski | Evolution: constantly avoiding mutation[END_REF][START_REF] Regoes | Viral mutation rates: modelling the roles of within-host viral dynamics and the trade-off between replication fidelity and speed[END_REF][START_REF] Sniegowski | The evolution of mutation rates: separating causes from consequences[END_REF][START_REF] Raynes | Experimental evolution and the dynamics of genomic mutation rate modifiers[END_REF] . However, these models are largely phenomenological and lack molecular realism. Where this balance occurs and what is the underlying molecular (or cellular) mechanism is still being debated both from theoretical and experimental standpoints because the structure of the fitness landscape, and consequently the distribution of fitness effects (DFE), is unclear. Assuming the DFE a priori is challenging because the DFE itself can change during evolution, as demonstrated in experiments [START_REF] Silander | Understanding the evolutionary fate of finite populations: the dynamics of mutational effects[END_REF] and in forward evolutionary simulations on biophysics-based fitness landscapes (in the context of protein folding and binding landscape [START_REF] Serohijos | Contribution of selection for protein folding stability in shaping the patterns of polymorphisms in coding regions[END_REF][START_REF] Heo | Interplay between pleiotropy and secondary selection determines rise and fall of mutators in stress response[END_REF] or in the context of transcription factor binding [START_REF] Mustonen | From fitness landscapes to seascapes: non-equilibrium dynamics of selection and adaptation[END_REF] ).

We present here an alternative approach by projecting the fitness landscape on the biophysical properties of the capsid. There is experimental support on what the mapping to fitness and protein properties could be in the context of viruses [START_REF] Gong | A biophysical protein folding model accounts for most mutational fitness effects in viruses[END_REF] . The biophysics-based fitness landscape has also been used to explain the measured DFE among viruses 21 and evolution of mutation rates [START_REF] Heo | Interplay between pleiotropy and secondary selection determines rise and fall of mutators in stress response[END_REF][START_REF] Heo | Emergence of species in evolutionary "simulated annealing[END_REF] . More importantly, the distribution of changes induced by random mutations on key biophysical properties such as folding stability are well-defined and established [START_REF] Tokuriki | The Stability Effects of Protein Mutations Appear to be Universally Distributed[END_REF][START_REF] Zeldovich | Protein stability imposes limits on organism complexity and speed of molecular evolution[END_REF] . Thus, the DFE does not need to be an assumption in the evolutionary model. Several groups including our own recently showed that a biophysical approach could account for the influences of mutations on proteins and could have an effect on the distribution of fitness effects and the dynamics of adaptation [START_REF] Chen | Lethal mutagenesis in viruses and bacteria[END_REF][START_REF] Zeldovich | Protein stability imposes limits on organism complexity and speed of molecular evolution[END_REF][START_REF] Serohijos | Contribution of selection for protein folding stability in shaping the patterns of polymorphisms in coding regions[END_REF]21,[START_REF] Goldstein | The evolution and evolutionary consequences of marginal thermostability in proteins[END_REF][START_REF] Bloom | Thermodynamics of neutral protein evolution[END_REF] . Our goal here is to provide the simplest model, although sophisticated enough at the level of biophysics, to determine the underlying molecular mechanism to investigate the contribution of biophysics to evolution. We investigate in this report the balance between two competing forces (the lethal mutagenesis and the maintenance of supply of beneficial mutations) during the evolution of viral populations under antibody stress by providing biophysical details into the molecular effects of mutations.

Methods

Biophysical model of viral evolution

To determine the interplay between the viral mutation rate, its demography, and the biophysical properties of viral proteins, we have decided to mimic in silico a typical "serial passaging" experiment (Figure 1-a). During these experiments, a fraction of an initial viral stock is added to a medium containing cells and other stressors (e.g. drugs or antibodies). The virus can then infect the cells under an external pressure (the drugs or the antibodies) and new viral particles are released, giving rise to a new stock. These steps constitute a single passage. After multiple passages, the virus titer can then be measured or the virus can be sequenced and compared to the initial viral stock. Our primary concern here is to model experimental evolution of viruses that use serial passaging; a premise is that the results from these kinds of studies are applicable to the evolution of natural populations. What is the appropriate demographic description of natural viral evolution can be debated. There is a claim that periodic bottlenecking resembles the interhost infection, whereby only a few viral particles are transmitted between hosts.

We start the computational model with a stock of N viral particles, from which 10 4 particles are randomly chosen to infect cells according to their fitness (sampling size) [START_REF] Wobus | Murine norovirus: a model system to study norovirus biology and pathogenesis[END_REF] . The fitness of a virus is considered to be proportional to its probability of infection Pinf which itself is assumed to be a function of the folding stability of all viral proteins and of the affinity to the antibody. Thus, this model applies to lytic cycles but not to lysogenic cycles where the integration of the viral genome in the host genome would require other parameters in the fitness function. We also assume that the multiplicity of infection (MOI) is equal to 1, meaning that each viral particle can encounter at least one host cell. If a cell is infected, M new particles are released. In our definition, the "burst size" M is the number of progenies after a successful infection of a mammalian cell. During replication, each new particle can be mutated at a specified per genome rate .

In our model, the fitness is expressed as the probability of infecting a cell (Pinf). We assumed that for a viral particle to infect a cell, the viral proteins must be folded (to be functional) and be free from the neutralizing antibody (while all proteins encoded by the viral genome must be folded, only one is considered to bind to the antibody). Considering a thermodynamic equilibrium between the capsid, the antibody and the complex (with an equilibrium constant of Kd), one can write if p viral proteins are encoded by the viral genome (with G=Gfolded-Gunfolded and =1/RT):

P inf = e bDG 1 1+ e bDG 1 + e bDG 1 [Ab] K d ae è ç ö ø ÷ ae è ç ç ç ç ö ø ÷ ÷ ÷ ÷ ´1 1+ e bDG k ae è ç ö ø ÷ k=2 p Õ (1)
assuming a two-state folding thermodynamics [START_REF] Privalov | Stability of proteins: small globular proteins[END_REF] . Here we distinguish between epitope-containing capsid protein, which we denote as the first viral protein, and the rest of viral proteins which have to be folded to function but do not interact with antibody. The first factor represents the probability that the epitope-containing capsid protein is folded and not sequestered by the antibody Ab. Assuming that the antibody can bind only to the folded capsid protein, there are three distinct species: folded and free capsid, folded and bound capsid, and unfolded capsid.

Thus, the fraction of folded and bound is determined by the ratio of the Boltzmann probability of being folded and the total probability of being in any of the three states. The second factor is the probability that the remaining p-1 viral proteins are in their native (functional) states under thermodynamic equilibrium. Thus, we do not rely on an assumption regarding the distribution of fitness effects as it has often been done before. Instead, we rely on a biophysical model of the fitness, for which data can be derived from theory or experiments. Here, the viral fitness depends on two biophysical traits of the viral proteins -the Gibbs free energy of folding of each protein (G) and the binding constant of the capsid protein to the antibody (Kd)-and on the external pressure (the concentration of antibody in the medium [Ab]).

In this model, we don't explicitly consider the capsid assembly, which would require additional assumptions on the strength of the interaction between individual capsid domains. The equation also assumes that the antibody could target only one of the capsid domains. The motivation for the latter assumption is that each antibody is large and it probably occludes other antibodies from binding.

Biophysical effects of mutations

When a mutation occurs in the viral genome, we estimate its effect on the folding Gibbs free energy in the following way: Gmutant = Gwildtype + G. The free energy of folding of the wild type (Gwildtype) for all proteins is set initially to -5.0 kcal/mol, which is close to the average folding stability of real proteins [START_REF] Kumar | ProTherm and ProNIT: thermodynamic databases for proteins and protein-nucleic acid interactions[END_REF] . This initial choice is not relevant to the overall results because the population will equilibrate to the folding stability dictated by mutation-selection balance [START_REF] Serohijos | Merging molecular mechanism and evolution: theory and computation at the interface of biophysics and evolutionary population genetics[END_REF] .

The biophysical effect G is drawn from a normal distribution with a mean value of DDG mean = -0.13 DG mean ( ) + 0.23 kcal/mol and a standard deviation of Gsd = 1.7 kcal/mol (ref. [START_REF] Zeldovich | Protein stability imposes limits on organism complexity and speed of molecular evolution[END_REF][START_REF] Tokuriki | The Stability Effects of Protein Mutations Appear to be Universally Distributed[END_REF] ). These values are derived from the distribution of experimentally measured G and G values collated in the ProTherm database [START_REF] Kumar | ProTherm and ProNIT: thermodynamic databases for proteins and protein-nucleic acid interactions[END_REF][START_REF] Adrian | Protein Biophysics Explains Why Highly Abundant Proteins Evolve Slowly[END_REF] . The linear dependence between Gmean and G is also derived from the ProTherm database [START_REF] Adrian | Protein Biophysics Explains Why Highly Abundant Proteins Evolve Slowly[END_REF] . It models the sequence depletion whereby arising mutations are on average more destabilizing when the wild-type protein is more stable, because there are fewer sequences that could lead to stabilization. That is, the more stable the protein is, the more likely it is to find destabilizing mutations [START_REF] Serohijos | Protein biophysics explains why highly abundant proteins evolve slowly[END_REF] .

Whereas all mutations can affect the free energy of folding, it is not the case for the binding affinity. The average number fraction of residues that participate in protein-protein interactions is approximately ~15% according to a bioinformatics study [START_REF] Hu | Conservation of polar residues as hot spots at protein interfaces[END_REF] . Thus, we assumed that the probability that a mutation hits a residue belonging to the epitope of the viral capsid epitopecontaining protein is 0.15. We also assumed that the perturbation of the protein-protein interaction is statistically analogous to mutational effect on the folding Gibbs free energy (i.e. the distribution of Gbinding is the same as the one for Gfolding):

( )

This equation reflects the assumption that the energetics that drive protein-protein interaction are largely similar to energetics that drive protein folding. Thus, the central part of our algorithm is an assumption regarding the distribution of biophysical effects of mutations. The average influence of G and of Kd on Pinf is shown in Supporting Information (see Figure S1).

Model system

Our model is not explicitly sequence-dependent because the effects of mutations are inferred from the distributions of Gfolding and Gbinding, but we account for the number of proteins and the genome length. We have used values for the murine norovirus for which three proteins are important during infection: VP1 (541 amino-acids), VP2 (208 amino-acids) and VPg (124 aminoacids) [START_REF] Hardy | Norovirus protein structure and function[END_REF] . The remaining proteins are only expressed once the virus has penetrated into the cell and were not considered in the model (thus, a value of p=3 was used, see Eq. 1). Realistic values of burst size for RNA viruses are in the range of 50 to ~1700 (ref. [START_REF] Sanjuan | Viral mutation rates[END_REF] ); however, to save computational time a value of M=100 was used unless otherwise noted.

Results

Evolution under antibody stressor

We first show the dynamics of the viral population aimed at mimicking the serial passaging experiments (Figure S2). Starting with 10 5 viral particles and no external pressure ([Ab]=0), the viral count first increases and quickly equilibrates. Since we start the simulations with very stable proteins (Gfolding=-5 kcal/mol), Pinf is close to 1 and all cells are infected; thus, the viral stock quickly increases. However, the dynamics of this initial equilibration is irrelevant because the following results are a consequence of the evolutionary dynamics already under equilibrium.

the distribution of selection coefficient recovers its shape prior to the subsequent change of antibody concentration.

To generalize this result, we explored the onset of extinction for various combinations of mutation rate and antibody concentration. For a given value of the mutation rate , we determined the highest concentration of antibody that still leads to viral escape. First, we performed evolutionary simulations to allow the population to reach mutation-selection balance at the chosen mutation rate but at zero antibody concentration. This equilibration phase was performed for 20 passages since we observed that this value is enough to equilibrate the virus stock (see for example Figure S2 where it can be seen that the steady-state has been reached at passage 20). Then, we stressed the population at the chosen antibody concentration. Extinction or survival is evaluated after 5 passages; we found that 5 passages are enough to observe the survival of a viral population at a given antibody concentration. If the virus has escaped, the viral stock increases at each step and there is no need to go beyond 5 passages.

The value of antibody concentration at which the virus escapes is an indication of how much pressure it can endure and still survive. By varying the mutation rate , we determined the phase diagram outlining the conditions at which we observed extinction and escape of the virus under antibody stress (see Figure 3). Under antibody pressure, increasing the mutation rate increases the likelihood of acquiring mutations that lower the binding free energy of the proteinantibody interaction, and then lead to escape. However, at sufficiently high mutation rates, further increase in mutation rate lead to lethal mutagenesis, because arising mutations decrease folding stability, and subsequently fitness by lowering the fraction of folded proteins Pf as given by Eq.1 (ref. [START_REF] Chen | Lethal mutagenesis in viruses and bacteria[END_REF][START_REF] Zeldovich | Protein stability imposes limits on organism complexity and speed of molecular evolution[END_REF][START_REF] Serohijos | Contribution of selection for protein folding stability in shaping the patterns of polymorphisms in coding regions[END_REF]21,[START_REF] Goldstein | The evolution and evolutionary consequences of marginal thermostability in proteins[END_REF] ). Under very low mutation rate, proteins are more stable, but the chance of acquiring escape mutation is decreased. Thus, there is a balance between lethal mutagenesis and maintenance of supply of escape mutations which is quantified by the phase diagram in Figure 3.

We find that under our model (which provides biophysically realistic estimates of the molecular effects of mutations) there exists a mutation rate at which virus can sustain highest stress. This mutation rate for an RNA virus in our model is found to be around log10(g)=0.92 or 8.3 mutations per genome per replication (for a burst size of M=100). This value is somewhat above the range of experimentally observed mutation rates for RNA virus (see Table S1 in Supporting

Information and Discussion below). It must be pointed out that the existence of a mutation rate allowing a maximal amount of stress does not prove that the population can reach this mutation rate, neither that this mutation rate is always optimal in a broad range of conditions; such a study will be the subject of a following paper. The lethal mutagenesis threshold is calculated to be log10(g)=1.69 or 48.8 mutations per genome per replication, which is higher than previously reported maximum value of lethal mutagenesis [START_REF] Chen | Lethal mutagenesis in viruses and bacteria[END_REF][START_REF] Zeldovich | Protein stability imposes limits on organism complexity and speed of molecular evolution[END_REF]21 : one must note that in these previous models, the burst size corresponding to conservative or semi-conservative replication was fixed at M=2.

Thus, the lethal mutagenesis occurs here at a higher value of the mutation rate due to much greater burst size compared to previous models.

Influence of the genome size

We next determined the effect of the genome size on the optimal mutation rate. In addition to the norovirus model with three genes detailed previously, we also performed simulations with only one of the genes. The phase diagrams of the four models are presented in Supporting Information (Figure S3). The peak in the phase diagram with respect to the mutation rate per genome is almost constant; however, for the mutation rate per nucleotide, an influence of the size of the genome is observed. For 124 amino acids, it is found to be log10(n)=-2.0; for 208 amino acids, it is log10(n)=-1.75; for 541 amino acids, it is log10(n)=-2.4; finally for 873 amino acids, it is log10(n)=-2.5. Thus, it appears that the longer the genome is, the lower the optimal mutation rate. This result reproduces experimental trend, as shown in Figure 4 where we report experimentally observed (ref. [START_REF] Sanjuan | Viral mutation rates[END_REF] and Table S1) and computed mutation rates allowing maximal stress for RNA virus. Two viruses reported in reference [START_REF] Sanjuan | Viral mutation rates[END_REF] are excluded in Table S1 and Figure 4 because these are bacteriophages and not interact with antibody. Nonetheless, there is an equivalent pressure for bacteriophages that manifests as an evolutionary "arms race model" [START_REF] Labrie | Bacteriophage resistance mechanisms[END_REF][START_REF] Stern | The phage-host arms race: shaping the evolution of microbes[END_REF] .

Thus, we envision that there is also an optimal mutation rate for these systems, but the underlying biophysics (such as protein-DNA interaction) could be different. Note that the longest genome we simulated is 873 amino-acids long, and we did not extend our study to longer genomes because of computational cost to perform simulations with longer genomes.

Influence of the burst size

Larger burst sizes increase the number of viral particles, which effectively lowers the chance of extinction. Thus, we also explored the effect of burst size on the phase diagram (Figure 5 and Figure S4). At higher burst sizes, the onset of lethal mutagenesis is delayed and the viral population is able to tolerate higher concentration of antibody: when a single viral particle gives rise to more progeny, the chances to find a stabilizing mutation increase, and a higher maximal mutation rate for the lethal mutagenesis is achieved. Interestingly, the effect of burst size on the onset of lethal mutagenesis (from 12.0 to 48.8 mutations per genome replication) is greater than its effect on the optimal mutation rate (Figure 5). Indeed, it appears that the burst size (M) has a small influence on the value of this mutation rate (g) -ranging from 4.2 to 8.3 mutations per genome replication.

Discussion and Conclusion

We presented here a simple biophysics-based model of dynamic of viral evolution that quantitatively explores the balance between lethal mutagenesis and the requirement for beneficial mutations that allows escape from antibody. Contrary to previous models of viral evolution that a priori assume a fixed distribution of fitness effects of mutations [START_REF] Bull | Theory of lethal mutagenesis for viruses[END_REF][START_REF] Badgett | Evolutionary dynamics of viral attenuation[END_REF] , we instead base our assumptions on the biophysical genotype-phenotype relationship and distribution of biophysical effects of mutations. As a consequence, the distribution of selection coefficient depends on the status of the system (is it equilibrated or not) and on the environment, reflecting epistasis in the viral fitness landscape. The distribution of biophysical effects can be explored computationally [START_REF] Tokuriki | The Stability Effects of Protein Mutations Appear to be Universally Distributed[END_REF] and experimentally [START_REF] Kumar | ProTherm and ProNIT: thermodynamic databases for proteins and protein-nucleic acid interactions[END_REF] for some proteins. First, using this biophysics-based model of evolution, we find that the survival and viability of the viral population under antibody stress is optimal between 4.2 and 8.3 mutations per genome per replication (for a burst size between 5 and 100). This value can indeed be explained as a balance between lowering the binding affinity with the antibody and keeping the viral proteins stable enough. At low mutation rate, there is only a small likelihood of finding a mutation that could strongly perturb the binding of the capsid protein to the antibody (strongly beneficial mutation). On the other hand, when the mutation rate is too high, some mutations could perturb the binding of the capsid to the antibody while other mutations could compromise the folding stability of any viral protein.

The mutation rate of RNA virus is commonly estimated to be approximately 1 mutation per genome per replication [START_REF] Drake | The distribution of rates of spontaneous mutation over viruses, prokaryotes, and eukaryotes[END_REF][START_REF] Drake | Rates of spontaneous mutation[END_REF] . We report in Table S1 the experimentally observed mutation rates of some RNA viruses [START_REF] Sanjuan | Viral mutation rates[END_REF] . The average mutation rate per genome for these RNA viruses is g=0.37, ranging from 0.024 to 1.2 mutations per genome per replication. Our model of viral evolution is simple and does not account for other biological factors that could influence mutation rates or their fitness effects (for example, cellular quality control response of host cells, the cost of protein production, selection for optimal codon or secondary RNA structures). Thus, while our model predicts a mutation rate that allows a maximal amount of stress, the exact values of these mutation rates may differ from the experimental values of actual viral mutation rates given in Table S1. There is a broad (approximately 2 orders of magnitude) variation in the measured and estimated values of mutation rates (Table S1) [START_REF] Sanjuan | Viral mutation rates[END_REF] , presumably due to confounding factors such as effective population and genome length. Indeed, in agreement with earlier experimental works [START_REF] Sanjuan | Viral mutation rates[END_REF][START_REF] Anderson | Viral error catastrophe by mutagenic nucleosides[END_REF] , we found that the mutation rate per nucleotide that allows a maximal amount of stress decreases when the genome length increases (Figure 4). We also showed the dependence of the mutation rate which allows maximal stress on burst size (Figure 5, see also Figure S4). The dependence of the optimal mutation rate on burst size is due to the increased likelihood at higher burst sizes of finding a variant that acquired a beneficial mutation. In a similar vein, the higher burst size results in an effective increase in fitness because of higher number of progenies.

Under mutation-selection balance, the evolution against the antibody should proceed towards the near-neutral regime [START_REF] Serohijos | Merging molecular mechanism and evolution: theory and computation at the interface of biophysics and evolutionary population genetics[END_REF] . In the review article of Wang & Bull [START_REF] Bull | REVIEW: Optimality models in the age of experimental evolution and genomics[END_REF] , they pointed out that fitness has three components: burst size, adsorption rate at which viruses could interact with cells, and lysis time (time interval between infect and burst). In our model, we account for adsorption rate (which is proportional to the probability of infection) and burst size, but not lysis time. Our main goal here is to immediately connect to biophysics, thus we mapped the fitness to the adsorption rate (Pinf). A possible compensation between lysis time and burst size will require additional assumptions [START_REF] Bull | REVIEW: Optimality models in the age of experimental evolution and genomics[END_REF] , which will be the subject of future works. It must also be pointed out that the current report assumes monoclonal antibodies. A very interesting extension of the study

FiguresFigure 1 .

 1 Figures

  Figure 1. Experimental and computational methodologies to study evolution of viral escape. The

Figure 2 .

 2 Figure 2. Viral evolution under an increasing pressure. In panel (a), the viral count first increases

Figure 3 .

 3 Figure 3. Phase diagram of virus survival under antibody stressor. At each mutation rate, the

Figure 4 .

 4 Figure 4. Influence of the genome length on the viral mutation rate. Computed mutation rates

Figure 5 .

 5 Figure 5. Influence of the burst size on the optimal mutation rate (blue open-triangle) and on the

Acknowledgments

The authors thank Andrew B. Feldman, Jeffrey S. Lin, Assaf Rotem, Thomas J. Smith and Christiane E. Wobus for fruitful discussions. They also thank Amy Gilson for comments on the manuscript. This work was supported by Defense Advanced Research Projects Agency (DARPA) contract HR0011-11-C-0093. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript. The computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.

Conflict of Interest Statement:

The authors declare there is no conflict of interest associated with the publication of this manuscript.

After some time, the population starts migrating towards mutation-selection balance by acquiring deleterious mutations and stays stable for the 100 first passages (Figure S2). Upon addition of antibody (at passage 101), the viral count dramatically drops, but then some particles find beneficial mutations that eventually lead to viral escape from antibody. After around 10 passages, the viral count reaches again a new steady state (Figure S2). When the external antibody pressure is too high, the population goes extinct (see Figure 2). We consider the population extinct when the titer is below the "bottleneck" count, which is the number of particles passaged to the next generation (10 4 particles). For a burst size of M=100, a population goes extinct when less than 100 cells are infected since it will give rise to less than 10 4 released virus particles. Expectedly, as the concentration of antibody is increased, the viral titer keeps dropping and equilibrating. However, the cycle of administration of the antibody affects the population dynamics of the viruses. Progressive increase of the antibody concentration allows the population to adapt before being stressed by another increase of the amount of antibody.

Thus the population can tolerate higher antibody concentration.

The novelty of our model lies in the non-fixed distribution of fitness effects. Thus, the distribution of the selection coefficient (defined as the probability of infection before and after mutation, s= P inf after -P inf before ) at different passages is different. We present in Figure 2-b the distribution of selection coefficient for the last passage with [Ab]=10 3 , the first passage with

[Ab]=10 [START_REF] Zeldovich | Protein stability imposes limits on organism complexity and speed of molecular evolution[END_REF] , and the last passage with [Ab]=10 [START_REF] Zeldovich | Protein stability imposes limits on organism complexity and speed of molecular evolution[END_REF] from Figure 2-a. When the concentration of antibody increases (from purple pentagons to blue triangles), the system is no more equilibrated.

Thus, the fraction of lethal mutations decreases (see insert on the left) and the fraction of beneficial mutations increases (see insert on the right). After 20 passages without changing the environment (from blue triangles to red circles), the system had enough opportunity to adapt and will be to include other critical players: host cells, viruses, and B-cells in the germinal centers [START_REF] Muyoung | Diversity Against Adversity: How Adaptive Immune System Evolves Potent Antibodies[END_REF] .

Indeed, a biophysics-based evolutionary model of the immune response could uncover its physical basis. In the future, other stressors could also be considered in the model presented here (such as inhibitors and pH) because effects of mutations on these quantities can be estimated.

Supplementary Materials

Table S1: Experimentally observed mutation rates and burst size of some RNA viruses.