Selection originating from protein stability/foldability: Relationships between protein folding free energy, sequence ensemble, and fitness

doi:10.1016/j.jtbi.2017.08.018

Highlights

•: A Boltzmann distribution with protein fitness is derived.
•: Relationships between folding free energy, inverse statistical potential and fitness.
•: Selective temperature, glass transition temperature and folding free energy are estimated.
•: Relationship between selective temperature and substitution rate (K_a/K_s).
•: Protein stability/foldability is kept in a balance of positive selection and random drift.

Abstract

Assuming that mutation and fixation processes are reversible Markov processes, we prove that the equilibrium ensemble of sequences obeys a Boltzmann distribution with $\exp (4 N_{e} m (1 - 1 / (2 N))),$ $\exp (4 N_{e} m (1 - 1 / (2 N))),$ where m is Malthusian fitness and N_e and N are effective and actual population sizes. On the other hand, the probability distribution of sequences with maximum entropy that satisfies a given amino acid composition at each site and a given pairwise amino acid frequency at each site pair is a Boltzmann distribution with $\exp (- ψ_{N}),$ $\exp (- ψ_{N}),$ where the evolutionary statistical energy ψ_N is represented as the sum of one body (h) (compositional) and pairwise (J) (covariational) interactions over all sites and site pairs. A protein folding theory based on the random energy model (REM) indicates that the equilibrium ensemble of natural protein sequences is well represented by a canonical ensemble characterized by $\exp (- Δ G_{N D} / k_{B} T_{s})$ $\exp (- Δ G_{N D} / k_{B} T_{s})$ or by $\exp (- G_{N} / k_{B} T_{s})$ $\exp (- G_{N} / k_{B} T_{s})$ if an amino acid composition is kept constant, where $Δ G_{N D} \equiv G_{N} - G_{D},$ $Δ G_{N D} \equiv G_{N} - G_{D},$ G_N and G_D are the native and denatured free energies, and T_s is the effective temperature representing the strength of selection pressure. Thus, $4 N_{e} m (1 - 1 / (2 N)),$ $4 N_{e} m (1 - 1 / (2 N)),$ $- Δ ψ_{N D} (\equiv - ψ_{N} + ψ_{D}),$ $- Δ ψ_{N D} (\equiv - ψ_{N} + ψ_{D}),$ and $- Δ G_{N D} / k_{B} T_{s}$ $- Δ G_{N D} / k_{B} T_{s}$ must be equivalent to each other. With h and J estimated by the DCA program, the changes (Δψ_N) of ψ_N due to single nucleotide nonsynonymous substitutions are analyzed. The results indicate that the standard deviation of $Δ G_{N} (= k_{B} T_{s} Δ ψ_{N})$ $Δ G_{N} (= k_{B} T_{s} Δ ψ_{N})$ is approximately constant irrespective of protein families, and therefore can be used to estimate the relative value of T_s. Glass transition temperature T_g and ΔG_ND are estimated from estimated T_s and experimental melting temperature (T_m) for 14 protein domains. The estimates of ΔG_ND agree with their experimental values for 5 proteins, and those of T_s and T_g are all within a reasonable range. In addition, approximating the probability density function (PDF) of Δψ_N by a log-normal distribution, PDFs of Δψ_N and K_a/K_s, which is the ratio of nonsynonymous to synonymous substitution rate per site, in all and in fixed mutants are estimated. The equilibrium values of ψ_N, at which the average of Δψ in fixed mutants is equal to zero, well match ψ_N averaged over homologous sequences, confirming that the present methods for a fixation process of mutations and for the equilibrium ensemble of ψ_N give a consistent result with each other. The PDFs of K_a/K_s at equilibrium confirm that T_s negatively correlates with the amino acid substitution rate (the mean of K_a/K_s) of protein. Interestingly, stabilizing mutations are significantly fixed by positive selection, and balance with destabilizing mutations fixed by random drift, although most of them are removed from population. Supporting the nearly neutral theory, neutral selection is not significant even in fixed mutants.

1. Introduction

Natural proteins can fold their sequences into unique structures. Protein’s stability and foldability result from natural selection and are not typical characteristics of random polymers (Bryngelson and Wolynes, 1987; Pande et al., 1997; Ramanathan and Shakhnovich, 1994; Shakhnovich and Gutin, 1993a; 1993b). Natural selection maintains protein’s stability and foldability over evolutionary timescales. On the basis of the random energy model (REM) for protein folding, it was discussed (Ramanathan and Shakhnovich, 1994; Shakhnovich and Gutin, 1993a; 1993b) that the equilibrium ensemble of natural protein sequences in sequence space is well represented by a canonical ensemble characterized by a Boltzmann factor $\exp (- Δ G_{N D} (σ) / k_{B} T_{s}),$ $\exp (- Δ G_{N D} (σ) / k_{B} T_{s}),$ where $Δ G_{N D} (σ) (\equiv G_{N} (σ) - G_{D} (σ))$ $Δ G_{N D} (σ) (\equiv G_{N} (σ) - G_{D} (σ))$ is the folding free energy of sequence σ, G_N and G_D are the free energies of the native and denatured states, k_B is the Boltzmann constant, and T_s is the effective temperature representing the strength of selection pressure and must satisfy T_s < T_g < T_m for natural proteins to fold into unique native structures; T_g is glass transition temperature and T_m is melting temperature. The REM also indicates that the free energy of denatured conformations (G_D) is a function of amino acid frequencies only and does not depend on amino acid order, and therefore the Boltzmann factor will be taken as $\exp (- G_{N} (σ) / k_{B} T_{s}),$ $\exp (- G_{N} (σ) / k_{B} T_{s}),$ if amino acid frequencies are kept constant. It was shown by lattice Monte Carlo simulations (Shakhnovich, 1994) that lattice protein sequences selected with this Boltzmann factor were not trapped by competing structures but could fold into unique native structures. Selective temperatures were also estimated (Dokholyan and Shakhnovich, 2001) for actual proteins to yield good correlations of sequence entropy between actual protein families and sequences designed with this type of Boltzmann factor.

On the other hand, the maximum entropy principle insists that the probability distribution of sequences in sequence space, which satisfies constraints on amino acid compositions at all sites and on amino acid pairwise frequencies for all site pairs, is a Boltzmann distribution with the Boltzmann factor $\exp (- ψ_{N} (σ)),$ $\exp (- ψ_{N} (σ)),$ where the total interaction ψ_N(σ) of a sequence σ is represented as the sum of one-body (h) (compositional) and pairwise (J) (covariational) interactions between residues in the sequence; ψ_N(σ) is called the evolutionary statistical energy by Hopf et al. (2017). The inverse statistical potentials, the one-body (h) and pairwise (J) interactions, that satisfy those constraints for homologous sequences have been estimated (Ekeberg et al., 2014; 2013; Marks et al., 2011; Morcos et al., 2011) as one of inverse Potts problems and successfully employed to predict contacting residue pairs in protein structures (Ekeberg et al., 2014; 2013; Hopf et al., 2012; Marks et al., 2011; Miyazawa, 2013; Morcos et al., 2011; Sułkowska et al., 2012). Morcos et al. (2014) noticed that the ψ_N in the Boltzmann factor is the dimensionless energy corresponding to G_N/k_BT_s, and estimated selective temperatures (T_s) for several protein families by comparing the difference (Δψ_ND) of ψ between the native and the molten globule states with folding free energies (ΔG_ND) estimated with associative-memory, water-mediated, structure, and energy model (AWSEM) (Davtyan et al., 2012).

A purpose of the present study is to establish relationships between protein foldability/stability, sequence distribution, and protein fitness. First, we prove that if mutation and fixation processes in protein evolution are reversible Markov processes, the equilibrium ensemble of genes will obey a Boltzmann distribution with the Boltzmann factor $\exp (4 N_{e} m (1 - 1 / (2 N))),$ $\exp (4 N_{e} m (1 - 1 / (2 N))),$ where N_e and N are effective and actual population sizes, and m is the Malthusian fitness of a gene. In other words, correspondences between $- Δ G_{N D} / k_{B} T_{s},$ $- Δ G_{N D} / k_{B} T_{s},$ $- Δ ψ_{N D} (\equiv ψ_{N} - ψ_{D})$ $- Δ ψ_{N D} (\equiv ψ_{N} - ψ_{D})$ and $4 N_{e} m (1 - 1 / (2 N))$ $4 N_{e} m (1 - 1 / (2 N))$ are obtained by equating these three Boltzmann distributions with each other; $ψ_{D} ≃ G_{D} / k_{B} T_{s} + constant$ $ψ_{D} ≃ G_{D} / k_{B} T_{s} + constant$ .

The second purpose is to analyze the effects (Δψ_N) of single amino acid substitutions on the evolutionary statistical energy of a protein, and to estimate from the distribution of Δψ_N the effective temperature of natural selection (T_s) and then glass transition temperature (T_g) and folding free energy (ΔG_ND) of protein. We estimate the one-body (h) and pairwise (J) interactions with the DCA program, which is available at “http://dca.rice.edu/portal/dca/home”, and then analyze the changes (Δψ_N) of the evolutionary statistical energy (ψ_N) of a natural sequence due to single amino acid substitutions caused by single nucleotide changes. The data of Δψ_N due to single nucleotide nonsynonymous substitutions for 14 protein domains show that the standard deviation of Δψ_N over all the substitutions at all sites hardly depends on the evolutionary statistical energy (ψ_N) of each homologous sequence and is nearly constant for each protein family, indicating that the standard deviation of ΔG_N ≃ k_BT_sΔψ_N is nearly constant irrespective of protein families. From this finding, T_s for each protein family has been estimated in relative to T_s for the PDZ family, which is determined by directly comparing $Δ Δ ψ_{N D} (\equiv Δ (ψ_{N} - ψ_{D}) ≃ Δ ψ_{N})$ $Δ Δ ψ_{N D} (\equiv Δ (ψ_{N} - ψ_{D}) ≃ Δ ψ_{N})$ with the experimental values of folding free energy changes, ΔΔG_ND, due to single amino acid substitutions. Also T_g and ΔG_ND for each protein family are estimated on the basis of the REM from the estimated T_s and an experimental melting temperature T_m. The estimates of T_s and T_g are all within a reasonable range, and those of ΔG_ND are well compared with experimental ΔG_ND for 5 protein families. The present method for estimating T_s is simpler than the method (Morcos et al., 2014) using AWSEM, and also is useful for the prediction of ΔG_ND, because the experimental data of ΔG_ND are limited in comparison with T_m, and also experimental conditions such as temperature and pH tend to be different among them. In addition, it has been revealed that Δψ_N averaged over all single nucleotide nonsynonymous substitutions is a linear function of ψ_N/L of each homologous sequence, where L is sequence length; the average of Δψ_N decreases as ψ_N/L increases. This characteristic is required for homologous proteins to stay at the equilibrium state of the native conformational energy G_N ≃ k_BT_sψ_N, and indicates a weak dependency (Miyazawa, 2016; Serohijos et al., 2012) of ΔΔG_ND on ΔG_ND/L of protein across protein families.

The third purpose is to study an amino acid substitution process in protein evolution, which is characterized by the fitness, $m = - Δ ψ_{N D} / (4 N_{e} (1 - 1 / (2 N)))$ $m = - Δ ψ_{N D} / (4 N_{e} (1 - 1 / (2 N)))$ . We employ a monoclonal approximation for mutation and fixation processes of genes, in which protein evolution proceeds with single amino acid substitutions fixed at a time in a population. In this approximation, ψ_N of a protein gene attains the equilibrium, $ψ_{N} = ψ_{N}^{e q},$ $ψ_{N} = ψ_{N}^{e q},$ when the average of Δψ_N( ≃ ΔΔψ_ND) over singe nucleotide nonsynonymous mutations fixed in a population is equal to zero. Approximating the distribution of Δψ_N due to singe nucleotide nonsynonymous mutations by a log-normal distribution, their distribution for fixed mutants is numerically calculated and used to calculate the averages of various quantities and also the probability density functions (PDF) of K_a/K_s in all arising mutants and also in fixed mutants only; K_a/K_s is defined as the ratio of nonsynonymous to synonymous substitution rate per site. There is a good agreement between the time average ( $ψ_{N}^{e q}$ $ψ_{N}^{e q}$ ) and ensemble average (⟨ψ_N⟩_σ), which is equal to the sample average, $\bar{ψ_{N}},$ $\bar{ψ_{N}},$ of ψ_N over homologous sequences, supporting the constancy of the standard deviation of Δψ_N assumed in the monoclonal approximation.

We also study protein evolution at equilibrium, $ψ_{N} = ψ_{N}^{e q}$ $ψ_{N} = ψ_{N}^{e q}$ . The common understanding of protein evolution has been that amino acid substitutions observed in homologous proteins are neutral (Kimura, 1968; 1969; Kimura and Ohta, 1971; 1974) or slightly deleterious (Ohta, 1973; 1992), and random drift is a primary force to fix amino acid substitutions in population. The PDFs of K_a/K_s in all arising mutations and in their fixed mutations are examined to see how significant each of positive, neutral, slightly negative,and negative selections is. Interestingly, stabilizing mutations are significantly fixed in population by positive selection, and balance with destabilizing mutations that are also significantly fixed by random drift, although most negative mutations are removed from population. Contrary to the neutral theory (Kimura, 1968; 1969; Kimura and Ohta, 1971; 1974) and supporting the nearly neutral theory (Ohta, 1973; 1992; 2002), the proportion of neutral selection is not large even in fixed mutants. It is also confirmed that the effective temperature (T_s) of selection negatively correlates with the amino acid substitution rate (K_a/K_s) of protein at equilibrium.

2. Methods

2.1. Knowledge of protein folding

A protein folding theory (Pande et al., 1997; Ramanathan and Shakhnovich, 1994; Shakhnovich and Gutin, 1993a; 1993b), which is based on a random energy model (REM), indicates that the equilibrium ensemble of amino acid sequences, $σ \equiv (σ_{1}, \dots, σ_{L})$ $σ \equiv (σ_{1}, \dots, σ_{L})$ where σ_i is the type of amino acid at site i and L is sequence length, can be well approximated by a canonical ensemble with a Boltzmann factor consisting of the folding free energy, ΔG_ND(σ, T) and an effective temperature T_s representing the strength of selection pressure.

(1)

\begin{matrix} P (σ) & \propto & P^{m u t} (σ) \exp (\frac{- Δ G_{N D} (σ, T)}{k_{B} T_{s}}) \end{matrix}

$\begin{matrix} P (σ) & \propto & P^{m u t} (σ) \exp (\frac{- Δ G_{N D} (σ, T)}{k_{B} T_{s}}) \end{matrix}$

(2)

\begin{matrix} \propto & \exp (\frac{- G_{N} (σ)}{k_{B} T_{s}}) if f (σ) = constant \end{matrix}

$\begin{matrix} \propto & \exp (\frac{- G_{N} (σ)}{k_{B} T_{s}}) if f (σ) = constant \end{matrix}$

(3)

\begin{matrix} Δ G_{N D} (σ, T) & \equiv & G_{N} (σ) - G_{D} (f (σ), T) \end{matrix}

$\begin{matrix} Δ G_{N D} (σ, T) & \equiv & G_{N} (σ) - G_{D} (f (σ), T) \end{matrix}$

where p^mut(σ) is the probability of a sequence (σ) randomly occurring in a mutational process and depends only on the amino acid frequencies f(σ), k_B is the Boltzmann constant, T is a growth temperature, and G_N and G_D are the free energies of the native conformation and denatured state, respectively. Selective temperature T_s quantifies how strong the folding constraints are in protein evolution, and is specific to the protein structure and function. The free energy G_D of the denatured state does not depend on the amino acid order but the amino acid composition, f(σ), in a sequence (Pande et al., 1997; Ramanathan and Shakhnovich, 1994; Shakhnovich and Gutin, 1993a; 1993b). It is reasonable to assume that mutations independently occur between sites, and therefore the equilibrium frequency of a sequence in the mutational process is equal to the product of the equilibrium frequencies over sites;

P^{m u t} (σ) = \prod_{i} p^{m u t} (σ_{i}),

$P^{m u t} (σ) = \prod_{i} p^{m u t} (σ_{i}),$ where p^mut(σ_i) is the equilibrium frequency of σ_i at site i in the mutational process.

The distribution of conformational energies in the denatured state (molten globule state), which consists of conformations as compact as the native conformation, is approximated in the random energy model (REM), particularly the independent interaction model (IIM) (Pande et al., 1997), to be equal to the energy distribution of randomized sequences, which is then approximated by a Gaussian distribution, in the native conformation. That is, the partition function Z for the denatured state is written as follows with the energy density n(E) of conformations that is approximated by a product of a Gaussian probability density and the total number of conformations whose logarithm is proportional to the chain length.

(4)

\begin{matrix} Z & = & \int \exp (\frac{- E}{k_{B} T}) n (E) d E \end{matrix}

$\begin{matrix} Z & = & \int \exp (\frac{- E}{k_{B} T}) n (E) d E \end{matrix}$

(5)

\begin{matrix} n (E) & \approx & \exp (ω L) N (\bar{E} (f (σ)), δ E^{2} (f (σ))) \end{matrix}

$\begin{matrix} n (E) & \approx & \exp (ω L) N (\bar{E} (f (σ)), δ E^{2} (f (σ))) \end{matrix}$

where ω is the conformational entropy per residue in the compact denatured state, and

N (\bar{E} (f (σ)), δ E^{2} (f (σ)))

$N (\bar{E} (f (σ)), δ E^{2} (f (σ)))$ is the Gaussian probability density with mean

\bar{E}

$\bar{E}$ and variance δE², which depend only on the amino acid composition of the protein sequence. The free energy of the denatured state is approximated as follows.

(6)

\begin{matrix} G_{D} (f (σ), T) & \approx & \bar{E} (f (σ)) - \frac{δ E^{2} (f (σ))}{2 k_{B} T} - k_{B} T ω L \end{matrix}

$\begin{matrix} G_{D} (f (σ), T) & \approx & \bar{E} (f (σ)) - \frac{δ E^{2} (f (σ))}{2 k_{B} T} - k_{B} T ω L \end{matrix}$

(7)

\begin{matrix} = & \bar{E} (f (σ)) - δ E^{2} (f (σ)) \frac{ϑ (T / T_{g})}{k_{B} T} \end{matrix}

$\begin{matrix} = & \bar{E} (f (σ)) - δ E^{2} (f (σ)) \frac{ϑ (T / T_{g})}{k_{B} T} \end{matrix}$

(8)

\begin{matrix} ϑ (T / T_{g}) & \equiv & {\begin{matrix} \frac{1}{2} (1 + \frac{T^{2}}{T_{g}^{2}}) & for T > T_{g} \\ \frac{T}{T_{g}} & for T \leq T_{g} \end{matrix} \end{matrix}

$\begin{matrix} ϑ (T / T_{g}) & \equiv & {\begin{matrix} \frac{1}{2} (1 + \frac{T^{2}}{T_{g}^{2}}) & for T > T_{g} \\ \frac{T}{T_{g}} & for T \leq T_{g} \end{matrix} \end{matrix}$

where

\bar{E}

$\bar{E}$ and δE² are estimated as the mean and variance of interaction energies of randomized sequences in the native conformation. T_g is the glass transition temperature of the protein at which entropy becomes zero (Pande et al., 1997; Ramanathan and Shakhnovich, 1994; Shakhnovich and Gutin, 1993a; 1993b);

- \partial G_{D} {/ \partial T |}_{T = T_{g}} = 0

$- \partial G_{D} {/ \partial T |}_{T = T_{g}} = 0$ . The conformational entropy per residue ω in the compact denatured state can be represented with T_g;

ω L = δ E^{2} / (2 {(k_{B} T_{g})}^{2})

$ω L = δ E^{2} / (2 {(k_{B} T_{g})}^{2})$ . Thus, unless T_g < T_m, a protein will be trapped at local minima on a rugged free energy landscape before it can fold into a unique native structure.

2.2. Probability distribution of homologous sequences with the same native fold in sequence space

The probability distribution P(σ) of homologous sequences with the same native fold, $σ = (σ_{1}, \dots, σ_{L})$ $σ = (σ_{1}, \dots, σ_{L})$ where σ_i ∈ {amino acids, deletion}, in sequence space with maximum entropy, which satisfies a given amino acid frequency at each site and a given pairwise amino acid frequency at each site pair, is a Boltzmann distribution (Marks et al., 2011; Morcos et al., 2011).

(9)

\begin{matrix} P (σ) & \propto & \exp (- ψ_{N} (σ)) \end{matrix}

$\begin{matrix} P (σ) & \propto & \exp (- ψ_{N} (σ)) \end{matrix}$

(10)

\begin{matrix} ψ_{N} (σ) & \equiv & - (\sum_{i}^{L} (h_{i} (σ_{i}) + \sum_{j > i} J_{i j} (σ_{i}, σ_{j}))) \end{matrix}

$\begin{matrix} ψ_{N} (σ) & \equiv & - (\sum_{i}^{L} (h_{i} (σ_{i}) + \sum_{j > i} J_{i j} (σ_{i}, σ_{j}))) \end{matrix}$

where h_i and J_ij are one-body (compositional) and two-body (covariational) interactions and must satisfy the following constraints.

(11)

\begin{matrix} \sum_{σ} P (σ) δ_{σ_{i} a_{k}} & = & P_{i} (a_{k}) \end{matrix}

$\begin{matrix} \sum_{σ} P (σ) δ_{σ_{i} a_{k}} & = & P_{i} (a_{k}) \end{matrix}$

(12)

\begin{matrix} \sum_{σ} P (σ) δ_{σ_{i} a_{k}} δ_{σ_{j} a_{l}} & = & P_{i j} (a_{k}, a_{l}) \end{matrix}

$\begin{matrix} \sum_{σ} P (σ) δ_{σ_{i} a_{k}} δ_{σ_{j} a_{l}} & = & P_{i j} (a_{k}, a_{l}) \end{matrix}$

where

δ_{σ_{i} a_{k}}

$δ_{σ_{i} a_{k}}$ is the Kronecker delta, P_i(a_k) is the frequency of amino acid a_k at site i, and P_ij(a_k, a_l) is the frequency of amino acid pair, a_k at i and a_l at j; a_k ∈ {amino acids, deletion}. The pairwise interaction matrix J satisfies

J_{i j} (a_{k}, a_{l}) = J_{j i} (a_{l}, a_{k})

$J_{i j} (a_{k}, a_{l}) = J_{j i} (a_{l}, a_{k})$ and

J_{i i} (a_{k}, a_{l}) = 0

$J_{i i} (a_{k}, a_{l}) = 0$ . Interactions h_i and J_ij can be well estimated from a multiple sequence alignment (MSA) in the mean field approximation (Marks et al., 2011; Morcos et al., 2011), or by maximizing a pseudo-likelihood (Ekeberg et al., 2014; 2013). Because ψ_N(σ) has been estimated under the constraints on amino acid compositions at all sites, only sequences with a given amino acid composition contribute significantly to the partition function, and other sequences may be ignored.

Hence, from Eqs. (2) and (9),

(13)

\begin{matrix} ψ_{N} (σ) & ≃ & G_{N} (σ) / (k_{B} T_{s}) + function of f (σ) \end{matrix}

$\begin{matrix} ψ_{N} (σ) & ≃ & G_{N} (σ) / (k_{B} T_{s}) + function of f (σ) \end{matrix}$

(14)

\begin{matrix} ψ_{D} (f (σ), T) & ≃ & G_{D} (f (σ), T) / (k_{B} T_{s}) + function of f (σ) \end{matrix}

$\begin{matrix} ψ_{D} (f (σ), T) & ≃ & G_{D} (f (σ), T) / (k_{B} T_{s}) + function of f (σ) \end{matrix}$

(15)

\begin{matrix} Δ ψ_{N D} (σ, T) & ≃ & Δ G_{N D} (σ, T) / (k_{B} T_{s}) \end{matrix}

$\begin{matrix} Δ ψ_{N D} (σ, T) & ≃ & Δ G_{N D} (σ, T) / (k_{B} T_{s}) \end{matrix}$

(16)

\begin{matrix} Δ ψ_{N D} (σ, T) & \equiv & ψ_{N} (σ) - ψ_{D} (f (σ), T) \end{matrix}

$\begin{matrix} Δ ψ_{N D} (σ, T) & \equiv & ψ_{N} (σ) - ψ_{D} (f (σ), T) \end{matrix}$

(17)

\begin{matrix} ψ_{D} (f (σ), T) & \approx & \bar{ψ} (f (σ)) - δ ψ^{2} (f (σ)) ϑ (T / T_{g}) T_{s} / T \end{matrix}

$\begin{matrix} ψ_{D} (f (σ), T) & \approx & \bar{ψ} (f (σ)) - δ ψ^{2} (f (σ)) ϑ (T / T_{g}) T_{s} / T \end{matrix}$

(18)

\begin{matrix} ω & = & {(T_{s} / T_{g})}^{2} δ ψ^{2} / (2 L) \end{matrix}

$\begin{matrix} ω & = & {(T_{s} / T_{g})}^{2} δ ψ^{2} / (2 L) \end{matrix}$

where the

\bar{ψ}

$\bar{ψ}$ and δψ² are estimated as the mean and variance of ψ_N over randomized sequences;

\bar{E} ≃ k_{B} T_{s} \bar{ψ}

$\bar{E} ≃ k_{B} T_{s} \bar{ψ}$ and δE² ≃ (k_BT_s)²δψ².

2.3. The equilibrium distribution of sequences in a mutation-fixation process

Here we assume that the mutational process is a reversible Markov process. That is, the mutation rate per gene, M_μν, from sequence $μ \equiv (μ_{1}, \dots, μ_{L})$ $μ \equiv (μ_{1}, \dots, μ_{L})$ to ν satisfies the detailed balance condition

(19)

\begin{matrix} P^{m u t} (μ) M_{μ ν} & = & P^{m u t} (ν) M_{ν μ} \end{matrix}

$\begin{matrix} P^{m u t} (μ) M_{μ ν} & = & P^{m u t} (ν) M_{ν μ} \end{matrix}$

where P^mut(ν) is the equilibrium frequency of sequence ν in a mutational process, M_μν. The mutation rate per population is equal to 2NM_μν for a diploid population, where N is the population size. The substitution rate R_μν from μ to ν is equal to the product of the mutation rate and the fixation probability with which a single mutant gene becomes to fully occupy the population (Crow and Kimura, 1970).

(20)

\begin{matrix} R_{μ ν} & = & 2 N M_{μ ν} u (s (μ \to ν)) \end{matrix}

$\begin{matrix} R_{μ ν} & = & 2 N M_{μ ν} u (s (μ \to ν)) \end{matrix}$

where u(s(μ → ν)) is the fixation probability of mutants from μ to ν the selective advantage of which is equal to s.

For genic selection (no dominance) or gametic selection in a Wright-Fisher population of diploid, the fixation probability, u, of a single mutant gene, the selective advantage of which is equal to s and the frequency of which in a population is equal to $q_{m} = 1 / (2 N),$ $q_{m} = 1 / (2 N),$ was estimated (Crow and Kimura, 1970) as

(21)

\begin{matrix} 2 N u (s) & = & 2 N \frac{1 - e^{- 4 N_{e} s q_{m}}}{1 - e^{- 4 N_{e} s}} \end{matrix}

$\begin{matrix} 2 N u (s) & = & 2 N \frac{1 - e^{- 4 N_{e} s q_{m}}}{1 - e^{- 4 N_{e} s}} \end{matrix}$

(22)

\begin{matrix} = & \frac{u (s)}{u (0)} with q_{m} = \frac{1}{2 N} \end{matrix}

$\begin{matrix} = & \frac{u (s)}{u (0)} with q_{m} = \frac{1}{2 N} \end{matrix}$

where N_e is effective population size. Eq. (21) will be also valid for haploid population if 2N_e and 2N are replaced by N_e and N, respectively. Also, for Moran population of haploid, 4N_e and 2N should be replaced by N_e and N, respectively. Fixation probabilities for various selection models, which are compiled from p. 192 and p. 424–427 of Crow and Kimura (1970) and from Moran (1958) and Ewens (1979), are listed in Table S.7. The selective advantage of a mutant sequence ν to a wildtype μ is equal to

(23)

\begin{matrix} s (μ \to ν) & = & m (ν) - m (μ) \end{matrix}

$\begin{matrix} s (μ \to ν) & = & m (ν) - m (μ) \end{matrix}$

where m(ν) is the Malthusian fitness of a mutant sequence, and m(μ) is for the wildtype.

This Markov process of substitutions in sequence is reversible, and the equilibrium frequency of sequence μ, P^eq(μ), in the total process consisting of mutation and fixation processes is represented by

(24)

\begin{matrix} P^{e q} (μ) & = & \frac{P^{m u t} (μ) \exp (4 N_{e} m (μ) (1 - q_{m}))}{\sum_{ν} P^{m u t} (ν) \exp (4 N_{e} m (ν) (1 - q_{m}))} \end{matrix}

$\begin{matrix} P^{e q} (μ) & = & \frac{P^{m u t} (μ) \exp (4 N_{e} m (μ) (1 - q_{m}))}{\sum_{ν} P^{m u t} (ν) \exp (4 N_{e} m (ν) (1 - q_{m}))} \end{matrix}$

because both the mutation and fixation processes satisfy the detailed balance conditions, Eq. (19) and the following equation, respectively.

(25)

\begin{matrix} \exp (4 N_{e} m (μ) (1 - q_{m})) u (s (μ \to ν)) \\ = \frac{\exp (- 4 N_{e} m (μ) q_{m}) - \exp (- 4 N_{e} m (ν) q_{m})}{\exp (- 4 N_{e} m (μ)) - \exp (- 4 N_{e} m (ν))} \end{matrix}

$\begin{matrix} \exp (4 N_{e} m (μ) (1 - q_{m})) u (s (μ \to ν)) \\ = \frac{\exp (- 4 N_{e} m (μ) q_{m}) - \exp (- 4 N_{e} m (ν) q_{m})}{\exp (- 4 N_{e} m (μ)) - \exp (- 4 N_{e} m (ν))} \end{matrix}$

(26)

\begin{matrix} = \exp (4 N_{e} m (ν) (1 - q_{m})) u (s (ν \to μ)) \end{matrix}

$\begin{matrix} = \exp (4 N_{e} m (ν) (1 - q_{m})) u (s (ν \to μ)) \end{matrix}$

As a result, the ensemble of homologous sequences in molecular evolution obeys a Boltzmann distribution.

2.4. Relationships between m(σ), ψ_N(σ), and ΔG_ND(σ) of protein sequence

From Eqs. (1), (9), and (24), we can get the following relationships among the Malthusian fitness m, the folding free energy change ΔG_ND and Δψ_ND of protein sequence.

(27)

\begin{matrix} P^{e q} (μ) & = & \frac{P^{m u t} (μ) \exp (4 N_{e} m (μ) (1 - q_{m}))}{\sum_{ν} P^{m u t} (ν) \exp (4 N_{e} m (ν) (1 - q_{m})))} \end{matrix}

$\begin{matrix} P^{e q} (μ) & = & \frac{P^{m u t} (μ) \exp (4 N_{e} m (μ) (1 - q_{m}))}{\sum_{ν} P^{m u t} (ν) \exp (4 N_{e} m (ν) (1 - q_{m})))} \end{matrix}$

(28)

\begin{matrix} = & \frac{P^{m u t} (\bar{μ}) \exp (- (ψ_{N} (μ) - ψ_{D} (\bar{f (μ)}, T)))}{\sum_{ν} P^{m u t} (\bar{ν}) \exp (- (ψ_{N} (ν) - ψ_{D} (\bar{f (ν)}, T)))} \end{matrix}

$\begin{matrix} = & \frac{P^{m u t} (\bar{μ}) \exp (- (ψ_{N} (μ) - ψ_{D} (\bar{f (μ)}, T)))}{\sum_{ν} P^{m u t} (\bar{ν}) \exp (- (ψ_{N} (ν) - ψ_{D} (\bar{f (ν)}, T)))} \end{matrix}$

(29)

\begin{matrix} ≃ & \frac{P^{m u t} (μ) \exp (- Δ G_{N D} (μ, T) / (k_{B} T_{s}))}{\sum_{ν} P^{m u t} (ν) \exp (- Δ G_{N D} (ν, T) / (k_{B} T_{s}))} \end{matrix}

$\begin{matrix} ≃ & \frac{P^{m u t} (μ) \exp (- Δ G_{N D} (μ, T) / (k_{B} T_{s}))}{\sum_{ν} P^{m u t} (ν) \exp (- Δ G_{N D} (ν, T) / (k_{B} T_{s}))} \end{matrix}$

where

\bar{f (σ)} \equiv \sum_{σ} f (σ) P (σ)

$\bar{f (σ)} \equiv \sum_{σ} f (σ) P (σ)$ and

\log P^{m u t} (\bar{σ}) \equiv \sum_{σ} P (σ) \log (\prod_{i} P^{m u t} (σ_{i}))

$\log P^{m u t} (\bar{σ}) \equiv \sum_{σ} P (σ) \log (\prod_{i} P^{m u t} (σ_{i}))$ . Then, the following relationships are derived for sequences for which

f (μ) = \bar{f (μ)}

$f (μ) = \bar{f (μ)}$ .

(30)

\begin{matrix} 4 N_{e} m (μ) (1 - q_{m}) & = & - Δ ψ_{N D} (μ, T) + constant \end{matrix}

$\begin{matrix} 4 N_{e} m (μ) (1 - q_{m}) & = & - Δ ψ_{N D} (μ, T) + constant \end{matrix}$

(31)

\begin{matrix} ≃ & \frac{- Δ G_{N D} (μ, T)}{k_{B} T_{s}} + constant \end{matrix}

$\begin{matrix} ≃ & \frac{- Δ G_{N D} (μ, T)}{k_{B} T_{s}} + constant \end{matrix}$

The selective advantage of ν to μ is represented as follows for

f (μ) = f (ν) = \bar{f (σ)}

$f (μ) = f (ν) = \bar{f (σ)}$ .

(32)

\begin{matrix} 4 N_{e} s (μ \to ν) (1 - q_{m}) \\ = (4 N_{e} m (ν) - 4 N_{e} m (μ)) (1 - q_{m}) \end{matrix}

$\begin{matrix} 4 N_{e} s (μ \to ν) (1 - q_{m}) \\ = (4 N_{e} m (ν) - 4 N_{e} m (μ)) (1 - q_{m}) \end{matrix}$

(33)

\begin{matrix} = - (Δ ψ_{N D} (ν, T) - Δ ψ_{N D} (μ, T)) = - (ψ_{N} (ν) - ψ_{N} (μ)) \end{matrix}

$\begin{matrix} = - (Δ ψ_{N D} (ν, T) - Δ ψ_{N D} (μ, T)) = - (ψ_{N} (ν) - ψ_{N} (μ)) \end{matrix}$

(34)

\begin{matrix} ≃ - (Δ G_{N D} (ν, T) - Δ G_{N D} (μ, T)) / (k_{B} T_{s}) \\ = - (G_{N} (ν) - G_{N} (μ)) / (k_{B} T_{s}) \end{matrix}

$\begin{matrix} ≃ - (Δ G_{N D} (ν, T) - Δ G_{N D} (μ, T)) / (k_{B} T_{s}) \\ = - (G_{N} (ν) - G_{N} (μ)) / (k_{B} T_{s}) \end{matrix}$

It should be noted here that only sequences for which

f (σ) = \bar{f (σ)}

$f (σ) = \bar{f (σ)}$ contribute significantly to the partition functions in Eq. (28), and other sequences may be ignored.

Eq. (33) indicates that evolutionary statistical energy ψ should be proportional to effective population size N_e, and therefore it is ideal to estimate one-body (h) and two-body (J) interactions from homologous sequences of species that do not significantly differ in effective population size. Also, Eq. (34) indicates that selective temperature T_s is inversely proportional to the effective population size N_e; T_s∝1/N_e, because free energy is a physical quantity and should not depend on effective population size.

2.5. The ensemble average of folding free energy, ΔG_ND(σ, T), over sequences

The ensemble average of ΔG_ND(σ, T) over sequences with Eq. (1) is

(35)

\begin{matrix} {〈 Δ G_{N D} (σ, T) 〉}_{σ} \end{matrix}

$\begin{matrix} {〈 Δ G_{N D} (σ, T) 〉}_{σ} \end{matrix}$

(36)

\begin{matrix} \equiv [\sum_{σ} Δ G_{N D} (σ, T) P^{m u t} (σ) \exp (- \frac{Δ G_{N D} (σ, T)}{k_{B} T_{s}})] / \\ [\sum_{σ} P^{m u t} (σ) \exp (- \frac{Δ G_{N D} (σ, T)}{k_{B} T_{s}})] \end{matrix}

$\begin{matrix} \equiv [\sum_{σ} Δ G_{N D} (σ, T) P^{m u t} (σ) \exp (- \frac{Δ G_{N D} (σ, T)}{k_{B} T_{s}})] / \\ [\sum_{σ} P^{m u t} (σ) \exp (- \frac{Δ G_{N D} (σ, T)}{k_{B} T_{s}})] \end{matrix}$

(37)

\begin{matrix} \approx [\sum_{σ | f (σ) = \bar{f (σ_{N})}} G_{N} (σ) \exp (- \frac{G_{N} (σ)}{k_{B} T_{s}})] / \\ [\sum_{σ | f (σ) = \bar{f (σ_{N})}} \exp (- \frac{G_{N} (σ)}{k_{B} T_{s}})] - G_{D} (\bar{f (σ_{N})}, T) \end{matrix}

$\begin{matrix} \approx [\sum_{σ | f (σ) = \bar{f (σ_{N})}} G_{N} (σ) \exp (- \frac{G_{N} (σ)}{k_{B} T_{s}})] / \\ [\sum_{σ | f (σ) = \bar{f (σ_{N})}} \exp (- \frac{G_{N} (σ)}{k_{B} T_{s}})] - G_{D} (\bar{f (σ_{N})}, T) \end{matrix}$

(38)

\begin{matrix} = {〈 G_{N} (σ) 〉}_{σ} - G_{D} (\bar{f (σ_{N})}, T) \end{matrix}

$\begin{matrix} = {〈 G_{N} (σ) 〉}_{σ} - G_{D} (\bar{f (σ_{N})}, T) \end{matrix}$

where σ_N denotes a natural sequence, and

\bar{f (σ_{N})}

$\bar{f (σ_{N})}$ denotes the average of amino acid frequencies f(σ_N) over homologous sequences. In Eq. (37), the sum over all sequences is approximated by the sum over sequences the amino acid composition of which is the same as that over the natural sequences.

The ensemble averages of G_N and ψ_N(σ) are estimated in the Gaussian approximation (Pande et al., 1997).

(39)

\begin{matrix} {〈 G_{N} (σ) 〉}_{σ} & \approx & \frac{\int E \exp (- E / (k_{B} T_{s})) n (E) d E}{\int \exp (- E / (k_{B} T_{s})) n (E) d E} \end{matrix}

$\begin{matrix} {〈 G_{N} (σ) 〉}_{σ} & \approx & \frac{\int E \exp (- E / (k_{B} T_{s})) n (E) d E}{\int \exp (- E / (k_{B} T_{s})) n (E) d E} \end{matrix}$

(40)

\begin{matrix} = & \bar{E} (\bar{f (σ_{N})}) - {δ E}^{2} (\bar{f (σ_{N})}) / (k_{B} T_{s}) \end{matrix}

$\begin{matrix} = & \bar{E} (\bar{f (σ_{N})}) - {δ E}^{2} (\bar{f (σ_{N})}) / (k_{B} T_{s}) \end{matrix}$

(41)

\begin{matrix} {〈 ψ_{N} (σ) 〉}_{σ} & \equiv & [\sum_{σ} ψ_{N D} (σ) \exp (- ψ_{N} (σ))] / \\ [\sum_{σ} \exp (- ψ_{N} (σ))] \end{matrix}

$\begin{matrix} {〈 ψ_{N} (σ) 〉}_{σ} & \equiv & [\sum_{σ} ψ_{N D} (σ) \exp (- ψ_{N} (σ))] / \\ [\sum_{σ} \exp (- ψ_{N} (σ))] \end{matrix}$

(42)

\begin{matrix} \approx & \bar{ψ} (\bar{f (σ_{N})}) - {δ ψ}^{2} (\bar{f (σ_{N})}) \end{matrix}

$\begin{matrix} \approx & \bar{ψ} (\bar{f (σ_{N})}) - {δ ψ}^{2} (\bar{f (σ_{N})}) \end{matrix}$

The ensemble averages of ΔG_ND(σ, T) and ψ_N(σ) over sequences are observable as the sample averages of ΔG_ND(σ_N, T) and ψ_N(σ_N) over homologous sequences fixed in protein evolution, respectively.

(43)

\begin{matrix} \bar{Δ G_{N D} (σ_{N}, T)} / (k_{B} T_{s}) & = & {〈 Δ G_{N D} (σ, T) 〉}_{σ} / (k_{B} T_{s}) \end{matrix}

$\begin{matrix} \bar{Δ G_{N D} (σ_{N}, T)} / (k_{B} T_{s}) & = & {〈 Δ G_{N D} (σ, T) 〉}_{σ} / (k_{B} T_{s}) \end{matrix}$

(44)

\begin{matrix} \approx & {δ ψ}^{2} (\bar{f (σ_{N})}) [ϑ (T / T_{g}) T_{s} / T - 1] \end{matrix}

$\begin{matrix} \approx & {δ ψ}^{2} (\bar{f (σ_{N})}) [ϑ (T / T_{g}) T_{s} / T - 1] \end{matrix}$

(45)

\begin{matrix} \bar{ψ_{N} (σ_{N})} & \equiv & \frac{\sum_{σ_{N}} w_{σ_{N}} ψ_{N} (σ_{N})}{\sum_{σ_{N}} w_{σ_{N}}} \end{matrix}

$\begin{matrix} \bar{ψ_{N} (σ_{N})} & \equiv & \frac{\sum_{σ_{N}} w_{σ_{N}} ψ_{N} (σ_{N})}{\sum_{σ_{N}} w_{σ_{N}}} \end{matrix}$

(46)

\begin{matrix} = & {〈 ψ_{N} (σ) 〉}_{σ} \end{matrix}

$\begin{matrix} = & {〈 ψ_{N} (σ) 〉}_{σ} \end{matrix}$

where the overline denotes a sample average with a sample weight

w_{σ_{N}}

$w_{σ_{N}}$ for each homologous sequence, which is used to reduce phylogenetic biases in the set of homologous sequences.

The folding free energy becomes equal to zero at the melting temperature T_m; ${〈 Δ G_{N D} (σ_{N}, T_{m}) 〉}_{σ} = 0$ ${〈 Δ G_{N D} (σ_{N}, T_{m}) 〉}_{σ} = 0$ . Thus, the following relationship must be satisfied (Pande et al., 1997; Ramanathan and Shakhnovich, 1994; Shakhnovich and Gutin, 1993a; 1993b).

(47)

\begin{matrix} ϑ (T_{m} / T_{g}) \frac{T_{s}}{T_{m}} & = & \frac{T_{s}}{2 T_{m}} (1 + \frac{T_{m}^{2}}{T_{g}^{2}}) = 1 with T_{s} \leq T_{g} \leq T_{m} \end{matrix}

$\begin{matrix} ϑ (T_{m} / T_{g}) \frac{T_{s}}{T_{m}} & = & \frac{T_{s}}{2 T_{m}} (1 + \frac{T_{m}^{2}}{T_{g}^{2}}) = 1 with T_{s} \leq T_{g} \leq T_{m} \end{matrix}$

2.6. Probability distributions of selective advantage, fixation rate and K_a/K_s

Let us consider the probability distributions of characteristic quantities that describe the evolution of genes. First of all, the probability density function (PDF) of selective advantage s, p(s), of mutant genes can be calculated from the PDF of the change of Δψ_ND due to a mutation from μ to ν, $Δ Δ ψ_{N D} (\equiv Δ ψ_{N D} (ν, T) - Δ ψ_{N D} (μ, T))$ $Δ Δ ψ_{N D} (\equiv Δ ψ_{N D} (ν, T) - Δ ψ_{N D} (μ, T))$ . The PDF of 4N_es, $p (4 N_{e} s) = p (s) / (4 N e),$ $p (4 N_{e} s) = p (s) / (4 N e),$ may be more useful than p(s).

(48)

\begin{matrix} p (4 N_{e} s) & = & p (Δ Δ ψ_{N D}) | \frac{d Δ Δ ψ_{N D}}{d 4 N_{e} s} | = p (Δ Δ ψ_{N D}) (1 - q_{m}) \end{matrix}

$\begin{matrix} p (4 N_{e} s) & = & p (Δ Δ ψ_{N D}) | \frac{d Δ Δ ψ_{N D}}{d 4 N_{e} s} | = p (Δ Δ ψ_{N D}) (1 - q_{m}) \end{matrix}$

where ΔΔψ_ND must be regarded as a function of 4N_es, that is,

Δ Δ ψ_{N D} = - 4 N_{e} s (1 - q_{m})

$Δ Δ ψ_{N D} = - 4 N_{e} s (1 - q_{m})$ ; see Eq. (33).

The PDF of fixation probability u can be represented by

(49)

\begin{matrix} p (u) = p (4 N_{e} s) \frac{d 4 N_{e} s}{d u} = p (4 N_{e} s) \frac{{(e^{4 N_{e} s} - 1)}^{2} e^{4 N_{e} s (q_{m} - 1)}}{q_{m} (e^{4 N_{e} s} - 1) - (e^{4 N_{e} s q_{m}} - 1)} \end{matrix}

$\begin{matrix} p (u) = p (4 N_{e} s) \frac{d 4 N_{e} s}{d u} = p (4 N_{e} s) \frac{{(e^{4 N_{e} s} - 1)}^{2} e^{4 N_{e} s (q_{m} - 1)}}{q_{m} (e^{4 N_{e} s} - 1) - (e^{4 N_{e} s q_{m}} - 1)} \end{matrix}$

where 4N_es must be regarded as a function of u.

The ratio of the substitution rate per nonsynonymous site (K_a) for nonsynonymous substitutions with selective advantage s to the substitution rate per synonymous site (K_s) for synonymous substitutions with s = 0 is

(50)

\begin{matrix} \frac{K_{a}}{K_{s}} & = & \frac{u (s)}{u (0)} = \frac{u (s)}{q_{m}} \end{matrix}

$\begin{matrix} \frac{K_{a}}{K_{s}} & = & \frac{u (s)}{u (0)} = \frac{u (s)}{q_{m}} \end{matrix}$

assuming that synonymous substitutions are completely neutral and mutation rates at both types of sites are the same. The PDF of K_a/K_s is

(51)

\begin{matrix} p (K_{a} / K_{s}) & = & p (u) \frac{d u}{d (K_{a} / K_{s})} = p (u) q_{m} \end{matrix}

$\begin{matrix} p (K_{a} / K_{s}) & = & p (u) \frac{d u}{d (K_{a} / K_{s})} = p (u) q_{m} \end{matrix}$

2.7. Probability distributions of ΔΔψ_ND, 4N_es, u, and K_a/K_s in fixed mutant genes

The PDF of ΔΔψ_ND in fixed mutants is proportional to that multiplied by the fixation probability.

(52)

\begin{matrix} p (Δ Δ ψ_{N D, f i x e d}) & = & p (Δ Δ ψ_{N D}) \frac{u (s (Δ Δ ψ_{N D}))}{〈 u (s (Δ Δ ψ_{N D})) 〉} \end{matrix}

$\begin{matrix} p (Δ Δ ψ_{N D, f i x e d}) & = & p (Δ Δ ψ_{N D}) \frac{u (s (Δ Δ ψ_{N D}))}{〈 u (s (Δ Δ ψ_{N D})) 〉} \end{matrix}$

(53)

\begin{matrix} 〈 u 〉 & \equiv & \int_{- \infty}^{\infty} u (s) p (Δ Δ ψ_{N D}) d Δ Δ ψ_{N D} \end{matrix}

$\begin{matrix} 〈 u 〉 & \equiv & \int_{- \infty}^{\infty} u (s) p (Δ Δ ψ_{N D}) d Δ Δ ψ_{N D} \end{matrix}$

Likewise, the PDF of selective advantage in fixed mutants is

(54)

\begin{matrix} p (4 N_{e} s_{f i x e d}) & = & p (4 N_{e} s) \frac{u (s)}{〈 u (s) 〉} \end{matrix}

$\begin{matrix} p (4 N_{e} s_{f i x e d}) & = & p (4 N_{e} s) \frac{u (s)}{〈 u (s) 〉} \end{matrix}$

and those of the u and K_a/K_s in fixed mutants are

(55)

\begin{matrix} p (u_{f i x e d}) & = & p (u) \frac{u}{〈 u 〉} \end{matrix}

$\begin{matrix} p (u_{f i x e d}) & = & p (u) \frac{u}{〈 u 〉} \end{matrix}$

(56)

\begin{matrix} p ({(\frac{K_{a}}{K_{s}})}_{f i x e d}) & = & p (\frac{K_{a}}{K_{s}}) \frac{u}{〈 u 〉} = p (\frac{K_{a}}{K_{s}}) \frac{\frac{K_{a}}{K_{s}}}{〈 \frac{K_{a}}{K_{s}} 〉} \end{matrix}

$\begin{matrix} p ({(\frac{K_{a}}{K_{s}})}_{f i x e d}) & = & p (\frac{K_{a}}{K_{s}}) \frac{u}{〈 u 〉} = p (\frac{K_{a}}{K_{s}}) \frac{\frac{K_{a}}{K_{s}}}{〈 \frac{K_{a}}{K_{s}} 〉} \end{matrix}$

The average of K_a/K_s in fixed mutants is equal to the ratio of the second moment to the first moment of K_a/K_s in all arising mutants;

{〈 K_{a} / K_{s} 〉}_{f i x e d} = 〈 {(K_{a} / K_{s})}^{2} 〉 / 〈 K_{a} / K_{s} 〉

${〈 K_{a} / K_{s} 〉}_{f i x e d} = 〈 {(K_{a} / K_{s})}^{2} 〉 / 〈 K_{a} / K_{s} 〉$ .

3. Materials

3.1. Sequence data

We study the single domains of 8 Pfam (Finn et al., 2016) families and both the single domains and multi-domains from 3 Pfam families. In Table 1, their Pfam ID for a multiple sequence alignment, and UniProt ID and PDB ID with the starting- and ending-residue positions of the domains are listed. The full alignments for their families at the Pfam are used to estimate one-body interactions h and pairwise interactions J with the DCA program from “http://dca.rice.edu/portal/dca/home ” (Marks et al., 2011; Morcos et al., 2011). To estimate the sample ( $\bar{ψ_{N}}$ $\bar{ψ_{N}}$ ) and ensemble (⟨ψ_N⟩_σ) averages of the evolutionary statistical energy, M unique sequences with no deletions are used. In order to reduce phylogenetic biases in the set of homologous sequences, we employ a sample weight ( $w_{σ_{N}}$ $w_{σ_{N}}$ ) for each sequence, which is equal to the inverse of the number of sequences that are less than 20% different from a given sequence in a given set of homologous sequences. Only representatives of unique sequences with no deletions, which are at least 20% different from each other, are used to calculate the changes of the evolutionary statistical energy (Δψ_N) due to single nucleotide nonsynonymous substitutions; the number of the representatives is almost equal to the effective number of sequences (M_eff) in Table 1.

Table 1. Protein families, and structures studied.

Pfam family	UniProt ID	N^a	N_eff^b,c	M^d	M_eff^c,e	L^f	PDB ID
HTH_3	RPC1_BP434/7-59	15315(15917)	11691.21	6286	4893.73	53	1R69-A:6-58
Nitroreductase	Q97IT9_CLOAB/4-76	6008(6084)	4912.96	1057	854.71	73	3E10-A/B:4-76^g
SBP_bac_3^h	GLNH_ECOLI/27-244	9874(9972)	7374.96	140	99.70	218	1WDN-A:5-222
SBP_bac_3	GLNH_ECOLI/111-204	9712(9898)	7442.85	829	689.64	94	1WDN-A:89-182
OmpA	PAL_ECOLI/73-167	6035(6070)	4920.44	2207	1761.24	95	1OAP-A:52-146
DnaB	DNAB_ECOLI/31-128	1929(1957)	1284.94	1187	697.30	98	1JWE-A:30-127
LysR_substrate^h	BENM_ACIAD/90-280	25138(25226)	20707.06	85(1)	67.00	191	2F6G-A/B:90-280^g
LysR_substrate	BENM_ACIAD/163-265	25032(25164)	21144.74	121(1)	99.27	103	2F6G-A/B:163-265^g
Methyltransf_5^h	RSMH_THEMA/8-292	1942(1953)	1286.67	578(2)	357.97	285	1N2X-A:8-292
Methyltransf_5	RSMH_THEMA/137-216	1877(1911)	1033.35	975(2)	465.53	80	1N2X-A:137-216
SH3_1	SRC_HUMAN:90-137	9716(16621)	3842.47	1191	458.31	48	1FMK-A:87-134
ACBP	ACBP_BOVIN/3-82	2130(2526)	1039.06	161	70.72	80	2ABD-A:2-81
PDZ	PTN13_MOUSE/1358-1438	13814(23726)	4748.76	1255	339.99	81	1GM1-A:16-96
Copper-bind	AZUR_PSEAE:24-148	1136(1169)	841.56	67(1)	45.23	125	5AZU-B/C:4-128^g

a: The number of unique sequences and the total number of sequences in parentheses; the full alignments in the Pfam (Finn et al., 2016) are used.
b: The effective number of sequences.
c: A sample weight ( $w_{σ_{N}}$ $w_{σ_{N}}$ ) for a given sequence is equal to the inverse of the number of sequences that are less than 20% different from the given sequence.
d: The number of unique sequences that include no deletion unless specified. The number in parentheses indicates the maximum number of deletions allowed.
e: The effective number of unique sequences that include no deletion or at most the specified number of deletions.
f: The number of residues.
g: Contacts are calculated in the homodimeric state for these protein.
h: These proteins consist of two domains, and other ones are single domains.

4. Results

First, We describe how one-body and pairwise interactions, h and J, are estimated. Then, the changes of evolutionary statistical energy (Δψ_N) due to single nucleotide nonsynonymous changes on natural sequences are analyzed with respect to dependences on the ψ_N of the wildtype sequences. The results indicate that the standard deviation of ΔG_N ≃ k_BT_sΔψ_N is almost constant over protein families. Hence, the selective temperatures, T_s, of various protein families can be estimated in a relative scale from the standard deviation of Δψ_N. The T_s of a reference protein is estimated by comparing the expected values of ΔΔG_ND with their experimental values. Folding free energies ΔG_ND are estimated from estimated T_s and experimental melting temperature T_m, and compared with their experimental values for 5 protein families. Glass transition temperatures T_g are also estimated from T_s and T_m.

Secondly, based on the distribution of Δψ_N, protein evolution is studied. Evolutionary statistical energy (ψ_N) attains the equilibrium when the average of Δψ_N over fixed mutations is equal to zero. The PDF of Δψ_N is approximated by log-normal distributions. The basic relationships are that 1) the standard deviation of Δψ_N is constant specific to a protein family, and 2) the mean of Δψ_N linearly depends on ψ_N. The equilibrium value of ψ_N is shown to agree with the mean of ψ_N over homologous proteins in each protein family. In the present approximation, the standard deviation of Δψ_N and selective temperature T_s at the equilibrium are simple functions of the equilibrium value of mean Δψ_N, ${\bar{Δ ψ_{N}}}^{e q}$ ${\bar{Δ ψ_{N}}}^{e q}$ . Lastly, the probability distribution of K_a/K_s, which is the ratio of nonsynonymous to synonymous substitution rate per site, is analyzed as a function of ${\bar{Δ ψ_{N}}}^{e q},$ ${\bar{Δ ψ_{N}}}^{e q},$ in order to examine how significant neutral selection is in the selection maintaining protein stability and foldability. Also, it is confirmed that selective temperature T_s negatively correlates with the mean of K_a/K_s, which represents the evolutionary rate of protein.

4.1. Important parameters in the estimations of one-body and pairwise interactions, h and J, and of the evolutionary statistical energy, ψ_N(σ)

The one-body (h) and pairwise (J) interactions for amino acid order in a protein sequence are estimated here by the DCA method (Marks et al., 2011; Morcos et al., 2011), although there are multiple methods for estimating them (Ekeberg et al., 2014; 2013). In the case of the DCA method, the ratio of pseudocount (0 ≤ p_c ≤ 1) defined in Eqs. (S.70) and (S.71) is a parameter and controls the values of the ensemble and sample averages of ψ_N in sequence space, ⟨ψ_N(σ)⟩_σ in Eq. (42) and $\bar{ψ_{N} (σ_{N})}$ $\bar{ψ_{N} (σ_{N})}$ in Eq. (45); a weight for observed counts is defined to be equal to $(1 - p_{c})$ $(1 - p_{c})$ . Sample average means the average over all homologous sequences with a weight for each sequence to reduce phylogenetic biases. An appropriate value must be chosen for the ratio of pseudocount in a reasonable manner.

Another problem is that the estimates of h and J(Marks et al., 2011; Morcos et al., 2011) may be noisy as a result of estimating many interaction parameters from a relatively small number of sequences. Therefore, only pairwise interactions within a certain distance are taken into account; the estimate of J is modified as follows, according to Morcos et al. (2014).

(57)

\begin{matrix} {\hat{J}}_{i j}^{q} (a_{k}, a_{l}) & = & J_{i j}^{q} (a_{k}, a_{l}) H (r_{c u t o f f} - r_{i j}) \end{matrix}

$\begin{matrix} {\hat{J}}_{i j}^{q} (a_{k}, a_{l}) & = & J_{i j}^{q} (a_{k}, a_{l}) H (r_{c u t o f f} - r_{i j}) \end{matrix}$

where J^q is the statistical estimate of J in the mean field approximation in which the amino acid a_q is the reference state, H is the Heaviside step function, and r_ij is the distance between the centers of amino acid side chains at sites i and j in a protein structure, and r_cutoff is a distance threshold for residue pairwise interactions. The one-body interactions h_i(a_k) are estimated in the isolated two-state model (Morcos et al., 2011) rather than the mean field approximation; see the Method section in the Supplement for details. The zero-sum gauge is employed to represent h and J;

\sum_{k} {\hat{h}}_{i}^{s} (a_{k}) = \sum_{k} \sum_{l} {\hat{J}}_{i j}^{s} (a_{k}, a_{l}) = 0

$\sum_{k} {\hat{h}}_{i}^{s} (a_{k}) = \sum_{k} \sum_{l} {\hat{J}}_{i j}^{s} (a_{k}, a_{l}) = 0$ in the zero-sum gauge.

Candidates for the cutoff distance may be about 8 Å for the first interaction shell and 15–16 Å for the second interaction shell between residues; distance between the centers of side chain atoms is employed for residue distance. Here both the distances are tested for the cutoff distance. Pseudocount in the Bayesian statistics is determined usually as a function of the number of samples (sequences), although the ratio of pseudocount $p_{c} = 0.5$ $p_{c} = 0.5$ was used for all proteins in the contact prediction (Morcos et al., 2011). Here, an appropriate value for the ratio of pseudocount for the certain cutoff distance, either about 8 Å or 15–16 Å, is chosen for each protein family in such a way that the sample average of the evolutionary statistical energies must be equal to the ensemble average, $\bar{ψ_{N}} = {〈 ψ_{N} 〉}_{σ}$ $\bar{ψ_{N}} = {〈 ψ_{N} 〉}_{σ}$ ; see Eqs. (42) and (46). As shown in Fig. S.1, the value of r_cutoff, where $\bar{ψ_{N}} = {〈 ψ_{N} 〉}_{σ}$ $\bar{ψ_{N}} = {〈 ψ_{N} 〉}_{σ}$ is satisfied, monotonously changes as a function of the ratio of pseudocount p_c. The values of p_c, where $\bar{ψ_{N}} = {〈 ψ_{N} 〉}_{σ}$ $\bar{ψ_{N}} = {〈 ψ_{N} 〉}_{σ}$ is satisfied near the specified values of r_cutoff, 8 Å and 15.5 Å, are employed for r_cutoff ≃ 8 Å and 15.5 Å, respectively. In the present multiple sequence alignment for the PDZ domain, with the ratios of pseudocount $p_{c} = 0.205$ $p_{c} = 0.205$ and $p_{c} = 0.33,$ $p_{c} = 0.33,$ the sample and ensemble averages agree with each other at the cutoff distances r_cutoff ∼ 8 Å and r_cutoff ∼ 15.5 Å, respectively; see Fig. S.1. In Fig. S.2, the reflective correlation and regression coefficients between the experimental ΔΔG_ND (Gianni et al., 2007) and Δψ_N due to single amino acid substitutions are plotted against the cutoff distance for pairwise interactions in the PDZ domain. The reflective correlation coefficient has the maximum at the r_cutoff ∼ 8 Å for $p_{c} = 0.205$ $p_{c} = 0.205$ and at r_cutoff ∼ 15.5 Å for

Details & settings

Outline

Figures (15)

Tables (3)

Journal of Theoretical Biology

Selection originating from protein stability/foldability: Relationships between protein folding free energy, sequence ensemble, and fitness

Highlights

Abstract

Keywords

1. Introduction

2. Methods

2.1. Knowledge of protein folding

2.2. Probability distribution of homologous sequences with the same native fold in sequence space

2.3. The equilibrium distribution of sequences in a mutation-fixation process

2.4. Relationships between m(σ), ψ_N(σ), and ΔG_ND(σ) of protein sequence

2.5. The ensemble average of folding free energy, ΔG_ND(σ, T), over sequences

2.6. Probability distributions of selective advantage, fixation rate and K_a/K_s

2.7. Probability distributions of ΔΔψ_ND, 4N_es, u, and K_a/K_s in fixed mutant genes

3. Materials

3.1. Sequence data

4. Results

4.1. Important parameters in the estimations of one-body and pairwise interactions, h and J, and of the evolutionary statistical energy, ψ_N(σ)

Details & settings

Outline

Figures (15)

Tables (3)

Journal of Theoretical Biology

Selection originating from protein stability/foldability: Relationships between protein folding free energy, sequence ensemble, and fitness

Highlights

Abstract

Keywords

1. Introduction

2. Methods

2.1. Knowledge of protein folding

2.2. Probability distribution of homologous sequences with the same native fold in sequence space

2.3. The equilibrium distribution of sequences in a mutation-fixation process

2.4. Relationships between m(σ), ψN(σ), and ΔGND(σ) of protein sequence

2.5. The ensemble average of folding free energy, ΔGND(σ, T), over sequences

2.6. Probability distributions of selective advantage, fixation rate and Ka/Ks

2.7. Probability distributions of ΔΔψND, 4Nes, u, and Ka/Ks in fixed mutant genes

3. Materials

3.1. Sequence data

4. Results

4.1. Important parameters in the estimations of one-body and pairwise interactions, h and J, and of the evolutionary statistical energy, ψN(σ)

2.4. Relationships between m(σ), ψ_N(σ), and ΔG_ND(σ) of protein sequence

2.5. The ensemble average of folding free energy, ΔG_ND(σ, T), over sequences

2.6. Probability distributions of selective advantage, fixation rate and K_a/K_s

2.7. Probability distributions of ΔΔψ_ND, 4N_es, u, and K_a/K_s in fixed mutant genes

4.1. Important parameters in the estimations of one-body and pairwise interactions, h and J, and of the evolutionary statistical energy, ψ_N(σ)