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.     

MCDOCK: A Monte Carlo simulation approach to the molecular docking problem

Prediction of the binding mode of a ligand (a drug molecule) to its macromolecular receptor, or molecular docking, is an important problem in rational drug design. We have developed a new docking method in which a non-conventional Monte Carlo (MC) simulation technique is employed. A computer program, MCDOCK, was developed to carry out the molecular docking operation automatically. The current version of the MCDOCK program (version 1.0) allows for the full flexibility of ligands in the docking calculations. The scoring function used in MCDOCK is the sum of the interaction energy between the ligand and its receptor, and the conformational energy of the ligand. To validate the MCDOCK method, 19 small ligands, the binding modes of which had been determined experimentally using X-ray diffraction, were docked into their receptor binding sites. To produce statistically significant results, 20 MCDOCK runs were performed for each protein–ligand complex. It was found that a significant percentage of these MCDOCK runs converge to the experimentally observed binding mode. The root-mean-square (rms) of all non-hydrogen atoms of the ligand between the predicted and experimental binding modes ranges from 0.25 to 1.84 Å for these 19 cases. The computational time for each run on an SGI Indigo2/R10000 varies from less than 1 min to 15 min, depending upon the size and the flexibility of the ligands. Thus MCDOCK may be used to predict the precise binding mode of ligands in lead optimization and to discover novel lead compounds through structure-based database searching.     

Near ab initio quality molecular electrostatic potential maps using hybridization displacement charges at PM3 level and effects of geometrical changes in amino groups

Charge distributions, dipole moments, and molecular electrostatic potentials (MEP) around several molecules consisting of carbon, nitrogen, oxygen, fluorine, sulfur, and chlorine atoms were studied using the PM3 semiempirical method and the results compared with those obtained using ab initio calculations at the RHF/6‐31G** level. Thus it is shown that relative MEP values near different atoms can be obtained using hybridization displacement charges (HDC) obtained by employing PM3 density matrices that usually agree quite satisfactorily with the ab initio ones. Further, positions of ab initio MEP minima are correctly located and the corresponding relative MEP values usually correctly predicted using the PM3(HDC) charges distributed continuously in three dimensions according to the forms of squares of valence s atomic orbitals. The necessary parameters for HDC calculations using the PM3 method were optimized. It is shown how within the frameworks of both PM3 and AM1 methods the π electrons or lone pairs associated with amino group nitrogen atoms and ring atoms can be satisfactorily treated in different situations. © 2001 John Wiley & Sons, Inc. Int J Quant Chem 82: 299–312, 2001     

Investigation of binding features: Effects on the interaction between CYP2A6 and inhibitors

A computational investigation has been carried out on CYP2A6 and its naphthalene inhibitors to explore the crucial molecular features contributing to binding specificity. The molecular bioactive orientations were obtained by docking (FlexX) these compounds into the active site of the enzyme. And the density functional theory method was further used to optimize the molecular structures with the subsequent analysis of molecular lipophilic potential (MLP) and molecular electrostatic potential (MEP). The minimal MLPs, minimal MEPs, and the band gap energies (the energy difference between the highest occupied molecular orbital and lowest unoccupied molecular orbital) showed high correlations with the inhibition activities (pIC₅₀s), illustrating their significant roles in driving the inhibitor to adopt an appropriate bioactive conformation oriented in the active site of CYP2A6 enzyme. The differences in MLPs, MEPs, and the orbital energies have been identified as key features in determining the binding specificity of this series of compounds to CYP2A6 and the consequent inhibitory effects. In addition, the combinational use of the docking, MLP and MEP analysis is also demonstrated as a good attempt to gain an insight into the interaction between CYP2A6 and its inhibitors. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Substituent Effect Parameters: Extending the Applications to Organometallic Chemistry

Typically, metal complexes are constituted of an acceptor metal ion and one or more Iigands containing the donor atoms. Accordingly, the properties of a metal complex are equally dependent on the nature of the metal ion and the ligands. Minute structural variations in the ligand will may result in linear changes in the respective energetic parameters and such linear relationships have paramount importance in organometallic chemistry. The variation in ligands is virtually limitless and substantial because of the extent of organic chemistry available for the modelling of desirable ligands, apart from the variation in metal ions. Anyhow, there is still a need for new parameters for the design and quantification of new ligands which in turn leads to the synthesis of metal complexes with possibly predictable chemical properties. Previous studies have demonstrated that quantum chemically derived molecular electrostatic potential (MESP) parameters can be listed as one of the superior quantifiers in this regard, which can act as an effective ligand electronic parameter. The interaction between the ligand part and metal-containing part will be crucial in assessing the reactivity of organometallic complexes. Here we are applying MESP based substituent constants derived from substituted benzenes to forecast the interaction energies in (pyr^*)W(CO)₅, (NHC^*)Mo(CO)₅ and (η⁶-arene^*)Cr(CO)₃ complexes. Ligands and metal ions are varied in each case for better understanding and transparency.     

A new method to determine electrostatic potential around a macromolecule in solution from molecular wave functions

The three-dimensional reference interaction site model integral equation theory (3D-RISM) combined with the ab initio molecular orbital method (3D-RISM-SCF) is applied to a solvated macromolecular system. The solvation structure around a solute molecule is obtained from the 3D-RISM integral equation under the electrostatic potential of the solute molecule, calculated by the ab initio molecular orbital theory. The electrostatic potential should be calculated on each grid point in the three-dimensional real space. Therefore, the calculation of the electrostatic potential is the most time consuming part in this method. In this article, we propose a new procedure to save the computational cost for calculating the electrostatic potential and the solvated fock matrix. The strategy of this procedure is to evaluate the electrostatic potential and the solvated fock matrix in different ways, depending on the distance between solute and solvent. Inside the repulsive cores of solute atoms, it is possible to avoid the calculation of electrostatic potential and solvated Fock matrix by assuming the potential to be infinity. In the region sufficiently far from solute, they are evaluated classically by putting the effective point charge on each atom. In the intermediate region, the electrostatic potential is evaluated directly by integrating the molecular orbitals of the solute molecule. The electronic structure and the energy gradient of Methionine-Enkephalin and solvation structure are estimated by using this procedure in aqueous solution, and are compared with the results from other procedures. The results are compared also with those from the continuum model.     

The atom assignment problem in automated de novo drug design. 3. Algorithms for optimization of fragment placement onto 3D molecular graphs

Summary Atom assignment onto 3D molecular graphs is a combinatoric problem in discrete space. If atoms are to be placed efficiently on molecular graphs produced in drug binding sites, the assignment must be optimized. An algorithm, based on simulated annealing, is presented for efficient optimization of fragment placement. Extensive tests of the method have been performed on five ligands taken from the Protein Data Bank. The algorithm is presented with the ligand graph and the electrostatic potential as input. Self placement of molecular fragments was monitored as an objective test. A hydrogen-bond option was also included, to enable the user to highlight specific needs. The algorithm performed well in the optimization, with successful replications. In some cases, a modification was necessary to reduce the tendency to give multiple halogenated structures. This optimization procedure should prove useful for automated de novo drug design.     

Development of a unique 3D interaction model of endogenous and synthetic peripheral benzodiazepine receptor ligands

Different classes of Peripheral-type Benzodiazepine Receptor (PBR) ligands were examined and common structural elements were detected and used to develop a rational binding model based on energetically allowed ligand conformations. Two lipophilic regions and one electrostatic interaction site are essential features for high affinity ligand binding, while a further lipophilic region plays an important modulator role. A comparative molecular field analysis, performed over 130 PBR ligands by means of the GRID/GOLPE methodology, led to a PLS model with both high fitting and predictive values (r² = 0.898, Q² = 0.761). The outcome from the 3D QSAR model and the GRID interaction fields computed on the putative endogenous PBR ligands DBI (Diazepam Binding Inhibitor) and TTN (Tetracontatetraneuropeptide) was used to identify the amino acids most probably involved in PBR binding. Three amino acids, bearing lipophilic side chains, were detected in DBI (Phe49, Leu47 and Met46) and in TTN (Phe33, Leu31 and Met30) as likely residues underlying receptor binding. Moreover, a qualitative comparison of the molecular electrostatic potentials of DBI, TTN and selected synthetic ligands indicated also similar electronic properties. Convergent results from the modeling studies of synthetic and endogenous ligands suggest a common binding mode to PBRs. This may help the rational design of new high affinity PBR ligands.     

Atomic radii: Incorporation of solvation effects

Atomic radii used to define the solute cavity in continuum-based methods are determined by reproducing the solvent-accessible surface defined as the loci of minima in a potential (solvent interaction potential) between the solute and a probe. This potential includes electrostatic interaction (ion–dipole, ion–quadrupole, and ion-induced dipole) terms as well as a Lennard–Jones energy term. The method alleviates the need to distinguish solute atoms in different chemical environments. These radii, when used in the calculation of solvation free energies, are shown to be superior to fixed atom-specific radii or to radii obtained from the electron isodensity surface from quantum-mechanical calculations. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1482–1493, 1998     

An anchor-dependent molecular docking process for docking small flexible molecules into rigid protein receptors

A molecular docking method designated as ADDock, anchor- dependent molecular docking process for docking small flexible molecules into rigid protein receptors, is presented in this article. ADDock makes the bond connection lists for atoms based on anchors chosen for building molecular structures for docking small flexible molecules or ligands into rigid active sites of protein receptors. ADDock employs an extended version of piecewise linear potential for scoring the docked structures. Since no translational motion for small molecules is implemented during the docking process, ADDock searches the best docking result by systematically changing the anchors chosen, which are usually the single-edge connected nodes or terminal hydrogen atoms of ligands. ADDock takes intact ligand structures generated during the docking process for computing the docked scores; therefore, no energy minimization is required in the evaluation phase of docking. The docking accuracy by ADDock for 92 receptor-ligand complexes docked is 91.3%. All these complexes have been docked by other groups using other docking methods. The receptor-ligand steric interaction energies computed by ADDock for some sets of active and inactive compounds selected and docked into the same receptor active sites are apparently separated. These results show that based on the steric interaction energies computed between the docked structures and receptor active sites, ADDock is able to separate active from inactive compounds for both being docked into the same receptor.     

Deterministic global optimization of molecular structures using interval analysis

The search for the global minimum of a molecular potential energy surface is a challenging problem. The molecular structure corresponding to the global minimum is of particular importance because it usually dictates both the physical and chemical properties of the molecule. The existence of an extremely large number of local minima, the number of which may increase exponentially with the size of the molecule, makes this global minimization problem extremely difficult. A new strategy is described here for solving such global minimization problems deterministically. The methodology is based on interval analysis, and provides a mathematical and computational guarantee that the molecular structure with the global minimum potential energy will be found. The technique is demonstrated using two sets of example problems. The first set involves a relatively simple potential model, and problems with up to 40 atoms. The second set involves a more realistic potential energy function, representative of those in current use, and problems with up to 11 atoms.     

Computational combinatorial ligand design: Application to human α-thrombin

Summary A new method is presented for computer-aided ligand design by combinatorial selection of fragments that bind favorably to a macromolecular target of known three-dimensional structure. Firstly, the multiple-copy simultaneous-search procedure (MCSS) is used to exhaustively search for optimal positions and orientations of functional groups on the surface of the macromolecule (enzyme or receptor fragment). The MCSS minima are then sorted according to an approximated binding free energy, whose solvation component is expressed as a sum of separate electrostatic and nonpolar contributions. The electrostatic solvation energy is calculated by the numerical solution of the linearized Poisson-Boltzmann equation, while the nonpolar contribution to the binding free energy is assumed to be proportional to the loss in solvent-accessible surface area. The program developed for computational combinatorial ligand design (CCLD) allows the fast and automatic generation of a multitude of highly diverse compounds, by connecting in a combinatorial fashion the functional groups in their minimized positions. The fragments are linked as two atoms may be either fused, or connected by a covalent bond or a small linker unit. To avoid the combinatorial explosion problem, pruning of the growing ligand is performed according to the average value of the approximated binding free energy of its fragments. The method is illustrated here by constructing candidate ligands for the active site of human -thrombin. The MCSS minima with favorable binding free energy reproduce the interaction patterns of known inhibitors. Starting from these fragments, CCLD generates a set of compounds that are closely related to high-affinity thrombin inhibitors. In addition, putative ligands with novel binding motifs are suggested. Probable implications of the MCSS-CCLD approach for the evolving scenario of drug discovery are discussed.     

New model potentials for sulfur-copper(I) and sulfur-mercury(II) interactions in proteins: from ab initio to molecular dynamics

We have developed new force field and parameters for copper(I) and mercury(II) to be used in molecular dynamics simulations of metalloproteins. Parameters have been derived from fitting of ab initio interaction potentials calculated at the MP2 level of theory, and results compared to experimental data when available. Nonbonded parameters for the metals have been calculated from ab initio interaction potentials with TIP3P water. Due to high charge transfer between Cu(I) or Hg(II) and their ligands, the model is restricted to a linear coordination of the metal bonded to two sulfur atoms. The experimentally observed asymmetric distribution of metal ligand bond lengths (r) is accounted for by the addition of an anharmonic (r3) term in the potential. Finally, the new parameters and potential, introduced into the CHARMM force field, are tested in short molecular dynamics simulations of two metal thiolates fragments in water. (Brooks BR et al. J Comput Chem 1983, 4, 1987.1).     

Estimation of important binding sites in compounds that interact with proteins

Proteins are one of the important substances in understanding biological activity, and many of them express the function by binding to other proteins or small molecules (ligands) on the molecular surface. This interaction often occurs in the hollows (pockets) on the molecular surface of the protein. It is known that when pockets are similar in structure and physical properties, they are likely to express similar functions and to bind similar ligands. Therefore, exploring the similarity of the structure and physical properties in pockets is very useful because it leads to the discovery of new ligands that are likely to bind. In addition, exploring the important structure when binding to the protein significant spot in the ligand will provide useful knowledge for the development of new ligands.In this study, we propose a method to search for proteins containing pockets that are structurally and physically similar to significant spot in the pocket of the analyzed protein, and to extract significant spots in the ligands that bind to them. We use feature points as data. Feature points are the 3-dimensional points that are extracted from 3D structure data of proteins with feature values quantifying hydrophobicity and electrostatic potential. The corresponding feature points are extracted by comparing structurally and physically the pockets of the search target proteins with the significant spot of the analyzed protein. By evaluating the similarity based on the comparison results of the feature values given to the extracted feature points, we search for proteins that are similar to the analyzed protein. From the ligands that bind to the searched proteins, atoms that are near the protein pocket and similar to the atoms in ligand binding to the analyzed protein are extracted. The site constituted by the extracted atoms is defined as a significant spot in the ligand.As a result of classifying ligands binding to the protein by using the extracted significant spot in the ligand, the effectiveness of the proposed method was confirmed.     

Nonexistence of local maxima in molecular electrostatic potential maps

It has been rigorously established by means of classical electrostatic arguments, that molecular electrostatic potential maps are devoid of local maxima. This forms a generalization of the earlier works of Politzer and co-workers which were restricted to the case of atoms.     

Examination of stabilty of molecular agglomerates in molecular crystals

We report an examination algorithm of stability of molecular aggregates based on the estimation of rigidity of intermolecular contacts in a crystal structure. The algorithm includes the intermolecular interaction energy calculation (in the atom-atom potential approximation) of a pair of molecules selected in the crystal structure. Further, the energy is minimized using a least-squares technique by varying the position and orientation of one of the molecules. The contact rigidity is quantitatively assessed by the minimized rms difference between the positions of the atoms in the original and optimized structures (Zorkii’s criterion). Every rigid contact revealed in the structure determines finite or infinite stable agglomerates. The paper presents the results of testing the computer program based on this algorithm with a number of real crystal structures previously determined by single crystal X-ray diffraction, and also the examples of the most common stable molecular agglomerates found with the aid of the program.     

The orientation of N-H...O=C and N-H...N hydrogen bonds in biological systems: How good is a point charge as a model for a hydrogen bonding atom?

In order to design new ligands for protein-binding sites of unknownstructure, it would be useful to predict the likely sites of hydrogenbonding of an unknown protein fragment to a known molecule. The positions ofmaxima and minima in the electrostatic potential at appropriate distancesfrom the van der Waals surface were calculated for various small molecules,nucleic acid bases, peptide units and amino acid side chains containinggroups which can form the biologically important N-H...O=C andN-H...N hydrogen bonds. Their ability to predict the positions of H andO/N in hydrogen bonded complexes, as predicted by optimising theelectrostatic interactions of pairs of such molecules constrained by themolecular shapes, was assessed. It is shown that extrema in theelectrostatic potential around the isolated molecules give worthwhilepredictions for the locations of hydrogen bonding partners. For moleculesbound by a single N-H...O=C hydrogen bond, the electrostatic maximumassociated with the H is usually less than 1 Å from an acceptor atom,while a C=O electrostatic minimum is generally less than 1.5 Å fromthe hydrogen bond proton. However, a significant number of hydrogen bondsform to the opposite lone pair from the electrostatic minimum, in which casethe separation is up to 3.3 Å. This reflects the broad electrostaticpotential well around a carbonyl oxygen between the lone pair directions.The model predicts when neighbouring atoms drastically change the hydrogenbonding characteristics of an N-H or C=O group. Although the geometries ofhydrogen bonded complexes are influenced by the other van der Waals contactsbetween the molecules, particularly multiple hydrogen bonds, theseinfluences are constant when considering hydrogen bonding to a specificuncharacterised binding site. Hence, the consideration of stericallyaccessible electrostatic extrema will be useful in the design of newligands.     

PEARLS: program for energetic analysis of receptor-ligand system

Analysis of the energetics of small molecule ligand-protein, ligand-nucleic acid, and protein-nucleic acid interactions facilitates the quantitative understanding of molecular interactions that regulate the function and conformation of proteins. It has also been extensively used for ranking potential new ligands in virtual drug screening. We developed a Web-based software, PEARLS (Program for Energetic Analysis of Ligand-Receptor Systems), for computing interaction energies of ligand-protein, ligand-nucleic acid, protein-nucleic acid, and ligand-protein-nucleic acid complexes from their 3D structures. AMBER molecular force field, Morse potential, and empirical energy functions are used to compute the van der Waals, electrostatic, hydrogen bond, metal-ligand bonding, and water-mediated hydrogen bond energies between the binding molecules. The change in the solvation free energy of molecular binding is estimated by using an empirical solvation free energy model. Contribution from ligand conformational entropy change is also estimated by a simple model. The computed free energy for a number of PDB ligand-receptor complexes were studied and compared to experimental binding affinity. A substantial degree of correlation between the computed free energy and experimental binding affinity was found, which suggests that PEARLS may be useful in facilitating energetic analysis of ligand-protein, ligand-nucleic acid, and protein-nucleic acid interactions. PEARLS can be accessed at http://ang.cz3.nus.edu.sg/cgi-bin/prog/rune.pl.     

Syntheses and crystal structures of three cesium salts: Cesium 5-sulfosalicylate, cesium 3,5-dinitrosalicylate and cesium 2,4-dinitrophenoxide monohydrate

In an attempt to probe a potential template role of the large alkali-metal cation cesium in organization of biorelevant ligands, 5-sulfosalicylate, 3,5-dinitrosalicylate and 2,4-dinitrophenol complexes of cesium were prepared and structurally investigated. The structures of cesium 5-sulfosalicylate, cesium 3,5-dinitrosalicylate and cesium 2,4-dinitrophenoxide monohydrate have been determined through X-ray diffraction analysis. The 5-sulfosalicylate anion has lost the proton at the −SO₃H group while the 3,5-dinitrosalicylate anion at −COOH group but both retains the usual intermolecular hydrogen bond between phenolic and carboxylic oxygen. In cesium 2,4-dinitrophenoxide monohydrate, the Cs⁺ cation is 12-coordinate by O atoms in anions and water molecules while the metal atoms in cesium 5-sulfosalicylate and cesium 3,5-dinitrosalicylate have coordination numbers 10 and 11, respectively, with an irregular coordination sphere made up exclusively of oxygen atoms. Even more in cesium 2,4-dinitrophenoxide monohydrate, the water molecules are in rare triply bridging positions between these cations. Both complexes have layer structures containing the cations and polar groups of the ligands in core domains sandwiched by the aromatic rings above and below. The organization of all layer structures appears to be governed mainly by steric effects and electrostatic forces with very little directional influence of the cations.     

Conformational sensitivity of polyether macrocycles to electrostatic potential: Partial atomic charges,molecular mechanics,and conformational prediction

Molecular recognition (whether by enzymes, the immune system, or chelating ligands) depends critically on molecular conformation. Molecular mechanics predicts energetically favorable molecular conformations by locating low energy conformations using an empirical fit of molecular potential energy as a function of internal coordinates. Molecular mechanics analysis of 18-crown-6 demonstrates that the nonbonded term (primarily the electrostatic part) is the largest contributor to the conformational energy. Nevertheless, common methods of treating the electrostatic interaction for 18-crown-6 yield inconsistent values for conformational energies partly because partial charges assigned to each atom can change with conformation due to through-space inductive effects which are not considered in most molecular mechanics programs. Similar findings from several other groups are reviewed to support our conclusions. We argue for care and caution in predicting conformational preferences of molecules with two or more highly polar atoms. We also discuss the desirability of using an empirical method of partial charge determination such as the charge equilibration algorithm of Rappé and Goddard (or a suitable generalization which includes polarization) as a method of including these effects in molecular mechanics and molecular dynamics calculations.