Article Information

• PDF (261 kb)
• Full Text with Large Figures
• Cited by in Scopus (9)

PubMed

• Articles by Claus O. Wilke
• Articles by Alpan Raval

Related Articles

• Kinetics of the Micelle-to-Vesicle Transition: Aqueous Lecit...

  Kinetics of the Micelle-to-Vesicle Transition: Aqueous Lecithin-Bile Salt Mixtures Biophysical Journal, Volume 85, Issue 3, 1 September 2003, Pages 1624-1646 J. Leng, S.U. Egelhaaf and M.E. Cates Abstract Important routes to lipid vesicles (liposomes) are detergent removal techniques, such as dialysis or dilution. Although they are widely applied, there has been only limited understanding about the structural evolution during the formation of vesicles and the parameters that determine their properties. We use time-resolved static and dynamic light scattering to study vesicle formation in aqueous lecithin-bile salt mixtures. The kinetic rates and vesicle sizes are found to strongly depend on total amphiphile concentration and, even more pronounced, on ionic strength. The observed trends contradict equilibrium calculations, but are in agreement with a kinetic model that we present. This model identifies the key kinetic steps during vesicle formation: rapid formation of disklike intermediate micelles, growth of these metastable micelles, and their closure to form vesicles once line tension dominates bending energy. A comparison of the rates of growth and closure provides a kinetic criterion for the critical size at which disks close and thus for the vesicle size. The model suggests that liposomes are nonequilibrium, kinetically trapped structures of very long lifetime. Their properties are hence controlled by kinetics rather than thermodynamics. Abstract | Full Text | PDF (398 kb)

• Dynamics of Cellular Focal Adhesions on Deformable Substrate...

  Dynamics of Cellular Focal Adhesions on Deformable Substrates: Consequences for Cell Force Microscopy Biophysical Journal, Volume 95, Issue 2, 15 July 2008, Pages 527-539 Alice Nicolas, Achim Besser and Samuel A. Safran Abstract Cell focal adhesions are micrometer-sized aggregates of proteins that anchor the cell to the extracellular matrix. Within the cell, these adhesions are connected to the contractile, actin cytoskeleton; this allows the adhesions to transmit forces to the surrounding matrix and makes the adhesion assembly sensitive to the rigidity of their environment. In this article, we predict the dynamics of focal adhesions as a function of the rigidity of the substrate. We generalize previous theories and include the fact that the dynamics of proteins that adsorb to adhesions are also driven by their coupling to cell contractility and the deformation of the matrix. We predict that adhesions reach a finite size that is proportional to the elastic compliance of the substrate, on a timescale that also scales with the compliance: focal adhesions quickly reach a relatively small, steady-state size on soft materials. However, their apparent sliding is not sensitive to the rigidity of the substrate. We also suggest some experimental probes of these ideas and discuss the nature of information that can be extracted from cell force microscopy on deformable substrates. Abstract | Full Text | PDF (424 kb)

• Flexible Charged Macromolecules on Mixed Fluid Lipid Membran...

  Flexible Charged Macromolecules on Mixed Fluid Lipid Membranes: Theory and Monte Carlo Simulations Biophysical Journal, Volume 89, Issue 5, 1 November 2005, Pages 2972-2987 Shelly Tzlil and Avinoam Ben-Shaul Abstract Fluid membranes containing charged lipids enhance binding of oppositely charged proteins by mobilizing these lipids into the interaction zone, overcoming the concomitant entropic losses due to lipid segregation and lower conformational freedom upon macromolecule adsorption. We study this energetic-entropic interplay using Monte Carlo simulations and theory. Our model system consists of a flexible cationic polyelectrolyte, interacting, via Debye-Hückel and short-ranged repulsive potentials, with membranes containing neutral lipids, 1% tetravalent, and 10% (or 1%) monovalent anionic lipids. Adsorption onto a fluid membrane is invariably stronger than to an equally charged frozen or uniform membrane. Although monovalent lipids may suffice for binding rigid macromolecules, polyvalent counter-lipids (e.g., phosphatidylinositol 4,5 bisphosphate), whose entropy loss upon localization is negligible, are crucial for binding flexible macromolecules, which lose conformational entropy upon adsorption. Extending Rosenbluth’s Monte Carlo scheme we directly simulate polymer adsorption on fluid membranes. Yet, we argue that similar information could be derived from a biased superposition of quenched membrane simulations. Using a simple cell model we account for surface concentration effects, and show that the average adsorption probabilities on annealed and quenched membranes coincide at vanishing surface concentrations. We discuss the relevance of our model to the electrostatic-switch mechanism of, e.g., the myristoylated alanine-rich C kinase substrate protein. Abstract | Full Text | PDF (427 kb)

• …more

Copyright © 2005 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 89, Issue 6, 3714-3720, 1 December 2005

doi:10.1529/biophysj.105.062125

Biophysical Theory and Modeling

Predicting the Tolerance of Proteins to Random Amino Acid Substitution

Claus O. Wilke^*, §, ,  , Jesse D. Bloom^†, §, D. Allan Drummond^‡, § and Alpan Raval^*, ¶

^* Keck Graduate Institute of Applied Life Sciences, Claremont, California
^† Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California
^‡ Program in Computation and Neural Systems, California Institute of Technology, Pasadena, California
^§ Digital Life Laboratory, California Institute of Technology, Pasadena, California
^¶ School of Mathematical Sciences, Claremont Graduate University, Claremont, California

Address reprint requests to C. O. Wilke at his present address, Section of Integrative Biology, University of Texas at Austin, Texas.

A protein’s tolerance to random amino acid substitutions is of fundamental importance both in protein engineering and molecular evolution. In molecular evolution, a protein’s neutrality, that is, the fraction of single amino acid substitutions that do not disrupt the protein’s function, has a substantial influence on how this protein evolves and accumulates mutations 1,2,3,4,5,6. In protein engineering, the knowledge of a protein’s tolerance to mutations helps one to optimize the mutagenesis conditions in directed protein evolution 7; several groups have characterized experimentally a protein’s loss of function under random mutations 8,9,10,11.

Protein mutagenesis studies suggest that a large fraction of deleterious amino acid substitutions disrupt a protein’s structure rather than specifically affecting functional residues 12,13,14. Therefore, the fraction of substitutions that disrupt a protein’s structure is a reasonable lower bound to the fraction of substitutions that will disrupt a protein’s function.

We 8 have recently proposed a thermodynamic model that allows one to calculate the probability P_f(m) with which a protein retains its structure after m amino acid substitutions. This model uses as input the distribution of free energy changes ΔΔG for individual amino acid substitutions. It is based on the idea that the free energy change caused by one amino acid substitution is independent of the change caused by another such substitution, and that the protein continues to fold correctly as long as its free energy of folding remains below some threshold level. If the protein’s free energy of folding is initially a distance C from the threshold, then the fraction of sequences with m substitutions that still fold correctly is given by the fraction of sums that are less than C, where the X_i are independent, identically distributed random variables taken from the ΔΔG distribution. For a small set of both simulated lattice proteins and real proteins, we 8 have shown that this model has excellent predictive power. Here, we are interested in three questions:

We implemented a maximally compact, 5×5 two-dimensional square lattice model, as previously described 15,5. In short, we folded simulated polypeptide chains of length L=25 residues into a maximally compact structure, representing one of the 1081 possible 16 self-avoiding compact walks of length 25 not related by rotational or reflection symmetry. (We neglected the vanishingly small fraction of palindromic sequences.) We used an alphabet of 20 amino acids, and calculated the contact energies between nonbonded neighboring residues according to Table 3 of Miyazawa and Jernigan 17. We calculated a lattice protein’s free energy of folding ΔG_f as described by Taverna and Goldstein 15, and considered the protein to be stably folded if ΔG_f was below a cutoff ΔG_cut. We carried out all analyses for three different cutoffs, ΔG_cut=−4.0kcal/mol, −5.0kcal/mol, and −6.0kcal/mol.

We first analyzed a dataset of 300 randomly chosen sequences, 100 at each cutoff. We generated these sequences in the following way: First, we generated random sequences and tried to fold them. We kept all those sequences whose free energy of folding was below ΔG_cut=−4.0kcal/mol, and whose native conformation was different from the native conformations of all stably folding sequences we had encountered so far. We repeated this procedure until we had 100 sequences that could stably fold into 100 unique conformations at ΔG_cut=−4.0kcal/mol. For the remaining two cutoffs, we used hill climbing and subsequent neutral evolution to obtain, at each cutoff, 100 additional sequences that could stably fold into the same 100 conformations as the original sequences. Under hill climbing, we repeatedly mutated a sequence, and accepted all mutations that increased the protein’s stability without changing the native conformation. Under neutral evolution, we repeatedly mutated a sequence, and accepted all mutations that did not destabilize the protein beyond the chosen cutoff and did not change the native conformation. We always repeated neutral evolution until we had accepted 1000 mutations.

For all 300 sequences, we estimated P_f(m), the fraction of mutant proteins that fold stably to the original native conformation after m amino-acid substitutions, by randomly sampling mutants according to the following procedure: We carried out all single-point mutations, and sampled 10⁴, 5×10⁴,10⁵,…10⁷ multiple-point mutations for m=2,3,4,…,8. We then calculated P_f(m) by dividing the number of correctly folded sequences that we found at the given mutational distance m by the total number of mutants we tried at that distance. We defined a protein as correctly folded if its minimum free energy was below the chosen cutoff ΔG_cut and if its native conformation was identical to that of the starting sequence. In the vast majority of these 300 replicates, we found between several hundred and several thousand correctly folded proteins at each mutational distance m. Consequently, our estimate for P_f(m) in lattice proteins is highly accurate.

We measured the ΔΔG distribution of each of the 300 sequences by carrying out all possible single-point mutations, and then calculating the differences between the minimum free energy of the original sequence and the mutated sequences.

We calculated the prediction for P_f(m) from the ΔΔG distribution as described 8. In short, we first binned the ΔΔG distribution into bins of width 0.01kcal/mol, and then calculated the m-fold convolution of this binned distribution using the fast Fourier transform of the software package R, version 1.9.1 18. Finally, we numerically integrated the convolved distribution from −∞ to C to obtain P_f(m).

We carried out a second set of simulations to determine the influence of the starting sequence on the neutrality 〈ν〉. We selected the sequences of 10 representative conformations (among the 100 unique conformations of the first data set), and generated, through neutral evolution as before, for each conformation at each cutoff nine additional sequences folding stably into this conformation. We measured then both P_f(m) and the ΔΔG distribution for these additional 270 sequences as described above.

P_f(m) decays approximately as 〈ν〉^m for large m. We estimated 〈ν〉 from the measured P_f(m) by carrying out a linear regression of ln P_f(m) versus m, where we restricted the range of m from 4 to 8 to capture the asymptotic behavior of P_f(m). The neutrality 〈ν〉 followed then as 〈ν〉=e^a, where a is the slope of the regression line.

We also calculated 〈ν〉 in the context of a number of approximation schemes, described in Appendix A Edgeworth expansion,Appendix B Cramér approximation,Appendix C Markov chain approximation,Appendix D Mean-field approximation, and summarized in Results, below. For the Cramér approximation (Appendix B ), we numerically minimized the moment-generating function ϕ(t) of the ΔΔG distribution. Let {ΔΔG_i} be the set of free energy changes caused by all single point mutations. Then, and its derivative We numerically found the value t^* at which ϕ′(t^*)=0, and then set 〈ν〉=ϕ(t^*).

For the Markov chain approximation (Appendix C ), we constructed the matrix W_ij using bins of width 0.015kcal/mol, and spanning a range of 25.0kcal/mol, from ΔG_cut to ΔG_cut −25.0kcal/mol. We calculated the largest eigenvalue of this matrix by repeatedly multiplying W_ij to a vector (with all components initially set to one), and then renormalizing the vector to unit length, until the vector had converged to the dominant eigenvector of W_ij. We then obtained the quantity 〈ν〉 from the change in length in the dominant eigenvector of W_ij after a single multiplication with W_ij.

First, we assess how well our method to predict P_f(m) works in a large data set. We 8 have previously studied only a handful of noncompact lattice proteins and three real proteins. Overall, we find that the method works very well for the compact lattice proteins we study here. Fig. 1 shows several typical examples. In many cases, we find that the prediction of P_f(m) is highly accurate up to m=8, which is the largest number of mutations we consider (Figure 1AD). In those cases that show some discrepancy between the predicted and the measured P_f(m), we typically find that the prediction works well up to m=3 or 4, but starts to deviate from the measured results for larger m. There is no clear tendency toward either over- or underestimation of the measured results by the prediction (Figure 1EH). Note that our data set covers a wide range of different conformations, as all 100 sequences at a given cutoff fold into a unique conformation.

Display large version of this figure

Figure 1

Fraction of correctly folded mutants P_f as a function of the number of mutations m, for eight lattice proteins that stably fold at ΔG_cut=−5.0kcal/mol. Points indicate measurement results, and solid lines indicate the prediction derived from the m-fold convolution of the ΔΔG distribution. A–D show cases for which the prediction works excellently, and E–H show cases for which there is some disagreement between the prediction and the measured values. The RMS values for these eight cases are (left to right and top to bottom): 0.127, 0.112, 0.157, 0.065, 1.423, 2.407, 1.406, and 0.902.

We can quantify the performance of our prediction using the root-mean-squared (RMS) deviation of the log-transformed P_f(m). Let be the predicted fraction of mutants that fold correctly, and P_f(m) the corresponding measured value. Then, we define the logarithmic RMS deviation ρ as

(1)

Display large version of this figure

Figure 2

Histogram of RMS values for 100 randomly chosen sequences each, at three cutoff levels.

Next, we are interested in asymptotic expressions of P_f(m) for small and large m. For small m, we can approximate P_f(m) using the Edgeworth expansion (Appendix A ). The Edgeworth expansion provides correction terms to the central limit theorem for finite sums of random variables. These correction terms take into account successively higher moments of the ΔΔG distribution. Fig. 3 shows how the Edgeworth expansion provides an increasingly more accurate approximation of P_f(m) as higher-order terms are included. However, whereas in some cases the Edgeworth expansion works very well with only three additional moments beyond mean and variance (Figure 3A), in other cases the Edgeworth expansion deviates significantly from P_f(m) in all orders we have considered (Figure 3B). Furthermore, because the Edgeworth expansion leads to a normal distribution function multiplied by a polynomial (Eq. (3)), it must inevitably break down as m becomes large.

Display large version of this figure

Figure 3

Prediction of P_f(m) according to the Edgeworth expansion. Points indicate the measured P_f(m), and lines indicate the Edgeworth expansion to various orders (dotted lines, normal term only; dot-dashed lines, normal term plus first-order corrections; dashed lines, normal term plus first- and second-order corrections; and solid lines, normal term plus first-, second-, and third-order corrections). (A) Example of a case where the expansion works well up to m=8. (B) Example of a case where the expansion works poorly.

For large m, empirical observations show that P_f(m) decays approximately as 〈ν〉^m (8,9,10 and Fig. 1). The value 〈ν〉 can vary substantially among sequences, but generally tends to increase with the cutoff (Fig. 4). We can interpret 〈ν〉 intuitively as the average neutrality of all sequences that stably fold into the given structure. We give a formal argument for this interpretation in Appendix C . An exponential decay of the form P_f(m) ≈ 〈ν〉^m follows from the Gaussian term in the Edgeworth expansion (Appendix A ). However, the value of 〈ν〉 predicted by this term is not very accurate (data not shown). The Gaussian approximation fails because, for large m, P_f(m) is extremely sensitive to small deviations from normality in the tail of the m-fold convolved ΔΔG distribution.

Display large version of this figure

Figure 4

Asymptotic neutralities 〈ν〉 (measured) at the three different cutoffs, sorted by magnitude and plotted against their rank.

Numerically, we can estimate 〈ν〉 by first calculating the prediction for P_f(m) using the m-fold convolution of the ΔΔG distribution, and then obtaining 〈ν〉 from a log-linear regression in the same way in which we estimate it from the measured P_f(m) (see Calculation of 〈ν〉, above). In the following, we refer to this method as the convolution method. The convolution method does not generate any new insight into what determines the value of 〈ν〉, but it serves as a useful test case. First, by comparing for a large set of proteins the measured 〈ν〉 to the 〈ν〉 predicted by the convolution method, we obtain an overall estimate of how well our model performs. Second, the convolution method is the correct benchmark for all other methods of estimating 〈ν〉: Because any deviation between the prediction from the convolution method and the measured 〈ν〉 is an inherent shortcoming of our model, we can only expect that any approximate method to estimate 〈ν〉 will work at most as well as the convolution method, and will generally perform worse. Figure 5A shows that the 〈ν〉 predicted by the convolution method correlates strongly with the measured (overall R² for all 300 data points R²=0.789, p<10⁻¹⁵), in agreement with our earlier observation that, overall, our model works very well.

Display large version of this figure

Figure 5

Prediction for 〈ν〉 according to three different methods. The solid lines indicate perfect agreement between predicted and measured 〈ν〉. Symbols indicate ΔG_cut=−4.0kcal/mol (solid circles), ΔG_cut=−5.0kcal/mol (open squares), and ΔG_cut=−6.0kcal/mol (shaded triangles). (A) Convolution method. Correlations with measured 〈ν〉, in order of increasing absolute cutoff value: R²=0.642, R²=0.540, and R²=0.757. Overall correlation, R²=0.789 (all p<10⁻¹⁵). (B) Cramér approximation. Correlations with measured 〈ν〉, in order of increasing absolute cutoff value: R²=0.584, R²=0.467, and R²=0.747. Overall correlation, R²=0.751 (all p<10⁻¹⁴). (C) Mean-field approximation. Correlations with measured 〈ν〉, in order of increasing absolute cutoff value: R²=0.590, R²=0.410, and R²=0.650. Overall correlation, R²=0.716 (all p<10⁻¹²).

A straightforward method to predict 〈ν〉 from the ΔΔG distribution follows from large-deviation probability theory. Cramér’s theorem implies that P_f(m) must decay exponentially, and implies that 〈ν〉 is approximately given by the unique minimum of the moment-generating function of the ΔΔG distribution (Appendix B ). In Figure 5B, we compare the 〈ν〉 predicted by the Cramér approximation to the measured 〈ν〉. We see that the Cramér approximation performs almost as well as the convolution method. The correlation between the 〈ν〉 values predicted according to the convolution method and the Cramér approximation is very strong (overall R² for all 300 data points R²=0.971, p<10⁻¹⁵).

The intuitive explanation for why P_f(m) decays approximately as 〈ν〉^m is that each correctly folded sequence has, on average, a fraction 〈ν〉 of correctly folded single-point neighbors, so that with each mutational step the total P_f(m) is reduced by a factor of 〈ν〉. We can make this reasoning more precise with the Markov chain approximation. The Markov chain approximation is based on the assumption that single-point mutants to sequences at distance m that do not fold correctly do not contribute to P_f(m+1). With this assumption, 〈ν〉 turns out to be the largest eigenvalue of a matrix W_ij that contains the transition probabilities from any stable protein to any other stable protein under single-point mutations (Appendix C ). We do not present results from the Markov chain approximation in Fig. 5, because they are very similar to those found with the Cramér approximation (overall R² for all 300 data points R²=0.9992, p<10⁻¹⁵). However, the 〈ν〉 values predicted by the Markov chain approximation tend to be slightly smaller than those predicted by the Cramér approximation, the reason being that the Markov chain approximation neglects mutations that stabilize previously unstable sequences (Appendix C ).

The last method we consider is the mean-field approximation. The mean-field approximation is based on the idea that we can replace the distribution of proteins with different neutralities by a single protein with an effective neutrality that equals 〈ν〉, and is extremely simple to calculate (Appendix D ). Figure 5C shows that the mean-field approximation performs only slightly worse than the Cramér approximation. The correlation between the 〈ν〉 values predicted from the convolution method and the mean-field approximation is also strong (overall R² for all 300 data points R²=0.939, p<10⁻¹⁵).

Finally, we have generated an additional data set of 10×10 sequences that fold into the same structure, to assess to what extent 〈ν〉 depends on the initial sequence or the structure. We find that although there is some spread in the estimated 〈ν〉 for different sequences folded into the same structure, the 〈ν〉 values for the different starting sequences clearly cluster around a mean value that is determined by the structure. Fig. 6 shows data for a representative five of the 10 structures we considered for this additional data set. We carried out a pairwise t-test for all 45 possible pairings of the 10 structures, at each cutoff, and found that (after applying the false-discovery-rate correction for multiple testing 19) only 12, 9, and 5 of the 45 pairs at cutoffs ΔG_cut=−4.0kcal/mol, −5.0kcal/mol, and −6.0kcal/mol do not have a statistically significant (at the 5% level) difference in

Display large version of this figure

Figure 6

Asymptotic neutralities 〈ν〉 (measured) for five different structures. Symbols indicate ΔG_cut=– 4.0kcal/mol (solid circles), ΔG_cut=−5.0kcal/mol (open squares), and ΔG_cut=−6.0 kcal/mol (shaded triangles).

We have extensively tested a model introduced earlier to describe and explain the tolerance of proteins to amino-acid substitutions 8. These tests were performed on an array of 100 structures and three cutoff levels. The model performs well across this data set, which gives strong support for the model’s central claims, its generality, and its theoretical underpinnings. The predicted emergence of an exponential decline in the P_f(m) that is parameterized by the mean neutrality 〈ν〉 is both observed and estimated by several independent methods, and the preliminary finding that 〈ν〉 is principally a structural property receives computational support through tests across 10 structures. Using a Markov chain method, we also explain why the rate of the asymptotic decay of P_f(m), as measured by 〈ν〉, is in fact related to the average neutrality of all sequences that can stably fold into the native conformation.

For computational efficiency, we have used maximally compact two-dimensional lattice proteins (with the full amino-acid alphabet). Compact lattice proteins have the drawback that the additional constraint of maximal compactness allows many more sequences to stably fold than otherwise would; also, noncompact lattice proteins rarely fold into maximally compact formations 20,21. However, in previous work 8, we had tested the model against a small set of two-dimensional noncompact lattice proteins, as well as two real proteins, and found the model to perform well in these cases. It therefore seems unlikely that the results that we report here are artifacts of the additional constraint of maximal compactness. Likewise, three-dimensional lattice proteins have substantially more conformations at the same sequence length than two-dimensional lattice proteins, and our model could, in principle, break down in three dimensions. We have no specific reason to believe that our model would perform substantially worse for three-dimensional lattice proteins than for two-dimensional lattice proteins, but this hypothesis remains to be tested.

A key advantage of our model is its extreme simplicity. Our finding that 〈ν〉 can be trivially computed with reasonable accuracy using either a mean-field approximation or a generating function approach that extends the model’s utility. Our finding that the Gaussian term in the Edgeworth expansion cannot accurately describe the data suggests that a Gaussian approximation for the initial ΔΔG distribution is simply not adequate for the estimation of 〈ν〉. Thus our model, although simple, is sensitive to the detailed form of the ΔΔG distribution, rather than just its mean and variance.

Whether these results extend to an equally broad class of naturally occurring proteins remains an open question. A useful feature of our model is that it depends, in a direct and relatively simple manner, on the distribution of the ΔΔG values, which are routinely measured in natural proteins and can be computationally estimated from crystal structures. In general, we do not know the difference C between the native stability of proteins and their minimum free energy cutoff. However, the existence of a cutoff is indicated by diverse observations such as the abundance of temperature-sensitive mutations and the steep (exponential) dependence on stability of the folded and unfolded protein concentrations at equilibrium. We do not know whether the cutoff is consistent across proteins or varies, like 〈ν〉, from structure to structure.

An important practical implication of our model is that the fraction of mutant proteins retaining fold can be increased in a predictable fashion by modest increases in wild-type protein stability. Mutagenesis experiments aimed at discovering functionally improved proteins may thus have stability-dependent optimal mutation rates 7 which, at least in principle, may be estimated using our model. Our results here offer strong support to the suggestion 8 that stability is a critical, but generally overlooked, parameter in directed evolution.

Appendix A

Edgeworth expansion

We wish to estimate the probability where X_i are independent, identically distributed random variables distributed according to the ΔΔG distribution, and C is the distance to the free-energy cutoff beyond which the protein does not stably fold. It is convenient to introduce the standardized random variable where and μ and σ are the mean and standard deviation of the ΔΔG distribution, respectively. Let κ_n be the n^th cumulant (see Appendix E ) of the ΔΔG distribution. We define and write the standard normal distribution function

(2)

Then, the Edgeworth expansion of the distribution function F(z) of the random variable Z is given by 22

(3)

An efficient algorithm to generate higher-order terms of the expansion has been presented by Blinnikov and Moessner 23. From F(z), we obtain P_f(m) via

(4)

The zero^th-order term of the Edgeworth expansion (the Gaussian term, which takes into account only the mean and variance of the ΔΔG distribution) is

(5)

where erfc(t) is the complementary error function. This expression predicts for large m that P_f(m) ≈ 〈ν〉^m with 〈ν〉=exp [−μ²/(2σ²)].

Appendix B

Cramér approximation

We can calculate the asymptotic behavior of P_f(m) for large m from large-deviation theory. According to the central limit theorem, for large m the sum (as introduced in Appendix A ) is approximately normally distributed with mean mμ and variance mσ², where μ and σ² are the mean and variance of the ΔΔG distribution. The probability is therefore a tail probability that becomes vanishingly small as m approaches infinity. Cramér’s theorem 24 for large deviation probability states that, for a<μ,

(6)

where ϕ(t) is the moment-generating function of the distribution of X_i, and t^* is the value of t at which ϕ(t) – at attains its minimum.

Cramér’s theorem can be used as a basis for approximating the asymptotic behavior of Prob(S_m/m≤a), namely, for large m,

(7)

The theorem therefore gives a theoretical justification for the exponential decay of tail probabilities. For the case of interest, we have a=C/m, which is small when m is large. For large m, we may therefore consider an approximation to P_f(m) of the form

(8)

and therefore estimate the average neutrality as 〈ν〉 ≈ ϕ(t^*).

Further refinements to Cramér’s theorem, especially in the context of placing bounds on tail probabilities for finite m, have been the subject of recent advances in large deviation probability theory (see, for example, Hahn and Klass 25 and references therein) and may be used to obtain more accurate estimates. For our purposes, Cramér’s theorem gives a simple and reasonably accurate estimate of P_f(m).

Appendix C

Markov chain approximation

An alternative method to estimate the asymptotic slope 〈ν〉 of P_f(m) is based on calculating the steady-state solution of a suitable Markov process. First, we subdivide the range of free energies of folding into discrete bins of width b. We number the bins consecutively and in such a way that all bins with index i≥0 represent stable proteins, and all other bins represent unstable proteins. Now, let p_i(m) be the fraction of proteins at mutation distance m in bin i. Clearly, we have Next, we introduce the matrix M_ij, which gives the probability that a single mutation to a protein in bin j moves that protein into bin i. (Note that under the assumptions of our theory, M_ij does not depend on m, and furthermore depends only on the difference i – j, but not on the specific values of i or j. The first assumption is necessary for the development of the Markov approximation; the second assumption of stationarity of the transition matrix could be, in principle, relaxed.) Then, we can write P_f(m+1) as

(9)

Our goal is to express the right-hand side in terms of P_f(m), so that we obtain a recursion relation for P_f(m). Unfortunately, the second sum on the right-hand side spans all values of j, positive as well as negative, whereas P_f(m) contains only information about p_j(m) with positive index j. Therefore, we now make the approximation that

(10)

This approximation is based on the assumption that mutations that stabilize an unstable protein are rare in comparison to mutations that do not destabilize a stable protein, and is the main difference between the Markov chain approximation and the m-fold convolution of the ΔΔG distribution. From here on, we will refer to the submatrix of M_ij for which the indices i and j are non-negative as W_ij. The distinction between M_ij and W_ij will become important when we discuss eigenvectors and eigenvalues of W_ij.

We can interpret as the neutrality of a protein in bin j (for j≥0), and the average neutrality of all proteins at distance m is given by

(11)

Since we can replace with P_f(m), we find that P_f(m+1) and P_f(m) are related to each other via

(12)

and, assuming that 〈ν(m)〉 approaches a limiting value 〈ν〉 for large m, we have

(13)

This equation implies that for large m, P_f(m) is proportional to 〈ν〉^m.

From we see that for large m, the p_i are proportional to the dominant eigenvector of W_ij, by virtue of the Frobenius-Perron theorem 26. (The Frobenius-Perron theorem holds if W_ij is primitive—the case whenever there is a path of mutations that leads from any bin i to any other bin j, and W_ii>0 for at least one i.) Furthermore, Eq. (11) implies

(14)

where p_i* is the dominant eigenvector of W_ij. Consequently, 〈ν〉 corresponds to the dominant eigenvalue of W_ij.

Appendix D

Mean-field approximation

A third method to calculate 〈ν〉 is the mean-field approximation. The idea of this approximation is that we can replace the distribution of proteins of different stabilities with a single protein of typical stability. The neutrality of this protein should correspond to the average neutrality of all stable proteins. We choose the stability of this protein such that its free energy of folding is identical to the average free energy of folding of all possible single-point mutants that fold correctly. In other words, the average change in free energy of a single mutation that does not destroy the protein’s ability to fold is zero. The neutrality of this protein is then the fraction of mutations that cause a change in free energy below a certain cutoff, where the cutoff is chosen such that the average change in free energy for all mutations below the cutoff is as close as possible to zero. We can formalize this condition as follows. Assume that the set {ΔΔG_i} contains the free-energy changes caused by all possible single-point mutations (of which there are n), and that the set is ordered such that ΔΔG_i<ΔΔG_i+1 for all i. Then, we have

(15)

Appendix E

Unbiased estimators of cumulants

Let be a set of n measurements, and define

(16)

According to Dressel 27, the following are unbiased estimators for the first five cumulants κ₁–κ₅ (note that κ₁ is the sample average, and κ₂ is the sample variance):

(17)

(18)

(19)

(20)

(21)