Transition path sampling for discrete master equations with absorbing states

Transition path sampling (TPS) algorithms have been implemented with deterministic dynamics, with thermostatted dynamics, with Brownian dynamics, and with simple spin flip dynamics. Missing from the TPS repertoire is an implementation with kinetic Monte Carlo (kMC), i.e., with the underlying dynamics coming from a discrete master equation. We present a new hybrid kMC-TPS algorithm and prove that it satisfies detailed balance in the transition path ensemble. The new algorithm is illustrated for a simplified Markov State Model of trp-cage folding. The transition path ensemble from kMC-TPS is consistent with that obtained from brute force kMC simulations. The committor probabilities and local fluxes for the simple model are consistent with those obtained from exact methods for simple master equations. The new kMC-TPS method should be useful for analysis of rare transitions in complex master equations where the individual states cannot be enumerated and therefore where exact solutions cannot be obtained.     

Isotropic-nematic phase transition in the Lebwohl-Lasher model from density of states simulations

Density of states Monte Carlo simulations have been performed to study the isotropic-nematic (IN) transition of the Lebwohl-Lasher model for liquid crystals. The IN transition temperature was calculated as a function of system size using expanded ensemble density of states simulations with histogram reweighting. The IN temperature for infinite system size was obtained by extrapolation of three independent measures. A subsequent analysis of the kinetics in the model showed that the transition occurs via spinodal decomposition through aggregation of clusters of liquid crystal molecules.     

Improved transition path sampling methods for simulation of rare events

The free energy surfaces of a wide variety of systems encountered in physics, chemistry, and biology are characterized by the existence of deep minima separated by numerous barriers. One of the central aims of recent research in computational chemistry and physics has been to determine how transitions occur between deep local minima on rugged free energy landscapes, and transition path sampling (TPS) Monte-Carlo methods have emerged as an effective means for numerical investigation of such transitions. Many of the shortcomings of TPS-like approaches generally stem from their high computational demands. Two new algorithms are presented in this work that improve the efficiency of TPS simulations. The first algorithm uses biased shooting moves to render the sampling of reactive trajectories more efficient. The second algorithm is shown to substantially improve the accuracy of the transition state ensemble by introducing a subset of local transition path simulations in the transition state. The system considered in this work consists of a two-dimensional rough energy surface that is representative of numerous systems encountered in applications. When taken together, these algorithms provide gains in efficiency of over two orders of magnitude when compared to traditional TPS simulations.     

Bottom-up view of water network-mediated CO2 reduction using cryogenic cluster ion spectroscopy and direct dynamics simulations

The transition states of a chemical reaction in solution are generally accessed through exchange of thermal energy between the solvent and the reactants. As such, an ensemble of reacting systems approaches the transition state configuration of reactant and surrounding solvent in an incoherent manner that does not lend itself to direct experimental observation. Here we describe how gas-phase cluster chemistry can provide a detailed picture of the microscopic mechanics at play when a network of six water molecules mediates the trapping of a highly reactive "hydrated electron" onto a neutral CO(2) molecule to form a radical anion. The exothermic reaction is triggered from a metastable intermediate by selective excitation of either the reactant CO(2) or the water network, which is evidenced by the evaporative decomposition of the product cluster. Ab initio molecular dynamics simulations of energized CO(2)·(H(2)O)(6)(-) clusters are used to elucidate the nature of the network deformations that mediate intracluster electron capture, thus revealing the detailed solvent fluctuations implicit in the Marcus theory for electron-transfer kinetics in solution.     

Atomistic simulations of the solid-liquid transition of 1-ethyl-3-methyl imidazolium bromide ionic liquid

Achieving melting point around room temperature is important for applications of ionic liquids. In this work, molecular dynamics simulations are carried out to investigate the solid-liquid transition of ionic liquid 1-ethyl-3-methyl imidazolium bromide ([emim]Br) by direct heating, hysteresis, void-nucleation, sandwich, and microcanonical ensemble approaches. Variations of the non-bonded energy, density, diffusion coefficient, and translational order parameter of [emim]Br are analyzed as a function of temperature, and a coexisting solid-liquid system is achieved in the microcanonical ensemble method. The melting points obtained from the first three methods are 547 ± 8 K, 429 ± 8 K, and 370 ± 6 K; while for the sandwich method, the melting points are 403 ± 4 K when merging along the x-axis by anisotropic isothermal-isobaric (NPT) ensemble, 393 ± 4 K when along the y-axis by anisotropic NPT ensemble, and 375 ± 4 K when along the y-axis by isotropic NPT ensemble. For microcanonical ensemble method, when the slabs are merging along different directions (x-axis, y-axis, and z-axis), the melting points are 364 ± 3 K, 365 ± 3 K, and 367 ± 3 K, respectively, the melting points we get by different methods are approximately 55.4%, 21.9%, 5.1%, 14.5%, 11.6%, 6.5%, 3.4%, 3.7%, and 4.3% higher than the experimental value of 352 K. The advantages and disadvantages of each method are discussed. The void-nucleation and microcanonical ensemble methods are most favorable for predicting the solid-liquid transition.     

Simulating rare events in equilibrium or nonequilibrium stochastic systems

We present three algorithms for calculating rate constants and sampling transition paths for rare events in simulations with stochastic dynamics. The methods do not require a priori knowledge of the phase-space density and are suitable for equilibrium or nonequilibrium systems in stationary state. All the methods use a series of interfaces in phase space, between the initial and final states, to generate transition paths as chains of connected partial paths, in a ratchetlike manner. No assumptions are made about the distribution of paths at the interfaces. The three methods differ in the way that the transition path ensemble is generated. We apply the algorithms to kinetic Monte Carlo simulations of a genetic switch and to Langevin dynamics simulations of intermittently driven polymer translocation through a pore. We find that the three methods are all of comparable efficiency, and that all the methods are much more efficient than brute-force simulation.     

Determination of the melting point of hard spheres from direct coexistence simulation methods

We consider the computation of the coexistence pressure of the liquid-solid transition of a system of hard spheres from direct simulation of the inhomogeneous system formed from liquid and solid phases separated by an interface. Monte Carlo simulations of the interfacial system are performed in three different ensembles. In a first approach, a series of simulations is carried out in the isothermal-isobaric ensemble, where the solid is allowed to relax to its equilibrium crystalline structure, thus avoiding the appearance of artificial stress in the system. Here, the total volume of the system fluctuates due to changes in the three dimensions of the simulation box. In a second approach, we consider simulations of the inhomogeneous system in an isothermal-isobaric ensemble where the normal pressure, as well as the area of the (planar) fluid-solid interface, are kept constant. Now, the total volume of the system fluctuates due to changes in the longitudinal dimension of the simulation box. In both approaches, the coexistence pressure is estimated by monitoring the evolution of the density along several simulations carried out at different pressures. Both routes are seen to provide consistent values of the fluid-solid coexistence pressure, p=11.54(4)k(B)T/sigma(3), which indicates that the error introduced by the use of the standard constant-pressure ensemble for this particular problem is small, provided the systems are sufficiently large. An additional simulation of the interfacial system is conducted in a canonical ensemble where the dimensions of the simulation box are allowed to change subject to the constraint that the total volume is kept fixed. In this approach, the coexistence pressure corresponds to the normal component of the pressure tensor, which can be computed as an appropriate ensemble average in a single simulation. This route yields a value of p=11.54(4)k(B)T/sigma(3). We conclude that the results obtained for the coexistence pressure from direct simulations of the liquid and solid phases in coexistence using different ensembles are mutually consistent and are in excellent agreement with the values obtained from free energy calculations.     

Folding transition-state and denatured-state ensembles of FSD-1 from folding and unfolding simulations

Characterization of the folding transition-state ensemble and the denatured-state ensemble is an important step toward a full elucidation of protein folding mechanisms. We report herein an investigation of the free-energy landscape of FSD-1 protein by a total of four sets of folding and unfolding molecular dynamics simulations with explicit solvent. The transition-state ensemble was initially identified from unfolding simulations at 500 K and was verified by simulations at 300 K starting from the ensemble structures. The denatured-state ensemble and the early-stage folding were studied by a combination of unfolding simulations at 500 K and folding simulations at 300 K starting from the extended conformation. A common feature of the transition-state ensemble was the substantial formation of the native secondary structures, including both the alpha-helix and beta-sheet, with partial exposure of the hydrophobic core in the solvent. Both the native and non-native secondary structures were observed in the denatured-state ensemble and early-stage folding, consistent with the smooth experimental melting curve. Interestingly, the contact orders of the transition-state ensemble structures were similar to that of the native structure and were notably lower than those of the compact structures found in early-stage folding, implying that chain and topological entropy might play significant roles in protein folding. Implications for FSD-1 folding mechanisms and the rate-limiting step are discussed. Analyses further revealed interesting non-native interactions in the denatured-state ensemble and early-stage folding and the possibility that destabilization of these interactions could help to enhance the stability and folding rate of the protein.     

Selective sampling of transition paths

In this short paper, we introduce an approximate method for the quick estimate of rate constants based on a simple sampling method of reactive transition paths over high energy barriers. It makes use of the previously introduced accelerated molecular dynamics (MD) simulation method to generate initial points for trajectory shooting. The accelerated MD simulations, although with the loss of real dynamics, lead to a quick calculation of thermodynamic properties and at the same time produce an ensemble of configurations with an enhanced sampling over the phase space that is more "reactive." The forward/backward trajectory shooting as that used in the transition path sampling method is then initiated from the configurations obtained from accelerated MD simulations to generate transition paths on the original unbiased potential. This method selectively enhances sampling of successful trajectories and at the same time accelerates significantly the calculation of rate constants.     

Determination of fluid--solid transitions in model protein solutions using the histogram reweighting method and expanded ensemble simulations

Protein crystallization conditions are usually identified by empirical screening methods because of the complexity of the process, such as the existence of nonequilibrium phases and the different crystal forms that may result from changes in solution conditions. Here the crystallization of a model protein is studied using computer simulation. The model consists of spheres that have both an isotropic interaction of short range and anisotropic interactions between patch-antipatch pairs. The free energy of a protein crystal is calculated using expanded ensemble simulations of the Einstein crystal, and NpT-Monte Carlo simulations with histogram reweighting are used to determine the fluid-solid coexistence. The histogram reweighting method is also used to trace out the complete coexistence curve, including multiple crystal phases, with varying reduced temperature, which corresponds to changing solution conditions. At a patch-antipatch interaction strength five times that of the isotropic interaction, the protein molecules form a stable simple cubic structure near room temperature, whereas an orientationally disordered face-centered-cubic structure is favored at higher temperatures. The anisotropic attractions also lead to a weak first-order transition between orientationally disordered and ordered face-centered-cubic structures at low temperature, although this transition is metastable. A complete phase diagram, including a fluid phase, three solid phases, and two triple points, is found for the six-patch protein model. A 12-patch protein model, consistent with the face-centered-cubic structure, leads to greater thermodynamic stability of the ordered phase. Metastable liquid-liquid phase equilibria for isotropic models with varying attraction tails are also predicted from Gibbs ensemble simulations.     

RedMD—Reduced molecular dynamics package

We developed a software package (RedMD) to perform molecular dynamics simulations and normal mode analysis of reduced models of proteins, nucleic acids, and their complexes. With RedMD one can perform molecular dynamics simulations in a microcanonical ensemble, with Berendsen and Langevin thermostats, and with Brownian dynamics. We provide force field and topology generators which are based on the one‐bead per residue/nucleotide elastic network model and its extensions. The user can change the force field parameters with the command line options that are passed to generators. Also, the generators can be modified, for example, to add new potential energy functions. Normal mode analysis tool is available for elastic or anisotropic network models. The program is written in C and C++ languages and the structure/topology of a molecule is based on an XML format. OpenMP technology for shared‐memory architectures was used for code parallelization. The code is distributed under GNU public licence and available at http://bionano.icm.edu.pl/software/ . © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Neural network studies. 4. Introduction to associative neural networks

Associative neural network (ASNN) represents a combination of an ensemble of feed-forward neural networks and the k-nearest neighbor technique. This method uses the correlation between ensemble responses as a measure of distance amid the analyzed cases for the nearest neighbor technique. This provides an improved prediction by the bias correction of the neural network ensemble. An associative neural network has a memory that can coincide with the training set. If new data becomes available, the network further improves its predictive ability and provides a reasonable approximation of the unknown function without a need to retrain the neural network ensemble. This feature of the method dramatically improves its predictive ability over traditional neural networks and k-nearest neighbor techniques, as demonstrated using several artificial data sets and a program to predict lipophilicity of chemical compounds. Another important feature of ASNN is the possibility to interpret neural network results by analysis of correlations between data cases in the space of models. It is shown that analysis of such correlations makes it possible to provide "property-targeted" clustering of data. The possible applications and importance of ASNN in drug design and medicinal and combinatorial chemistry are discussed. The method is available on-line at http://www.vcclab.org/lab/asnn.     

Probing molecular kinetics with Markov models: metastable states, transition pathways and spectroscopic observables

Markov (state) models (MSMs) have attracted a lot of interest recently as they (1) can probe long-term molecular kinetics based on short-time simulations, (2) offer a way to analyze great amounts of simulation data with relatively little subjectivity of the analyst, (3) provide insight into microscopic quantities such as the ensemble of transition pathways, and (4) allow simulation data to be reconciled with measurement data in a rigorous and explicit way. Here we sketch our current perspective of Markov models and explain in short their theoretical basis and assumptions. We describe transition path theory which allows the entire ensemble of protein folding pathways to be investigated and that combines naturally with Markov models. Experimental observations can be naturally linked to Markov models with the dynamical fingerprint theory, by which experimentally observable timescales can be equipped with an understanding of the structural rearrangement processes that take place at these timescales. The concepts of this paper are illustrated by a simple kinetic model of protein folding.     

Microsecond timescale MD simulations at the transition state of PmHMGR predict remote allosteric residues

Understanding the mechanisms of enzymatic catalysis requires a detailed understanding of the complex interplay of structure and dynamics of large systems that is a challenge for both experimental and computational approaches. More importantly, the computational demands of QM/MM simulations mean that the dynamics of the reaction can only be considered on a timescale of nanoseconds even though the conformational changes needed to reach the catalytically active state happen on a much slower timescale. Here we demonstrate an alternative approach that uses transition state force fields (TSFFs) derived by the quantum-guided molecular mechanics (Q2MM) method that provides a consistent treatment of the entire system at the classical molecular mechanics level and allows simulations at the microsecond timescale. Application of this approach to the second hydride transfer transition state of HMG-CoA reductase from Pseudomonas mevalonii (PmHMGR) identified three remote residues, R396, E399 and L407, (15–27 Å away from the active site) that have a remote dynamic effect on enzyme activity. The predictions were subsequently validated experimentally via site-directed mutagenesis. These results show that microsecond timescale MD simulations of transition states are possible and can predict rather than just rationalize remote allosteric residues.

Transition state force fields enable MD simulations at the transition state of HMGCoA reductase that sample the transition state ensemble on the μs timescale to identify remote residues that affect the reaction rate.     

Simulation studies of the fidelity of biomolecular structure ensemble recreation

We examine the ability of Bayesian methods to recreate structural ensembles for partially folded molecules from averaged data. Specifically we test the ability of various algorithms to recreate different transition state ensembles for folding proteins using a multiple replica simulation algorithm using input from "gold standard" reference ensembles that were first generated with a Go-like Hamiltonian having nonpairwise additive terms. A set of low resolution data, which function as the "experimental" phi values, were first constructed from this reference ensemble. The resulting phi values were then treated as one would treat laboratory experimental data and were used as input in the replica reconstruction algorithm. The resulting ensembles of structures obtained by the replica algorithm were compared to the gold standard reference ensemble, from which those "data" were, in fact, obtained. It is found that for a unimodal transition state ensemble with a low barrier, the multiple replica algorithm does recreate the reference ensemble fairly successfully when no experimental error is assumed. The Kolmogorov-Smirnov test as well as principal component analysis show that the overlap of the recovered and reference ensembles is significantly enhanced when multiple replicas are used. Reduction of the multiple replica ensembles by clustering successfully yields subensembles with close similarity to the reference ensembles. On the other hand, for a high barrier transition state with two distinct transition state ensembles, the single replica algorithm only samples a few structures of one of the reference ensemble basins. This is due to the fact that the phi values are intrinsically ensemble averaged quantities. The replica algorithm with multiple copies does sample both reference ensemble basins. In contrast to the single replica case, the multiple replicas are constrained to reproduce the average phi values, but allow fluctuations in phi for each individual copy. These fluctuations facilitate a more faithful sampling of the reference ensemble basins. Finally, we test how robustly the reconstruction algorithm can function by introducing errors in phi comparable in magnitude to those suggested by some authors. In this circumstance we observe that the chances of ensemble recovery with the replica algorithm are poor using a single replica, but are improved when multiple copies are used. A multimodal transition state ensemble, however, turns out to be more sensitive to large errors in phi (if appropriately gauged) and attempts at successful recreation of the reference ensemble with simple replica algorithms can fall short.     

A Smoluchowski model of crystallization dynamics of small colloidal clusters

We investigate the dynamics of colloidal crystallization in a 32-particle system at a fixed value of interparticle depletion attraction that produces coexisting fluid and solid phases. Free energy landscapes (FELs) and diffusivity landscapes (DLs) are obtained as coefficients of 1D Smoluchowski equations using as order parameters either the radius of gyration or the average crystallinity. FELs and DLs are estimated by fitting the Smoluchowski equations to Brownian dynamics (BD) simulations using either linear fits to locally initiated trajectories or global fits to unbiased trajectories using Bayesian inference. The resulting FELs are compared to Monte Carlo Umbrella Sampling results. The accuracy of the FELs and DLs for modeling colloidal crystallization dynamics is evaluated by comparing mean first-passage times from BD simulations with analytical predictions using the FEL and DL models. While the 1D models accurately capture dynamics near the free energy minimum fluid and crystal configurations, predictions near the transition region are not quantitatively accurate. A preliminary investigation of ensemble averaged 2D order parameter trajectories suggests that 2D models are required to capture crystallization dynamics in the transition region.     

Link between allosteric signal transduction and functional dynamics in a multisubunit enzyme: S-adenosylhomocysteine hydrolase

S-adenosylhomocysteine hydrolase (SAHH), a cellular enzyme that plays a key role in methylation reactions including those required for maturation of viral mRNA, is an important drug target in the discovery of antiviral agents. While targeting the active site is a straightforward strategy of enzyme inhibition, evidence of allosteric modulation of active site in many enzymes underscores the molecular origin of signal transduction. Information of co-evolving sequences in SAHH family and the key residues for functional dynamics that can be identified using native topology of the enzyme provide glimpses into how the allosteric signaling network, dispersed over the molecular structure, coordinates intra- and intersubunit conformational dynamics. To study the link between the allosteric communication and functional dynamics of SAHHs, we performed Brownian dynamics simulations by building a coarse-grained model based on the holo and ligand-bound structures. The simulations of ligand-induced transition revealed that the signal of intrasubunit closure dynamics is transmitted to form intersubunit contacts, which in turn invoke a precise alignment of active site, followed by the dimer-dimer rotation that compacts the whole tetrameric structure. Further analyses of SAHH dynamics associated with ligand binding provided evidence of both induced fit and population shift mechanisms and also showed that the transition-state ensemble is akin to the ligand-bound state. Besides the formation of enzyme-ligand contacts at the active site, the allosteric couplings from the residues distal to the active site are vital to the enzymatic function.     

Calculating the free energy of self-assembled structures by thermodynamic integration

We discuss a method for calculating free energy differences between disordered and ordered phases of self-assembling systems utilizing computer simulations. Applying an external, ordering field, we impose a predefined structure onto the fluid in the disordered phase. The structure in the presence of the external, ordering field closely mimics the structure of the ordered phase (in the absence of an ordering field). Self-consistent field theory or density functional theory provides an accurate estimate for choosing the strength of the ordering field. Subsequently, we gradually switch off the external, ordering field and, in turn, increase the control parameter that drives the self-assembly. The free energy difference along this reversible path connecting the disordered and the ordered state is obtained via thermodynamic integration or expanded ensemble simulation techniques. Utilizing Single-Chain-in-Mean-Field simulations of a symmetric diblock copolymer melt we illustrate the method and calculate the free energy difference between the disordered phase and the lamellar structure at an intermediate incompatibility chiN=20. Evidence for the first-order character of the order-disorder transition at fixed volume is presented. The transition is located at chi(ODT)N=13.65+/-0.10 for an invariant degree of polymerization of N=14 884. The magnitude of the shift of the transition from the mean field prediction qualitatively agrees with other simulations.     

Structural transitions and melting in LJ(74-78) Lennard-Jones clusters from adaptive exchange Monte Carlo simulations

Phase changes in Lennard-Jones (LJ) clusters containing between 74 and 78 atoms are investigated by means of exchange Monte Carlo simulations in the canonical ensemble. The replica temperatures are self-adapted to facilitate the convergence. Although the 74- and 78-atom clusters have icosahedral global minima, the clusters with 75-77 atoms have decahedral ground-state structures and they undergo a structural transition to icosahedral minima before melting. The structural transitions are characterized by quenching and by looking at the Q4 and Q6 orientational bond order parameters. The transition temperatures are estimated to be 0.114, 0.065, and 0.074 reduced units for LJ75, LJ76, and LJ77, respectively. These values, their ordering and the associated latent heats are compared with other estimates based on the harmonic superposition approach.