Many tools for error detection are currently available (Marti-Renom et al. 2000), most of which use statistical potentials or molecular mechanics force fields to assess the accuracy of a protein model. They include contact potentials (Miyazawa and Jernigan 1985; DeBolt and Skolnick 1996; Park and Levitt 1996; Park et al. 1997; Melo et al. 2002), residue–residue distance-dependent potentials (Sippl 1993a; Jones 1999; Xia et al. 2000; Melo et al. 2002), residue-based solvent-accessibility potentials (Jones et al. 1992; Sippl 1993b; Melo et al. 2002), atomic solvent accessibility, and pairwise interaction potentials (Melo and Feytmans 1997, 1998; Lazaridis and Karplus 1998; Samudrala and Moult 1998; Lu and Skolnick 2001; Zhou and Zhou 2002; Wang et al. 2004). All these scoring functions have proven to be able to discriminate between native and misfolded proteins.

The accurate detection of structural errors constitutes a challenging problem, and it is clear that different scoring functions exhibit varying degrees of performance (Park and Levitt 1996). The success of a given energy or scoring function at detecting structural errors will depend on the nature of the model being assessed and on the parameters of the energy function itself. In this respect, most of the protein decoys available for the testing of new scoring functions are not complete and contain many models of low resolution, with incomplete atomic coordinates and significant structural errors (Sanchez and Sali 1998; Jones 1999; Tsai et al. 2003). Therefore, if the ability of scoring functions in the detection of errors of small magnitude is going to be tested, new decoys containing highly accurate and complete protein models are needed.

A scoring function capable of detecting structural errors in highly accurate protein models should have some specific requirements. First, the scoring function must be able to detect a few specific errors of small magnitude that would probably be scattered along the protein sequence and mapping onto different locations in the three-dimensional (3D) structure. Then, this scoring function should be highly sensitive and therefore parameterized at the atomic level with a high-resolution description of the interactions. Second, it is also probable that few localized errors of small magnitude within a structure will not interact with each other in 3D space at a short distance range. Thus, a successful scoring function at detecting errors in highly accurate protein models should not require a strong nonlocal component, nor have a long interacting distance range to describe the atomic interactions.

The requirements described above point toward a central role of a local and short distance range component in a scoring function that is capable of detecting small and localized errors in protein structures. To date, all knowledge-based potentials described have a local term with a lower limit that accounts for the interactions of atoms belonging to two consecutive amino acids. In contrast, force fields like CHARMM (Brooks et al. 1983; MacKerell et al. 1998), AMBER (Weiner and Kollman 1981), and CVFF (Hagler et al. 1974) go beyond this limit and calculate the energy from the 1–4 nonbonded interactions for the Lennard-Jones and the Coulomb terms of the energy function. Therefore, local interactions within a residue are currently considered and their energy calculated when using force fields, but not when using knowledge-based potentials.

Knowledge-based potentials do not incorporate nonbonded terms within a residue because of the lack of an appropriate reference system to accurately calculate them. These energy functions are derived from a set of native protein conformations that have been experimentally solved. The reference state is also calculated from the same data set, because an unbiased and representative experimental data source for the unfolded state of proteins is not yet available. Both the uniform density (Sippl 1990; Samudrala and Moult 1998; Lu and Skolnick 2001) and the distance-scaled finite ideal-gas (Zhou and Zhou 2002) reference states cannot explicitly deal with atom–atom connectivity issues, which become extremely important for nonbonded interactions of short distance range. Additionally, the assumption of independence among interatomic distances for these type of interactions is clearly invalid (Dill 1997). The use of current reference systems to derive the pseudo-energy of nonbonded interactions at a short distance range leads to functions that are either flat (not informative) or nearly identical for different interacting atom pairs (not discriminative).

We have already described the calculation of an atomic knowledge-based potential that considered local and nonlocal interactions (Melo and Feytmans 1997). As in other potentials, the nonbonded interactions within residues were not included, and sequence separation was used to discriminate between local and nonlocal interactions. The shape of each pairwise energy curve of this potential exhibited a high variability for very short sequence separations, which resulted from its local energy component (Melo and Feytmans 1997). On the other hand, when the sequence separation was increased to values larger than ~10 amino acids, a monotonic curve shape was observed for a given pairwise interaction (Melo and Feytmans 1997). Based on those observations, we concluded that the splitting of nonlocal interactions into different functions is unnecessary, and therefore two atoms that belong to sequentially distant amino acids could be considered as free particles. Then, in a subsequent study, we derived a new atomic potential that included only the nonlocal interactions. In this nonlocal potential, all interacting atom pairs belonging to residues with sequence separations above 10 residues were collapsed and represented by a single energy function (Melo and Feytmans 1998). This grouping substantially increased the number of observations and allowed us to increase the number of classes of interatomic distances, resulting in a more accurate potential. Additionally, a short distance range that allows a maximum of about two atomic shells was used, thus minimizing the impact of side effects due to atom connectivity. Since then, the resulting potential, which is called ANOLEA, has been used to assess the quality of protein structure models (Melo and Feytmans 1998; Gonzalez et al. 2004; Feliciangeli et al. 2006; Pizzo et al. 2006), and it is currently incorporated on the SwissModel Web server and repository (Schwede et al. 2003; Kopp and Schwede 2004).

The ANOLEA potential was also incorporated into MODELLER (Sali and Blundell 1993), a comparative protein structure modeling software, and used for the prediction of the conformation of loops through molecular dynamics simulations (Fiser et al. 2000). In that work, the ANOLEA potential (Melo and Feytmans 1998) was combined with the CHARMM-22 molecular mechanics force field (Brooks et al. 1983; MacKerell et al. 1998) by replacing Lennard-Jones and Coulomb nonbonded terms from CHARMM with the complete set of terms in ANOLEA. Therefore, the 1–2 bond, 1–3 angle, and 1–4 dihedral and improper dihedral energies were obtained from CHARMM; and the 1–5 and above nonbonded pseudo-energies were obtained by cubic spline interpolation from the discrete ANOLEA functions. Thus, all 1–5 or higher nonbonded interactions were obtained from a single pseudo-energy function that was derived nonlocally for each atomic pair. Contributions from both potentials, as well as residue side-chain dihedral angle pseudo-energies derived from the observed statistical preferences in experimental data (Sali and Blundell 1993), were equally weighted and combined to get the total pseudo-energy of the system. It turned out that this approach resulted in the most accurate predictions of loop conformations (Fiser et al. 2000).

In this study, we have evaluated the impact of incorporating 1–5 nonbonded terms on the ability to detect small and localized errors in a benchmark of protein structure models built by comparative modeling. To avoid the connectivity issues and the problems with the reference system mentioned above, and motivated by the results obtained in the modeling of loops (Fiser et al. 2000), the distance-dependent 1–5 nonbonded terms of each atom–atom interaction were directly extrapolated from our previous nonlocal knowledge-based potential (Melo and Feytmans 1998).

We begin the Results section by describing the building of a new benchmark set of accurate comparative models, which was used to test the potentials in error detection. We continue by assessing the performance in error detection of different potentials with and without the calculation of nonbonded terms. After that, we compare the performance of the potential that includes nonbonded terms described here against the performance of other existing and commonly used energy functions. We conclude by summarizing and discussing the main lessons learned from this study. In the Materials and Methods section, we describe in detail the building of the benchmark set of models and define the derivation and calculation scheme for the different knowledge-based potentials, as well as the criteria used to assess their performance.

Figure 1. Flowchart of the procedure used to build the sets of comparative protein models. The source protein model set was originally constituted of 3.375 models with the correct fold, which were built in a previous study using MODPIPE (Sanchez and Sali 1998). (more ...)

The first set, the class A set, contains 55 highly accurate protein models with >95% of their Cα atoms within 3.5 Å from their corresponding native structures, after the optimal superposition of the structures in 3D space. The total Cα root mean square deviation (RMSD) of all residues in each model is <1.1 Å. The second set, the class B set, contains 57 accurate protein models with >90% of their Cα atoms within 3.5 Å from their corresponding native structures. The total Cα RMSD of all residues in each model is >1.5 Å and <2.6 Å.

The structural modeling accuracy of each residue in a model was assessed by comparison to the corresponding native structure. The 55 models in the class A set contain a total of 10,295 residues, out of which only 201 contain structural errors. The 57 models in the class B set contain 10,714 residues, out of which a total of 1257 were erroneously modeled. For details about the building of these two sets of models and the definition of structural errors, see Materials and Methods.

Table 1. Variants of knowledge-based potentials assessed in this work

To evaluate the effect of the type of interactions considered on the ability to detect small and localized errors in protein structure models, ROC curves were built for each potential based on the energy profiles obtained with each protein model in the class A and class B sets (Fig. 2). The results clearly show that the consideration of those interactions among pairs of atoms with a short connectivity range is important to accurately detect the wrongly modeled residues, as it is demonstrated by the performance achieved with the NL-NB potential that include the calculation of all 1–5 and above nonbonded terms. The NL-NB potential is the best classifier of residue accuracy in the tested sets, because it has a higher sensitivity and specificity than the other potentials for all possible classification thresholds. The observed performance of the NL-NB potential is higher for the class A set, which contains more accurate models than those in the class B set. Interestingly, the NB-NB potential exhibits a worse performance than the NL-NB and L-L (local-local) potentials, in both sets of models. The performance of the NB-NB potential is similar to that observed for the All-All potential, although to the fact that the former potential includes a significantly larger number of terms for the evaluation of the protein models. The observed effect is similar in both sets of models.

Figure 2. Global performance on error detection of local, nonlocal, and nonbonded potentials. The different potentials calculated in this study are assessed globally by means of ROC curve analysis. (A) Results obtained for the highly accurate models or class A (more ...)

The second best classifier is the L-L potential, but the overall difference in performance from the NL-NB potential is very small for the class B set. These two potentials are equivalent in performance when a classifier of high specificity is required for the class B set. The NL-NL potential exhibits a similar performance to that of the NB-NB and All-All potentials for the class A set, and it is the worst classifier of residue modeling accuracy for the class B set.

The comparative performance of optimal classifiers based on these potentials was also assessed in the class A and class B sets of models (Table 2). The results show that for the class A set, the NL-NB potential is clearly the best, achieving an accuracy of 80%, a specificity of 80%, and a sensitivity of 70%. The second best classifier, the L-L potential, achieves an accuracy of 72%, a specificity of 73%, and a sensitivity of 67%. The other potentials perform clearly worse on this set upon the evaluation of these measures. For the class B set, although a similar trend is observed, the NL-NB potential is not clearly better than the L-L potential.

Table 2. Optimal performance on error detection of local, nonlocal, and nonbonded potentials

The results of this comparison show that the NL-NB potential clearly outperforms the other potentials in both sets of models, for almost all possible classification thresholds (Fig. 3). The most important observation is the large difference in performance for the class A set when a high specificity is required. For example, for a fixed specificity of 95% (i.e., a false-positive rate of 5%), the NL-NB potential exhibits a sensitivity of 45%. For the same specificity value, DFIRE-C_α and DFIRE-C_β have a sensitivity of ~18%, and ProSa exhibits a sensitivity of ~11%. Analogously, for a fixed specificity of 80%, the NL-NB potential exhibits a sensitivity of 71%. This performance compares favorably against that obtained for DFIRE-C_α, DFIRE-C_β, and ProSa, which exhibit sensitivities in the range between 30% and 40%. For the class B set, the same trend is observed, although the differences are not as large as those observed for the class A set.

Figure 3. Global performance on error detection of other energy functions. The global performance of DFIRE and ProSa potentials and the semiempirical CHARMM force field are compared, by means of ROC curve analysis, against that obtained for the NL-NB potential (more ...)

As a negative control, the semiempirical potential CHARMM was also assessed. This potential was used in molecular dynamics simulations followed by energy minimization when the models of the two sets were built. Therefore, it was initially expected that the ability of this force field to detect errors in these models would be limited. The results show that the CHARMM force field exhibits a very poor performance, similar to that of a random classifier, at identifying those residues that were wrongly modeled (Fig. 3).

The comparative performance of optimal classifiers based on the NL-NB and the external potentials was also assessed in the class A and class B sets of models (Table 3). In accordance with the observations described above, the NL-NB potential is clearly the best classifier of residue modeling accuracy in both sets of models.

Table 3. Optimal performance on error detection of other energy functions

In this study, we have assessed the impact of using close nonbonded interactions from a knowledge-based potential on the detection of small and localized errors in protein structure models. Since data sets of accurate protein structure models were not available, we have first built two sets containing accurate and highly accurate protein models. Although the residues in these protein models were arbitrarily categorized as either correctly modeled or in error, it is important to mention that at least this definition was adopted by considering the visual inspection of the models after the optimal superposition with their corresponding native structures. We believe that this is an important feature that should be incorporated in other data sets of models. Most of the existing data sets have been automatically generated, and the models are often arbitrarily classified upon some strict, but not necessarily well established, criteria. We acknowledge, however, that this task is not always feasible or straightforward to achieve. The two data sets of models generated in this study have been publicly released through our Web site (http://protein.bio.puc.cl/sup-mat.html) thus they are freely available for everyone willing to test the performance of different energy or scoring functions in the detection of small and localized errors in protein structures.

The central point addressed in this study involves the incorporation of close nonbonded terms for the evaluation of protein structure conformations. By close nonbonded terms, we refer to those nonbonded interactions occurring between pair of atoms that are closely connected through a few intermediate covalent bonds (i.e., 1–5 interactions and above). This type of interaction is clearly not independent of the protein backbone and side-chain configurations of the amino acids since some pairs of atoms are restricted to be close in 3D space because of the existing connectivity. Additionally, some interactions such as those occurring between atoms that belong to a ring are very restricted and do not adopt alternative conformations. This poses a challenge to current methodologies used to derive the knowledge-based potentials, if these terms are to be included. The problem resides in the reference system used to derive the potential. The derivation of knowledge-based potentials for protein structure prediction is carried out by using the inverse Boltzmann law, which in its general form states:

where δE(s) represents the change of energy associated with the transition between the unfolded and the folded states defined by the variable s; p^F(s) represents the probability of occurrence of the subsystem defined by s in the folded state F; p^U(s) represents the probability of occurrence of the subsystem defined by s in the unfolded state U; R is the gas constant; and T is the absolute temperature measured in kelvins. Unfortunately, the term p^U(s), which describes the reference state, cannot be directly calculated because a homogeneous and unbiased experimental sample of the unfolded state of proteins is not available. As an attempt to circumvent this problem, the observed interactions among any pair of atoms by assuming a null interaction model are often used as a reference system to derive the potential. However, these probability estimations are also derived from the same experimental data set as the p^F(s) term, this is the folded state F, and it could not necessarily correspond to the real probabilities that exist in the unfolded state U.

When a reference system based on the folded state is used to derive a knowledge-based potential, most of the close nonbonded terms result in noninformative or nondiscriminate energy functions. A possibility to circumvent this problem, as it has been proposed in this work, is to derive a nonlocal potential with a short distance range and then to extrapolate the resulting energy values of nonlocal atom pairs to the close nonbonded terms. The calculation of a nonlocal potential with a short distance range is a good approximation to capture atom–atom interactions as if they were free particles, because there should not exist any connectivity constraints for two nonlocal atoms to be close in space. In addition, the short distance range also allows us to minimize some side effects or dependency on bonded atoms when deriving the potential. For example, the existence of a real interaction between two atoms i and j should not lead to the observation of an artificial favorable interaction between atoms bonded to i and atoms bonded to j.

Based on a benchmark of accurate protein models, we have demonstrated that this approximation exhibited the best performance at classifying residue modeling accuracy. The incorporation of close nonbonded terms by direct extrapolation improved the potential's performance as compared to other potentials that either included or not more distant local terms. On the other hand, the direct derivation of a potential including close nonbonded terms resulted in a poor classifier of residue modeling accuracy. These results suggest that close nonbonded terms are important to describe native protein properties and should be included in the potentials, but their direct derivation by using the current methodologies and experimental data sets does not lead to an accurate description of this type of interaction.

The accurate detection of small and localized errors in the benchmark of protein structures used in this study is not an easy task. This is demonstrated by the maximum accuracy, specificity, and sensitivity achieved by all the potentials tested here. The best performance was achieved by the NL-NB potential, with a maximum accuracy of ~80%, a sensitivity of 70%, and a specificity of 80%. This means a rate of 20% false positives and a rate of 30% false negatives, which certainly does not constitute an ideal classifier. However, it is important to mention that this performance was achieved by considering only distance-dependent pairwise terms. It is expected that the incorporation of solvation terms will further improve the current performance.

Knowledge-based potentials are informatics functions (Solis and Rackovsky 2006). Their capacity to properly describe the atomic interactions that are recurrent in native protein conformations depends on many parameters and on how the data are compressed and classified. In addition to this, the performance of a potential will not only depend on how the information was extracted, compressed, and classified, but also on how the information is used. In this study, we have demonstrated that a potential that was derived by considering only the nonlocal interactions can be successfully used to infer the close nonbonded interactions by direct extrapolation. We are not claiming that these extrapolated terms would reflect the true existing nonbonded interactions, but only that this methodology allows to improve the performance for the task of classifying residue modeling accuracy. This improvement could only be explained by the fact that there is a given amount of information gain for this classification task when using this methodology. This information gain was certainly not achieved when using the NB-NB potential, although the total number of terms considered to estimate the accuracy of a given residue was exactly the same as that for the NL-NB potential.

The limited performance of the ProSa and DFIRE potentials in the benchmark of models used in this study is partially due to the fact that these are residue-based potentials, and, as such, they only include a description of the C_α, C_β, or backbone atom interactions. Therefore, a fair comparison of their performance in error detection against full-atom potentials is simply not possible. However, on the other hand, these potentials are distributed in this form, and therefore it is the only way that they can be used and compared against other potentials in a particular benchmark and error detection task. It is evident that, for the detection of errors of small magnitude like the ones present in our benchmark of models, residue-based potentials exhibit a low sensitivity and specificity.

Future directions to further improve the performance of knowledge-based potentials for the task of error detection include the incorporation of implicit solvation models to account for the interactions between protein and solvent atoms, the consideration of atom shielding to describe more accurately the effective atomic interactions, and the exploration of different reference systems to derive the potentials. We are presently working on these topics and expect to incorporate these elements in a new generation of knowledge-based potentials that should be more accurate for the task of protein structure assessment.

The distance-dependent potentials were calculated as described (Sippl 1993a; Melo and Feytmans 1997, 1998):

where M_ijk is the number of occurrences for the interaction of atom types i and j at sequence separation k:

r is the number of classes of distance; and σ is the weight given to each observation. σ = 0.02 was used (Sippl 1990), so that with 50 observations f_k^ij(d) and f_k^xx(d) have equal weights for the calculation of δE_k^ij(d) · f_k^ij(d) is the relative frequency of occurrence for the interaction of atom types i and j at sequence separation k in the class of distance d:

f_k^xx(d) is the relative frequency of occurrence for all the interactions of any two atom types at sequence separation k in the class of distance d:

where n is the number of different atom types; and r is the number of distance classes. The temperature T was set to 293 K, so that RT is equivalent to 0.582 kcal/mol.

The common parameters of the potentials tested in this work are the same as those previously optimized for the ANOLEA potential (Melo and Feytmans 1997). The potentials were derived for a maximal distance range of 7.0 Å using 35 discrete classes of 0.2 Å each. A total of 40 atom types were defined for all heavy atoms in the 20 standard amino acids.

The detection of residues with structural errors was carried out by means of smoothed and normalized energy profiles. The normalized total energy score per residue (E^R) is defined as follows:

where N is the total number of i atoms belonging to a given residue, and M is the total number of atoms j that interact with atom i at sequence separation k and below the distance range defined in the potential. E_k^ij(d) corresponds to the energy score assigned by the statistical potential to the interaction between atoms i and j at sequence separation k and at distance d. T^R is the total sum of i–j interactions recorded for residue R and lies in the range 0 ≤ T^R ≤ N × M. The normalized energy profiles were smoothed by a sliding window with a length of seven residues. These normalized and smoothed energy profiles were finally used to assess the structural errors in our benchmark of protein models. Thus, for each protein model, the three amino acids at the C and N termini were not assessed. We initially calculated smoothing windows of three, five, seven, and nine residues. Because the windows of five, seven, and nine residues gave similar results, but slightly better than those obtained when using a window of three residues, we decided to fix this parameter to the middle value of seven. Different potentials showed the same trend upon changing the smoothing window size parameter value.

The models in the class A set represent 35 distinct folds. In the class B set, 34 distinct folds are represented. Therefore, even though not all the models in these two sets were built for proteins representing different folds, they are not strongly biased to any particular fold since a large fraction of the models has a distinct fold. In terms of a more general classification, as the one defined by the composition and arrangement of secondary structure elements, a fairly good and poorly biased representation is achieved. For the class A set, 26% of the models contain only α-helices as secondary structure elements, 24% contain only β-sheets, 18% contain α and β, and 33% of the models contain α + β secondary structure elements in their structures. In the class B set of models, 36% of the models contain only α-helices, 31% contain only β-sheets, 13% contain α and β, and 20% contain α + β.

The 57 accurate and 55 highly accurate protein models were selected from a existing set of 3375 models with a correct fold that has been described previously (Sanchez and Sali 1998; Melo et al. 2002). This original set of 3375 models was built by the comparative modeling of representative chains of the Protein Data Bank (PDB) (Berman et al. 2002). The models were built based on the correct templates and mostly correct alignments between the target sequences and the template structures. The models were obtained by applying MODPIPE to 1.085 representative chains of the PDB (Sanchez and Sali 1998). These representative sequences corresponded to the protein chains in PDB that shared <30% sequence identity or were >30 residues different in length. The templates for comparative modeling, selected by MODPIPE, were 1637 PDB chains with <80% sequence identity to each other or >30 residues different in length. Each target sequence was aligned separately with each one of the 1637 known structures using the program ALIGN, which implements local sequence alignment by dynamic programming (Altschul 1998). Only the target–template alignments with a significance score higher than 22 nats (corresponding approximately to the PSI-BLAST e-value of 10⁻⁴) were used, resulting in 3993 models. Models with <30% structural overlap with the actual experimentally determined structure were eliminated. Structural overlap was defined as the fraction of the equivalent C_α atoms upon least squares superposition of the two structures with the 3.5 Å cutoff. This procedure also removed models based on correct templates that had a poor alignment and models based on templates that had large domain or rigid-body movements with respect to the target structure. The final set contained 3375 correct models (Melo et al. 2002).

The set of 3375 correct models was initially filtered by updating the target and template structures currently available at the PDB and checking the sequence alignments originally used to build the models. A total of 132 models presented inconsistencies between the target sequence in the original alignment used to build the model and the current target sequence available at the PDB. These 132 models were removed, thus resulting in a total of 3243 entries (Fig. 1). Then a second filter was applied, and we selected only those protein models of a length >100 residues for which >90% of their residues were possible to model. Finally, as explained above, two independent filters were applied to select those models belonging to the class A and class B sets (Fig. 1). All models in both sets were built for target monomeric proteins. The 3D coordinates of these models in PDB format are available as supplemental material at http://protein.bio.puc.cl/sup-mat.html.

where tp is the true-positive rate, TP is the count of true-positive instances, P is the total number of positive instances (i.e., the sum of TP and FN counts), fp is the false-positive rate, FP is the count of false-positive instances, and N is the total number of negative instances (i.e., the sum of TN and FP counts).

Based on these two rates, receiver operating characteristic (ROC) curves were calculated for each potential as previously described (Fawcett 2004) and used to assess its performance. Receiver operating characteristic (ROC) curves are two-dimensional graphs in which the TP rate is plotted on the Y-axis and the FP rate is plotted on the X-axis (Swets 1988; Swets et al. 2000; Fawcett 2004). An ROC graph depicts relative tradeoffs between benefits (true positives) and costs (false positives) for all possible decision thresholds.

In addition to the ROC curves, which constitute the best way to compare two classifiers, we also calculated four overall measures to assess and compare the performance of the potentials in the detection of errors in accurate and highly accurate protein structure models. The first measure was the area under the ROC curve (AUC), which ranges from 0.5 to 1.0. The AUC represents somehow an overall measure or summary of the ROC itself, and it has an important statistical property: The AUC of a classifier is equivalent to the probability that the classifier will rank a randomly chosen positive instance higher than a randomly chosen negative instance (Fawcett 2004). The second measure was the accuracy (ACC), which is defined as:

Finally we define the sensitivity (Sn) and specificity (Sp) measures as follows:

Since almost all the metrics mentioned above depend on the particular decision threshold that is used to classify the instances (with AUC as the only exception), we report these metrics at a single and fixed value that is called the optimal threshold (OT). The OT is uniquely defined by the point in the ROC curve that has the minimal distance to the upper-left corner of the ROC graph, which would correspond to a perfect classifier (i.e., fp = 0.0 and tp = 1.0).

This work was funded by grant 1051112 from FONDECYT.