Committors, first-passage times, fluxes, Markov states, milestones, and all that

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.

High Nucleotide Skew Palindromic DNA Sequences Function as Potential Replication Origins due to their Unzipping Propensity

Locations of DNA replication initiation in prokaryotes, called "origins of replication", are well-characterized. However, a mechanistic understanding of the sequence dependence of the local unzipping of double-stranded DNA, the first step towards replication initiation, is lacking. Here, utilizing a Markov chain model that was created to address the directional nature of DNA unzipping and replication, we model the sequence dependence of local melting of double-stranded linear DNA segments. We show that generalized palindromic sequences with high nucleotide skews have a low kinetic barrier for local melting near melting temperatures. This allows for such sequences to function as potential replication origins. We support our claim with evidence for high-skew palindromic sequences within the replication origins of mitochondrial DNA, bacteria, archaea and plasmids.

Kinetic Network in Milestoning: Clustering, Reduction, and Transition Path Analysis

We present a reduction of the Milestoning (ReM) algorithm to analyze the high-dimensional Milestoning kinetic network. The algorithm reduces the Milestoning network to low dimensions but preserves essential kinetic information, such as local residence time, exit time, and mean first passage time between any two states. This is achieved in three steps. First, nodes (milestones) in the high-dimensional Milestoning network are grouped into clusters based on the metastability identified by an auxiliary continuous-time Markov chain. Our clustering method is applicable not only to time-reversible networks but also to nonreversible networks generated from practical simulations with statistical fluctuations. Second, a reduced network is established via network transformation, containing only the core sets of clusters as nodes. Finally, transition pathways are analyzed in the reduced network based on the transition path theory. The algorithm is illustrated using a toy model and a solvated alanine dipeptide in two and four dihedral angles.

QM/MM Free Energy Calculations of Long-Range Biological Protonation Dynamics by Adaptive and Focused Sampling

Water-mediated proton transfer reactions are central for catalytic processes in a wide range of biochemical systems, ranging from biological energy conversion to chemical transformations in the metabolism. Yet, the accurate computational treatment of such complex biochemical reactions is highly challenging and requires the application of multiscale methods, in particular hybrid quantum/classical (QM/MM) approaches combined with free energy simulations. Here, we combine the unique exploration power of new advanced sampling methods with density functional theory (DFT)-based QM/MM free energy methods for multiscale simulations of long-range protonation dynamics in biological systems. In this regard, we show that combining multiple walkers/well-tempered metadynamics with an extended system adaptive biasing force method (MWE) provides a powerful approach for exploration of water-mediated proton transfer reactions in complex biochemical systems. We compare and combine the MWE method also with QM/MM umbrella sampling and explore the sampling of the free energy landscape with both geometric (linear combination of proton transfer distances) and physical (center of excess charge) reaction coordinates and show how these affect the convergence of the potential of mean force (PMF) and the activation free energy. We find that the QM/MM-MWE method can efficiently explore both direct and water-mediated proton transfer pathways together with forward and reverse hole transfer mechanisms in the highly complex proton channel of respiratory Complex I, while the QM/MM-US approach shows a systematic convergence of selected long-range proton transfer pathways. In this regard, we show that the PMF along multiple proton transfer pathways is recovered by combining the strengths of both approaches in a QM/MM-MWE/focused US (FUS) scheme and reveals new mechanistic insight into the proton transfer principles of Complex I. Our findings provide a promising basis for the quantitative multiscale simulations of long-range proton transfer reactions in biological systems.

Information thermodynamics of transition paths between multiple mesostates

A central concern across the natural sciences is a quantitative understanding of the mechanism governing rare transitions between two metastable states. Recent research has uncovered a fundamental equality between the time-reversal asymmetry of the ensemble of such transition paths and the informativeness of system dynamics about the reactivity of a given trajectory, immediately leading to quantitative criteria for judging the importance of distinct system coordinates for the transition. Here we generalize this framework to multiple mesostates. We find that the main system-wide and coordinate-specific results generalize intuitively, while the combinatorial diversity of pairwise transitions raises new questions and points to new concepts. This work increases the previous framework's generality and applicability and forges connections to enhanced-sampling and coarse-grained dynamical approaches such as milestoning and Markov-state models.

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).

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.

Controlling Uncertainty of Empirical First-Passage Times in the Small-Sample Regime

We derive general bounds on the probability that the empirical first-passage time τ[over ¯]_{n}≡∑_{i=1}^{n}τ_{i}/n of a reversible ergodic Markov process inferred from a sample of n independent realizations deviates from the true mean first-passage time by more than any given amount in either direction. We construct nonasymptotic confidence intervals that hold in the elusive small-sample regime and thus fill the gap between asymptotic methods and the Bayesian approach that is known to be sensitive to prior belief and tends to underestimate uncertainty in the small-sample setting. We prove sharp bounds on extreme first-passage times that control uncertainty even in cases where the mean alone does not sufficiently characterize the statistics. Our concentration-of-measure-based results allow for model-free error control and reliable error estimation in kinetic inference, and are thus important for the analysis of experimental and simulation data in the presence of limited sampling.

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.

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.

The effect of time resolution on the observed first passage times in diffusive dynamics

Single-molecule and single-particle tracking experiments are typically unable to resolve fine details of thermal motion at short timescales where trajectories are continuous. We show that, when a diffusive trajectory xt is sampled at finite time intervals δt, the resulting error in measuring the first passage time to a given domain can exceed the time resolution of the measurement by more than an order of magnitude. Such surprisingly large errors originate from the fact that the trajectory may enter and exit the domain while being unobserved, thereby lengthening the apparent first passage time by an amount that is larger than δt. Such systematic errors are particularly important in single-molecule studies of barrier crossing dynamics. We show that the correct first passage times, as well as other properties of the trajectories such as splitting probabilities, can be recovered via a stochastic algorithm that reintroduces unobserved first passage events probabilistically.

Variational kinetic clustering of complex networks

Efficiently identifying the most important communities and key transition nodes in weighted and unweighted networks is a prevalent problem in a wide range of disciplines. Here, we focus on the optimal clustering using variational kinetic parameters, linked to Markov processes defined on the underlying networks, namely, the slowest relaxation time and the Kemeny constant. We derive novel relations in terms of mean first passage times for optimizing clustering via the Kemeny constant and show that the optimal clustering boundaries have equal round-trip times to the clusters they separate. We also propose an efficient method that first projects the network nodes onto a 1D reaction coordinate and subsequently performs a variational boundary search using a parallel tempering algorithm, where the variational kinetic parameters act as an energy function to be extremized. We find that maximization of the Kemeny constant is effective in detecting communities, while the slowest relaxation time allows for detection of transition nodes. We demonstrate the validity of our method on several test systems, including synthetic networks generated from the stochastic block model and real world networks (Santa Fe Institute collaboration network, a network of co-purchased political books, and a street network of multiple cities in Luxembourg). Our approach is compared with existing clustering algorithms based on modularity and the robust Perron cluster analysis, and the identified transition nodes are compared with different notions of node centrality.

Chasing collective variables using temporal data-driven strategies

The convergence of free-energy calculations based on importance sampling depends heavily on the choice of collective variables (CVs), which in principle, should include the slow degrees of freedom of the biological processes to be investigated. Autoencoders (AEs), as emerging data-driven dimension reduction tools, have been utilised for discovering CVs. AEs, however, are often treated as black boxes, and what AEs actually encode during training, and whether the latent variables from encoders are suitable as CVs for further free-energy calculations remains unknown. In this contribution, we review AEs and their time-series-based variants, including time-lagged AEs (TAEs) and modified TAEs, as well as the closely related model variational approach for Markov processes networks (VAMPnets). We then show through numerical examples that AEs learn the high-variance modes instead of the slow modes. In stark contrast, time series-based models are able to capture the slow modes. Moreover, both modified TAEs with extensions from slow feature analysis and the state-free reversible VAMPnets (SRVs) can yield orthogonal multidimensional CVs. As an illustration, we employ SRVs to discover the CVs of the isomerizations of N-acetyl-N'-methylalanylamide and trialanine by iterative learning with trajectories from biased simulations. Last, through numerical experiments with anisotropic diffusion, we investigate the potential relationship of time-series-based models and committor probabilities.

Modeling Ensembles of Enzyme Reaction Pathways with Hi-MSM Reveals the Importance of Accounting for Pathway Diversity

Conventional quantum mechanical-molecular mechanics (QM/MM) simulation approaches for modeling enzyme reactions often assume that there is one dominant reaction pathway and that this pathway can be sampled starting from an X-ray structure of the enzyme. These assumptions reduce computational cost; however, their validity has not been extensively tested. This is due in part to the lack of a rigorous formalism for integrating disparate pathway information from dynamical QM/MM calculations. Here, we present a way to model ensembles of reaction pathways efficiently using a divide-and-conquer strategy through Hierarchical Markov State Modeling (Hi-MSM). This approach allows information on multiple, distinct pathways to be incorporated into a chemical kinetic model, and it allows us to test these two assumptions. Applying Hi-MSM to the reaction carried out by dihydrofolate reductase (DHFR) we find (i) there are multiple, distinct pathways significantly contributing to the overall flux of the reaction that the conventional approach does not identify and (ii) that the conventional approach does not identify the dominant reaction pathway. Thus, both assumptions underpinning the conventional approach are violated. Since DHFR is a relatively small enzyme, and configuration space scales exponentially with protein size, accounting for multiple reaction pathways is likely to be necessary for most enzymes.

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.

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.

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.

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.

Information Thermodynamics of the Transition-Path Ensemble

The reaction coordinate describing a transition between reactant and product is a fundamental concept in the theory of chemical reactions. Within transition-path theory, a quantitative definition of the reaction coordinate is found in the committor, which is the probability that a trajectory initiated from a given microstate first reaches the product before the reactant. Here we develop an information-theoretic origin for the committor and show how selecting transition paths from a long ergodic equilibrium trajectory induces entropy production which exactly equals the information that system dynamics provide about the reactivity of trajectories. This equality of entropy production and dynamical information generation also holds at the level of arbitrary individual coordinates, providing parallel measures of the coordinate's relevance to the reaction, each of which is maximized by the committor.

Stochastic Model of Solvent Exchange in the First Coordination Shell of Aqua Ions

Ion microsolvation is a basic, yet fundamental, process of ionic solutions underlying many relevant phenomena in either biological or nanotechnological applications, such as solvent reorganization energy, ion transport, catalytic activity, and so on. As a consequence, it is a topic of extensive investigations by various experimental techniques, ranging from X-ray diffraction to NMR relaxation and from calorimetry to vibrational spectroscopy, and theoretical approaches, especially those based on molecular dynamics (MD) simulations. The conventional microscopic view of ion solvation is usually provided by a "static" cluster model representing the first ion-solvent coordination shell. Despite the merits of such a simple model, however, ion coordination in solution should be better regarded as a complex population of dynamically interchanging molecular configurations. Such a more comprehensive view is more subtle to characterize and often elusive to standard approaches. In this work, we report on an effective computational strategy aiming at providing a detailed picture of solvent coordination and exchange around aqua ions, thus including the main structural, thermodynamic, and dynamic properties of ion microsolvation, such as the most probable first-shell complex structures, the corresponding free energies, the interchanging energy barriers, and the solvent-exchange rates. Assuming the solvent coordination number as an effective reaction coordinate and combining MD simulations with enhanced sampling and master-equation approaches, we propose a stochastic model suitable for properly describing, at the same time, the thermodynamics and kinetics of ion-water coordination. The model is successfully tested toward various divalent ions (Ca²⁺, Zn²⁺, Hg²⁺, and Cd²⁺) in aqueous solution, considering also the case of a high ionic concentration. Results show a very good agreement with those issuing from brute-force MD simulations, when available, and support the reliable prediction of rare ion-water complexes and slow water exchange rates not easily accessible to usual computational methods.

Transition rate theory, spectral analysis, and reactive paths

The kinetics of a dynamical system dominated by two metastable states is examined from the perspective of the activated-dynamics reactive flux formalism, Markov state eigenvalue spectral decomposition, and committor-based transition path theory. Analysis shows that the different theoretical formulations are consistent, clarifying the significance of the inherent microscopic lag-times that are implicated, and that the most meaningful one-dimensional reaction coordinate in the region of the transition state is along the gradient of the committor in the multidimensional subspace of collective variables. It is shown that the familiar reactive flux activated dynamics formalism provides an effective route to calculate the transition rate in the case of a narrow sharp barrier but much less so in the case of a broad flat barrier. In this case, the standard reactive flux correlation function decays very slowly to the plateau value that corresponds to the transmission coefficient. Treating the committor function as a reaction coordinate does not alleviate all issues caused by the slow relaxation of the reactive flux correlation function. A more efficient activated dynamics simulation algorithm may be achieved from a modified reactive flux weighted by the committor. Simulation results on simple systems are used to illustrate the various conceptual points.

A Companion Guide to the String Method with Swarms of Trajectories: Characterization, Performance, and Pitfalls

The string method with swarms of trajectories (SMwST) is an algorithm that identifies a physically meaningful transition pathway─a one-dimensional curve, embedded within a high-dimensional space of selected collective variables. The SMwST algorithm leans on a series of short, unbiased molecular dynamics simulations spawned at different locations of the discretized path, from whence an average dynamic drift is determined to evolve the string toward an optimal pathway. However conceptually simple in both its theoretical formulation and practical implementation, the SMwST algorithm is computationally intensive and requires a careful choice of parameters for optimal cost-effectiveness in applications to challenging problems in chemistry and biology. In this contribution, the SMwST algorithm is presented in a self-contained manner, discussing with a critical eye its theoretical underpinnings, applicability, inherent limitations, and use in the context of path-following free-energy calculations and their possible extension to kinetics modeling. Through multiple simulations of a prototypical polypeptide, combining the search of the transition pathway and the computation of the potential of mean force along it, several practical aspects of the methodology are examined with the objective of optimizing the computational effort, yet without sacrificing accuracy. In light of the results reported here, we propose some general guidelines aimed at improving the efficiency and reliability of the computed pathways and free-energy profiles underlying the conformational transitions at hand.

Nonlinear Reaction Coordinate of an Enzyme Catalyzed Proton Transfer Reaction

We present an in-depth study on the theoretical calculation of an optimum reaction coordinate as a linear or nonlinear combination of important collective variables (CVs) sampled from an ensemble of reactive transition paths for an intramolecular proton transfer reaction catalyzed by the enzyme human carbonic anhydrase (HCA) II. The linear models are optimized by likelihood maximization for a given number of CVs. The nonlinear models are based on an artificial neural network with the same number of CVs and optimized by minimizing the root-mean-square error in comparison to a training set of committor estimators generated for the given transition. The nonlinear reaction coordinate thus obtained yields the free energy of activation and rate constant as 9.46 kcal mol^-1 and 1.25 × 10⁶ s^-1, respectively. These estimates are found to be in quantitative agreement with the known experimental results. We have also used an extended autoencoder model to show that a similar analysis can be carried out using a single CV only. The resultant free energies and kinetics of the reaction slightly overestimate the experimental data. The implications of these results are discussed using a detailed microkinetic scheme of the proton transfer reaction catalyzed by HCA II.

Markovian Weighted Ensemble Milestoning (M-WEM): Long-Time Kinetics from Short Trajectories

We introduce a rare-event sampling scheme, named Markovian Weighted Ensemble Milestoning (M-WEM), which inlays a weighted ensemble framework within a Markovian milestoning theory to efficiently calculate thermodynamic and kinetic properties of long-time-scale biomolecular processes from short atomistic molecular dynamics simulations. M-WEM is tested on the Müller-Brown potential model, the conformational switching in alanine dipeptide, and the millisecond time-scale protein-ligand unbinding in a trypsin-benzamidine complex. Not only can M-WEM predict the kinetics of these processes with quantitative accuracy but it also allows for a scheme to reconstruct a multidimensional free-energy landscape along additional degrees of freedom, which are not part of the milestoning progress coordinate. For the ligand-receptor system, the experimental residence time, association and dissociation kinetics, and binding free energy could be reproduced using M-WEM within a simulation time of a few hundreds of nanoseconds, which is a fraction of the computational cost of other currently available methods, and close to 4 orders of magnitude less than the experimental residence time. Due to the high accuracy and low computational cost, the M-WEM approach can find potential applications in kinetics and free-energy-based computational drug design.

Nearly reducible finite Markov chains: Theory and algorithms

Finite Markov chains, memoryless random walks on complex networks, appear commonly as models for stochastic dynamics in condensed matter physics, biophysics, ecology, epidemiology, economics, and elsewhere. Here, we review exact numerical methods for the analysis of arbitrary discrete- and continuous-time Markovian networks. We focus on numerically stable methods that are required to treat nearly reducible Markov chains, which exhibit a separation of characteristic timescales and are therefore ill-conditioned. In this metastable regime, dense linear algebra methods are afflicted by propagation of error in the finite precision arithmetic, and the kinetic Monte Carlo algorithm to simulate paths is unfeasibly inefficient. Furthermore, iterative eigendecomposition methods fail to converge without the use of nontrivial and system-specific preconditioning techniques. An alternative approach is provided by state reduction procedures, which do not require additional a priori knowledge of the Markov chain. Macroscopic dynamical quantities, such as moments of the first passage time distribution for a transition to an absorbing state, and microscopic properties, such as the stationary, committor, and visitation probabilities for nodes, can be computed robustly using state reduction algorithms. The related kinetic path sampling algorithm allows for efficient sampling of trajectories on a nearly reducible Markov chain. Thus, all of the information required to determine the kinetically relevant transition mechanisms, and to identify the states that have a dominant effect on the global dynamics, can be computed reliably even for computationally challenging models. Rare events are a ubiquitous feature of realistic dynamical systems, and so the methods described herein are valuable in many practical applications.

String Method with Swarms-of-Trajectories, Mean Drifts, Lag Time, and Committor

The kinetics of a dynamical system comprising two metastable states is formulated in terms of a finite-time propagator in phase space (position and velocity) adapted to the underdamped Langevin equation. Dimensionality reduction to a subspace of collective variables yields familiar expressions for the propagator, committor, and steady-state flux. A quadratic expression for the steady-state flux between the two metastable states can serve as a robust variational principle to determine an optimal approximate committor expressed in terms of a set of collective variables. The theoretical formulation is exploited to clarify the foundation of the string method with swarms-of-trajectories, which relies on the mean drift of short trajectories to determine the optimal transition pathway. It is argued that the conditions for Markovity within a subspace of collective variables may not be satisfied with an arbitrary short time-step and that proper kinetic behaviors appear only when considering the effective propagator for longer lag times. The effective propagator with finite lag time is amenable to an eigenvalue-eigenvector spectral analysis, as elaborated previously in the context of position-based Markov models. The time-correlation functions calculated by swarms-of-trajectories along the string pathway constitutes a natural extension of these developments. The present formulation provides a powerful theoretical framework to characterize the optimal pathway between two metastable states of a system.

Numerical analysis of first-passage processes in finite Markov chains exhibiting metastability

We describe state-reduction algorithms for the analysis of first-passage processes in discrete- and continuous-time finite Markov chains. We present a formulation of the graph transformation algorithm that allows for the evaluation of exact mean first-passage times, stationary probabilities, and committor probabilities for all nonabsorbing nodes of a Markov chain in a single computation. Calculation of the committor probabilities within the state-reduction formalism is readily generalizable to the first hitting problem for any number of alternative target states. We then show that a state-reduction algorithm can be formulated to compute the expected number of times that each node is visited along a first-passage path. Hence, all properties required to analyze the first-passage path ensemble (FPPE) at both a microscopic and macroscopic level of detail, including the mean and variance of the first-passage time distribution, can be computed using state-reduction methods. In particular, we derive expressions for the probability that a node is visited along a direct transition path, which proceeds without returning to the initial state, considering both the nonequilibrium and equilibrium (steady-state) FPPEs. The reactive visitation probability provides a rigorous metric to quantify the dynamical importance of a node for the productive transition between two endpoint states and thus allows the local states that facilitate the dominant transition mechanisms to be readily identified. The state-reduction procedures remain numerically stable even for Markov chains exhibiting metastability, which can be severely ill-conditioned. The rare event regime is frequently encountered in realistic models of dynamical processes, and our methodology therefore provides valuable tools for the analysis of Markov chains in practical applications. We illustrate our approach with numerical results for a kinetic network representing a structural transition in an atomic cluster.

Diffusive barrier crossing rates from variationally determined eigenvalues

Kramers' procedure for calculating the rate of activated processes involves partitioning space into reactant, barrier, and product regions by introducing two dividing surfaces. Then, a nonequilibrium steady state is established by injecting particles on one surface and removing them when they reach the other. The rate is obtained as the ratio of the steady-state flux between the surfaces and the population of the initial well. An alternative procedure that seems less artificial is to estimate the first non-zero eigenvalue of the operator that describes the dynamics and then equate its magnitude to the sum of the forward and backward rate constants. Here, we establish the relationship between these approaches for diffusive dynamics, starting with the variational principle for the eigenvalue of interest and then using a trial function involving two adjustable surfaces. We show how Kramers' flux-over-population expression for the rate constant can be obtained from our variationally determined eigenvalue in the special case where the reactant and product regions are separated by a high barrier. This work exploits the modern theory of activated rate processes where the committor (the probability of reaching one dividing surface before the other) plays a central role. Surprisingly, our upper bound for the eigenvalue can be expressed solely in terms of mean first-passage times and the mean transition-path time between the two dividing surfaces.

Graph transformation and shortest paths algorithms for finite Markov chains

The graph transformation (GT) algorithm robustly computes the mean first-passage time to an absorbing state in a finite Markov chain. Here we present a concise overview of the iterative and block formulations of the GT procedure and generalize the GT formalism to the case of any path property that is a sum of contributions from individual transitions. In particular, we examine the path action, which directly relates to the path probability, and analyze the first-passage path ensemble for a model Markov chain that is metastable and therefore numerically challenging. We compare the mean first-passage path action, obtained using GT, with the full path action probability distribution simulated efficiently using kinetic path sampling, and with values for the highest-probability paths determined by the recursive enumeration algorithm (REA). In Markov chains representing realistic dynamical processes, the probability distributions of first-passage path properties are typically fat-tailed and therefore difficult to converge by sampling, which motivates the use of exact and numerically stable approaches to compute the expectation. We find that the kinetic relevance of the set of highest-probability paths depends strongly on the metastability of the Markov chain, and so the properties of the dominant first-passage paths may be unrepresentative of the global dynamics. Use of a global measure for edge costs in the REA, based on net productive fluxes, allows the total reactive flux to be decomposed into a finite set of contributions from simple flux paths. By considering transition flux paths, a detailed quantitative analysis of the relative importance of competing dynamical processes is possible even in the metastable regime.

Electrophysiological Properties from Computations at a Single Voltage: Testing Theory with Stochastic Simulations

We use stochastic simulations to investigate the performance of two recently developed methods for calculating the free energy profiles of ion channels and their electrophysiological properties, such as current–voltage dependence and reversal potential, from molecular dynamics simulations at a single applied voltage. These methods require neither knowledge of the diffusivity nor simulations at multiple voltages, which greatly reduces the computational effort required to probe the electrophysiological properties of ion channels. They can be used to determine the free energy profiles from either forward or backward one-sided properties of ions in the channel, such as ion fluxes, density profiles, committor probabilities, or from their two-sided combination. By generating large sets of stochastic trajectories, which are individually designed to mimic the molecular dynamics crossing statistics of models of channels of trichotoxin, p7 from hepatitis C and a bacterial homolog of the pentameric ligand-gated ion channel, GLIC, we find that the free energy profiles obtained from stochastic simulations corresponding to molecular dynamics simulations of even a modest length are burdened with statistical errors of only 0.3 kcal/mol. Even with many crossing events, applying two-sided formulas substantially reduces statistical errors compared to one-sided formulas. With a properly chosen reference voltage, the current–voltage curves can be reproduced with good accuracy from simulations at a single voltage in a range extending for over 200 mV. If possible, the reference voltages should be chosen not simply to drive a large current in one direction, but to observe crossing events in both directions.

An effective one-dimensional approach to calculating mean first passage time in multi-dimensional potentials

Thermally activated escape processes in multi-dimensional potentials are of interest to a variety of fields, so being able to calculate the rate of escape-or the mean first-passage time (MFPT)-is important. Unlike in one dimension, there is no general, exact formula for the MFPT. However, Langer's formula, a multi-dimensional generalization of Kramers's one-dimensional formula, provides an approximate result when the barrier to escape is large. Kramers's and Langer's formulas are related to one another by the potential of mean force (PMF): when calculated along a particular direction (the unstable mode at the saddle point) and substituted into Kramers's formula, the result is Langer's formula. We build on this result by using the PMF in the exact, one-dimensional expression for the MFPT. Our model offers better agreement with Brownian dynamics simulations than Langer's formula, although discrepancies arise when the potential becomes less confining along the direction of escape. When the energy barrier is small our model offers significant improvements upon Langer's theory. Finally, the optimal direction along which to evaluate the PMF no longer corresponds to the unstable mode at the saddle point.

Optimal dimensionality reduction of Markov chains using graph transformation

Markov chains can accurately model the state-to-state dynamics of a wide range of complex systems, but the underlying transition matrix is ill-conditioned when the dynamics feature a separation of timescales. Graph transformation (GT) provides a numerically stable method to compute exact mean first passage times (MFPTs) between states, which are the usual dynamical observables in continuous-time Markov chains (CTMCs). Here, we generalize the GT algorithm to discrete-time Markov chains (DTMCs), which are commonly estimated from simulation data, for example, in the Markov state model approach. We then consider the dimensionality reduction of CTMCs and DTMCs, which aids model interpretation and facilitates more expensive computations, including sampling of pathways. We perform a detailed numerical analysis of existing methods to compute the optimal reduced CTMC, given a partitioning of the network into metastable communities (macrostates) of nodes (microstates). We show that approaches based on linear algebra encounter numerical problems that arise from the requisite metastability. We propose an alternative approach using GT to compute the matrix of intermicrostate MFPTs in the original Markov chain, from which a matrix of weighted intermacrostate MFPTs can be obtained. We also propose an approximation to the weighted-MFPT matrix in the strongly metastable limit. Inversion of the weighted-MFPT matrix, which is better conditioned than the matrices that must be inverted in alternative dimensionality reduction schemes, then yields the optimal reduced Markov chain. The superior numerical stability of the GT approach therefore enables us to realize optimal Markovian coarse-graining of systems with rare event dynamics.

Membrane permeability of small molecules from unbiased molecular dynamics simulations

Permeation of many small molecules through lipid bilayers can be directly observed in molecular dynamics simulations on the nano- and microsecond timescale. While unbiased simulations provide an unobstructed view of the permeation process, their feasibility for computing permeability coefficients depends on various factors that differ for each permeant. The present work studies three small molecules for which unbiased simulations of permeation are feasible within less than a microsecond, one hydrophobic (oxygen), one hydrophilic (water), and one amphiphilic (ethanol). Permeabilities are computed using two approaches: counting methods and a maximum-likelihood estimation for the inhomogeneous solubility diffusion (ISD) model. Counting methods yield nearly model-free estimates of the permeability for all three permeants. While the ISD-based approach is reasonable for oxygen, it lacks precision for water due to insufficient sampling and results in misleading estimates for ethanol due to invalid model assumptions. It is also demonstrated that simulations using a Langevin thermostat with collision frequencies of 1/ps and 5/ps yield oxygen permeabilities and diffusion constants that are lower than those using Nosé-Hoover by statistically significant margins. In contrast, permeabilities from trajectories generated with Nosé-Hoover and the microcanonical ensemble do not show statistically significant differences. As molecular simulations become more affordable and accurate, calculation of permeability for an expanding range of molecules will be feasible using unbiased simulations. The present work summarizes theoretical underpinnings, identifies pitfalls, and develops best practices for such simulations.

Dynamical matrix propagator scheme for large-scale proton dynamics simulations

We derive a matrix formalism for the simulation of long range proton dynamics for extended systems and timescales. On the basis of an ab initio molecular dynamics simulation, we construct a Markov chain, which allows us to store the entire proton dynamics in an M × M transition matrix (where M is the number of oxygen atoms). In this article, we start from common topology features of the hydrogen bond network of good proton conductors and utilize them as constituent constraints of our dynamic model. We present a thorough mathematical derivation of our approach and verify its uniqueness and correct asymptotic behavior. We propagate the proton distribution by means of transition matrices, which contain kinetic data from both ultra-short (sub-ps) and intermediate (ps) timescales. This concept allows us to keep the most relevant features from the microscopic level while effectively reaching larger time and length scales. We demonstrate the applicability of the transition matrices for the description of proton conduction trends in proton exchange membrane materials.

Defining, Calculating, and Converging Observables of a Kinetic Transition Network

Kinetic transition networks (KTNs) of local minima and transition states are able to capture the dynamics of numerous systems in chemistry, biology, and materials science. However, extracting observables is numerically challenging for large networks and generally will be sensitive to additional computational discovery. To have any measure of convergence for observables, these sensitivities must be regularly calculated. We present a matrix formulation of the discrete path sampling framework for KTNs, deriving expressions for branching probabilities, transition rates, and waiting times. Using the concept of the quasi-stationary distribution, a clear hierarchy of expressions for network observables is established, from exact results to steady-state approximations. We use these results in combination with the graph transformation method to derive the sensitivity, with respect to perturbations of the known KTN, giving explicit terms for the pairwise sensitivity and discussing the pathwise sensitivity. These results provide guidelines for converging the network, with respect to additional sampling, focusing on the estimates obtained for the overall rate coefficients between product and reactant states. We demonstrate this procedure for transitions in the double-funnel landscape of the 38-atom Lennard-Jones cluster.

Correlation functions, mean first passage times, and the Kemeny constant

Markov processes are widely used models for investigating kinetic networks. Here, we collate and present a variety of results pertaining to kinetic network models in a unified framework. The aim is to lay out explicit links between several important quantities commonly studied in the field, including mean first passage times (MFPTs), correlation functions, and the Kemeny constant. We provide new insights into (i) a simple physical interpretation of the Kemeny constant, (ii) a relationship to infer equilibrium distributions and rate matrices from measurements of MFPTs, and (iii) a protocol to reduce the dimensionality of kinetic networks based on specific requirements that the MFPTs in the coarse-grained system should satisfy. We prove that this protocol coincides with the one proposed by Hummer and Szabo [J. Phys. Chem. B 119, 9029 (2014)], and it leads to a variational principle for the Kemeny constant. Finally, we introduce a modification of this protocol, which preserves the Kemeny constant. Our work underpinning the theoretical aspects of kinetic networks will be useful in applications including milestoning and path sampling algorithms in molecular simulations.

Hierarchical Markov State Model Building to Describe Molecular Processes

Markov state models can describe ensembles of pathways via kinetic networks but are difficult to create when large free-energy barriers limit unbiased sampling. Chain-of-states simulations allow sampling over large free-energy barriers but are often constructed using a single pathway that is unlikely to thermodynamically average over orthogonal degrees of freedom in complex systems. Here, we combine the advantages of these two approaches in the form of a Markov state model of Markov state models, which we call a Hierarchical Markov state model. In this approach, independent Markov models are constructed in regions of configuration space that are locally well sampled but are separated by large free-energy barriers from other regions. A string method is used to construct an ensemble of pathways connecting the states of these different local Markov models, and the rate through each pathway is then estimated. These rates are then combined with the rate information from the local Markov models in a master equation to predict global rates, fluxes, and populations. By applying this hierarchical approach to tractable systems, a toy potential and dipeptides, we demonstrate that it is more accurate than the conventional single-pathway description. The advantages of this approach are that it (i) is more realistic than the conventional chain-of-states approach, as an ensemble of pathways rather than a single pathway is used to describe processes in high-dimensional systems, and (ii) it resolves the issue of poor sampling in Markov State model building when large free-energy barriers are present. The divide-and-conquer strategy inherent to this approach should make this procedure straightforward to apply to more complex systems.

Long-time methods for molecular dynamics simulations: Markov State Models and Milestoning

Molecular dynamics (MD) studies of biomolecules require the ability to simulate complex biochemical systems with an increasingly larger number of particles and for longer time scales, a problem that cannot be overcome by computational hardware advances alone. A main problem springs from the intrinsically high-dimensional and complex nature of the underlying free energy landscape of most systems, and from the necessity to sample accurately such landscapes for identifying kinetic and thermodynamic states in the configurations space, and for accurate calculations of both free energy differences and of the corresponding transition rates between states. Here, we review and present applications of two increasingly popular methods that allow long-time MD simulations of biomolecular systems that can open a broad spectrum of new studies. A first approach, Markov State Models (MSMs), relies on identifying a set of configuration states in which the system resides sufficiently long to relax and loose the memory of previous transitions, and on using simulations for mapping the underlying complex energy landscape and for extracting accurate thermodynamic and kinetic information. The Markovian independence of the underlying transition probabilities creates the opportunity to increase the sampling efficiency by using sets of appropriately initialized short simulations rather than typically long MD trajectories, which also enhances sampling. This allows MSM-based studies to unveil bio-molecular mechanisms and to estimate free energy barriers with high accuracy, in a manner that is both systematic and relatively automatic, which accounts for their increasing popularity. The second approach presented, Milestoning, targets accurate studies of the ensemble of pathways connecting specific end-states (e.g., reactants and products) in a similarly systematic, accurate and highly automatic manner. Applications presented range from studies of conformational dynamics and binding of amyloid-forming peptides, cell-penetrating peptides and the DFG-flip dynamics in Abl kinase. As highlighted by the increasing number of studies using both methods, we anticipate that they will open new avenues for the investigation of systematic sampling of reactions pathways and mechanisms occurring on longer time scales than currently accessible by purely computational hardware developments.