Vol.:(0123456789)1 3 Journal of Computer-Aided Molecular Design (2022) 36:1–9 https://doi.org/10.1007/s10822-021-00439-w Binding free energies for the SAMPL8 CB8 “Drugs of Abuse” challenge from umbrella sampling combined with Hamiltonian replica exchange Daniel Markthaler1 · Hamzeh Kraus1 · Niels Hansen1 Received: 8 September 2021 / Accepted: 11 December 2021 / Published online: 3 January 2022 © The Author(s) 2021 Abstract Umbrella sampling along a one-dimensional order parameter in combination with Hamiltonian replica exchange was employed to calculate the binding free energy of five guest molecules with known affinity to cucurbit[8]uril. A simple empirical approach correcting for the overestimation of the affinity by the GAFF force field was proposed and subsequently applied to the seven guest molecules of the “Drugs of Abuse” SAMPL8 challenge. Compared to the uncorrected binding free energies, the systematic error decreased but quantitative agreement with experiment was only reached for a few compounds. From a retrospective analysis a weak point of the correction term was identified. Keywords SAMPL8 · Umbrella sampling · Host–guest complex Introduction Among the many possible methodological variants to com- pute binding free energies for host-guest systems within a rigorous thermodynamic framework [1], umbrella sampling (US) [2] has emerged as an easy-to-use approach that is sup- ported by many publicly available pre- and post-analysis tools [3]. If certain artefacts related to insufficient sampling are accounted for [4], the method yields robust estimates. Based on our previous experience with this approach in the context of cyclodextrin host-guest systems [4], we aimed for the evaluation of the simulation protocol in the SAMPL8 challenge, featuring more complex guest molecules than considered before. We hope that the comparison to other free-energy methods applied to the same host-guest systems within the SAMPL challenge is of use for the continuous evaluation of methods and force fields pushing forward the field of binding free energy calculations, which, for complex systems, often faces a combination of force-field and sam- pling issues [5, 6]. Computational approach Molecular model Coordinates of host and guest molecules (see Fig. 1) for the SAMPL8 CB8 challenge were obtained from the SAMPL8 GitHub repository [7] in form of mol2 files. In order to incorporate pH effects in the simulation protocol, for every ligand a deprotonated and corresponding protonated form was parametrized. The deprotonated species were created by manually removing the respective hydrogen atoms from the structure files. To generate the topology files, general Amber force field (GAFF) [8] atom types were assigned using the antechamber module of Ambertools20 [9] with the molecule structures as an input. For the protonated types, the partial charges in the provided mol2 file were retained while in the deprotonated cases these were obtained from the semiem- pirical AM1-BCC model [10, 11]. Due to slight rounding errors in the calculated values, the resulting excess charge had to be distributed among the different atoms in a weighted manner, which was done using the Python community package ParmEd [12]. These molecule files were then used to generate topology files in Amber format using the leap module of Ambertools before converting them into the GROMACS format using the Par- mEd package. The same pipeline was used to parametrize the CB7 host and a further set of five guest molecules used for training purposes (see Fig. 2). * Niels Hansen hansen@itt.uni-stuttgart.de 1 Institute of Thermodynamics and Thermal Process Engineering, University of Stuttgart, 70569 Stuttgart, Germany http://orcid.org/0000-0003-4366-6120 http://crossmark.crossref.org/dialog/?doi=10.1007/s10822-021-00439-w&domain=pdf 2 Journal of Computer-Aided Molecular Design (2022) 36:1–9 1 3 Fig. 1 Structures of the CB8 host and the seven guest mol- ecules forming the SAMPL8 “Drugs of Abuse” challenge G1 (Methamphetamine) G2 (Fentanyl) G3 (Morphine) G4 (Hydromorphone) G5 (Ketamine) G6 (Phencyclidine) G7 (Cocaine) Cucurbit[8]uril (CB8) GT1 (3,5-Dimethyl-1,1-adamantanamine) GT2 (Cyclododecanamine) GT3 (3-Amino-1-adamantanol) GT4 (Cycloheptanamine) GT5 (Cyclooctanamine) Fig. 2 Structures of the training set molecules GT1 to GT5 3Journal of Computer-Aided Molecular Design (2022) 36:1–9 1 3 The mol2 file of CB7 was obtained from the Amber tutorials website [13]. The additional five guest molecules were 3,5-dimethyl-1,1-adamantanamine (GT1), cyclodode- canamine (GT2), 3-amino-1-adamantan-1-ol (GT3), cyclo- heptanamine (GT4) and cyclooctanamine (GT5). Except for GT1, these molecules were part of the SAMPL6 chal- lenge and the mol2 files were obtained from the correspond- ing GitHub repository. For GT1, initial coordinates were obtained by DFT optimization [14–16] and the molecular model was parametrized as described above. Water was rep- resented by the TIP3P model [17] while the chloride ion used to neutralize the simulation box was described with the parameters of Joung and Cheatham [18]. Alternative partial charges for the guest molecules were calculated with the DDEC6 approach [19]. The effect on binding free energies was marginal for GT1 and more sig- nificant for G1 (compare Fig. S1), but far below the dis- crepancy between simulation results and experiment. There- fore, all further results correspond to the AM1-BCC partial charges. System preparation For the given concentration of sodium phosphate buffer of 20 mM [20] we estimated the concentration of H2PO − 4 and HPO2− 4 at pH 7.4 by means of the Henderson-Hasselbalch equation [21] using an available online tool [22], resulting in concentrations of 5.40 mM and 14.60 mM, respectively. Although AMBER/GAFF-compatible parameters have been reported for these ions [23], the low concentration (approxi- mately 1 or 2.5 molecules, respectively, per 10000 water molecules) suggests that not capturing the ionic strength is a rather secondary effect for this dataset. Moreover, low ion concentrations may cause sampling issues [24, 25]. There- fore, we decided to include only one Cl− counterion for the protonated guest molecules and lumped a possible error caused by not accounting for the ionic strength into a cor- rection anyway applied to the raw MD data (see below). Available p Ka-values for the seven ligands suggest that the protonated form of the molecule is the dominant one in a solution of pH 7.4, except for guest 5, which has a p Ka-value of 7.5, and therefore the concentration of the protonated form is only marginally larger than that of the unprotonated form. For all guest molecules both the protonated and depro- tonated form were considered in the binding free energy calculations. Reported is the value of the form with higher binding affinity, which was always the protonated form. Molecular dynamics simulations All simulations were conducted with the GROMACS 2016.4 program [26] patched to the free-energy library PLUMED 2.4.2 [27] for restraints definition. Free-energy profiles were constructed from the time series of a single order parameter sampled via umbrella sampling [2]. The order parameter is given by the projection of the instantaneous separation vector between the centers of mass (COM) of the binding partners onto the host’s instantaneous symmetry axis and was sampled between −2 and +2 nm in steps of 0.1 nm and using a force constant of 4000 kJmol−1nm−2 for keeping the ligand in the respective umbrella window. Lateral movement of the ligand at every umbrella window was restricted with the aid of a flat-bottom potential acting on the orthogonal displacement of the ligand’s COM from the host’s molecular axis. The flat-bottom width of 0.5 nm was chosen such that no force was acting on the bound ligand. The host’s orien- tation was aligned along the z-axis of the simulation box alongside with a translational restraint ( 500 kJmol−1nm−2 ) to keep its COM close to the box center. The order parameter was written to file every 100 steps. To facilitate configu- rational sampling, a Hamiltonian replica exchange scheme [28, 29] was used that attempts an exchange move between neighboring umbrella windows every 1000 steps. Each sampling window included a minimum of 40 ns simulation time. Each system was solvated with ∼ 2000 TIP3P waters in an orthorhombic box whose dimensions were approximately 36 × 36 × 52Å3 . Starting configurations for each umbrella window were generated by displacing the ligand out of the host’s binding pocket through application of the set of restraining potentials described above and simulat- ing for 100 ps in order to relax the system. The bound state of the ligand was determined in a preceding step by allow- ing the ligand to be captured into the CB cavity within a 2 ns NPT simulation via a harmonic distance restraint (force constant: 2000 kJmol−1nm−2 ) applied for the COM-COM distance. The stability of the complex configuration found this way was then further examined within another 2 ns NPT simulation without any restraining potential between ligand and host. This procedure was conducted 3 times within 3 independent simulations and the final snapshots of the bound state configurations were visually compared. In none of the studied cases involving CB8, significant differences could be determined in the 3 respective complex configura- tions and all bound states were stable within 2 ns unbiased simulations. Production simulations were conducted in the NPT ensemble at 300 K and ambient pressure, with temperature control using a Langevin thermostat [30] with an inverse friction constant of �T = 1.0 ps and Parrinello-Rahman barostat [31, 32] with coupling constant �p = 2.0 ps . A Ver- let-buffered neighbor list [33] which was updated every 10 steps, was applied for the treatment of short-range electro- static and van der Waals interactions with potentials shifted to zero at 0.9 nm. The latter were modeled by the Lennard- Jones potential. Analytic dispersion corrections were applied for energy and pressure calculation. Long-range electrostatic 4 Journal of Computer-Aided Molecular Design (2022) 36:1–9 1 3 interactions were treated with the smooth particle-mesh Ewald (PME) method [34, 35] using a real-space cut-off of 0.9 nm with a cubic splines interpolation scheme and a grid spacing of 0.1 nm. The center of mass translation of the computational box was removed every 100 steps. All bond lengths involving hydrogens were constrained using the SHAKE algorithm [36] with a relative tolerance of 0.0001. The hydration free energy and octanol water partition coefficient were calculated for the neutral form of G1 from alchemical decoupling of the molecule from the surround- ing bulk phase consisting of cubic boxes with one guest molecule in 1000 water molecules or in a mixture of 200 water and 800 1-octanol molecules, respectively, the latter described by the GAFF-DC parameters [37]. The scaling of the non-bonded interactions between the ligand and its environment was controlled via a coupling parameter � , such that � = 0 and � = 1 represents the fully interacting and fully decoupled ligand, respectively, while retaining the intramolecular interactions. The decoupling was conducted in a sequence of 20 discrete steps, using simulations times of 10 ns per �-state. In the applied perturbation scheme, elec- trostatic interactions were deactivated first within 5 steps, followed by the deactivation of the Lennard-Jones interac- tions. To avoid numerical problems close to the end states, soft-core (sc) potentials were used with parameters �sc = 0.5 , �sc = 0.3 nm and a power for the soft-core scaling function of psc = 1 [26]. The free energy changes were estimated from the sampled potential energy differences between all �-states using the MBAR estimator [38] as implemented in a freely available Python program [3]. Analysis Free energy profiles were analyzed using the umbrella inte- gration (UI) method [39–41] and compared to results from the weighted histogram analysis method (WHAM) [42–44] which showed no substantial differences. If a free energy offset between the two flat bulk water regions of the free energy profiles was present, it was found to be relatively small (below 8 kJmol−1 ) compared to the global well depth in all cases, except for the largest ligand G2 ( ≈ 35 kJmol−1 ), probably pointing towards non-sufficient sampling. In cases for which the offset exceeded the thermal noise threshold ( = RT ), binding free energies were calcu- lated from free energy profiles estimated with the modified WHAM algorithm of Hub et al. [45], using an additional constraint to suppress such offsets, noting that this is an ad- hoc solution that may lead to deviation from the standard binding free energy obtained from an offset-free free energy profile [4]. The standard binding free energy was obtained from the well depth of the free energy profile after correcting for the orthogonal restraint and the standard state concentra- tion as described previously [4, 46, 47]. Empirical correction From the results of previous SAMPL-challenges [48] and from own test calculations using guest molecules of known bind- ing affinity, it became evident that the used force field over- estimates binding substantially. This was addressed in some SAMPL6 contributions by a correction term of the form where the slope and offset coefficients (i.e., a and b respec- tively) were trained on data generated for previous rounds of the challenge [48]. While this correction could reduce the RMSE, it did not appreciably impact correlation statistics [48]. Therefore, an alternative approach was considered in the present work. The raw binding free energy values were corrected by an additive term such that The functional form of the (possibly very approximate) cor- rection term �Gcorr was chosen to be similar to the expres- sion for the binding free energy in the Linear Interaction Energy approach [49], comprising a contribution from van der Waals interactions and one associated with electrostatic interactions, approximated by a simple descriptor, In this equation, Nheavy denotes the number of heavy atoms that interact with the host molecule, VvdW hg represents the average short-range LJ interactions between host and guest in the minimum of the free-energy profile and �SASA is the change in solvent accessible surface area upon binding. For the set of five training molecules the ratio VvdW hg ∕�SASA turned out to be almost constant. TPSA describes the topo- logical polar surface area obtained from PubChem while � and � are fit parameters, which were determined based on a set of five guest molecules with known binding affinity towards CB8. Since the training set comprised relatively small molecules that fit well into the CB8 cavity, the above expression was slightly altered for the set of seven drug mol- ecules such that an effective number of heavy atoms was defined by where the denominator represents the average short-range LJ interactions per heavy atom for the training set and equals −2.94 kJmol−1 atom−1. (1)�Gcorr = a�Gcalc + b (2)�Gcorr pred = �Graw sim − �Gcorr (3)�Gcorr = Nheavy ( VvdW hg �SASA ) ⋅ � + TPSA ⋅ � (4)Neff heavy = VvdW hg ⟨VvdW hg ∕Nheavy⟩training 5Journal of Computer-Aided Molecular Design (2022) 36:1–9 1 3 For the set of five training molecules the contribution of the first term to the total correction in Eq. (3) was 97.6%, 81.2%, 70.7%, 71.6% and 75.3%. Results and discussion An initial assessment of force field and methodology to cal- culate binding free energies was attempted based on experi- mental data from the literature regarding the hydration free energy of G1 as well as on binding free energies of G1, G2, G5, G6 and G7 to CB7. Subsequently, the binding free ener- gies of guests GT1 to GT5 to CB8 was calculated and used to parametrize an empirical correction model that was then applied to the SAMPL8 set of guest molecules. Solvation free energies For the neutral form of G1 (methamphetamine) a hydra- tion free energy of −20.2 kJmol−1 was calculated. Because experimental values were not available directly, estimates were obtained from the relation [50] where cs G1 is the aqueous solubility expressed in molar con- centration, R is the universal gas constant, T is the tempera- ture, psat G1 is the vapor pressure of G1 in equilibrium with pure condensed G1 and po is the pressure (24.77 bar) of an ideal gas at 1 molar concentration and 298.15 K. Using a solu- bility from PubChem of 0.928 g l−1 and different estimates of the vapor pressure ranging from 0.72 Pa (PubChem) to 38 Pa [51], a range of possible hydration free energies was obtained with values between −14.9 and −24.7 kJmol−1 . Other theoretical estimates reported in the literature are (5)csG1 = ( psat G1 po ) exp [ −�Ghyd G1 RT ] −6.65 kJmol−1 and −32.5 kJmol−1 coming from classical density functional calculation and the estimation program interface (EPI) SuiteTM, respectively [52]. For the solva- tion free energy in the hydrated octanol phase a value of −34.95 kJmol−1 was calculated in the present work, lead- ing to an octanol-water partition coefficient of 2.6, which compares reasonably well with the value of 2.07 reported on PubChem, which is, however, not an experimental esti- mate. Due to a lack of experimental data for the other com- pounds, additional solvation free energy simulations were not conducted. Binding free energies of SAMPL8 guests to CB7 Binding constants for some guest molecules to CB7 were measured by direct 1 H NMR titration (G5, G6, G7), 1 H NMR competition assay (G1) or ITC (G2), respectively [53]. The binding free energy calculations showed some sampling problems at the cavity entrance, illustrated by the gaps in the time series of replica exchanges, see Fig. 3. For the more bulky guest molecules these sampling issues became more severe such that the simulated binding free energies are likely contaminated by insufficient sampling. The results are reported in Table 1 and show a rather het- erogeneous picture that does not allow for an unambiguous assessment of force-field adequacy. The free-energy profiles are shown in the Supplementary Information (see Fig. S2). Binding free energies of training set to CB8 In contrast to CB7, no significant sampling problems at the cavity entrance were observed for the CB8 host. The free-energy profiles of the five training guest molecules are shown in the Supplementary Information (see Fig. S3). In 0 10 20 30 40 50 time / ns 0 10 20 30 40 re pl ic a in de x -2 -1 0 1 2 z-coordinate / nm 0 20 40 60 80 fr ee e ne rg y / k J m ol -1 Fig. 3 Replica exchanges over simulation time (left) and corresponding free energy profile (right) for the system G1 (protonated form) binding to CB7. Differently colored traces represent different replicas 6 Journal of Computer-Aided Molecular Design (2022) 36:1–9 1 3 some cases, e.g. GT1, some offset was observed between the two bulk regions. However, compared to the profile well depth this effect is of minor importance. Fig. 4 shows that the calculated binding free energies correlate well with experiment data but are consistently too negative. By means of Eq. (3), these raw binding free ener- gies were corrected, resulting in an almost perfect agree- ment with the experimental data. For GT3, GT4 and GT5 the advantage of individualized contributions to the correction term are visible as the corrected data are much less scattered compared to the raw data. Binding free energies of SAMPL8 guests to CB8 Free-energy profiles of the seven challenge guest molecules are shown in the Supplementary Information (see Fig. S4). Fig. 5 shows the correlation of raw and corrected binding free energies with experimental values. Note that for the corrected values, Eq. (4) was used to calculate an effective number of heavy atoms, while the parameters � and � enter- ing Eq. (3) were transferred from the training set discussed in the previous paragraph. In contrast to the training mol- ecules, the slope of the linear regression curve is larger than one. Unfortunately this slope even gets larger when applying the corrections. Therefore, the systematic error decreased but the statistical correlation did not improve. If the correction model according to Eq. (3) is re-evalu- ated using the actual number of heavy atoms instead of the effective number according to Eq. (4), the correlation with experiment improves for those molecules which fit well into the CB8 cavity, i.e. G3, G4 and G5. For G1, the difference between the actual and effective number of heavy atoms is only 0.2. Only for the molecules that do not fit well into the CB8 cavity, i.e. G2 and G7, the evaluation of the correction model with the effective number of heavy atoms leads to bet- ter agreement with experiment. For G6, the predicted value gets closer to experiment when using the actual number of heavy atoms but is still far off. The results of this retrospec- tive analysis are displayed in Fig. 6 and may point to an improved correction model for future applications. We also note, that the present results, together with other SAMPL contributions could be re-analysed by efficient reweight- ing schemes [54] to study the sensitivity of the binding free energy with respect to the various force field parameters. Conclusions Based on previous rounds of the SAMPL challenge in which the GAFF force field overestimated binding affini- ties to the CB8 host [48], quantitative agreement between simulation and experiment from the raw simulation results was not expected in the present study. An empirical cor- rection model that goes beyond a simple linear correc- tion by incorporating specific properties of the guest mol- ecules as well as of their interactions with the host was proposed. The model worked rather well for the small set of five training molecules. For the seven guest molecules Table 1 Experimental [53] and simulated binding free energies ( kJmol−1 ) of SAMPL8 guest molecules to CB7. Simulated values are presented for the neutral/protonated form, respectively The uncertainty in the simulated values resulting from free-energy profile offsets and statistical errors estimated from fluctuations of the sampled order parameter via error propagation according to the UI method [41] are on the order of 4.2 kJmol−1 . Cases for which no stable bound state (i.e. with the ligand binding either to the inner or outer host surface) could be determined within a series of 3 independ- ent unbiased simulations of 2 ns simulation time, are referred to as as “no binding” Guest molecule �Gexp �Gsim G1 (Methamphetamine) −46.1 −35.7 / −55.6 G2 (Fentanyl) −41.4 No binding G3 (Morphine) No binding No binding/−16.1 G4 (Hydromorphone) No binding No binding/−17.1 G5 (Ketamine) −16.0 No binding/−5.8(R); −7.4(S) G6 (Phencyclidine) −20.8 No binding/−22.3 G7 (Cocaine) −19.2 −7.5 / −44.0 Fig. 4 Correlation between calculated and experimental binding free energies for the set of five training molecules with and without an empirical correction. The model parameters required to evalu- ate Eq.  (3) are provided in the Supplementary Information. The fit parameters were � = 17.8Å2 and � = −0.21 kJmol−1 Å−2 . The uncer- tainty in the simulated values resulting from free-energy profile off- sets and statistical errors estimated from fluctuations of the sampled order parameter via error propagation according to the UI method [41] are on the order of 4.2 kJmol−1 7Journal of Computer-Aided Molecular Design (2022) 36:1–9 1 3 of the SAMPL8 challenge it reduced the systematic error but did not improve the statistical correlation. The good agreement between simulation and experiment obtained by Henchman and co-workers in the SAMPL8 challenge [55] using the GAFF2 force field [9] suggest to reconsider the US simulations in future work to explore the effect of the re-optimized covalent and Lennard-Jones interactions (possibly combined with a revised charge model [56]) rela- tive to the GAFF force field. Umbrella sampling combined with Hamiltonian replica exchange turned out to be a sim- ple and robust simulation protocol, as was also concluded by other participants of the SAMPL8 challenge [57]. Supplementary Information The online version of this article con- tains supplementary material available https:// doi. org/ 10. 1007/ s10822- 021- 00439-w. Acknowledgements This work was funded by Deutsche Forschungsge- meinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC 2075 - 390740016 and by Ministerium für Wissenschaft, Forschung und Kunst Baden-Württemberg. We acknowledge the support by the Stuttgart Center for Simulation Sci- ence (SimTech), the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the State of Baden-Württemberg through bwHPC and the DFG through grant no INST 37/935-1 FUGG. We appreciate the National Institutes of Health for its support of the SAMPL project via R01GM124270 to David L. Mobley (UC Irvine). Funding Open Access funding enabled and organized by Projekt DEAL. Data availability The datasets generated and analysed in this work are available in the data repository of the University of Stuttgart (DaRUS): https:// doi. org/ 10. 18419/ darus- 2109. Open Access This article is licensed under a Creative Commons Attri- bution 4.0 International License, which permits use, sharing, adapta- tion, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. Fig. 5 Correlation between calculated and experimental binding free energies for the set of seven challenge molecules with and without an empirical correction. The model parameters required to evalu- ate Eq.  (3) are provided in the Supplementary Information. The fit parameters were � = 17.8Å2 and � = −0.21 kJmol−1 Å−2 . The uncer- tainty in the simulated values resulting from free-energy profile off- sets and statistical errors estimated from fluctuations of the sampled order parameter via error propagation according to the UI method [41] are on the order of 4.2 kJmol−1 Fig. 6 Correlation between calculated and experimental binding free energies for the set of five training molecules and seven challenge molecules after applying Eq.  (3) with the actual number of heavy atoms, except for molecules G2 and G7 for which the same effective number was used as in Fig.  5. The fit parameters were � = 17.8Å2 and � = −0.21 kJmol−1 Å−2 . The uncertainty in the simulated val- ues resulting from free-energy profile offsets and statistical errors estimated from fluctuations of the sampled order parameter via error propagation according to the UI method [41] are on the order of 4.2 kJmol−1 https://doi.org/10.1007/s10822-021-00439-w https://doi.org/10.1007/s10822-021-00439-w https://doi.org/10.18419/darus-2109 http://creativecommons.org/licenses/by/4.0/ 8 Journal of Computer-Aided Molecular Design (2022) 36:1–9 1 3 References 1. Chipot C, Pohorille A (eds) (2007) Free energy calculations. Springer, Berlin 2. Torrie GM, Valleau JP (1977) J. Comput. Phys. 23:187. https:// doi. org/ 10. 1016/ 0021- 9991(77) 90121-8 3. Klimovich PV, Shirts MR, Mobley DL (2015) J Comput Aided Mol Des 29:397. https:// doi. org/ 10. 1007/ s10822- 015- 9840-9 4. Markthaler D, Jakobtorweihen S, Hansen N (2019) Living J Comp Mol Sci 1(2):11073. https:// doi. org/ 10. 33011/ livec oms.1. 2. 11073 5. Schindler CEM, Baumann H, Blum A, Böse D, Buchstaller HP, Burgdorf L, Cappel D, Chekler E, Czodrowski P, Dorsch D, Eguida MKI, Follows B, Fuchß T, Grädler U, Gunera J, Johnson T, Jorand Lebrun C, Karra S, Klein M, Knehans T, Koetzner L, Krier M, Leiendecker M, Leuthner B, Li L, Mochalkin I, Musil D, Neagu C, Rippmann F, Schiemann K, Schulz R, Steinbrecher T, Tanzer EM, Unzue Lopez A, Viacava Follis A, Wegener A, Kuhn D (2020) J Chem Inf Model 60(11):5457. https:// doi. org/ 10. 1021/ acs. jcim. 0c009 00 6. King E, Aitchison E, Li H, Luo R (2021) Front Mol Biosci 8:712085. https:// doi. org/ 10. 3389/ fmolb. 2021. 712085 7. Mobley DL, Amezcua M (2020) SAMPL8 CB8 “drugs of abuse” challenge inputs (2020). https:// doi. org/ 10. 5281/ zenodo. 40295 60 8. Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA (2004) J Comput Chem 25(9):1157. https:// doi. org/ 10. 1002/ jcc. 20035 9. Case DA, Belfon K, Ben-Shalom IY, Brozell SR, Cerutti DS, Cheatham TE III, Cruzeiro VWD, Darden TA, Duke RE, Giambasu G, Gilson MK, Gohlke H, Goetz AW, Harris R, Izadi S, Izmailov SA, Kasavajhala K, Kovalenko A, Krasny R, Kurtzman T, Lee TS, LeGrand S, Li P, Lin C, Liu J, Luchko T, Luo R, Man V, Merz KM, Miao Y, Mikhailovskii O, Monard G, Nguyen H, Onufriev A, Pan F, Pantano S, Qi R, Roe DR, Roitberg A, Sagui C, Schott-Verdugo S, Shen J, Simmerling CL, Skrynnikov NR, Smith J, Swails J, Walker RC, Wang J, Wilson L, Wolf RM, Wu X, Xiong Y, Xue Y, York DM, Kollman PA (2020) AMBER 2020 10. Jakalian A, Bush BL, Jack DB, Bayly CI (2000) J Comput Chem 21(2):132 11. Jakalian A, Jack DB, Bayly CI (2002) J Comput Chem 23(16):1623. https:// doi. org/ 10. 1002/ jcc. 10128 12. Shirts MR, Klein C, Swails JM, Yin J, Gilson MK, Mobley DL, Case DA, Zhong ED (2017) J Comput Aided Mol Des 31(1):147. https:// doi. org/ 10. 1007/ s10822- 016- 9977-1 13. https:// amber md. org/ tutor ials/ advan ced/ tutor ial21/ files/ cb7_ am1- bcc. mol2 14. Ahlrichs R, Bär M, Häser M, Horn H, Kölmel C (1989) Chem Phys Lett 162:165 15. TURBOMOLE v.7.2.1 2017, A development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2007, TURBOMOLE GmbH, since 2007. http:// www. turbo mole. com 16. Balasubramani SG, Chen GP, Coriani S, Diedenhofen M, Frank MS, Franzke YJ, Furche F, Grotjahn R, Harding ME, Hättig C, Hell- weg A, Helmich-Paris B, Holzer C, Huniar U, Kaupp M, Marefat Khah A, Karbalaei Khani S, Müller T, Mack F, Nguyen BD, Parker SM, Perlt E, Rappoport D, Reiter K, Roy S, Rückert M, Schmitz G, Sierka M, Tapavicza E, Tew DP, van Wüllen C, Voora VK, Weigend F, Wodyński A, Yu JM (2020) J Chem Phys 152(18):184107. https:// doi. org/ 10. 1063/5. 00046 35 17. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML (1983) J Chem Phys 79(2):926. https:// doi. org/ 10. 1063/1. 445869 18. Joung IS, Cheatham TE (2008) J Phys Chem B 112(30):9020. https:// doi. org/ 10. 1021/ jp800 1614 19. Manz TA, Limas NG (2016) RSC Adv 6:47771 20. Murkli S, Klemm J, Brockett AT, Shuster M, Briken V, Roesch MR, Isaacs L (2021) Chem Eur J 27(9):3098. https:// doi. org/ 10. 1002/ chem. 20200 4380 21. Hippler M, Metcalfe GD (2020) Bunsen-Magazin 22(5), 102 (2020). https:// doi. org/ 10. 26125/ y7p7- an56 22. Barton SC, https:// www. egr. msu. edu/ ~scb- group- web/ buffe rs/ buffe rs. html 23. Kashefolgheta S, Vila Verde A (2017) Phys Chem Chem Phys 19:20593. https:// doi. org/ 10. 1039/ C7CP0 2557B 24. Donnini S, Mark AE, Juffer AH, Villa A (2005) J Comput Chem 26(2):115. https:// doi. org/ 10. 1002/ jcc. 20156 25. König G, Brooks BR (2012) J Comput Aided Mol Des 26:543. https:// doi. org/ 10. 1007/ s10822- 011- 9525-y 26. Abraham MJ, Murtola T, Schulz R, Páll S, Smith JC, Hess B, Lin- dahl E (2015) SoftwareX 1–2:19. https:// doi. org/ 10. 1016/j. softx. 2015. 06. 001 27. Tribello GA, Bonomi M, Branduardi D, Camilloni C, Bussi G (2014) Comput Phys Commun 185(2):604. https:// doi. org/ 10. 1016/j. cpc. 2013. 09. 018 28. Sugita Y, Kitao A, Okamoto Y (2000) J Chem Phys 113(15):6042. https:// doi. org/ 10. 1063/1. 13085 16 29. Bussi G (2014) Mol Phys 112(3–4):379. https:// doi. org/ 10. 1080/ 00268 976. 2013. 824126 30. Goga N, Rzepiela AJ, de Vries AH, Marrink SJ, Berendsen HJC (2012) J Chem Theory Comput 8(10):3637. https:// doi. org/ 10. 1021/ ct300 0876 31. Parrinello M, Rahman A (1981) J Appl Phys 52(12):7182. https:// doi. org/ 10. 1063/1. 328693 32. Nosé S, Klein M (1983) Mol Phys 50(5):1055. https:// doi. org/ 10. 1080/ 00268 97830 01028 51 33. Páll S, Hess B (2013) Comput Phys Commun 184(12):2641. https:// doi. org/ 10. 1016/j. cpc. 2013. 06. 003 34. Darden T, York D, Pedersen L (1993) J Chem Phys 98(12):10089. https:// doi. org/ 10. 1063/1. 464397 35. Essmann U, Perera L, Berkowitz M, Darden T, Lee H, Pedersen L (1995) J Chem Phys 103(19):8577. https:// doi. org/ 10. 1063/1. 470117 36. Ryckaert JP, Ciccotti G, Berendsen HJC (1977) J Comput Phys 23:327. https:// doi. org/ 10. 1016/ 0021- 9991(77) 90098-5 37. Fennel CJ, Wymer KL, Mobley DL (2014) J Phys Chem B 118:6438. https:// doi. org/ 10. 1021/ jp411 529h 38. Shirts MR, Chodera JD (2008) J Chem Phys 129:124105. https:// doi. org/ 10. 1063/1. 29781 77 39. Billeter S, van Gunsteren W (2000) J Phys Chem A 104(15):3276. https:// doi. org/ 10. 1021/ jp994 127q 40. Kästner J, Thiel W (2005) J Chem Phys 123(14):144104. https:// doi. org/ 10. 1063/1. 20526 48 41. Kästner J, Thiel W (2006) J Chem Phys 124(23):234106. https:// doi. org/ 10. 1063/1. 22067 75 42. Ferrenberg AM, Swendsen RH (1989) Phys Rev Lett 63(12):1195. https:// doi. org/ 10. 1103/ PhysR evLett. 63. 1195 43. Kumar S, Rosenberg JM, Bouzida D, Swendsen RH, Kollman PA (1992) J Comput Chem 13(8):1011. https:// doi. org/ 10. 1002/ jcc. 54013 0812 44. Roux B (1995) Comput Phys Commun 91(1–3):275. https:// doi. org/ 10. 1016/ 0010- 4655(95) 00053-I 45. Hub JS, De Groot BL, Van Der Spoel D (2010) J Chem Theory Comput 6(12):3713. https:// doi. org/ 10. 1021/ ct100 494z 46. Doudou S, Burton NA, Henchman RH (2009) J Chem Theory Com- put 5(4):909. https:// doi. org/ 10. 1021/ ct800 2354 47. Markthaler D, Gebhardt J, Jakobtorweihen S, Hansen N (2017) Chem Ing Tech 89(10):1306. https:// doi. org/ 10. 1002/ cite. 20170 0057 48. Rizzi A, Murkli S, McNeill JN, Yao W, Sullivan M, Gilson MK, Chiu MW, Isaacs L, Gibb BC, Mobley DL, Chodera JD (2018) J Comput Aided Mol Des 32:937. https:// doi. org/ 10. 1007/ s10822- 018- 0170-6 49. Åqvist J, Medina C, Samuelsson JE (1994) Protein Eng 7(3):385. https:// doi. org/ 10. 1093/ prote in/7. 3. 385 https://doi.org/10.1016/0021-9991(77)90121-8 https://doi.org/10.1016/0021-9991(77)90121-8 https://doi.org/10.1007/s10822-015-9840-9 https://doi.org/10.33011/livecoms.1.2.11073 https://doi.org/10.1021/acs.jcim.0c00900 https://doi.org/10.1021/acs.jcim.0c00900 https://doi.org/10.3389/fmolb.2021.712085 https://doi.org/10.5281/zenodo.4029560 https://doi.org/10.1002/jcc.20035 https://doi.org/10.1002/jcc.10128 https://doi.org/10.1007/s10822-016-9977-1 https://doi.org/10.1007/s10822-016-9977-1 https://ambermd.org/tutorials/advanced/tutorial21/files/cb7_am1-bcc.mol2 https://ambermd.org/tutorials/advanced/tutorial21/files/cb7_am1-bcc.mol2 http://www.turbomole.com https://doi.org/10.1063/5.0004635 https://doi.org/10.1063/5.0004635 https://doi.org/10.1063/1.445869 https://doi.org/10.1021/jp8001614 https://doi.org/10.1002/chem.202004380 https://doi.org/10.1002/chem.202004380 https://doi.org/10.26125/y7p7-an56 https://www.egr.msu.edu/%7escb-group-web/buffers/buffers.html https://www.egr.msu.edu/%7escb-group-web/buffers/buffers.html https://doi.org/10.1039/C7CP02557B https://doi.org/10.1002/jcc.20156 https://doi.org/10.1007/s10822-011-9525-y https://doi.org/10.1016/j.softx.2015.06.001 https://doi.org/10.1016/j.softx.2015.06.001 https://doi.org/10.1016/j.cpc.2013.09.018 https://doi.org/10.1016/j.cpc.2013.09.018 https://doi.org/10.1063/1.1308516 https://doi.org/10.1080/00268976.2013.824126 https://doi.org/10.1080/00268976.2013.824126 https://doi.org/10.1021/ct3000876 https://doi.org/10.1021/ct3000876 https://doi.org/10.1063/1.328693 https://doi.org/10.1063/1.328693 https://doi.org/10.1080/00268978300102851 https://doi.org/10.1080/00268978300102851 https://doi.org/10.1016/j.cpc.2013.06.003 https://doi.org/10.1016/j.cpc.2013.06.003 https://doi.org/10.1063/1.464397 https://doi.org/10.1063/1.470117 https://doi.org/10.1016/0021-9991(77)90098-5 https://doi.org/10.1021/jp411529h https://doi.org/10.1063/1.2978177 https://doi.org/10.1063/1.2978177 https://doi.org/10.1021/jp994127q https://doi.org/10.1063/1.2052648 https://doi.org/10.1063/1.2052648 https://doi.org/10.1063/1.2206775 https://doi.org/10.1063/1.2206775 https://doi.org/10.1103/PhysRevLett.63.1195 https://doi.org/10.1002/jcc.540130812 https://doi.org/10.1002/jcc.540130812 https://doi.org/10.1016/0010-4655(95)00053-I https://doi.org/10.1016/0010-4655(95)00053-I https://doi.org/10.1021/ct100494z https://doi.org/10.1021/ct8002354 https://doi.org/10.1002/cite.201700057 https://doi.org/10.1002/cite.201700057 https://doi.org/10.1007/s10822-018-0170-6 https://doi.org/10.1007/s10822-018-0170-6 https://doi.org/10.1093/protein/7.3.385 9Journal of Computer-Aided Molecular Design (2022) 36:1–9 1 3 50. Thompson JD, Cramer CJ, Truhlar DG (2003) J Chem Phys 119(3):1661. https:// doi. org/ 10. 1063/1. 15794 74 51. Thornton M, Gobble C, Chickos J (2014) J Chem Thermodyn 73:51. https:// doi. org/ 10. 1016/j. jct. 2013. 08. 005 52. Li J, Fu J, Huang X, Lu D, Wu J (2016) J. Phys. Condens. Matter 28:344001. https:// doi. org/ 10. 1088/ 0953- 8984/ 28/ 34/ 344001 53. Ganapati S, Grabitz SD, Murkli S, Scheffenbichler F, Rudolph MI, Zavalij PY, Eikermann M, Isaacs L (2017) ChemBioChem 18(16):1583. https:// doi. org/ 10. 1002/ cbic. 20170 0289 54. Messerly RA, Razavi SM, Shirts MR (2018) J Chem Theory Com- put 14(6):3144. https:// doi. org/ 10. 1021/ acs. jctc. 8b002 23 55. Ali HS, Chakravorty A, Kalayan J, de Visser SP, Henchman RH (2021) J Comput Aided Mol Des 35:911. https:// doi. org/ 10. 1007/ s10822- 021- 00406-5 56. He X, Man VH, Yang W, Lee TS, Wang J (2020) J Chem Phys 153:114502. https:// doi. org/ 10. 1063/5. 00190 56 57. Ghorbani M, Hudson PS, Jones MR, Aviat F, Meana-Pañeda R, Klauda JB, Brooks BR (2021) J Comput Aided Mol Des 35:667. https:// doi. org/ 10. 1007/ s10822- 021- 00385-7 Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. https://doi.org/10.1063/1.1579474 https://doi.org/10.1016/j.jct.2013.08.005 https://doi.org/10.1088/0953-8984/28/34/344001 https://doi.org/10.1002/cbic.201700289 https://doi.org/10.1021/acs.jctc.8b00223 https://doi.org/10.1007/s10822-021-00406-5 https://doi.org/10.1007/s10822-021-00406-5 https://doi.org/10.1063/5.0019056 https://doi.org/10.1007/s10822-021-00385-7 Binding free energies for the SAMPL8 CB8 “Drugs of Abuse” challenge from umbrella sampling combined with Hamiltonian replica exchange Abstract Introduction Computational approach Molecular model System preparation Molecular dynamics simulations Analysis Empirical correction Results and discussion Solvation free energies Binding free energies of SAMPL8 guests to CB7 Binding free energies of training set to CB8 Binding free energies of SAMPL8 guests to CB8 Conclusions Acknowledgements References