Keywords

• Neutral theory;
• Positive selection;
• Evolutionary rate;
• Structural constraints;
• Protein abundance

1. Introduction

The common understanding of protein evolution has been that amino acid substitutions observed in homologous proteins are selectively neutral (Kimura, 1968, Kimura, 1969, Kimura and Ohta, 1971 and Kimura and Ohta, 1974) or slightly deleterious (Ohta, 1973 and Ohta, 1992), and random drift is a primary force to fix amino acid substitutions in population. The rate of protein evolution has been understood to be determined primarily by the proportion of neutral mutations, which may be measured by the ratio of nonsynonymous to synonymous substitution rate per site (K_a/K_s$K_a/K_s$) (Miyata and Yasunaga, 1980) and determined by functional density (Zuckerkandl, 1976) weighted by the relative variability at specific-function sites of a protein (Go and Miyazawa, 1980). Since then, these theories have been widely accepted, however, recently a question has been raised on whether the diversity of protein evolutionary rate among genes can be explained only by the proportion and the variability of specific-function sites, and molecular and population-genetic constraints on protein evolutionary rate have been explored.

Recent works have revealed that protein evolutionary rate is correlated with gene expression level; highly expressed genes evolve slowly, accounting for as much as 34% of rate variation in yeast (Pál et al., 2001). Of course, there are many reports that support a principle of lower evolution rate for stronger functional density. Broadly expressed proteins in many tissues tend to evolve slower than tissue-specific ones (Kuma et al., 1995 and Duret and Mouchiro, 2000). The connectivity of well-conserved proteins in a network is shown (Fraser et al., 2002) to be negatively correlated with their rate of evolution, because a greater proportion of the protein sequence is directly involved in its function. A fitness cost due to protein–protein misinteraction affects the evolutionary rate of surface residues (Yang et al., 2012). Protein dispensability in yeast is correlated with the rate of evolution (Hirsh and Fraser, 2001 and Hirsh and Fraser, 2003), although there is a report insisting on no correlation between them (Pál et al., 2003). Other reports indicate that the correlation between gene dispensability and evolutionary rate, although low, is significant (Zhang and He, 2005, Wall et al., 2005 and Jordan et al., 2002).

It was proposed (Drummond et al., 2005, Drummond and Wilke, 2008 and Geiler-Samerotte et al., 2011) that low substitution rates of highly expressed genes could be explained by fitness costs due to functional loss and toxicity (Stoebel et al., 2008 and Geiler-Samerotte et al., 2011) of misfolded proteins. Misfolding reduces the concentration of functional proteins, and wastes cellular time and energy on production of useless proteins. Also misfolded proteins form insoluble aggregates (Geiler-Samerotte et al., 2011). Fitness cost due to misfolded proteins is larger for highly expressed genes than for less expressed ones.

Fitness cost due to misfolded proteins was formulated (Drummond and Wilke, 2008 and Geiler-Samerotte et al., 2011) to be related to the proportion of misfolded proteins. Knowledge of protein folding indicates that protein folding primarily occurs in two-state transition (Miyazawa and Jernigan, 1982a and Miyazawa and Jernigan, 1982b), which means that the ensemble of protein conformations are a mixture of completely folded and unfolded conformations. Free energy (ΔG$\mathrm{\Delta }G$) of protein stability, which is equal to the free energy of the denatured state subtracted from that of the native state, and stability change (ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$) due to amino acid substitutions are collected in the ProTherm database (Kumar et al., 2006), although the data are not sufficient. Prediction methods, however, for ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ are improved enough to reproduce real distributions of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ (Schymkowitz et al., 2005 and Yin et al., 2007). Therefore, on the biophysical basis, the distribution of fitness can be estimated and protein evolution can be studied. Shakhnovich group studied protein evolution on the basis of knowledge of protein folding (Serohijos and Shakhnovich, 2014 and Dasmeh et al., 2014) and showed (Serohijos et al., 2012) that the negative correlation between protein abundance and K_a/K_s$K_a/K_s$ was caused by the distribution of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ that negatively correlates with the ΔG$\mathrm{\Delta }G$ of a wild type. Also, it was shown (Serohijos et al., 2013) that highly abundant proteins had to be more stable than low abundant ones. Relationship between evolutionary rate and protein stability is studied from various points of view (Echave et al., 2015 and Faure and Koonin, 2015).

Here we study relationship between evolutionary rate and selection on protein stability in a monoclonal approximation. A fitness assumed here for a protein is a generic form to which all formulations (Drummond and Wilke, 2008, Geiler-Samerotte et al., 2011, Serohijos et al., 2012, Serohijos et al., 2013, Serohijos and Shakhnovich, 2014 and Dasmeh et al., 2014) previously employed for protein fitness are reduced in the condition of exp(βΔG)⪡1$\exp \left(\beta \mathrm{\Delta }G\right)⪡1$, which is satisfied in the typical range of folding free energies shown in Fig. 1; β=1/(kT)$\beta =1/\left(\mathit{kT}\right)$, k is the Boltzmann constant and T   is absolute temperature. The generic form of Malthusian fitness of a protein-coding gene is $m\equiv -\kappa \exp \left(\beta \mathrm{\Delta }G\right)$, where κ   is a parameter, which may be a function of protein abundance and dispensability; see Methods for details. The distribution of stability change ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ due to single amino acid substitutions is approximated as a weighted sum of two Gaussian functions that was shown (Tokuriki et al., 2007) to well reproduce actual distributions of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$. One of the two Gaussian functions describes substitutions at structurally less-constrained surface sites, and the other at more-constrained core sites of proteins. The proportion of less-constrained surface sites is a parameter (θ).

The fixation probability of a mutant with ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ can be calculated for a duploid population with effective population size N_e (Crow and Kimura, 1970). In the population of genes with such a fitness protein stability is evolutionarily maintained at equilibrium, and equilibrium stability (ΔG_e$\mathrm{\Delta }{G}_e$) negatively correlates with protein abundance/dispensability (κ  ). The range of ΔG_e$\mathrm{\Delta }{G}_e$ is consistent with the observed range of folding free energies shown in Fig. 1.

The probability density functions (PDF) of K_a/K_s$K_a/K_s$, the ratio of nonsynonymous to synonymous substitution rate per site (Miyata and Yasunaga, 1980), at equilibrium and also in the vicinity of equilibrium are numerically examined over a whole domain of the parameters, $0\le \log 4N_e\kappa \le 20$ and 0≤θ≤1$0\le \theta \le 1$. The dependences of evolutionary rate on protein abundance/dispensability and on structural constraint are quantitatively described, and it is shown that both factors cannot be ignored on protein evolutionary rate, although protein abundance/indispensability more affect evolutionary rate for less constrained proteins, and structural constraint for less abundant, less essential proteins. Like protein abundance, protein indispensability must correlate with evolutionary rate, but a correlation between them may be hidden by the variation of protein abundance as well as effective population size, and detected only in low-abundant proteins. It has also become clear that nearly neutral selection is predominant only in low-abundant, non-essential proteins with $\log 4N_e\kappa <2$ or ΔG_e>−2.5$\mathrm{\Delta }{G}_e>-2.5$ kcal/mol, and in the other proteins positive selection is significant to more stabilize a less-stable wild type. Also, a significant amount of slightly negative mutants are fixed in population by random drift. This view of protein evolution is contrary to the previous understanding. The present model based on a biophysical knowledge of protein stability also indicates that protein stability (−βΔG_e)$(-\beta \mathrm{\Delta }{G}_e)$ and the average of K_a/K_s$K_a/K_s$ decrease as growth temperature increases.

2. Methods

2.1. Fitness costs due to misfolded proteins

Misfolding can impose costs in three distinct ways (Geiler-Samerotte et al., 2011); loss of function, diversion of protein synthesis resources away from essential proteins, and toxicity of the misfolded molecules. Fitness cost due to functional loss was formulated (Drummond and Wilke, 2008) by taking account of protein dispensability. Assuming that fitness cost of each gene is additive in the Malthusian fitness scale, the total Malthusian fitness of a genome was estimated as

equation(1)

$m_{\mathrm{dispensability}}\equiv -\sum _i{\gamma }_i(1-f_i^{\mathrm{native}})$

Turn MathJax on

$f_i^{\mathrm{native}}$

Protein folding primarily occurs in the two-state transition, which means that protein conformations are a mixture of completely folded and unfolded conformations (Miyazawa and Jernigan, 1982a and Miyazawa and Jernigan, 1982b). Therefore, if the completely folded (native) state is more stable by a free energy difference ΔG$\mathrm{\Delta }G$ than the unfolded (denatured) state, then the native fraction in the conformational ensemble will be equal to

equation(2)

$f^{\mathrm{native}}=\frac{e^{-\beta \mathrm{\Delta }G}}{1+e^{-\beta \mathrm{\Delta }G}}$

Turn MathJax on

Thus, Eq. (1) for the Malthusian fitness of a genome can be transformed as follows in terms of the folding free energy ΔG$\mathrm{\Delta }G$ of the native conformation:

equation(3)

$m_{\mathrm{dispensability}}=-\sum _i{\gamma }_i\frac{e^{\beta \mathrm{\Delta }{G}_i}}{e^{\beta \mathrm{\Delta }{G}_i}+1}$

Turn MathJax on

$\exp \left(\beta \mathrm{\Delta }G\right)⪡1$

equation(4)

$m_{\mathrm{dispensability}}=-\sum _i{\gamma }_i\left[e^{\beta \mathrm{\Delta }{G}_i}-O\left(e^{2\beta \mathrm{\Delta }{G}_i}\right)\right]$

Turn MathJax on

Drummond and Wilke (2008) took notice of toxicity of misfolded proteins as well as diversion of protein synthesis resources, and formulated the Malthusian fitness (m_misfolds$m_{\mathrm{misfolds}}$) of a genome to be negatively proportional to the total amount of misfolded proteins, which must be produced to obtain the necessary amount of folded proteins (Serohijos et al., 2012):

equation(5)

$m_{\mathrm{misfolds}}=-c\sum _i{A}_i\frac{1-f_i^{\mathrm{native}}}{f_i^{\mathrm{native}}}$

Turn MathJax on

equation(6)

$m_{\mathrm{misfolds}}=-c\sum _i{A}_i{e}^{\beta \mathrm{\Delta }{G}_i}$

Turn MathJax on

2.2. Fitness of a linear metabolic pathway

Serohijos and Shakhnovich (2014) examined the evolution of a linear metabolic pathway whose Wrightian fitness was defined as

equation(7)

$w_{\mathrm{linear}\mathrm{pathway}}\equiv w_{\mathrm{flux}}+w_{\mathrm{misfolds}}$

Turn MathJax on

equation(8)

$w_{\mathrm{flux}}\equiv \frac{{\sum }_i{\epsilon }_i{A}_i^{-1}}{{\sum }_i{\epsilon }_i{(A_i{f}_i^{\mathrm{native}})}^{-1}}$

Turn MathJax on

equation(9)

$w_{\mathrm{misfolds}}\equiv -c\sum _i{A}_i(1-f_i^{\mathrm{native}})$

Turn MathJax on

equation(10)

m_{linear_pathway}$m_{\mathrm{linear\_pathway}}$

Turn MathJax on

equation(11)

$=\log \left[{\left\{1+\sum _i\frac{{\epsilon }_i{A}_i^{-1}}{{\sum }_i{\epsilon }_i{A}_i^{-1}}{\mathit{se}}^{\beta \mathrm{\Delta }{G}_i}\right\}}^{-1}-c\sum _i{A}_i{e}^{\beta \mathrm{\Delta }{G}_i}{(1+e^{\beta \mathrm{\Delta }{G}_i})}^{-1}\right]=-\sum _i\left\{\frac{{\epsilon }_i{A}_i^{-1}}{\left({\sum }_i{\epsilon }_i{A}_i^{-1}\right)}+{\mathit{cA}}_i\right\}{e}^{\beta \mathrm{\Delta }{G}_i}+O\left({\left(\sum _i\left\{\frac{{\epsilon }_i{A}_i^{-1}}{\left({\sum }_i{\epsilon }_i{A}_i^{-1}\right)}+{\mathit{cA}}_i\right\}{e}^{\beta \mathrm{\Delta }{G}_i}\right)}^2\right)$

Turn MathJax on

${\sum }_{i=1}^{10}{\mathit{cA}}_i\exp \left(\beta \mathrm{\Delta }{G}_i\right)<0.03$

equation(12)

$m_{\mathrm{linear\_pathway}}\equiv -\sum _i\left\{\frac{{\epsilon }_i{A}_i^{-1}}{{\sum }_i{\epsilon }_i{A}_i^{-1}}+{\mathit{cA}}_i\right\}{e}^{\beta \mathrm{\Delta }{G}_i}$

Turn MathJax on

2.3. Other formulations of protein fitness

Also, the following simple definition for fitness to maintain protein stability was used (Dasmeh et al., 2014):

equation(13)

w∝f_native$w\propto f_{\mathrm{native}}$

Turn MathJax on

equation(14)

m=−e^βΔG+O(e^2βΔG)+constant$m=-e^{\beta \mathrm{\Delta }G}+O\left(e^{2\beta \mathrm{\Delta }G}\right)+\mathrm{constant}$

Turn MathJax on

In addition, Eq. (3) for functional loss was employed with γ_i⇒cA_i${\gamma }_i\Rightarrow {\mathit{cA}}_i$ to represent toxicity of misfolded proteins in (Serohijos et al., 2012 and Serohijos et al., 2013).

2.4. A generic form of protein fitness

Thus, all expressions above for Malthusian fitness of protein can be well approximated by the following expression, because of exp(βΔG)⪡1$\exp \left(\beta \mathrm{\Delta }G\right)⪡1$ in the typical range of folding free energies shown in Fig. 1:

equation(15)

$m\equiv -\sum _i{\kappa }_i{e}^{\beta \mathrm{\Delta }{G}_i}\mathrm{with}{\kappa }_i\ge 0$

Turn MathJax on

equation(16)

κ_i=cA_i+γ_i≥0${\kappa }_i={\mathit{cA}}_i+{\gamma }_i\ge 0$

Turn MathJax on

The selective advantage of a mutant, in which each protein is destabilized by ΔΔG_i$\mathrm{\Delta }\mathrm{\Delta }{G}_i$, to the wild type can be represented by

equation(17)

$s\equiv m^{\mathrm{mutant}}-m^{\mathrm{wildtype}}=\sum _i{s}_i$

Turn MathJax on

equation(18)

$s_i={\kappa }_i{e}^{\beta \mathrm{\Delta }{G}_i}(1-e^{\beta \mathrm{\Delta }\mathrm{\Delta }{G}_i})\mathrm{with}{\kappa }_i\ge 0$

Turn MathJax on

3. Results

3.1. Protein stability and fitness

Here, we consider the evolution of a single protein-coding gene in which the selective advantage of mutant proteins in Malthusian parameters is assumed to be

equation(19)

$s=\kappa e^{\beta \mathrm{\Delta }G}\left(1-e^{\beta \mathrm{\Delta }\mathrm{\Delta }G}\right)\mathrm{with}\kappa \ge 0$

Turn MathJax on

equation(20)

s≤κe^βΔG$s\le \kappa e^{\beta \mathrm{\Delta }G}$

Turn MathJax on

$T=298°$

equation(21)

κ=cA+γ$\kappa =\mathit{cA}+\gamma$

Turn MathJax on

$\kappa \exp \left(\beta \mathrm{\Delta }G\right)$

$0\le \log 4N_e\kappa \le 20$

Based on measurements of stability changes due to single amino acid substitutions in proteins, which are collected in the ProTherm database (Kumar et al., 2006), Serohijos et al. (2012) reported that the distribution of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ is approximately a Gaussian distribution with mean=1 kcal/mol and standard deviation=1.7 kcal/mol. In addition, it was shown (Serohijos et al., 2012) that the mean of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ is negatively proportional to ΔG$\mathrm{\Delta }G$, and that this dependence of the mean of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ on ΔG$\mathrm{\Delta }G$ is not large but still important to cause the observed negative correlation between protein abundance and evolutionary rate. On the other hand, (Tokuriki et al., 2007) computationally predicted ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ for all possible single amino acid substitutions in 21 different globular, single domain proteins, and showed that the predicted distributions of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ were strikingly similar despite a range of protein sizes and folds and largely follow a bi-Gaussian distribution: one of the two Gaussian distributions results from substitutions on protein surfaces and is a narrow distribution with a mildly destabilizing mean ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$, whereas the other due to substitutions in protein cores is a wider distribution with a stronger destabilizing mean (Tokuriki et al., 2007).

Here, according to Tokuriki et al. (2007), the distribution of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ due to single amino acid substitutions is approximated as a bi-Gaussian function with the dependence of mean ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ on ΔG$\mathrm{\Delta }G$, in order to examine the effects of structural constraint on evolutionary rate. The probability density function (PDF) of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$, p(ΔΔG)$p\left(\mathrm{\Delta }\mathrm{\Delta }G\right)$, for nonsynonymous substitutions is assumed to be

equation(22)

p(ΔΔG)=θN(μ_s,σ_s)+(1−θ)N(μ_c,σ_c)$p\left(\mathrm{\Delta }\mathrm{\Delta }G\right)=\theta \mathcal{N}({\mu }_s,{\sigma }_s)+(1-\theta )\mathcal{N}({\mu }_c,{\sigma }_c)$

Turn MathJax on

equation(23)

${\mu }_s=-0.14\mathrm{\Delta }G-0.17,{\sigma }_s=0.90$

Turn MathJax on

equation(24)

${\mu }_c=-0.14\mathrm{\Delta }G+1.23,{\sigma }_c=1.93$

Turn MathJax on

1.

Equations (23) and (24) satisfy μ_s(ΔG₀)=0.56${\mu }_s\left(\mathrm{\Delta }{G}_0\right)=0.56$ and μ_c(ΔG₀)=1.96${\mu }_c\left(\mathrm{\Delta }{G}_0\right)=1.96$, which were estimated for single nucleotide substitutions in Tokuriki et al. (2007), at a certain value (ΔG₀$\mathrm{\Delta }{G}_0$) of ΔG$\mathrm{\Delta }G$.

2.

The total mean of the two Gaussian functions agrees with the regression line, μ=−0.14ΔG+0.49$\mu =-0.14\mathrm{\Delta }G+0.49$. The value of θ is taken to be 0.53, which is equal to the average of θ over proteins used in Tokuriki et al. (2007).

$\log 4N_e\kappa$

$\theta =1.27-0.33·{\log }_{10}\left(\mathrm{protein}\mathrm{length}\right)\mathrm{for}50\le \mathrm{length}\le 330$

The dependence of the PDF, p(ΔΔG)$p\left(\mathrm{\Delta }\mathrm{\Delta }G\right)$, of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ on θ is shown in Fig. 3. Also, the PDF of selective advantage, p(4N_es)=−p(ΔΔG)dΔΔG/d4N_es$p\left(4N_{e}s\right)=-p\left(\mathrm{\Delta }\mathrm{\Delta }G\right)d\mathrm{\Delta }\mathrm{\Delta }G/d4N_{e}s$, is shown in Fig. S.4 to have a peak at a small, positive value of selective advantage, which moves toward more positive values as θ   and/or 4N_eκ$4N_e\kappa$ increase.

3.2. Equilibrium state of protein stability in protein evolution

The fixation probability u for a mutant gene with selective advantage s and gene frequency q in a duploid system of effective population size N_e was given as a function of 4N_es$4N_{e}s$ and q by ( Crow and Kimura, 1970)

equation(25)

$u\left(4N_{e}s\right)=\frac{1-e^{-4N_e\mathit{sq}}}{1-e^{-4N_{e}s}}$

Turn MathJax on

equation(26)

$\frac{K_a}{K_s}=\frac{u\left(4N_{e}s\right)}{u\left(0\right)}=\frac{u\left(4N_{e}s\right)}{q}\mathrm{with}q=\frac{1}{2N}$

Turn MathJax on

equation(27)

$\simeq \frac{4N_{e}s}{1-e^{-4N_{e}s}}\mathrm{for}\frac{\left|4N_e\mathit{sq}\right|}{2}⪡1$

Turn MathJax on

The PDF of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ of fixed mutant genes, p(ΔΔG_fixed)$p\left(\mathrm{\Delta }\mathrm{\Delta }{G}_{\mathrm{fixed}}\right)$, is

equation(28)

$p\left(\mathrm{\Delta }\mathrm{\Delta }{G}_{\mathrm{fixed}}\right)\equiv p\left(\mathrm{\Delta }\mathrm{\Delta }G\right)\frac{u\left(4N_{e}s\right)}{〈u〉}$

Turn MathJax on

equation(29)

$〈u〉\equiv {\int }_{-\infty }^{\infty }u\left(4N_{e}s\right)p\left(\mathrm{\Delta }\mathrm{\Delta }G\right)d\mathrm{\Delta }\mathrm{\Delta }G$

Turn MathJax on

${〈\mathrm{\Delta }\mathrm{\Delta }G〉}_{\mathrm{fixed}}\equiv {\int }_{-\infty }^{\infty }\mathrm{\Delta }\mathrm{\Delta }Gp\left(\mathrm{\Delta }\mathrm{\Delta }{G}_{\mathrm{fixed}}\right)d\mathrm{\Delta }\mathrm{\Delta }G$

Fig. 4 shows the average of the ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ over fixed mutant genes, 〈ΔΔG〉_fixed${〈\mathrm{\Delta }\mathrm{\Delta }G〉}_{\mathrm{fixed}}$, to monotonically decrease with ΔG$\mathrm{\Delta }G$, indicating that the temporal process of ΔG$\mathrm{\Delta }G$ is stable at 〈ΔΔG〉_fixed(ΔG_e)=0${〈\mathrm{\Delta }\mathrm{\Delta }G〉}_{\mathrm{fixed}}\left(\mathrm{\Delta }{G}_e\right)=0$ due to the balance between random drift on destabilizing mutations and positive selection on stabilizing mutations; ΔG_e$\mathrm{\Delta }{G}_e$ is the folding free energy at the equilibrium state. If a wild-type protein becomes less stable than the equilibrium, ΔG>ΔG_e$\mathrm{\Delta }G>\mathrm{\Delta }{G}_e$, more stabilizing mutants will fix due to primarily positive selection and secondarily random drift, because stabilizing mutants will increase due to negative shifts of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ and also the effect of stability change on selective advantage will be more amplified; see Eqs. (23) and (24) for the dependence of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ on ΔG$\mathrm{\Delta }G$, and Eq. (19) for the fitness of stability change. As shown in Fig. S.6, the probability of K_a/K_s>1.0$K_a/K_s>1.0$, that is, positive selection, significantly increases as ΔG$\mathrm{\Delta }G$ becomes more positive than the equilibrium stability ΔG_e$\mathrm{\Delta }{G}_e$. On the other hand, if a wild-type protein becomes more stable than the equilibrium, ΔG<ΔG_e$\mathrm{\Delta }G<\mathrm{\Delta }{G}_e$, more destabilizing mutants will fix due to random drift, because destabilizing mutants will increase due to positive shifts of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ and also more destabilizing mutants become nearly neutral due to the less-amplified effect of stability change on selective advantage. As shown later, the PDF of K_a/K_s$K_a/K_s$ in the vicinity of equilibrium confirms this mechanism for maintaining protein stability at equilibrium.

It was claimed (Serohijos et al., 2012 and Serohijos et al., 2013) that the equilibrium point would correspond to the minimum of the average fixation probability. However, in Fig. 4 for $\log 4N_e\kappa =7.550$ and θ=0.53$\theta =0.53$, the average 〈s〉_fixed${〈s〉}_{\mathrm{fixed}}$ of selective advantage in fixed mutants has a minimum at ΔG=−5.50$\mathrm{\Delta }G=-5.50$ kcal/mol and changes its sign at ΔG=−4.58$\mathrm{\Delta }G=-4.58$ kcal/mol, where the average 〈K_a/K_s〉=〈u〉/q$〈K_a/K_s〉=〈u〉/q$ has a minimum and which is more positive than the equilibrium stability ΔG_e=−5.24$\mathrm{\Delta }{G}_e=-5.24$ kcal/mol. In other words, Fig. 4 and Fig. S.16 show that the values of ΔG$\mathrm{\Delta }G$ at 〈ΔΔG〉_fixed=0${〈\mathrm{\Delta }\mathrm{\Delta }G〉}_{\mathrm{fixed}}=0$ and at the minimum of 〈K_a/K_s〉$〈K_a/K_s〉$ may be close but differ from each other, and indicate that the value of ΔG$\mathrm{\Delta }G$ corresponding to the minimum of 〈K_a/K_s〉$〈K_a/K_s〉$ is not a good approximation for the equilibrium stability, because 〈K_a/K_s〉$〈K_a/K_s〉$ gently changes in the vicinity of the equilibrium stability as shown in Fig. S.16.

3.3. Equilibrium stability, ΔG_e$\mathrm{\Delta }{G}_e$

The equilibrium value, ΔG_e$\mathrm{\Delta }{G}_e$, of ΔG$\mathrm{\Delta }G$ that satisfies 〈ΔΔG_fixed〉$〈\mathrm{\Delta }\mathrm{\Delta }{G}_{\mathrm{fixed}}〉$=0 in fixed mutants depends directly on θ   and indirectly on 4N_eκ$4N_e\kappa$ through fixation probability; see Eqs. (19) and (25). As shown in Fig. 5, ΔG_e$\mathrm{\Delta }{G}_e$ depends weakly on θ  . On the other hand, ΔG_e$\mathrm{\Delta }{G}_e$ depends more strongly on and is almost negatively proportional to $\log 4N_e\kappa$, as also shown in real proteins (Serohijos et al., 2013). If the dependence of the means, μ_s and μ_c in Eqs. (23) and (24), of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ in all mutants on ΔG$\mathrm{\Delta }G$ could be neglected, $4N_e\kappa \exp \left(\beta \mathrm{\Delta }G\right)$ could be regarded as a single parameter, and so $4N_e\kappa \exp \left(\beta \mathrm{\Delta }{G}_e\right)$ would be constant, irrespective of 4N_eκ$4N_e\kappa$. Thus, the dependence of $\log 4N_e\kappa \exp \left(\beta \mathrm{\Delta }{G}_e\right)$ on $\log 4N_e\kappa$ shown in Fig. 5 is caused solely by the linear dependence of the means μ_s and μ_c of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ on ΔG$\mathrm{\Delta }G$ (Serohijos et al., 2012). It is interesting to know that as $\log 4N_e\kappa$ varies from 0 to 20,ΔG_e$20,\mathrm{\Delta }{G}_e$ changes from $-1.5\mathrm{to}-12.5\mathrm{kcal}/\mathrm{mol}$, the range of which is consistent with experimental values of protein folding free energies shown in Fig. 1.

3.4. K_a/K_s$K_a/K_s$ at equilibrium, ΔG=ΔG_e$\mathrm{\Delta }G=\mathrm{\Delta }{G}_e$

Equations (23) and (24) indicate that the distribution of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ shifts toward the positive direction as ΔG$\mathrm{\Delta }G$ becomes more negative. Hence, increasing 4N_eκ$4N_e\kappa$ that makes ΔG_e$\mathrm{\Delta }{G}_e$ more negative results in positive shifts of the distribution of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$, which increase destabilizing mutations. In addition, as indicated by Eq. (19), the upper bound of $4N_{e}s,4N_e\kappa exp\left(\beta \Delta G_e\right)$, scales the effect of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ on protein fitness. The larger $4N_e\kappa \exp (\beta \mathrm{\Delta }{G}_e)$ is, the larger the effect of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ on selective advantage becomes. Thus, the increase of 4N_eκexp(βΔG_e)$4N_e\kappa \exp \left(\beta \mathrm{\Delta }{G}_e\right)$ caused by the increase of κ$\kappa$ and/or N_e$N_e$ increases both destabilizing mutations and their fitness costs, and results in slow evolutionary rates for proteins with large κ$\kappa$ and/or N_e$N_e$. In other words, highly expressed and indispensable genes, and genes with a large effective population size must evolve slowly. On the other hand, the decrease of θ$\theta$, that is, the increase of highly constrained residues directly shifts the average of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ in all arising mutants toward the positive direction, and causes slow evolutionary rates.

The average of K_a/K_s$K_a/K_s$ over all mutants, which can be observed as the ratio of average nonsynonymous substitution rate per nonsynonymous site to average synonymous substitution rate per synonymous site, and also that over fixed mutants only are shown in Fig. 6. At any value of θ,〈K_a/K_s〉$\theta ,〈K_a/K_s〉$ decreases as $\log 4N_e\kappa$ increases, explaining the observed relationship that highly expressed and indispensable genes evolve slowly (Drummond et al., 2005, Drummond and Wilke, 2008 and Serohijos et al., 2012). Likewise, at any value of $\log 4N_e\kappa$, 〈K_a/K_s〉$〈K_a/K_s〉$ increases as θ$\theta$ increases. In other words, the more structurally constrained a protein is, the more slowly it evolves. The effect of protein abundance/indispensability on evolutionary rate is more remarkable for less constrained proteins and the effect of structural constraint is more remarkable for less abundant, less essential proteins.

The average of K_a/K_s$K_a/K_s$ over all mutants, 〈K_a/K_s〉$〈K_a/K_s〉$, is less than 1.0 on the whole domain shown in the figure, indicating that the average of K_a/K_s$K_a/K_s$ over a long time interval and over many sites should not show any positive selection. Even the average of K_a/K_s$K_a/K_s$ over fixed mutants is less than 1.0, and falls into a narrow range of 0.97–0.85, which is much narrower than a range of 0.96–0.15 for that over all mutants; the average of K_a/K_s$K_a/K_s$ over fixed mutants is equal to 〈(K_a/K_s)²〉/〈K_a/K_s〉$〈{(K_a/K_s)}^2〉/〈K_a/K_s〉$, and as a matter of course must be equal to or larger than the averages of K_a/K_s$K_a/K_s$ over all mutants. However, the average of K_a/K_s$K_a/K_s$ over a short time interval and over a small number of sites may exhibit values larger than one. In Fig. 7, the PDFs of K_a/K_s$K_a/K_s$ for all mutants and also for fixed mutants only are shown; p(K_a/K_s)=p(4N_es)d(4N_es)/d(K_a/K_s)$p(K_a/K_s)=p\left(4N_{e}s\right)d\left(4N_{e}s\right)/d(K_a/K_s)$. A significant fraction of fixed mutants fix with K_a/K_s>1$K_a/K_s>1$.

Arbitrarily, the value of K_a/K_s$K_a/K_s$ is categorized into four classes; negative, slightly negative, nearly neutral, and positive selection categories whose K_a/K_s$K_a/K_s$ are within the range of K_a/K_s≤0.5,0.5<K_a/K_s≤0.95,0.95<K_a/K_s≤1.05$K_a/K_s\le 0.5,0.5<K_a/K_s\le 0.95,0.95<K_a/K_s\le 1.05$, and 1.05<K_a/K_s$1.05<K_a/K_s$, respectively. Then, the probabilities of each selection category in all mutants and in fixed mutants are calculated and shown in Fig. S.10 and Fig. 8, respectively. At the largest abundance ($\log 4N_e\kappa =20$) most arising mutations are negative mutations whose K_a/K_s$K_a/K_s$ are less than 0.5. This is reasonable, because at this condition the wild-type protein is very stable with low equilibrium values ΔG_e$\mathrm{\Delta }{G}_e$ as shown in Fig. 5, and therefore most mutations destabilize the wild-type and tend to be removed from population. Most fixed mutants are positive mutants or slightly negative mutants fixed by random drift. Nearly neutral mutants are less than 3% of all mutants, and less than 15% of fixed mutants.

On the other hand, at the other extreme of $\log 4N_e\kappa <2$, there are no mutations of the positive selection category, this is because the upper bound of K_a/K_s$K_a/K_s$, which corresponds to the upper bound ($4N_e\kappa \exp \left(\beta \mathrm{\Delta }G\right)$) of 4N_es$4N_{e}s$, at the equilibrium stability ΔG_e$\mathrm{\Delta }{G}_e$ becomes less than 1.05 that is the lower bound for the positive selection category; see Eq. (19). The significant amount of mutations become nearly neutral. As θ$\theta$ changes from 1 to 0, that is, structural constraints increase, the proportion of nearly neutral mutations changes from 0.75(0.56) to 0.31(0.22), and instead negative mutations increase and most of them are removed from population. Thus, the selection mechanism for fixing stabilizing mutants in little expressed, non-essential genes ($\log 4N_e\kappa <2$) is not positive selection but nearly neutral selection, that is, random drift.

The probability of each selection category in fixed mutants depends strongly on 4N_eκ$4N_e\kappa$, but much less on θ$\theta$. Current common understanding is that amino acid substitutions in protein evolution are either neutral (Kimura, 1968) or lethal, at most slightly deleterious (Ohta, 1973) or lethal, unless functional selection operates and functional changes occur. On the contrary, nearly neutral fixations are predominant only in proteins with $\log 4N_e\kappa <2$ or ΔG_e>−2.5$\mathrm{\Delta }{G}_e>-2.5$ kcal/mol, and positive selection is significant in the other proteins. On the other hand, slightly negative selection is always significant. An interesting result is that the effects of structural constraint on K_a/K_s$K_a/K_s$ are the most remarkable in proteins with $\log 4N_e\kappa <2$ or instead ΔG_e>−2.5$\mathrm{\Delta }{G}_e>-2.5$ kcal/mol in which nearly neutral fixations are predominant.

3.5. K_a/K_s$K_a/K_s$ in the vicinity of equilibrium

In Fig. 4, the 〈ΔΔG〉_fixed±${〈\mathrm{\Delta }\mathrm{\Delta }G〉}_{\mathrm{fixed}}\pm$ standard deviation of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ of fixed mutants are also drawn. The standard deviation of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ of fixed mutants is equal to 0.84 kcal/mol at the equilibrium, ΔG_e$\mathrm{\Delta }{G}_e$, indicating that protein stability ΔG$\mathrm{\Delta }G$ fluctuates more or less within ΔG_e±0.84$\mathrm{\Delta }{G}_e\pm 0.84$ kcal/mol instantaneously. Such a deviation from the equilibrium must be canceled by compensatory substitutions that consecutively occur, otherwise the protein stability would far depart from its equilibrium point.

In Fig. 9 and Fig. 10 and Figs. S.12 and S.14, the probabilities of each selection category in fixed mutants and in all arising mutants are shown as a function of ΔG$\mathrm{\Delta }G$ and 4N_eκ$4N_e\kappa$ or θ$\theta$, respectively. The range of ΔG$\mathrm{\Delta }G$ around ΔG_e$\mathrm{\Delta }{G}_e$ shown by a blue line on the surface grid is within two times of the standard deviation of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ over fixed mutants at ΔG_e$\mathrm{\Delta }{G}_e$.

As indicated by Eqs. (23) and (24), it is shown in Figs. S.12 and S14 that stabilizing mutations increase due to negative shifts of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ as the wild type becomes less stable than the equilibrium, ΔG>ΔG_e$\mathrm{\Delta }G>\mathrm{\Delta }{G}_e$, and that destabilizing mutations increase due to positive shifts of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ as ΔG<ΔG_e$\mathrm{\Delta }G<\mathrm{\Delta }{G}_e$. In addition, as indicated by Eq. (19), it is shown in Fig. 9 and Fig. 10 that positive selection on stabilizing mutants is more amplified as ΔG>ΔG_e$\mathrm{\Delta }G>\mathrm{\Delta }{G}_e$, and that more destabilizing mutations become nearly neutral due to the less-amplified effect of stability change on selective advantage as ΔG<ΔG_e$\mathrm{\Delta }G<\mathrm{\Delta }{G}_e$. This is a mechanism of maintaining protein stability at equilibrium.

3.6. Lower bound of K_a/K_s$K_a/K_s$ for adaptive substitutions on protein function

The observed value of K_a/K_s>1$K_a/K_s>1$ is often used to indicate functional selection. The averages of K_a/K_s$K_a/K_s$ over all mutants and even over fixed mutants are less than 1 as shown in Fig. 6. Therefore, the average of K_a/K_s$K_a/K_s$ over long time or many sites does not indicate positive selection. However, the probability of K_a/K_s>1$K_a/K_s>1$ is not negligible as shown in Fig. 7 and Fig. 8. Then, a question is how large K_a/K_s$K_a/K_s$ due to selection on protein stability can be.

The distribution of K_a/K_s$K_a/K_s$ significantly changes with ΔG$\mathrm{\Delta }G$, as shown in Fig. S.6 and Fig. 11. It may be appropriate to see the average of K_a/K_s,〈K_a/K_s〉_fixed$K_a/K_s,{〈K_a/K_s〉}_{\mathrm{fixed}}$, in mutants fixed at ΔG>ΔG_e$\mathrm{\Delta }G>\mathrm{\Delta }{G}_e$, because the upper bound of K_a/K_s$K_a/K_s$ becomes larger for ΔG>ΔG_e$\mathrm{\Delta }G>\mathrm{\Delta }{G}_e$ than at the equilibrium, and also positive mutants must fix to improve the protein stability of the wild type. Fig. 11 shows that 〈K_a/K_s〉_fixed${〈K_a/K_s〉}_{\mathrm{fixed}}$ can be very large for proteins with low equilibrium stabilities (large 4N_eκ$4N_e\kappa$ and small θ$\theta$), although 〈K_a/K_s〉_fixed ~ 1in ΔG < ΔG_e in which nearly neutral and slightly negative selections are predominant 1.7(1.2) at $\mathrm{\Delta }{G}_e+\mathrm{\Delta }\mathrm{\Delta }{G}_{\mathrm{fixed}}^{\mathrm{sd}}$ and 6.1(5.6)$6.1\left(5.6\right)$ at $\mathrm{\Delta }{G}_e+2·\mathrm{\Delta }\mathrm{\Delta }{G}_{\mathrm{fixed}}^{\mathrm{sd}}$ for $\log 4N_e\kappa =20\left(\theta =0.0\right)$, where $\mathrm{\Delta }\mathrm{\Delta }{G}_{\mathrm{fixed}}^{\mathrm{sd}}$ means the standard deviation of ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ in fixed mutants at ΔG_e$\mathrm{\Delta }{G}_e$. The 85% of fixed mutants have ΔΔG$\mathrm{\Delta }\mathrm{\Delta }G$ within the standard deviation. Therefore, a lower bound for adaptive substitutions may be taken to be about 1.7, which is almost equal to the upper bound of K_a/K_s$K_a/K_s$ at the equilibrium for $\log 4N_e\kappa =20$ and θ=1$\theta =1$; see Fig. S.17. However, as shown in Fig. S.17, the more genes are expressed and/or the stronger structural constraints are, the larger the upper bound of K_a/K_s$K_a/K_s$ at the equilibrium is. Judging of adaptive changes may need not only K_a/K_s>1$K_a/K_s>1$ but also other supporting evidences; such that substitutions are localized at specific sites.

3.7. Dependences of protein stability (βΔG_e)$\left(\beta \mathrm{\Delta }{G}_e\right)$ and evolutionary rate (〈K_a/K_s〉)$(〈K_a/K_s〉)$ on growth temperature

It is natural that the folding free energies, ΔG_e$\mathrm{\Delta }{G}_e$, of proteins in organisms growing at higher temperatures must be lower than those at the normal temperature, in order to attain the same stabilities and fitnesses as in the normal temperature. Equations (2) and (15) indicate that the same stability and fitness will be attained if βΔG_e$\beta \mathrm{\Delta }{G}_e$ is constant. It means that it is sufficient for ΔG_e$\mathrm{\Delta }{G}_e$ at the $100°\mathrm{C}$ to decrease 373/298=1.25$373/298=1.25$ times of that at the normal temperature $(25°\mathrm{C})$. Is this a figure expected for folding free energies of thermophilic proteins at high growth temperature? It is not enough data of ΔG$\mathrm{\Delta }G$ at high temperature in the ProTherm (Kumar et al., 2006) to answer this question; $\mathrm{\Delta }G(T=75°\mathrm{C})=10.76$ kcal/mol for oxidized and 4.3 for reduced CuA domain of cytochrome oxidase from Thermus thermophilus (Wittung-Stafshede et al., 1998) and $\mathrm{\Delta }G(T=60°\mathrm{C})=13.01$ kcal/mol for pyrrolidone carboxyl peptidase from Pyrococcus (Ogasahara et al., 1998). The present model indicates that βΔG_e$\beta \mathrm{\Delta }{G}_e$ slightly increases as growth temperature increases.