Maximum-Likelihood Phylogenetic Inference with Selection on Protein Folding Stability

Abstract

Despite intense work, incorporating constraints on protein native structures into the mathematical models of molecular evolution remains difficult, because most models and programs assume that protein sites evolve independently, whereas protein stability is maintained by interactions between sites. Here, we address this problem by developing a new mean-field substitution model that generates independent site-specific amino acid distributions with constraints on the stability of the native state against both unfolding and misfolding. The model depends on a background distribution of amino acids and one selection parameter that we fix maximizing the likelihood of the observed protein sequence. The analytic solution of the model shows that the main determinant of the site-specific distributions is the number of native contacts of the site and that the most variable sites are those with an intermediate number of native contacts. The mean-field models obtained, taking into account misfolded conformations, yield larger likelihood than models that only consider the native state, because their average hydrophobicity is more realistic, and they produce on the average stable sequences for most proteins. We evaluated the mean-field model with respect to empirical substitution models on 12 test data sets of different protein families. In all cases, the observed site-specific sequence profiles presented smaller Kullback-Leibler divergence from the mean-field distributions than from the empirical substitution model. Next, we obtained substitution rates combining the mean-field frequencies with an empirical substitution model. The resulting mean-field substitution model assigns larger likelihood than the empirical model to all studied families when we consider sequences with identity larger than 0.35, plausibly a condition that enforces conservation of the native structure across the family. We found that the mean-field model performs better than other structurally constrained models with similar or higher complexity. With respect to the much more complex model recently developed by Bordner and Mittelmann, which takes into account pairwise terms in the amino acid distributions and also optimizes the exchangeability matrix, our model performed worse for data with small sequence divergence but better for data with larger sequence divergence. The mean-field model has been implemented into the computer program Prot_Evol that is freely available at http://ub.cbm.uam.es/software/Prot_Evol.php.

Keywords: folding stability; maximum-likelihood estimate; misfolded state; structurally constrained substitution models.

Fig. 1.

Site-specific average hydrophobicity (left) and entropy (right) of the MF distributions as a function of the number of native contacts for the protein with PDB code 153l. As expected, there is a very strong correlation between hydrophobicity and number of contacts and the entropy reaches a maximum at an intermediate number of contacts.

Fig. 2.

Left: Log likelihood of various MF models as a function of the log likelihood of the purely mutation model. Each point represents a protein. Right: Mean log likelihood of the five types of MF models. The plotted statistical errors show that differences are significant except for the two rightmost bars.

Fig. 3.

Left: Average hydrophobicity of the MF models versus the average hydrophobicity of the PDB sequence. Each point represents a protein. Right: Average folding free energy (native minus misfolded) $\overline{\Delta G}$ of the MF models versus the average folding free energy of the PDB sequence. $\overline{\Delta G}<0$ means that the MF model describes on the average stable proteins. Each point represents a protein.

Fig. 4.

Mean log-likelihood of the proteins in the test set with respect to the model $P_a^{\text{MF}2,i}$ versus the temperature in arbitrary units set by our contact interaction energy function.

Fig. 5.

Difference of KL divergence from the observed amino acid profile between the empirical model and the MF model (KLDobs_emp–KLDobs_mf) for the 12 studied protein families, under different conditions on the minimum sequence identity allowed. Positive differences mean that the observed profile agrees better with the MF model than with the empirical model.