A simple and efficient dispersion correction to the Hartree-Fock theory (3): A comprehensive performance comparison of HF-Dtq with MP2 and DFT-Ds.

Accurate prediction of the intermolecular interaction energy (ΔEbind) has been a challenging and serious problem. Current in silico drug screening demands efficient and accurate evaluation of ΔEbind for ligands and their target proteins. It is desirable that ΔEbind including the dispersion interaction energy (Edisp) is calculated using a post-Hartree-Fock (HF) theory, such as the high-order coupled-cluster one, with a larger basis set. However, it remains computationally too expensive to apply such a one to large molecular systems. As another problem, it is necessary to consider the contribution of the basis set superposition error (BSSE) in calculation of ΔEbind. In Bioorg. Med. Chem. Lett. 2014 and 2015, we proposed simple and efficient corrections of dispersion and BSSE for the HF theory, which is not able to express the dispersion interaction energy correctly. The current Letter, as the final one in the series, aims to verify the HF theory enhanced by the dispersion correction (HF-Dtq) in the light of reproducibility of 'accurate' intermolecular ligand-protein interaction energy values, with comprehensive comparison with the MP2 and recently proposed various DFT-D theories. Taking ΔEbind calculated with the coupled-cluster theory coupled with a complete basis set as a reference, ΔEbind of over a hundred small sized noncovalent complexes as well as real ligand-protein complexes models was systematically examined in terms of accuracy and computational cost. The comprehensive comparison in the current work showed that HF-Dtq is a practical and reliable approach for in silico drug screening and quantitative structure-activity relationships.

Molecular dynamics and density functional studies on the metabolic selectivity of antipsychotic thioridazine by cytochrome P450 2D6: Connection with crystallographic and metabolic results.

CYP2D6, a cytochrome P450 isoform, significantly contributes to the metabolism of many clinically important drugs. Thioridazine (THD) is one of the phenothiazine-type antipsychotics, which exhibit dopamine D2 antagonistic activity. THD shows characteristic metabolic profiles compared to other phenothiazine-type antipsychotics such as chlorpromazine. The sulfur atom attached to the phenothiazine ring is preferentially oxidized mainly by CYP2D6, that is, the 2-sulfoxide is a major metabolite, and interestingly this metabolite shows more potent activity against dopamine D2 receptors than THD. On the other hand, the formation of this metabolite causes many serious problems for its clinical use. Wójcikowski et al. (Drug Metab. Dispos. 2006, 34, 471) reported a kinetic study of THD formed by CYP2D6. Recently, Wang et al. (J. Biol. Chem. 2012, 287, 10834 and J. Biol. Chem. 2015, 290, 5092) revealed the crystallographic structure of THD with CYP2D6. In the current study, the binding and reaction mechanisms at the atomic and electronic levels were computationally examined based on the assumption as to whether or not the different crystallographic binding poses correspond to the different metabolites. The binding and oxidative reaction steps in the whole metabolic process were investigated using molecular dynamics and density functional theory calculations, respectively. The current study demonstrated the essential importance of the orientation of the substrate in the reaction center of CYP2D6 for the metabolic reaction.

Binding interaction of SGLT with sugar and thiosugar by the molecular dynamics simulation.

The human sodium-glucose co-transporter 2 (hSGLT2) is a transporter responsible for reabsorption of glucose in the proximal convoluted tubule of the kidney. hSGLT2 inhibitors, including luseogliflozin, have been developed as drugs for type 2 diabetes mellitus. Only luseogliflozin contains a thiosugar ring in its chemical structure, while other hSGLT2 inhibitors contain glucose rings. Consequently, we focused on the binding interactions of hSGLT2 with sugars and thiosugars. We first revealed that the binding affinities of thiosugars are stronger than those of sugars through molecular dynamics simulations of Vibrio parahaemolyticus, sodium-galactose co-transporter, and human hSGLT2. We then demonstrated that Na(+) dissociates from the protein to the cytoplasmic solution more slowly in the thiosugar system than in the sugar system. These differences between sugars and thiosugars are discussed on the basis of the different binding modes due to the atom at the 5-position of the sugar and thiosugar rings. Finally, as a result of Na(+) dissociation, we suggest that the dissociation of thiosugars is slower than that of sugars.

A simple and efficient dispersion correction to the Hartree-Fock theory (2): Incorporation of a geometrical correction for the basis set superposition error.

One of the most challenging problems in computer-aided drug discovery is the accurate prediction of the binding energy between a ligand and a protein. For accurate estimation of net binding energy ΔEbind in the framework of the Hartree-Fock (HF) theory, it is necessary to estimate two additional energy terms; the dispersion interaction energy (Edisp) and the basis set superposition error (BSSE). We previously reported a simple and efficient dispersion correction, Edisp, to the Hartree-Fock theory (HF-Dtq). In the present study, an approximation procedure for estimating BSSE proposed by Kruse and Grimme, a geometrical counterpoise correction (gCP), was incorporated into HF-Dtq (HF-Dtq-gCP). The relative weights of the Edisp (Dtq) and BSSE (gCP) terms were determined to reproduce ΔEbind calculated with CCSD(T)/CBS or /aug-cc-pVTZ (HF-Dtq-gCP (scaled)). The performance of HF-Dtq-gCP (scaled) was compared with that of B3LYP-D3(BJ)-bCP (dispersion corrected B3LYP with the Boys and Bernadi counterpoise correction (bCP)), by taking ΔEbind (CCSD(T)-bCP) of small non-covalent complexes as 'a golden standard'. As a critical test, HF-Dtq-gCP (scaled)/6-31G(d) and B3LYP-D3(BJ)-bCP/6-31G(d) were applied to the complex model for HIV-1 protease and its potent inhibitor, KNI-10033. The present results demonstrate that HF-Dtq-gCP (scaled) is a useful and powerful remedy for accurately and promptly predicting ΔEbind between a ligand and a protein, albeit it is a simple correction procedure.

Novel pH-dependent regulation of human cytosolic sialidase 2 (NEU2) activities by siastatin B and structural prediction of NEU2/siastatin B complex.

Human cytosolic sialidase (Neuraminidase 2, NEU2) catalyzes the removal of terminal sialic acid residues from glycoconjugates. The effect of siastatin B, known as a sialidase inhibitor, has not been evaluated toward human NEU2 yet. We studied the regulation of NEU2 activity by siastatin B in vitro and predicted the interaction in silico. Inhibitory and stabilizing effects of siastatin B were analyzed in comparison with DANA (2-deoxy-2,3-dehydro-N-acetylneuraminic acid) toward 4-umbelliferyl N-acetylneuraminic acid (4-MU-NANA)- and α2,3-sialyllactose-degrading activities of recombinant NEU2 produced by E. coli GST-fusion gene expression. Siastatin B exhibited to have higher competitive inhibitory activity toward NEU2 than DANA at pH 4.0. We also revealed the stabilizing effect of siastatin B toward NEU2 activity at acidic pH. Docking model was constructed on the basis of the crystal structure of NEU2/DANA complex (PDB code: 1VCU). Molecular docking predicted that electrostatic neutralization of E111 and E218 residues of the active pocket should not prevent siastatin B from binding at pH 4.0. The imino group (1NH) of siastatin B can also interact with D46, neutralized at pH 4.0. Siastatin B was suggested to have higher affinity to the active pocket of NEU2 than DANA, although it has no C7-9 fragment corresponding to that of DANA. We demonstrated here the pH-dependent affinity of siastatin B toward NEU2 to exhibit potent inhibitory and stabilizing activities. Molecular interaction between siastatin B and NEU2 will be utilized to develop specific inhibitors and stabilizers (chemical chaperones) not only for NEU2 but also the other human sialidases, including NEU1, NEU3 and NEU4, based on homology modeling.

A QSAR study on the inhibition mechanism of matrix metalloproteinase-12 by arylsulfone analogs based on molecular orbital calculations.

A binding mechanism between human matrix metalloproteinase-12 (MMP-12) and eight arylsulfone analogs having two types of carboxylic and hydroxamic acids as the most representative zinc binding group is investigated using a quantitative structure-activity relationship (QSAR) analysis based on a linear expression by representative energy terms (LERE). The LERE-QSAR analysis quantitatively reveals that the variation in the observed (experimental) inhibitory potency among the arylsulfone analogs is decisively governed by those in the intrinsic binding and dispersion interaction energies. The results show that the LERE-QSAR analysis not only can excellently reproduce the observed overall free-energy change but also can determine the contributions of representative free-energy changes. An inter-fragment interaction energy difference (IFIED) analysis based on the fragment molecular orbital (FMO) method (FMO-IFIED) leads to the identification of key residues governing the variation in the inhibitory potency as well as to the understanding of the difference between the interactions of the carboxylic and hydroxamic acid zinc binding groups. The current results that have led to the optimization of the inhibitory potency of arylsulfone analogs toward MMP-12 to be used in the treatment of chronic obstructive pulmonary disease may be useful for the development of a new potent MMP-12 inhibitor.

Binding free-energy calculation is a powerful tool for drug optimization: calculation and measurement of binding free energy for 7-azaindole derivatives to glycogen synthase kinase-3ß.

Present computational lead (drug)-optimization is lacking in thermodynamic tactics. To examine whether calculation of binding free-energy change (ΔG) is effective for the lead-optimization process, binding ΔGs of 7-azaindole derivatives to the ATP binding site of glycogen synthase kinase-3ß (GSK-3ß) were calculated. The result was a significant correlation coefficient of r = 0.895 between calculated and observed ΔGs. This indicates that calculated ΔG reflects the inhibitory activities of 7-azaindole derivatives. In addition to quantitative estimation of activity, ΔG calculation characterizes the thermodynamic behavior of 7-azaindole derivatives, providing also useful information for inhibitor optimization on affinity to water molecules.

A simple and efficient dispersion correction to the Hartree-Fock theory.

One of the most challenging problems in computational chemistry and in drug discovery is the accurate prediction of the binding energy between a ligand and a protein receptor. It is well known that the binding energy calculated with the Hartree-Fock molecular orbital theory (HF) lacks the dispersion interaction energy that significantly affects the accuracy of the total binding energy of a large molecular system. We propose a simple and efficient dispersion energy correction to the HF theory (HF-Dtq). The performance of HF-Dtq was compared with those of several recently proposed dispersion corrected density functional theory methods (DFT-Ds) as to the binding energies of 68 small non-covalent complexes. The overall performance of HF-Dtq was found to be nearly equivalent to that of more sophisticated B3LYP-D3. HF-Dtq will thus be a useful and powerful method for accurately predicting the binding energy between a ligand and a protein, albeit it is a simple correction procedure based on HF.

Connecting Classical QSAR and LERE Analyses Using Modern Molecular Calculations, LERE-QSAR (VI): Hydrolysis of Substituted Hippuric Acid Phenyl Esters by Trypsin.

The reaction mechanism of trypsin was studied by applying DFT and ab initio molecular orbital (MO) calculations to complexes of trypsin with a congeneric series of eight para-substituted hippuric acid phenyl esters, for which a previous quantitative structureactivity relationship (QSAR) study revealed nice linearity of Hammett substitution constant σ(-) with logarithmic values of the MichaelisMenten and catalytic rate constants. Based on the LERE procedure, we performed QSAR analyses on each elementary reaction step during the acylation process. The present calculations showed that the rate-determining step during the acylation process is the transition state (TS) between the enzymesubstrate complex (ES) and tetrahedral intermediate (TET), and that the proton transfer occurs from Ser195 to His57, not between His57 and Asp102. The LERE-QSAR analysis statistically suggested that the variation of overall free-energy changes leading to formation of TS is governed mostly by that of activation energies required to form TS from ES. In spite of a very limited number of congeneric ligands in the current work, it is critically essential to clarify and verify physicochemical meanings of a typical QSAR/Chemoinformatics parameter, Hammett σ(-) based on quantum chemical calculations on the proteinligand kinetics; how Hammett σ(-) behaves in terms of proteinligand interaction energies.

Combined QM/MM (ONIOM) and QSAR approach to the study of complex formation of matrix metalloproteinase­9 with a series of biphenylsulfonamides−LERE-QSAR analysis (V).

We previously proposed a novel QSAR (quantitative structure-activity relationship) procedure called LERE (linear expression by representative energy terms)-QSAR involving molecular calculations such as an ab initio fragment molecular orbital ones. In the present work, we applied LERE-QSAR to complex formation of matrix metalloproteinase-9 (MMP-9) with a series of substituted biphenylsulfonamides. The results shows that the overall free-energy change accompanying complex formation is due to predominantly the contribution from the electrostatic interaction with the zinc atom in the active site of MMP-9. Carbonic anhydrase (CA) belongs to the zinc-containing protease family. In contrast to the current case of MMP-9, the overall free-energy change during complex formation of CA with a series of benzenesulfonamides is due to the contributions from the solvation and dissociation free-energy changes, as previously reported. Comparison of the two sets of results indicates quantitative differences in the relative contributions of free-energy components to the overall free-energy change between the two data sets, corresponding with those in the respective classical QSAR equations. The LERE-QSAR procedure was demonstrated to quantitatively reveal differences in the binding mechanisms between the two cases involving similar but different zinc-containing proteins at the electronic and atomic levels.

Reassessment of Hammett σ as an effective parameter representing intermolecular interaction energy-links between traditional and modern QSAR approaches.

The Hammett σ constant has for a long time been known to be one of most important linear free-energy related parameters that correlate with biological activity. It is a conventionally used electronic parameter in studies of enzymatic quantitative structure-activity relationships (QSAR). However, it is not necessarily obvious why σ represents variations in the free-energy change associated with the complex formation between a congeneric series of ligands with their target protein. So far, several powerful molecular calculations, such as the ab initio fragment molecular orbital (FMO) one, that are directly applicable to ligand-protein complexes have emerged. In this study, we comprehensively reevaluate experimentally derived parameter σ confirming it represents intermolecular interaction energy terms, by applying molecular orbital (MO) calculations to a simple ligand-protein complex model. The current results provide a rational and quantitative basis for bridging the gap between the traditional QSAR approach and 'the modern QSAR one', which involves the molecular calculations to evaluate the overall free-energy change for complex formation.

Correlation analyses on binding affinity of sialic acid analogues and anti-influenza drugs with human neuraminidase using ab initio MO calculations on their complex structures--LERE-QSAR analysis (IV).

We carried out full ab initio fragment molecular orbital (FMO) calculations for complexes comprising human neuraminidase-2 (hNEU2) and sialic acid analogues including anti-influenza drugs zanamivir (Relenza) and oseltamivir (Tamiflu) in order to examine the variation in the observed inhibitory activity toward hNEU2 at the atomic and electronic levels. We recently proposed the LERE (linear expression by representative energy terms)-QSAR (quantitative structure-activity relationship) procedure. LERE-QSAR analysis quantitatively revealed that the complex formation is driven by hydrogen-bonding and electrostatic interaction of hNEU2 with sialic acid analogues. The most potent inhibitory activity, that of zanamivir, is attributable to the strong electrostatic interaction of a positively charged guanidino group in zanamivir with negatively charged amino acid residues in hNEU2. After we confirmed that the variation in the observed inhibitory activity among sialic acid analogues is excellently reproducible with the LERE-QSAR equation, we examined the reason for the remarkable difference between the inhibitory potencies of oseltamivir as to hNEU2 and influenza A virus neuraminidase-1 (N1-NA). Several amino acid residues in close contact with a positively charged amino group in oseltamivir are different between hNEU2 and N1-NA. FMO-IFIE (interfragment interaction energy) analysis showed that the difference in amino acid residues causes a remarkably large difference between the overall interaction energies of oseltamivir with hNEU2 and N1-NA. The current results will be useful for the development of new anti-influenza drugs with high selectivity and without the risk of adverse side effects.

Density functional theory study of hydrogen atom abstraction from a series of para-substituted phenols: why is the Hammett σ(p)+ constant able to represent radical reaction rates?

The rate of hydrogen atom abstraction from phenolic compounds by a radical is known to be often linear with the Hammett substitution constant σ(+), defined using the S(N)1 solvolysis rates of substituted cumyl chlorides. Nevertheless, a physicochemical reason for the above "empirical fact" has not been fully revealed. The transition states of complexes between the 2,2-diphenyl-1-picrylhydrazyl radical (dpph·) and a series of para-substituted phenols were determined by DFT (Density Functional Theory) calculations, and then the activation energy as well as the homolytic bond dissociation energy of the O-H bond and charge distribution in the transition state were calculated. The heterolytic bond dissociation energy of the C-Cl bond and charge distribution in the corresponding para-substituted cumyl chlorides were calculated in parallel. Excellent correlations among σ(+), charge distribution, and activation and bond dissociation energies revealed quantitatively that there is a strong similarity between the two reactions, showing that the electron-deficiency of the π-electron system conjugated with a substituent plays a crucial role in determining rates of the two reactions. The results provide a new insight into and physicochemical understanding of σ(+) in the hydrogen abstraction from substituted phenols by a radical.

Correlation analyses on binding affinity of substituted benzenesulfonamides with carbonic anhydrase using ab initio MO calculations on their complex structures (II).

We proposed a novel QSAR (quantitative structure-activity relationship) procedure called LERE (linear expression by representative energy terms)-QSAR involving molecular calculations such as ab initio fragment molecular orbital and generalized Born/surface area ones. We applied LERE-QSAR to two datasets for the free-energy changes during complex formation between carbonic anhydrase and a series of substituted benzenesulfonamides. The first compound set (Set I) and the second one (Set II) include relatively small substituents and alkyl chains of different lengths in the benzene ring, respectively. Variation of the inhibitory activity in Set I is expressed as the combination of Hammett σ and the hydrophobic substituent constant π in classical QSAR, and variation in Set II only by π. LERE-QSAR analyses clearly revealed that effects of σ and π on the activity variations in Sets I and II are consistently explainable with the energy terms in the LERE formulation, and provide more detailed and direct information as to the binding mechanism. The proposed procedure was demonstrated to provide a quantitative basis for understanding ligand-protein interactions at the electronic and atomic levels.

Correlation analyses on binding affinity of sialic acid analogues with influenza virus neuraminidase-1 using ab initio MO calculations on their complex structures.

We carried out full ab initio molecular orbital calculations on complexes between neuraminidase-1 (N1-NA) in the influenza A virus and a series of eight sialic acid analogues including oseltamivir (Tamiflu) in order to quantitatively examine the binding mechanism and variation in the inhibitory potency at the atomic and electronic levels. FMO-MP2-IFIE (interfragment interaction energy at the MP2 level of ab initio fragment molecular orbital calculations) analyses quantitatively revealed (1) that the complex formation is driven by strong electrostatic interactions of charged functional groups in the analogues with ionized amino acid residues and water molecules in the active site of N1-NA, and (2) that the variation in the inhibitory potency among the eight analogues is determined by the dispersion and/or hydrophobic interaction energies of the 3-pentyl ether and charged amino moieties in oseltamivir with certain residues and water molecules in the active site of N1-NA. The current results will be useful for the development of new antiinfluenza drugs with high potency against various subtypes of wild-type and drug-resistant NAs.

Correlation analyses on binding affinity of substituted benzenesulfonamides with carbonic anhydrase using ab initio MO calculations on their complex structures.

Quantitative structure-activity relationship analyses on the free energy change during complex formation between substituted benzenesulfonamides (BSAs) and bovine carbonic anhydrase II (bCA II) were performed using generilized Born/surface area (GB/SA) and ab initio fragment molecular orbital (FMO) calculations for the whole complex structures. The result shows that the overall free energy change is governed by the contribution from solvation and dissociation free energy changes accompanying by complex formation. The FMO-IFIE (interfragment interaction energy) analysis quantitatively revealed that the intrinsic interaction energy of bCA II with BSAs is mostly from interactions with amino acid residues in the active site of bCA II. The "Zn block" (Zn(2+) and three histidine residues coordinated to Zn(2+)) in the active site shows the lowest interaction energy and the greatest variance of interaction energy with BSAs through their coordination interaction. The proposed procedure was demonstrated to provide a quantitative basis for understanding a ligand-protein interaction at electronic and atomic levels.

Expression and molecular dynamics studies on effect of amino acid substitutions at Arg344 in human cathepsin A on the protein local conformation.

Human lysosomal protective protein/cathepsin A (CathA) is a multifunctional protein that exhibits not only protective functions as to lysosomal glycosidases, i.e., neuraminidase 1 (NEU1) and beta-galactosidase (GLB), but also its own serine carboxypeptidase activity, and exhibits conserved structural similarity to yeast and wheat homologs (CPY and CPW). Our previous study revealed that the R344 (Arg344) residue in CathA could contribute to the binding and recognition of the serine peptidase inhibitor chymostatin. We examined here the effects of substitution of R344 with other amino acids, including A, D, E, G, I, K, M, N, P, Q, S, and V, denoted as R344X, including the wild-type CathA, on expression of CathA activity and intracellular processing. Among the mutant gene products, the 54-kDa precursor/zymogen with the R344D substitution was not processed to the 32/20-kDa mature form with CathA activity in a fibroblastic cell line derived from a galactosialidosis patient. Molecular dynamics (MD) simulations on the total twelve R344X mutants and the wild-type revealed that only R344D takes on a significantly different conformation of S293-D295 in the excision peptide (M285-R298) compared to the other R344X mutants; the side chains of S293 and D295 in R344D are exposed on the molecular surface, although those in the other twelve R344X mutants are buried inside the protein. The results of the current work strongly suggest that the distinct conformational change of the S293-D295 region in the R344D protein causes the processing defect of the 54-kDa precursor of the R344D mutant gene product in cultured cells.

Flavonoids as substrates and inhibitors of myeloperoxidase: molecular actions of aglycone and metabolites.

Myeloperoxidase (MPO), secreted by activated neutrophils and macrophages at the site of inflammation, may be implicated in the oxidation of protein/lipoprotein during the development of cardiovascular diseases. Flavonoids have been suggested to act as antioxidative and anti-inflammatory agents in vivo; however, their molecular actions have not yet been fully understood. In this study, we examined the molecular basis of the inhibitory effects of dietary flavonoids, such as quercetin, and their metabolites on the catalytic reaction of MPO using a combination of biological assays and theoretical calculation studies. Immunohistochemical staining showed that a quercetin metabolite was colocalized with macrophages, MPO, and dityrosine, an MPO-derived oxidation product of tyrosine, in human atherosclerotic aorta. Quercetin and the plasma metabolites inhibited the formation of dityrosine catalyzed by the MPO enzyme and HL-60 cells in a dose-dependent manner. Spectrometric analysis indicated that quercetin might act as a cosubstrate of MPO resulting in the formation of the oxidized quercetin. Quantitative structure-activity relationship studies showed that the inhibitory actions of flavonoids strongly depended not only on radical scavenging activity but also on hydrophobicity (log P). The requirement of a set of hydroxyl groups at the 3, 5, and 4'-positions and C2-C3 double bond was suggested for the inhibitory effect. The binding of quercetin and the metabolites to a hydrophobic region at the entrance to the distal heme pocket of MPO was also proposed by a computer docking simulation. The current study provides the structure-activity relationships for flavonoids as the anti-inflammatory dietary constituents targeting the MPO-derived oxidative reactions in vivo.

Simple computational models of type I/type II cells in Fas signaling-induced apoptosis.

Distinct apoptotic response of the type I/type II cells against Fas-ligand stimulation is considered to arise from the difference in dominant signaling pathways involved. In the type I cells, apoptotic signaling predominantly takes place via the direct activation of caspase-3 by activated caspase-8 (D channel) while mitochondrial pathway (M channel) plays a major role in the type II cells. To elucidate the selection mechanism of dominant pathway, we carried out systematic model analysis of the Fas signaling-induced apoptosis network. An increase in the expression level of caspase-8 induced a switch of dominant pathway from M- to D-channel (M-D transition), showing a phenotypic change from type II to type I cells. With the aid of sensitivity analysis and kinetic considerations, we succeeded in constructing a minimal network model relevant for the M-D transition, which revealed that mechanistic origin of the transition lies in the competition between the activated forms of caspase-8 and caspase-9 for their common substrate caspase-3. The pathway dominance was found to be primarily controlled by the balance between the activation rate of caspase-8 and the initial level of caspase-9. In the full network model, we showed that differential formation ability of the death-inducing signaling complex (DISC) can also induce M-D transition, in accordance with the experimental observations.