Transition Path Ensemble Research Articles

Splitting probabilities as optimal controllers of rare reactive events.

The committor constitutes the primary quantity of interest within chemical kinetics as it is understood to encode the ideal reaction coordinate for a rare reactive event. We show the generative utility of the committor in that it can be used explicitly to produce a reactive trajectory ensemble that exhibits numerically exact statistics as that of the original transition path ensemble. This is done by relating a time-dependent analog of the committor that solves a generalized bridge problem to the splitting probability that solves a boundary value problem under a bistable assumption. By invoking stochastic optimal control and spectral theory, we derive a general form for the optimal controller of a bridge process that connects two metastable states expressed in terms of the splitting probability. This formalism offers an alternative perspective into the role of the committor and its gradients in that they encode force fields that guarantee reactivity, generating trajectories that are statistically identical to the way that a system would react autonomously.

Probing reaction channels via reinforcement learning

Chemical reactions are dynamical processes involving the correlated reorganization of atomic configurations, driving the conversion of an initial reactant into a result product. By virtue of the metastability of both the reactants and products, chemical reactions are rare events, proceeding fleetingly. Reaction pathways can be modelled probabilistically by using the notion of reactive density in the phase space of the molecular system. Such density is related to a function known as the committor function, which describes the likelihood of a configuration evolving to one of the nearby metastable regions. In theory, the committor function can be obtained by solving the backward Kolmogorov equation (BKE), which is a partial differential equation (PDE) defined in the full dimensional phase space. However, using traditional methods to solve this problem is not practical for high dimensional systems. In this work, we propose a reinforcement learning based method to identify important configurations that connect reactant and product states along chemical reaction paths. By shooting multiple trajectories from these configurations, we can generate an ensemble of states that concentrate on the transition path ensemble. This configuration ensemble can be effectively employed in a neural network-based PDE solver to obtain an approximation solution of a restricted BKE, even when the dimension of the problem is very high. The resulting solution provides an approximation for the committor function that encodes mechanistic information for the reaction, paving a new way for understanding of complex chemical reactions and evaluation of reaction rates.

Variational deep learning of equilibrium transition path ensembles.

We present a time-dependent variational method to learn the mechanisms of equilibrium reactive processes and efficiently evaluate their rates within a transition path ensemble. This approach builds off of the variational path sampling methodology by approximating the time-dependent commitment probability within a neural network ansatz. The reaction mechanisms inferred through this approach are elucidated by a novel decomposition of the rate in terms of the components of a stochastic path action conditioned on a transition. This decomposition affords an ability to resolve the typical contribution of each reactive mode and their couplings to the rare event. The associated rate evaluation is variational and systematically improvable through the development of a cumulant expansion. We demonstrate this method in both over- and under-damped stochastic equations of motion, in low-dimensional model systems, and in the isomerization of a solvated alanine dipeptide. In all examples, we find that we can obtain quantitatively accurate estimates of the rates of the reactive events with minimal trajectory statistics and gain unique insights into transitions through the analysis of their commitment probability.

Deep learning collective variables from transition path ensemble.

The study of the rare transitions that take place between long lived metastable states is a major challenge in molecular dynamics simulations. Many of the methods suggested to address this problem rely on the identification of the slow modes of the system, which are referred to as collective variables. Recently, machine learning methods have been used to learn the collective variables as functions of a large number of physical descriptors. Among many such methods, Deep Targeted Discriminant Analysis has proven to be useful. This collective variable is built from data harvested from short unbiased simulations in the metastable basins. Here, we enrich the set of data on which the Deep Targeted Discriminant Analysis collective variable is built by adding data from the transition path ensemble. These are collected from a number of reactive trajectories obtained using the On-the-fly Probability Enhanced Sampling flooding method. The collective variables thus trained lead to more accurate sampling and faster convergence. The performance of these new collective variables is tested on a number of representative examples.

Analysis of transition path ensemble in the exactly solvable models via overdamped langevin equation

Transition of a system between two states is an important but difficult problem in natural science. In this article we study the transition problem in the framework of transition path ensemble. Using the overdamped Langevin method, we introduce the path integral formulation of the transition probability and obtain the equation for the minimum action path in the transition path space. For the effective sampling in the transition path ensemble, we derive a conditional overdamped Langevin equation. In two exactly solvable models, the free particle system and the harmonic system, we present the expression of the conditional probability density and the explicit solutions for the conditional Langevin equation and the minimum action path. The analytic results demonstrate the consistence of the conditional Langevin equation with the desired probability distribution in the transition. It is confirmed that the conditional Langevin equation is an effective tool to sample the transition path ensemble, and the stationary action principle actually leads to the most probable path.

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.

Energy Decomposition along 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.

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.

Measurement of irreversibility and entropy production via the tubular ensemble

The appealing theoretical measure of irreversibility in a stochastic process, as the ratio of the probabilities of a trajectory and its time reversal, cannot be accessed directly in experiments since the probability of a single trajectory is zero. We regularize this definition by considering, instead, the limiting ratio of probabilities for trajectories to remain in the tubular neighborhood of a smooth path and its time reversal. The resulting pathwise medium entropy production agrees with the formal expression from stochastic thermodynamics and can be obtained from measurable tube probabilities. Estimating the latter from numerically sampled trajectories for Langevin dynamics yields excellent agreement with theory. By combining our measurement of pathwise entropy production with a Markov chain Monte Carlo algorithm, we infer the entropy-production distribution for a transition path ensemble directly from short recorded trajectories. Our work enables the measurement of irreversibility along individual paths and path ensembles in a model-free manner.

Optimizing reaction coordinate by flux maximization in the transition path ensemble.

The transition path ensemble is a collection of reactive trajectories, all of which largely keep going forward along the transition channel from the reactant state to the product one, and is believed to possess the information necessary for the identification of reaction coordinates. Previously, the full coordinates (both position and momentum) of the snapshots in the transition path ensemble were utilized to obtain the reaction coordinate [W. Li and A. Ma, J. Chem. Phys. 144, 114103 (2016) and W. Li, J. Chem. Phys. 148, 084105 (2018)]. Here, with the conformational (or position) coordinates alone, it is suggested that a meaningful one-dimensional representation of the reaction coordinate can be obtained by maximizing the flux (the net number of transitions through the dividing surface) of a given coordinate in the transition path ensemble. In the application to alanine dipeptide in vacuum, dihedral angles ϕ and θ were identified to be the two best reaction coordinates, which was consistent with the results in the existing studies. A linear combination of these two coordinates gave a better reaction coordinate, which is highly correlated with the committor. Most importantly, the method obtained a linear combination of pairwise distances between heavy atoms, which was highly correlated with the committor as well. The standard deviation of the committor at the transition region defined by the optimized reaction coordinate is as small as 0.08. In addition, the effects of practical factors, such as the choice of transition path sub-ensembles and saving interval between frames in transition paths, on reaction coordinate optimization were also considered.

Efficient and exact sampling of transition path ensembles on Markovian networks.

The problem of flickering trajectories in standard kinetic Monte Carlo (kMC) simulations prohibits sampling of the transition path ensembles (TPEs) on Markovian networks representing many slow dynamical processes of interest. In the present contribution, we overcome this problem using knowledge of the metastable macrostates, determined by an unsupervised community detection algorithm, to perform enhanced sampling kMC simulations. We implement two accelerated kMC methods to simulate the nonequilibrium stochastic dynamics on arbitrary Markovian networks, namely, weighted ensemble (WE) sampling and kinetic path sampling (kPS). WE-kMC utilizes resampling in pathway space to maintain an ensemble of representative trajectories covering the state space, and kPS utilizes graph transformation to simplify the description of an escape trajectory from a trapping energy basin. Both methods sample individual trajectories governed by the linear master equation with the correct statistical frequency. We demonstrate that they allow for efficient estimation of the time-dependent occupation probability distributions for the metastable macrostates, and of TPE statistics, such as committor functions and first passage time distributions. kPS is particularly attractive, since its efficiency is essentially independent of the degree of metastability, and we suggest how the algorithm could be coupled with other enhanced sampling methodologies. We illustrate our approach with results for a network representing the folding transition of a tryptophan zipper peptide, which exhibits a separation of characteristic timescales. We highlight some salient features of the dynamics, most notably, strong deviations from two-state behavior, and the existence of multiple competing mechanisms.

On Calculating Free Energy Differences Using Ensembles of Transition Paths.

The free energy of a process is the fundamental quantity that determines its spontaneity or propensity at a given temperature. In particular, the binding free energy of a drug candidate to its biomolecular target is used as an objective quantity in drug design. Recently, binding kinetics—rates of association (kon) and dissociation (koff)—have also demonstrated utility for their ability to predict efficacy and in some cases have been shown to be more predictive than the binding free energy alone. Some methods exist to calculate binding kinetics from molecular simulations, although these are typically more difficult to calculate than the binding affinity as they depend on details of the transition path ensemble. Assessing these rate constants can be difficult, due to uncertainty in the definition of the bound and unbound states, large error bars and the lack of experimental data. As an additional consistency check, rate constants from simulation can be used to calculate free energies (using the log of their ratio) which can then be compared to free energies obtained experimentally or using alchemical free energy perturbation. However, in this calculation it is not straightforward to account for common, practical details such as the finite simulation volume or the particular definition of the “bound” and “unbound” states. Here we derive a set of correction terms that can be applied to calculations of binding free energies using full reactive trajectories. We apply these correction terms to revisit the calculation of binding free energies from rate constants for a host-guest system that was part of a blind prediction challenge, where significant deviations were observed between free energies calculated with rate ratios and those calculated from alchemical perturbation. The correction terms combine to significantly decrease the error with respect to computational benchmarks, from 3.4 to 0.76 kcal/mol. Although these terms were derived with weighted ensemble simulations in mind, some of the correction terms are generally applicable to free energies calculated using physical pathways via methods such as Markov state modeling, metadynamics, milestoning, or umbrella sampling.

Residue-Residue Mutual Work Analysis of Retinal-Opsin Interaction in Rhodopsin: Implications for Protein-Ligand Binding.

Energetic contributions at the single-residue level for retinal-opsin interactions in rhodopsin were studied by combining molecular dynamics simulations, transition path sampling, and a newly developed energy decomposition approach. The virtual work at an infinitesimal time interval was decomposed into the work components on one residue due to its interaction with another residue, which were then averaged over the transition path ensemble along a proposed reaction coordinate. Such residue-residue mutual work analysis on 62 residues within the active center of rhodopsin resulted in a very sparse interaction matrix, which is generally not symmetric but antisymmetric to some extent. Fourteen residues were identified to be major players in retinal relaxation along a plausible pathway from bathorhodopsin to the blue-shifted intermediate, which is in good agreement with an existing NMR study. Based on the matrix of mutual work, a comprehensive network was constructed to provide detailed insights into the chromophore-protein interaction from a viewpoint of energy flow.

Interlayer hopping dynamics of bilayer water confined between graphene sheets

Abstract The ubiquitous existence of water confined in nano-capillaries is key to fundamental biological and technological applications. Using molecular dynamics simulations, we analyzed the hopping-like interlayer relocation dynamics of bilayer water confined between two parallel graphene sheets. In contrary to the common scheme that relocation is driven by density fluctuations, analysis of the transition path ensemble revealed that interlayer hopping is induced by local hydrogen bond configuration fluctuations coupled with activated consecutive transient angular reorientations. Our results consolidated the anisotropic nature of water relocation under strong ordering, which provided a mechanistic insight into the relaxation behavior at glass-forming water’s fragile-to-strong crossover.

Transition path theory from biased simulations.

Transition Path Theory (TPT) provides a rigorous framework to investigate the dynamics of rare thermally activated transitions. In this theory, a central role is played by the forward committor function q+(x), which provides the ideal reaction coordinate. Furthermore, the reactive dynamics and kinetics are fully characterized in terms of two time-independent scalar and vector distributions. In this work, we develop a scheme which enables all these ingredients of TPT to be efficiently computed using the short non-equilibrium trajectories generated by means of a specific combination of enhanced path sampling techniques. In particular, first we further extend the recently introduced self-consistent path sampling algorithm in order to compute the committor q+(x). Next, we show how this result can be exploited in order to define efficient algorithms which enable us to directly sample the transition path ensemble.

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.

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.

Reconciling transition path time and rate measurements in reactions with large entropic barriers.

Recent experiments and simulation studies showed that protein/DNA folding barriers inferred from folding rates or from potentials of mean force are often much higher than the barriers estimated from the distributions of transition path times. Here a toy model is used to explain a possible origin of this effect: It is shown that when the transition in question involves an entropic barrier, the one-dimensional Langevin model commonly used to interpret experimental data, while adequately predicting the transition rate, fails to describe the properties of the subset of the trajectories that form the transition path ensemble; the latter may still be describable in terms of a one-dimensional model, but with a different potential, just as observed experimentally.

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.

Reaction mechanism and reaction coordinates from the viewpoint of energy flow

Reaction coordinates are of central importance for correct understanding of reaction dynamics in complex systems, but their counter-intuitive nature made it a daunting challenge to identify them. Starting from an energetic view of a reaction process as stochastic energy flows biased towards preferred channels, which we deemed the reaction coordinates, we developed a rigorous scheme for decomposing energy changes of a system, both potential and kinetic, into pairwise components. The pairwise energy flows between different coordinates provide a concrete statistical mechanical language for depicting reaction mechanisms. Application of this scheme to the C7eq → C7ax transition of the alanine dipeptide in vacuum revealed novel and intriguing mechanisms that eluded previous investigations of this well studied prototype system for biomolecular conformational dynamics. Using a cost function developed from the energy decomposition components by proper averaging over the transition path ensemble, we were able to identify signatures of the reaction coordinates of this system without requiring any input from human intuition.