Received June 11, 2007

Introduction

A protein–ligand docking program has been used to efficiently discover lead compounds for a target protein from a huge compound database. Since the pioneering work by Kuntz et al.,¹ numerous docking programs have been developed.^2–15 There are two standpoints for assessing docking programs: conformation search or scoring. For the former, it is assessed whether the pose or binding mode of a crystal structure can be reproduced. For the latter, it is assessed whether the predicted score is correlated with experimentally measured binding affinity. Many reports assessing the performance of a docking program have been published.^16–29 Of particular interest is the recent critical assessment reported by Warren et al.²¹ These reports generally reached a conclusion that docking algorithms are highly successful at generating good binding modes, while scoring functions are less successful at correctly identifying a binding mode. More specifically, some of the docking programs can generate a conformation very close to the crystal binding structure, but one cannot know which conformation created by which program is close to the true structure. Under such situation, it is desirable to further improve scoring method to rerank conformations generated. To supplement the weakness of a single scoring function, a consensus scoring (CS) method combining multiple scoring functions was proposed.^22–30 Although CS typically does not perform better than the best of the scoring functions combined, it has the advantage of generally providing stable results. Since which scoring function is the best can only be known after the analysis of results, not before, CS is a reliable method for a blind trial. It has been demonstrated^22–30 that hit rate or false positive is improved by CS.

In this study, we developed a bootstrap-based consensus scoring (BBCS) method. BBCS combines multiple scorings, each of which has the same function form but different energy-parameter set. BBCS is derived from the idea of bootstrap aggregating (bagging),³¹ which is one kind of ensemble learning. In general, it is not easy to uniquely determine the energy-parameter set that provides the global minimum of binding free energy over all complexes when the number of data in a training set is not enough. To circumvent this problem, BBCS based on ensemble learning is a rational method.

Multiple energy-parameter sets are determined in two steps: (1) generation of training sets by the bootstrap method and (2) optimization of an energy-parameter set by a Z-score approach against each training set. The Z-score approach is based on energy landscape theory, which leads to describes a funnel-shaped landscape of protein folding³² and has been applied to folding simulation, structure prediction^33–38 and de novo design³⁹ of proteins and protein-peptide docking.⁴⁰ We consider that protein–ligand docking also has a funnel-shaped landscape, as discussed by Camacho and Vajda,⁴¹ and the Z-score approach works well for protein–ligand docking.

Machine leaning has been applied in several ways to overcome difficulties in structure-based consensus scoring.^{27, 42, 43} For example, Jacobsson et al. used the ensemble method (bagging) to create rule-based prediction models based on the scoring matrices.²⁷ Another example is in a recent paper of Antes et al., who applied an ensemble method to neural networks and k-nearest neighbor models.⁴² Although the aim of theirs and our work is to find a more plausible binding mode among numerous modes, there are marked differences. They used an ensemble method to optimize a parameter set efficiently, while we use it to improve the rank of the crystal structure relative to decoy structures. In other words, they determined a single parameter set used for a scoring function, not multiple parameter sets as ours does. Their parameter set is optimized so that the binding mode with the lowest rmsd becomes close to the crystal structure, while our parameter set is optimized so that crystal structure is maximally separated from decoys with respect to score.

Methods

Preparation of training and test sets. In this study, 98 protein–ligand complexes were used. These complexes are classified into three types²⁶ according to the dominant interactions between the protein and ligand: electrostatic interaction for type 1, hydrophobic for type 2, and mixed electrostatic and hydrophobic for type 3. The complexes in each type were randomly divided into “training set” and “test set,” as listed in Table 1. For each complex, the crystal structure and 100 docked structures (decoys) generated by Wang et al.²⁶ were used. Their docked structures are carefully generated in conformational space of ligand as completely as possible so that their distribution is not biased by the docking algorithm employed. Such structures allow various scoring functions to be compared on the same ground. Wang et al. used AutoDock,⁷ because it can generate a wide variety of conformations during long multiple stochastic searches on potential energy surface. Wang’s data set is useful for testing the performance of scoring independently of the performance of docking algorithm. Our intention of using Wang’s data set is to proceed to refining scoring functions beyond comparing their performances. To show an example of clear advantage of Wang’s data set over other data sets such as the data set generated by FlexX program with default protocol, we compared the distributions of conformations in the two data sets, in Supporting Information S-1. There is seen a remarkable difference for some of the complexes; FlexX fails to generate conformations sufficiently thoroughly in a low range of rmsd for the complex shown. If such insufficient data set of conformations is used for test of scoring, the performance of scoring should be coupled with the performance of docking algorithm. On the other hand, Wang’s data set is more complete, containing conformations searched in wide range. Therefore, we consider that Wang’s data set is useful as benchmark set for scoring.

Wang’s data set is published online at http://sw16.im.med.umich.edu/software/xtool. In the published structures, bound water molecules are eliminated but metal atoms are included. In this study, metal atoms are also eliminated.

Scoring Function. In BBCS, any scoring function proposed can be employed. We employ the FlexX scoring function, which is given as²⁴

For convenience in optimizing energy-parameter, we rewrite eq 1 in the form

Determining Multiple Energy-Parameter Sets. Multiple energy-parameter sets are determined in two steps, as shown in Figure 1.

Figure 1. Process for determining multiple parameter sets. Click to Enlarge

Step 1: Generation of multiple training sets.

Each observation in our sampling is taken as one crystal structure plus N_decoy decoy structures: C ={ r^crystal, r₁^decoy, r₂^decoy, r₃^decoy, …, r_Ndecoy^decoy}, for a given protein–ligand system (complex), where r represents the Cartesian coordinates of the structure. The size of the original sample is the number N_complex of complexes. Then, from the original sample, a bootstrap sample T^(B) = {C⁽¹⁾, C⁽²⁾, C⁽³⁾, …, C^(N_s)} of size N_s is generated by randomly sampling with replacement, with the superscript on C indicating the sampled complex. This means that some observation C in the original sample may be included several times in the bootstrap sample T^(B), i.e., C⁽ⁱ⁾ = C^(j) is allowed for i ≠ j. By repeating the above procedure N_train times, multiple training sets denoted by T_Mul = {T⁽¹⁾, T⁽²⁾, T⁽³⁾, …,T^(N_train)} are generated. In this study, N_decoy, N_s, and N_train are all 100. In various applications of bagging, N_train of order 10–10³ is empirically found to be suitable.³¹

Step 2: Optimization of the parameter set.

We want to determine the parameter set p^(B) optimized for the training set T^(B). For a set of structures C = {r^crystal, r₁^decoy, r₂^decoy, r₃^decoy, …, r_Ndecoy^decoy} of a given complex, the Z-score is defined as a function of parameter set p as follows,

Combining Multiple Predictors (Consensus Ranking Method). Multiple predictors can be combined in various ways to make a final prediction. Here, we consider two methods, rank-by-wcs (weighted consensus score) and rank-by-number. Multiple predictors under consideration are p_Mul = {p⁽¹⁾, p⁽²⁾, p⁽³⁾, …, p^(N_pred)} for BBCS and P_cs = {FlexX, DOCK, GOLD, PMF, ChemScore} for CS.

(i) Rank-by-wcs. We introduce a weighted consensus score as an extension of rank-by-rank.³⁰ First, the protein–ligand docking structures including crystal structure and decoys are ranked according to the ascending order of scores calculated by each predictor. A matrix of rank is then obtained as shown in Figure 2a. Next, for each structure, population P_R = n_R/N_pred is calculated, where n_R is the number of the predictors that rank the structure as Rth among N_pred predictors. As a result, a matrix of P_R is obtained, as shown in Figure 2b. For example, DS_1 gets first, second, and third rank by 25%, 19%, and 12% of all predictors, respectively. A weighted consensus score, S_WCS, is defined by using P_R as a weight factor:

Figure 2. Illustration of determining the weight factor of rank-by-wcs. Pred and DS represent the predictor and the docking structure, respectively. Table a of rank is converted into Table b of the population PR = nR/Npred, where nR is the number of the predictors that rank the structure as Rth, and Npred is the total number of the predictors. Click to Enlarge

(ii) Rank-by-Number. As an alternative, consensus rank by rank-by-number³⁰ is also calculated. First, an average score is calculated for every docking structure as follows:

Why BBCS Works. A semi-complete scoring function can be constructed if all possible conformations of all complexes existing in nature are used for optimizing a parameter set. Note that it is complete only within a parameter space (e.g., adjustment factor p = {α_i}) of a given scoring function. However, because it is impossible to use all complexes as a training set, some error between actually constructed scoring function and semi-complete one is inevitable. The error becomes large, especially when the number of data in the training set is small. In this case, BBCS, whose key idea is based on bootstrap aggregating (bagging) in ensemble learning, reduces the error of scores between the actual scoring function and the semi-complete one. On the basis of bagging theory,³¹ we can explain why BBCS works, for example, in case of rank-by-number. Here, e_A is defined as the error of the combined predictor:

Results and Discussion

Validation of Z-score Approach (Optimization Method). Figure 3 shows the correlation between binding scores and rmsd values of the heavy atom positions of a ligand against the crystal structure as a reference, where scores are obtained with the original FlexX parameter set⁴⁵ and a typical optimized parameter set in p_Mul. For the plot of 1d3p (Figure 3a), the original parameter set succeeds in ranking the crystal structure at the top but leaves only a marginal score gap between the crystal structure and decoys. In contrast, the optimized parameter set makes the score gap quite clear. For the plot of 1buc (Figure 3b), the optimized parameter set makes the crystal structure top-ranked (low score), but the original parameter set does not. These results are examples of how the Z-score approach successfully discriminates crystal structure from decoys. The plots with the optimized parameter set appear to describe a funnel shaped energy landscape.

Figure 3. Typical plots of binding score vs rmsd value: (a) 1d3p and (b) 1buc . The triangles represent values obtained from eq 1 by using the original FlexX parameter set and the circles represent values obtained from eq 2 by using a typical optimized parameter set in pMul. Click to Enlarge

Figure 4 shows the average Z-scores calculated by

Figure 4. Average of Z-score, ZBBCSave, for each complex in the training and test sets. The average is taken over all parameter sets contained in multiple parameter sets pMul. Click to Enlarge

It is interesting to compare Z_BBCS^ave with Z_FlexX calculated with the original parameter set⁴⁵ of the FlexX scoring function. Table 2 lists the distribution of Z_BBCS^ave − Z_FlexX. Negative values of Z_BBCS^ave − Z_FlexX mean that p_Mul provides a larger score gap than the original parameter set does, i.e., higher performance in discriminating the crystal structure from decoys. High percentages for the training set (76%) and test set (68.75%) were obtained.

Validation of BBCS. To determine whether multiple parameter sets used in BBCS are more effective than a single parameter set optimized using all complexes in the training set without bootstrapping, we compared the results of BBCS with those of one reparameterization of the FlexX scoring function, referred to as “single scoring”. To make a fair comparison, the parameter set for single scoring was also optimized with the Z-score approach.

We compared their performances against 48 complexes in the test set in two ways. The first comparison was with respect to the success rate of recognizing the crystal structure, which was calculated as the percentage of complexes whose crystal structure or low rmsd structures are top-ranked relative to decoys. The results are summarized in Table 3, where results mentioned in the next subsections are also included. BBCS attains the highest success rate of 52%, compared with 18–33% of the others in recognizing crystal structure. It is clear that both BBCSs, especially BBCS by rank-by-wcs, perform better than single scoring. This conclusion remains valid, even if we relax the success criteria by including low rmsd structures with rmsd < 1.0 or 2.0 Å as hits.

The second comparison was of the rank R of the crystal structure relative to 100 decoys. Table 4 lists the percentage of complexes giving positive, zero, and negative D_S−B, where D_S−B is an indicator defined by

Figure 5. Average of DX−B values, DX−B, for values over 48 complexes. Note that the averages of Rsingle, RFlexX and RCS are 10.4, 10.8, and 11.5, respectively. Click to Enlarge

Comparison of BBCS and the Original FlexX Scoring Function. Now that we have established the advantage of BBCS over single scoring on parameter sets determined by the same Z-score approach, let us proceed to a comparison of BBCS with the original FlexX scoring function.⁴⁵ Table 3 shows that BBCS has a higher success rate for the recognizing crystal structure or low rmsd structures than the original FlexX, especially for rank-by-wcs. Figure 6 shows the D_F−B value of each complex in the test set, where

Figure 6. DF−B = RFlexX − RBBCS for each complex in the test set, where RFlexX and RBBCS are the rank of crystal structure by the original FlexX and by BBCS, respectively. The first and second values in parentheses are RBBCS and RFlexX, respectively. (a) Rank-by-wcs. (b) Rank-by-number. Click to Enlarge

Comparing single scoring with the Z-score approach and original FlexX, one sees from Tables 3 and 4 that they are quite similar in performance in spite of the difference in optimization method.

In this paper, we have restricted ourselves to rescore conformations generated by AutoDock. When one wants to achieve best performance BBCS for the conformations generated by FlexX docking program, it may be advantageous to newly train scoring function only against a large number of conformations created by FlexX program. In practice, generated conformations are biased to some extent, because conformation search and scoring function are coupled. Thus it is not efficient to train against a set including conformations particularly disfavored by FlexX docking program. However, its determination is beyond the scope of this paper.

Comparison of BBCS and CS. The scoring functions used for CS were FlexX and D-score, G-score, PMF, and ChemScore, the latter four of which are included in CScore of the Sybyl module (Tripos Inc.). The number of predictors giving the best or better performance for each method should be used when different methods are compared. Thus, the number of predictors is 100 for BBCS and 5 for CS in this study.

Table 3 shows that compared with CS, BBCS gives a higher success rate for recognizing crystal structure, especially for rank-by-wcs. Figure 7 shows D_C−B for each complex in the test set, where

Figure 7. DC−B = RCS − RBBCS for each complex in test set, where RCS and RBBCS are the consensus rank of crystal structure by CS and by BBCS, respectively. The first and second values in parentheses are RBBCS and RCS, respectively. (a) Rank-by-wcs. (b) Rank-by-number. Click to Enlarge

CS gave the worst performance among all scorings in Tables 3 and 4; it was even worse than the predictor FlexX contained in CS. The difference between BBCS and CS lies in whether predictors are constructed based on a statistically rational method (bootstrapping) or constructed by simply gathering scoring functions. For CS, there is no guarantee that predictors will supplement each other’s weak points if ones that give quite different results or almost the same results are gathered.

CPU time required is another important standpoint. Note that as long as docked conformations are generated by single procedure and these conformations are rescored by scoring functions, CPU time is not proportional to the number of predictors, i.e. CPU time of BBCS is not simply 20 times longer than that of CS in this study. Actually, there was little difference of CPU time between them. This is because, for BBCS, once values of F_i in eq 1, are calculated as a function of coordinates, scores for any parameter sets are simply obtained as Σc_iF_i, whereas for CS, all scoring functions must be separately calculated because their function forms are different.

Worst Case by BBCS. There remains some complexes for which BBCS still fails to highly rank crystal structure among decoys; namely, the worst case is 1tet for which R_BBCS was 75th in 101 protein–ligand docking structures for rank-by-wcs. Hence, it is interesting to examine why BBCS failed for 1tet . As shown in Figure 8a, the ligand is a citrate. Figure 8b shows the difference between ligand binding sites of the crystal structure and of the top-ranked (incorrect) structure. For the crystal structure, the citrate binds to the receptor on the surface. To our knowledge, many docking programs regard such structure as energetically and structurally unstable. On the other hand, for the top-ranked structure, the citrate binds to the receptor in the pocket, which is mostly regarded as a more plausible site. A correct binding site can not be predicted. The reason for this is explained as follows. The binding mode of the crystal structure in Figure 8c shows two important histidine residues. It is likely that these imidazole rings are positively charged and make strong ion–ion interactions with two carboxylate ions of citrate, as discussed by Shoham et al.⁴⁶ In our study, however, the histidine residues were treated as neutral. That is a main reason why the crystal structure was top-ranked. However, it must be remarked that one generally assumes pH = 7.0 for a blind trial, without considering delicate protonation state of histidines and that, in this case, histidines are neutral. On the condition that histidines are neutral, it is no wonder that the incorrect structure obtained top rank, because hydrogen bonds between ligand and four residues are plausibly formed for that structure. Thus, not only the performance of the scoring function but also adequate assignment of protonation state is very important.

Figure 8. Structure of 1tet . (a) Ligand: citrate. (b) Binding site of crystal structure and the top-ranked (incorrect) structure obtained by BBCS. (c) Binding mode of crystal structure. (d) Binding mode of (incorrect) top-ranked structure. Click to Enlarge

Conclusion

BBCS was examined for pose prediction of protein–ligand docking. BBCS comes from two key ideas: the bagging to reduce the statistical error and the Z-score approach to distinguish the crystal structure from decoys. BBCS was shown to be an improvement over single scoring in a comparison of results obtained with the same Z-score approach. BBCS also performed better than original FlexX and CS in ranking the crystal structure. In particular, BBCS’s success rate of recognizing the crystal structure at the top was higher. As for the strategies of combining multiple scoring functions, rank-by-wcs performed better than rank-by-number. Once the parameter sets optimized in BBCS are obtained, one can rescore placements generated by arbitrary docking algorithm in the following steps.

(1) Perform only scoring of their placements by FlexX module to obtain individual energy components.

(2) Calculate rescores of the placements by using the parameter sets optimized in BBCS.

(3) Perform consensus scoring based on the rescores obtained.

It would be interesting to apply BBCS to other scoring functions different in form from the FlexX function.

Supporting Information Available

Sections S1 and S2. This material is available free of charge via the Internet at http://pubs.acs.org.

* Corresponding author. E-mail address: h-fukunishi@bu.jp.nec.com. Phone: +81 298 856 6155 . Fax: +81 298 856 6136.

† Nano Electronics Research Laboratories, NEC Corporation.

‡ Bio-IT Center, NEC Corporation.

§ Riken, Next-Generation Supercomputer R&D Center.

Table 2. Distribution of ZBBCSave − ZFlexX