Discrete path sampling

Decoding the energy landscape: extracting structure, dynamics and thermodynamics

Describing a potential energy surface in terms of local minima and the transition states that connect them provides a conceptual and computational framework for understanding and predicting observable properties. Visualizing the potential energy landscape using disconnectivity graphs supplies a graphical connection between different structure-seeking systems, which can relax efficiently to a particular morphology. Landscapes involving competing morphologies support multiple potential energy funnels, which may exhibit characteristic heat capacity features and relaxation time scales. These connections between the organization of the potential energy landscape and structure, dynamics and thermodynamics are common to all the examples presented, ranging from atomic and molecular clusters to biomolecules and soft and condensed matter. Further connections between motifs in the energy landscape and the interparticle forces can be developed using symmetry considerations and results from catastrophe theory.

Evolution of the potential energy landscape with static pulling force for two model proteins

The energy landscape is analyzed for off-lattice bead models of protein L and protein G as a function of a static pulling force. Two different pairs of attachment points (pulling directions) are compared in each case, namely, residues 1/56 and 10/32. For the terminal residue pulling direction 1/56, the distinct global minimum structures are all extended, aside from the compact geometry that correlates with zero force. The helical turns finally disappear at the highest pulling forces considered. For the 10/32 pulling direction, the changes are more complicated, with a variety of competing arrangements for beads outside the region where the force is directly applied. These alternatives produce frustrated energy landscapes, with low-lying minima separated by high barriers. The calculated folding pathways in the absence of force are in good agreement with previous work. The N-terminal hairpin folds first for protein L and the C-terminal hairpin for protein G, which exhibits an intermediate. However, for a relatively low static force, where the global minimum retains its structure, the folding mechanisms change, sometimes dramatically, depending on the protein and the attachment points. The scaling relations predicted by catastrophe theory are found to hold in the limit of short path lengths.

The Activation-Relaxation Technique: ART Nouveau and Kinetic ART

The evolution of many systems is dominated by rare activated events that occur on timescale ranging from nanoseconds to the hour or more. For such systems, simulations must leave aside the full thermal description to focus specifically on mechanisms that generate a configurational change. We present here the activation relaxation technique (ART), an open-ended saddle point search algorithm, and a series of recent improvements to ART nouveau and kinetic ART, an ART-based on-the-fly off-lattice self-learning kinetic Monte Carlo method.

Energy landscapes of colloidal clusters: thermodynamics and rearrangement mechanisms

New experiments involving direct observation of colloidal clusters by optical microscopy promise to deliver a wealth of new information about such systems. Calculations suggest that some of the observable properties may be predicted using a simple pairwise potential to represent the interparticle forces, but in a range of parameter space that is distinctly different from previous representations of atomic clusters. The present contribution provides some benchmark calculations and predictions of structure, thermodynamics and rearrangement mechanisms for colloidal clusters containing up to 80 particles. The results suggest that distinct features characteristic of short-ranged interactions should be observable in terms of the structure, thermodynamics and dynamical properties. Analysis of a kinetic transition network for the 19-particle cluster reveals super-Arrhenius behaviour in the dynamics, analogous to a 'fragile' glass-former.

Efficient Construction of Mesostate Networks from Molecular Dynamics Trajectories

The coarse-graining of data from molecular simulations yields conformational space networks that may be used for predicting the system's long time scale behavior, to discover structural pathways connecting free energy basins in the system, or simply to represent accessible phase space regions of interest and their connectivities in a two-dimensional plot. In this contribution, we present a tree-based algorithm to partition conformations of biomolecules into sets of similar microstates, i.e., to coarse-grain trajectory data into mesostates. On account of utilizing an architecture similar to that of established tree-based algorithms, the proposed scheme operates in near-linear time with data set size. We derive expressions needed for the fast evaluation of mesostate properties and distances when employing typical choices for measures of similarity between microstates. Using both a pedagogically useful and a real-word application, the algorithm is shown to be robust with respect to tree height, which in addition to mesostate threshold size is the main adjustable parameter. It is demonstrated that the derived mesostate networks can preserve information regarding the free energy basins and barriers by which the system is characterized.

A survey of the potential energy surface for the (benzene)13 cluster

We report on a survey of the potential energy surface for the 13-molecule benzene cluster, (C(6)H(6))(13), bound by an atom-atom intermolecular potential developed from first principles. The potential, which has an anisotropic repulsion term, is found to support distinct pairs of structures of C(3), C(i), and S(6) symmetry as low-lying minima, including a C(3) global minimum. The organisation of the low-lying region of the potential energy surface suggests that one of the S(6) structures is likely to act as a kinetic trap, hindering efficient relaxation to the global minimum, in agreement with the hypothesis that two isomers coexist in this system, suggested by spectroscopic experiments.

Quantitative analysis of ligand migration from transition networks

In this work we use transition network analysis for the first time to investigate ligand migration in truncated hemoglobin (trHbN) and obtain kinetic information about the docking-site dynamics in the protein. A comparison with explicit water molecular dynamics simulations (100 ns in total) shows that the rate constants derived from the network analysis are realistic. The transition network analysis provides 1) The time-resolved connectivity network in the protein; 2) The half-lives of the docking sites; 3) The transition timescales between two given docking sites; and 4) The extent of population transfer among different docking sites of the protein as a function of lag time. We investigate the role of the Tyr33 and Gln58 residues in ligand migration by studying ligand migration in four mutants of trHbN. The mutation study suggests that residues Tyr33 and Gln58 stabilize the NO ligand in the Xe2 docking site of trHbN, thus facilitating the efficiency of the NO detoxification reaction.

A Bayesian method for construction of Markov models to describe dynamics on various time-scales

The dynamics of many biological processes of interest, such as the folding of a protein, are slow and complicated enough that a single molecular dynamics simulation trajectory of the entire process is difficult to obtain in any reasonable amount of time. Moreover, one such simulation may not be sufficient to develop an understanding of the mechanism of the process, and multiple simulations may be necessary. One approach to circumvent this computational barrier is the use of Markov state models. These models are useful because they can be constructed using data from a large number of shorter simulations instead of a single long simulation. This paper presents a new Bayesian method for the construction of Markov models from simulation data. A Markov model is specified by (τ,P,T), where τ is the mesoscopic time step, P is a partition of configuration space into mesostates, and T is an N(P)×N(P) transition rate matrix for transitions between the mesostates in one mesoscopic time step, where N(P) is the number of mesostates in P. The method presented here is different from previous Bayesian methods in several ways. (1) The method uses Bayesian analysis to determine the partition as well as the transition probabilities. (2) The method allows the construction of a Markov model for any chosen mesoscopic time-scale τ. (3) It constructs Markov models for which the diagonal elements of T are all equal to or greater than 0.5. Such a model will be called a "consistent mesoscopic Markov model" (CMMM). Such models have important advantages for providing an understanding of the dynamics on a mesoscopic time-scale. The Bayesian method uses simulation data to find a posterior probability distribution for (P,T) for any chosen τ. This distribution can be regarded as the Bayesian probability that the kinetics observed in the atomistic simulation data on the mesoscopic time-scale τ was generated by the CMMM specified by (P,T). An optimization algorithm is used to find the most probable CMMM for the chosen mesoscopic time step. We applied this method of Markov model construction to several toy systems (random walks in one and two dimensions) as well as the dynamics of alanine dipeptide in water. The resulting Markov state models were indeed successful in capturing the dynamics of our test systems on a variety of mesoscopic time-scales.

Symmetrization of the AMBER and CHARMM force fields

The AMBER and CHARMM force fields are analyzed from the viewpoint of the permutational symmetry of the potential for feasible exchanges of identical atoms and chemical groups in amino and nucleic acids. In each case, we propose schemes for symmetrizing the potentials, which greatly facilitate the bookkeeping associated with constructing kinetic transition networks via geometry optimization.

Free-energy landscapes from adaptively biased methods: application to quantum systems

Several parallel adaptive biasing methods are applied to the calculation of free-energy pathways along reaction coordinates, choosing as a difficult example the double-funnel landscape of the 38-atom Lennard-Jones cluster. In the case of classical statistics, the Wang-Landau and adaptively biased molecular-dynamics (ABMD) methods are both found efficient if multiple walkers and replication and deletion schemes are used. An extension of the ABMD technique to quantum systems, implemented through the path-integral MD framework, is presented and tested on Ne38 against the quantum superposition method.

A self-learning algorithm for biased molecular dynamics

A new self-learning algorithm for accelerated dynamics, reconnaissance metadynamics, is proposed that is able to work with a very large number of collective coordinates. Acceleration of the dynamics is achieved by constructing a bias potential in terms of a patchwork of one-dimensional, locally valid collective coordinates. These collective coordinates are obtained from trajectory analyses so that they adapt to any new features encountered during the simulation. We show how this methodology can be used to enhance sampling in real chemical systems citing examples both from the physics of clusters and from the biological sciences.

Modelling proteins: conformational sampling and reconstruction of folding kinetics

In the last decades biomolecular simulation has made tremendous inroads to help elucidate biomolecular processes in-silico. Despite enormous advances in molecular dynamics techniques and the available computational power, many problems involve long time scales and large-scale molecular rearrangements that are still difficult to sample adequately. In this review we therefore summarise recent efforts to fundamentally improve this situation by decoupling the sampling of the energy landscape from the description of the kinetics of the process. Recent years have seen the emergence of many advanced sampling techniques, which permit efficient characterisation of the relevant family of molecular conformations by dispensing with the details of the short-term kinetics of the process. Because these methods generate thermodynamic information at best, they must be complemented by techniques to reconstruct the kinetics of the process using the ensemble of relevant conformations. Here we review recent advances for both types of methods and discuss their perspectives to permit efficient and accurate modelling of large-scale conformational changes in biomolecules. This article is part of a Special Issue entitled: Protein Dynamics: Experimental and Computational Approaches.

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.

The energy landscape, folding pathways and the kinetics of a knotted protein

The folding pathway and rate coefficients of the folding of a knotted protein are calculated for a potential energy function with minimal energetic frustration. A kinetic transition network is constructed using the discrete path sampling approach, and the resulting potential energy surface is visualized by constructing disconnectivity graphs. Owing to topological constraints, the low-lying portion of the landscape consists of three distinct regions, corresponding to the native knotted state and to configurations where either the N or C terminus is not yet folded into the knot. The fastest folding pathways from denatured states exhibit early formation of the N terminus portion of the knot and a rate-determining step where the C terminus is incorporated. The low-lying minima with the N terminus knotted and the C terminus free therefore constitute an off-pathway intermediate for this model. The insertion of both the N and C termini into the knot occurs late in the folding process, creating large energy barriers that are the rate limiting steps in the folding process. When compared to other protein folding proteins of a similar length, this system folds over six orders of magnitude more slowly.

The "weighted ensemble" path sampling method is statistically exact for a broad class of stochastic processes and binning procedures

The "weighted ensemble" method, introduced by Huber and Kim [Biophys. J. 70, 97 (1996)], is one of a handful of rigorous approaches to path sampling of rare events. Expanding earlier discussions, we show that the technique is statistically exact for a wide class of Markovian and non-Markovian dynamics. The derivation is based on standard path-integral (path probability) ideas, but recasts the weighted-ensemble approach as simple "resampling" in path space. Similar reasoning indicates that arbitrary nonstatic binning procedures, which merely guide the resampling process, are also valid. Numerical examples confirm the claims, including the use of bins which can adaptively find the target state in a simple model.

Emergent complexity from simple anisotropic building blocks: shells, tubes, and spirals

We describe a remarkably simple, generic, coarse-grained model involving anisotropic interactions, and characterize the global minima for clusters as a function of various parameters. Appropriate choices for the anisotropic interactions can reproduce a wide variety of complex morphologies as global minima, including spheroidal shells, tubular, helical and even head-tail morphologies, elucidating the physical principles that drive the assembly of these mesoscopic structures. Our model captures several experimental observations, such as the existence of competing morphologies, capsid polymorphism, and the effect of scaffolding proteins on capsid assembly.

Principal component analysis of potential energy surfaces of large clusters: allowing the practical calculation of the master equation

The number of variables in many-particle systems is typically unmanageably large; some way to reduce that number and still retain access to the important information about the system of interest is one of the great challenges in the broad topic of complexity. Principal components and principal coordinates provide a powerful means to extract--from unwieldy, large data sets--a reduced collection of variables that provide the information one needs, in a relatively efficient way and useful form. We investigate the application of principal components to the analysis of kinetics of the atomic motions in atomic clusters, particularly of clusters that are large enough so that a full description in terms of the entire high-dimensional potential surface is entirely impractical. A specific application is the use of principal components linking minima with their adjacent saddles, permitting the evaluation of rate coefficients (in the context of transition state theory) as ratios of partition functions of only one or two key variables.

Graph-based analysis of kinetics on multidimensional potential-energy surfaces

The aim of this paper is twofold: one is to give a detailed description of an alternative graph-based analysis method, which we call saddle connectivity graph, for analyzing the global topography and the dynamical properties of many-dimensional potential-energy landscapes and the other is to give examples of applications of this method in the analysis of the kinetics of realistic systems. A Dijkstra-type shortest path algorithm is proposed to extract dynamically dominant transition pathways by kinetically defining transition costs. The applicability of this approach is first confirmed by an illustrative example of a low-dimensional random potential. We then show that a coarse-graining procedure tailored for saddle connectivity graphs can be used to obtain the kinetic properties of 13- and 38-atom Lennard-Jones clusters. The coarse-graining method not only reduces the complexity of the graphs, but also, with iterative use, reveals a self-similar hierarchical structure in these clusters. We also propose that the self-similarity is common to many-atom Lennard-Jones clusters.

Recovering kinetics from a simplified protein folding model using replica exchange simulations: a kinetic network and effective stochastic dynamics

We present an approach to recover kinetics from a simplified protein folding model at different temperatures using the combined power of replica exchange (RE), a kinetic network, and effective stochastic dynamics. While RE simulations generate a large set of discrete states with the correct thermodynamics, kinetic information is lost due to the random exchange of temperatures. We show how we can recover the kinetics of a 2D continuous potential with an entropic barrier by using RE-generated discrete states as nodes of a kinetic network. By choosing the neighbors and the microscopic rates between the neighbors appropriately, the correct kinetics of the system can be recovered by running a kinetic simulation on the network. We fine-tune the parameters of the network by comparison with the effective drift velocities and diffusion coefficients of the system determined from short-time stochastic trajectories. One of the advantages of the kinetic network model is that the network can be built on a high-dimensional discretized state space, which can consist of multiple paths not consistent with a single reaction coordinate.

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.

Energy landscapes for shells assembled from pentagonal and hexagonal pyramids

We present new rigid body potentials that should favour efficient self-assembly of pentagonal and hexagonal pyramids into icosahedral shells over a wide range of temperature. By adding an extra repulsive site opposite the existing apex sites of the pyramids considered in a previously published model, frustrated energy landscapes are transformed into systems identified with self-assembling properties. The extra interaction may be considered analogous to a hydrophobic-hydrophilic repulsion, as in micelle formation.

Transition networks for modeling the kinetics of conformational change in macromolecules

The kinetics and thermodynamics of complex transitions in biomolecules can be modeled in terms of a network of transitions between the relevant conformational substates. Such a transition network, which overcomes the fundamental limitations of reaction-coordinate-based methods, can be constructed either based on the features of the energy landscape, or from molecular dynamics simulations. Energy-landscape-based networks are generated with the aid of automated path-optimization methods, and, using graph-theoretical adaptive methods, can now be constructed for large molecules such as proteins. Dynamics-based networks, also called Markov State Models, can be interpreted and adaptively improved using statistical concepts, such as the mean first passage time, reactive flux and sampling error analysis. This makes transition networks powerful tools for understanding large-scale conformational changes.

Folding pathways and rates for the three-stranded beta-sheet peptide Beta3s using discrete path sampling

The discrete path sampling method was used to investigate the folding of a three-stranded antiparallel beta-sheet peptide, Beta3s, described by an empirical potential and implicit solvent model. After application of a coarse-graining scheme that groups together sets of minima in local equilibrium, the calculated folding time was in reasonable agreement with other simulations and consistent with the experimental upper bound. The folding mechanism exhibited by the most significant discrete paths involves early formation of the C-terminal hairpin followed by docking of the N-terminal strand.

Pathways for conformational change in nitrogen regulatory protein C from discrete path sampling

Pathways corresponding to the conformational change in nitrogen regulatory protein C are calculated using the CHARMM19 force field with an implicit solvation model. Our analysis employs the discrete path sampling approach to grow a database of local minima and transition states from the potential energy surface that contains kinetically relevant pathways. The pathways with the largest contribution to the phenomenological two-state rate constants are found to exhibit a number of structural features that agree with experimental observations. Further details of the calculated pathways for conformational change may therefore provide useful predictions of how this large-scale motion is achieved.

Thermodynamics and kinetics of aggregation for the GNNQQNY peptide

The energy landscape of the monomer and dimer are explored for the amyloidogenic heptapeptide GNNQQNY from the N-terminal prion-determining domain of the yeast protein Sup35. The peptide is modeled by a united-atom potential and an implicit solvent representation. Replica exchange molecular dynamics is used to explore the conformational space, and discrete path sampling is employed to investigate the pathways that interconvert the most populated minima on the free energy surfaces. For the monomer, we find a rapid fluctuation between four different conformations, where a geometry intermediate between compact and extended structures is the most thermodynamically favorable. The GNNQQNY dimer forms three stable sheet structures, namely in-register parallel, off-register parallel, and antiparallel. The antiparallel dimer is stabilized by strong electrostatic interactions resulting from interpeptide hydrogen bonds, which restrict its conformational flexibility. The in-register parallel dimer, which is close to the amyloid beta-sheet structure, has fewer interpeptide hydrogen bonds, making hydrophobic interactions more important and increasing the conformational entropy compared to the antiparallel sheet. The estimated two-state rate constants indicate that the formation of dimers from monomers is fast and that the dimers are kinetically stable against dissociation at room temperature. Interconversions between the different dimers are feasible processes and are more likely than dissociation.

Path ensembles and path sampling in nonequilibrium stochastic systems

Markovian models based on the stochastic master equation are often encountered in single molecule dynamics, reaction networks, and nonequilibrium problems in chemistry, physics, and biology. An efficient and convenient method to simulate these systems is the kinetic Monte Carlo algorithm which generates continuous-time stochastic trajectories. We discuss an alternative simulation method based on sampling of stochastic paths. Utilizing known probabilities of stochastic paths, it is possible to apply Metropolis Monte Carlo in path space to generate a desired ensemble of stochastic paths. The method is a generalization of the path sampling idea to stochastic dynamics, and is especially suited for the analysis of rare paths which are not often produced in the standard kinetic Monte Carlo procedure. Two generic examples are presented to illustrate the methodology.

Connectivity graph: multiple connectivity on potential energy surface does affect the dynamics

A mapping called a connectivity graph is proposed to visualize the complicated networks between local minima via saddles on many-dimensional potential energy surfaces. With this mapping, the effect of the topography of a complex potential energy surface on the dynamics is studied in model funnel potentials and a Lennard-Jones cluster. We present strong evidence that energetically expensive but dynamically relevant saddles are indispensable for kinetic dynamics. The results show that the connectivity graphs provide a sound basis for understanding nonequilibrium dynamics.

Potential energy and free energy landscapes

Familiar concepts for small molecules may require reinterpretation for larger systems. For example, rearrangements between geometrical isomers are usually considered in terms of transitions between the corresponding local minima on the underlying potential energy surface, V. However, transitions between bulk phases such as solid and liquid, or between the denatured and native states of a protein, are normally addressed in terms of free energy minima. To reestablish a connection with the potential energy surface we must think in terms of representative samples of local minima of V, from which a free energy surface is projected by averaging over most of the coordinates. The present contribution outlines how this connection can be developed into a tool for quantitative calculations. In particular, stepping between the local minima of V provides powerful methods for locating the global potential energy minimum, and for calculating global thermodynamic properties. When the transition states that link local minima are also sampled we can exploit statistical rate theory to obtain insight into global dynamics and rare events. Visualizing the potential energy landscape helps to explain how the network of local minima and transition states determines properties such as heat capacity features, which signify transitions between free energy minima. The organization of the landscape also reveals how certain systems can reliably locate particular structures on the experimental time scale from among an exponentially large number of local minima. Such directed searches not only enable proteins to overcome Levinthal's paradox but may also underlie the formation of "magic numbers" in molecular beams, the self-assembly of macromolecular structures, and crystallization.