Integrating water exclusion theory into β contacts to predict binding free energy changes and binding hot spots

Abstract

Background

Binding free energy and binding hot spots at protein-protein interfaces are two important research areas for understanding protein interactions. Computational methods have been developed previously for accurate prediction of binding free energy change upon mutation for interfacial residues. However, a large number of interrupted and unimportant atomic contacts are used in the training phase which caused accuracy loss.

Results

This work proposes a new method, β ACV _ASA , to predict the change of binding free energy after alanine mutations. β ACV _ASA integrates accessible surface area (ASA) and our newly defined β contacts together into an atomic contact vector (ACV). A β contact between two atoms is a direct contact without being interrupted by any other atom between them. A β contact’s potential contribution to protein binding is also supposed to be inversely proportional to its ASA to follow the water exclusion hypothesis of binding hot spots. Tested on a dataset of 396 alanine mutations, our method is found to be superior in classification performance to many other methods, including Robetta, FoldX, HotPOINT, an ACV method of β contacts without ASA integration, and ACV _ASA methods (similar to β ACV _ASA but based on distance-cutoff contacts). Based on our data analysis and results, we can draw conclusions that: (i) our method is powerful in the prediction of binding free energy change after alanine mutation; (ii) β contacts are better than distance-cutoff contacts for modeling the well-organized protein-binding interfaces; (iii) β contacts usually are only a small fraction number of the distance-based contacts; and (iv) water exclusion is a necessary condition for a residue to become a binding hot spot.

Conclusions

β ACV _ASA is designed using the advantages of both β contacts and water exclusion. It is an excellent tool to predict binding free energy changes and binding hot spots after alanine mutation.

Background

A binding hot spot is a small area in a protein binding interface whose mutation can lead to a big change in binding free energy. The determination of its accurate location in the interface is a fundamental problem in structural biology, and is useful for applications such as rational drug design and protein engineering [1]. In wet labs, a residue’s contribution to binding free energy can be determined through mutation experiments. For example, alanine scanning mutagenesis [2] mutates interfacial residues individually into an alanine, and then measures the change of binding free energy (Δ Δ G) to quantify the contribution of the side chain of the mutated residue. Based on these wet-lab experimental outcomes and databases [3–6], it has been reported that binding free energy is unevenly distributed in protein interfaces [7]. In fact, there are always a small fraction of interfacial residues—the binding hot spot—which make major contribution to the binding [7, 8] with Δ Δ G≥2 kcal/mol [3]. But wet-lab experiments are both time and cost expensive. Reliable computational methods are thus needed for accurate prediction of binding free energy change.

FoldX [9, 10], Robetta [11, 12] and CC/PBSA [13] are some well-known physics-based methods for this prediction problem. These methods use empirical terms (such as hydrogen bonds), the van der Waals terms and Coulomb electrostatics to learn a linear function for estimating the effect on the change of binding free energy after residue mutations. However, the predicted energy by these methods has a large discrepancy from experimentally measured Δ Δ G[14]. Thus, other methods have been proposed to qualitatively identify binding hot spots. For example, protein sequences are used by [15] and ISIS [16], while protein tertiary structures are used together with docking techniques by [17]. Protein quaternary structures have been also widely used [18]. For example, Hotsprint [19] and HotPOINT [20] generate rules to identify binding hot spots from features such as conservation, accessible surface area (ASA), residue propensity and/or residue pairwise potentials. Machine learning models are also widely used for predicting binding hot spots. Decision trees are used in MINERVA [14] to induce rules at different levels of protein information including structure, sequence and molecular interactions. Later, machine learning algorithms SVM and its ensemble are employed to combine energetic terms such as van der Waals potentials, solvation energy, hydrogen bonds and Coulomb electrostatics, and/or other protein sequences and structure information for a better hot spot prediction performance. Recently, Bayesian Networks are used to combine three main sources of information related to conservation, FoldX-calculated Δ Δ G and atomic contacts for a novel probabilistic model of binding hot spots prediction [21]. Very recently, random forests have been proposed to predict hot spots [22] by using structural neighborhood properties of mutated residues and other conventional physicochemical features [23, 24]. Besides alanine mutations, hot spots after mutations to any other type of residues are also investigated [6] and their binding free energy changes can be predicted [13, 25] with good performance. Several of these methods are also assessed in a community-wide test for predicting mutation effects on protein-protein interaction affinity [26].

In spite of intensive research, the prediction still needs a big improvement. The existing methods usually used those atomic contacts based on Voronoi diagram or simply defined by a distance threshold with little consideration on the local atomic organization of the contacts. If the distance threshold is too large, e.g., larger than 6 Å, an atomic contact between two atoms i and j may have no direct contact area, because the space between i and j can accommodate other atoms. Such interrupted contacts constitute a large proportion of the traditionally used contacts. It is highly questionable whether they are really important to protein binding. In fact, important contacts in hot spot prediction [10, 11] or those closely related to binding hot spots [14] are generally not interrupted, such as hydrogen bonds, salt bridges and π−π contacts.

To overcome these drawbacks, we propose a novel classifier βACV _ASA for predicting Δ Δ G and binding hot spots. The main idea of βACV _ASA is to use atomic contact vector (ACV) of β contacts (that’s why our classifier is named βACV for short) instead of distance-cutoff contacts. β contact, found on β-skeletons [27], is our newly defined contact [28]. A β contact between two atoms restricts that there is no other atoms between these two atoms, and requires that the two atoms should have enough direct contact area to form an interaction. The definition of β contacts can filter out a lot of unimportant and interrupted distance-cutoff contacts. Our analysis has found that β contacts are only a small fraction number of those contacts based on a distance threshold [28], but they are effective to distinguish crystal packing from homodimers [28] and to predict protein-ligand binding affinity [29].

Another important idea is that the relative ASA properties are integrated by our βACV classifier based on the water exclusion hypothesis of binding hot spots. The water exclusion hypothesis states that the topological shape of a binding hot spot and its surrounding residues can be characterized as an O-ring structure [3]. Few residues on the O-ring, which are largely exposed to bulk solvent water, can contribute significantly to the protein binding. Thus in βACV, the energy contribution of a β contact to protein binding is required to be inversely proportional to its ASA.

Our βACV _ASA was tested on a dataset of 396 alanine mutations to show its superior performance. We compared βACV _ASA with the following methods: (i) ACV methods using distance-cutoff contacts to reveal the importance of β contacts to protein binding; (ii) a βACV method without ASA integration to confirm whether the water exclusion theory is necessary for binding hot spots; and (iii) several widely-used state-of-the-art methods such as Robetta, FoldX, HotPOINT and KFC to show the overall better prediction capability of βACV _ASA .

Methods

Dataset

The data stored in the ASEdb database [3] and the mutations in BID [4] having Δ Δ G measurements are both used for evaluating our method. In total, our dataset contains 22 protein-protein complexes (detailed in Additional file 1: Table S2). All of them have quaternary structures in PDB and meet the following three requirements. First, no redundancy exists among these protein complexes. Given two protein complexes (e.g., interacted pair A and B, and interacted pair C and D), a sequence identity is calculated through BLAST with the default setting for A and C, A and D, B and C, and B and D, denoted by S(A,C), S(A,D), S(B,C) and S(B,D) respectively. These two protein complexes are redundant if S(A,C)≥40% and S(B,D)≥40%, or S(A,D)≥40% and S(B,C)≥40%. According to this criterion, most of the protein complexes in our dataset are non-redundant. For those redundant complexes, the mutations in the similar proteins must be in different positions. Our requirement on this sequence identity is reasonable, since atomic contacts and ASA used in this work are derived from complexes only, rather than from sequences (such as required by conservation scores). Secondly, only alanine mutations are considered. Thirdly, mutated atoms before mutation must have at least one distance-cutoff atomic contact with the partner proteins. The mutated atoms are those atoms except N, CA, C, O and CB. Under these requirements, our dataset has 396 alanine mutations (detailed in Additional file 1: Table S3). Of these mutations, 86 are binding hot spot residues having Δ Δ G≥2 kcal/mol.

Atomic βcontacts in protein binding interfaces

Atomic β contact is a recently proposed notion of atomic contacts for modeling the well-organized protein 3D structures [28]. Its detail can be found in [28]. For easy reference, we give a brief description of β contacts and a simple method to produce β contacts from a protein quaternary structure.

Atomic βcontacts: a definition

Given a quaternary structure of a protein complex p, a β contact between two atoms i and j in p requires that (i) the spatial distance between i and j is less than a threshold T _d plus the sum of their van der Waals radii defined by [30] (distance-cutoff contacts for short), (ii) i and j share a Voronoi facet in p’s Voronoi diagram, and (iii) the contact cannot break p’s β-skeleton. The β-skeleton [27] of a discrete set p is an undirected graph in computational geometry. In this graph, two points i and j have an edge if angle ikj is sharper than a threshold determined by β, ∀k∈p,k≠i,j. This angle threshold is denoted as ∠ β, which actually defines a forbidden region fr of the contact between i and j, e.g., the gray regions in Figure 1. When β=1, namely ∠ β=90°, fr is the sphere with the midpoint of i and j as the center and with c’s length as the diameter as shown in Figure 1(b). This sphere is similar to van der Waals radii of atoms. The forbidden region fr of a β contact usually does not cover any other atoms. Otherwise, if there is an atom k in fr, for example as shown in Figure 1(a) when ∠ β=90° or in Figure 1(c) when ∠ β=75°, the contact between i and j is not a β contact. A β contact suggests that its two atoms should have enough direct contact area to form an important interaction. The number of atomic β contacts in protein binding interfaces is only a small fraction number of distance-based contacts or less than half the number of contacts in the Voronoi diagrams when T _d =3.3 as found by [28]. Interestingly, the use of β contacts can achieve better prediction performance for distinguishing false binding of crystal packing from homodimers.

Figure 1

Examples of β contacts and non- β contacts. Three points, denoted by i, j and k, represent atoms. The dashed circles represent the van der Waals spheres in 2D space. The lines in yellow are of interest.

A method to produce βcontacts

A protein complex p can be represented as an atomic β contact graph b(p), if all of the heavy atoms are represented by nodes, and the β contacts are represented by edges. To produce b(p) for p, Qhull is first used to obtain the Delaunay triangulation [31] for all nodes. After that, the distance threshold T _d is used to remove those atomic contacts whose distances are too large. T _d is set as 3.3 Å (the diameter of a water molecule 2.8 Å plus 0.5 Å). This threshold is an insensitive factor to β contacts when it is large enough. Please refer to the Additional file 1 for an analysis of β contacts under several different T _d s. Thirdly, each atomic contact is checked to guarantee that it satisfies the β skeleton requirements. To sharpen the difference of those mutations with higher Δ Δ G and those with lower Δ Δ G, the angle threshold ∠ β is set as 75° in this work, whose forbidden region fr is larger than that of ∠ β=90° as shown in Figure 1(b) and (c). That is, ∠ β=75° is a stricter condition than using ∠ β=90° to produce β contacts. The rationale to choose the stricter condition ∠ β=75° is illustrated in the following situation. Assume (i) A, B and C are the center points of three atoms with van der Waals radii 1.8 Å (for example, the radius of common Carbon atoms in protein structures), (ii) there is a covalent bond between A and B (AB for short) with spatial distance 1.5 Å, and a non-covalent bond between A and C (AC for short) and one between B and C (BC for short), and (iii) the van der Waals sphere of A and that of C are circumscribed each other with the spatial distance 3.6 Å, and the same for the van der Waals sphere of A and that of C. Then, the angle ∠ A B C=78°. A stricter threshold than 78° is 75°, which was chosen by this work.

Our prediction methods

This section describes how to construct our βACV _ASA classifier, including how to define βACV and how to integrate ASA into βACV.

βACV: a vector representation for interfacial alanine mutations

We use an atomic contact vector (ACV) [32] of β contacts to represent an interfacial alanine mutation. We also use β contacts of the interfacial bound water molecules to update the vector, and use the atomic environment of the mutation to expand the basic vector.

Constructing a basic βACV vector

To produce a basic βACV for an alanine mutation m u t r ala of residue r in a protein complex p, we build the β contact graph b(p) for p. We then remove the coordinates of the mutated atoms in r to get a quaternary structure resulted from the mutation, denoted by p r mut. Then, we produce another β contact graph b(p r mut) for p r mut.

By the mutation m u t r ala, some contacts in b(p) may disappear in b(p r mut), namely those mutated contacts, while some new contacts can be formed in b(p r mut), called new contacts. Both of those mutated and new contacts are represented in a βACV vector. As the heavy atoms from the 20 standard residues are grouped into 8 atomic types (shown in Additional file 1: Table S4) by this work, our βACV vector has 36 pairs of atomic types as elements. The value v(T i,T j) of each element in βACV with atomic types Ti and Tj is calculated, using Equation 1.

$$\begin{array}{ll}v(\mathit{\text{Ti}},\mathit{\text{Tj}})=& \sum _{(x,y)\in M(\mathit{\text{Ti}},\mathit{\text{Tj}})}\frac{1}{d_{(x,y)}^2}-\sum _{(x^{\prime },y^{\prime })\in N(\mathit{\text{Ti}},\mathit{\text{Tj}})}\frac{1}{d_{(x^{\prime },y^{\prime })}^2}\end{array}$$

(1)

where x and x^′ are of the atomic type Ti, y and y^′ are of the atomic type Tj, (x,y) and (x^′,y^′) are two atomic pairs, d(∗,∗) is the spatial distance between a pair of atoms, and M(T i,T j)(or N(T i,T j)) is the set of all those mutated (or the set of all those new) contacts whose atomic types are Ti and Tj. Here term d²(∗,∗) is specially used to follow the same idea as Coulomb’s law which also uses the inverse of squared distance. Note that the other common contacts between b(p) and b(p r mut) are not used in βACV. Alanine mutations of Ala are assumed to have insignificant Δ Δ G and alanine mutations of Gly are not considered.

It can be seen that a basic βACV considers all mutated contacts and new contacts, including both acrossinterface contacts and those contacts from the same proteins or same biological units. However, atomic contacts between covalent-bond nearby atoms are not used in Equation 1. The covalent-bond nearby atoms of a given atom i are those atoms that have not more than three covalent-bond steps from i. For example, suppose i−j−k−l−m, where − indicates a covalent bond. From i, the covalent-bond step is 0 to i, is 1 to j, is 2 to k, is 3 to l and is 4 to m, respectively. Thus, i,j,k and l are covalent-bond nearby atoms of i, while m is not. In βACV, the contacts between i and its covalent-bond nearby atoms are excluded from M or N in Equation 1. This is reasonable, because spatially close distances between i and its covalent-bond nearby atoms are more likely due to the rigidity of their covalent bonds.

Bound water molecules in protein interfaces

Protein folding and binding occur in a solvent environment in vivo. Water molecules are heavily involved in protein binding and sometimes they can form a compulsory part of the protein interfaces. In this work, a water molecule in PDB is considered as a part of a binding interface if (i) it has at least 3 potential hydrogen-bonds contacts, or (ii) it has 2 potential hydrogen-bond contacts and also has at least 2 other contacts with spatial distances less than 4 Å. A potential hydrogen-bond contact is required to have a spatial distance less than 3.2 Å between a hydrogen donor (such as a nitrogen atom) and a hydrogen acceptor (such as an oxygen atom). Water molecules under this requirement, named bound water molecules, are such closely involved in protein folding and binding that they can play an integral part. Bound water molecules are then grouped into the Oxygen atomic type with more than one hydrogen atom (shown in Additional file 1: Table S4) to update the values of the elements in the basic βACV vector. We did not consider the contacts between any two water molecules.

The neighbourhood atoms of mutated residues

Information of neighbourhood atoms of m u t r ala is used to expand the basic βACV vector. Assume that S is a set of atoms which have β contacts with the mutated atoms under ∠ β=90°. For each atom in S including the mutated atoms, its nearby atoms are added into S. Then, an atomic vector with the above 8 atomic types is also used to represent those atoms in S in the bound state. The value of its element Tk is calculated using v(T k)^b in Equation 2. Similarly, v(T k)^u in Equation 2 is used to calculate another atomic vector for representing the atoms in S in the unbound state.

$$v{\left(\mathit{\text{Tk}}\right)}^b=\sum _{j\in S,t_j=\mathit{\text{Tk}}}{E}_j^{\mathit{\text{loc}}};v{\left(\mathit{\text{Tk}}\right)}^u=\sum _{j\in S,t_j=\mathit{\text{Tk}}}{E}_j^{\mathit{\text{loc}}}\left(u\right)$$

(2)

where E j loc (or E j loc(u)) is the relative local ASA of atom j in the bound (or unbound) state calculated via Equation 6 below. Water molecules were not considered here. Thus, each basic βACV vector is now expanded by another 16 atomic types for representing surrounding information of mutated atoms. So, the expanded βACV is a vector representation with 52 elements.

βACV _ASA : integrating the water exclusion hypothesis into βACV

Solvent water is compulsory for protein binding, but water exclusion—small accessible surface area (ASA)—is a necessary condition for a residue to become binding hot spot [3, 33, 34]. Few highly exposed residues can make significant contribution to protein binding strength [34]. Thus, we integrate ASA information into each atomic pair of βACV in Equation 1, and name the method βACV _ASA . We note that except Equation 1, the other definitions in βACV _ASA are the same as those in βACV.

Given a protein complex p, we take the following steps to integrate the water exclusion theory into Equation 1. The first step is to use NACCESS [30] to produce ASA for all of the atoms and residues in both bound and unbound states. For p in the bound state, we then define special ASA terms for an atom i using Equation 3, and for a residue R _i using Equation 4 and Equation 5.

$$\begin{array}{llll}{E}_i& =\sqrt{\frac{{\mathit{\text{ASA}}}_i}{50}}& B_i& =\mathit{\text{max}}(0,1-E_i)\end{array}$$

(3)

$$\begin{array}{llll}{E}_{R_i}^{\mathit{\text{bb}}}& =\sqrt{\frac{{\mathit{\text{ASA}}}_{R_i}^{\mathit{\text{bb}}}}{\mathit{\text{max}}\left({\mathit{\text{ASA}}}_{R_i}^{\mathit{\text{bb}}}\right)}}& B_{R_i}^{\mathit{\text{bb}}}& =\mathit{\text{max}}(0,1-E_{R_i}^{\mathit{\text{bb}}})\end{array}$$

(4)

$$\begin{array}{llll}{E}_{R_i}^{\mathit{\text{sc}}}& =\sqrt{\frac{{\mathit{\text{ASA}}}_{R_i}^{\mathit{\text{sc}}}}{\mathit{\text{max}}\left({\mathit{\text{ASA}}}_{R_i}^{\mathit{\text{sc}}}\right)}}& B_{R_i}^{\mathit{\text{sc}}}& =\mathit{\text{max}}(0,1-E_{R_i}^{\mathit{\text{sc}}})\end{array}$$

(5)

In Equation 3, ASA _i is accessible surface area of atom i, while E _i is its relative ASA, and B _i is its relative ASA burial compared to the maximum ASA, where number 50 is roughly half of NACCESS-calculated ASA of a single water molecule without any neighbor atoms (the ASA of a water molecule is 98.47=4×3.14×2.8² Ų). In Equations 4 and 5, ${\mathit{\text{ASA}}}_{R_i}^{\mathit{\text{bb}}}$ and ${\mathit{\text{ASA}}}_{R_i}^{\mathit{\text{sc}}}$ are accessible surface area of backbone atoms (i.e., bb) and of side-chain atoms (i.e., sc) for a residue R _i , while $E_{R_i}^{\ast }$ and $B_{R_i}^{\ast }$ are the relative ASA and the relative ASA burial of ∗∈{b b,s c}. $\mathit{\text{max}}\left({\mathit{\text{ASA}}}_{R_i}^{\mathit{\text{bb}}}\right)$ and $\mathit{\text{max}}\left({\mathit{\text{ASA}}}_{R_i}^{\mathit{\text{sc}}}\right)$ are the maximum ASA of backbone atoms and of side-chain atoms for the residue type of R _i , which are calculated in a triplet of ALA- R _i -ALA by NACCESS. These backbone atoms and side-chain atoms are defined in the same way as those in [30].

We compute the local ASA E i loc and local ASA burial B i loc of an atom i via Equation 6.

$$\begin{array}{ll}{E}_i^{\mathit{\text{loc}}}& =\left\{\begin{array}{c}{E}_i\times E_{R_i}^{\mathit{\text{bb}}}\text{if}i\text{is a backbone atom of}{R}_i\\ E_i\times E_{R_i}^{\mathit{\text{sc}}}\text{if}i\text{is a side-chain atom of}{R}_i\end{array}\right.\\ B_i^{\mathit{\text{loc}}}& =\left\{\begin{array}{c}{B}_i\times B_{R_i}^{\mathit{\text{bb}}}\text{if}i\text{is a backbone atom of}{R}_i\\ B_i\times B_{R_i}^{\mathit{\text{sc}}}\text{if}i\text{is a side-chain atom of}{R}_i\end{array}\right.\end{array}$$

(6)

where the multiplication of relative ASA burial of both atom i and its residue is used to calculate local ASA burial B i loc. This is because relative ASA of both an atom and its residue are critical in describing the accessibility of an atom. For example, an atom may be buried with small ASA but its covalent-bond atoms might be exposed. When relative ASA of atoms or residues are used individually, the performance was worse (data not shown).

To integrate water exclusion theory into Equation 1, we determine the value v(T i,T j) of each element in βACV _ASA through Equation 7 instead of Equation 1.

$$\begin{array}{c}v(\mathit{\text{Ti}},\mathit{\text{Tj}})=\sum _{(x,y)\in M(\mathit{\text{Ti}},\mathit{\text{Tj}})}\frac{B_x^{\mathit{\text{loc}}}\times B_y^{\mathit{\text{loc}}}}{d_{(x,y)}^2}\\ -\sum _{(x^{\prime },y^{\prime })\in N(\mathit{\text{Ti}},\mathit{\text{Tj}})}\frac{B_{x^{\prime }}^{\mathit{\text{loc}}}\times B_{y^{\prime }}^{\mathit{\text{loc}}}}{d_{(x^{\prime },y^{\prime })}^2}\end{array}$$

(7)

where T_∗, x, y, x^′, y^′, M and N have the same meaning as those in Equation 1.

Comparison of βcontacts with distance-cutoff contacts

To compare the performance of β contacts with distance-cutoff contacts for predicting Δ Δ G, ACV _ASA based on distance-cutoff contacts is constructed in a similar way to constructing βACV _ASA . To further show the importance of β contacts in protein binding interfaces, non βACV _ASA is also constructed for alanine mutations at the setting of ∠ β=90°. In non βACV _ASA , the values of its elements are the difference of the values of the 52 elements between ACV _ASA and βACV _ASA . To highlight the advantage of β contacts, ACV _ASA is also evaluated with different spatial distance thresholds (from 2.9 Å to 5 Å) for defining atomic contacts across interfaces and within binding sites.

Our βACV _ASA classifier and its variants described above are summarized in Table 1.

Table 1 Description of β ACV _ASA and its variant methods

Ridge regression: predict Δ Δ Gand binding hot spots

Ridge regression in Matlab is used here to learn a relation between atomic contact vectors and Δ Δ G. By this regression, values in each column are normalized for the training dataset. Ridge regression minimizes average square error SE between the experimental ($\mathrm{\Delta \Delta }{G}_i^e$) and the predicted $\mathrm{\Delta \Delta }{G}_i^p$ in the training data with N mutations where $\mathit{\text{SE}}=\frac{\sum _i{(\mathrm{\Delta \Delta }{G}_i^e-\mathrm{\Delta \Delta }{G}_i^p)}^2}{N-1}$.

In our evaluation, leave-one-out cross-validation is used for all of the 396 mutations, and then the correlation coefficient R and average standard deviation $\delta =\sqrt{\mathit{\text{SE}}}$ are calculated. Under this evaluation framework, there is one outlier prediction by β ACV _ASA and one outlier by βACV for the whole dataset with 396 mutations. These outliers have less than -3 kcal/mol predicted Δ Δ G, or more than 11 kcal/mol predicted Δ Δ G, as shown in Additional file 1: Table S3. This may be due to limited alanine mutations of a high Δ Δ G in the dataset.

In this work, predicted hot spot residues are those residue mutations with a predicted Δ Δ G≥2 kcal/mol, same as the true hot spot definition.

Hot spot prediction and evaluation measures

βACV _ASA is also assessed by applying to the classification problem of binding hot spots. Classification performance is measured by p r e c i s i o n(p.), r e c a l l(r.), a c c u r a c y(a c c.) and F1 whose definitions are given in Equation 8.

$$\begin{array}{c}\mathit{\text{precision}}\left(\mathit{p.}\right)=\frac{\mathit{\text{TP}}}{\mathit{\text{TP}}+\mathit{\text{FP}}}\\ \mathit{\text{recall}}\left(\mathit{r.}\right)=\frac{\mathit{\text{TP}}}{\mathit{\text{TP}}+\mathit{\text{FN}}}\\ \mathit{\text{accuracy}}\left(\mathit{acc.}\right)=\frac{\mathit{\text{TP}}+\mathit{\text{TN}}}{\mathit{\text{TP}}+\mathit{\text{TN}}+\mathit{\text{FP}}+\mathit{\text{FN}}}\\ F1=\frac{2\times \mathit{\text{precision}}\times \mathit{\text{recall}}}{\mathit{\text{precision}}+\mathit{\text{recall}}}\end{array}$$

(8)

where binding hot spots are considered as the true cases, while non-hot spots as the false cases; TP, FP, TN and FN are true positives, false positives, true negatives and false negatives, respectively. Hence, precision is the number of correct hot spot predictions divided by the number of positive predictions, recall is the fraction of correct hot spot predictions over all hot spots, while accuracy is the number of correctly predicted hot spots and non-hot spots divided by the number of all mutations. These measures are also used in [14, 20, 35] with the same definitions.

Results and discussion

β contacts are better than distance-cutoff contacts for predicting Δ Δ G

Our β ACV _ASA classifier is compared with ACV _ASA and with non βACV _ASA to show the importance of β contacts in the prediction of Δ Δ G under alanine mutations. The prediction results are presented in Figure 2(a), (b) and (c).

Figure 2

Δ Δ G predicted by different methods. In (a) - (d), ‘o’ represents non-polar residues, while ‘x’ represents polar residues. R is the Pearson correlation coefficient and δ is the average standard deviation. R is specially calculated for (a) and (d) based on the data set after removing the one outlier prediction. R is slightly changed to 0.543 or 0.515 respectively when the outlier is not removed. The two diagonal red lines represent y=x±1.5.

It can be seen from Figure 2(a) and (b) that βACV _ASA has a better Δ Δ G prediction performance according to both correlation coefficient R and average standard deviation δ. The number of β contacts used by βACV _ASA is only a small fraction of the number of distance-cutoff contacts used by ACV _ASA . For example, there are 54,286 distance-cutoff contacts across binding interfaces for the 22 protein complexes, but there are only 9,830 β contacts across the binding interfaces (∠ β=90°), and 4,096 β contacts under the setting of ∠ β=75° which is actually used by βACV _ASA . So, βACV _ASA uses only 7.55% number of distance-cutoff atomic contacts but it achieves a better prediction performance.

The comparison between βACV _ASA and non βACV _ASA (Figure 2(a) and (c)) further suggests the importance of β contacts in Δ Δ G prediction. In Figure 2(c), non βACV _ASA has much lower R (0.469) and a higher δ (1.486) than βACV _ASA , but non βACV _ASA uses all non- β contacts of ∠ β=90°, that is, 81.8% number of distance-cutoff atomic contacts.

A lot of alanine mutations are not binding hot spot residues, having a small Δ Δ G, i.e., <2 kcal/mol. These mutations heavily affect the calculation of R and δ. On the other hand, the prediction of residue mutations with a high Δ Δ G is more important. Thus, the classification performance for these binding hot spots is also assessed. The results are shown in Table 2. It is noted that F1 is not the objective function to be optimized in the regression process.

Table 2 Prediction performances by different methods for the same set of binding hot spots

In Table 2 among βACV _ASA , ACV _ASA and non βACV _ASA , βACV _ASA has the highest precision, recall, F1 and accuracy, while non βACV _ASA has the lowest. For example, βACV _ASA ’s F1 is 0.604, 0.122 higher than non βACV _ASA ’s F1. However, ACV _ASA and non βACV _ASA show quite similar performances. The reason would be that non- β contacts (used by non βACV _ASA ) are often dominant in distance-cutoff contacts (used by ACV _ASA ).

The performances of ACV _ASA with different spatial distance thresholds are shown in Table 3. We can see that the performance has a growing tendency when the threshold increases. Nevertheless, the best performance of ACV _ASA (F1 = 0.5 in Table 2) is much lower than that of βACV _ASA (F1 = 0.604). Of special interest, when the spatial distance threshold is set at 3.6 Å, ACV _ASA has a number of distance-cutoff contacts nearly the same as the number of our used β contacts. In this special case of almost the same number of contacts used, ACV _ASA has much worse performance than βACV _ASA , and only about half of the distance-based contacts are β contacts across the 22 protein-protein binding interfaces. These results affirm that β contacts are advantageous over distance-based contacts for predicting binding hot spot residues.

Table 3 Prediction performance and the numbers of used contacts by β ACV _ASA and ACV _ASA

Water exclusion is a necessary condition of hot spot binding

Literature work has reported that water exclusion is a necessary condition for an interfacial residue to become a hot spot residue [3, 33]. To confirm the importance of water exclusion in the prediction of Δ Δ G, the performance by βACV when ASA is not integrated is assessed. This prediction performance is shown in Figure 2(d) and Table 2. Comparing Figure 2(d) with Figure 2(a), βACV _ASA has a better regression performance than βACV indeed. It also has a better F1 performance as seen in Table 2. These results confirm that water exclusion plays an important role in hot spot prediction, and it should be a necessary condition for an interfacial residue to become a binding hot spot residue in protein-protein complexes.

Our method βACV _ASA is superior to several widely-used methods

Our βACV _ASA classifier is also assessed against the state-of-the-art methods FoldX [9, 10], Robetta [11], HotPOINT [20] and KFC [36]. The prediction performances of these previous methods were obtained through their web servers (Robetta, HotPOINT and KFC) or the standalone executable program (FoldX with default settings).

Comparison results

Figures 2(e) and 2(f) show that the prediction performance of FoldX and Robetta are much worse than our βACV _ASA . These two methods have a R of 0.324 or 0.485, much smaller than βACV _ASA ’s 0.569; their δ is 1.788 or 1.554, much larger than βACV _ASA ’s 1.349. Table 2 also shows their classification performance on the 396 mutations: FoldX’s F1 is 0.44, while Robetta’s F1 is 0.494, both worse than βACV _ASA ’s 0.604. From Table 2, our method also has better classification performance than HotPOINT and KFC. Other performance comparison results are provided in the Additional file 1 when tested on BID (including protein-peptide interfacial residues) or under the leave-one-complex-out cross-validation framework.

An example of hot spot predictions

We use 3HFM as a case study to illustrate the difference of the binding hot spot prediction results by βACV _ASA , FoldX and Robetta. The 3HFM complex is an antibody-antigen binding between HyHEL-10 and hen egg white lysozyme. According to ASEdb, a total of 25 alanine mutations were experimented, 11 of which have Δ Δ G more than 2 kcal/mol.

Our βACV _ASA correctly identified 9 binding hot spot residues with a recall of 0.818, but made 3 false positive predictions with a precision of 0.75 (Figure 3). This gives an F1 of 0.783. However, FoldX made only one hot spot prediction which is correct with a recall of 0.091, and Robetta has a recall of 0.455 (5 out of 11) and a precision of 0.833 (5 out of 6), namely an F1 of 0.588. Both of these methods have a lower prediction performance than our βACV _ASA . What is more important is that the four positive predictions correctly made only by our βACV _ASA , not by FoldX or Robetta, have a high Δ Δ G, such as Trp in position 98 of Chain H (Δ Δ G=5.5 kcal/mol) and Tyr in position 50 of Chain L (Δ Δ G=4.6 kcal/mol). Please refer to Additional file 1: Table S3 for detail.

Figure 3

Prediction results by β ACV _ASA for the residues in the interface between Chain Y and Chain HL(together) in 3HFM. In (a) and (b), the true predicted hot spot residues are in magenta, the false predicted non-hot spot residues are in yellow, the false predicted hot spot residues are in orange, while the true predicted non-hot spots are in cyan; X-YZZ stands for residue Y in position ZZ of Chain X.

Conclusion

A new classifier βACV _ASA has been proposed to predict Δ Δ G and binding hot spot residues. The novel idea of this classifier is to integrate the water exclusion theory into β contacts. Tested on a data set of 396 alanine mutations, βACV _ASA has been found to outperform many other methods. This confirms that β contacts are truly better than traditional distance-cutoff contacts to capture the energetic characteristics of protein folding and binding. This also confirms that water exclusion is a necessary condition for a residue to become a binding hot spot residue.

Availability of supporting data

All the supporting data are included as additional files.