article: Similarly to protein folding, the association of two proteins is driven by a free energy funnel, determined by favorable interactions in some neighborhood of the native state. We describe a docking method based on stochastic global minimization of funnel-shaped energy functions in the space of rigid body motions (SE (3) ) while accounting for flexibility of the interface side chains. The method, called semi-definite programming-based underestimation (SDU), employs a general quadratic function to underestimate a set of local energy minima and uses the resulting underestimator to bias further sampling. While SDU effectively minimizes functions with funnel-shaped basins, its application to docking in the rotational and translational space SE (3) is not straightforward due to the geometry of that space. We introduce a strategy that uses separate independent variables for side-chain optimization, center-to-center distance of the two proteins, and five angular descriptors of the relative orientations of the molecules. The removal of the center-to-center distance turns out to vastly improve the efficiency of the search, because the five-dimensional space now exhibits a well-behaved energy surface suitable for underestimation. This algorithm explores the free energy surface spanned by encounter complexes that correspond to local free energy minima and shows similarity to the model of macromolecular association that proceeds through a series of collisions. Results for standard protein docking benchmarks establish that in this space the free energy landscape is a funnel in a reasonably broad neighborhood of the native state and that the SDU strategy can generate docking predictions with less than 5 Å ligand interface Cα root-mean-square deviation while achieving an approximately 20-fold efficiency gain compared to Monte Carlo methods.
  A minimization approach which is specific to funnel-like functions can be based on the concept of underestimation. The existence of a funnel implies that the free energy can be locally underestimated by a convex function (Figure 1). The original free energy function is extremely rugged with a huge number of local minima even in a small region of conformational space. Yet its convex underestimator is much smoother and still captures the overall funnel-like landscape, which provides a handle to free energy minimization. The quality of minimization through underestimation depends on the choice of underestimator functions, the way they are constructed and utilized to locate the global minimum, as well as how structured the free energy funnels are in conformational space. The Convex Global Underestimation (CGU) method [26] employed canonical quadratic functions as underestimators without any cross-terms. In that case the underestimator, based on a set of local minima, can be constructed by solving a Linear Programming (LP) problem. Uniformly distributed samples in the neighborhood of the underestimator' s global minimum were then used to bias further sampling. The process was iterated with the set of local minima being updated, and the search region being reduced until certain convergence criteria are satisfied. CGU has been a very promising method with various applications in molecular structure prediction, including protein folding [27] and docking small molecules to proteins [28]. However, its restriction of using canonical quadratic functions limits its success in some cases [29], since the principal axes of the free energy surface are not necessarily aligned with the canonical coordinates. The performance further deteriorates as the dimensionality of the search space increases. We have used theoretical analysis to show and simple test problems to demonstrate that this restriction can lead to incorrect convergence [30]. Motivated by the potential advantages of underestimation and practical shortcomings of the CGU algorithm, we have recently extended the method using general quadratic underestimators and introduced biased sampling guided by the underestimator [30], [31]. Since the tightest underestimator in this class is obtained by solving a semi-definite programming problem, the method is termed SDU (Semi-Definite programming-based Underestimation). Semi-definite programming is computationally more demanding than the linear programming (still solvable in polynomial-time though) used in the CGU method. However, SDU typically requires fewer iterations and substantially improves optimization performance. The SDU method starts from a set of K local minima x1, …, xK of a funnel-like function f (x): within a given region of the search space. Throughout the course of the algorithm we maintain a set of local minima in the search region; initially. To capture the global funnel-like structure of f (x) within we construct a smooth (convex) quadratic function U (x) = x′Qx+b′x+c, where Q is a positive semi-definite matrix, , c is a scalar, and prime denotes transpose, such that U (x) underestimates f (x) at all local minima in, i. e. , U (xi) ≤f (xi) for all i = {1,2, …, K}. The tightest possible underestimator (with an L1 norm metric) can be found by solving a semi-definite programming problem [30]. U (x) is a general convex quadratic function. The underestimator U is used to guide further sampling. The minimum of U, denoted by xP and referred to as the predictive conformation, is in an energetically favorable region, and hence a new conformation can be generated by local minimization starting from xP. Additional conformations are obtained by local minimization with randomly generated starting points such that points in the vicinity of xP have a higher probability of being selected than points further away. To that end, we simply sample within using a density function shaped as −U. The set is being updated by adding these newly obtained conformations while removing unfavorable (i. e. , higher energy) conformations, and the search area is being reduced to a neighborhood of xP. Using the updated conformations in we repeat the underestimation step and the whole process is being iterated until a convergence criterion is met. The SDU algorithm has a number of favorable properties when applied to funnel-shaped functions. Theoretical analysis shows [30] that xP converges in probability to the global minimum x* of the funnel-like function f as the number of samples K grows. When applied to test functions resembling the funnel-like free energy functions, SDU has been shown outperforming CGU and a simulated annealing algorithm which adaptively tunes its parameters, with much less required function evaluations and much higher success rates [30]. The main goal of this paper is to develop and test docking methods that use the SDU algorithm in order to find the global minimum of a funnel-like function approximating the free energy over regions of the conformational space. Over the last few years we have developed a multistage docking method that starts with rigid body search based on Fast Fourier Transform (FFT), selects and clusters 1000 to 2000 low energy docked structures, and retains the 10 to 30 largest clusters for further processing [32]. The conformational space is decomposed into separate regions defined by the clusters, where free energy attraction basins are believed to exist and the free energy landscape is assumed to be funnel-like. Then SDU is called upon to locate the global minimum within each region, by utilizing such funnel-like behavior. As will be shown, in spite of the success of SDU as an optimizer for functions with funnel-shaped basins, its application to docking turned out to be far from straightforward. Although the method yields meaningful moves in either translational or rotational subspaces [30], minimization in the full 6D space of rigid body motions poses a challenge. This difficulty is well known in the robotics literature [33]. The Euclidean group SE (3), which is the semidirect product of (translations) and the special orthogonal Euclidean group SO (3) (rotations), is a nonlinear manifold and its parametrization is critical to any optimization procedure [33]. In particular, we will show that the funnel-like shape of the free energy surface is affected by the parameterization of the search space, and that underestimation tends to fail unless appropriate parameters are selected. We will describe two implementations of the SDU method: (i) SDU1, a cyclic coordinate descent strategy where rotational and translational moves alternate; and (ii) SDU2, a 5D strategy in which the distance of the two proteins is separately optimized and the SDU-driven search is performed in a lower dimensional space defined by 5 angular coordinates. Since the most energetically favorable distances occur when the two proteins are in contact but do not overlap, this strategy explores the free energy surface spanned by encounter complexes [15]. Results will show that the methods discover broad energy funnels, generate high quality docking predictions, and produce a substantial efficiency gain compared to Monte Carlo methods. The comparison with Monte Carlo methods is based on our earlier preliminary work [34] where we have tested a method similar to SDU1 against a 10-protein set using reduced Gō-type potentials. Even with these relatively smooth potentials, underestimation in the full 6D space has not been effective, and we had to rely on optimizing in each (rotational or translational) subspace separately. In the present paper, we show that the SDU1 strategy is effective against a larger benchmark set using much more refined energy potentials. However, our main contribution is the introduction of SDU2 which is more effective and improves efficiency by a factor of 10 compared to SDU1. In addition, SDU2 provides interesting biophysical insights with its resemblance to docking by repeated micro-collisions.
  To define the docking problem we fix the position and orientation of the first (receptor) protein. The six-dimensional vector x specifies the position and orientation of the second (ligand) protein. The variables s account for the side chain conformations in both proteins, and ΔG (x, s) denotes the free energy function. The minimization of ΔG (x, s) with respect to s is restricted to the side chains in the interface and is carried out by local methods. We define ΔG* (x) = ΔG (x, s*), where s* is an optimal solution; then the protein docking problem is finding the lowest minimum of ΔG* (x) in the region defined by the cluster, where x belongs to the space of rigid body motions, i. e. , the Euclidean group SE (3). As described before, the parameterization of SE (3) is critical for optimization purposes, and hence we first describe our results concerning the parameterization of the search space. The space SE (3) is the semidirect product of (translations) and SO (3) (rotations). The rotation group SO (3) is a 3-dimensional Lie group consisting of rotation matrices, i. e, , RR′ = I, det (R) = 1. The Lie algebra of SO (3), denoted by so (3), may be represented by the real skew-symmetric matriceswhere. It is well-known that the one-parameter subgroups of SO (3), i. e. , are geodesics [33], [35], [36], i. e. , the shortest paths between two points on SO (3). Moreover, for R0, R1∈SO (3), Ω = log (R0′R1) ∈so (3) and the distance between R0 and R1 can be defined by ρ (R0, R1) = ∥log (R0′R1) ∥ = ∥ω1−ω0∥. This distance is a natural Riemannian metric on SO (3), i. e. , it is bi-invariant with respect to the actions of the group (rotations). The exponential map from so (3) to SO (3) defined by Ω → R0eΩ is a local diffeomorphism, i. e. , there exist an open neighborhood of 0∈so (3) and an open neighborhood of R0∈SO (3) and an invertible, surjective, and smooth map from one neighborhood to the other whose inverse is also smooth. The local diffeomorphism induces a coordinate chart in a neighborhood of R0 that is known as exponential coordinate system. Given the definition of so (3) this coordinate system can be parameterized by. Given the favorable properties of SO (3) (more generally SO (n) ) such as the existence of a natural bi-invariant metric and, in particular, the simplicity of determining the geodesics of the manifold, the standard gradient based optimization algorithms on can be generalized for optimization on SO (n) [35]. However, generalization to the entire SE (3) is more difficult [33], [36]. Although SE (3) is also Lie group, its one-parameter subgroups are no longer geodesics. Moreover, SE (3) does not admit a natural (bi-invariant) Riemannian metric while its subspaces SO (3) and do [33], [36]. Attempting to implement SDU to SE (3) by simultaneous translational and rotational optimization, we have generally failed to construct useful underestimators. We have, however, found two strategies that were able to overcome this problem. The SDU1 and SDU2 algorithms were tested on the protein pairs given in protein docking benchmark sets [38], [39] that contain enzyme–inhibitor, antigen–antibody, and “other” types of complexes, using the independently determined (unbound) structures of the component proteins. The algorithms were used to refine the top (most populated) clusters of docked structures generated by the rigid body docking program PIPER [22]. We consider the 10 largest clusters for enzyme–inhibitor complexes, and 30 clusters for antigen–antibody and “other” complexes. Table 1 shows results both from the rigid body docking and the SDU-based refinement procedures for all three types of complexes, each defined by its Protein Data Bank (PDB) code [5] in Column 1 of Table 1. We emphasize that for most complexes we docked the unbound (separately crystallized) protein structures rather than their bound conformations obtained by separating the complex. The exceptions are mostly a few antigen–antibody complexes for which no separate antibody structures were available and hence were taken from the complex. However, even in these cases we used the separately crystallized structure of the antigen as given in the benchmark sets [38], [39]. PDB codes for these “semi-bound” complexes are shown in bold italic fonts in Table 1. Columns 2–4 of Table 1 describe the docked complex structures generated by the rigid body docking [22] before any refinement. As described in the Methods, the PIPER docking program evaluates the energy for billions of docked conformations. We retain the 1000 best scoring structures, and cluster them using the pairwise RMSD as the distance measure and an optimally selected clustering radius. We have observed that the near-native structures tend to be in one of the largest clusters, and hence rank the clusters on the basis of their size. In fact, calculated for the rigid protein structures, the energy function is approximate, and better discrimination of the native structures can be achieved by focusing on the large clusters. The properties shown in Table 1 are the number of near-native conformations (or “hits” with less than 10 Å ligand interface Cα RMSD) among the 1000 best scoring structures retained from the PIPER results, the Root Mean Square deviation (RMSD) between the native ligand structure and the docked structure at the center of the first cluster that includes a near-native conformation, and the rank of the particular cluster based on cluster size. Results are fairly good for enzyme–inhibitor complexes, as the 1000 structures, on the average, include over 200 hits, and the average RMSD is less than 5 Å between the native structure and the center of one of the three largest clusters. Although PIPER yields almost 200 hits for the “other” types of complexes, discrimination by cluster size is more difficult, and retaining the top 6 clusters results in 6. 16 Å average RMSD. For antibody-antigen complexes PIPER generates much fewer hits, and on the average we have to retain 15 clusters to have a near-native structure in them. The next 6 columns in Table 1 show the results for the SDU1 and SDU2 algorithms, in each case listing the RMSD for the lowest energy conformation in the first cluster that includes a near-native structure, the rank of the cluster, and the average number n̅ of function evaluations. It is important to note that after the SDU refinement the clusters are ranked based on the energies of the SDU solutions rather than on cluster size. For comparison we also generated 2000 conformations (based on the original protocol) for each complex using version 2. 3 of the rigid body docking program ZDOCK [21], and refined the structures using RDOCK [40], which performs local energy minimization. The last two columns of Table 1 show the RMSD of the first near-native structure found and its rank based on the RDOCK energy function. For enzyme–inhibitor complexes both SDU1 and SDU2 give excellent results. On the average, the RMSD is reduced by almost 1 Å and the average rank is around 2. As the overall performance measure of the refinement we consider the number of complexes that have at least one prediction with less than 5 Å ligand interface Cα RMSD from the native structure in the top 5 clusters. Such predictions would be termed “5 Å models” later for simplicity. As noted, the clusters of PIPER-generated structures are ranked based on their size (i. e. , the number of structures). However, after the refinement the clusters are ranked based on the energy of their lowest energy structures. Among the PIPER-generated structures the top 5 clusters include less than 5 Å ligand interface RMSD predictions for only 11 of the 25 enzyme–inhibitor complexes. As shown in Table 1, both SDU1 and SDU2 increase this number to 17, i. e. , a more than 50% improvement relative to the PIPER results. Notice that the 5 lowest energy predictions obtained by the ZDOCK/RDOCK procedure include 5 Å models only for 7 complexes. However, from the ZDOCK/RDOCK runs we retain low energy models, whereas the results provided by SDU1 and SDU2 are low energy clusters, and ranking clusters rather than individual structures generally removes some false positives, i. e. , conformations that have low energy but are far from the native. It is well known that antigen–antibody and “other” complexes are more difficult to predict than enzyme–inhibitor complexes [20]. For antibody-antigen pairs only the refinement by SDU2 improves the PIPER results. Although both the RMSD and the average rank of the first near-native cluster are reduced, at 12. 7 the latter remains high (Table 1). Similarly, for the complexes in the “other” category only SDU2 improves both the RMSD and the rank. Even SDU2 yields only a total of seven 5 Å models in the top five clusters for antibody–antigen and “other” types of complexes. This result is somewhat disappointing, but note that docking by ZDOCK and refining by RDOCK leads to 5 Å models only for three complexes, and no near-native solution is found in four cases. Figure 2 shows the RMSD vs. calculated free energy for the 25 enzyme–inhibitor complexes. Again we show PDB codes for the “semi-bound” complexes in bold italic fonts. Each point represents a structure sampled in the process of refining the 10 largest clusters using the SDU2 algorithm. The encircled blue asterisk indicates the native structure and the first hit is shown as a red square. In spite of the one-dimensional representation of the energy function defined in the 6D rotational and translational space, the figure demonstrates the multi-funnel behavior over a relatively broad region (within 20 Å RMSD) of the native state. For most complexes the figure shows a well defined deep funnel within 5 Å RMSD from the native structure. According to Table 1, for 12 of the 25 complexes (including 5 “semi-bound” cases), this funnel is deepest among the 10 clusters sampled. For example, for 4HTC the near-native cluster is the 7th largest, but it is energetically the most favorable after refinement. For the remaining 13 complexes, clusters farther from the native structure yield the deepest funnels, resulting in false positive predictions. One of the worst behaviors can be observed for 2PTC with the first near-native cluster ranked 4, and Figure 2 shows a number of deep non-native funnels. However, on the average, one of the 3 deepest funnels is near-native (Table 1). This also shows the power of heuristic combination of entropic and enthalpic measurement, i. e. , cluster size as a filter and refined cluster depth as ranking parameter. Figure 3 shows the RMSD vs. calculated free energy for the antigen–antibody and “other” complexes, sampled in the refinement of the 30 largest clusters by SDU2. All but 1SPB and 2BTF have multiple funnels, and the funnels are further from the native state than for enzyme–inhibitor pairs, demonstrating the well-known difficulty of estimating free energy using simple models, particularly for antigen–antibody complexes. Adjusting the conformation of interface side chains only by local minimization and keeping the backbone rigid also limits the accuracy of free energy calculation. As shown in Table 1, the SDU2 strategy is much more efficient than SDU1, and also provides substantially better results for the more difficult problem of docking antigen–antibody and “other” complexes. These differences are demonstrated in Figure 4A and 4B that show, for the near-native cluster of the complex 4HTC, the conformations sampled and minimized by each algorithm. The horizontal and vertical axes, respectively, represent translational and rotational distances between each sampled ligand conformation and the one in the native structure, the latter placed at the origin of this coordinate system. The rotational distance here is defined as the length of the minimum geodesic (i. e. , the minimum rotation in radian) between two rotations [36]. As discussed earlier, this distance between R0 and R1 can be defined by ρ (R0, R1) = ∥log (R0′R1) ∥ = ∥ω1−ω0∥. The points are color-coded according to their energies, from low (blue) to high (red) energy structures. As shown in Figure 4A, the separate treatment of translational and rotational subspaces in the SDU1 algorithm is highly inefficient, resulting in the sampling of many conformations with relatively high energies. Although all energy values shown are obtained by local minimization, the latter is unable to reduce the energy if the two proteins are too far from each other, which frequently occurs with the SDU1 algorithm. In addition, as indicated by the parallel lines, very similar translational regions are re-sampled at slightly different rotational coordinates, and vice versa. In contrast, as shown in Figure 4B, the SDU2 algorithm smoothly and efficiently descends toward the bottom of the free energy funnel. Based on the ranges of sampled free energy values, SDU2 sampled much lower energy regions compared to SDU1. Since the search is restricted to biophysically meaningful encounter complexes, the more consistent energy values facilitate the construction of better underestimators during the search. According to Figure 4B, for SDU2 there is a clear trend that structures closer to the native complex generally have lower energies, resulting in a deep and broad free energy funnel in this space. The existence of such free energy funnels is much less obvious, even in each separate subspace, when sampled by the SDU1 algorithm (Figure 4A). In test docking problems even the slower SDU1 algorithm outperformed a standard Monte Carlo method by reducing the number of function evaluations by a factor of two [34]. As shown in Table 1, SDU2 further reduces the computational costs by a factor of 11 to 15, depending on the type of complex. Since most of the computational time is spent in energy evaluations, the computational gain of SDU2 over SDU1 is more than a factor of 10, and we estimate that SDU2 achieves more than 20-fold efficiency gain compared to Monte Carlo methods. We tested our algorithms on a 128-node biowulf cluster (IBM eServer×Series). Each node contains dual 1 GHz PIII processors with 2 GB memory. A typical refinement by SDU2 for each PIPER-generated cluster would take 2 to 6 processor-hours. The running time varies with the protein complex size (especially the interface size) and the number of iterations before convergence. Notice that on the average SDU2 samples only about 100 encounter structures for each cluster. No particular efforts have been made to accelerate either the interface side-chain search or the line search to determine the center-to-center distance. The CPU times can be compared to those reported for a server [41] based on the RosettaDock algorithm [42]. The server performs 1000 independent Monte Carlo simulations within 30 Å Cα RMSD of a starting structure, as described in the RosettaDock protocol [42]. A typical run requires about 65 processor-hours. SDU based docking algorithms are stochastic in nature. Although we gave a theoretical guarantee of probabilistic convergence to the global minimum of funnel-like functions under some fairly general conditions [30], practical protein docking problems do not necessarily satisfy all these conditions. To check the variations in our results we have repeatedly run SDU2 to refine the near-native cluster from PIPER (columns 3 and 4 in Table 1) for 10 randomly selected complexes, including 4 enzyme–inhibitor, 4 antigen–antibody and 2 “other” complexes. SDU2 was run 10 times independently for each complex with the same set of parameters. The results (ligand interface Cα RMSD versus corresponding free energy values) are shown in Figure 5, where black circles represent the rigid body predictions from PIPER before refinement and blue asterisks represent the independent SDU2 predictions. Some of the latter overlap, resulting in less than 10 distinct solutions. SDU2 is able to lower the free energy values in all cases and to improve RMSD in most of them. Table 2 shows the mean (indicated by overline) and the standard deviation (indicated by σ) for both the ligand interface RMSD and the free energy for each complex. Although the average standard deviation is relatively large (about 4. 4 Kcal/mol) for free energies, it is less than 0. 5 Å for RMSDs. This result indicates a good level of robustness, considering that our free energy model does not discriminate among structures within 1 Å RMSD from each other. Note that the SDU-based algorithms rely more on the collective distribution of a set of encounter structures in the free energy funnel, rather than on a single low energy structure, which reduces the sensitivity of the results to the variations in the starting structures. The successful application of the Semi-Definite programming-based Underestimation (SDU) search algorithm to protein–protein docking further validates the assumption that the free energy landscape of the complex is a funnel in some neighborhood of the native state. However, the direct application of SDU in the space SE (3) of rotations and translations fails to yield useful underestimators. Alternating searches in rotational and translational subspaces yields a feasible but inefficient algorithm. We have substantially improved performance by separately optimizing the center-to-center distance and describing SE (3) in terms of five angles. It is potentially important that this strategy samples encounter complexes, and hence it is reminiscent of the model of molecular association through a series of micro-collisions [37]. Results emphasize that the funnel-like shape of the free energy surface seen in this parameterization of SE (3) is largely lost when changing to the straightforward description of the space in terms of rotational and translational coordinates. The underestimation approach has been used in the latest rounds of CAPRI with considerable success [19], [43], and it provides a promising platform for improving docking methods. We note that Marcia et al. [44] recently reported the application of SDU to the docking problem using the general quadratic underestimation method we have earlier developed [30], [31]. However, the central problem of parameterizing the search space was not discussed and the method was applied only to five bound docking problems using co-crystallized structures, which is much easier than docking two separately crystallized proteins. In fact, we believe that the separate minimization along the center-to-center vector and the use of five angular descriptors can improve the performance of any minimization algorithm used for docking.
  The SDU1 and SDU2 algorithms were tested on complexes from the protein docking benchmark sets [38], [39]. These sets contain enzyme–inhibitor, antigen–antibody, and “other” types of complexes. As described, the algorithms were used to refine the separate clusters generated by the rigid body docking program PIPER. For enzyme–inhibitor complexes we refined only the 10 largest clusters and hence restricted consideration to complexes from the two benchmark sets for which these clusters included at least one near-native conformation. In addition, we disregarded four enzyme–inhibitor pairs that form oligomeric rather than binary complexes, resulting in the 25 test problems. The refinement algorithms were also applied to 11 antigen–antibody and 8 “other” type of complexes from the benchmark set 1 [38] that had at least one near-native structure in the 30 largest clusters. We emphasize that tests for most complexes involved separately determined protein structures as given in the benchmark sets [38], [39]. The exceptions are a few of the antigen–antibody complexes in which the antigen was separately solved but the antibody structure was taken from the complex. The rigid body docking program PIPER [22], based on the FFT correlation approach, systematically samples billions of docked conformations on a grid. Compared with other FFT-based approaches that use only shape complementarity and electrostatics for scoring, the scoring function in PIPER also includes the statistical pairwise potential DARS (Decoy As Reference States) [22]. Since the potential is represented as the sum of a few correlation functions through the eigenvalue-eigenvector decomposition of the matrix of the DARS interactions energy coefficients, the energy can be very efficiently evaluated using Fourier transforms. In conjunction with this higher accuracy scoring function, the PIPER program significantly enriches the hit rates among top ranked predictions for the benchmark sets described above. For each complex, we retained the 1000 lowest energy predictions and clustered them using pairwise RMSD as the distance metric [45]. The resulting clusters were ranked based on their size, reflecting a preference for local minima with broad regions of attraction [32]. We have retained at most 30 clusters for each complex, each one being roughly 10 Å RMSD wide. Docking by the SDU algorithms involves the use of two different free energy models. In the rigid body global search the scoring function iswhere the desolvation free energy ΔGdes is estimated by the Atomic Contact Potential [49], an atom-level extension of the Miyazawa-Jernigan potential [50], the electrostatic energy ΔEelec is based on the Coulombic formula with distance-dependent dielectrics ϵ = 4r, and the Van der Waals term ΔEvdw is adopted from the Charmm potential [51]. The scaling factor λ∈[0,1] is dynamically adjusted during the course of the refinement to improve the quality of the underestimator by dampening the effect of the van der Waals term ΔEvdw and thereby smoothing ΔG. Specifically, whenever a “flat” underestimator U (this can be determined when the minimum eigenvalue of Q is close to zero) is being computed with λ = 1, we gradually reduce λ with a stepsize of 0. 1 until we obtain a more informative underestimator or λ reaches 0. In the flexible local minimization we use the Charmm potential with ϵ = 4r, including the internal energy terms, and perform 100 steps of adopted base Newton-Raphson (ABNR) minimization allowing for side-chain flexibility in the interface [51]. The distance between two components is separately optimized. Specifically, the ligand is pushed towards or pulled apart from the receptor at a stepsize of 0. 5 Å along the line segment connecting the two centers of mass. The maximum shift in distance (which defines locality) is 2 Å in this study. Each position is followed by a Charmm minimization described above and energy evaluation. Only non-clashing structures are accepted, judged by the condition ΔEvdw<0. In fact, we found that PIPER had a higher tolerance toward clashing structures and pushing component proteins closer generally resulted in increased positive ΔEvdw. To reduce calls to Charmm minimization, in practice we only pull them apart if necessary, i. e. , when ΔEvdw>0 is found. For simplicity we will call the work involved in evaluating ΔG* (x) for each conformation x a function evaluation, although it generally involves several evaluations of ΔG (x, s). As a measure of prediction quality, we select the ligand Cα atoms in the binding site, and calculate the RMSD between their predicted and observed positions. A ligand residue is considered to be in the binding site if any of its atom is within 10 Å of an atom on the receptor. We refer to a structure as near-native (or a “hit”) if its ligand binding site Cα RMSD is less than 10 Å. Although such structures are not really close to the native complex, by rigid body docking it is generally difficult to obtain better results. In fact, since the near-native binding region is selected by cluster size using a clustering radius on the order of 10 Å, the goal of this first step is to generate as many such 10 Å RMSD structures as possible. To show the improvements due to the SDU method, we list the rank of the first cluster that includes a hit, as well as the RMSD between the native structure and the center of the cluster. However, as a more appropriate overall performance measure of the refinement, we also note the number of complexes that have at least one prediction with less than 5 Å RMSD from the native structure in the top 5 clusters.
summary: Protein–protein interactions play a central role in various aspects of the structural and functional organization of the cell, and their elucidation is crucial for a better understanding of processes such as metabolic control, signal transduction, and gene regulation. Genomewide proteomics studies, primarily yeast two-hybrid assays, will provide an increasing list of interacting proteins, but only a small fraction of the potential complexes will be amenable to direct experimental analysis. Thus, it is important to develop computational docking methods that can elucidate the details of specific interactions at the atomic level. Protein–protein docking generally starts with a rigid body search that generates a large number of docked conformations with good shape, electrostatic, and chemical complementarity. The conformations are clustered to obtain a manageable number of models, but the current methods are unable to select the most likely structure among these models. Here we describe a refinement algorithm that, applied to the individual clusters, improves the quality of the models. The better models are suitable for higher-accuracy energy calculation, thereby increasing the chances that near-native structures can be identified, and thus the refinement increases the reliability of the entire docking algorithm.
keywords: computational biology/macromolecular structure analysis
  computer science/applications
  biophysics/structural genomics
  biophysics/theory and simulation
  biophysics/biomacromolecule-ligand interactions
  computer science/numerical analysis and theoretical computing
title: Protein Docking by the Underestimation of Free Energy Funnels in the Space of Encounter Complexes