A biophysical protein folding model accounts for most mutational fitness effects in viruses

Nearly Neutral Thermodynamic Fitness Landscape.

Model genomes comprise Γ genes, each encoding an essential protein completely described by a free energy of folding (ΔG_i). We do not explicitly represent nucleotides. Many real proteins, or their domains, fold two state (12) and fluctuate in thermal equilibrium between a unique native conformation and a set of unfolded decoy conformations. By elementary statistical mechanics, the fraction of time spent in the native state (P^nat) is given by (Eq. 1)

where k_b is Boltzman's constant and T is temperature. We assume that native structures are perfectly functional, whereas all unfolded conformations are completely nonfunctional. Thus, P^nat is proportional to the concentration of functional protein. This paradigm assumes that function requires well-defined structure, which is not true for intrinsically disordered proteins (13). Although viral proteins often have short disordered regions, longer stretches (e.g., disordered domains) exist mostly in eukaryotes (14).

Viability requires that all essential proteins fold and have their native activity. This suggests the AND-like fitness function, (Eq. 2)

where Γ is the number of essential genes and b is the birth rate (i.e., fitness). b_o summarizes the activity of all proteins and equals zero (lethal phenotype) if a mutation abolishes activity of any protein (see below). Truly beneficial mutations would correspond to increases in b_o, but we focus on short evolutionary timescales during which these mutations are unlikely to occur. Thus, without loss of generality, we henceforth scale time such that b_o = 1 for all viable viruses. Eq. 2 states that fitness is reduced by the partial absence of folded protein. In reality, fitness is also compromised by the presence of unfolded proteins, which tend to aggregate and poison the organism. We heuristically capture this effect by also assuming that lethal phenotypes occur whenever P_i^nat < 0.5 (i.e., ΔG_i > 0) for any gene. The resulting fitness function is illustrated in Fig. 1A. Note that we do not assume any tradeoff between fitness and stability. Such tradeoffs often result from mutations to an enzyme's active site (15–17) when evolving a new function (i.e., when adapting, which we do not consider).

We consider two types of nonsynonymous mutations. The first type abolishes activity by introducing STOP codons or disrupting critical residues (e.g., the active site), leading to b_o = 0. We conservatively estimate that, together, these unconditionally lethal mutations comprise 10% of all nonsynonymous mutations (SI Text). The remaining 90% alter thermodynamic stability (and thus, b) by an amount ΔΔG drawn from a Gaussian distribution [p(ΔΔG)] with mean + 1 kcal/mole and SD = 1.7 kcal/mole (18). Our Gaussian form closely approximates p(ΔΔG) obtained computationally (19) as well as obtained from thousands of biophysical measurements taken from the ProTherm database (20) (Fig. 1B). Because (ΔΔG)_mean ∼ 1 kcal/mole, a protein with stability ΔG can tolerate ΔG mutations, on average, before unfolding. Mutations that push ΔG > 0 are considered lethal (see above). Note that relatively stable proteins have more viable single mutant neighbors than do less stable proteins (21). Additionally, we assume that p(ΔΔG) is (i) independent of ΔG and (ii) the same for all proteins. Both i and ii are roughly supported by computational studies (19, 22). Of course, i must fail for some proteins (e.g., the most stable sequence folding to a particular structure), but hyper-stable sequences are sufficiently rare that they are not visited by evolution in simulations or observed in reality (20).

We refer to Eqs. 1 and 2, along with our Gaussian form of p(ΔΔG), as the nearly neutral thermodynamic landscape, which features continuous fitness effects (Fig. 1C) and complex epistatic patterns (Figs. 1 and 2). A natural measure of mutations’ strength is the selection coefficient s ≡ (b_after − b_before)/b_before. Note that s is independent of bystander proteins not involved with the mutation, because those factors of P^nat cancel in Eq. 2. Using Eq. 1, it is easy to show (SI Text) that (3)

which is plotted in Fig. 1C. Expression 3 has two important consequences. First, because ΔΔG is Gaussian, s is log-normally distributed (for a single gene), featuring long tails of deleterious mutations (Fig. S1) when ΔG << 0.

A second important consequence of expression 3 is that the landscape is epistatic, because s depends not only on mutations’ biophysical impact (ΔΔG) but also on protein stability (ΔG) before mutation. We will fully explore epistasis in future work—for now, we merely point out important features. First, although destabilizing mutations are energetically independent, their fitness effects interact synergistically when they occur in the same gene. This negative epistasis for fitness follows from the sigmoidal shape of Eq. 1, which follows from cooperative protein unfolding. Recently, negative epistasis was also observed in yeast regulatory sites using an energy-based approach (23). Second, our fitness landscape appears smooth when defined over the space of protein stability (Fig. 1A), but it is actually rugged when viewed in sequence space because of the random statistical way that we assign mutational effects through p(ΔΔG). This ruggedness is evident in Fig. 2, where colored curves depict random walks through sequence space that terminate on lethal mutations. Compensatory mutations occur frequently (Fig. 2 has additional commentary).

When applied to a single gene, our landscape is verified experimentally with no adjustable parameters. Refs. 24 and 25 introduced mutations by in vitro mutagenesis of an antibiotic resistance enzyme (TEM-1 β-lactamase), which was then expressed in Escherichia coli. Fitness was measured as the fraction of clones surviving antibiotic treatment, averaged over all clones carrying a common number of mutations. Using Eqs. 1 and 2 and Γ = 1, our model quantitatively matches the observed decline in enzyme function (fitness) as the number of mutations (k) increases (Fig. 2). For small k, fitness declines exponentially solely because of unconditionally lethal mutations. For larger k, after thermodynamic stability is exhausted, the decline becomes steeper. An approach by Bloom et al. (24) that approximates protein function by a binary variable (zero or one) also matches the averaged data. However, such landscapes (which we refer to as strictly neutral) preclude both compensatory mutations and continuously varying fitness effects (Fig. 2).

Distribution of Fitness Effects.

The collection of first mutational steps in Fig. 2 encodes the DFE for a particular gene: TEM-1 β-lactamase. Recently, data from ref. 25 were reused to estimate the DFE for this enzyme (26). The resulting DFE was mostly unimodal, with ∼60% of mutations neutral, 8% lethal, and 0% beneficial/compensatory. Although that result is valuable, the question remains as to how mutations impact the fitness of arbitrary proteins and the organism (virus) as a whole. To address this question in the context of our nearly neutral thermodynamic landscape, we must first possess the distribution of stabilities [p(ΔG)] among all of the virus's (essential) proteins. One approach is to use the global, empirical [p(ΔG)] from ProTherm (20) (SI Text). However, in general, we predict that p(ΔG) depends on the evolutionary forces of mutation, selection, and stochastic drift (see below). To this end, we prepared the desired distribution p(ΔG) by simulating populations of stochastically evolving asexual viruses in our near-neutral model (Methods), each containing Γ = 20 essential proteins, under a range of population sizes (N) and realistic mutation rates (m).

After evolving for a long time (50,000 generations), populations reached a steady state of mutation/selection/drift balance independent of initial conditions (Movies S1 and S2). Starting from steady state, we separately mutated genes in each virus and measured the resulting distribution (DFE) of selection coefficients (s). We then averaged these distributions over the population and independent runs. Fig. 3 compares our DFE with an experiment on the ssRNA-based vesicular stomatitis virus (VSV) in which 40 random nonsynonymous mutations were introduced by site-directed mutagenesis (7). Experiments on other viruses yielded qualitatively similar results (10) (Fig. S2). The experimental and model DFE both feature a bimodal shape with a peak near neutrality (s = 0) and another peak at lethality (s = −1). Interestingly, both model and experiment show some beneficial mutations (s > 0), although the population is not being challenged with a new environment. Our model interprets these mutations as stabilizing mutations that compensate for partially destabilized proteins rather than true novel adaptations. Fig. 4 shows that our model can account for intermediate to large deleterious and compensatory fitness effects (e.g., nonlethal deleterious mutations decrease fitness by 2–12% on average). By contrast, on strictly neutral landscapes (18, 24, 27), all mutations must be either lethal or neutral. Because the DFE for a single gene is log-normal, the overall model DFE is a linear combination of log-normals, with weights given by p(ΔG). This helps explain why the experimental DFE is easily fit by log-normal distributions (10, 11).

A striking feature of Figs. 3 and 4 is that the DFE depends on population size (N) and mutation rate (m). In other words, the DFE depends not only on molecular constraints and the fitness landscape implied by them but also evolutionary forces that pull the population to the flat region of Fig. 1A or push it to the cliff. In particular, we find that mutations are more potent (i.e., genomes are less robust) in small populations with high mutation rates. We now analyze the origins of these effects by considering how a population's location on the fitness landscape, measured by p(ΔG), depends on N and m.

Small Population Size Decreases Evolved Thermodynamic Stability.

Fig. 5 illustrates the distribution of protein stability p(ΔG) for several values of N and m. Each instance of p(ΔG) vanishes at ΔG = 0, peaks when −ΔG equals a few kilocalories per mole, and decays for more negative ΔG. These qualitative features agree with previous work (analytic and simulation) (18, 27, 28) as well as the empirical distribution sampled from different genes and organisms in all kingdoms of life obtained from the ProTherm database (20). For appropriately chosen N and m, p(ΔG) from simulations agrees quantitatively with the global distribution from ProTherm (22) (Figs. S3 and S4). We note that, because fitness increases monotonically (although exponentially weakly) with −ΔG in our model (Eq. 1), perfectly fit proteins would have ΔG = −∞, which is obviously not observed. Therefore, the distribution's peak near ΔG ∼ −2 to −5 kcal/mole does not reflect optimality of this stability range but rather, a balance between mutation, selection, and random drift. Specifically, deleterious mutations with |s| < 1/N have a significant probability of taking over the population (by random drift), making small populations especially poor optimizers (4). Fig. 5 shows that, in mutation/selection/drift steady state, partially destabilized (i.e., less fit) proteins are more heavily represented in small populations. Mutations that occur within these destabilized proteins have large fitness effects (Fig. 1C) that disperse the DFE (Figs. 3 and 4).

We can semiquantitatively understand how p(ΔG) depends on N. In an ideal (mostly monoclonal) population, destabilizing mutations will cease to fix when Ns ∼ −Ne ^ΔG/kT ∼ 1 (4) (i.e., when ΔG ∼ −k_bT ln N). This expression matches our data reasonably well for Nm < 1 (Fig. 5B). For larger Nm values (polyclonal populations), N should be replaced by some effective population size N_e (which is less than N), reflecting the additional stochasticity conferred by “background selection” (29). Evidently, populations as large as N = 10⁵ cannot be treated deterministically (i.e., N→ ∞), because ΔG continues to decrease with N (Fig. 5B). We also note that deep population bottlenecks (Methods) abet protein destabilization in our simulations (Fig. 5A), and therefore, this effect is not specific to constant N population genetic models.

High Mutation Rates Decrease Evolved Thermodynamic Stability.

Deleterious mutations, most of which are not completely lethal, are incessantly produced with rate ∼ m. Meanwhile, before being purged from the population, deleterious mutations linger for ∼b_avg/(b_avg − b) ∼ 1/s generations, independent of m. Therefore, as m increases, more deleterious mutations, in the form of slightly destabilized proteins, accumulate in the population. This reduction in fitness, known as the mutation load, (30) is evident in Fig. 5A and can be quite strong (e.g., ∼3.5 kcal/mole for N = 10⁵) (Fig. S5). Mutations that arise in the background of excess destabilized proteins have large fitness effects, thereby dispersing the DFE (Figs. 3 and 4) and decreasing the robustness of genomes to mutations.

Although this phenomenon is grounded in simple population genetics principles, it clashes with several existing models of neutral networks (18, 21, 31, 32) in which a mutationally connected set of genomes is perfectly fit, whereas another set is completely lethal. Those studies predict that, at large m, sequences connected to many viable mutational neighbors will predominate over less connected sequences. In the language of our model, the prediction from neutral network theory is that greater stability will occur at large m, which contradicts Fig. 5A. In SI Text, we show that our model reproduces the neutral network prediction if we artificially set temperature (T) to zero (i.e., if we eliminate all curvature from our fitness landscape). This point is very important in understanding the general relationship between mutation rate and mutational robustness (33, 34).