Transition rate theory, spectral analysis, and reactive paths

Reaction Coordinates Are Optimal Channels of Energy Flow

Reaction coordinates (RCs) are the few essential coordinates of a protein that control its functional processes, such as allostery, enzymatic reaction, and conformational change. They are critical for understanding protein function and provide optimal enhanced sampling of protein conformational changes and states. Since the pioneering work in the late 1990s, identifying the correct and objectively provable RCs has been a central topic in molecular biophysics and chemical physics. This review summarizes the major advances in identifying RCs over the past 25 years, focusing on methods aimed at finding RCs that meet the rigorous committor criterion, widely accepted as the true RCs. Notably, the newly developed physics-based energy flow theory and generalized work functional method provide a general and rigorous approach for identifying true RCs, revealing their physical nature as the optimal channels of energy flow in biomolecules.

Slow dynamical modes from static averages

In recent times, efforts have been made to describe the evolution of a complex system not through long trajectories but via the study of probability distribution evolution. This more collective approach can be made practical using the transfer operator formalism and its associated dynamics generator. Here, we reformulate in a more transparent way the result of Devergne et al. [Adv. Neural Inform. Process. Syst. 37, 75495-75521 (2024)] and show that the lowest eigenfunctions and eigenvalues of the dynamics generator can be efficiently computed using data easily obtainable from biased simulations. We also show explicitly that the long time dynamics can be reconstructed by using the spectral decomposition of the dynamics operator.

Mechanism of Fe(II) Chemisorption on Hematite(001) Revealed by Reactive Neural Network Potential Molecular Dynamics

Atomic-scale understanding of important geochemical processes including sorption, dissolution, nucleation, and crystal growth is difficult to obtain from experimental measurements alone and would benefit from strong continuous progress in molecular simulation. To this end, we present a reactive neural network potential-based molecular dynamics approach to simulate the interaction of aqueous ions on mineral surfaces in contact with liquid water, taking Fe(II) on hematite(001) as a model system. We show that a single neural network potential predicts rate constants for water exchange for aqueous Fe(II) and for the exergonic chemisorption of aqueous Fe(II) on hematite(001) in good agreement with experimental observations. The neural network potential developed herein allows one to converge free energy profiles and transmission coefficients at density functional theory-level accuracy outperforming state-of-the-art classical force field potentials. This suggests that machine learning potential molecular dynamics should become the method of choice for atomistic studies of geochemical processes.

Statistical Mechanical Theories of Membrane Permeability

A popular theoretical framework to compute the permeability coefficient of a molecule is provided by the classic Smoluchowski-Kramers treatment of the steady-state diffusive flux across a free-energy barrier. Within this framework, commonly termed "inhomogeneous solubility-diffusion" (ISD), the permeability, P, is expressed in closed form in terms of the potential of mean force and position-dependent diffusivity of the molecule of interest along the membrane normal. In principle, both quantities can be calculated from all-atom MD simulations. Although several methods exist for calculating the position-dependent diffusivity, each of these is at best an estimate. In addition, the ISD model does not account for memory effects along the chosen reaction coordinate. For these reasons, it is important to seek alternative theoretical formulations to determine the permeability coefficient that are able to account for the factors ignored by the ISD approximation. Using Green-Kubo linear response theory, we establish the familiar constitutive relation between the flux density across the membrane and the difference in the concentration of a permeant molecule, j = PΔC. On this basis, we derive a time-correlation function expression for the nonequilibrium flux across a membrane that is reminiscent of the transmission coefficient in the reactive flux formalism treatment of transition rates. An analysis based on the transition path theory framework is exploited to derive alternative expressions for the permeability coefficient. The different strategies are illustrated with stochastic simulations based on the generalized Langevin equation in addition to unbiased molecular dynamics simulations of water permeation of a lipid bilayer.

Spectral Map for Slow Collective Variables, Markovian Dynamics, and Transition State Ensembles

Understanding the behavior of complex molecular systems is a fundamental problem in physical chemistry. To describe the long-time dynamics of such systems, which is responsible for their most informative characteristics, we can identify a few slow collective variables (CVs) while treating the remaining fast variables as thermal noise. This enables us to simplify the dynamics and treat it as diffusion in a free-energy landscape spanned by slow CVs, effectively rendering the dynamics Markovian. Our recent statistical learning technique, spectral map [Rydzewski, J. J. Phys. Chem. Lett. 2023, 14(22), 5216-5220], explores this strategy to learn slow CVs by maximizing a spectral gap of a transition matrix. In this work, we introduce several advancements into our framework, using a high-dimensional reversible folding process of a protein as an example. We implement an algorithm for coarse-graining Markov transition matrices to partition the reduced space of slow CVs kinetically and use it to define a transition state ensemble. We show that slow CVs learned by spectral map closely approach the Markovian limit for an overdamped diffusion. We demonstrate that coordinate-dependent diffusion coefficients only slightly affect the constructed free-energy landscapes. Finally, we present how spectral maps can be used to quantify the importance of features and compare slow CVs with structural descriptors commonly used in protein folding. Overall, we demonstrate that a single slow CV learned by spectral map can be used as a physical reaction coordinate to capture essential characteristics of protein folding.

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.

Computing the committor with the committor to study the transition state ensemble

The study of the kinetic bottlenecks that hinder the rare transitions between long-lived metastable states is a major challenge in atomistic simulations. Here we propose a method to explore the transition state ensemble, which is the distribution of configurations that the system passes through as it translocates from one metastable basin to another. We base our method on the committor function and the variational principle that it obeys. We find its minimum through a self-consistent procedure that starts from information limited to the initial and final states. Right from the start, our procedure allows the sampling of very many transition state configurations. With the help of the variational principle, we perform a detailed analysis of the transition state ensemble, ranking quantitatively the degrees of freedom mostly involved in the transition and enabling a systematic approach for the interpretation of simulation results and the construction of efficient physics-informed collective variables.

Collective Variable-Based Enhanced Sampling: From Human Learning to Machine Learning

Enhanced-sampling algorithms relying on collective variables (CVs) are extensively employed to study complex (bio)chemical processes that are not amenable to brute-force molecular simulations. The selection of appropriate CVs characterizing the slow movement modes is of paramount importance for reliable and efficient enhanced-sampling simulations. In this Perspective, we first review the application and limitations of CVs obtained from chemical and geometrical intuition. We also introduce path-sampling algorithms, which can identify path-like CVs in a high-dimensional free-energy space. Machine-learning algorithms offer a viable approach to finding suitable CVs by analyzing trajectories from preliminary simulations. We discuss both the performance of machine-learning-derived CVs in enhanced-sampling simulations of experimental models and the challenges involved in applying these CVs to realistic, complex molecular assemblies. Moreover, we provide a prospective view of the potential advancements of machine-learning algorithms for the development of CVs in the field of enhanced-sampling simulations.

Dynamics of activation in the voltage-sensing domain of Ciona intestinalis phosphatase Ci-VSP

The Ciona intestinalis voltage-sensing phosphatase (Ci-VSP) is a membrane protein containing a voltage-sensing domain (VSD) that is homologous to VSDs from voltage-gated ion channels responsible for cellular excitability. Previously published crystal structures of Ci-VSD in putative resting and active conformations suggested a helical-screw voltage sensing mechanism in which the S4 helix translocates and rotates to enable exchange of salt-bridge partners, but the microscopic details of the transition between the resting and active conformations remained unknown. Here, by combining extensive molecular dynamics simulations with a recently developed computational framework based on dynamical operators, we elucidate the microscopic mechanism of the resting-active transition at physiological membrane potential. Sparse regression reveals a small set of coordinates that distinguish intermediates that are hidden from electrophysiological measurements. The intermediates arise from a noncanonical helical-screw mechanism in which translocation, rotation, and side-chain movement of the S4 helix are only loosely coupled. These results provide insights into existing experimental and computational findings on voltage sensing and suggest ways of further probing its mechanism.

How is Membrane Permeation of Small Ionizable Molecules Affected by Protonation Kinetics?

According to the pH-partition hypothesis, the aqueous solution adjacent to a membrane is a mixture of the ionization states of the permeating molecule at fixed Henderson-Hasselbalch concentrations, such that each state passes through the membrane in parallel with its own specific permeability. An alternative view, based on the assumption that the rate of switching ionization states is instantaneous, represents the permeation of ionizable molecules via an effective Boltzmann-weighted average potential (BWAP). Such an assumption is used in constant-pH molecular dynamics simulations. The inhomogeneous solubility-diffusion framework can be used to compute the pH-dependent membrane permeability for each of these two limiting treatments. With biased WTM-eABF molecular dynamics simulations, we computed the potential of mean force and diffusivity of each ionization state of two weakly basic small molecules: nicotine, an addictive drug, and varenicline, a therapeutic for treating nicotine addiction. At pH = 7, the BWAP effective permeability is greater than that determined by pH-partitioning by a factor of 2.5 for nicotine and 5 for varenicline. To assess the importance of ionization kinetics, we present a Smoluchowski master equation that includes explicitly the protonation and deprotonation processes coupled with the diffusive motion across the membrane. At pH = 7, the increase in permeability due to the explicit ionization kinetics is negligible for both nicotine and varenicline. This finding is reaffirmed by combined Brownian dynamics and Markov state model simulations for estimating the permeability of nicotine while allowing changes in its ionization state. We conclude that for these molecules the pH-partition hypothesis correctly captures the physics of the permeation process. The small free energy barriers for the permeation of nicotine and varenicline in their deprotonated neutral forms play a crucial role in establishing the validity of the pH-partitioning mechanism. Essentially, BWAP fails because ionization kinetics are too slow on the time scale of membrane crossing to affect the permeation of small ionizable molecules such as nicotine and varenicline. For the singly protonated state of nicotine, the computational results agree well with experimental measurements (P1 = 1.29 × 10-7 cm/s), but the agreement for neutral (P0 = 6.12 cm/s) and doubly protonated nicotine (P2 = 3.70 × 10-13 cm/s) is slightly worse, likely due to factors associated with the aqueous boundary layer (neutral form) or leaks through paracellular pathways (doubly protonated form).

Robust fragment-based method of calculating hydrogen atom transfer activation barrier in complex molecules

Dynamic processes driven by non-covalent interactions (NCI), such as conformational exchange, molecular binding, and solvation, can strongly influence the rate constants of reactions with low activation barriers, especially at low temperatures. Examples of this may include hydrogen-atom-transfer (HAT) reactions involved in the oxidative stress of an active pharmaceutical ingredient (API). Here, we develop an automated workflow to generate HAT transition-state (TS) geometries for complex and flexible APIs and then systematically evaluate the influences of NCI on the free activation energies, based on the multi-conformational transition-state theory (MC-TST) within the framework of a multi-step reaction path. The two APIs studied: fesoterodine and imipramine, display considerable conformational complexity and have multiple ways of forming hydrogen bonds with the abstracting radical-a hydroxymethyl peroxyl radical. Our results underscore the significance of considering conformational exchange and multiple activation pathways in activation calculations. We also show that structural elements and NCIs outside the reaction site minimally influence TS core geometry and covalent activation barrier, although they more strongly affect reactant binding and consequently the overall activation barrier. We further propose a robust and economical fragment-based method to obtain overall activation barriers, by combining the covalent activation barrier calculated for a small molecular fragment with the binding free energy calculated for the whole molecule.

Molecular Free Energies, Rates, and Mechanisms from Data-Efficient Path Sampling Simulations

Molecular dynamics is a powerful tool for studying the thermodynamics and kinetics of complex molecular events. However, these simulations can rarely sample the required time scales in practice. Transition path sampling overcomes this limitation by collecting unbiased trajectories and capturing the relevant events. Moreover, the integration of machine learning can boost the sampling while simultaneously learning a quantitative representation of the mechanism. Still, the resulting trajectories are by construction non-Boltzmann-distributed, preventing the calculation of free energies and rates. We developed an algorithm to approximate the equilibrium path ensemble from machine-learning-guided path sampling data. At the same time, our algorithm provides efficient sampling, mechanism, free energy, and rates of rare molecular events at a very moderate computational cost. We tested the method on the folding of the mini-protein chignolin. Our algorithm is straightforward and data-efficient, opening the door to applications in many challenging molecular systems.

Perspectives on Ligand/Protein Binding Kinetics Simulations: Force Fields, Machine Learning, Sampling, and User-Friendliness

Computational techniques applied to drug discovery have gained considerable popularity for their ability to filter potentially active drugs from inactive ones, reducing the time scale and costs of preclinical investigations. The main focus of these studies has historically been the search for compounds endowed with high affinity for a specific molecular target to ensure the formation of stable and long-lasting complexes. Recent evidence has also correlated the in vivo drug efficacy with its binding kinetics, thus opening new fascinating scenarios for ligand/protein binding kinetic simulations in drug discovery. The present article examines the state of the art in the field, providing a brief summary of the most popular and advanced ligand/protein binding kinetics techniques and evaluating their current limitations and the potential solutions to reach more accurate kinetic models. Particular emphasis is put on the need for a paradigm change in the present methodologies toward ligand and protein parametrization, the force field problem, characterization of the transition states, the sampling issue, and algorithms' performance, user-friendliness, and data openness.

Predicting rare events using neural networks and short-trajectory data

Estimating the likelihood, timing, and nature of events is a major goal of modeling stochastic dynamical systems. When the event is rare in comparison with the timescales of simulation and/or measurement needed to resolve the elemental dynamics, accurate prediction from direct observations becomes challenging. In such cases a more effective approach is to cast statistics of interest as solutions to Feynman-Kac equations (partial differential equations). Here, we develop an approach to solve Feynman-Kac equations by training neural networks on short-trajectory data. Our approach is based on a Markov approximation but otherwise avoids assumptions about the underlying model and dynamics. This makes it applicable to treating complex computational models and observational data. We illustrate the advantages of our method using a low-dimensional model that facilitates visualization, and this analysis motivates an adaptive sampling strategy that allows on-the-fly identification of and addition of data to regions important for predicting the statistics of interest. Finally, we demonstrate that we can compute accurate statistics for a 75-dimensional model of sudden stratospheric warming. This system provides a stringent test bed for our method.

Discovering Reaction Pathways, Slow Variables, and Committor Probabilities with Machine Learning

A significant challenge faced by atomistic simulations is the difficulty, and often impossibility, to sample the transitions between metastable states of the free-energy landscape associated with slow molecular processes. Importance-sampling schemes represent an appealing option to accelerate the underlying dynamics by smoothing out the relevant free-energy barriers, but require the definition of suitable reaction-coordinate (RC) models expressed in terms of compact low-dimensional sets of collective variables (CVs). While most computational studies of slow molecular processes have traditionally relied on educated guesses based on human intuition to reduce the dimensionality of the problem at hand, a variety of machine-learning (ML) algorithms have recently emerged as powerful alternatives to discover meaningful CVs capable of capturing the dynamics of the slowest degrees of freedom. Considering a simple paradigmatic situation in which the long-time dynamics is dominated by the transition between two known metastable states, we compare two variational data-driven ML methods based on Siamese neural networks aimed at discovering a meaningful RC model─the slowest decorrelating CV of the molecular process, and the committor probability to first reach one of the two metastable states. One method is the state-free reversible variational approach for Markov processes networks (VAMPnets), or SRVs─the other, inspired by the transition path theory framework, is the variational committor-based neural networks, or VCNs. The relationship and the ability of these methodologies to discover the relevant descriptors of the slow molecular process of interest are illustrated with a series of simple model systems. We also show that both strategies are amenable to importance-sampling schemes through an appropriate reweighting algorithm that approximates the kinetic properties of the transition.

Inexact iterative numerical linear algebra for neural network-based spectral estimation and rare-event prediction

Understanding dynamics in complex systems is challenging because there are many degrees of freedom, and those that are most important for describing events of interest are often not obvious. The leading eigenfunctions of the transition operator are useful for visualization, and they can provide an efficient basis for computing statistics, such as the likelihood and average time of events (predictions). Here, we develop inexact iterative linear algebra methods for computing these eigenfunctions (spectral estimation) and making predictions from a dataset of short trajectories sampled at finite intervals. We demonstrate the methods on a low-dimensional model that facilitates visualization and a high-dimensional model of a biomolecular system. Implications for the prediction problem in reinforcement learning are discussed.

Density functional theory and ab initio molecular dynamics reveal atomistic mechanisms for carbonate clumped isotope reordering

Carbon (¹³C) and oxygen (¹⁸O) isotopes in carbonates form clumped isotope species inversely correlated with temperature, providing a valuable paleothermometer for sedimentary carbonates and fossils. However, this signal resets ("reorders") with increasing temperature after burial. Research on reordering kinetics has characterized reordering rates and hypothesized the effects of impurities and trapped water, but the atomistic mechanism remains obscure. This work studies carbonate-clumped isotope reordering in calcite via first-principles simulations. We developed an atomistic view of the isotope exchange reaction between carbonate pairs in calcite, discovering a preferred configuration and elucidating how Mg²⁺ substitution and Ca²⁺ vacancies lower the free energy of activation (ΔA^‡) compared to pristine calcite. Regarding water-assisted isotopic exchange, the H⁺-O coordination distorts the transition state configuration and reduces ΔA^‡. We proposed a water-mediated exchange mechanism showing the lowest ΔA^‡ involving a reaction pathway with a hydroxylated four-coordinated carbon atom, confirming that internal water facilitates clumped isotope reordering.

Large barrier behavior of the rate constant from the diffusion equation

Many processes in chemistry, physics, and biology depend on thermally activated events in which the system changes its state by surmounting an activation barrier. Examples range from chemical reactions to protein folding and nucleation events. Parameterized forms of the mean field potential are often employed in the stochastic modeling of activated processes. In this contribution, we explore the alternative of employing parameterized forms of the equilibrium distribution by means of symmetric linear combination of two Gaussian functions. Such a procedure leads to flexible and convenient models for the landscape and the energy barrier whose features are controlled by the second moments of these Gaussian functions. The rate constants are examined through the solution of the corresponding diffusion problem, that is, the Fokker-Planck-Smoluchowski equation specified according to the parameterized equilibrium distribution. Numerical calculations clearly show that the asymptotic limit of large barriers does not agree with the results of the Kramers theory. The underlying reason is that the linear scaling of the potential, the procedure justifying the Kramers theory, cannot be applied when dealing with parameterized forms of the equilibrium distribution. A different kind of asymptotic analysis is then required and we introduce the appropriate theory when the equilibrium distribution is represented as a symmetric linear combination of two Gaussian functions: first in the one-dimensional case and afterward in the multidimensional diffusion model.

Generalized Transition State Theory Treatment of Water-Assisted Proton Transport Processes in Proteins

Transition state theory (TST) is widely employed for estimating the transition rate of a reaction when combined with free energy sampling techniques. A derivation of the transition theory rate expression for a general n-dimensional case is presented in this work which specifically focuses on water-assisted proton transfer/transport reactions, especially for protein systems. Our work evaluates the TST prefactor calculated at the transition state dividing surface compared to one sampled, as an approximation, in the reactant state in four case studies of water-assisted proton transport inside membrane proteins and highlights the significant impact of the prefactor position dependence in proton transport processes.

Committor-Consistent Variational String Method

The treatment of slow and rare transitions in the simulation of complex systems poses a great computational challenge. A powerful approach to tackle this challenge is the string method, which represents the transition path as a one-dimensional curve in a multidimensional space of collective variables. Commonly used strategies for pathway optimization include aligning the tangent of the string to the local mean force or to the mean drift determined from swarms of short trajectories. Here, a novel strategy is proposed, allowing the string to be optimized based on a variational principle involving the unidirectional reactive flux expressed in terms of the time-correlation function of the committor. The method is illustrated with model systems and then probed with the alanine dipeptide and a coarse-grained model of the barstar-barnase protein complex. Successive iterations variationally refine the string toward an optimal transition pathway following the gradient of the committor between two metastable states.

Relations among Unidirectional Fluxes at Equilibrium, Committors, and First Passage and Transition Path Times

For multidimensional diffusive dynamics, we algebraically derive remarkable analytical expressions that relate the mean first passage and transition path times between two dividing surfaces with the number of unidirectional transitions per unit time (fluxes) at equilibrium between the two surfaces and the committor (the probability of reaching one surface before the other). In one dimension, such relationships can be easily derived because analytical expressions for all the above-mentioned quantities can be found. This is not possible in higher dimensions, and at first sight, the problem seems much harder. We circumvent the difficulty that the equations determining the mean first passage and transition path times cannot be solved analytically by multiplying these equations by the committor, integrating both sides and finally using the divergence theorem. A byproduct of our derivation is an analytical expression for the starting point distribution over which individual first passage and transition path times must be averaged. It turns out that this distribution is not the Boltzmann one, but it has a simple physical interpretation.

Effects of Hydrodynamic Backflow on the Transmission Coefficient of a Barrier-Crossing Brownian Particle

The slow power law decay of the velocity autocorrelation function of a particle moving stochastically in a condensed-phase fluid is widely attributed to the momentum that fluid molecules displaced by the particle transfer back to it during the course of its motion. The forces created by this backflow effect are known as Basset forces, and they have been found in recent analytical work and numerical simulations to be implicated in a number of interesting dynamical phenomena, including boosted particle mobility in tilted washboard potentials. Motivated by these findings, the present paper is an investigation of the role of backflow in thermally activated barrier crossing, the governing process in essentially all condensed-phase chemical reactions. More specifically, it is an exact analytical calculation, carried out within the framework of the reactive-flux formalism, of the transmission coefficient κ(t) of a Brownian particle that crosses an inverted parabola under the influence of a colored noise process originating in the Basset force and a Markovian time-local friction. The calculation establishes that κ(t) is significantly enhanced over its backflow-free limit.