JComputAidedMolDes(2017)31:29–44 DOI10.1007/s10822-016-9956-6 A combined treatment of hydration and dynamical effects for the modeling of host–guest binding thermodynamics: the SAMPL5 blinded challenge Rajat Kumar Pal1,5 • Kamran Haider2 • Divya Kaur6 • William Flynn4,7 • Junchao Xia4 • Ronald M Levy4 • Tetiana Taran3 • Lauren Wickstrom3 • Tom Kurtzman2,5,6 • Emilio Gallicchio1,5,6 Received:17June2016/Accepted:25August2016/Publishedonline:30September2016 (cid:2)SpringerInternationalPublishingSwitzerland2016 Abstract As part of the SAMPL5 blinded experiment, we charged guests. Binding free energies obtained with computed the absolute binding free energies of 22 host– Debye–Hu¨ckel treatment of salt effects were in good guestcomplexesemployinganovelapproachbasedonthe agreement with experimental measurements. Water dis- BEDAM single-decoupling alchemical free energy proto- placement effects contributed favorably and very signifi- colwithparallelreplicaexchangeconformationalsampling cantly to the observed binding affinities; without it, the and the AGBNP2 implicit solvation model specifically modeling predictions would have grossly underestimated customized to treat the effect of water displacement as binding. The work validates the implicit/explicit solvation modeled by the Hydration Site Analysis method with approach employed here and it shows that comprehensive explicit solvation. Initial predictions were affected by the physical models can be effective at predicting binding lack of treatment of ionic charge screening, which is very affinities of molecular complexes requiring accurate treat- significant for these highly charged hosts, and resulted in ment of conformational dynamics and hydration. poor relative ranking of negatively versus positively Keywords SAMPL5 (cid:2) Hydration Site Analysis (HSA) (cid:2) Debye–Hu¨ckel (cid:2) Salt effects (cid:2) AGBNP2 (cid:2) BEDAM Theoriginalversionofthisarticlewasrevised:Correctionsdonein theoriginalarticlehasbeenpublishedintheerratum. Introduction & EmilioGallicchio [email protected] Accurate modeling of hydration is important for under- 1 DepartmentofChemistry,BrooklynCollege,2900Bedford standing the thermodynamics of binding and molecular Avenue,Brooklyn,NewYork11210,USA recognition[1–9].Bindingfreeenergymodelswithexplicit 2 DepartmentofChemistry,LehmanCollege,TheCity solvation are considered as the ‘‘gold standard’’ for mod- UniversityofNewYork,250BedfordParkBlvd.West, eling receptor-ligand complexation because of the critical Bronx,NewYork,NY10468,USA role of desolvation, the hydrophobic effect, and bridging 3 BoroughofManhattanCommunityCollege,Departmentof waters. In recent years, there has been significant progress Science,TheCityUniversityofNewYork,199Chambers in the development of computational tools to quantify the Street,NewYork,NY10007,USA structure and thermodynamics of water molecules on host 4 CenterforBiophysicsandComputationalBiology,Institute and protein surfaces [10–15]. Many of these approaches ofComputationalMolecularScienceandDepartmentof Chemistry,TempleUniversity,Philadelphia,PA,USA have highlighted the importance of displacing water 5 Ph.D.PrograminBiochemistry,TheGraduateCenterofthe molecules from enclosed regions of protein active sites in CityUniversityofNewYork,NewYork,NY10016,USA boosting the binding affinity [16–20]. 6 Ph.D.PrograminChemistry,TheGraduateCenteroftheCity Inthiswork,weemployHydrationSiteAnalysis(HSA) UniversityofNewYork,NewYork,NY10016,USA [10, 18, 20], a methodology based on Inhomogeneous 7 DepartmentofPhysicsandAstronomy,RutgersUniversity, SolvationTheory(IST)[21]tomapoutsolvationstructural Piscataway,NJ08854,USA and thermodynamic features in protein binding sites. High 123 30 JComputAidedMolDes(2017)31:29–44 density regions of water, termed hydration sites, are clas- AGBNP hydration site spheres to mimic the thermody- sified based on a number of energetic and structural mea- namicsofenclosedwatermolecules[24,32].Theresulting sures such as the water-protein interaction energy, water hybrid explicit/implicit model was then employed to pre- entropy and the strength of local water-water and water- dict the binding affinities of the SAMPL5 datasets. The protein hydrogen bonding. These measures are used to SAMPL5 blinded experiment, based on the octa-acid, [33] identifyregionsfromwhichthedisplacementofwaterleads tetra-endo-methyl octa-acid [34] and CB-clip host–guest to either significant gains or significant penalties to the systems [35], offers a unique opportunity to validate the predicted binding affinity. Ligands that displace high free approach outlinedhere.The tetra-endo-methylocta-acid is energy water molecules are morelikelyto bind strongly. referred to as methyl octa-acid throughout the article. We Moleculesareinherentlyflexible.Hence,computational have studied the octa-acid system before as part of the models ofmolecular recognition can not ignore dynamical previous SAMPL challenge [36], where we trained the aspects of binding. In this work, we employ the Binding AGBNP water sites empirically to reproduce available Energy Distribution Analysis Method (BEDAM) [22–24], experimental data. Here instead, we intend to train the a well established alchemical protocol for the quantitative model systematically using data from the HSA method, prediction of absolute binding free energies. BEDAM has grounded in well established hydration theories. Further- been specifically designed to capture conformational reor- more, the SAMPL5 datasets give us the opportunity to ganization and entropic contributions to the binding free validate our approach in an unbiased manner on hosts for energies. To overcome the slow time-scales of conforma- which experimental data is not available. tional reorganization, the method employs parallel multi- dimensional replica exchange conformational sampling algorithms[25,26]coupledwithanimplicitrepresentation Methods of the solvent [27–29]. The implicit solvent representation enables the formu- Hydration Site Analysis of the host binding cavity lationofbindingintermsofthedirecttransferoftheligand from the implicit solvent continuum to the receptor site We investigated the hydration properties of the three host [22]. This single-decoupling feature circumvents some of molecules using Hydration Site Analysis (HSA) approach the convergence difficulties encountered with double-de- [10,18].HSAutilizesanexplicitsolventMDsimulationof coupling approaches with explicit solvation [30, 31]. Key the solute molecule to identify spherical sites where water to the method is the sampling of the effective binding is present at high density. The thermodynamic quantities energy function u(x), which represents the ligand-receptor forthesesites,suchasenthalpy,first-orderentropyandfree interaction energy of a specific conformation x of the energy are calculated using inhomogeneous solvation the- complex averaged over the conformations of the solvent. ory (IST) [21]. A detailed description of the simulation The benefits of the single-decoupling approach is particu- setup used for HSA is given in the Computational Details larly significant for large and charged ligands, as in the section. In summary, the explicit solvent simulation for present application. The same approach is simply not fea- each of the host system in a restrained configuration is sible withexplicitsolvation since itwouldentailalengthy processed to obtain hydration site locations where water potential of mean force calculation for each evaluation of molecules are present with at least twice the bulk density. the effective binding energy function [30]. The clustering procedure to obtain hydration site locations The advantages of implicit solvation modeling in terms is described in detail elsewhere [20]. Briefly, for each of conformational sampling and the single decoupling hydration site, the total energy of the water, E , is cal- total pathway are, however, counterbalanced by the lack of a culatedasthesumofitsmeanwater-waterE ,andsolute- ww detailed representation of hydration effects. A way around water interaction energies E , where E is one half the sw ww these conflicting requirements is to train the implicit sol- mean interaction energy of the water in the site with all ventenergyfunctionbasedonexplicitsolvationdataaswe other waters, and E is one half the mean interaction sw set out to do in this work. BEDAM employs the Analytic energy of the water in the site with the solute. The factors Generalized Born plus Non-Polar (AGBNP) model which of one half follow the convention that half of the interac- has been specifically tailored for protein-drug binding tion energy in a pairwise interaction is assigned to each applications[29].Oneimportantaspectofthemodelisthe molecule of the pair. The locations and average solvation treatment of short ranged water-solute interactions by energies of the hydration sites obtained for each host means of hydration site spheres [24]. molecule are shown in Fig. 1 and Table 1. The solvation In this work, structural and thermodynamic data energies were compared to TIP3P neat liquid simulations obtained from explicit solvent Hydration Site Analysis are todeterminewhichhydrationsitescontainedfavorableand used to inform the placement and energetic strength of unfavorable water molecules relative to bulk solvent. The 123 JComputAidedMolDes(2017)31:29–44 31 Fig.1 Locationofhydrationsiteswithinthebindingcavityofthehostsasidentifiedfromhydrationsiteanalysis.aOcta-acid,bmethylocta- acid,cCB-clip Table1 SummaryofHydrationSiteAnalysis(HSA)resultsforthreehostsystemsinvestigatedinthisstudy Siteid Location Occupancy(p) Esw Eww Etotal ½Etotal(cid:3)Ebulk(cid:4)p Nnbrs fenc A. Octa-acid 0 Bottom 0.62 -4.57 -2.42 (cid:3)6:99 1.58 1.27 0.76 1 Top 0.40 -0.64 -7.59 (cid:3)8:23 0.52 3.06 0.42 2 Top 0.41 -0.61 -7.41 (cid:3)8:02 0.62 2.89 0.45 3 Top 0.40 -0.58 -7.62 (cid:3)8:19 0.53 3.12 0.41 4 Top 0.40 -0.58 -7.55 (cid:3)8:13 0.56 2.96 0.44 5 Middle 0.29 -1.48 -7.50 (cid:3)8:98 0.16 2.89 0.45 B. Methylocta-acid 0 Bottom 0.60 -4.52 -2.14 (cid:3)6:66 1.72 1.32 0.75 1 Top 0.36 -0.48 (cid:3)7:00 (cid:3)7:48 0.74 2.27 0.57 2 Top 0.31 -0.73 -6.84 (cid:3)7:57 0.61 2.54 0.52 3 Rim 0.31 -0.02 (cid:3)8:49 (cid:3)8:51 0.31 3.19 0.39 4 Top 0.32 -0.85 -6.83 (cid:3)7:68 0.60 2.50 0.52 5 Top 0.30 -0.62 -6.68 (cid:3)7:30 0.66 2.24 0.57 6 Rim 0.27 0.60 -9.67 (cid:3)9:06 0.66 3.64 0.31 C. CB-clip 0 Center 0.92 -5.76 -2.27 (cid:3)8:03 1.38 1.75 0.67 3 Naphthalenerings 0.42 -1.55 (cid:3)6:97 (cid:3)8:52 0.42 3.43 0.35 Allenergeticquantitiesareexpressedinkcal/mol.ForCB-clip,atotalofninehydrationsitewereobtainedinthecavity,thetwositeusedfor AGBNP2parametrizationarereportedhere.E =-9.53kcal/mol bulk relative solvation energies and positions of the explicit below.TheHSAscores(listedinTable 1)aregivenbythe hydration sites were used to position and parameterize the expression ½E ðiÞ(cid:3)E (cid:4)pðiÞ where i is the index of the total bulk strength of the AGBNP2 hydration site adjustments which HSAhydrationsites,p(i)isthewateroccupancyofthesite, account for enclosed hydration effects. We chose to only E ðiÞ is the total energy of the site and E is the cor- total bulk useHSAenergiestoparameterizethestrengthofhydrogen responding reference value obtained from TIP3P neat bond correction parameters because HSA entropies were water.Twowatermoleculesareconsideredtobefirst-shell previouslyshowntobenon-predictiveinascoringfunction neighbors in a given MD frame if their oxygen atoms are usedtopredictrelativebindingaffinitiesofligandsbinding separatedbylessthan3.5A˚ Thisdistancecriterionisbased to Factor Xa [20]. For each system, we derive scores for onthe location ofthe minima ofthe oxygen-oxygen radial hydration sites and apply them as corrections to AGBNP2 distributionfunctionforbulkwater[37].Themeannumber sites for the corresponding system (Fig. 2a) as described of these first shell neighbors for water in each hydration 123 32 JComputAidedMolDes(2017)31:29–44 Fig. 2 a Position of the hydration spheres added to AGBNP2 between two naphthalene rings; the center of masses of carbons parameters for each hosts, b the carbon atoms used to position the markedingreenwereusedtopositionwater-sitesinthemiddlecavity hydration sphere in all three hosts; the carbon atoms marked in of the octa-acids; carbons marked in orange were used to model a magentaareusedforpositioningthewatersitesonthetopcavityof watersiteatthebottomcavityoftheocta-acids;inCB-clip,awater octa-acids;inCB-clip,thosecarbonsareusedtopositionwatersitein siteispositionedatthecenterofmassofthosefourcarbons site is termed Nnbrs, and this first-shell neighbor count is DG ¼DG þDG þDG ð2Þ h elec np hs used to define the fractional enclosure of a hydration site, Theelectrostaticcomponentofthehydrationfreeenergyis byreferencingthenumberofitsfirst-shellneighborstothe represented by a pair-wise descreening implementation of mean number of first-shell neighbors of a TIP3P water molecule in bulk, Nnbrs: the Generalized Born model [28, 40, 41]. The non-polar bulk component of the hydration free energy is modeled as the Nnbrs sum of cavity and solute-solvent dispersions interactions f ¼1(cid:3) ð1Þ enc Nnbrs [28]. AGBNP2 includes algorithms to place and score bulk hydration spheres on the surface of the solute to describe This quantity indicates the degree to which the water in a short-range solute-solvent interactions such as hydrogen hydration site is blocked from contact with other water bonding [29]. The functional form of the short-ranged molecules. The value of Nnbrs is 5.25 as calculated from hydrogen bonding component of the solvation free energy bulk TIP3P pure water. in Eq. (2) above is X DG ¼ h Sðw Þ hs s s ð3Þ The AGBNP2 implicit solvation model s wherew isthewateroccupancyfactorofthesitemeasured s Here we employed the OPLS-AA/AGBNP2 effective as the fraction of the volume of the hydration site sphere potential,whereOPLS-AA[38,39]representsthecovalent actually accessible to water defined as and non-bonded inter-atomic interactions and solvation is Vfree modeled implicitly by the Analytic Generalized Born plus w ¼ s ð4Þ s V non-polar (AGBNP2) model [29]. The hydration free s energyDGh iscomputedasthesumofelectrostaticDGelec, whereVs isthevolumeofthewatersphereandVsfree isthe non-polar DGnp, and short-range solute-water hydrogen free volume of the water site, that is the volume not bonding interaction DGhs: occludedbysoluteatoms,obtainedbysummingallthetwo 123 JComputAidedMolDes(2017)31:29–44 33 body,threebody,etc.overlapvolumesofthewatersphere sodium phosphate solution used for the measurements of with the solute atoms i, j, k, etc. as [29] CB-clip guest binding, was used throughout. X X X Vfree ¼V (cid:3) V þ V (cid:3) V s s si sij sijk ð5Þ i i\j i\j\k Binding free energy protocol In Eq. (3), S is a switching function, and the h are s adjustable parameters which measure the strength of the Binding free energies were computed using the Binding hydrogen bonding energy (more precisely the portion of it EnergyDistributionAnalysisMethod (BEDAM)[22].The not captured by the Generalized Born model) [29]. The methodusesasingle-decouplingalchemicalapproachwith specification of the solute atoms carrying such hydration implicit solvation. An alchemical progress parameter k is spheresisaccomplishedbymeansofSMARTSpatterns.A introducedrangingfrom0,correspondingtotheuncoupled parameterdatabaseisusedtomapSMARTSpatternstothe state of the complex, to 1, corresponding to the coupled energeticparametersandplacementgeometryofhydration state of the complex. The k-dependent effective soft-core spheres [29]. The h parameters are normally negative and hybrid potential energy function is: s adjusted to describe favorable solute-solvent hydrogen U ðrÞ¼U ðrÞþkuðrÞ ð7Þ k 0 bonding interactions. In this context, the DG term disfa- hs vors binding by increasing the ligand and receptor desol- where r ¼ðr ;r Þ are the atomic coordinates of the A B vationpenalties.Inaddition,inthiswork,weusethesame receptor-ligand complex and r and r denote the coordi- A B functionalformtomodelenclosedwatermolecules,which nates of the receptor and ligand, respectively. U ðrÞ¼ 0 contribute favorably to binding, when displaced by the Uðr ÞþUðr Þisthepotentialenergyofthecomplexwhen A B ligand.Thekeydistinctionbetweenhydrogenbondingsites the receptor and ligand are uncoupled, i.e. as if they are andenclosedhydrationsitesisthesignoftheh parameter, separatedatinfinitedistancefromeachother.Thequantity s whichispositiveforthelatterandnegativefortheformer. u(r) is the binding energy and is defined as the change in The specific values of h assigned to enclosed hydration the effective potential energy of the complex for bringing s sites are derived from explicit solvent HSA analysis and the receptor and ligand from infinite separation to the are given in the Results section below. conformation r of the complex: uðrÞ¼Uðr ;r Þ(cid:3)ðUðr ÞþUðr ÞÞ ð8Þ A B A B Incorporation of Debye–Hu¨ckel parameters In this work we adopted a soft-core form of the binding in the implicit solvation model energy mentioned above and as described [23, 43]. UWHAM analysis of the binding energy samples was In this work the Generalized Born components of the carried as described [43] to obtain binding free-energy AGBNP2 model included modifications based on the estimates and their corresponding statistical uncertainties. Debye–Hu¨ckel screening aimed at modeling the effects of Average interaction energies, DE were obtained by aver- b electrolytes present in the aqueous solution. Specifically, aging the binding energy values collected at k=1 (the we adopted the following functional form for the GB pair bound state). The uncertainty on interaction energies were potential [42] measured as the standard error of the mean. The entropic u ði;jÞ¼(cid:3)(cid:2) 1 (cid:3)e(cid:3)jrij(cid:3) qiqj contributiontothefreeenergy,theeffectofreorganization GB (cid:2) (cid:2) rffihffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiiffiffi energy and reorganization free energy to binding were out in ri2jþBiBje(cid:3)ðri2j=4BiBjÞ derived [26, 44] from binding energies and other thermo- dynamic parameters as described in the computational ð6Þ details section. The uncertainties of all these quantities where r is the distance between atoms i and j, (cid:2) ¼80 were computed by error propagation. ij out and (cid:2) ¼1 are the dielectric constants of the solvent and In this work, BEDAM conformational sampling is in the interior of the solute respectively; q and q are the acceleratedbymeansoftwo-dimensionalreplicaexchange i j partialchargesofatomsiandj.ThetermsB andB arethe along alchemical and temperature parameters [26, 45]. k i j Born radii of atoms i and j. The Debye–Hu¨ckel screening exchanges facilitate the mixing of intermolecular degrees p parameter is defined as j¼ ffi8ffiffipffiffiffikffiffiffiffiIffiffi where I is the ionic of freedom, while temperature exchanges activate B strength of the solution and k is the Bjerrum length intramolecular degrees of freedom thus allowing confor- B [40, 42]. For water at room temperature, k =7.0 A˚. An mational transitions to explore the conformational space B ionic strength of 0.12 M, corresponding to the 20 mM more efficiently [46]. 123 34 JComputAidedMolDes(2017)31:29–44 Computational details the water molecules in those locations. For each host, 190, 000 frames of the trajectory were processed with Explicit solvent calculations and Hydration Site Analysis HSA analysis. Hydration site locations were identified (HSA) using a subset of 10,000 frames which were evenly spaced in time. High density spherical regions (hydration For each host molecule, a molecular dynamics simulation sites) of 1A˚ radius were identified using a clustering was performed using the OPLS 2005 force field (Schro- procedure [10] on the water molecules that were found dinger Inc.) in a TIP3P [47] cubic water box, using the within 5 A˚ of the host molecule. The resulting hydration DESMOND package [48]. The initial startingstructure for sites were each populated by retrieving all water mole- eachsimulationwasarepresentativeholostructurewithout cules, which had oxygen atoms within 1.0 A˚ from the the ligand. These representative structures were obtained corresponding hydration site center. The hydration sites from BEDAM simulations of the octa-acid and methyl were then enumerated according to their occupancies, octa-acid hosts with guest 4 and the CB-Clip host with with the highest populated site given the index 0. The guest1.Eachwaterboxwasbuiltwithaminimumdistance average solvation energies were calculated based on of 10 A˚ between the solute and the edge of the box. The energies of the water molecules that were within 1.0 A˚ default Desmond equilibration protocol was used to mini- from the coordinates of each hydration site center in the mize and thermalize the systems for the production run. original 190,000 frames. The system was further relaxed for 25 ns under NPT conditions at 1 atm and 300 K. An NVT production run was started from a configuration in which the box volume Systempreparationforthebindingfreeenergycalculations was close to the average of the previous NPT simulation. The production MD simulations were run for 100 ns and Thestructures ofthe octa-acid and CB-clip hostsandtheir snapshotswerecollectedevery0.5ps.Thefirst1nsofthe guests were prepared using Maestro (Schrodinger Inc.) production run was discarded. In the octa-acid and methyl starting from structure files provided by the SAMPL5 octa-acidsimulations,positionalrestraintswereusedonall organizers. Bond orders and formal charges were adjusted oftheheavyatomswithaforceconstantof5.0kcal/mol/A˚2 asappropriate.Alltheeightcarboxylategroupsoftheocta- which allowed for free rotation of hydrogens. In the CB- acids were modeled as deprotonated resulting in a net clip simulations, position restraints were used on all of the negativechargeof(cid:3)8.Similarly,thefoursulfonategroups heavy atoms of the host with a force constant of 5.0 CB-cliphostweremodeledasallionized,resultinginanet kcal/mol/A˚2, excluding the sulfonate groups. The use of overall charge of (cid:3)4. Alternate protonation and tau- one rigid receptor should not have a significant effect on tomerization states of the CB-clip guests were generated determining the hydration sites and the following free using the LigPrep facility (Schrodinger Inc.) using Epik energycalculationsduetothenatureofthescalingfunction [49] with a pH range of 7±2. Ionization and tautomeriza- used to calculate enclosure energies and the rigidity of the tion free energy penalties were recorded and added to the SAMPL5 hosts. If the host structures are unrestrained BEDAM binding free energy estimates to compute the duringHSAanalysis,theoccupanciesofthehydrationsites predictedbindingfreeenergies.LigPreppredictedmultiple are expected to decrease due to the effect of the confor- protonation and tautomerization states for the CB-clip mationalfluctuationsofthehostonthesolvent.Incontrast, guests. Multiple protonation states of guests 1, 4, 5, 6, 8 the energies of these hydration sites may become more and 9 were tested individually. Guests 3 and 5 which favorable due to the ability of the solvent to optimize its contain alkylammonium groups were modeled as posi- interactionswitheithertheneighboringwatermoleculesor tively charged and the guests containing carboxyl groups, nearby solute atoms. These changes to the HSA energies guests 1, 2, and 4 were modeled as deprotonated each and occupancies will have an impact on the global having a single negative charge. Guest 6, which is parameter c and should not impact the enclosure energies nitrobenzoic acid, holds an overall single negative charge usedinthefreeenergycalculations.Ontheotherhand,this contributed by the carboxyl moiety. To match the experi- treatmentisinappropriateformoreflexiblehoststhatadopt mental pH, the acidic guests of the octa-acid set were multiple bound states because the hydration site locations modeledasdeprotonatedandionizationpenalties werenot will vary for each bound representative structure. SHAKE applied. wasusedtoconstrainanybondsinvolvinghydrogenatoms. All individual guests were manually docked to the Hydration Site Analysis involved two steps: (1) identi- bindingcavityofthehostsintwooftheocta-acidsandthe fying regions of high water density around each host sur- in-betweenspaceofthetwonaphthaleneringsonCB-clip. face and (2) calculating the average solvation energies of Joint Hamiltonian and temperature 2-D Asynchronous 123 JComputAidedMolDes(2017)31:29–44 35 Replica-exchange Molecular dynamics simulations [26] Results employed18intermediate kstepsasfollows: k=0,0.002, 0.004, 0.008, 0.01, 0.02, 0.04, 0.07, 0.1, 0.17, 0.25, 0.35, The locations of hydration sites obtained for each host 0.5, 0.6, 0.7, 0.8, 0.9, and 1.0 and 8 temperatures were moleculeareshowninFig. 1.Fortheocta-acidhost,atotal distributed between 300 K and 534 K (T = 300, 328, 358 of six hydration sites are identified, one in the bottom ,390,424,458,494,534).Atotalof144replicaswereused cavity,oneinthemiddlecavityandfourintheuppercavity corresponding to all possible pairing of temperature and k (Fig. 1,left).Forthemethylocta-acidhost,atotalofseven values. The replica-exchange (AsyncRE) simulations were hydrationsitesareidentified,oneinthebottomcavity,four started from energy-minimized and thermalized structures in the upper cavity and two at the level of the rim of the fromthemanuallydockedmodels.Aflat-bottomharmonic host (Fig. 1, middle). The introduction of the methyl sub- restraint with a tolerance of 5 A˚ between the centers of stituentsmainlycausespartialorderingofwatermolecules massofthehostandtheguestwasapplied.Eachcycleofa near the rim of the octa-ocid host, which are otherwise replica lastedfor100picosecondswith 1fstime-step. The bulk-like. In the particular conformation examined, two average sampling time for a replica was approximately 10 partially ordered sites are found at the center of the mouth ns. Calculations were performed on the campus computa- of the cavity(Fig. 1, middle). For the CB-clip host, a total tional grids at Brooklyn College and at Temple University of nine hydration sites were identified, one of which is and on the XSEDE SuperMIC cluster, the latter using 3 located at the center of the host and the remaining eight compute nodes utilizing all CPU’s and both MIC devices. were symmetrically distributed, spanning the region from The binding energies obtained from all replicas were ana- naphthalene rings to the outer edges of the cavity. In the lyzed using UWHAM [43] method and the R-statistical bindingfreeenergycalculations(seebelow)weconsidered package to compute the binding free energy DG(cid:5) and thecentralsiteplusthesitesandwichedinbetweenthetwo b average interaction energies DE of individual complexes. naphthalene rings (Fig. 1, right). b Energy values from the first 1 ns of each replica were In this work, as in the previous SAMPL4 host–guest excluded from the free energy calculation. Calculation at edition [36] , we used custom AGBNP2 hydration spheres multiple temperature allows estimation of the conforma- to model the effect of enclosed water molecules displaced tional entropy of binding DS(cid:5) from the temperature upon ligand binding (see Methods). The idea is that these conf derivative of the binding free energy (here we used a sites have large water occupancy ðws ’1Þ in the unbound simple finite difference estimator using adjacent binding state and small water occupancy ðws ’0Þ in the presence free energy values). Reorganization free energies DG(cid:5) oftheboundligand.Hencepositivevaluesofthewatersite reorg were measured as the difference between the binding free energy parameter ðhsÞ can mimic the favorable effect of energy and the average interaction energy at the bound water expulsion towards binding. The approach we fol- state as lowed here for SAMPL5 is to place and score enclosed hydration sites which are distinct from the empirical DG(cid:5)reorg ¼DG(cid:5)b(cid:3)DEb ð9Þ approach adopted for the octa-acid host in SAMPL4. In SAMPL4 hydration site parameterization was conducted Finally, the reorganization energy for each complex was essentially by trial and error in an attempt to reproduce measured by subtracting the entropic component from the known experimental data. In this work instead the place- reorganization free energy: ment and scoring of AGBNP2 hydration sites is based on HSA analysis of explicit solvent data. Furthermore, DE(cid:5) ¼DG(cid:5) þTDS(cid:5) ð10Þ reorg reorg conf because of limitations of the AGBNP2 hydration site Statistical uncertainties for the reorganization free energy, placement algorithm, in SAMPL4 some of the effects of conformational entropy and reorganization energy were water expulsion were captured by empirical modifications obtained from error propagation of the uncertainties of the of surface tension and solute-solvent van der Waals binding free energy and the average binding energy. The parameters.Forthiswork,wedevelopednovelgeometrical data submitted to SAMPL5 challenge were obtained from algorithms to place hydration sites based on geometrical the 2D replica exchange binding free energy calculations centersofarbitrarygroupsofatomsandmoresignificantly, without inclusion of salt effects. Uncertainties were thanks to HSA analysis, we were able to reduce substan- reported as twice the standard error. The binding free tially, the number of empirically-derived parameters as energy predictions were submitted to the SAMPL5 host– described next. guest challenge on January 29, 2016. For our submissions, The locations of enclosed hydration sites identified by the octa-acid, methyl octa-acid and CB-clip were assigned HSA analysis were reproduced as best as possible using prediction IDs # 18, 20 and 10 respectively. AGBNP2 water spheres and the anchoring methods 123 36 JComputAidedMolDes(2017)31:29–44 available in AGBNP2 [29] and those developed for this octa-acid set [36]. Deviations from the experiments sug- work. Solute anchoring points correspond to hydrogen gested that the HSA analysis over-emphasized the benefit bondingsitessuchasalongpolarhydrogenbondsandlone of displacing water from the bottom site. To correct this pairs in various geometries [29], in addition to anchoring deficiency we decided to assign h parameters by dis- s points based on geometrical centers of groups of atoms. tributing the total enclosure energy equally over the six Waterspherelocationsdefinedinthiswaynaturallyfollow sites. This led to results in reasonable agreement with solutemovements.Thefourtopsitesoftheocta-acidhost, experimental affinities with a scaling factor of c¼2:88 sites (1, 2, 3, and 4) in Table 1, were modeled with eight withoutionic screeningandc¼1:85withionicscreening. AGBNP2 hydration spheres placed perpendicularly to The predictions for the octa-acid and methyl octa-acid inner phenyl rings of the cavity (Fig. 2b, left). The middle hosts were produced using this scheme. For the CB-clip hydration site of the octa-acid host (site 5, Table 1A) was host, lacking experimental evidence supporting one strat- modeled by an AGBNP2 site placed at the geometrical egy over another, we employed the more physically rea- center of a set of aromatic carbons as shown in Fig. 2b, sonableapproachofderivingh parametersfromindividual s left.Thebottomwatersite(site0inTable 1A)wasplaced enclosure energies as described above. The same scaling similarly.Thesix hydration sitesidentifiedforthe methyl- factors derived for the octa-acid system were used for the octa-acid host (Table 1B), were modeled using three CB-clip system. AGBNP2hydrationspheresasindicatedinFig 2a,middle. The blinded binding free energy predictions DG(cid:5) calc The one at the rim ofthe cavity corresponds to sites 3 and submitted to the SAMPL5 challenge, obtained without the 6,thesecondinthetopregionofthecavitycorrespondsto inclusionofsalteffects,areshowninTables 2, 3and 4for sites 1, 2, 4 and 5, and finally the bottom site. All these the octa-acid, methyl octa-acid and CB-clip sets, respec- individual water sites were placed at the center of mass of tively and Fig. 3, together with the corresponding experi- specific set of carbons surroundingthe three regionsin the mentalvalues.Thesepredictionsdidnotreproducewellthe pocket (Fig. 2b, middle). For the CB-clip host, one water experimental measurements, with the most noticeable sitewasplacedatthecenterofmassof4aromaticcarbons shortcoming being the significant overestimation of the and4 alkyl carbonsinthe center ofthe cavity as shownin affinitiesofpositivelychargedguests(guests3and5ofthe (Fig. 2b, right). The second water site was placed in octa-acid guests, and guests 6 and 10 of the CB-clip host). betweenthetwonaphthaleneringsatthecenterofmassof Theresultswiththeinclusionofionicchargescreeningare 4 aromatic carbon atoms indicated (Fig. 2b, right). We inmuchbetteragreementwiththeexperiments.Rootmean have not attempted to model hydration sites located above squareerrors(RMSE)relativetoexperiments(Tables 2, 3; and below the CB-clip host cavity partly because of their Fig. 3a, b) are reduced by more than a factor of 2 for the polar character already captured by AGBNP2 hydrogen- octa-acidhostandbymorethanafactorof3forthemethyl bonding sites [29]. octa-acid host. The RMSE for the CB-clip set (Table 4; The solvation parameters h associated with each Fig. 3c) is reduced from 4.7 to 3.6 kcal/mol. Introduction s AGBNP2 water site of the octa-acid host were derived of salt effects also improved the correlation between from HSA analysis using the water enclosure energies as computational predictions and the experiments. Pearson correlationcoefficients (r)improvedfrom0.66to0.89and E ðiÞ¼c½E ðiÞ(cid:3)E (cid:4)pðiÞ ð11Þ enclosure total bulk from0.67to0.90forthemethylocta-acidandCB-clipsets, wherep(i)isthewateroccupancyofthesiteand½E ðiÞ(cid:3) respectively. total E (cid:4)istheaverageenergyofthesiterelativetotheenergy The correlation coefficient for the octa-acid set bulk of bulk water (Table 1). Enclosure energies were rescaled improvedsignificantlywiththeinclusionofsalteffectsbut by a global parameter, c, to obtain AGBNP2 hydration remained low. This can be considered an acceptable result parameters for enclosed water spheres. The global scaling given that, with the exception of guest 4, the range of parameterwasadjustedsoastoreproducetheexperimental experimental affinities of the octa-acid set is small (less binding free energies of the SAMPL4 octa-acid set as than1.8kcal/mol).Thestrongestbinderinthissetisguest described below. Scaled enclosure energies of hydration 4. The calculations do not reproduce this feature well. sites corresponding to a single AGBNP2 hydration sphere Instead,guest3,apositivelychargedligand,ispredictedto were aggregated into the h parameter of that sphere. bindthestrongestwith-8.8kcal/molwithionicscreening. s Conversely, the h parameters for AGBNP2 sites corre- The trend of overprediction of the affinity of positively s spondingtomorethanonehydrationsiteswereobtainedby charged ligands, while significantly reduced, remain even distributing enclosure energies equally. after inclusion of ionic screening. The origin of the rela- Using the scheme above, we were unable to reproduce tivelylargedeviationbetweenexperimentalandcalculated accurately the experimental affinities for the SAMPL4 affinities for guest 4 is not clearly evident. This complex 123 JComputAidedMolDes(2017)31:29–44 37 Table2 Calculatedandexperimentalbindingfreeenergieswithandwithoutsalteffectsforocta-acidsetandthermodynamicdecompositionof bindingfreeenergieswiththeinclusionofsalteffects Guests Structure DG(cid:5)a Withoutsalteffects Withsalteffects exp DG(cid:5)b DG(cid:5)c DEd DG(cid:5)e DEf (cid:3)TDS(cid:5)g calc calcDH b reorg reorg conf G1 O (cid:3)5:4 (cid:3)2:2(cid:6)0:06 (cid:3)3:8(cid:6)0:08 (cid:3)15:5(cid:6)0:01 11.7±0.09 1.0±1.4 10.7±1.4 O– G2 O (cid:3)4:7 (cid:3)6:3(cid:6)0:06 (cid:3)6:6(cid:6)0:09 (cid:3)15:9(cid:6)0:01 9.2±0.09 0.9±1.6 8.4±1.5 O– N G3 (cid:3)4:5 (cid:3)15:4(cid:6)0:06 (cid:3)8:8(cid:6)0:09 (cid:3)20:4(cid:6)0:01 11.6±0.09 -0.9±1.5 12.5±1.5 N+ G4 Br (cid:3)9:4 (cid:3)6:3(cid:6)0:08 (cid:3)6:8(cid:6)0:2 (cid:3)20:0(cid:6)0:02 13.2±0.1 1.3±2.6 11.9±2.6 O –O G5 (cid:3)3:7 (cid:3)11:3(cid:6)0:06 (cid:3)6:3(cid:6)0:1 (cid:3)16:6(cid:6)0:02 10.3±0.1 1.7±1.7 8.6±1.7 N+ G6 O O (cid:3)5:3 (cid:3)7:4(cid:6)0:06 (cid:3)7:0(cid:6)0:09 (cid:3)18:8(cid:6)0:02 11.8±0.1 0.2±1.6 11.6±1.6 N+ –O O– RMSE 5.8 2.6 Correlationcoefficient(r) (cid:3)0:39 (cid:3)0:04 Allvaluesinkcal/mol.aExperimentalITCbindingfreeenergy[35].bPredictedbindingfreeenergywithoutadditionofioniceffects.cBinding freeenergypredictionwiththeadditionofsalteffects.dAverageinteractionenergyoftheboundcomplexat=1.Theuncertaintyisestimatedas thestandard errorofthemean. eReorganization free energyascalculatedusingEq.9, fReorganization energyorintra-molecular strainas calculatedusingEq.10, gBindingentropycomputedbyfinitedifferencemethodandevaluatedat300K;theuncertaintyiscomputedbyerror propagationforallthedecompositions has the largest reorganization penalty in the set, which is negatively-charged counterpart, guest 7, is correctly pre- surprising, given the rigidity of the guest. We ascribe this dicted to be a relatively weak binder. Guest 2 is correctly effect to the side-ways orientation of the bromine atom predictedastheweakestbinderoftheset.Thebindingfree relative to the carboxylate substituent. This causes con- energies of the guests in the middle of the pack are formationalstraininthe hostandtiltingofthe carboxylate reproduced reasonably well and within the expected to accomodate the bromine atom. We speculate that our accuracy limits of the model. The inclusion of ionic model may be over-estimating the steric hindrance expe- screening has the general effect of weakening binding and rienced by the bromine substituent. shifting the predictions closer to the experimental values. Binding trends in the methyl octa-acid are better repro- For some guests the shift is very substantial (as much as 9 duced, partlythankstolargervariationsinaffinities inthis kcal/molforguest3),whiletheaffinitiesofotherguests(5, set. Focusing on the more accurate results with ionic 8and9forexample)arebarelyaffectedbyionicscreening. screening, guest 4, the strongest binding for the octa-acid These differences appear to be related to the number of host, is here the weakest binder both experimentally and charged groups and the degree of solvent screening affor- computationally.Thecalculationsreproducethelargelossof ded by alkyl substituents. affinityofguest4uponmethylationofthehostreasonably In Tables 2, 3 and 4, we also report the results of free well. Furthermore, in agreement with the experiments, the energy decomposition analyses [44] (values without ionic calculationspredictthatthestrongestbinderisguest3. screeningarenotshown).AveragebindingenergiesðDE Þ, b The two strongest binders of the CB-clip host, guests 6 obtainedfrom theaverage ofthe binding energyvalues,u, and 10, are correctly ranked the best by the computational in the bound state (k¼1), reflect the strength of ligand- model, although their relative rankings are reversed. The host interactions, including desolvation. Subtracting these complexwithguest6ispredictedtohaveanunreasonably from the binding free energies yields reorganization bind- favorable binding free energy (-16.8 kcal/mol). Its ing free energies (DG(cid:5) ). Using data from higher reorg 123 38 JComputAidedMolDes(2017)31:29–44 Table 3 Experimental and computed binding free energies with and without salt effects for methyl octa-acid hosts and thermodynamic decompositionsofbindingfreeenergiescalculatedwithinclusionofsalteffects Guests Structure DG(cid:5)a Withoutsalteffects Withsalteffects exp DG(cid:5)b DG(cid:5)c DEd DG(cid:5)e DEf (cid:3)TDS(cid:5)g calc calcDH b reorg reorg conf G1 O (cid:3)5:3 (cid:3)7:5(cid:6)0:06 -4.8±0.08 -17.2±0.01 12.4±0.08 0.8±1.4 11.5±1.4 O– G2 O (cid:3)5:1 (cid:3)10:6(cid:6)0:07 -7.7±0.09 -17.6±0.01 9.9±0.1 1.4±1.7 8.5±1.6 O– N G3 (cid:3)6 (cid:3)16:2(cid:6)0:07 -9.4±0.1 -21.4±0.01 12±0.1 0.7±1.7 11.3±1.7 N+ G4 Br (cid:3)2:4 (cid:3)2:4(cid:6)0:1 -0.7±0.2 -21.1±0.02 20.4±0.2 6.7±2.7 13.7±2.7 O –O G5 (cid:3)3:9 (cid:3)14:5(cid:6)0:07 -5.7±0.1 -16.5±0.02 10.8±0.1 1.2±1.8 9.6±1.8 N+ G6 O O (cid:3)4:5 (cid:3)8:9(cid:6)0:06 -5.8±0.09 -18.0±0.02 12.2±0.1 1.4±1.6 10.8±1.6 N+ –O O– RMSE 6.7 2.1 Correlationcoefficient(r) 0.66 0.89 Allvaluesinkcal/mol.aExperimentalITCbindingfreeenergy,exceptforG4andG5forwhichonlyNMRmeasurementsareavailable[35]. bPredictedbindingfreeenergywithoutadditionofioniceffects. cBindingfreeenergypredictionwiththeadditionofsalteffects. dAverage interactionenergyoftheboundcomplexat=1.Theuncertaintyisestimatedasthestandarderrorofthemean. eReorganizationfreeenergyas calculated using Eq.9, fReorganization energy or intra-molecular strain as calculated using Eq.10, gBinding entropy computed by finite differencemethodandevaluatedat300K;theuncertaintyiscomputedbyerrorpropagationforallthedecompositions simulated temperatures, the latter is further decomposed predictors of affinity. For example, guest 3, which is the intoconformationalentropy((cid:3)TDS(cid:5) )andreorganization second weakest binder in the set, is predicted to have the conf energy (DE ) components. The conformational entropy strongest interactions with the host (-47.2 kcal/mol). reorg component measures the entropy loss for the formation of Conversely, guest 5, which is one of the top binders, has the complex. The reorganization energy measures the relatively weak interactions with the host (-26.2 kcal/- intramolecular energy change of the guest and the host as mol). Clearly, as observed previously [32, 36], reorgani- theychangeconformationstobindtoeachother.Thelatter zation(entropiclossandintramolecularenergystrain)isan is commonly referred as intramolecular energy strain important determinant of binding. because it often (but not always) opposes binding. As further discussed below, decompositions of binding free energiesprovideusefulphysicalinterpretationsofobserved Discussion bindingaffinitytrends.Asoftenobserved[32],thebinding free energy is the result of a large compensation between Water molecules solvating hydrophobic enclosures exhibit theaveragebindingenergyandreorganizationfreeenergy. unique thermodynamic and structural signatures. Because Forexample,thepooraffinityofguest4tothemethylocta- they are in concave regions, they have fewer proximal acid host can be ascribed primarily to a large and unfa- water neighbors with which they can form favorable vorable reorganization energy (20.4 kcal/mol) rather than hydrogen bonding interactions. Simultaneously, due to the the strength of the host–guest interaction energy (-21.1 hydrophobicityofthesolute,theyareunabletocompensate kcal/mol), which is among the most favorable in the set. for these lost interactions by forming favorable hydrogen The average binding energies for the CB-clip complexes bond contacts with the solute surface. As a result, these have a particular wide range, from -20.2 kcal/mol to watermoleculesareenergeticallyunfavorablewithrespect -47.2 kcal/mol. Interestingly, these are imperfect to bulk and their displacement from the surface into the 123