Electrostatic complementarity as a fast and ... - ACS Publications

Subscriber access provided by Washington University | Libraries

Article

Electrostatic complementarity as a fast and effective tool to optimize binding and selectivity of protein-ligand complexes Matthias Rolf Bauer, and Mark Denis Mackey J. Med. Chem., Just Accepted Manuscript • DOI: 10.1021/acs.jmedchem.8b01925 • Publication Date (Web): 26 Feb 2019 Downloaded from http://pubs.acs.org on February 27, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Medicinal Chemistry

Electrostatic complementarity as a fast and effective tool to optimize binding and selectivity of proteinligand complexes Matthias R. Bauer‡,*, Mark D. Mackey‡ Cresset, New Cambridge House, Bassingbourn Road, Litlington, Cambridgeshire, SG8 0SS, UK Electrostatic scoring, electrostatic complementarity, molecular recognition

ABSTRACT

Electrostatic interactions between small molecules and their respective receptors are essential for molecular recognition and are also key contributors to the binding free energy. Assessing the electrostatic match of protein-ligand complexes therefore provides important insights into why ligands bind and what can be changed to improve binding. Ideally, ligand and protein electrostatic potentials at the protein-ligand interaction interface should maximize their complementarity while minimizing desolvation penalties. In this work, we present a fast and efficient tool to calculate and visualize the electrostatic complementarity (EC) of protein-ligand complexes. Using several benchmark sets compiled from mainly electrostatically driven SAR, including data of the PPI target XIAP and the GPCR mGLU5, we demonstrate that the EC method can visualize, rationalize,

ACS Paragon Plus Environment

1

Journal of Medicinal Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 66

and predict electrostatically driven ligand affinity changes and help to predict compound selectivity. The methodology presented here for analysis of electrostatic complementarity is a powerful and versatile tool for drug design.

INTRODUCTION Typical electrostatic interactions between small molecules and their receptors include hydrogen bonding, ionic/coulombic, cation-π, π-π, lone-pair sigma-hole (e.g., halogen bonding), and orthogonal multipolar interactions (e.g., fluorine bonding).1–6 These interactions are essential for molecular recognition and are also key contributors to the binding free energy (particularly the enthalpic term) of protein-ligand complexes. Assessing the electrostatic match of protein-ligand complexes therefore provides important insights into why ligands bind and what can be changed to improve binding. Ideally, corresponding ligand and protein electrostatic potential (ESP) values on the protein-ligand interaction interface should be complementary (i.e. the same magnitude but with opposite sign). Tools for assessing the degree of electrostatic complementarity of host-guest complexes should inform the design of polar, enthalpic binders, which have typically better selectivity and pharmacokinetic parameters than entropic binders and have even been suggested to be ‘better’ drugs.7,8 Electrostatic potential (ESP) calculations on small molecules have been employed since the 1970s to predict molecular properties such as molecular reactivity or characterize biological recognition processes.9,10 Both ab initio methods and approaches that rely on atom-centered partial charges have been used.7,8 Davis et al. described the use of molecular fields11 and in particular the pattern of local ESP minima to analyze and guide the design of selective PDE


2

Page 3 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


inhibitors.12 A review by Kukic and Nielsen summarized different molecular-mechanics-based approaches for modelling electrostatic interactions, including implicit solvation models such as Poisson-Boltzmann and Generalized-Born as well as all-atom models.13 Although the computationally inexpensive atom-centered charge (ACC) approaches for ESP computation are frequently used, they cannot properly describe directional electrostatic interactions such as halogen bonding or hydrogen bonding, which require atomic charge anisotropy (e.g., halogen sigma hole or lone electron pairs of carbonyl groups). The concept and usefulness of analyzing protein-ligand electrostatic complementarity (EC) has been described and investigated in many studies. In the 1980s, Weiner et al. visualized and qualitatively analyzed the EC in ligand-macromolecule complexes using Mulliken net atomic charges obtained from STO-3G and CNDO/2 calculations.14 Naray-Szabo et al. defined the main qualitative features of electrostatic complementarity and analyzed trypsin-ligand complexes.15,16 Nakamura et al. were the first to attempt EC quantification of protein-ligand complexes.17 In a series of consecutive articles, Chau and Dean progressed the concept of EC from qualitative to quantitative analysis by computing various parameters of diverse PDB structures using Pearson’s R and Spearman’s Rho rank correlation coefficients. Their studies concluded that EC has to be optimized as a global property instead of analyzing isolated fragment ligand moieties and that in many cases principal charge patterns exist that are largely responsible for the observed EC.18–20 Later studies by Kangas and Tidor investigated the optimization of binding electrostatics in terms of balancing favorable electrostatic interactions and desolvation contributions in a novel theoretical framework using continuum electrostatic theory (Poisson-Boltzmann), and applied the concept to analyzing the barnase-barstar protein-protein interaction (PPI) and to optimizing putative transition state analogues of chorismate mutase.21–24 Naray-Szabo used EC to analyze


3


Page 4 of 66

sulfur macrocycle packing and showed correlation of EC with experimental binding affinities of a ligand against five serine proteases.25 Muzet et al. derived ESPs experimentally using an ultrahigh resolution X-ray structure of an inhibitor-aldose-reductase complex, and compared their results with quantum-mechanical (QM) density functional theory (DFT) and molecular mechanics (MM) AMBER calculations.26 The study concluded that the co-factor shows an especially high EC to the aldose-reductase binding pocket, that experimental and DFT-computed ESPs agree reasonably well, and demonstrated the relevance of lone electron pairs and atomic polarization effects to the description of host-guest complexes.26 Huggins et al. discussed EC in the context of improving ligand selectivity, describing examples of designed electrostatic ligand selectivity for inhibitors of trypsin/thrombin, thrombin/FXA, and PTP isoforms. Their conclusion was that EC provides a direct means for gaining selectivity and that optimizing EC can be a particularly effective strategy for charged or highly polar target binding sites.27 Several recent studies have made use of electrostatic analysis to study the molecular adhesion between collagen and proteins belonging to the integrin family with respect to the EC principle,28 to screen for small molecules that mimic the electrostatic properties of protein ligands at PPI interfaces,29 to rescore docking poses with an EC term,30 and to include a protein-ligand EC parameter in the newly developed scoring functions.31,32 Despite the development of several electrostatic potential/complementarity metrics with application to protein-ligand complexes,21–23,33–40 the available tools either lack the description of directional electrostatic interactions due to use of atom-centered charges (molecular mechanics Poisson-Boltzmann methods) or are routinely employed only for ligand electrostatics because of prohibitive computational costs and difficult setup for large macromolecules when ab initio methods are used.41


4

Page 5 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In this work, we present a fast and efficient method to visualize and calculate the electrostatic complementarity of protein-ligand complexes. The electrostatic calculations use the polarizable eXtended Electron Distribution (XED) force field,42,43 which, unlike classical force fields, employs off-center charges in addition to atom-centered charges. The XED charge model provides a better description of systems involving electron anisotropy, such as the carbonyl oxygen atom lone pair directionality or the sigma hole of heavier halogens.44 A distancedependent dielectric is used to account for the highly charged environment typically found in proteins. The EC values can be plotted on the solvent-accessible surface of the ligand, giving an intuitive visualization of where the ligand is well-matched to the protein and where possibilities exist for improvement. This can be used to guide ligand optimization, especially in later stages where small tweaks are desired to fine-tune affinity. In order to validate the method, we apply it to a range of literature data sets characterized by electrostatically driven SAR, including many fragment-based drug discovery studies such as for the PPI target XIAP or the GPCR mGLU5, 45,46

to test whether our EC calculations correlate with and can predict reported bioactivity

differences. In further examples, we apply EC analysis to kinase selectivity prediction and streptavidin mutant analysis. Our results show that EC assessment can be a powerful tool for analysis and optimization of electrostatic protein-ligand interactions, making it possible to quantify electrostatically driven SAR and predict electrostatic target selectivity. MATERIALS AND METHODS ESP calculation


5


Page 6 of 66

The base of electrostatic potential calculations is the polarizable XED force-field (version 3).42,43 Polarization effects and description of atomic charge anisotropy are especially useful for computing electrostatic properties of aromatic or unsaturated hydrocarbons, sp2-hybridized oxygen atoms, nitrogen atoms (sp, sp2, and sp3), and halogens (sigma hole of Cl, Br, and I).42,43,47 As previously described,44 the electrostatic interaction (potential) energy Vc for a probe p can be calculated using equation 1: 𝑖

(1) 𝑉𝑐 =

∑ 1

1 𝑞𝑝𝑞𝑖 332.17 4𝜋𝜀0 𝑟𝑖,𝑝 𝐷(𝑟)

where i is the number of point charges, ε0 the electrical field constant (permittivity in vacuum), q the signed magnitude of charges, ri,p the distance between probe p and charge i, and D(r) the following distance dependent dielectric function described by Mehler et al.48:

(2) 𝐷(𝑟) = 𝐴 +

𝐷0 ― 𝐴 (

)

1 + 𝑘𝑒 ―𝜆 𝐷0 ― 𝐴 𝑟

Parameters in equation 2 are: A = -8.5525, k = 7.7839, λ = 0.003627 and D0 = 78.4 (water relative permittivity at 25 ˚C).48 For ESP calculation protein or ligand molecules are selfpolarized, but not with respect to each other. EC calculation The protein-ligand EC is calculated from comparison of protein and ligand ESP values at all vertices of a generated ligand or protein solvent accessible surface (SAS).49,50 Perfect electrostatic complementarity would mean that at each vertex the ligand ESP value would be paired with a protein ESP value of the same magnitude with opposite sign. Figure 1 visualizes


6

Page 7 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ESP and EC of the biotin-streptavidin complex on both ligand and protein SAS, showing how a close matching of positive and negative electrostatic potential areas of protein and ligand leads to a good EC. As previously described by Dean et al.,18 one option to assess EC is to employ the Pearson’s R correlation test on the raw protein and ligand ESP values, yielding a coefficient of 1 for completely complementary and -1 for anti-complementary protein-ligand complexes:

(3) 𝐸𝐶𝑅 = 𝑟𝐸𝑆𝑃𝐿,

―1 𝐸𝑆𝑃𝑃 = 𝑛―1

𝑛

𝐸𝑆𝑃𝐿𝑖 ― 𝐸𝑆𝑃𝐿

∑(

𝑖=1

𝑠𝐸𝑆𝑃𝐿

𝐸𝑆𝑃𝑃𝑖 ― 𝐸𝑆𝑃𝑃

) ( ×

𝑠𝐸𝑆𝑃𝑃

)

With n being the number of vertices on the ligand surface, s the standard deviation, ESPL and ESPP the ligand and protein ESP values, ESPthe mean value of all protein or all ligand ESP values. The Pearson’s R coefficient is a global property calculated over the whole ligand surface, which is especially robust against strong electrical background fields. However, it is not well suited for scoring and coloring of small surface areas as a meaningful correlation test can only be performed if the respective area contains sufficiently large differences in both protein and ligand ESP values. Alternatively, we propose a score that provides an approximate correction for some desolvation effects and allows local visualization of EC on a protein or ligand solvent-accessible surface:

(4) 𝐸𝐶 =

∬(

1―

𝑆

𝐸𝑆𝑃𝐿 + 𝐸𝑆𝑃𝑃

)

𝑚𝑎𝑥(𝐸𝑆𝑃𝐿, 𝐸𝑆𝑃𝑃,𝑘)

𝑑𝑆

where the integral is over the ligand SAS, ESPL and ESPP the ligand and protein ESP values, and max(ESPL,ESPP,k) the protein or ligand ESP value with the largest deviation from zero, or a constant k if that is larger. A k value of 5 was chosen heuristically. Both ESPL and ESPP were capped to a maximum deviation from zero of 12 for the EC score. The capping value was


7


Page 8 of 66

derived from XED ESP computation of a water molecule and taking the 95th percentile of its ESP values. The rationale behind this approach is that any ligand comprising a charge/ESP area that is higher than the maximum ESP for water cannot be perfectly electrostatically complementary to solvent anymore and also pays a desolvation penalty upon complex formation. The ESP cap accentuates electrostatic clashes with at least one large ESP value, leading to a decreased EC score because of the desolvation contribution, and rewards the complementarity of large ESP values, which lead to large, attractive electrostatic interaction energies even if they are not identical in absolute values. EC scores range from 1 (perfect EC) to -1 (complete electrostatic clash). As solvent-exposed portions of the ligand contribute less information about the electrostatic complementarity of protein-ligand complexes, regions of the ligand SAS that are more than 3Å away from any protein atom are scaled down by a distance-dependent factor. This factor downweighs the respective EC value by the quotient of 3 divided by distance in angstrom (e.g., at a distance of 6Å the local EC value would be reduced to half of its original value). Unlike the ECR score, which is a statistical property of the entire surface, the EC score can be decomposed into perregion contributions. This allows coloring of the solvent-accessible surface by the EC contribution. We adopt a convention where regions with positive EC are colored green, regions with low or zero EC are colored white, and regions with negative EC are colored red.


8

Page 9 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Electrostatic complementarity (EC) of biotin-streptavidin complex (PDB 3RY2). The electrostatic potentials (ESP) of protein and ligand (blue = negative and red = positive ESP) and the protein-ligand EC (green = complementary, red = electrostatic clash) are projected on the protein (left panel) or ligand (right panel) solvent-accessible surface (SAS), respectively. Data sets and experimental setup An overview of data sets extracted from literature is given in Table S1. Reported SAR data (e.g., fragment-based drug discovery SAR as listed in fragment-to-lead medicinal chemistry reviews45,46) were pre-screened using the following criteria: (1) availability of at least one highquality X-ray structure that enables modelling of compound binding modes, (2) minimum 1 log


9


Page 10 of 66

unit difference in compound bioactivities, (3) neglect of SAR that is mainly determined by hydrophobic space filling in binding pocket or conformational effects, (4) absence of significant changes in binding mode due to ligand reorientation or protein flexibility. Protein-ligand complex structures were prepared with BuildModel51,52 in the structure-based modelling suite Flare53 (version 2.0) using default settings. Binding sites were visually inspected to correct protonation states of ligands and amino acid side chains and optimize water hydrogen bonding networks where necessary. As relative orientations and positions of water molecules are mostly dynamic, we chose to keep only water molecules in or close to the binding site that have at least 2 hydrogen bonding contacts to the protein or at least 1 hydrogen bond to both ligand and protein for EC calculations. If bridging water hydrogen bonding interactions were only possible for a subset of the reported binders, an alternative receptor without the bridging water was prepared for compounds that cannot interact with the respective bridging water molecule. Ligands were either prepared by manually editing X-ray binding complexes in Flare (including soft local minimization of edited moieties as required) or aligned to the original X-ray ligand using molecular field and shape-guided substructure alignment in Forge.44,54 If orientation of substituents was not clear from the reference X-ray binding mode (e.g. rotamers of ortho or meta substituents on aromatic rings), the orientation with the higher EC value was chosen. More details of the respective ligand and protein preparation procedure for each data set can be found in the supporting information and Table S1. The final coordinates of ligand and protein structures used for EC calculations can be found as associated content. All pictures were rendered in Flare 2.0.53 RESULTS


10

Page 11 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


To test our EC method, we applied it to several data sets that contain electrostatically driven SAR. Table 1 summarizes the coefficients obtained for diverse targets using EC and ECR scoring. We present two scores as the ECR score has been used by previous researchers,18–20 while the EC score allows a direct validation of the coloring the surface by electrostatic complementarity. In the following sections, we will discuss results and findings for each target in detail. Table 1. Correlation of EC scores with bioactivity

Target

Receptor

Data set size

Affinity range

mGLU5

5CGC, 5CGD

22

4 – 9.3 (pKi)

XIAP

5C7D

11

XIAP

5C7Ab

PIM1

5VUC

Mcl-1

6B4La

Mcl-1

6B4La

17d

MEK1

3DY7

imatinib-kinase biotin-SA

R2 coefficient

3.0 – 6.8 (pIC50)

EC 0.68, 0.61a 0.50

ECR 0.67, 0.52a 0.60

11

3.0 – 6.8 (pIC50)

0.70

0.65

6

6.3 – 8.7 (pIC50)

0.68

0.84

40

)c

0.18

0.16

)c

0.21

0.36

4

5.6 – 7.3 (pKD 6.9 – 8.7 (pIC50)

0.99

0.94

Various

9e

4.4 – 9.1 (pKD)

0.32

0.52

3RY2

4e

5.7 – 13.3 (pKD)

0.96

0.86

4.2 – 7.3 (pKD

a

Forge substructure alignment using the respective X-ray ligands as reference, b Relaxation of Lys297 side chain for each ligand (see text), c pKD values calculated from ΔG values at T=298K using ΔG° = -RT ln K, d Chlorine scan ligand subset, e Number of different receptor structures for same ligand.

mGLU5 Starting from a fragment hit against a thermostabilized metabotropic glutamate mGLU5 GPCR, Christopher et al. developed negative allosteric modulator lead compounds using X-ray crystallography and structure-based design.55 X-ray structures of the most active compounds (1


11


Page 12 of 66

and 2) in complex with mGLU5 served as starting points for model building (Figure 2). The published SAR comprises 22 compounds and contains several changes in bioactivity that appear to be mainly driven by changes in ligand electrostatics.

Figure 2. Structures of mGLU5 inhibitors and X-ray binding modes of compounds 1 and 2. Both 1 (left panel; PDB 5CGC) and 2 (right panel; PDB 5CGD) are negative allosteric modulators of mGLU5 and bind at the mavoglurant binding site, which is formed by residues of TM2, TM3, TM5, TM6, and TM7. The largest changes in pIC50 values were observed upon introduction of electron donating or withdrawing groups onto the ligand phenyl moiety, which is in close proximity to the π-electron surface of the Trp785 side chain and the thioether group of Met802. Figure 3 (and Figure S1 in


12

Page 13 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


an alternative view) visualize the electrostatic potential and complementarity for reported mGLU5 binders with different electrostatic properties at the ligand phenyl group. Electron donating groups such as methoxy (compound 3) increase the negative potential of the π-cloud and lead to a stronger electrostatic clash between the two aromatic ring planes. Additionally, the close proximity between the oxygen of the methoxy group and the Met802 sulfur atom (4.3Å) is also electrostatically unfavorable (Figure S1). Successive introduction of strongly electron withdrawing groups such as Cl, F, and nitrile in compounds 5, 6, and 1 gradually decrease the negative potential of the ligand phenyl ring (Figure 3) and consequently minimize the electrostatic clash between the phenyl ring and the indole ring while preserving the electrostatic complementarity to the rest of the binding site. The electrostatic changes between compound 5 and 6 are relatively small, and do not fully explain their large pIC50 difference. The compounds only differ in the chloro substituent of 6, which fills a hydrophobic subpocket of the binding site. It is worth noting that compounds 1a and 1 likewise differ only in the addition of a single chloro group, and have a similar activity delta to 5 and 6. Bioactivity changes caused by space filling and shape complementarity are not included in the EC calculations, which are size independent. As a result, compounds 5 and particularly 1a fall away from the regression lines in Figure 3.


13


Page 14 of 66

Figure 3. Comparison of ESP and EC surfaces of 5 representative mGLU5 ligands with increasing pIC50 values and correlation of mGLU5 bioactivities with EC scores. Positive and negative ESP values are colored in red and blue and electrostatically complementary and clashing regions are colored in green and red, respectively. The ESP ligand surfaces in the top panel show that the negative potential (blue color) on the face of the phenyl ring gradually decreases with introduction of electron withdrawing substituents. The EC surface in the middle


14

Page 15 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


panel shows that the electrostatic clash (red color) between the phenyl ring and the Trp785 indole (about 4Å distance to ligand phenyl ring) is minimized with decreasing negative electrostatic potential of the phenyl ring, which is in line with the gradually improving pIC50 values. The visually increasing electrostatic complementarity is also reflected in the improving EC scores. The correlation plots show that both EC and ECR scores correlate well with pIC50 values of the whole data set. Outliers 1a and 10 are highlighted and discussed in the text. Fluorination of the pyrazole moiety at position 4 of compound 6 lead to a seven-fold decrease in bioactivity (compound 7; Figure 4A). The fluorine substituent is in close proximity to the side chain oxygen atom of Ser658 and the nitrogen atom of the Ile625 backbone amide (about 3Å). The backbone carbonyl groups of Ile621 and Ser654 are also close to the fluorine (about 4Å). This electrostatically unfavorable alignment is reflected by the red ligand EC surface and was quantified by the respective EC scores (Figure 4A). Conversely, fluorination of the pyridine-2-yl substituent at position 5 of 8 increased the bioactivity by a factor of eight (compound 2). Here, the fluorine substituent reduces the electrostatic clash of the pyridine π-system with the Tyr659 phenol ring and the Ser654 carbonyl group and improves the electrostatic interactions with the Cα proton of Ala810 and the methyl side chain of Ala813 (Figure 4B). These two ligand pairs demonstrate nicely that similar chemical transformations can have very different impacts on EC and ligand bioactivities, strongly depending on their relative orientation and alignment within the protein pocket.


15


Page 16 of 66

Figure 4. Comparison of ESP and EC surfaces of compound 6 (5CGC model) and 8 (5CGD model) with their fluorinated analogs. Positive and negative ESP values are colored in red and blue and electrostatically complementary and clashing regions are colored in green and red, respectively. (A) Upon 4-fluorination of the pyrazole ring, the negative electrostatic potential of the fluorine substituent electrostatically clashes with oxygen and nitrogen atoms of the surrounding residues. (B) Fluorination of the 5 position of the pyridine ring of compound 8 decreases the electrostatic clash with Ser654 and Tyr659 and improves the complementarity with Ala810 and Ala813. Replacing the nitrogen atoms in the pyrimidine ring of 9 with carbon atoms to give 10 and 11 lead to an approximately seven-fold decrease in IC50 values in both cases (Figure S2). Interestingly, only compound 10 showed a significantly decreased EC, mainly due to the loss of a weak hydrogen bonding interaction between the heterocyclic nitrogen at position 3 of the


16

Page 17 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


scaffold and the polarized Cα protons of Gly628 and methylene side chain proton of Ser654 (Figure S2). Compound 10 also appears to be an outlier in the EC-bioactivity correlation plots (Figure 2). The reference compound 1, which was used for modelling the mGLU5 SAR, contains a nitrogen atom at the modified ring position of 10 that may induce binding site conformations that enable weak hydrogen bonding with the ligand. Upon binding of 10 the polarized proton positions of Gly628 and Ser654 may adapt to decrease the electrostatic clash, giving a potential explanation for the relatively low EC scores of 10 when using the reference binding site conformation. Replacing the other nitrogen atom with carbon (compound 11) did not strongly change the EC, suggesting that other factors such as conformational preferences (e.g. nitrogen is at ortho position between two aromatic 6 rings) may contribute to the observed bioactivity change. Overall, the reported pIC50 values for the mGLU5 compound data set show a good correlation with computed EC scores (Figure 2). Although the correlation seems to rely on very few data points at lower pIC50 values (e.g., 3 and 4), there is still a reasonable EC-bioactivity correlation for data points above a pIC50 of 7 (R2 values of 0.42 and 0.46 for EC and ECR, respectively), indicating that EC scoring is still predictive for this narrower but more densely populated bioactivity space. Using field and shape-guided substructure alignment as implemented in Forge45 with the 6B4L ligand as reference yielded similar results to manual model building starting from the X-ray ligand, although with slightly reduced correlation coefficients (see Table 1).

XIAP


17


Page 18 of 66

Chessari et al. recently reported small molecule inhibitors of the XIAP / c-IAP-caspase proteinprotein interaction.41 Both proteins are members of the inhibitor of apoptosis protein (IAP) family, which are frequently deregulated in cancer causing tumor growth and treatment resistance and thus constitute promising anticancer targets. The researchers used fragment-based drug discovery and structure-based design to develop low nanomolar lead compounds. A particularly interesting feature of this data set was the strong correlation of compound potency with Hammett σp values of substituents on the aromatic part of the indoline ring. This finding was rationalized with the electrostatic clash between the aromatic indoline system with a backbone carbonyl group and a phenolic tyrosine oxygen using DFT ab initio calculations and XED-based molecular fields.44


18

Page 19 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. Structures of XIAP inhibitors, binding mode of inhibitors 12 and 13 bound to XIAP (PDB 5C7A and 5C7D) and their EC projected on the protein SAS. Lead compound 13, which is about 65-fold more potent than 12, shows an improved EC, especially in proximity to the Lys297 side chain and the Gly306 backbone carbonyl group. Binding modes of the starting compound 12 and an optimized lead compound 13 are shown in Figure 5. The depicted EC maps show two areas of improved complementarity: (1) the area around the positively charged Lys297 side chain and (2) the area close to the carbonyl group of


19


Page 20 of 66

the Gly306 backbone. The X-ray structure (PDB 5C7D) of the optimized reference compound was used for modelling the data set. Its Lys297 side chain, which adapted its conformation to the chloro substituent of 13, can accommodate most of the ligand substituents without any clashes. However, for compound 14, the Lys297 side chain was optimized with the XED force field to avoid a clash with the methylsulfonyl substituent. With increasing pIC50 values, the XIAP inhibitors tend to have more electron-withdrawing substituents with negative potential in proximity to the Lys297 side chain (Figure 6). Compound 15, a 6-amino analog of 12, shows the worst EC because of the positive ESP and electrondonating properties of its amino group. Conversely, compound 13 shows good EC to Lys297 and a decreased electrostatic clash with Gly306 due to its negative, electron-withdrawing chloro substituent as well as the 5-nitrogen heteroatom that additionally decreases the negative potential of the π-system. Comparison of compounds 16 and 17, which possess similar EC scores, reveals that the space filling effect of the larger chlorine substituent of 17 also contributes to increased bioactivity. For the data set as a whole, there was generally a good correlation between pIC50 values and EC scores (Figure 6), demonstrating that ligand electrostatics play an essential role for XIAP-ligand binding. However, the most potent compound of the data set, 14, appears to be an outlier particularly in the EC score correlation plot. This is because the methylsulfonyl substituent of 14 is significantly larger than the corresponding substituent of the other ligands. As seen in the mGLU5 data set, the EC calculations do not take into account changes in the overall size of the ligand, and hence are generally only quantitative for ligands of similar heavy atom count. In addition, the protein conformation had to be modified to accommodate this substituent, and any energetic effects of this were not included in the EC calculation.


20

Page 21 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


An additional experiment to analyze the data set in a more prospective manner was performed by modelling all ligands using just the X-ray structure of XIAP in complex with the original fragment 12 (PDB 5C7A). For each modelled ligand, we allowed relaxation of Lys297 side chain atoms to avoid steric clashes and optimize protein-ligand interactions. EC scores were then calculated using the respective optimized protein structure for each ligand. This protocol yielded even better correlation between EC and activity (R2 = 0.70 and 0.65 for EC and ECR scores, respectively, data not shown).


21


Page 22 of 66

Figure 6. Electrostatic complementarity of XIAP inhibitors. With increasing bioactivity, the degree of electrostatic clash (red areas) with Gly306 carbonyl and Lys297 side chain, as exemplified by compound 12 and 15, is lessened or even turned into electrostatic complementarity (green areas). Electron-withdrawing 6-substituents decrease the negative potential of the aromatic π-plane and change the electrostatic potential at this position from positive (15, amine group) to a negative potential (17, chloro substituent).

PIM1 Recently, Watanabe et al. described a data set of serine/threonine PIM1 kinase benzofuranone inhibitors, which show large affinity differences (225-fold from best to worst) upon introduction of a single nitrogen atom into an aromatic indole system at different positions.56 These activity cliffs were studied in depth using the fragment molecular orbital method with molecular mechanics Poisson−Boltzmann surface area (FMO+MM-PBSA) approach. Using QM/MMoptimized structures, the authors reported an excellent correlation (R2 = 0.85) between experimental pIC50 values and the calculated FMO+MM-PBSA energies, which include MM correction terms for ligand strain (deformation energy upon binding compared to lowest energy conformer in solution) and desolvation free energy calculated from Poisson-Boltzmann electrostatics. The X-ray binding mode of reference compound 18 and its EC to the receptor is shown in Figure 7A. Interestingly, the EC surface highlights that there are electrostatic clashes in proximity to position 2 and 7 of the indole scaffold which may be amenable to further optimization. Aza-analogs of 18 and their respective bioactivities are shown in Figure 7B. EC scoring of this data set against the 5VUC receptor yielded good correlations with the reported


22

Page 23 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


pIC50 values (Fig. 7C), agreeing with Watanabe et al.’s results but using a much simpler protocol and orders of magnitude less calculation time.

Figure 7. Binding mode of benzofuranone inhibitor 18 of PIM1 and bioactivities of derived aza analogs. (A) X-ray structure of 18 in complex with PIM1 (PDB 5VUC) and its electrostatic complementarity to the binding pocket. (B) Chemical structures and respective bioactivities of aza analogs of 18. (C) Correlation of pIC50 values with EC scores.


23


Page 24 of 66

The ligand ESP plots show that introduction of nitrogen atoms into the indole system slightly reduces the negative potential of the π-plane, which is in close contact with alkyl side chains of Leu44, Val52, Ala65, Leu174, and Ile185 to form CH-π interactions (Figure S3). A decreased negative potential of the azaindole π-plane should also decrease the electrostatic interactions with the surrounding alkyl side chain protons, however, both compounds 19 and 20 show similar activity to 18. These compounds compensate the weaker CH-π interactions with improved electrostatics in the previously highlighted areas with electrostatic clashes, which are in proximity to position 2 and 7 of the indole (Figure 7). Compounds 21-23 show relatively similar ECR scores, although 23 is about 7-10 fold less active than 21 and 22. A contributing factor could be the internal electrostatic clash between the introduced nitrogen atom and the oxygen of the benzofuran-3(2H)-one system (about 3Å distance), which may impact the stability of the binding confirmation of 23 and therefore lead to increased strain energy upon binding. The EC scores for this data set are generally noisier and do not distinguish as clearly between the two groups, which is also reflected by the lower R2 coefficient.

Mcl-1 Friberg et al. reported the development of selective myeloid cell leukemia 1 (Mcl-1) inhibitors using structure-based design and a fragment-merging approach.57 Mcl-1 is a member of the Bcl-2 protein family and also a promising anticancer target, as Mcl-1 overexpression prevents cancer cell apoptosis. The reported Mcl-1 SAR was recently compiled for a free energy benchmark data set and tested with free binding energy simulations that correlated well with compound activities (R2 = 0.60).58 Ligand structures were taken from the benchmark data set and were aligned to the


24

Page 25 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


co-crystallized ligand (PDB 6B4L) in Forge,54 and the Build Model-prepared51,52 receptor was then used to perform EC calculations. A detailed look at the ligand SAR showed that the majority of ligand activity variations are due to the size and shape of a hydrophobic anchor moiety (e.g. naphth-1-yl with oxypropyl linker), which fills a hydrophobic pocket (Figure 8).

Figure 8. Upper and left panel: Mcl-1 inhibitor structures and binding mode of Mcl-1 inhibitor 24 (PDB 6B4L). Right panel: a subset of the compound data shows systematic variation of chlorine substitution on the core scaffold.

The EC scores for the full data set showed low correlation with the reported bioactivity (R2 values of 0.18 and 0.16). This outcome is not surprising, as EC metrics are unable to capture activity differences based on hydrophobic space-filling and shape-based differences. To focus on


25


Page 26 of 66

electrostatic differences between ligands, we compiled a subset of all compounds that contain a chlorine substituent on the core scaffold. This experiment should test whether EC maps and scores can predict the most favorable position for chloro substitution. Where available, the corresponding des-chloro compound was added to provide a baseline. This ‘chlorine scan’ subset showed improved correlations of EC scores with the binding free energies (R2 values of 0.21 and 0.36 for EC and ECR) and ranked the most favorable chlorine substitution among the top EC scores. EC analysis of the chlorine scan subset is illustrated in Figure 9. Compounds with 6chloro substitution (e.g., 29) are clearly the most potent binders in this subset and show also the highest EC scores. Compound 29 shows a mostly green to white EC map on the chlorine atom, which can be rationalized by the presence of a favorable halogen bond between the 6-chlorine and the Ala227 backbone carbonyl group and polar interactions with polarized protons of surrounding residue side chains. Chloro substitution at position 5 (e.g., 28) shows the least favorable EC as indicated by the red map in proximity to the side chain of Phe270 (Figure 9). This arrangement leads to an electrostatic clash of the chlorine atom with the π-plane of the side chain phenyl ring Phe270. The EC results for 28 also agree well with its high ΔG value. The 7chloro substituted 30 and its corresponding des-chloro analog 26 showed very similar EC maps and scores, although 30 was slightly more active by a factor of 4. The larger ligand-protein contact area for 30 may be a contributing factor. Chloro substitution at position 4 of compounds 24 and 31 either decreased or increased ΔG values, depending on the hydrophobic anchor: for the 1-naphthyl moiety it decreased binding by 5-fold (33), whereas for the 3,5-dimethyl-4chlorophenyl moiety it marginally increased ΔG by 3-fold (32; see Figure S4). As the EC scores slightly decreased in both cases, conformational effects of the 4-chloro substitution may be the key factor for the observed ΔG change. Interestingly, for the equipotent compounds 24 and 31


26

Page 27 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the difference in EC maps and scores is not reflected in ΔG. The naphthyl anchor of 24 fills more hydrophobic volume of the pocket whereas the 3,5-dimethyl-4-chlorophenyl anchor of 31 shows a better EC, which suggests that the compensation of shape matching versus EC results in similar bioactivities here. Taken together, our results demonstrate that the EC surface maps in this system do provide valuable guidance about the electrostatic fit of the ligand to the receptor and can be quantitatively related to activity.


27


Page 28 of 66

Figure 9. Correlation of EC scores with bioactivities for the Mcl-1 chlorine scan subset. EC maps show that 6-chloro substitution (29) shows the most favorable electrostatics due to formation of a halogen bond with the Ala227 backbone carbonyl group. In both EC and ECR plots, compounds with 6-Cl substitution (data points shown as black ■) are enriched at high EC scores.


28

Page 29 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Halogen Bonding Halogen bonds are mainly electrostatic interactions. The so-called sigma hole, an electron deficient area on the tip of the heavier halogens (Cl, Br, I), can interact with electron-donating Lewis bases (e.g. backbone carbonyl atoms) to form halogen bonds. Due to the anisotropic halogen electron distribution, the directional nature of this interaction cannot be properly described with atom-centered charges. Wilcken et al. investigated this interaction type in detail using ab initio calculations and provided literature examples that showed increasing halogen bonding strength with increasing halogen weight and sigma hole size (from Cl to I).1 However, the different halogen substituents not only differ in the electrostatic properties of the halogen itself, but also change the electrostatics of aromatic scaffolds (which can alter binding affinity by tuning other interactions such as CH-π) and have different vdW volumes in the pocket (which can contribute to the measured binding affinity via dispersion effects and entropic contributions). We used EC to analyze a series of MEK1 ligands.59 This example shows the increasing electrostatic complementarity of the halogen sigma hole with the carbonyl atom (Figure 10) in both maps and scores. Replacing the heavier halogens with fluorine, which does not typically have a positive sigma hole that allows formation of halogen bonds, resulted in an electrostatic clash with the carbonyl oxygen atom. The correlation of EC scores with experimental pIC50 values yielded R2 values of 0.99 and 0.94 for EC and ECR, respectively.


29


Page 30 of 66

Figure 10. EC analysis of MEK1 inhibitors that form halogen bond to the backbone (BB) carbonyl group of Val127. The depicted structures show the complex structure of the iodine compound or modelled analogs with MEK1 (PDB 3DY7). The electrostatic complementarity


30

Page 31 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


decreases (both EC maps and scores) from iodine to fluorine (no halogen bonding interaction), as expected from theoretical halogen bonding strengths.1

Correlation of EC values to changes in protein sequence The anticancer drug imatinib inhibits tyrosine kinases such as ABL/BCR-ABL, c-Kit, and PDGF-R and is mainly used for treating chronic myelogenous leukemia and acute lymphocytic leukemia. Due to the high conservation of the ATP binding site within the kinase family, kinase inhibitors such as imatinib are usually not completely selective for their targets and can inhibit other kinases to a varying degree, which poses a major challenge for the development of kinase inhibitor drugs. As there is extensive selectivity and structural data available for imatinib,60,61 we investigated whether differences in the EC of imatinib-kinase complexes could be an important factor for imatinib kinase selectivity. The data set consists of all imatinib-kinase complex structures that were available in the Protein Data Bank (September 2018, Figure 11A).62 Structures for EC experiments were prepared as described in the Methods section.


31


Page 32 of 66

Figure 11. Electrostatic complementarity analysis of imatinib kinase selectivity. (A) Imatinib structure and its respective pKD values for imatinib-kinase complexes available in the PDB. (B) Correlation of EC with imatinib pKD values for kinase subset. The c-SRC data point was not included in the Pearson’s R2 correlation scores as weak imatinib binding to c-SRC is mainly due


32

Page 33 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


to an energetic penalty of the imatinib bioactive conformation.63 (C) EC maps highlight differences between imatinib binding to ABL and the ancestral ANC-AS kinases.

Apart from the SYK kinase (PDB 1XBB), where imatinib binds in a compact cis-conformation,64 imatinib is bound in similar extended conformations to the respective kinase ATP pockets. There was a reasonable correlation between experimental pKD values and EC scores, suggesting that electrostatic features of kinase binding pockets and their spatial arrangement around imatinib can be an important factor for kinase selectivity (Figure 11B). Comparison of imatinib-kinase complexes with good EC (ABL) and bad EC (ancestral SRC/ABL kinase, ANC-AS) revealed several factors that decrease the EC of imatinib in complex with ANC-AS (see Figure 11C): (1) electrostatic clash of the 3-pyridyl moiety with the polarized Cα and backbone NH protons (Gly87 and Ser88), (2) missing CH-π interactions of the 4-(3-pyridyl)pyrimidin-2-yl substructure with a tyrosine side chain (Tyr272 in ABL), (3) position of water network around Lys37, Glu52, Asp147, and Phe148 is closer to the toluene substructure causing an electrostatic clash, and (4) a less favorable EC of the piperazine ring with the ANC-AS receptor due to a different orientation and conformation. The EC-pKD correlation plots show that both the proto-oncogene tyrosineprotein kinase Src (c-SRC) and the mitogen-activated protein kinase 14 (p38a) are outliers. Seeliger et al. reported that imatinib binds to an inactive Abl/c-Kit-like conformation of c-SRC, which leads to a distributed thermodynamic penalty due to conformational effects outside the ATP binding pocket upon imatinib binding. This finding rationalizes the low pKD value despite good EC with the c-SRC binding pocket. Additionally, Namboodiri et al. reported for the imatinib-p38a and c-SRC complexes that their higher solvent accessible surface areas may contribute to the relatively weak binding of imatinib.65 When excluding both c-SRC and p38a


33


Page 34 of 66

from the data set, the R2 correlation coefficients for EC and ECR scores further improve to 0.77 and 0.71, respectively. These results show that EC analysis can reveal electrostatic features that are important for binding selectivity and score their contributions accordingly. The biotin-streptavidin complex shows one of the strongest non-covalent affinities known to date and has been subject of many studies. Its high degree of EC has already been shown in Figure 1 in the Methods section. Hyre et al. tested the effect of S45A, D128A, and S45A/D128A mutants, which form hydrogen bonds with the 2-oxoimidazolidine moiety of biotin, to investigate the thermodynamic and structural cooperativity of biotin binding.66 The authors reported that the single mutants show a 1000-fold and the double-mutant a 2107-fold reduced binding, indicating that the two mutations interact cooperatively. Starting from a high-resolution wild-type biotinstreptavidin X-ray structure (PDB 3RY2), we introduced the respective alanine mutations and applied EC analysis. The maps clearly show the gradual decrease in EC upon mutation of the binding site and the EC scores correlate with the determined binding affinities (Figure 12). This example highlights that the electrostatic analysis of protein mutations on protein-ligand binding is feasible and predictive.


34

Page 35 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 12. EC analysis of wild-type and mutant biotin-streptavidin complexes. EC scores can rank the experimentally determined binding affinities. The ligand EC maps highlight the decreased electrostatic complementarity in proximity to S45A and D128A mutations due to loss of hydrogen bonding acceptor motifs.


35


Page 36 of 66

DISCUSSION It is well known that hydrophobic, entropically-driven protein-ligand interactions can significantly increase ligand binding affinity.67 However, analysis of a large number of highquality protein-ligand X-ray structures with binding data (taken from Binding MOAD 68,69) showed that there is no clear correlation between ligand size or protein-ligand contact area and binding affinity.70 Different conformational, steric, and entropic contributions to the free energy of binding may explain this lack of correlation; however, the highest ligand efficiencies of the data set were interestingly observed for ligands that bind to highly charged protein pockets, comprising either several charged protein side chains or metal ions. Smith et al. concluded that electrostatic interactions define the maximum ligand efficiency, which highlights the importance of assessing electrostatic interactions and complementarity of protein-ligand complexes more routinely during drug development.70 The electrostatic complementarity methodology presented here is a powerful and versatile tool for drug design. In particular, the visualization of the EC values across the surface of a ligand and/or protein is a powerful aide to lead optimization. It highlights which parts of the ligandprotein interface are electrostatically sub-optimal and provides a simple way of assessing whether proposed changes to the structure are likely to improve binding. Such guidance is, of course, only useful if it actually does correlate with activity. In multiple test cases where the protein-ligand affinity differences could reasonably be expected to be driven by electrostatics the EC scores provided a good correlation with activity. The method is computationally inexpensive, allowing it to be applied interactively. EC scoring of hundreds of molecules can be done within minutes (less than 60s for 100 molecules) on a


36

Page 37 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


desktop computer (see Table S2 for detailed benchmark results). The presented EC method can easily be combined with experiments such as molecular docking or field and shape-based bioisosteric replacement or screening procedures, provided structural information is available, to enrich molecules with high electrostatic complementarity to the respective receptor. In context of virtual screening (VS) especially the size-independent property of EC scoring has to be considered, as it neglects the varying degrees of ligand size and the resulting protein-ligand contact surface area that also contribute to the binding free energy. The use of ‘normalized’ EC scores using descriptors for ligand space filling or buried protein-ligand surface area may be an avenue to improve ranking of compound bioactivity. However, the potential benefit of EC scoring for virtual screening (VS) has to be still investigated in detail. Target classes like PPIs (e.g., XIAP) may be particularly well suited for EC analysis, due to the intrinsically high EC of protein-protein interaction complexes71 and less pronounced contributions from factors such as shape effects, hydrophobic space filling, and ligand desolvation to the binding free energy. It is worth noting that EC analysis alone is not going to provide a complete prediction of the binding free energies of protein-ligand complexes due to missing key factors such as conformational contributions for both proteins and ligands, desolvation effects, van der Waals and entropic contributions, and energetics of binding site water molecules. However, it is very common during the lead optimization process to want to improve the electrostatic interactions between the ligand and the protein by making relatively small changes. The ability to visualize the EC allows the identification of the regions of the ligand where the electrostatics are suboptimal, and the EC score in this situation provides excellent correlations with activity. It may be possible to improve predictions further by combining the EC scores with computed estimates of ligand strain energy calculations,72,73 correction for desolvation contributions,4,73,74 and EC score


37


Page 38 of 66

normalization with protein-ligand interaction surface area to account for vdW contacts (EC scores are independent from size of ligands and vdW contact area). In this work the well-defined ligand SAS has been used for EC computation. Alternatives, such as the buried protein-ligand surface, may help to minimize the contribution of unwanted solvent exposed ligand areas and give a better estimate of vdW interactions. The quality of EC analysis strongly relies on a reasonable model of the correct ligand binding conformation, which preferably should be experimentally determined or modelled from a suitable structure of a reference protein-ligand complex. For our data sets, manual building from reference X-ray structures or field-based substructure alignment in Forge both provided good starting structures for EC calculations, with further guidance from the EC method itself to determine the best suited orientation of substituents in case experimental guidance is not available. Ligand docking is a further option to generate ligand binding poses, however, in our hands free docking caused some noise in the EC scores due to small changes in ligand scaffold placement of congeneric SAR (data not shown). Benchmarking electrostatic complementarity metrics is not trivial as the electrostatic potential and protein-ligand complementarity cannot be directly determined. Agreement of EC methods with experiments can only be indirectly assessed by analyzing electrostatically driven SAR. Unlike most scoring functions EC scores are independent of ligand size. As most reported SAR includes changes that are caused by previously mentioned factors (e.g. large differences in ligand size, conformational contributions, desolvation, hydrophobic space filling, reordering of water, potential changes in binding mode or bioactive conformation), we made an effort to compile focused SAR data sets where affinity changes were clearly electrostatically driven.


38

Page 39 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In this paper we presented two different EC scores: EC and ECR. The reason behind developing the EC score was the primary interest in visualization of protein-ligand EC, which requires a scoring metric that provides local values for surface coloring. However, the coloring should provide a meaningful assessment of EC and, thus, had to be validated against electrostaticallydriven bioactivity data. Previous EC studies employed the Pearson’s R correlation coefficient as EC metric. We included this approach (ECR) for direct comparison against earlier work. Its advantage for numeric scoring is its robustness in context of strong electric background fields. However, it cannot be locally displayed in a meaningful way due to its global nature of assessing the correlation between protein and ligand ESPs. Local ECR coloring would require to define surface areas with sufficiently large differences in protein and ligand ESP values to check for linear correlation of ESP data points, but such differences are not always available in small surface regions. The benchmark results of the EC scores show that the ECR score appears to correlate slightly better with experimental data than the EC score. For four of the data sets listed in Table 1 the ECR correlation coefficient is higher than for EC and the other six data sets show comparable results for EC and ECR (assuming a meaningful R2 difference must be at least 0.1). The imatinib selectivity correlation ECR also yields slightly better results than EC. A possible reason for this trend is that truncation of electrostatic potentials in EC scores, to correct for desolvation effects, may smooth strongly complementary or clashing areas too much. Alternatively, the improved performance of the ECR score may be due to its insensitivity to uniform electrostatic fields over the active site (e.g., where the protein has a large net charge). As the correlation coefficients have large error bars for small data sets (N < 25), more and larger data sets would be needed to confirm this trend with statistical significance. Although most of our analyzed data sets


39


Page 40 of 66

contained several formal charges within the binding pocket (see Table S1), which can be an issue for electrostatic calculations due to overemphasis of formal charges, EC analysis worked reliably for both uncharged (mGLU5) and highly-charged binding sites (XIAP).

CONCLUSION AND OUTLOOK Analysis of electrostatic complementarity is a powerful tool for drug design. It can be employed to highlight ligand substructures that can be electrostatically optimized to improve protein binding, rank ligand affinity changes due to differences in EC, and analyze selectivity changes for different proteins or protein mutations. The presented method is computationally inexpensive, can be applied to large data sets, and includes atomic electrostatic anisotropy and polarization effects thanks to its XED force-field foundations. As part of the validation of this method we compiled benchmark sets focused on electrostatically driven SAR and provide them as associated content. Future improvements of the presented method will include a more rigorous exclusion or down weighting of solvent exposed areas of ligands for EC scoring and improved description of desolvation effects.

ASSOCIATED CONTENT Supporting Information. The following files are available free of charge. Detailed protein and ligand preparation procedures for electrostatic benchmark sets, tables with


40

Page 41 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


detailed data set statistics and EC performance benchmark, additional EC analysis figures for several data sets (PDF) SD and PDB files of ligands and protein structures used for EC analysis (ZIP) AUTHOR INFORMATION Corresponding Author *E-mail: [email protected]. Phone: +44 (0)1223 858890

ORCID Matthias R. Bauer: https://orcid.org/0000-0003-4015-6483 Mark D. Mackey: https://orcid.org/0000-0001-5131-7583

Author Contributions The manuscript was written through contributions of all authors. M.R.B. and M.D.M. designed research; M.D.M developed and wrote code for EC method; M.R.B. performed research; M.R.B. and M.D.M. analyzed data; and M.R.B. and M.D.M. wrote the paper. All authors have given approval to the final version of the manuscript. ‡These authors contributed equally. Funding Sources

Notes


41


Page 42 of 66

The authors declare the following competing financial interest(s): Mark D. Mackey is a shareholder in Cresset.

ACKNOWLEDGMENT We acknowledge and thank Paolo Tosco, Giovanna Tedesco, Tim Cheeseright, and Andy Vinter for helpful discussions and comments on both the method and this paper. ABBREVIATIONS ABL, tyrosine kinase ABL1; c-Kit, tyrosine kinase KIT; c-SRC, proto-oncogene tyrosine-protein kinase Src; ESP, electrostatic potential; EC, electrostatic complementarity; FMO, fragment molecular orbital; Mcl-1, induced myeloid leukemia cell differentiation protein Mcl-1; MEK1, dual specificity mitogen-activated protein kinase kinase 1; mGLU5, metabotropic glutamate receptor 5; MM, molecular mechanics; p38a, mitogen-activated protein kinase 14; PDGF-R, Platelet-derived growth factor receptors; PIM1, proto-oncogene serine/threonine-protein kinase Pim-1; QM, quantum mechanical; SAS, solvent accessible surface; XIAP, X-linked inhibitor of apoptosis protein REFERENCES (1)

Wilcken, R.; Zimmermann, M. O.; Lange, A.; Joerger, A. C.; Boeckler, F. M. Principles

and Applications of Halogen Bonding in Medicinal Chemistry and Chemical Biology. J. Med. Chem. 2013, 56 (4), 1363–1388. https://doi.org/10.1021/jm3012068. (2)

Zhou, P.; Zou, J.; Tian, F.; Shang, Z. Fluorine Bonding — How Does It Work In

Protein−Ligand

Interactions?

J.

Chem.

Inf.

Model.

2009,

49

(10),

2344–2355.

https://doi.org/10.1021/ci9002393.


42

Page 43 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(3)

Meanwell, N. A. Fluorine and Fluorinated Motifs in the Design and Application of

Bioisosteres

for

Drug

Design.

J.

Med.

Chem.

2018,

61

(14),

5822–5880.

https://doi.org/10.1021/acs.jmedchem.7b01788. (4)

Bauer, M. R.; Jones, R. N.; Baud, M. G. J.; Wilcken, R.; Boeckler, F. M.; Fersht, A. R.;

Joerger, A. C.; Spencer, J. Harnessing Fluorine–Sulfur Contacts and Multipolar Interactions for the Design of P53 Mutant Y220C Rescue Drugs. ACS Chem. Biol. 2016, 11 (8), 2265–2274. https://doi.org/10.1021/acschembio.6b00315. (5)

Ma, J. C.; Dougherty, D. A. The Cation−π Interaction. Chem. Rev. 1997, 97 (5), 1303–

1324. https://doi.org/10.1021/cr9603744. (6)

Meyer, E. A.; Castellano, R. K.; Diederich, F. Interactions with Aromatic Rings in

Chemical and Biological Recognition. Angew. Chem. Int. Ed. Engl. 2003, 42 (11), 1210–1250. (7)

Hann, M. M. Molecular Obesity, Potency and Other Addictions in Drug Discovery. Med.

Chem. Commun. 2011, 2 (5), 349–355. https://doi.org/10.1039/C1MD00017A. (8)

Hann, M. M.; Keserü, G. M. Finding the Sweet Spot: The Role of Nature and Nurture in

Medicinal

Chemistry.

Nat.

Rev.

Drug

Discov.

2012,

11

(5),

355–365.

https://doi.org/10.1038/nrd3701. (9)

Scrocco, E.; Tomasi, J. The Electrostatic Molecular Potential as a Tool for the

Interpretation of Molecular Properties. New Concepts II 1973, 42, 95–170. (10) Politzer, P.; Laurence, P. R. Molecular Electrostatic Potentials: An Effective Tool for the Elucidation of Biochemical Phenomena. Environ. Health Perspect. 1985, 61, 191–202.


43


Page 44 of 66

(11) Tosco, P.; Mackey, M. Lessons and Successes in the Use of Molecular Fields. In Comprehensive Medicinal Chemistry III; Elsevier, 2017; pp 253–296. (12) Davis, A.; Warrington, B. H.; Vinter, J. G. Strategic Approaches to Drug Design. II. Modelling Studies on Phosphodiesterase Substrates and Inhibitors. J. Comput. Aided Mol. Des. 1987, 1 (2), 97–119. https://doi.org/10.1007/BF01676955. (13) Kukić, P.; Nielsen, J. E. Electrostatics in Proteins and Protein-Ligand Complexes. Future Med Chem 2010, 2 (4), 647–666. https://doi.org/10.4155/fmc.10.6. (14) Weiner, P. K.; Langridge, R.; Blaney, J. M.; Schaefer, R.; Kollman, P. A. Electrostatic Potential Molecular Surfaces. Proc. Natl. Acad. Sci. U.S.A. 1982, 79 (12), 3754–3758. https://doi.org/10.1073/pnas.79.12.3754. (15) Náray‐Szabó, G.; Nagy, P. Electrostatic Complementarity in Molecular Aggregates. 9: Protein‐ligand

Complexes.

Int.

J.

Quantum

Chem.

1989,

35

(1),

215–221.

https://doi.org/10.1002/qua.560350114. (16) Náray-Szabó, G. Electrostatic Complementarity in Molecular Associations. J. Mol. Graph. 1989, 7 (2), 76–81. https://doi.org/10.1016/S0263-7855(89)80003-7. (17) Nakamura, H.; Komatsu, K.; Nakagawa, S.; Umeyama, H. Visualization of Electrostatic Recognition by Enzymes for Their Ligands and Cofactors. J. Mol. Graph. 1985, 3 (1), 2–11. https://doi.org/10.1016/0263-7855(85)80007-2. (18) Chau, P.-L.; Dean, P. M. Electrostatic Complementarity between Proteins and Ligands. 1. Charge Disposition, Dielectric and Interface Effects. J. Comput. Aided Mol. Des. 1994, 8 (5), 513– 525. https://doi.org/10.1007/BF00123663.


44

Page 45 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(19) Chau, P.-L.; Dean, P. M. Electrostatic Complementarity between Proteins and Ligands. 3. Structural

Basis.

J.

Comput.

Aided

Mol.

Des.

1994,

8

(5),

545–564.

https://doi.org/10.1007/BF00123665. (20) Chau, P.-L.; Dean, P. M. Electrostatic Complementarity between Proteins and Ligands. 2. Ligand

Moieties.

J.

Comput.

Aided

Mol.

Des.

1994,

8

(5),

527–544.

https://doi.org/10.1007/BF00123664. (21) Kangas, E.; Tidor, B. Optimizing Electrostatic Affinity in Ligand–Receptor Binding: Theory, Computation, and Ligand Properties. J. Chem. Phys. 1998, 109 (17), 7522–7545. https://doi.org/10.1063/1.477375. (22) Kangas, E.; Tidor, B. Charge Optimization Leads to Favorable Electrostatic Binding Free Energy. Phys. Rev. E 1999, 59 (5), 5958–5961. https://doi.org/10.1103/PhysRevE.59.5958. (23) Kangas, E.; Tidor, B. Electrostatic Complementarity at Ligand Binding Sites: Application to

Chorismate

Mutase.

J.

Phys.

Chem.

B

2001,

105

(4),

880–888.

https://doi.org/10.1021/jp003449n. (24) Lee, L.-P.; Tidor, B. Optimization of Electrostatic Binding Free Energy. J. Chem. Phys. 1997, 106 (21), 8681–8690. https://doi.org/10.1063/1.473929. (25) Naray-Szabo, G. Quantitative Estimation of Electrostatic Complementarity. MATCH Commun. Math. Comput. Chem. 2001, 44 (44), 297–304. (26) Muzet, N.; Guillot, B.; Jelsch, C.; Howard, E.; Lecomte, C. Electrostatic Complementarity in an Aldose Reductase Complex from Ultra-High-Resolution Crystallography and First-


45


Page 46 of 66

Principles Calculations. Proc. Natl. Acad. Sci. U.S.A. 2003, 100 (15), 8742–8747. https://doi.org/10.1073/pnas.1432955100. (27) Huggins, D. J.; Sherman, W.; Tidor, B. Rational Approaches to Improving Selectivity in Drug Design. J. Med. Chem. 2012, 55 (4), 1424–1444. https://doi.org/10.1021/jm2010332. (28) Ansalone, P.; Bottauscio, O.; Manzin, A. Cation‐mediated Electrostatic Interaction in Collagen–Integrin Complex. Surf. Interface Anal. 2014, 46 (10–11), 693–697. (29) Voet, A.; Berenger, F.; Zhang, K. Y. J. Electrostatic Similarities between Protein and Small Molecule Ligands Facilitate the Design of Protein-Protein Interaction Inhibitors. PLOS ONE 2013, 8 (10), e75762. https://doi.org/10.1371/journal.pone.0075762. (30) Kurkinen, S. T.; Postila, P. A. Improving Docking Performance Using Negative ImageBased

Rescoring.

Front.

Pharmacol.

2018.

https://doi.org/https://doi.org/10.3389/fphar.2018.00260. (31) Axenopoulos, A.; Daras, P.; Papadopoulos, G. E.; Houstis, E. N. SP-Dock: Protein-Protein Docking Using Shape and Physicochemical Complementarity. IEEE/ACM Trans. Comput. Biol. Bioinform. 2013, 10 (1), 135–150. https://doi.org/10.1109/TCBB.2012.149. (32) Gowthaman, R.; Lyskov, S.; Karanicolas, J. DARC 2.0: Improved Docking and Virtual Screening

at

Protein

Interaction

Sites.

PLOS

ONE

2015,

10

(7),

e0131612.

https://doi.org/10.1371/journal.pone.0131612. (33) Rocchia, W.; Sridharan, S.; Nicholls, A.; Alexov, E.; Chiabrera, A.; Honig, B. Rapid GridBased Construction of the Molecular Surface and the Use of Induced Surface Charge to Calculate


46

Page 47 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Reaction Field Energies: Applications to the Molecular Systems and Geometric Objects. J Comput Chem 2002, 23 (1), 128–137. https://doi.org/10.1002/jcc.1161. (34) Rocchia, W.; Alexov, E.; Honig, B. Extending the Applicability of the Nonlinear Poisson−Boltzmann Equation: Multiple Dielectric Constants and Multivalent Ions. J. Phys. Chem. B 2001, 105 (28), 6507–6514. https://doi.org/10.1021/jp010454y. (35) Holst, M.; Baker, N.; Wang, F. Adaptive Multilevel Finite Element Solution of the Poisson–Boltzmann Equation I. Algorithms and Examples. J. Comput. Chem. 2000, 21 (15), 1319– 1342. https://doi.org/10.1002/1096-987X(20001130)21:153.0.CO;2-8. (36) Baker, N.; Holst, M.; Wang, F. Adaptive Multilevel Finite Element Solution of the Poisson–Boltzmann Equation II. Refinement at Solvent-Accessible Surfaces in Biomolecular Systems. J. Comput. Chem. 2000, 21 (15), 1343–1352. https://doi.org/10.1002/1096987X(20001130)21:153.0.CO;2-K. (37) Nicholls, A.; Bharadwaj, R.; Honig, B. GRASP - Graphical Representation and Analysis of Surface Properties. Biophys J. 1993, 64 (part 2), A166. (38) Klamt, A.; Thormann, M.; Wichmann, K.; Tosco, P. COSMOsar3D: Molecular Field Analysis Based on Local COSMO σ-Profiles. J. Chem. Inf. Model. 2012, 52 (8), 2157–2164. https://doi.org/10.1021/ci300231t. (39) Thormann, M.; Klamt, A.; Wichmann, K. COSMOsim3D: 3D-Similarity and Alignment Based on COSMO Polarization Charge Densities. J. Chem. Inf. Model. 2012, 52 (8), 2149–2156. https://doi.org/10.1021/ci300205p.


47


Page 48 of 66

(40) Thormann, M.; Klamt, A.; Hornig, M.; Almstetter, M. COSMOsim: Bioisosteric Similarity Based on COSMO-RS σ Profiles. J. Chem. Inf. Model. 2006, 46 (3), 1040–1053. https://doi.org/10.1021/ci050464m. (41) Chessari, G.; Buck, I. M.; Day, J. E. H.; Day, P. J.; Iqbal, A.; Johnson, C. N.; Lewis, E. J.; Martins, V.; Miller, D.; Reader, M.; Rees, D. C.; Rich, S. J.; Tamanini, E.; Vitorino, M.; Ward, G. A.; Williams, P. A.; Williams, G.; Wilsher, N. E.; Woolford, A. J.-A. Fragment-Based Drug Discovery Targeting Inhibitor of Apoptosis Proteins: Discovery of a Non-Alanine Lead Series with Dual Activity Against CIAP1 and XIAP. J. Med. Chem. 2015, 58 (16), 6574–6588. https://doi.org/10.1021/acs.jmedchem.5b00706. (42) Vinter, J. G. Extended Electron Distributions Applied to the Molecular Mechanics of Some Intermolecular

Interactions.

J

Comput

Aided

Mol

Des

1994,

8

(6),

653–668.

https://doi.org/https://doi.org/10.1007/BF00124013. (43) Vinter, J. G. Extended Electron Distributions Applied to the Molecular Mechanics of Some Intermolecular Interactions. II. Organic Complexes. J Comput Aided Mol Des 1996, 10 (5), 417– 426. https://doi.org/https://doi.org/10.1007/BF00124473. (44) Cheeseright, T.; Mackey, M.; Rose, S.; Vinter, A. Molecular Field Extrema as Descriptors of Biological Activity: Definition and Validation. J. Chem. Inf. Model. 2006, 46 (2), 665–676. https://doi.org/10.1021/ci050357s. (45) Johnson, C. N.; Erlanson, D. A.; Murray, C. W.; Rees, D. C. Fragment-to-Lead Medicinal Chemistry

Publications

in

2015.

J.

Med.

Chem.

pp

89–99,

60

(1).

https://doi.org/10.1021/acs.jmedchem.6b01123.


48

Page 49 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(46) Johnson, C. N.; Erlanson, D. A.; Jahnke, W.; Mortenson, P. N.; Rees, D. C. Fragment-toLead Medicinal Chemistry Publications in 2016. J. Med. Chem. 2018, 61 (5), 1774–1784. https://doi.org/10.1021/acs.jmedchem.7b01298. (47) Chessari, G.; Hunter, C. A.; Low, C. M.; Packer, M. J.; Vinter, J. G.; Zonta, C. An Evaluation of Force-Field Treatments of Aromatic Interactions. Chem. Eur. J. 2002, 8 (13), 2860– 2867. (48) Mehler, E. L. The Lorentz-Debye-Sack Theory and Dielectric Screening of Electrostatic Effects in Proteins and Nucleic Acids. In Molecular Electrostatic Potentials: Concepts and Applications, Theoretical and Computational Chemistry; Elsevier Science B.V.: Amsterdam, Netherlands, 1996; Vol. 3. (49) Xu, D.; Li, H.; Zhang, Y. Protein Depth Calculation and the Use for Improving Accuracy of

Protein

Fold

Recognition.

J

Comput

Biol

2013,

20

(10),

805–816.

https://doi.org/10.1089/cmb.2013.0071. (50) Xu, D.; Zhang, Y. Generating Triangulated Macromolecular Surfaces by Euclidean Distance

Transform.

PLoS

ONE

2009,

4

(12),

e8140.

https://doi.org/10.1371/journal.pone.0008140. (51) BioMolTech. www.biomoltech.com. (52) Stroganov, O. V.; Novikov, F. N.; Zeifman, A. A.; Stroylov, V. S.; Chilov, G. G. TSAR a New Graph Theoretical Approach to Computational Modeling of Protein Side Chain Flexibility Modeling of Ionization Properties of Proteins. Proteins 2011, 79 (9), 2693–2710. https://doi.org/10.1002/prot.23099.


49


Page 50 of 66

(53) Flare Version 2.0. http://www.cresset-group.com/flare/. (54) Forge Version 10.5. http://www.cresset-group.com/forge/. (55) Christopher, J. A.; Aves, S. J.; Bennett, K. A.; Doré, A. S.; Errey, J. C.; Jazayeri, A.; Marshall, F. H.; Okrasa, K.; Serrano-Vega, M. J.; Tehan, B. G.; Wiggin, G. R.; Congreve, M. Fragment and Structure-Based Drug Discovery for a Class C GPCR: Discovery of the MGlu5 Negative Allosteric Modulator HTL14242 (3-Chloro-5-[6-(5-Fluoropyridin-2-Yl)Pyrimidin-4Yl]Benzonitrile).

J.

Med.

Chem.

2015,

58

(16),

6653–6664.

https://doi.org/10.1021/acs.jmedchem.5b00892. (56) Watanabe, C.; Watanabe, H.; Fukuzawa, K.; Parker, L. J.; Okiyama, Y.; Yuki, H.; Yokoyama, S.; Nakano, H.; Tanaka, S.; Honma, T. Theoretical Analysis of Activity Cliffs among Benzofuranone-Class Pim1 Inhibitors Using the Fragment Molecular Orbital Method with Molecular Mechanics Poisson–Boltzmann Surface Area (FMO+MM-PBSA) Approach. J. Chem. Inf. Model. 2017, 57 (12), 2996–3010. https://doi.org/10.1021/acs.jcim.7b00110. (57) Friberg, A.; Vigil, D.; Zhao, B.; Daniels, R. N.; Burke, J. P.; Garcia-Barrantes, P. M.; Camper, D.; Chauder, B. A.; Lee, T.; Olejniczak, E. T.; Fesik, S. W. Discovery of Potent Myeloid Cell Leukemia 1 (Mcl-1) Inhibitors Using Fragment-Based Methods and Structure-Based Design. J. Med. Chem. 2013, 56 (1), 15–30. https://doi.org/10.1021/jm301448p. (58) Wang, L.; Wu, Y.; Deng, Y.; Kim, B.; Pierce, L.; Krilov, G.; Lupyan, D.; Robinson, S.; Dahlgren, M. K.; Greenwood, J.; Romero, D. L.; Masse, C.; Knight, J. L.; Steinbrecher, T.; Beuming, T.; Damm, W.; Harder, E.; Sherman, W.; Brewer, M.; Wester, R.; Murcko, M.; Frye, L.; Farid, R.; Lin, T.; Mobley, D. L.; Jorgensen, W. L.; Berne, B. J.; Friesner, R. A.; Abel, R. Accurate and Reliable Prediction of Relative Ligand Binding Potency in Prospective Drug


50

Page 51 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Discovery by Way of a Modern Free-Energy Calculation Protocol and Force Field. J. Am. Chem. Soc. 2015, 137 (7), 2695–2703. https://doi.org/10.1021/ja512751q. (59) Hardegger, L. A.; Kuhn, B.; Spinnler, B.; Anselm, L.; Ecabert, R.; Stihle, M.; Gsell, B.; Thoma, R.; Diez, J.; Benz, J.; Plancher, J. M.; Hartmann, G.; Isshiki, Y.; Morikami, K.; Shimma, N.; Haap, W.; Banner, D. W.; Diederich, F. Halogen Bonding at the Active Sites of Human Cathepsin L and MEK1 Kinase: Efficient Interactions in Different Environments. ChemMedChem 2011, 6 (11), 2048–2054. https://doi.org/10.1002/cmdc.201100353. (60) Wodicka, L. M.; Ciceri, P.; Davis, M. I.; Hunt, J. P.; Floyd, M.; Salerno, S.; Hua, X. H.; Ford, J. M.; Armstrong, R. C.; Zarrinkar, P. P.; Treiber, D. K. Activation State-Dependent Binding of Small Molecule Kinase Inhibitors: Structural Insights from Biochemistry. Chem. Biol. 2010, 17 (11), 1241–1249. https://doi.org/10.1016/j.chembiol.2010.09.010. (61) Davis, M. I.; Hunt, J. P.; Herrgard, S.; Ciceri, P.; Wodicka, L. M.; Pallares, G.; Hocker, M.; Treiber, D. K.; Zarrinkar, P. P. Comprehensive Analysis of Kinase Inhibitor Selectivity. Nat. Biotechnol. 2011, 29 (11), 1046–1051. https://doi.org/10.1038/nbt.1990. (62) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. The Protein Data Bank. Nucleic Acids Res 2000, 28 (1), 235–242. (63) Seeliger, M. A.; Nagar, B.; Frank, F.; Cao, X.; Henderson, M. N.; Kuriyan, J. C-Src Binds to the Cancer Drug Imatinib with an Inactive Abl/c-Kit Conformation and a Distributed Thermodynamic

Penalty.

Structure

2007,

15

(3),

299–311.

https://doi.org/10.1016/j.str.2007.01.015.


51


Page 52 of 66

(64) Atwell, S.; Adams, J. M.; Badger, J.; Buchanan, M. D.; Feil, I. K.; Froning, K. J.; Gao, X.; Hendle, J.; Keegan, K.; Leon, B. C.; Müller-Dieckmann, H. J.; Nienaber, V. L.; Noland, B. W.; Post, K.; Rajashankar, K. R.; Ramos, A.; Russell, M.; Burley, S. K.; Buchanan, S. G. A Novel Mode of Gleevec Binding Is Revealed by the Structure of Spleen Tyrosine Kinase. J. Biol. Chem. 2004, 279 (53), 55827–55832. https://doi.org/10.1074/jbc.M409792200. (65) Namboodiri, H. V.; Bukhtiyarova, M.; Ramcharan, J.; Karpusas, M.; Lee, Y.; Springman, E. B. Analysis of Imatinib and Sorafenib Binding to P38α Compared with C-Abl and b-Raf Provides Structural Insights for Understanding the Selectivity of Inhibitors Targeting the DFGOut

Form

of

Protein

Kinases.

Biochemistry

2010,

49

(17),

3611–3618.

https://doi.org/10.1021/bi100070r. (66) Hyre, D. E.; Le Trong, I.; Merritt, E. A.; Eccleston, J. F.; Green, N. M.; Stenkamp, R. E.; Stayton, P. S. Cooperative Hydrogen Bond Interactions in the Streptavidin–Biotin System. Protein Sci 2006, 15 (3), 459–467. https://doi.org/10.1110/ps.051970306. (67) Gohlke, H.; Klebe, G. Approaches to the Description and Prediction of the Binding Affinity of Small-Molecule Ligands to Macromolecular Receptors. Angewandte Chemie International Edition

2002,

41

(15),

2644–2676.

https://doi.org/10.1002/1521-

3773(20020802)41:153.0.CO;2-O. (68) Hu, L.; Benson, M. L.; Smith, R. D.; Lerner, M. G.; Carlson, H. A. Binding MOAD (Mother Of All Databases). Proteins: Structure, Function, and Bioinformatics 2005, 60 (3), 333– 340. https://doi.org/10.1002/prot.20512.


52

Page 53 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(69) Ahmed, A.; Smith, R. D.; Clark, J. J.; Dunbar, J. B.; Carlson, H. A. Recent Improvements to Binding MOAD: A Resource for Protein–Ligand Binding Affinities and Structures. Nucleic Acids Res 2015, 43 (D1), D465–D469. https://doi.org/10.1093/nar/gku1088. (70) Smith, R. D.; Engdahl, A. L.; Dunbar, J. B.; Carlson, H. A. Biophysical Limits of ProteinLigand Binding. J Chem Inf Model 2012, 52 (8), 2098–2106. https://doi.org/10.1021/ci200612f. (71) McCoy, A. J.; Chandana Epa, V.; Colman, P. M. Electrostatic Complementarity at Protein/Protein Interfaces. Journal of Molecular Biology 1997, 268 (2), 570–584. https://doi.org/10.1006/jmbi.1997.0987. (72) Ehrlich, S.; Göller, A. H.; Grimme, S. Towards Full Quantum-Mechanics-Based Protein– Ligand

Binding

Affinities.

ChemPhysChem

2017,

18

(8),

898–905.

https://doi.org/10.1002/cphc.201700082. (73) Mobley, D. L.; Dill, K. A. Binding of Small-Molecule Ligands to Proteins: “What You See”

Is

Not

Always

“What

You

Get.”

Structure

2009,

17

(4),

489–498.

https://doi.org/10.1016/j.str.2009.02.010. (74) Klamt, A.; Jonas, V.; Bürger, T.; Lohrenz, J. C. W. Refinement and Parametrization of COSMO-RS. J. Phys. Chem. A 1998, 102 (26), 5074–5085. https://doi.org/10.1021/jp980017s.


53


Page 54 of 66

Table of Contents graphic


54

Page 55 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Electrostatic complementarity (EC) of biotin-streptavidin complex (PDB 3RY2). The electrostatic potentials (ESP) of protein and ligand (blue = negative and red = positive ESP) and the protein-ligand EC (green = complementary, red = electrostatic clash) are projected on the protein (left panel) or ligand (right panel) solvent-accessible surface (SAS), respectively.



Figure 2. Structures of mGLU5 inhibitors and X-ray binding modes of compounds 1 and 2. Both 1 (left panel; PDB 5CGC) and 2 (right panel; PDB 5CGD) are negative allosteric modulators of mGLU5 and bind at the mavoglurant binding site, which is formed by residues of TM2, TM3, TM5, TM6, and TM7.


Page 56 of 66

Page 57 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Comparison of ESP and EC surfaces of 5 representative mGLU5 ligands with increasing pIC50 values and correlation of mGLU5 bioactivities with EC scores. Positive and negative ESP values are colored in red and blue and electrostatically complementary and clashing regions are colored in green and red, respectively. The ESP ligand surfaces in the top panel show that the negative potential (blue color) on the face of the phenyl ring gradually decreases with introduction of electron withdrawing substituents. The EC surface in the middle panel shows that the electrostatic clash (red color) between the phenyl ring and the Trp785 indole (about 4Å distance to ligand phenyl ring) is minimized with decreasing negative electrostatic potential of the phenyl ring, which is in line with the gradually improving pIC50 values. The visually increasing electrostatic complementarity is also reflected in the improving EC scores. The correlation plots show that both EC and ECR scores correlate well with pIC50 values of the whole data set. Outliers 1a and 10 are highlighted and discussed in the text.



Figure 4. Comparison of ESP and EC surfaces of compound 6 (5CGC model) and 8 (5CGD model) with their fluorinated analogs. Positive and negative ESP values are colored in red and blue and electrostatically complementary and clashing regions are colored in green and red, respectively. (A) Upon 4-fluorination of the pyrazole ring, the negative electrostatic potential of the fluorine substituent electrostatically clashes with oxygen and nitrogen atoms of the surrounding residues. (B) Fluorination of the 5 position of the pyridine ring of compound 8 decreases the electrostatic clash with Ser654 and Tyr659 and improves the complementarity with Ala810 and Ala813.


Page 58 of 66

Page 59 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. Structures of XIAP inhibitors, binding mode of inhibitors 12 and 13 bound to XIAP (PDB 5C7A and 5C7D) and their EC projected on the protein SAS. Lead compound 13, which is about 65-fold more potent than 12, shows an improved EC, especially in proximity to the Lys297 side chain and the Gly306 backbone carbonyl group.



Figure 6. Electrostatic complementarity of XIAP inhibitors. With increasing bioactivity, the degree of electrostatic clash (red areas) with Gly306 carbonyl and Lys297 side chain, as exemplified by compound 12 and 15, is lessened or even turned into electrostatic complementarity (green areas). Electron-withdrawing 6-substituents decrease the negative potential of the aromatic π-plane and change the electrostatic potential at this position from positive (15, amine group) to a negative potential (17, chloro substituent).


Page 60 of 66

Page 61 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7. Binding mode of benzofuranone inhibitor 18 of PIM1 and bioactivities of derived aza analogs. (A) X-ray structure of 18 in complex with PIM1 (PDB 5VUC) and its electrostatic complementarity to the binding pocket. (B) Chemical structures and respective bioactivities of aza analogs of 18. (C) Correlation of pIC50 values with EC scores.



Figure 8. Upper and left panel: Mcl-1 inhibitor structures and binding mode of Mcl-1 inhibitor 24 (PDB 6B4L). Right panel: a subset of the compound data shows systematic variation of chlorine substitution on the core scaffold.


Page 62 of 66

Page 63 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 9. Correlation of EC scores with bioactivities for the Mcl-1 chlorine scan subset. EC maps show that 6-chloro substitution (29) shows the most favorable electrostatics due to formation of a halogen bond with the Ala227 backbone carbonyl group. In both EC and ECR plots, compounds with 6-Cl substitution (data points shown as black ■) are enriched at high EC scores.



Figure 10. EC analysis of MEK1 inhibitors that form halogen bond to the backbone (BB) carbonyl group of Val127. The depicted structures show the complex structure of the iodine compound or modelled analogs with MEK1 (PDB 3DY7). The electrostatic complementarity decreases (both EC maps and scores) from iodine to fluorine (no halogen bonding interaction), as expected from theoretical halogen bonding strengths.1


Page 64 of 66

Page 65 of 66 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 11. Electrostatic complementarity analysis of imatinib kinase selectivity. (A) Imatinib structure and its respective pKD values for imatinib-kinase complexes available in the PDB. (B) Correlation of EC with

imatinib pKD values for kinase subset. The c-SRC data point was not included in the Pearson’s R2 correlation scores as weak imatinib binding to c-SRC is mainly due to an energetic penalty of the imatinib bioactive conformation.63 (C) EC maps highlight differences between imatinib binding to ABL and the ancestral ANCAS kinases.



Figure 12. EC analysis of wild-type and mutant biotin-streptavidin complexes. EC scores can rank the experimentally determined binding affinities. The ligand EC maps highlight the decreased electrostatic complementarity in proximity to S45A and D128A mutations due to loss of hydrogen bonding acceptor motifs.


Page 66 of 66

Electrostatic complementarity as a fast and ... - ACS Publications

Recommend Documents