A benchmark for reaction coordinates in the transition path ensemble

Reactive Vortexes in a Naturally Activated Process: Non-Diffusive Rotational Fluxes at Transition State Uncovered by Persistent Homology

The dynamics of reaction coordinates during barrier-crossing are key to understanding activated processes in complex systems such as proteins. The default assumption from Kramers' physical intuition is that of a diffusion process. However, the dynamics of barrier-crossing in natural complex molecules are largely unexplored. Here we investigate the transition dynamics of alanine dipeptide isomerization, the simplest complex system with a large number of non-reaction coordinates that can serve as an adequate thermal bath feeding energy into the reaction coordinates. We separate conformations along the time axis and construct the dynamic probability surface of reaction. We quantify its topological structure and rotational flux using persistent homology and differential form. Our results uncovered a region with a strong reactive vortex in the configuration-time space, where the highest probability peak and the transition state ensemble are located. This reactive region contains strong rotational fluxes: Most reactive trajectories swirl multiple times around this region in the subspace of the two most important reaction coordinates. Furthermore, the rotational fluxes result from cooperative movement along the isocommitter surfaces and orthogonal barrier-crossing. Overall, our findings offer a first glimpse into the reactive vortex regions that characterize the non-diffusive dynamics of barrier-crossing of a naturally occurring activation process.

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.

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.

Nonparametric Analysis of Nonequilibrium Simulations

We extend the nonparametric framework of reaction coordinate optimization to nonequilibrium ensembles of (short) trajectories. For example, we show how, starting from such an ensemble, one can obtain an equilibrium free-energy profile along the committor, which can be used to determine important properties of the dynamics exactly. A new adaptive sampling approach, the transition-state ensemble enrichment, is suggested, which samples the configuration space by "growing" committor segments toward each other starting from the boundary states. This framework is suggested as a general tool, alternative to the Markov state models, for a rigorous and accurate analysis of simulations of large biomolecular systems, as it has the following attractive properties. It is immune to the curse of dimensionality, does not require system-specific information, can approximate arbitrary reaction coordinates with high accuracy, and has sensitive and rigorous criteria to test optimality and convergence. The approaches are illustrated on a 50-dimensional model system and a realistic protein folding trajectory.

Exact Topology of the Dynamic Probability Surface of an Activated Process by Persistent Homology

To gain insight into the reaction mechanism of activated processes, we introduce an exact approach for quantifying the topology of high-dimensional probability surfaces of the underlying dynamic processes. Instead of Morse indexes, we study the homology groups of a sequence of superlevel sets of the probability surface over high-dimensional configuration spaces using persistent homology. For alanine-dipeptide isomerization, a prototype of activated processes, we identify locations of probability peaks and connecting ridges, along with measures of their global prominence. Instead of a saddle point, the transition state ensemble (TSE) of conformations is at the most prominent probability peak after reactants/products, when proper reaction coordinates are included. Intuition-based models, even those exhibiting a double-well, fail to capture the dynamics of the activated process. Peak occurrence, prominence, and locations can be distorted upon subspace projection. While principal component analysis accounts for conformational variance, it inflates the complexity of the surface topology and destroys the dynamic properties of the topological features. In contrast, TSE emerges naturally as the most prominent peak beyond the reactant/product basins, when projected to a subspace of minimum dimension containing the reaction coordinates. Our approach is general and can be applied to investigate the topology of high-dimensional probability surfaces of other activated processes.

Dissecting the Dynamics during Enzyme Catalysis: A Case Study of Pin1 Peptidyl-Prolyl Isomerase

Free energy surfaces have played a central role in studying protein conformational changes and enzymatic reactions over decades. Yet, free energy barriers and kinetics are highly dependent on the coordinates chosen to define the surface, and furthermore, the dynamics during the reactions are often overlooked. Our recent study on the Pin1-catalyzed isomerization reaction has indicated that the isomerization transition events remarkably deviate from the free energy path, highlighting the need to understand the reaction dynamics in more detail. To this end, here we investigate the reaction coordinates that describe the transition states of the free energy and transition pathways by minimizing the cross-entropy function. We show that the isomerization transition events can be expressed by the concerted changes in the improper torsion angle ζ and nearby backbone torsional angles of the ligand, whereas the transition state of the free energy surface involves changes in a broad range of coordinates including multiple protein-ligand interactions. The current result supports the previous finding that the isomerization transitions occur quickly from the conformational excited states, which is in sharp contrast to the slow and collective changes suggested from the free energy path. Our results further indicate that the coordinates derived from the transition trajectories are not sufficient for finding the transition states on the free energy surfaces due to the lack of information from conformational excited states.

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.

Orthogonal order parameters to model the reaction coordinate of an enzyme catalyzed reaction

The choice of suitable collective variables in formulating an optimal reaction coordinate is a challenging task for activated transitions between a pair of stable states especially when dealing with biochemical changes such as enzyme catalyzed reactions. A detailed benchmarking study is carried out on the choice of collective variables that can distinguish between the stable states unambiguously. We specifically address the issue if these variables may be directly used to model the optimal reaction coordinate, or if it would be better to use their orthogonalized counterparts. The proposed computational scheme is applied to the rate determining intramolecular proton transfer step in the enzyme human carbonic anhydrase II. The optimum reaction coordinate is determined with and without orthogonalization of the collective variables pertinent to a key conformational fluctuation and the actual proton transfer step at the active site of the enzyme. Suitability of the predicted reaction coordinates in different processes is examined in terms of the free energy profile projected along the reaction coordinate, the rate constant of transition and the underlying molecular mechanism of barrier crossing. Our results indicate that a better agreement with earlier simulation and experimental data is obtained when the orthogonalized collective variables are used to model the reaction coordinate.

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.