ProPose: a docking engine based on a fully configurable protein–ligand interaction model

Virtual high-throughput screening of molecular databases and in particular high-throughput protein–ligand docking are both common methodologies that identify and enrich hits in the early stages of the drug design process. Current protein–ligand docking algorithms often implement a program-specific model for protein–ligand interaction geometries. However, in order to create a platform for arbitrary queries in molecular databases, a new program is desirable that allows more manual control of the modeling of molecular interactions.For that reason, ProPose, an advanced incremental construction docking engine, is presented here that implements a fast and fully configurable molecular interaction and scoring model. This program uses user-defined, discrete, pharmacophore-like representations of molecular interactions that are transformed on-the-fly into a continuous potential energy surface, allowing for the incorporation of target specific interaction mechanisms into docking protocols in a straightforward manner. A torsion angle library, based on semi-empirical quantum chemistry calculations, is used to provide minimum energy torsion angles for the incremental construction algorithm. Docking results of a diverse set of protein–ligand complexes from the Protein Data Bank demonstrate the feasibility of this new approach.As a result, the seamless integration of pharmacophore-like interaction types into the docking and scoring scheme implemented in ProPose opens new opportunities for efficient, receptor-specific screening protocols. Figure ProPose — a fully configurable protein-ligand docking program — transforms pharmacophores into a smooth potential energy surface.This revised version was published online in October 2004 with corrections to the Graphical Abstract.     

Molecular modeling of the von Willebrand factor A2 Domain and the effects of associated type 2A von Willebrand disease mutations

A homology model for the A2 domain of von Willebrand factor (VWF) is presented. A large number of target–template alignments were combined into a consensus alignment and used for constructing the model from the structures of six template proteins. Molecular dynamics simulation was used to study the structural and dynamic effects of eight mutations introduced into the model, all associated with type 2A von Willebrand disease. It was found that the group I mutations G1505R, L1540P and S1506L cause significant deviations over multiple regions of the protein, coupled to significant thermal fluctuations for G1505R and L1540P. This suggests that protein instability may be responsible for their intracellular retention. The group II mutations R1597W, E1638K and G1505E caused single loop displacements near the physiologic VWF proteolysis site between Y1605–M1606. These modest structural changes may affect interactions between VWF and the ADAMTS13 protease. The group II mutations I1628T and L1503Q caused no significant structural change in the protein, suggesting that inclusion of the protease in this model is necessary for understanding their effect.Figure Homology model of the von Willebrand factor A2 domainElectronic Supplementary Material Supplementary material is available in the online version of this article at http://dx.doi.org/10.1007/s00894-004-0194-9     

Analysis of the interactions of ribonuclease inhibitor with kanamycin

The interaction of ribonuclease inhibitor (RI) with kanamycin was studied by molecular modeling. The preliminary binding model was constructed using the Affinity module of the Insight II molecular modeling program and the key residues involved in the combination of RI binding to kanamycin were determined. Meanwhile, we determined relevant surface characteristics determining the interaction behavior. The modeling results indicated that electrostatic interactions and H-bond forces may work as major factors for the molecular interaction between kanamycin and RI. The above results are useful for elucidating the molecular principles upon which the selectivity of a kanamycin is based. The quartz-crystal microbalance (QCM) is a new method usually used to monitor the binding function of macromolecules with samples online in a flow-injection analysis (FIA) system. The experimental results demonstrate that kanamycin has an extraordinary affinity to the basic protein RI, and our result is consistent with the molecular modeling results. These principles can in turn be used to study the molecular recognition mechanism and design a mimic of kanamycin for the development of new RI binders.Figure Proposed binding model of kanamycin to RI obtained by the computer-aided docking method and optimized with molecular mechanics with the CVFF forcefield.     

Prototropic tautomerism of imidazolone in aqueous solution: a density functional approach using the combined discrete/self-consistent reaction field (SCRF) models

A systematic investigation of the proton transfer in the keto-amino/enol tautomerization of imidazolone was undertaken. Calculations in aqueous solution were performed using both combined discrete/self-consistent reaction field (SCRF) and SCRF methods. Complexes containing one to three water molecules around the hydrophilic site of imidazolone were used for the combined discrete/SCRF calculations. The DFT results predict that the barrier height for non-water-assisted intramolecular proton transfer is very high (214.8 kJmol^–1). Hydrogen bonding between imidazolone and the water molecule(s) will dramatically lower the barrier by a concerted multiple proton transfer mechanism. The proton transfer process through a eight-member ring formed by imidazolone and two water molecules is found to be more efficient and the calculated barrier height is ca. 61 kJmol^–1.Figure DFT calculations in aqueous solution predict the H-bonding between imidazolone(IZ) and the water molecule(s) will dramatically lower the tautomeric barrier by a concerted multiple proton transfer mechanism, in which an eight-member ring structure formed by IZ and 2H₂O is found to be more efficient and the barrier is 60.8 kJ mol^–1, much less than 214.8 kJ mol^–1 in the non-water-assisted mechanism.     

Decomposition mechanism of 3-<Emphasis Type="Italic">N</Emphasis>-morpholinosydnonimine (<Emphasis Type="Bold">SIN-1</Emphasis>)—a density functional study on intrinsic structures and reactivities

Selected intrinsic aspects of the mode of action of 3-N-morpholinosydnonimine (SIN-1) and its follow-up products are investigated by means of density functional theory. Besides the well known radical-cationic Feelisch–Schoenafinger pathway, an alternative anionic route via a neutral radical is proposed and included in the study. The individual reaction pathways are followed. Most notably, the overall activation barrier for the cationic route is calculated to be 184.04 kcal mol^–1, while the one for the anionic route is predicted to be more favorable with 14.09 kcal mol^–1.Figure A mere conformational change is predicted by DFT theory to elicit NO ejection from the SIN-1a radical cation, formed after one electron oxidation of the N-nitrosohydrazine SIN-1a.

Homology modeling and PAPS ligand (cofactor) binding study of bovine phenol sulfotransferase

In order to understand the mechanisms of ligand binding and the interaction between the ligand and the bovine phenol sulfotransferase, (bSULT1A1, EC 2.8.2.1) a three-dimensional (3D) model of the bSULT1A1 is generated based on the crystal structure of the estrogen sulfotransferase (PDB code 1AQU) by using the InsightII/Homology module. With the aid of the molecular mechanics and molecular dynamics methods, the final refined model is obtained and is further assessed by Profile-3D and ProStat, which show that the refined model is reliable. With this model, a flexible docking study is performed and the results indicate that 3-phosphoadenosine-5- phosphosulfate (PAPS) is a more preferred ligand than coenzyme A (CoA), and that His108 forms hydrogen bond with PAPS, which is in good agreement with the experimental results. From these docking studies, we also suggest that Phe255, Phe24 and Tyr169 in bSULT1A1 are three important determinant residues in binding as they have strong van-der-Waals contacts with the ligand. The hydrogen–bonding interactions also play an important role for the stability of the complex. Our results may be helpful for further experimental investigations.Figure The final 3D-structure of bSULT1A1. The structure is obtained by energy minimizing an average conformation over the last 100 ps of MD simulation. The -helix is represented in red and the -sheet in yellow.     

Evolutionary trace analysis of ionotropic glutamate receptor sequences and modeling the interactions of agonists with different NMDA receptor subunits

The ionotropic N-methyl-d-aspartate (NMDA) receptor is of importance in neuronal development, functioning, and degeneration in the mammalian central nervous system. The functional NMDA receptor is a heterotetramer comprising two NR1 and two NR2 or NR3 subunits. We have carried out evolutionary trace (ET) analysis of forty ionotropic glutamate receptor (IGRs) sequences to identify and characterize the residues forming the binding socket. We have also modeled the ligand binding core (S1S2) of NMDA receptor subunits using the recently available crystal structure of NR1 subunit ligand binding core which shares ~40% homology with other NMDA receptor subunits. A short molecular dynamics simulation of the glycine-bound form of wild-type and double-mutated (D481N; K483Q) NR1 subunit structure shows considerable RMSD at the hinge region of S1S2 segment, where pore forming transmembrane helices are located in the native receptor. It is suggested that the disruption of domain closure could affect ion-channel activation and thereby lead to perturbations in normal animal behavior. In conclusion, we identified the amino acids that form the ligand-binding pocket in many ionotropic glutamate receptors and studied their hydrogen bonded and nonbonded interaction patterns. Finally, the disruption in the S1S2 domain conformation (of NR1 subunit- crystal structure) has been studied with a short molecular dynamics simulation and correlated with some experimental observations.Figure The figure shows the binding mechanism of glutamate with NR2B subunit of the NMDA receptor. Glutamate is shown in cpk, hydrogen bonds in dotted lines and amino acids in blue. The amino acids shown here are within a 4-Å radius of the ligand (glutamate)     

Chiral recognition of aromatic compounds by β-cyclodextrin based on bimodal complexation

The chiral recognition of the selected aromatic chiral compounds by native -cyclodextrin (-CD) based on bimodal complexation was studied using a flexible docking algorithm FDOCK. A quantitative empirical free energy relationship model was employed to predict the complex stability constants and the preferred binding modes. The results showed that the calculated complex stability constants are in good agreement with the experimental data. Furthermore, the main force responsible for host-guest complexation is the van der Waals force and the chiral molecules are completely included into the -CD cavity. The chiral recognition for the selected aromatic chiral compounds is the result of the van der Waals force counterbalancing with the other effects, such as the electrostatic interaction and the hydrophobic effect.Figure The favorable structures for the inclusion complexes of (S)_phenylbutyric with -CD. View in the plane normal to the Z-axis. -CD is shown in surface and (S)_phenylbutyric in CPK representation.     

Modeling the E. coli 4-hydroxybenzoic acid oligoprenyltransferase (ubiA transferase) and characterization of potential active sites

4-Hydroxybenzoate oligoprenyltransferase of E. coli, encoded in the gene ubiA, is an important key enzyme in the biosynthetic pathway to ubiquinone. It catalyzes the prenylation of 4-hydroxybenzoic acid in position 3 using an oligoprenyl diphosphate as a second substrate. Up to now, no X-ray structure of this oligoprenyltransferase or any structurally related enzyme is known. Knowledge of the tertiary structure and possible active sites is, however, essential for understanding the catalysis mechanism and the substrate specificity.With homology modeling techniques, secondary structure prediction tools, molecular dynamics simulations, and energy optimizations, a model with two putative active sites could be created and refined. One active site selected to be the most likely one for the docking of oligoprenyl diphosphate and 4-hydroxybenzoic acid is located near the N-terminus of the enzyme. It is widely accepted that residues forming an active site are usually evolutionary conserved within a family of enzymes. Multiple alignments of a multitude of related proteins clearly showed 100% conservation of the amino acid residues that form the first putative active site and therefore strongly support this hypothesis. However, an additional highly conserved region in the amino acid sequence of the ubiA enzyme could be detected, which also can be considered a putative (or rudimentary) active site. This site is characterized by a high sequence similarity to the aforementioned site and may give some hints regarding the evolutionary origin of the ubiA enzyme.Semiempirical quantum mechanical PM3 calculations have been performed to investigate the thermodynamics and kinetics of the catalysis mechanism. These results suggest a near S_N1 mechanism for the cleavage of the diphosphate ion from the isoprenyl unit. The 4-hydroxybenzoic acid interestingly appears not to be activated as benzoate anion but rather as phenolate anion to allow attack of the isoprenyl cation to the phenolate, which appeared to be the rate limiting step of the whole process according to our quantum chemical calculations. Our models are a basis for developing inhibitors of this enzyme, which is crucial for bacterial aerobic metabolism. Figure Structure of the model of ubiA oligoprenyltransferase derived from the photosynthetic reaction center (1PRC). Putative active amino acid residues and substrates are shown as capped sticks to describe their location and geometry in the putative active sites. The violet spheres identify Mg²⁺.This revised version was published online in April 2005 with corrections to Table 3 and the page make-up.     

3D QSAR studies on peroxisome proliferator-activated receptor γ agonists using CoMFA and CoMSIA

The peroxisome proliferator-activated receptors (PPARs) have increasingly become attractive targets for developing novel anti-type 2 diabetic drugs. We employed comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) to study three-dimensional quantitative structure–activity relationship (3D QSAR) based on existing agonists of PPAR (including five thiazolidinediones and 74 tyrosine-based compounds). Predictive 3D QSAR models with conventional r² and cross-validated coefficient (q²) values up to 0.974 and 0.642 for CoMFA and 0.979 and 0.686 for COMSIA were established using the SYBYL package. These models were validated by a test set containing 18 compounds. The CoMFA and CoMSIA field distributions are in general agreement with the structural characteristics of the binding pockets of PPAR, which demonstrates that the 3D QSAR models built here are very useful in predicting activities of novel compounds for activating PPAR.      

Structure and dynamics of Candida rugosa lipase: the role of organic solvent

The effect of organic solvent on the structure and dynamics of proteins was investigated by multiple molecular dynamics simulations (1 ns each) of Candida rugosa lipase in water and in carbon tetrachloride. The choice of solvent had only a minor structural effect. For both solvents the open and the closed conformation of the lipase were near to their experimental X-ray structures (C rms deviation 1–1.3 Å). However, the solvents had a highly specific effect on the flexibility of solvent-exposed side chains: polar side chains were more flexible in water, but less flexible in organic solvent. In contrast, hydrophobic residues were more flexible in organic solvent, but less flexible in water. As a major effect solvent changed the dynamics of the lid, a mobile element involved in activation of the lipase, which fluctuated as a rigid body about its average position. While in water the deviations were about 1.6 Å, organic solvent reduced flexibility to 0.9 Å. This increase rigidity was caused by two salt bridges (Lys85–Asp284, Lys75–Asp79) and a stable hydrogen bond (Lys75–Asn 292) in organic solvent. Thus, organic solvents stabilize the lid but render the side chains in the hydrophobic substrate-binding site more mobile. Figure Superimposition of open (black, PDB entry 1CRL) and closed (gray, PDB entry 1TRH) conformers of C. rugosa lipase. The mobile lid is indicatedThis revised version was published online in October 2004 with corrections to the Graphical Abstract.     

Combination of anti-hypertensive drugs: a molecular dynamics simulation study

The anti-hypertensive drugs amlodipine, atenolol and lisinopril, in ordinary and PEGylated forms, with different combined-ratios, were studied by molecular dynamics simulations using GROMACS software. Twenty simulation systems were designed to evaluate the interactions of drug mixtures with a dimyristoylphosphatidylcholine (DMPC) lipid bilayer membrane, in the presence of water molecules. In the course of simulations, various properties of the systems were investigated, including drug location, diffusion and mass distribution in the membrane; drug orientation; the lipid chain disorder as a result of drug penetration into the DMPC membrane; the number of hydrogen bonds; and drug surface area. According to the results obtained, combined drugs penetrate deeper into the DMPC lipid bilayer membrane, and the lipid chains remain ordered. Also, the combined PEGylated drugs, at a combination ratio of 1:1:1, enhance drug penetration into the DMPC membrane, reduce drug agglomeration, orient the drug in a proper angle for easy penetration into the membrane, and decrease undesirable lipotoxicity due to distorted membrane self-assembly and thickness.

Open image in new window

Graphical abstract ?

     

Molecular modeling of dipeptide and its analogous systems with water

Intermolecular hydrogen-bond interactions in the monohydrated complexes of formamide, N-methylacetamide and glycylglycine have been studied using ab initio and DFT methods. The geometries were optimized using second-order Møller–Plesset perturbation theory and the B3LYP DFT functional with the 6-311++G** basis set. It is observed that hydrogen-bond interactions at the carbonyl group of the peptide moiety are stronger than those at the amino group of the formamide and N-methylacetamide molecules. Because of the presence of cyclic hydrogen-bonding interactions in glycylglycine, the interaction at the amino group is higher than at the carbonyl. The ¹³C and ¹⁵N NMR shielding values were calculated for the non-hydrated and monohydrated complexes. Condensed Fukui functions have also been calculated for non-hydrated formamide, N-methylacetamide and glycylglycine molecules at the B3LYP/6-311++G** level of theory, and the results are discussed.Figure Structure of hydrated glycylglycine dipeptide     

Molecular dynamics simulation of hepatitis C virus IRES IIId domain: structural behavior,electrostatic and energetic analysis

The dynamic behavior of the HCV IRES IIId domain is analyzed by means of a 2.6-ns molecular dynamics simulation, starting from an NMR structure. The simulation is carried out in explicit water with Na⁺ counterions, and particle-mesh Ewald summation is used for the electrostatic interactions. In this work, we analyze selected patterns of the helix that are crucial for IRES activity and that could be considered as targets for the intervention of inhibitors, such as the hexanucleotide terminal loop (more particularly its three consecutive guanines) and the loop-E motif. The simulation has allowed us to analyze the dynamics of the loop substructure and has revealed a behavior among the guanine bases that might explain the different role of the third guanine of the GGG triplet upon molecular recognition. The accessibility of the loop-E motif and the loop major and minor groove is also examined, as well as the effect of Na⁺ or Mg²⁺ counterion within the simulation. The electrostatic analysis reveals several ion pockets, not discussed in the experimental structure. The positions of these ions are useful for locating specific electrostatic recognition sites for potential inhibitor binding. Figure Superposition of 14 structures representative of the evolution of IRES IIId RNA along 2.6-ns MD simulation     

Hydrogen bonding interactions in E- or Z-2-phenyl-3-( X'-pyridyl)propenoic acid ( X=2, 3 or 4) assemblies--a molecular modeling study

The agglomeration properties of stereoisomeric 2-phenyl-3-(X-pyridyl)propenoic acids (X=2, 3 or 4) were studied by the PM3 semiempirical quantum chemical method. Calculations revealed that dimers kept together by the intermolecular hydrogen bonding interactions of the carboxylic groups could be built from both stereoisomers irrespective of the position of the nitrogen heteroatom. The dimers of the Z-isomers could also be built through (aromatic)C–H...N hydrogen bonds between the dimer units. The longest agglomerate was the pentamer of the dimers when the nitrogen was in the 2 position. Longer hydrogen-bonded agglomerates than dimers could only be constructed from the E-isomer with the nitrogen in position 4. Here, the trimer of the dimers proved to be the longest hydrogen-bonded entity and similarly to the Z-isomers, the dimer units are kept together by (aromatic)C–H...N intermolecular hydrogen bonds.Figure The helical structure of the pentamer of dimers for the Z-2-phenyl-3-(2-pyridyl)propenoic acid molecule, when in the initial geometry (aromatic)C(4)–H...N close contact was enforced.     

Effect of Pt and Sn on the adsorption of n-heptane in γ-Al2O3 catalyst models

The Grand Canonical Monte-Carlo (GCMC) method has been used to carry out simulations of the adsorption of n-heptane in models of naphtha-reforming catalysts. Models used in the study differed in the number and distribution of metal atoms—Pt and Sn. The number of adsorbed n-heptane molecules grows linearly with increasing number of metal atoms. The effect of Pt content on the adsorption of n-heptane molecules is most distinct at approximately 100 kPa and within the lower range of the temperatures investigated. In the models of bimetallic catalysts, the effect of the two metals is additive.Figure Effect of Pt and Sn on number of n-heptane molecules adsorbed in Al₂O₃ catalyst in 773 K and 1000 kPa.      

Quantitative structure-activity relationship by CoMFA for cyclic urea and nonpeptide-cyclic cyanoguanidine derivatives on wild type and mutant HIV-1 protease

Abstact 3D-QSAR studies using the Comparative Molecular Field Analysis (CoMFA) methodology were conducted to predict the inhibition constants, K_i, and the inhibitor concentrations, IC₉₀ of 127 symmetrical and unsymmetrical cyclic urea and cyclic cyanoguanidine derivatives containing different substituent groups such as: benzyl, isopropyl, 4-hydroxybenzyl, ketone, oxime, pyrazole, imidazole, triazole and having anti-HIV-1 protease activities. A significant cross-validated correlation coefficient (q²) of 0.63 and a fitted correlation coefficient r² of 0.70 were obtained, indicating that the models can predict the anti-protease activity from poorly to highly active compounds reliably. The best predictions were obtained for: XV643 (predicted log 1/K_i=9.86), a 3,5-dimethoxy-benzyl cyclic urea derivate (molec60, predicted log 1/K_i=8.57) and a benzyl cyclic urea derivate (molec 61, predicted log 1/IC₉₀=6.87). Using the CoMFA method, we also predicted the biological activity of 14 cyclic urea derivatives that inhibit the HIV-1 protease mutants V82A, V82I and V82F. The predicted biological activities of the: (i) XNO63 (inhibitory activity on the mutant HIV-1 PR V82A), (ii) SB570 (inhibiting the mutant HIV-1 PR V82I) and also (iii) XV652 (during the interaction with the mutant HIV-1 PR V82F) were in good agreement with the experimental values.Figure Stereoview of the contour plots of the CoMFA steric and electrostatic fields. The favorable (indicated by blue polyhedra) and unfavorable (represented by red polyhedra) electrostatic areas and also the favorable (shown by green polyhedra) and unfavorable (shown by yellow polyhedra) steric areas formed around the most active molecule, 6a.     

Local sequential minimization of double stranded B-DNA using Monte Carlo annealing

A software algorithm has been developed to investigate the folding process in B-DNA structures in vacuum under a simple and accurate force field. This algorithm models linear double stranded B-DNA sequences based on a local, sequential minimization procedure. The original B-DNA structures were modeled using initial nucleotide structures taken from the Brookhaven database. The models contain information at the atomic level allowing one to investigate as accurately as possible the structure and characteristics of the resulting DNA structures. A variety of DNA sequences and sizes were investigated containing coding and non-coding, random and real, homogeneous or heterogeneous sequences in the range of 2 to 40 base pairs. The force field contains terms such as angle bend, Lennard-Jones, electrostatic interactions and hydrogen bonding which are set up using the Dreiding II force field and defined to account for the helical parameters such as twist, tilt and rise. A close comparison was made between this local minimization algorithm and a global one (previously published) in order to find out advantages and disadvantages of the different methods. From the comparison, this algorithm gives better and faster results than the previous method, allowing one to minimize larger DNA segments. DNA segments with a length of 40 bases need approximately 4 h, while 2.5 weeks are needed with the previous method. After each minimization the angles between phosphate–oxygen-carbon A₁, the oxygen–phosphate–oxygen A₂ and the average helical twists were calculated. From the generated fragments it was found that the bond angles are A₁=150°±2°and A₂=130°±10°, while the helical twist is 36.6°±2° in the A strand and A₁=150°±6° and A₂=130±6° with helical twist 39.6°±2° in the B strand for the DNA segment with the same sequence as the Dickerson dodecamer.Figure The final minimized DNA segment of the Dickerson dodecamer sequence represented by ball drawings and viewed (left) perpendicular and (right) down the helical axis     

Prediction of lower critical solution temperature of N-isopropylacrylamide–acrylic acid copolymer by an artificial neural network model

In this paper, we have investigated the lower critical solution temperature (LCST) of N-isopropylacrylamide–acrylic acid (NIPAAm-AAc) copolymer as a function of chain-transfer agent/initiator mole ratio, acrylic acid content of copolymer, concentration, pH and ionic strength of aqueous copolymer solution. Aqueous solutions with the desired properties were prepared from previously purified polymers, synthesized at 65 °C by solution polymerization using ethanol. The effects of each parameter on the LCST were examined experimentally.In addition, an artificial neural network model that is able to predict the lower cretical solution temperature was develeped. The predictions from this model compare well against both training and test data sets with an average error less than 2.53%.Figure Cross plot of predicted and experimental LCST values for the testing data set.     

Organization of rhodopsin molecules in native membranes of rod cells–an old theoretical model compared to new experimental data

It has been shown that rhodopsin forms an oligomer in the shape of long double rows of monomers. Because of the importance of rhodopsin as a template for all G protein-coupled receptors, its dimeric, tetrameric and higher-oligomeric structures also provide a useful pattern for similar structures in GPCRs. New experimental data published recently are discussed in the context of a proposed model of the rhodopsin oligomer 1N3M deposited in the protein data bank. The new rhodopsin structure at 2.2 Å resolution with all residues resolved as well as an electron cryomicroscopy structure from 2D crystals of rhodopsin are in agreement with the 1N3M model. Accommodation of movement of transmembrane helix VI, regarded as a major event during the activation of rhodopsin, in a steady structure of the oligomer is also discussed.Figure Superimposition of the 1U19 (red wire), 1GZM (purple wire) and 1N3M (blue wire) rhodopsin structures. Size of the wires is proportional to thermal factors of backbone C atoms, view parallel to the membrane.