Systematic exploration of multiple drug binding sites
BackgroundTargets with multiple (prerequisite or allosteric) binding sites have an increasing importance in drug design. Experimental determination of atomic resolution structures of ligands weakly bound to multiple binding sites is often challenging. Blind docking has been widely used for fast mapping of the entire target surface for multiple binding sites. Reliability of blind docking is limited by approximations of hydration models, simplified handling of molecular flexibility, and imperfect search algorithms.ResultsTo overcome such limitations, the present study introduces Wrap ‘n’ Shake (WnS), an atomic resolution method that systematically “wraps” the entire target into a monolayer of ligand molecules. Functional binding sites are extracted by a rapid molecular dynamics shaker. WnS is tested on biologically important systems such as mitogen-activated protein, tyrosine-protein kinases, key players of cellular signaling, and farnesyl pyrophosphate synthase, a target of antitumor agents. Electronic supplementary materialThe online version of this article (10.1186/s13321-017-0255-6) contains supplementary material, which is available to authorized users.
Background<!><!>Wrapper<!>Shaker<!><!>Systems and test metrics<!><!>Association or dissociation?<!>A systematic approach<!><!>A systematic approach<!><!>Cases with large ligands<!>Conclusions<!>
<p>Molecular docking complements experimental structure determination and it has become a standard tool of drug discovery for the determination of protein–ligand complex structures [1]. The technique in practice is a compromise between computational cost and accuracy. Its high speed necessitates the use of severe approximations such as (i) restriction of the search space to the surroundings of the binding site, (ii) no or inadequate explicit hydration of the ligand-target interface, (iii) partial or complete neglect of target flexibility [2–5] during ligand binding, (iv) and non-deterministic search algorithms [1, 6] based on random number generation. Approximations i–iv seriously limit the applicability of docking methods for the following reasons. Restriction of the search to a primary binding site requires knowledge of its location and also neglects multiple sites such as allosteric ones [7, 8]. Water molecules often play a role in ligand binding [9–11] and ignoring interfacial water positions during docking may drive the ligands into pockets which are or should be filled with water molecules, resulting in incorrectly docked ligand poses [12]. Potential water release is also important during ligand binding especially through its entropic contributions [13, 14]. Neglecting or limiting the flexibility of target molecules is obviously incorrect at binding situations with induced fit [15]. Eventuality of random number generation in search engines such as Monte-Carlo or genetic algorithms [1, 5, 6] is a natural barrier of the reproducibility and reliability of the results.</p><p>The blind docking (BD) approach was introduced [16, 17] to extend the docking search to the entire target surface. In BD, previous knowledge and restriction of the search to a primary binding site are not necessary, and therefore, it can be used in search of multiple binding sites, as well. Indeed, BD has gained popularity [18–20] and has been used for finding allosteric [21–23], or multiple [24–28] binding sites. Thus, BD addresses the above first challenge and performs a global search instead of a focused one at an increased computational cost. However, approximations ii–iv cannot be remediated as simply as the first one. Promising approaches using explicit water molecules in the binding pocket [10] (approximation ii) and treating target flexibility (approximation iii) have been reported for focused docking [29]. However, such approaches have not been implemented in conjunction with solving the global search problem of BD on the entire target surface. Statistical evaluation of multiple docking trials has been shown to increase reproducibility of a BD search [17]; by using multiple randomized (approximation iv) initial ligand positions. Thus, it has become common to perform several docking trials with different initial positions in a BD search to ensure that the largest possible part of the target surface is scanned. However, even such a statistical evaluation cannot guarantee systematic and reproducible exploration of the entire target surface during BD.</p><p>Molecular dynamics (MD) simulations have an increasing impact on drug development [30–32]. A series of pioneering studies have reported the use of MD for tracking the ligand binding process [33–37], at atomic resolution. MD calculations also allow the use of explicit water molecules and flexible targets overcoming the above limitations from approximations ii and iii [38–40] potentially opening a new avenue for improvement of BD. MD simulations typically use random starting conformations for the ligands, likewise to BD. Generally, long MD calculation times are required for successful navigation of the ligand into the binding site such that the computational time necessary for accurate docking of a ligand may be prohibitive in practice. Pocket search methods were also developed, exploiting the above-mentioned advantages of MD [41]. A recent review [30] also concludes that "Improper preparation of the initial structure or insufficient equilibration of the initial structure(s) can impact the quality of the MD results". The present study is aimed at overcoming the above uncertainties of present fast BD and molecular dynamics techniques, by combination of their advantages into a new strategy. Test applications are presented with successful identification of multiple binding sites on biologically important systems such as MAP and tyrosine-protein kinases, key players of cellular signaling as well as farnesyl pyrophosphate synthase, a target of antitumor agents.</p><!><p>Wrap 'n' Shake flowchart featuring the main steps of the method. A quick overview is also presented in Additional file 1: Supporting Movie 1</p><p>Systematic calibration of εX and RX. a A section of the VXO(r, εX, RX) LJ potential function at r = 2 Å. Scenarios Sc1-Sc3 are shown with the magnitude of the RX values corresponding to a short range repulsion of VXO ≈ 1 kcal/mol (dashed lines) b VXO LJ potential functions of scenarios Sc1-Sc3. VOO of an oxygen atom pair is also shown for comparison. c An example of excluded atoms X (red, Cycle 1, Rank 2, System 9). Docked ligand conformation, is presented with sticks and the binding pocket is shown as surface. d Ligand of System 2 (dark blue sticks) bound to its target farnesyl pyrophosphate synthase (grey lines and surface) after Shaker (MDF step). Explicit water molecules surrounding the ligand within 7 Å are shown with sticks and light blue surface</p><!><p>Finally, instead of the above-mentioned, oversimplified attempt of switching off all intermolecular terms of Einter we elaborated a new protocol which produced the desired ligand monolayer by introduction of an excluded atom type (X). In this protocol, all ligand copies docked in a cycle and their surrounding target atoms are excluded from the next cycle (red in Fig. 2c), and only the unbound target surface (grey) is used for a next BD cycle. The neighboring target atoms are selected by an interface tolerance of 3.5 Å, the maximal distance between a target heavy atom and the closest docked ligand heavy atom. The above exclusion of certain atoms during docking is physically achieved by modification of the non-bonding terms of Einter. For this, the new atom type X is assigned for excluded atoms (red in Fig. 2c) by a C program Wrp developed for this study. Wrp switches off ECoulomb by setting the partial charge of X to zero and also assigns new LJ parameters.</p><p>The new LJ parameters were fine-tuned for atom type X in order to produce the necessary weak repulsions described above. Briefly, the LJ parameters of X were calibrated considering the pairwise LJ potential between atom types X and Y (VXY) at three common atom types (Y=O, C and H). A systematic search of both equilibrium potential well-depth (εX, Fig. 2a) and inter-nuclear distance (RX) was conducted. Numerous docking runs were performed to evaluate the effect of the selected LJ parameters. A pre-defined value of r = 2 Å (ca. a covalent bond + 0.5 Å) was used as a minimal distance where short-range repulsion should act at a desired maximal value not exceeding a VXY≈1 kcal/mol. Three scenarios (Sc1-Sc3) were evaluated as shown in the r = 2 Å section of VXO (r, εX, RX) function (Fig. 2a) calculated for the XO atom type pair. Sc2 (green line, Fig. 2a, b) was identified as an optimal scenario with an εX = 10−4 kcal/mol and an RX of 3.2 Å (approximate distance between heavy atoms in an H-bond). In this case, available target surface is optimally used without large ligand-free zones in the monolayer. A short-range repulsion was achieved (green line in Fig. 2b) with a zero value beyond the repulsion zone. If RX was too large (Sc1, red in Fig. 2a, b) then the repulsion zone around the docked ligand copies would also increase with a VXO curve shifted to the right if compared to the green curve of Sc2 (Fig. 2b) resulting in large ligand-free zones, i.e. a non-optimal arrangement of the ligand copies in the monolayer. Importantly, the repulsion zone in the optimal VXO curve of Sc2 starts at lower distances (r) than in the VOO curve. VOO is shifted to the right of the red curve (Sc1), which would result in even larger ligand-free regions than Sc1. Thus, using only a repulsion term of VOO would have not been adequate for exclusions of atoms in wrapping. On the other hand, if RX was too small (Sc3, blue in Fig. 2a, b), then unwanted attractive effects such as aggregation between docked ligand copies would still happen similar to Trial 1, in Additional file 2: Table S2. Accordingly, in Sc3 the corresponding blue curve is shifted to the left from the green Sc2 curve (Fig. 2b). The same procedure was repeated for atom types Y = C and H and an average RX value of 3.6 Å was concluded (Additional file 2: Table S3) and used in Wrapper along with the above εX = 10−4 kcal/mol.</p><p>These calibrated LJ parameters of X allowed elimination of the above-mentioned unwanted interactions between the newly docked ligand copies and the already filled binding pockets (Fig. 2c). As the introduced repulsive potential acts on a short range, the ligands can still dock to other, unbound parts of the target surface. The new atom type and parameters also maximize target-ligand interactions adding the maximal number of ligand copies to the mono-layer during a BD cycle.</p><p>Wrapper cycles are terminated by either the drop of uncovered surface area of the target below one percent of its total (ligand-free, initial) surface area, or positive target-ligand interaction energy in every cluster representative (ECW in Fig. 1). As a last step, a trimming is performed to remove all ligand copies situated more than 3.5 Å from the target. Wrapper results in a target wrapped in N ligand copies (target-ligandN complex) provided as a single Protein Databank (PDB) file. Wrapper is implemented in a new open source package WnS as shell scripts and a C program Wrp available for download together with a User's Manual at www.wnsdock.xyz.</p><!><p>Shaker selects functional binding sites by removing non-specific, loosely bound ligand copies from the target surface. The target-ligandN complex is placed in a box filled with water and subjected to MD simulations in consecutive cycles. The cycles are performed until a 75% of the ligand copies are eliminated (Exit Criterion of Shaker, ECS Fig. 1). In each Shaker cycle, distance and energy metrics are calculated describing target-ligand interactions at each time step (frame) of a trajectory. The metrics include the closest distances between the target and the ligand as well as ELJ, calculated using Amber parameters. Based on these metrics, filtering (Additional file 2: Table S4) and subsequent removal of the corresponding ligand co-ordinates (Washing, Fig. 1) are applied to exclude ligand positions dissociated from their starting binding positions. The filtering involves two distance-based steps and two final steps based on ELJ.</p><!><p>Test systems</p><p>aPDB ID of the holo X-ray structure</p><p>bMolecular weight of the ligand</p><!><p>A diverse set of ten target-ligand systems were selected (Table 1) and prepared (Additional file 2: Table S1) as test cases of WnS. Challenging systems with multiple (prerequisite or allosteric) binding sites were included (Table 1). Our selection contains both small ligands and bulky, flexible ones. Apo protein structures were used as targets except System 8. In the case of System 5 another protein tyrosine-protein kinase was used as apo structure similar to a previous study [33].</p><!><p>Results for the test systems</p><p>aTotal count of ligand copies after Wrapper</p><p>bCount of ligands surviving the Shaker, after MDBSA</p><p>cRank serial number of the structure with the best RMSD value, after MDBSA and after MDF</p><p>dCount of cluster representatives (final solutions) Shaker</p><p>eShaker Rate</p><p>fTotal computational time required for MDB, MDBSA and MDF, as explained in Additional file 2: Table S12</p><p>gFor System 1, WnS was performed with different seeds for data reproduction purposes</p><p>hFinal clustering was done using van der Waals and Coulomb interactions due to interactions of zinc ion with the ligand</p><p>iWrapper process was done, using the LJ interaction as a scoring function, instead of AD4 (Additional file 2: Table S13)</p><p>jFinal clustering was done with 6 Å distance limit between clusters</p><p>Pilot molecular dynamics simulations. Benzamidine ligand (sticks) started the MD simulations from three positions at different distances (as indicated in the legend) from the native binding site on the trypsin target (grey cartoon). Arrows in a point from starting (t = 0 ns) to final (t = 1000 ns) ligand positions. Only two of the three 1000 ns-long simulations with the closest starting position succeeded in finding the reference binding pose (*) known from the crystallographic structure (3ptb). b Time-dependence of root mean squared deviation (RMSD) of the ligand measured from its reference pose</p><!><p>Analysis of the trajectories shows that the crystallographic binding position was found in two out of the three simulations after 81 and 690 ns simulation time (drop of red and green lines in Fig. 3b), respectively. In the 3rd case with the largest starting distance, 1 µs was not enough to dock to the native site by association (blue line). Thus, the usefulness of association MD runs for docking strongly depends on the starting ligand position even in the easy case of benzamidine. MD needs a simulation time comparable to the real association time of the ligand (Fig. 3b). This can be considerable, as migration of the ligand is hindered by friction in the surrounding water. Previous studies [33, 36, 45], have also reported simulations of several hundreds of nanoseconds for navigation of the ligand to the desired binding pocket.</p><p>All-in-all, the necessary time for successful docking by association MD depends on the actual starting position of the ligand, the size and shape of the target, ligand etc. To overcome such uncertainties on simulation length and still use the benefits of MD we elaborated a new strategy, the Wrap 'n' Shake (WnS, Fig. 1). Instead of simulating the association process, WnS is based on the dissociation of the ligand. Dissociation is fast and reproducible at binding sites of low stability.</p><!><p>Naturally, a dissociation approach requires a set of ligand copies bound to the target. Systematic mapping of all possible ligand positions (sites) cannot be guaranteed in a single BD cycle (Introduction) even if it contains hundreds of fast BD trials [17]. A truly systematic algorithm should completely wrap the entire surface of the target in a monolayer of copies of the ligand molecule. Our initial guess of such a Wrapper algorithm was based on a previous finding [17] that the coverage of the target can be increased with several, successive fast BD cycles where accumulated docked ligand copies from the previous cycle are considered as part of the target in the next cycle. However, additional experiments with such successive BD cycles showed that previously and newly docked ligand copies can easily form multi-layer aggregates with each-other instead of the target (Additional file 2: Table S2). The formation of such aggregates hinders wrapping of the target surface into the desired monolayer of ligand copies.</p><!><p>Wrapping tyrosine-protein kinase Src target into a mono-layer of ligand copies (System 5). a Unbound (ligand free) accessible surface area (ASA) of the target and the lowest Einter of the cluster representatives in consecutive wrapping cycles. Target-ligand interaction energy (Einter) increases with increasing number of cycles finding strong binding sites in the first few cycles, before the final, saturation region. ASA finally decreases below 1%. Structural images show the wrapping of the target (grey surface) with ligands (red). b The monolayer arrangement of the ligands (red sticks) wrapping the entire target surface (grey) after the final cycles. c A close-up of a section of the monolayer showing that the ligand copies are evenly arranged without overlap</p><p>Structural fits quantified as root mean squared deviation (RMSD) with values given in Å. Ligand conformations after Shaker (grey) compared to the crystallographic references (red sticks). System# is bold</p><!><p>For example, B-factors of measured atomic positions of ligand MSo (System 2) vary in a range between 54 and 95 Å2 (Additional file 2: Table S1). During MDF simulations we found that the RMSD varied between 2.5 and 5.1 Å (Additional file 2: Table S8), and a final #Rank of 4 and an RMSD of 3.1 Å were obtained. Considering the above high B-factor values, it is realistic to assume that ligand MSo adopts various conformations when bound to farnesyl phosphate synthase (System 2) including the one close to the assigned position found with an RMSD of 2.5 Å. This conformational variability of the bound MSo is probably due to its carboxylate group with the highest B-factor of 95 Å2. This group is hydrated by bulk water molecules, helping the dissociation of MSo from the target (Fig. 2d). At the same time, MD simulations with explicit water molecules also account for a hydrophobic, anchoring interaction between the benzofuran part of MSo (no waters present, Fig. 2d) and the target. This example shows the necessity of use of explicit water model during the shaking process in order to account for all, even antagonistic interactions.</p><p>In our pilot study ("Association or dissociation?" section) it was demonstrated that MD methods following the association pathways often need large amount of computational time and/or a fortunate starting conformation in order to find the primary site correctly for System 1. WnS yielded the correct solution for this system (Additional file 2: Table S8) in a 5-ns-long MDB simulation which is at least one order of magnitude shorter than the lengthy association times discussed in "Association or dissociation?" section. Elimination of ligand excess (dissociation of ligand copies) (Tables S14 and S15) at an SR of 11 was facilitated by hydrogen bonding with explicit water molecules [46, 47]. Thermal motion of water molecules also contributed to fast "shake off" of the ligand copies especially in the cases of Systems with small ligands with the application of the simulated annealing protocol (MDBSA, see an SR of 23 in case of System 2 in Table 2).</p><!><p>Haematopoetic cell kinase (HCK, System 5) with ligand copies remaining after Shaker. Ligand copies are colored by the calculated target-ligand interaction energy E, and the #Rank assigned. The previously reported pockets 1(ATP), 2(A-loop), 3(PIF site), 4(G-loop) and 5(MYR) are numbered by their increasing ELJ</p><p>During Shaker, conformational changes of the pentapeptide KQTSV are observed, while remains bound to its pocket on the PDZ domain (System 10). Red sticks represent the native ligand conformation from PDB (1be9). Teal sticks represent ligand conformations at different Shaker stages starting with the conformation right after Wrapper (1), and continuing with conformation after MDBSA (2), and after MDF (3). The changes of target-ligand interaction energy (ELJ) and the RMSD during the MD stages in the Shaker protocol are plotted below the structural snapshots. See also Additional file 3: Supporting Movie 2 for further details of conformational changes</p><!><p>A re-ranking (Table 2) from Rank 2 to Rank 1 was also observed after MDF. For comparison, the wrapped target-ligandN complex of System 10 was subjected directly to an MDF simulation skipping the MDB and MDBSA steps of Shaker. In this case, an RMSD of 11.3 Å (Line 10b in Additional file 2: Table S8) was obtained which was worse than the RMSD obtained with the complete Shaker protocol (1.7 Å, Fig. 5). This demonstrates that both MDB and MDBSA steps of Shaker are necessary to find the correct position. After Wrapper, the pentapeptide was in a closed, cyclic conformation (Fig. 7, Snapshot 1). This unrealistic arrangement was opened up (Snapshots 2 and 3) by interacting water molecules. It can be also observed that limited protein flexibility during MDB and MDBSA allowed only moderate reduction of the ligand RMSD by improvement of the target-ligand interactions. Most of the RMSD and interaction energy improvement was achieved after MDF, and rearrangement of K380 inside the pocket was necessary, to improve the conformation of the simulated ligand (Fig. 7). All-in-all, MD steps including target flexibility have a significant influence on the results of WnS for large ligands. Introduction of MDF considerably improved structural precision, in the above case studies of large ligands (Systems 9 and 10).</p><!><p>In the present study, a systematic strategy, the Wrap 'n' Shake was introduced for exploration of multiple binding sites and modes of drugs on their macromolecular targets. Wrap 'n' Shake systematically wraps the target into a monolayer of ligand copies using a modified blind docking approach and selects stable positions by shaking off loose binders. The method offers a computationally feasible solution for the present problems of the field (Introduction). Wrapper requires only fast blind docking cycles with a program package such as AutoDock 4.2.3. The Shaker process is fairly short and can be performed by available MD packages. Shaker is further accelerated by simulated annealing and uses all benefits of explicit water model and target flexibility. Wrap 'n' Shake is suitable to study interactions of protein targets with even large peptide ligands. We have started the extension of the method towards protein ligands using a fragment-based approach with post hoc reconstruction of the ligand. In future applications, Wrap 'n' Shake could be also used for general pocket search, besides docking of individual ligands. We envision that Wrap 'n' Shake can become the tool of choice for systematic exploration of multiple binding sites and modes of ligands in drug design and structural biology.</p><!><p>Additional file 1. Supporting Movie 1 featuring the processes of Wrapper and Shaker in the case of System 5. The first part presents the results of 15 wrapping cycles. The second part contains MDB and two MDBSA cycles of Shaker. Final cluster representatives are the outputs of WnS. Additional refinement steps are shown in Supporting Movie 2 (Additional file 3).</p><p>Additional file 2. Supporting Tables S1–S17 and Appendix 1–4 with detailed methods and results.</p><p>Additional file 3. Supporting Movie 2 featuring conformational changes of pentapeptide KQTSV, bound to PDZ-domain (System 10) during 65 ns simulations performed Shaker. The binding pocket of KQTSV on the PDZ domain is presented with grey surface. The simulated and crystallographic reference structures of KQTSV are presented as teal and red sticks.</p><p>blind docking</p><p>exit criterion of Wrapper</p><p>exit criterion of Shaker</p><p>intermolecular interaction energy</p><p>Lennard-Jones potential</p><p>molecular dynamics with backbone position restraints</p><p>molecular dynamics with backbone position restraints and simulated annealing</p><p>molecular dynamics without restraints (flexible simulation</p><p>protein data bank</p><p>root mean squared deviation</p><p>Shaker Rate</p><p>Wrap 'n' Shake</p><p>Electronic supplementary material</p><p>The online version of this article (10.1186/s13321-017-0255-6) contains supplementary material, which is available to authorized users.</p>