<i>AutoDockFR</i>: Advances in Protein-Ligand Docking with Explicitly Specified Binding Site Flexibility

<div><p>Automated docking of drug-like molecules into receptors is an essential tool in structure-based drug design. While modeling receptor flexibility is important for correctly predicting ligand binding, it still remains challenging. This work focuses on an approach in which receptor flexibility is modeled by explicitly specifying a set of receptor side-chains a-priori. The challenges of this approach include the: 1) exponential growth of the search space, demanding more efficient search methods; and 2) increased number of false positives, calling for scoring functions tailored for flexible receptor docking. We present <i>AutoDockFR</i>–<i>AutoDock</i> for Flexible Receptors (<i>ADFR</i>), a new docking engine based on the <i>AutoDock4</i> scoring function, which addresses the aforementioned challenges with a new Genetic Algorithm (GA) and customized scoring function. We validate <i>ADFR</i> using the Astex Diverse Set, demonstrating an increase in efficiency and reliability of its GA over the one implemented in <i>AutoDock4</i>. We demonstrate greatly increased success rates when cross-docking ligands into <i>apo</i> receptors that require side-chain conformational changes for ligand binding. These cross-docking experiments are based on two datasets: 1) SEQ17 –a receptor diversity set containing 17 pairs of <i>apo-holo</i> structures; and 2) CDK2 –a ligand diversity set composed of one CDK2 <i>apo</i> structure and 52 known bound inhibitors. We show that, when cross-docking ligands into the <i>apo</i> conformation of the receptors with up to 14 flexible side-chains, <i>ADFR</i> reports more correctly cross-docked ligands than <i>AutoDock Vina</i> on both datasets with solutions found for 70.6% vs. 35.3% systems on SEQ17, and 76.9% vs. 61.5% on CDK2. <i>ADFR</i> also outperforms <i>AutoDock Vina</i> in number of top ranking solutions on both datasets. Furthermore, we show that correctly docked CDK2 complexes re-create on average 79.8% of all pairwise atomic interactions between the ligand and moving receptor atoms in the <i>holo</i> complexes. Finally, we show that down-weighting the receptor internal energy improves the ranking of correctly docked poses and that runtime for <i>AutoDockFR</i> scales linearly when side-chain flexibility is added.</p></div

Scaling of docking runtimes as function of the number of flexible receptor side-chains.

<p>The Y-axis represents multiples of the rigid cross-docking runtimes. The times used in this graph are averages taken over all docking runs for the 52 complexes of the CDK2 cross-docking experiments. For <i>AutoDock Vina</i> the times corresponding to the default exhaustiveness 8 are used. The X-axis indicates the number of flexible receptor side-chains. <i>ADFR</i> scales by a factor of 2, while <i>Vina</i>8 scales by a factor of 62, when 12 protein side-chains are made flexible.</p

Overall flowchart of ADFR.

<p>The flexibility information of the ligand (i.e. rotatable bonds) and receptor (i.e. flexible side-chains) is used to assemble the genome from which an initial list of solutions (i.e. population) is created. The population is scored, sorted, and top-ranking solutions are clustered. The GA seeds the next generation with the best solution of each cluster and completes it by crossing-over, mutating, and minimizing individuals from the mating population. The optimization stops when one of the termination criteria (maximum number of generations or evaluations) is reached or the search converges, at which point the solutions within 1 kcal/mol of the best solution are written out.</p

Comparison of side-chain conformations between <i>apo</i>, <i>holo</i>, and successfully docked solution.

<p>This figure provides a pairwise comparison of the conformations of the <i>apo</i> (4EK3), <i>holo</i> complex (1YKR), and the 1YKR ligand docked solution with the 12 flexible receptor side-chains displayed as ball-and-sticks. A) <i>Apo</i> vs. <i>holo</i>: The native bound ligand is displayed as sticks with green carbon atoms along with a partially transparent green molecular surface. The 2 lysine side-chains in the <i>apo</i> conformation severely overlap with the space occupied by the ligand. B) Docked vs. <i>apo</i>. The docked solution is shown with purple carbon atoms and partially transparent ligand molecular surface. The <i>apo</i> structure is shown with orange carbon atoms. All 12 side-chains in the docked solution adopt conformations different from the initial <i>apo</i> conformation. Most of them settle for conformations corresponding to small adjustments while others adopt substantially different conformations to resolve steric clashes (Lys33 and Lys89). C) The docked solution (purple carbon atoms) is shown with the <i>holo</i> receptor (green carbon atoms). The ligand is docked perfectly (RMSD from the crystallographic structure is 0.34Å) and the receptor side-chains changed their conformations to accommodate the ligand binding in the correct binding mode.</p

Affinity maps processing.

<p>A) A cross-section of the <i>AutoDock</i> carbon affinity map. B) The same cross-section after processing the map to create a gradient inside the protein. Besides creating a potential gradient inside the receptor, this processing also removes the local minima inside the receptor volume. The color gradient outside the protein surface indicates favorable interactions going from weak (green) to strong (blue). Inside the protein surface the color gradient indicates unfavorable interactions going from low (yellow) to highly unfavorable (red).</p

Impact of making 12 receptor side-chains flexible when docking ligands into the native <i>holo</i> receptor and the <i>apo</i> receptor.

<p>An expected loss of accuracy is observed when making the native <i>holo</i> receptor flexible, reflecting shortcomings in the scoring function and search method. Adding flexibility to the <i>apo</i> receptor, however, improves the docking success rate. <i>Holo</i> docking success rates are shown for ligand RMSD < 2Å. The success rate for <i>apo</i> cross-docking increases from 17.3% to 36.5% with a 2.0 Å RMSD cutoff. This success rate increases from 23.1% to 44.2% when using a 2.5Å RMSD cutoff (darker shade bars).</p

Heat map of ligand-flexible receptor atomic contacts reproduced in docked poses.

<p>The 43 systems reported in this table are the ones for which ADFR correctly reports the docked solution (i.e. ligand RMSD < 2.5Å). The rank of the solution among 50 GA runs is reported. White cells correspond to flexible side-chains not interacting with the ligand in either the <i>holo</i> or the docked complex. Grey cells indicate interactions formed in the docked solution, which do not exist in the <i>holo</i> complex. The remainder of the cells is colored using a red to green color scale indicating the percentage of <i>holo</i> interacting atomic pairs reproduced by the docked solution. A green cell (rate of 100%) indicates that every pairwise atomic interaction between ligand atoms and the side-chain atoms of the residue corresponding to that cell are reproduced in the docked solution. The histogram displays the percentage of <i>holo</i> interactions that are reproduced across all 12 side-chains for every ligand. The ligand reproduced at least 57.1% of all the interacting pairs in the <i>holo</i> complex, with an average of 79.8% interactions.</p

Astex Diverse Set re-docking.

<p>A) The bars depict the energy differences between lowest energy solution found by <i>ADFR</i> and <i>AD</i>2.5M (dark), and <i>ADFR</i> and <i>AD</i>25M (light). Negative values indicate a lower energy for the <i>ADFR</i> solution. Only complexes with at least one of the two differences larger than 0.5 kcal/mol are shown. 1R1H is the only complex where <i>ADFR</i> finds a significantly better solution than <i>AutoDock</i> (i.e. difference > 2 kcal/mol). B) This histogram shows the distribution of number of evaluations of the scoring function performed by <i>ADFR</i> in the GA evolution leading to lowest energy solution. C) Each docking consists of 50 GA evolutions, each producing a solution. The 50 solutions are clustered with an RMSD cutoff of 2Å. In this diagram the 85 complexes are binned based on the cluster size of the lowest energy solution indicating how many of the 50 GA runs identified the pose corresponding to the lowest energy pose found across the 50 runs, i.e. the <i>reliability</i> of the GA.</p

Cross-docking results comparison between <i>ADFR</i> and <i>AutoDock Vina</i> with 0, 4, 10, 12 flexible receptor side-chains.

<p>The left section (“Flexible side-chains”) show the number of side-chains considered flexible in each docking and the corresponding number of rotatable χ angles. The right section reports the number of systems (and percentage) for which the ligand-RMSD of the lowest energy solutions found by <i>ADFR</i> and <i>AutoDock Vina</i> are less than 2.5Å (rank 1) or a correctly docked solution is within the 10 top solutions (rank < 10).</p><p>Cross-docking results comparison between <i>ADFR</i> and <i>AutoDock Vina</i> with 0, 4, 10, 12 flexible receptor side-chains.</p

Frequency of receptor side-chain changes in the GA population during a successful docking of the 4EK6 ligand docked into the corresponding 4EK3 <i>apo</i> receptor with 12 flexible side-chains.

<p>The figure plots the evolution of the average number of receptor side-chains with a modified conformation over successive generations of the GA optimization. In the initial population all receptor side-chains are in the <i>apo</i> conformation. The number of side-chains changing rotameric state in individuals of the optimized population quickly increases in the first few generations and reaches a plateau. This profile is typical and observed in all runs for all system.</p