Finite Temperature String Method for the Study of Rare Events†

Path finding on high-dimensional free energy landscapes

We present a method for determining the average transition path and the free energy along this path in the space of selected collective variables. The formalism is based upon a history-dependent bias along a flexible path variable within the metadynamics framework but with a trivial scaling of the cost with the number of collective variables. Controlling the sampling of the orthogonal modes recovers the average path and the minimum free energy path as the limiting cases. The method is applied to resolve the path and the free energy of a conformational transition in alanine dipeptide.

Mapping L1 ligase ribozyme conformational switch

L1 ligase (L1L) molecular switch is an in vitro optimized synthetic allosteric ribozyme that catalyzes the regioselective formation of a 5'-to-3' phosphodiester bond, a reaction for which there is no known naturally occurring RNA catalyst. L1L serves as a proof of principle that RNA can catalyze a critical reaction for prebiotic RNA self-replication according to the RNA world hypothesis. L1L crystal structure captures two distinct conformations that differ by a reorientation of one of the stems by around 80Å and are presumed to correspond to the active and inactive state, respectively. It is of great interest to understand the nature of these two states in solution and the pathway for their interconversion. In this study, we use explicit solvent molecular simulation together with a novel enhanced sampling method that utilizes concepts from network theory to map out the conformational transition between active and inactive states of L1L. We find that the overall switching mechanism can be described as a three-state/two-step process. The first step involves a large-amplitude swing that reorients stem C. The second step involves the allosteric activation of the catalytic site through distant contacts with stem C. Using a conformational space network representation of the L1L switch transition, it is shown that the connection between the three states follows different topographical patterns: the stem C swing step passes through a narrow region of the conformational space network, whereas the allosteric activation step covers a much wider region and a more diverse set of pathways through the network.

Minimum free energy path of ligand-induced transition in adenylate kinase

Large-scale conformational changes in proteins involve barrier-crossing transitions on the complex free energy surfaces of high-dimensional space. Such rare events cannot be efficiently captured by conventional molecular dynamics simulations. Here we show that, by combining the on-the-fly string method and the multi-state Bennett acceptance ratio (MBAR) method, the free energy profile of a conformational transition pathway in Escherichia coli adenylate kinase can be characterized in a high-dimensional space. The minimum free energy paths of the conformational transitions in adenylate kinase were explored by the on-the-fly string method in 20-dimensional space spanned by the 20 largest-amplitude principal modes, and the free energy and various kinds of average physical quantities along the pathways were successfully evaluated by the MBAR method. The influence of ligand binding on the pathways was characterized in terms of rigid-body motions of the lid-shaped ATP-binding domain (LID) and the AMP-binding (AMPbd) domains. It was found that the LID domain was able to partially close without the ligand, while the closure of the AMPbd domain required the ligand binding. The transition state ensemble of the ligand bound form was identified as those structures characterized by highly specific binding of the ligand to the AMPbd domain, and was validated by unrestrained MD simulations. It was also found that complete closure of the LID domain required the dehydration of solvents around the P-loop. These findings suggest that the interplay of the two different types of domain motion is an essential feature in the conformational transition of the enzyme.

Transition path from two apposed membranes to a stalk obtained by a combination of particle simulations and string method

The formation of an hourglass-shaped passage (stalk) connecting two apposed membranes is an essential initial step in membrane fusion. The most probable transition path from two separate membranes to a stalk, i.e., the minimum free-energy path (MFEP), is constructed using a combination of particle simulations and string method. For the reversible transition path in the coarse-grained membrane model, a collective order parameter, m, can be identified as the local difference of hydrophilic and hydrophobic densities. In particle simulations, the free energy F[m] as a functional of m is not readily available. This difficulty is overcome by an equation-free approach, where the morphology and the excess free energy along the MFEP are obtained by an on-the-fly string method. The transition state is confirmed by diagonalization of order-parameter fluctuations and by the probability of reaching either stalk or bilayer morphology from different positions along the MFEP.

Transition path sampling algorithm for discrete many-body systems

We propose a Monte Carlo method for efficiently sampling trajectories with fixed initial and final conditions in a system with discrete degrees of freedom. The method can be applied to any stochastic process with local interactions, including systems that are out of equilibrium. We combine the proposed path sampling algorithm with thermodynamic integration to calculate transition rates. We demonstrate our method on the well-studied two-dimensional Ising model with periodic boundary conditions, and show agreement with other results for both large and small system sizes. The method scales well with the system size, allowing one to simulate systems with many degrees of freedom, and providing complementary information with respect to other algorithms.

Atomistic molecular simulations of protein folding

Theory and experiment have provided answers to many of the fundamental questions of protein folding; a remaining challenge is an accurate, high-resolution picture of folding mechanism. Atomistic molecular simulations with explicit solvent are the most promising method for providing this information, by accounting more directly for the physical interactions that stabilize proteins. Although simulations of folding with such force fields are extremely challenging, they have become feasible as a result of recent advances in computational power, accuracy of the energy functions or 'force fields', and methods for improving sampling of folding events. I review the recent progress in these areas, and highlight future challenges and questions that we may hope to address with these methods. I also attempt to place atomistic models into the context of the energy landscape view of protein folding, and coarse-grained simulations.

Multiple replica repulsion technique for efficient conformational sampling of biological systems

Here, we propose a technique for sampling complex molecular systems with many degrees of freedom. The technique, termed "multiple replica repulsion" (MRR), does not suffer from poor scaling with the number of degrees of freedom associated with common replica exchange procedures and does not require sampling at high temperatures. The algorithm involves creation of multiple copies (replicas) of the system, which interact with one another through a repulsive potential that can be applied to the system as a whole or to portions of it. The proposed scheme prevents oversampling of the most populated states and provides accurate descriptions of conformational perturbations typically associated with sampling ground-state energy wells. The performance of MRR is illustrated for three systems of increasing complexity. A two-dimensional toy potential surface is used to probe the sampling efficiency as a function of key parameters of the procedure. MRR simulations of the Met-enkephalin pentapeptide, and the 76-residue protein ubiquitin, performed in presence of explicit water molecules and totaling 32 ns each, investigate the ability of MRR to characterize the conformational landscape of the peptide, and the protein native basin, respectively. Results obtained for the enkephalin peptide reflect more closely the extensive conformational flexibility of this peptide than previously reported simulations. Those obtained for ubiquitin show that conformational ensembles sampled by MRR largely encompass structural fluctuations relevant to biological recognition, which occur on the microsecond timescale, or are observed in crystal structures of ubiquitin complexes with other proteins. MRR thus emerges as a very promising simple and versatile technique for modeling the structural plasticity of complex biological systems.

Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation

Metropolis Monte Carlo simulation is a powerful tool for studying the equilibrium properties of matter. In complex condensed-phase systems, however, it is difficult to design Monte Carlo moves with high acceptance probabilities that also rapidly sample uncorrelated configurations. Here, we introduce a new class of moves based on nonequilibrium dynamics: Candidate configurations are generated through a finite-time process in which a system is actively driven out of equilibrium, and accepted with criteria that preserve the equilibrium distribution. The acceptance rule is similar to the Metropolis acceptance probability, but related to the nonequilibrium work rather than the instantaneous energy difference. Our method is applicable to sampling from both a single thermodynamic state or a mixture of thermodynamic states, and allows both coordinates and thermodynamic parameters to be driven in nonequilibrium proposals. Whereas generating finite-time switching trajectories incurs an additional cost, driving some degrees of freedom while allowing others to evolve naturally can lead to large enhancements in acceptance probabilities, greatly reducing structural correlation times. Using nonequilibrium driven processes vastly expands the repertoire of useful Monte Carlo proposals in simulations of dense solvated systems.

Catalytic mechanism of RNA backbone cleavage by ribonuclease H from quantum mechanics/molecular mechanics simulations

We use quantum mechanics/molecular mechanics simulations to study the cleavage of the ribonucleic acid (RNA) backbone catalyzed by ribonuclease H. This protein is a prototypical member of a large family of enzymes that use two-metal catalysis to process nucleic acids. By combining Hamiltonian replica exchange with a finite-temperature string method, we calculate the free energy surface underlying the RNA-cleavage reaction and characterize its mechanism. We find that the reaction proceeds in two steps. In a first step, catalyzed primarily by magnesium ion A and its ligands, a water molecule attacks the scissile phosphate. Consistent with thiol-substitution experiments, a water proton is transferred to the downstream phosphate group. The transient phosphorane formed as a result of this nucleophilic attack decays by breaking the bond between the phosphate and the ribose oxygen. In the resulting intermediate, the dissociated but unprotonated leaving group forms an alkoxide coordinated to magnesium ion B. In a second step, the reaction is completed by protonation of the leaving group, with a neutral Asp132 as a likely proton donor. The overall reaction barrier of ∼15 kcal mol(-1), encountered in the first step, together with the cost of protonating Asp132, is consistent with the slow measured rate of ∼1-100/min. The two-step mechanism is also consistent with the bell-shaped pH dependence of the reaction rate. The nonmonotonic relative motion of the magnesium ions along the reaction pathway agrees with X-ray crystal structures. Proton-transfer reactions and changes in the metal ion coordination emerge as central factors in the RNA-cleavage reaction.

Free energy of conformational transition paths in biomolecules: the string method and its application to myosin VI

A set of techniques developed under the umbrella of the string method is used in combination with all-atom molecular dynamics simulations to analyze the conformation change between the prepowerstroke (PPS) and rigor (R) structures of the converter domain of myosin VI. The challenges specific to the application of these techniques to such a large and complex biomolecule are addressed in detail. These challenges include (i) identifying a proper set of collective variables to apply the string method, (ii) finding a suitable initial string, (iii) obtaining converged profiles of the free energy along the transition path, (iv) validating and interpreting the free energy profiles, and (v) computing the mean first passage time of the transition. A detailed description of the PPS↔R transition in the converter domain of myosin VI is obtained, including the transition path, the free energy along the path, and the rates of interconversion. The methodology developed here is expected to be useful more generally in studies of conformational transitions in complex biomolecules.

Entropic effect on the rate of dislocation nucleation

Dislocation nucleation is essential to our understanding of plastic deformation, ductility, and mechanical strength of crystalline materials. Molecular dynamics simulation has played an important role in uncovering the fundamental mechanisms of dislocation nucleation, but its limited timescale remains a significant challenge for studying nucleation at experimentally relevant conditions. Here we show that dislocation nucleation rates can be accurately predicted over a wide range of conditions by determining the activation free energy from umbrella sampling. Our data reveal very large activation entropies, which contribute a multiplicative factor of many orders of magnitude to the nucleation rate. The activation entropy at constant strain is caused by thermal expansion, with negligible contribution from the vibrational entropy. The activation entropy at constant stress is significantly larger than that at constant strain, as a result of thermal softening. The large activation entropies are caused by anharmonic effects, showing the limitations of the harmonic approximation widely used for rate estimation in solids. Similar behaviors are expected to occur in other nucleation processes in solids.

Mobility of a 5|7 defect in carbon nanotubes

The movement of a 5|7 defect in the plastic deformations of a carbon nanotube (CNT) plays the role of dislocation glide in the plastic deformations of metals. This work is concerned with the atomic shift mechanism of the 5|7 defect and the energy barrier when the CNT is subjected to tensile loading. Action-derived molecular dynamics (ADMD) is applied to find the minimum energy path and the energy barrier. It is found that the tensile loads make it easy for the 5|7 defect to glide, and lower the energy barrier. The minimum load level that makes a 5|7 defect glide freely with no energy barrier in the presence of an adatom is obtained.

A computational study of the mechanism of the selective crystallization of α- and β-glycine from water and methanol-water mixture

Understanding the control of polymorphism in organic crystals is of paramount importance to the pharmaceutical, chemical, and food industries. In this work, we investigated two mechanisms described in the literature about the selective crystallization of α- and β-glycine from water and from mixtures of water and methanol using molecular simulations. The link hypothesis (J. Phys. Chem. B 2008, 112, 7794; Cryst. Growth Des. 2006, 6, 1788; J. Inclusion Phenom. Mol. Recognit. Chem. 1990, 8, 395; J. Am. Chem. Soc. 1986, 108, 5871.), which tries to relate the structure of the polymorph obtained from crystallization to the structure of the prenucleation aggregates in the solutions, says the abundance of glycine cyclic dimers in aqueous solutions leads to the crystallization of α-glycine, the polymorph using cyclic dimers as the packing units. This hypothesis was studied first. We revisited the self-assembly of glycine molecules in solution using molecular dynamics to address the debate (Phys. Rev. Lett. 2007, 99, 115702; J. Phys. Chem. B 2008, 112, 7280; J. Am. Chem. Soc. 2008, 130, 13973.) about which is the dominating species in the glycine aqueous solutions and whether there is a link between the solution chemistry and the polymorphic outcome of crystallization. The structures of the glycine clusters were characterized using a structural parameter called cyclic dimer fraction. The glycine clusters in methanol-water mixtures have higher cyclic dimer compositions than those in the pure aqueous solutions. Moreover, the glycine open-chain dimer is more stable than the cyclic dimer regardless of the presence of methanol. All these suggest that the link hypothesis does not work for the polymorphic system of glycine, and the selective crystallization of α- and β-glycine from water and methanol-water mixture, respectively, is not due to the abundance of glycine aggregates in the solution phase with a similar structure to the crystallizing solid form. The hypothesis of the methanol inhibition on the growth of α-glycine {010} and {010} faces, proposed by Weissbuch (Angew. Chem., Int. Ed. 2005, 44, 3226.), was also studied. The interfaces between the {010} and {010} faces of both crystal forms (α and β) and both solvents (water and methanol-water 3:7 mixture) were studied using molecular simulation. No strong binding of methanol onto the {010} and {010} faces of both crystal forms was observed, and the addition of methanol dilutes the crystal-solvent interactions on all faces. Therefore, the selective crystallization of β and α-glycine with and without methanol does not follow either of the two mechanisms in the literature.

Simulations of allosteric transitions

Allosteric transitions are one of the subtler mechanisms used by nature to fine tune protein activity. Effector binding to a specific site on the protein surface induces significant activity change, and initiates a conformational transition that frequently includes domain motions and is very large. From a theoretical and biophysical perspective two problems are particularly intriguing. The first is the way in which a launching signal, which is spatially confined and includes only a few interacting atoms, is propagated to a large-scale conformational transition we frequently see in allosteric transitions. Hence, there is the question of how a small perturbation is magnified to yield motions of thousands of atoms. The second puzzle is of focus, coherence, and efficiency. The impact of the binding of the effector is sometimes extended over tens of angstroms. How the signal is transmitted and kept significant over such large distances in the 'noisy' molecular environment is another major direction of investigation. In the present review we examined different theoretical and computational attempts to solve the questions.

Overcoming barriers in trajectory space: mechanism and kinetics of rare events via Wang-Landau enhanced transition path sampling

Within the framework of transition path sampling (TPS), activation energies can be computed as path ensemble averages without a priori information about the reaction mechanism [C. Dellago and P. G. Bolhuis, Mol. Simul. 30, 795 (2004)]. Activation energies computed for different conditions can then be used to determine by numerical integration the rate constant for a system of interest from the rate constant known for a reference system. However, in systems with complex potential energy surfaces, multiple reaction pathways may exist making ergodic sampling of trajectory space difficult. Here, we present a combination of TPS with the Wang-Landau (WL) flat-histogram algorithm for an efficient sampling of the transition path ensemble. This method, denoted by WL-TPS, has the advantage that from one single simulation, activation energies at different temperatures can be determined even for systems with multiple reaction mechanisms. The proposed methodology for rate constant calculations does not require the knowledge of the reaction coordinate and is generally applicable to Arrhenius and non-Arrhenius processes. We illustrate the applicability of this technique by studying a two-dimensional toy system consisting of a triatomic molecule immersed in a fluid of repulsive soft disks. We also provide an expression for the calculation of activation volumes from path averages such that the pressure dependence of the rate constant can be obtained by numerical integration.

Generating stereochemically acceptable protein pathways

We describe a new method for rapidly generating stereochemically acceptable pathways in proteins. The method, called geometric targeting, is publicly available at the webserver http://pathways.asu.edu, and includes tools for visualization of the pathway and creating movie files for use in presentations. The user submits an initial structure and a target structure, and a pathway between the two input states is generated automatically. Besides visualization, the structural quality of the pathways makes them useful as input pathways into pathway refinement techniques and further computations. The approach in geometric targeting is to gradually change the system's RMSD relative to the target structure while enforcing a set of geometric constraints. The generated pathways are not minimum free energy pathways, but they are geometrically plausible pathways that maintain good covalent bond distances and angles, keep backbone dihedral angles in allowed Ramachandran regions, avoid eclipsed side-chain torsion angles, avoid non-bonded overlap, and maintain a set of hydrogen bonds and hydrophobic contacts. Resulting pathways for over 20 proteins featuring a wide variety of conformational changes are reported here, including the very large GroEL complex.

Slow conformational motions that favor sub-picosecond motions important for catalysis

It has been accepted for many years that functionally important motions are crucial to binding properties of ligands in such molecules as hemoglobin and myoglobin. In enzymatic reactions, theory and now experiment are beginning to confirm the importance of motions on a fast (ps) time scale in the chemical step of the catalytic process. What is missing is a clear physical picture of how slow conformational fluctuations are related to the fast motions that have been identified as crucial. This paper presents a theoretical analysis of this issue for human heart lactate dehydrogenase. We will examine how slow conformational motions bring the system to conformations that are distinguished as catalytically competent because they favor specific fast motions.

Separating forward and backward pathways in nonequilibrium umbrella sampling

Umbrella sampling enforces uniform sampling of steady-state distributions that are functions of arbitrary numbers of order parameters. The key to applying such methods to nonequilibrium processes is the accumulation of fluxes between regions. A significant difference between microscopically reversible and irreversible systems is that, in the latter case, the transition path ensemble for a reaction can be significantly different for "forward" and "backward" trajectories. Here, we show how to separately treat forward and backward pathways in nonequilibrium umbrella sampling simulations by working in an extended space. In this extended space, the exact rate (for equilibrium or nonequilibrium processes) can be calculated "for free" as a flux in phase space. We compare the efficiency of this rate calculation with forward flux sampling for a two-dimensional potential and show that nonequilibrium umbrella sampling is more efficient when an intermediate is present. We show that this technique can also be used to describe steady-state limit cycles by examining a simulation of circadian oscillations. We obtain the path of the limit cycle in a space of 22 order parameters, as well as the oscillation period. The relation of our method to others is discussed.

Onsager-Machlup action-based path sampling and its combination with replica exchange for diffusive and multiple pathways

For sampling multiple pathways in a rugged energy landscape, we propose a novel action-based path sampling method using the Onsager-Machlup action functional. Inspired by the Fourier-path integral simulation of a quantum mechanical system, a path in Cartesian space is transformed into that in Fourier space, and an overdamped Langevin equation is derived for the Fourier components to achieve a canonical ensemble of the path at a finite temperature. To avoid "path trapping" around an initially guessed path, the path sampling method is further combined with a powerful sampling technique, the replica exchange method. The principle and algorithm of our method is numerically demonstrated for a model two-dimensional system with a bifurcated potential landscape. The results are compared with those of conventional transition path sampling and the equilibrium theory, and the error due to path discretization is also discussed.

On the Origin of the Chemical Barrier and Tunneling in Enzymes

This paper presents both a review of some recent results from our group and experimental groups, and some new theoretical results all of which are helping to form a more physically rigorous picture of the process of enzymatic catalysis. A common classical picture of enzymatic catalysis is the transition state tight binding model. Schwartz and Schramm1 have recently argued from both theoretical and experimental results that this picture is incorrect. We now investigate what the nature of barriers might be in enzymatic reactions, and what this viewpoint might imply for tunneling in a hydrogen transfer enzyme. For lactate dehydrogenase we conclude that the enzymes role in catalysis is at least partially to hunt through configuration space for those configurations that minimize chemical free energy barriers. Those configurations do not seem to be stable basins on the free energy surface, and in fact the overall free energy barrier to reaction may well largely be due to this stochastic hunt - both probabilistically and energetically. We suggest further computations to test this hypothesis.

Exploration of effective potential landscapes using coarse reverse integration

We describe a reverse integration approach for the exploration of low-dimensional effective potential landscapes. Coarse reverse integration initialized on a ring of coarse states enables efficient navigation on the landscape terrain: Escape from local effective potential wells, detection of saddle points, and identification of significant transition paths between wells. We consider several distinct ring evolution modes: Backward stepping in time, solution arc length, and effective potential. The performance of these approaches is illustrated for a deterministic problem where the energy landscape is known explicitly. Reverse ring integration is then applied to noisy problems where the ring integration routine serves as an outer wrapper around a forward-in-time inner simulator. Two versions of such inner simulators are considered: A Gillespie-type stochastic simulator and a molecular dynamics simulator. In these "equation-free" computational illustrations, estimation techniques are applied to the results of short bursts of inner simulation to obtain the unavailable (in closed-form) quantities (local drift and diffusion coefficient estimates) required for reverse ring integration; this naturally leads to approximations of the effective landscape.

Forward flux sampling for rare event simulations

Rare events are ubiquitous in many different fields, yet they are notoriously difficult to simulate because few, if any, events are observed in a conventional simulation run. Over the past several decades, specialized simulation methods have been developed to overcome this problem. We review one recently developed class of such methods, known as forward flux sampling. Forward flux sampling uses a series of interfaces between the initial and final states to calculate rate constants and generate transition paths for rare events in equilibrium or nonequilibrium systems with stochastic dynamics. This review draws together a number of recent advances, summarizes several applications of the method and highlights challenges that remain to be overcome.

Minimum action method for the Kardar-Parisi-Zhang equation

We apply a numerical minimum action method derived from the Wentzell-Freidlin theory of large deviations to the Kardar-Parisi-Zhang equation for the height profile of a growing interface. In one dimension we find that the transition pathway between different height configurations is determined by the nucleation and subsequent propagation of facets or steps, corresponding to moving domain walls or growth modes in the underlying noise-driven Burgers equation. This transition scenario is in accordance with recent analytical studies of the one-dimensional Kardar-Parisi-Zhang equation in the asymptotic weak noise limit. We also briefly discuss transitions in two dimensions.

Some recent techniques for free energy calculations

A few recent techniques to calculate free energies in the context of molecular dynamics simulations are discussed: temperature-accelerated molecular dynamics, which is a method to explore fast the important regions in the free energy landscape associated with a set of continuous collective variables without having to know where these regions are beforehand; the single sweep method, which is a variational method to interpolate the free energy globally given a set of mean forces (i.e., a set of gradients of the free energy) calculated at specific points, or centers, on the free energy landscape; and a Voronoi-based free energy method for the calculation of the free energy of the Voronoi tessellation associated with a set of centers. We also discuss how this last technique can be used in conjunction with the string method, and how kinetic information such as reaction rates can be calculated by milestoning using the edges of a Voronoi tessellation as milestones.

On-the-path random walk sampling for efficient optimization of minimum free-energy path

A sampling method is proposed for the efficient acquisition of minimum free-energy path (MFEP). Here, the MFEP optimization is realized based on the sampling via single on-the-path random walk simulation. The present strategy naturally ensures the on-the-path structural continuity so that the MFEP optimization can be simplified and its efficiency can be improved as well.

Free energy and kinetics of conformational transitions from Voronoi tessellated milestoning with restraining potentials

Milestoning is a method aimed at reconstructing the statistical properties of the long-time dynamics of a system by exploiting the crossing statistics of a set of hypersurfaces, called the "milestones", placed along the reaction coordinate [Faradjian and Elber, J. Chem. Phys.2004, 120, 10880]. Recently, Vanden-Eijnden and Venturoli [J. Chem. Phys.2009, 130, 194101] showed that when a complete Voronoi tessellation of the configurational space is available, milestoning can be reformulated to utilise the statistics from a series of independent simulations, each confined within a single cell via strict reflections at the boundaries. As a byproduct, this "Voronoi tessellated milestoning" method also permits to compute the free energy of the tessellation. Here, the method is extended to support the usage of differentiable restraining potentials to confine the trajectories within each cell.

Transition path sampling and forward flux sampling. Applications to biological systems

The last decade has seen a rapid growth in the number of simulation methods and applications dealing with the sampling of transition pathways of rare nanoscale events. Such studies are crucial, for example, for understanding the mechanism and kinetics of conformational transitions and enzymatic events associated with the function of biomolecules. In this review, a broad account of transition path sampling approaches is provided, starting from the general concepts, progressing to the specific principles that underlie some of the most important methods, and eventually singling out the so-called forward flux sampling method for a more detailed description. This is done because forward flux sampling, despite its appealing simplicity and potential efficiency, has thus far received limited attention from practitioners. While path sampling methods have a widespread application to many types of rare transitional events, here only recent applications involving biomolecules are reviewed, including isomerization, protein folding, and enzyme catalysis.

Structuring and sampling complex conformation space: Weighted ensemble dynamics simulations

Based on multiple simulation trajectories, which started from dispersively selected initial conformations, the weighted ensemble dynamics method is designed to robustly and systematically explore the hierarchical structure of complex conformational space through the spectral analysis of the variance-covariance matrix of trajectory-mapped vectors. The nondegenerate ground state of the matrix directly predicts the ergodicity of simulation data. The ground state could be adopted as statistical weights of trajectories to correctly reconstruct the equilibrium properties, even though each trajectory only explores part of the conformational space. Otherwise, the degree of degeneracy simply gives the number of metastable states of the system under the time scale of individual trajectory. Manipulation on the eigenvectors leads to the classification of trajectories into nontransition ones within the states and transition ones between them. The transition states may also be predicted without a priori knowledge of the system. We demonstrate the application of the general method both to the system with a one-dimensional glassy potential and with the one of alanine dipeptide in explicit solvent.

A kinetic model of trp-cage folding from multiple biased molecular dynamics simulations

Trp-cage is a designed 20-residue polypeptide that, in spite of its size, shares several features with larger globular proteins.Although the system has been intensively investigated experimentally and theoretically, its folding mechanism is not yet fully understood. Indeed, some experiments suggest a two-state behavior, while others point to the presence of intermediates. In this work we show that the results of a bias-exchange metadynamics simulation can be used for constructing a detailed thermodynamic and kinetic model of the system. The model, although constructed from a biased simulation, has a quality similar to those extracted from the analysis of long unbiased molecular dynamics trajectories. This is demonstrated by a careful benchmark of the approach on a smaller system, the solvated Ace-Ala3-Nme peptide. For theTrp-cage folding, the model predicts that the relaxation time of 3100 ns observed experimentally is due to the presence of a compact molten globule-like conformation. This state has an occupancy of only 3% at 300 K, but acts as a kinetic trap.Instead, non-compact structures relax to the folded state on the sub-microsecond timescale. The model also predicts the presence of a state at Calpha-RMSD of 4.4 A from the NMR structure in which the Trp strongly interacts with Pro12. This state can explain the abnormal temperature dependence of the Pro12-delta3 and Gly11-alpha3 chemical shifts. The structures of the two most stable misfolded intermediates are in agreement with NMR experiments on the unfolded protein. Our work shows that, using biased molecular dynamics trajectories, it is possible to construct a model describing in detail the Trp-cage folding kinetics and thermodynamics in agreement with experimental data.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

The stochastic separatrix and the reaction coordinate for complex systems

We present a new approach to the identification of degrees of freedom which comprise a reaction coordinate in a complex system. The method begins with the generation of an ensemble of reactive trajectories. Each trajectory is analyzed for its equicommittor position or transition state; then the transition state ensemble is identified as the stochastic separatrix. Numerical analysis of the points along the separatrix for variability of coordinate location correctly identifies the components of the reaction coordinate in a test system of a double well coupled to a promoting vibration and a bath of linearly coupled oscillators.