## Abstract

Recent work with *Saccharomyces cerevisiae* shows a linear relationship between the evolutionary rate of sites and the relative solvent accessibility (RSA) of the corresponding residues in the folded protein. Here, we aim to develop a mathematical model that can reproduce this linear relationship. We first demonstrate that two models that both seem reasonable choices (a simple model in which selection strength correlates with RSA and a more complex model based on RSA-dependent amino acid distributions) fail to reproduce the observed relationship. We then develop a model on the basis of observed site-specific amino acid distributions and show that this model behaves appropriately. We conclude that evolutionary rates are directly linked to the distribution of amino acids at individual sites. Because of this link, any future insight into the biophysical mechanisms that determine amino acid distributions will improve our understanding of evolutionary rates.

THE requirement for successful and efficient protein folding imposes significant biophysical constraints on coding sequences. These constraints shape how sequences evolve. Mutations that interfere with correct folding will generally be removed by purifying selection. Likewise, mutations that do not interfere with folding are often neutral, or nearly so, and accumulate over time. As a consequence of this interaction between protein biophysics and molecular evolution, signatures of protein structure can be found in the divergence patterns of coding sequences (Franzosa and Xia 2008; Lobkovsky*et al.* 2009; Wilke and Drummond 2010).

Mutagenesis experiments have shown that different positions in proteins have widely differing tolerances to amino acid substitutions (Reidhaar-Olson and Sauer 1988; Bowie*et al.* 1990; Lau and Dill 1990; Guo*et al.* 2004; Campbell-Valois*et al.* 2005; Smith and Raines 2006). On average, however, mutations introduced at solvent-exposed sites are less likely to disrupt protein structure and function than mutations introduced at buried sites. The latter tend to destabilize proteins, through steric hindrance and introduction of strained conformations (Chothia and Finkelstein 1990).

The higher tolerance of solvent-exposed sites to amino acid substitutions results in a correlation between the rate at which individual sites in coding sequences accumulate mutations over evolutionary time and the solvent exposure that these sites have in the expressed protein. Studies that have linked evolutionary rate with solvent exposure have consistently found that buried sites are more conserved and evolve slower than exposed sites (Overington*et al.* 1992; Goldman*et al.* 1998; Mirny and Shakhnovich 1999; Bustamante*et al.* 2000; Bloom*et al.* 2006; Conant and Stadler 2009; Franzosa and Xia 2009). At the same time, however, proteins with a larger core (more buried residues) evolve faster than proteins with a smaller core (Bloom*et al.* 2006; Ferrada and Wagner 2008; Zhou*et al.* 2008; Franzosa and Xia 2009). This apparent paradox can be resolved by observing that a larger core allows surface residues to vary more (Shakhnovich*et al.* 2005; Bloom*et al.* 2006; Franzosa and Xia 2009).

Recently, Franzosa and Xia (2009) developed a novel approach to analyze the relationship between evolutionary rate and solvent accessibility. They first mapped a large fraction of the genome of the yeast *Saccharomyces cerevisiae* onto homologous crystal structures from the Protein Data Bank (PDB). On the basis of this mapping, Franzosa and Xia (2009) determined relative solvent accessibility (RSA) for ∼300,000 sites in the yeast genome. They then grouped these sites into bins of similar RSA value and calculated for each bin the average evolutionary rate *d*_{N}/*d*_{S} in a phylogeny of four yeast species. They found a strikingly linear relationship between evolutionary rate and RSA. Every 1% increase in RSA was associated with an increase in *d*_{N}/*d*_{S} of 0.001 (Franzosa and Xia 2009). Why the relationship between evolutionary rate and RSA is linear remains unknown.

Here, we employ mathematical modeling and bioinformatics analysis to explore what mechanism could be responsible for the linear relationship. We first show that a two-allele model in which selection strength correlates with RSA fails to reproduce this relationship. We then develop a more sophisticated model on the basis of amino acid frequencies and show that this model fails as well. The second model fails because amino acid frequencies averaged over many sites with comparable RSA differ dramatically from the distributions of allowed amino acids at individual sites. By building a model on the basis of the latter distributions, we can reproduce an approximately linear relationship between evolutionary rate and RSA.

## METHODS

#### Evolutionary rate as a function of RSA:

To verify the linear relationship between evolutionary rate and RSA at the amino acid level, we reproduced Franzosa and Xia’s (2009) results using amino acid distance instead of *d*_{N}/*d*_{S}. First, we obtained orthologs between *S. bayanus* and *S. cerevisiae* from the Saccharomyces Genome Database as in Zhou*et al.* (2008) and aligned sequences with MUSCLE (Edgar 2004). We mapped the *S. cerevisiae* sequences to structures using three iterations of PSI-BLAST against the PDB, requiring a minimum of 80% sequence identity for a match. We ended up with 525 matching structures. For these matched structures we used the program DSSP (Kabsch and Sander 1983) to calculate solvent accessibility at each site. To obtain RSA, we normalized the solvent accessibilities calculated by DSSP with respect to an extended Gly-X-Gly peptide (Creighton 1992). We binned sites by RSA and then calculated evolutionary rate *K* with the PAML package codeml (Yang 2007), using the Whelan and Goldman (WAG) model for amino acid distance (Whelan and Goldman 2001).

#### Amino acid distribution over many yeast proteins:

To calculate amino acid distributions, we used the same set of *S. cerevisiae* ORFs mapped to protein structures. We binned all sites by RSA as above. (A few residues had RSA > 1 and we treated them as if they had RSA = 1.) We then calculated the relative frequency of each amino acid in each RSA bin. For visualization, we ordered amino acids by hydrophobicity using the Fauchere–Pliska octanol scale (Fauchere and Pliska 1983).

#### Coordination number and RSA correlation:

We computed the correlation between normalized coordination number and RSA using the same set of *S. cerevisiae* proteins as above. The coordination number of a site is the number of sites it is in contact with, and we considered two sites to be in contact if any two heavy atoms are within 4.5 Å of each other (excluding sequence neighbors). We used the BioPython module Bio.PDB (Hamelryck and Manderick 2003) to compute coordination numbers, which for each site we normalized by the average over the entire protein.

#### Variation at individual sites over structural homologs:

To compute distributions at individual sites across structurally similar proteins, we employed a PSI-BLAST search of the NCBI nonredundant database (NR) to construct alignments from various seed proteins. As seed proteins, we used the proteins obtained by mapping the yeast genome to the PDB (as described above). We then filtered the alignment given by PSI-BLAST such that the remaining sequences all had between 40% and 80% pairwise sequence similarity with all other sequences in the alignment. This filtering procedure excluded redundant sequences while still ensuring structural similarity (Chothia and Lesk 1986; Holm*et al.* 1992). We retained only filtered alignments that contained at least 50 sequences. Our final data set consisted of 162 distinct alignments. In the filtered alignments, we classified each site by RSA of the seed protein at this site and placed sites into bins of similar RSA. We then calculated the alignment-wide amino acid distribution for every site. At each site, we ranked residues by declining frequency at that site. We then averaged the frequency-sorted amino acid distributions over all sites within each bin.

To characterize these averaged distributions with a single parameter, we fitted the one-parameter exponential function *e*^{−λk} to the average amino acid frequency as a function of the amino acid rank *k*.

Analysis scripts and data to reproduce this analysis are provided in supporting information, File S1.

#### Parameter choices:

To study numerically the behavior of our mathematical models of protein evolution, we had to choose suitable values for the parameters *N*_{e} (effective population size) and μ (mutation rate). We chose values that are approximately correct for yeast, namely *N*_{e} = 5 × 10^{6} individuals and μ = 3.3 × 10^{−10} mutations per site per generation (Lynch*et al.* 2008; Lancaster*et al.* 2010).

## RESULTS

Franzosa and Xia (2009) found a strong linear relationship between *d*_{N}/*d*_{S} and RSA. While their result was likely driven by selection on the amino acid level, their use of *d*_{N}/*d*_{S} does not allow us to draw this conclusion *a priori*. Their result could be confounded by varying levels of selection on synonymous sites; synonymous codon usage is not uniform across genes and covaries with protein structure (Akashi 1994; Drummond and Wilke 2008; Zhou*et al.* 2009; Lee*et al.* 2010).

Therefore, to verify that Franzosa and Xia (2009) had indeed identified an effect occurring at the amino acid level, we repeated their analysis with amino acid sequences. We aligned orthologous genes from the two yeasts *S. cerevisiae* and *S. bayanus* and classified sites into bins of similar RSA values. We then concatenated all sites within each bin and calculated the amino acid distance *K* between the *S. cerevisiae* and the *S. bayanus* sequence in each bin. Amino acid distance is a measure of evolutionary rate on the amino acid level (Whelan and Goldman 2001).

We found a near-perfect linear relationship between evolutionary rate *K* and RSA (Figure 1). We interpret this result as a signal of purifying selection acting on the amino acid sequence. On average, buried sites experience stronger purifying selection than exposed sites and thus evolve slower. The increased selective constraints on buried amino acids presumably reflect the requirement for proteins to fold and function properly.

That buried sites are more constrained than exposed sites is well known. Much existing theory, experiments, and sequence data support the notion that substitutions in the core of a protein are more likely to be disruptive than substitutions in solvent-exposed regions. Yet the perfectly linear relationship between evolutionary rate and RSA is surprising and deserves an explanation. We thus proceeded to explore what kind of evolutionary models could potentially reproduce this observation.

#### A simple two-allele model:

The simplest model we can consider is a multiplicative multisite, two-allele model; in this model, an organism's genome consists of a finite number of sites, each of which can exist in two alleles. All sites contribute multiplicatively to the overall fitness of the organism. At each site *i*, one of the two alleles is preferred, and we assume it has fitness 1. The second allele is selected against and has fitness 1−*s _{i}*. We assume that all sites mutate with the same rate μ. In such a model, in equilibrium, sites with larger

*s*will evolve slower than sites with smaller

_{i}*s*. For sufficiently small

_{i}*s*, sites will evolve neutrally at rate μ.

_{i}Here and throughout, we consider haploid, asexual organisms and assume that the product of mutation rate and effective population size *N*_{e} is small, μ*N*_{e} ≪ 1. In this case, and because we consider a multiplicative model, the evolutionary rate of a genome of length *L* is the average of the evolutionary rates of *L* single-site models with identical selection coefficients. Therefore, in what follows, we consider only the evolutionary rate at a single site and ask how it changes with selection coefficient *s*. For simplicity, we drop the site index *i*.

We refer to the two alleles at a site as *A* and *a*. Allele *A* has fitness 1 and allele *a* has fitness 1 – *s*. The probability that allele *a* goes to fixation in a background of allele *A* is given by Kimura (1962):(1)Likewise, the probability that allele *A* goes to fixation in a background of allele *a* is given by(2)

In equilibrium, and averaged over long periods of time, both alleles will be present at the site some fraction of time. We denote these fractions as *F*(*A*) and *F*(*a*), with *F*(*A*) + *F*(*a*) = 1. We have(3)Evolutionary rate *K* is the rate with which mutations originate and go to fixation. Thus, *K* is given by(4)

The evolutionary rate *K* is of course a function of *s*. Thus, we can now ask how *K* changes as *s* changes. Assuming *s* ≪ *N*_{e} and using standard approximations for the fixation probabilities, we obtain(5)For *s* 1/*N*_{e}, evolution is neutral, and *K*(*s*) ≈ μ. For larger *s*, the evolutionary rate declines exponentially in *s*, .

We now assume that the selection coefficient *s* is a function of RSA. We denote RSA by *r* in mathematical expressions. An increase in *K* as *r* increases corresponds to a greater tolerance to mutation; hence, the selection coefficient *s*(*r*) should be a decreasing function of *r*. We assume that *r* can take on any value in the interval [0, 1]. The function *s*(*r*) maps this interval into some interval of *s* values. Thus, we have to ask: Is there a reasonable mapping from *r* to *s*(*r*) such that *K*(*s*(*r*)) is approximately a linear, increasing function in *r*? We found that generally, such a mapping does not exist. Because *K* decreases exponentially with *s*, any function *s*(*r*) that might result in approximately linear behavior of *K* will necessarily have an exponentially small range of possible *s* values. To illustrate this result, we defined three functions with parameters that give similar ranges in [0, 1]: a linearly decaying function whose image is [0, 1.5 × 10^{−5}], an exponential function with image [7.3 × 10^{−6}, 2 × 10^{−5}], and a logarithmic function spanning ∼[0, 1.8 × 10^{−5}]. For all three definitions of *s*(*r*), Equation 4 still produces exponentially fast growth of *K* as a function of *r* (Figure 2). More generally, we can show that even if the difference in *s* corresponding to fully buried (*r* = 0) and fully exposed (*r* = 1) sites is only on the order of 1/*N*_{e}, the deviation from linearity is larger than the magnitude of the evolutionary rate *K* itself (see appendix). We conclude that the two-allele model does not seem to be an appropriate model to describe the effect of relative solvent accessibility on evolutionary rate.

#### A model based on amino acid frequencies:

We believe that the main reason why the two-allele model gives unsatisfactory results is that it replaces 20 different amino acids by only two different states, preferred and unpreferred. In real proteins, it may well be that at one site 3 amino acids are preferred and 17 unpreferred, while at a different site 5 are preferred and 15 unpreferred. All else being equal, the second site will evolve faster than the first. This reasoning suggested to us that we should aim to develop a model on the basis of amino acid frequencies. The sites with the broadest distributions of amino acids should evolve the fastest, and the sites with the narrowest distributions the slowest.

Amino acid distributions in proteins have been studied extensively. The general consensus is that amino acid frequencies follow a Boltzmann distribution. The individual frequencies at sites can be calculated either from stability effects [ΔΔ*G* values (Dokholyan and Shakhnovich 2001; Dokholyan*et al.* 2002; Godoy-Ruiz*et al.* 2004; Bloom and Glassman 2009; Schmidt am Busch *et al.* 2010] or from the protein's connectivity matrix (Porto*et al.* 2004; Bastolla*et al.* 2005; Pokarowski*et al.* 2005; Wolff*et al.* 2008; Bastolla*et al.* 2008). In particular, Porto*et al.* (2004) showed that the frequency of amino acid *a* is proportional to *e*^{−βh(a)}, where β measures properties of the site under consideration and *h*(*a*) measures properties of the amino acid. The quantity β can be derived from the protein structure's contact matrix. It varies almost linearly with the site's coordination number normalized by the protein's average. The quantity *h*(*a*) is the *interactivity* of amino acid *a*, a quantity highly correlated with hydrophobicity (Bastolla*et al.* 2005).

Because solvent occlusion happens through interresidue contacts, we hypothesized that the normalized coordination number should correlate strongly with RSA and that the theory of Porto*et al.* (2004) should provide at least a qualitatively correct description of the amino acid distribution in different RSA bins. We found both to be the case in yeast. The normalized coordination number correlated well with RSA (Pearson's *r* = 0.66, *P* < 2.2 × 10^{−16}). Amino acid distributions were strongly skewed toward hydrophobic residues at low RSA and toward hydrophilic residues at high RSA. For intermediate RSA, corresponding to β = 0, both hydrophobic and hydrophilic residues had comparable frequencies (Figure S1). Having found this correspondence, we proceeded to obtain the evolutionary rates predicted by the theory of Porto*et al.* (2004).

The amino acid distribution at a site, combined with effective population size *N*_{e} and mutation rate μ, fully specifies the evolutionary rate at the site, under the assumption that sites evolve independently. The link between amino acid distribution and evolutionary rate is established by Sella–Hirsh theory (Sella and Hirsh 2005). This theory demonstrates that equilibrium frequencies of alleles follow a Boltzmann distribution just like the one found by Porto*et al.* (2004). Thus, from the equilibrium frequencies of alleles we can infer the relative fitness of alleles and their fixation probabilities in various backgrounds.

According to Porto*et al.* (2004), the distribution of amino acids is given by(6)where *a* is a specific amino acid as before and the sum in the denominator runs over all 20 amino acids. Fixation probabilities follow as(7)(Sella and Hirsh 2005). (These fixation probabilities are equivalent to the Kimura probabilities used in the previous subsection; see Sella and Hirsh 2005 for details.) We can now express evolutionary rate in terms of amino acid distribution and fixation probabilities as(8)Remember that *K* is a function of β, and β is an approximate measure of solvent accessibility. Highly buried sites will have a large negative β, highly exposed sites will have a large positive β, and intermediate sites will have a β close to zero. Thus, to be consistent with data (*e.g.*, Figure 1), Equation 8 should be an increasing function of β. Instead, however, we found that Equation 8 predicts *K* to be maximal at β = 0 (Figure 3A) and to decline in both directions as the absolute value of β increases. This result makes intuitive sense, as the distribution defined by Equation 6 is the broadest for β = 0. However, we have to conclude that the theory of Porto*et al.* (2004) cannot be used to explain the linear relationship between evolutionary rate and RSA.

We emphasize that the failure of Equation 8 does not imply that the amino acid distributions calculated by Porto*et al.* (2004) and given by Equation 6 are incorrect. In fact, we used Equations 7 and 8 to predict evolutionary rates from the observed amino acid frequencies in yeast and found similarly that the predicted evolutionary rate peaked at intermediate RSA (Figure 3B).

#### An alternative model based on amino acid frequencies:

The failure of the previous model implies that the model missed some important aspect of protein evolution. We hypothesized that the model failed because Equation 6 was valid only for the entire class of sites with similar β, but not for any individual site in this class. It is entirely possible that the distribution of amino acids at a specific site, when observed over evolutionarily long periods of time, does not agree with Equation 6, even though the average distribution of all sites with similar β or RSA does. Both previously published tests of Equation 6 (Porto*et al.* 2004) and our amino acid distributions as a function of RSA (Figure S1) were obtained by averaging over many sites and thus would not reveal any deviation from Equation 6 at individual sites.

To determine the distribution of amino acids at individual sites, we built large alignments of structurally similar proteins (see methods). We found that the distributions at individual sites were highly variable and looked nothing like Equation 6. In general, at any given site, only a small number of different amino acids were actually present, and there was often no obvious relationship between which amino acids were present and what their hydrophobicity was. However, when averaging over many sites with similar RSA, we could recover distributions comparable to Equation 6.

Even though the specific amino acids preferred at individual sites were highly variable, we found that the frequency distributions at different sites were similar. When we ordered amino acids by their relative frequency at each site, we found that the frequencies were proportional to an exponential, exp(−λ*k*), where *k* counts amino acids in descending order of frequency, *k* = 0, 1, … , 19. We averaged the reordered amino acid distributions over all sites within bins of similar RSA (Figure 4A) and fitted exp(−λ*k*) to these averaged distributions. We thus obtained λ as a function of RSA and found that λ decayed approximately linearly with RSA (Figure 4B and Figure S2).

We carried out this analysis on 162 yeast proteins and found that generally (i) λ was approximately a linear function of RSA and (ii) λ decayed with increasing RSA (Figure 5). For each protein, we fitted a linear function λ(*r*) = *c*_{1} + *c*_{2}*r* to the data and generally found a negative slope *c*_{2} and a good model fit. The few cases with an apparent positive slope *c*_{2} could be traced back to a single outlying λ-value at the highest RSA bin (see Figure S3 for an example). This bin generally encompassed the fewest number of sites (see also discussion) and thus its λ-value was not always reliable.

On the basis of these findings, we can model the evolutionary process at individual sites such that it produces steady-state amino acid frequencies(9)where *a* and *b* index amino acids, in the appropriate order, and run from 0 to 19. The parameter λ declines with RSA.

As in the previous subsection, we can use the Sella and Hirsh (2005) method to map these steady-state frequencies onto a unique evolutionary process. The fitness values for individual amino acids are given by(10)It might seem disconcerting that we measure fitness here in units of the effective population size *N*_{e}. After all, the fitness contribution of a particular amino acid in a particular protein of an organism should not depend on the size of the population of that organism. However, this scaling by population size is merely a mathematical convenience to keep the actually observable quantities (amino acid distributions, evolutionary rates) free of any explicit dependency on *N*_{e}. For real organisms, we expect that *w*(*a*) is independent of *N*_{e} but that λ, *F*(*a*), and evolutionary rate *K* all depend on *N*_{e}.

Fixation probabilities follow from Equation 7. Making the approximation *N*_{e} − 1 ≈ *N*_{e}, we find(11)We obtain the average evolutionary rate for this model by substituting Equations 9 and 11 into Equation 8. We find(12)where is the partition function.

For large *N*_{e}, we can approximate so that(13)The absence of *N*_{e} from this equation shows that if we scale *w*(*a*) with *N*_{e}, as in (10), then *K* is approximately independent of *N*_{e}.

To obtain evolutionary rate as a function of RSA, we substitute λ = *c*_{1} + *c*_{2}*r* into Equation 13. Figure 6 shows resulting evolutionary rates for three representative proteins. The curves *K*(*r*) are roughly linear and *K* is approximately of the correct order of magnitude. However, *K*(*r*) is not perfectly linear; there is some clear upward curvature. The curvature tends to increase with the absolute magnitude of *c*_{2}. We comment on this issue in the discussion. Also, note that the units for *K* are not the same in Figure 1 as they are in the other figures. In Figure 1, *K* is estimated as the number of substitutions per site per unit time. The time unit is the total divergence time between the species that are being compared. By contrast, our mathematical models predict *K* in units of substitutions per site per generation. We estimate that ∼10^{11} generations separate *S. cerevisiae* and *S. bayanus*, 40 million yr × 4000 generations/yr.

## DISCUSSION

We have shown that the linear relationship between evolutionary rate and RSA reflects a selection pressure on the amino acid level. Further, we have demonstrated that a simple two-allele model and a more elaborate model based on observed mean amino acid frequencies for sites with similar RSA cannot reproduce this linear relationship. The first model fails because it is too simplistic; individual sites in proteins can, at least in principle, assume 1 of 20 different states. The second model fails because amino acid frequencies averaged over many sites are not representative of amino acid frequencies at individual sites. We have found that the latter frequencies follow a Boltzmann distribution that becomes increasingly broad as RSA increases. Finally, we have shown that a mathematical model based on this observation can reproduce the linear relationship between evolutionary rate and relative solvent accessibility.

Our analysis highlights how important it is to distinguish between amino acid frequencies averaged over a large class of sites with similar property (such as RSA) and amino acid frequencies at individual sites. In both cases, frequencies are Boltzmann distributed, and thus it is easy to mistake one for the other. However, the properties of these two distributions are very different. For example, in yeast, at sites with RSA close to 0.2 nearly all amino acids occur at comparable frequencies. Yet at any given site, only a small number of amino acids are actually permissible. Evolutionary rate, which measures the rate at which mutations at individual sites arise and go to fixation, is governed by the amino acid distribution of individual sites, not the average distribution over a broad class of sites.

However, averaging distributions of similarly exposed sites from many proteins seems to agree qualitatively with distributions predicted by Porto*et al.* (2004). This agreement suggests that any future theory attempting to predict site-specific distributions should also be able to predict average distributions of sites with similar β (or RSA). These average distributions should reduce to something similar to the theory of Porto*et al.* (2004) and the data shown in Figure S1.

Our model describes the variation in steady-state distribution at sites using the exponential parameter λ, which we defined above as a linear function of RSA: λ(*r*) = *c*_{1} + *c*_{2}*r*. In this way λ(*r*) describes the level of variation in the distribution function (Equation 9) for a given RSA. The intercept of λ(*r*), the largest value it takes, corresponds to the strongest selective pressure and the minimal level of variation for the most buried sites, at *r* = 0. This maximal selective pressure in turn determines the value of the minimal evolutionary rate. Likewise, the slope of λ(*r*) determines the rate of increase of *K*(*r*): a steeper slope (more negative *c*_{2}) signifies a greater tolerance of alternative residues as *r* increases compared to a shallower slope (less negative *c*_{2}), and greater tolerance of alternative residues implies a greater increase in *K* as *r* grows.

We emphasize that different RSA bins contain different numbers of sites (see also Figure 2 from Franzosa and Xia 2009). Bins below RSA values of 0.1 tend to contain more than twice as many sites as bins for RSA values between 0.1 and 0.6. Bins for higher RSA values are even less occupied. In our data set of all yeast genes, we have 69,521 sites in the lowest-RSA bin but only 1452 sites in the highest-RSA bin. Because of the comparative scarcity of high-RSA sites, our estimates for amino acid distributions at these sites are not always reliable, as exemplified in Figure S3. In our experience, the amino acid distributions at high RSA are reliable when a linear model produces a negative slope for λ(*r*) and they are unreliable otherwise.

In our analysis of amino acid distributions at individual sites in individual proteins, we generally observed only a few (on the order of 5) different amino acids at each site. This outcome was expected for Boltzmann-distributed amino acid frequencies. For the proteins we investigated, we found that λ typically fell somewhere between 0.3 and 1.2. Even for the smallest λ in this range, λ = 0.3, the expected frequency of the 10th most abundant amino acid under a Boltzmann distribution is only ∼0.02, and the expected frequency of the 20th most abundant amino acid is ∼0.001. For larger λ, the expected frequencies are much smaller. In our alignments, which mostly ranged from 50 sequences to ∼200 sequences, with a few cases going up to 360 sequences, we could not properly sample amino acids that have such low expected abundances. In fact, in our distributions, the least abundant amino acid generally has absolute frequencies in low single digits, and thus we cannot expect to see any other amino acids that should, according to theory and the overall pattern we see, arise at less than single-digit frequencies.

By measuring fitness in units of *N*_{e} for ease of analysis, we have implicitly made where What we have then is a relation linking the original λ(*r*) to *N*_{e} and the parameters and . Note that the original λ(*r*) is a statistically measurable function describing variation at sites by RSA. If we could obtain estimates of and independently of *K*, say from an *ab initio* model of protein folding, and then the relationships were formally attached to biophysical quantities that proved reliably measurable, this relationship could provide a novel method by which to estimate effective population size.

While our final model produces an approximately linear relationship between evolutionary rate and RSA, the model predictions are not perfectly linear. In particular for proteins with larger absolute *c*_{2} values, we see a clear upward curvature in evolutionary rate as a function of RSA (Figure 6). In our modeling approach, we made several approximating assumptions, and each of them could potentially be the source of the curvature. First, we assumed that amino acid distributions are Boltzmann distributed. This assumption may not be entirely correct. In fact, if amino acid distributions were perfectly Boltzmann distributed, then the data in Figure 4A should be perfectly linear. Instead, they seem to display a moderate amount of curvature. Second, λ may not be a linear function of RSA. We did see a fair amount of noise in λ for some proteins (*e.g.*, Figure S2D), but we did not see any systematic deviation from the linear trend. Third, when modeling how amino acid distributions relate to evolutionary rate, we completely neglected any interactions among sites. While models without interactions have been successful in related studies, *e.g.*, in predicting the effect of multiple mutations on protein stability (Bloom*et al.* 2005) and in linking mutation frequencies to stability effects (ΔΔ*G* values) (Godoy-Ruiz*et al.* 2004; Zeldovich*et al.* 2007; Bloom and Glassman 2009), epistatic interactions among sites in proteins are well documented and may be important for precise prediction of evolutionary rates.

## Acknowledgments

We thank Markus Porto and Eugene Shakhnovich for helpful comments on this work. This work was supported by National Institutes of Health grant R01 GM088344 and by the National Science Foundation under Cooperative Agreement DBI-0939454.

## APPENDIX

In the main body of this article, we have shown that in the two-allele model the evolutionary rate declines exponentially in *s* for *s* > 1/*N*_{e}. We may ask whether it is possible for *s*(*r*) to map *r* to a sufficiently small range [*s*_{1}, *s*_{2}] ⊂ [0, 1] so that *K*(*r*) is approximately linear over that range. To this end, take *s*_{1}, *s*_{2} with 1/*N*_{e} < *s*_{1} < *s*_{2}.

We judge linearity in the range [*s*_{1}, *x*] by the magnitude of the function(A1)where *L*(*x*)=*K*′(*s*_{1})(*x* − *s*_{1}) + *K*(*s*_{1}) is the line tangent to *K* at *s*_{1}. We examine the behavior of *D*(*s*_{2}) for a fixed distance ε = *s*_{2} − *s*_{1}. Substituting and *s*_{2} = *s*_{1} + ε, we find for Equation A1,(A2)(A3)This function decreases with both ε and *s*_{1}. Setting *s*_{2} = *s*_{1} + 1/*N*_{e} gives us(A4)(A5)Wherever the approximations based on *s*_{1} > 1/*N*_{e} are valid, a value of ε on the order of 1/*N*_{e} gives a difference between *K*(*s*_{2}) and *L*(*s*_{2}) larger than the magnitude of *K*(*s*_{1}) itself. Thus, in the two-allele model, the relationship between *K* and RSA remains highly nonlinear even if the difference in selection pressure at fully exposed and fully buried sites becomes minute.

## Footnotes

Communicating editor: L. M. Wahl

- Received February 21, 2011.
- Accepted March 16, 2011.

- Copyright © 2011 by the Genetics Society of America