Path Sampling Simulations Research Articles

Advantages of a distant cellulase catalytic base

The inverting glycoside hydrolase Trichoderma reesei (Hypocrea jecorina) Cel6A is a promising candidate for protein engineering for more economical production of biofuels. Until recently, its catalytic mechanism had been uncertain: The best candidate residue to serve as a catalytic base, Asp-175, is farther from the glycosidic cleavage site than in other glycoside hydrolase enzymes. Recent unbiased transition path sampling simulations revealed the hydrolytic mechanism for this more distant base, employing a water wire; however, it is not clear why the enzyme employs a more distant catalytic base, a highly conserved feature among homologs across different kingdoms. In this work, we describe molecular dynamics simulations designed to uncover how a base with a longer side chain, as in a D175E mutant, affects procession and active site alignment in the Michaelis complex. We show that the hydrogen bond network is tuned to the shorter aspartate side chain, and that a longer glutamate side chain inhibits procession as well as being less likely to adopt a catalytically productive conformation. Furthermore, we draw comparisons between the active site in Trichoderma reesei Cel6A and another inverting, processive cellulase to deduce the contribution of the water wire to the overall enzyme function, revealing that the more distant catalytic base enhances product release. Our results can inform efforts in the study and design of enzymes by demonstrating how counterintuitive sacrifices in chemical reactivity can have worthwhile benefits for other steps in the catalytic cycle.

Reaction Coordinate, Free Energy, and Rate of Intramolecular Proton Transfer in Human Carbonic Anhydrase II.

The role of structure and dynamics of an enzyme has been investigated at three different stages of its function including the chemical event it catalyzes. A one-pot computational method has been designed for each of these stages on the basis of classical and/or quantum mechanical-molecular mechanical molecular dynamics and transition path sampling simulations. For a pair of initial and final states A and B separated by a high free-energy barrier, using a two-stage selection process, several collective variables (CVs) are identified that can delineate A and B. However, these CVs are found to exhibit strong cross-coupling over the transition paths. A set of mutually orthogonal order parameters is then derived from these CVs and an optimal reaction coordinate, r, determined applying half-trajectory likelihood maximization along with a Bayesian information criterion. The transition paths are also used to project the multidimensional free energy surface and barrier crossing dynamics along r. The proposed scheme has been applied to the rate-determining intramolecular proton transfer reaction of the well-known enzyme human carbonic anhydrase II. The potential of mean force, F( r), in the absence of the chemical step is found to reproduce earlier results on the equilibrium population of two side-chain orientations of key residue His-64. Estimation of rate constants, k, from mean first passage times for the three different stages of catalysis shows that the rate-determining step of intramolecular proton transfer occurs with k ≃ 1.0 × 106 s-1, in close agreement with known experimental results.

Atomistic insight into the non-classical nucleation mechanism during solidification in Ni.

Nucleation is a key step during crystallization, but a complete understanding of the fundamental atomistic processes remains elusive. We investigate the mechanism of nucleation during solidification in nickel for various undercoolings using transition path sampling simulations. The temperature dependence of the free energy barriers and rate constants that we obtain is consistent with the predictions of classical nucleation theory and experiments. However, our analysis of the transition path ensemble reveals a mechanism that deviates from the classical picture of nucleation: the growing solid clusters have predominantly non-spherical shapes and consist of face-centered-cubic and random hexagonal-close-packed coordinated atoms surrounded by a cloud of prestructured liquid. The nucleation initiates in regions of supercooled liquid that are characterized by a high orientational order with structural features that predetermine the polymorph selection. These results provide atomistic insight not only into the nucleation mechanism of nickel but also into the role of the preordered liquid regions as precursors for crystallization.

Path Sampling Simulations of the Mechanisms and Rates of Transitions between Watson-Crick and Hoogsteen Base Pairing in DNA

DNA duplexes predominantly contain Watson-Crick (WC) base pairs. However, at any time, a non-negligible fraction of base pairs exhibit a different interaction pattern, called Hoogsteen (HG) where the purine base is rotated approximately 180o relative to its orientation in the WC pairing. The conversion from WC to HG alters the conformation of DNA, and may play a role in several processes, including recognition and replication. The transient nature of these transitions hamper thorough experimental investigation. Molecular dynamics simulations complement experimental work by providing insights at very high spatial and temporal resolution. Path sampling methods focus the molecular dynamics effort on the dynamics during a transition, thus avoiding the long waiting times in stable states. We apply two different path sampling techniques to study the transitions between the WC and HG base pairing motifs. By using transition path sampling we observed that the WC to HG conversion can proceed along several mechanistic routes with varying degrees of exposure of the purine. We computed the rate for each process by employing transition interface sampling.

Thermodynamic and kinetic solid-liquid interface properties from transition path sampling.

We perform transition path sampling simulations to determine two of the key quantities in solidification, the solid-liquid interface energy and velocity, in a Lennard-Jones system. Our approach is applicable to a wide range of temperature and pressure conditions, at the melting temperature and out-of-equilibrium. We show that small system sizes are sufficient for good values of interface energies and velocities. The transition path sampling method thus offers an attractive and robust alternative for the evaluation of solid-liquid interface properties.

Triple Isotope Effects Support Concerted Hydride and Proton Transfer and Promoting Vibrations in Human Heart Lactate Dehydrogenase.

Transition path sampling simulations have proposed that human heart lactate dehydrogenase (LDH) employs protein promoting vibrations (PPVs) on the femtosecond (fs) to picosecond (ps) time scale to promote crossing of the chemical barrier. This chemical barrier involves both hydride and proton transfers to pyruvate to form l-lactate, using reduced nicotinamide adenine dinucleotide (NADH) as the cofactor. Here we report experimental evidence from three types of isotope effect experiments that support coupling of the promoting vibrations to barrier crossing and the coincidence of hydride and proton transfer. We prepared the native (light) LDH and a heavy LDH labeled with 13C, 15N, and nonexchangeable 2H (D) to perturb the predicted PPVs. Heavy LDH has slowed chemistry in single turnover experiments, supporting a contribution of PPVs to transition state formation. Both the [4-2H]NADH (NADD) kinetic isotope effect and the D2O solvent isotope effect were increased in dual-label experiments combining both NADD and D2O, a pattern maintained with both light and heavy LDHs. These isotope effects support concerted hydride and proton transfer for both light and heavy LDHs. Although the transition state barrier-crossing probability is reduced in heavy LDH, the concerted mechanism of the hydride-proton transfer reaction is not altered. This study takes advantage of triple isotope effects to resolve the chemical mechanism of LDH and establish the coupling of fs-ps protein dynamics to barrier crossing.

Allosteric Gating Pathways for the Pentameric Ligand-Gated Ion Channel Glic

Pentameric ligand-gated ion channels (pLGICs) play essential roles in the transmission and interpretation of synaptic signals. The recent solution of multiple X-ray structures of pLGICs has allowed detailed examination of the molecular mechanisms of synaptic transduction at the atomic level. We have described the mechanism of activation for the proton-gated ion channel GLIC, shedding light on the allosteric mechanisms for the entire family of pLGICs. Using a large-scale string method to explore pathways connecting open and the closed state structures, we have described the conformational changes in the agonist-binding extracellular domain, the ion-conducting trans-membrane domain, and the gating interface that communicates the change. We have observed switching of key interactions, allowing for the sensing of agonist/proton binding, leading to a discrete change in the channel pore. These changes involve several residues conserved across the pLGIC family. We have employed transition analysis of data, collected during millions of short path sampling simulations, to compute free energy surfaces that suggest multiple possible pathways, quantitatively reproduce allosteric stabilization of the open state at low pH, and reveal intermediates that are promising as target states for new drugs and anesthetics.

Using Bayes formula to estimate rates of rare events in transition path sampling simulations.

Transition path sampling is a method for estimating the rates of rare events in molecular systems based on the gradual transformation of a path distribution containing a small fraction of reactive trajectories into a biased distribution in which these rare trajectories have become frequent. Then, a multistate reweighting scheme is implemented to postprocess data collected from the staged simulations. Herein, we show how Bayes formula allows to directly construct a biased sample containing an enhanced fraction of reactive trajectories and to concomitantly estimate the transition rate from this sample. The approach can remediate the convergence issues encountered in free energy perturbation or umbrella sampling simulations when the transformed distribution insufficiently overlaps with the reference distribution.

Enzyme homologues have distinct reaction paths through their transition states.

Recent studies of the bacterial enzymes EcMTAN and VcMTAN showed that they have different binding affinities for the same transition state analogue. This was surprising given the similarity of their active sites. We performed transition path sampling simulations of both enzymes to reveal the atomic details of the catalytic chemical step, which may be the key for explaining the inhibitor affinity differences. Even though all experimental data would suggest the two enzymes are almost identical, subtle dynamic differences manifest in differences of reaction coordinate, transition state structure, and eventually significant differences in inhibitor binding. Unlike EcMTAN, VcMTAN has multiple distinct transition states, which is an indication that multiple sets of coordinated protein motions can reach a transition state. Reaction coordinate information is only accessible from transition path sampling approaches, since all experimental approaches report averages. Detailed knowledge could have a significant impact on pharmaceutical design.

The water-catalyzed mechanism of the ring-opening reaction of glucose.

The hexopyranose mutarotation is an important focus for carbohydrate chemistry for more than 150 years. The paper describes the results of advanced computational studies aimed at elucidating the ring-opening reaction of glucose. Molecular simulations based on the combination of the DFT method with the molecular dynamics formalism allowed for a detailed insight into the mechanism of the process accompanied by the information of the kinetic and dynamic nature. The results indicate that the process is initiated by deprotonation of the anomeric hydroxyl group by water molecules and the subsequent proton transfer to the ring oxygen atom. The latter event has been identified as a 'bottleneck' of the process triggering the ring cleavage. The most time-consuming steps of the ring-opening reaction are the orientational rearrangements of water molecule(s) participating in the proton transfer(s) and the final extension of the newly-formed aldehyde chain. The orientational preferences of the aldehyde group present in the acyclic form of D-hexopyranoses are responsible for the anomeric equilibrium characteristics.

Reaction coordinate of incipient methane clathrate hydrate nucleation.

Nucleation from solution is a ubiquitous phenomenon with relevance to myriad scientific disciplines, including pharmaceuticals, biomineralization, and disease. One prominent example is the nucleation of clathrate hydrates, multicomponent crystalline inclusion compounds relevant to the energy industry where they block pipelines and also constitute a potential vast energy resource. Despite their importance, the molecular mechanism of incipient hydrate formation remains unknown. Herein, we employ advanced molecular simulation tools (pB histogram, equilibrium path sampling) to provide a statistical-mechanical basis for extracting physical insight into the molecular steps by which clathrates form. Through testing the Mutually Coordinated Guest (MCG) order parameter, we demonstrate that both guest (methane) and host (water) structuring are crucial to accurately describe the nucleation of hydrates and determine a critical nucleus size of MCG-1 = 16 at 255 K and 500 bar. Equipped with a validated (and novel) reaction coordinate, subsequent equilibrium path sampling simulations yield the free energy barrier and nucleation rate. The resulting quantitative nucleation process is described by the MCG clustering mechanism. This constitutes a significant advance in the field of hydrates research, as the fitness of a molecular descriptor has never been statistically verified. More broadly, this work has significance to a wide range of multicomponent nucleation contexts wherein the formation mechanism depends on contributions from both solute and solvent.

Structure, Kinetics, and Thermodynamics of the Aqueous Uranyl(VI) Cation

In this work, molecular simulation techniques were employed to gain insight into the structural, kinetic, and thermodynamic properties of the uranyl(VI) cation (UO2(2+)) in aqueous solution. The simulations made use of an atomistic potential model (force field) derived in this work and based on the model of Guilbaud and Wipff [J. Mol. Struct. (THEOCHEM) 1996 , 366 , 55 - 63]. Reactive flux and thermodynamic integration calculations show that the derived potential model yields predictions for the water exchange rate and free energy of hydration, respectively, that are in agreement with experimental data. The water binding energies, hydration shell structure, and self-diffusion coefficient were also calculated and analyzed. Finally, a combination of metadynamics and transition path sampling simulations was employed to probe the mechanisms of water exchange reactions in the first hydration shell of the uranyl ion. These atomistic simulations indicate, based on two-dimensional free energy surfaces, that water exchanges follow an associative interchange mechanism. The nature and structure of the water exchange transition states were also determined. The improved potential model is expected to lead to more accurate predictions of uranyl adsorption energies at mineral surfaces using potential-based molecular dynamics simulations.

α-helix to β-hairpin transition of human amylin monomer.

The human islet amylin polypeptide is produced along with insulin by pancreatic islets. Under some circumstances, amylin can aggregate to form amyloid fibrils, whose presence in pancreatic cells is a common pathological feature of Type II diabetes. A growing body of evidence indicates that small, early stage aggregates of amylin are cytotoxic. A better understanding of the early stages of the amylin aggregation process and, in particular, of the nucleation events leading to fibril growth could help identify therapeutic strategies. Recent studies have shown that, in dilute solution, human amylin can adopt an α-helical conformation, a β-hairpin conformation, or an unstructured coil conformation. While such states have comparable free energies, the β-hairpin state exhibits a large propensity towards aggregation. In this work, we present a detailed computational analysis of the folding pathways that arise between the various conformational states of human amylin in water. A free energy surface for amylin in explicit water is first constructed by resorting to advanced sampling techniques. Extensive transition path sampling simulations are then employed to identify the preferred folding mechanisms between distinct minima on that surface. Our results reveal that the α-helical conformer of amylin undergoes a transformation into the β-hairpin monomer through one of two mechanisms. In the first, misfolding begins through formation of specific contacts near the turn region, and proceeds via a zipping mechanism. In the second, misfolding occurs through an unstructured coil intermediate. The transition states for these processes are identified. Taken together, the findings presented in this work suggest that the inter-conversion of amylin between an α-helix and a β-hairpin is an activated process and could constitute the nucleation event for fibril growth.

Understanding the nature of “superhard graphite”

Numerous experiments showed that on cold compression graphite transforms into a new superhard and transparent allotrope. Several structures with different topologies have been proposed for this phase. While experimental data are compatible with most of these models, the only way to solve this puzzle is to find which structure is kinetically easiest to form. Using state-of-the-art molecular-dynamics transition path sampling simulations, we investigate kinetic pathways of the pressure-induced transformation of graphite to various superhard candidate structures. Unlike hitherto applied methods for elucidating nature of superhard graphite, transition path sampling realistically models nucleation events necessary for physically meaningful transformation kinetics. We demonstrate that nucleation mechanism and kinetics lead to M-carbon as the final product. W-carbon, initially competitor to M-carbon, is ruled out by phase growth. Bct-C4 structure is not expected to be produced by cold compression due to less probable nucleation and higher barrier of formation.

Insight into the Alternating Access Mechanism of the Sodium Symporter Mhp1 using Path Sampling Simulations

Sodium coupled co-transporters of the LeuT superfamily use an alternating access mechanism to move small molecules across the cell membrane. A key step in the transport cycle involves the transporter undergoing a conformational transition from a state poised to bind extracellular substrates to one that is competent to deliver substrate to the cytoplasm. To investigate this step in semi-atomistic detail, we construct a coarse-grained model of the benzylhydantoin transporter Mhp1 that incorporates experimental structures of the outward- and inward-open states. We obtain a rigorous sampling of the unbiased outward- to inward-facing transition path ensemble (TPE) using the weighted ensemble path-sampling method. Analysis of the TPE reveals two distinct sets of pathways connecting the inward- and outward-facing conformations. The first is consistent with a strict alternating access mechanism, while in the other, decoupling of the inner and outer gates transiently forms a putative permeation pathway through the transporter. Our simulations additionally show that transitions between the outward- and inward-open states result from rigid body motions of the hash motif relative to the substrate bundle, supporting the rocking bundle hypothesis. This methodology provides the groundwork for more chemically detailed investigations of the transport mechanism of superfamily members using fully atomistic simulations.

Waste-recycling Monte Carlo with optimal estimates: Application to free energy calculations in alloys

The estimator proposed recently by Delmas and Jourdain for waste-recycling Monte Carlo achieves variance reduction optimally with respect to a control variate that is evaluated directly using the simulation data. Here, the performance of this estimator is assessed numerically for free energy calculations in generic binary alloys and is compared to those of other estimators taken from the literature. A systematic investigation with varying simulation parameters of a simplified system, the anti-ferromagnetic Ising model, is first carried out in the transmutation ensemble using path-sampling. We observe numerically that (i) the variance of the Delmas-Jourdain estimator is indeed reduced compared to that of other estimators; and that (ii) the resulting reduction is close to the maximal possible one, despite the inaccuracy in the estimated control variate. More extensive path-sampling simulations involving an FeCr alloy system described by a many-body potential additionally show that (iii) gradual transmutations accommodate the atomic frustrations; thus, alleviating the numerical ergodicity issue present in numerous alloy systems and eventually enabling the determination of phase coexistence conditions.

Theoretical Study of the Decomposition of BCl3 Induced by a H Radical

We report on a theoretical study of the gas-phase decomposition of boron trichloride in the presence of hydrogen radicals using ab initio energetic calculations coupled to TST, RRKM, and VTST-VRC kinetic calculations. In particular, we present an addition-elimination mechanism (BCl(3) + H → BHCl(2) + Cl) allowing for a much more rapid consumption of BCl(3) than the direct abstraction reaction (BCl(3) + H → BCl(2) + HCl) considered up to now. At low temperatures, T ≤ 800 K, our results show that a weakly stabilized complex BHCl(3) is formed with a kinetic law compatible with the consumption rate measured in the former experiments. At higher temperatures, this complex is not stable and then easily eliminates a chlorine atom. Our work also shows that a very similar mechanism, involving the same intermediate and sharing the same transition state, allows for the elimination of HCl. A dividing coefficient between these two elimination pathways is obtained from both a potential energy surface based statistical analysis and an ab initio molecular dynamics transition path sampling simulation. It finally allows partitioning of the global consumption rate of BCl(3) in terms of the formation of (i) BHCl(3), (ii) BHCl(2) + Cl through a H addition/Cl elimination mechanism, (iii) BCl(2) + HCl through a H addition/HCl elimination mechanism, and (iv) BCl(2) + HCl through direct abstraction.

Molecular pathways in DNA-DNA hybridization of surface-bound oligonucleotides

A molecular model of DNA, that has been shown to reproduce hybridization, is used here to study the transition from two single-stranded molecules to a double-helical complex in situations where one strand is attached to a surface. Transition Path Sampling simulations reveal that hybridization is initiated by forming specific contacts between distinct bases of the molecule. Molecules with repetitive sequences form double helices by sliding along each other after having formed initial contacts. For random sequences, it is shown that the pathway that molecules follow to hybridization depends on the linker length as well as the temperature, and is significantly different from that observed in the bulk.

Recent advances in transition path sampling: accurate reaction coordinates, likelihood maximisation and diffusive barrier-crossing dynamics

Because transition states are too rare and transient to observe in experiments, simulations are among the best sources of molecular-level insight on transition states and reaction mechanisms. A first step in many rare-event simulations is to identify a reaction coordinate. For reactions that break and make strong bonds, the reaction coordinate is the unstable eigenmode at a saddle point on the potential energy surface. However, processes such as nucleation and protein folding may disrupt and reorganise thousands of interactions, so identifying a reaction coordinate is a major challenge. An accurate reaction coordinate should result in a free energy profile that is consistent with the projected dynamics. Ten years ago, this intuitive requirement led to committor analysis, a trial-and-error procedure for testing putative reaction coordinates. Since then, our ability to identify accurate reaction coordinates has dramatically improved. First, new versions of committor analysis are quantitative and more efficient. Second, likelihood maximisation can systematically and efficiently identify accurate reaction coordinates from thousands of candidates using only the shooting point data from a path sampling simulation. Finally, the new aimless shooting version of transition path sampling retains a high sampling efficiency even for systems with highly diffusive barrier-crossing dynamics.

Nonlinear reaction coordinate analysis in the reweighted path ensemble

We present a flexible nonlinear reaction coordinate analysis method for the transition path ensemble based on the likelihood maximization approach developed by Peters and Trout [J. Chem. Phys. 125, 054108 (2006)]. By parametrizing the reaction coordinate by a string of images in a collective variable space, we can optimize the likelihood that the string correctly models the committor data obtained from a path sampling simulation. The collective variable space with the maximum likelihood is considered to contain the best description of the reaction. The use of the reweighted path ensemble [J. Rogal et al., J. Chem. Phys. 133, 174109 (2010)] allows a complete reaction coordinate description from the initial to the final state. We illustrate the method on a z-shaped two-dimensional potential. While developed for use with path sampling, this analysis method can also be applied to regular molecular dynamics trajectories.