A thermodynamic model of protein structure evolution explains empirical amino acid substitution matrices

Abstract Proteins evolve under a myriad of biophysical selection pressures that collectively control the patterns of amino acid substitutions. These evolutionary pressures are sufficiently consistent over time and across protein families to produce substitution patterns, summarized in global amino acid substitution matrices such as BLOSUM, JTT, WAG, and LG, which can be used to successfully detect homologs, infer phylogenies, and reconstruct ancestral sequences. Although the factors that govern the variation of amino acid substitution rates have received much attention, the influence of thermodynamic stability constraints remains unresolved. Here we develop a simple model to calculate amino acid substitution matrices from evolutionary dynamics controlled by a fitness function that reports on the thermodynamic effects of amino acid mutations in protein structures. This hybrid biophysical and evolutionary model accounts for nucleotide transition/transversion rate bias, multi‐nucleotide codon changes, the number of codons per amino acid, and thermodynamic protein stability. We find that our theoretical model accurately recapitulates the complex yet universal pattern observed in common global amino acid substitution matrices used in phylogenetics. These results suggest that selection for thermodynamically stable proteins, coupled with nucleotide mutation bias filtered by the structure of the genetic code, is the primary driver behind the global amino acid substitution patterns observed in proteins throughout the tree of life.


= 2.57
The protein misfolding-fitness function behaves as the sigmoidal stability-fitness function with a modified offset Another fitness model was proposed by Drummond and Wilke 3 . They found that a major selection pressure is selection against cytotoxicity caused by protein misfolding. In their model the fraction of unfolded protein has an exponentially decreasing effect on fitness, which depends on a toxicity parameter c and the protein abundance parameter A.
When sΔG is low (sΔG<-3) the misfolding fitness function can be simplified to the fraction folded fitness function although offset by log(cA): Where ωfolding is the fitness function used throughout the paper:

Estimation of parameter uncertainty by bootstrap analysis
To estimate the uncertainty of the TSM parameters we performed a bootstrap analysis.
Specifically, we compute the phylogenetic log-likelihood (IQ-TREE) for each parameter combination explored in Fig. S1 for each of the 52 sequence alignments individually. Next, using random sampling with replacement (sample size = 52) over the alignment dataset, we took 10,000 bootstrap samples, finding the argmax parameter set for each: From these samples, we computed the standard deviation of each parameter.

Q-matrix estimation error upon noise injection
Deviations between inferred rates and rates predicted with the TSM-model could originate from ΔΔG prediction errors. The magnitudes of ΔΔG prediction errors were estimated from the correlation between predicted ΔΔG values from Rosetta and known experimental values for the same 52 proteins described in the main text (σerr=1.7 kcal/mol based on 590 measurements) 4 . To estimate the impact of such errors, we predicted rates with the TSM model as described in the main text but modified the ΔΔG values by adding errors sampled from the empirical error distribution.
The expected correlation between Q-matrix parameters with and without sampled error is rπ+R 2 =0.79. Thus, ΔΔG prediction error alone likely contributes at least 21% of the variance in the LG matrix that is unexplained by our model. We note that the true error could be greater, as half of the substitutions in the experimental ΔΔG dataset are to alanine, and alanine substitutions are more accurately predicted than other mutations (σerr,X->A=1.4 vs σerr,X->!A=1.9 kcal/mol).

Fig. S4. The TMS model recapitulates the mean substitution behavior between amino acids.
A) Comparison of the amino acid exchangeability matrix predicted by TSM (black circles) and the WAG global exchangeability matrix (red circles). Inset, correlation between exchangeabilities from TSM (x-axis) vs WAG (y-axis). B) Correlation between amino acid equilibrium probabilities predicted by TSM (x-axis) and values from WAG (y-axis). C) The same as A, but without selection. D) the same as B but without selection.   Starting from the optimal parameters of the TMS matrix, we expanded the model by replacing the K80 nucleotide proposal model with a T92 model. The figure shows that the K80 assumption of no GC-bias (i.e., %-GC=0.5) results in exchangeabilities and stationary frequencies that correlate optimally with the LG matrix parameters. Other values for GC-bias correlate less well. shows likelihoods for the MSAs with a delta log likehood of less than 5. Out of 500 MSAs, our Q-matrix achieved a higher likelihood than LG for 8.