Journal of Computational Chemistry

Illingworth, Christopher J. R.; Scott, Paul D.; Parkes, Kevin E. B.; Snell, Christopher R.; Campbell, Matthew P.; Reynolds, Christopher A.

doi: 10.1002/jcc.21561pmid: 20839295

Here, we describe a family of methods based on residue–residue connectivity for characterizing binding sites and apply variants of the method to various types of protein–ligand complexes including proteases, allosteric‐binding sites, correctly and incorrectly docked poses, and inhibitors of protein–protein interactions. Residues within ligand‐binding sites have about 25% more contact neighbors than surface residues in general; high‐connectivity residues are found in contact with the ligand in 84% of all complexes studied. In addition, a k‐means algorithm was developed that may be useful for identifying potential binding sites with no obvious geometric or connectivity features. The analysis was primarily carried out on 61 protein–ligand structures from the MEROPS protease database, 250 protein–ligand structures from the PDBSelect (25%), and 30 protein–protein complexes. Analysis of four proteases with crystal structures for multiple bound ligands has shown that residues with high connectivity tend to have less variable side‐chain conformation. The relevance to drug design is discussed in terms of identifying allosteric‐binding sites, distinguishing between alternative docked poses and designing protein interface inhibitors. Taken together, this data indicate that residue–residue connectivity is highly relevant to medicinal chemistry. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010

Chintapalli, Sree V.; Yew, Boon K.; Illingworth, Christopher J. R.; Upton, Graham J. G.; Reeves, Philip J.; Parkes, Kevin E. B.; Snell, Christopher R.; Reynolds, Christopher A.

doi: 10.1002/jcc.21562pmid: 20839296

Nonoverlapping closed loops of around 25–35 amino acids formed via nonlocal interactions at the loop ends have been proposed as an important unit of protein structure. This hypothesis is significant as such short loops can fold quickly and so would not be bound by the Leventhal paradox, giving insight into the possible nature of the funnel in protein folding. Previously, these closed loops have been identified either by sequence analysis (conservation and autocorrelation) or studies of the geometry of individual proteins. Given the potential significance of the closed loop hypothesis, we have explored a new strategy for determining closed loops from the insertions identified by the structural alignment of proteins sharing the same overall fold. We determined the locations of the closed loops in 37 pairs of proteins and obtained excellent agreement with previously published closed loops. The relevance of NMR structures to closed loop determination is briefly discussed. For cytochrome c, cytochrome b562 and triosephophate isomerase, independent folding units have been determined on the basis of hydrogen exchange experiments and misincorporation proton‐alkyl exchange experiments. The correspondence between these experimentally derived foldons and the theoretically derived closed loops indicates that the closed loop hypothesis may provide a useful framework for analyzing such experimental data. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010

Elking, Dennis M.; Perera, Lalith; Duke, Robert; Darden, Thomas; Pedersen, Lee G.

doi: 10.1002/jcc.21563pmid: 20839297

In standard treatments of atomic multipole models, interaction energies, total molecular forces, and total molecular torques are given for multipolar interactions between rigid molecules. However, if the molecules are assumed to be flexible, two additional multipolar atomic forces arise because of (1) the transfer of torque between neighboring atoms and (2) the dependence of multipole moment on internal geometry (bond lengths, bond angles, etc.) for geometry‐dependent multipole models. In this study, atomic force expressions for geometry‐dependent multipoles are presented for use in simulations of flexible molecules. The atomic forces are derived by first proposing a new general expression for Wigner function derivatives $\partial D_{m'm}^l /\partial \Omega$. The force equations can be applied to electrostatic models based on atomic point multipoles or Gaussian multipole charge density. Hydrogen‐bonded dimers are used to test the intermolecular electrostatic energies and atomic forces calculated by geometry‐dependent multipoles fit to the ab initio electrostatic potential. The electrostatic energies and forces are compared with their reference ab initio values. It is shown that both static and geometry‐dependent multipole models are able to reproduce total molecular forces and torques with respect to ab initio, whereas geometry‐dependent multipoles are needed to reproduce ab initio atomic forces. The expressions for atomic force can be used in simulations of flexible molecules with atomic multipoles. In addition, the results presented in this work should lead to further development of next generation force fields composed of geometry‐dependent multipole models. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010

Murata, Katsumi; Nagata, Naoya; Nakanishi, Isao; Kitaura, Kazuo

doi: 10.1002/jcc.21565pmid: 20839298

We previously reported that solvent dipole ordering (SDO) at the ligand binding site of a protein indicates an outline of the preferred shape and binding pose of the ligands. We suggested that SDO‐mimetic pseudo‐molecules that mimic the 3D shape of the SDO region could be used as molecular queries with a shape similarity matching method in virtual screening. In this work, a virtual screening method based on SDO, named SDOVS, was proposed. This method was applied to virtual screening of ligands for four typical drug target proteins and the performance compared with that of FRED (well‐known rigid docking method); the efficiency of SDOVS was demonstrated to be better than FRED. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010

Valiente, Pedro A.; Gil L., Alejandro; Batista, Paulo R.; Caffarena, Ernesto R.; Pons, Tirso; Pascutti, Pedro G

doi: 10.1002/jcc.21566pmid: 20839299

The standard parameterization of the Linear Interaction Energy (LIE) method has been applied with quite good results to reproduce the experimental absolute binding free energies for several protein–ligand systems. However, we found that this parameterization failed to reproduce the experimental binding free energy of Plasmepsin II (PlmII) in complexes with inhibitors belonging to four dissimilar scaffolds. To overcome this fact, we developed three approaches of LIE, which combine systematic approaches to predict the inhibitor‐specific values of α, β, and γ parameters, to gauge their ability to calculate the absolute binding free energies for these PlmII‐Inhibitor complexes. Specifically: (i) we modified the linear relationship between the weighted nonpolar desolvation ratio (WNDR) and the α parameter, by introducing two models of the β parameter determined by the free energy perturbation (FEP) method in the absence of the constant term γ, and (ii) we developed a new parameterization model to investigate the linear correlation between WNDR and the correction term γ. Using these parameterizations, we were able to reproduce the experimental binding free energy from these systems with mean absolute errors lower than 1.5 kcal/mol. © 2010 Wiley Periodicals, Inc. J Comput Chem 2010

Mandado, Marcos; Cordeiro, M. Natália D. S.

doi: 10.1002/jcc.21567pmid: 20839300

Thermodynamic stability of metal–aminoacid complexes in water is discussed in terms of the Gibbs free energy of water–ligand exchange processes, and the electronic stabilizing factors thoroughly investigated by means of 1‐electron and 2‐electron density properties. Hexacoordinated complexes formed between iron cations and glycine molecules acting as monodentate or bidentate ligands have been chosen as targets for the current study. Results agree with experimental findings, and complexes formed with bidentate ligands are found to be more stable than those formed with monodentate ones. The larger the number of the coordinated glycine molecules the more stable is the complex. Fe(III) complexes are more stable than Fe(II) ones, but differences are small and the Fe3+/Fe2+ exchange process appears to be energetically feasible for these complexes. Formation of the second glycine–iron interaction involving the amino nitrogen in the bidentate ligands is enthalpycally unfavorable but takes place due to the large entropy rise of the process. The larger stability of Fe(III) complexes is due however to the balance between energetic and solvation terms, which is favorable to these complexes. Electron density properties account satisfactorily for the electronic energy changes along the complex formation in terms of ligand–metal electron transfer and covalent bond orders. © 2010 Wiley Periodicals, Inc. J Comput Chem 2010

Sano, Eri; Li, Weihua; Yuki, Hitomi; Liu, Xinli; Furihata, Tomomi; Kobayashi, Kaoru; Chiba, Kan; Neya, Saburo; Hoshino, Tyuji

doi: 10.1002/jcc.21568pmid: 20839301

Cytochrome P450 (CYP) is deeply involved in the metabolism of chemicals including pharmaceuticals. Therefore, polymorphisms of this enzyme have been widely studied to avoid unfavorable side effects of drugs in chemotherapy. In this work, we performed computational analysis of the mechanism of the decrease in enzymatic activity for three typical polymorphisms in CYP 2C9 species: *2, *3, and *5. Based on the equilibrated structure obtained by molecular dynamics simulation, the volume of the binding pocket and the fluctuation of amino residues responsible for substrate holding were compared between the wild type and the three variants. Further docking simulation was carried out to evaluate the appropriateness of the binding pocket to accommodate substrate chemicals. Every polymorphic variant was suggested to be inferior to the wild type in enzymatic ability from the structural viewpoint. F‐G helices were obviously displaced outward in CYP2C9*2. Expansion of the binding pocket, especially the space near F′ helix, was remarkable in CYP2C9*3. Disappearance of the hydrogen bond between K helix and β4 loop was observed in CYP2C9*5. The reduction of catalytic activity of those variants can be explained from the deformation of the binding pocket and the consequent change in binding mode of substrate chemicals. The computational approach is effective for predicting the enzymatic activity of polymorphic variants of CYP. This prediction will be helpful for advanced drug design because calculations forecast unexpected change in drug efficacy for individuals. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010

Lee, Tzong‐Yi; Hsu, Justin Bo‐Kai; Lin, Feng‐Mao; Chang, Wen‐Chi; Hsu, Po‐Chiang; Huang, Hsien‐Da

doi: 10.1002/jcc.21569pmid: 20839302

Protein acetylation, which is catalyzed by acetyltransferases, is a type of post‐translational modification and crucial to numerous essential biological processes, including transcriptional regulation, apoptosis, and cytokine signaling. As the experimental identification of protein acetylation sites is time consuming and laboratory intensive, several computational approaches have been developed for identifying the candidates of experimental validation. In this work, solvent accessibility and the physicochemical properties of proteins are utilized to identify acetylated alanine, glycine, lysine, methionine, serine, and threonine. A two‐stage support vector machine was applied to learn the computational models with combinations of amino acid sequences, and the accessible surface area and physicochemical properties of proteins. The predictive accuracy thus achieved is 5% to 14% higher than that of models trained using only amino acid sequences. Additionally, the substrate specificity of the acetylated site was investigated in detail with reference to the subcellular colocalization of acetyltransferases and acetylated proteins. The proposed method, N‐Ace, is evaluated using independent test sets in various acetylated residues and predictive accuracies of 90% were achieved, indicating that the performance of N‐Ace is comparable with that of other acetylation prediction methods. N‐Ace not only provides a user‐friendly input/output interface but also is a creative method for predicting protein acetylation sites. This novel analytical resource is now freely available at http://N‐Ace.mbc.NCTU.edu.tw/. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010

Zhang, Hao; Zhang, Yun‐Hong

doi: 10.1002/jcc.21570pmid: 20839303

In the present article, two focal subjects, i.e., hydration of the NO 3− and associated ion species in the Mg(NO3)2 solution are researched by using the ab initio method. Nitrate ions with the hydration number of 1–6 are optimized at the HF/6‐31+G* level. Their relative energies, binding energies, and v1‐NO 3− frequencies are also presented. The investigation of the binding energies shows the hydration number is 3–6 in the solvent abundant environment. The associated species, including ion pairings, triple‐ and multiple‐ion clusters, are also optimized at the same level and their v1‐NO 3− frequencies are calculated for comparing with the results in experiments. From the comparison, the new associated process via aqueous free ions → solvent‐shared ion pairings → solvent‐shared triple and multiple ion clusters → contact multiple ion clusters → amorphous crystal is proposed. © 2010 Wiley Periodicals, Inc. J Comput Chem 2010

Du, Juan; Xi, Lili; Lei, Beilei; Lu, Jing; Li, Jiazhong; Liu, Huanxiang; Yao, Xiaojun

doi: 10.1002/jcc.21571pmid: 20839304

Structure‐based quantitative structure‐activity relationship (QSAR) studies on a series of checkpoint kinase 1 (Chk1) inhibitors were performed to find the key structural features responsible for their inhibitory activity. Molecular docking was employed to explore the binding mode of all inhibitors at the active site of Chk1 and determine the active conformation for the QSAR studies. Ligand and structure‐based descriptors incorporating the ligand‐receptor interaction were generated based on the docked complex. Genetic Algorithm‐Multiple Linear Regression (GA‐MLR) method was used to build 2D QSAR model. The 2D QSAR model gave a squared correlation coefficient R2 of 0.887, cross‐validated Q2 of 0.837 and the prediction squared correlation coefficient R 2pred of 0.849, respectively. Furthermore, three‐dimensional quantitative structure‐activity relationship (3D QSAR) model using comparative molecular field analysis (CoMFA) with R2 of 0.983, Q2 of 0.550 and R 2pred of 0.720 was also developed. The obtained results are helpful for the design of novel Chk1 inhibitors with improved activities. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010