Abstract
Background
Proteinprotein interactions are fundamental for the majority of cellular processes and their study is of enormous biotechnological and therapeutic interest. In recent years, a variety of computational approaches to the proteinprotein docking problem have been reported, with encouraging results. Most of the currently available proteinprotein docking algorithms are composed of two clearly defined parts: the sampling of the rotational and translational space of the interacting molecules, and the scoring and clustering of the resulting orientations. Although this kind of strategy has shown some of the most successful results in the CAPRI blind test http://www.ebi.ac.uk/msdsrv/capri webcite, more efforts need to be applied. Thus, the sampling protocol should generate a pool of conformations that include a sufficient number of nearnative ones, while the scoring function should discriminate between nearnative and nonnearnative proposed conformations. On the other hand, protocols to efficiently include full flexibility on the protein structures are increasingly needed.
Results
In these work we present new computational tools for proteinprotein docking. We describe here the RotBUS (RotationBased Uniform Sampling) method to generate uniformly distributed sets of rigidbody docking poses, with a new fast calculation of the optimal contacting distance between molecules. We have tested the method on a standard benchmark of unbound structures and we can find nearnative solutions in 100% of the cases. After applying a new fast filtering scheme based on residuebased desolvation, in combination with FTDock plus pyDock scoring, nearnative solutions are found with rank ≤ 50 in 39% of the cases. Knowledgebased experimental restraints can be easily included to reduce computational times during sampling and improve success rates, and the method can be extended in the future to include flexibility of the sidechains.
Conclusions
This new sampling algorithm has the advantage of its high speed achieved by fast computing of the intermolecular distance based on a coarse representation of the interacting surfaces. In addition, a fast desolvation scoring permits the screening of millions of conformations at low computational cost, without compromising accuracy. The protocol presented here can be used as a framework to include restraints, flexibility and ensemble docking approaches.
Background
Proteinprotein interactions are essential for living organisms. They are involved in most of the key biological and biochemical processes, such as signal transduction, redox reactions, immune response and protein transport, among many others. Thus, understanding proteinprotein association is the object of an increasing interest, not only from a basic physicochemical point of view, but also for biotechnological and therapeutic reasons, with promising applications for drug design. However, experimental data on complex formation is scarce. Although in recent years techniques like NMR or Xray crystallography have fuelled the field of structural biology, the number of known proteinprotein complex structures remains low. Therefore, the prediction of proteinprotein interactions have become one of the most active and creative fields in computational physicochemistry and biology. Strictly the native complex consists of an ensemble of structures belonging to the free energy minimum of the system formed by the receptor, the ligand, and the water and ions surrounding them. Thus, in principle, using molecular simulations via molecular dynamics or Monte Carlo strategies, one should calculate the free energy of the system, and discriminate between native and non native conformations (for a recent review of free energy calculations we refer the reader to [1]). However, the amount of calculation that must be carried in a system with a pair of mediumsized proteins is prohibitive. In this way, docking strategies follow two main simplifications. On the one hand the native state of the complex is not seen as the ensemble of structures referred below, but as a single structure. On the other hand, free energy calculations are replaced by a scoring function which should distinguish the native structures in a pool of conformations.
Most of the docking algorithms are clearly divided in two fundamental parts: the search of conformations, and the scoring function. The first stage should generate a pool of conformations with a sufficient number of them similar to the native complex; the latter should discriminate between the nearnative and nonnearnative proposed conformations. One of the limitations of current docking simulations is precisely the conformational search phase, which is not able to produce sufficient number of nearnative docking solutions within the docking pool that can be correctly identified in subsequent scoring or refinement steps.
There are many available search methods, but most of them can be classified in two basic types: systematic and heuristic. On the one side, systematic sampling methods try to explore the whole or a subspace of the available conformations of the system. The large number of degrees of freedom makes unfeasible a complete conformational search in a reasonable time, so most of the reported docking algorithms treat both proteins of the complex as rigid bodies. This dramatic approximation has worked reasonably well for a significant number of proteinprotein complexes, especially those with neither hingemotions nor disorderorder transitions and where the conformational changes are limited to the sidechains. On the basis of this approach, protein docking has benefited from Fast Fourier Transform (FFT) [27], Spherical Polar Correlations [8] or GeometricHashing [9] algorithms for a fast search of the position of the interacting molecules, so these have always been the most popular docking methods. However, challenging cases involve major conformational changes upon binding, in which rigidbody docking methods often struggle to find any reasonable solution. Another major limitation is that these methods have difficulties in including sophisticated scoring functions to evaluate the docking orientations and thus are mostly based on geometrical criteria such as surface complementarity. A last (but not least) problem is that, in most of FFTbased methods, one monomer is fixed while the other is translated and rotated. This sampling strategy is obviously not the one that uses nature to face up the interacting proteins, but was chosen for practical issues of implementation. However, as recently reported [10], other sampling techniques based on orientational changes of both molecules can provide a more natural sampling and lead to an improvement of computational efficiency.
On the other side, heuristic searching methods do not need to exhaustively sample the whole conformational space, which presents some advantages, such as being able to include flexibility during the conformational search and use more sophisticated scoring functions. For instance, the ICMDISCO docking method was based on a MonteCarlo search of the ligand molecule using the binding potentials precalculated in a 3 D grid for the receptor, followed by the optimisation of the ligand interface sidechains in the torsional space [11]. These methods can provide more accurate geometries, but they are usually much slower than the systematic FFTbased search approaches.
In this paper we will try to overcome the limitations of current docking search methods, and propose an algorithm to uniformly generate docking orientations for a two rigid body system in a fast manner, with potential inclusion of flexibility. The method, to which we refer as RotBUS (
 Rot
 B
 U
 S
Methods
We describe here new procedures to generate different configurations in a twoprotein system, based on uniform sampling of the rotational space of the two interacting proteins and fast calculation of their optimal distance. The two proteins are considered as rigid bodies. In this work, the larger protein is defined as the receptor, and the smaller one the ligand. First we will define uniform distributions of points over the spheres around the interacting molecules. This will be used later to uniformly sample orientations by three different methods, inspired by early work on uniform random rotations [14]. Finally, for each orientation the optimal distance between these molecules will be calculated.
Uniform sampling of molecular orientations
Uniform distribution of points around the molecule
In order to define the different molecule orientations (as described later), we first needed to generate a uniform distribution of N points around each molecule. The total number of points N depended on the desired sampling resolution, and was calculated as follows. A triangular mesh (with equilateral faces of size ρ) was generated on the expanded molecular surface (figure 1A). In this way, the number of vertices of the mesh defines the number of points N that we need to distribute to get a sampling resolution ρ. The expanded molecular surface helped to overcome details on the molecular surface that could otherwise introduce error in the mesh generation, and was calculated by using MMTK module MolecularSurface [15] with a probe of radius 14 Å. This value is roughly the radius of gyration for a mediumsized protein, so the resolution ρ can be visualized as the distance between the center of mass of different positions of a mediumsized protein in contact with the given molecule. Once we calculated the number of points N that we needed to generate in order to do the sampling at a given resolution ρ and according to the geometrical characteristics of the protein, the next step was to distribute the N points over a sphere around the molecule, so that the orientations generated from these points could be uniformly distributed in the rotational space (on the contrary, the rotations directly based on the mesh points would not be uniformly distributed).
Figure 1. Initial distribution of points. A) A triangular mesh of points on the expanded surface of the barnase protein (PDB code 1bgs). B) On the left, a naïve distribution of points over a sphere based on uniform sampling of polar coordinates. On the right, a uniform distribution based on the method used in this work [18].
The problem of uniformly distributing a given number of points over a sphere has been studied over centuries, and constitutes one of the mathematical challenges of the XXIst century [16]. In spite of being of critical importance in many scientific disciplines, such as physics, chemistry, and biology, no analytic solution is possible. A naïve approach, based on a uniform exploration of the spherical coordinates, generates a distribution strongly biased to the poles, as can be seen in figure 1B. However, there are geometric algorithms that produce asymptotically correct solutions (i.e. they provide uniform distributions when the number of points tends to infinite) [17]. For the sake of simplicity, we have used here previously reported tables containing the distribution of point unit charges that minimizes the potential energy over a unit sphere (tables are available for N < 133, and N = 192, 212, 272, or 282) [18]. When these tables were not available, we used an algorithm based on geometrical considerations [19]. The method supposes that when the number of points is high, there is one set of points uniformly distributed that defines a tiling of the sphere by identical squares. Considering that each square edge has a certain Euclidean length, the method is able to distribute almost uniformly any number of points. This can be seen in figure 1B (on the right).
We used these uniform distributions of points around receptor (P_{rec}) and ligand (P_{lig}) to generate uniformly sets of rigidbody docking orientations by using three different methods, as explained below.
Rotational sampling method RRT ()
The first method we devised to generate uniform orientations for receptor and ligand molecules is described as follows. We first generated a series of rotations for the ligand from a distribution of points over a sphere P_{lig }using a method based on Arvo's work [20]. Each ligand rotation R_{lig }was generated by i) first performing a rotation of the molecule by a specific angle ψ around the axis that went through the center of mass of the molecule, and ii) then rotating the molecule so that the north pole pointed towards a point p ∈ P_{lig }(figure 2A). In order to ensure a uniform distribution of such defined rotations, the distributions of the angle ψ ∈ [0, 2π] and the points p ∈ P_{lig }must be uniform, and the number of sampled angle values must correspond to the square root of the number of points in P_{lig}. We have generated uniform orientations for receptor and ligand molecules as follows. Once the rotations for the ligand R_{lig }were generated as above described, then, for each of them the ligand was translated around the receptor using a uniform set of spherical coordinates (r, ϕ, θ) computed from the uniform set of points P_{rec }around the receptor (figure 2A). Thus, the direction of the translation was defined by the pair of angles ϕ, θ (i.e. the center of mass of the ligand was placed in the axis defined by each point in P_{rec }and the center of mass of the receptor), and then the value of the radial coordinate r that brings both molecules in contact was computed by a new method described later (see section "Optimal distance between interacting molecules").
Figure 2. Orientational sampling strategies used in this work. (A) RRT method; (B) TRR method; (C) RRR method. See text for details.
Rotational sampling method ()
In this second method, the ligand was first translated using the pairs of polar angles ϕ, θ defined from the uniform set of points P_{rec }around the receptor (i.e. the ligand was translated so that its center of mass was placed in every axis defined by the receptor center of mass and each receptor point p ∈ P_{rec}). For each ligand position, the ligand was rotated around its center of mass so that every ligand point p ∈ P_{lig }ended in the axis that joined the centers of mass of the two molecules (figure 2B). Finally, the ligand was rotated around the axis that joined the centers of mass of both molecules. In order to ensure a uniform sampling, as in the first method, the distribution of the angle ψ ∈ [0, 2π] and the number of sampled angle values must correspond to the square root of the number of points in P_{lig }(in a previously described method based on ICM software, http://www.molsoft.com webcite, surface points around receptor and ligand were used in a similar manner to define the translation and orientation of the ligand, but a fixed angle of 60° was used there for the final set of rotations [21]). Finally, for each orientation, the value of the radial coordinate r that brings both molecules in contact was computed as described later (see section "Optimal distance between interacting molecules").
Rotational sampling method RRR ()
In this last method, both receptor and ligand were rotated so that each one of their surface points p_{rec }∈ P_{rec }and P_{lig }∈ P_{lig }were placed in the axis that joined their centers of mass ( axis), being thus each molecule facing towards the other one. Then, the ligand was rotated around the mentioned axis (figure 2C). In order to ensure a uniform sampling, as in the other two methods, the distribution of angle ψ ∈ [0, 2π] must be uniform, and the number of sampled angle values ψ must correspond to the square root of the number of visited ligand points (those in P_{lig}). Finally, the distance r between the centers of mass of the two molecules was computed as described in section "Optimal distance between interacting molecules". This is in practice very similar to method TRR, but it reproduces in a more "natural" way the actual rotational movement of both molecules when they are interacting.
We have to note that there are other reported methods to sample rigidbody docking orientations [8,22]. For instance, Mitchell's method [22] is implemented in wellknown FFTbased docking programs like ZDOCK [6]. However, these programs usually keep a fix number of rotations for all proteins, which implies that large interfaces are sampled with lower efficiency in Euclidean distance. With our approach (i.e. defining uniform points over a sphere and then using them to generate rotations) we wanted to explore the possibility of uniform rotational sampling with fixed Euclidean resolution at the sampled interfaces independently on the size of the molecules. Our RRR rotational sampling is very similar to the method of ref [8], but they do not include the mutual receptorligand twist angle (ψ) in the formalism for uniform sampling. For instance, in their manuscript they perform a test using 492 vertices and 72 twist increments of angle ψ (about 5° angular resolution). However, for a uniform rotation plan given 492 vertices (around 10° angular resolution) only twist increments would be required.
Optimal distance between interacting molecules
For each receptor/ligand orientation, we had to compute the optimal distance r between the centers of mass of the molecules. For this, we projected the molecular surfaces of the molecules as 2 D grids on the planes defined through the geometric center of receptor and ligand, respectively, both perpendicular to the axis that joined the geometric centers of the two molecules. The distance from each grid point to the molecular surface was calculated, thus generating a distancetosurface matrix that represents the gridprojected surface for each molecule. Then, the resulting distance r was obtained by the maximum element of the sum of the distance tosurface matrices of receptor and ligand (figure 3). In order to speed up calculations, we computed these simplified gridprojected surfaces without transforming the molecular coordinates. This made the practical implementation for the last rotational sampling method slightly different from that of the other two methods. The details are explained below.
Figure 3. New method for determining the optimal separation between the interacting molecules, for each sampled orientation. (A) A triangular grid (orthogonal to the direction of the interaction v_{r}) is used to generate a 2 D projection of the receptor (in blue ribbon) surface, by computing the distance from each grid point to the furthest surface backbone atoms (grey spheres). B) The grid points are shown here coloured according to the distance to the surface backbone atoms from black to red (red indicates maximum distance to the surface). C) This process is repeated for the ligand molecule, and the two grids are summed up. The value of the maximum element of the sum matrix (plus some fixed extra space ϵ to account for the volume of the side chains) will define the final separation distance r between the centers of mass of the molecules.
For the first two rotational sampling methods, RRT and TRR, rotations and translations were applied only to the ligand, so the receptor was kept fixed. We first generated a grid representation of the molecular surface for the receptor, as follows. Let be the direction that points from the center of mass of the receptor to a given point p ∈ P_{rec }(i.e. it is the direction that joins the centers of mass of both molecules when ligand is translated to that given point p as described before). We defined and to form a basis of the perpendicular plane to , so that we generated an equilateral triangular lattice of step λ with them (fig. 3A). The equilateral triangular lattice provided a closer packing than a square one, and hence a better resolution. We used a lattice of size 17 × 15 cells, with λ = 3.5 Å. Then, for each grid point we took the farthest surface backbone atom (defined as those with a nonzero accessible surface area) inside the cylinder of radius λ/2, centered at that point with axis parallel to (fig. 3A) (the mathematical condition is that the cross product between and the vector that points from the grid point to the atom has a value smaller than ). The distance between the grid point and the selected atom is calculated for each grid point, so that we obtained a distance matrix that mimicked the shape of the receptor from each grid point (fig. 3B).
Applying the algorithm to the ligand is analogous. Given the rotation matrix R_{lig}, the ligand has its interaction on the direction given by:
To find a basis and of its perpendicular plane, we used the transformation in equation 1 over and . Then we generated a grid defined through the center of masses of the ligand and the corresponding distance matrix for the ligand (as explained before for the receptor).
Finally we summed both matrices, selected the maximum element that had both summands no null and added ϵ = 4.7 Å (arbitrary value selected from previous tests; data not shown) to this quantity in order to account for the extra volume required to include the sidechains of both proteins. The resulting value was assigned as to r, accounting for the distance between centers of mass that kept the surfaces in contact (figure 3C).
For the third rotational sampling method, RRR, both receptor and ligand molecules were rotated so that the interaction was always on the axis. Then, in order to generate the gridprojected surface of the receptor we took the direction v_{rec }(the direction from the geometric center of the receptor to the given point p ∈ P_{rec}), and the difference with the other two methods for generating rotations is that now we need to generate the vector basis depending on R_{rec }as:
Similarly, for the ligand, we had as the direction corresponding to the given point p ∈ P_{lig}, and the vector basis:
The rest was computed as above described.
Scoring the conformations
We have used pyDock [13] for the final scoring of the orientations generated by RotBUS. The function, based on electrostatics and desolvation, with weighted van der Waals term, is specially suited to study rigidbody proteinprotein docking. In the original pyDock benchmark [13], a weight of 0.1 for van der Waals term was optimal to tolerate the small number of atomic clashes found in the FFTgenerated rigidbody docking poses. However, in RotBUS docking orientations might not have such good geometric complementarity as in FFTbased approaches. Thus, in order not to filter out acceptable orientations, we had to consider lowering van der Waals weight (see Results for further details).
Given the large amount of docking orientations generated by RotBUS, pyDock evaluation became computationally too expensive at high resolution (e.g. ρ = 9 Å). Thus we devised a fast residuebased desolvation scoring. Using the atomic parameters σ_{i }in a previous work [21], we have first calculated the contribution of each surface residue to the desolvation energy in the individual proteins, :
Then, for each docking orientation the total desolvation energy was computed by summing up the residue contributions in all intermolecular contact residue pairs (i.e. residues from each molecule that have at least one atom within 8 Å distance), as in the following equation:
The contribution of each residue thus depended on its intrinsic contribution (precomputed in the individual proteins) and on the number of contact pairs formed with residues from the partner molecule. The concern with this approach was that the residues involved in several pairs might be overrepresented and thus could contribute excessively to the desolvation energy. However, we checked that in this method (as opposed to that of considering each interface residue only once) the resulting energy correlated better with the originally described atombased desolvation energy. The main advantage is that accessible surface area was calculated only once (on the individual proteins), and thus during the scoring process only distances needed to be computed. This scoring method showed to be much faster than pyDock and therefore was used as a prepyDock scoring filter.
Data set for benchmarking
We analyzed the predictive capabilities of RotBUS, using a standard benchmark for docking formed by 84 proteinprotein cases with known Xray structures both for the unbound and the bound subunits [12]. In order to perform a test as much realistic as possible, the orientation of the initial structures was randomized (the structures in the benchmark are provided as superimposed to the xray structure of the complexes, which could yield artificial better results). The success rates of the predictions were defined by the percentage of test cases in which at least one acceptable docking solution was found within the top N solutions. A nearnative or acceptable docking solution was defined as the one that had below 10 Å of RMSD of the ligand Calpha atoms with respect to the corresponding ones in the reference structure, after superimposing the receptor molecules. In order to focus our test onto the rigidbody results, for the first test regarding the number of nearnative solutions at different sampling conditions, the complex reference was formed by the unbound subunits superimposed onto the complex structure (thus avoiding the high RMSD values that cases with large unboundbound differences could have even with correct docking orientations).
Results and Discussion
Sampling efficiency by the different rotation methods
We have devised three different methods for rigidbody sampling (see Methods), based on uniform rotations, and have evaluated the quality of the generated docking orientations on a widelyused benchmark set formed by 84 proteinprotein complexes [12].
First we analyzed the percentage of cases with nearnative docking solutions (defined with respect to a reference structure formed by superimposing the unbound subunits onto the complex structure) generated by each method, at different levels of resolution (see table 1). The resolution value ρ indicated the minimal distance between the grid points on the expanded surface that were used to define the number of rotations, as described in Methods. The resolution value can be intuitively related to the distance between the geometrical centers of the different ligand orientations sampled during the procedure. For instance, a 14 Å resolution means that, for standardsized molecules, the docking orientations would be uniformly sampled so that the minimal distance between the geometrical centers of the ligands (after superimposing the receptors) would be 14 Å, while the contacting surfaces would be approximately 9 Å separated (on the other side, at 9 Å resolution, the contacting surface would be approximately 6 Å separated). In order to compare with standard angular resolution units, a resolution of 14 Å would be equivalent to 14.0° for a big protein (1de4 ligand, 1278 residues) and to 28.4° for a small one (1eaw ligand, 58 residues). On the other side, a resolution of 9 Å is equivalent to 8.5° for a large protein (1de4 ligand, 1278 residues) and 17.5° for a small one (1eaw ligand, 58 residues). In this way, for the largest cases, at 9 Å resolution we have better surface sampling than with the standard 12° resolution used in FFTbased approaches.
Table 1. Docking results with different rotational sampling methods.
The three methods had similar results. In all of them, the number of cases with nearnative docking solutions increased at higher resolution values (i.e. those with smaller grid cell size), but obviously the number of total docking poses also increased. Given that the RRR method gave slightly best results for high +resolution (below 14 Å), and since it was more intuitive (as it directly described the docking orientations based solely on molecular rotations without needing any translation), we decided to use it in this work.
Because of the systematic sampling, the number of nearnative solutions was proportional to the total number of docking poses, and this in turns depended on the resolution. At 9 Å, all cases had nearnative conformations, although the total number of generated docking poses was quite high in average (circa 1.2 million). On the other side, at 12 Å resolution, the average total number of docking poses was more reasonable (below 100,000), while the percentage of cases with nearnative solutions was still high (77 out of 84 cases). Moreover, we can see in figure 4A that the concentration of nearnative solutions at 14 Å at 12.5 Å and at 9 Å resolution are very similar. This is important, as we have recently reported that the concentration of nearnative solutions is a key factor for pyDock scoring success [23].
Figure 4. Analysis of the RotBUS orientational sampling. A) Percentage of cases with concentration of nearnative solutions (9 Å ligand RMSD with respect to a reference formed by the unbound molecules superimposed on the complex structure) above the value indicated in abscissas for the rigidbody docking sets generated by RRR method at 9 Å (blue), 12.5 Å (green) and 14 Å resolution (red). B) Success rates (i.e. percentage of cases with nearnative solutions 10 Å ligand RMSD with respect to the complex structure within the top ranked docking poses indicated in abscissas) after scoring by pyDock the docking sets generated at 14 Å resolution, using weighting factors for the van der Waals term 0 (red), 0.02 (green), and 0.1 (blue). C) Same as in B but in logarithmic scale.
For the only purposes of exploring the number of nearnative solutions found at the different sampling conditions, as described in this subsection, the ligand RMSD of each docking pose was calculated with respect to a reference formed by the unbound structures superimposed onto the complex structure (although we have to note that for all the success rate results given later, we have calculated the ligand RMSD with respect to the Xray structure of the complex, as standard). In this way, conformational changes in the subunits were ignored so that we could focus on the capability of the sampling algorithm to find correct orientations. There are 77 cases that had nearnative solutions at 14 Å resolution, but we needed to increase resolution up to 12.5 Å in order to have nearnative solutions for all the remaining 7 cases. When RMSD was computed with respect the complex structure, the only case in the combined set (14 Å/12.5 Å) that did not have any nearnative solution (i.e. RMSD < 10 Å) was 1h1v (one of the set of cases run at 12.5 Å resolution), in which the unbound ligand was 14 Å from the bound one. For these 83 complexes with nearnative solution, there was an average of 15.3 nearnative conformations (with a standard deviation of 13.2). Moreover, 53 of these cases had 10 or more acceptable solutions (defined with respect to the complex structure). For comparison, FTDock (12° angular resolution, 0.7 Å grid size, with electrostatics, generating 10,000 conformations per complex) failed to find a nearnative solution in 17 cases, found only one nearnative solution in 11 other cases, and found 10 or more acceptable solutions in only 15 cases. All these data indicate that the method presented here can properly sample the 2rigid body conformational space, with optimal parameters for the number of rotations, radial distance r and resolution.
Scoring with pyDock
Having proven the efficiency of the new sampling method in generating rigidbody docking orientations, we next focused on the scoring, with the goal of accurately identifying the nearnative docking solutions (based on ligand RMSD computed with respect to the complex structures, as standard). Figures 4B,C show the performance of pyDock scoring on the solutions generated by the RRR method at 14 Å resolution (except for the 7 cases with no nearnative solution, in which we used 12.5 Å resolution). We checked the use of van der Waals and there was a little improvement with 0.02 weight. As can be seen in figure 4B,C a weight of 0.02 gave the best success rate for top 10 solutions, and it gave better results for all rank values. This small van der Waals contribution was sufficient to remove false positives in which electrostatics was artificially high. Moreover, imposing a van der Waals cutoff to avoid cases with high van der Waals values did not change the results (data not shown). On the other side, higher van der Waals weights were not helping either, probably because of the noise derived from the rigidbody assumption (as discussed in previous work [13]). The main problem is that at 14 Å resolution, the number of docking poses was too high for the practical application of pyDock (the average computational time was more than 130 hours per case). For this reason, we further explored the use of a fast residuebased desolvation (see Methods) as a first filtering step in order to reduce the number of docking poses to be scored by pyDock.
Filter with fast solvation energy
First we studied the most efficient filtering protocol with the new residuebased desolvation in terms of recovery of nearnative solutions. When we applied the fast desolvation to the sets derived at 9 Å and 12.5 Å/14 Å resolution, we could keep a similar number of total docking poses by selecting the top 1% in the 9 Å resolution set (yielding an average of about 11800 docking poses) or the top 10% in the 12.5 Å/14 Å resolution set (yielding an average of about 11300 docking poses). However, although the final number of selected conformations was similar, we saw an important difference: the concentration of nearnative solutions (defined with respect to the unbound subunits superimposed on the complex structure) was higher after selecting the top 1% lowestdesolvation orientations in the 9 Å resolution set (figure 5A). In these conditions, the number of nearnative solutions found in the top 2000 conformations generated by RotBUS is in many cases significantly higher than when generated by FTDock (see [Additional file 1]).
Additional file 1. Number of nearnative structures in the top 2000 conformations generated by RotBUS and FTDock. We show the number of nearnative structures in the top 2000 conformations generated by the method presented here (RotBUS 9A resolution, 1% lowest residuebased desolvation), as compared to when generated by the wellknown FFTbased method FTDock. For the sake of completeness, we have also shown in brackets the number of nearnative structures in the top 2000 conformations generated by each method after scoring by pyDock.
Format: PDF Size: 23KB Download file
This file can be viewed with: Adobe Acrobat Reader
Figure 5. Analysis of the residuebased desolvation filtered subsets. A) Percentage of cases with concentration of nearnative solutions (10 Å ligand RMSD with respect to a reference formed by the unbound molecules superimposed on the complex structure) above the value indicated in abscissas at 9 Å (green) and 12.5 Å/14 Å resolution (red), when only the top 1% and 10% docking poses with the best residuebased desolvation values are considered. B) Success rates (nearnative solutions defined with respect to the complex structure) after scoring by pyDock (w_{vdw }= 0.02) the top 1% docking poses with the best residuebased desolvation of the docking sets generated at 9 Å and 12.5 Å/14 Å resolution (in red and green, respectively), and the top 10% docking poses with the best residuebased desolvation of the docking sets generated at 9 Å resolution (blue). C) Same as in B but in logarithmic scale.
In figure 5B,C we can see that filtering the set at 12.5 Å/14 Å resolution with the new residuebased desolvation permits a fast removal of 90% of the docking poses, without losing too much efficiency in the pyDock scoring. Computational times for a small case like PDB 1eaw were less than 2 minutes for the orientation search and 40 seconds for the residuebased solvation calculation on a 2.4 GHz Dual Core AMD Opteron CPU. For one of the largest cases (PDB code 1de4) these times increased up to 32 and 52 minutes, respectively. In any case, given that most of the time is spent in pyDock scoring, removal of 90% of the poses makes the total procedure around ten times faster.
Besides, given that selecting the top 1% at 9 Å resolution seemed more efficient in terms of recovery of nearnative solutions, we have analyzed pyDock success rates using this strategy (from now on, nearnative solutions are defined based on ligand RMSD computed with respect to the complex structures, as standard). In figure 5B,C (and in [Additional file 2]) we can see the results for 9 Å resolution after selecting the 1% docking poses with the lowest residuebased desolvation plus final scoring by pyDock (with van der Waals 0.02). For low ranks, 9 Å resolution (1% filtering) was clearly better than 14 Å resolution (10% filtering).
Additional file 2. Benchmark results for the RotBUS+PyDock protocol. Docking results after scoring with PyDock the set of rigid body poses generated by RotBUS at 9 Å resolution and filtered up to 1% with best residuebased solvation. RMSD is calculated for ligand Calpha atoms with respect to the complex structure.
Format: PDF Size: 23KB Download file
This file can be viewed with: Adobe Acrobat Reader
When selecting instead the top 10% (9 Å resolution), the results were slightly better at high rank values, but they were worse at low rank values (figure 5). At 9 Å resolution, computational times for a small case like PDB 1eaw were 10 minutes for the orientation search and 3 minutes for the residuebased solvation calculation on a 2.4 GHz Dual Core AMD Opteron CPU. For the exceptionally large case PDB 1de4 these times dramatically increased up to 5 and 10 hours, respectively. In any case, these times were smaller than those of FTDock (for 1eaw and 1de4, 6 and 60 hours respectively) and pyDock (the latter depending on the number of solutions to score). Considering that pyDock scoring is 10 times faster with 1% filtering than with 10% filtering, this is clearly the strategy of choice.
Integration of RotBUS and FTDock sets
In figure 6 and [Additional file 1] we have compared the results of RotBUS plus pyDock with those obtained from FTDock plus pyDock. For rank values of 10 or below, the results are very similar. However, for larger rank values the success rates are better in FTDock. We have explored whether both methods could be complementary and thus generate good solutions in cases in which the other one could have difficulties, and vice versa. We have checked that in many cases the number of nearnative solutions in the top 2000 conformations is high when generated by RotBUS or by FTDock, but not by both simultaneously (see [Additional file 1]). Figure 6A,B shows the success rates of the combined docking sets from i) 9 Å resolution plus 1% filtering, scored by pyDock with 0.02 van der Waals, with the final value weighted by a factor of 0.5, and ii) FTDock, scored by pyDock with 0.1 van der Waals. Other weighting factors for the pyDock scores of the RotBUS poses did not improve the results (see [Additional file 3]). With 9 Å filtered at 10%, the results did not improve (figure 6C,D). We checked that considering a fix number of top solutions (e.g. 10,000) after solvation filtering gave similar results (data not shown).
Additional file 3. Extra weight values for the analysis of the combined RotBUS and FTDock sets. Success rates for RotBUS (generated at 9 Å resolution and scored by pyDock with 0.02 van der Waals weight; in blue diamonds), FTDock (scored by pyDock with 0.1 van der Waals weight; in green triangles), and the combined sets from RotBUS and FTDock, in which the scoring values of the RotBUS docking poses have been weighted by different factors from 0.1 to 1. A) For RotBUS, only the top 1% docking poses with the best residuebased desolvation values were considered. B) As in A but in logarithmic scale.
Format: EPS Size: 1.1MB Download file
Figure 6. Analysis of the combined RotBUS and FTDock sets. Success rates for RotBUS (generated at 9 Å resolution and scored by pyDock with 0.02 van der Waals weight; in red), FTDock (scored by pyDock with 0.1 van der Waals weight; in green), and the combined sets from RotBUS and FTDock, in which the scoring values of the RotBUS docking poses have been weighted by a factor of 0.5 (in blue). A) For RotBUS, only the top 1% docking poses with the best residuebased desolvation values were considered. B) As in A but in logarithmic scale. C) For RotBUS, only the top 10% docking poses with the best residuebased desolvation values were considered. D) As in C but in logarithmic scale.
RotBUS with residuebased desolvation filtering in combination with FTDock gave better general results than each method on their own, especially for rank values over 50. However, we saw specific examples that were very dependent on the choice of the method. For instance, there were cases that had poor results with RotBUS (9 Å resolution) filtering at 1% residuebased desolvation, but had good results with FTDock, such as those with PDB code 1ay7, 1b6c, 1buh, 1ml0, 2btf, and 2jel. These cases can improve the results if we apply 10% filtering, which indicates that for them, geometrical complementarity rather than desolvation is the important factor to detect the nearnative solutions. On the contrary, there were cases that had good results with RotBUS at 1% filtering, but not with FTDock, such as those with PDB code 1de4, 1n2c, 2hmi, and 2qfw. Moreover, the results did not further improve with 10% filtering. In these cases, desolvation seems to be determinant (indeed, in all of these cases, except the ligand of 1n2c, the interacting molecules had highly significant ODA values [24]), and thus our residuebased desolvation is helping to identify the near native solutions even if they do not have optimal geometrical complementarity or are too large for FTDock sampling. For these cases, the FTDock results could improve by inclusion of some solvation descriptor during the FFTbased search. Indeed, we have checked that the program ZDOCK 3.0 [25], which includes solvation, gives for these cases much better results than FTDock (rank 21, 72, 544, 1 for 1de4, 1n2c, 2hmi and 2qfw, respectively), which suggests that FFTbased approaches can cover the same sets of orientations of RotBUS with appropriate inclusion of solvation. In general, in order for RotBUS sampling to have the same success rates as FTDock (especially at high rank values), we would need to use 9 Å resolution at 10% filtering. However, it seems much more efficient to combine RotBUS at 1% filtering with FTDock if we need to produce more nearnative solutions in a selection of docking poses for further refinement, for instance.
Success rates by size...does it improve FTDock?
We have studied the success rates of FTDock and RotBUS according to the size of the proteins. To be consistent with previous studies we have used the size of the grid, s, generated by FTDock [23], according to:
where is the maximum radius of the receptor (distance from the center of coordinates of the receptor protein to its farthest atom), maximum radius of the ligand, and δ the resolution of the grid used by FTDock (0.7 Å in our case). As we recently saw, pyDock had very good results on the docking poses generated by FTDock for small cases (s < 150), but very bad results for the large cases (s > 250). The reason was that for these cases, FTDock had problems in generating a sufficient number of acceptable docking solutions. We have found here that FTDock is giving slightly better pyDock success rates than RotBUS for small cases (s < 200). But for large cases (s > 250), top 100 success rates were by far much better when docking poses were generated by RotBUS (33.3%) than by FTDock (0%). Some interesting examples are the cases with PDB code 1de4, 1n2c and 2hmi, which had pyDock ranks between 12 and 32 when docking poses were generated by RotBUS, although they did not have any acceptable solution generated by FTDOCK. It is clear that the limitations of FTDock with respect to the largesized cases can be overcome with RotBUS efficient sampling.
RotBUS and beyond
When external information about the proteinprotein interaction is provided, e.g., computationally predicted hot spots, residue conservation or experimental data, the search can be restricted to defined subzones of the subunits. With RotBUS, this can be easily done by checking whether the vectors point to the expected interacting surfaces for receptor and/or ligand, respectively. If they do not, the distance r is not calculated and the orientation is disregarded. This dramatically reduces computational costs, making possible the study of larger cases.
As an example of practical application, we have performed this restricted search on two CAPRI http://www.ebi.ac.uk/msdsrv/capri webcite[26] targets, T26 (TolB/Pal) and T27 (Ubc9/E225K). We first projected the spherical uniform set of points (which represents all the directions to be explored a priori) on the surface of each subunit. For each residue expected to be at the interface, we have selected the direction corresponding to the closest projected point (figure 7A). We selected also the nearest neighbours of the corresponding direction in order to ensure a wider sampling. The expected interface residues (easily obtained from previous literature knowledge) for CAPRI target T26 receptor TolB were His246, Ala249 and Thr292, while for ligand Pal were Ala88 and Glu102 [27]. For target T27, the expected interface residues in receptor Ubc9 were Gln126, Gln130, Ala131, Glu132, Tyr134, and T135 [28], while the ones for ligand E225K were Arg8, Lys10, Arg11, Phe13 and Lys14 [29]. For target T26, the residue restraints reduced the number of explored docking poses from about 780,000 to only 1,300 orientations, at a sampling resolution of 9 Å (no filtering with residuebased desolvation). For target T27, the original sampling of about 820,000 orientations was reduced to only 6,000 ones with the restricted search. When pyDock was applied to these restraintfiltered sets of docking poses, a nearnative solution was detected at rank 1 for T26, as can be seen in figure 8A (for comparison, in our CAPRI submission, unrestricted docking obtained a nearnative solution at rank 19, which was improved to rank 6 with pyDockRST). For T27, a nearnative solution was found at rank 25 using the same restrictedsearch method (for comparison, in our CAPRI submission, unrestricted docking obtained a nearnative solution at rank 977, which was improved to rank 5 with pyDockRST) [30].
Figure 7. Restricted sampling with RotBUS. A) Diagram of the chosen points to be explored, when there is available external information on the expected binding regions. Red dot is the expected interface residue. Red lines represent the directions selected for sampling (the one corresponding to the expected interface residue, and the neighbour directions). B) CAPRI target 26 receptor and ligand respectively. The residues depicted in red are expected to be in the interface, the arrows in blue account for the selected directions to be explored in both receptor and ligand molecules.
Figure 8. Results of restraintfiltered RotBUS for the CAPRI targets T26 and T27. A) Nearnative solution (rank 1) for the complex T26, with RMSD_{lig }= 6.9 Å. Receptor in magenta, ligand in cyan and reference ligand in green. B) Nearnative solution (rank 25) for the complex T27, with RMSD_{lig }= 7.6 Å. Same colour code as in A.
Similarly, our sampling method could be easily extended or adapted to other algorithms. For instance, flexibility can be implemented by loading several structures of the proteins when the radial coordinate is being computed. These structures could be generated, for example, from a molecular dynamics ensemble. In this way, future implementations of this method may provide an efficient approach to the problem of flexibility in proteinprotein docking, especially for cases in which large conformational changes invalidate other more standard algorithms such as FFTbased methods. Moreover, this method could be implemented in combination with any other conformational search protocol, for example as part of a Monte Carlo strategy to perform simple minimizations or thermodynamic calculations.
Conclusions
We have presented here a new systematic approach to generate rigidbody orientations of a receptorligand system, based on three novel algorithms: i) uniform definition of rotations of receptor and ligand; ii) fast computing of optimal distance between proteins; and iii) fast filtering with residuebased desolvation. The method generates docking orientations at low computational cost and good efficiency, and can be complemented with those generated by other methods (e.g. FFTbased). The final scoring with the previously developed pyDock function yields competitive success rates and opens the door to efficient treatment of flexibility by using presampled ensembles or onthefly conformational search methods.
Authors' contributions
JFR devised the concept, directed the research and finalized the draft. AS performed the calculations and drafted the paper. All authors read and approved the final manuscript.
Acknowledgements
This work was supported by Plan Nacional I+D+i Grant BIO200802882.
References

van Gunsteren WF, Bakowies D, Baron R, Chandrasekhar I, Markus C, Daura X, Gee P, Geerke DP, Glattli A, Hunenberger PH, Kastenholz MA, Oostenbrink C, Schenk M, Trzesniak D, van der Vegt NFA, Yu HB: Biomolecular modelling: goals, problems, perspectives.
Angew Chem Int Ed 2006, 45:40644092. Publisher Full Text

KatchalskiKatzir E, Shariv I, Eisenstein M, Friesem AA, Aflalo C, Vakser IA: Molecular surface recognition: Determination of geometric fit between proteins and their ligands by corrrelation techniques.
Proc Nat Acad Sci 1992, 89:21952199. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gabb H, Jackson R, Sternberg M: Modelling protein docking using shape complementarity, electrostatics and biochemical information.
J Mol Biol 1997, 272:106120. PubMed Abstract  Publisher Full Text

Kozakov D, Brenke R, Comeau SR, Vajda S: PIPER: an FFTbased protein docking program with pairwise potentials.
Proteins 2006, 65:392406. PubMed Abstract  Publisher Full Text

Carter P, Lesk VI, Islam SA, Sternberg MJE: Proteinprotein docking using 3DDock in rounds 3,4, and 5 of CAPRI.
Proteins 2005, 60:281288. PubMed Abstract  Publisher Full Text

Chen R, Li L, Weng Z: ZDOCK: an initialstage proteindocking algorithm.
Proteins 2003, 52:8087. PubMed Abstract  Publisher Full Text

Vakser IA: Protein docking for lowresolution structures.
Prot Eng 1995, 8:371377. Publisher Full Text

Ritchie D, Kemp G: Protein docking using spherical polar Fourier correlations.
Proteins 2000, 39(2):178194. PubMed Abstract  Publisher Full Text

Fischer D, Lin SL, Wolfson HL, Nussinov R: A geometrybased suite of molecular docking processes.

Shen Y, Paschalidis JC, Vakili P, Vajda S: Protein docking by the underestimation of free energy funnels in the space of encounter complexes.
PLoS Comput Biol 2008, 4:e1000191. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

FernandezRecio J, Totrov M, Abagyan R: Soft proteinprotein docking in internal coordinates.
Prot Sci 2002, 11:280291. Publisher Full Text

Mintseris J, Wiehe K, Pierce B, Anderson R, Chen R, Janin J, Weng Z: Proteinprotein docking benchmark 2.0: an update.
Proteins 2005, 60:214216. PubMed Abstract  Publisher Full Text

Cheng TK, Blundell TL, FernandezRecio J: byDock: Electrostatics and desolvation for effective scoring of rigidbody proteinprotein docking.
Proteins 2007, 68:503515. PubMed Abstract  Publisher Full Text

Hinsen K: The molecular modelling toolkit: a new approach to molecular simulations.
J Comp Chem 2000, 21:7985. Publisher Full Text

Smale S: Mathematical problems for the next century.
Math Intelligencer 1998, 20:715. Publisher Full Text

Saaf EB, Kuijlaarss ABJ: Distributing many points on a sphere.
Math Intelligencer 1997, 19:511. Publisher Full Text

Sloane NJA, Hardin RH, Smith WD: [http://www.research.att.com/~njas/packings/] webcite

Lovisolo L, da Silva EAB: Uniform distribution of points on a hypersphere with applications to vector bitplane encoding.
IEE Proc  Vis., Image Signal Process 2001, 148:187193. Publisher Full Text

FernandezRecio J, Totrov M, Abagyan R: Identification of proteinprotein interaction sites from docking energy landscapes.
J Mol Biol 2004, 335:843865. PubMed Abstract  Publisher Full Text

Mitchell JC: Sampling Rotation Groups by Successive Orthogonal Images.
SIAM J Sci Comp 2008, 30:525547. Publisher Full Text

Pons C, Grosdidier S, Solernou A, PerezCano L, FernandezRecio J: Present and future challenges and limitations in proteinprotein docking.

FernandezRecio J, Totrov M, Skorodumov C, Abagyan R: Optimal docking area: A new method for predicting proteinprotein interaction sites.
Proteins 2005, 58:134143. PubMed Abstract  Publisher Full Text

Mintseris J, Pierce B, Wiehe K, Anderson R, Chen R, Weng Z: Integrating statistical pair potentials into protein complex prediction.
Proteins 2007, 69(3):511. PubMed Abstract  Publisher Full Text

Lensink MF, Mendez R, Wodak SJ: Docking and scoring protein complexes: CAPRI 3rd Edition.
Proteins 2007, 69:704718. PubMed Abstract  Publisher Full Text

Ray MC, Germon P, Vianney A, Portalier R, Lazzaroni JC: Identification by genetic suppression of Escherichia coli TolB residues important for TolBPal interaction.
J Bacteriol 2000, 182:821824. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lin D, Tatham MH, Yu B, Kim S, Hay RT, Chen Y: Identification of a substrate recognition site on Ubc9.
J Biol Chem 2002, 277:2174021748. PubMed Abstract  Publisher Full Text

Pichler A, Knipscheer P, Oberhofer E, van Dijk WJ, Korner R, Olsen JV, Jentsch S, Melchior F, Sixma TK: SUMO modification of the ubiquitinconjugating enzyme E225K.
Nat Struct Mol Biol 2005, 12:264269. PubMed Abstract  Publisher Full Text

Grosdidier S, Pons C, Solernou A, FernandezRecio J: Prediction and scoring of docking poses with pyDock.
Proteins 2007, 69:852858. PubMed Abstract  Publisher Full Text