Automatic discovery of metastable states for the construction of Markov models of macromolecular conformational dynamics

To meet the challenge of modeling the conformational dynamics of biological macromolecules over long time scales, much recent effort has been devoted to constructing stochastic kinetic models, often in the form of discrete-state Markov models, from short molecular dynamics simulations. To construct useful models that faithfully represent dynamics at the time scales of interest, it is necessary to decompose configuration space into a set of kinetically metastable states. Previous attempts to define these states have relied upon either prior knowledge of the slow degrees of freedom or on the application of conformational clustering techniques which assume that conformationally distinct clusters are also kinetically distinct. Here, we present a first version of an automatic algorithm for the discovery of kinetically metastable states that is generally applicable to solvated macromolecules. Given molecular dynamics trajectories initiated from a well-defined starting distribution, the algorithm discovers long lived, kinetically metastable states through successive iterations of partitioning and aggregating conformation space into kinetically related regions. The authors apply this method to three peptides in explicit solvent-terminally blocked alanine, the 21-residue helical F(s) peptide, and the engineered 12-residue beta-hairpin trpzip2-to assess its ability to generate physically meaningful states and faithful kinetic models.     

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.     

Algorithms for clustering molecular dynamics configurations

Two traditional clustering algorithms are applied to configurations from a long molecular dynamics trajectory and compared using two sets of test data. First, a subset of atoms was chosen to present conformations which naturally fall into a number of clusters. Second, a subset of atoms was selected to span a relatively continuous region of conformational space rather than form discrete conformational classes. Of the two algorithms used, the single linkage method is inappropriate for this kind of data. The divisive hierarchical method, based on minimizing the difference between cluster centroids and extrema, is successful but also prone to imposing clustering hierarchy where none can be justified. © 1994 by John Wiley & Sons, Inc.     

Efficient evaluation of sampling quality of molecular dynamics simulations by clustering of dihedral torsion angles and Sammon mapping

A direct conformational clustering and mapping approach for peptide conformations based on backbone dihedral angles has been developed and applied to compare conformational sampling of Met-enkephalin using two molecular dynamics (MD) methods. Efficient clustering in dihedrals has been achieved by evaluating all combinations resulting from independent clustering of each dihedral angle distribution, thus resolving all conformational substates. In contrast, Cartesian clustering was unable to accurately distinguish between all substates. Projection of clusters on dihedral principal component (PCA) subspaces did not result in efficient separation of highly populated clusters. However, representation in a nonlinear metric by Sammon mapping was able to separate well the 48 highest populated clusters in just two dimensions. In addition, this approach also allowed us to visualize the transition frequencies between clusters efficiently. Significantly, higher transition frequencies between more distinct conformational substates were found for a recently developed biasing-potential replica exchange MD simulation method allowing faster sampling of possible substates compared to conventional MD simulations. Although the number of theoretically possible clusters grows exponentially with peptide length, in practice, the number of clusters is only limited by the sampling size (typically much smaller), and therefore the method is well suited also for large systems. The approach could be useful to rapidly and accurately evaluate conformational sampling during MD simulations, to compare different sampling strategies and eventually to detect kinetic bottlenecks in folding pathways.     

Hierarchical analysis of conformational dynamics in biomolecules: transition networks of metastable states

Molecular dynamics simulation generates large quantities of data that must be interpreted using physically meaningful analysis. A common approach is to describe the system dynamics in terms of transitions between coarse partitions of conformational space. In contrast to previous work that partitions the space according to geometric proximity, the authors examine here clustering based on kinetics, merging configurational microstates together so as to identify long-lived, i.e., dynamically metastable, states. As test systems microsecond molecular dynamics simulations of the polyalanines Ala(8) and Ala(12) are analyzed. Both systems clearly exhibit metastability, with some kinetically distinct metastable states being geometrically very similar. Using the backbone torsion rotamer pattern to define the microstates, a definition is obtained of metastable states whose lifetimes considerably exceed the memory associated with interstate dynamics, thus allowing the kinetics to be described by a Markov model. This model is shown to be valid by comparison of its predictions with the kinetics obtained directly from the molecular dynamics simulations. In contrast, clustering based on the hydrogen-bonding pattern fails to identify long-lived metastable states or a reliable Markov model. Finally, an approach is proposed to generate a hierarchical model of networks, each having a different number of metastable states. The model hierarchy yields a qualitative understanding of the multiple time and length scales in the dynamics of biomolecules.     

Semiclassical nonadiabatic dynamics based on quantum trajectories for the O(3P,1D) + H2 system

The O(3P,1D) + H2 --> OH + H reaction is studied using trajectory dynamics within the approximate quantum potential approach. Calculations of the wave-packet reaction probabilities are performed for four coupled electronic states for total angular momentum J = 0 using a mixed coordinate/polar representation of the wave function. Semiclassical dynamics is based on a single set of trajectories evolving on an effective potential-energy surface and in the presence of the approximate quantum potential. Population functions associated with each trajectory are computed for each electronic state. The effective surface is a linear combination of the electronic states with the contributions of individual components defined by their time-dependent average populations. The wave-packet reaction probabilities are in good agreement with the quantum-mechanical results. Intersystem crossing is found to have negligible effect on reaction probabilities summed over final electronic states.     

Validation of Markov state models using Shannon's entropy

Markov state models are kinetic models built from the dynamics of molecular simulation trajectories by grouping similar configurations into states and examining the transition probabilities between states. Here we present a procedure for validating the underlying Markov assumption in Markov state models based on information theory using Shannon's entropy. This entropy method is applied to a simple system and is compared with the previous eigenvalue method. The entropy method also provides a way to identify states that are least Markovian, which can then be divided into finer states to improve the model.     

Identification of metastable states in peptide's dynamics

A recently developed spectral method for identifying metastable states in Markov chains is used to analyze the conformational dynamics of a four-residue peptide valine-proline-alanine-leucine. We compare our results to empirically defined conformational states and show that the found metastable states correctly reproduce the conformational dynamics of the system.     

Milestoning with transition memory

Milestoning is a method used to calculate the kinetics and thermodynamics of molecular processes occurring on time scales that are not accessible to brute force molecular dynamics (MD). In milestoning, the conformation space of the system is sectioned by hypersurfaces (milestones), an ensemble of trajectories is initialized on each milestone, and MD simulations are performed to calculate transitions between milestones. The transition probabilities and transition time distributions are then used to model the dynamics of the system with a Markov renewal process, wherein a long trajectory of the system is approximated as a succession of independent transitions between milestones. This approximation is justified if the transition probabilities and transition times are statistically independent. In practice, this amounts to a requirement that milestones are spaced such that trajectories lose position and velocity memory between subsequent transitions. Unfortunately, limiting the number of milestones limits both the resolution at which a system's properties can be analyzed, and the computational speedup achieved by the method. We propose a generalized milestoning procedure, milestoning with transition memory (MTM), which accounts for memory of previous transitions made by the system. When a reaction coordinate is used to define the milestones, the MTM procedure can be carried out at no significant additional expense as compared to conventional milestoning. To test MTM, we have applied its version that allows for the memory of the previous step to the toy model of a polymer chain undergoing Langevin dynamics in solution. We have computed the mean first passage time for the chain to attain a cyclic conformation and found that the number of milestones that can be used, without incurring significant errors in the first passage time is at least 8 times that permitted by conventional milestoning. We further demonstrate that, unlike conventional milestoning, MTM permits milestones to be spaced such that trajectories do not have enough time to lose their velocity memory between successively crossed milestones.     

Effects of solute-solvent proton exchange on polypeptide chain dynamics: a constant-pH molecular dynamics study

A method for performing implicit-solvent molecular dynamics simulations at constant pH was applied to a pentapeptide acetyl-Ala-Asp-Ala-Lys-Ala-amide at pH 4. As a reference, molecular dynamics simulations were done for the same peptide with two variants of its fixed protonation patterns expected to dominate at pH 4, i.e., with a protonated and a deprotonated side chain of the Asp residue and the protonated Lys residue in both cases. The dynamic trajectories of the peptide were used to discuss the problem of the significance of the solute-solvent proton exchange phenomena for the dynamics and structural distributions of the polypeptide chain. The Asp-Lys distance was used as a probe of the overall molecular structure of the investigated pentapeptide. To characterize the dynamics, distributions of the "waiting" times for a transition from a "short" distance conformation to a "long" distance conformation were constructed, based on the generated molecular dynamics trajectories. We show that the relaxation time for the transitions, derived from the constant-pH simulations, is very close to the relaxation time characterizing a permanently protonated molecule, although the average protonation probability of the short-distance conformation is close to zero. However, the distribution of the Asp-Lys distances obtained from constant-pH simulations cannot be reproduced as a linear combination of the distributions resulting from the simulations with fixed protonation states.     

On the calculation of velocity-dependent properties in molecular dynamics simulations using the leapfrog integration algorithm

Widely used programs for molecular dynamics simulation of (bio)molecular systems are the Verlet and leapfrog algorithms. In these algorithms, the particle velocities are less accurately propagated than the positions. Important quantities for the simulation such as the temperature and the pressure involve the squared velocities at full time steps. Here, we derive an expression for the squared particle velocity at full time step in the leapfrog scheme, which is more accurate than the standardly used one. In particular, this allows us to show that the full time step kinetic energy of a particle is more accurately computed as the average of the kinetic energies at previous and following half steps than as the square of the average velocity as implemented in various molecular dynamics codes. Use of the square of the average velocity introduces a systematic bias in the calculation of the instantaneous temperature and pressure of a molecular dynamics system. We show the consequences when the system is coupled to a thermostat and a barostat.     

Calculation of the distribution of eigenvalues and eigenvectors in Markovian state models for molecular dynamics

Markovian state models (MSMs) are a convenient and efficient means to compactly describe the kinetics of a molecular system as well as a formalism for using many short simulations to predict long time scale behavior. Building a MSM consists of grouping the conformations into states and estimating the transition probabilities between these states. In a previous paper, we described an efficient method for calculating the uncertainty due to finite sampling in the mean first passage time between two states. In this paper, we extend the uncertainty analysis to derive similar closed-form solutions for the distributions of the eigenvalues and eigenvectors of the transition matrix, quantities that have numerous applications when using the model. We demonstrate the accuracy of the distributions on a six-state model of the terminally blocked alanine peptide. We also show how to significantly reduce the total number of simulations necessary to build a model with a given precision using these uncertainty estimates for the blocked alanine system and for a 2454-state MSM for the dynamics of the villin headpiece.     

Optimal use of data in parallel tempering simulations for the construction of discrete-state Markov models of biomolecular dynamics

Parallel tempering (PT) molecular dynamics simulations have been extensively investigated as a means of efficient sampling of the configurations of biomolecular systems. Recent work has demonstrated how the short physical trajectories generated in PT simulations of biomolecules can be used to construct the Markov models describing biomolecular dynamics at each simulated temperature. While this approach describes the temperature-dependent kinetics, it does not make optimal use of all available PT data, instead estimating the rates at a given temperature using only data from that temperature. This can be problematic, as some relevant transitions or states may not be sufficiently sampled at the temperature of interest, but might be readily sampled at nearby temperatures. Further, the comparison of temperature-dependent properties can suffer from the false assumption that data collected from different temperatures are uncorrelated. We propose here a strategy in which, by a simple modification of the PT protocol, the harvested trajectories can be reweighted, permitting data from all temperatures to contribute to the estimated kinetic model. The method reduces the statistical uncertainty in the kinetic model relative to the single temperature approach and provides estimates of transition probabilities even for transitions not observed at the temperature of interest. Further, the method allows the kinetics to be estimated at temperatures other than those at which simulations were run. We illustrate this method by applying it to the generation of a Markov model of the conformational dynamics of the solvated terminally blocked alanine peptide.     

Coherent dynamics of meta-toluidine investigated by quasielastic neutron scattering

The coherent dynamics of a typical fragile glass former, meta-toluidine, was investigated at the molecular level using quasielastic neutron scattering, with time-of-flight and neutron spin echo spectrometers. It is well known that the static structure factor of meta-toluidine shows a prepeak originating from clustering of the molecules through hydrogen bonding between the amine groups. The dynamics of meta-toluidine was measured for several values of the wavevector transfer Q, which is equivalent to an inverse length scale, in a range encompassing the prepeak and the structure factor peak. Data were collected in the temperature range corresponding to the liquid and supercooled states, down to the glass transition. At least two dynamical processes were identified. This paper focuses on the slowest relaxation process in the system, the α-relaxation, which was found to scale with the macroscopic shear viscosity at all the investigated Q values. No evidence of "de Gennes" narrowing associated with the prepeak was observed, in contrast with what happens at the Q value corresponding to the interparticle distance. Moreover, using partially deuterated samples, the dynamics of the clusters was found to be correlated to the single-particle dynamics of the meta-toluidine molecules.     

Using molecular dynamics and quantum mechanics calculations to model fluorescence observables

We provide a critical examination of two different methods for generating a donor-acceptor electronic coupling trajectory from a molecular dynamics (MD) trajectory and three methods for sampling that coupling trajectory, allowing the modeling of experimental observables directly from the MD simulation. In the first coupling method we perform a single quantum-mechanical (QM) calculation to characterize the excited state behavior, specifically the transition dipole moment, of the fluorescent probe, which is then mapped onto the configuration space sampled by MD. We then utilize these transition dipoles within the ideal dipole approximation (IDA) to determine the electronic coupling between the probes that mediates the transfer of energy. In the second method we perform a QM calculation on each snapshot and use the complete transition densities to calculate the electronic coupling without need for the IDA. The resulting coupling trajectories are then sampled using three methods ranging from an independent sampling of each trajectory point (the independent snapshot method) to a Markov chain treatment that accounts for the dynamics of the coupling in determining effective rates. The results show that the IDA significantly overestimates the energy transfer rate (by a factor of 2.6) during the portions of the trajectory in which the probes are close to each other. Comparison of the sampling methods shows that the Markov chain approach yields more realistic observables at both high and low FRET efficiencies. Differences between the three sampling methods are discussed in terms of the different mechanisms for averaging over structural dynamics in the system. Convergence of the Markov chain method is carefully examined. Together, the methods for estimating coupling and for sampling the coupling provide a mechanism for directly connecting the structural dynamics modeled by MD with fluorescence observables determined through FRET experiments.     

Adaptive partitioning by local density‐peaks: An efficient density‐based clustering algorithm for analyzing molecular dynamics trajectories

We present an efficient density‐based adaptive‐resolution clustering method APLoD for analyzing large‐scale molecular dynamics (MD) trajectories. APLoD performs the k‐nearest‐neighbors search to estimate the density of MD conformations in a local fashion, which can group MD conformations in the same high‐density region into a cluster. APLoD greatly improves the popular density peaks algorithm by reducing the running time and the memory usage by 2–3 orders of magnitude for systems ranging from alanine dipeptide to a 370‐residue Maltose‐binding protein. In addition, we demonstrate that APLoD can produce clusters with various sizes that are adaptive to the underlying density (i.e., larger clusters at low‐density regions, while smaller clusters at high‐density regions), which is a clear advantage over other popular clustering algorithms including k‐centers and k‐medoids. We anticipate that APLoD can be widely applied to split ultra‐large MD datasets containing millions of conformations for subsequent construction of Markov State Models. © 2016 Wiley Periodicals, Inc.     

Ligand-to-metal charge-transfer dynamics in a blue copper protein plastocyanin: a molecular dynamics study

Equilibrium and nonequilibrium dynamics of a blue copper protein plastocyanin in an oxidized state are studied by molecular dynamics (MD) simulation. Potential energy functions of the lowest seven electronic states, including ligand-to-metal charge-transfer (LMCT) and copper d --> d excited states, were taken from our previous work (Ando, K. J. Phys. Chem. B 2004, 108, 3940), which employed ab initio molecular orbital and density functional calculations on the active-site model. The equilibrium MD simulations in the ground state indicate that ligand motions coupled to transition from the ground state to the LMCT state are mostly represented by stretching and bending vibrations of the Cu-S(Cys) distance, Ndelta(His)-Cu-Ndelta(His) angle, and S(Cys)-Cu-[Ndelta(His)]2 trigonal pyramid structure. The nonequilibrium dynamics on the LMCT potential exhibit rapid decays in which surface crossings to the d --> d and the first excited states occur in 70-80 fs. The crossing dynamics mostly correlate with cleavage of the Cu-S(Cys) bond and the associated response in the Ndelta(His)-Cu-Ndelta(His) moiety. The average dynamics of the vertical energy gap coordinates exhibit an overdamped decay with a recurrence oscillation in 500 fs, which shows clear coherence surviving after the ensemble averaging. This oscillation stems mostly from the recoiling motion of the Ndelta(His)-Cu-Ndelta(His) part. The dynamics of the energy gaps after this coherent oscillation are randomized such that the ensemble average yields flat profiles along time, although each single trajectory exhibits fluctuations with amplitudes large enough to reach surface crossings. These indicate that the relaxation from the LMCT state first occurs via ballistic and coherent potential crossings in 70-80 and 500 fs, followed by thermally activated random transitions.     

Communication: κ-dynamics--an exact method for accelerating rare event classical molecular dynamics

κ-dynamics is an accelerated molecular dynamics method for systems with slow transitions between stable states. Short trajectories are integrated from a transition state separating a reactant state from products. The first trajectory found that leads directly to a product without recrossing the transition state and starts in the reactant state is followed. The transition time is drawn from a distribution given by the transition state theory rate and the number of attempted trajectories. Repeating this procedure from each state visited gives a statistically exact state-to-state trajectory.     

Generation of initial trajectories for transition path sampling of chemical reactions with ab initio molecular dynamics

Transition path sampling is an innovative method for focusing a molecular dynamics simulation on a reactive event. Although transition path sampling methods can generate an ensemble of reactive trajectories, an initial reactive trajectory must be generated by some other means. In this paper, the authors have evaluated three methods for generating initial reactive trajectories for transition path sampling with ab initio molecular dynamics. The authors have tested each of these methods on a set of chemical reactions involving the breaking and making of covalent bonds: the 1,2-hydrogen elimination in the borane-ammonia adduct, a tautomerization, and the Claisen rearrangement. The first method is to initiate trajectories from the potential energy transition state, which was effective for all reactions in the test set. Assigning atomic velocities found using normal mode analysis greatly improved the success of this method. The second method uses a high temperature molecular dynamics simulation and then iteratively reduces the total energy of the simulation until a low temperature reactive trajectory is found. This was effective in generating a low temperature trajectory from an initial trajectory run at 3000 K of the tautomerization reaction, although it failed for the other two. The third uses an orbital based bias potential to find a reactive trajectory and uses this trajectory to initiate an unbiased trajectory. The authors found that a highest occupied molecular orbital-lowest unoccupied molecular orbital bias could be used to find a reactive trajectory for the Claisen rearrangement, although it failed for the other two reactions. These techniques will help make it practical to use transition path sampling to study chemical reaction mechanisms that involve bond breaking and forming.