Normal mode analysis (NMA) methods are widely used to study dynamic aspects of protein structures. Two critical components of NMA methods are coarse-graining in the level of simplification used to represent protein structures and the choice of potential energy functional form. There is a trade-off between speed and accuracy in different choices. In one extreme one finds accurate but slow molecular-dynamics based methods with all-atom representations and detailed atom potentials. On the other extreme, fast elastic network model (ENM) methods with Cα−only representations and simplified potentials that based on geometry alone, thus oblivious to protein sequence. Here we present ENCoM, an Elastic Network Contact Model that employs a potential energy function that includes a pairwise atom-type non-bonded interaction term and thus makes it possible to consider the effect of the specific nature of amino-acids on dynamics within the context of NMA. ENCoM is as fast as existing ENM methods and outperforms such methods in the generation of conformational ensembles. Here we introduce a new application for NMA methods with the use of ENCoM in the prediction of the effect of mutations on protein stability. While existing methods are based on machine learning or enthalpic considerations, the use of ENCoM, based on vibrational normal modes, is based on entropic considerations. This represents a novel area of application for NMA methods and a novel approach for the prediction of the effect of mutations. We compare ENCoM to a large number of methods in terms of accuracy and self-consistency. We show that the accuracy of ENCoM is comparable to that of the best existing methods. We show that existing methods are biased towards the prediction of destabilizing mutations and that ENCoM is less biased at predicting stabilizing mutations.
Biological macromolecules are dynamic objects. In the case of proteins, such movements form a continuum ranging from bond and angle vibrations, sub-rotameric and rotameric side-chain rearrangements [1], loop or domain movements through to folding. Such movements are closely related to function and play important roles in most processes such as enzyme catalysis [2], signal transduction [3] and molecular recognition [4] among others. While the number of proteins with known structure is vast with around 85K structures for over 35K protein chains (at 90% sequence identity) in the PDB database [5], our view of protein structure tends to be somewhat biased, even if unconsciously, towards considering such macromolecules as rigid objects. This is due in part to the static nature of images used in publications to guide our interpretations of how structural details influence protein function. However, the main reason is that most known structures were solved using X-ray crystallography [6] where dynamic proprieties are limited to b-factors and the observation of alternative locations. Despite this, it is common to analyze larger conformational changes using X-ray structures with the comparison of different crystal structures for the same protein obtained in different conditions or bound to different partners (protein, ligand, nucleic acid). It is particularly necessary to consider the potential effect of crystal packing [7], [8] when studying dynamic properties using X-ray structures. Nuclear magnetic resonance (NMR) is a powerful technique that gives more direct information regarding protein dynamics [9], [10]. Different NMR methodologies probe distinct timescales covering 15 orders of magnitude from 10−12 s side chain rotations via nuclear spin relaxation to 103 s using real time NMR [10]. In practice, there is a limitation on the size of proteins that can be studied (between 50–100 kDa) although this boundary is being continuously pushed [11] providing at least partial dynamic information on extremely large systems [12]. However, only a small portion (around 10%) of the available proteins structures in the Protein Data Bank (PDB) are the result of NMR experiments [5]. Molecular dynamic simulations numerically solve the classical equations of motion for an ensemble of atoms whose interactions are modeled using empirical potential energy functions [13]–[15]. At each time step the positions and velocities of each atom are calculated based on their current position and velocity as a result of the forces exerted by the rest of the system. The first MD simulation of a protein (Bovine Pancreatic Trypsin Inhibitor, BPTI) ran for a total 8. 8 ps [16] followed by slightly longer simulations (up to 56 ps) [17]. A film of the latter can be seen online (http: //youtu. be/_hMa6G0ZoPQ). Despite using simplified potentials and structure representations (implicit hydrogen atoms) as well as ignoring the solvent, these first simulations showed large oscillations around the equilibrium structure, concerted loop motions and hydrogen bond fluctuations that correlate with experimental observations. Nowadays, the latest breakthroughs in molecular dynamics simulations deal with biological processes that take place over longer timescales. For example, protein folding [18], transmembrane receptor activation [19], [20] and ligand binding [21]. These simulations require substantial computer power or purpose built hardware such as Anton that pushes the current limit of MD simulations to the millisecond range [22]. Despite powerful freely available programs like NAMD [23] and GROMACS [24] and the raise of computational power over the last decade, longer simulations reaching timescales where most biological processes take place are still state-of-the-art. Normal modes are long established in the analysis of the vibrational properties of molecules in physical chemistry [25]. Their application to the study of proteins dates back to just over 30 years [26]–[30]. These earlier Normal Mode Analysis (NMA) methods utilized either internal or Cartesian coordinates and complex potentials (at times the same ones used in MD). As with earlier MD methods their application was restricted to relatively small proteins. Size limitations notwithstanding, these early studies were sufficient to demonstrate the existence of modes representing concerted delocalized motions, showing a facet of protein dynamics that is difficult to access with MD methods. Some simplifications were later introduced and shown to have little effect on the slowest vibrational modes and their utility to predict certain molecular properties such as crystallographic b-factors. These simplifications included the use of a single-parameter potential [31], blocks of consecutive amino acids considered as units (nodes) [32] and the assumptions of isotropic [33] fluctuations in the Gaussian Network Model (GNM) or anisotropic fluctuations [34]. These approximations have drastically reduced the computational time required, thus permitting a much broader exploration of conformational space using conventional desktop computers in a matter of minutes. Of these, the most amply used method is the Anisotropic Network Model (ANM) [35], [36]. ANM is often referred simply as an elastic network model; one should however bear in mind that all normal mode analysis methods are examples of elastic network models. ANM uses a simple Hook potential that connects every node (a point mass defined at the position of an alpha carbon), within a predetermined cut-off distance (usually 18 Å). More recently, a simplified model, called Spring generalized Tensor Model (STeM), that uses a potential function with four terms (covalent bond stretching, angle bending, dihedral angle torsion and non-bonded interaction) has been proposed [37]. The normal mode analysis of a macromolecule produces a set of modes (eigenvectors and their respective eigenvalues) that represent possible movements. Any conformation of the macromolecule can in principle be reached from any other using a linear combination of amplitudes associated to eigenvectors. It is essential however to not loose sight of the limitations of normal mode analysis methods. Namely, normal modes tell us absolutely nothing about the actual dynamics of a protein in the sense of the evolution in time of atomic coordinates. Plainly speaking, normal mode analysis is informative about the possible movements but not actual movements. Additionally, normal modes tell us of the possible movements around equilibrium. These two caveats clearly place normal mode analysis and molecular dynamics apart. First, molecular dynamics gives an actual dynamics (insofar as the potential is realistic and quantum effects can be ignored). Second, while the equilibrium state (or the starting conformation) affects the dynamics, one can explore biologically relevant timescales given sufficient computational resources to perform long simulations. The vast majority of coarse-grained NMA models only use the geometry of the protein backbone (via Cα Cartesian position) disregarding the nature of the corresponding amino acid, in doing so a lot of information is lost. To our knowledge, there have been three independent attempts at expanding coarse-grained NMA models over the years to include extra information based on backbone and side-chain atoms. Micheletti et al. [38] developed the βGM model in which the protein is represented by Cβ atoms for all residues except Glycine in addition to the Cα atoms. The Hamiltonian is a function exclusively of Cα and Cβ distances. As a Gaussian model, βGM does not give information about directions of movement but only their magnitude and can be used solely to predict b-factors. As Cβ atoms do not change position, this model cannot be used by definition to predict the effect of mutations, either on dynamics or stability. The authors report results on b-factor prediction comparable to GNM. Lopéz-Blanco et al. [39] developed a NMA model in internal coordinates with three different levels of representation: 1. Heavy-atoms, 2. five pseudo-atoms (backbone: NH, Cα, CO and side-chain: Cβ and one at the center of mass of the remaining side chain atoms) and 3. Cα representation. While the potential is customizable, the default potential uses a force constant that is distance dependent but atom type independent. The method is validated through overlap analysis on a dataset of 23 cases. The authors report no significant differences in overlap for the different representations. Lastly, Kurkcuoglu et al. [40] developed a method that mixes different levels of coarse-graining. The authors test the method on a single protein, Triosephosphate Isomerase with higher atomic representation for a loop and achieve better overlap values compared to Cα only representation. All the methods above, while adding more detail to the representation utilize force constants that are not atom-type dependent. Therefore, while less coarse-grained, all the methods above are still atom-type and amino-acid type agnostic. By definition, irrespective of the level of coarse-graining, such models cannot account for the effect of mutations on protein dynamics or stability. It has been shown that different amino acids interact differently and that single mutations can have a high impact on protein function and stability [41]–[43]. Mutations on non-catalytic residues that participate into concerted (correlated) movements have been shown to disrupt protein function in NMR relaxation experiments [44]–[46]. Several cases have been documented of mutations that don' t affect the global fold of the protein, but affect protein dynamics and disrupt enzyme function [47]. To overcome this limitation of coarse-grained NMA methods while maintaining the advantages of simplified elastic network models, we developed a new mixed coarse-grained NMA model called Elastic Network Contact Model (ENCoM). ENCoM employs a potential function based on the four bodies potential of STeM with an addition to take in consideration the nature and the orientation of side chains. Side-chain atomic contacts are used to modulate the long range interaction term with a factor based on the surface area in contact [48] and the type of each atom in contact. Additionally, we introduce a non-specific version of ENCoM (ENCoMns) where all interactions between atom types are the same. ENCoM and ENCoMns were validated trough comparison to ANM, GNM and STeM with respect to the prediction of crystallographic b-factors and conformational changes, two properties conventionally used to test ENM methods. Moreover, we test the ability of ENCoM and ENCoMns to predict the effect of mutations with respect to protein stability and compare the ability of ENCoM and ENCoMns to a large number of existing methods specifically designed for the prediction of the effect of single point mutations on protein stability. Finally, we use ENCoM to predict the effect of mutations on protein function in the absence of any effects on protein stability.
We utilized a dataset of 113 non-redundant high-resolution crystal structures [49] to predict b-factors using the calculated ENCoM eigenvectors and eigenvalues as described previously [35] (Equation 4). We compared the predicted b-factors using ENCoM, ENCoMns, ANM, STeM and GNM to the experimental Cα b-factors for the above dataset (Supplementary Table S1). For each protein we calculate the Pearson correlation between experimental and predicted values. The results in Figure 1 represent the bootstrapping average of 10000 iterations. We observe that while comparable, ENCoM (median = 0. 54) and ENCoMns (median = 0. 56) have lower median values than STeM (median = 0. 60) and GNM (median = 0. 59) but similar or higher than ANM (median = 0. 54). It should be noted that it is possible to find specific parameter sets that maximize b-factor correlations beyond the values obtained with STeM and GNM (see methods). However we observe a trade-off between the prediction of b-factors on one side and overlap and the effect of mutations on the other (see methods). Ultimately we opted for a parameter set that maximizes overlap and the prediction of mutations with complete disregard to b-factor predictions. Nonetheless, as shown below, even the lower correlations obtained with ENCoM are sufficiently high to detect functionally relevant local variations in b-factors as a result of mutations. As GNM does not provide information on the direction of movements or the effect of mutations, it is not considered further in the present study. By definition, any conformation of a protein can be described as a linear combination of amplitudes associated to the eigenvectors representing normal modes. It should be stressed that such conformations are as precise as the choice of structure representation used and correct within the quadratic approximation of the potential around equilibrium. Those limitations notwithstanding, one application of NMA is to explore the conformational space of macromolecules using such linear combinations of amplitudes. Pairs of distinct protein conformations, often obtained by X-ray crystallography are used to assess the extent to which the eigenvectors calculated from a starting conformation could generate movements that could lead to conformational changes in the direction of a target conformation. Rather than an optimization to determine the amplitudes for a linear combination of eigenvectors, this is often simplified to the analysis of the overlap (Equation 5), i. e. , the determination of the single largest contribution from a single eigenvector towards the target conformation. In a sense the overlap represents a lower bound on the ability to predict conformational changes without requiring the use of an optimization process. The analysis of overlap for ANM, STeM, ENCoM and ENCoMns was performed using the Protein Structural Change Database (PSCDB) [50], which contains 839 pairs of protein structures undergoing conformational change upon ligand binding. The authors classify those changes into seven types: coupled domain motions (59 entries), independent domain motions (70 entries), coupled local motions (125 entries), independent local motions (135 entries), burying ligand motions (104 entries), no significant motion (311 entries) and other type of motions (35 entries). The independent movements are movements that don' t affect the binding pocket, while dependent movements are necessary to accommodate ligands in the pose found in the bound (holo) form. Burying movements are associated with a significant change of the solvent accessible surfaces of the ligand, but with small structural changes (backbone RMSD variation lower than 1 Å). Despite differentiating between types of movements based on the ligands, the ligands were not used as part of the normal mode analysis. Since side-chain movements associated to the burying movements cannot be predicted with coarse-grained NMA methods, we restrict the analysis to domain and loop movements [51] as these involve backbone movements amenable to analysis using coarse grained NMA methods. For practical purposes, in order to simplify the calculations in this large-scale analysis, NMR structures were not considered. It is worth stressing however that all NMA methods presented here don' t have any restriction with respect to the structure determination method and can also be used with modeled structures. A total of 736 conformational changes, half representing apo to holo changes and the other half holo to apo (in total 368 entries from PSCDB) are used in this study (Supplementary Table S2). Overlap calculations were performed from the unbound (apo) form to the bound form (holo) and from the bound form to the unbound form. Bootstrapped results based onto the best overlap found within the first 10 slowest modes [52], [53] for the different types of conformational changes, domain or loop are shown in figures 2 and 3 respectively. In each case a set of box-plots represent the performance of the four methods being compared, namely STeM, ANM, ENCoM and ENCoMns. The left-most set of box-plots represents the average over all data while subsequent sets represent distinct subsets of the dataset as labeled. The first observation (comparing Figures 2 and 3) is that all tested NMA models show higher average overlaps for domain movements (Figure 2) than loop movements (Figure 3). This confirms earlier observations that NMA methods capture essential cooperative global (delocalized) movements associated with domain movements [51]. Loop movements on the other hand are likely to come about from a more fine tuned combination of normal mode amplitudes than what can be adequately described with a single eigenvector as measured by the overlap. The second observation is that STeM performs quite poorly compared to other methods irrespective of the type of movement (domain or loop). This is somewhat surprising when one compares with ENCoM or ENCoMns considering how similar the potentials are. This suggests that the modulation of interactions by the surface area in contact (the βij terms in Equation 1) of the corresponding side-chains as well as the specific parameters used are crucial. Focusing for a moment on domain movements (Figure 2), ENCoM/ENCoMns outperform all other methods for domain movements in general as well as for every sub category of types of motions therein. Independent movements show lower overlaps than coupled ones, a fraction of those movements may not be biologically relevant due to crystal packing. Interestingly, while there are no differences between the overlap for independent movements starting from the apo or holo forms, this is not the case for coupled movements. In this case (right-most two sets in Figure 2), it is easier to use the apo (unbound) form to predict the holo (bound) form, suggesting that the lower packing in the apo form (as this are frequently more open) generates eigenvectors that favor a more comprehensive exploration of conformational space. Lastly, with respect to loop movements (Figure 3), while it is more difficult to obtain good overlaps irrespective of the method or type of structure used, overall ENCoM/ENCoMns again outperforms ANM. Some of the same patterns observed for domain movements are repeated here. For example, the higher overlap for coupled apo versus holo movements. We observe that ENCoMns consistently performs almost as well as ENCoM irrespective of the type of motion used (all sets in Figures 2 and 3). As side-chain conformations in crystal structures tent to minimize unfavourable interactions, the modulation of interactions by atom types that differentiate ENCoM from ENCoMns plays as minor but still positive role. Normal mode resonance frequencies (eigenvalues) are related to vibrational entropy [54], [55] (see methods). Therefore, it is reasonable to assume that the information contained in the eigenvectors can be used to infer differences in protein stability between two structures differing by a mutation under the assumption that the mutation does not drastically affect the equilibrium structure. A mutation may affect stability due to an increase in the entropy of the folded state by lowering its resonance frequencies, thus making more microstates accessible around the equilibrium. We utilize experimental data from the ProTherm database [56] on the thermodynamic effect of mutations to validate the use of ENCoM to predict protein stability. Here we benefit from the manual curation efforts previously performed to generate a non-redundant subset of ProTherm comprising 303 mutations used for the validation of the PoPMuSiC-2. 0 [57]. The dataset contains 45 stabilizing mutations (ΔΔG<−0. 5 kcal/mol), 84 neutral mutations (ΔΔG [−0. 5,0. 5] kcal/mol) and 174 destabilizing mutations (ΔΔG>0. 5 kcal/mol) (Supplementary Table S2). Each protein in the dataset have at least one structure in the PDB database [58]. As we calculate the eigenvectors in the mutated form we require model structures of the mutants. We generate such models using Modeller [59] and are thus assume that the mutation does not drastically affect the structure. Mutations were generated using the mutated. py script from the standard Modeller software distribution. Modeller utilizes a two-pass minimization. The first one optimizes only the mutated residue, with the rest of the protein fixed. The second pass optimizes the non-mutated neighboring atoms. It is important to stress that our goal is to model the mutated protein as accurately as possible and thus using any method that unrealistically holds the backbone fixed to model the mutant form would be an unnecessary simplification. We observe a backbone RMSD for the whole protein of 0. 01+/−0. 01 Å on average. Considering that RMSD is a global measure that could mask more drastic local backbone rearrangements, we also calculated the average maximum Cα displacement but with a value of 0. 13+/−0. 12 Å we are confident that while not fixed, backbone rearrangements are indeed minimal. In the present work we predict the effect of mutations (Equation 6) for ENCoM, ENCoMns, ANM and STeM and compare the results to existing methods for the prediction of the effect of mutations using the PoPMuSiC-2. 0 dataset above. We compare our results to those reported by Dehouck et al. [57] for different existing techniques: CUPSAT, a Boltzman mean-force potential [60]; DMutant, an atom-based distance potential [61]; PoPMuSiC-2. 0, a neural network based method [57]; Eris, a force field based approach [62]; I-Mutant 2. 0, a support vector machine method [63]; and AUTO-MUTE, a knowledge based four body potential [64]. We used the same dataset to generate the data for FoldX 3. 0, an empirical full atom force-field [65] and Rosetta [66], based on the knowledge based Rosetta energy function. A negative control model was build with a randomized reshuffling of the experimental data. Figure 4 presents RMSE results for each model. The raw data for the 303 mutations is available in Supplementary Table S3. ANM and STeM are as good as the random model when considering all types of mutations together (Figure 4). This is not surprising as the potentials used ANM and STeM are exclusively geometry-based and are thus agnostic to sequence. ENCoMns, ENCoM, Eris, CUPSAT, DMutant, I-Mutant 2. 0 and give similar results and predict significantly better than the random model. AUTO-MUTE, FoldX 3. 0, Rosetta and in particular PoPMuSiC-2. 0 outperform all of the other models. The RMSE values for the subset of 174 destabilizing mutations (Figure 5) shows similar trends as the whole dataset with the exceptions of DMutant losing performance and PoPMuSiC-2. 0 as well as AUTO-MUTE gaining performance compared to the others. It is important to stress that the low RMSE of PoPMuSiC-2. 0 on the overall dataset is to a great extent due to its ability to predict destabilizing mutations. The subset of 45 stabilizing mutations (Figure 6) gives completely different results as those obtained for destabilizing mutations. AUTO-MUTE, Rosetta, FoldX 3. 0 and PoPMuSiC-2. 0 that outperformed all of the models on the whole dataset or the destabilizing mutations dataset cannot predict better than the random model. This is also true for CUPSAT, I-Mutant 2. 0 and Eris. ENCoM and DMutant are the only models with significantly better than random RMSE values for the prediction of stabilizing mutations. ANM and STeM outperform all models on the neutral mutations (Figure 7). All other models fail to predict neutral mutations any better than random. While the accuracy of ANM and STeM to predict neutral mutations may seem surprising at first, it is in fact an artifact of the methodology. As the wild type or mutated structures are assumed to maintain the same general backbone structure, the eigenvectors/eigenvalues calculated with ANM or STeM will always be extremely similar for wild type and mutant forms. Any differences will arise as a result of small variations in backbone conformation produced by Modeller. As such, ANM and STeM predict almost every mutation as neutral, explaining their high success in this case. At first glance, the comparison of ENCoMns and ENCoM could suggest that a large part of the effect observed come from a consideration of the total area in contact and not the specific types of amino acids in contact. However, the side chains in contact are already in conformations that minimize unfavourable contacts to the extent that is acceptable in reality (in the experimental structure) or as a result of the energy minimization performed by Modeller for the mutant form given the local environments. The fact that ENCoM is able to improve on ENCoMns is the actual surprising result and points to the existence of frustration in molecular interactions [67]. Considering that none of the existing models can reasonably predict neutral mutations, the only models that achieve a certain balance in predicting both destabilizing as well as stabilizing mutations better than random and with low bias are ENCoM and DMutant. The analysis of the performance of ENCoM in the prediction of different types of mutations in terms of amino acid properties shows that mutations from small (ANDCGPSV) to big (others) residues are the most accurately predicted followed by mutations between non-polar or aromatic residues (ACGILMFPWV). ENCoM performs poorly on exposed residues (defined as having more than 30% of the surface area exposed to solvent) (Figure 8). It may in principle be possible to find particular linear combinations of ENCoM and other methods that further improve predictions given the widely different (and potentially complementary) nature of the various approaches with respect to ENCoM. We performed linear regressions to find parameters involving ENCoM and each of the other methods in turn that maximize the RMSE difference between the combined models and random predictions (Eq. 8). When considering all types of mutations together, all mixed models perform better than either model individually (left-most column in Figure 8). By definition, a mixed model cannot perform worst than the better of the two models individually. The contribution of ENCoM to the improved performance of the combined model varies according to the model. The ratio of the relative contributions (in parenthesis), broadly classifies the methods into three categories: 1. Methods where ENCoM contributes highly, including I-Mutant (1. 18±0. 25), DMutant (1. 07±0. 31), CUPSAT (1. 02±0. 09) and Eris (1. 02±0. 11); 2. Methods where the contribution of ENCoM is smaller than that of the other method but still significant, including Rosetta (0. 90±0. 11) and FoldX3 (0. 89±0. 08); and finally methods where the addition of ENCoM have a small beneficial effect, including in this class Automute (0. 69±0. 13) and especially PoPMuSIC (0. 18±0. 10). The left-hand side dendrogram in Figure 8 clusters the methods according to their overall accuracy relative to random based on the entire profile of predictions (Eq. 8) for different subsets of the data (columns) according to the type of mutations being predicted. This clustering of methods shows that the relative position of the methods is maintained throughout except for a small rearrangement due to changes in the predictions for CUPSAT. This result suggests that the contribution from the combination of ENCoM to other methods is uniform irrespective of the type of mutations studied. One basic requirement for a system that predicts the effect of mutations on stability is that it should be self-consistent, both unbiased and with small error with respect to the prediction of the forward or back mutations as reported by Thiltgen et al. [68]. The authors built a non-redundant set of 65 pairs of PDB structures containing single mutations (called form A and form B) and utilized different models to predict the effect of each mutation going from the form A to form B and back. From a thermodynamic point of view, the predicted variation in free energy variation should be of the same magnitude for the forward or back mutations, ΔΔGA→B = −ΔΔGB→A. Using the Thiltgen dataset we performed a similar analysis for ENCoM, ENCoMns, ANM, STeM, CUPSAT, DMutant, PoPMuSiC-2. 0 and a random model (Gaussian prediction with unitary standard deviation). For the remaining methods (Rosetta, Eris and I-Mutant) we utilize the data provided by Thiltgen. We removed three cases involving prolines as such cases produce backbone alterations. Furthermore, PoPMuSiC-2. 0 failed to return results for five cases. The final dataset therefore contains 57 pairs (Supplementary Table S4). The CUPSAT and AUTO-MUTE servers failed to predict 25 and 32 cases respectively. As these failure rates are significant considering the size of the dataset, we prefer to not include these two methods in figures 9 and 10 (the remaining cases appear however in Supplementary Table S4). The results in figure 9 show that compared to the random model (a positive control in this case), Rosetta and FoldX 3. 0 show moderate bias while PoPMuSiC-2. 0 and I-Mutant show significant bias. All biased methods are biased toward the prediction of destabilizing mutations (data not shown) in agreement with the results in Figure 3. DMutant, Eris, ENCoM and ENCoMns are the only models with bias comparable to that of the random model (the positive control in this experiment). ENCoM, ENCoMns, and to a lesser extent Rosetta and DMutant have lower errors than the random model (Figure 10). All other methods display an error equal or higher than that of the random model. ENCoM and ENCoMns vastly outperform all the others models in terms of error. Lastly, STeM and ANM show low and moderate biases respectively and errors equivalent to random (data not shown) but as mentioned, these methods cannot be used for the prediction of mutations (other than neutral mutations as an artefact). Mutations may not only affect protein stability but also protein function. While experimental data is less abundant, one protein in particular, dihydrofolate reductase (DHFR) from E. coli, has been widely used experimentally to understand this relationship [69], [70]. Recently, Boher et al. [47] have analyzed the effect of the G121V mutation on protein dynamics in DHFR by NMR spectroscopy. This mutation is located 15 Å away from the binding site but reduces enzyme catalysis by 200 fold with negligible effect on protein stability (0. 70 kcal/mol). The authors evaluated, among many other parameters, the S2 parameter of the folic acid bound form for the wild type and mutated forms and identify the regions where the mutation affects flexibility. We calculated b-factor differences (Equation 4) between the folate-bound wild type (PDB ID 1RX7) and the G121V mutant (modeled with Modeller) forms of DHFR (Supplementary Table S5). We obtain a good agreement (Pearson correlation = 0. 61) between our predicted b-factor difference and S2 differences (Figure 11). As mentioned earlier, the overall correlation of 0. 54 in the prediction of b-factors (Figure 1) appears at least in this case to be sufficient to capture essential functional information.
Our results show that a small modification of the long-range interaction term in the potential energy function of STeM had an important positive impact on the model. This small change improves the method in comparison to existing NMA methods in the traditional areas such as the prediction of b-factors and conformational sampling (overlap) where coarse-grained normal mode analysis are applied. More importantly however, it opens an entire new area of application to coarse-grained normal mode analysis methods. Specifically ENCoM is the first coarse-grained normal-mode analysis method that permits to take in consideration the specific sequence of the protein in addition to the geometry. This is introduced through a modification in the long-range interactions to account for types of atoms in contact modulated by their surface in contact. As a validation of the approach we explored the ability of the method to predict the effect of mutations in protein stability. In doing so we created the first entropy-based methodology to predict the effect of mutations on the thermodynamic stability of proteins. This methodology is entirely orthogonal to existing methods that are either machine learning or enthalpy based. Not only the approach is novel but also the method performs extremely favourably compared to other methods when viewed in terms of both error and bias. As the approach taken in ENCoM is completely different from existing methods for the prediction of the effect of mutations on protein stability, a new opportunity arises to combine ENCoM with enthalpy and machine-learning methods. Unfortunately, we tried to create a naïve method based on linear combinations of the predictions of ENCoM and the different methods presented without success, perhaps due to the large bias characteristic to the different methods. To assess the relative importance of contact area and the modulation of interactions with atom types, we tested a model that has non-specific atom-type interactions (ENCoMns), this model is atom type insensitive, but is sensitive the orientation of side-chain atoms. While a large fraction of the observed effect can be attributed to surfaces in contact only, ENCoM is consistently better than ENCoMns, particularly at predicting destabilizing mutations where the possibility to accommodate unfavourable interactions is more restricted. We cannot however exclude the effect of the intrinsic difficulty in modeling destabilizing mutations. For stabilizing mutations, the near equivalence of ENCoM and ENCoMns may be explained in part by the successful energy minimization of the mutated side-chain performed by Modeller. ANM and STeM failed to predict the effect of mutations on the whole dataset. They were not expected to perform well because their respective potentials only take in account the position of alpha carbons (backbone geometry). As such ANM and STeM tend to predict mutations as neutral, explaining their excellent performance onto the neutral subset and failure otherwise. Our results suggest that surfaces in contact are essential in a coarse-grained NMA model to predict the effect of mutation and that the specific interactions between atom types is necessary to get more subtle results, particularly stabilizing mutations. ENCoM is consistently better than ENCoMns in the prediction of loop or domain movements irrespective of the dependency of the coupling of this movement to ligand binding or the starting structure (apo or holo form) and both outperform ANM and STeM. Our results corroborate previous work on a mix coarse-grained method adding a atomistic resolution to loops capable of improving the prediction of loop movements [40]. ENCoM performs considerably better than STeM throughout despite having very similar potentials, showing the importance of surfaces in contact in the prediction of movements. There is little difference between ENCoMns and ENCoM in the prediction of b-factors, but both perform worst than ANM, STeM and GNM. At least in the case of DHFR b-factor differences capture some essential characteristics of the system as calculated by NMR. However, one should be careful in placing too much emphasis on the validation of b-factor predictions using experimental data derived from crystals as these are affected to a great extent by rigid body motions within the crystal [71]. PoPMuSiC-2. 0, AUTO-MUTE, FoldX 3. 0 and Rosetta perform better than other models in the whole test dataset of mutations. However, the dataset consists of 15% stabilizing mutation, 57% of destabilizing and 28% of neutral mutations. When looking at each subset, machine learning or enthalpy based models failed to predict better than random on the stabilizing mutations subset. Biases in the dataset may have affected the training of machine-learning methods. For example the training set of PoPMuSiC-2. 0 contains 2648 mutations in proportions that are similar to those in the testing set with 60%, 29% and 11% destabilizing, neutral and stabilizing mutations respectively. While it is true that most mutations tend to be destabilizing, if one is interested in detecting stabilizing mutations, a method over trained on destabilizing mutations will not meet expectations. Indeed, PoPMuSiC-2. 0 and I-mutant the two machine learning based methods, have larger biases and errors than other methods in their predictions. Our method relies on a model structure of the mutant. As the modeling may fail to find the most stable side-chain conformation, it could have a bias toward giving slightly higher energies to the mutant. Notwithstanding this potential bias, ENCoM have the lowest error and bias. This may be a case where less is more as the coarse-grained nature of the method makes it also less sensitive to errors in modeling that may affect enthalpy-based methods to a greater extent. Finally, there is one more advantage in the approach taken in ENCoM. As the network model is a global connected model it considers indirectly the entire protein, while in existing enthalpy or machine-learning methods the effect of a mutation is calculated mostly from a local point of view. The prediction of the thermodynamic effect of mutations is very important to understand disease-causing mutations as well as in protein engineering. With respect to human diseases, and particularly speaking of cancer mutations, one of the factors that may lead to tumour suppressor or oncogenic mutations is their effect on stability (the authors thank Gaddy Getz from the Broad Institute for first introducing us to this hypothesis). Specifically, destabilizing mutations in tumour suppressor genes or alternatively stabilizing mutations in oncogenes may be driver mutations in cancer. Therefore the prediction of stabilizing mutations may be very important to predict driver mutations in oncogenes. Likewise, in protein engineering, one major goal is that of improving protein stability with the prediction of stabilizing mutations. Such mutations may be useful not as the final goal (for purification or industrial purposes) but also to create a ‘stability buffer’ that permits the introduction of potentially destabilizing additional mutations that may be relevant to create the intended new function. The work presented here is to our knowledge also the most extensive test of existing methods for the prediction of the effect of mutations in protein stability. The majority of methods tested in the present work fail to predict stabilizing mutations. However, we are aware that the random reshuffled model used may be too stringent given the excessive number of destabilizing mutations in the dataset. The only models that predict stabilizing as well as destabilizing mutations are ENCoM and DMutant, however ENCoM is the only method with low self-consistency bias and error. While the contribution of side chain entropy to stability is well established [72], [73], here we use backbone normal modes to predict stability. As a consequence of the relationship between normal modes and entropy, our results attest to the importance of backbone entropy to stability and increase our understanding of the overall importance of entropy to stability. The strong trend observed on the behaviour of different parameters sets with respect to the α4 parameter is very interesting. Lower values are associated with better predictions of conformational changes while higher values are associated with better b-factor predictions. One way to rationalize this observation is to consider that higher α4 values lead to a rigidification of the structure, adding constraints and restricting overall motion. Likewise, lower α4 values remove constraints and thus lead to higher overlap. We used ENCoM to predict the functional effects of the G121V mutant of the E. coli DHFR compared to NMR data. This position is part of a correlated network of residues that play a role in enzyme catalysis but with little effect on stability. The mutation affects this network by disrupting the movement of residues that are far from the binding site. We can predict the local changes in S2 with ENCoM. As these predictions are based on b-factor calculations, this result shows that at least in this case, even with b-factor prediction correlation lower than ANM, STeM and GNM we can detect functionally relevant variations. Clearly, despite the greater performance of GNM, ANM or STeM in the calculation of b-factors, these methods cannot predict b-factor differences as a consequence of mutations, as their predictions are the same for the two forms. While a more extensive study is necessary involving S2 NMR parameters, our results serve as an example against relying too heavily on crystallographic b-factors for the evaluation of normal mode analysis methods.
In order to obtain a set of parameters to be used with ENCoM we performed a sparse exhaustive integer search of the logarithm of parameters with for to maximize the prediction ability of the algorithm in terms of overlap and prediction of mutations. In other words, we searched all combinations of 13 distinct relative orders of magnitude for the set of 4 parameters. For each parameter set, we calculated the bootstrapped median RMSE (see below) Z-score sum for the prediction of stabilizing and destabilizing mutations, . Keeping in mind that lower RMSE values represent better predictions, the 2000 parameter sets (out of 28561 combinations) with highest were then used to calculate Z-scores for overlap in domain and loop movements. As our goal is to obtain a parameter set that combines low RMSE and high overlap, we ranked the 1000 parameter sets according to. The parameter set with highest is (solid black line in Figure 13). The optimization of the bootstrapped median is equivalent to a training procedure with leave-many-out testing. Given the dichotomy in predicting the effect of mutations and overlap on the one hand and b-factors on the other, we provide, the following is the best parameter set observed for the prediction of b-factors with average b-factor correlation of. The exploration of parameter space shows that there is a clear trade-off between the prediction of mutations (low RMSE), conformational sampling (high overlap) and b-factors (high correlations). Parameter sets that improve the prediction of b-factors are invariably associated with poor conformational changes (low overlap) associated to both domain and loop movements and variable RMSE for the prediction of mutations (red lines in Figure 13). On the other hand, parameter sets that predict poorly b-factors, perform better in the prediction of conformational changes and the effect of mutations (blue lines in Figure 13). The parameters used in STeM, are arbitrary, taken without modifications from a previous study focusing on folding [77]. As expected, this set of parameters can be considerably improved upon as can be observed in Figure 13 (dashed line). The four right-most variables in the parallel coordinates plot in Figure 13 show the logarithm of the α parameters for each parameter set. Either class of parameter sets, better for b-factors (in red) or better for overlap/RMSE (in blue) come about from widely diverging values for each parameter across several orders of magnitude. There are however some patterns. Most notably for α4, where there is an almost perfect separation of parameter sets around α4 = 1. Interestingly, higher values of α1 and α2, associated with stronger constraints on distances and angles tend also to be associated to better overlap values. While it is likely that a better-performing set of parameters can be found, the wide variation of values across many orders of magnitude show that within certain limits, the method is robust with respect to the choice of parameters. This result justifies the sparse search employed. Bootstrapping is a simple and general statistical technique to estimate standard errors, p-values, and other quantities associated with finite samples of unknown distributions. In particular, bootstrapping help mitigate the effect of outliers and offers better estimates in small samples. Bootstrapping is a process by which the replicates (here 10000 replicates) of the sample points are stochastically generated (with repetitions) and used to measure statistical quantities. In particular, bootstrapping allows the quantification of error of the mean [78]–[82]. Explained in simple terms, two extreme bootstrapping samples would be one in which the estimations of the real distribution of values is entirely made of replicates of the best case and another entirely of the worst case. Some more realistic combination of cases in fact better describes the real distribution. Thus, bootstrapping, while still affected by any biases present in the sample of cases, helps alleviate them to some extent. One of the most common types of experimental data used to validate normal mode models is the calculation of predicted b-factors and their correlation to experimentally determined b-factors. B-factors measure how much each atom oscillates around its equilibrium position [6]. Predicted b-factors are calculated as previously described [35]. Namely, for a given Cα node (i), one calculates the sum over all eigenvectors representing internal movements (n = 7 to 3N) of the sum of the squared ith component of each eigenvector in the spatial coordinates x, y and z normalized by the corresponding eigenvalues: (4) We calculate the Pearson correlation between predicted and experimental b-factors for each protein and average random samples according to the bootstrapping protocol described above. The overlap is a measure that quantifies the similarity between the direction of movements described by eigenvectors calculated from a starting structure and the variations in coordinates observed between that conformation and a target conformation [52], [53]. In other words, the goal of overlap is to quantify to what extent movements based on particular eigenvectors can describe another conformation. The overlap between the nth mode, , described by the eigenvector is given by (5) where represent the vector of displacements of coordinates between the starting and target conformations. The larger the overlap, the closer one can get to the target conformation from the starting one though the movements defined entirely and by the nth eigenvector. We calculate the best overlap among the first 10 slowest modes representing internal motions. Insofar as the simplified elastic network model captures essential characteristics of the dynamics of proteins around their equilibrium structures, the eigenvalues obtained from the normal mode analysis can be directly used to define entropy differences around equilibrium. Following earlier work [54], [55], the vibrational entropy difference between two conformations in terms of their respective sets of eigenvalues is given by: (6) In the present work the enthalpic contributions to the free energy are completely ignored. Therefore, in the present work we directly compare experimental values of ΔG to predicted ΔS values. In order to use the same nomenclature as the existing published methods, we utilize ΔΔG to calculate the variation of free energy variation as a measure of conferred stability of a mutation. A linear regression going through the origin is build between predicted ΔΔG and experimental ΔΔG values to evaluate the prediction ability of the different models. The use of this type of regression is justified by the fact that a comparison of a protein to itself (in the absence of any mutation) should not have any impact on the energy of the model and the model should always predict an experimental variation of zero. However, a linear regression that is not going through the origin would predict a value different from zero equal to the intercept term. In other words, the effect of two consecutive mutations, going from the wild type to a mutated form back to the wild type form (WT→M→WT) would not end with the expected net null change. The accuracy of the different methods was evaluated using a bootstrapped average root mean square error of a linear regression going through the origin between the predicted and experimental values. We refer to this as RMSE for short and use it to describe the strength of the relationship between experimental and predicted data. If one was to plot the predicted energies variation of ΔΔGA→B versus ΔΔGB→A and trace a line y = −x, the bias would represent a tendency of a model to have points not equally distributed above or below that line while the error would represent how far away a point is from this line. In other words, considering a dataset of forward and back predictions, the error is a measure of how the predicted ΔΔG differ and the bias how skewed the predictions are towards the forward or back predictions [68]. A perfect model, both self-consistent and unbiased, would have all the points in the line. Statistically, the measures of bias and error are positively correlated. The higher the error for a particular method, the higher the chance of bias. We determined the efficacy of linear combinations involving ENCoM and any of the other models for the prediction of the effect of mutations on thermostability as follows. For a given bootstrap sample of the data, we rescaled the predictions of each model as follows: (7) where the vector notation signifies all data points in the particular bootstrap sample and the index represents each model as well as the random model. We then use singular value decomposition to determine the best parameter the normalized predicted values to calculate the parameters that maximize the RMSE difference between the linear combination model and the random predictions as follows (8) whereand (9) The bootstrapped average is then calculated from the 10000 bootstrap iterations. The relative contribution of ENCoM and the model under consideration is given by the ratio of the parameters. It is interesting to note that this ratio could be seen as an effective temperature factor, particularly considering that predicted values are primarily enthalpic in nature for certain methods and entropy based in ENCoM.
Normal mode analysis (NMA) methods can be used to explore potential movements around an equilibrium conformation by mean of calculating the eigenvectors and eigenvalues associated to different normal modes. Each normal mode represents a global collective, correlated and complex, form of motion of the entire protein. Any conformation around equilibrium can be represented as a weighted combination of normal modes. Differences in the magnitudes of the set of eigenvalues between two structures can be used to calculate differences in entropy. We introduce ENCoM the first coarse-grained NMA method to consider atom-specific side-chain interactions and thus account for the effect of mutations on eigenvectors and eigenvalues. ENCoM performs better than existing NMA methods with respect to traditional applications of NMA methods but is the first to predict the effect of mutations on protein stability and function. Comparing ENCoM to a large set of dedicated methods for the prediction of the effect of mutations on protein stability shows that ENCoM performs better than existing methods particularly on stabilizing mutations. ENCoM is the first entropy-based method developed to predict the effect of mutations on protein stability.
biochemistry
proteins
protein structure
biology and life sciences
computational biology
biophysics
molecular biology
macromolecular structure analysis
biophysical simulations
A Coarse-Grained Elastic Network Atom Contact Model and Its Use in the Simulation of Protein Dynamics and the Prediction of the Effect of Mutations