Trajectory-Based Rare Event Simulations. In: Lipkowitz KB, editor. Reviews in Computational Chemistry. Hoboken: John Wiley & Sons, Inc.; 2010. pp. 111-210. [DOI: 10.1002/9780470890905.ch3]

Unified perspective on exponential tilt and bridge algorithms for rare trajectories of discrete Markov processes

This article analyzes and compares two general techniques of rare event simulation for generating paths of Markov processes over fixed time horizons: exponential tilting and stochastic bridge. These two methods allow us to accurately compute the probability that a Markov process ends within a rare region which is unlikely to be attained. Exponential tilting is a general technique for obtaining an alternative or tilted sampling probability measure, under which the Markov process becomes likely to hit the rare region at terminal time. The stochastic bridge technique involves conditioning paths towards two endpoints: the terminal point and the initial one. The terminal point is generated from some appropriately chosen probability distribution that covers well the rare region. We show that both methods belong to the class of importance sampling procedures by providing a common mathematical framework of these two conceptually different methods of sampling rare trajectories. We also conduct a numerical comparison of these two methods, revealing distinct areas of application for each Monte Carlo method, where they exhibit superior efficiency. Detailed simulation algorithms are provided.

Optimizing molecular potential models by imposing kinetic constraints with path reweighting

Empirical force fields employed in molecular dynamics simulations of complex systems are often optimized to reproduce experimentally determined structural and thermodynamic properties. In contrast, experimental knowledge about the interconversion rates between metastable states in such systems is hardly ever incorporated in a force field due to a lack of an efficient approach. Here, we introduce such a framework based on the relationship between dynamical observables, such as rate constants, and the underlying molecular model parameters using the statistical mechanics of trajectories. Given a prior ensemble of molecular dynamics trajectories produced with imperfect force field parameters, the approach allows for the optimal adaption of these parameters such that the imposed constraint of equally predicted and experimental rate constant is obeyed. To do so, the method combines the continuum path ensemble maximum caliber approach with path reweighting methods for stochastic dynamics. When multiple solutions are found, the method selects automatically the combination that corresponds to the smallest perturbation of the entire path ensemble, as required by the maximum entropy principle. To show the validity of the approach, we illustrate the method on simple test systems undergoing rare event dynamics. Next to simple 2D potentials, we explore particle models representing molecular isomerization reactions and protein-ligand unbinding. Besides optimal interaction parameters, the methodology gives physical insights into what parts of the model are most sensitive to the kinetics. We discuss the generality and broad implications of the methodology.

Free Energy Landscapes, Diffusion Coefficients, and Kinetic Rates from Transition Paths

We address the problem of constructing accurate mathematical models of the dynamics of complex systems projected on a collective variable. To this aim we introduce a conceptually simple yet effective algorithm for estimating the parameters of Langevin and Fokker-Planck equations from a set of short, possibly out-of-equilibrium molecular dynamics trajectories, obtained for instance from transition path sampling or as relaxation from high free-energy configurations. The approach maximizes the model likelihood based on any explicit expression of the short-time propagator, hence it can be applied to different evolution equations. We demonstrate the numerical efficiency and robustness of the algorithm on model systems, and we apply it to reconstruct the projected dynamics of pairs of C₆₀ and C₂₄₀ fullerene molecules in explicit water. Our methodology allows reconstructing the accurate thermodynamics and kinetics of activated processes, namely free energy landscapes, diffusion coefficients, and kinetic rates. Compared to existing enhanced sampling methods, we directly exploit short unbiased trajectories at a competitive computational cost.

Practical guide to replica exchange transition interface sampling and forward flux sampling

Path sampling approaches have become invaluable tools to explore the mechanisms and dynamics of the so-called rare events that are characterized by transitions between metastable states separated by sizable free energy barriers. Their practical application, in particular to ever more complex molecular systems, is, however, not entirely trivial. Focusing on replica exchange transition interface sampling (RETIS) and forward flux sampling (FFS), we discuss a range of analysis tools that can be used to assess the quality and convergence of such simulations, which is crucial to obtain reliable results. The basic ideas of a step-wise evaluation are exemplified for the study of nucleation in several systems with different complexities, providing a general guide for the critical assessment of RETIS and FFS simulations.

Multitask Machine Learning of Collective Variables for Enhanced Sampling of Rare Events

Computing accurate reaction rates is a central challenge in computational chemistry and biology because of the high cost of free energy estimation with unbiased molecular dynamics. In this work, a data-driven machine learning algorithm is devised to learn collective variables with a multitask neural network, where a common upstream part reduces the high dimensionality of atomic configurations to a low dimensional latent space and separate downstream parts map the latent space to predictions of basin class labels and potential energies. The resulting latent space is shown to be an effective low-dimensional representation, capturing the reaction progress and guiding effective umbrella sampling to obtain accurate free energy landscapes. This approach is successfully applied to model systems including a 5D Müller Brown model, a 5D three-well model, the alanine dipeptide in vacuum, and an Au(110) surface reconstruction unit reaction. It enables automated dimensionality reduction for energy controlled reactions in complex systems, offers a unified and data-efficient framework that can be trained with limited data, and outperforms single-task learning approaches, including autoencoders.

Artificial Intelligence Resolves Kinetic Pathways of Magnesium Binding to RNA

Magnesium is an indispensable cofactor in countless vital processes. In order to understand its functional role, the characterization of the binding pathways to biomolecules such as RNA is crucial. Despite the importance, a molecular description is still lacking since the transition from the water-mediated outer-sphere to the direct inner-sphere coordination is on the millisecond time scale and therefore out of reach for conventional simulation techniques. To fill this gap, we use transition path sampling to resolve the binding pathways and to elucidate the role of the solvent in the binding process. The results reveal that the molecular void provoked by the leaving phosphate oxygen of the RNA is immediately filled by an entering water molecule. In addition, water molecules from the first and second hydration shell couple to the concerted exchange. To capture the intimate solute-solvent coupling, we perform a committor analysis as the basis for a machine learning algorithm that derives the optimal deep learning model from thousands of scanned architectures using hyperparameter tuning. The results reveal that the properly optimized deep network architecture recognizes the important solvent structures, extracts the relevant information, and predicts the commitment probability with high accuracy. Our results provide detailed insights into the solute-solvent coupling which is ubiquitous for kosmotropic ions and governs a large variety of biochemical reactions in aqueous solutions.

An extended autoencoder model for reaction coordinate discovery in rare event molecular dynamics datasets

The reaction coordinate (RC) is the principal collective variable or feature that determines the progress along an activated or reactive process. In a molecular simulation using enhanced sampling, a good description of the RC is crucial for generating sufficient statistics. Moreover, the RC provides invaluable atomistic insight into the process under study. The optimal RC is the committor, which represents the likelihood of a system to evolve toward a given state based on the coordinates of all its particles. As the interpretability of such a high dimensional function is low, a more practical approach is to describe the RC by some low-dimensional molecular collective variables or order parameters. While several methods can perform this dimensionality reduction, they usually require a preselection of these low-dimension collective variables (CVs). Here, we propose to automate this dimensionality reduction using an extended autoencoder, which maps the input (many CVs) onto a lower-dimensional latent space, which is subsequently used for the reconstruction of the input as well as the prediction of the committor function. As a consequence, the latent space is optimized for both reconstruction and committor prediction and is likely to yield the best non-linear low-dimensional representation of the committor. We test our extended autoencoder model on simple but nontrivial toy systems, as well as extensive molecular simulation data of methane hydrate nucleation. The extended autoencoder model can effectively extract the underlying mechanism of a reaction, make reliable predictions about the committor of a given configuration, and potentially even generate new paths representative for a reaction.

Homogenous nucleation rate of CO2 hydrates using transition interface sampling

Carbon dioxide and water can form solid clathrate structures in which water cages encapsulate the gas molecules. Such hydrates have sparked much interest due to their possible application in CO₂ sequestration. How the solid structure forms exactly from the liquid phase via a homogenous nucleation process is still poorly understood. This nucleation event is rare on the molecular timescale even under moderate undercooling or supersaturation conditions because of the large free energy barrier toward crystallization, rendering a brute force simulation of hydrate nucleation unfeasible for moderate undercooling or supersaturation. Here, we perform transition interface sampling simulations to quantify the homogenous nucleation rate for CO₂ hydrate formation using accurate atomistic force fields at 500 bars for three different temperatures between 260 and 273 K. Collecting more than 100 000 pathways comprising roughly two milliseconds of simulation time, we computed a nucleation rate in the amorphous phase of ∼10²¹ nuclei s^-1 cm^-3 for a temperature of 260 K and a rate of ∼10¹² nuclei s^-1 cm^-3 for a temperature of 265 K. For a temperature of 273 K, we find that the hydrate forms an sI crystalline phase with a rate of order of ∼10¹ nuclei s^-1 cm^-3. We compare these rates to classical nucleation theory estimates as well as experiments, and to nucleation rate estimates for methane hydrates and discuss possible causes of the observed differences. Our findings shed light on the kinetics of this important clathrate and should assist in future hydrate formation investigation.

Comparison of Carbohydrate Force Fields in Molecular Dynamics Simulations of Protein-Carbohydrate Complexes

In this paper, we present the results of molecular dynamics simulations aimed at critical comparison of classical, biomolecular force fields (FFs) in the context of their capabilities to describe the structural and thermodynamic features of carbohydrate-protein interactions. We have considered the three main families of FFs (CHARMM, GROMOS, and GLYCAM/AMBER) by applying them to investigate the seven different carbohydrate-protein complexes. The results indicate that although the qualitative pattern of several structural descriptors (intermolecular hydrogen bonding, ligand dynamic location, etc.) is conserved among the compared FFs, there also exists a number of significant divergences (mainly the patterns of contacts between particular amino acid residues and bound carbohydrate). The carbohydrate-protein unbinding free energies also vary from one FF to another, displaying diversified trends in deviations from the experimental data. The magnitude of those deviations is not negligible and indicates the need for refinement in the currently existing combinations of carbohydrate- and protein-dedicated biomolecular force fields. In spite of the lack of explicit functional terms responsible for the corresponding intermolecular forces, all tested FFs are capable of adequately reproducing the CH-π interactions, crucial for carbohydrate-protein binding.

Rate Prediction for Homogeneous Nucleation of Methane Hydrate at Moderate Supersaturation Using Transition Interface Sampling

The crystallization of methane hydrates via homogeneous nucleation under natural, moderate conditions is of both industrial and scientific relevance, yet still poorly understood. Predicting the nucleation rates at such conditions is notoriously difficult due to high nucleation barriers, and requires, besides an accurate molecular model, enhanced sampling. Here, we apply the transition interface sampling technique, which efficiently computes the exact rate of nucleation by generating ensembles of unbiased dynamical trajectories crossing predefined interfaces located between the stable states. Using an accurate atomistic force field and focusing on specific conditions of 280 K and 500 bar, we compute for nucleation directly into the sI crystal phase at a rate of ∼10^-17 nuclei per nanosecond per simulation volume or ∼10² nuclei per second per cm³, in agreement with consensus estimates for nearby conditions. As this is most likely fortuitous, we discuss the causes of the large differences between our results and previous simulation studies. Our work shows that it is now possible to compute rates for methane hydrates at moderate supersaturation, without relying on any assumptions other than the force field.

Kinetics of Drug Binding and Residence Time

The kinetics of drug binding and unbinding is assuming an increasingly crucial role in the long, costly process of bringing a new medicine to patients. For example, the time a drug spends in contact with its biological target is known as residence time (the inverse of the kinetic constant of the drug-target unbinding, 1/ koff). Recent reports suggest that residence time could predict drug efficacy in vivo, perhaps even more effectively than conventional thermodynamic parameters (free energy, enthalpy, entropy). There are many experimental and computational methods for predicting drug-target residence time at an early stage of drug discovery programs. Here, we review and discuss the methodological approaches to estimating drug binding kinetics and residence time. We first introduce the theoretical background of drug binding kinetics from a physicochemical standpoint. We then analyze the recent literature in the field, starting from the experimental methodologies and applications thereof and moving to theoretical and computational approaches to the kinetics of drug binding and unbinding. We acknowledge the central role of molecular dynamics and related methods, which comprise a great number of the computational methods and applications reviewed here. However, we also consider kinetic Monte Carlo. We conclude with the outlook that drug (un)binding kinetics may soon become a go/no go step in the discovery and development of new medicines.

Contour forward flux sampling: Sampling rare events along multiple collective variables

Many rare event transitions involve multiple collective variables (CVs), and the most appropriate combination of CVs is generally unknown a priori. We thus introduce a new method, contour forward flux sampling (cFFS), to study rare events with multiple CVs simultaneously. cFFS places nonlinear interfaces on-the-fly from the collective progress of the simulations, without any prior knowledge of the energy landscape or appropriate combination of CVs. We demonstrate cFFS on analytical potential energy surfaces and a conformational change in alanine dipeptide.

Sampling the equilibrium kinetic network of Trp-cage in explicit solvent

We employed the single replica multiple state transition interface sampling (MSTIS) approach to sample the kinetic (un)folding network of Trp-cage mini-protein in explicit water. Cluster analysis yielded 14 important metastable states in the network. The MSTIS simulation thus resulted in a full 14 × 14 rate matrix. Analysis of the kinetic rate matrix indicates the presence of a near native intermediate state characterized by a fully formed alpha helix, a slightly disordered proline tail, a broken salt-bridge, and a rotated arginine residue. This intermediate was also found in recent IR experiments. Moreover, the predicted rate constants and timescales are in agreement with previous experiments and simulations.