Influence of mutation bias and hydrophobicity on the substitution rates and sequence entropies of protein evolution

Stability constrained fitness model

Stability constrained models of protein evolution assume that the fitness of a protein with sequence A is proportional to the fraction of protein that is in the native state, which can be computed from the folding free energy as (Goldstein, 2011; Serohijos & Shakhnovich, 2014) (1)${\displaystyle f\left(\mathbf{A}\right)={\mathrm{e}}^{-\Delta G\left(\mathbf{A}\right)∕kT}∕\left(1+{\mathrm{e}}^{-\Delta G\left(\mathbf{A}\right)∕kT}\right).}$ The computation is performed assuming that the native contact matrix C^nat does not change in evolution. Upon single mutation, the free energy change ΔG_mut = ΔG_wt + ΔΔG is predicted adopting some models of protein stability (see below).

Equilibrium distribution

Another approximation that is often used in these models is that the mutation rate is extremely slow (Nμ ≪ 1) so that at every time there is only one mutant gene that “competes” with the wild-type gene for fixation in the population with effective population size of N individuals. Under this scenario, the probability that the mutation gets fixed in the population can be computed with Kimura’s formula (Kimura, 1962) as (2)${\displaystyle P_{\mathrm{fix}}\left({\mathbf{A}}^{\mathrm{wt}}\to {\mathbf{A}}^{\mathrm{mut}}\right)=\frac{{\mathrm{e}}^{-\left(\phi \left({\mathbf{A}}^{\mathrm{mut}}\right)-\phi \left({\mathbf{A}}^{\mathrm{wt}}\right)\right)}-1}{{\mathrm{e}}^{-N\left(\phi \left({\mathbf{A}}^{\mathrm{mut}}\right)-\phi \left({\mathbf{A}}^{\mathrm{wt}}\right)\right)}-1}}$ where φ(A) = log(f(A)) is the logarithmic fitness associated with the amino acid sequence A, Eq. (1). As it is well known, the fixation probability tends to the neutral limit P_fix = 1∕N when Δφ tends to zero, it tends exponentially to zero when Δφ is negative and large, and it tends to 1 − e^−Δφ when Δφ is positive. Nearly neutral mutations with selective effect |Δφ| ≈ 1∕N are likely to be fixed even when their effect is deleterious (Ohta, 1976). Importantly, the above fixation probability defines a Monte Carlo process in sequence space that fulfils detailed balance, so that its stationary distribution can be computed exactly (Sella & Hirsh, 2005, Mustonen & Lässig, 2005), except for the normalization constant, which would require a sum over 20^L possible sequences A = A₁⋯A_L: (3)${\displaystyle P\left(A_1\cdots A_L\right)\propto exp\left(\left(N-1\right)\phi \left(A_1\cdots A_L\right)\right)}$ Note the analogy between this formula and the Boltzmann distribution with energy equal to −φ and temperature equal to 1∕(N − 1). This explicit formula holds when the mutation process is unbiased, so that all sequences are equally probable under the mutation model. In the presence of mutation bias, the stationary distribution can be determined as the distribution with minimal Kullback–Leibler divergence from the mutational distribution, $d_{\mathrm{KL}}={\sum }_A{P}^{\mathrm{mut}}\left(\mathbf{A}\right)\left[log\left(P^{\mathrm{mut}}\left(\mathbf{A}\right)\right)-log\left(P\left(\mathbf{A}\right)\right)\right]$, with a constraint on the average fitness ∑_AP(A)φ(A). This condition generalizes the Boltzmann principle, and it was adopted for developing the mean-field model of protein evolution (Arenas, Sanchez-Cobos & Bastolla, 2015).

Mean-field model of protein evolution

The mean-field (MF) model assumes that the equilibrium amino acid distribution is the product of independent distributions at each protein site, (4)${\displaystyle P\left(A_1,\dots ,A_L\right)=\prod _{i=1}^L{P}^i\left(A_i\right).}$ Of course this assumption is not realistic, since different sites determine protein stability through their interactions, but it is needed for performing likelihood computations in an efficient way. Our strategy consists in determining the effect of a mutation at site i self-consistently, with respect to the MF distribution at all other sites. For simplicity, we shall sometimes use the vectorial notation $P_a^i$ for indicating Pⁱ(a), where a denotes one of the twenty amino acid types.

The mean-field distribution is determined by minimizing the Kullback–Leibler divergence (distance between distribution) with respect to a global mutational distribution $P_a^{\mathrm{mut}}$, i.e.,  ${\sum }_{ia}{P}_a^{i}log\left(P_a^i∕P_a^{\mathrm{mut}}\right)$. We impose a constraint on the average fitness, which is transformed into a constraint on the folding free energy ΔG. This condition on stability is imposed through the Lagrange multiplier Λ that represents the strength of selection and is related with the effective population size. Furthermore, we impose the normalization constraints ${\sum }_a{P}_a^i=1$ at all sites.

Since the parameters that determine the folding free energy are fixed for all proteins (see below), the only free parameters of the model are Λ and $P_a^{\mathrm{mut}}$. The frequencies are generally determined from the observed sequences in the protein of known structures and the other sequences of the protein family, while Λ is determined by maximizing the log-likelihood of the PDB sequence, ${\sum }_{i}log\left(P^i\left(A_i^{\mathrm{PDB}}\right)\right)$, which yields a well-defined single maximum. The pre-computation of the moments of the contacts makes the computation very fast, it runs in a few minutes even for proteins of several hundreds of amino acids. For further computational details see (Arenas, Sanchez-Cobos & Bastolla, 2015).

Wild-type model of protein evolution

In the wild-type model (Jimenez, Arenas & Bastolla, 2018), we also assume that sites evolve independently. We further assume that the site-specific distribution $P_a^i$ of amino acid a at position i is proportional to the background distribution $P_a^{\mathrm{mut}}$ multiplied by the exponential of the logarithmic fitness of the corresponding mutation in which the wild-type amino acid in the PDB $A_i^{WT}$ is substituted by the new amino acid a: (5)${\displaystyle P_a^{\mathrm{WT},i}\propto P_a^{\mathrm{mut}}exp\left(\Lambda \phi \left(\mathrm{mut}\left(A_i^{WT}\to a\right)\right)\right).}$ The fitness of a sequence is computed as in Eq. (1). The parameter Λ is again determined by maximizing the likelihood of the wild-type sequence, ${\sum }_{i}log\left(P^{\mathrm{WT},i}\left(A_i^{\mathrm{WT}}\right)\right)$.

Sequence entropy

The sequence entropy at position i measures the variability of this position as (6)${\displaystyle S_i=-\sum _{a=1}^{20}{P}_a^{i}log\left(P_a^i\right),}$ where $P_a^i$ is obtained either from the evolutionary model (mean-field or wild-type) or from a MSA or from pooled amino acids at equivalent structural positions with the same number of contacts.

Halpern-Bruno exchangeability matrices

To fully specify the site-specific substitution processes, besides the site-specific frequencies $P_a^i$ we need to compute consistent exchangeability matrices with the Halpern-Bruno formulas (Halpern & Bruno, 1998).

Given a site-specific amino acid distribution that reflects selective constraints, the Halpern-Bruno method allows computing the rate matrices of the associated site-specific substitution processes $Q_{ab}^i=E_{ab}^i{P}_b^i$ that are produced by a global (not site-specific) mutation process consistent together with Kimura’s fixation probability, Eq. (2).

Without loss of generality, we parametrize the rate matrix of the global mutation process as $Q_{ab}^{\mathrm{mut}}=E_{ab}^{\mathrm{mut}}{P}_b^{\mathrm{mut}}$, where $P_a^{\mathrm{mut}}$ is the stationary matrix of the mutation process and $E_{ab}^{\mathrm{mut}}$ is its exchangeability matrix. To simplify formulas, here we assume detailed balance, i.e., we assume that $E_{ab}^{\mathrm{mut}}$ is a symmetric matrix (this condition can be easily relaxed). We write the rate matrices as $Q_{ab}^i=Q_{ab}^{\mathrm{mut}}{P}_{\mathrm{fix}}\left(f_a^i,f_b^i\right)$, where $f_a^i$ is the “fitness” of amino acid a at site i. We impose that P_fix is the fixation probability Eq. (2). Halpern and Bruno showed that the site-specific fitness can be inferred from the stationary distribution from $P_a^i=P_a^{\mathrm{mut}}{\left(f_a^i\right)}^N$, yielding the following site-specific substitution process (7)${\displaystyle Q_{ab}^i=E_{ab}^i{P}_b^i}$ (8)${\displaystyle E_{ab}^i=E_{ab}^{\mathrm{mut}}\left(\frac{ln\left(F_b^{\mathrm{sel},i}\right)-ln\left(F_a^{\mathrm{sel},i}\right)}{F_b^{\mathrm{sel},i}-F_a^{\mathrm{sel},i}}\right)}$ (9)${\displaystyle \mathrm{with}{F}_a^{\mathrm{sel},i}=\frac{P_a^i}{P_a^{\mathrm{mut}}}}$ The selective factors $F_a^{\mathrm{sel},i}$ quantify how much the site-specific distribution $P_a^i$ deviates from the background distribution $P_a^{\mathrm{mut}}$ induced by mutation alone.

It can be immediately seen that the exchangeability matrices $E_{ab}^i$ are symmetric, which implies that detailed balance holds and $P_a^i$ is the stationary distribution.

Evolutionary rates

For neutral substitutions with $F_a^{\mathrm{sel},i}=F_b^{\mathrm{sel},i}$, in particular synonymous substitutions a = b, applying l’Hopital’s rule we find $E_{ab}^i=E_{ab}^{\mathrm{mut}}∕F_b^{\mathrm{sel},i}$ and $Q_{ab}^i=Q_{ab}^{\mathrm{mut}}$, i.e., the rate of synonymous substitutions equals the mutation rate, in agreement with Kimura’s theory. If the amino acid b is favored by selection with respect to amino acid a, $F_b^{\mathrm{sel},i}>F_a^{\mathrm{sel},i}$, then the substitution rate is enhanced with respect to the neutral rate, and it is decreased in the opposite case. Because of detailed balance, the flux in one direction and the other are equal, $R_{ab}^i=P_a^i{P}_b^i{E}_{ab}^i=R_{ba}^i$, with (10)${\displaystyle R_{ab}^i=\left(P_a^{\mathrm{mut}}{P}_b^{\mathrm{mut}}{E}_{ab}^{\mathrm{mut}}\right)F_a^{\mathrm{sel},i}{F}_b^{\mathrm{sel},i}\frac{ln\left(F_b^{\mathrm{sel},i}\right)-ln\left(F_a^{\mathrm{sel},i}\right)}{F_b^{\mathrm{sel},i}-F_a^{\mathrm{sel},i}}}$ In the above equation, the flux is partitioned into a global component that is attributed to the mutation process (superscript mut) and a site-specific component that is attributed to selection (superscript sel), which allows analysing the contributions of mutation and selection separately. The flux is maximal for substitutions ab that have large and almost equal selective factors $F_a^{\mathrm{sel},i}\approx F_b^{\mathrm{sel},i}$ and have large mutational flux $P_a^{\mathrm{mut}}{P}_b^{\mathrm{mut}}{E}_{ab}^{\mathrm{mut}}$. The site-specific substitution rates are computed as the weighted average of the substitution rate matrix $Q_{ab}=E_{ab}^i{P}_b^i$, (11)${\displaystyle R^i=\sum _{a\ne b}{P}_a^i{E}_{ab}^i{P}_b^i=\sum _{a\ne b}{R}_{ab}^i.}$ Since the flux between any pair of amino acids a and b decreases when their difference of fitness increases, Halpern & Bruno, (1998) argued that the substitution rate Rⁱ is higher at position with higher sequence entropy. However, this argument is not rigorously valid, since it neglects the fact that the substitution rate is enhanced at sites where the site-specific selective favtors $F_a^{\mathrm{sel},i}{F}_b^{\mathrm{sel},i}$ are correlated with the global mutational flux $\left(P_a^{\mathrm{mut}}{P}_b^{\mathrm{mut}}{E}_{ab}^{\mathrm{mut}}\right)$. Consistent with this argument, we observed that the substitution rate is not a strictly increasing function of sequence entropy but, for the same sequence entropy, it tends to increase at sites that favor polar amino acids (see ‘Results’), which are characterized by higher mutational fluxes than hydrophobic amino acids.

Mutation process

Finally, we have to define the global exchangeability matrix $E_{ab}^{\mathrm{mut}}$ that characterizes the mutation process. For this, we consider four types of mutational models. To compare the resulting substitution rates, in all cases we fix the scale of the exchangeability matrix equating the substitution rate under mutation alone, ${\sum }_{a\ne b}{P}_a^{\mathrm{mut}}{P}_b^{\mathrm{mut}}{E}_{ab}^{\mathrm{mut}}=1$.

1. In the first model, the global exchangeability matrix is equal to the empirical exchangeability matrix (WAG, Whelan & Goldman, 2001; or JTT, Jones, Taylor & Thornton, 1992), i.e., $E_{ab}^{\mathrm{mut}}\equiv E_{ab}^{\mathrm{emp}}$. We call this model the empirical (emp) exchangeability matrix. Since empirical substitution processes include information both on mutation and selection, we expect that they strongly correlate with the selection process.

2. In the second model, we remove the effect of selection from the empirical substitution model by imposing that for each pair of amino acids, the flux predicted by the global model and averaged over all positions is equal to the empirical flux $P_a^{\mathrm{emp}}{P}_b^{\mathrm{emp}}{E}_{ab}^{\mathrm{emp}}$, which is the observational data from which empirical models are deduced: (12)${\displaystyle \left(P_a^{\mathrm{mut}}{P}_b^{\mathrm{mut}}{E}_{ab}^{\mathrm{flux}}\right)\frac{1}{L}\sum _i{F}_a^{\mathrm{sel},i}{F}_b^{\mathrm{sel},i}\frac{ln\left(F_b^{\mathrm{sel},i}\right)-ln\left(F_a^{\mathrm{sel},i}\right)}{F_b^{\mathrm{sel},i}-F_a^{\mathrm{sel},i}}=P_a^{\mathrm{emp}}{P}_b^{\mathrm{emp}}{E}_{ab}^{\mathrm{emp}}}$ where we use more compact matricial notation. We call the corresponding exchangebility matrix $E_{ab}^{\mathrm{flux}}$ the flux matrix (flux).This mutation model yields optimal results in phylogenetic inference (Arenas, Sanchez-Cobos & Bastolla, 2015).

3. Thirdly, we model the mutational process at the nucleotide level, using the genetic code and parameterizing the process through the nucleotide frequencies and the transition-transversion ratio κ. The four free parameters are fixed by imposing that the resulting background distribution $P_a^{\mathrm{mut}}$ yields amino acid frequencies as close as possible to those observed in the data, $P_a^{\mathrm{obs}}$ (Arenas, Sanchez-Cobos & Bastolla, 2015), as detailed below. We call the corresponding exchangeability matrix the optimized nucleotide (nuc_opt) matrix.

4. The last model is identical to the nuc_opt model, except that the nucleotide frequencies are not optimized but they are input parameters. In this way, we can vary the average hydrophobicity of the complete model by varying the Thymine content, since hydrophobic amino acids are enriched in the T base at second codon position. We call this model the nuc_var model.

In the nuc models, for any set of nucleotide frequencies and transition-tranversion rate we combine the substitution process at the nucleotide level with a selection process that assigns fitness one to sense codons and fitness zero to stop codons. Detailed balance is fulfilled at the nucleotide level, but it is only approximated at the codon level because of this selection against stop codons, therefore the transition to transversion rate can influence the stationary frequencies and we have to compute the stationary distribution of the 61 sense codons numerically.

More precisely, we model the mutation rate between two codons differing at one position, say the third one (n₁n₂n₃ and n₁n₂n₃′) as $\mu \kappa \left(n_3,n_3^{\prime }\right)f^{\mathrm{nuc}}\left(n_3^{\prime }\right)S\left(n_1{n}_2{n}_3\prime \right)$, where μ is a global rate parameter, $\kappa \left(n_3,n_3^{\prime }\right)$ is one if $n_3,n_3^{\prime }$ are related through a transversion and is the transition-tranversion rate otherwise, $f^{\mathrm{nuc}}\left(n_3^{\prime }\right)$ is the stationary frequency of the new nucleotide and S(n₁n₂n₃′) is zero if n₁n₂n₃′ is a stop codon, one otherwise. After the frequencies of the 61 sense codons evolve to their equilibrium state, the stationary frequencies of amino acids $P_a^{\mathrm{mut}}$ are computed summing over codons and the exchangeability matrix is computed from the equilibrium fluxes between pairs of codons that code for any pair of amino acids. In the nuc_opt model, the score of each set of mutation parameters is computed as the likelihood of the observed number of amino acids, ${\sum }_a{n}^{\mathrm{obs}}\left(a\right)log\left(P_a^{\mathrm{mut}}\right)$, and the parameters that maximize the likelihood are chosen.

Data and observed substitution rates

We performed our computations on 213 proteins that were examined in a previous study (Echave, Jackson & Wilke, 2015). The results were qualitatively identical from one protein to the other.

The observed substitution rates of 213 proteins that we show for comparison were estimated in (Echave, Jackson & Wilke, 2015) from the MSA of homologous sequences through the program Rate4Site (Pupko et al., 2002), which builds the phylogenetic tree using a neighbour-joining algorithm (Saitou & Nei, 1987) and estimates rates with an empirical Bayesian approach adopting the JTT model of sequence evolution (Jones, Taylor & Thornton, 1992). The multiple sequence alignments were generously provided by Julián Echave and are publicly available at the url https://github.com/wilkelab/therm_constraints_rate_variation.

Modelling stability against unfolded and misfolded states

Finally, for completeness we descibe here how we estimate the folding free energy ΔG of the experimentally known native state of a protein.

For this purpose, we adopt the contact matrix representation of the protein structure, consisting in the following: for each pair of residues at positions i and j along the polypeptidic chain, C_ij equals one if the residues are in contact and zero otherwise. We define two residues to be in contact if any pair of their heavy atoms are closer than 4.5 Å. Since contacts with |i − j| ≤ 2 are formed in almost all structures, they do not contribute to the free energy difference between the native and the misfolded ensemble, and we set C_ij = 0 if |i − j| ≤ 2. The free energy of a protein in the mesoscopic structure described by C_ij is modelled as a sum of contact interactions, E(C, A) = ∑_i<jC_ijU(A_i, A_j), which depends on the type of amino acids in contact A_i and A_j and on 210 contact interaction parameters U(a, b), for which we adopt the parameters determined in (Bastolla, Vendruscolo & Knapp, 2000).

For simplicity, we neglect the conformational entropy of the folded native state and estimate its free energy as $G_{\mathrm{nat}}\left(C^{\mathrm{nat}},A\right)\approx {\sum }_{i<j}{C}_{ij}^{\mathrm{nat}}U\left(A_i,A_j\right)$. Regarding the unfolded state, we neglect their contact interactions and estimate its free energy as G_U ≈  − TLS_U, where T is the temperature in units in which k_B = 1, L is chain length and S_U is the conformational entropy per residue of an unfolded chain. We compute the free energy of the misfolded state from the partition function of the contact energy E(C, A) over a set of compact contact matrices C of L residues that are obtained from the PDB. In agreement with previous studies (Garel & Orland, 1988; Shakhnovich & Gutin, 1989; Bryngelson et al., 1995), the resulting free energy is approximately described by the Random Energy Model (REM) (Derrida, 1981), with the addition of the third moment of the contact energy (Minning, Porto & Bastolla, 2013): (13)${\displaystyle G_{\mathrm{misf}}\equiv -Tlog\left(\sum _C{\mathrm{e}}^{-\sum _{i<j}{C}_{ij}U\left(A_i,A_j\right)∕T+S\left(C\right)}\right)\approx 〈E〉-\frac{〈{\left(E-〈E〉\right)}^2〉}{2T}+\frac{〈{\left(E-〈E〉\right)}^3〉}{6T^2}-L{S}_{C}T}$ where LS_C is the logarithm of the number of compact contact matrices, $〈.〉$ represents the average over the set of alternative compact contact matrices of L residues. This estimate only holds above the freezing temperature, while the free energy is kept constant below the freezing temperature (Derrida, 1981). We assume for simplicity that the conformational entropy, S(C_ij), is approximately the same for all compact structures including the native one, and it can be neglected for computing free energy differences. The mean values of the energy can be computed from the mean values of the contacts, which are computed at the beginning and tabulated to accelerate the computation: $〈E〉={\sum }_{i<j}〈C_{ij}〉U_{ij}$, $〈\left(E-〈E〉\right)〉={\sum }_{i<j,k<l}\left(〈C_{ij}{C}_{kl}〉-〈C_{ij}〉〈C_{kl}〉\right)U_{ij}{U}_{kl}$ with U_ij = U(A_i, A_j). We also adopt the approximation that $〈C_{ij}〉$ only depends on |i − j| (Minning, Porto & Bastolla, 2013).

Putting together these free energy estimates, we obtain the free energy difference between the native and the non-native states as (14)${\displaystyle \Delta G\left(C^{\mathrm{nat}},A\right)=G_{\mathrm{nat}}-kTlog\left({\mathrm{e}}^{-G_{\mathrm{misf}}∕kT}+{\mathrm{e}}^{-G_U∕kT}\right),}$ where the free energy of the non-native state is computed as a Boltzmann average, which is essentially equal to G_misf when the sequence is hydrophobic (G_misf − G_U∕kT ≪  − kT) and is essentially equal to G_U when the sequence is hydrophylic (G_misf − G_U∕kT ≫ kT). For neglecting stability against misfolding, we compute ΔG = G_nat(C^nat, A) + LS_U.

Dependence of the substitution rate on the global exchangeability model

In this work we studied two measures of the evolutionary variability of protein sites, sequence entropy and substitution rate, predicted through the site-specific stability constrained substitution models that we introduced and studied recently (Arenas, Sanchez-Cobos & Bastolla, 2015; Jimenez, Arenas & Bastolla, 2018). The results that we present arise from the predictions of our computational models, not from the analysis of natural protein sequences, thus they ignore important aspects of protein biology such as active sites and protein-protein interactions. We have to pay this price in order to address general questions such as the comparison between the two measures of evolutionary variability and how they are affected by the mutational process and by selection on protein stability.

The models adopted in this study are described in Methods and schematically represented in Fig. 1. They simulate an evolutionary process where mutations follow a global mutational process modelled through a global amino acid distribution $P_a^{\mathrm{mut}}$ and exchangeability matrix $E_{ab}^{\mathrm{mut}}$ (Fig. 1A), and they modify the stability of the protein ΔG, which determines fitness through Eq. (1) (Fig. 1B). We adopt the approximation that sites of the protein evolve independently under site-specific substitution processes globally governed by the fitness function of Eq. (1) (Fig. 1C). Since the stability ΔG depends on pairwise interactions, the fitness effect of a mutation at a site depends on the amino acids present at all other sites. In order to compute independent amino acid frequencies at each site, we perform two type of approximations: the mean-field approximation (Fig. 1D; Arenas, Sanchez-Cobos & Bastolla, 2015), which evaluates the effect of the mutation on stability considering the mean-field of the other sites, and the wild-type approximation (Fig. 1E; Jimenez, Arenas & Bastolla, 2018), which evaluates the effect of the mutation on stability when the other sites are occupied by the same amino acid as in the PDB sequence. We then compute site-specific substitution processes obtained by combining the global mutational process and the site-specific fixation probabilities computed from the site-specific amino acid frequencies through the Halpern-Bruno formula, Eq. (10).

We briefly report here previous results obtained with the stability-constrained models. Since the main component of contact interaction matrices like the one that we adopt here is hydrophobicity (i.e., U(a, b) ≈ ϵh(a)h(b), where h(a) is related with the hydrophobicity of amino acid a), the site-specific amino acid frequencies obtained under the stability-constrained model yield high average hydrophobicity at buried sites at which the native contact matrix has many contacts. More precisely, sites constrained to be highly hydrophobic have large components of the principal eigenvector of the contact matrix (Bastolla et al., 2005) or, almost equivalently, large effective connectivity (Bastolla et al., 2008). These structural descriptors are strongly related with the number of contacts but do not coincide with them. These stability constraints in turn influence the variability of the amino acid distribution (Porto et al., 2005) in such a way that sites that are constrained to have either very high or very low average hydrophobic (averaged over the site-specific amino acid distribution) are characterized by low entropy, while sites that can accomodate both polar and hydrophobic amino acids are characterized by high entropy. Because of this reason, the plot of the sequence entropy versus hydrophobicity and versus number of contacts has a bell shape.

These stability constraints that influence the sequence entropy also influence the site-specific substitution rates. However, we found that the substitution rates depend not only on the selective forces that act specifically at each protein site, but also on the global exchangeability matrix that represents the mutation process.

We considered three models of global exchangeability matrices (see ‘Materials and Methods’): (1) Empirical (emp) exchangeability matrices, such as the familiar WAG (Whelan & Goldman, 2001) or JTT (Jones, Taylor & Thornton, 1992) matrices; (2) Flux exchangeability matrices (flux), which are obtained from empirical exchangeability matrices removing the selective factors represented in the stability-constrained mean-field model, so that the average flux predicted by the model between any pair of amino acids coincides with the observed empirical flux, see Eq. (12); (3) Exchangeability matrices between amino acids obtained from a mutational process at the nucleotide level with parameters optimized by maximizing the likelihood of the observed amino acid composition (nuc_opt); (4) Exchangeability matrices obtained from a mutational process at the nucleotide level with varying parameters, that allows studying the effect of varying hydrophobicity (nuc_var).

We found that empirical exchangeability matrices (emp) produce the larger substitution rates (Figs. 2B and 2C ). These matrices take into account both the mutation process and the selection process, since they have been obtained from substitutions that have been fixed through natural selection. From Eq. (10) we can see that the substitution rate is enhanced when the site-specific selective factors $F_a^{\mathrm{sel},i}{F}_b^{\mathrm{sel},i}$ are correlated with the global mutational flux $\left(P_a^{\mathrm{mut}}{P}_b^{\mathrm{mut}}{E}_{ab}^{\mathrm{mut}}\right)$. Since the empirical substitution models were determined in such a way that their flux equals the flux observed in real data, which accounts both for the mutational process and for selection, we expect and find that the empirical flux is strongly correlated with the selective factors $F_a^{\mathrm{sel},i}{F}_b^{\mathrm{sel},i}$ averaged across all protein sites. This argument explains why the empirical exchangeability matrices yield the highest substitution  rates.

The flux exchangeability matrices remove from the empirical exchangeability matrix the effect of natural selection that is represented in the mean-field model averaged across sites. Consistently, we find that the substitution rates determined through the flux model are smaller than those determined with the emp model (Figs. 2B and 2C). We also found in previous work that the flux model yields larger likelihood in phylogenetic inference (Arenas, Sanchez-Cobos & Bastolla, 2015). Because of these results, the flux model is our default exchangeability model.

Figure 2C shows that the nuc_opt model with mutations at the nucleotide level and optimized parameters produces lower substitution rates than the flux model when associated with the MF model of selective constraints, which indicates that the flux model combined with the MF model may still include some selection. However, when the WT model of selection is applied, the nuc_opt model again produces lower substitution rates than the flux model for exposed sites with few contacts and high entropy, but the flux model produces lower substitution rates for buried sites with many contacts and low entropy (Figs. 2B and 2E). Note that the WT model represents stronger selective constraints than the MF model, since it generally predicts lower sequence entropies and substitution rates. Thus, these results suggest that the flux model associated with the WT model is effective in removing selective constraints for sites with many contacts, but less effective for sites with few contacts and high entropy.

Note that the curves that represent the substitution rates versus the sequence entropy tend to collapse for very high rates. This is due to the fact that all the global exchangeability matrices are normalized in such a way that their average flux equals one, ${\sum }_{ab}{P}_a^{\mathrm{mut}}{P}_b^{\mathrm{mut}}{E}_{ab}^{\mathrm{mut}}=1$. This flux is achieved at neutral sites where the selective factors $F_a^{\mathrm{sel},i}$ are equal and the entropy is maximal.

Substitution rates are different for hydrophobic and hydrophylic sites with the same entropy

Next, we investigated more in detail the relationship between site-specific sequence entropies and substitution rates. Since the flux between any pair of amino acids a and b, Eq. (10), decreases when their difference of fitness increases, Halpern and Bruno argued that the substitution rate Rⁱ is higher at position with higher sequence entropy (Halpern & Bruno, 1998). However, this argument is not rigorously valid, since it neglects the fact that the substitution rate is enhanced at sites where the site-specific selective factors $F_a^{\mathrm{sel},i}{F}_b^{\mathrm{sel},i}$ are correlated with the global mutational flux $\left(P_a^{\mathrm{mut}}{P}_b^{\mathrm{mut}}{E}_{ab}^{\mathrm{mut}}\right)$. One can see in Figs. 2B–2C that sites with larger entropy tend to have on the average larger substitution rates, as predicted by Halpern & Bruno (1998), but the substitution rate is not a strictly increasing function of sequence entropy, not even when it is averaged over different sites.

We then show in Fig. 3 the detailed plot of the substitution rate versus the sequence entropy for all sites of a small protein, chosen in such a way that we can spot all of the sites. We can clearly see in Fig. 3 two branches that correspond to different numbers of native contacts. Sites with few contacts, which tend to be occupied by polar amino acids, evolve faster than sites with many contacts, occupied by hydrophobic amino acids, even if their sequence entropy is equal. With the flux model of the mutation process, which we consider the most reliable model since it reproduces the empirically observed flux between all pairs of amino acids, this happens for both the MF and the WT model of natural selection (plots C and D). All other studied proteins present the same trend (see Figs. S1 and S2), but for large proteins the representation is less clear. When we model the mutation process at the codon level through the model nuc_opt, whose parameters are separately optimized for each protein, the differences between the two branches decrease considerably, in particular when we apply the WT model of selection (Fig. 3B) and the trend may change from one protein to the other, since different proteins evolve under different mutation processes, in such a way that the buried branch may evolve faster than the exposed branch for some protein, see Figs. S3 and S4.

Since sequence entropy is a measure of the selective constraints, this difference between sites with equal sequence entropy must be attributed to the mutation process embodied in the exchangeability matrix, not to natural selection. The explanation of this finding is based on Eq. (10): polar amino acids tend to have higher mutational fluxes $\left(P_a^{\mathrm{mut}}{P}_b^{\mathrm{mut}}{E}_{ab}^{\mathrm{mut}}\right)$ than hydrophobic amino acids, therefore exposed sites in which polar amino acids have larger site-specific selective factors $F_a^{\mathrm{sel},i}{F}_b^{\mathrm{sel},i}$ tend to evolve faster than buried sites.

This result contradicts the expectation that the site-specifi substitution rates are monotonic functions of the sequence entropy (Halpern & Bruno, 1998): one can see from Fig. 3 that sites characterized by lower entropy can substitute faster if they are exposed sites occupied by polar amino acids.

One may note that in Fig. 3 some white points with small number of contacts overlap with the black points with large number of contacts. This is due to the fact that the number of contacts is only an approximate predictor of the selective factors favoring hydrophobicity, while better correlation with the average hydrophobicity is achieved using the principal eigenvector of the contact matrix (Bastolla et al., 2005) or the effective connectivity profile that generalizes it for multidomain proteins (Bastolla et al., 2008), which are not merely local descriptors as the number of contacts but also represent the global topology of the native contact matrix.

More hydrophobic proteins substitute more slowly, but mutation bias towards hydrophobicity increases the substitution rates

After investigating the relationship between hydrophobicity and substitution rates comparing individual sites, we perform the same analysis comparing different proteins. For this purpose, we group the 213 proteins in our data set according to their predicted average hydrophobicity under the same mutational process and compare the substitution rates of groups characterized by different hydophobicity. In Fig. 4, each point represents a group of proteins with similar mean hydrophobicity. Each curve is obtained for a different mutation process with its background distribution P^mut and exchangeability matrix E^mut. One can see that, for the flux mutation process (black circles in Fig. 4), more polar proteins tend to evolve more rapidly, consistent with what we observed for polar sites in Fig. 3. Once again, this behavior may be explained by considering that more polar amino acids have higher mutational fluxes.

The above also holds for the nuc_opt model, in which the mutation process is separately optimized for each protein, but in this case when the WT selection model is used the maximum of the substitution rate is achieved for proteins of intermediate hydrophobicity (see red squares in Fig. 4B), consistent with the observation that with the nuc_opt model the hydrophobic sites may evolve slightly faster than the polar sites in some proteins characterized by a mutation process that favors hydrophobic residues, see Figs. S3 and S4).

We then consider the same mutation process for all proteins, parameterized by the G+C content at the mutational equilibrium (nuc_var model). Since Thymine at second codon position almost always codes for hydrophobic amino acids, there is a negative correlation between G+C content of the mutation model and the average hydrophobicity of the protein sequence. Varying the mutation bias we construct different sets of model proteins that present varying hydrophobicity and are characterized by different mutational fluxes. In this way, we can investigate how the mutation bias influcence the biophysical properties (hydrophobicity) and the evolutionary properties (substitution rate, sequence entropy) of an evolving protein. When the GC content is high and the hydrophobicity is low (GC ≥0.35 with the MF model, Fig. 4A and GC ≥0.5 with the WT model, Fig. 4B) more polar proteins tend to evolve faster, as we observed with the flux model. However, when the mutation process induces high hydrophobicity (GC = 0.2 in Fig. 4A and GC <0.5 in Fig. 4B), the substitution rate becomes an increasing function of hydrophobicity. This is easily rationalized by the fact that, when the background distribution $P_a^{\mathrm{mut}}$ is biased towards hydrophobic amino acids, the mutational flux $P_a^{\mathrm{mut}}{P}_b^{\mathrm{mut}}{E}_{ab}^{\mathrm{mut}}$ is higher between pairs of hydrophobic residues, and the other way round.

Different G+C content in Fig. 4 represent different mutational processes, which may be interpreted as bacterial species characterized by different GC bias. For the MF model (plot A), one can see that mutational processes with extreme bias (very high or very low G+C content and hydrophobicity) tend to increase the average substitution rate, while for the WT model (plot B) mutational processes biased towards hydrophobic residues such as GC = 0.2 have higher substitution rates. Therefore, although hydrophobicity is nega- tively correlated with the substitution rate when we compare different proteins evolving with the flux mutational process (black points in Figs. 4 and Figs. 3C and 3D) or with mutational processes with high G+C, the correlation becomes positive if we compare different mutational processes.

Influence of the mutation bias on sequence entropies

Next, we study how the shape of the entropy-hydrophobicity curve depends on the mutation bias. In Fig. 5 each point represents a set of protein sites with similar hydrophobicity in the stability-constrained evolutionary model. The sequence entropy has an almost universal bell shape as a function of the average hydrophobicity of the site, with a maximum when the average hydrophobicity is approximately 0.14, which is the average hydrophobicity of the uniform distribution of amino acids. This result is of course not surprising, since sites with very high or low average hydrophobicity have distributions that favor only the most hydrophobic or polar amino acids, while all amino acids are allowed when the average hydrophobicity equals the one of the uniform distribution.

Changes in the mutation bias do not change the position of the maximum, but they strongly shift the sequence entropy curves in the vertical direction. The largest entropies are obtained for the mutation model nuc_opt optimized separately from each PDB sequence (thick black line) and for the mutation bias with G+C content equal to 0.40, which has a small bias towards slightly hydrophobic sequences. Extreme mutation bias both towards hydrophobic (low GC) and polar (high GC) amino acids yield very reduced sequence entropies, which means that selection must impose stronger constraints in order to preserve the average hydrophobicity needed for stable proteins. This result is consistent with the finding that, for equal population size, the average fitness achieved in evolution has a maximum as a function of the mutation bias, and it is low for extreme mutation bias either toward hydrophobic or towards hydrophylic sequences (Mendez et al., 2010; Bastolla, Dehouck & Echave, 2017).

Influence of hydrophobicity on substitution rates

We now study the relationship between site-specific hydrophobicity and site-specific substitution rate. As in the previous figure, also in Fig. 6 each point represents a set of protein sites with similar average hydrophobicity in the stability-constrained evolutionary model, and we plot the average substitution rate versus the average hydrophobicity. Different from the shape of the sequence entropy, the shape of the substitution rate curve clearly depends on the mutation bias. The average hydropobicity at which the maximum substitution rate is achieved decreases with the G+C content or, equivalently, it increases with the average hydrophobicity of the mutation process. In other words, when the mutation process favors the exchange of hydrophobic amino acids, the maximum of the substitution rate is achieved at sites that are more hydrophobic, as expected on the basis of Eq. (10) that suggests that sites where the site-specific selective factors are correlated with the mutational flux have higher substitution rates. This result confirms that the mutation process has a strong influence on the substitution rates.

Consistent with Fig. 4, the substitution rate at the maximum tends to increase for mutation processes that favor higher hydrophobicity (lower G+C bias), but for the MF model (plot A) they also increase for extremely polar mutation bias (G+C content 0.8). Consistent with the results reported in Fig. 2, the flux model of the exchangeability matrix (thick black line) predicts higher substitution rates than the nuc_opt model (red) when applied together with the MF model (Fig. 4A), but when it is applied together with the WT model it predicts lower substitution rates at hydrophobic sites with many contacts (Fig. 4B), suggesting that the WT model is effective at removing the effect of selection from the flux model at these buried sites.