Conformational search using a molecular dynamics–minimization procedure: Applications to clusters of coulombic charges,Lennard–Jones particles,and waters

A new conformational search method, molecular dynamics–minimization (MDM), is proposed, which combines a molecular dynamics sampling strategy with energy minimizations in the search for low-energy molecular structures. This new method is applied to the search for low energy configurations of clusters of coulombic charges on a unit sphere, Lennard–Jones clusters, and water clusters. The MDM method is shown to be efficient in finding the lowest energy configurations of these clusters. A closer comparison of MDM with alternative conformational search methods on Lennard–Jones clusters shows that, although MDM is not as efficient as the Monte Carlo–minimization method in locating the global energy minima, it is more efficient than the diffusion equation method and the method of minimization from randomly generated structures. Given the versatility of the molecular dynamics sampling strategy in comparison to Monte Carlo in treating molecular complexes or molecules in explicit solution, one anticipates that the MDM method could be profitably applied to conformational search problems where the number of degrees of freedom is much greater. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 60–70, 1998     

An exploration of a novel strategy for superposing several flexible molecules

Summary This paper describes a computational strategy for the superposition of a set of flexible molecules. The combinatorial problems of searching conformational space and molecular matching are reduced drastically by the combined use of simulated annealing methods and cluster analysis. For each molecule, the global minimum of the conformational energy is determined by annealing and the search trajectory is retained in a history file. All the significantly different low-energy conformations are extracted by cluster analysis of data in the history file. Each pair of molecules, in each of their significantly different conformations, is then matched by simulated annealing, using the difference-distance matrix as the objective function. A set of match statistics is then obtained, from which the best match taken from all different conformations can be found. The molecules are then superposed either by reference to a base molecule or by a consensus method. This strategy ensures that as wide a range of conformations as possible is considered, but at the same time that the smallest number of significantly different conformations is used. The method has been tested on a set of six angiotensin II antagonists with between 7–11 rotatable bonds.     

Efficiency of simulated annealing for peptides with increasing geometrical restrictions

Simulated annealing (SA) is a popular global minimizer that can conveniently be applied to complex macromolecular systems. Thus, a molecular dynamics or a Monte Carlo simulation starts at high temperature, which is decreased gradually, and the system is expected to reach the low-energy region on the potential energy surface of the molecule. However, in many cases this process is not efficient. Alternatively, the low-energy region can be reached more effectively by minimizing the energy of selected molecular structures generated along the simulation pathway. The efficiency of SA to locate energy-minimized structures within 5 kcal/mol above the global energy minimum is studied as applied to three peptide models with increasing geometrical restrictions: (1) The linear pentapeptide Leu-enkephalin described by the ECEPP potential, (2) a cyclic hexapeptide described by the GROMOS force field energy E_GRO alone, and (3) the same cyclic peptide with E_GRO combined with a restraining potential based on 31 proton–proton restraints obtained from nuclear magnetic resonance (NMR) experiments. The efficiency of SA is compared to that of the Monte Carlo minimization (MCM) method of Li and Scheraga, and to our local torsional deformations (LTD) method for the conformational search of cyclic molecules. The results for the linear peptide show that SA provides a relatively weak guidance towards the most stable energy region; as expected, this guidance increases for the cyclic peptide and the cyclic peptide with NMR restraints. However, in general, MCM and LTD are significantly more efficient than SA as generators of low-energy minimized structures. This suggests that LTD might provide a better search tool than SA in structure determination of protein regions for which a relatively small number of restraints are provided by NMR. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1659–1670, 1999     

Conformational search by potential energy annealing: Algorithm and application to cyclosporin A

Summary A major problem in modelling (biological) macromolecules is the search for low-energy conformations. The complexity of a conformational search problem increases exponentially with the number of degrees of freedom which means that a systematic search can only be performed for very small structures. Here we introduce a new method (PEACS) which has a far better performance than conventional search methods.To show the advantages of PEACS we applied it to the refinement of Cyclosporin A and compared the results with normal molecular dynamics (MD) refinement. The structures obtained with PEACS were lower in energy and agreed with the NMR parameters much better than those obtained with MD. From the results it is further clear that PEACS samples a much larger part of the available conformational space than MD does.     

Toward absolute density of states calculations for proteins

The density of states (DOS), which gives the number of conformations with a particular energy E, is a prerequisite in computing most thermodynamic quantities and in elucidating important biological processes such as the mechanism of protein folding. However, current methods for computing DOS of large systems such as proteins generally yield only the ratios of microstate counts for different energies, which could yield absolute conformation counts if the total number of conformations in phase space is known, thus motivating this work. Here, the total number of energy minima of 50-mer polyalanine, whose size corresponds to naturally occurring small proteins, was estimated under an all-atom potential energy function based on the cumulative distribution function (CDF) of conformational differences to be approximately 10(38). This estimate can place any DOS function, such as the Gaussian DOS distribution in the random energy model, on an absolute scale. Comparing the absolute conformational counts from a Gaussian DOS function with those from the CDF derived from quenched molecular dynamics ensembles shows that the former are far greater than the latter, indicating far fewer low-energy minima actually exist. In addition to showing how CDF and relative DOS calculations can yield absolute DOS for a discrete system, we also show how they can yield absolute DOS for continuous variable systems to a specified atomic variance. In the context of protein folding, knowing this phase-space "volume" of conformations in a DOS function, as well as characteristic transition times, constrains the set of possible folding mechanisms.     

Cages, baskets, ladders, and tubes: conformational studies of polyhedral oligomeric silsesquioxanes

The conformational flexibility of a series of cage, basket, ladder, and tube polyhedral oligomeric silsesquioxanes (POSS) has been examined using the Low Mode:Monte Carlo conformational search method in conjunction with the MM3/GBSA(CHCl3) surface. An ensemble of low energy structures was generated and used to explore the molecular shape and flexibility of each system. The results indicate that, except for the ladder molecule, the incompletely condensed systems that are studied are relatively rigid. Even in cases where the molecule is able to adopt numerous low energy conformations, the overall shape remains cage-like and the conformations differ only by small angles or substituent orientations. The ladder molecule is the most flexible and this ensemble clusters into two families: one that is cage-like and the other that is more open and ladder-like. The conformational flexibilities in the gas and solvent phases, as approximated using the GBSA continuum solvent model, are very similar.     

Conformational analysis and flexibility of carbohydrates using the CICADA approach with MM3

The potential energy hypersurfaces (PES) of several carbohydrate molecules were studied with a new algorithm for conformational searches, CICADA (Channels in Conformational Space Analyzed by Driver Approach) interfaced with the molecular mechanics program MM3(92). The method requires (1) one or a few low-energy conformations as starting points; and (2) designation of the torsion angles important for understanding the conformational behavior of the molecule. The PES is explored by driving separately each selected torsion angle (in both directions) with a concomitant full-geometry optimization at each increment (except for the driven angle). When a minimum has been detected, the molecule is freely optimized, and the minima so detected are then stored if not encountered previously. The new minima serve as starting structures for further explorations. The results from CICADA permit prediction of relative and absolute flexibility and conformational softness for both the entire molecule as well as for individual group rotations and local minima. The carbohydrates analyzed were Me-α-D -glucopyranoside, β-D -GlcNAc(1-2)α-D -Man, and α-D -GalNAc(1-3)[α-L -Fuc(1-2)]Gal-O-Me. All the low-energy conformers along with the transition states and flexibilities features were characterized. CICADA found all minima and low-energy conversion pathways for the disaccharide that were found by a traditional grid search. In contrast to the grid search method, CICADA concentrates mostly on the exploration of the low-energy regions of the PES, thereby saving a significant amount of computational time. The performance of the method opens new routes for exploring conformational space of larger molecules, such as oligosaccharides. © 1995 by John Wiley & Sons, Inc.     

Protein structure predictions by parallel simulated annealing molecular dynamics using genetic crossover

We propose a conformational search method to find a global minimum energy structure for protein systems. The simulated annealing is a powerful method for local conformational search. On the other hand, the genetic crossover can search the global conformational space. Our method incorporates these attractive features of the simulated annealing and genetic crossover. In the previous works, we have been using the Monte Carlo algorithm for simulated annealing. In the present work, we use the molecular dynamics algorithm instead. To examine the effectiveness of our method, we compared our results with those of the normal simulated annealing molecular dynamics simulations by using an α-helical miniprotein. We used genetic two-point crossover here. The conformations, which have lower energy than those obtained from the conventional simulated annealing, were obtained.     

Models for the binding of amiodarone to the thyroid hormone receptor

Summary The antiarrhythmic drug amiodarone has recently been characterized as the first known thyroid hormone antagonist. Its mode of interaction with the thyroid hormone receptor is therefore of interest. A computational analysis of the conformational flexibility of amiodarone using molecular mechanics and the semiempirical molecular orbital method AM1 has been performed. The molecular mechanics studies show that the low-energy conformations of the benzoylbenzofuran portion of amiodarone can be grouped into 4 distinct classes, while the diethylaminoethoxy side chain is extremely flexible. Conformers representative of the 4 low-energy classes were fitted to an extended thyroid hormone receptor model. Four independent modes in which amiodarone could bind to the thyroid hormone receptor site were evaluated.     

Free energy landscape of receptor-mediated cell adhesion

Receptor-mediated cell adhesion plays a critical role in cell migration, proliferation, signaling, and survival. A number of diseases, including cancer, show a strong correlation between integrin activation and metastasis. A better understanding of cell adhesion is highly desirable for not only therapeutic but also a number of tissue engineering applications. While a number of computational models and experimental studies have addressed the issue of cell adhesion to surfaces, no model or theory has adequately addressed cell adhesion at the molecular level. In this paper, the authors present a thermodynamic model that addresses receptor-mediated cell adhesion at the molecular level. By incorporating the entropic, conformational, solvation, and long- and short-range interactive components of receptors and the extracellular matrix molecules, they are able to predict adhesive free energy as a function of a number of key variables such as surface coverage, interaction distance, molecule size, and solvent conditions. Their method allows them to compute the free energy of adhesion in a multicomponent system where they can simultaneously study adhesion receptors and ligands of different sizes, chemical identities, and conformational properties. The authors' results not only provide a fundamental understanding of adhesion at the molecular level but also suggest possible strategies for designing novel biomaterials.     

New optimization method for conformational energy calculations on polypeptides: Conformational space annealing

A new optimization method is presented to search for the global minimum-energy conformations of polypeptides. The method combines essential aspects of the build-up procedure and the genetic algorithm, and it introduces the important concept of “conformational space annealing.” Instead of considering a single conformation, attention is focused on a population of conformations while new conformations are obtained by modifying a “seed conformation.” The annealing is carried out by introducing a distance cutoff, D_cut, which is defined in the conformational space; D_cut effectively divides the whole conformational space of local minima into subdivisions. The value of D_cut is set to a large number at the beginning of the algorithm to cover the whole conformational space, and annealing is achieved by slowly reducing it. Many distinct local minima designed to be distributed as far apart as possible in conformational space are investigated simultaneously. Therefore, the new method finds not only the global minimum-energy conformation but also many other distinct local minima as by-products. The method is tested on Met-enkephalin, a 24-dihedral angle problem. For all 100 independent runs, the accepted global minimum-energy conformation was obtained after about 2600 minimizations on average. © 1997 John Wiley & Sons, Inc. J Comput Chem 18: 1222–1232     

Conformational sampling and ensemble generation by molecular dynamics simulations: 18-Crown-6 as a test case

Using the crown ether 18-crown-6 as a test system, molecular dynamics has been evaluated as a technique for conformational searching and thermodynamic ensemble generation. By running a series of 200 ps and 2 ns simulations, an “optimum” temperature range for conformational searching, i.e., the temperature at which one finds the largest number of low energy structures, was demonstrated to be dependent on the time interval at which one examines the structure. By considering conformational degeneracy and entropy with the rigid rotor harmonic oscillator approximation we have been able to demonstrate that the ensemble generated approaches thermodynamic equilibrium in about 6 ns of simulation. To our knowledge this is the first time this has been demonstrated for a complex organic molecule and it highlights the power and usefulness of molecular dynamics as a method for thermodynamic ensemble generation and conformational searching.     

A fast and Space-efficient boundary element method for computing electrostatic and hydration effects in large molecules

At present, there are two widely used approaches for computing molecular hydration and electrostatic effects within the continuum approximation: the finite difference method, in which the electric potential is directly computed on a cubic grid, and the induced polarization charge or boundary element method, in which an induced charge distribution is first computed on the molecular surface and in which solvation effects are then calculated by reference to the reaction field arising from this induced surface charge. While the induced surface charge approach has a number of advantages over finite differences, especially in the computation of hydration forces and solvent stabilization, the applications of this technique have been largely restricted to small molecules. This is primarily due to the very large system of equations that results when the surface of a macromolecule is discretized into elements small enough to ensure an acceptable level of numerical accuracy within the continuum model. This article describes a new algorithm for implementing boundary element calculations within the continuum model. The essence of our approach is only to compute explicitly those interactions between surface elements that are relatively close together and to approximate long-range interactions by grid-based multipole expansion. The resulting system of equations has a relatively sparse coefficient matrix and requires disk storage that increases linearly with molecular surface area. The technique has numerous applications in the analysis of solvation effects in large molecules, especially in the area of conformational analysis, where it is critical to accurately estimate the global hydration energy for the entire structure. © 1996 by John Wiley & Sons, Inc.     

Conformational searching methods for small molecules. II. Genetic algorithm approach

We demonstrate the use of a genetic algorithm (GA) search procedure for finding low-energy conformations of small to medium organic molecules (1–12 rotatable bonds). GAS are in a class of biologically motivated optimization methods that evolve a population of individuals where individuals who are more “fit” have a higher probability of surviving into subsequent generations. Here, an individual is a conformation of a given molecule and the fitness is the molecule's conformational energy. In the course of a simulated evolution, the population produces conformations having increasingly lower energy. We test the GA method on a suite of 72 molecules and compare the performance against the CSEARCH algorithm in Sybyl. For molecules with more than eight rotatable bonds, the GA method is more efficient computationally and as the number of rotatable bonds increases the relative efficiency of the GA method grows. The GA method also found energies equal to or lower than the energy of the relaxed crystal structure in the large majority of cases. © John Wiley & Sons, Inc.     

Efficiency of tabu‐search‐based conformational search algorithms

Efficient conformational search or sampling approaches play an integral role in molecular modeling, leading to a strong demand for even faster and more reliable conformer search algorithms. This article compares the efficiency of a molecular dynamics method, a simulated annealing method, and the basin hopping (BH) approach (which are widely used in this field) with a previously suggested tabu‐search‐based approach called gradient only tabu search (GOTS). The study emphasizes the success of the GOTS procedure and, more importantly, shows that an approach which combines BH and GOTS outperforms the single methods in efficiency and speed. We also show that ring structures built by a hydrogen bond are useful as starting points for conformational search investigations of peptides and organic ligands with biological activities, especially in structures that contain multiple rings. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

Assessment of Two Surface Monte Carlo (TSMC) method to find stationary points of (H2O)15 and (H2O)20 clusters

The Two Surface Monte Carlo (TSMC) technique reduces computational cost by using a computationally cheap biasing potential, which guides the molecular system to explore the potential energy surface of interest. It was shown earlier that the Effective Fragment Potential (EFP) can be a good choice for this biasing potential (Bandyopadhyay, J Chem Phys 122:091102, 2005) when the potential energy surface of interest is quantum mechanical. This may help in expanding the applicability of TSMC, since finding a good biasing potential is a major challenge. In the present work, the viability of TSMC method in finding stationary points of large molecular system is investigated using EFP as the biasing potential and RHF theory as the potential of interest. TSMC is applied to find the stationary points of water clusters of size 15 and 20. A semi-automated method starting from random geometries, without using any chemical intuition, found several stationary points. The simulated annealing method was used to refine the structures obtained from TSMC. Among the several low-energy structures obtained for 15 water cluster, one minimum, about 1 kcal/mol higher than the global minimum, was found. However, for 20 water cluster, no structure very close to the global minimum was obtained. Several strategies, learned from the experience of the present work, are discussed for improving the TSMC method, including the acceptance between the two energy surfaces. Contribution to the Mark S. Gordon 65th Birthday Festschrift Issue.     

Structures of protonated arginine dimer and bradykinin investigated by density functional theory: further support for stable gas-phase salt bridges

The gas-phase structures and energetics of both protonated arginine dimer and protonated bradykinin were investigated using a combination of molecular mechanics with conformational searching to identify candidate low-energy structures, and density functional theory for subsequent minimization and energy calculations. For protonated arginine dimer, a good correlation (R = 0.88) was obtained between the molecular mechanics and EDF1 6-31+G* energies, indicating that mechanics with MMFF is suitable for finding low-energy conformers. For this ion, the salt-bridge or ion-zwitterion form was found to be 5.7 and 7.2 kcal/mol more stable than the simple protonated or ion-molecule form at the EDF1 6-31++G** and B3LYP 6-311++G** levels. For bradykinin, the correlation between the molecular mechanics and DFT energies was poor (R = 0.28), indicating that many low-energy structures are likely passed over in the mechanics conformational searching. This result suggests that structures of this larger peptide ion obtained using mechanics calculations alone are not necessarily reliable. The lowest energy structure of the salt-bridge form of bradykinin is 10.6 kcal/mol lower in energy (EDF1) than the lowest energy simple protonated form at the 6-311G* level. Similarly, the average energy of all salt-bridge structures investigated is 13.6 kcal/mol lower than the average of all the protonated forms investigated. To the extent that a sufficient number of structures are investigated, these results provide some additional support for the salt-bridge form of bradykinin in the gas phase.     

Accelerating rare reactive events by means of a finite electronic temperature

The range of chemical problems that are directly accessible to first-principles molecular dynamics simulations based on density functional theory is extended with a novel method apt to accelerate rare reactive events. The introduction of a finite electronic temperature within the Mermin formalism leads to a lowering of chemical activation barriers and thus to an exponential enhancement of the rate at which these reactions are observed during a first-principles molecular dynamics simulation. The method presented here makes direct use of the intrinsic chemical information encoded in the electronic structure, and is therefore able to lower selectively chemically relevant activation energies even in systems where many competing low-energy pathways for conformational transitions or diffusive motions are present. The performance of this new approach is demonstrated for a series of prototypical chemical reactions in gas and in condensed phase. A typical acceleration that can be achieved is, for example, a factor of 10(5) for the cis-trans isomerization of peroxynitrous acid in aqueous solution at room temperature.     

Photochemistry and vibrational spectra of matrix-isolated 5-ethoxy-1-phenyl-1H-tetrazole

A combined matrix isolation FT-IR and theoretical DFT(B3LYP)/6-311++G(d,p) study of the molecular structure and photochemistry of 5-ethoxy-1-phenyl-1H-tetrazole (5EPT) was performed. A new method of synthesis of the compound is described. Calculations show three minima, very close in energy and separated by low-energy barriers (less than 4 kJ mol-1), in the ground-state potential energy profile of the molecule. The method of matrix isolation enabled the reduction of the number of populated conformational states in the experiment at low temperature due to the effect known as conformational cooling. As a result, the spectrum of the as-deposited matrix of 5EPT closely matches that of the most stable conformer predicted theoretically, pointing to the existence of only this conformer in the low-temperature matrixes. In this structure, the dihedral angle between the two rings, phenyl and tetrazole, is ca. 30 degrees, whereas the ethyl group stays nearly in the plane of the tetrazole ring and is as far as possible from the phenyl group. In situ UV irradiation (lambda > 235 nm) of the matrix-isolated 5EPT induced unimolecular decomposition, which led mainly to production of ethylcyanate and phenylazide, this later compound further reacting to yield, as final product, 1-aza-1,2,4,6-cycloheptatetraene. Anti-aromatic 3-ethoxy-1-phenyl-1H-diazirene was also observed experimentally as minor photoproduct, resulting from direct extrusion of molecular nitrogen from 5EPT. This species has not been described before and is now characterized by infrared spectroscopy for the first time.