Computed Interaction Energies Research Articles

Polarization effects in protein–ligand calculations extend farther than the actual induction energy

The various roles that polarizabilities play in the calculation of protein–ligand interaction energies with a polarizable force field are investigated, and the importance and distance dependence of some common approximations is determined for each of these roles separately, using quantum-mechanical calculations as the reference. It is found that the pure induction energy, if defined as the energetic gain from the charge redistribution upon interaction between the protein and ligand, is a rather short-ranged effect that becomes independent of the exact implementation at distances above ∼4 A. On the other hand, the polarization between the protein residues in the assembly of a protein from separately computed fragments (as is routinely done in force field development) has a significant effect on the computed interaction energies, even for residues as far as 15 A from the ligand. Finally, polarization improves the transferability of partial charges, but the simple polarization model used in, for example, the Amber force field explains only 14–19% of the conformational variation of the charges. In all cases, more advanced polarization models, especially involving anisotropic polarizabilities, seem to give significantly better descriptions of these effects. The study suggests that an accurate treatment of polarization can be important even in systems where the actual induction energy is small in magnitude.

Double kink mechanisms for discrete dislocations in BCC crystals

We present an application of the discrete dislocation theory to the characterization of the energetics of kinks in Mo, Ta and W body-centered cubic (BCC) crystals. The discrete dislocation calculations supply detailed predictions of formation and interaction energies for various double-kink formation and spreading mechanisms as a function of the geometry of the double kinks, including: the dependence of the formation energy of a double kink on its width; the energy of formation of a double kink on a screw dislocation containing a pre-existing double kink; and energy of formation of a double kink on a screw dislocation containing a pre-existing single kink. The computed interaction energies are expected to facilitates the nucleation of double kinks in close proximity to each other and to pre-existing kinks, thus promoting clustering of double kinks on screw segments and of 'daughter’ double kinks ahead of ‘mother’ kinks. The predictions of the discrete dislocation theory are found to be in good agreement with the full atomistic calculations based on empirical interatomic potentials available in the literature.

Novel high‐affinity binders of human interferon gamma derived from albumin‐binding domain of protein G

Recombinant ligands derived from small protein scaffolds show promise as robust research and diagnostic reagents and next generation protein therapeutics. Here, we derived high-affinity binders of human interferon gamma (hIFNγ) from the three helix bundle scaffold of the albumin-binding domain (ABD) of protein G from Streptococcus G148. Computational interaction energy mapping, solvent accessibility assessment, and in silico alanine scanning identified 11 residues from the albumin-binding surface of ABD as suitable for randomization. A corresponding combinatorial ABD scaffold library was synthesized and screened for hIFNγ binders using in vitro ribosome display selection, to yield recombinant ligands that exhibited K(d) values for hIFNγ from 0.2 to 10 nM. Molecular modeling, computational docking onto hIFNγ, and in vitro competition for hIFNγ binding revealed that four of the best ABD-derived ligands shared a common binding surface on hIFNγ, which differed from the site of human IFNγ receptor 1 binding. Thus, these hIFNγ ligands provide a proof of concept for design of novel recombinant binding proteins derived from the ABD scaffold.

Protein–Ligand Interaction Energies with Dispersion Corrected Density Functional Theory and High-Level Wave Function Based Methods

With dispersion-corrected density functional theory (DFT-D3) intermolecular interaction energies for a diverse set of noncovalently bound protein-ligand complexes from the Protein Data Bank are calculated. The focus is on major contacts occurring between the drug molecule and the binding site. Generalized gradient approximation (GGA), meta-GGA, and hybrid functionals are used. DFT-D3 interaction energies are benchmarked against the best available wave function based results that are provided by the estimated complete basis set (CBS) limit of the local pair natural orbital coupled-electron pair approximation (LPNO-CEPA/1) and compared to MP2 and semiempirical data. The size of the complexes and their interaction energies (ΔE(PL)) varies between 50 and 300 atoms and from -1 to -65 kcal/mol, respectively. Basis set effects are considered by applying extended sets of triple- to quadruple-ζ quality. Computed total ΔE(PL) values show a good correlation with the dispersion contribution despite the fact that the protein-ligand complexes contain many hydrogen bonds. It is concluded that an adequate, for example, asymptotically correct, treatment of dispersion interactions is necessary for the realistic modeling of protein-ligand binding. Inclusion of the dispersion correction drastically reduces the dependence of the computed interaction energies on the density functional compared to uncorrected DFT results. DFT-D3 methods provide results that are consistent with LPNO-CEPA/1 and MP2, the differences of about 1-2 kcal/mol on average (<5% of ΔE(PL)) being on the order of their accuracy, while dispersion-corrected semiempirical AM1 and PM3 approaches show a deviating behavior. The DFT-D3 results are found to depend insignificantly on the choice of the short-range damping model. We propose to use DFT-D3 as an essential ingredient in a QM/MM approach for advanced virtual screening approaches of protein-ligand interactions to be combined with similarly "first-principle" accounts for the estimation of solvation and entropic effects.

Unusually Long, Multicenter, Cationδ+⋅⋅⋅Anionδ−Bonding Observed for Several Polymorphs of [TTF][TCNE

The α, β, and δ polymorphs of [TTF][TCNE] (TTF=tetrathiafulvalene; TCNE=tetracyanoethylene) exhibit a new type of long, multicenter bonding between the [TTF](δ+) and [TCNE](δ-) moieties, demonstrating the existence of long, hetero-multicenter bonding with a cationic(δ+)···anionic(δ-) zwitterionic-like structure. These diamagnetic π-[TTF](δ+) [TCNE](δ-) heterodimers exhibit a transfer of about 0.5 e(-) from the TTF to the TCNE fragments, as observed from experimental studies, in accord with theoretical predictions, that is, [TTF(δ+)···TCNE(δ-)] (δ≅0.5). They have several interfragment distances <3.4 Å, and a computed interaction energy of -21.2 kcal mol(-1), which is typical of long, multicenter bonds. The lower stability of [TTF](δ+) [TCNE](δ-) with respect to typical ionic bonds is due, in part, to the partial electron transfer that reduces the electrostatic bonding component. This reduced electrostatic interaction, and the large interfragment dispersion stabilize the long, heterocationic/anionic multicenter interaction, which in [TTF(δ+)···TCNE(δ-)] always involves two electrons, but have ten, eight, and eight bond critical points (bcps) involving C-C, N-S, and sometimes C-S and C-N components for the α, β, and δ polymorphs, respectively. In contrast, γ-[TTF][TCNE] possesses [TTF](2)(2+) and [TCNE](2)(2-) dimers, each with long, homo-multicenter 2e(-)/12c (c=center, 2 C+4 S) [TTF](2)(2+) cationic(+)···cationic(+) bonds, as well as long, homo-multicenter 2e(-)/4c [TCNE](2)(2-) anionic(-)···anionic(-) bonding. The MO diagrams for the α, β, and δ polymorphs have all of the features found for conventional covalent C-C bonds, and for all of the previously studied multicenter long bonds, for example, π-[TTF](2)(2+) and π-[TCNE](2)(2-). The HOMOs for α-, β-, and δ-[TTF][TCNE] have 2c C-S and 3c C-C-C orbital-overlap contributions between the [TTF](δ+-) and [TCNE](δ-) moieties; these are the shortest intra [TTF···TCNE] separations. Thus, from an orbital-overlap perspective, the bonding has 2c and 3c components residing over one S and four C atoms.

Structures of protonated benzene dimer and intermolecular interaction decomposition via symmetry-adapted perturbation theory

Intermolecular interactions of the protonated benzene dimer are examined through computed interaction energies of three distinct dimer configurations and the decomposition of the interaction into electrostatic, induction, exchange-repulsion, and dispersion components. Protonated benzene dimer is an effective prototype for the study of positively charged, non-covalent interactions commonly encountered throughout aromatic regions of complex biological systems. Three distinct configurations of the dimer are identified: parallel-displaced, T-shaped, and canted. Equilibrium interaction energies computed with CCSD(T) range from −11.4 to −10.0 kcal mol −1 where the lowest energy structures involve direct interaction of the protonation site with the top face of benzene. All three conformers contain unorthodox, cationic, CH–π hydrogen bonds. Comparison of the parallel-displaced structures of C 6 H 6 – C 6 H 7 + and ( C 6 H 6 ) 2 + reveals that the interaction between benzene and protonated benzene is weaker than the interaction between benzene and benzene radical cation. Symmetry adapted perturbation theory provides evidence that the benzene–benzenium attractive interaction consists primarily of equal contributions from dispersion and electrostatics.

Interaction energies of histidine with cations (H+, Li+, Na+, K+, Mg2+, Ca2+)

The imidazol side group of histidine has two nitrogen atoms capable of being protonated or participating in metal binding. Hence, histidine can take on various metal-bound and protonated forms in proteins. Because of its variable structural state, histidine often functions as a key amino acid residue in enzymatic reactions. Ab initio (HF and MP2) calculations were done in modeling the cation (H+, Li+, Na+, K+, Mg2+, Ca2+) interaction with side chain of histidine. The region selectivity of metal ion complexation is controlled by the affinity of the side of attack. In the imidazol unite of histidine the ring nitrogen has much higher metal ion (as well as proton) affinity. The complexation energies with the model systems decrease in the following order: Mg2+ > Ca2+ > Li+ > Na+ > K+. The variation of the bond lengths and the extent of charge transfer upon complexation correlate well with the computed interaction energies.

Effect of Solvation on Ion Binding to Imidazole and Methylimidazole

Quantum chemical [MP2(FULL)/6-311++G-(d,p)] calculations are done on the binding of hydrated Li(+), Na(+), K(+), Mg(2+), Cu(+), and Zn(2+) metal ions with biologically relevant heteroaromatics such as imidazole and methylimidazole. The computed interaction energies are found to be in good agreement with the available experimental data. The effect of hydration on hydrogen bonding has been studied in detail and it shows that the hydrogen bond strength between H(2)O···H-N(1) substantially increases in the presence of metal ions. The present study quantifies the cooperativity between M···imidazole (M = Li(+), Na(+), K(+), Mg(2+), Cu(+), and Zn(2+)) and N(1)-H···OH(2) interactions. Topological atoms in molecules (AIM) analysis and charge analysis support the variation in hydrogen-bonding strength and the variation in M···imidazole binding strength. Effect of hydration on N(1)-H stretching frequency is studied, and it shows a clear shift in the stretching frequency after sequential hydration of metal ion as well as the N(1) of imidazole. The present study provides a detailed account on the biologically important M-histidine motif interaction with metal ions, where histidine is modeled by imidazole and methylimidazole.

Formal Estimation of Errors in Computed Absolute Interaction Energies of Protein-ligand Complexes.

A largely unsolved problem in computational biochemistry is the accurate prediction of binding affinities of small ligands to protein receptors. We present a detailed analysis of the systematic and random errors present in computational methods through the use of error probability density functions, specifically for computed interaction energies between chemical fragments comprising a protein-ligand complex. An HIV-II protease crystal structure with a bound ligand (indinavir) was chosen as a model protein-ligand complex. The complex was decomposed into twenty-one (21) interacting fragment pairs, which were studied using a number of computational methods. The chemically accurate complete basis set coupled cluster theory (CCSD(T)/CBS) interaction energies were used as reference values to generate our error estimates. In our analysis we observed significant systematic and random errors in most methods, which was surprising especially for parameterized classical and semiempirical quantum mechanical calculations. After propagating these fragment-based error estimates over the entire protein-ligand complex, our total error estimates for many methods are large compared to the experimentally determined free energy of binding. Thus, we conclude that statistical error analysis is a necessary addition to any scoring function attempting to produce reliable binding affinity predictions.

Optimal π-Stacking Interaction Energies in Parallel-Displaced Aryl/Aryl Dimers are Predicted by the Dimer Heavy Atom Count

It is generally accepted that large basis sets are required for the accurate calculation of interaction energies for weakly bound aryl dimers, but it has recently been reported that MP2(full)/6-31G* energies, although inaccurate in absolute terms, are well-correlated with estimated CCSD(T)/CBS results. It is now shown that this correlation holds for MP2/aug-cc-pVDZ and SCS-MP2/aug-cc-pVDZ values. Linear regression of published CCSD(T)/CBS results with MP2 or SCS-MP2 results has been used to correct systematic errors observed with both MP2 theories, and these corrections are applied to 27 parallel-displaced aromatic dimers of interest in medicinal chemistry. The optimal computed interaction energies are found to be strongly correlated with the heavy atom counts of the aryl/aryl dimers. This relationship between heavy atom count and interaction energy also applies to a series of 12 aryl/non-aryl dimers such that a single linear regression equation accounts for all of the dimers studied.

Antitumor activity of bent metallocenes: electronic structure analysis using DFT computations

The antitumor activities of bent metallocenes [Cp-M-Cp](2+) (M = Ti, V, Nb, Mo) and complexes of them with guanine, adenine, thymine and cytosine nucleotides have been probed using electronic structure calculations. DFT/BP86 calculations have revealed that the bent metallocene-nucleotide interaction strongly depends on the stability of the hydrolyzed form of the bent metallocene dichloride [Cp(2)M](2+) species, and in turn the stability of the [Cp(2)M](2+) species strongly depends on the electronic structure of [Cp(2)M](2+). Detailed electronic structure and Walsh energy analyses have been carried out for the hydrolyzed forms of four [Cp-M-Cp](2+) (M = Ti, V, Nb, Mo) species to find out why the bent structure is unusually stable. Energy changes that occur during the bending process in frontier molecular orbitals as well as the p(π)-d(π) overlap have been invoked to account for the anticipated antitumor activities of these species. The bonding situation and the interactions in bent metallocene-nucleotide adducts were elucidated by fragment analysis. Of the four nucleotides complexed with the four bent metallocenes, adenine and guanine show better binding abilities than the other two nucleotides. Metallocenes of second-row transition metals exhibit better binding with pyrimidine-base nucleotides. In particular, the Lewis acidic bent metallocenes interact strongly with nucleotides. The antitumor activity is directly related to the binding strength of the bent metallocene with nucleotide adducts, and the computed interaction energy values correlate very well with the experimentally observed antitumor activities.

Natural Lignans as Adhesives for Cellulose: Computational Interaction Energy vs Experimental Results

Comparison between the molecular mechanics calculated energy of interaction of lignan dimers and trimers with a cellulose I crystallite and the experimental values of Young's modulus obtained by thermomechanical analysis (TMA) of cellulose paper impregnated with low molecular weight lignins showed good correlation between calculated and experimental results. The oligomer composition of the four low molecular weight lignins tested was obtained by MALDI-TOF mass spectrometry. This showed that these lignins were predominantly composed of dimers and trimers rendering them ideal for correlation testing. The lignan/cellulose crystallite interaction energy is determined by the oligomer molecular weight as well as the type of linkage within the lignan oligomers. Lignans with higher molecular weight in which the units are linked as β–O–4 give interaction energy values indicating stronger attraction with cellulose.

Explicit polarization (X-Pol) potential using ab initio molecular orbital theory and density functional theory.

The explicit polarization (X-Pol) method has been examined using ab initio molecular orbital theory and density functional theory. The X-Pol potential was designed to provide a novel theoretical framework for developing next-generation force fields for biomolecular simulations. Importantly, the X-Pol potential is a general method, which can be employed with any level of electronic structure theory. The present study illustrates the implementation of the X-Pol method using ab initio Hartree-Fock theory and hybrid density functional theory. The computational results are illustrated by considering a set of bimolecular complexes of small organic molecules and ions with water. The computed interaction energies and hydrogen bond geometries are in good accord with CCSD(T) calculations and B3LYP/aug-cc-pVDZ optimizations.

A simple model of hydrogen bonding with particular application to trends in hydrogen-bonded dimers

A simple model has been proposed to explain trends in the computed interaction energy, bond length changes, frequency shifts and infrared intensities for the chlorofluoromethanes CFnClmH, FH and FArH on complexation with the isoelectronic diatomics BF, CO, N2 and the rare gas atoms Kr, Ar, Ne to form a series of linear or nearly linear hydrogen-bonded complexes. The dipole moment derivative of the proton donor (with respect to the stretching coordinate) and the chemical hardness of the hydrogen-bonded atom of the proton acceptor are identified as two useful parameters for rationalizing the changes in some of the molecular properties of the proton donor when the hydrogen bond is formed. © 2009 Wiley Periodicals, Inc. Int J Quantum Chem, 2009

Model study on sorption of polycyclic aromatic hydrocarbons to goethite

A systematic DFT study of interactions between a set of mono- and polyaromatic hydrocarbons (PAHs) and the (110) goethite surface have been performed in this work. It was found that PAHs form relatively weak surface complexes having their molecular plane practically parallel to the surface plane. The origin of the interactions is in the polarization of the π-system by polar OH groups and in the formation of weak hydrogen bonds where the π-system acts as a proton acceptor. The computed perpendicular distances of the molecular plane of PAHs to the hydrogen atoms of the surface OH groups range from 2.3 to 2.7 Å. Computed interaction energies regularly increase for the linearly shaped molecules from benzene to anthracene. Two other PAHs with a nonlinear shape, phenanthrene and pyrene, are less strongly bound to the surface although they have a similar (phenanthrene) or even larger size (pyrene) than anthracene. These differences were explained by the specific configuration of the surface hydroxyl groups of goethite. The three types of OH groups, μ-OH, μ 3 -O IIH, and OH, form a valley, the width of it fits very well the molecular shape of the linear PAHs. It was found that with anthracene as example the linear PAHs can easily slide along the valley of OH groups with practically no barrier. In summary it is concluded from our results that the (110) goethite surface will withhold linear PAHs significantly better than nonlinear ones.

DFT study of urea interaction with potassium chloride surfaces

In this work, interaction of urea on the important surfaces of potassium chloride has been investigated using density functional theory (DFT) approach. We have performed the DFT calculations using slab models for the interactions of urea with {100}, {110} and {111} surfaces of KCl. Furthermore, the influence of solvent effect was computed with the COSMO continuum model. While sodium chloride crystals are known to change their habit in presence of urea impurity, the higher homologue KCl has not been extensively studied. In one of our earlier articles, we have shown that the preferential interaction of urea with the specific surface of sodium chloride can lead to change the habit from cubes to octahedrons [A. Singh, S. Chakraborty and B. Ganguly, Computational study of urea and its homologue glycinamide: conformations, rotational barriers, and relative interactions with sodium chloride, Langmuir 23 (2007), pp. 5406–5411]. The similar surface docking approach for the interaction of urea with the KCl surfaces suggests that urea would not influence the morphology of KCl. The computed interaction energies were found to be comparable for all the important surfaces of KCl with urea. We have performed experimental studies to examine our computational analysis. The experimental results suggest that the growing KCl crystals in urea solution are cubic in nature. The combined computational and experimental studies corroborate the model of effective interaction of additives with the certain face of crystals and responsible for the change in the morphology [R. Speidel, Neues Jahrbuch fur Mineralogie, Part 4, 81m (1961); N. Cabrera, D.A. Vermilyea. Growth and Perfection of Crystals; John Wiley & Sons Inc. 393. (1958); C.W. Bunn, Adsorption oriented overgrowth and the mixed crystal formation, Proc. R. Soc. A 141 (1933), pp. 567–593].

Substituent Effects in the Benzene Dimer are Due to Direct Interactions of the Substituents with the Unsubstituted Benzene

The prevailing views of substituent effects in the sandwich configuration of the benzene dimer are flawed. For example, in the polar/pi model of Cozzi and co-workers (J. Am. Chem. Soc. 1992, 114, 5729), electron-withdrawing substituents enhance binding in the benzene dimer by withdrawing electron density from the pi-cloud of the substituted ring, reducing the repulsive electrostatic interaction with the nonsubstituted benzene. Conversely, electron-donating substituents donate excess electrons into the pi-system and diminish the pi-stacking interaction. We present computed interaction energies for the sandwich configuration of the benzene dimer and 24 substituted dimers, as well as sandwich complexes of substituted benzenes with perfluorobenzene. While the computed interaction energies correlate well with sigmam values for the substituents, interaction energies for related model systems demonstrate that this trend is independent of the substituted ring. Instead, the observed trends are consistent with direct electrostatic and dispersive interactions of the substituents with the unsubstituted ring.

Experimental and Computational Study of HXeY···HX Complexes (X, Y = Cl and Br): An Example of Exceptionally Large Complexation Effect

The complexes of xenon hydrides HXeY (Y = Cl and Br) with hydrogen halides HX (X = Cl and Br) have been studied both computationally and experimentally in a xenon matrix. The experiments revealed three new complexes: HXeBr...HBr, HXeBr...HCl, and HXeCl...HCl. The experimental assignments were done on the basis of the strong H-Xe stretching absorption of HXeY (Y = Cl and Br) molecules and supported by theoretical results. We experimentally obtained monomer-to-complex blue-shifts of this vibrational mode for all the studied systems (up to approximately 150 cm (-1)). The electronic structure calculations revealed three local structures for each HNgY...HX complexes and their computed interaction energies varied between -460 and -2800 cm (-1). The computational estimates of the vibrational shifts were in agreement with the experimental values. We also found possible experimental absorption belonging to HXeBr...(HBr) 2 trimer and its vibrational shift (+245 cm (-1)) is similar to the computational estimate of a cyclic ternary complex (+252 cm (-1)).

An accurate three-dimensional potential energy surface for the He-Na2 complex

An accurate three-dimensional potential energy surface (PES) for the He-Na2 van der Waals complex was calculated at the coupled cluster singles-and-doubles with noniterative inclusion of connected triple (CCSD(T)) level of theory. A mixed basis set, aug-cc-pVQZ for the He atom and cc-pCVQZ for the sodium atom, and an additional (3s3p2d1f) set of midbond functions were used. The computed interaction energies in 819 configurations were fitted to a 96-parameter analytic potential model by least squares fitting. The PES has two shallow wells corresponding to the T-shaped structure and the linear configuration, which are located at 12.5a0 and 14 a0 with depths of 1.769 and 1.684 cm−1, respectively. The who potential energy surface exhibits weak anisotropy. Based on the fitted PES, state-to-state differential cross sections were calculated.

Interaction of CO with Kr: Potential energy surface and bound states

The first ab initio potential energy surface of the Kr-CO complex is developed using single and double excitation coupled-cluster theory with noniterative treatment of triple excitations. Mixed basis sets, aug-cc-pVQZ for the C and O atoms and aug-cc-pVQZ-PP for the Kr atom, with an additional (3s3p2d2f1g) set of midbond functions are used. The computed interaction energies in 336 configurations are analytically fitted to a two-dimensional potential model by a least squares fit. The potential has a minimum of -119.68 cm(-1) with Re=7.35a 0 at an approximate T-shaped geometry (theta e=98.5 degrees ). Bound state energies are calculated up to J=12, thus enabling a comprehensive comparison between theory and available experimental data as well as providing detailed guidance for future spectroscopic investigations of higher-lying states. The predicted transition frequencies and spectroscopic constants are in good agreement with the experimental results.