On the calculation of reaction rate constants in the transition path ensemble

Substrate Turnover Dynamics Guide Ketol-Acid Reductoisomerase Redesign for Increased Specific Activity

The task of adapting enzymes for specific applications is often hampered by our incomplete ability to tune and tailor catalytic functions, particularly when seeking increased activity. Here, we develop and demonstrate a rational approach to address this challenge, applied to ketol-acid reductoisomerase (KARI), which has uses in industrial-scale isobutanol production. While traditional structure-based computational enzyme redesign strategies typically focus on the enzyme-bound ground state (GS) and transition state (TS), we postulated that additionally treating the underlying dynamics of complete turnover events that connect and pass through both states could further elucidate the structural properties affecting catalysis and help identify mutations that lead to increased catalytic activity. To examine the dynamics of substrate conversion with atomistic detail, we adapted and applied computational methods based on path sampling techniques to gather thousands of QM/MM simulations of attempted substrate turnover events by KARI: both productive (reactive) and unproductive (nonreactive) attempts. From these data, machine learning models were constructed and used to identify specific conformational features (interatomic distances, angles, and torsions) associated with successful, productive catalysis. Multistate protein redesign techniques were then used to select mutations that stabilized reactive-like structures over nonreactive-like ones while also meeting additional criteria consistent with enhanced specific activity. This procedure resulted in eight high-confidence enzyme mutants with a significant improvement in calculated specific activity relative to wild type (WT), with the fastest variant's increase in calculated k cat being (2 ± 1) × 104-fold. Collectively, these results suggest that introducing mutations designed to increase the population of reaction-promoting conformations of the enzyme-substrate complex before it reaches the barrier can provide an effective approach to engineering improved enzyme catalysts.

Sampling a Rare Protein Transition Using Quantum Annealing

Simulating spontaneous structural rearrangements in macromolecules with classical molecular dynamics is an outstanding challenge. Conventional supercomputers can access time intervals of up to tens of μs, while many key events occur on exponentially longer time scales. Path sampling techniques have the advantage of focusing the computational power on barrier-crossing trajectories, but generating uncorrelated transition paths that explore diverse conformational regions remains a problem. We employ a hybrid path-sampling paradigm that addresses this issue by generating trial transition paths using a quantum annealing (QA) machine. We first employ a classical computer to perform an uncharted exploration of the conformational space. The data set generated in this exploration is then postprocessed using a path integral-based method to yield a coarse-grained network representation of the reactive kinetics. By resorting to a quantum annealer, quantum superposition can be exploited to encode all of the transition pathways in the initial quantum state, thus potentially solving the path exploration problem. Furthermore, each QA cycle yields a completely uncorrelated trial trajectory. We previously validated this scheme on a prototypically simple transition, which could be extensively characterized on a desktop computer. Here, we scale up in complexity and perform an all-atom simulation of a protein conformational transition that occurs on the millisecond time scale, obtaining results that match those of the Anton special-purpose supercomputer. Despite limitations due to the available quantum annealers, our study highlights how realistic biomolecular simulations provide potentially impactful new ground for applying, testing, and advancing quantum technologies.

Diffusion in a two-dimensional energy landscape in the presence of dynamical correlations and validity of random walk model

Diffusion in a multidimensional energy surface with minima and barriers is a problem of importance in statistical mechanics and it also has wide applications, such as protein folding. To understand it in such a system, we carry out theory and simulations of a tagged particle moving on a two-dimensional periodic potential energy surface, both in the presence and absence of noise. Langevin dynamics simulations at multiple temperatures are carried out to obtain the diffusion coefficient of a solute particle. Friction is varied from zero to large values. Diffusive motion emerges in the limit of a long time, even in the absence of noise. Noise destroys the correlations and increases diffusion at small friction. Diffusion thus exhibits a nonmonotonic friction dependence at the intermediate value of the damping, ultimately converging to our theoretically predicted value. The latter is obtained using the well-established relationship between diffusion and random walk. An excellent agreement is obtained between theory and simulations in the high-friction limit but not so in the intermediate regime. We explain the deviation in the low- to intermediate-friction regime using the modified random walk theory. The rate of escape from one cell to another is obtained from the multidimensional rate theory of Langer. We find that enhanced dimensionality plays an important role. To quantify the effects of noise on the potential-imposed coherence on the trajectories, we calculate the Lyapunov exponent. At small friction values, the Lyapunov exponent mimics the friction dependence of the rate.

Supervised learning and the finite-temperature string method for computing committor functions and reaction rates

A central object in the computational studies of rare events is the committor function. Though costly to compute, the committor function encodes complete mechanistic information of the processes involving rare events, including reaction rates and transition-state ensembles. Under the framework of transition path theory, Rotskoff et al. [ Proceedings of the 2nd Mathematical and Scientific Machine Learning Conference, Proceedings of Machine Learning Research (PLMR, 2022), Vol. 145, pp. 757–780] proposes an algorithm where a feedback loop couples a neural network that models the committor function with importance sampling, mainly umbrella sampling, which collects data needed for adaptive training. In this work, we show additional modifications are needed to improve the accuracy of the algorithm. The first modification adds elements of supervised learning, which allows the neural network to improve its prediction by fitting to sample-mean estimates of committor values obtained from short molecular dynamics trajectories. The second modification replaces the committor-based umbrella sampling with the finite-temperature string (FTS) method, which enables homogeneous sampling in regions where transition pathways are located. We test our modifications on low-dimensional systems with non-convex potential energy where reference solutions can be found via analytical or finite element methods, and show how combining supervised learning and the FTS method yields accurate computation of committor functions and reaction rates. We also provide an error analysis for algorithms that use the FTS method, using which reaction rates can be accurately estimated during training with a small number of samples. The methods are then applied to a molecular system in which no reference solution is known, where accurate computations of committor functions and reaction rates can still be obtained.

Nonadiabatic transition paths from quantum jump trajectories

We present a means of studying rare reactive pathways in open quantum systems using transition path theory and ensembles of quantum jump trajectories. This approach allows for the elucidation of reactive paths for dissipative, nonadiabatic dynamics when the system is embedded in a Markovian environment. We detail the dominant pathways and rates of thermally activated processes and the relaxation pathways and photoyields following vertical excitation in a minimal model of a conical intersection. We find that the geometry of the conical intersection affects the electronic character of the transition state as defined through a generalization of a committor function for a thermal barrier crossing event. Similarly, the geometry changes the mechanism of relaxation following a vertical excitation. Relaxation in models resulting from small diabatic coupling proceeds through pathways dominated by pure dephasing, while those with large diabatic coupling proceed through pathways limited by dissipation. The perspective introduced here for the nonadiabatic dynamics of open quantum systems generalizes classical notions of reactive paths to fundamentally quantum mechanical processes.

Potential Energy Weighted Reactive Flux and Total Rate of Change of Potential Energy: Theory and Illustrative Applications

Reactive flux can be largely nonzero in a nonequilibrium ensemble of trajectories and provide insightful information for reactive transitions from the reactant state to the product state. Based on the reactive flux, a theoretical framework is proposed here for two quantities, the potential energy weighted reactive flux and the total rate of change of potential energy, which are useful for the identification of the mechanism from a nonequilibrium ensemble. From such quantities, two multidimensional free-energy analogues can be derived in the subspace of collective variables and they are equivalent in the regions where the reactive flux is divergence-free. These free-energy analogues are assumed to be closely related to the free energy in the subspace of collective variables, and they are reduced in the one-dimensional case to be the ensemble average of the potential energy weighted with reactive flux intensity, which was proposed recently [Li, W. J. Phys. Chem. A 2022, DOI: 10.1021/acs.jpca.2c04130] and could be decomposed into energy components at the per-coordinate level. In the subspace of collective variables, the decomposition of the multidimensional free-energy analogues at the per-coordinate level is theoretically possible and is numerically difficult to be calculated. Interestingly, the total rate of change of potential energy is able to identify the location of the transition state ensemble or the stochastic separatrix, in addition to the locations of the reactant and product states. The total rate of change of potential energy can be decomposed at the per-coordinate level, and its components can quantify the contribution of a coordinate to the reactive transition in the subspace of collective variables. We then illustrated the main insights and objects that can be provided by the approach in the applications to a two-dimensional system with various diffusion anisotropies and the alanine peptide in vacuum in various nonequilibrium ensembles of short trajectories, from which the results were found to be consistent.

Energy Decomposition along the Reaction Coordinate: Theory and Applications to Nonequilibrium Ensembles of Trajectories

A theoretical framework is proposed for an energy decomposition scheme along the reaction coordinate, in which the ensemble average of the potential energy weighted with reactive flux intensity is decomposed into energy components at the per-coordinate level. The decomposed energy quantity is demonstrated to be closely related to the free energy along the reaction coordinate, and its connection to the emergent potential energy is provided. In the application to alanine dipeptide under vacuum, illustrative calculations were performed in three nonequilibrium ensembles of trajectories: (1) transition path ensemble sampled with transition path sampling; (2) ensemble of short trajectories initiated from configurations around the transition-state region; and (3) ensemble of short trajectories shooting from configurations in several transition paths. The energy components on each coordinate were found to be consistent among the three ensembles of trajectories, indicating a broad applicability of the approach in biomolecular studies. In addition, the free energies along an optimized reaction coordinate obtained with these nonequilibrium ensembles were largely overlapped with a reference free energy calculated from a long equilibrium trajectory.

Sampling rare conformational transitions with a quantum computer

Structural rearrangements play a central role in the organization and function of complex biomolecular systems. In principle, Molecular Dynamics (MD) simulations enable us to investigate these thermally activated processes with an atomic level of resolution. In practice, an exponentially large fraction of computational resources must be invested to simulate thermal fluctuations in metastable states. Path sampling methods focus the computational power on sampling the rare transitions between states. One of their outstanding limitations is to efficiently generate paths that visit significantly different regions of the conformational space. To overcome this issue, we introduce a new algorithm for MD simulations that integrates machine learning and quantum computing. First, using functional integral methods, we derive a rigorous low-resolution spatially coarse-grained representation of the system's dynamics, based on a small set of molecular configurations explored with machine learning. Then, we use a quantum annealer to sample the transition paths of this low-resolution theory. We provide a proof-of-concept application by simulating a benchmark conformational transition with all-atom resolution on the D-Wave quantum computer. By exploiting the unique features of quantum annealing, we generate uncorrelated trajectories at every iteration, thus addressing one of the challenges of path sampling. Once larger quantum machines will be available, the interplay between quantum and classical resources may emerge as a new paradigm of high-performance scientific computing. In this work, we provide a platform to implement this integrated scheme in the field of molecular simulations.

Time-Lagged Flux in the Transition Path Ensemble: Flux Maximization and Relation to Transition Path Theory

The transition path ensemble is of special interest in reaction coordinate identification as it consists of reactive trajectories that start from the reactant state and end in the product one. As a theoretical framework for describing the transition path ensemble, the transition path theory has been introduced more than 10 years ago, and so far, its applications have only been illustrated in several low-dimensional systems. Given the transition path ensemble, expressions for calculating flux, current (a vector field), and principal curves are derived here in the space of collective variables from the transition path theory, and they are applicable to time series obtained from molecular dynamics simulations of high-dimensional systems, i.e., the position coordinates as a function of time in the transition path ensemble. The connection of the transition path theory is made to a density-weighted average flux, a quantity proposed in a previous work to appraise the relevance of a coordinate to the reaction coordinate [Li, W. J. Chem. Phys. 2022, 156, 054117]. Most importantly, as an extension of the existing quantities, time-lagged quantities such as flux and current are also proposed. The main insights and objects provided by these time-lagged quantities are illustrated in the application to the alanine peptide in vacuum.

Direct Evaluation of Rare Events in Active Matter from Variational Path Sampling

Active matter represents a broad class of systems that evolve far from equilibrium due to the local injection of energy. Like their passive analogs, transformations between distinct metastable states in active matter proceed through rare fluctuations; however, their detailed balance violating dynamics renders these events difficult to study. Here, we present a simulation method for evaluating the rate and mechanism of rare events in generic nonequilibrium systems and apply it to study the conformational changes of a passive solute in an active fluid. The method employs a variational optimization of a control force that renders the rare event a typical one, supplying an exact estimate of its rate as a ratio of path partition functions. Using this method we find that increasing activity in the active bath can enhance the rate of conformational switching of the passive solute in a manner consistent with recent bounds from stochastic thermodynamics.

Collective variable-based enhanced sampling and machine learning

ABSTRACT

Collective variable-based enhanced sampling methods have been widely used to study thermodynamic properties of complex systems. Efficiency and accuracy of these enhanced sampling methods are affected by two factors: constructing appropriate collective variables for enhanced sampling and generating accurate free energy surfaces. Recently, many machine learning techniques have been developed to improve the quality of collective variables and the accuracy of free energy surfaces. Although machine learning has achieved great successes in improving enhanced sampling methods, there are still many challenges and open questions. In this perspective, we shall review recent developments on integrating machine learning techniques and collective variable-based enhanced sampling approaches. We also discuss challenges and future research directions including generating kinetic information, exploring high-dimensional free energy surfaces, and efficiently sampling all-atom configurations.

GRAPHIC ABSTRACT

Homogenous nucleation rate of CO2 hydrates using transition interface sampling

Carbon dioxide and water can form solid clathrate structures in which water cages encapsulate the gas molecules. Such hydrates have sparked much interest due to their possible application in CO₂ sequestration. How the solid structure forms exactly from the liquid phase via a homogenous nucleation process is still poorly understood. This nucleation event is rare on the molecular timescale even under moderate undercooling or supersaturation conditions because of the large free energy barrier toward crystallization, rendering a brute force simulation of hydrate nucleation unfeasible for moderate undercooling or supersaturation. Here, we perform transition interface sampling simulations to quantify the homogenous nucleation rate for CO₂ hydrate formation using accurate atomistic force fields at 500 bars for three different temperatures between 260 and 273 K. Collecting more than 100 000 pathways comprising roughly two milliseconds of simulation time, we computed a nucleation rate in the amorphous phase of ∼10²¹ nuclei s^-1 cm^-3 for a temperature of 260 K and a rate of ∼10¹² nuclei s^-1 cm^-3 for a temperature of 265 K. For a temperature of 273 K, we find that the hydrate forms an sI crystalline phase with a rate of order of ∼10¹ nuclei s^-1 cm^-3. We compare these rates to classical nucleation theory estimates as well as experiments, and to nucleation rate estimates for methane hydrates and discuss possible causes of the observed differences. Our findings shed light on the kinetics of this important clathrate and should assist in future hydrate formation investigation.

Homogeneous nucleation of sheared liquids: advances and insights from simulations and theory

One of the most ubiquitous and technologically important phenomena in nature is the nucleation of homogeneous flowing systems. The microscopic effects of shear on a nucleating system are still imperfectly understood, although in recent years a consistent picture has emerged. The opposing effects of shear can be split into two major contributions for simple atomic and molecular liquids: increase of the energetic cost of nucleation, and enhancement of the kinetics. In this perspective, we describe the latest computational and theoretical techniques which have been developed over the past two decades. We collate and unify the overarching influences of shear, temperature, and supersaturation on the process of homogeneous nucleation. Experimental techniques and capabilities are discussed, against the backdrop of results from simulations and theory. Although we primarily focus on simple systems, we also touch upon the sheared nucleation of more complex systems, including glasses and polymer melts. We speculate on the promising directions and possible advances that could come to fruition in the future.

Estimating reciprocal partition functions to enable design space sampling

Reaction rates are a complicated function of molecular interactions, which can be selected from vast chemical design spaces. Seeking the design that optimizes a rate is a particularly challenging problem since the rate calculation for any one design is itself a difficult computation. Toward this end, we demonstrate a strategy based on transition path sampling to generate an ensemble of designs and reactive trajectories with a preference for fast reaction rates. Each step of the Monte Carlo procedure requires a measure of how a design constrains molecular configurations, expressed via the reciprocal of the partition function for the design. Although the reciprocal of the partition function would be prohibitively expensive to compute, we apply Booth's method for generating unbiased estimates of a reciprocal of an integral to sample designs without bias. A generalization with multiple trajectories introduces a stronger preference for fast rates, pushing the sampled designs closer to the optimal design. We illustrate the methodology on two toy models of increasing complexity: escape of a single particle from a Lennard-Jones potential well of tunable depth and escape from a metastable tetrahedral cluster with tunable pair potentials.

Improving double-ended transition state searches for soft-matter systems

Transitions between different stable configurations of biomolecules are important in understanding disease mechanisms, structure-function relations, and novel molecular-scale engineering. The corresponding pathways can be characterized efficiently using geometry optimization schemes based on double-ended transition state searches. An interpolation is first constructed between the known states and then refined, yielding a band that contains transition state candidates. Here, we analyze an example where various interpolation schemes lead to bands with a single step transition, but the correct pathway actually proceeds via an intervening, low-energy minimum. We compare a number of different interpolation schemes for this problem. We systematically alter the number of discrete images in the interpolations and the spring constants used in the optimization and test two schemes for adjusting the spring constants and image distribution, resulting in a total of 2760 different connection attempts. Our results confirm that optimized bands are not necessarily a good description of the transition pathways in themselves, and further refinement to actually converge transition states and establish their connectivity is required. We see an improvement in the optimized bands if we employ the adjustment of spring constants with doubly-nudged elastic band and a smaller improvement from the image redistribution. The example we consider is representative of numerous cases we have encountered in a wide variety of molecular and condensed matter systems.

Revisiting the Volmer-Heyrovský mechanism of hydrogen evolution on a nitrogen doped carbon nanotube: constrained molecular dynamics versus the nudged elastic band method

Density functional theory (DFT) based computational electrochemistry has the potential to serve as a tool with predictive power in the rational development and screening of electrocatalysts for renewable energy technologies. It is, however, of paramount importance that simulations are conducted rigorously at a level of theory that is sufficiently accurate in order to obtain physicochemically sensible results. Herein, we present a comparative study of the performance of the static climbing image nudged elastic band method (CI-NEB) vs. DFT based constrained molecular dynamics simulations with thermodynamic integration in estimating activation and reaction (free) energies of the Volmer-Heyrovský mechanism on a nitrogen doped carbon nanotube. Due to cancellation of errors within the CI-NEB calculations, static and dynamic activation barriers are observed to be surprisingly similar, while a substantial decrease in reaction energies is seen upon incorporation of solvent dynamics. This finding is attributed to two competing effects; (1) solvent reorganization that stabilizes the transition and, in particular, the product states with respect to the reactant state and (2) destabilizing entropic contributions due to solvent fluctuations. Our results highlight the importance of explicitly sampling the interfacial solvent dynamics when studying hydrogen evolution at solid-liquid interfaces.

Simulating ion channel activation mechanisms using swarms of trajectories

Atomic-level studies of protein activity represent a significant challenge as a result of the complexity of conformational changes occurring on wide-ranging timescales, often greatly exceeding that of even the longest simulations. A prime example is the elucidation of protein allosteric mechanisms, where localized perturbations transmit throughout a large macromolecule to generate a response signal. For example, the conversion of chemical to electrical signals during synaptic neurotransmission in the brain is achieved by specialized membrane proteins called pentameric ligand-gated ion channels. Here, the binding of a neurotransmitter results in a global conformational change to open an ion-conducting pore across the nerve cell membrane. X-ray crystallography has produced static structures of the open and closed states of the proton-gated GLIC pentameric ligand-gated ion channel protein, allowing for atomistic simulations that can uncover changes related to activation. We discuss a range of enhanced sampling approaches that could be used to explore activation mechanisms. In particular, we describe recent application of an atomistic string method, based on Roux's "swarms of trajectories" approach, to elucidate the sequence and interdependence of conformational changes during activation. We illustrate how this can be combined with transition analysis and Brownian dynamics to extract thermodynamic and kinetic information, leading to understanding of what controls ion channel function. © 2019 Wiley Periodicals, Inc.

Transition path sampling for non-equilibrium dynamics without predefined reaction coordinates

We develop two novel transition path sampling (TPS) algorithms for harvesting ensembles of rare event trajectories using non-equilibrium dynamics. These methods have the advantage that no predefined reaction coordinate is needed. Instead, an instantaneous reaction coordinate is based on the current path. Constituting a Monte Carlo random walk in trajectory space, the algorithms can be viewed as bridging between the original TPS methodology and the Rosenbluth based forward flux sampling methodology. We illustrate the new methods on toy models undergoing equilibrium and non-equilibrium dynamics, including an active Brownian particle system. For the latter, we find that transitions between steady states occur via states that are locally ordered but globally disordered.

Electrode polarization effects on interfacial kinetics of ionic liquid at graphite surface: An extended lagrangian-based constant potential molecular dynamics simulation study

Computational models including electrode polarization can be essential to study electrode/electrolyte interfacial phenomena more realistically. We present here a constant-potential classical molecular dynamics simulation method based on the extended Lagrangian formulation where the fluctuating electrode atomic charges are treated as independent dynamical variables. The method is applied to a graphite/ionic liquid system for the validation and the interfacial kinetics study. While the correct adiabatic dynamics is achieved with a sufficiently small fictitious mass of charge, static properties have been shown to be almost insensitive to the fictitious mass. As for the kinetics study, electrical double layer (EDL) relaxation and ion desorption from the electrode surface are considered. We found that the polarization slows EDL relaxation greatly whereas it has little impact on the ion desorption kinetics. The findings suggest that the polarization is essential to estimate the kinetics in nonequilibrium processes, not in equilibrium. © 2019 Wiley Periodicals, Inc.

Past-future information bottleneck for sampling molecular reaction coordinate simultaneously with thermodynamics and kinetics

The ability to rapidly learn from high-dimensional data to make reliable bets about the future is crucial in many contexts. This could be a fly avoiding predators, or the retina processing gigabytes of data to guide human actions. In this work we draw parallels between these and the efficient sampling of biomolecules with hundreds of thousands of atoms. For this we use the Predictive Information Bottleneck framework used for the first two problems, and re-formulate it for the sampling of biomolecules, especially when plagued with rare events. Our method uses a deep neural network to learn the minimally complex yet most predictive aspects of a given biomolecular trajectory. This information is used to perform iteratively biased simulations that enhance the sampling and directly obtain associated thermodynamic and kinetic information. We demonstrate the method on two test-pieces, studying processes slower than milliseconds, calculating free energies, kinetics and critical mutations.

Efficient sampling of rare events in all-atom molecular dynamics simulations remains a challenge. Here, the authors adapt the Predictive Information Bottleneck framework to sample biomolecular structure and dynamics through iterative rounds of biased simulations and deep learning.

OpenPathSampling: A Python Framework for Path Sampling Simulations. 1. Basics

Transition path sampling techniques allow molecular dynamics simulations of complex systems to focus on rare dynamical events, providing insight into mechanisms and the ability to calculate rates inaccessible by ordinary dynamics simulations. While path sampling algorithms are conceptually as simple as importance sampling Monte Carlo, the technical complexity of their implementation has kept these techniques out of reach of the broad community. Here, we introduce an easy-to-use Python framework called OpenPathSampling (OPS) that facilitates path sampling for (bio)molecular systems with minimal effort and yet is still extensible. Interfaces to OpenMM and an internal dynamics engine for simple models are provided in the initial release, but new molecular simulation packages can easily be added. Multiple ready-to-use transition path sampling methodologies are implemented, including standard transition path sampling (TPS) between reactant and product states and transition interface sampling (TIS) and its replica exchange variant (RETIS), as well as recent multistate and multiset extensions of transition interface sampling (MSTIS, MISTIS). In addition, tools are provided to facilitate the implementation of new path sampling schemes built on basic path sampling components. In this paper, we give an overview of the design of this framework and illustrate the simplicity of applying the available path sampling algorithms to a variety of benchmark problems.

Advances in milestoning. II. Calculating time-correlation functions from milestoning using stochastic path integrals

In the milestoning framework, and more generally in related transition interface sampling schemes, one significantly enhances the calculation of relaxation rates for complex equilibrium kinetics from molecular dynamics simulations between the milestones or interfaces. The goal of the present paper is to advance milestoning applications into the realm of non-equilibrium statistical mechanics, in particular, to calculate entire time correlation functions. In order to accomplish this, we introduce a novel methodology for obtaining the flux through a given milestone configuration as a function of both time and initial configuration and build upon it with a novel formalism describing autocorrelation for Langevin motion in a discrete configuration space. The method is then applied to three different test systems: a harmonic oscillator, which we solve analytically, a two-well potential, which is solved numerically, and an atomistic molecular dynamics simulation of alanine dipeptide.

Solid-phase nucleation free-energy barriers in truncated cubes: interplay of localized orientational order and facet alignment

The nucleation of ordered phases from the bulk isotropic phase of octahedron-like particles has been studied via Monte Carlo simulations and umbrella sampling. In particular, selected shapes that form ordered (plastic) phases with various symmetries (cubic and tetragonal) are chosen to unveil trends in the free-energy barrier heights (ΔG*'s) associated with disorder to order transitions. The shapes studied in this work have truncation parameter (s) values of 0.58, 0.75, 0.8 and 1. The case of octahedra (s = 1.0) is studied to provide a counter-example where the isotropic phase nucleates directly into a (Minkowski) crystal phase rather than a rotator phase. The simulated ΔG*'s for these systems are compared with those previously reported for hard spheres and truncated cubes with s = 0.5 (cuboctahedra, CO) and s = 2/3 (truncated octahedra, TO). The comparison shows that, for comparable degrees of supersaturation, all rotator phases nucleate with smaller ΔG*'s than that of the hard sphere crystal, whereas the octahedral crystal nucleates with a larger ΔG*. Our analysis of near-critical translationally ordered nuclei of octahedra shows a strong bias towards an orientational alignment which is incompatible with the tendency to form facet-to-facet contacts in the disordered phase, thus creating an additional entropic penalty for crystallization. For rotator phases of octahedra-like particles, we observe that the strength of the localized orientational order correlates inversely with ΔG*. We also observe that for s > 0.66 shapes and similar to octahedra, configurations with high facet alignment do not favor high orientational order, and thus ΔG*'s increase with truncation.

Equipartition terms in transition path ensemble: Insights from molecular dynamics simulations of alanine dipeptide

Transition path ensemble consists of reactive trajectories and possesses all the information necessary for the understanding of the mechanism and dynamics of important condensed phase processes. However, quantitative description of the properties of the transition path ensemble is far from being established. Here, with numerical calculations on a model system, the equipartition terms defined in thermal equilibrium were for the first time estimated in the transition path ensemble. It was not surprising to observe that the energy was not equally distributed among all the coordinates. However, the energies distributed on a pair of conjugated coordinates remained equal. Higher energies were observed to be distributed on several coordinates, which are highly coupled to the reaction coordinate, while the rest were almost equally distributed. In addition, the ensemble-averaged energy on each coordinate as a function of time was also quantified. These quantitative analyses on energy distributions provided new insights into the transition path ensemble.

Fast Decorrelating Monte Carlo Moves for Efficient Path Sampling

Many relevant processes in chemistry, physics, and biology are rare events from a computational perspective as they take place beyond the accessible time scale of molecular dynamics (MD). Examples are chemical reactions, nucleation, and conformational changes of biomolecules. Path sampling is an approach to break this time scale limit via a Monte Carlo (MC) sampling of MD trajectories. Still, many trajectories are needed for accurately predicting rate constants. To improve the speed of convergence, we propose two new MC moves, stone skipping and web throwing. In these moves, trajectories are constructed via a sequence of subpaths obeying superdetailed balance. By a reweighting procedure, almost all paths can be accepted. Whereas the generation of a single trajectory becomes more expensive, the reduced correlation results in a significant speedup. For a study on DNA denaturation, the increase was found to be a factor 12.

Structurally Linked Dynamics in Lactate Dehydrogenases of Evolutionarily Distinct Species

We present new findings about how primary and secondary structure affects the role of fast protein motions in the reaction coordinates of enzymatic reactions. Using transition path sampling and committor distribution analysis, we examined the difference in the role of these fast protein motions in the reaction coordinate of lactate dehydrogenases (LDHs) of Apicomplexa organisms Plasmodium falciparum and Cryptosporidium parvum. Having evolved separately from a common malate dehydrogenase ancestor, the two enzymes exhibit several important structural differences, notably a five-amino acid insertion in the active site loop of P. falciparum LDH. We find that these active site differences between the two organisms' LDHs likely cause a decrease in the contribution of the previously determined LDH rate-promoting vibration to the reaction coordinate of P. falciparum LDH compared to that of C. parvum LDH, specifically in the coupling of the rate-promoting vibration and the hydride transfer. This effect, while subtle, directly shows how changes in structure near the active site of LDH alter catalytically important motions. Insights provided by studying these alterations would prove to be useful in identifying LDH inhibitors that specifically target the isozymes of these parasitic organisms.

Jump Markov models and transition state theory: the quasi-stationary distribution approach

We are interested in the connection between a metastable continuous state space Markov process (satisfying e.g. the Langevin or overdamped Langevin equation) and a jump Markov process in a discrete state space. More precisely, we use the notion of quasi-stationary distribution within a metastable state for the continuous state space Markov process to parametrize the exit event from the state. This approach is useful to analyze and justify methods which use the jump Markov process underlying a metastable dynamics as a support to efficiently sample the state-to-state dynamics (accelerated dynamics techniques). Moreover, it is possible by this approach to quantify the error on the exit event when the parametrization of the jump Markov model is based on the Eyring-Kramers formula. This therefore provides a mathematical framework to justify the use of transition state theory and the Eyring-Kramers formula to build kinetic Monte Carlo or Markov state models.

A benchmark for reaction coordinates in the transition path ensemble

The molecular mechanism of a reaction is embedded in its transition path ensemble, the complete collection of reactive trajectories. Utilizing the information in the transition path ensemble alone, we developed a novel metric, which we termed the emergent potential energy, for distinguishing reaction coordinates from the bath modes. The emergent potential energy can be understood as the average energy cost for making a displacement of a coordinate in the transition path ensemble. Where displacing a bath mode invokes essentially no cost, it costs significantly to move the reaction coordinate. Based on some general assumptions of the behaviors of reaction and bath coordinates in the transition path ensemble, we proved theoretically with statistical mechanics that the emergent potential energy could serve as a benchmark of reaction coordinates and demonstrated its effectiveness by applying it to a prototypical system of biomolecular dynamics. Using the emergent potential energy as guidance, we developed a committor-free and intuition-independent method for identifying reaction coordinates in complex systems. We expect this method to be applicable to a wide range of reaction processes in complex biomolecular systems.

Primary Fibril Nucleation of Aggregation Prone Tau Fragments PHF6 and PHF6

We performed replica exchange molecular dynamics and forward flux sampling simulations of hexapeptide VQIINK and VQIVYK systems, also known as, respectively, fragments PHF6* and PHF6 from the tau protein. Being a part of the microtubule binding region, these fragments are known to be aggregation prone, and at least one of them is a prerequisite for fibril formation of the tau protein. Using a coarse-grained force field, we establish the phase behavior of both fragments, and investigate the nucleation kinetics for the conversion into a β-sheet fibril. As the conversion is, in principle, a reversible process, we predict the rate constants for both the fibril formation and melting, and examine the corresponding mechanisms. Our simulations indicate that, while both fragments form disordered aggregates, only PHF6 is able to form β-sheet fibrils. This observation provides a possible explanation for the lack of available steric zipper crystal structures for PHF6*.

Pathways to self-organization: Crystallization via nucleation and growth

Crystallization, a prototypical self-organization process during which a disordered state spontaneously transforms into a crystal characterized by a regular arrangement of its building blocks, usually proceeds by nucleation and growth. In the initial stages of the transformation, a localized nucleus of the new phase forms in the old one due to a random fluctuation. Most of these nuclei disappear after a short time, but rarely a crystalline embryo may reach a critical size after which further growth becomes thermodynamically favorable and the entire system is converted into the new phase. In this article, we will discuss several theoretical concepts and computational methods to study crystallization. More specifically, we will address the rare event problem arising in the simulation of nucleation processes and explain how to calculate nucleation rates accurately. Particular attention is directed towards discussing statistical tools to analyze crystallization trajectories and identify the transition mechanism.

Path Sampling Methods for Enzymatic Quantum Particle Transfer Reactions

The mechanisms of enzymatic reactions are studied via a host of computational techniques. While previous methods have been used successfully, many fail to incorporate the full dynamical properties of enzymatic systems. This can lead to misleading results in cases where enzyme motion plays a significant role in the reaction coordinate, which is especially relevant in particle transfer reactions where nuclear tunneling may occur. In this chapter, we outline previous methods, as well as discuss newly developed dynamical methods to interrogate mechanisms of enzymatic particle transfer reactions. These new methods allow for the calculation of free energy barriers and kinetic isotope effects (KIEs) with the incorporation of quantum effects through centroid molecular dynamics (CMD) and the full complement of enzyme dynamics through transition path sampling (TPS). Recent work, summarized in this chapter, applied the method for calculation of free energy barriers to reaction in lactate dehydrogenase (LDH) and yeast alcohol dehydrogenase (YADH). We found that tunneling plays an insignificant role in YADH but plays a more significant role in LDH, though not dominant over classical transfer. Additionally, we summarize the application of a TPS algorithm for the calculation of reaction rates in tandem with CMD to calculate the primary H/D KIE of YADH from first principles. We found that the computationally obtained KIE is within the margin of error of experimentally determined KIEs and corresponds to the KIE of particle transfer in the enzyme. These methods provide new ways to investigate enzyme mechanism with the inclusion of protein and quantum dynamics.