Constrained search of conformational hyperspace

Summary We introduce a new method for determining pharmacophore or active site geometries by analysis of the structures of a series of active compounds. The method, constrained search, and the key concepts on which it is based, is described and illustrated by its application to 28 potent inhibitors of angiotensin-converting enzyme (ACE). The data set is one utilized by Mayer et al. [J. Comput.-Aided Mol. Design, 1 (1987) 3–16] to determine a unique geometry for the active site. Our experiment validated the previously reported results, obtained by a systematic search, while reducing the computer time requirement by more than two orders of magnitude. The experiment also identified a previously unrecognized alternative active site geometry for the ACE series.     

Strategies for the determination of pharmacophoric 3D database queries

Strategies are described for constructing pharmacophoric 3D database queries, based on aseries of active and inactive analogs. The results are highly selective database queries, whichare consistent with the generally accepted pharmacophore for a number of systems. Thefoundation of these strategies is the method of Mayer, Naylor, Motoc and Marshall [J.Comput.-Aided Mol. Design, 1 (1987) 3] for inferring a unique binding geometry forangiotensin-converting enzyme (ACE) inhibitors. The strategies described here generalize theirapproach to cases where the chemical features responsible for binding are not a prioriapparent, and to cases where the binding geometry deduced by that method is not unique. Thekey new insight, the selectivity principle, is to rank the multiple solutions produced by themethod of Mayer et al. by their selectivity, a value that is related to the proportion of adatabase that is returned as a database hit list. Retrospective analyses are described for D2-antagonists, ACE inhibitors, fibrinogen antagonists, and 2-antagonists.     

Alchemical prediction of hydration free energies for SAMPL

Hydration free energy calculations have become important tests of force fields. Alchemical free energy calculations based on molecular dynamics simulations provide a rigorous way to calculate these free energies for a particular force field, given sufficient sampling. Here, we report results of alchemical hydration free energy calculations for the set of small molecules comprising the 2011 Statistical Assessment of Modeling of Proteins and Ligands challenge. Our calculations are largely based on the Generalized Amber Force Field with several different charge models, and we achieved RMS errors in the 1.4-2.2 kcal/mol range depending on charge model, marginally higher than what we typically observed in previous studies (Mobley et al. in J Phys Chem B 111(9):2242-2254, 2007, J Chem Theory Comput 5(2):350-358, 2009, J Phys Chem B 115:1329-1332, 2011; Nicholls et al. in J Med Chem 51:769-779, 2008; Klimovich and Mobley in J Comput Aided Mol Design 24(4):307-316, 2010). The test set consists of ethane, biphenyl, and a dibenzyl dioxin, as well as a series of chlorinated derivatives of each. We found that, for this set, using high-quality partial charges from MP2/cc-PVTZ SCRF RESP fits provided marginally improved agreement with experiment over using AM1-BCC partial charges as we have more typically done, in keeping with our recent findings (Mobley et al. in J Phys Chem B 115:1329-1332, 2011). Switching to OPLS Lennard-Jones parameters with AM1-BCC charges also improves agreement with experiment. We also find a number of chemical trends within each molecular series which we can explain, but there are also some surprises, including some that are captured by the calculations and some that are not.     

On the efficiency of VBSCF algorithms,a comment on “An efficient algorithm for energy gradients and orbital optimization in valence bond theory”

We comment on the paper [Song et al., J. Comput. Chem. 2009, 30, 399]. and discuss the efficiency of the orbital optimization and gradient evaluation in the Valence Bond Self Consistent Field (VBSCF) method. We note that Song et al. neglect to properly reference Broer et al., who published an algorithm [Broer and Nieuwpoort, Theor. Chim. Acta 1988, 73, 405] to use a Fock matrix to compute a matrix element between two different determinants, which can be used for an orbital optimization. Further, Song et al. publish a misleading comparison with our VBSCF algorithm [Dijkstra and van Lenthe, J. Chem. Phys. 2000, 113, 2100; van Lenthe et al., Mol. Phys. 1991, 73, 1159] to enable them to favorably compare their algorithm with ours. We give detail timings in terms of different orbital types in the calculation and actual timings for the example cases. © 2012 Wiley Periodicals, Inc.     

A distance geometry heuristic for expanding the range of geometries sampled during conformational search

A recent study of crystal structures of protein-ligand complexes has shown that bioactive conformations tend to be more extended than random ones (Diller and Merz, J. Comput. Aid. Mol. Des. 2002, 16, 105-112). Existing conformational sampling techniques produce molecular conformations with a distribution of geometric sizes that may not cover that of the bioactive conformations. Here, we describe a simple heuristic for biasing the conformational search toward more extended or compact conformations, while maintaining excellent sampling. The method uses a boosting strategy to generate a series of conformations, each of which is at least as extended (or compact) as the previous one. We demonstrate that this method significantly expands the range of geometric sizes generated during the search and thus increases the efficiency of sampling bioactive conformations.     

A virtual library of constrained cyclic tetrapeptides that mimics all four side-chain orientations for over half the reverse turns in the protein data bank

Reverse turns are often recognition sites for protein/protein interactions and, therefore, valuable potential targets for determining recognition motifs in development of potential therapeutics. A virtual combinatorial library of cyclic tetrapeptides (CTPs) was generated and the bonds in the low-energy structures were overlapped with canonical reverse-turn Calpha-Cbeta bonds (Tran et al., J Comput Aided Mol Des 19(8):551-566, 2005) to determine the utility of CTPs as reverse-turn peptidomimetics. All reverse turns in the Protein Data Bank (PDB) with a crystal structures resolution < or = 3.0 A were classified into the same known canonical reverse-turn Calpha-Cbeta bond clusters (Tran et al., J Comput Aided Mol Des 19(8):551-566, 2005). CTP reverse-turn mimics were compiled that mimicked both the relative orientations of three of the four as well as all four Calpha-Cbeta bonds in the reverse turns of the PDB. 54% of reverse turns represented in the PDB had eight or more CTPs structures that mimicked the orientation of all four of the Calpha-Cbeta bonds in the reverse turn.     

Quantitative structure-activity relationships by neural networks and inductive logic programming. I. The inhibition of dihydrofolate reductase by pyrimidines

Summary Neural networks and inductive logic programming (ILP) have been compared to linear regression for modelling the QSAR of the inhibition of E. coli dihydrofolate reductase (DHFR) by 2,4-diamino-5-(substitured benzyl)pyrimidines, and, in the subsequent paper [Hirst, J.D., King, R.D. and Sternberg, M.J.E., J. Comput.-Aided Mol. Design, 8 (1994) 421], the inhibition of rodent DHFR by 2,4-diamino-6,6-dimethyl-5-phenyl-dihydrotriazines. Cross-validation trials provide a statistically rigorous assessment of the predictive capabilities of the methods, with training and testing data selected randomly and all the methods developed using identical training data. For the ILP analysis, molecules are represented by attributes other than Hansch parameters. Neural networks and ILP perform better than linear regression using the attribute representation, but the difference is not statistically significant. The major benefit from the ILP analysis is the formulation of understandable rules relating the activity of the inhibitors to their chemical structure.     

Do benzodiazepines mimic reverse-turn structures?

The role of benzodiazepine derivatives (BZD) as a privileged scaffold that mimics beta-turn structures (Ripka et al. (1993) Tetrahedron 49:3593-3608) in peptide/protein recognition was reexamined in detail. Stable BZD ring conformers were determined with MM3, and experimental reverse-turn structures were extracted from the basis set of protein crystal structures previously defined by Ripka et al. Ideal beta-turns were also modeled and similarly compared with BZD conformers. Huge numbers of conformers were generated by systematically scanning the torsional degrees of freedom for BZDs, as well as those of ideal beta-turns for comparison. Using these structures, conformers of BZDs were fit to experimental structures as suggested by Ripka et al., or modeled classical beta-turn conformers, and the root-mean-square deviation (RMSD) values were calculated for each pairwise comparison. Pairs of conformers with the smallest RMSD values for overlap of the four alpha-beta side-chain orientations were selected. All overlaps of BZD conformers with experimental beta-turns yielded one or more comparisons where the least RMSD was significantly small, 0.48-0.86 angstroms, as previously suggested. Utilizing a different methodology, the overall conclusion that benzodiazepines could serve as reverse-turn mimetics of Ripka et al. is justified. The least RMSD values for the overlap of BZDs and modeled classical beta-turns were also less than 1 angstrom. When comparing BZDs with experimental or classical beta-turns, the set of experimental beta-turns selected by Ripka et al. fit the BZD scaffolds better than modeled classical beta-turns; however, all the experimental beta-turns did not fit a particular BZD scaffold better. A single BZD ring conformation, and/or chiral orientation, can mimic some, but not all, of the experimental beta-turn structures. BZD has two central ring conformations and one chiral center that explains why the four variations of the BZD scaffold can mimic all types of beta-turn structure examined. It was found, moreover, that the BZD scaffold also mimics each of the nine clusters of experimental orientations of side chains of reverse turns in the Protein Data Bank, when the new classification scheme for the four side-chain directions (the relative orientations of alpha-beta vectors of residues i through i+3) was considered (Tran et al. (2005) J Comput-Aided Mol Des 19:551-566).     

Designing novel nicotinic agonists by searching a database of molecular shapes

Summary We introduce an approach by which novel ligands can be designed for a receptor if a pharmacophore geometry has been established and the receptor-bound conformations of other ligands are known. We use the shape-matching method of Kuntz et al. [J. Mol. Biol., 161 (1982) 269–288] to search a database of molecular shapes for those molecules which can fit inside the combined volume of the known ligands and which have interatomic distances compatible with the pharmacophore geometry. Some of these molecules are then modified by interactive modeling techniques to better match the chemical properties of the known ligands. Our shape database (about 5000 candidate molecules) is derived from a subset of the Cambridge Crystallographic Database [Allen et al., Acta Crystallogr., Sect. B,35 (1979) 2331–2339]. We show, as an example, how several novel designs for nicotinic agonists can be derived by this approach, given a pharmacophore model derived from known agonists [Sheridan et al., J. Med. Chem., 29 (1986) 889–906]. This report complements our previous report [DesJarlais et al., J. Med. Chem., in press], which introduced a similar method for designing ligands when the structure of the receptor is known.     

Partition coefficients for the SAMPL5 challenge using transfer free energies

SAMPL challenges (Mobley et al. in J Comput Aided Mol Des 28:135–150, 2014; Skillman in J Comput Aided Mol Des 26:473–474, 2012; Geballe in J Comput Aided Mol Des 24:259–279, 2010; Guthrie in J Phys Chem B 113:4501–4507, 2009) provide excellent opportunities to assess theoretical approaches on new data sets with a goal of gaining greater insight towards protein and ligand modeling. In the SAMPL5 experiment, cyclohexane–water partition coefficients were determined using a vertical solvation scheme in conjunction with the SMD continuum solvent model. Several DFT functionals partnered with correlation consistent basis sets were evaluated for the prediction of the partition coefficients. The approach chosen for the competition, a B3PW91 vertical solvation scheme, yields a mean absolute deviation of 1.9 logP units and performs well at estimating the correct hydrophilicity and hydrophobicity for the full SAMPL5 molecule set.     

A theoretical study of P4O10: vibrational analysis, infrared and Raman spectra

The normal mode frequencies and the corresponding vibrational assignments of tetraphosphorus decaoxide (P4O10) in tetrahedral (Td) symmetry are examined theoretically and experimentally. The Gaussian 98 set of quantum chemistry codes at the HF/6-311G*, MP2/6-311G*, and DFT/B3LYP/6-311G* levels of theory are used. By comparison to experimental normal mode frequencies deduced by Gilliam et al. [J. Phys. Chem. B 107 (2003) 2892], Chapman [Spectrochim. Acta A, 24 (1968) 1687], Beattie et al. [J. Chem. Soc. A (1970) 449], Konings et al. [J. Mol. Spectrosc. 152 (1992) 29] and the present work, correction factors for predominant vibrational motions are determined and compared. Normal modes were decomposed into five non-redundant motions (P-O stretch, P=O stretch, P-O-P bend, P-O-P wag, and P=O wag). Standard deviations found for the HF, MP2, and DFT corrected frequencies compared to experiment are 9, 5, and 4 cm(-1), respectively. Electron distribution for selected molecular orbitals is considered.     

Pharmacophore and receptor models for neurokinin receptors

Three neurokinin (NK) antagonist pharmacophore models (Models 1-3) accounting for hydrogen bonding groups in the 'head' and 'tail' of NK receptor ligands have been developed by use of a new procedure for treatment of hydrogen bonds during superimposition. Instead of modelling the hydrogen bond acceptor vector in the strict direction of the lone pair, an angle is allowed between the hydrogen bond acceptor direction and the ideal lone pair direction. This approach adds flexibility to hydrogen bond directions and produces more realistic RMS values. By using this approach, two novel pharmacophore models were derived (Models 2 and 3) and a hydrogen bond acceptor was added to a previously published NK2 pharmacophore model [Poulsen et al., J. Comput.-Aided Mol. Design, 16 (2002) 273] (Model 1). Model 2 as well as Model 3 are described by seven pharmacophore elements: three hydrophobic groups, three hydrogen bond acceptors and a hydrogen bond donor. Model 1 contains the same hydrophobic groups and hydrogen bond donor as Models 2 and 3, but only one hydrogen bond acceptor. The hydrogen bond acceptors and donor are represented as vectors. Two of the hydrophobic groups are always aromatic rings whereas the other hydrophobic group can be either aromatic or aliphatic. In Model 1 the antagonists bind in an extended conformation with two aromatic rings in a parallel displaced and tilted conformation. Model 2 has the same two aromatic rings in a parallel displaced conformation whereas Model 3 has the rings in an edge to face conformation. The pharmacophore models were evaluated using both a structure (NK receptor homology models) and a ligand based approach. By use of exhaustive conformational analysis (MMFFs force field and the GB/SA hydration model) and least-squares molecular superimposition studies, 21 non-peptide antagonists from several structurally diverse classes were fitted to the pharmacophore models. More antagonists could be fitted to Model 2 with a low RMS and a low conformational energy penalty than to Models 1 and 3. Pharmacophore Model 2 was also able to explain the NK1, NK2 and NK3 subtype selectivity of the compounds fitted to the model. Three NK 7TM receptor models were constructed, one for each receptor subtype. The location of the antagonist binding site in the three NK receptor models is identical. Compounds fitted to pharmacophore Model 2 could be docked into the NK1, NK2 and NK3 receptor models after adjustment of the conformation of the flexible linker connecting the head and tail. Models I and 3 are not compatible with the receptor models.     

Constrained proper sampling of conformations of transition state ensemble of protein folding

Characterizing the conformations of protein in the transition state ensemble (TSE) is important for studying protein folding. A promising approach pioneered by Vendruscolo et al. [Nature (London) 409, 641 (2001)] to study TSE is to generate conformations that satisfy all constraints imposed by the experimentally measured φ values that provide information about the native likeness of the transition states. Fai?sca et al. [J. Chem. Phys. 129, 095108 (2008)] generated conformations of TSE based on the criterion that, starting from a TS conformation, the probabilities of folding and unfolding are about equal through Markov Chain Monte Carlo (MCMC) simulations. In this study, we use the technique of constrained sequential Monte Carlo method [Lin et al., J. Chem. Phys. 129, 094101 (2008); Zhang et al. Proteins 66, 61 (2007)] to generate TSE conformations of acylphosphatase of 98 residues that satisfy the φ-value constraints, as well as the criterion that each conformation has a folding probability of 0.5 by Monte Carlo simulations. We adopt a two stage process and first generate 5000 contact maps satisfying the φ-value constraints. Each contact map is then used to generate 1000 properly weighted conformations. After clustering similar conformations, we obtain a set of properly weighted samples of 4185 candidate clusters. Representative conformation of each of these cluster is then selected and 50 runs of Markov chain Monte Carlo (MCMC) simulation are carried using a regrowth move set. We then select a subset of 1501 conformations that have equal probabilities to fold and to unfold as the set of TSE. These 1501 samples characterize well the distribution of transition state ensemble conformations of acylphosphatase. Compared with previous studies, our approach can access much wider conformational space and can objectively generate conformations that satisfy the φ-value constraints and the criterion of 0.5 folding probability without bias. In contrast to previous studies, our results show that transition state conformations are very diverse and are far from nativelike when measured in cartesian root-mean-square deviation (cRMSD): the average cRMSD between TSE conformations and the native structure is 9.4 A? for this short protein, instead of 6 A? reported in previous studies. In addition, we found that the average fraction of native contacts in the TSE is 0.37, with enrichment in native-like β-sheets and a shortage of long range contacts, suggesting such contacts form at a later stage of folding. We further calculate the first passage time of folding of TSE conformations through calculation of physical time associated with the regrowth moves in MCMC simulation through mapping such moves to a Markovian state model, whose transition time was obtained by Langevin dynamics simulations. Our results indicate that despite the large structural diversity of the TSE, they are characterized by similar folding time. Our approach is general and can be used to study TSE in other macromolecules.     

Molecular Mechanics Investigation of an Adenine-Adenine Non-Canonical Pair Conformational Change

Conformational changes are important in RNA for binding and catalysis and understanding these changes is important for understanding how RNA functions. Computational techniques using all-atom molecular models can be used to characterize conformational changes in RNA. These techniques are applied to an RNA conformational change involving a single base pair within a nine base pair RNA duplex. The Adenine-Adenine (AA) non-canonical pair in the sequence 5'GGUGAAGGCU3' paired with 3'PCCGAAGCCG5', where P is Purine, undergoes conformational exchange between two conformations on the timescale of tens of microseconds, as demonstrated in a previous NMR solution structure [Chen, G., et al., Biochemistry, 2006. 45: 6889-903]. The more populated, major, conformation was estimated to be 0.5 to 1.3 kcal/mol more stable at 30 °C than the less populated, minor, conformation. Both conformations are trans-Hoogsteen/sugar edge pairs, where the interacting edges on the adenines change with the conformational change. Targeted Molecular Dynamics (TMD) and Nudged Elastic Band (NEB) were used to model the pathway between the major and minor conformations using the AMBER software package. The adenines were predicted to change conformation via intermediates in which they are stacked as opposed to hydrogen-bonded. The predicted pathways can be described by an improper dihedral angle reaction coordinate. Umbrella sampling along the reaction coordinate was performed to model the free energy profile for the conformational change using a total of 1800 ns of sampling. Although the barrier height between the major and minor conformations was reasonable, the free energy difference between the major and minor conformations was the opposite of that expected based on the NMR experiments. Variations in the force field applied did not improve the misrepresentation of the free energies of the major and minor conformations. As an alternative, the Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) approximation was applied to predict free energy differences between the two conformations using a total of 800 ns of sampling. MM-PBSA also incorrectly predicted the major conformation to be higher in free energy than the minor conformation.     

Conformational variety for the ansa chain of rifamycins: Comparison of observed crystal structures and molecular dynamics simulations

The antibiotic activity (via inhibition of DNA-dependent RNA polymerase, DDRP) of rifamycins has been correlated to the conformation of the ansa chain, which can be described by means of 17 torsion angles defined along the ansa backbone. It has been shown that favourable or unfavourable conformations of the ansa chain in rifamycin crystals are generally diagnostic of activity or inactivity against isolated DDRP. The principles of structure correlation suggest that the torsional variety observed in rifamycin crystals should mimic the dynamic flexibility of the ansa chain in solution. Twenty-six crystal structures of rifamycins are grouped into two classes (active and non-active). For each class the variance of the 17 ansa backbone torsion angles is analysed. Active compounds show a well-defined common pattern, while non-active molecules are more scattered, mainly due to steric constraints forcing the molecules into unfavourable conformations. The experimental distributions of torsion angles are compared to the torsional freedom of the ansa chain simulated by molecular dynamics calculations performed at different temperatures and conditions on rifamycin S and rifamycin O, which represent a typical active and a typical sterically constrained molecule, respectively. It is shown that the torsional variety found in the crystalline state samples the dynamic behaviour of the ansa chain for active compounds. The methods of circular statistics are illustrated to describe torsion angle distributions.     

Systematic search in conformational analysis

The coupling of conformation to activity and reactivity is a widely accepted concept, and as such has driven the development of tools which execute conformational searches in rapid and robust fashion [T.F. Havel, Prog. Biophys. Molec. Biol., 56 (1991) 43–78; A.R. Leach, In Rev. Comput. Chem.; K.B. Lipkowitz and D.B. Boyd, Ed.; VCH Publishers, Inc.: New York, N.Y., 1991, Vol. II, pp. 1–55]. Among the aims of these methods are the determination of a complete set of local minima from which the global energy minimum can be identified, or the generation of conformations consistent with constraints derived from SAR or structural studies. Most methods fall into two broad categories: those which are random or stochastic, and those which are systematic. Yet another group consists of those which are based on heuristics and artificial intelligence [A.R. Leach, K. Prout, D.P. Dolata, J. Comput. Chem. 11 (1990) 680–693]. The first category is typified by molecular dynamics [W.F. van Gunsteren and H.J.C. Berendsen, Angew. Chem. Int. Ed. Eng., 29 (1990) 992–1023], Monte Carlo [M.P. Alien and D.J. Tildesley, Computer Simulation of Liquids, Oxford Science Publications, 1989], distance geometry [J.M. Blaney and J.S. Dixon, in K.B. Lipkowitz and D.B. Boyd (Eds.), Reviews in Computational Chemistry, VCH, New York, Vol. 5, pp. 299–335, 1994], and other approaches [M. Saunders, J. Comput. Chem., 10 (1989) 203–208] in which the path by which conformational space is examined is ideally completely random, but bounded by the geometries of covalent bond lengths and angles. In traditional systematic searches, the variable to be examined, e.g. torsion angles, is divided into a regular grid. Each and every grid point is evaluated in a systematic fashion to determine its validity. The path through the grid points is regular and defined. In principle, systematic search can, within the resolution of the grid, identify all sterically allowed conformations of a molecule. Consequently, systematic search is an ideal tool for conformational analysis because it is not path dependent and cannot become entrapped in local minima. In this article we review some of the basics of systematic search, algorithmic improvements that have enhanced its speed, and new developments that have increased its accuracy by moving away from the limitations of a fixed torsional grid.     

Alignment of molecules by the Monte Carlo optimization of molecular similarity indices

3D-QSAR uses statistical techniques to correlate calculated structural properties with target properties like biological activity. The comparison of calculated structural properties is dependent upon the relative orientations of molecules in a given data set. Typically molecules are aligned by performing an overlap of common structural units. This “alignment rule” is adequate for a data set, that is closely related structurally, but is far more difficult to apply to either a diverse data set or on the basis of some structural property other than shape, even for sterically similar molecules. In this work we describe a new algorithm for molecular alignment based upon optimization of molecular similarity indices. We show that this Monte Carlo based algorithm is more effective and robust than other optimizers applied previously to the similarity based alignment problem. We show that QSARs derived using the alignments generated by our algorithm are superior to QSARs derived using the more common alignment of fitting of common structural units. © 1997 by John Wiley & Sons, Inc. J Comput Chem 18 : 1344–1353, 1997