Protein‐specific force field derived from the fragment molecular orbital method can improve protein–ligand binding interactions

Accurate computational estimate of the protein–ligand binding affinity is of central importance in rational drug design. To improve accuracy of the molecular mechanics (MM) force field (FF) for protein–ligand simulations, we use a protein‐specific FF derived by the fragment molecular orbital (FMO) method and by the restrained electrostatic potential (RESP) method. Applying this FMO‐RESP method to two proteins, dodecin, and lysozyme, we found that protein‐specific partial charges tend to differ more significantly from the standard AMBER charges for isolated charged atoms. We did not see the dependence of partial charges on the secondary structure. Computing the binding affinities of dodecin with five ligands by MM PBSA protocol with the FMO‐RESP charge set as well as with the standard AMBER charges, we found that the former gives better correlation with experimental affinities than the latter. While, for lysozyme with five ligands, both charge sets gave similar and relatively accurate estimates of binding affinities. © 2013 Wiley Periodicals, Inc.     

Fragment‐based quantum mechanical calculation of protein–protein binding affinities

The electrostatically embedded generalized molecular fractionation with conjugate caps (EE‐GMFCC) method has been successfully utilized for efficient linear‐scaling quantum mechanical (QM) calculation of protein energies. In this work, we applied the EE‐GMFCC method for calculation of binding affinity of Endonuclease colicin–immunity protein complex. The binding free energy changes between the wild‐type and mutants of the complex calculated by EE‐GMFCC are in good agreement with experimental results. The correlation coefficient (R) between the predicted binding energy changes and experimental values is 0.906 at the B3LYP/6‐31G*‐D level, based on the snapshot whose binding affinity is closest to the average result from the molecular mechanics/Poisson–Boltzmann surface area (MM/PBSA) calculation. The inclusion of the QM effects is important for accurate prediction of protein–protein binding affinities. Moreover, the self‐consistent calculation of PB solvation energy is required for accurate calculations of protein–protein binding free energies. This study demonstrates that the EE‐GMFCC method is capable of providing reliable prediction of relative binding affinities for protein–protein complexes. © 2018 Wiley Periodicals, Inc.     

Importance of CH/π hydrogen bonds in recognition of the core motif in proline-recognition domains: an ab initio fragment molecular orbital study

We examined CH/π hydrogen bonds in protein/ligand complexes involving at least one proline residue using the ab initio fragment molecular orbital (FMO) method and the program CHPI. FMO calculations were carried out at the Hartree–Fock (HF)/6‐31G*, HF/6‐31G**, second‐order Møller–Plesset perturbation (MP2)/6‐31G*, and MP2/6‐31G** levels for three Src homology 3 (SH3) domains and five proline‐recognition domains (PRDs) complexed with their corresponding ligand peptides. PRDs use a conserved set of aromatic residues to recognize proline‐rich sequences of specific ligands. Many CH/π hydrogen bonds were identified in these complexes. CH/π hydrogen bonds occurred, in particular, in the central part of the proline‐rich motifs. Our results suggest that CH/π hydrogen bonds are important in the recognition of SH3 and PRDs by their ligand peptides and play a vital role in the signal transduction system. Combined use of the FMO method and CHPI analysis is a valuable tool for the study of protein/protein and protein/ligand interactions and may be useful in rational drug design. © 2011 Wiley Periodicals, Inc. J Comput Chem 2011     

Accuracy of the three-body fragment molecular orbital method applied to Møller-Plesset perturbation theory

The three‐body energy expansion in the fragment molecular orbital method (FMO) was applied to the 2nd order Møller–Plesset theory (MP2). The accuracy of both the two and three‐body expansions was determined for water clusters, alanine n‐mers (α‐helices and β‐strands) and one synthetic protein, using the 6‐31G* and 6‐311G* basis sets. At the best level of theory (three‐body, two molecules/residues per fragment), the absolute errors in energy relative to ab initio MP2 were at most 1.2 and 5.0 mhartree, for the 6‐31G* and 6‐311G* basis sets, respectively. The relative accuracy was at worst 99.996% and 99.96%, for 6‐31G* and 6‐311G*, respectively. A three‐body approximation was introduced and the optimum threshold value was determined. The protein calculation (6‐31G*) at the production level (FMO2/2) took 3 h on 36 3.2‐GHz Pentium 4 nodes and had the absolute error in the MP2 correlation energy of only 2 kcal/mol. © 2007 Wiley Periodicals, Inc. J Comput Chem, 2007     

Theoretical studies on opioid receptors and ligands. I. Molecular modeling and QSAR studies on the interaction mechanism of fentanyl analogs binding to μ‐opioid receptor

Based on our previous result of the three‐dimensional model of the μ‐opioid receptor, binding conformations of 13 fentanyl analogs and three‐dimensional structures for the complexs of these analogs with μ‐opioid receptor were constructed employing the molecular modeling method and our binding conformation search program for ligands (BCSPL). Energetic calculation and quantitative structure–activity relationship (QSAR) analysis indicated a good correlation between the calculated binding energies of fentanyl analogs and their binding affinities, pK_i's and pK's, and analgesic activities, − log ED₅₀'s. Based on the three‐dimensional models, the possible interaction mechanism of fentanyl analogs with μ‐opioid receptor can be illustrated and the available structure–activity relationship of these analgesic agents can be explained reasonably. © 2000 John Wiley & Sons, Inc. Int J Quant Chem 78: 285–293, 2000     

Hydrogen bonding interaction in sarcosine–water complex using ab initio and DFT method

The hydrogen bonding interaction in the Sarcosine (N‐methylglycine)–water complex is studied using ab initio, MP2, and density functional theory (DFT/B3LYP). For this complex, binding energies, dipole–dipole interactions, chemical hardness, and chemical potential have been calculated. Three different basis sets, viz. 6‐311+G, 6‐311++G, and 6‐311++G*, have been used to optimize the geometries by all three methods. The basis set superposition errors are also calculated, and the corrected binding energies are reported for this complex. © 2004 Wiley Periodicals, Inc. Int J Quantum Chem, 2005     

Calculating protein–ligand binding affinities with MMPBSA: Method and error analysis

Molecular Mechanics Poisson–Boltzmann Surface Area (MMPBSA) methods have become widely adopted in estimating protein–ligand binding affinities due to their efficiency and high correlation with experiment. Here different computational alternatives were investigated to assess their impact to the agreement of MMPBSA calculations with experiment. Seven receptor families with both high‐quality crystal structures and binding affinities were selected. First the performance of nonpolar solvation models was studied and it was found that the modern approach that separately models hydrophobic and dispersion interactions dramatically reduces RMSD's of computed relative binding affinities. The numerical setup of the Poisson–Boltzmann methods was analyzed next. The data shows that the impact of grid spacing to the quality of MMPBSA calculations is small: the numerical error at the grid spacing of 0.5 Å is already small enough to be negligible. The impact of different atomic radius sets and different molecular surface definitions was further analyzed and weak influences were found on the agreement with experiment. The influence of solute dielectric constant was also analyzed: a higher dielectric constant generally improves the overall agreement with experiment, especially for highly charged binding pockets. The data also showed that the converged simulations caused slight reduction in the agreement with experiment. Finally the direction of estimating absolute binding free energies was briefly explored. Upon correction of the binding‐induced rearrangement free energy and the binding entropy lost, the errors in absolute binding affinities were also reduced dramatically when the modern nonpolar solvent model was used, although further developments were apparently necessary to further improve the MMPBSA methods. © 2016 Wiley Periodicals, Inc.     

BSSE‐corrected geometry and harmonic and anharmonic vibrational frequencies of formamide–water and formamide–formamide dimers

The basis set superposition error (BSSE) influence in the geometry structure, interaction energies, and intermolecular harmonic and anharmonic vibrational frequencies of cyclic formamide–formamide and formamide–water dimers have been studied using different basis sets (6‐31G, 6‐31G**, 6‐31++G**, D95V, D95V**, and D95V++**). The a posteriori “counterpoise” (CP) correction scheme has been compared with the a priori “chemical Hamiltonian approach” (CHA) both at the Hartree–Fock (HF) and second‐order Møller–Plesset many‐body perturbation (MP2) levels of theory. The effect of BSSE on geometrical parameters, interaction energies, and intermolecular harmonic vibrational frequencies are discussed and compared with the existing experimental data. As expected, the BSSE‐free CP and CHA interaction energies usually show less deep minima than those obtained from the uncorrected methods at both the HF and MP2 levels. Focusing on the correlated level, the amount of BSSE in the intermolecular interaction energies is much larger than that at the HF level, and this effect is also conserved in the values of the force constants and harmonic vibrational frequencies. All these results clearly indicate the importance of the proper BSSE‐free correlation treatment with the well‐defined basis functions. At the same time, the results show a good agreement between the a priori CHA and a posteriori CP correction scheme; this agreement is remarkable in the case of large and well‐balanced basis sets. The anharmonic frequency correction values also show an important BSSE dependence, especially for hydrogen bond stretching and for low frequencies belonging to the intermolecular normal modes. © 2005 Wiley Periodicals, Inc. Int J Quantum Chem, 2005     

Ab initio quantum chemical calculation of electron density,electrostatic potential,and electric field of biomolecule based on fragment molecular orbital method

Efficient quantum chemical calculations of electrostatic properties, namely, the electron density (EDN), electrostatic potential (ESP), and electric field (EFL), were performed using the fragment molecular orbital (FMO) method. The numerical errors associated with the FMO scheme were examined at the HF, MP2, and RI‐MP2 levels of theory using 4 small peptides. As a result, the FMO errors in the EDN, ESP, and EFL were significantly smaller than the magnitude of the electron correlation effects, which indicated that the FMO method provides sufficiently accurate values of electrostatic properties. In addition, an attempt to reduce the computational effort was proposed by combining the FMO scheme and a point charge approximation. The error due to this approximation was examined using 2 proteins, prion protein and human immunodeficiency virus type 1 protease. As illustrative examples, the ESP values at the molecular surface of these proteins were calculated at the MP2 level of theory.     

Affinity Screening Using Competitive Binding with Fluorine‐19 Hyperpolarized Ligands

Fluorine‐19 NMR and hyperpolarization form a powerful combination for drug screening. Under a competitive equilibrium with a selected fluorinated reporter ligand, the dissociation constant (K_D) of other ligands of interest is measurable using a single‐scan Carr–Purcell–Meiboom–Gill (CPMG) experiment, without the need for a titration. This method is demonstrated by characterizing the binding of three ligands with different affinities for the serine protease trypsin. Monte Carlo simulations show that the highest accuracy is obtained when about one‐half of the bound reporter ligand is displaced in the binding competition. Such conditions can be achieved over a wide range of affinities, allowing for rapid screening of non‐fluorinated compounds when a single fluorinated ligand for the binding pocket of interest is known.     

DFT,CBS‐Q,W1BD and G4MP2 calculation of the proton and electron affinities,gas phase basicities and ionization energies of saturated and unsaturated carboxylic acids (C1–C4)

The proton affinities, gas phase basicities and ionization energies of formic acid, acetic acid, propanoic acid, 2‐propenoic acid, propiolic acid, butanoic acid, 2‐butenoic acid, 3‐botenoic acid, 2‐methyl‐propanoic acid and 2‐methyl‐2‐propenoic acid were calculated using the computational methods including B3LYP/6‐311++G(2df,p), CBS‐Q and G4MP2. Also, the considered properties were calculated using W1BD method only for formic and acetic acids. In addition, the electron affinities of the acids were calculated using B3LYP, CBS‐Q, G4MP2 and G2MP2 methods, separately. The calculations showed that the PA and gas phase basicity increase with the increase in the number of carbon atoms. The calculated Ionization energies of the unsaturated carboxylic acids are less than the corresponding saturated acids, which are in good agreement with the experimental results. © 2013 Wiley Periodicals, Inc.     

Rapid evaluation of the binding energies in hydrogen-bonded amide-thymine and amide-uracil dimers in gas phase

The binding energies and the equilibrium hydrogen bond distances as well as the potential energy curves of 48 hydrogen‐bonded amide–thymine and amide–uracil dimers are evaluated from the analytic potential energy function established in our lab recently. The calculation results show that the potential energy curves obtained from the analytic potential energy function are in good agreement with those obtained from MP2/6‐311+G** calculations by including the BSSE correction. For all the 48 dimers, the analytic potential energy function yields the binding energies of the MP2/6‐311+G** with BSSE correction within the error limits of 0.50 kcal/mol for 46 dimers, only two differences are larger than 0.50 kcal/mol and the largest one is only 0.60 kcal/mol. The analytic potential energy function produces the equilibrium hydrogen bond distances of the MP2/6‐311+G** with BSSE correction within the error limits of 0.050 Å for all the 48 dimers. The analytic potential energy function is further applied to four more complicated hydrogen‐bonded amide–base systems involving amino acid side chain and β‐sheet. The values of the binding energies and equilibrium hydrogen bond distances obtained from the analytic potential energy function are also in good agreement with those obtained from MP2 calculations with the BSSE correction. These results demonstrate that the analytic potential energy function can be used to evaluate the binding energies in hydrogen‐bonded amide–base dimers quickly and accurately. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011     

Toward quantitative estimates of binding affinities for protein–ligand systems involving large inhibitor compounds: A steered molecular dynamics simulation route

Understanding binding mechanisms between enzymes and potential inhibitors and quantifying protein – ligand affinities in terms of binding free energy is of primary importance in drug design studies. In this respect, several approaches based on molecular dynamics simulations, often combined with docking techniques, have been exploited to investigate the physicochemical properties of complexes of pharmaceutical interest. Even if the geometric properties of a modeled protein – ligand complex can be well predicted by computational methods, it is still challenging to rank with chemical accuracy a series of ligand analogues in a consistent way. In this article, we face this issue calculating relative binding free energies of a focal adhesion kinase, an important target for the development of anticancer drugs, with pyrrolopyrimidine‐based ligands having different inhibitory power. To this aim, we employ steered molecular dynamics simulations combined with nonequilibrium work theorems for free energy calculations. This technique proves very powerful when a series of ligand analogues is considered, allowing one to tackle estimation of protein – ligand relative binding free energies in a reasonable time. In our cases, the calculated binding affinities are comparable with those recovered from experiments by exploiting the Michaelis – Menten mechanism with a competitive inhibitor.     

Geometry prediction of bridged zirconocene dichlorides by quantum chemical methods

The ab initio Hartree–Fock theory has been demonstrated to give accurate geometry predictions for bridged zirconocene dichlorides. Equilibrium geometries of crystallographically characterized bridged zirconocene dichlorides were optimized by Hartree–Fock, MP2, BLYP, and B3LYP methods, with basis sets ranging from 3‐21G* to 6‐311G**. Selected geometrical parameters were compared with experimental crystal structures. The least expensive HF/3‐21G* method proved to be notably accurate. The accuracy of HF/3‐21G* was verified by a comprehensive data set of 62 bridged zirconocene dichlorides. Furthermore, experimental corrections were applied to the optimized geometry parameters to eliminate systematic deviations. Corrections resulted in considerably improved accuracy for systematically overestimated metal–ligand distances, with maximum deviation falling from 0.081 to 0.039 Å, and absolute average deviations from 0.048 to 0.012 Å. Ligand–metal–ligand angles were predicted accurately with absolute average deviations of 0.7–1.3°. Zirconium–chlorine distances and chlorine–zirconium–chlorine angles are relatively constant in the studied molecules. Zirconium–cyclopentadienyl distances can be influenced mainly by modifying the ligand structure, whereas cyclopentadienyl–zirconium– cyclopentadienyl angles and cyclopentadienyl–cyclopentadienyl plane angles can be controlled by bridge modifications. The HF/3‐21G* method can be applied for the estimation of steric effects in zirconocene catalyzed polymerization reactions, therefore being suitable for the construction of structure–polymerization property correlations. © 2000 John Wiley & Sons, Inc. J Comput Chem 22: 51–64, 2001     

Rapid and accurate assessment of GPCR–ligand interactions Using the fragment molecular orbital‐based density‐functional tight‐binding method

The reliable and precise evaluation of receptor–ligand interactions and pair‐interaction energy is an essential element of rational drug design. While quantum mechanical (QM) methods have been a promising means by which to achieve this, traditional QM is not applicable for large biological systems due to its high computational cost. Here, the fragment molecular orbital (FMO) method has been used to accelerate QM calculations, and by combining FMO with the density‐functional tight‐binding (DFTB) method we are able to decrease computational cost 1000 times, achieving results in seconds, instead of hours. We have applied FMO‐DFTB to three different GPCR–ligand systems. Our results correlate well with site directed mutagenesis data and findings presented in the published literature, demonstrating that FMO‐DFTB is a rapid and accurate means of GPCR–ligand interactions. © 2017 Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Rapid evaluation of the binding energies between peptide amide and DNA base

The binding energies and the equilibrium hydrogen bond distances as well as the potential energy curves of 20 hydrogen‐bonded amide–base dimers are evaluated from the analytic potential energy function established in our laboratory recently. The analytic potential energy function is used to calculate the N? H···N, N? H···O?C, C? H···N, and C? H···O?C dipole–dipole attractive interaction energies and C?O···O?C, N? H···H? N, and N? H···H? C dipole–dipole repulsive interaction energies in the 20 dimers composed of DNA bases adenine, guanine, cytosine, or thymine and peptide amide. The calculation results show that the potential energy curves obtained from the analytic potential energy function are in good agreement with those obtained from MP2/6‐311+G** calculations by including the basis set superposition error (BSSE) correction. For all the 20 dimers, the analytic potential energy function yields the binding energies of the MP2/6‐311+G** with BSSE correction within the error limits of 0.50 kcal/mol for 19 dimers, only one difference is larger than 0.50 kcal/mol and the difference is only 0.61 kcal/mol. The analytic potential energy function produces the equilibrium hydrogen bond distances of the MP2/6‐311+G** with BSSE correction within the error limits of 0.030 Å for all the 20 dimers. The analytic potential energy function is further applied to four more complicated DNA base‐peptide amide systems involving amino acid side chain and β‐sheet. The values of the binding energies and equilibrium hydrogen bond distances obtained from the analytic potential energy function are also in good agreement with those obtained from MP2 calculations with the BSSE correction. These results demonstrate that the analytic potential energy function can be used to evaluate the binding energies in hydrogen‐bonded peptide amide–DNA base dimers quickly and accurately. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

Fast and accurate predictions of binding free energies using MM‐PBSA and MM‐GBSA

In the drug discovery process, accurate methods of computing the affinity of small molecules with a biological target are strongly needed. This is particularly true for molecular docking and virtual screening methods, which use approximated scoring functions and struggle in estimating binding energies in correlation with experimental values. Among the various methods, MM‐PBSA and MM‐GBSA are emerging as useful and effective approaches. Although these methods are typically applied to large collections of equilibrated structures of protein‐ligand complexes sampled during molecular dynamics in water, the possibility to reliably estimate ligand affinity using a single energy‐minimized structure and implicit solvation models has not been explored in sufficient detail. Herein, we thoroughly investigate this hypothesis by comparing different methods for the generation of protein‐ligand complexes and diverse methods for free energy prediction for their ability to correlate with experimental values. The methods were tested on a series of structurally diverse inhibitors of Plasmodium falciparum DHFR with known binding mode and measured affinities. The results showed that correlations between MM‐PBSA or MM‐GBSA binding free energies with experimental affinities were in most cases excellent. Importantly, we found that correlations obtained with the use of a single protein‐ligand minimized structure and with implicit solvation models were similar to those obtained after averaging over multiple MD snapshots with explicit water molecules, with consequent save of computing time without loss of accuracy. When applied to a virtual screening experiment, such an approach proved to discriminate between true binders and decoy molecules and yielded significantly better enrichment curves. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

DFT‐based QSAR study and molecular design of AHMA derivatives as potent anticancer agents

A quantitative structure–activity relationship (QSAR) of 3‐(9‐acridinylamino)‐5‐hydroxymethylaniline (AHMA) derivatives and their alkylcarbamates as potent anticancer agents has been studied using density functional theory (DFT), molecular mechanics (MM+), and statistical methods. In the best established QSAR equation, the energy (E_NL) of the next lowest unoccupied molecular orbital (NLUMO) and the net charges (Q_FR) of the first atom of the substituent R, as well as the steric parameter (MR₂) of subsituent R₂ are the main independent factors contributing to the anticancer activity of the compounds. A new scheme determining outliers by “leave‐one‐out” (LOO) cross‐validation coefficient (q) was suggested and successfully used. The fitting correlation coefficient (R²) and the “LOO” cross‐validation coefficient (q²) values for the training set of 25 compounds are 0.881 and 0.829, respectively. The predicted activities of 5 compounds in the test set using this QSAR model are in good agreement with their experimental values, indicating that this model has excellent predictive ability. Based on the established QSAR equation, 10 new compounds with rather high anticancer activity much greater than that of 34 compounds have been designed and await experimental verification. © 2006 Wiley Periodicals, Inc. Int J Quantum Chem, 2007