Binding Free Energy Calculations for Lead Optimization: Assessment of Their Accuracy in an Industrial Drug Design Context

Correctly ranking compounds according to their computed relative binding affinities will be of great value for decision making in the lead optimization phase of industrial drug discovery. However, the performance of existing computationally demanding binding free energy calculation methods in this context is largely unknown. We analyzed the performance of the molecular mechanics continuum solvent, the linear interaction energy (LIE), and the thermodynamic integration (TI) approach for three sets of compounds from industrial lead optimization projects. The data sets pose challenges typical for this early stage of drug discovery. None of the methods was sufficiently predictive when applied out of the box without considering these challenges. Detailed investigations of failures revealed critical points that are essential for good binding free energy predictions. When data set-specific features were considered accordingly, predictions valuable for lead optimization could be obtained for all approaches but LIE. Our findings lead to clear recommendations for when to use which of the above approaches. Our findings also stress the important role of expert knowledge in this process, not least for estimating the accuracy of prediction results by TI, using indicators such as the size and chemical structure of exchanged groups and the statistical error in the predictions. Such knowledge will be invaluable when it comes to the question which of the TI results can be trusted for decision making

Structure of human GS, the dimeric model system, and the binding site in the crystal structure and during MD simulations.

<p><b>(A):</b> Schematic representation of the GS decamer in top (top, left) and side (top, right) view. Subunits are labelled A to J. Below, the crystal structure of GS (PDB entry 2QC8 [<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004693#pcbi.1004693.ref022" target="_blank">22</a>]) is shown in cartoon representation. Subunits A (beige) and B (grey) used for the dimeric model system are highlighted, as in the schematic representation. (<b>B):</b> Close-up view of the dimeric model system. Subunits A (beige) and B (grey) extracted from the GS decamer are shown in cartoon representation. ADP (orange) and the GS inhibitor L-methionine-<i>S</i>-sulfoximine phosphate (MSO-P, cyan) are depicted in ball-and-stick representation bound to the bifunnel-shaped catalytic site in the interface between two subunits; manganese ions (Mn²⁺) ions are shown as black spheres. ATP binding promotes a shift of helix 8 (H8; magenta from PDB entry 2UU7 of canine GS in the <i>apo</i> form [<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004693#pcbi.1004693.ref022" target="_blank">22</a>]; violet from PDB entry 2QC8 of human GS bound to ATP and MSO-P [<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004693#pcbi.1004693.ref022" target="_blank">22</a>]) that enables glutamate binding. (<b>C):</b> Close-up view of the binding site of GS in the crystal structure with ADP (orange), MSO-P (cyan), and both mutated residues [<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004693#pcbi.1004693.ref017" target="_blank">17</a>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004693#pcbi.1004693.ref018" target="_blank">18</a>] R324 (green) and R341 (blue) in ball-and-stick representation. Mn²⁺ ions are shown as black spheres. Residue R341 is separated by ~ 10 Å from the center of the binding site (dashed line). (<b>D):</b> Backbone RMSD relative to the starting structure during 100 ns of MD simulations of the GS decamer (Decamer) and the dimeric model including all residues (Dimer) or only residues of the core region (Dimer₉₀); the GS_ADP+GGP state was simulated. The core region comprises 90% of the residues with the lowest RMSF. Respective mean RMSD values are listed in brackets; SEM < 0.1 Å in all cases. <b>(E):</b> Residue wise RMSF for subunits A and B in the GS decamer and the dimeric model system during 100 ns of MD simulations of the GS_ADP+GGP state. The table lists residues that are separated by ≤ 4 Å from ADP or GGP; regions encompassing such residues are highlighted with an arrow and labeled in the figure. (<b>F):</b> Backbone RMSD of residues listed in the table in panel E relative to the starting structure during 100 ns of MD simulations for ten dimeric pairs in the GS decamer and the dimeric model. For the decamer, the backbone RMSD was plotted as smoothed cubic spline. Respective mean RMSD values are listed in brackets; SEM < 0.1 Å in all cases. (<b>G):</b> RMSD of ADP relative to the starting structure after superimpositioning of the protein atoms during 100 ns of MD simulations for ten dimeric pairs in the GS decamer and the dimeric model. For the decamer the RMSD was plotted as smoothed cubic spline. Respective mean RMSD values are listed in brackets; SEM < 0.1 Å in all cases. (<b>H):</b> RMSD of GGP relative to the starting structure after superimpositioning of the protein atoms during 100 ns of MD simulations for ten dimeric pairs in the GS decamer and the dimeric model. For the decamer the RMSD was plotted as smoothed cubic spline. Respective mean RMSD values are listed in brackets; SEM < 0.1 Å in all cases.</p

Structural and stability changes in the R341C mutant.

<p><b>(A):</b> Close-up view of the crystal structure of human GS (PDB entry 2QC8 [<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004693#pcbi.1004693.ref022" target="_blank">22</a>]) around R341. The triad composed of residues D339, R340, and R341, and residues H281, H284, and Y288 on helix 8 (H8; raspberry) are shown in ball-and-stick representation. The salt bridge between D339 and R341 (SB¹ and SB²)) and the interaction between R340 and L-methionine-<i>S</i>-sulfoximine phosphate (MSO-P) are indicated by black dashed lines. Interactions between R341 and H281, H284, and Y288, respectively, are indicated by red dashed lines. ADP (orange) and MSO-P (cyan) are depicted in ball-and-stick representation, and Mn²⁺ ions are shown as black spheres. (<b>B):</b> Mean distances between terminal guanidino nitrogens in R340 and the oxygens of the γ-phosphate group of ATP oriented towards R340, the center of oxygens of the γ-carboxylic group of glutamate, and the carbonylic oxygen in GGP. SEM < 0.1 Å in all cases. GS_ATP, GS_ATP+GLU, and GS_ADP+GGP were considered. (<b>C):</b> Mean distances of interactions SB¹ and SB² (see panel A) for wild type GS and when considering the thiol group of C341 in the R341C mutant. SEM < 0.1 Å in all cases. Stars indicate a significant difference (<i>p</i> < 0.05) between wild type and mutant. (<b>D):</b> Mean occupancy of interactions SB¹ and SB² (see panel A) for wild type GS and when considering the thiol group of C341 in the R341C mutant. Error bars denote the SEM; stars indicate a significant difference (<i>p</i> < 0.05) between wild type and mutant. <b>(E):</b> All-atom RMSF of residue R340 in wild type GS and the R341C mutant. Error bars denote the SEM; stars indicate a significant difference (<i>p</i> < 0.05) between wild type and mutant. (<b>F):</b> Stability map depicting significant differences (<i>p</i> < 0.05) in the structural stability as computed by CNA between wild type GS and the R341A mutant. Protein structures were extracted from the GS_ADP+GGP state: Blue colors indicate that two residues are less stably connected in wild type, red colors that two residues are less stably connected in the R341A mutant. The secondary structure of GS is depicted on the top, with orange bars representing β-strands and blue bars representing α-helices; H8 is labelled. Subunits are indicated by arrows. (<b>G):</b> Probability for residues 277 to 288 of H8 to be in a loop conformation during MD simulations of wild type GS, the R341C mutant, and the HHY mutant in the GS_ATP state. Error bars denote the SEM; stars indicate significant differences (<i>p</i> < 0.05) with respect to the wild type. Results in panels B-G are based on snapshots recorded during the 20–100 ns interval of the respective MD simulations.</p

Structural changes and water structure in the binding sites of the R324S and R324C mutants.

<p><b>(A):</b> Backbone RMSD of catalytic site residues (for definition see <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004693#pcbi.1004693.g001" target="_blank">Fig 1E</a>) of the R324S (dark grey) and R324C (light grey) mutants in the GS_APO state during 100 ns of MD simulations (each subpanel shows MD simulations initiated from a different starting structure (see section “Experimental procedures” above)). Respective mean RMSD values are listed in brackets; SEM < 0.1 Å in all cases. (<b>B):</b> Mean distances between R324 (wild type GS), or S324 and C324 in GS mutants, and the β-phosphate group of ATP in states GS_ATP and GS_ATP+GLU or ADP in state GS_ADP+GGP, respectively. Stars indicate significant differences (<i>p</i> < 0.05) with respect to the wild type. In all cases, SEM < 0.1 Å. (<b>C, D):</b> Density distribution of water around ATP in the binding site during MD simulations of R324S (C) and R324C (D) in the GS_ATP+GLU state. Regions where water is most present are indicated by water density grids for three MD simulations (cyan, light blue, and dark blue; isopleths were plotted such that they encompass 80% of the maximum occupancy). ATP (orange) and the mutated amino acid 324 are shown in ball-and-stick representation. The red oval indicates an area of pronounced difference in the water density between the R324S and R342C mutants. (<b>E, F):</b> Radial distribution function (RDF) of water oxygens around the side chain oxygen or sulfur, respectively, of S324 (E) and C324 (F) in the GS_ATP+GLU state. The solid line shows the mean RDF, and dashed lines indicate ± SEM. (<b>G, H):</b> Mean relative occurrence of water-mediated hydrogen bonds between the β-phosphate group (G) or the γ-phosphate group (H) of ATP and residues S324 (gray) or C324 (white), respectively, in the GS_ATP+GLU state. The distance cutoff for the hydrogen bonds was set to 2.8 Å for strong hydrogen bonds and 3.2 Å for weak hydrogen bonds. Error bars denote the SEM; stars indicate a significant difference (<i>p</i> < 0.05) between both mutants. For panels B—H, data from the 20–100 ns intervals of the respective MD simulations was taken.</p

Mean relative effective binding energies of ATP or glutamate.

<p><b>(A):</b> Time-series of effective binding energies calculated for 4000 snapshots extracted in 20 ps intervals from the last 80 ns of MD simulations of glutamate bound to wild type GS in the GS_ATP+GLU state (black line) and least-squares line of best fit from a correlation analysis (grey line). The mean of the effective binding energies and the slope of the least-squares line of best fit are given in the legend. (<b>B):</b> Mean effective binding energies with respect to wild type GS (ΔΔ<i>G</i>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004693#pcbi.1004693.e001" target="_blank">eq 1</a>). ΔΔ<i>G</i> values were calculated by the MM-PBSA approach for ATP in the GS_ATP state and for glutamate in the GS_ATP+GLU state for GS mutants R324C, R324S, and R341C. Error bars indicate <i>SEM</i>_total (<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004693#pcbi.1004693.e003" target="_blank">eq 3</a>); stars indicate a significant difference (<i>p</i> < 0.05) between wild type and mutant.</p

Interpreting Thermodynamic Profiles of Aminoadamantane Compounds Inhibiting the M2 Proton Channel of Influenza A by Free Energy Calculations

The development of novel anti-influenza drugs is of great importance because of the capability of influenza viruses to occasionally cross interspecies barriers and to rapidly mutate. One class of anti-influenza agents, aminoadamantanes, including the drugs amantadine and rimantadine now widely abandoned due to virus resistance, bind to and block the pore of the transmembrane domain of the M2 proton channel (M2TM) of influenza A. Here, we present one of the still rare studies that interprets thermodynamic profiles from isothermal titration calorimetry (ITC) experiments in terms of individual energy contributions to binding, calculated by the computationally inexpensive implicit solvent/implicit membrane molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) approach, for aminoadamantane compounds binding to M2TM of the avian “Weybridge” strain. For all eight pairs of aminoadamantane compounds considered, the trend of the predicted relative binding free energies and their individual components, effective binding energies and changes in the configurational entropy, agrees with experimental measures (ΔΔ<i>G</i>, ΔΔ<i>H</i>, <i>T</i>ΔΔ<i>S</i>) in 88, 88, and 50% of the cases. In addition, information yielded by the MM-PBSA approach about determinants of binding goes beyond that available in component quantities (Δ<i>H</i>, Δ<i>S</i>) from ITC measurements. We demonstrate how one can make use of such information to link thermodynamic profiles from ITC with structural causes on the ligand side and, ultimately, to guide decision making in lead optimization in a prospective manner, which results in an aminoadamantane derivative with improved binding affinity against M2TM_Weybridge