From residue coevolution to protein conformational ensembles and functional dynamics

Contact Predictions via Boltzmann Learning.

Entropy maximization provides a simple procedure for building a probabilistic model that is consistent with a set of available measures. If we know the average values {Fi} of a set of variables {xi}, the associated maximum entropy distribution is given by P∝exp(∑iλixi), with a Lagrange multiplier λi for each variable xi (40). Fixing the frequencies for single and pairs of amino acids in a multiple sequence alignment (MSA), P takes the form (11): P({a})=Z−1exp[∑ihi(ai)+∑i<jJi,j(ai,aj)] over protein sequences. The parameters hi(α) and Ji,j(α,β) are the Lagrange multipliers that fix the averages of the model, fi(α) and fi,j(α,β), to the empirical frequencies Fi(α) and Fi,j(α,β) computed from the MSA, where α and β denote two particular amino acids and i, j two particular residues along the protein sequence. Due to the Boltzmann-like form of the previous equation, the parameter Ji,j(α,β) can be interpreted as the direct interaction between amino acids α and β at positions i and j, after the contributions from the interaction with other positions through indirect pathways have been disentangled (11, 12, 22).

To our knowledge, since the early work of Lapedes et al. (11), the numerical route of likelihood maximization via importance sampling, or Boltzmann learning (27), has not been explored, probably due to its computational complexity. As outlined in the introduction, this approach [see, e.g., Roudi et al. (41) for an extensive comparison between Boltzmann learning and various approximated schemes] has the advantage of having unbounded precision in retrieving the parameters of a maximum-entropy model. We tested Boltzmann learning on an assorted set of 18 proteins with varying length, from 63 to 216 residues, and different secondary structure composition (Table 1). For each alignment, we inferred a pairwise model by maximizing a regularized version of the log-likelihood of the sequences with respect to parameters {h} and {J} (see Materials and Methods for details). In all 18 cases, we were able to reproduce the empirical frequencies {F} within a reasonable error, as we checked through extensive sampling from the final models (SI Appendix, Fig. S1). Depending on the size of the protein sequence, we obtained mean absolute relative errors 〈|Fi,j−fi,j|/Fi,j〉 between the model pair frequencies fi,j and the empirical Fi,j ranging from 1% to 3% (0.2–2% for Fi,j>0.01) for the different protein families. Being that our approach is based on importance sampling, the presence of many isolated modes in the distribution of sequences could lead to poor mixing of the Markov chain and, consequently, to a large error in the estimate of the gradient of the likelihood function. Indeed, clusters of sequences in multiple sequence alignments are common and are known to reflect potential functional subfamilies among the members of a single protein family (7). As a cross-check, we verified that the fitted models capture the organization of the original alignment in clusters of sequences (Fig. 1 and SI Appendix). We point out that “external” sources of variation—such as changes in functional requirements within the same protein family—should be explicitly taken into account in future, improved models to discriminate between intrinsic and extrinsic correlations in the sequence alignment.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

(A and B) Energy surfaces obtained projecting sequences from the ADH_zinc_N domain family (A) and sequences sampled from the model (B) over the first two principal components of the MSA (7), and taking the negative logarithm of the resulting distribution. High-probability regions in sequence space are in dark blue. The cluster structure of the alignment is clearly reproduced by the simulated trajectories. (C and D) Mean precision of the top-ranked predictions for different values of the scaled rank (rank/total number of contacts), and the mean true positive rate for different values of the false-positive rate (ROC curve). The color bands show the SD and the interval between the minimum and maximum values.

For each protein, from the estimated set of {J} we computed a matrix of coevolutionary couplings using a protocol first proposed by Ekeberg et al. (25). Assuming that a strong, direct coevolution is an evidence of physical contact between two residues, we finally ranked the pairs of residues using the value of coevolutionary coupling as a score. The accurate determination of couplings through the Boltzmann learning algorithm resulted in high-quality predictions. Fig. 1 summarizes the performance in terms of mean precision against the top-scoring predictions rank, and as mean true positive rate as a function of false positive rate. We defined a pair of residues to be in contact when the distance between their C_β atoms (C_α in case of GLY) is smaller than 8 Å (42, 43), and included in the analysis all of the pairs with a sequence separation larger than four. For the top N predictions, where N denotes the number of amino acids in a protein structure, the algorithm obtained a mean precision of 0.67±0.03, with a maximum of 0.89 for the cNMP_binding domain (PDB ID code 3FHI). A comparison with predictions obtained through a more standard mean-field solution (18, 19) is included in SI Appendix.

A Coarse-Grained Model for Structure Prediction and Conformational Sampling.

To take full advantage of the accurate contact predictions resulting from the Boltzmann learning algorithm, we propose a coarse-grained model that combines an all-heavy-atom description of the protein backbone with a Cβ description of the side chains (Materials and Methods). Similar coarse-grained approaches have already been successfully applied to study the thermodynamics of model proteins (30). Because we expect the residues to be evolutionary coupled through their side chains, we set the predicted interactions to act between the Cβ atoms—i.e., the first atom to branch out of the main chain. More precisely, in our model we used the N best-predicted contacts, where N is the number of residues (Materials and Methods). Moreover, the inclusion of all of the heavy atoms of the main chain allows use of a transferable potential that acts on the actual degrees of freedom of the backbone enabling the correct reproduction of the experimental population of Ramachandran backbone angles. To complement long-range predictions, we estimated the helical propensity of each residue solely from the predicted contacts and translated it in a stabilizing hydrogen bond-mimicking interaction between residues i,i+4 (see SI Appendix for details). For all of the α and α/β proteins, all of the helices are correctly predicted with the exception of h1 for 3D7I, h5 for 2GJ3, and h2 for 5P21 structures (SI Appendix, Fig. S2).

Using this simplified model, we investigated the conformational space accessible to the protein domains in Table 1. The nine best-predicted structures obtained in the folding runs are shown in Fig. 2 superimposed to their respective native folds. Those structures correspond to the conformations with minimum distance RMSD of the α carbon atoms (Cα-dRMSD; Materials and Methods) to the native conformation. For all of them, not only is the global fold correctly predicted, but also most of the secondary structure elements are present and correctly packed.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

The nine best-predicted protein structures of the 18 simulated are shown in blue, overlaid to the native conformation in cyan (transparent). The proteins are disposed with increasing size from left to right and from top to bottom. Their Pfam identifier is shown below.

For domains of similar size, we found a clear correlation between the precision of the contact predictions and the final dRMSD values (e.g., domains Trans_Reg_C, CMD, and PAS have lower precision in contact prediction and higher dRMSD values than other domains with a similar number of amino acids; Table 1), confirming that the precision of contacts prediction is the main determinant of the quality of the final fold reconstruction. However, all 18 proteins fold within 3 Å to the native protein, except for the 2GJ3 structure (PAS domain), which performs worse (Table 1 and SI Appendix, Fig. S6). We found that the model is tolerant to both a conspicuous presence of nonnative contacts among the set of predicted contacts as well as to the absence of a large proportion of native contacts among the predicted ones.

To have a reference baseline to compare with, we also simulated the proteins using the same conditions but replacing the set of predicted contacts with the set of native contacts (see details in SI Appendix). The best dRMSDs obtained in this case are always below ∼2 Å (SI Appendix, Fig. S6, black curve). These values quantify the theoretical limit of the present model with perfectly predicted contacts. The quality of the predicted structures is as good as those predicted by the structure-based reference simulation in several cases (HTH_31, Sigma70_r2, RRM_1, Ras). Interestingly, the dRMSD of the minimum energy structure (dashed curve in SI Appendix, Fig. S6) does not deviate much from the absolute minimum dRMSD structure sampled in the trajectory. Indeed, for 17 of 18 proteins, the minimum energy structure is below 4 Å to the native structure, which indicates that the model and the contact predictions are sound and lead to a funnel-shaped energy landscape whose minimum corresponds to the crystallographic structure.

In the single case of the PAS domain, we observe a larger deviation, 7 vs. 3.6 Å for the dRMSD of the minimum energy structure and the absolute minimum, respectively (Table 1), which can be only partially explained by the presence of an associated cofactor in the experimental structure. Such a large difference indicates that though the minimum dRMSD structure still belongs to the native basin, it does not coincide with the minimum energy structure, which is the defining property of the native fold. Inspecting the set of predicted contacts, we found that five of them are in the dimeric interface of the corresponding homodimer (SI Appendix, Fig. S7). After removing these interdomain contacts, the minimum energy dRMSD of the monomeric structure decreases to 4.7 Å (SI Appendix, Fig. S6), showing that these few contacts were responsible for a large displacement from the native conformation. This result is also corroborated by the fact that when five randomly chosen false-positive contacts are removed, the dRMSD of the minimum energy structure does not improve, which we verified running 20 independent simulations (SI Appendix). As far as we know, the effect on fold reconstruction of strongly coupled pairs of residues at homodimeric (or homooligomeric) interfaces, when incorrectly classified as contacts in the monomeric structure, have not been described before. Here we show that this effect can be large, depending on the relative position of the pairs of residues at the dimer interfaces. Reasonably, similar but weaker effects are present for other cases, being that 10 over 18 proteins in our set were homodimeric in the corresponding PDB structure. The analysis demonstrates that the ability of PAS domains to form homodimers (44) is subjected to evolutionary pressure, and identifies a set of five pairs of positions involving the N-terminal helix (I40-F27; Y49-F27; A117-P35; L130-Q29; M132-A34 in the 2GJ3 PDB numbering) that emerge as crucial for the stabilization of the PAS homodimer across the protein family, as supported by their proximity in the crystal structure, strong coevolutionary coupling, and simulation.

Conformational Heterogeneity and Residue Coevolution.

To investigate the ability of coevolutionary couplings to also encode for dynamical and functional information, we analyzed the conformational space close to the native fold in the case of the catalytic domain (CD) of SRC tyrosine kinase and characterized the full folding reaction of the Ras domain from a thermodynamical point of view.

Protein kinases are known to undergo a large conformational rearrangement of a centrally located loop during activation, called activation loop or A-loop. In the case of the CD of the SRC tyrosine kinase, the A-loop spans more than 20 residues, and the structural rearrangement moves the backbone across ≈25 Å from the inactive conformation where the A-loop is folded, to the fully active structure where the A-loop is in an extended conformation. The A-loop is known to be flexible to the point of being invisible in many crystal structures. It is thus interesting to see if traces of this conformational transition can be observed by an exhaustive sampling of the native fold basin with our model. We performed a 500-ns-long parallel tempering (PT) simulation of the 250 residues CD starting from the inactive conformation with the A-loop in the so-called “half-closed conformation,” corresponding to the structure with PDB code 2SRC. We used the 250 best-predicted contacts calculated over 3,812 effective sequences from PFAM Pkinase_Tyr family, without adding any structure-based bias. Indeed, only 130 of the 250 contacts are native contacts (where a native contact is defined between CB atoms within 8 Å in the native structure). Consistent with the fact that no contacts have been predicted for the A-loop, we observe an extensive sampling of the conformational space available to the loop while the rest of the protein correctly maintains the native fold.

In Fig. 3 we compare the experimental and simulated structural flexibility of different positions along the chain. Even though we expect a sequence-dependent effect on structural order/disorder, we note that flexible regions (A-loop, 298–303 loop and αG-helix) are captured by our model with few exceptions (αC-helix, G-rich loop), suggesting that chain flexibility itself is partially encoded in the sequences of the kinase protein family (SI Appendix, Fig. S8). Indeed, these regions are known to be involved in the activation or in protein–protein interaction (45). We note that flexibility cannot be trivially deducted from the number of predicted contacts involving each residue, as shown in SI Appendix, Fig. S8, but is rather a consequence of the whole fold and network of interactions. Moreover, both the active and inactive conformations of the A-loop are repeatedly sampled within 5 Å Cα-RMSD to the crystal structure (SI Appendix, Fig. S8). The ability to reach both endpoints of the complex conformational transition in the catalytic domain is very encouraging.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Comparison of the fluctuations of the SRC catalytic domain. The experimental X-ray structure (PDB ID code 2SRC) is shown (Left) with the normalized B-factor color and thickness coded. The same structure with the B-factor calculated from the simulation of the SRC catalytic domain (Right). For easier comparison, the scales have been normalized. Higher values correspond to larger fluctuations.

The Ras protein is a α/β globular protein and is a crucial mediator of cellular proliferation and differentiation involved in several signaling pathways (46, 47). We performed a 180-ns-long PT simulation to explore a wide range of temperatures and to have a solid sampling. The simulated folding transition is highly cooperative with a clear peak in the heat capacity at constant volume at the folding temperature Tf (SI Appendix, Fig. S9). In Fig. 4, we show the free energy as a function of dRMSD at three different temperatures: Tlow<Tf, T≈Tf, and Thigh>Tf. In the unfolded state (U), the β-strands β2 and β3 are unfolded and lead to an extended tail departing from a rather structured core around the partially formed α3 and α4 helices. The collapse of this unfolded tail and its correct positioning lead to an intermediate, partially folded state (I, dRMSD = 6.5 Å), where helices α4 and α3 are formed. Eventually, the native basin (N, dRMSD = 3.8 Å) is reached with the formation of helices α3 and α5 and their correct packing by crossing a free-energy barrier of 4 kJ/mol. These features are not simply encoded in the native fold, because a coarse-grained Cα structure-based model (48) is unable to capture them (SI Appendix). It is interesting to note that the existence of a folding intermediate has been reported at high pressure and in denaturing conditions (49). Though we cannot directly compare our result with the experimental structures and the different experimental conditions, it is promising to observe the emergence of typical folding features, such as intermediate folding states, from a genomic-derived coarse-grained model.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Folding free energy of Ras protein at three different temperatures around the folding temperature Tf=317K as a function of the Cα-dRMSD. Representative structures of the native (N), intermediate (I), and unfolded (U) states are shown above with the secondary structure elements labeled.