A pairwise residue contact area-based mean force potential for discrimination of native protein structure

Considering energy function to detect a correct protein fold from incorrect ones is very important for protein structure prediction and protein folding. Mainly, two different types of potential energy function are currently in use, either on the identification of native protein models from a large set of decoys or protein fold recognition and threading studies [1–10]. The first class of potentials, named physical-based potential, is based on the fundamental analysis of the forces between the particles referred to as physical energy function. The second type is knowledge-based energy function based on information from known protein structures. In physical energy function, a molecular mechanics force field is used. Molecular mechanics force fields are parameterized from ab-initio calculation and small molecule structural data. They are essentially the sum of pairwise electrostatic and Van der Waals interaction energies, bonds, angles and dihedral angle terms [11–14]. In addition, terms that are not included such as entropy and the solvent effect are implicitly considered. Although, physical energy function is widely used in molecular dynamic simulation of proteins in their native and denatured states which can be used to distinguish the decoy/native structures, but these functions have not been efficient in protein structure prediction because of their greater computational cost. To reduce the computational complexity of the protein folding problem, knowledge-based or empirical mean-force potential is widely used. Since the structure of folded proteins reflects the free energy of the interaction of all their components, including all enthalpic and entropic contributions, as well as solvent effects, such potentials provide an excellent shortcut towards a powerful objective function. It can be used to force the system to obtain potential between groups of atoms by the use of experimentally determined structures. In this approach, statistical thermodynamics is used in an analysis of the frequency of observed states to estimate the underlying free energy [15]. Most often, the distribution of pairwise distances are used to extract a set of effective potential between residues or atoms. The distribution of pairwise distances can be compiled from the protein structure database and by defining a reference state, Boltzmann equation is used to calculate the interaction energy of a particular pair. The total potential energy of a protein is simply taken as a sum over all pairwise interactions. In most cases, one or two points for each residue are used to represent a protein [16–18]. These points are usually C(alpha), C(beta) or the centre of mass of each side chain. Each interaction can be distance-dependent. A large variety of knowledge-based potential of mean force have been developed by introducing additional interactions such as surface area terms, the main chain and side chain dihedral angles, three and four body terms and heavy atoms [6, 19–23].

In the contact potential, either distance-dependent or contact based, the distance between the centres of two C(alpha), C(beta) or centre of mass of two residues or the all heavy atoms of two residues are calculated and the observed frequency of contacts between residues converts to free energy using Boltzmann equation. In this way, two problems may be encountered. First, when an atom or centre of mass is selected for each residue, calculated potential is independent of orientation of the side chains and when the distance between two atoms of two residues are equal to the distance of two atoms of other residues in other positions, the same potentials are assigned to them although the orientation of two residue side chains may be quite different. Second, the atoms of two residues may not have direct contact with each other and some atoms may be located in an interval close to them.

In this study, we develop a new approach to calculate a knowledge-based potential energy using pairwise residue contact area. We calculate the parts of each pairwise residue area that are in contact in Å2 by rolling a probe ball of different sizes around the atoms of a residue to determine the contacts area of each pair. This pairwise contact area is used to determine statistical contact area preference between each residue pairs, when a contact area preference estimates a sum of energetic interactions and structural constraints.

A good energy function at its minimum should discriminate native structures from decoys. So, to test the effectiveness of this new potential, we calculate it on several decoy sets to measure its ability to discriminate native structure from decoys. Several decoy sets that contain one to hundreds of decoy structures generated in different ways are used and in the most cases this potential has been able to distinguish native structures from the decoys. Calculated Z-score and Pe, which are useful measures of the validity of the computed potential, show high value for all protein datasets.

One of the best ways to show the performance of a force field is its ability to find the native structure in a large set of decoys. Different decoys sets have been used to evaluate how well knowledge-based potentials and physical potentials discriminate native structures. In this study, the performance of our model based on pairwise contact area was tested on models from different decoy sets containing misfold, DecoyForMMPBSA, fisa, hgstructal, semfold, vhp-mcmd, 4state_reduced, lmds, ig_structal, ig_structal_hires, HRDecoy and Rosetta_Tsai. The quality of the models in decoy sets and the members of decoy structures are very different.

From the principle of statistical mechanics, we suppose that the energy of the native structure has to be minimum energy among all conformations and much lower than the average energy of all possible conformations. Then, in addition to finding the fold with lower energy, the Z-Score has been calculated. Energy profiles have been made for each residue pairwise separated by d residues in sequence (10 distinct values for sequence separation have been considered) using different probe sizes (r = 0.25, 0.5, 0.75, 1, 1.5, 2, 2.5 Å). So, we have 70 different energy profiles for energy calculation. The total energy of particular structure of a protein can be calculated as the sum of all the pairwise interactions, but the best discrimination has been achieved when only energy profile for d = 1 has been considered. This shows that contact area of consequent neighbors has more important role in distinguishing native folds from the decoys. In this situation, with r = 0.25Å we could discriminate native folds from decoys in almost all decoy sets. The increase in the probe size has had slightly improvement on the discrimination power shown by Z-Score. Since the increase of the probe radius resulted in an increase in the amount of calculations, so choosing the great probe radius was not efficient.

Tables in Additional file 1 show the details of results for proteins in each decoy set. Although these decoy sets have been produced in different ways andfinding native structures in a set of low quality models would not be difficult, but our pairwise contact area based potential ranked the native structures first in the most cases. There are some exceptions. First, in semfold, protein 1nkl has the rank five. Second in hg_structal, 1gdm has the rank 17, but in the last version of PDB this protein has been replaced by 2gdm and it is notable that when energy has been calculated for 2gdm, this protein ranked first in decoy structures. In the Rosetta_Tsai dataset, our model could not find any native structures in the first rank.

The results obtained on different decoy sets show good performance of our methodology to discriminate native folds. However, an experiment to evaluate the performance of an energy model when performing ab-inito folding is to discriminate between the native-like and non-native structures. In Additional file 1, contact area energy is plotted against the RMSD from native structure for native and all decoy structures. Different datasets have different distributions of RMSD for non-native proteins. Usually RMSD's are calculated using C _α distances and residue side chain atoms are not considered. Since our method is based on content area of side chain atoms, then it is very sensitive to orientation of side chains atoms although the change in the position of C _α may not be very large. As shown in plots in Additional file 1, in the most cases decoy structures have RMSD more than 2 Å and in these cases the contact area of side chain atoms may be far from native structure and there is a cosiderable distance between the energy of native and decoy structures in most datasets.

The aim of this study was to evaluate a mean force potential based on contact area of residues instead of contact or distance to separate correct from incorrect folds. This was done by calculation of contact area of all atoms of residues considering the Van der Walls spheres of atoms and obtaining a coefficient from training dataset used to quantify pairwise potential in a protein fold. The new potential not only is residue orientation-dependent, but also gives residue contact area in angstrom square for each pair of atoms in adjacent residues that provides a better quantification of atomic interaction than distance-based methods.

The analysis in this work showed that the best definition is the one involving the contact area between Van der Walls spheres of atoms of any two consecutive residues with employing a cut off distance around 0.5 Å. Considering atomic radii, those distances around 5 Å has been considered. This, in fact, is close to the cut off distance considered in the contact-based potential methods. Contact area-based potential was able to recognize the native structures on different decoy sets with a high degree of accuracy, as evident from the Z-score. Only one of 1386 native fold in 11 decoy sets was not ranked first, however this protein, 1nlk in semfold decoy set, was ranked five among 11600 models.

These results show that in addition to the important role of contact area between two atoms in improvement of potential function that reflects the orientation of residue side chain, short range contact between neighbor residues play more important role than long range contact.