BioMed Research International
Volume 2017 (2017), Article ID 5760612, 17 pages
https://doi.org/10.1155/2017/5760612

Research Article

All-Atom Four-Body Knowledge-Based Statistical Potentials to Distinguish Native Protein Structures from Nonnative Folds

Majid Masso

School of Systems Biology, George Mason University, 10900 University Blvd. MS 5B3, Manassas, VA 20110, USA

Correspondence should be addressed to Majid Masso; mmasso@gmu.edu

Received 27 June 2017; Revised 13 August 2017; Accepted 23 August 2017; Published 8 October 2017

Academic Editor: Rita Casadio

Copyright © 2017 Majid Masso. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Recent advances in understanding protein folding have benefitted from coarse-grained representations of protein structures. Empirical energy functions derived from these techniques occasionally succeed in distinguishing native structures from their corresponding ensembles of nonnative folds or decoys which display varying degrees of structural dissimilarity to the native proteins. Here we utilized atomic coordinates of single protein chains, comprising a large diverse training set, to develop and evaluate twelve all-atom four-body statistical potentials obtained by exploring alternative values for a pair of inherent parameters. Delaunay tessellation was performed on the atomic coordinates of each protein to objectively identify all quadruplets of interacting atoms, and atomic potentials were generated via statistical analysis of the data and implementation of the inverted Boltzmann principle. Our potentials were evaluated using benchmarking datasets from Decoys-‘R’-Us, and comparisons were made with twelve other physics- and knowledge-based potentials. Ranking 3rd, our best potential tied CHARMM19 and surpassed AMBER force field potentials. We illustrate how a generalized version of our potential can be used to empirically calculate binding energies for target-ligand complexes, using HIV-1 protease-inhibitor complexes for a practical application. The combined results suggest an accurate and efficient atomic four-body statistical potential for protein structure prediction and assessment.

1. Introduction

Over recent years, exponential growth of the Protein Data Bank (PDB) [1] has facilitated the selection of larger, nonredundant subsets of experimentally solved protein structures at higher resolutions, which in turn have provided the data used in developing more effective knowledge-based statistical potentials for improved structure prediction. In contrast to physics-based energy functions, statistical potentials generally perform better and are more computationally efficient at identifying the native structure as a global minimum [2, 3]. Distance-dependent statistical potentials often focus on pairwise atomic contacts within macromolecular structures [4, 5]; however, such energy functions fail to take into consideration important higher-order contributions based on multibody interactions [6, 7]. Indeed, use of an “atomic environment potential” for which neighborhood sizes vary by atom previously demonstrated improved performance at discriminating between native and near-native protein structures [3]. In the present work we employed the well-established computational geometry tiling technique of Delaunay tessellation [8], for objectively identifying all quadruplets of nearest neighbor atoms in order to develop, evaluate, and apply all-atom four-body statistical potentials for protein structure prediction.

Four-body statistical potentials were derived based on PDB atomic coordinate file data corresponding to single chains selected from over 1400 diverse protein structures. Delaunay tessellation was applied to the three-dimensional (3D) atomic coordinates of each protein chain, whereby atoms were treated as vertices to generate a convex hull encompassing thousands of space-filling, nonoverlapping, irregular tetrahedra (Figure 1). For assurance that each tetrahedron identifies at its four vertices a quadruplet of atoms that are pairwise all within a prescribed distance from one another, a subsequent edge-length cutoff parameter may be introduced; removal of a tetrahedral edge between a pair of atoms longer than this cutoff eliminates from the tessellation all tetrahedra sharing that edge. Depending on the size of the atomic alphabet used for labeling points, the four atoms appearing at vertices of any particular tetrahedron in these tessellations represent one of 35 ( letters), 330 (), or 8855 () possible distinct atomic quadruplet types. For each cutoff (if any) and alphabet size, statistical data obtained from the protein chains and their tessellations included the following: () observed relative frequencies of interaction for each type of atomic quadruplet, based on their rates of occurrence as tetrahedral vertices; and () rates expected by chance for each atomic quadruplet type, based on relative frequencies of individual atom types in the protein chains and use of a multinomial reference distribution. Through application of the inverted Boltzmann principle [9, 10], the negative logarithm of the ratio of observed to expected rates of occurrence was used to calculate an empirical energy of interaction for each atomic quadruplet, which collectively form an atomic four-body statistical potential.

Figure 1: HIV-1 protease (a) ribbon and (b) atomic ball-and-stick diagrams. The atomic coordinates are used as tetrahedral vertices to generate (c) the Delaunay tessellation of the protein chain, a convex hull consisting of thousands of space-filling and nonoverlapping tetrahedra, each of whose vertices objectively identifies a quadruplet of nearest neighbor atoms. The modified tessellation in (d) is obtained by removing all edges longer than 12 Å between pairs of atoms, thereby eliminating all tetrahedra that share those edges and excluding their corresponding atomic quadruplets from consideration as nearest neighbors.

The approach implemented here at the atomic level was motivated by its prior successful application at the residue level [11–16]. All atomic 3D coordinates in proteins are considered in this work to generate the all-atom four-body statistical potentials, while previously developed residue-based four-body potentials used only a single point per amino acid (e.g., or residue center of mass). Clearly, there is a degree of information loss with the coarser-grained residue representation of proteins relative to the finer all-atom representation. Both approaches implement the Delaunay tessellation algorithm, which uses the respective point-sets to serve as vertices for generating a tetrahedral tiling of the protein structure that objectively identifies quadruplets of nearest neighbors (i.e., either residues or atoms). Given its significantly more sparse point-set, a residue-based (i.e., one point per residue) tessellation typically yields a few hundred tetrahedra, whereas tessellation applied to all the atoms in the same protein structure has on the order of a few thousand tetrahedra.

Upon selecting an atomic alphabet and edge-length cutoff as parameters, the energy of any folded protein chain would subsequently be calculated with the atomic four-body potential as follows: label and tessellate the 3D atomic coordinates of the structure according to the same parameters, refer to the previously derived atomic four-body potential under those parameters to assign a score to each tetrahedron in the tessellation equal to the interaction energy of the atomic quadruplet found at its four vertices, and add up the scores of all the tetrahedra in the tessellation. The all-atom four-body statistical potentials that we developed were each evaluated by scoring multiple decoy directories in the Decoys-‘R’-Us benchmarking database [17]. We compared these four-body potentials to one another, based on standard performance metrics, as well as to the knowledge-based potentials of Fogolari et al. [18] and Summa et al. [3]; the latter study detailed performance results for 10 diverse physics- and knowledge-based potentials to conduct their own comparisons, hence providing us an opportunity to assess our four-body potentials relative to a dozen other methods in total. Lastly, we report on a practical application, related to predicting target-inhibitor binding energy, by implementing a modification of our best performing four-body potential.

2. Methods

2.1. Protein Training Set

A nonredundant set of 1417 high-resolution (≤2.2 Å) crystallographic structures, with atomic coordinate files deposited in the PDB, were culled using the PISCES server [28] with the constraint that the single protein chains selected from the structures shared low (<30%) sequence identity (http://binf2.gmu.edu/automute/tessellatable1417.txt). The ensemble of structure files is diverse, consisting of single- and multichain proteins, the vast majority of which are additionally complexed to small molecular or peptide ligands. Coordinates of hydrogen atoms and water molecules were removed from all files prior to proceeding with the analyses.

2.2. Designation of Atoms

For each of the 1417 protein chains, three alphabets were explored for defining atom types and labeling points corresponding to their 3D atomic coordinates. In the first instance, a simple four-letter alphabet (C, N, O, and S) accounts for all atoms and ensures sufficient frequency data are collected for all possible atomic quadruplets observed at the four vertices of tetrahedra in Delaunay tessellations (Table 1). Clearly, the same atom type may appear at more than one of the four vertices of any tetrahedron in a protein tessellation, and given that those vertices are unordered, all permutations of the four atomic letters at the vertices of a tetrahedron refer to the same quadruplet, so that an alphabetical ordering (e.g., COON) of the atoms can be used as a singular representation. In this case, a combinatorial argument [29] shows that the number of distinct subsets of size letters that can be formed from an atomic alphabet of size is given by

Hence, letters admit distinct atomic quadruplets. Next, an atomic alphabet consisting of letters (amino acid backbone: , , , and ; side-chain: , , , and ) differentiates between backbone alpha- and carbonyl-carbon atoms, distinguishes residue backbone atoms from those in side-chains, and can form distinct atomic quadruplets. Lastly, we explored a maximum diversity of quadruplet atomic interactions with letters as described in Summa et al. [3], which groups atoms based on common traits, including bonding pattern, partial charge, and hydrophobicity, and generates distinct atomic quadruplets.

Table 1: Summary data for the protein structure training set (1417 single chains).

Atomic alphabet
Count
Proportion

C
1572222
0.634149

N
425874
0.171774

O
469869
0.189520

S
11299
0.004557

Total atom count
2479264

Total tetrahedron counts

    No edge-length cutoff
16152638

    12 Å edge-length cutoff
15497203

    8 Å edge-length cutoff
14567713

    4.8 Å edge-length cutoff
9569503

Same counts regardless of atomic alphabet size.

2.3. Derivation of the Atomic Four-Body Statistical Potentials

Delaunay tessellations for the 1417 single protein chains were generated by submitting their respective atomic coordinates as input to the Qhull program [30], visualizations of the tessellated structures were obtained by utilizing the output data from Qhull to create plots within Matlab, and molecular graphics were produced with Chimera [31] (Figure 1). An in-house suite of Perl programs was used for all data formatting and analyses related to the tessellated structures (Table 1). In particular, for each atomic alphabet of size the relative frequencies of occurrence for all types of atomic quadruplets (, , , ) were calculated as the proportion of tetrahedra among all the tessellations for which the four atoms appear on the vertices. Four separate sets of relative frequencies were calculated for each of the three atomic alphabets explored, based on the original protein tessellations (no edge-length cutoff applied), as well as tessellations modified by introducing cutoffs of length 12 Å, 8 Å, and 4.8 Å. The use of an 8 Å cutoff is consistent with that used by other researchers to generate atomic pair potentials [32], while the other two cutoffs were also selected to identify the appropriate choice for an atomic four-body potential.

For each of the three atomic alphabets, we additionally computed relative frequencies of occurrence , for the atom types in all 1417 single protein chains. These frequencies, in turn, were needed for calculating the rate expected by chance for all types of atomic quadruplets (, , , ) obtained using a multinomial reference distribution, given by

In the formula above, represents the number of occurrences of atom type in the quadruplet. For each pair of parameters selected (i.e., alphabet size and cutoff), we applied the inverted Boltzmann principle to calculate a score for quantifying the interaction energy for all types of atomic quadruplets (, , , ), as described by Sippl [9, 10], thus defining a particular atomic four-body statistical potential function (Table 2 and Figure 2). A total of 12 four-body potentials were generated and evaluated in this study (3 atomic alphabets × 4 edge-length cutoffs, which includes the case where no cutoff is applied to the original tessellations).

Table 2: Atomic four-body statistical potentials employing a 4-letter alphabet.

Quad
No cutoff
12 Å cutoff
8 Å cutoff
4.8 Å cutoff

Count

Count

Count

Count

CCCC
1711740
0.183570
1692278
0.170562
1643281
0.156444
922742
0.224573

CCCN
1823746
0.190871
1790819
0.180796
1726843
0.169730
1229054
0.134910

CCCO
2807489
0.046213
2750767
0.037094
2616669
0.031930
1533444
0.081509

CCCS
119435
−0.201538
117868
−0.213776
114416
−0.227745
53175
−0.077468

CCNN
832442
0.140340
803373
0.137776
764343
0.132553
612799
0.046021

CCNO
3838549
−0.179748
3746112
−0.187158
3593186
−0.195913
2695736
−0.253612

CCNS
53655
0.055872
52228
0.049592
50119
0.040628
28817
0.098480

CCOO
1643096
−0.069578
1559797
−0.064976
1408171
−0.047423
796121
0.017752

CCOS
86638
−0.109530
84188
−0.115055
79331
−0.116117
36077
0.043594

CCSS
6408
−0.898511
6343
−0.912057
6183
−0.927840
3685
−0.885580

CNNN
64504
0.507783
51089
0.591028
39927
0.671251
18568
0.821250

CNNO
961282
−0.145664
903031
−0.136526
842880
−0.133431
647785
−0.201598

CNNS
7628
0.335839
6916
0.360394
6293
0.374540
3499
0.446952

CNOO
1380693
−0.260215
1283449
−0.246502
1168615
−0.232641
757532
−0.226873

CNOS
44097
−0.082434
41962
−0.078877
39453
−0.078957
22706
−0.021519

CNSS
2153
−0.691035
2131
−0.704562
2085
−0.721950
1312
−0.703279

COOO
336824
−0.081946
278989
−0.018132
207773
0.083015
65646
0.400894

COOS
17883
0.051201
16090
0.079093
13492
0.128713
4711
0.403174

COSS
2068
−0.630846
2009
−0.636259
1897
−0.638215
1001
−0.543084

CSSS
214
−1.741768
214
−1.759742
207
−1.772176
125
−1.735618

NNNN
4632
0.482308
2068
0.814499
1118
1.054783
133
1.796871

NNNO
36223
0.233848
22370
0.425144
15881
0.547103
5948
0.791107

NNNS
407
0.564301
263
0.735926
203
0.821548
58
1.183114

NNOO
190771
−0.268893
158110
−0.205359
137906
−0.172816
97822
−0.206171

NNOS
3088
0.204035
2530
0.272582
2171
0.312201
1061
0.440643

NNSS
236
−0.599166
230
−0.605982
224
−0.621354
79
−0.351235

NOOO
129494
−0.234026
92243
−0.104725
67495
0.004100
26079
0.234579

NOOS
5426
0.001929
4626
0.053197
4007
0.088737
1832
0.246129

NOSS
436
−0.522015
418
−0.521701
405
−0.534836
220
−0.452305

NSSS
138
−2.118466
138
−2.136453
137
−2.160159
48
−1.887182

OOOO
39551
−0.278298
23332
−0.067103
12170
0.188717
1542
0.903421

OOOS
1462
0.137035
1007
0.280952
636
0.453674
78
1.182534

OOSS
158
−0.339520
143
−0.314191
125
−0.282625
32
0.126633

OSSS
22
−1.278314
22
−1.296297
21
−1.302962
5
−0.862215

SSSS
50
−3.855858
50
−3.873832
50
−3.900710
31
−3.875603

(a)

(b)

Figure 2: Graphical representations for two four-body potentials, based on an eight-letter alphabet with a 12 Å edge-length cutoff, and a twenty-letter alphabet with a 4.8 Å edge-length cutoff. Here = alpha-carbon, = backbone carbonyl-carbon, and = side-chain sulfur (from either cysteine or methionine) represent the same atom types in both alphabets, with quadruplets and appearing at the same extremes of both potentials. Despite millions of tetrahedra generated by the 1417 protein tessellations irrespective of the cutoff length (see Table 1), note that 3 of 330 atomic quadruplet types (, and ) did not appear at all as tetrahedral vertices based on an 8-letter atomic alphabet with a 12 Å cutoff ( = side-chain nitrogen atom), while 1935 of 8855 quadruplets types were not observed under a 20-letter alphabet with a 4.8 Å cutoff.

In any atomic tessellation of a protein structure, two adjacent tetrahedra may share a common vertex (1 atom), a common edge (2 atoms), or a common triangular face (3 atoms). Although the two adjacent tetrahedra represent two sets of atomic quadruplets that may share up to 3 atoms in common, those quadruplets are distinct by virtue of the atom(s) that the two tetrahedra do not share. The collective interaction of a quadruplet of atoms as a fundamental unit in this four-body scenario is analogous to the interaction of two atoms in the development of a pair potential, whereby a given atom may be considered to interact with each of several neighboring atoms by virtue of satisfying a prescribed distance cutoff between itself and each of the neighbors, and therefore the atom is shared by all of those pairs; likewise, two atomic quadruplets from adjacent tetrahedra in a tessellation may share up to 3 atoms and yet remain fundamentally distinct quadruplets. Moreover, since Delaunay tessellation does not distinguish types of bonds and generates a tetrahedral tiling by objectively identifying quadruplets of nearest neighbor atoms based solely on their six collective pairwise distances from each other, all covalent bonds as well as noncovalent interactions between particular pairs of atoms are included together in these tetrahedral atomic quadruplets without the need to explicitly identify and segregate them. Recent studies suggest that covalent interactions are informative when combined with nonbonded interactions [33, 34].

2.4. Decoy Database

A significant collection of models provided in the Decoys-‘R’-Us database (http://compbio.buffalo.edu/dd/) form a well-established and challenging standard for benchmarking the performance of energy functions. Several categories are located under the heading “The multiple decoy sets,” each containing a number of decoy model directories. Each such directory is named after the PDB accession code of the native crystallographic protein structure and contains coordinate files for that native structure as well as for numerous decoy model structures (i.e., alternative conformations for a given native structure); additionally, the directory includes a file that provides the root mean square deviations (rmsds) for all the alternative models relative to the native structure. For this work, we focused on the following decoy set categories: 4_state_reduced, fisa, fisa_casp3, hg_structal, ig_structal, ig_structal_hires, lattice_ssfit, and lmds.

3. Results

3.1. Energy Calculations and Benchmark Evaluation Measures

Energy calculations were made for 145 native protein structures as well as for all of their respective decoy models downloaded from the 8 decoy set categories in the Decoys-‘R’-Us database. To this end, all native and decoy structures were tessellated, and their energies were repeatedly computed using all twelve four-body potentials under their respective parameters of atomic alphabet size and tessellation edge-length cutoff. Given the energy scores for a native structure and its collection of decoys, all calculated using the same four-body potential, the following measures of performance were evaluated.

(1) Native Rank. Among the native protein and all decoys, the structures are ranked in ascending order according to increasing energy (i.e., lowest energy structure has rank 1).

(2) Z-Score. This measurement is defined as

where is energy of the native structure, is mean energy over all decoy models, and is standard deviation of the distribution of decoy energies [3]. A large positive -score indicates a wide gap between the energy of the native protein and the mean decoy energy.

(3) Correlation Coefficient (r). It is the linear correlation between calculated energy and rmsd. For decoys with low rmsds relative to the native structure, good correlation is preferable; however, this is unlikely if decoys are significantly misfolded with high rmsds.

(4) Fractional Enrichment (FE). It is the proportion of decoy structures corresponding to the lowest 10% of rmsds that are also found among those corresponding to the lowest 10% of calculated energy scores.

The raw performance data obtained with the four-body potential derived using a 4-letter atomic alphabet and a 12 Å tessellation edge-length cutoff as parameters (Table 2) are presented in Table 3. As such, we employed the same 4-letter alphabet to label the atomic vertices in the tessellations of the 145 native structures and all of their respective decoys, and edges longer than 12 Å were removed from all tessellations prior to calculating total energies as described in the last paragraph of the Introduction. Data analogous to that of Table 3 were obtained using each of the 11 other four-body potentials generated for this study under alternative parameter-pair values for atomic alphabet size and tessellation edge-length cutoff (raw data not shown).

Table 3: Performance evaluation on 145 benchmarks in 8 decoy sets from Decoys-‘R’-Us, based on energies obtained with the four-body potential derived using a 4-letter alphabet and 12 Å cutoff as parameters.

Decoy set
PDB ID
Native rank
-score

FE

4state_reduced
1ctf
156/631
0.7
0.27
11.1

1r69
6/676
2.5
0.18
14.9

1sn3
23/660
1.8
0.42
48.5

2cro
5/674
2.4
0.40
20.9

3icb
137/654
0.8
0.34
21.5

4pti
1/687
3.5
0.48
45.6

4rxn
83/677
1.1
0.31
19.4

fisa
1fc2
497/501
−2.2
−0.22
4.0

1hddC
442/501
−1.2
−0.02
8.0

2cro
32/501
1.5
0.07
10.0

4icb
465/500
−1.4
−0.02
6.0

fisa_casp3
1bg8A
10/1201
2.4
0.06
14.2

1bl0
883/972
−1.4
−0.19
5.2

1eh2
1624/2414
−0.5
0.09
16.2

1jwe
539/1408
0.3
0.04
9.3

smd3
1/1201
2.8
0.07
10.0

hg_structal
1ash
2/30
2.2
0.47
33.3

1babB
1/30
3.7
0.17
66.7

1colA
1/30
3.9
0.46
33.3

1cpcA
4/30
1.4
0.19
0.0

1ecd
11/30
0.3
0.00
0.0

1emy
1/30
3.3
0.38
66.7

1flp
1/30
2.3
0.21
66.7

1gdm
3/30
1.5
0.12
33.3

1hbg
1/30
2.8
0.03
33.3

1hbhA
1/30
3.4
0.18
33.3

1hbhB
3/30
1.7
−0.27
33.3

1hdaA
1/30
5.2
0.27
33.3

1hdaB
1/30
4.2
0.42
33.3

1hlb
1/30
2.0
0.25
33.3

1hlm
5/30
0.9
0.12
0.0

1hsy
12/30
0.3
0.49
33.3

1ithA
1/30
3.4
0.17
66.7

1lht
13/30
0.1
0.12
0.0

1mba
1/30
5.7
0.31
33.3

1mbs
30/30
−2.8
0.16
33.3

1mygA
3/30
1.3
0.25
66.7

1myjA
14/30
0.1
0.29
0.0

1myt
20/30
−0.8
0.06
33.3

2dhbA
30/30
−2.0
−0.23
33.3

2dhbB
24/30
−0.9
0.45
33.3

2lhb
1/30
4.7
0.49
66.7

2pghA
6/30
1.1
0.24
33.3

2pghB
1/30
2.7
0.29
33.3

4sdhA
1/30
2.7
−0.01
33.3

lattice_ssfit
1beo
1/1999
4.9
0.04
7.0

1ctf
88/2001
1.7
−0.04
11.0

1dktA
13/1999
2.5
0.00
8.5

1fca
3/2001
3.1
0.03
16.0

1nkl
92/1998
1.7
0.01
17.1

1pgb
506/2000
0.7
−0.04
13.0

1trlA
68/2000
1.8
0.03
10.5

4icb
223/2000
1.2
−0.03
13.0

ig_structal_hires
1dvf
1/20
2.1
0.13
50.0

1fgv
1/20
3.9
0.20
50.0

1flr
16/20
−0.3
−0.08
0.0

1fvc
1/20
3.5
0.56
50.0

1gaf
1/20
3.8
−0.01
50.0

1hil
1/20
2.6
0.55
50.0

1ind
7/20
0.2
0.36
0.0

1kem
1/20
2.4
0.46
50.0

1mfa
1/20
2.0
0.10
50.0

1mlb
9/20
0.2
−0.51
0.0

1nbv
11/20
−0.2
0.20
0.0

1opg
20/20
−3.8
−0.45
0.0

1vfa
1/20
2.2
0.39
50.0

1vge
2/20
1.9
0.42
50.0

2cgr
1/20
2.6
0.39
50.0

2fb4
7/20
0.3
−0.09
0.0

2fbj
19/20
−1.1
−0.12
0.0

6fab
3/20
1.5
0.25
0.0

7fab
10/20
0.2
0.27
0.0

8fab
1/20
2.6
0.14
50.0

ig_structal
1acy
56/61
−1.5
−0.01
0.0

1baf
3/61
2.1
0.12
16.7

1bbd
61/61
−3.1
−0.15
0.0

1bbj
23/61
0.1
−0.02
0.0

1dbb
24/61
0.1
−0.31
0.0

1dfb
8/61
1.1
0.20
0.0

1dvf
1/61
2.7
0.11
16.7

1eap
1/61
2.6
0.16
16.7

1fai
2/61
1.9
0.29
33.3

1fbi
12/61
0.8
−0.06
33.3

1fgv
1/61
2.9
0.03
16.7

1fig
2/61
2.0
0.11
16.7

1flr
40/61
−0.2
−0.08
16.7

1for
1/61
3.5
0.16
33.3

1fpt
1/61
3.7
0.12
16.7

1frg
4/61
1.5
0.39
16.7

1fvc
1/61
2.8
0.16
33.3

1fvd
9/61
1.1
0.09
33.3

1gaf
1/61
3.5
0.04
16.7

1ggi
9/61
1.1
−0.03
0.0

1gig
10/61
1.0
0.16
0.0

1hil
2/61
2.5
0.41
16.7

1hkl
1/61
4.3
0.01
16.7

1iai
1/61
3.0
0.15
33.3

1ibg
30/61
0.0
−0.08
16.7

1igc
55/61
−1.3
0.02
0.0

1igf
6/61
1.3
0.10
16.7

1igi
2/61
1.9
0.13
16.7

1igm
43/61
−0.5
0.09
0.0

1ikf
1/61
2.7
0.35
50.0

1ind
20/61
0.5
0.31
0.0

1jel
1/61
3.5
0.02
16.7

1jhl
1/61
3.1
0.16
33.3

1kem
1/61
2.8
0.26
33.3

1mam
1/61
2.0
0.31
33.3

1mcp
60/61
−1.8
0.23
0.0

1mfa
4/61
1.6
−0.01
16.7

1mlb
23/61
0.4
−0.19
0.0

1mrd
2/61
2.6
0.28
33.3

1nbv
36/61
−0.2
0.03
0.0

1ncb
12/61
0.9
−0.02
16.7

1ngq
13/61
0.8
0.00
0.0

1nmb
1/61
2.9
0.38
16.7

1nsn
8/61
1.0
−0.11
0.0

1opg
61/61
−2.7
−0.01
0.0

1plg
2/61
1.8
−0.02
16.7

1rmf
3/61
1.8
0.01
16.7

1tet
2/61
2.6
−0.15
16.7

1ucb
1/61
4.3
0.33
16.7

1vfa
1/61
2.5
0.20
16.7

1vge
3/61
2.1
0.03
16.7

1yuh
22/61
0.6
−0.02
16.7

2cgr
1/60
2.5
0.17
16.7

2fb4
19/61
0.4
−0.12
0.0

2fbj
53/61
−0.9
0.12
16.7

2gfb
1/61
2.9
0.28
16.7

3hfl
12/61
1.0
0.20
16.7

3hfm
61/61
−4.3
−0.16
0.0

6fab
7/61
1.4
0.02
0.0

7fab
23/61
0.3
0.01
0.0

8fab
1/61
3.1
0.04
16.7

lmds
1b0nB
1/498
3.9
0.05
10.2

1bba
496/501
−2.1
0.03
18.0

1ctf
32/498
1.6
0.03
14.3

1dtk
206/216
−1.7
0.04
14.3

1fc2
247/501
0.0
0.05
16.0

1igd
342/501
−0.4
0.14
10.0

1shfA
223/438
0.1
−0.01
4.7

2cro
2/501
2.3
0.22
14.0

2ovo
110/348
0.4
0.05
11.8

4pti
26/344
1.5
0.02
14.7

smd3
1/501
3.5
0.01
10.0

Native rank = (rank of native structure with given PDB ID)/(total number of decoys); a rank of 1 is optimal and means the calculated energy of the native structure is lower than that of all its decoys; = correlation coefficient; FE = fractional enrichment.

The plots of energy versus rmsd in Figure 3, based on 4 native proteins and their collections of decoys evaluated with a varied selection of four-body potentials that we generated, are illustrative of the strengths and weaknesses of the performance measures defined above. In particular, since 4state_reduced is known to contain native-like alternative conformations for each protein in the set, reasonably good correlation () and fractional enrichment (FE > 10%) are expected from a reliable energy function [18, 35], and this is illustrated by the plot for 4pti. Next, ig_structal_hires and hg-structal contain decoys built by homology modeling for immunoglobulin (ig) and globin (hg) proteins, all of which are native-like structures with very low rmsds relative to native [18]. The plots for 1fvc and 1hdaB reflect the expected strong correlation and fractional enrichment; additionally, despite the fact that the native protein and very low rmsd decoys all have a good chance of achieving the lowest energy conformation, both native proteins rank 1 for these examples. Finally, the set lattice_ssfit consists of decoys selected with an all-atom energy function and refined using coarse lattice models, and rmsd > 4 Å for all decoys in this set relative to their native proteins [18, 36]. The plot for 1beo shows that, as expected in such cases of significantly misfolded decoys, there is no correlation between energy and rmsd relative to the native structure; furthermore, the fractional enrichment is low and the -score is relatively large, as commonly encountered by such decoys and suggested by the plot.

Figure 3: Sampling of calculated energy versus rmsd plots for four decoy sets. A different atomic four-body statistical potential energy function (i.e., distinct pairs of atomic alphabet size and tessellation edge-length cutoff parameters) was selected to compute the energy values for each plot. The plots reveal wide variability in the number of alternative conformations for a given native structure based on decoy category, and they highlight the relative strengths and weaknesses of native rank, correlation coefficient , -score, and fractional enrichment (FE) as performance measures under a range of conditions, hence reinforcing their collective importance for evaluating energy functions.

3.2. Four-Body Potentials: Relative Performance

To effectively rank all twelve four-body potentials generated for this study, first we identified, for each of the 145 native proteins and their respective decoys, the best native rank and largest -score, correlation coefficient, and fractional enrichment values obtained, without regard to which potential yielded those optimal values of the performance measures. Next, for each potential separately, we counted the number of times (out of 145) that the potential either matched or singularly provided each optimal value recorded for a performance measure concerning a native protein and its set of decoys (Table 4, numbers above parentheses). For each performance measure, we then ranked these counts across all the potentials (Table 4, numbers in parentheses); subsequently for each potential separately, we averaged its rankings across the four performance measures (Table 4, next to bottom row). Finally, those averaged ranks were used to generate an overall ranking of the twelve four-body potentials (Table 4, bottom row). The ranking approach based on relative performance employed here and in the subsequent section was inspired by the technique described in Summa et al. [3] for comparing the performance of their potential to other related methods.

Table 4: Relative performance among twelve atomic four-body statistical potentials.

Alphabet size
4
4
4
4
8
8
8
8
20
20
20
20

Cutoff (Å)
4.8
8
12
None
4.8
8
12
None
4.8
8
12
None

Native rank
42
48
55
61
23
29
27
25
33
41
34
33

(4)
(3)
(2)
(1)
(12)
(9)
(10)
(11)
(7)
(5)
(6)
(7)

-score
25
9
18
37
1
4
1
1
15
21
2
11

(2)
(7)
(4)
(1)
(10)
(8)
(10)
(10)
(5)
(3)
(9)
(6)

10
9
6
24
4
6
7
9
5
28
21
16

(5)
(6)
(9)
(2)
(12)
(9)
(8)
(6)
(11)
(1)
(3)
(4)

FE
38
46
49
47
23
34
32
38
37
56
43
42

(7)
(4)
(2)
(3)
(12)
(10)
(11)
(7)
(9)
(1)
(5)
(6)

Average of ranks
4.5
5
4.25
1.75
11.5
9
9.75
8.5
8
2.5
5.75
5.75

Overall ranking
4
5
3
1
12
10
11
9
8
2
6
6

Numbers above parentheses in each row reflect how many decoy sets (out of 145) for which the given potential matches the best performance value achieved among all 12 potentials tested; numbers in parentheses are the rankings of the counts in that row; = correlation coefficient; FE = fractional enrichment.

In general, four-body potentials derived using a 4-letter atomic alphabet ranked highest, followed by those based on letters, while potentials generated using 8 atom types ranked poorly over all four choices of tessellation edge-length cutoff parameter values. Among the 4-letter alphabet potentials, the one based on full structure tessellations (i.e., no edge-length cutoff) outperformed that using a 12 Å cutoff; however, the latter case is preferable since, without a fixed cutoff, false-positive atomic quadruplet interactions are admitted into the analyses based on those tessellations. A satisfying solution to this dilemma is revealed in the subsequent section as these four-body potentials are compared to those developed by other research groups.

3.3. Relative Performance: Comparisons with Related Methods

Next, an approach similar to that described in the previous section is used to individually compare each of our 12 four-body potentials to those of a dozen related methods. Using the Decoys-‘R’-Us database, Summa et al. [3] compared their “atomic environment potential” to ten other well-known physics- and knowledge-based potentials, providing us with valuable raw data to make comparisons with our four-body potentials. In addition, Fogolari et al. [18] developed an energy function employing two centers of interactions per amino acid. They also used the Decoys-‘R’-Us database to evaluate its performance, and these data are also included in our evaluations. Out of the 145 decoy sets that we used for benchmarking our four-body potentials relative to one another, 129 sets overlap with those used by both of those studies and form the basis of comparisons reported here. Lastly, the ten related methods investigated by Summa et al. [3] for comparing relative performance and used by us for a similar purpose include the following: three taken from the AMBER force field (a simple van der Waals potential, a pairwise electrostatic potential term, and the sum of these two terms, the latter representing the entire nonbonded contact energy of a typical molecular mechanics force field without either an explicit or implicit solvent model) [19]; two taken from CHARMM19 (both a van der Waals and a coulombic term) [20]; the and potentials of Delarue and Koehl [21], and the potential of Koehl and Delarue [22]; and distance-dependent atomic potentials RAPDF [23] and DFIRE [24].

For each of the 129 decoy sets common to the studies, we obtained the raw performance data (i.e., native rank, -score, correlation, coefficient, and fractional enrichment as presented in Table 3 for one of the four-body potentials) generated by each of the twelve methods described above. Next, we selected one of our four-body potentials and included its raw performance data, for a total of 13 methods to be compared. With every decoy set, we identified the best native rank achieved and the largest values obtained for -score, correlation coefficient, and fractional enrichment, without regard to which of the 13 methods was responsible for each optimal measurement. For each of these 13 methods separately, we counted the number of times (out of 129) that the method either matched or singularly provided each optimal value recorded for a performance measure concerning a native protein and its set of decoys (Table 5, numbers to the left of those in parentheses). For each performance measure, we then ranked these counts across all 13 methods (Table 5, numbers in parentheses); subsequently for each method separately, we averaged its rankings across the four performance measures (Table 5, next to last column). Finally, those averaged ranks were used to generate an overall ranking of the 13 methods (Table 5, last column).

Table 5: Relative performance among the four-body potential derived using a 4-letter alphabet and 12 Å cutoff as parameters and twelve other state-of-the-art methods.

Native rank
-score

FE
Average of ranks
Overall ranking

4 letters/12 Å cutoff
47 (5)
17 (3)
8 (5)
27 (5)
4.5
3

Summa et al. [3]
97 (1)
44 (1)
7 (6)
19 (8)
4
2

United-atom vdW
(AMBER) [19]
25 (12)
0 (12)
4 (8)
10 (12)
11
12

Coulombic (AMBER) [19]
33 (8)
8 (6)
2 (11)
18 (9)
8.5
10

United-atom vdW + coulombic (AMBER) [19]
26 (11)
4 (9)
4 (8)
11 (11)
9.75
11

United-atom vdW (CHARM19) [20]
31 (10)
9 (4)
4 (8)
36 (3)
6.25
7

Coulombic (CHARM19) [20]
76 (2)
22 (2)
5 (7)
22 (7)
4.5
3

(Delarue and Koehl) [21]
70 (3)
7 (7)
0 (13)
17 (10)
8.25
9

(Koehl and Delarue) [22]
33 (8)
0 (12)
32 (2)
26 (6)
7
8

(Delarue and Koehl) [21]
20 (13)
1 (11)
1 (12)
9 (13)
12.25
13

RAPDF [23]
53 (4)
3 (10)
28 (3)
36 (3)
5
5

DFIRE [24]
38 (7)
9 (4)
33 (1)
51 (1)
3.25
1

Fogolari et al. [18]
40 (6)
5 (8)
15 (4)
39 (2)
5
5

Numbers not in parentheses in each column reflect how many decoy sets (out of 129) for which the given method matches the best performance value achieved among all 13 methods tested; numbers in parentheses are the rankings of the counts in that column; = correlation coefficient; FE = fractional enrichment.

The overall rankings in Table 5 reveal that our four-body potential, derived using a 4-letter alphabet and 12 Å cutoff as parameters, outperformed 9 other methods, tied in overall ranking (3rd) with the coulombic term from CHARMM19, and was outperformed by DFIRE and the “atomic environment potential” of Summa et al. The methodology described in the previous paragraph was repeated separately for each of the twelve four-body potentials that we investigated, and the overall rankings in each case are reported in the columns of Table 6. Four-body potentials employing a 4-letter alphabet again appear to be the most competitive, and in particular those derived using unmodified tessellations (i.e., no edge-length cutoff) and a 12 Å cutoff achieved the highest overall rankings (3rd) among all of the four-body potentials, each in comparison to the 12 related state-of-the-art methods. As mentioned in the prior section, the introduction of false-positive atomic quadruplet interactions into the analyses is a concern when edge-length cutoffs are not considered after tessellation. Given that both of these potentials are equally competitive when compared with the other related methods, we conclude that a 4-letter alphabet and 12 Å tessellation edge-length cutoff provide the best pair of parameters with which to derive a four-body potential for calculating protein structure energies and effectively distinguishing native folds from nonnative decoy structures.

Table 6: Overall rankings by separately comparing each four-body potential with twelve other methods.

Alphabet size
4
4
4
4
8
8
8
8
20
20
20
20

Cutoff (Å)
4.8
8
12
None
4.8
8
12
None
4.8
8
12
None

Four-body potential
10
4
3
3
12
11
12
11
11
9
10
11

Summa et al. [3]
2
2
2
2
2
2
2
2
2
2
2
2

United-atom vdW (AMBER) [19]
12
12
12
12
11
12
11
12
12
12
12
12

Coulombic (AMBER) [19]
9
10
10
9
9
9
9
9
9
9
9
9

United-atom vdW + coulombic (AMBER) [19]
11
11
11
11
10
10
10
10
10
11
11
10

United-atom vdW (CHARM19) [20]
6
7
7
7
6
6
6
6
6
6
6
6

Coulombic (CHARM19) [20]
3
3
3
3
3
3
3
2
3
2
3
2

ΔE (Delarue and Koehl) [21]
8
9
9
10
8
8
8
8
8
8
8
8

(Koehl and Delarue) [22]
7
8
8
8
7
7
7
7
7
7
7
7

(Delarue and Koehl) [21]
13
13
13
13
12
13
13
13
13
13
13
13

RAPDF [23]
5
5
5
5
5
5
5
5
5
5
5
5

DFIRE [24]
1
1
1
1
1
1
1
1
1
1
1
1

Fogolari et al. [18]
4
5
5
5
4
4
4
4
4
4
4
4

Note that the overall rankings in this column correspond to those in the final column of Table 5; the remaining columns in this table were obtained by repeating the data analyses that generated Table 5 with respect to each of the other 11 four-body potentials.

4. Discussion

Energy calculations for single protein chains have been the sole focus up to this point, so it has been appropriate to consider atomic alphabets based only on the four heavy atom types found in proteins: carbon, nitrogen, oxygen, and sulfur. As mentioned earlier in Methods, a training set of 1417 diverse structures of single chain proteins were used for deriving the four-body potentials; however, the atomic coordinate for these protein chains was each obtained from a distinct PDB coordinate file, and the vast majority of these files are for structures of proteins complexed to small molecular or peptide ligands. Therefore, in order to tessellate the entirety of each of these PDB files, an expansion of the atomic alphabet is necessary to accommodate all atom types. The fact that such tessellations have an important function will become apparent as an application is introduced for predicting target-ligand binding affinities.

4.1. Generalized Four-Body Potential: An Alphabet Incorporating All Atom Types

Given the impressive performance on protein structures by the four-body potential derived using a 4-letter atomic alphabet and 12 Å tessellation edge-length cutoff, we simply expanded the alphabet to 6 letters in order to include atoms found exclusively in molecular ligands: M = all metals and X = all nonmetals other than (, , , ). The atomic frequency data and total number of tetrahedra generated by tessellating the totality of the atomic coordinate data in the 1417 PDB structure files (hydrogen atoms and water molecules excluded, as discussed in Methods), after filtering out edges longer than 12 Å, are provided in Table 7. Since we are now working with a = 6 letter alphabet, the atoms at the four vertices of each tetrahedron of a tessellation represent one of = 126 atomic quadruplets. A retracing of the steps described in Methods yields the all-atom four-body statistical potential presented in Table 8. Note that 11 of the 126 atomic quadruplet types are not represented by any of the 36,406,467 tetrahedra obtained from the tessellations.

Table 7: Summary data for the 1417 PDB coordinate files.

Atom types
Count
Proportion

C (carbon)
3612988
0.633193

N (nitrogen)
969253
0.169866

O (oxygen)
1088410
0.190749

S (sulfur)
28502
0.004995

M (all metals)
2529
0.000443

X (all other nonmetals)
4299
0.000754

Total atom count
5705981

Total tetrahedron count
36406467

Table 8: All-atom four-body statistical potential derived using a 6-letter alphabet and a 12 Å cutoff.

Quad
Count

CCCC
4107297

CCCM
1924

CCCN
4142684

CCCO
6462239

CCCS
297980

CCCX
2996

CCMM
157

CCMN
3758

CCMO
6511

CCMS
2320

CCMX
15

CCNN
1871781

CCNO
8544461

CCNS
128008

CCNX
2159

CCOO
3686844

CCOS
205846

CCOX
4995

CCSS
15467

CCSX
148

CCXX
161

CMMM
29

CMMN
164

CMMO
293

CMMS
665

CMMX
1

CMNN
2643

CMNO
7243

CMNS
2610

CMNX
30

CMOO
9551

CMOS
1041

CMOX
77

CMSS
2052

CMSX
13

CMXX
6

CNNN
122810

CNNO
2117811

CNNS
16884

CNNX
631

CNOO
2981894

CNOS
99630

CNOX
2400

CNSS
4318

CNSX
38

CNXX
68

COOO
683049

COOS
38976

COOX
24064

COSS
4524

COSX
64

COXX
84

CSSS
320

CSSX
5

CSXX
4

CXXX
12

MMMM
83

MMMN
42

MMMO
31

MMMS
379

MMMX
0

—

MMNN
85

1.836638

MMNO
113

1.608913

MMNS
364

3.698853

MMNX
0

—

MMOO
320

2.311659

MMOS
104

3.104429

MMOX
3

2.386025

MMSS
254

5.375177

MMSX
2

3.791851

MMXX
0

—

MNNN
1048

0.520184

MNNO
1323

0.093906

MNNS
562

1.303999

MNNX
6

0.153922

MNOO
4193

0.544515

MNOS
352

0.74942

MNOX
31

0.515747

MNSS
793

2.985098

MNSX
5

1.305273

MNXX
9

2.683083

MOOO
5790

1.111435

MOOS
167

0.676269

MOOX
171

1.508056

MOSS
211

2.359752

MOSX
4

1.158007

MOXX
55

3.418848

MSSS
62

3.8869

MSSX
2

2.739925

MSXX
0

—

MXXX
16

NNNN
5639

NNNO
60175

NNNS
538

NNNX
39

NNOO
384854

NNOS
6209

NNOX
354

NNSS
319

NNSX
6

NNXX
7

NOOO
227156

NOOS
11871

NOOX
3214

NOSS
951

NOSX
13

NOXX
66

NSSS
35

NSSX
0

—

NSXX
0

—

NXXX
3

OOOO
61473

OOOS
5019

OOOX
9614

OOSS
331

OOSX
45

OOXX
144

OSSS
38

OSSX
3

OSXX
0

—

OXXX
5

2.624158

SSSS
11

2.686034

SSSX
0
0

—

SSXX
0
0

—

SXXX
0
0

—

XXXX
0
0

—

4.2. Application: Target-Ligand Binding Affinity Prediction

The four-body potential derived in the previous section can be used to calculate the energy of any macromolecular structure. First, the 3D atomic coordinates of the structure are each labeled using the 6-letter alphabet and those points are tessellated subject to a 12 Å cutoff, then each tetrahedron in the tessellation is scored according to the atomic quadruplet identified at its four vertices by referring to the four-body potential previously derived in Table 8, and finally, the scores of all the tetrahedra are added up to determine the energy of the structure. Using the notation (i.e., total potential) to refer to the energy of a structure calculated in this way, we empirically calculate target-ligand binding affinity in the following manner (Figure 4):

(1)Tessellate the entire macromolecular complex and calculate .(2)Tessellate only atomic coordinates for the target protein and calculate .(3)The calculated target-ligand binding affinity is given by the difference

The above formula is a simplified model that is valid in the case of small ligands for which tetrahedra formed at the protein interface dominate any purely internal quadruplet atomic interactions within the ligand; hence, the relative energy contribution of the ligand is negligible [37].

Figure 4: Visualization of a procedure based on a simplified model to calculate target-ligand binding affinity () with the four-body potential.

4.3. Example: Predicting HIV-1 Protease-Inhibitor Binding Energy

To validate the approach for empirically calculating binding affinity, PDB accession codes and experimental binding energies were obtained from Jenwitheesuk and Samudrala [25] for twenty-five HIV-1 protease-inhibitor complexes (Table 9); they converted experimental inhibition constants () to experimental binding energies (, Gibbs free energy of binding, in units of kcal/mol) by applying the equation , where is the gas constant (1.987 cal K⁻¹ mol⁻¹) and is the absolute temperature (room temperature, 300 K). By following the steps outlined in the previous section, we determined for each of these complexes and used it as the calculated binding energy. As shown in Figure 5, the experimental and calculated binding energies for these complexes were highly correlated (). From a subsequent search of the Binding MOAD database [26, 27], we identified 115 additional HIV-1 protease-inhibitor complexes with experimental structures in PDB, for which experimental inhibition constants are available (Table 9). As before, values were obtained and used for representing the calculated binding energy, and values were converted to experimental binding energies. The correlation remained robust () when these data were combined with those of the initial plot (Figure 5).

Table 9: PDB accession codes for 140 HIV-1 protease-inhibitor complexes.

Twenty-five complexes culled from Jenwitheesuk and Samudrala [25]

1gno
1hbv
1hef
1heg
1hih
1hiv
1hps
1hpv
1hte
1htf
1htg
1hvi
1hvj

1hvk
1hvl
1hvr
1hvs
1pro
1sbg
2upj
4hvp
4phv
5hvp
8hvp
9hvp

An additional 115 complexes obtained from the Binding MOAD database [26, 27]

1a30
1a8g
1a8k
1a94
1a9m
1ajv
1b6j
1b6k
1b6m
1b6p
1c70
1d4h
1d4i

1d4k
1d4l
1daz
1dmp
1dw6
1ebk
1ebw
1eby
1ebz
1ec1
1ec2
1ec3
1fej

1ff0
1fff
1ffi
1fg8
1fgc
1fqx
1g2k
1g35
1gnm
1gnn
1hpo
1hsg
1hwr

1hxb
1iiq
1izh
1k1t
1k1u
1k2b
1k2c
1k6p
1k6v
1lzq
1m0b
1met
1mrx

1msn
1mtr
1nh0
1odw
1ody
1ohr
1qbr
1qbs
1qbt
1qbu
1sdu
1t7k
1tcx

1vij
1vik
1z1h
1zj7
1zlf
1zp8
1zpa
2aod
2aog
2aqu
2avm
2avo
2bpv

2bqv
2cem
2cen
2f3k
2f80
2fgu
2idw
2ien
2ieo
2nxd
2nxl
2nxm
2o4s

2p3b
2psu
2psv
2pym
2pyn
2q54
2q55
2q5k
2q63
2qd6
2qd8
2qi1
2qi4

2qi7
2r38
2wkz
2wl0
3bgb
3cyw
3cyx
3d1x
3d1y
3d1z
3d20

(a)

(b)

Figure 5: Scatter plots of experimental versus calculated binding energy for (a) twenty-five HIV-1 protease-inhibitor complexes culled by Jenwitheesuk and Samudrala [25] and (b) a larger set of 140 such complexes which include the initial twenty-five, with the remainder obtained by searching the Binding MOAD database [26, 27]. Both experimental binding energies and crystallographic structures are available for these complexes, and the latter were required for calculating binding energy () as outlined in the text and in Figure 4.

5. Conclusion

In this study, we derived and evaluated twelve distinct atomic four-body knowledge-based statistical potentials for protein structure prediction, by altering two parameter values: atom type (4-, 8-, or 20-letter alphabets) and distance cutoff for atomic interactions (none, 12 Å, 8 Å, or 4.8 Å). The best potential employed a simple 4-letter atomic alphabet and considered any quadruplet of atoms to be interacting when they were all pairwise within 12 Å of each other. In a head-to-head comparison of methods using 129 benchmarks from the Decoys-‘R’-Us database, our potential ranked 3rd and was outperformed by only two out of twelve other state-of-the-art methods. In addition to its simplicity and relative accuracy, our method is faster and more efficient in general, with some of the other physics- and knowledge-based potentials used for comparison employing well over one hundred different atom types. Future plans for improvement include combining this four-body potential together with other knowledge-based potentials, as well as subsequently implementing them together in conjunction with statistical machine learning tools.

Conflicts of Interest

The author declares that there are no conflicts of interest regarding the publication of this paper.

References

H. Berman, K. Henrick, H. Nakamura, and J. L. Markley, “The worldwide protein data bank (wwPDB): ensuring a single, uniform archive of PDB data,” Nucleic Acids Research, vol. 35, supplement 1, pp. D301–D303, 2007.

D. Rykunov and A. Fiser, “New statistical potential for quality assessment of protein models and a survey of energy functions,” BMC Bioinformatics, vol. 11, article 128, 2010.

C. M. Summa, M. Levitt, and W. F. DeGrado, “An atomic environment potential for use in protein structure prediction,” Journal of Molecular Biology, vol. 352, no. 4, pp. 986–1001, 2005.

C. Zhang, S. Liu, H. Zhou, and Y. Zhou, “An accurate, residue-level, pair potential of mean force for folding and binding based on the distance-scaled, ideal-gas reference state,” Protein Science, vol. 13, no. 2, pp. 400–411, 2004.

D. Rykunov and A. Fiser, “Effects of amino acid composition, finite size of proteins, and sparse statistics on distance-dependent statistical pair potentials,” Proteins: Structure, Function and Genetics, vol. 67, no. 3, pp. 559–568, 2007.

E. Ferrada and F. Melo, “Effective knowledge-based potentials,” Protein Science, vol. 18, no. 7, pp. 1469–1485, 2009.

K. E. Johansson and T. Hamelryck, “A simple probabilistic model of multibody interactions in proteins,” Proteins: Structure, Function and Bioinformatics, vol. 81, no. 8, pp. 1340–1350, 2013.

M. de Berg, O. Cheong, M. van Kreveld, and M. Overmars, Computational Geometry: Algorithms and Applications, Springer-Verlag, Berlin, Germany, Third edition, 2008.

M. J. Sippl, “Boltzmann's principle, knowledge-based mean fields and protein folding. n approach to the computational determination of protein structures,” Journal of Computer-Aided Molecular Design, vol. 7, no. 4, pp. 473–501, 1993.

M. J. Sippl, “Knowledge-based potentials for proteins,” Current Opinion in Structural Biology, vol. 5, no. 2, pp. 229–235, 1995.

C. W. Carter Jr., B. C. Lefebvre, S. A. Cammer, A. Tropsha, and M. H. Edgell, “Four-body potentials reveal protein-specific correlations to stability changes caused by hydrophobic core mutations,” Journal of Molecular Biology, vol. 311, no. 4, pp. 625–638, 2001.

B. Krishnamoorthy and A. Tropsha, “Development of a four-body statistical pseudo-potential to discriminate native from non-native protein conformations,” Bioinformatics, vol. 19, no. 12, pp. 1540–1548, 2003.

M. Masso and I. I. Vaisman, “Accurate prediction of enzyme mutant activity based on a multibody statistical potential,” Bioinformatics, vol. 23, no. 23, pp. 3155–3161, 2007.

M. Masso and I. I. Vaisman, “AUTO-MUTE: Web-based tools for predicting stability changes in proteins due to single amino acid replacements,” Protein Engineering, Design and Selection, vol. 23, no. 8, pp. 683–687, 2010.

M. Masso and I. I. Vaisman, “AUTO-MUTE 2.0: a portable framework with enhanced capabilities for predicting protein functional consequences upon mutation,” Advances in Bioinformatics, vol. 2014, Article ID 278385, 2014.

R. K. Singh, “Delaunay tessellation of proteins: Four body nearest-neighbor propensities of amino acid residues,” Journal of Computational Biology, vol. 3, no. 2, pp. 213–221, 1996.

R. Samudrala and M. Levitt, “Decoys 'R' Us: a database of incorrect conformations to improve protein structure prediction,” Protein Science, vol. 9, no. 7, pp. 1399–1401, 2000.

F. Fogolari, L. Pieri, A. Dovier et al., “Scoring predictive models using a reduced representation of proteins: Model and energy definition,” BMC Structural Biology, vol. 7, article no. 15, 2007.

S. J. Weiner, P. A. Kollman, D. A. Case et al., “A new force field for molecular mechanical simulation of nucleic acids and proteins,” Journal of the American Chemical Society, vol. 106, no. 3, pp. 765–784, 1984.

B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M. Karplus, “CHARMM: a program for macromolecular energy, minimization, and dynamics calculations,” Journal of Computational Chemistry, vol. 4, no. 2, pp. 187–217, 1983.

M. Delarue and P. Koehl, “Atomic environment energies in proteins defined from statistics of accessible and contact surface areas,” Journal of Molecular Biology, vol. 249, no. 3, pp. 675–690, 1995.

P. Koehl and M. Delarue, “Polar and nonpolar atomic environments in the protein core: Implications for folding and binding,” Proteins: Structure, Function, and Bioinformatics, vol. 20, no. 3, pp. 264–278, 1994.

R. Samudrala and J. Moult, “An all-atom distance-dependent conditional probability discriminatory function for protein structure prediction,” Journal of Molecular Biology, vol. 275, no. 5, pp. 895–916, 1998.

H. Zhou and Y. Zhou, “Distance-scaled, finite ideal-gas reference state improves structure-derived potentials of mean force for structure selection and stability prediction,” Protein Science, vol. 11, no. 11, pp. 2714–2726, 2002.

E. Jenwitheesuk and R. Samudrala, “Improved prediction of HIV-1 protease-inhibitor binding energies by molecular dynamics simulations,” BMC Structural Biology, vol. 3, no. 2, 2003.

M. L. Benson, R. D. Smith, N. A. Khazanov et al., “Binding MOAD, a high-quality protein-ligand database,” Nucleic Acids Research, vol. 36, no. 1, pp. D674–D678, 2008.

L. Hu, M. L. Benson, R. D. Smith, M. G. Lerner, and H. A. Carlson, “Binding MOAD (Mother of All Databases),” Proteins: Structure, Function and Genetics, vol. 60, no. 3, pp. 333–340, 2005.

G. Wang and R. L. Dunbrack Jr., “PISCES: a protein sequence culling server,” Bioinformatics, vol. 19, no. 12, pp. 1589–1591, 2003.

R. A. Brualdi, Introductory Combinatorics, Prentice Hall, New Jersey, NJ, USA, 3rd edition, 1999.

C. B. Barber, D. P. Dobkin, and H. Huhdanpaa, “The quickhull algorithm for convex hulls,” Association for Computing Machinery. Transactions on Mathematical Software, vol. 22, no. 4, pp. 469–483, 1996.

E. F. Pettersen, T. D. Goddard, C. C. Huang et al., “UCSF Chimera—a visualization system for exploratory research and analysis,” Journal of Computational Chemistry, vol. 25, no. 13, pp. 1605–1612, 2004.

J. B. O. Mitchell, R. A. Laskowski, A. Alex, and J. M. Thornton, “BLEEP - potential of mean force describing protein-ligand interactions: I. Generating potential,” Journal of Computational Chemistry, vol. 20, no. 11, pp. 1165–1176, 1999.

G. Q. Dong, H. Fan, D. Schneidman-Duhovny, B. Webb, A. Sali, and A. Tramontano, “Optimized atomic statistical potentials: assessment of protein interfaces and loops,” Bioinformatics, vol. 29, no. 24, pp. 3158–3166, 2013.

M. J. O'Meara, A. Leaver-Fay, M. D. Tyka et al., “Combined covalent-electrostatic model of hydrogen bonding improves structure prediction with Rosetta,” Journal of Chemical Theory and Computation, vol. 11, no. 2, pp. 609–622, 2015.

B. Park and M. Levitt, “Energy functions that discriminate X-ray and near-native folds from well-constructed decoys,” Journal of Molecular Biology, vol. 258, no. 2, pp. 367–392, 1996.

Y. Xia and M. Levitt, “Extracting knowledge-based energy functions from protein structures by error rate minimization: Comparison of methods using lattice model,” Journal of Chemical Physics, vol. 113, no. 20, pp. 9318–9330, 2000.

M. Masso, “Four-body atomic potential for modeling protein-ligand binding affinity: Application to enzyme-inhibitor binding energy prediction,” BMC Structural Biology, vol. 13, no. 1, article no. S1, 2013.