| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||



* Keck Graduate Institute of Applied Life Sciences, Claremont, California;
Division of Chemistry and Chemical Engineering,
Program in Computation and Neural Systems,
Digital Life Laboratory, California Institute of Technology, Pasadena, California; and ¶ School of Mathematical Sciences, Claremont Graduate University, Claremont, California
Correspondence: Address reprint requests to C. O. Wilke at his present address, Section of Integrative Biology, University of Texas at Austin, Texas. E-mail: cwilke{at}mail.utexas.edu.
| ABSTRACT |
|---|
| INTRODUCTION |
|---|
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
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 Pf(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 Xi 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:

G distribution? 
G values? | METHODS |
|---|
Gf as described by Taverna and Goldstein (15
Gf was below a cutoff
Gcut. We carried out all analyses for three different cutoffs,
Gcut = 4.0 kcal/mol, 5.0 kcal/mol, and 6.0 kcal/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
Gcut = 4.0 kcal/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
Gcut = 4.0 kcal/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 Pf(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 104, 5 x 104,105,...107 multiple-point mutations for m = 2,3,4,...,8. We then calculated Pf(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
Gcut 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 Pf(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 Pf(m) from the 
G distribution as described (8
). In short, we first binned the 
G distribution into bins of width 0.01 kcal/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 Pf(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 Pf(m) and the 
G distribution for these additional 270 sequences as described above.
Calculation of 


Pf(m) decays approximately as 

m for large m. We estimated 

from the measured Pf(m) by carrying out a linear regression of ln Pf(m) versus m, where we restricted the range of m from 4 to 8 to capture the asymptotic behavior of Pf(m). The neutrality 

followed then as 

= ea, where a is the slope of the regression line.
We also calculated 

in the context of a number of approximation schemes, described in Appendices AD, 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 {
Gi} 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 Wij using bins of width 0.015 kcal/mol, and spanning a range of 25.0 kcal/mol, from
Gcut to
Gcut 25.0 kcal/mol. We calculated the largest eigenvalue of this matrix by repeatedly multiplying Wij 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 Wij. We then obtained the quantity 

from the change in length in the dominant eigenvector of Wij after a single multiplication with Wij.
| RESULTS |
|---|
|
be the predicted fraction of mutants that fold correctly, and Pf(m) the corresponding measured value. Then, we define the logarithmic RMS deviation
as
![]() | (1) |
1.0 or higher (Fig. 1). When we consider all 300 replicates, we find that the majority of the cases have an RMS value below 1.0, and only rarely does the RMS value exceed 2.0 (Fig. 2). There seems to be a slight tendency for the RMS value to increase as the cutoff value becomes more stringent (i.e., from
Gcut = 4 kcal/mol to 6 kcal/mol).
|

G distribution. Fig. 3 shows how the Edgeworth expansion provides an increasingly more accurate approximation of Pf(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 (Fig. 3 A), in other cases the Edgeworth expansion deviates significantly from Pf(m) in all orders we have considered (Fig. 3 B). 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.
|


m ((8

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 Pf(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, Pf(m) is extremely sensitive to small deviations from normality in the tail of the m-fold convolved 
G distribution.
|


by first calculating the prediction for Pf(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 Pf(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. Fig. 5 A shows that the 

predicted by the convolution method correlates strongly with the measured (overall R2 for all 300 data points R2 = 0.789, p < 1015), in agreement with our earlier observation that, overall, our model works very well.
|


from the 
G distribution follows from large-deviation probability theory. Cramér's theorem implies that Pf(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 Fig. 5 B, 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 R2 for all 300 data points R2 = 0.971, p < 1015).
The intuitive explanation for why Pf(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 Pf(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 Pf(m + 1). With this assumption, 

turns out to be the largest eigenvalue of a matrix Wij 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 R2 for all 300 data points R2 = 0.9992, p < 1015). 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). Fig. 5 C 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 R2 for all 300 data points R2 = 0.939, p < 1015).
Finally, we have generated an additional data set of 10 x 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
Gcut = 4.0 kcal/mol, 5.0 kcal/mol, and 6.0 kcal/mol do not have a statistically significant (at the 5% level) difference in
|
| DISCUSSION |
|---|


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 Pf(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 |
|---|
where Xi 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 nth cumulant (see Appendix E) of the 
G distribution. We define
and write the standard normal distribution function
![]() | (2) |
![]() | (3) |
![]() | (4) |

G distribution) is
![]() | (5) |


m with 

= exp [µ2/(2
2)]. | APPENDIX B: CRAMÉR APPROXIMATION |
|---|
(as introduced in Appendix A) is approximately normally distributed with mean mµ and variance m
2, where µ and
2 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
![]() | (6) |
(t) is the moment-generating function of the distribution of Xi, 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(Sm/m
a), namely, for large m,
![]() | (7) |
![]() | (8) |


(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 Pf(m).
| APPENDIX C: MARKOV CHAIN APPROXIMATION |
|---|


of Pf(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 pi(m) be the fraction of proteins at mutation distance m in bin i. Clearly, we have
Next, we introduce the matrix Mij, 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, Mij 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 Pf(m + 1) as
![]() | (9) |
![]() | (10) |

G distribution. From here on, we will refer to the submatrix of Mij for which the indices i and j are non-negative as Wij. The distinction between Mij and Wij will become important when we discuss eigenvectors and eigenvalues of Wij.
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) |
with Pf(m), we find that Pf(m + 1) and Pf(m) are related to each other via
![]() | (12) |

(m)
approaches a limiting value 

for large m, we have
![]() | (13) |


m.
From
we see that for large m, the pi are proportional to the dominant eigenvector of Wij, by virtue of the Frobenius-Perron theorem (26
). (The Frobenius-Perron theorem holds if Wij is primitivethe case whenever there is a path of mutations that leads from any bin i to any other bin j, and Wii > 0 for at least one i.) Furthermore, Eq. 11 implies
![]() | (14) |


corresponds to the dominant eigenvalue of Wij. | APPENDIX D: MEAN-FIELD APPROXIMATION |
|---|


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 {
Gi} 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 
Gi < 
Gi+1 for all i. Then, we have
![]() | (15) |
| APPENDIX E: UNBIASED ESTIMATORS OF CUMULANTS |
|---|
be a set of n measurements, and define
![]() | (16) |
1
5 (note that
1 is the sample average, and
2 is the sample variance):
![]() | (17) |
![]() | (18) |
![]() | (19) |
![]() | (20) |
![]() | (21) |
| ACKNOWLEDGEMENTS |
|---|
Submitted on February 28, 2005; accepted for publication August 16, 2005.
| REFERENCES |
|---|
2. Bornberg-Bauer, E., and H. S. Chan. 1999. Modeling evolutionary landscapes: mutational stability, topology, and superfunnels in sequence space. Proc. Natl. Acad. Sci. USA. 96:1068910694.
3. Broglia, R. A., G. Tiana, H. E. Roman, E. Vigezzi, and E. Shakhnovich. 1999. Stability of designed proteins against mutations. Phys. Rev. Lett. 82:47274730.[CrossRef]
4. Chan, H. S., and E. Bornberg-Bauer. 2002. Perspectives on protein evolution from simple exact models. Appl. Bioinformat. 1:121144.
5. Wilke, C. O. 2004. Molecular clock in neutral protein evolution. BMC Genet. 5:25.[CrossRef][Medline]
6. Xia, Y., and M. Levitt. 2004. Simulating protein evolution in sequence and structure space. Curr. Opin. Struct. Biol. 14:202207.[CrossRef][Medline]
7. Drummond, D. A., B. L. Iverson, G. Georgiou, and F. H. Arnold. 2005. Why high-error-rate random mutagenesis libraries are enriched in functional and improved proteins. J. Mol. Biol. 350:806816.[CrossRef][Medline]
8. Bloom, J. D., J. J. Silberg, C. O. Wilke, D. A. Drummond, C. Adami, and F. H. Arnold. 2005. Thermodynamic prediction of protein neutrality. Proc. Natl. Acad. Sci. USA. 102:606611.
9. Daugherty, P. S., G. Chen, B. L. Iverson, and G. Georgiou. 1999. Quantitative analysis of the effect of the mutation frequency on the affinity maturation of single chain Fv antibodies. Proc. Natl. Acad. Sci. USA. 97:20292034.[CrossRef]
10. Guo, H. H., J. Choe, and L. A. Loeb. 2004. Protein tolerance to random amino acid change. Proc. Natl. Acad. Sci. USA. 101:92059210.
11. Shafikhani, S., R. A. Siegel, E. Ferrari, and V. Schnellenberger. 1997. Generation of large libraries of random mutants in Bacillus subtilis by PCR-based plasmid multimerization. Biotechniques. 23:304310.[Medline]
12. Loeb, D. D., R. Swanstrom, L. Everitt, M. Manchester, S. E. Stamper, and C. A. Hutchison. 1989. Complete mutagenesis of the HIV-1 protease. Nature. 340:397400.[CrossRef][Medline]
13. Pakula, A. A., V. B. Young, and R. T. Sauer. 1986. Bacteriophage
cro mutations: effects on activity and intracellular degradation. Proc. Natl. Acad. Sci. USA. 83:88298833.
14. Shortle, D., and B. Lin. 1985. Genetic analysis of staphylococcal nuclease: identification of three intragenic "global" suppressors of nuclease-minus mutations. Genetics. 110:539555.
15. Taverna, D. M., and R. A. Goldstein. 2002. Why are proteins so robust to site mutations? J. Mol. Biol. 315:479484.[CrossRef][Medline]
16. Kloczkowski, A., and R. L. Jernigan. 1997. Computer generation and enumeration of compact self-avoiding walks within simple geometries on lattices. Comput. Theor. Polym. Sci. 7:163173.[CrossRef]
17. Miyazawa, S., and R. L. Jernigan. 1996. Residue-residue potentials with a favorable contact pair term and an unfavorable high packing density term, for simulation and threading. J. Mol. Biol. 256:623644.[CrossRef][Medline]
18. Venables, W. N., and D. M. Smith. 2002. The R Development Core Team. An Introduction to R. Network Theory Ltd., Bristol, UK.
19. Benjamini, Y., and Y. Hochberg. 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. Roy. Stat. Soc. B. 57:289300.
20. Chan, H. S., and K. A. Dill. 1996. Comparing folding codes for proteins and polymers. Proteins Struct. Funct. Genet. 24:335344.[CrossRef][Medline]
21. Irbäck, A., and C. Troein. 2002. Enumerating designing sequences in the HP model. J. Biol. Phys. 28:115.[CrossRef]
22. Cramér, H. 1946. Mathematical Methods of Statistics. Princeton University Press, Princeton, NJ.
23. Blinnikov, S., and R. Moessner. 1998. Expansions for nearly Gaussian distributions. Astron. Astrophys. Suppl. Ser. 130:193205.[CrossRef]
24. Cramér, H. 1938. On a new limit theorem in the theory of probability. In Colloquium on the Theory of Probability. Hermann, Paris, France.
25. Hahn, M. G., and M. J. Klass. 1997. Approximation of partial sums of arbitrary i.i.d. random variables and the precision of the usual exponential upper bound. Annals Prob. 25:14511470.[CrossRef]
26. Varga, R. S. 2000. Matrix Iterative Analysis, 2nd Ed. Springer-Verlag, New York.
27. Dressel, P. L. 1940. Statistical semivariants and their estimates with particular emphasis on their relation to algebraic invariants. Annals Math. Stat. 11:3357.[CrossRef]
This article has been cited by other articles:
![]() |
C. O. Wilke and D. A. Drummond Population Genetics of Translational Robustness Genetics, May 1, 2006; 173(1): 473 - 481. [Abstract] [Full Text] [PDF] |
||||
| |||||||||||||||||||||||||||||||||||||