p(TP|q) peak maximization: Necessary but not sufficient for reaction coordinate accuracy

Optimal Reaction Coordinates and Kinetic Rates from the Projected Dynamics of Transition Paths

Finding optimal reaction coordinates and predicting accurate kinetic rates for activated processes are two of the foremost challenges of molecular simulations. We introduce an algorithm that tackles the two problems at once: starting from a limited number of reactive molecular dynamics trajectories (transition paths), we automatically generate with a Monte Carlo approach a sequence of different reaction coordinates that progressively reduce the kinetic rate of their projected effective dynamics. Based on a variational principle, the minimal rate accurately approximates the exact one, and it corresponds to the optimal reaction coordinate. After benchmarking the method on an analytic double-well system, we apply it to complex atomistic systems: the interaction of carbon nanoparticles of different sizes in water.

Explaining reaction coordinates of alanine dipeptide isomerization obtained from deep neural networks using Explainable Artificial Intelligence (XAI)

A method for obtaining appropriate reaction coordinates is required to identify transition states distinguishing product and reactant in complex molecular systems. Recently, abundant research has been devoted to obtaining reaction coordinates using artificial neural networks from deep learning literature, where many collective variables are typically utilized in the input layer. However, it is difficult to explain the details of which collective variables contribute to the predicted reaction coordinates owing to the complexity of the nonlinear functions in deep neural networks. To overcome this limitation, we used Explainable Artificial Intelligence (XAI) methods of the Local Interpretable Model-agnostic Explanation (LIME) and the game theory-based framework known as Shapley Additive exPlanations (SHAP). We demonstrated that XAI enables us to obtain the degree of contribution of each collective variable to reaction coordinates that is determined by nonlinear regressions with deep learning for the committor of the alanine dipeptide isomerization in vacuum. In particular, both LIME and SHAP provide important features to the predicted reaction coordinates, which are characterized by appropriate dihedral angles consistent with those previously reported from the committor test analysis. The present study offers an AI-aided framework to explain the appropriate reaction coordinates, which acquires considerable significance when the number of degrees of freedom increases.

Two distinct catalytic pathways for GH43 xylanolytic enzymes unveiled by X-ray and QM/MM simulations

Xylanolytic enzymes from glycoside hydrolase family 43 (GH43) are involved in the breakdown of hemicellulose, the second most abundant carbohydrate in plants. Here, we kinetically and mechanistically describe the non-reducing-end xylose-releasing exo-oligoxylanase activity and report the crystal structure of a native GH43 Michaelis complex with its substrate prior to hydrolysis. Two distinct calcium-stabilized conformations of the active site xylosyl unit are found, suggesting two alternative catalytic routes. These results are confirmed by QM/MM simulations that unveil the complete hydrolysis mechanism and identify two possible reaction pathways, involving different transition state conformations for the cleavage of xylooligosaccharides. Such catalytic conformational promiscuity in glycosidases is related to the open architecture of the active site and thus might be extended to other exo-acting enzymes. These findings expand the current general model of catalytic mechanism of glycosidases, a main reaction in nature, and impact on our understanding about their interaction with substrates and inhibitors.

Identification of simple reaction coordinates from complex dynamics

Reaction coordinates are widely used throughout chemical physics to model and understand complex chemical transformations. We introduce a definition of the natural reaction coordinate, suitable for condensed phase and biomolecular systems, as a maximally predictive one-dimensional projection. We then show that this criterion is uniquely satisfied by a dominant eigenfunction of an integral operator associated with the ensemble dynamics. We present a new sparse estimator for these eigenfunctions which can search through a large candidate pool of structural order parameters and build simple, interpretable approximations that employ only a small number of these order parameters. Example applications with a small molecule's rotational dynamics and simulations of protein conformational change and folding show that this approach can filter through statistical noise to identify simple reaction coordinates from complex dynamics.

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.

Reducing the cost of evaluating the committor by a fitting procedure

Correct identification of reaction coordinates in complex systems is essential for understanding the mechanisms of their reaction dynamics. Existing methods for identifying reaction coordinates typically require knowledge of the committor--the probability of a given configuration to reach the product basin. The high computational cost of evaluating committors has limited applications of methods for identifying reaction coordinates. We proposed a fitting procedure that can reduce the cost of evaluating committors by an order of magnitude or more. The method only requires evaluating the committors of a few configurations in a transition path by the standard and costly shooting procedure. The committors of the other configurations are then estimated with great accuracy by a sigmoid function derived from fitting the few numerically evaluated committors. The method has been systematically tested on a model system of a Brownian particle moving in a one-dimensional double-well potential, and a small biomolecular system--the isomerization of alanine dipeptide in vacuum and in explicit water.

Recent developments in methods for identifying reaction coordinates

In the study of rare events in complex systems with many degrees of freedom, a key element is to identify the reaction coordinates of a given process. Over recent years, a number of methods and protocols have been developed to extract the reaction coordinates based on limited information from molecular dynamics simulations. In this review, we provide a brief survey over a number of major methods developed in the past decade, some of which are discussed in greater detail, to provide an overview of the problems that are partially solved and challenges that still remain. A particular emphasis has been placed on methods for identifying reaction coordinates that are related to the committor.

Reaction coordinates, one-dimensional Smoluchowski equations, and a test for dynamical self-consistency

We propose a method for identifying accurate reaction coordinates among a set of trial coordinates. The method applies to special cases where motion along the reaction coordinate follows a one-dimensional Smoluchowski equation. In these cases the reaction coordinate can predict its own short-time dynamical evolution, i.e., the dynamics projected from multiple dimensions onto the reaction coordinate depend only on the reaction coordinate itself. To test whether this property holds, we project an ensemble of short trajectory swarms onto trial coordinates and compare projections of individual swarms to projections of the ensemble of swarms. The comparison, quantified by the Kullback-Leibler divergence, is numerically performed for each isosurface of each trial coordinate. The ensemble of short dynamical trajectories is generated only once by sampling along an initial order parameter. The initial order parameter should separate the reactants and products with a free energy barrier, and distributions on isosurfaces of the initial parameter should be unimodal. The method is illustrated for three model free energy landscapes with anisotropic diffusion. Where exact coordinates can be obtained from Kramers-Langer-Berezhkovskii-Szabo theory, results from the new method agree with the exact results. We also examine characteristics of systems where the proposed method fails. We show how dynamical self-consistency is related (through the Chapman-Kolmogorov equation) to the earlier isocommittor criterion, which is based on longer paths.

Friction contribution to water-bond breakage kinetics

Based on the trajectories of the separation between water molecule pairs from MD simulations, we investigate the bond breakage dynamics in bulk water. From the spectrum of mean first-passage times, the Fokker-Planck equation allows us to derive the diffusivity profile along the separation coordinate and thus to unambiguously disentangle the effects of free-energy and local friction on the separation kinetics. For tightly coordinated water, the friction is six times higher than in bulk, which can be interpreted in terms of a dominant reaction path that involves additional orthogonal coordinates.

Splitting probabilities as a test of reaction coordinate choice in single-molecule experiments

To explain the observed dynamics in equilibrium single-molecule measurements of biomolecules, the experimental observable is often chosen as a putative reaction coordinate along which kinetic behavior is presumed to be governed by diffusive dynamics. Here, we invoke the splitting probability as a test of the suitability of such a proposed reaction coordinate. Comparison of the observed splitting probability with that computed from the kinetic model provides a simple test to reject poor reaction coordinates. We demonstrate this test for a force spectroscopy measurement of a DNA hairpin.

Optimizing Nucleus Size Metrics for Liquid-Solid Nucleation from Transition Paths of Near-Nanosecond Duration

We determine the mechanism for the liquid-solid phase transition in the Lennard-Jones fluid close to coexistence with aimless shooting and likelihood maximization. The reaction coordinate for this process is a product of a structural descriptor and the size of the nascent solid nucleus and is quantitatively verified with the committor probability histogram test. This study identifies the first accurate scalar reaction coordinate for the liquid-solid nucleation process in Lennard-Jonesium, which will likely extend to nucleation processes in other spherically symmetric fluids. On the basis of our results, we propose a structural correction factor for the commonly cited nucleus size reaction coordinate from classical nucleation theory that enables connection of simulation data to stochastic models of nucleation kinetics. In addition, we show that aimless shooting is able to obtain reasonable acceptance rates for transitions with highly diffusive characteristics, which has been problematic for transition path sampling methods for diffusive processes such as nucleation and macromolecular transitions.