0JCTG
[ournal of Chemical Theoiy and Computation
THC ll«l*tU4TV Of
NEWCASTLE
AU3TRMJA
Subscriber access provided by University of Newcastle, Australia
Article
Rapid and Reliable Binding Affinity Prediction of Bromodomain Inhibitors: a Computational Study
Shunzhou Wan, Agastya P. Bhati, Stefan J. Zasada, Ian Wall, Darren Green, Paul Bamborough, and Peter Vivian Coveney
J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00794 • Publication Date (Web): 22 Dec 2016
Downloaded from http://pubs.acs.org on December 26, 2016
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 free 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 accessible to all readers and 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.
Journal of Chemical Theory and Computation 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.
10 11 12
Rapid and Reliable Binding Affinity Prediction of Bromodomain Inhibitors: a Computational Study
Centre for Computational Science, Department of Chemistry, University College London,
GlaxoSmithKline, Gunnels Wood Road, Stevenage, Hertfordshire, SG1 2NY, United Kingdom
17 Shunzhou Wan,a Agastya P. Bhati,a Stefan J. Zasada, a Ian Wall,b Darren Green,b Paul
19 Bamborough,b and Peter V. Coveneya*
20 21 22
25 London WC1H 0AJ, United Kingdom.
32 Corresponding to: p.v.coveney@ucl.ac.uk
36 ABSTRACT. Binding free energies of bromodomain inhibitors are calculated with recently
38 formulated approaches, namely ESMACS (enhanced sampling of molecular dynamics with
40 approximation of continuum solvent) and TIES (thermodynamic integration with enhanced
43 sampling). A set of compounds is provided by GlaxoSmithKline, which represents a range of
45 chemical functionality and binding affinities. The predicted binding free energies exhibit a good
50 excellent correlation of 0.92 from the TIES approach where applicable. Given access to suitable
52 high end computing resources and a high degree of automation, we can compute individual
57 larger than 0.34 kcal/mol and 1.71 kcal/mol for the 1- and 3-trajectory ESMACS approaches.
Spearman correlation of 0.78 with the experimental data from the 3-trajectory ESMACS, and an
binding affinities in a few hours with precisions no greater than 0.2 kcal/mol for TIES, and no
1. INTRODUCTION
discovery programmes. In particular, the FEP+ implementation of Free Energy Perturbation
6 Computational chemistry has been an established tool in drug discovery for a number of years.
8 The number of crystal structures available in the public domain and within pharmaceutical
11 companies on which to base computational studies continues to rise rapidly. Despite the increase
13 in resources applied and in experimental data on which to base the studies, industrial structure-
15 based design approaches have evolved very little in recent times.1 In particular, the approaches
18 are largely qualitative and largely dependent on the experience and knowledge of the
20 practitioner. - Attempts to quantify protein-ligand binding affinities are rare. Expert practitioners
23 have little confidence in existing tools to make robust predictions and certainly not to do so on a
25 timescale that can substantially impact medicinal chemistry programmes.
27 Recently, there has been a renewed interest in the use of free energy calculations in drug
32 (FEP) has shown potential to improve the ability to predict protein-ligand binding affinities on an
34 industrially relevant timescale.4 Research is ongoing to understand how broadly applicable the
37 method is, and how accurate its predictions are when applied to active drug discovery
39 programmes.
41 Although FEP+ applies replica exchange solute tempering (REST) where exchange moves are
44 made between different X windows, its predictions, like those from many other approaches, are
46 generated from a single output for each pair of mutations. Advances in high-end computing
49 capabilities offer the opportunity to run vast numbers of calculations in parallel. The application
51 of these computational capabilities to free energy calculations allows results to be returned
53 quickly and multiple replicas of simulations5 to be run, leading to tighter control of standard
56 errors. If such an approach could be validated and implemented in an industrial setting it would
represent a major step forward in structure-based design capabilities. The first step in this
antagonise the binding of acetylated histone tails to these modules have been shown to exert
6 process is to validate the performance on an industrially relevant data set.
8 Depending on the reliability, rapidity and throughput of these calculations, they might find
11 application at various stages of the drug discovery and development process across the wider
13 pharmaceutical industry. As these methods require significant compute resources beyond
15 existing in-house industrial capacity, assessment (and any subsequent adoption) of the
18 methodology requires access to high performance computing resources.
20 Research into epigenetic proteins is currently a major and rapidly evolving focus for the
22 6 8
23 pharmaceutical industry. " Bromodomain-containing proteins, and in particular the four
25 members of the BET (bromodomain and extra terminal domain) family, each containing two
27 bromodomains, have been widely studied. Small molecule inhibitors able to competitively
32 profound effects on gene expression and have shown promising preclinical efficacy in
34 pathologies ranging from cancer to inflammation. Indeed, several compounds are progressing
37 through early-stage clinical trials and are showing exciting early results.9 Most inhibitors
39 reported show similar binding potencies to all BET family bromodomains. A representative
41 inhibitor-protein structure is given in Figure 1, showing the key elements for the inhibitor
44 binding. This study will concentrate on the first bromodomain of bromodomain-containing
46 protein 4 (BRD4-BD1) for which extensive crystallographic and ligand binding data are
48 available.10-12
51 The purpose of the present study is to assess the potential for rapid, accurate, precise and
53 reproducible binding affinity calculations based on the use of a Binding Affinity Calculator
56 (BAC) software tool and associated services. The approach is based on the use of high
10 11 12
20 21 22
performance computing in an automated workflow which builds models, runs large numbers of replica calculations, and analyses the output data in order to place reliable standard error bounds on predicted binding affinities.
Figure 1. The Bromodomain inhibitor I-BET726 and its binding mode in BRD4-BD1. Two views are displayed for the binding mode (PDB id: 4BJX14), in which I-BET72615 is represented as stick in cyan/blue/red/green, the protein by cartoon in silver, the crystallographic water molecules in red ball, and clipped protein surfaces are shown in orange.
2. COMPUTATIONAL SECTION
Models. In this study, chemical structures of 16 BRD4-BD1 ligands based on a single tetrahydroquinoline (THQ) template15 were provided by GSK (Table 1). The compound set was designed to represent a range of chemical functionality and binding affinities, but also to contain
4 sets of closely related compounds with key SAR trends. Specifically, there are two growth
6 vectors which cause a drop off in potency, a growth vector where substantial structural variation
8 is tolerated and an enantiomeric pair where one isomer is significantly more potent than the
11 other. These will be described below. The calculations were then performed by the UCL group,
13 initially blind, to investigate the ability of BAC to reproduce the SAR trends. The experimental
15 binding affinities were not released to the UCL group until all of the computations had been
18 completed and the estimation of binding free energies reported (see Figure 3 below). These
20 predictions were subsequently compared with the independently measured experimental free
23 energies of binding (Table 1) as part of the assessment of the reliability of the method.
25 The X-ray crystal structure used to represent BRD4-BD1 (published coordinates PDB id:
27 4BJX14) is a complex with a THQ (I-BET72615) chosen in order to reproduce the likely
32 binding mode of I-BET726 is shown in Figure 1. Both of the significant protein-ligand
34 hydrogen-bonding interactions take place through the THQ N1-acetamide carbonyl group, which
37 interacts directly with the Asn140 sidechain of the acetyl-lysine binding site and also through the
39 W1 water to the sidechain of Tyr97. The W1 water forms part of a chain of water molecules
41 buried within the site (W2-W4, Figure 1). Substituents larger than acetamide at the N1-position
44 extend into the water-filled part of the site. Our ligand data-set included a small series of
46 increasing size at N1 (1, 8 and 9, see Table 1) to probe whether the computer methodology can
conditions where free-energy calculations might be used in a real drug discovery program. The
accurately predict the outcome of growing into this region.
10 11 12
20 21 22
Table 1. Compounds used in this study, ordered according to their R4 group, and their experimental IC50S with standard deviations converted into binding free energies AG (kcal/mol).
Compounda R1 R2 R3 R4 pIC50b PIC50 n AGb
2 H Me Me H 5.6 ± 0.06 6 -7.7 ± 0.08
3 H Me Me 6.8 ± 0.04 2 -9.3 ± 0.05
4 H Me Me 6.8 ± 0.09 3 -9.3 ± 0.12
5 H Me Me 7.9 ± 0.03 2 -10.8 ± 0.04
6 H Me Me 5.6 ± 0.01 2 -7.7 ± 0.01
7 H Me Me h=:rN 5.8 ± 0.07 4 -8.0 ± 0.10
1 H Me Me 7.0 ± 0.05 3 -9.6 ± 0.07
8 H Me Et 6.5 ± 0.01 2 -8.9 ± 0.01
9 H Me i-Pr < 4.3c 4 > -5.9
10 H Me Me 0 7.6 ± 0.1 16 -10.4 ± 0.14
16d H Me Me -CK 0 5.4 ± 0.28 5 -7.4 ± 0.38
11 H Et Me -CK 0 6.8 ± 0.3 4 -9.3 ± 0.41
12 H Pr Me 0 5.5 ± 0.01 4 -7.5 ± 0.01
13 H Pr Me ^CM, 5.4 ± 0.01 2 -7.4 ± 0.01
14 H Et Me 6.7 ± 0.2 4 -9.2 ± 0.27
15 Cl Me Me 7.8 ± 0.1 44 -10.7 ± 0.14
a Compounds 1-9 are electrostatically neutral, 10-12 and 16 are positively charged, and 13-15 are negatively charged.
b Statistical errors were calculated from repeated IC50 measurements. c There was no activity at the highest concentration (50 (iM) tested.
d All compounds are the 2-(S) 4-(R) isomers (Figure 1) expect compound 16 which is 2-(R) 4-(S).
The THQ 4-anilino substituent projects onto the "WPF-shelf" subsite outside the acetyl-lysine
4 The THQ R2-position (S)-methyl group occupies a small lipophilic site between the sidechains
6 of three residues of the BRD4-BD1 ZA-loop (Val87, Leu92 and Tyr97). Carbon-carbon contacts
8 for all three lie within 4.25A, apparently offering little room for extension of this group without
11 some structural reorganization. Our ligand data-set includes a small number of analogues
13 exploring larger substituents at the R2-position (10-14) to see whether the simulation can
15 accurately predict the consequences of pushing on these residues.
18 The THQ 6-phenyl substituent fills the narrow "ZA channel" between Trp81 and Leu92. This
20 ring can be substituted at the R3- and R4- positions, but R2-substitution is detrimental because
22 this causes an increase in the inter-ring torsion. The ZA channel does not seem to accommodate
25 the resulting wider 6-substituents. The remainder of our data-set includes a variety of
27 substituents probing electrostatic and steric changes at this position.
32 pocket, close to Trp81. SAR around this ring has been published previously , showing that
34 substitution had small effects on potency, probably because one edge of the aniline ring is
37 solvent-exposed. Therefore we did not attempt to vary substituents at this position in our data-
39 set, preferring to keep it as constant as possible while exploring changes elsewhere.
41 In this series the 2-(S) 4-(R) isomers (see Figure 1) are significantly more potent than their
44 alternative trans enantiomers. We included an enantiomeric pair of 2-(S) 4-(R) (compound 10)
46 and 2-(R) 4-(S) (compound 16) isomers to explore whether the simulations were capable of
51 The ligand-protein complexes were constructed by replacing the ligand in the PDB file with
53 the ligands of interest (see structures in Supporting Information). For the congeneric compounds
56 studied here, it is plausible to assume that they bind in the same mode in which all key
distinguishing between them.
4 compound-protein interactions are preserved (Figure 1). For compounds with significant
6 structural differences, computational docking may be required to generate reasonable complex
8 structures16. Preparation and setup of the simulations were implemented using BAC13, including
11 parameterization of the compounds, solvation of the complexes, electrostatic neutralization of
13 the systems by adding counterions and generation of configurations files for the simulations. The
15 AMBER ff99SB-ILDN17 force field was used for the protein, and TIP3P for water molecules.
18 Standard protonation states were assigned to all titratable residues at pH 7, with histidines
20 protonated on the s position (HIE). Compound parameters were produced using the general
22 AMBER force field (GAFF)18 with Gaussian 0319 to optimise compound geometries and to
25 determine electrostatic potentials at the Hartree-Fock level with 6-31G** basis functions. The
27 restrained electrostatic potential (RESP) module in the AMBER package20 was used to calculate
30 the partial atomic charges for the compounds. All systems were solvated in orthorhombic water
32 boxes with a minimum extension from the protein of 14 A.
34 Some of the ligands could adopt several rotamers for the R4 group (Table 1) in solution, which
37 due to the constrained environment of the active site are unlikely to be sampled in a single bound
39 simulation run with the protocol used (see Simulations sub-section). To decide which rotamer(s)
42 to start with, metadynamics simulations21 were employed which used a history-dependent
44 biasing potential to explore the conformational space of the chosen degrees of freedom, here the
46 rotatable bond(s) involving groups R2 and R4 (see details in the Supporting Information). Five
49 replicas were used for each metadynamics study to have a reasonable convergence of the
51 potential profile of the chosen dihedral angle(s) of a ligand in complexed form. The calculated
53 potential of mean force (PMF) was used to determine the most favourable rotamer(s) from which
56 ESMACS simulations initiate. When unambiguous result ensued, the energetically most
11 and 12 which are positively charged; and set 3 including ligands 13, 14 and 15 which are
4 favourable rotamer was chosen. For some ligands, more than one rotamer showed similar free
6 energies. In these cases, multiple initial structures were prepared using the rotamers suggested by
8 the metadynamics study. The ligands with multiple rotamers were: 4, 7, 10, 11, 12 and 16. There
11 are two rotamers for each of these, except for 7 for which there are three generated by the flip
13 and twist between the two ring planes. This resulted in a total number of 23 complexes being
15 simulated in the ESMACS study.
18 For each ligand with multiple rotamers, the energetic properties were analysed and the most
20 favourable rotamer was chosen as the final result of the ESMACS study (see Results section),
23 and was used as the initial structure in the TIES study. For the latter, three subsets of the ligands
25 were selected, within which relative binding free energies were calculated with the TIES
27 approach. The subsets of the ligands were: set 1 including ligands 1-9; set 2 including ligands 10,
32 negatively charged. The division of the full set of the cognate ligands into subsets is necessary
34 because TIES, just as any other TI and FEP based method, encounters specific difficulties owing
37 to major adjustments in long range electrostatic interactions for the "congeners" where the net
39 charges change.22 For each subset, ligands were paired based on their similarities, and TIES
41 calculations were performed to alchemically mutate one ligand to another (see Theoretical Basis
44 below).
46 Theoretical Basis. The UCL group have recently introduced new protocols for binding free
51 continuum solvent" (ESMACS) " and "thermodynamic integration with enhanced sampling"
53 (TIES).23,25 ESMACS is based on the molecular mechanics Poisson-Boltzmann surface area
56 method (MMPBSA)26 which makes a continuum approximation for the aqueous solvent, while
energy calculations, termed "enhanced sampling of molecular dynamics with approximation of
4 TIES centres, as the abbreviation indicates, on thermodynamic integration (TI). Although the
6 approaches are built around the standard MMPBSA and TI methodologies, our abbreviations are
8 used to emphasise the central importance of the ensemble based nature of the protocols
11 employed as well as, in the case of ESMACS the wider generality and flexibility of the
13 protocol.5'23-24 The size of the statistical mechanical ensembles is determined systematically so
j c 23 24 27
15 that predictions are accurate, precise and reliable. - ' Moreover, the term "MMPBSA" is used
18 to mean very different things in numerous j ournal articles and textbooks, including calculations
20 based on single docked structures4 or on simulation trajectories, calculations with or without the
23 inclusion of conformational entropies, and almost wholly using the so-called 1-trajectory
25 approach. In the ESMACS protocol, we always mean a fully converged, ensemble-based
27 determination of the free energy of binding from either a one-, two- or three-trajectory method:
29 27 29
30 this includes both the configurational entropy and the association free energy, and (where
31 24 30
32 appropriate) the relative or absolute adaptation energy. Statistical analyses were performed
34 throughout the study for all of the quantities obtained. A Binding Affinity Calculator (BAC)13
37 was used to perform ESMACS and TIES studies. BAC constitutes a computational pipeline built
39 from selected software tools and services (see Supporting Information), and relies on access to a
41 range of computing resources. It automates much of the complexity of building, running and
44 marshalling the molecular dynamics simulations, as well as collecting and analysing data.
46 Integration and automation are central to the reliability of the method, ensuring that the results
51 In ESMACS, the free energy is evaluated approximately based on the extended MMPBSA
24 27 29
53 method, , including configurational entropy, and the free energy of association. Free energy
are reproducible and can be delivered rapidly.
changes (AGbinding) are determined for the molecules in their solvated states. The binding free
6 energy change is then calculated as:
9 binding = Gcomplex — Greceptor — Gligand. (-0
10 where Gcompiex, Greceptor and Giigand are free energies of the complex, the receptor and the ligands, 12
13 respectively. The Possion-Boltzmann calculation is performed here using the Amber built-in
15 PBSA solver31 with radii taken from the parameter/topology files, and the configurational
18 entropy is calculated by a normal mode (NMODE) method.
20 The three terms on the right hand side of Equation 1 can be generated from single simulations
23 of the complexes or from separated simulations of complexes, receptor(s) and ligand(s). The
25 former is the so-called 1-trajectory (1-traj) method of which the trajectories of the receptor(s) and
27 the ligand(s) are extracted from those of the complexes. The latter is the so-called 3-trajectory (3-
32 components. In the drug development field, binding is usually investigated for a set of ligands
34 bound to the same protein target. The free energy of the receptor Greceptor will then be the same
37 for all ligands in the 3-traj method. In the current study, we employ three approaches to calculate
39 the binding free energies AGbinding: (A) using single simulations of the complexes (1-traj); (B)
44 (denoted as 2-traj); and (C) the same as (B) but invoking separated ligand simulations for the
46 derivation of Gligand (denoted as 3-traj).
traj) method of which the trajectories are generated from separate simulations of the three
same as (A) but using the average of the receptor free energies (Greceptor) from the 1-traj method
In TIES, an alchemical transformation for the mutated entity, either the ligand or the protein, is
51 used in both aqueous solution and within the ligand-protein complex. The free energy changes of
53 the alchemical mutation processes, AGalfh and AGaOmplex, are calculated by:
8 and the binding free energy difference is calculated from
1 0 AAGMnding = AGalcn - LG^lex (3)
13 Here X (0 < X < 1) is a coupling parameter such that X=0 and X=1 correspond to the initial and
15 final thermodynamic states, and V(X) is the potential energy of an intermediate state X.
18 denotes an ensemble average over configurations representative of the state X. To avoid the well-
20 known "end-point catastrophe"32, a soft-core potential is used, along with an approach to
23 decouple at different rates the electrostatic and van der Waals interactions of the perturbed atoms
25 with their environment (for more details see the Supporting Information). The foregoing is a
27 textbook account of thermodynamic integration; TIES differs by virtue of performing ensemble
29 25 33
30 based calculations, ' in the same manner as for ESMACS, thereby providing control of
32 standard errors based on length of simulation and number of ensembles used. Indeed, it is
37 described by Gaussian random processes, which makes it possible to draw on the theory of
39 stochastic calculus to compute relevant properties reliably.5
42 Simulations. The ESMACS and TIES studies follow the protocol developed recently23-25 (see
44 Supporting Information). The Amber package20 was used for the set-up of the systems and
47 analyses of the results, and the MD package NAMD2.934 was used throughout the equilibration
49 and production runs of all simulations, including the metadynamics. An ensemble simulation was
51 performed for each molecular system with identical atomic coordinates for all replicas. Energy
56 The initial velocities were then generated independently from a Maxwell-Boltzmann distribution
important to note that the energy derivatives (Equation 2) calculated for all X windows are well
minimizations were first performed with heavy protein atoms restrained at their initial positions.
at 50 K, and the systems were heated up to and kept at 300K within 60 ps. A series of
run on the BlueWonder2 supercomputer, an IBM NextScale Cluster (8,640 cores) located in
6 equilibration runs, totalling 2 ns, were conducted, while the restraints on heavy atoms were
8 gradually reduced. Finally, 4 ns production simulations were run for each replica for all
11 ESMACS and TIES simulations. In ESMACS, we used 25 replicas in an ensemble calculation
13 for each ligand. In TIES, 5 replicas were used for each pair of the ligands. Simulation of each
24,27,30,35
15 replica consists of 2 ns for equilibration and 4 ns for production run. Previous studies ' ' '
18 have shown that the combination of the simulation length and the size of the ensemble provides a
20 trade off between computational cost and precision.
23 Two independent ESMACS runs were performed to assess the reproducibility of the
25 calculations. One was run on ARCHER, a Cray XC30 supercomputer (equipped with ca 118,000
27 cores), the UK's National High Performance Computing Service located in Edinburgh; the other
32 Science and Technology Facilities Council's (STFC's) Hartree Centre. A comparison between
34 the two runs is very instructive, both in terms of performance and time to solution as well as
37 providing an opportunity to conduct a valuable reproducibility study. Thanks to the large number
39 of cores on ARCHER, we ran the entire set of ESMACS binding affinity calculations
41 concurrently and had the findings available within 7 hours. On BlueWonder2, we were hampered
44 by a number of issues, resulting in the calculations taking significantly longer than on ARCHER.
46 On BlueWonder2 and ARCHER, about 440,000 and 335,000 core hours were consumed in the
51 dynamics simulations and 140,000 and 110,000 core hours for free energy calculations,
53 respectively.
course of this study, including 300,000 and 225,000 core hours for production molecular
10 11 12
20 21 22
The TIES study was performed entirely on ARCHER. TIES simulations of seven out of twelve pairs were packed together (the currently enforced limit of maximum tasks for one job in ARCHER) and submitted as one single job, requiring 43,680 cores, which has high priority and executes rapidly on the machine (waiting time in the ARCHER queue was around 8 hours). This made it possible to complete the entire TIES study of seven pairs within one day including queuing time. Further speed up is feasible on supercomputers with hybrid architectures which would permit for example GPU acceleration for some parts of these calculations. 3. RESULTS
To assess the accuracy and precision of the method, we evaluated the binding affinities of the ligands (Table 1) to BRD4, and compared the computed results with experimental data. ESMACS was used for the full set of the ligands, including the stereoisomer which cannot be considered as a perturbation of any others. TIES was applied to three subsets of the ligands, of which each includes congeners with the same net charge.
1 1 _ slope: 1 1 1.00 + 0.05 1 1 1 1 y-
intercept: -0.16 ±0.40
correlation: 0.98 + 0.0! • • /
- • ; • / / * »/ y * - s* m
• / , 1 y • < • H 1 i , I
o Archer o BW2
AGk . (kcal/mol)
binding
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Ligands
Figure 2. Correlation and standard errors of the calculated binding free energies from two independent studies of the BRD4-ligand models performed on BlueWonder2 and ARCHER. a)
4 Correlation of the predictions, including all rotamers, from 1-traj calculations performed on
6 BlueWonder2 (BW2, horizontal axis) and ARCHER (vertical axis). Solid line: regression of the
8 data using means of the calculated free energies; dotted line: 1: 1 ideal regression. b) the averages
11 and their standard errors from the two separate calculations. One rotamer is used for each ligand.
15 Reproducibility. It is well known that many complex systems exhibit sensitive dependence on
18 initial conditions.5'36 The differences of the initial velocities among individual simulations lead to
20 rapid divergence of trajectories. The calculated thermodynamic properties from individual
23 simulations will therefore inevitably differ. Two sets of ESMACS simulations were performed
25 for the complexes independently on ARCHER and BlueWonder2 (see the Computational Section
27 below). The results are compared by linear regression. Figure 2 shows the variances and
32 correlation coefficient r is 0.98±0.01, with no statistical differences between slopes and
34 intercepts of the calculated regression line and the ideal line (x=y) (Figure 2a). The two studies
37 produce consistent results, with an average difference of 0.07±0.10 kcal/mol, and an average
39 unsigned difference of 0.39±0.06 kcal/mol. Only 3 out of 16 predictions do not overlap within
41 their variances (Figure 2b). For the purpose of ranking compound selectivity, the Spearman
44 correlation coefficient (rS) is also calculated, which shows the rankings from the two studies are
46 very close to one another, with a highly significant Spearman correlation of 0.98±0.02. In the
51 based calculations produced very similar results.
correlation of the calculated binding free energies from the two sets of simulations. The
following analyses, only the results from BlueWonder2 simulations are reported. The ARCHER
53 Choosing Rotamers. Some of the ligands (Table 1) would be able to adopt more than one
rotamer in solution. When the ligands are complexed with the protein, they are usually trapped in
a specific rotameric state. The energy barriers between rotamers are sufficiently high, especially for groups occupying the narrow ZA channel (Figure 1), that it is not feasible to achieve equilibrated populations of rotamers using the current protocol. We therefore initiated molecular dynamics simulations of complexes from all possible rotamers for six out of the sixteen ligands. Simulations including different rotamers are only applied to the complexes, as the free ligands in water are able to fully sample all rotamers. The most favourable rotamer is selected (Figure 3) for each ligand in its complex form, based on the lowest ligand energy of the rotamer, which is approximated by (Gligand + AGbinding).
Figure 3. The calculated binding free energies from simulations on BlueWonder2 and ARCHER. The ligands are numbered as per Table 1. Circles with red/blue colours are the results based on studies with different rotamers. The circles with crosses are the final results with selected rotamers which are chosen based on the sum of energies Gligand and AGbinding (see Equation 1). All of the calculated binding free energies are associated with standard errors of less than 1.7 kcal/mol, and are not shown in the figures for reasons of clarity.
Comparison between ESMACS Calculations and Experiments. The predicted binding free energies from the 1-traj approach exhibit a weak Spearman correlation with the experimental data (Figure 4a). In the current study, as commonly employed in pharmaceutical drug
development projects, a series of ligands are investigated for the same protein target. The energy
whole set (Figure 4). It should be noted, however, that although the 2- and 3-traj methods
6 of the protein therefore does not affect the ranking of the binding affinities (Equation 1). Treating
8 it as a constant, the 2-traj approach significantly improves the correlation between the predictions
11 and experimental data (Figure 4b). Incorporating the free energies of the ligands from their
13 individual simulations in free form, the 3-traj approach, further improves the correlation (Figure
15 4c). The most significant improvement appears in the compounds involving modifications at the
18 R2 position of the tetrahydroquinoline (Table 1), for which the 1-traj approach predicts an
20 increase (compound 11 to 12) or no change (compound 14 to 13) in the binding affinities when
23 an ethyl is replaced by an n-propyl. The 2- and 3-traj approaches correctly rank the binding
25 potencies of the two pairs, and indeed of these four compounds. Although the improvement is
27 most evident for the listed four compounds, it is indeed applied to the overall correlation for the
32 significantly improve the correlations between predictions and experimental data, relatively large
34 standard deviations, up to 1.7 kcal/mol, are associated with the calculated free energies. This
37 indicates that the binding free energies must differ by no less than 3.4 kcal/mol to be statistically
39 significant to a 95% confidence level. With the experimental free energy differences ranging
41 from 0.0 to 4.9 kcal/mol, a larger ensemble would be required to reduce the relatively large
44 uncertainties in order to render all the predictions statistically significant. It is important to
46 mention here that experimental binding affinities usually contain errors to a greater extent37 than
51 important as that from simulations to make the comparisons statistically significant.
are implied by Table 1. Generating reproducible binding free energies from experiments is as
10 11 12
20 21 22
. „ 1 -traj _,com „com „com ACj = Cj - CJ - Cj,.
com rec hg
. „2-traj „com , „corru „com
ACj = Cj - (Cj > - Cj,.
com rec hg
AG3-traJ = Gc„m com lig com rec hg
1 rs = 0.29 ± 0.26 1 1
^xr" -
o ®LI' ® 1 ■
L14 1 ® L12 1 1
-'-1-1-
r = 0.78 ±0.14
-4 -12
o Àê /
/ ® / o ®
/ « / ®
AGexp(kcal/mol)
Figure 4. Spearman ranking correlations of the calculated binding free energies and the experimental data from 1 -traj (left panel), 2-traj (centre) and 3-traj (right panel) ESMACS approaches. The equations on the sub-figures indicate the calculations used in each case. The subscripts (com/rec/lig) and the superscripts (com/lig) in the equations indicate the components (complexes, receptor and ligands) and the simulations (complexes and free ligands), respectively. The ligands with modifications at the R2-position of the tetrahydroquinoline are marked with crosses; they are all significantly improved in the 2- and 3-trajectory version. The standard errors, which are 0.19 - 0.34 kcal/mol for the 1 -traj and 1.02 - 1.71 kcal/mol for the 2- and 3-traj approaches, are not shown for reasons of clarity. They are calculated using a bootstrapping method (see Supporting Information). The 2- and 3-traj approaches have similar errors because the energy of the receptor is treated as a constant and hence the uncertainties are dominated by the energies of the complexes.
The improvements are achieved by including the adaptation energies of the receptor and the ligands (Figure 5). While the adaptation energies for the ligands are calculated as the differences between the free energies in their bound and free states, those for the receptor are calculated relative to the average of the receptor free energies in the bound states. The absolute adaptation energies of protein for each ligand binding would require a converged ensemble simulation of
10 11 12
20 21 22
protein in solvent, preferably initiated from structures of the protein in its free state. As the ranking of the ligand binding is the main concern of the study, and the relative adaptation energies of the protein are energetically as informative as the absolute ones in the binding affinity comparison, no attempt has been made here to compute accurate absolute adaptation energies for the protein (all adaptation energies for the protein are relative ones in the current paper). Unfavourable adaptation energies are usually induced within the protein by the introduction of larger functional groups at R2 and R3 positions of the compounds (Figure 5c). It is interesting to note that, for the ligands with higher binding affinities (AGbinding < -8 kcal/mol), only small adaptation energies arise for both the receptor and the ligands, which makes the binding apparently closer to a "lock and key" recognition mechanism; while for the ligands with lower binding affinities, non-negligible adaptation energies are found for both the receptor and the ligands (Figure 5b). The one exception is the stereoisomer ligand 16 (Table 1) whose binding induces the smallest relative adaptation energy of the receptor but the largest adaptation energy for the ligand (Figure 5c). This compound was included in the set as a negative control to investigate the ability of the computational procedure to highlight the weak activity of the enantiomer of compound 10. Compound 16 has a very high ligand strain, suggesting that this isomer has the wrong shape for the binding site and must adopt a high energy conformation to interact with the protein.
AO"'"- >MV'*'
AG^^AG3"^
o o g o
AG J => AG"
AG" ' => AG
o 3 -
o o ° O O 0 o oo o o
--¿5-0S
°° O o
-6 -4 -12 -10 -8 -6
AG (kcal/mol)
-12 -10
-4 -12 -10 AG (kcal/mol)
10 11 12
20 21 22
(c) «
AG J => AG
AG 1 => AG
5 6 7 8 9 10 II 12 13 14 15 16 Ligands
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Ligands
Figure 5. Improvement of the predictions by inclusion of the adaptation free energies of the receptor and the ligands. a) the binding free energy changes between the 1-traj (black circles) and 2-traj (magenta circles) indicate the relative adaptation energies of the receptor; those between the 2-traj (magenta circles) and 3-traj (orange circles) show the adaptation energies of the ligands. The adaptation energies can be seen more clearly in b) as a function of binding affinities, and in c) for each ligand.
The ESMACS predictions are precise and accurate in terms of ranking the binding free energies, and are so in a reproducible way. However, ESMACS does not provide accurate absolute free energies. In our study, the 3-traj ESMACS protocol yields a much better correlation than the 1-traj ESMACS approach, while they have similar mean absolute deviations (1.97 ± 1.33 kcal/mol and 1.74 ± 0.29 kcal/mol, respectively). The mean absolute deviations can be significantly larger for more flexible ligands such as peptides binding to major histocompatibility complex (MHC)24 Some studies shown that improved binding free energy prediction could be obtained by including important water molecules between ligands and protein38. Inclusion of the bridging water molecules (Figure 1) in the current ESMCA study, however, does not improve the correlation between the calculated binding free energies and the experimental measurements (see Figure S3 in the Supporting Information). Indeed, ours is a "generic" methodology which treats solvent implicitly and works well as such. Embarking on any approaches for treating water
10 11 12
20 21 22
molecules explicitly would lead to a loss of simplicity in the method without conferring any benefit.
Component analyses of the binding free energies could provide insight into the mechanism of compound binding. In our study of peptide-MHC binding24, for example, the van der Waals interaction was shown to be the dominant component and to manifest a good correlation with the experimental binding free energies in the 3-trajectory ESMACS study. In the current molecular systems, however, there are no significant correlations between any energy component and the experimental data, except the bonded (including bond stretching, angle bending, torsions and improper torsions) and non-bonded (van der Waals, electrostatic and solvation) interactions as shown in Figure 6. The sum of the two interactions improves the correlation further, with a similar correlation coefficient as shown in Figure 4 for the 3-trajectory approach. It can therefore be deduced that the conformational entropy component does not contribute to the quality of the ranking prediction. As for the impact of the configurational entropy on the ranking of calculated
27 39 24
binding affinities, different conclusions - improving, worsening or having no effect - have been drawn for diverse protein-ligand complexes. The differences may indeed reflect the mechanism of protein-ligand binding which can be driven predominantly by enthalpy, by entropy, or by both.40
measurements. The TIES approach is both accurate and precise, yielding a Spearman ranking
4 Figure 6. Correlations of free energy components and the experimental data from 3-traj
6 approaches. Both bonded and non-bonded energy terms contribute to the ranking of binding
8 affinities. Their combination (the MMPBSA energy) exhibits a better correlation with
11 experimental data than the components themselves.
14 Comparison between TIES Calculations and Experiments. Good agreement is found for
16 the binding free energy differences between the TIES calculations and the experimental
21 coefficient of 0.92 for the means of calculated and experimental binding free energy differences.
23 The bootstrap analyses give a Spearman coefficient of 0.86±0.15, drawn from a non-normal
26 distribution of resampling correlation coefficients (see Figure S4 in SI). An average difference of
28 0.06±0.26 kcal/mol, and an average unsigned difference of 0.75±0.14 kcal/mol, are found
31 between the predictions and experimental measurements (Figure 7). It should be noted that
33 experimental binding affinities may have statistical errors of around 24% from isothermal
35 titration calorimetry (ITC)37 which is a more reliable method to measure binding affinities than
40 attributed to these experimental errors. The regression line is close to an ideal 1:1 regression,
42 with a slope of 0.93±0.24 and an intercept of -0.10±0.33 kcal/mol. A similar level of accuracy
45 was recently reported16 in an absolute binding free energy calculation for a set of drug-like
47 molecules binding to BRD4 albeit little attention was paid to reproducibility in that work.
the IC50 measurements used in this study. The variances observed in the study can be partially
10 11 12
20 21 22
Figure 7. Correlations of the calculated binding free energy differences from the TIES study and from experimental measurement. The standard error bars from the TIES calculations are all no greater than 0.2 kcal/mol.
4. DISCUSSION
In order to analyse the results it is necessary to consider the differences between the 1-, 2- and 3-trajectory methods. The 1-trajectory simulation calculates the AGbinding assuming there is no energetic penalty for the adaptation of the protein or ligand to the binding conformation, the 2-trajectory method takes the change in protein energy into consideration and the 3-trajectory method accounts for the changes in both protein and ligand energy.
Figure 8. Calculated vs experimental binding free energies for ligands 1, 8 and 9 which are
The test set for this evaluation was designed to investigate the ability of ESMACS to predict
that, in the 1-traj simulation, compound 1 is predicted to be more potent than compounds 8 and
6 labelled in the 1-traj sub-figure.
10 11 12
14 specific SAR trends. The first trend is the loss of potency observed as the hydrophobic part of the
16 acetyl lysine mimetics is grown from methyl (1) to ethyl (8) to isopropyl (9). Figure 8 shows
21 9, but no difference is predicted between 8 and 9. This is in contrast to the experimental results
23 which show a modest loss of potency going from methyl to ethyl, but a substantial loss going
26 from ethyl to isopropyl. On comparing the same 3 compounds in the 2-traj simulation, the
28 predicted affinity of compound 9 is substantially reduced giving the correct ranking of the three
31 compounds. The predictions for the three compounds are qualitatively the same for the 3-traj
33 calculation. These results indicate that the low activity of compound 9 is the result of the
35 increased bulk of the isopropyl group causing strain in the system, which is manifested in an
40 energy of the whole dataset (Figure 6c). The TIES study gives the same ranking for the binding
42 affinities as those from the 2- and 3-traj ESMACS, but with an enhanced correlation coefficient
47 is rigid and so the strain caused by the increased bulk cannot be accommodated by the ligand and
49 hence is transmitted to the protein.
52 The second SAR trend of interest is the effect of increasing bulk at the R2-position of the
54 tetrahydroquinoline (THQ). Similar to the previous instance, increasing the size of the
56 substituent at this position from methyl (10) to ethyl (11) to n-propyl (12) results in a sequential
increased internal energy of the protein. In fact, compound 9 has the highest protein internal
compared to that of ESMACS. Structurally, this result makes sense because the isopropyl group
loss of activity. In this case the 1-traj method inverts the rank order of the three compounds
ESMACS methods, indicating the loss in potency for this bulky group is due to a combination of
6 (Figure 9), relative to experiment, predicting the n-propyl to be the tightest binder. The 2-traj
8 calculation predicts a decrease of the potency of compound 12, consequently predicting it to be
11 substantially weaker than compound 11, as observed experimentally. In the 3-traj method the
13 potency of compound 12 drops further, relative to the other compounds. In all three ESMACS
15 methods, the potency of compound 10 is under-predicted relative to the other compounds; the
18 reason for this is unclear. The TIES study also predicts the binding of compound 12 to be the
20 weakest in this sub-group, but could not distinguish the affinities of the compounds 10 and 11. A
23 second pair of compounds also has the ethyl (14) to n-propyl (13) modification at the same
25 position, and qualitatively the results parallel those described above. In this case, the n-propyl
27 compound shows a correctly predicted loss of activity in both the TIES and 2- and 3-traj
32 protein and ligand strain (Figure 6c). Again, this is consistent with the molecular structure of
34 compound 12, because the n-propyl group is more flexible than the isopropyl group in the first
37 example and consequently is more able to adopt a slightly strained ligand conformation.
39 The carbonyl group in all compounds forms a hydrogen bond with a conserved asparagine
41 Asn140 in the BRD4. This is a key interaction observed in bromodomain-ligand complexes as it
44 mimics the interaction made by the carbonyl of the acetyl lysine of the substrate. The
46 occupancies of this hydrogen bond in all of the compounds are very similar, and do not appear to
51 motif, as described in the previous cases, the loss of potency does not result from the disruption
53 of this hydrogen bond. Conversely, the interaction is maintained at the expense of creating
correlate with their binding affinities. Hence, when bulk is increased adjacent to this key binding
internal strain in the system. This is consistent with the observation that this is a ubiquitous interaction in this protein and a key binding pharmacophore.
l-traj ESMACS 2-traj ESMACS 3-traj ESMACS TIES
AGei (kcal/mol) AGex (kcal/mol) AGex (kcal/mol) AGex (kcal/mol)
Figure 9. Calculated vs experimental binding free energies for ligands 10, 11 and 12 which are labelled in the 1-traj sub-figure.
For compounds 13, 14 and 15, there is an extra hydrogen bond formed between the carboxylate group of the ligands and the Lys 91 of the protein. The replacement of n-propyl (13) with ethyl (14) and methyl (15) slightly changes the orientations of the compounds, making the carboxylate group better aligned for the formation of hydrogen bonds. The occupancies of this hydrogen bond are 16.6%, 31.7% and 33.3% for the compounds 13, 14 and 15, respectively (the cut-off values are 3Á and 120° for the donor-acceptor distance and donor-hydrogen-acceptor angle). The occupancies of the hydrogen bonds for these 3 compounds appear to correlate with their binding affinities (-7.41, -9.20 and -10.51 kcal/mol). However, the energetic contribution of this solvent exposed interaction is difficult to quantify, and the energetic analysis elsewhere in this section suggests differences in intermolecular interactions in the complex contribute modestly to the differences in binding free energy compared to differences in internal energies. Further evidence for this is seen in the pairs of compounds 11 and 14, 12 and 13, which are identical except for a carboxylic acid to piperidine (base) modification in this region. These pairs
10 11 12
20 21 22
of compounds have very similar binding affinities both in terms of experimental and predicted values.
These data suggest that the internal energy of the protein and ligand are major contributors to the relative binding affinities of this series of compounds. The internal energy contributions are exemplified by plotting the difference between the binding free energies predicted using the different ESMACS methods (Figure 10).
Figure 10. Correlations of the internal energy contributions to the calculated binding free energies and experimental measurement. The internal energy changes are calculated as the differences of the binding free energies between those from the 1-traj and 3-traj approaches:
_ A¡-.3-traj 1binding
AAGcalc = AGbind;ng - AGbndZg.
-,1-traj 1 binding.
Figure 10 plots the difference between the AGbinding calculated by the 3-traj method and the 1-traj method against the experimental AGbinding. This is a plot of the internal energy of the protein and ligand against the experimental AGbinding; it shows a strong correlation (r = 0.87±0.13). If this is contrasted with the 1 trajectory method - which essentially contains the intermolecular contributions to AGbinding, Figure 4a (r = 0.29±0.26) - one can conclude that, for this protein and
10 11 12
20 21 22
set of ligands, the internal energy of the ligand and the protein differentiate between the ligands to a far greater extent than protein-ligand interaction energies.
(a)»r H (b)r (c)"
Figure 11. (a) The binding free energies obtained from 3-traj ESMACS calculation (AGesmacs), (b) the binding free energy differences from TIES calculation (AAGTIES) and the binding free energies derived from IC50 data (AGexp). AAGTIES for the subset 1 (ligands 1, 2, 3, 5, 6, 7, 8 and 9), which are relative to ligand 3 (black cross), are shown in black circles; those for the subset 2 (ligands 10, 11 and 12), which are relative to ligand 11 (red cross), are shown in red circles; and those for subset 3 (ligands 13, 14 and 15), which are relative to ligand 14 (blue cross), are shown in blue circles.
5. CONCLUSIONS
The Binding Affinity Calculator (BAC) software environment has been used to run ESMACS and TIES calculations against a series of inhibitors of BRD4 from the THQ chemical series. High performance computing was used in an automated workflow to build models, run multiple replicas of calculations and analyse the output, placing reliable standard error bounds on predicted binding affinities.
Despite the challenging dataset used for this evaluation, good agreement is found for the binding free energy differences between the experimental measurements and the theoretical
4 calculations, for both ESMACS and TIES. ESMACS is good as an "absolute" method for
6 ranking an arbitrary set of ligands, which may be very diverse in terms of structures and
8 electronic properties. To produce good rankings it is necessary in this case to use a 3-trajectory
11 version, which has the further benefit of clearly distinguishing between ligands that adhere more
13 closely to the "lock and key" or to the "induced fit" binding hypotheses. TIES performs well in
15 determining relative binding free energies of congeners when there is no change of net charge,
18 even when there are differences in structural features between pairs of compounds which one
20 would not intuitively regard as minor.
23 Overall, ESMACS provides reliable binding affinity rankings and clear mechanistic insight
25 into what factors drive binding processes in individual ligand-protein systems. TIES offers more
27 quantitative accuracy in its predictions but, owing to its alchemical basis which necessarily keeps
32 mechanistic insights into binding. The standard errors in ESMACS and TIES are fully controlled
34 through simulation length and number of replicas used in the ensembles chosen24'27.
37 The results also offer insight into possible design strategies for BRD4 ligands. In particular it
39 has been noted that, within a chemical series, differences in binding affinity do not seem to result
41 from differences in protein-ligand interaction energies. On the contrary, to avoid unfavourable
44 (in this case, mainly steric) interactions both ligands and protein adopt strained states resulting in
46 large adaptation energies. Hence, during the optimization of a chemical series for BRD4 it would
51 cavity. It should be noted, however, that this study did not include molecules in which
53 unfavourable electrostatic protein-ligand interactions were introduced.
two congeneric ligands in play at all times, is less able to provide similar structural and
be prudent to carefully consider the shape complementarity of the ligand and the active site
1083-1-191), the MRC Medical Bioinformatics project (MR/L016311/1), the EU H2020
4 The results offer encouragement that BAC, and its underlying approach of running multiple
6 replicas of each simulation to improve predictive power , can accurately rank ligands by their
8 binding free energy. Further evaluations are planned to confirm this potential. As large and
11 secure computing resources become more routinely available, for example through cloud
13 computing, it will become increasingly easy for industrial groups to access approaches like the
15 one outlined in this study. Consequently, the robust prediction of protein-ligand binding affinities
18 in an industrial setting should become more routine and offer a long awaited development in the
20 field of structure based design.
22 ACKNOWLEDGEMENT
25 The authors would like to acknowledge the support of EPSRC via the 2020 Science
27 programme (http://www.2020science.net/, EP/I017909/1), the Qatar National Research Fund (7-
32 CompBioMed (675451) and ComPat (671564) projects, and funding from the UCL Provost. We
34 acknowledge use of Hartree Centre computer BlueWonder2 in this work and assistance from its
37 scientific support staff. The STFC Hartree Centre is a research collaboratory in association with
39 IBM providing High Performance Computing platforms funded by the UK's investment in e-
41 Infrastructure. We made use of ARCHER, the UK's national High Performance Computing
44 Service, funded by the Office of Science and Technology through EPSRC's High-End
46 Computing Programme. Access to ARCHER was provided through the 2020 Science
51 access to Prometheus, the fastest Polish supercomputer run by ACK Cyfronet AGH in Krakow,
53 was provided.
57 Supporting Information
programme. This research was partially supported by the PLGrid Infrastructure through which
4 Detailed description of the methods used, additional energetic analyses of the metadynamics,
6 free energy calculation with inclusion of explicit water molecules, alongside the atomic
8 coordinates of the compound-protein complexes and experimental data on compound binding.
11 This information is available free of charge via the Internet at http://pubs.acs.org.
14 REFERENCES
16 1. Lounnas, V.; Ritschel, T.; Kelder, J.; McGuire, R.; Bywater, R. P.; Foloppe, N., Current 18 progress in structure-based rational drug design marks a new mindset in drug discovery. Comput. -19 Struct. Biotechnol. J. 2013, 5, e201302011.
20 2. Green, D. V.; Leach, A. R.; Head, M. S., Computer-aided molecular design under the
21 SWOTlight. J. Comput. AidedMol. Des. 2012, 26 (1), 51-56.
22 3. Kuhn, B.; Guba, W.; Hert, J.; Banner, D.; Bissantz, C.; Ceccarelli, S.; Haap, W.; Korner,
24 M.; Kuglstatter, A.; Lerner, C.; Mattei, P.; Neidhart, W.; Pinard, E.; Rudolph, M. G.; Schulz-
25 Gasch, T.; Woltering, T.; Stahl, M., A real-world perspective on molecular design. J. Med.
26 Chem. 2016, 59 (9), 4087-4102.
27 4. Wang, L.; Wu, Y.; Deng, Y.; Kim, B.; Pierce, L.; Krilov, G.; Lupyan, D.; Robinson, S.;
28 Dahlgren, M. K.; Greenwood, J.; Romero, D. L.; Masse, C.; Knight, J. L.; Steinbrecher, T.;
29 Beuming, T.; Damm, W.; Harder, E.; Sherman, W.; Brewer, M.; Wester, R.; Murcko, M.; Frye,
31 L.; Farid, R.; Lin, T.; Mobley, D. L.; Jorgensen, W. L.; Berne, B. J.; Friesner, R. A.; Abel, R.,
32 Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery
33 by way of a modern free-energy calculation protocol and force field. J. Am. Chem. Soc. 2015,
34 137 (7), 2695-2703.
35 5. Coveney, P. V.; Wan, S., On the calculation of equilibrium thermodynamic properties
37 from molecular dynamics. Phys. Chem. Chem. Phys. 2016, 18 (44), 30236-30240.
38 6. Arrowsmith, C. H.; Bountra, C.; Fish, P. V.; Lee, K.; Schapira, M., Epigenetic protein
39 families: a new frontier for drug discovery. Nat Rev Drug Discov 2012, 11 (5), 384-400.
40 7. Copeland, R. A., Epigenetic medicinal chemistry. ACS Med. Chem. Lett. 2016, 7 (2),
41 124-127.
42 8. Copeland, R. A.; Olhava, E. J.; Scott, M. P., Targeting epigenetic enzymes for drug
44 discovery. Curr. Opin. Chem. Biol. 2010, 14 (4), 505-510.
45 9. Theodoulou, N. H.; Tomkinson, N. C.; Prinjha, R. K.; Humphreys, P. G., Clinical
46 progress and pharmacology of small molecule bromodomain inhibitors. Curr. Opin. Chem. Biol.
47 2016, 33, 58-66.
48 10. Bamborough, P.; Diallo, H.; Goodacre, J. D.; Gordon, L.; Lewis, A.; Seal, J. T.; Wilson, 40 D. M.; Woodrow, M. D.; Chung, C. W., Fragment-based discovery of bromodomain inhibitors
51 part 2: optimization of phenylisoxazole sulfonamides. J. Med. Chem. 2012, 55 (2), 587-596.
52 11. Chung, C. W.; Coste, H.; White, J. H.; Mirguet, O.; Wilde, J.; Gosmini, R. L.; Delves, C.;
53 Magny, S. M.; Woodward, R.; Hughes, S. A.; Boursier, E. V.; Flynn, H.; Bouillot, A. M.;
54 Bamborough, P.; Brusq, J. M.; Gellibert, F. J.; Jones, E. J.; Riou, A. M.; Homes, P.; Martin, S.
55 L.; Uings, I. J.; Toum, J.; Clement, C. A.; Boullay, A. B.; Grimley, R. L.; Blandel, F. M.;
57 Prinjha, R. K.; Lee, K.; Kirilovsky, J.; Nicodeme, E., Discovery and characterization of small
58 molecule inhibitors of the BET family bromodomains. J. Med. Chem. 2011, 54 (11), 3827-3838.
4 12. Chung, C. W.; Dean, A. W.; Woolven, J. M.; Bamborough, P., Fragment-based discovery
5 of bromodomain inhibitors part 1: inhibitor binding modes and implications for lead discovery.
6 J. Med. Chem. 2012, 55 (2), 576-586.
7 13. Sadiq, S. K.; Wright, D.; Watson, S. J.; Zasada, S. J.; Stoica, I.; Coveney, P. V.,
8 Automated molecular simulation based binding affinity calculator for ligand-bound HIV-1
9 proteases. J. Chem. Inf. Model. 2008, 48 (9), 1909-1919.
11 14. Wyce, A.; Ganji, G.; Smitheman, K. N.; Chung, C. W.; Korenchuk, S.; Bai, Y.; Barbash,
-|2 O.; Le, B.; Craggs, P. D.; McCabe, M. T.; Kennedy-Wilson, K. M.; Sanchez, L. V.; Gosmini, R.
13 L.; Parr, N.; McHugh, C. F.; Dhanak, D.; Prinjha, R. K.; Auger, K. R.; Tummino, P. J., BET
14 inhibition silences expression of MYCN and BCL2 and induces cytotoxicity in neuroblastoma
15 tumor models. PLoS One 2013, 8 (8), e72967.
17 15. Gosmini, R.; Nguyen, V. L.; Toum, J.; Simon, C.; Brusq, J. M.; Krysa, G.; Mirguet, O.;
18 Riou-Eymard, A. M.; Boursier, E. V.; Trottet, L.; Bamborough, P.; Clark, H.; Chung, C. W.;
19 Cutler, L.; Demont, E. H.; Kaur, R.; Lewis, A. J.; Schilling, M. B.; Soden, P. E.; Taylor, S.;
20 Walker, A. L.; Walker, M. D.; Prinjha, R. K.; Nicodeme, E., The discovery of I-BET726
21 (GSK1324726A), a potent tetrahydroquinoline ApoA1 up-regulator and selective BET
23 bromodomain inhibitor. J. Med. Chem. 2014, 57 (19), 8111-8131.
24 16. Aldeghi, M.; Heifetz, A.; Bodkin, M. J.; Knapp, S.; Biggin, P. C., Accurate calculation of
25 the absolute free energy of binding for drug molecules. Chem. Sci. 2016, 7 (1), 207-218.
26 17. Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.;
27 Shaw, D. E., Improved side-chain torsion potentials for the Amber ff99SB protein force field.
28 Proteins 2010, 78 (8), 1950-1958.
20 18. Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A., Development and
31 testing of a general amber force field. J. Comput. Chem. 2004, 25 (9), 1157-1174.
32 19. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman,
33 J. R.; Montgomery, J. J.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.;
34 Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.;
36 Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima,
37 T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.;
38 Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.;
39 Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador,
40 P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.;
41 Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul,
43 A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.;
44 Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara,
45 A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.;
46 Pople, J. A., Gaussian 03. Gaussian, Inc.: 2004.
47 20. Case, D. A.; Cheatham, T. E., 3rd; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M., Jr.;
49 Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J., The Amber biomolecular simulation
50 programs. J. Comput. Chem. 2005, 26 (16), 1668-1688.
51 21. Laio, A.; Gervasio, F. L., Metadynamics: a method to simulate rare events and
52 reconstruct the free energy in biophysics, chemistry and material science. Rep. Prog. Phys. 2008,
53 71 (12), 126601.
54 22. Lin, Y. L.; Aleksandrov, A.; Simonson, T.; Roux, B., An overview of electrostatic free
56 energy computations for solutions and proteins. J. Chem. Theory Comput. 2014, 10 (7), 269057 2709.
4 23. Wan, S.; Bhati, A. P.; Zasada, S. J.; Coveney, P. V., Rapid, accurate, precise and
5 reproducible ligand-protein binding free energy prediction. 2016, preprint.
6 24. Wan, S.; Knapp, B.; Wright, D. W.; Deane, C. M.; Coveney, P. V., Rapid, precise, and
7 reproducible prediction of peptide-MHC binding affinities from molecular dynamics that
8 correlate well with experiment. J. Chem. Theory Comput. 2015, 11 (7), 3346-3356.
90 25. Bhati, A. P.; Wan, S.; Wright, D. W.; Coveney, P. V., Rapid, accurate, precise and
,|1 reliable relative free energy prediction using ensemble based thermodynamic integration. J.
-|2 Chem. Theory Comput. 2016, DOI: 10.1021/acs.jctc.6b00979
13 26. Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.;
14 Duan, Y.; Wang, W.; Donini, O.; Cieplak, P.; Srinivasan, J.; Case, D. A.; Cheatham III, T. E.,
15 Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 2000, 33 (12), 889-897.
27. Wright, D. W.; Hall, B. A.; Kenway, O. A.; Jha, S.; Coveney, P. V., Computing
19 Clinically Relevant Binding Free Energies of HIV-1 Protease Inhibitors. J. Chem. Theory
20 Comput. 2014, 10 (3), 1228-1241.
21 28. Genheden, S.; Ryde, U., The MM/PBSA and MM/GBSA methods to estimate ligand-
22 binding affinities. Expert Opin. Drug Discov. 2015, 10 (5), 449-461.
24 29. Swanson, J. M.; Henchman, R. H.; McCammon, J. A., Revisiting free energy
25 calculations: a theoretical connection to MM/PBSA and direct calculation of the association free
26 energy. Biophys. J. 2004, 86, 67-74.
27 30. Genheden, S.; Ryde, U., How to obtain statistically converged MM/GBSA results. J.
28 Comput. Chem. 2010, 31 (4), 837-846.
20 31. Luo, R.; David, L.; Gilson, M. K., Accelerated Poisson-Boltzmann calculations for static
31 and dynamic systems. J. Comput. Chem. 2002, 23 (13), 1244-1253.
32 32. Beveridge, D. L.; Dicapua, F. M., Free-energy via molecular simulation - applications to
33 chemical and biomolecular systems. Annu. Rev. Biophys. Biophys. Chem. 1989, 18, 431-492.
34 33. Bunney, T. D.; Wan, S.; Thiyagarajan, N.; Sutto, L.; Williams, S. V.; Ashford, P.; Koss, 36 H.; Knowles, M. A.; Gervasio, F. L.; Coveney, P. V.; Katan, M., The effect of mutations on drug 3367 sensitivity and kinase activity of fibroblast growth factor receptors: A combined experimental
38 and theoretical study. EBioMedicine 2015, 2 (3), 194-204.
39 34. Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.;
40 Skeel, R. D.; Kale, L.; Schulten, K., Scalable molecular dynamics with NAMD. J. Comput.
41 Chem. 2005, 26 (16), 1781-1802.
43 35. Sadiq, S. K.; Wright, D. W.; Kenway, O. A.; Coveney, P. V., Accurate ensemble
44 molecular dynamics binding free energy ranking of multidrug-resistant HIV-1 proteases. J.
45 Chem. Inf. Model. 2010, 50 (5), 890-905.
46 36. Norman, G. E.; Stegailov, V. V., Stochastic theory of the classical molecular dynamics
47 method. Math. Models Comput. Simul. 2013, 5 (4), 305-333.
49 37. Chodera, J. D.; Mobley, D. L., Entropy-enthalpy compensation: role and ramifications in
50 biomolecular ligand recognition and design. Annu. Rev. Biophys. 2013, 42, 121-142.
51 38. Zhu, Y. L.; Beroza, P.; Artis, D. R., Including explicit water molecules as part of the
52 protein structure in MM/PBSA calculations. J. Chem. Inf. Model. 2014, 54 (2), 462-469.
53 39. Rastelli, G.; Del Rio, A.; Degliesposti, G.; Sgobba, M., Fast and accurate predictions of
54 binding free energies using MM-PBSA and MM-GBSA. J. Comput. Chem. 2010, 31 (4), 79755 810.
4 40. Reynolds, C. H.; Holloway, M. K., Thermodynamics of ligand binding and efficiency.
5 ACS Med. Chem. Lett. 2011, 2 (6), 433-437.
10 11 12
20 21 22
10 11 12
20 21 22
Table of Contents graph
Rapid and Reliable Binding Affinity Prediction for Analysis of Bromodomain Inhibitors: a Computational Study
Shunzhou Wan,a Agastya P. Bhati,a Stefan J. Zasadaa, Ian Wall,b Darren Green,b Paul Bamborough,b and Peter V. Coveneya*
