Mixed implicit/explicit solvation models in quantum mechanical calculations of binding enthalpy for protein–ligand complexes

An approach to quantum mechanical investigation of interactions in protein–ligand complexes has been developed that treats the solvation effect in a mixed scheme combining implicit and explicit solvent models. In this approach, the first solvation shell of the solvent around the solute is modeled with a limited number of hydrogen bonded explicit solvent molecules. The influence of the remaining bulk solvent is treated as a surrounding continuum in the conductor‐like screening model (COSMO). The enthalpy term of the binding free energy for the protein–ligand complexes was calculated using the semiempirical PM3 method implemented in the MOPAC package, applied to a trimmed model of the protein–ligand complex constructed with special rules. The dependence of the accuracy of binding enthalpy calculations on size of the trimmed model and number of optimized parameters was evaluated. Testing of the approach was performed for 12 complexes of different ligands with trypsin, thrombin, and ribonuclease with experimentally known binding enthalpies. The root‐mean‐square deviation (RMSD) of the calculated binding enthalpies from experimental data was found as ～1 kcal/mol over a large range. © 2006 Wiley Periodicals, Inc. Int J Quantum Chem, 2006     

TANGO: A high through-put conformation generation and semiempirical method-based optimization tool for ligand molecules

Lead optimization is one of the crucial steps in the drug discovery pipeline. After identifying the lead molecule and obtaining its 2D geometry, understanding the best conformation it would attain in 3D still remains one of the most challenging steps in drug discovery. There have been multiple methods and algorithms that are directed toward achieving best conformation for the lead molecules. TANGO focuses on conformation generation and its optimization using semiempirical energy calculations. The conformation generation is based on torsion angle rotation of the exocyclic bonds. The energy calculations are performed using MOPAC. The unique feature of this tool lies in the implementation of Message Passing Interface (MPI) for conformation generation and semiempirical-based optimization. A well-defined architecture handling the input and output generation has been used. The master and slave approach to handle operations involved in torsion angle rotation and energy calculations has helped in load balancing the process of conformation generation. The benchmarking results suggest that TANGO scales significantly well across eight nodes with each node utilizing 16 cores. This tool may prove to very useful in high throughput generation of semiempirically optimized small molecule conformations. The use of semiempirical methods for optimization generates a conformational ensemble thereby helping to obtain stable and alternate stable conformers for a given ligand molecule. © 2018 Wiley Periodicals, Inc.     

Cyclophosphamides as hypoxia-activated diffusible cytotoxins: A theoretical study

Cyclophosphamides have been in clinical use as anti-cancer drugs for a long time and much research has been directed towards reducing their side effects. Here we have performed a theoretical investigation into the possibility of designing bioreductive analogues of cyclophosphamides. Our calculations have employed semiempirical molecular orbital AM1-SM2 and PM3-SM3 calculations, as implemented in MOPAC 93, which include a modified Born method for the treatment of solvation. We have investigated the effect of bioreductive activation on the -elimination reaction that is central to the activation of cyclophosphamides. The approach was tested on two known bioreductive agents, including CB1954, and gave results in agreement with experiment. Non-local density functional calculations on CB1954 and its metabolites, including the radical anion, were in agreement with the semiempirical calculations. The calculations have identified a number of potentially novel bioreductive cyclophosphamides. In particular, our calculations identified compounds in which the initial one-electron reduction was not activating. Such compounds are likely to be more effective bioreductive agents, as the -elimination will not compete under oxic conditions with the important re-oxidation required for the protection of oxic tissue.     

Calculations of the enthalpies of formation of oxygen-containing adamantane derivatives by quantum-chemical methods

The enthalpies of formation of oxygen-containing adamantane derivatives were calculated using the semiempirical quantum-chemical PM3, MINDO, AM1, and MNDO methods implemented in the MOPAC package. The calculated and experimental values were compared. The best correlation was obtained for AM1 calculations. This method was therefore used to calculate the enthalpies of formation of 20 oxygen-containing adamantane derivatives.     

QMQSAR: utilization of a semiempirical probe potential in a field-based QSAR method

A semiempirical quantum mechanical approach is described for the creation of molecular field-based QSAR models from a set of aligned ligand structures. Each ligand is characterized by a set of probe interaction energy (PIE) values computed at various grid points located near the surface of the ligand. Single-point PM3 calculations afford these PIE values, which represents a pool of independent variables from which multilinear regression models of activity are built. The best n-variable fit is determined by constructing an initial regression using standard forward stepwise selection, followed by refinement using a simulated annealing technique. The resulting fit provides an easily interpreted 3D physical model of ligand binding affinity. Validation against three literature datasets demonstrates the ability of the semiempirical potential to model critical binding interactions in diverse systems.     

Pose scoring by NMR

Recently, we have developed a fast approach to calculate NMR chemical shifts using the divide and conquer method at the semiempirical level. To demonstrate the utility of this approach for characterizing protein-ligand interactions, we used the deviation of calculated chemical shift perturbations from experiment to determine the orientation of a ligand (GPI-1046) in the binding pocket of the FK506 binding protein (FKBP12). Moreover, we were able to select the native state of the ligand from a collection of decoy poses. A key hydrogen bond between O1 and HN in Ile56 was also identified. Our results suggest that ligand-induced chemical shift perturbations can be used to refine protein/ligand structures.     

Asymmetric synthesis of alpha,alpha-disubstituted alpha-amino acids by diastereoselective alkylation of camphor-based tricyclic iminolactone

A novel and convenient route for the preparation of chiral tricyclic iminolactones 9 and 10 from camphorquinone has been developed. Alkylation of iminolactones 9 and 10 provided iminolactones 16 and 17 in high yields which were, in turn, alkylated again to afford the alpha,alpha-disubstituted products in good yields (70-90%) and excellent diastereoselectivities (>98%). Hydrolysis of the alkylated iminolactones furnished the desired alpha,alpha-disubstituted alpha-amino acids in good yields and high enantiomeric excesses with good recovery yields of the chiral auxiliary 12 and 13. The extremely high endo-face selectivity for alkylation is discussed using semiempirical (MOPAC 93) calculations.     

Enthalpies of formation of nitro-, nitroxy-, and nitrosoadamantanes

The enthalpies of formation of nitro-, nitroso-, and nitroxyadamantanes were calculated using semiempirical quantum-chemical PM3, MINDO, AM1, and MNDO methods incorporated into MOPAC software package. The best correlation with the experimental data was obtained with the results of PM3 calculations. Using the corresponding linear regression equation, the enthalpies of formation of 22 adamantane derivatives having nitro-, nitroso-, and nitroxy groups were calculated.     

XAFS study of the complex of an acetylacetonate-based ligand and copper ion

An X-ray absorption fine structure (XAFS) study has been conducted to reveal the local structure and chemical state of the copper in the complex of an acetylacetonate-based ligand (L1) and copper ion in acetonitrile solution. The copper ion in the complex was found to be divalent from the Cu K-edge X-ray absorption near-edge structure (XANES) spectrum. The FEFF (ab initio multiple scattering calculations of XAFS) were performed with the model compounds, whose structures were optimized by using MOPAC program with AM1 Hamiltonian. The comparative study of the experimental XAFS spectra and theoretical calculations from FEFF gave the perspectives for clarifying the coordination structure of the complex of L1 and copper ion.     

Large scale free energy calculations for blind predictions of protein–ligand binding: the D3R Grand Challenge 2015

We describe binding free energy calculations in the D3R Grand Challenge 2015 for blind prediction of the binding affinities of 180 ligands to Hsp90. The present D3R challenge was built around experimental datasets involving Heat shock protein (Hsp) 90, an ATP-dependent molecular chaperone which is an important anticancer drug target. The Hsp90 ATP binding site is known to be a challenging target for accurate calculations of ligand binding affinities because of the ligand-dependent conformational changes in the binding site, the presence of ordered waters and the broad chemical diversity of ligands that can bind at this site. Our primary focus here is to distinguish binders from nonbinders. Large scale absolute binding free energy calculations that cover over 3000 protein–ligand complexes were performed using the BEDAM method starting from docked structures generated by Glide docking. Although the ligand dataset in this study resembles an intermediate to late stage lead optimization project while the BEDAM method is mainly developed for early stage virtual screening of hit molecules, the BEDAM binding free energy scoring has resulted in a moderate enrichment of ligand screening against this challenging drug target. Results show that, using a statistical mechanics based free energy method like BEDAM starting from docked poses offers better enrichment than classical docking scoring functions and rescoring methods like Prime MM-GBSA for the Hsp90 data set in this blind challenge. Importantly, among the three methods tested here, only the mean value of the BEDAM binding free energy scores is able to separate the large group of binders from the small group of nonbinders with a gap of 2.4 kcal/mol. None of the three methods that we have tested provided accurate ranking of the affinities of the 147 active compounds. We discuss the possible sources of errors in the binding free energy calculations. The study suggests that BEDAM can be used strategically to discriminate binders from nonbinders in virtual screening and to more accurately predict the ligand binding modes prior to the more computationally expensive FEP calculations of binding affinity.     

A molecular dynamics investigation of CDK8/CycC and ligand binding: conformational flexibility and implication in drug discovery

Abnormal activity of cyclin-dependent kinase 8 (CDK8) along with its partner protein cyclin C (CycC) is a common feature of many diseases including colorectal cancer. Using molecular dynamics (MD) simulations, this study determined the dynamics of the CDK8-CycC system and we obtained detailed breakdowns of binding energy contributions for four type-I and five type-II CDK8 inhibitors. We revealed system motions and conformational changes that will affect ligand binding, confirmed the essentialness of CycC for inclusion in future computational studies, and provide guidance in development of CDK8 binders. We employed unbiased all-atom MD simulations for 500 ns on twelve CDK8-CycC systems, including apoproteins and protein–ligand complexes, then performed principal component analysis (PCA) and measured the RMSF of key regions to identify protein dynamics. Binding pocket volume analysis identified conformational changes that accompany ligand binding. Next, H-bond analysis, residue-wise interaction calculations, and MM/PBSA were performed to characterize protein–ligand interactions and find the binding energy. We discovered that CycC is vital for maintaining a proper conformation of CDK8 to facilitate ligand binding and that the system exhibits motion that should be carefully considered in future computational work. Surprisingly, we found that motion of the activation loop did not affect ligand binding. Type-I and type-II ligand binding is driven by van der Waals interactions, but electrostatic energy and entropic penalties affect type-II binding as well. Binding of both ligand types affects protein flexibility. Based on this we provide suggestions for development of tighter-binding CDK8 inhibitors and offer insight that can aid future computational studies.     

Determination of the thermodynamic functions of lead carbonate on the basis of the calculation of vibrational states

The force constants of bond angles and bonds and parameters of the interatomic potential for the natural carbonate cerussite were determined from the valence force field calculation of the vibrational states of its crystal structure. The initial force constants were calculated by the semiempirical PM5 method using the MOPAC quantum-chemical program package. As the criterion of adequacy of calculations, the consistency between the simulated IR and Raman spectra and the experimental spectra of the compound was used. The heat capacity of lead carbonate as a function of temperature was calculated based on the theory of crystal lattice dynamics and by quantum-chemical methods. The best fit to the experimental data was provided by the semiempirical PM5 method. From the calculated heat capacities, the entropy values of the compound were found.     

X-ray diffraction and quantum-chemical analysis of a single crystal of 2,5-dimethyl-3,4-dihydro-2h-pyran-2-carboxylic acid

The structure of synthesized 2,5-dimethyl-3,4-dihydro-2H-pyran-2-carboxylic acid was investigated by the single crystal X-ray diffraction analysis method. It was established that the molecule of the acid exists in the form of the endo isomer while the single crystal exists as a racemate of the two enantiomeric endo stereomers. Quantum-chemical calculations of a model of the macrocell of the acid by means of the semiempirical MOPAC2009 program agree well with the X-ray diffraction data.     

Semiempirical configuration interaction calculations for ru‐centered dyes*

Computational investigation of the photochemical properties of transition‐metal‐centered dyes typically involves optimization of the molecular structure followed by calculation of the UV/visible spectrum. At present, these steps are usually carried out using density functional theory (DFT) and time‐dependent DFT calculations. Recently, we demonstrated that semiempirical methods with appropriate parameterization could yield geometries that were in very good agreement with DFT calculations, allowing large sets of molecules to be screened quickly and efficiently. In this article, we modify a configuration interaction (CI) method based on a semiempirical PM6 Hamiltonian to determine the UV/visible absorption spectra of Ru‐centered complexes. Our modification to the CI method is based on a scaling of the two‐center, two‐electron Coulomb integrals. This modified, PM6‐based method shows a significantly better match to the experimental absorption spectra versus the default configuration interaction method (in MOPAC) on a training set of 13 molecules. In particular, the modified PM6 method blue‐shifts the location of the metal‐to‐ligand charge‐transfer (MLCT) peaks, in better agreement with experimental and DFT‐based computational results, correcting a significant deficiency of the unmodified method. Published 2018. This article is a U.S. Government work and is in the public domain in the USA     

Spin–orbit coupling in the intersystem crossing of the ring-opened oxirane biradical

The intersystem crossing (ISC ) between the lowest triplet and singlet states occurring in the reaction of atomic oxygen with ethylene was studied. The importance of spin–orbit coupling (SOC ) in oxirane biradicals (?R′R″—CRR*—?) is stressed through calculations where the spin–orbit matrix elements over the full Breit–Pauli SOC operator has been obtained in the singlet–triplet crossing region. The calculations are performed with a multiconfigurational linear response approach, in which the spin–orbit couplings are obtained from triplet response functions using differently correlated singlet-reference-state wave functions. Computational results confirm earlier semiempirical predictions of the spin–orbit coupling as an important mechanism behind the ring opening of oxiranes and addition of oxygen O(³P) atoms to alkenes. © 1995 John Wiley & Sons, Inc.     

How to obtain statistically converged MM/GBSA results

The molecular mechanics/generalized Born surface area (MM/GBSA) method has been investigated with the aim of achieving a statistical precision of 1 kJ/mol for the results. We studied the binding of seven biotin analogues to avidin, taking advantage of the fact that the protein is a tetramer with four independent binding sites, which should give the same estimated binding affinities. We show that it is not enough to use a single long simulation (10 ns), because the standard error of such a calculation underestimates the difference between the four binding sites. Instead, it is better to run several independent simulations and average the results. With such an approach, we obtain the same results for the four binding sites, and any desired precision can be obtained by running a proper number of simulations. We discuss how the simulations should be performed to optimize the use of computer time. The correlation time between the MM/GBSA energies is ～5 ps and an equilibration time of 100 ps is needed. For MM/GBSA, we recommend a sampling time of 20–200 ps for each separate simulation, depending on the protein. With 200 ps production time, 5–50 separate simulations are required to reach a statistical precision of 1 kJ/mol (800–8000 energy calculations or 1.5–15 ns total simulation time per ligand) for the seven avidin ligands. This is an order of magnitude more than what is normally used, but such a number of simulations is needed to obtain statistically valid results for the MM/GBSA method. © 2009 Wiley Periodicals, Inc. J Comput Chem 2010     

Theoretical study of aminoalkylation in the Mannich reaction of furan with methyleneimminium salt

The potential energy surface for the reaction of furan and methyleneimminium cation with formation of a Mannich base has been studied using AM1 and PM3 semiempirical calculations. Nonspecific solvent effects were taken account of in the framework of the multicavity self-consistent reaction field approach. Characteristics of the reaction path elucidated for various media are discussed. © 1998 John Wiley & Sons, Inc. Int J Quant Chem 67: 359–366, 1998     

Binding affinities by alchemical perturbation using QM/MM with a large QM system and polarizable MM model

The most general way to improve the accuracy of binding‐affinity calculations for protein–ligand systems is to use quantum‐mechanical (QM) methods together with rigorous alchemical‐perturbation (AP) methods. We explore this approach by calculating the relative binding free energy of two synthetic disaccharides binding to galectin‐3 at a reasonably high QM level (dispersion‐corrected density functional theory with a triple‐zeta basis set) and with a sufficiently large QM system to include all short‐range interactions with the ligand (744–748 atoms). The rest of the protein is treated as a collection of atomic multipoles (up to quadrupoles) and polarizabilities. Several methods for evaluating the binding free energy from the 3600 QM calculations are investigated in terms of stability and accuracy. In particular, methods using QM calculations only at the endpoints of the transformation are compared with the recently proposed non‐Boltzmann Bennett acceptance ratio (NBB) method that uses QM calculations at several stages of the transformation. Unfortunately, none of the rigorous approaches give sufficient statistical precision. However, a novel approximate method, involving the direct use of QM energies in the Bennett acceptance ratio method, gives similar results as NBB but with better precision, ～3 kJ/mol. The statistical error can be further reduced by performing a greater number of QM calculations. © 2015 Wiley Periodicals, Inc.