computational 
 chemistry has been an established tool in drug discovery 
 for a number of years .
the number of crystal structures available 
 in the public domain and within pharmaceutical companies on which 
 to base computational studies continues to rise rapidly . despite the 
 increase in resources applied and in experimental data on which to 
 base the studies ,
industrial structure - based design approaches have 
 evolved very little in recent times . in 
 particular , the approaches are largely qualitative and largely dependent 
 on the experience and knowledge of the practitioner . attempts to quantify protein  ligand binding affinities are 
 rare .
expert practitioners have little confidence in existing tools 
 to make robust predictions and certainly not to do so on a time scale 
 that can substantially impact medicinal chemistry programmes .
recently , there has been a renewed interest in the use of free 
 energy calculations in drug discovery programmes . in particular , the 
 fep+ implementation of free energy perturbation ( fep ) has shown potential 
 to improve the ability to predict protein  ligand binding affinities 
 on an industrially relevant time scale .
research is ongoing to understand how broadly applicable the method 
 is , and how accurate its predictions are when applied to active drug 
 discovery programmes .
although fep+ applies replica exchange 
 solute tempering ( rest ) 
 in which exchange moves are made between different  windows , 
 its predictions , like those from many other approaches , are generated 
 from a single output for each pair of mutations .
advances in high - end 
 computing capabilities offer the opportunity to run vast numbers of 
 calculations in parallel .
the application of these computational capabilities 
 to free energy calculations allows results to be returned quickly 
 and multiple replicas of simulations to 
 be run , leading to tighter control of standard 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 process is to validate the performance 
 on an industrially relevant data set .
depending on the reliability , 
 rapidity , and throughput of these 
 calculations , they might find application at various stages of the 
 drug discovery and development process across the wider pharmaceutical 
 industry . as these methods require significant compute resources beyond 
 existing in - house industrial capacity , assessment ( and any subsequent 
 adoption ) of the methodology
research into epigenetic proteins is currently a 
 major and rapidly 
 evolving focus for the pharmaceutical industry .
bromodomain - containing 
 proteins , and in particular the four members of the bet ( bromodomain 
 and extra terminal domain ) family , each containing two bromodomains , 
 have been widely studied .
small molecule inhibitors able to competitively 
 antagonize the binding of acetylated histone tails to these modules 
 have been shown to exert profound effects on gene expression and have 
 shown promising preclinical efficacy in pathologies ranging from cancer 
 to inflammation .
indeed , several compounds are progressing through 
 early stage clinical trials and are showing exciting early results .
most inhibitors reported show similar binding 
 potencies to all bet family bromodomains . a representative inhibitor - protein 
 structure is given in figure 1 , showing the key elements for the inhibitor binding .
this 
 study will concentrate on the first bromodomain of bromodomain - containing 
 protein 4 ( brd4-bd1 ) for which extensive crystallographic and ligand 
 binding data are available .
two views are displayed for the binding mode ( pdb i d : 4bjx(14 ) ) , in which i - bet726 is represented 
 as stick in cyan / blue / red / green , the protein is shown as cartoon in 
 silver , the crystallographic water molecules are shown as red balls , 
 and clipped protein surfaces are shown in orange .
the purpose of the present study is to assess the potential 
 for 
 rapid , accurate , precise , and reproducible binding affinity calculations 
 based on the use of a binding affinity calculator ( bac ) software tool and associated services including 
 a python - based toolkit , fabsim , to automate 
 data transfer and job submission .
the approach is based on the use 
 of high performance computing in an automated workflow which builds 
 models , runs large numbers of replica calculations , and analyzes the 
 output data in order to place reliable standard error bounds on predicted 
 binding affinities .
in this study , chemical structures of 16 brd4-bd1 
 ligands based on a single tetrahydroquinoline ( thq ) template 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 sets of closely related compounds with key sar trends . 

specifically , there are two growth vectors which cause a drop off 
 in potency , a growth vector where substantial structural variation 
 is tolerated and an enantiomeric pair where one isomer is significantly 
 more potent than the other .
the calculations 
 were then performed by the ucl group , initially blind , to investigate 
 the ability of bac to reproduce the sar trends .
the experimental binding 
 affinities were not released to the ucl group until all of the computations 
 had been completed and the estimation of binding free energies were 
 reported ( see figure 3 below ) .
these predictions were subsequently compared with the independently 
 measured experimental free energies of binding ( table 1 ) as part of the assessment of the reliability 
 of the method . compounds 
 19 are electrostatically 
 neutral ,
compounds 1012 and 16 are positively charged , and 
 compounds 1315 are negatively charged .
all compounds are the 2-(s ) 4-(r ) isomers ( figure 1 ) except compound 16 which is 2-(r ) 4-(s ) .
the x - ray 
 crystal structure used to represent brd4-bd1 ( published 
 coordinates pdb i d : 4bjx(14 ) ) is a complex with a thq ( i - bet726 ) chosen in order to reproduce the likely conditions 
 for which free - energy calculations might be used in a real drug discovery 
 program .
the binding mode of i - bet726 is shown in figure 1 . both of the significant protein 
ligand 
 hydrogen - bonding interactions take place through the thq n1-acetamide 
 carbonyl group , which interacts directly with the asn140 side chain 
 of the acetyl - lysine binding site and also through the w1 water to 
 the side chain of tyr97 .
the w1 water forms part of a chain of water 
 molecules buried within the site ( w2w4 , figure 1 ) .
substituents larger than acetamide at 
 the n1-position extend into the water - filled part of the site .
our 
 ligand data set included a small series of increasing size at n1 ( 1 , 
 8 , and 9 , see table 1 ) to probe whether the computer methodology can accurately predict 
 the outcome of growing into this region . the thq r2-position 
 ( s)-methyl 
 group occupies a small lipophilic site between the side chains of 
 three residues of the brd4-bd1 za - loop ( val87 , leu92 , and tyr97 ) . 

carbon  carbon contacts for all three lie within 4.25  , 
 apparently offering little room for extension of this group without 
 some structural reorganization .
our ligand data set includes a small 
 number of analogues exploring larger substituents at the r2-position 
 ( 1014 ) to see whether the simulation can accurately predict 
 the consequences of pushing on these residues .
the thq 6-phenyl 
 substituent fills the narrow  za channel  
 between trp81 and leu92 .
this ring can be substituted at the r3- and 
 r4-positions , but r2-substitution is detrimental because this causes 
 an increase in the inter - ring torsion .
the remainder of 
 our data set includes a variety of substituents probing electrostatic 
 and steric changes at this position .
the thq 4-anilino substituent 
 projects onto the  wpf - shelf  
 subsite outside the acetyl - lysine pocket , close to trp81 .
sar around 
 this ring has been published previously , showing that substitution had small effects on potency , probably 
 because one edge of the aniline ring is solvent - exposed . therefore
, 
 we did not attempt to vary substituents at this position in our data 
 set , preferring to keep it as constant as possible while exploring 
 changes elsewhere . in this series the 2-(s ) 

4-(r ) isomers ( see figure 1 ) are significantly more potent than their 
 alternative trans enantiomers .
we included an enantiomeric 
 pair of 2-(s ) 4-(r ) ( compound 10 ) 
 and 2-(r ) 
 4-(s ) ( compound 16 ) isomers to explore whether the 
 simulations were capable of distinguishing between them .
the 
 ligand - protein complexes were constructed by replacing the 
 ligand in the pdb file with the ligands of interest ( see structures 
 in supporting information ) . for the congeneric 
 compounds studied here , it is plausible to assume that they bind in 
 the same mode in which all key compound - protein interactions are preserved 
 ( figure 1 ) .
for compounds 
 with significant structural differences , computational docking may 
 be required to generate reasonable complex structures .
preparation and setup of the simulations were 
 implemented using bac , including parametrization 
 of the compounds , solvation of the complexes , electrostatic neutralization 
 of the systems by adding counterions and generation of configurations 
 files for the simulations .
the amber ff99sb - ildn force field was used for the protein , and tip3p was used 
 for water molecules .
standard protonation states were assigned to 
 all titratable residues at ph 7 , with histidines protonated on the 
  position ( hie ) .
compound parameters were produced using the 
 general amber force field ( gaff ) with 
 gaussian 03 to optimize compound geometries 
 and to determine electrostatic potentials at the hartree 
the restrained electrostatic potential 
 ( resp ) module in the amber package was 
 used to calculate the partial atomic charges for the compounds .
all 
 systems were solvated in orthorhombic water boxes with a minimum extension 
 from the protein of 14 
. some of the ligands could adopt 
 several rotamers for the r4 group 
 ( table 1 ) in solution , 
 which due to the constrained environment of the active site are unlikely 
 to be sampled in a single bound simulation run with the protocol used 
 ( see simulations subsection ) . to decide 
 which rotamer(s ) with which to start ,
metadynamics simulations were employed which used a history - dependent 
 biasing potential to explore the conformational space of the chosen 
 degrees of freedom , here the rotatable bond(s ) involving groups r2 
 and r4 ( see details in the supporting information ) .
five replicas were used for each metadynamics study to have a 
 reasonable convergence of the potential profile of the chosen dihedral 
 angle(s ) of a ligand in complexed form .
the calculated potential of 
 mean force ( pmf ) was used to determine the most favorable rotamer(s ) 
 from which esmacs simulations initiate . when unambiguous results ensued , 
 the energetically most favorable rotamer was chosen . for some ligands , 

the ligands with multiple rotamers were : 
 4 , 7 , 10 , 11 , 12 , and 16 .
there are two rotamers for each of these , 
 except for 7 for which there are three generated by the flip and twist 
 between the two ring planes .
this resulted in a total number of 23 
 complexes being simulated in the esmacs study . for each ligand 
 with multiple rotamers ,
the energetic properties 
 were analyzed and the most favorable rotamer was chosen as the final 
 result of the esmacs study ( see results section ) 
 and was used as the initial structure in the ties study . for the latter , 
 three subsets of the ligands were selected , within which relative 
 binding free energies were calculated with the ties approach .
the 
 subsets of the ligands were as follows : set 1 including ligands 19 ; 
 set 2 including ligands 10 , 11 , and 12 which are positively charged ; 
 and set 3 including ligands 13 , 14 , and 15 which are negatively charged . 

the division of the full set of the cognate ligands into subsets is 
 necessary because ties , just as any other ti and fep based method , 
 encounters specific difficulties owing to major adjustments in long - range 
 electrostatic interactions for the  congeners  where 
 the net charges change . for each subset , 

ligands were paired based on their similarities , and ties calculations 
 were performed to alchemically mutate one ligand to another ( see theoretical basis below ) .
the ucl group have recently introduced 
 new protocols for binding free energy calculations , termed  enhanced 
 sampling of molecular dynamics with approximation of continuum solvent  

( esmacs ) and  thermodynamic integration 
 with enhanced sampling  ( ties ) .
boltzmann 
 surface area method ( mmpbsa ) which makes 
 a continuum approximation for the aqueous solvent , while ties centers , 
 as the abbreviation indicates , on thermodynamic integration ( ti ) . 

although the approaches are built around the standard mmpbsa and ti 
 methodologies , our abbreviations are used to emphasize the central 
 importance of the ensemble - based nature of the protocols employed 
 as well as , in the case of esmacs the wider generality and flexibility 
 of the protocol .
the size of the statistical mechanical 
 ensembles is determined systematically so that predictions are accurate , 
 precise , and reliable . moreover ,
the term  mmpbsa  
 is used to mean very different things in numerous journal articles and textbooks , including calculations 
 based on single docked structures or on 
 simulation trajectories , calculations with or without the inclusion 
 of configurational entropies , and almost wholly using the so - called 
 1-trajectory approach . in the esmacs protocol , we always mean a fully 
 converged , ensemble - based determination of the free energy of binding 
 from either a one- , two- , or three - trajectory method : this includes 
 both the configurational entropy and the association free energy , and ( where appropriate ) the relative or absolute adaptation energy .
a binding affinity calculator ( bac ) was 
 used to perform esmacs and ties studies .
bac constitutes a computational 
 pipeline built from selected software tools and services ( see supporting information ) , and relies on access 
 to a range of computing resources .
it automates much of the complexity 
 of building , running , and marshalling the molecular dynamics simulations , 
 as well as collecting and analyzing data .
integration and automation 
 are central to the reliability of the method , ensuring that the results 
 are reproducible and can be delivered rapidly . in esmacs
, the 
 free energy is evaluated approximately based on the extended mmpbsa 
 method , including configurational entropy , and the 
 free energy of association .
free energy 
 changes ( gbinding ) are determined 
 for the molecules in their solvated states .
the binding free energy 
 change is then calculated as1where gcomplex , greceptor , and gligand are free 
 energies of the complex , the receptor and the 
 ligands , respectively .
boltzmann calculation is 
 performed here using the amber built - in pbsa solver with radii taken from the parameter / topology files , and 
 the configurational entropy is calculated by a normal mode ( nmode ) 
 method .
the three terms on the right - hand side of eq 1 can be generated from single simulations 
 of the complexes or from separated simulations of complexes , receptor(s ) , 
 and ligand(s ) .
the former is the so - called 1-trajectory ( 1-traj ) method 
 of which the trajectories of the receptor(s ) and the ligand(s ) are 
 extracted from those of the complexes .
the latter is the so - called 
 3-trajectory ( 3-traj ) method of which the trajectories are generated 
 from separate simulations of the three components . in the drug development 
 field , binding is usually investigated for a set of ligands bound 
 to the same protein target .
the free energy of the receptor greceptor will then be the same for all ligands 
 in the 3-traj method . in the current study
, we employ three approaches 
 to calculate the binding free energies gbinding : ( a ) using single simulations of the complexes ( 1-traj ) ; 
 ( b ) same as ( a ) but using the average of the receptor free energies 
 greceptor from the 1-traj 
 method ( denoted as 2-traj ) ; and ( c ) the same as ( b ) but invoking separated 
 ligand simulations for the derivation of gligand ( denoted as 3-traj ) . in ties , an alchemical transformation 
 for the mutated entity , either 
 the ligand or the protein ,
the free energy changes of the 
 alchemical mutation processes , gaqalch and gcomplexalch , are calculated by2and the binding free energy difference 
 is 
 calculated from3here  ( 0 
  
 1 ) is a coupling parameter such that  = 0 and  = 1 correspond 
 to the initial and final thermodynamic states , and v( ) is the potential energy of an intermediate state . 
  denotes an ensemble 
 average over configurations representative of the state . to 
 avoid the well - known  end - point catastrophe  , a soft - core potential is used , along with an 
 approach to decouple at different rates the electrostatic and van 
 der waals interactions of the perturbed atoms with their environment 
 ( for more details see the supporting information ) .
the foregoing is a textbook account of thermodynamic integration ; 
 ties differs by virtue of performing ensemble based calculations , in the same manner as for esmacs , thereby providing control of standard 
 errors based on length of simulation and number of ensembles used . 

indeed , it is important to note that the energy derivatives ( eq 2 ) calculated for all  
 windows are well described by gaussian random processes , which makes 
 it possible to draw on the theory of stochastic calculus to compute 
 relevant properties reliably .
the esmacs and ties studies follow the 
 protocol developed recently ( see supporting information ) .
the amber package was used for the setup of the systems and analyses of the 
 results , and the md package namd2.9 was 
 used throughout the equilibration and production runs of all simulations , 
 including the metadynamics .
an ensemble simulation was performed for 
 each molecular system with identical atomic coordinates for all replicas . 

energy minimizations were first performed with heavy protein atoms 
 restrained at their initial positions .
the initial velocities were 
 then generated independently from a maxwell  boltzmann distribution 
 at 50 k , and the systems were heated up to and kept at 300 k within 
 60 ps .
a series of equilibration runs , totalling 2 ns , were conducted , 
 while the restraints on heavy atoms were gradually reduced .
finally , 
 4 ns production simulations were run for each replica for all esmacs 
 and ties simulations . in esmacs
, we used 25 replicas in an ensemble 
 calculation for each ligand . in ties ,
simulation of each replica consists of 2 ns for 
 equilibration and 4 ns for production run .
previous studies have shown that the combination of the simulation length and the 
 size of the ensemble provides a trade off between computational cost 
 and precision .
one was run on archer , 
 a cray xc30 supercomputer ( equipped with ca 118  000 
 cores ) , the uk s national high performance computing service 
 located in edinburgh ; the other run on the bluewonder2 supercomputer , 
 an ibm nextscale cluster ( 8640 cores ) located in science and technology 
 facilities council s ( stfc s ) hartree centre .
a comparison 
 between the two runs is very instructive , both in terms of performance 
 and time to solution as well as providing an opportunity to conduct 
 a valuable reproducibility study .
thanks to the large number of cores 
 on archer , we ran the entire set of esmacs binding affinity calculations 
 concurrently and had the findings available within 7 h. on bluewonder2 , 
 we were hampered by a number of issues , resulting in the calculations 
 taking significantly longer than on archer .
on bluewonder2 and archer , 
 about 440  000 and 335 000 core hours were consumed in the course 
 of this study , including 300  000 and 225  000 core hours 
 for production molecular dynamics simulations and 140  000 and 
 110  000 core hours for free energy calculations , respectively .
ties simulations 
 of 7 out of 12 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 h ) .
this made it possible to complete the entire ties 
 study of seven pairs within 1 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 .
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 can not 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 .
it is well - known that many complex 
 systems exhibit sensitive dependence on initial conditions .
the differences of the initial velocities among individual simulations 
 lead to rapid divergence of 
 trajectories .
two sets of esmacs simulations 
 were performed for the complexes independently on archer and bluewonder2 
 ( see the computational section above ) .
figure 2 shows the variances and correlation of the 
 calculated binding free energies from the two sets of simulations . 

the correlation coefficient r is 0.98  0.01 , 
 with no statistical differences between slopes and intercepts of the 
 calculated regression line and the ideal line ( x = y ) ( figure 2a ) .
the two studies produce consistent results , with an average difference 
 of 0.07  0.10 kcal / mol , and an average unsigned difference of 
 0.39  0.06 kcal / mol .
only 3 out of 16 predictions do not overlap 
 within their variances ( figure 2b ) . for the purpose of ranking compound selectivity ,
the spearman 
 correlation coefficient ( rs ) is also calculated , 
 which shows that the rankings from the two studies are very close 
 to one another , with a highly significant spearman correlation of 
 0.98  0.02 . in the following analyses ,
the archer - based calculations 
 produced very similar results . 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 ) correlation of the predictions , including 
 all rotamers , from 1-traj calculations performed on bluewonder2 ( bw2 , 
 horizontal axis ) and archer ( vertical axis ) .
solid line , regression 
 of the data using means of the calculated free energies ; dotted line , 
 1:1 ideal regression .
( b ) the averages and their standard errors from 
 the two separate calculations .
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 16 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 favorable 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 + gbinding ) .
circles with red / blue colors are the results 
 based on studies with different rotamers .
the circles with crosses 
 are the final results with selected rotamers which are chosen on the 
 basis of the sum of energies gligand and gbinding ( see eq 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 .
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 
 of the protein therefore does not affect the ranking of the binding 
 affinities ( eq 1 ) .
treating 
 it as a constant , the 2-traj approach significantly improves the correlation 
 between the predictions and experimental data ( figure 4b ) . incorporating the free energies of the 
 ligands from their individual simulations in free form , the 3-traj 
 approach ,
the most significant improvement appears 
 in the compounds involving modifications at the r2 position of the 
 tetrahydroquinoline ( table 1 ) , for which the 1-traj approach predicts an increase ( compound 
 11 to 12 ) or no change ( compound 14 to 13 ) in the binding affinities 
 when an ethyl is replaced by an n - propyl .
the 2- 
 and 3-traj approaches correctly rank the binding potencies of the 
 two pairs , and indeed of these four compounds .
although the improvement 
 is most evident for the listed four compounds , it is indeed applied 
 to the overall correlation for the whole set ( figure 4 ) .
it should be noted , however , that although 
 the 2- and 3-traj methods significantly improve the correlations between 
 predictions and experimental data , relatively large standard deviations , 
 up to 1.7 kcal / mol , are associated with the calculated free energies . 

this indicates that the binding free energies must differ by no less 
 than 3.4 kcal / mol to be statistically significant to a 95% confidence 
 level . with the experimental free energy differences ranging from 
 0.0 to 4.9 kcal / mol
, a larger ensemble would be required to reduce 
 the relatively large uncertainties in order to render all the predictions 
 statistically significant .
it is important to mention here that experimental 
 binding affinities usually contain errors to a greater extent than are implied by table 1 .
generating reproducible binding free energies 
 from experiments is as important as that from simulations to make 
 the comparisons statistically significant .
spearman ranking correlations 
 of the calculated binding free energies 
 and the experimental data from 1-traj ( left panel ) , 2-traj ( center ) , 
 and 3-traj ( right panel ) esmacs approaches .
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.190.34 kcal / mol for the 1-traj 
 and 1.021.71 kcal / mol for the 2- and 3-traj approaches , are 
 not shown for reasons of clarity .
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 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 ) .
unfavorable 
 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 ( gbinding < 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 .
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 panels b as a function of binding affinities , 
 and in panels 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
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 a major histocompatibility 
 complex ( mhc ) .
some studies show that 
 improved binding free energy prediction could be obtained by including 
 important water molecules between ligands and protein .
inclusion of the bridging water molecules ( figure 1 ) in the current 
 esmacs study , however , does not improve the correlation between the 
 calculated binding free energies and the experimental measurements 
 ( see figure s3 in the supporting information ) .
generic  methodology which treats 
 solvent implicitly and works well as such . embarking on any approaches 
 for treating water 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 
 binding , 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 nonbonded ( 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 configurational 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 binding affinities , 
 different conclusions  improving , worsening , or having no effect  have been drawn for diverse protein 
ligand 
 binding which can be driven predominantly by enthalpy , by entropy , 
 or by both .
correlations of free 
 energy components and the experimental data 
 from 3-traj approaches . both bonded and nonbonded
their combination ( the mmpbsa 
 energy ) exhibits a better correlation with experimental data than 
 the components themselves .
good 
 agreement is found for the binding free energy differences between 
 the ties calculations and the experimental measurements .
the ties 
 approach is both accurate and precise , yielding a spearman ranking 
 coefficient of 0.92 for the means of calculated and experimental binding 
 free energy differences .
the bootstrap analyses give a spearman coefficient 
 of 0.86  0.15 , drawn from a non - normal distribution of resampling 
 correlation coefficients ( see figure s4 in the supporting information ) .
an average difference of 0.06  
 0.26 kcal / mol , and an average unsigned difference of 0.75  0.14 
 kcal / mol , are found between the predictions and experimental measurements 
 ( figure 7 ) .
it should 
 be noted that experimental binding affinities may have statistical 
 errors of around 24% from isothermal titration calorimetry ( itc ) which is a more reliable method to measure binding 
 affinities than the ic50 measurements used in this study . 

the variances observed in the study can be partially attributed to 
 these experimental errors .
the regression line is close to an ideal 
 1:1 regression , with a slope of 0.93  0.24 and an intercept 
 of 0.10  0.33 kcal / mol .
a similar level of accuracy 
 was recently reported in an absolute 
 binding free energy calculation for a set of drug - like molecules binding 
 to brd4 , albeit little attention was paid to reproducibility in that 
 work .
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 .
to analyze the results it is necessary to consider the differences 
 between the 1- , 2- and 3-trajectory methods .
the 1-trajectory simulation 
 calculates the gbinding 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 .
the test set for this evaluation was designed to investigate the 
 ability of esmacs to predict specific sar trends .
the first trend 
 is the loss of potency observed as the hydrophobic part of the acetyl 
 lysine mimetics is grown from methyl ( 1 ) to ethyl ( 8) to isopropyl 
 ( 9 ) .
figure 8 shows 
 that , in the 1-traj simulation , compound 1 is predicted to be more 
 potent than compounds 8 and 9 , but no difference is predicted between 
 8 and 9 .
this is in contrast to the experimental results which show 
 a modest loss of potency going from methyl to ethyl , but a substantial 
 loss going from ethyl to isopropyl . on comparing the same three compounds 
 in the 2-traj simulation ,
the predicted affinity of compound 9 is 
 substantially reduced giving the correct ranking of the three compounds . 

the predictions for the three compounds are qualitatively the same 
 for the 3-traj calculation .
these results indicate that the low activity 
 of compound 9 is the result of the increased bulk of the isopropyl 
 group causing strain in the system , which is manifested in an increased 
 internal energy of the protein .
in fact , compound 9 has the highest 
 protein internal energy of the whole data set ( figure 6c ) .
the ties study gives the same ranking 
 for the binding affinities as those from the 2- and 3-traj esmacs , 
 but with an enhanced correlation coefficient compared to that of esmacs . 

structurally , this result makes sense because the isopropyl group 
 is rigid and so the strain caused by the increased bulk can not be 
 accommodated by the ligand and hence is transmitted to the protein .
calculated 
 vs experimental binding free energies for ligands 1 , 
 8 , and 9 which are labeled in the 1-traj subfigure .
the second sar trend of interest is the effect 
 of increasing bulk 
 at the r2-position of the tetrahydroquinoline ( thq ) .
similar to the 
 previous instance , increasing the size of the substituent at this 
 position from methyl ( 10 ) to ethyl ( 11 ) to n - propyl 
 ( 12 ) results in a sequential loss of activity . in this case
the 1-traj 
 method inverts the rank order of the three compounds ( figure 9 ) , relative to experiment , 
 predicting the n - propyl to be the tightest binder . 

the 2-traj calculation predicts a decrease of the potency of compound 
 12 , consequently predicting it to be substantially weaker than compound 
 11 , as observed experimentally . in the 3-traj method
the potency of 
 compound 12 drops further , relative to the other compounds . in all 
 three esmacs methods ,
the potency of compound 10 is under - predicted 
 relative to the other compounds ; the reason for this is unclear .
the 
 ties study also predicts the binding of compound 12 to be the weakest 
 in this subgroup , but could not distinguish the affinities of the 
 compounds 10 and 11 .
a second pair of compounds also has the ethyl 
 ( 14 ) to n - propyl ( 13 ) modification at the same position , 
 and qualitatively the results parallel those described above . in this 
 case , the n - propyl compound shows a correctly predicted 
 loss of activity in both the ties and 2- and 3-traj esmacs methods , 
 indicating the loss in potency for this bulky group is due to a combination 
 of protein and ligand strain ( figure 6c ) .
again , this is consistent with the molecular structure 
 of compound 12 , because the n - propyl group is more 
 flexible than the isopropyl group in the first example and consequently 
 is more able to adopt a slightly strained ligand conformation .
calculated 
 vs experimental binding free energies for ligands 10 , 
 11 , and 12 which are labeled in the 1-traj subfigure .
the carbonyl group in all compounds forms a hydrogen 
 bond with 
 a conserved asparagine asn140 in the brd4 .
this is a key interaction 
 observed in bromodomain - ligand complexes as it mimics the interaction 
 made by the carbonyl of the acetyl lysine of the substrate .
the occupancies 
 of this hydrogen bond in all of the compounds are very similar , and 
 do not appear to correlate with their binding affinities .
hence , when 
 bulk is increased adjacent to this key binding motif , as described 
 in the previous cases , the loss of potency does not result from the 
 disruption of this hydrogen bond .
conversely , the interaction is maintained 
 at the expense of creating 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 . 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 cutoff values 
 are 3  and 120 for the donor 
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 that 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,14 and 12,13 , which 
 are identical except for a carboxylic acid to piperidine ( base ) modification 
 in this region .
these pairs 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 ) .
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 : gcalc = gbinding3-traj  gbinding1-traj .
figure 10 plots 
 the difference between the gbinding calculated by the 3-traj method and the 1-traj method against the 
 experimental gbinding .
this is 
 a plot of the internal energy of the protein and ligand against the 
 experimental gbinding ; 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 gbinding , figure 4a ( r = 0.29  0.26)one can conclude 
 that , for this protein and 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 .
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 analyze the output , placing reliable 
 standard error bounds on predicted binding affinities . despite 
 the challenging data set used for this evaluation ,
good 
 agreement is found for the binding free energy differences between 
 the experimental measurements and the theoretical calculations , for 
 both esmacs and ties .
 
 method for ranking an arbitrary set of ligands , which may be very 
 diverse in terms of structures and electronic properties . to produce 
 good rankings
it is necessary in this case to use a 3-trajectory version , 
 which has the further benefit of clearly distinguishing between ligands 
 that adhere more closely to the  lock and key  or to 
 the  induced fit  binding hypotheses .
ties performs 
 well in determining relative binding free energies of congeners when 
 there is no change of net charge , even when there are differences 
 in structural features between pairs of compounds which one would 
 not intuitively regard as minor .
overall , esmacs provides reliable 
 binding affinity rankings and 
 clear mechanistic insight into what factors drive binding processes 
 in individual ligand  protein systems .
ties offers more quantitative 
 accuracy in its predictions but , owing to its alchemical basis which 
 necessarily keeps two congeneric ligands in play at all times , is 
 less able to provide similar structural and mechanistic insights into 
 binding .
the standard errors in esmacs and ties are fully controlled 
 through simulation length and number of replicas used in the ensembles 
 chosen .
in particular it has been noted that , within a chemical 
 series , differences in binding affinity do not seem to result from 
 differences in protein 
ligand interaction energies . on the 
 contrary , to avoid unfavorable ( in this case , mainly steric ) interactions , 
 both ligands and protein adopt strained states resulting in large 
 adaptation energies .
hence , during the optimization of a chemical 
 series for brd4 it would be prudent to carefully consider the shape 
 complementarity of the ligand and the active site cavity .
it should 
 be noted , however , that this study did not include molecules in which 
 unfavorable electrostatic protein  ligand interactions were 
 introduced .
the results offer encouragement that bac , and its 
 underlying approach 
 of running multiple replicas of each simulation to improve predictive 
 power , can accurately rank ligands by 
 their binding free energy .
as large and secure computing resources become more 
 routinely available , for example through cloud computing , it will 
 become increasingly easy for industrial groups to access approaches 
 like the one outlined in this study .
consequently , the robust prediction 
 of protein  ligand binding affinities in an industrial setting 
 should become more routine and offer a long awaited development in 
 the field of structure - based design .
 binding 
 free energies of bromodomain inhibitors are calculated 
 with recently formulated approaches , namely esmacs ( enhanced sampling 
 of molecular dynamics with approximation of continuum solvent ) and 
 ties ( thermodynamic integration with enhanced sampling ) . 
 a set of 
 compounds is provided by glaxosmithkline , which represents a range 
 of chemical functionality and binding affinities . the predicted binding 
 free energies exhibit a good spearman correlation of 0.78 with the 
 experimental data from the 3-trajectory esmacs , and an excellent correlation 
 of 0.92 from the ties approach where applicable . 
 given access to suitable 
 high end computing resources and a high degree of automation , we can 
 compute individual binding affinities in a few hours with precisions 
 no greater than 0.2 kcal / mol for ties , and no larger than 0.34 and 
 1.71 kcal / mol for the 1- and 3-trajectory esmacs approaches . 
Introduction
Computational Section
Results
Discussion
Conclusions