Search

Computational methods to calculate binding affinity in protein−ligand interaction are of immense interest because of obvious practical applications in structure-based drug design. Scoring functions attempt to calculate the variation in binding affinity of ligands−inhibitors bound to protein targets at various levels of theory. In this study we use semiempirical quantum mechanics to design a scoring function that can calculate the electrostatic interactions and solvation free energy expected during complexation. This physically based approach has the ability to capture binding affinity trends in a diverse range of protein−ligand complexes. We also show the predictive power of this scoring function within protein targets and its ability to score ligand poses docked to a protein target. We also demonstrate the ability of this scoring function to discriminate between native and decoy poses and highlight the crucial role played by electrostatic interactions in molecular recognition. Finally we compare the performance of our scoring function with other available scoring functions in the literature.

On the basis of molecular dynamics simulations, we investigate the dynamic properties of the carbapenem- and cephamycin-resistant dinuclear zinc metallo-β-lactamase from Bacteroides fragilis and its complex with a biphenyl tetrazole inhibitor, 2-butyl-6-hydroxy-3-[2‘-(1H-tetrazol-5-yl)biphenyl-4-ylmethyl]-3H-quinazolin-4-one 1 (L-159061). The results obtained with the newly developed force field parameters for the coordination environment of the catalytic zinc ions show that the active site gorge comprising major and minor loops gets deeper and narrower upon binding of the inhibitor, which supports the previous experimental implication that the structural flexibility of the loop structures plays a significant role in enzymatic action. In the presence of the inhibitor, the Trp32 side chain at the apex of the major loop covers the entrance of active site channel, thereby contributing to the stabilization of the enzyme−inhibitor complex. In addition to a direct coordination of the inhibitor tetrazole ring to the second zinc ion in the active site, the hydrogen bonding of Lys167 to the inhibitor carbonyl group and hydrophobic interactions between the inhibitor and side chains of loop residues prove to be significant binding forces of the enzyme−inhibitor complex.

Herein, we present theoretical results aimed at elucidating the origin of the kinetic preference for penicillins over cephalosporins characteristic of the TEM/SHV subgroup of class A β-lactamases. First, we study the conformational properties of cephalothin showing that the C2-down conformer of the dihydrothiazine ring is preferred over the C2-up one by ~2 kcal/mol in solution (0.4−1.4 kcal/mol in the gas phase). Second, the TEM-1 β-lactamase complexed with cephalothin is investigated by carrying out a molecular dynamics simulation. The ΔGbinding energy is then estimated using molecular mechanics Poisson−Boltzmann surface area (MM-PBSA) and quantum chemical PBSA (QM-PBSA) computational schemes. The preferential binding of benzylpenicillin over cephalothin is reproduced by the different energetic calculations, which predict relative ΔΔGbinding energies ranging from 1.8 to 5.7 kcal/mol. The benzylpenicillin/cephalothin ΔΔGbinding energy is most likely due to the lower efficacy of cephalosporins than that of penicillins in order to simultaneously bind the “carboxylate pocket” and the “oxyanion hole” in the TEM-1 active site.

“Fail early and fail fast” is the current paradigm that the pharmaceutical industry has adopted widely. Removing non-drug-like compounds from the drug discovery lifecycle in the early stages can lead to tremendous savings of resources. Thus, fast screening methods are needed to profile the large collection of synthesized and virtual libraries involved in the early stage. Solubility is one of the filters that are applied extensively to ensure that the compounds are reasonably soluble so that synthesis of the compounds and assay studies of pharmacokinetics and toxicity are feasible. To address this need, we have developed a fast quantitative structure−property relationship (QSPR) model for the prediction of aqueous solubility (at 298 K, unbuffered solution) from the molecular structures. Multiple linear regressions and genetic algorithms were used to develop the models. The model was based on a set of diverse compounds including small organic molecules and drug and drug-like species. The predicted solubility for the training and test sets agrees well with the experimental values. The coefficient of determination is R² = 0.84 for the training set of 775 compounds and the RMS error = 0.87. This model was validated on four sets of compounds. The RMS error for the 1665 compounds from the four validation data sets (including compounds from the Physician's Desk References and Comprehensive Medicinal Chemistry databases) is 1 log unit and the unsigned error is 0.77. This model does not require 3-D structure generation which is rather time-consuming. Using 2-D structure as input, this model is able to compute solubility for 90 000−700 000 compounds/h on a SGI Origin 2000 workstation. This kind of fast calculation allows the model to be used in data mining and screening of large synthesized or virtual libraries.

Sodium is very important as a counterion in biology. However, when used with the most common semiempirical Hamiltonians, such as AM1 or PM3, sodium is modeled as a point charge that can accept no electron density, called a “sparkle”. To better model sodium, we derived two sets of sodium parameters, which treat sodium on the same footing as other atoms parametrized in semiempirical methods. One set is compatible with the AM1 parameter set, while the second is compatible with PM3. These parameters were derived using a modified genetic algorithm with a diverse set of 71 compounds. The average unsigned error for the heats of formation was 10.3 kcal/mol for AM1 and 10.5 kcal/mol for PM3.

Herein we present results of a computational study on the benzylpenicillin acylation of the class A TEM1 β-lactamase via hydroxyl-only and hydroxyl and carboxylate assisted processes. These mechanisms correspond to a one-step Ser130-assisted process and a second route in which the β-lactam carboxylate and the Ser130 hydroxyl group help the proton transfer from the hydroxyl group of Ser70 to the β-lactam leaving N atom. The internal geometry of the reactive part of the TEM1-benzylpenicillin system is taken from a B3LYP/6-31+G* computational study on the methanol-assisted methanolysis reaction of a penicillin model compound (3α-carboxypenam). The 6-acylamino side chain and the 2-methyl groups of benzylpenicillin, together with the closer residues around the essential Ser70, are relaxed by carrying out geometry optimizations with a hybrid QM/MM method. The corresponding relative energies in the protein combine the B3LYP/6-31+G* electronic energies of the reactive subsystem with semiempirical PM3 energies of the TEM1-benzylpenicillin systems both in vacuo and in solution. The PM3 calculations on the TEM1-benzylpenicillin systems are performed with a Divide and Conquer linear-scaling method. The hydroxyl and carboxylate assisted pathway, which is the most favored one, is in agreement with the experimentally observed kinetic isotope effects and is also compatible with the effects of mutagenesis experiments on the Ser130 residue. These results suggest that a similar mechanism for the formation of acylenzyme intermediates could be relevant to other active-site serine penicillin-recognizing enzymes.

Drug discovery research is increasingly dedicated to biological screening on a massive scale, which seems to imply a basic rejection of many computer-assisted techniques originally designed to add rationality to the early stages of discovery. While ever-faster and more clever 3D methodologies continue to be developed and rejected as alternatives to indiscriminant screening, simpler tools based on 2D structure have carved a stable niche in the high-throughput paradigm of drug discovery. Their staying power is due in no small part to simplicity, ease of use, and demonstrated ability to explain structure−activity data. This observation led us to wonder whether an even simpler view of structure might offer an advantage over existing 2D and 3D methods. Accordingly, we introduce 1D representations of chemical structure, which are generated by collapsing a 3D molecular model or a 2D chemical graph onto a single coordinate of atomic positions. Atoms along this coordinate are differentiated according to elemental type, hybridization, and connectivity. By aligning 1D representations to match up identical atom types, a measure of overall structural similarity is afforded. In extensive structure−activity validation tests, 1D similarities consistently outperform both Daylight 2D fingerprints and Cerius² pharmacophore fingerprints, suggesting that this new, simple means of representing and comparing structures may offer a significant advantage over existing tried-and-true methods.

The IIA binding site of human serum albumin (HSA) preferentially binds hydrophobic organic anions of medium size (e.g., aspirin, benzylpenicillin, warfarin, etc.) and bilirubin. This binding ability is particularly important for the distribution, metabolism, and efficacy of drugs. In addition, HSA can also covalently link to different IIA substrates owing to the presence of a highly reactive residue, Lys199, which is strategically located in the IIA site. Herein, we present results of three restrained molecular dynamics (MD) simulations of the IIA binding site on the HSA protein. From these simulations, we have determined the influence that the ionization state of the key residue, Lys199, and the nearby Lys195 has on the structure and dynamics of the IIA binding site. When Lys199 is neutral the computed average distances for the most significant interresidue contacts are in good agreement with those estimated from the X-ray coordinates. The analysis of the solvent structure and dynamics indicates that the basic form of Lys199 is likely connected to the acid form of Lys195 through a network of H-bonding water molecules with a donor → acceptor character. The presence of these water bridges can be important for stabilizing the configuration of the IIA binding site and/or promoting a potential Lys195 → Lys199 proton-transfer process. These results suggest that both lysine residues located in the IIA binding site of HSA, Lys195 and Lys199, could play a combined and comparable chemical role. Our simulations also give insight into the binding of bilirubin to HSA.

Literature data on compounds both well- and poorly-absorbed in humans were used to build a statistical pattern recognition model of passive intestinal absorption. Robust outlier detection was utilized to analyze the well-absorbed compounds, some of which were intermingled with the poorly-absorbed compounds in the model space. Outliers were identified as being actively transported. The descriptors chosen for inclusion in the model were PSA and AlogP98, based on consideration of the physical processes involved in membrane permeability and the interrelationships and redundancies between available descriptors. These descriptors are quite straightforward for a medicinal chemist to interpret, enhancing the utility of the model. Molecular weight, while often used in passive absorption models, was shown to be superfluous, as it is already a component of both PSA and AlogP98. Extensive validation of the model on hundreds of known orally delivered drugs, “drug-like” molecules, and Pharmacopeia, Inc. compounds, which had been assayed for Caco-2 cell permeability, demonstrated a good rate of successful predictions (74−92%, depending on the dataset and exact criterion used).

In this article we investigate the importance of many-body and nonclassical effects, such as polarization and charge transfer, on the folding of the betanova protein. Our calculations show that these effects are crucial in stabilization of the system. Moreover, both polarization and charge transfer significantly alter the charge distribution of the system. Our detailed study shows that these fluctuations in charge are solely dependent on the local environment and not on the overall fold of the protein. Moreover, the contributions of polarization and charge transfer are roughly constant during the protein folding process. This means that the folding driving force is largely determined by the electrostatic energy. Our findings indicate that the folding of betanova can be accurately described by effective two-body potentials, despite the absence of explicit polarization and charge transfer in these models.

The paper presents a variational method to distribute the charge transferred (CT) between an ion and a dielectric continuum onto a boundary element (BE) grid used to represent the dielectric continuum media. The CT charge is transferred into the dielectric using a methodology that combines a semiempirical effective Hamiltonian with a modified Poisson−Boltzmann equation, which includes CT in the form of a surface charge density positioned at the dielectric interface. The variational method eliminates the bias of our earlier method used to distribute CT charge into the dielectric, i.e., the need to specify the solute atoms used in constructing the BE grid. The surface “pattern” of the magnitudes of the CT charges placed onto the BE grid is obtained by minimizing an energy functional defined as the sum of the CT charge self-interaction energy and the energies due to CT charge interaction with the reaction field (RF) and solute partial atomic charges. The total amount of charge, distributed on the BE grid, is constrained to the amount of charge transferred from/to the solute HOMO/LUMO molecular orbital, by using the method of Lagrange multipliers.

Molecular dynamics (MD) free energy perturbation (FEP) simulations were carried out in order to obtain insights into the structures and dynamics of a series of small organic solutes in water and water-saturated 1-octanol. Relative free energies of solvation were computed for each solute in both solvents, and these results were used to estimate the relative octanol/water partition coefficients (log Pow). The relative octanol/water partition coefficients were in good agreement with experimental log Pow values (average unsigned error = 0.74 log units), if one omits the acetamide−acetone simulation which proved problematic. Partition coefficients were also calculated using the newly developed GB/SA octanol continuum solvation model in order to compare the MD-FEP and continuum model results. Interestingly, the computationally much more efficient GB/SA calculations proved to be more accurate (average unsigned error in log Pow = 0.50 log units) than FEP for this set of 12 solutes.

In this paper we report molecular dynamics (MD) and free energy perturbation (FEP) studies carried out on enzyme−inhibitor (two hydroxamates that only differ by a carbon−carbon double bond) complexes of human fibroblast collagenase to obtain insights into the structural and energetic preferences of these inhibitors. We have developed a bonded model for the catalytic and structural zinc centers (Hoops, S. C.; et al. J. Am. Chem. Soc. 1991, 113, 8262−8270) where the electrostatic representation for this model was derived using a novel quantum-mechanical/molecular-mechanical (QM/MM) minimization procedure followed by electrostatic potential fitting. The resulting bonded model for the zinc ions was then used to generate MD trajectories for structural analysis and FEP studies. This model has satisfactorily reproduced the structural features of the active site, and furthermore, the FEP simulations gave relative free energies of binding in good agreement with experimental results. MD simulations in conjunction with the FEP are able to provide a structural explanation regarding why one hydroxamate inhibitor is favored over the other, and we are also able to make predictions about changes in the inhibitor that would enhance protein−inhibitor interactions.

Polyphosphazene single-substituent polymers were synthesized with the general formula [NP(OCH2CH2OCH2CH2XCH3)2] where X = oxygen for polymer 5 or X = sulfur for polymer 6. Characterization of these materials made use of ¹H, ¹³C, ¹⁵N, and ³¹P nuclear magnetic resonance (NMR) spectroscopy, differential scanning calorimetry, gel permeation chromatography, and elemental microanalysis. The polymers were complexed with LiSO3CF3 and AgSO3CF3 and examined both as solid electrolyte media and in the presence of dimethylformamide solvent. The ionic conductivities of these materials were determined at 25 °C through the use of complex impedance analysis. The mechanism of ionic conduction in the polymer−salt complexes was probed through an examination of ¹³C, ³¹P, and ¹⁵N NMR shifts and ¹³C NMR spin−lattice relaxation times (T1) for d7-DMF solutions. Molecular dynamics simulations were also carried out in order to investigate the interactions within the polymer−salt−DMF complexes.

Herein we describe the implementation of the multiple time step method, in conjunction with a reversible integrator and the Nosé−Hoover chain method for temperature and pressure control, into ROAR 1.0. We have extensively tested the MTS method on three systems:  an antifreeze peptide in water, the organic solvent dimethylformamide (DMF), and a dimyristoylphosphatidylcholine (DMPC)-based lipid bilayer. From these test simulations, we observe that the MTS method was capable of producing stable trajectories even when a long time step (e.g., 8 fs) is used, while the SHAKE method was unable to do so. The SHAKE method also disturbs the bond vibrational motion while MTS algorithm does not when the time step was smaller than 5 fs. We also observe that we can conservatively obtain a 2.5-fold speed-up using the MTS method over a SHAKE simulation using a 1 fs time step. Overall, the MTS method gives a solid speed-up over the traditional SHAKE method, while simultaneously giving much more stable trajectories.

In this paper we report a method for solving the Schrödinger equation for large molecules in solution which involved merging a linear scaling divide and conquer (D&C) semiempirical algorithm with the Poisson−Boltzmann (PB) equation. We then assess the performance of our self-consistent reaction field (SCRF) approach by comparing our D&C-PB calculations for a set of 29 neutral and 36 charged molecules with those obtained by ab initio GVB and DFT (B3LYP) methods, Cramer and Truhlar's semiempirical generalized-Born SM5 model, and with the experimental solvation free energies. Furthermore, we show that our SCRF method can be used to perform fully quantum mechanical calculations of proteins in solution in a reasonable amount of time on a modern workstation. We believe that all electrostatic interactions in biological systems require a quantum mechanical description in order to obtain an accurate representation. Thus, our new SCRF method should have an impact on the computational study of physical and chemical phenomena occurring in proteins and nucleic acids, which are, in general, strongly influenced by electrostatic interactions. Moreover, this may lead to novel insights into classic problems like protein folding or drug design.

We present a new energy decomposition scheme, which decomposes interaction energies into electrostatic, polarization, and charge transfer contributions. The method is founded on the divide and conquer approach but is only prone to errors introduced by the basis set and the applied Hamiltonian. We illustrate the method by AM1/PM3 interaction energy decompositions on bimolecular systems and 64 water systems. Our decomposition scheme also allows for screening of charge transfer or charge transfer and polarization from intermolecular interactions. Screening of these interactions is illustrated by the charge analysis of a screened 64 water system. Our calculations indicate the importance of charge transfer, even for intermolecular separations 0.5−1.0 Å from equilibrium.

We have implemented and tested a frozen density matrix (FDM) approximation to the basic divide and conquer (DC) semiempirical algorithm. Molecular dynamics and Monte Carlo simulations were performed to estimate the advantages of the method. Results were compared to those obtained from the original DC method and the combined quantum mechanical/molecular mechanical (QM/MM) method. We found that the FDM approximation speeds DC calculations up significantly, while only introducing small errors. We also found that the FDM DC scheme performs better that the standard QM/MM approach in terms of defining the electronic and structural properties of the systems studied herein.

A continuum solvation model has been developed for octanol on the basis of the generalized Born/surface area (GB/SA) formalism. This model was parameterized to fit experimental octanol free energies of solvation for a training set of 66 diverse organic solutes. The resulting model is able to reproduce experimental octanol free energies of solution with an unsigned average error of 0.50 kcal/mol when using the recently developed Merck molecular force field (MMFF). Results are poorer when the model is applied to other force fields such as OPLS, where the atomic charges vary significantly from MMFF. The GB/SA octanol model can be used in conjunction with computed or experimental aqueous free energies of solvation to calculate log Pow directly. The average unsigned error in calculated log Pow values for the same set of compounds is 0.37.