Correcting electrostatic artifacts due to net-charge changes in the calculation of ligand binding free energies

Alchemically derived free energies are artifacted when the perturbed moiety has a nonzero net charge. The source of the artifacts lies in the effective treatment of the electrostatic interactions within and between the perturbed atoms and remaining (partial) charges in the simulated system. To treat the electrostatic interactions effectively, lattice-summation (LS) methods or cutoff schemes in combination with a reaction-field contribution are usually employed. Both methods render the charging component of the calculated free energies sensitive to essential parameters of the system like the cutoff radius or the box side lengths. Here, we discuss the results of three previously published studies of ligand binding. These studies presented estimates of binding free energies that were artifacted due to the charged nature of the ligands. We show that the size of the artifacts can be efficiently calculated and raw simulation data can be corrected. We compare the corrected results with experimental estimates and nonartifacted estimates from path-sampling methods. Although the employed correction scheme involves computationally demanding continuum-electrostatics calculations, we show that the correction estimate can be deduced from a small sample of configurations rather than from the entire ensemble. This observation makes the calculations of correction terms feasible for complex biological systems. To show the general applicability of the proposed procedure, we also present results where the correction scheme was used to correct independent free energies obtained from simulations employing a cutoff scheme or LS electrostatics. In this work, we give practical guidelines on how to apply the appropriate corrections easily.     

Modeling competitive guest binding to beta-cyclodextrin molecular printboards

Anchoring of functionalized guest molecules to self-assembled monolayers (SAMs) is key to the development of molecular printboards for nanopatterning. One very promising system involves guest binding to immobilized beta-cyclodextrin (beta-CD) hosts, with guest:host recognition facilitated by a hydrophobic interaction between uncharged anchor groups on the guest molecule and beta-CD hosts self-assembled at gold surfaces. We use molecular dynamics free energy (MDFE) simulations to describe the specificity of guest:beta-CD association. We find good agreement with experimental thermodynamic measurements for binding enthalpy differences between three commonly used phenyl guests: benzene, toluene, and t-butylbenzene. van der Waals interaction with the inside of the host cavity accounts for almost all of the net stabilization of the larger phenyl guests in beta-CD. Partial and full methylation of the secondary rim of beta-CD decreases host rigidity and significantly impairs binding of both phenyl and larger adamantane guest molecules. The beta-CD cavity is also very intolerant of guest charging, penalizing the oxidized state of ferrocene by at least 7 kcal/mol. beta-CD hence expresses moderate specificity toward uncharged organic guest molecules by van der Waals recognition, with a much higher specificity calculated for electrostatic recognition of organometallic guests.     

SAMPL6 host–guest challenge: binding free energies via a multistep approach

In this effort in the SAMPL6 host–guest binding challenge, a combination of molecular dynamics and quantum mechanical methods were used to blindly predict the host–guest binding free energies of a series of cucurbit[8]uril (CB8), octa-acid (OA), and tetramethyl octa-acid (TEMOA) hosts bound to various guest molecules in aqueous solution. Poses for host–guest systems were generated via molecular dynamics (MD) simulations and clustering analyses. The binding free energies for the structures obtained via cluster analyses of MD trajectories were calculated using the MMPBSA method and density functional theory (DFT) with the inclusion of Grimme’s dispersion correction, an implicit solvation model to model the aqueous solution, and the resolution-of-the-identity (RI) approximation (MMPBSA, RI-B3PW91-D3, and RI-B3PW91, respectively). Among these three methods tested, the results for OA and TEMOA systems showed MMPBSA and RI-B3PW91-D3 methods can be used to qualitatively rank binding energies of small molecules with an overbinding by 7 and 37 kcal/mol respectively, and RI-B3PW91 gave the poorest quality results, indicating the importance of dispersion correction for the binding free energy calculations. Due to the complexity of the CB8 systems, all of the methods tested show poor correlation with the experimental results. Other quantum mechanical approaches used for the calculation of binding free energies included DFT without the RI approximation, utilizing truncated basis sets to reduce the computational cost (memory, disk space, CPU time), and a corrected dielectric constant to account for ionic strength within the implicit solvation model.     

Molecular dynamics simulation study of the structural features and inclusion capacities of cucurbit[6]uril derivatives in aqueous solutions

Molecular dynamics (MD) simulations were performed for cucurbit[6]uril (CB6) methyl and cyclohexyl derivatives in aqueous solutions. Furthermore, MD simulations have been conducted to study the inclusion complexes between each CB6 derivative with α,ω-pentane diammonium ion (NH₃⁺(CH₂)₅NH₃⁺) to estimate the binding free energies, the complex geometries and the intermolecular forces responsible for complex formation. Results show a complete inclusion of the guest molecule in the cavity of the host for all complexes. Results also indicate that the guest dynamics inside the cavity of the substituted host is similar to that for the unsubstituted host. This demonstrates that the molecular recognition of the host is not affected by the alkyl substitution at the equator. Also, there is an insignificant conformational change of the macrocyclic structure upon inclusion of the guest. Molecular mechanics/Poisson Boltzmann surface area method was used to estimate the binding free energy of each complex. Results indicate that host–guest electrostatic interactions make the largest contribution to the complex binding free energy. Moreover, van der Waals interactions add significantly to the complex stability. The guest molecules show more or less similar binding free energies with the substituted CB6 that exhibits slightly more negative values than unsubstituted CB6 which is proved also by umbrella sampling.     

Overcoming dissipation in the calculation of standard binding free energies by ligand extraction

This article addresses calculations of the standard free energy of binding from molecular simulations in which a bound ligand is extracted from its binding site by steered molecular dynamics (MD) simulations or equilibrium umbrella sampling (US). Host–guest systems are used as test beds to examine the requirements for obtaining the reversible work of ligand extraction. We find that, for both steered MD and US, marked irreversibilities can occur when the guest molecule crosses an energy barrier and suddenly jumps to a new position, causing dissipation of energy stored in the stretched molecule(s). For flexible molecules, this occurs even when a stiff pulling spring is used, and it is difficult to suppress in calculations where the spring is attached to the molecules by single, fixed attachment points. We, therefore, introduce and test a method, fluctuation‐guided pulling, which adaptively adjusts the spring's attachment points based on the guest's atomic fluctuations relative to the host. This adaptive approach is found to substantially improve the reversibility of both steered MD and US calculations for the present systems. The results are then used to estimate standard binding free energies within a comprehensive framework, termed attach‐pull‐release, which recognizes that the standard free energy of binding must include not only the pulling work itself, but also the work of attaching and then releasing the spring, where the release work includes an accounting of the standard concentration to which the ligand is discharged. © 2013 Wiley Periodicals, Inc.     

Host Perturbation in a β‐Hydroquinone Clathrate Studied by Combined X‐ray/Neutron Charge‐Density Analysis: Implications for Molecular Inclusion in Supramolecular Entities

X‐ray/neutron (X/N) diffraction data measured at very low temperature (15 K) in conjunction with ab initio theoretical calculations were used to model the crystal charge density (CD) of the host–guest complex of hydroquinone (HQ) and acetonitrile. Due to pseudosymmetry, information about the ordering of the acetonitrile molecules within the HQ cavities is present only in almost extinct, very weak diffraction data, which cannot be measured with sufficient accuracy even by using the brightest X‐ray and neutron sources available, and the CD model of the guest molecule was ultimately based on theoretical calculations. On the other hand, the CD of the HQ host structure is well determined by the experimental data. The neutron diffraction data provide hydrogen anisotropic thermal parameters and positions, which are important to obtain a reliable CD for this light‐atom‐only crystal. Atomic displacement parameters obtained independently from the X‐ray and neutron diffraction data show excellent agreement with a |ΔU| value of 0.00058 Ų indicating outstanding data quality. The CD and especially the derived electrostatic properties clearly reveal increased polarization of the HQ molecules in the host–guest complex compared with the HQ molecules in the empty HQ apohost crystal structure. It was found that the origin of the increased polarization is inclusion of the acetonitrile molecule, whereas the change in geometry of the HQ host structure following inclusion of the guest has very little effect on the electrostatic potential. The fact that guest inclusion has a profound effect on the electrostatic potential suggests that nonpolarizable force fields may be unsuitable for molecular dynamics simulations of host–guest interaction (e.g., in protein–drug complexes), at least for polar molecules.     

Effects of system net charge and electrostatic truncation on all‐atom constant pH molecular dynamics

Constant pH molecular dynamics offers a means to rigorously study the effects of solution pH on dynamical processes. Here, we address two critical questions arising from the most recent developments of the all‐atom continuous constant pH molecular dynamics (CpHMD) method: (1) What is the effect of spatial electrostatic truncation on the sampling of protonation states? (2) Is the enforcement of electrical neutrality necessary for constant pH simulations? We first examined how the generalized reaction field and force‐shifting schemes modify the electrostatic forces on the titration coordinates. Free energy simulations of model compounds were then carried out to delineate the errors in the deprotonation free energy and salt‐bridge stability due to electrostatic truncation and system net charge. Finally, CpHMD titration of a mini‐protein HP36 was used to understand the manifestation of the two types of errors in the calculated pK_a values. The major finding is that enforcing charge neutrality under all pH conditions and at all time via cotitrating ions significantly improves the accuracy of protonation‐state sampling. We suggest that such finding is also relevant for simulations with particle mesh Ewald, considering the known artifacts due to charge‐compensating background plasma. © 2014 Wiley Periodicals, Inc.     

The F130L mutation in streptavidin reduces its binding affinity to biotin through electronic polarization effect

Recently, Baugh et al. discovered that a distal point mutation (F130L) in streptavidin causes no distinct variation to the structure of the binding pocket but a 1000‐fold reduction in biotin binding affinity. In this work, we carry out molecular dynamics simulations and apply an end‐state free energy method to calculate the binding free energies of biotin to wild type streptavidin and its F130L mutant. The absolute binding affinities based on AMBER charge are repulsive, and the mutation induced binding loss is underestimated. When using the polarized protein‐specific charge, the absolute binding affinities are significantly enhanced. In particular, both the absolute and relative binding affinities are in line with the experimental measurements. Further investigation indicates that polarization effect is indispensable in both the generation of structural ensembles and the calculation of interaction energies. This work verifies Baugh's conjecture that electrostatic polarization effect plays an essential role in modulating the binding affinity of biotin to the streptavidin through F130L mutation. © 2013 Wiley Periodicals, Inc.     

Coupled semiempirical quantum mechanics and molecular mechanics (QM/MM) calculations on the aqueous solvation free energies of ionized molecules

The aqueous solvation free energies of ionized molecules were computed using a coupled quantum mechanical and molecular mechanical (QM/MM) model based on the AM1, MNDO, and PM3 semiempirical molecular orbital methods for the solute molecule and the TIP3P molecular mechanics model for liquid water. The present work is an extension of our model for neutral solutes where we assumed that the total free energy is the sum of components derived from the electrostatic/polarization terms in the Hamiltonian plus an empirical “nonpolar” term. The electrostatic/polarization contributions to the solvation free energies were computed using molecular dynamics (MD) simulation and thermodynamic integration techniques, while the nonpolar contributions were taken from the literature. The contribution to the electrostatic/polarization component of the free energy due to nonbonded interactions outside the cutoff radii used in the MD simulations was approximated by a Born solvation term. The experimental free energies were reproduced satisfactorily using variational parameters from the vdW terms as in the original model, in addition to a parameter from the one-electron integral terms. The new one-electron parameter was required to account for the short-range effects of overlapping atomic charge densities. The radial distribution functions obtained from the MD simulations showed the expected H-bonded structures between the ionized solute molecule and solvent molecules. We also obtained satisfactory results by neglecting both the empirical nonpolar term and the electronic polarization of the solute, i.e., by implementing a nonpolarization model. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1028–1038, 1999     

Computation of methodology-independent ionic solvation free energies from molecular simulations. I. The electrostatic potential in molecular liquids

The computation of ionic solvation free energies from atomistic simulations is a surprisingly difficult problem that has found no satisfactory solution for more than 15 years. The reason is that the charging free energies evaluated from such simulations are affected by very large errors. One of these is related to the choice of a specific convention for summing up the contributions of solvent charges to the electrostatic potential in the ionic cavity, namely, on the basis of point charges within entire solvent molecules (M scheme) or on the basis of individual point charges (P scheme). The use of an inappropriate convention may lead to a charge-independent offset in the calculated potential, which depends on the details of the summation scheme, on the quadrupole-moment trace of the solvent molecule, and on the approximate form used to represent electrostatic interactions in the system. However, whether the M or P scheme (if any) represents the appropriate convention is still a matter of on-going debate. The goal of the present article is to settle this long-standing controversy by carefully analyzing (both analytically and numerically) the properties of the electrostatic potential in molecular liquids (and inside cavities within them). Restricting the discussion to real liquids of "spherical" solvent molecules (represented by a classical solvent model with a single van der Waals interaction site), it is concluded that (i) for Coulombic (or straight-cutoff truncated) electrostatic interactions, the M scheme is the appropriate way of calculating the electrostatic potential; (ii) for non-Coulombic interactions deriving from a continuously differentiable function, both M and P schemes generally deliver an incorrect result (for which an analytical correction must be applied); and (iii) finite-temperature effects, including intermolecular orientation correlations and a preferential orientational structure in the neighborhood of a liquid-vacuum interface, must be taken into account. Applications of these results to the computation methodology-independent ionic solvation free energies from molecular simulations will be the scope of a forthcoming article.     

Theoretical studies on the effects of methods and parameterization on the calculated free energy of hydration for small molecules

Free energies of hydration (FEH) have been computed for 13 neutral and nine ionic species as a difference of theoretically calculated Gibbs free energies in solution and in the gas phase. In‐solution calculations have been performed using both SCIPCM and PCM polarizable continuum models at the density functional theory (DFT)/B3LYP and ab initio Hartree–Fock levels with two basis sets (6‐31G* and 6‐311++G**). Good linear correlation has been obtained for calculated and experimental gas‐phase dipole moments, with an increase by ～30% upon solvation due to solute polarization. The geometry distortion in solution turns out to be small, whereas solute polarization energies are up to 3 kcal/mol for neutral molecules. Calculation of free energies of hydration with PCM provides a balanced set of values with 6‐31G* and 6‐311++G** basis sets for neutral molecules and ionic species, respectively. Explicit solvent calculations within Monte Carlo simulations applying free energy perturbation methods have been considered for 12 neutral molecules. Four different partial atomic charge sets have been studied, obtained by a fit to the gas‐phase and in‐solution molecular electrostatic potentials at in‐solution optimized geometries. Calculated FEH values depend on the charge set and the atom model used. Results indicate a preference for the all‐atom model and partial charges obtained by a fit to the molecular electrostatic potential of the solute computed at the SCIPCM/B3LYP/6‐31G* level. © 2004 Wiley Periodicals, Inc. Int J Quantum Chem, 2004     

环双(对-蒽基-对草快)的分子识别与谱学性质

环双(对-蒽基-对草快)是一种新型的缺电子大环仿生主体, 分子识别是其最重要的应用之一. 考察主体对一系列客体分子如水、氨、醇及杂环等的识别能力, 用密度泛函理论(DFT)中的B3LYP/3-21G基组对主客体复合物的结构进行优化. 在B3LYP/6-31G(d)水平上进行单点能计算, 校正后得到复合物的结合能. 用B3LYP/3-21G方法计算13C和3He化学位移. 结果表明, 主体对客体分子的识别主要靠客体上的杂原子与主体上的氢原子之间的氢键进行. 复合物的稳定化能受氢键的数目和距离影响. 氢键的形成导致部分复合物LUMO与HOMO能隙增大, 同时导致与氢键相连的C—H键上C原子的化学位移向低场移动. 复合物的芳香性与其结合能的大小及结合方式有关. 主体的芳香性因其与客体之间的弱相互作用而提高, 但太强的相互作用及客体在主体空腔内都将影响主体的环电流, 从而削弱其芳香性.     

Guest Binding with Sulfated Cyclodextrins: Does the Size of Cavity Matter?

Sulfated cyclodextrins have recently emerged as potential candidates for producing host–induced guest aggregation with properties better than p-sulfonatocalixarenes that have previously shown numerous applications involving the phenomena of host-induced guest aggregation. In the class of sulfated cyclodextrins (SCD), sulfated β-cyclodextrin (β-SCD) remains the most extensively investigated host molecule. Although it is assumed that the host-induced guest aggregation is predominantly an outcome of interaction of the guest molecule with the charges on the exterior of SCD cavity, it has not been deciphered whether the variation in the cavity size will make a difference in the efficiency of host-induced guest-aggregation process. In this investigation, we present a systematic study of host–induced guest aggregation of a cationic molecular rotor dye, Thioflavin T (ThT) with three different sulfated cyclodextrin molecules, α-SCD, β-SCD and γ-SCD, which differ in their cavity size, using steady-state emission, ground-state absorption and time-resolved emission measurements. The obtained photophysical properties of ThT, upon interaction with different SCD molecules, indicate that the binding strength of ThT with different SCD molecules correlate with the cavity size of the host molecule, giving rise to the strongest complexation of ThT with the largest host molecule (γ-SCD). The binding affinity of ThT towards different host molecules has been supported by molecular docking calculations. The results obtained are further supported with the temperature and ionic strength dependent studies performed on the host-guest complex. Our results indicate that for host–induced guest aggregation, involving oppositely charged molecules, the size of the cavity also plays a crucial role beside the charge density on the exterior of host cavity.     

Theoretical study of the inclusion properties of 1,2,4,5-tetra(morpholinocarbonyl)benzene

Molecular inclusion by the new amide host molecule (TMB) has been reconsidered by calculating the crystal stabilization energies for the guest molecules in the TMB + guest system from the simple intermolecular potential functions of Caillet and Claverie. Water, ethylene glycol, methanol, and ethanol have been employed as guest molecules and their relative stabilities have been considered. Water has been found to be the most suitable guest molecule in the TMB + guest system. It also has been found that the guest host interaction is the most important contributor in determining the relative stabilities of the guest molecules in the TMB + guest system, but the guest guest interaction is very important, too. Moreover, the electrostatic interaction has been found to be the most important contributor to the total interaction energy in the TMB + guest system.     

Converging ligand‐binding free energies obtained with free‐energy perturbations at the quantum mechanical level

In this article, the convergence of quantum mechanical (QM) free‐energy simulations based on molecular dynamics simulations at the molecular mechanics (MM) level has been investigated. We have estimated relative free energies for the binding of nine cyclic carboxylate ligands to the octa‐acid deep‐cavity host, including the host, the ligand, and all water molecules within 4.5 Å of the ligand in the QM calculations (158–224 atoms). We use single‐step exponential averaging (ssEA) and the non‐Boltzmann Bennett acceptance ratio (NBB) methods to estimate QM/MM free energy with the semi‐empirical PM6‐DH2X method, both based on interaction energies. We show that ssEA with cumulant expansion gives a better convergence and uses half as many QM calculations as NBB, although the two methods give consistent results. With 720,000 QM calculations per transformation, QM/MM free‐energy estimates with a precision of 1 kJ/mol can be obtained for all eight relative energies with ssEA, showing that this approach can be used to calculate converged QM/MM binding free energies for realistic systems and large QM partitions. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Converging free energies of binding in cucurbit[7]uril and octa-acid host–guest systems from SAMPL4 using expanded ensemble simulations

Molecular containers such as cucurbit[7]uril (CB7) and the octa-acid (OA) host are ideal simplified model test systems for optimizing and analyzing methods for computing free energies of binding intended for use with biologically relevant protein–ligand complexes. To this end, we have performed initially blind free energy calculations to determine the free energies of binding for ligands of both the CB7 and OA hosts. A subset of the selected guest molecules were those included in the SAMPL4 prediction challenge. Using expanded ensemble simulations in the dimension of coupling host–guest intermolecular interactions, we are able to show that our estimates in most cases can be demonstrated to fully converge and that the errors in our estimates are due almost entirely to the assigned force field parameters and the choice of environmental conditions used to model experiment. We confirm the convergence through the use of alternative simulation methodologies and thermodynamic pathways, analyzing sampled conformations, and directly observing changes of the free energy with respect to simulation time. Our results demonstrate the benefits of enhanced sampling of multiple local free energy minima made possible by the use of expanded ensemble molecular dynamics and may indicate the presence of significant problems with current transferable force fields for organic molecules when used for calculating binding affinities, especially in non-protein chemistries.     

Finite difference Poisson-Boltzmann electrostatic calculations: Increased accuracy achieved by harmonic dielectric smoothing and charge antialiasing

A common problem in the calculation of electrostatic potentials with the Poisson-Boltzmann equation using finite difference methods is the effect of molecular position relative to the grid. Previously a uniform charging method was shown to reduce the grid dependence substantially over the point charge model used in commercially available codes. In this article we demonstrate that smoothing the charge and dielectric values on the grid can improve the grid independence, as measured by the spread of calculated values, by another order of magnitude. Calculations of Born ion solvation energies, small molecule solvation energies, the electrostatic field of superoxide dismutase, and protein-protein binding energies are used to demonstrate that this method yields the same results as the point charge model while reducing the positional errors by several orders of magnitude. © 1997 by John Wiley & Sons, Inc.     

Squeezing,Then Stacking: From Breathing Pores to Three‐Dimensional Ionic Self‐Assembly under Electrochemical Control

We demonstrate the spontaneous and reversible transition between the two‐ and three‐dimensional self‐assembly of a supramolecular system at the solid–liquid interface under electrochemical conditions, using in situ scanning tunneling microscopy. By tuning the interfacial potential, we can selectively organize our target molecules in an open porous pattern, fill these pores to form an auto‐host–guest structure, or stack the building blocks in a stratified bilayer. Using a simple electrostatic model, we rationalize which charge density is required to enable bilayer formation, and conversely, which molecular size/charge ratio is necessary in the design of new building blocks. Our results may lead to a new class of electrochemically controlled dynamic host–guest systems, artificial receptors, and smart materials.     

Theoretical Characterization of Charge Transport in One‐Dimensional Collinear Arrays of Organic Conjugated Molecules

A great deal of interest has recently focused on host–guest systems consisting of one‐dimensional collinear arrays of conjugated molecules encapsulated in the channels of organic or inorganic matrices. Such architectures allow for controlled charge and energy migration processes between the interacting guest molecules and are thus attractive in the field of organic electronics. In this context, we characterize here at a quantum‐chemical level the molecular parameters governing charge transport in the hopping regime in 1D arrays built with different types of molecules. We investigate the influence of several parameters (such as the symmetry of the molecule, the presence of terminal substituents, and the molecular size) and define on that basis the molecular features required to maximize the charge carrier mobility within the channels. In particular, we demonstrate that a strong localization of the molecular orbitals in push–pull compounds is generally detrimental to the charge transport properties.