Scholarly article on topic 'Rapid and Reliable Binding Affinity Prediction of Bromodomain Inhibitors: a Computational Study'

Rapid and Reliable Binding Affinity Prediction of Bromodomain Inhibitors: a Computational Study Academic research paper on "Chemical sciences"

Share paper
Academic journal
J. Chem. Theory Comput.
OECD Field of science

Academic research paper on topic "Rapid and Reliable Binding Affinity Prediction of Bromodomain Inhibitors: a Computational Study"


[ournal of Chemical Theoiy and Computation

THC ll«l*tU4TV Of



Subscriber access provided by University of Newcastle, Australia


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 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:

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


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.


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


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)


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


r = 0.78 ±0.14

-4 -12

o Àê /

/ ® / o ®

/ « / ®


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'*'


o o g o

AG J => AG"

AG" ' => AG

o 3 -

o o ° O O 0 o oo o o


°° 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.


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.


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.


25 The authors would like to acknowledge the support of EPSRC via the 2020 Science

27 programme (, 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


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*