MC(JBW): Simple but smart Monte Carlo algorithm for free energy simulations of multiconformational molecules

Many of the most common molecular simulation methods, including Monte Carlo (MC) and molecular or stochastic dynamics (MD or SD), have significant difficulties in sampling the space of molecular potential energy surfaces characterized by multiple conformational minima and significant energy barriers. In such cases improved sampling can be obtained by special techniques that lower such barriers or somehow direct search steps toward different low energy regions of space. We recently described a hybrid MC/SD algorithm [MC(JBW)/SD] incorporating such a technique that directed MC moves of selected torsion and bond angles toward known low energy regions of conformational space. Exploration of other degrees of freedom was left to the SD part of the hybrid algorithm. In the work described here, we develop a related but simpler simulation algorithm that uses only MC to sample all degrees of freedom (e.g., stretch, bend, and torsion). We term this algorithm MC(JBW). Using simulations on various model potential energy surfaces and on simple molecular systems (n-pentane, n-butane, and cyclohexane), MC(JBW) is shown to generate ensembles of states that are indistinguishable from the canonical ensembles generated by classical Metropolis MC in the limit of very long simulations. We further demonstrate the utility of MC(JBW) by evaluating the room temperature free energy differences between conformers of various substituted cyclohexanes and the larger ring hydrocarbons cycloheptane, cyclooctane, cyclononane, and cyclodecane. The results compare favorably with available experimental data and results from previously reported MC(JBW)/SD conformational free energy calculations. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1736–1745, 1998     

Torsional flexing: Conformational searching of cyclic molecules in biased internal coordinate space

A new stochastic (Monte Carlo) procedure, termed torsional flexing, has been devised for searching the conformational space of cyclic molecules. Torsional flexing causes a local, torsion angle-biased, distortion of a ring bond in a cyclic molecule. Because torsional flexing does not cause large atomic movements, even when it is applied to several bonds simultaneously, subsequent energy minimization generally proceeds rapidly. Nevertheless, the torsional flexing method is prone to generate structures that cross energy barriers so that the structure resulting after energy minimization is frequently a different conformer of the cyclic molecule. Conformational searches on cycloheptadecane, oxobrefeldin A, cyclopenta-L -alanine, and rifamycin SV based upon torsional flexing indicated that torsional flexing is among the best methods yet devised for searching the conformational space of flexible cyclic molecules. © 1993 John Wiley & Sons, Inc.     

Evaluating changes of writhe in computer simulations of supercoiled DNA

We compute changes in the writhe of a polygonal space curve when one of the vertices is displaced. The resulting expressions can be used in simulations of supercoiled DNA. For Brownian dynamics simulations, the expressions can be used to eliminate the explicit twisting degree of freedom. For Monte Carlo simulations, they can be used in fast local moves. Preliminary Monte Carlo simulations using only such fast local moves show that these can be used to efficiently simulate small DNA supercoils.     

Stochastic generator of chemical structure. 4. Building polymeric systems with specified properties

A novel computational technique to generate close‐to‐equilibrium crosslinked polymeric systems is proposed. Compared to the current state‐of‐the‐art equilibration methods, the new technique appears to be faster by several orders of magnitude. The main advantage of the technique is that one can circumvent the bottlenecks in configuration space that inhibit relaxation in molecular dynamics or Monte Carlo simulations. The problem of polymer equilibration described by continuous equations in molecular dynamics is reduced to a discrete representation where solutions are approximated by simple algorithms. In the current study, a series of coarse‐grained, united‐atom, and fully atomistic crosslinked networks has been generated. Network statistics and topology, X‐ray scattering intensities, and elastic properties are tested vs. experimental results and similar models generated using molecular dynamics and Monte Carlo simulations. The results demonstrate the efficiency of this new method for generating large realistic polymeric systems up to 1.4 M atoms. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 580–590, 2001     

Application of combined Monte Carlo and molecular dynamics method to simulation of dipalmitoyl phosphatidylcholine lipid bilayer

We describe a new equilibration procedure for the atomic level simulation of a hydrated lipid bilayer. The procedure consists of alternating molecular dynamics trajectory calculations in a constant surface tension and temperature ensemble with configurational bias Monte Carlo moves to different regions of the configuration space of the bilayer, in a constant volume and temperature ensemble. The procedure is described in detail and is applied to a bilayer of 100 molecules of dipalmitoyl phosphatidylcholine (DPPC) and 3205 water molecules. We find that the hybrid simulation procedure enhances the equilibration of the bilayer as measured by the convergence of the area per molecule and the segmental order parameters, as compared with a simulation using only molecular dynamics (MD). Progress toward equilibration is almost three times as fast in CPU time, compared with a purely MD simulation. Equilibration is complete, as judged by the lack of energy drift in three separate 200-ps runs of continuous MD started from different initial states. Results of the simulation are presented and compared with experimental data and with other recent simulations of DPPC. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1153–1164, 1999     

Optimization of Monte Carlo trial moves for protein simulations

Closed rigid-body rotations of residue segments under bond-angle restraints are simple and effective Monte Carlo moves for searching the conformational space of proteins. The efficiency of these moves is examined here as a function of the number of moving residues and the magnitude of their displacement. It is found that the efficiency of folding and equilibrium simulations can be significantly improved by tailoring the distribution of the number of moving residues to the simulation temperature. In general, simulations exploring compact conformations are more efficient when the average number of moving residues is smaller. It is also demonstrated that the moves do not require additional restrictions on the magnitude of the rotation displacements and perform much better than other rotation moves that do not restrict the bond angles a priori. As an example, these results are applied to the replica exchange method. By assigning distributions that generate a smaller number of moving residues to lower temperature replicas, the simulation times are decreased as long as the higher temperature replicas are effective.     

RNA conformational sampling. I. Single-nucleotide loop closure

We consider the loop-closure problem for nucleic acids and describe an efficient numerical algorithm for closing single-nucleotide loops in nucleic acids. Using six new internal coordinates to represent the nucleotide conformation, which we call the R-representation, the original closure problem with six free torsion angles in each nucleotide can be reduced to one with only four degrees of freedom. Simple numerical techniques have been used to solve the resulting loop-closure equations, and a test of the closure algorithm on a set of RNAs consisting of more than 7000 nucleotides was able to regenerate the native torsion angles in every nucleotide in the test set without exception. We show how the conformational probability density transforms when the original torsion angle representation is mapped onto the new R-representation. We also present statistical evidence showing that the delta and nu(2) torsion angles are coupled, and how this coupling affects the conformation probability density in the R-representation. In addition to the backbone, the same loop-closure algorithm can also be applied to close the ribose ring. The algorithm is freely available at http://tyrosine.use.edu/closure.     

A kinematic view of loop closure

We consider the problem of loop closure, i.e., of finding the ensemble of possible backbone structures of a chain segment of a protein molecule that is geometrically consistent with preceding and following parts of the chain whose structures are given. We reduce this problem of determining the loop conformations of six torsions to finding the real roots of a 16th degree polynomial in one variable, based on the robotics literature on the kinematics of the equivalent rotator linkage in the most general case of oblique rotators. We provide a simple intuitive view and derivation of the polynomial for the case in which each of the three pair of torsional axes has a common point. Our method generalizes previous work on analytical loop closure in that the torsion angles need not be consecutive, and any rigid intervening segments are allowed between the free torsions. Our approach also allows for a small degree of flexibility in the bond angles and the peptide torsion angles; this substantially enlarges the space of solvable configurations as is demonstrated by an application of the method to the modeling of cyclic pentapeptides. We give further applications to two important problems. First, we show that this analytical loop closure algorithm can be efficiently combined with an existing loop-construction algorithm to sample loops longer than three residues. Second, we show that Monte Carlo minimization is made severalfold more efficient by employing the local moves generated by the loop closure algorithm, when applied to the global minimization of an eight-residue loop. Our loop closure algorithm is freely available at http://dillgroup. ucsf.edu/loop_closure/.     

Fast protein structure prediction using Monte Carlo simulations with modal moves

Using normal modes to generate torsion space moves in Monte Carlo simulations of peptides and proteins is not a new idea; nevertheless, despite its power it has not received widespread application. We show that such a "Modal Monte Carlo" approach is an efficient tool for ab initio predictions of small-protein structures. We apply this method to the Trp cage, a 20-residue polypeptide designed to fold rapidly into a structure that includes tertiary contacts, despite its short length. We achieve a high-quality ab initio structure prediction in about 2 orders of magnitude less computation time than state of the art molecular dynamics techniques.     

Ab initio Calculations on the thermal ring closure of acrolein

The electrocyclic reaction which interconverts acrolein and oxetene is investigated by ab initio calculations, using a 4-31G basis set. Like the transition state of the thermal ring closure of cisbutadiene, the transition state of this reaction shows a torsion of about 18 ° around the central bond, in spite of the absence of H-H repulsion. The resulting oxetene is 43 kcal/mol endothermic with regard to acrolein. The calculated activation energies for ring closure and -opening are 3.4 and 1.0 eV. The conformational result of the ring opening reaction is discussed.     

Conformational memories with variable bond angles

Conformational Memories (CM) is a simulated annealing/Monte Carlo method that explores peptide and protein dihedral conformational space completely and efficiently, independent of the original conformation. Here we extend the CM method to include the variation of a randomly chosen bond angle, in addition to the standard variation of two or three randomly chosen dihedral angles, in each Monte Carlo trial of the CM exploratory and biased phases. We test the hypothesis that the inclusion of variable bond angles in CM leads to an improved sampling of conformational space. We compare the results with variable bond angles to CM with no bond angle variation for the following systems: (1) the pentapeptide Met-enkephalin, which is a standard test case for conformational search methods; (2) the proline ring pucker in a 17mer model peptide, (Ala)(8)Pro(Ala)(8); and (3) the conformations of the Ser 7.39 chi(1) in transmembrane helix 7 (TMH7) of the cannabinoid CB1 receptor, a 25-residue system. In each case, analysis of the CM results shows that the inclusion of variable bond angles results in sampling of regions of conformational space that are inaccessible to CM calculations with only variable dihedral angles, and/or a shift in conformational populations from those calculated when variable bond angles are not included. The incorporation of variable bond angles leads to an improved sampling of conformational space without loss of efficiency. Our examples show that this improved sampling leads to better exploration of biologically relevant conformations that have been experimentally validated.     

Grand canonical Monte Carlo simulation of a model colloid-polymer mixture: coexistence line, critical behavior, and interfacial tension

Grand canonical Monte Carlo simulations are used to study phase separation in a simple colloid-polymer model, the so-called Asakura-Oosawa model. To overcome the problem of small acceptance rates of the grand-canonical moves, cluster moves are introduced. Successive umbrella sampling, recently introduced by Virnau and Muller [J. Chem. Phys. 120, 10925 (2004)], is used to access the phase-separated regime. The unmixing binodal and the interfacial tension are measured and compared to theoretical predictions. By means of finite-size scaling, the behavior close to the critical point is also investigated. Close to criticality, we observe substantial deviations from mean-field behavior.     

Assessing protein loop flexibility by hierarchical Monte Carlo sampling

Loop flexibility is often crucial to protein biological function in solution. We report a new Monte Carlo method for generating conformational ensembles for protein loops and cyclic peptides. The approach incorporates the triaxial loop closure method which addresses the inverse kinematic problem for generating backbone move sets that do not break the loop. Sidechains are sampled together with the backbone in a hierarchical way, making it possible to make large moves that cross energy barriers. As an initial application, we apply the method to the flexible loop in triosephosphate isomerase that caps the active site, and demonstrate that the resulting loop ensembles agree well with key observations from previous structural studies. We also demonstrate, with 3 other test cases, the ability to distinguish relatively flexible and rigid loops within the same protein.     

Simulations of rotational isomeric state models for poly(propylene) melts on a high coordination lattice

In the last two years, a method was developed for the Monte Carlo simulation of coarse-grained representations of the chains in polyethylene (PE) melts, under conditions where individual snapshots can be“reverse-mapped” back to continuous space, with all atoms present. In its original form, the symmetry of the torsion potential energy function, E(ϕ) = E(-ϕ), was exploited in the mapping of a coarse-grained version of the rotational isomeric state (RIS) model of PE onto a high coordination lattice. Recently the symmetry restriction was relaxed, so that the simulation could treat isolated RIS poly(propylene) (PP) chains of any stereochemical sequence. Here that simulation of isolated PP chains is extended to PP melts, by introducing the intermolecular interactions required for maintenance of cohesion with the proper cohesive energy. The method is applied to melts of isotactic PP and syndiotactic PP, and satisfactory cohesion is achieved. Satisfactory equilibration of the PP melts requires utilization of a reptation move, in addition to the single bead moves employed previously in simulations of PE melts.     

Monte Carlo sampling algorithm for searching a scale-transformed energy space of polypeptides

A Monte Carlo sampling algorithm for searching a scale-transformed conformational energy space of polypeptides is presented. This algorithm is based on the assumption that energy barriers can be overcome by a uniform sampling of the logarithmically transformed energy space. This algorithm is tested with Met-enkephalin. For comparison, the entropy sampling Monte Carlo (ESMC) simulation is performed. First, the global minimum is easily found by the optimization of a scale-transformed energy space. With a new Monte Carlo sampling, energy barriers of 3000 kcal/mol are frequently overcome, and low-energy conformations are sampled more efficiently than with ESMC simulations. Several thermodynamic quantities are calculated with good accuracy.     

Hybrid Monte Carlo implementation of the Fourier path integral algorithm

This paper formulates a hybrid Monte Carlo implementation of the Fourier path integral (FPI-HMC) approach with partial averaging. Such a hybrid Monte Carlo approach allows one to generate collective moves through configuration space using molecular dynamics while retaining the computational advantages associated with the Fourier path integral Monte Carlo method. In comparison with the earlier Metropolis Monte Carlo implementations of the FPI algorithm, the present HMC method is shown to be significantly more efficient for quantum Lennard-Jones solids and suggests that such algorithms may prove useful for efficient simulations of a range of atomic and molecular systems.     

Conformational analysis of furanose rings with PSEUROT: parametrization for rings possessing the arabino,lyxo, ribo,and xylo stereochemistry and application to arabinofuranosides

The solution conformation of a furanose ring can be assessed through PSEUROT analysis of three-bond (1)H-(1)H coupling constants ((3)J(HH)) of the ring hydrogens. For each coupling constant, PSEUROT requires two parameters, A and B, which are used to translate the H[bond]C[bond]C[bond]H dihedral angle predicted from the (3)J(HH) into an endocyclic torsion angle from which the identity of the conformers can be determined. In this paper, we have used density functional theory methods to generate a family of envelope conformers for methyl furanosides 1-8. From these structures, A and B were calculated for each H[bond]C[bond]C[bond]H fragment. In turn, the values of these parameters for the arabinofuranose ring were used in PSEUROT calculations to determine the conformers populated by monosaccharides 1 and 2 as well as the furanose rings in oligosaccharides 9-15. The results of these analyses are consistent with the low-energy conformers identified from previous computational and X-ray crystallographic studies of 1 and 2.     

Conformational memories and a simulated annealing program that learns: Application to LTB4

The use of computer simulations in all areas of chemistry is growing rapidly because of the powerful insights that they have provided into many interesting phenomena. As investigators continuously examine more sophisticated problems, they need increasingly more powerful tools. Hence, much effort has gone into the development of algorithms which might extend the scope and power of standard dynamic and Monte Carlo techniques. In the Monte Carlo regime, the most common area subject to improvement is the choice of a trial move. In the ordinary case, trial moves are generated uniformly at random. In the extended and hopefully improved case, trial moves are generated randomly but not uniformly. In this article we present a new and totally general method of biased sampling which is applicable to any flexible molecule. In our method, multiple simulated annealing runs are performed to reveal populated and unpopulated regions of the multidimensional conformation space. The second phase of the simulation is done at a fixed temperature with sampling only from populated regions found in the first phase. Because the simulated annealing runs quickly reveal unpopulated regions of the conformation space, the volume of conformation space that needs to be sampled in the second phase of the algorithm is reduced by many orders of magnitude. Additionally, because no energy minimization is used, these populations represent a canonical ensemble which may be used to estimate conformational free energies. © 1995 by John Wiley & Sons, Inc.     

Monte Carlo simulation of polarizable systems: Early rejection scheme for improving the performance of adiabatic nuclear and electronic sampling Monte Carlo simulations

An early rejection scheme for trial moves in adiabatic nuclear and electronic sampling Monte Carlo simulation (ANES-MC) of polarizable intermolecular potential models is presented. The proposed algorithm is based on Swendsen–Wang filter functions for prediction of success or failure of trial moves in Monte Carlo simulations. The goal was to reduce the amount of calculations involved in ANES-MC electronic moves, by foreseeing the success of an attempt before making those moves. The new method was employed in Gibbs ensemble Monte Carlo (GEMC) simulations of the polarizable simple point charge-fluctuating charge (SPC-FQ) model of water. The overall improvement in GEMC depends on the number of swap attempts (transfer molecules between phases) in one Monte Carlo cycle. The proposed method allows this number to increase, enhancing the chemical potential equalization. For a system with 300 SPC-FQ water molecules, for example, the fractions of early rejected transfers were about 0.9998 and 0.9994 at 373 and 423 K, respectively. This means that the transfer moves consume only a very small part of the overall computing effort, making GEMC almost equivalent to a simulation in the canonical ensemble.