Protein-protein binding site identification by enumerating the configurations

Abstract

Background

The ability to predict protein-protein binding sites has a wide range of applications, including signal transduction studies, de novo drug design, structure identification and comparison of functional sites. The interface in a complex involves two structurally matched protein subunits, and the binding sites can be predicted by identifying structural matches at protein surfaces.

Results

We propose a method which enumerates “all” the configurations (or poses) between two proteins (3D coordinates of the two subunits in a complex) and evaluates each configuration by the interaction between its components using the Atomic Contact Energy function. The enumeration is achieved efficiently by exploring a set of rigid transformations. Our approach incorporates a surface identification technique and a method for avoiding clashes of two subunits when computing rigid transformations. When the optimal transformations according to the Atomic Contact Energy function are identified, the corresponding binding sites are given as predictions. Our results show that this approach consistently performs better than other methods in binding site identification.

Conclusions

Our method achieved a success rate higher than other methods, with the prediction quality improved in terms of both accuracy and coverage. Moreover, our method is being able to predict the configurations of two binding proteins, where most of other methods predict only the binding sites. The software package is available at http://sites.google.com/site/guofeics/dobi for non-commercial use.

Background

Most of the existing efforts to identify the binding sites in protein-protein interaction are based on analyzing the differences between interface residues and non-interface residues, often through the use of machine learning or statistical methods. These methods differ in the features analyzed, that is, the sequence and structural or physical attributes. Chung et al.[1] used multiple structure alignments of the individual components in known complexes to derive structurally conserved residues. Sequence profile and accessible surface area information are combined with the conservation score to predict protein-protein binding sites by using a Support Vector Machine. Ofran et al.[2] employed neural networks to predict binding sites, using the sequence environment, the profile and the structural features as input. The random forest algorithm is used to utilize these features from sequences or 3D structures for the binding site prediction [3, 4]. PSIVER [5] uses sequence features for training a Naïve Bayes classifier to predict binding sites. In PSIVER, conditional probabilities of each sequence feature are estimated using a kernel density estimation method.

Besides the machine learning and statistical approaches, 3D structural algorithms and other methods have also been used to identify binding sites through investigating protein surface structures. ProBiS [6] predicts binding sites by local surface structure alignment. It compares the query protein to 3D protein structures in a database to detect proteins with structurally similar sites on the surfaces. Burgoyne et al.[7] analyzed clefts in protein surfaces that are likely to correspond to the binding sites. They ranked them according to sequence conservation and simple measures of physical properties including hydrophobicity, desolvation, electrostatic and van der Waals potentials. Ortuso et al.[8] defined most relevant interaction areas in complexes deriving pharmacophore models from 3D structure information. It is based on 3D maps computed by the GRID program on structurally known molecular complexes.

ProMate [9] is based on the idea of interface and non-interface circles. A circle is first created around each residue. Then, features are extracted from these circles. Statistics are performed and histograms are created for each feature. Thereafter, the probability for each circle of a test protein to be an interface is estimated. The interface circles are clustered for each test protein to identify the binding patch.

Bradford et al.[10] proposed an approach (PPI-Pred) which uses SVM (Support Vector Machine) on surface patch features to predict binding sites. PPI-Pred generates an interacting patch and a non-interacting patch for each protein. Seven features are extracted for each patch to build an SVM model, which is then used to predict if a given test patch is an interacting patch.

In PINUP [11], an empirical scoring function is presented to predict binding sites. The function is a linear combination of energy score, interface propensity and residue conservation score. A patch is formed by a residue and its spatial neighbors within the protein subunit. PINUP takes the top 5% scoring patches and ranks residues based on their occurrences in these patches. The top 15 ranked residues are predicted as the interface residues.

Li et al.[12] proposed another SVM approach (core-SVM). The residues of the proteins are divided into four classes: the interior residues, the core interface residues, the rim interface residues, and the non-interface residues. The core interface and rim interface residues are distinguished by the percentage of their neighboring residues which are interface residues. An SVM is built over eight features extracted from the interface residues, and used to compute the probability of whether a residue is a core interface residue.

Meta-servers have also been constructed to combine the strengths of existing approaches. The program called meta-PPISP [13] combines three individual servers, namely cons-PPISP, ProMate and PINUP; another program called metaPPI [14] combines five prediction methods, namely PPI-Pred, PINUP, PPISP, ProMate, and SPPIDER [15].

Another approach in binding site prediction is to examine the possible structural configurations, or referred to as poses, of protein subunits, that is, how the subunits may dock. Docking methods based on fast Fourier transformation (FFT) [16, 17], geometric surface matching [18], as well as intermolecular energy [19–21] have been proposed. Fernández-Recio et al.[22] simulated protein docking and analyzed the interaction energy landscapes. Their method uses a global docking method based on multi-start global energy optimization of the ligand. It explores the conformational space around the whole receptor, and uses the rigid-body docking configurations to project the docking energy landscapes onto the surfaces. The low-energy regions are predicted as the binding sites.

In this paper, we propose a method which enumerates the configurations of two binding proteins (that is, the possible positions of the two subunits in a complex), and identify binding sites by evaluating the interaction between the components using the Atomic Contact Energy (ACE) function [23]. We perform rigid transformation to enumerate the configurations of two binding proteins. The enumeration is performed in conjunction with a surface identification technique for avoiding clashes between protein subunits when computing rigid transformations. The transformations which result in the minimum score according to the Atomic Contact Energy function are found; the corresponding interacting residues are reported as binding sites. Our method is implemented in a program called DoBi^a.

We perform experiment to compare DoBi with the existing methods using commonly used measures for assessments. The program outperforms the other methods on these measures. DoBi achieved a success rate higher than all the other methods, improving prediction quality in terms of both accuracy and coverage. In addition, it predicts the configurations of two binding proteins, as opposed to giving only the binding sites.

Methods

The main idea of our method is to enumerate “all” configurations between two proteins, where a configuration refers to the 3D coordinates representing the relative position and orientation of two protein subunits in a complex. We use the Atomic Contact Energy (ACE) function to compute the score for a configuration. The configurations with the lowest score are chosen, and the corresponding interacting residues are predicted as binding sites. We use rigid transformation to enumerate the configurations. The key techniques required here contain (1) an efficient algorithm to enumerate “all” configurations (rigid transformations) and (2) a good energy score.

Atomic contact energy

Atomic Contact Energy (ACE) is an atomic desolvation energy measure developed in [24]. It is defined over the energy of replacing a protein-atom/water contact, with a protein-atom/protein-atom contact. The ACE score takes into account 18 atom types, hence resulting in 18×18 possible atom pairs. The score for each atom pair has been determined, based on a statistical analysis of atom-pairing frequencies in known proteins. These pre-determined scores are given as log likelihood values in [24], thus allowing the summation of these values. The pre-determined score of effective contact energy between atom type i and type j is defined as

$$T[i,j]=-ln\frac{N_{i,j}/C_{i,j}}{(N_{i,0}/C_{i,0})\times (N_{j,0}/C_{j,0})}$$

(1)

where type 0 corresponds to the solvent. The number of i-j contact (N_i,j) and the number of i-0 contact (N_j,0) are estimates of the actual contact numbers of known complexes. In addition, C_i,j and C_i,0 are defined as the expected numbers of i j contact and i-0 contact.

For a given configuration, the ACE score is a summation of each of the atom pairs (one from each subunit) within threshold distance d, and d = 6Å is used in this paper. Denote the sets of atoms from the two subunits as S₁ and S₂, respectively, then the ACE is computed as

$$E_{\mathit{\text{ACE}}}=\sum _{s\in S_1,t\in S_2,\left|\right|s-t\left|\right|\le d}T[s,t]$$

(2)

where |s−t| is the Euclidean distance between s and t, and T[s,t] is the pre-determined score of the atom pair s and t.

The ACE score can be considered an estimate of the change in desolvation energy of the two proteins in going from the unbound state to the complex. A lower ACE value implies a lower (and hence more favorable) desolvation free energy.

Enumeration of the configurations

In this paper, we assume that subunits are rigid. A protein structure consists of a sequence of residues. Each residue consists of a set of atoms. We assume that the atoms in a residue are ordered as a sequence. Hence, the whole protein structure can be represented by a sequence of atoms. In the rest of this subsection, we let A and B denote two protein structures (subunit), and write A = (a₁,a₂,…,b _m ), and B = (b₁,b₂,…,b _n ), where a _i , and b _j are atoms of structure A and B. Without loss of generality, we assume that n ≥ m. We also assume that we know the 3D coordinates of each atom in both input proteins. We use A[i:j] to denote the subsequence (a _i ,…,a _j ), and refer to a subsequence of atoms as a structural fragment.

To enumerate all the configurations, we assume B is fixed, and we perform rotations and translations (referred to as rigid transformations, and simply, transformations, in the rest of the paper) on A. The method proposed here is modified from the algorithms for structure comparison [25].

Assume that two points a _i and a _j of A interact with two points b _i′ and b _j′ of B, then we know that ||a _i − b _i′ || ≤ d and ||a _j − b _j ′ || ≤ d. To enumerate the configurations, we enumerate the positions for atoms a _i and a _j first, and for each fixed positions of a _i and a _j , we rotate A about the line formed by a _i and a _j . Let the d-ball of an atom a be the ball with radius d centered at a. We discretize the d-ball of b _i′ with step size εd, where ε is a small constant (and we choose ε = 0.1 for this paper). Each grid point in the d-ball of b _i′ is used as a candidate position for atom a _i for the binding. When a _i is fixed at one of the grid points, the possible positions for a _j form a sphere cap, where the sphere is centered at a _i with radius |a _i −a _j |, and the cap is the portion of the spheres enclosed in the d-ball of b _j′ . Again, we discretize the sphere cap with step size εd. Each grid point on the sphere cap is a candidate position for a _j . This gives us a total of $O\left({\left(\frac{1}{\epsilon }\right)}^5\right)$ possible positions for the pair of a _i and a _j . After a _i and a _j are fixed on their respective grid points, the only degree of freedom to move A[i,j] is to rotate it around the axis through a _i and a _j . We use a 1° step size; that is, we explore 360 different positions for the remaining atoms through 360 rotations. Figure 1 illustrates the steps to compute a transformation.

Figure 1

Steps to obtain a transformation. (1) put a _i at one of the $O\left({\left(\frac{1}{\epsilon }\right)}^3\right)$ grid points d-ball of b _i′ . (2) put a _j at a grid point on the intersection of the sphere centered at a _i with radius |a _i a _j | and d-ball of b _j′ . There are at most $O\left({\left(\frac{1}{\epsilon }\right)}^2\right)$ grid points on the intersection. (3) use a _i and a _j as the rotation axis.

The method will work well if we know two interaction pairs (a _i ,b _i′ ) and (a _j ,b _j′ ). We can simply enumerate all the atoms pairs as the interaction pair candidate. However, there will be O(n⁴) such cases, which makes the computer program too slow in practice. This is perhaps one of the reasons that such a method has not been tried. The focus of the following subsection is to identify two pairs (a _i ,b _i′ ) and (a _j ,b _j′ ) which are more likely to be interaction pairs.

When enumerating “all” configurations, we also want to make sure that (1) only surface fragments can be candidate binding sites for a configuration and (2) there is no clash between the two proteins in such a configuration. Before presenting the details of the method, we define the surface atoms and clashes of two subunits first.

Surface atoms

The interface residues of two proteins are necessarily surface residues. Inspired by the work in LIGSITE _csc [26, 27], we propose a method to identify the surface atoms of a protein.

First, we build a 3D grid with step size 1Å around the protein. Then, each grid point is labeled as a protein point if it is within distance 2Å of any atom, and labeled as empty otherwise. We further subdivide the protein grid points into two types: interior or surface. A protein grid point is labeled as surface if at least one of its six neighboring grid points is empty, otherwise it is labeled as interior. With the grid points labeled, we can label the atoms. an atom is labeled as a surface atom if it is within distance 1.5Å of a surface grid point, otherwise it is labeled as an interior atom.

Figure 2 gives an example in 2D, where a protein grid point is labeled as interior if it has all four neighbors as protein points. In 3D, a protein grid point should be labeled as interior if all of its six neighbors are labeled as protein.

Figure 2

The surface atoms are indicated in 2D. (A) the grid is created, and grid points are labeled as either empty or protein; (B) the grid points labeled as protein are relabeled as surface or interior; (C) an atom is labeled either as surface or as interior. We use 2D as an illustration.

Clashes of two subunits

A configuration cannot result in two subunits to have clashes. The following method is used to capture if a configuration resulted in clashes. Given a configuration, we build a 3D grid as in the previous subsection. For each of the structures A and B, we mark the grid points as interior, surface, or empty. We use a threshold θ to identify whether two subunits clash, by calculating the proportion of interior points for both of them. We say that the two subunits clash if they share more than θ × 100% of their interior points; that is, if X is the number of interior grid points which are shared by both proteins, and X _A and X _B are the number of interior grid points of each subunit, respectively, then we require that X ≤ θ × min{X _A ,X _B } if the subunits do not clash.

Finding the two interaction pairs

In the following subsections, we present the details to explore the potential interaction pairs.

Identify candidate fragment pairs

We first select fragment pairs that are potential binding sites. As discussed in Section “Enumeration of the configurations”, there are O(n⁴) possible fragment pairs (a _i , a_i′) and (b _j , b_j′) for each binding site. To reduce the computational complexity, we adopt a local alignment algorithm to accelerate this selection. This is a raw estimation and we hope that the actual binding sites are not discarded by this process.

We first use a heuristic to quickly discard fragments pairs that are unlikely to bind. The heuristic simplifies the problem, as follows: (1) every atom is within the threshold value required in the ACE computation (that is, we ignore the geometry of the structure); (2) each atom interacts with at most one atom; (3) interacting pairs follows a sequential order. That is, for any two pairs of interacted atoms (a _i , b _i′ ) and (a _j , b _j′ ), we have either i < i^′and j < j^′, or i^′< i and j^′< j. With these three simplifications, the standard Smith-Waterman local alignment algorithm [28] can be employed, with the ACE scores used as the penalty (negation of the score) for alignment. We use a penalty of 1 for aligning an atom to a space. Each local aligned segment gives us two fragments, where each atom in the fragment is either aligned to another atom from the partner, or aligned to nothing (i.e., aligned to space).

We present details here. For two sequences P₁and P₂, an alignment of P₁ and P₂ can be obtained by (1) inserting spaces into the two sequences P₁ and P₂ such that the two resulting sequences with inserted spaces P′₁ and P′₂ have the same length and (2) overlap the two resulting sequences P′₁ and P′₂. The score of the alignment is the sum of the scores for all the columns, where each column has a pair of letters (including spaces) and for each pair of letters there is a pre-defined score. A subsequence α of P₁ and a subsequence β of P₂ can be formed as a local aligned segment such that the score between α and β is minimum. Here we want to find all (non-overlapping) pairs of subsequences with a score of at most x. For our purpose, we set x = 0 throughout the paper.

Due to the simplifications, there are many false positive results, and some of the interaction pairs can be filtered. The latter issue can be handled to some extend by raising the threshold. The former issue is tackled by further refinement in the next subsection. In practice, our program outputs 70 to 120 fragment pairs as potential binding sites, which is much smaller than O(n⁴), where the number of atoms n in a protein is from 500 to a few thousands.

Since a binding site is necessarily on the surface of a subunit, we filter out fragments with only very few atoms on the surface. To achieve this, we use a sliding window of length 15 to parse the aligned fragment pair. For each window, if the surface atoms are at least 2/3 (that is, ten atoms) for both fragments, the fragment pair of this window is kept for further processing and this fragment pair is extracted from the alignment. We continue this process on the un-extracted portion of the alignment. If the window does not contain sufficient surface atoms, we continue at the next window. Our choice of 2/3 comes from observations with a docking decoy set from the Dockground [29], where 94% of the binding sites have more than 2/3 of surface atoms.

Identify configurations of fragment pairs

From the fragment pairs obtained in the previous step, a second step is used to further filter out fragment pairs of ACE scores below a threshold. Given two structural fragments A[i,j] = (a _i ,…,a _j ), and B[i′,j′] = (b i′,…,b j′), we assume that a _i interacts with b _i′ , and a _j interacts with b _j′ . Using the enumeration method described earlier, we enumerate different configurations for A and B and compute the corresponding ACE score for the atom sets A[i,j] and B[i^′,j^′]. We do not consider any configuration which causes A and B to clash. In this step, a pair of structural fragment which does not give any configuration with an ACE score below a specified threshold is discarded. In this paper, we define the threshold value as 400, since the ACE scores of actual interface in the docking decoy set from Dockground are all less than 400. After this step, it is unlikely for two protein structures which cannot be bound to have an unfiltered fragment pair.

Identify the configuration for the two subunits

In the third step, for each pair of protein structures with at least one remaining fragment pair, we enumerate all the potential configurations for the structures. We want to use the begin and end atoms of the identified fragments for our choice of (a _i , b _i′ ) and (a _j , b _j′ ) in the enumeration, since these are the atoms that are likely to be interacting. Assuming that there are k fragment pairs from the same two proteins left after the filtration of the second step, we will have a maximum of 2k distinct atom pairs to choose. Thus, there is a total of at most $\left(\genfrac{}{}{0ex}{}{2k}{2}\right)$ combinations to consider for the choice of (a _i , b _i′ ) and (a _j , b _j′ ).

When the best configuration is obtained, two residues, one from each subunit, are reported as the interface residues if they can be connected with a pair of atoms within distance 4.5Å. In our search for the best configuration, we also require the configurations to be free from clashes.

Results and discussion

Three commonly used measures are utilized to assess the performance of DoBi. Accuracy and Coverage are two common measures to assess the quality of the binding sites adopted by a method [11]. The accuracy of the predicted interface is the fraction of correctly predicted residues over the total number of predicted interface residues; the coverage of the predicted interface is the fraction of correctly predicted interface residues over the total number of actual interface residues. F-score ($F=2\times \frac{\text{Accuracy}\times \text{Coverage}}{\text{Accuracy}+\text{Coverage}}$) is a weighted average of the accuracy and coverage, where an F-score reaches its best score at 1 and worst score at 0. Another common measure is success rate, which is defined in [9]. A reported result is claimed as a success if at least half of the predicted residues are actual interface residues; that is, the accuracy is no less than 50%. The success rate is the fraction of successful predicted cases in the total number of predicted proteins.

A protein complex may contain several subunits, and multiple binding sites. Each binding site in a protein complex consists of a pair of subunits. Two residues in a pair of subunits are called interface residues if any two atoms, one from each residue, interact. By interact, we mean the distance between the two atoms is less than the sum of the van der Waals radius of the two atoms plus 1Å. The number of residues on interface is referred to as the interface size.

Training set

We use the unbound protein structures from Dockground [29] as the training set to calculate the parameters of DoBi. The docking decoys from Dockground were generated by GRAMM-X scan. The GRAMM-X docking scan was used to generate 102 unbound-unbound complexes and 131 unbound-bound complexes. By excluding the proteins used in the comparison, 36 unbound-unbound complexes and 80 unbound-bound complexes can be used to calculate the value of the threshold θ. When we set θ = 0.17, the overall F-score of DoBi on the training set is 60.5%, which is the best score that DoBi achieves under different threshold values. The details on the training set are shown in Table 1.

Table 1 Details of DoBi on the training set

Comparison to the existing methods

We divide our comparisons into four separate groups, where in each group we compare a different set of methods. The reason that we cannot compare all the methods with the same data set is due to the unavailability of some methods, in which case the only comparison possible is with the results in the respective publications.

Comparison to Fernández-Recio et al.’s method

DoBi is compared to the method introduced by Fernández-Recio et al. in [22], using the test data therein, which consists of 43 complexes. The results are reported in Table 2. The overall accuracy and coverage for DoBi are 44.3% and 70.5%. Fernández-Recio et al. ’s method achieved the overall accuracy and coverage of 39.3% and 72.7%, respectively. The success rate for DoBi is 39.6%, improving over the success rate of 37.2% reported by Fernández-Recio et al.. The F-score is 0.54 for DoBi, and 0.51 for Fernández-Recio et al.’s method.

Table 2 Comparison of DoBi and Fernández-Recio et al. ’s method

The average predicted sizes for DoBi and Fernández-Recio et al.’s method are 37.5 residues and 46.3 residues respectively, while the average actual size is 21.1 residues. The standard deviation of the sizes predicted by DoBi is 29.0, while that of the sizes predicted by Fernández-Recio et al.’s method is 40.0.

Table 3 displays the detailed results for all unbound structures of 43 complexes. Each row corresponds to a pair of proteins. We can observe from the table that the binding sites are identified accurately for the complexes 2sni(E:I), 2sic(E:I), 1ay7(A:B) and 1wq1(G:R).

Table 3 Detailed Results of DoBi and Fernández-Recio et al. ’s method

Comparison to metaPPI, meta-PPISP and PPI-Pred

In this group of our comparisons, the test set in [14] is used. It consists of 41 complexes from the benchmark v2.0 [30] and 27 targets from the CAPRI experiment[31]. The 41 complexes are divided into two categories, enzyme-inhibitor (EI) and others. We compare our method to metaPPI, meta-PPISP and PPI-Pred with this group of data. The overall accuracy and coverage of each prediction method are shown in Table 4. DoBi has an F-score of 0.55, where in contrast, metaPPI, meta-PPISP and PPI-Pred have the F-scores 0.35, 0.43 and 0.32 respectively. DoBi has a success rate of 53.7%, as well as overall accuracy and coverage of 50.0% and 60.0% respectively.

Table 4 Comparisons of DoBi, metaPPI, meta-PPISP and PPI-Pred

The detailed results on all the unbound structures of the 41 complexes are displayed in Table 5. The detailed results on 27 CAPRI targets are displayed in Table 6. Each row displays the results of the methods tested on the two corresponding binding partners.

Table 5 Detailed Results of DoBi, metaPPI, meta-PPISP and PPI-Pred on 41 complexes

Table 6 Detailed Results of DoBi, metaPPI, meta-PPISP and PPI-Pred on 27 targets

Besides the identification of binding sites, our program also estimates the orientations and positions of the proteins after binding. Figure 3 displays the orientation and position discovered by our program for 1qa9(A:B). The C _α interface RMSD (root mean squared deviation) (iRMSD) between the experimental structure and the predicted complex is 2.36Å.

Figure 3

Configuration discovered by DoBi for 1qa9(A:B). (A) is the figuration by DoBi; and (B) is the experimental structure. The C _α iRMSD between two complexes is 2.36Å.

Comparison to ProMate and PINUP

In this experiment, DoBi is compared to ProMate and PINUP. The test data is originally used by ProMate, and consists of 57 non-homologous proteins. The results are reported in Table 7. DoBi has an F-score of 0.56, while PINUP and ProMate have the F-scores 0.43 and 0.21 respectively. The overall accuracy and coverage of DoBi are 54.2% and 59.1%. The success rate of DoBI is 64.9%. Hence the success rate is improved by at least 1.8%, while the overall accuracy and coverage are improved by at least 1.7% and 16.6% respectively.

Table 7 Comparison to PINUP and ProMate

The average of the sizes predicted by DoBi, PINUP and ProMate are 23.5 residues, 19.0 residues and 5.4 residues respectively, while the actual average size (average size of actual interface residues) is 21.0 residues. The number of residues correctly predicted to be on interface by DoBi, PINUP and ProMate are 12.3 residues, 8.3 residues and 2.7 residues respectively.

Table 8 shows the detailed results of 57 unbound proteins. DoBi performed better for most of the cases. However, for some cases where all three methods do not perform well, DoBi is usually the worst, e.g. 1avu_, 1aye_, 1qqrA and 1b1eA.

Table 8 Detailed Comparison to PINUP and ProMate

Comparison to core-SVM

In this study, we compare DoBi to core-SVM using the same data set of 50 dimers which core-SVM was tested against [12]. The results are shown in Table 9. The overall accuracy and coverage for our method are 59.0% and 61.1%, while those for core-SVM are 53.4% and 60.6%. The success rate of DoBi is 70.0% on 50 pairs of proteins in those binary complexes. The F-score is 0.60 for DoBi, and 0.56 for core-SVM. The average of the size predicted by DoBi is 39.0 residues (with standard deviation 19.1), while the average actual size is 40.3 residues. The number of residues correctly predicted by DoBi to be on the interface is 22.5.

Table 9 Comparison to core-SVM

Table 10 shows the details for DoBi on the data set used by core-SVM. The performance of DoBi is particularly good on several proteins such as 1aym2 and 1rzhM.

Table 10 Detailed Results for DoBi on the data set used by core-SVM

Evaluation on benchmark v4.0

To further evaluate our method, we perform tests on the protein-protein docking benchmark v4.0 [32, 33]. This benchmark consists of 176 complexes. Proteins dynamically change their conformations upon binding with other proteins [34]. A single protein without binding with any other structure is referred to as unbound, whereas a protein with a binding partner in a complex is referred to as bound. We test our method in both the bound and the unbound cases.

" Running time "

We used a Pentium(R) 4 (CPU of 3.40GHz) to run DoBi. The computation for each of the 176 complexes took 100 seconds on average.

Results on bound states

The complexes are classified into broad biochemical categories: Enzyme-Inhibitor (52), Antibody-Antigen (25) and Others (99). The average accuracy and coverage of DoBi are 61.8% and 67.9% respectively on the 52 complexes in Enzyme-Inhibitor, 51.6% and 70.1% on the 25 complexes in Antibody-Antigen, and 58.2% and 69.1% on the 99 complexes in Others. A success rate of 77.6% is achieved for the Enzyme-Inhibitor complexes. The details are shown in Table 11.

Table 11 DoBi’s performance for proteins of benchmark v4.0 in bound states

Results on unbound states

The pairs of unbound proteins are classified into three categories: 121 rigid-body (easy) cases, 30 medium difficult cases, and 25 difficult cases, according to the magnitude of conformational change after binding [30]. The average accuracy and coverage of DoBi are 43.6% and 65.4% on the 121 rigid-body cases, 34.1% and 56.7% on the 30 medium difficult cases, and 32.4% and 53.4% on the 25 difficult cases. The success rate of DoBi is 41.7% for the rigid-body cases, which is significantly better than for the other two categories. In general, the accuracy and coverage decrease as the magnitude of conformational increases. The details are shown in Table 12.

Table 12 DoBi’s performance for proteins of benchmark v4.0 in unbound states

DoBi discovered several good configurations for the medium difficult cases. One of the instances is 1wq1(R:G). Its configuration discovered by DoBi is shown in Figure 4. The C _α iRMSD between the experimental structure and the predicted complex is 4.12Å.

Figure 4

Configuration discovered by DoBi for 1wq1(R:G). (A) is the configuration by DoBi; and (B) is the experimental structure. The C _α iRMSD between two complexes is 4.12Å.

Docking result of DoBi

DoBi is optimized for binding site prediction, but it also can be used to dock two protein structures. We compare DoBi’s poses to the best configurations obtained by ZDOCK and 3D-Dock. ZDOCK [35] uses a fast Fourier transform to search all possible binding modes for the proteins, and evaluates them based on shape complementarity, desolvation energy, and electrostatics. It can produce structures with the smallest iRMSD values in top 1000 predictions with minimum energy. 3D-Dock[36, 37] uses an initial grid-based shape complementarity search to produce lots of potential interacting conformations. They can be ranked by using interface residue propensities and interaction energies. It reports structures with the smallest iRMSD values in top ten predictions.

Table 13 The Docking Results of DoBi, ZDOCK and 3D-Dock on CAPRI

We calculate the predicted structures by different methods on the complexes in benchmark v4.0 and the targets in CAPRI. CAPRI is a community-wide experiment to assess the capacity of protein docking methods to predict protein-protein interactions [31]. The C _α iRMSD, F-score and the fraction of native contacts are used to evaluate the results by different methods. The fraction of native contacts is used by 3D-Dock[37]. It is calculated as the total number of native contacts for the predicted configuration divided by the total number of contacts in the native structure. A native contact exists between residues i and j if distances between them in native structure and in predicted configuration are both less than 4.5Å.

We compare the docking results of DoBi, ZDOCK and 3D-DOCK on the CAPRI targets. The results are shown in Table 14. The top 1,000 configurations predicted by DoBi and ZDOCK are used for comparison. Among the 1,000 predictions, we choose the configuration of the best iRMSD value to evaluate the methods. The average iRMSD values for DoBi and ZDOCK are 7.5Å and 6.9Å, respectively. However, the average fractions of native contacts for DoBi and ZDOCK are 40.6% and 35.2%, respectively. DoBi improves the F-score of binding site prediction by at least 1.3%. DoBi’s performance on docking is worse than ZDOCK, but its performance on binding site prediction is more accurate than ZDOCK.

Table 14 The Docking Results of DoBi and ZDOCK on Benchmark v4.0

Each of DoBi and 3D-Dock produced ten results for each target, and the configurations with smallest iRMSD values among those ten predictions are used for comparison. The average iRMSD values for DoBi and 3D-Dock are 9.2Å and 9.1Å. However, the overall fractions of native contacts for DoBi and 3D-Dock are 29.1% and 26.8%. DoBi’s performance on binding site prediction is better than that of 3D-Dock.

The docking results obtained by DoBi and ZDOCK on Benchmark v4.0 are shown in Table 15. Similarly, we compare the best configurations in the top 1000 predictions from each method of DoBi and ZDOCK for each target. The average iRMSD values of DoBi and ZDOCK are 6.1Å and 4.9Å, respectively. For the binding site prediction, the overall F-score values of ligand proteins by DoBi and ZDOCK are 69.5% and 69.4%, and those of receptor proteins by DoBi and ZDOCK are 68.2% and 66.1%, respectively. These results indicate that DoBi’s performance on binding site prediction is better than ZDOCK. The docking quality of DoBi requires further efforts to improve.

Table 15 Comparison of Atomic Contact Energy for the Predicted Complexes and the Experimental Structures on Benchmark v4.0

We calculate the docking results of 1i4d. The C _α iRMSD values between the experimental structure and the configurations by DoBi and ZDOCK are 2.97Å and 1.97Å, respectively. DoBi improves F-score value of ligand protein by 2.7%, and that of receptor protein by 0.4%. The configurations produced by methods are shown in Figure 5.

Figure 5

Configuration discovered by DoBi and ZDOCK for 1i4d. (A) is the configuration by DoBi; (B) is the configuration by ZDOCK; (C) is the experimental structure.

Factors affecting the performance of DoBi

We notice that DoBi performed badly on a few specific instances. We analyze this performance issue with Table 13, which compares the ACE scores for the experimental structures and predicted complexes, for the bound states of proteins in the benchmark v4.0. Among the 176 complexes, only 43 of them have an ACE score for experimental structures lower than that of the predicted complexes. This implies that in 133 cases, DoBi is able to find a configuration of a lower score than the experimental structures. These anomalies suggest that the score function currently used in DoBi may be inaccurate, and this inaccuracy may have contributed to the poorly performed cases of DoBi. We also note that the search space currently explored by our method is incomplete, and this may have contributed as well to the inaccuracy of DoBi in some cases.

Figure 6 shows the protein complex incorrectly predicted by DoBi as well as the experimental structure for 1kxq(H:A). The iRMSD between the two complexes is 18.87Å. The ACE score of the docking structure predicted by DoBi, -497.6, is lower than the ACE score of the experimental structure, 63.7. The binding sites predicted by DoBi are incorrect as well.

Figure 6

DoBi fails to solve the instance 1kxq(H:A). (A) is the predicted complex; and (B) is the experimental structure.

Conclusions

In this work, we proposed an approach to identify binding sites in protein complexes by docking protein subunits. The method is implemented in a program called DoBi. DoBi consistently and significantly performed better than existing techniques in predicting binding sites in experimental results.

We identify a few potential areas for future improvements to our method. The first area to work on is in the energy function used. Currently, DoBi uses a simple score function. As suggested by the experiment results, a better energy function is able to improve the performance of DoBi.

A second area for improvement is in our current assumption that protein structures are rigid when binding. In reality, protein structures may vary sightly or even dramatically when they bind. Hence, further studies on this issue are very much in demand.

Although our method shows better overall performance, there are some protein complexes where other methods outperformed DoBi. It will be beneficial if we could combine the strengths of these existing programs with DoBi, to come up with a more reliable method.

Endnote

^aThe initial two letters from each of the two words, Docking and Binding, were taken.