Conformation‐Family Monte Carlo (CFMC): An Efficient Computational Method for Identifying the Low‐Energy States of a Macromolecule

A highly efficient method, Conformation‐Family Monte Carlo (CFMC), has been developed for searching the conformational space of a macromolecule and identifying its low‐energy conformations. This method maintains a database of low‐energy conformations that are clustered into families. The conformations in this database are improved iteratively by a Metropolis‐type Monte Carlo procedure, together with energy minimization, in which the search is biased towards investigating the regions of the lowest‐energy families. The CFMC method has the advantages of our earlier potential‐smoothing methods (in that it `coarse‐grains' the conformational space and exploits information about nearby low‐energy states), but avoids their disadvantages (such as the displacement of the global minimum at large smoothings). The CFMC method is applied to a test protein, domain B of Staphylococcal protein A. Independent CFMC runs yielded the same low‐energy families of conformations from random starts, indicating that the thermodynamically relevant conformational space of this protein has been explored thoroughly. The CFMC method is highly efficient, performing as well as or better than competing methods, such as Monte Carlo with minimization, conformational‐space annealing, and the self‐consistent basin‐to‐deformed‐basin method.     

Bond dissocation and conformational energetics of tetrasulfur: a quantum Monte Carlo study

Variational Monte Carlo (VMC) and fixed-node diffusion Monte Carlo (DMC) calculations are performed for S4. The effect of single- and multireference trial functions, as well as choice of orbitals, is investigated for its effect on the quality of the Monte Carlo estimates. Estimates of symmetric (two S2 molecules) and asymmetric (S atom and S3 molecule) bond dissociation are reported. The conformational change of S4 from C2v to D2h defines a double-well potential and is also estimated. Multireference DMC with natural orbitals (DMC/NO) estimates the energy of the conformational change as 1.20(20) kcal/mol; the dissociation of the long S-S single bond is estimated at 21.1(1.3) kcal/mol, and the asymmetric bond energy is estimated as 53.2(2.4) kcal/mol. An estimate of the total atomization energy using multireference DMC/NO gives a value of 219.5(2.2) kcal/mol. The relative quality of result and implications for simplified trial function design are discussed.     

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.     

Kinetic Mechanism of Conformational Switch between Bistable RNA Hairpins

Transitions between the different conformational states play a critical role in many RNA catalytic and regulatory functions. In this study, we use the Kinetic Monte Carlo method to investigate the kinetic mechanism for the conformational switches between bistable RNA hairpins. We find three types of conformational switch pathways for RNA hairpins: refolding after complete unfolding, folding through basepair-exchange pathways and through pseudoknot-assisted pathways, respectively. The result of the competition between the three types of pathways depends mainly on the location of the rate-limiting base stacks (such as the GC base stacks) in the structures. Depending on the structural relationships between the two bistable hairpins, the conformational switch can follow single or multiple dominant pathways. The predicted folding pathways are supported by the activation energy results derived from the Arrhenius plot as well as the NMR spectroscopy data.     

The conformational statistics of polymer chains as obtained from their energy landscape: a conceptual view of the rotational isomeric state approach

In the rotational isomeric state (RIS) approach the conformational statistics of a polymer chain can be conveniently described in terms of the statistical matrix. In this work it is shown from a conceptual point of view to which degree the projection onto a few states, inherently present in any RIS approach, may lead to systematic deviations for estimation of the a priori probabilities. It is shown that these deviations scale with some power of the inverse temperature. Furthermore an algorithm is presented, deriving the statistical matrix from computer simulations. This algorithm is applied to polyethylene. Comparison with the conformational statistics as obtained from Monte Carlo simulations allows the quality of the statistical matrix elements to be estimated. It turns out that the systematic deviation does not hamper practical applications. Questions like the relevance of the anharmonicity of the energy landscape or consideration of more than nearest-neighbor correlations in the RIS approach are discussed.     

Conformational search of peptides and proteins: Monte Carlo minimization with an adaptive bias method applied to the heptapeptide deltorphin

The energy function of a protein consists of a tremendous number of minima. Locating the global energy minimum (GEM) structure, which corresponds approximately to the native structure, is a severe problem in global optimization. Recently we have proposed a conformational search technique based on the Monte Carlo minimization (MCM) method of Li and Scheraga, where trial dihedral angles are not selected at random within the range [-180 degrees,180 degrees ] (as with MCM) but with biased probabilities depending on the increased structure-energy correlations as the GEM is approached during the search. This method, called the Monte Carlo minimization with an adaptive bias (MCMAB), was applied initially to the pentapeptide Leu-enkephalin. Here we study its properties further by applying it to the larger peptide with bulky side chains, deltorphin (H-Tyr-D-Met-Phe-His-Leu-Met-Asp-NH(2)). We find that on average the number of energy minimizations required by MCMAB to locate the GEM for the first time is smaller by a factor of approximately three than the number required by MCM-in accord with results obtained for Leu-enkephalin.     

Conformational transition free energy profiles of an adsorbed, lattice model protein by multicanonical Monte Carlo simulation

Proteins often undergo changes in internal conformation upon interacting with a surface. We investigate the thermodynamics of surface induced conformational change in a lattice model protein using a multicanonical Monte Carlo method. The protein is a linear heteropolymer of 27 segments (of types A and B) confined to a cubic lattice. The segmental order and nearest neighbor contact energies are chosen to yield, in the absence of an adsorbing surface, a unique 3x3x3 folded structure. The surface is a plane of sites interacting either equally with A and B segments (equal affinity surface) or more strongly with the A segments (A affinity surface). We use a multicanonical Monte Carlo algorithm, with configuration bias and jump walking moves, featuring an iteratively updated sampling function that converges to the reciprocal of the density of states 1/Omega(E), E being the potential energy. We find inflection points in the configurational entropy, S(E)=k ln Omega(E), for all but a strongly adsorbing equal affinity surface, indicating the presence of free energy barriers to transition. When protein-surface interactions are weak, the free energy profiles F(E)=E-TS(E) qualitatively resemble those of a protein in the absence of a surface: a free energy barrier separates a folded, lowest energy state from globular, higher energy states. The surface acts in this case to stabilize the globular states relative to the folded state. When the protein surface interactions are stronger, the situation differs markedly: the folded state no longer occurs at the lowest energy and free energy barriers may be absent altogether.     

Order parameters and molecular shapes of phosphonic acid dialkyl esters in a uniaxial hydrophobic environment. A theoretical study

The conformational behaviour of phosphonic acid dialkyl ester (PAE) molecules in a uniaxial hydrophobic environment was studied using the Monte Carlo method with the well-known Metropolis algorithm. The influence of the surrounding molecules on ordering processes of the intrinsic PAE molecules was taken into account by a molecular field. The intramolecular energy was calculated using 6–12 atom-atom and torsion potential functions. Conformations were analysed using torsion angle distributions, segment and bond order parameters. The order parameters of the C-H bonds are compared with experimental results of²H NMR. Further a general method of determining the effective shape of a molecule is presented. The shape is defined by probability distribution for finding atomic coordinates within volume elements during a Monte Carlo run. Thus the asymmetry of a molecule can be visualized. PAE molecules studied show a positive asymmetry in all cases.     

Monte Carlo method for computing density of states and quench probability of potential energy and enthalpy landscapes

The thermodynamics and kinetics of a many-body system can be described in terms of a potential energy landscape in multidimensional configuration space. The partition function of such a landscape can be written in terms of a density of states, which can be computed using a variety of Monte Carlo techniques. In this paper, a new self-consistent Monte Carlo method for computing density of states is described that uses importance sampling and a multiplicative update factor to achieve rapid convergence. The technique is then applied to compute the equilibrium quench probability of the various inherent structures (minima) in the landscape. The quench probability depends on both the potential energy of the inherent structure and the volume of its corresponding basin in configuration space. Finally, the methodology is extended to the isothermal-isobaric ensemble in order to compute inherent structure quench probabilities in an enthalpy landscape.     

Molecular modeling of porous carbons using the hybrid reverse Monte Carlo method

We apply a simulation protocol based on the reverse Monte Carlo (RMC) method, which incorporates an energy constraint, to model porous carbons. This method is called hybrid reverse Monte Carlo (HRMC), since it combines the features of the Monte Carlo and reverse Monte Carlo methods. The use of the energy constraint term helps alleviate the problem of the presence of unrealistic features (such as three- and four-membered carbon rings), reported in previous RMC studies of carbons, and also correctly describes the local environment of carbon atoms. The HRMC protocol is used to develop molecular models of saccharose-based porous carbons in which hydrogen atoms are taken into account explicitly in addition to the carbon atoms. We find that the model reproduces the experimental pair correlation function with good accuracy. The local structure differs from that obtained with a previous model (Pikunic, J.; Clinard, C.; Cohaut, N.; Gubbins, K. E.; Guet, J. M.; Pellenq, R. J.-M.; Rannou, I.; Rouzaud, J. N. Langmuir 2003, 19 (20), 8565). We study the local structure by calculating the nearest neighbor distribution, bond angle distribution, and ring statistics.     

Path integral Monte Carlo method for ab initio calculation

The Feynman path integral method is applied to the many-electron problem. We first give new closure relations in terms of ordinary complex and real numbers, which could be derived from an arbitrary complete set of state vectors. Then, in the path integral form, the partition function of the system and the ensemble average of energy are explicitly expressed in terms of these closure relations. It is impossible to evaluate the path integral by direct numerical integrations because of its huge amount of integration variables. Therefore, we develop an algorithm by the Monte Carlo method with constraints corresponding to the normalization condition of states to calculate the required integral. Finally, the ensemble average of energy for the hydrogen molecule is explicitly evaluated by the quantum Monte Carlo method and results are compared with the result obtained by the ordinary full configuration interaction (CI) method. © 1996 John Wiley & Sons, Inc.     

Efficiency of the multicanonical simulation method as applied to peptides of increasing size: the heptapeptide deltorphin

The advantage of the multicanonical (MUCA) simulation method of Berg and coworkers over the conventional Metropolis method is in its ability to move a system effectively across energy barriers thereby providing results for a wide range of temperatures. However, a MUCA simulation is based on weights (related to the density of states) that should be determined prior to a production run and their calculation is not straightforward. To overcome this difficulty a procedure has been developed by Berg that calculates the MUCA weights automatically. In a previous article (Ya?ar et al. J Comput Chem 2000, 14, 1251-1261) we extended this procedure to continuous systems and applied it successfully to the small pentapeptide Leu-enkephalin. To investigate the performance of the automated MUCA procedure for larger peptides, we apply it here to deltorphin, a linear heptapeptide with bulky side chains (H-Tyr(1)-D-Met(2)-Phe(3)-His(4)-Leu(5)-Met(6)-Asp(7)-NH(2)). As for Leu-enkephalin, deltorphin is modeled in vacuum by the potential energy function ECEPP. MUCA is found to perform well. A weak second peak is seen for the specific heat, which is given a special attention. By minimizing the energy of structures along the trajectory it is found that MUCA provides a good conformational coverage of the low energy region of the molecule. These latter results are compared with conformational coverage obtained by the Monte Carlo minimization method of Li and Scheraga.     

Estimation of protein folding rate from Monte Carlo simulations and entropy capacity

The problem of protein self-organization is one of the most important problems of molecular biology nowadays. Despite the recent success in the understanding of general principles of protein folding, details of this process are yet to be elucidated. Moreover, the prediction of protein folding rates has its own practical value due to the fact that aggregation directly depends on the rate of protein folding. The time of folding has been calculated for 67 proteins with known experimental data at the point of thermodynamic equilibrium between unfolded and native states using a Monte Carlo model where each residue is considered to be either folded as in the native state or completely disordered. The times of folding for 67 proteins which reach the native state within the limit of 10(8) Monte Carlo steps are in a good correlation with the experimentally measured folding rate at the mid-transition point (the correlation coefficient is -0.82). Theoretical consideration of a capillarity model for the process of protein folding demonstrates that the difference in the folding rate for proteins sharing more spherical and less spherical folds is the result of differences in the conformational entropy due to a larger surface of the boundary between folded and unfolded phases in the transition state for proteins with more spherical fold. The capillarity model allows us to predict the folding rate at the same level of correlation as by Monte Carlo simulations. The calculated model entropy capacity (conformational entropy per residue divided by the average contact energy per residue) for 67 proteins correlates by about 78% with the experimentally measured folding rate at the mid-transition point.     

Quantum Monte Carlo study of first-row atoms using transcorrelated variational Monte Carlo trial functions

The effect of using the transcorrelated variational Monte Carlo (TC-VMC) approach to construct a trial function for fixed node diffusion Monte Carlo (DMC) energy calculations has been investigated for the first-row atoms, Li to Ne. The computed energies are compared with fixed node DMC energies obtained using trial functions constructed from Hartree-Fock and density functional levels of theory. Despite major VMC energy improvement with TC-VMC trial functions, no improvement in DMC energy was observed using these trial functions for the first-row atoms studied. The implications of these results on the nodes of the trial wave functions are discussed.     

On the correspondence between classical and quantum mechanics in macromolecular systems: Too much classical chaos

In a recent study we found the classical dynamics of a polyethylene (PE) chain to exhibit low dimensional chaos at temperatures as low as a few Kelvin. These results strongly suggest that classical molecular dynamic simulations in polymer systems can grossly overestimate vibrational motion, which consequently results in disordered structures. In contrast, quantum mechanical calculations using Internal Coordinate Quantum Monte Carlo (an improved method with an initial conjecture for the correct wave function) indicate that the quantum ground state for a three-dimensional model PE chain is far more rigid than determined from molecular dynamics (MD) simulations, even at energies as low as a small fraction of the ground state energy. This result casts uncertainty on the reliability of MD estimates of dynamical or structural quantities relevant to the study of some macromolecular systems.     

Flexible docking simulations: Scaled collective variable Monte Carlo minimization approach using Bezier splines,and comparison with a standard Monte Carlo algorithm

An algorithm for docking a flexible ligand onto a flexible or rigid receptor, using the scaled‐collective‐variables Monte Carlo with energy minimization approach, is presented. Energy minimization is shown to be one of the best techniques for distinguishing between native‐ and nonnative‐generated conformations. Incorporation of this technique into a Monte Carlo procedure enables one to distinguish the native conformation directly during the conformational search. It avoids the generation of a large number of ligand conformers for which more sophisticated energy evaluation tools would have had to be applied to identify the nativelike conformations. The efficiency of the Monte Carlo minimization was greatly improved by incorporating a new grid‐based energy evaluation technique using Bezier splines for which the energy function, as well as all of its derivatives, can be deduced from the values at the grid points. Comparison between our ECEPP/3‐based algorithm and the Monte Carlo algorithm presented elsewhere (Hart, T. N.; Read, R. J. Prot Struct Funct Genet 1992, 13, 206–222) has been made for docking NH₂ D Phe Pro Arg COOH, the noncovalent analog of NH₂ D Phe Pro Arg chloromethylketone (PPACK), onto the active site of human α‐thrombin. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 244–252, 1999     

Self-consistent field equations in the distribution function theory of polymeric liquids

In the distribution function approach to the conformational and thermodynamic properties of polymeric liquids site-site (pair) distribution functions are essential components of the theory. These site-site pair distribution functions are basically mean fields obeying integral equations. In our recent works, a set of self-consistent field equations has been proposed for site-site pair correlation functions which allow us to study conformational and thermodynamic properties of polymeric liquids. In this article, we present a short review of the theory and its applications to a number of aspects of polymeric liquids we have made until now. We also present a self-consistent version of the polymer reference interaction site model where the integral equations for the intramolecular site-site correlation functions are obtained from the Kirkwood hierarchy on the basis of the present theory. The present theory is shown to predict correctly the scaling properties associated with swollen and collapsed polymers in good and poor solvents, respectively. At finite densities, self-consistent solutions of the intra- and intermolecular equations yield the structures and thermodynamics of polymer melts which are favorably compared with Monte Carlo simulation results. Self-consistent theory results are found to be more accurate than the non-self-consistent approaches that use an ideal Gaussian chain conformation distribution function. © 1995 John Wiley & Sons, Inc.     

A version of diffusion Monte Carlo method based on random grids of coherent states. II. Six-dimensional simulation of electronic states of H2

We report a new version of the diffusion Monte Carlo (DMC) method, based on coherent-state quantum mechanics. Randomly selected grids of coherent states in phase space are used to obtain numerical imaginary time solutions of the Schrodinger equation, with an iterative refinement technique to improve the quality of the Monte Carlo grid. Accurate results were obtained, for the appropriately symmetrized two lowest states of the hydrogen molecule, by Monte Carlo sampling and six-dimensional propagation in the full phase space.     

Molecular vibrational energy flow and dilution factors in an anharmonic state space

A fourth-order resonance Hamiltonian is derived from the experimental normal-mode Hamiltonian of SCCl2. The anharmonic vibrational state space constructed from the effective Hamiltonian provides a realistic model for vibrational energy flow from bright states accessible by pulsed laser excitation. We study the experimentally derived distribution PE(sigma) of dilution factors sigma as a function of energy. This distribution characterizes the dynamics in the long-time limit. State space models predict that PE(sigma) should be bimodal, with some states undergoing facile intramolecular vibrational energy redistribution (small sigma), while others at the same total energy remain "protected" (sigma approximately 1). The bimodal distribution is in qualitative agreement with analytical and numerical local density of states models. However, there are fewer states protected from energy flow, and the protected states begin to fragment at higher energy, shifting from sigma approximately 1 to sigma approximately 0.5. We also examine how dilution factors are distributed in the vibrational state space of SCCl2 and how the power law specifying the survival probability of harmonic initial states correlates with the dilution factor distribution of anharmonic initial states.     

The SAAP force field. A simple approach to a new all-atom protein force field by using single amino acid potential (SAAP) functions in various solvents

A simple strategy to compose a new all-atom protein force field (named as the SAAP force field), which utilizes the single amino acid potential (SAAP) functions obtained in various solvents by ab initio molecular orbital calculation applying the isodensity polarizable continuum model (IPCM), is presented. We considered that the total energy function of a protein force field (E(TOTAL)) is divided into three components; a single amino acid potential term (E(SAAP)), an interamino acid nonbonded interaction term (E(INTER)), and a miscellaneous term (E(OTHERS)), which is ignored (or considered to be constant) at the current version of the force field. The E(INTER) term consists of electrostatic interactions (E(ES')) and van der Waals interactions (E(LJ')). Despite simplicity, the SAAP force field implicitly involves the correlation among individual terms of the Lifson's potential function within a single amino acid unit and can treat solvent effects unambiguously by choosing the SAAP function in an appropriate solvent and the dielectric constant (D) of medium. Application of the SAAP force field to the Monte Carlo simulation of For-Ala(2)-NH(2) in vacuo reasonably reproduced the results of the extensive conformational search by ab initio molecular orbital calculation. In addition, the preliminary Monte Carlo simulations for For-Gly(10)-NH(2) and For-Ala(10)-NH(2) showed reversible transitions from the extended to the pseudosecondary structures in water (D = 78.39) as well as in ether (D = 4.335). The result suggested that the new approach is efficient for fast modeling of protein structures in various environments. Decomposition analysis of the total energy function (E(TOTAL)) by using the SAAP force field suggested that conformational propensities of single amino acids (i.e., the E(SAAP) term) may play definitive roles on the topologies of protein secondary structures.