On the calculation of reaction rate constants in the transition path ensemble

Stochastic transition states: Reaction geometry amidst noise

Classical transition state theory (TST) is the cornerstone of reaction-rate theory. It postulates a partition of phase space into reactant and product regions, which are separated by a dividing surface that reactive trajectories must cross. In order not to overestimate the reaction rate, the dynamics must be free of recrossings of the dividing surface. This no-recrossing rule is difficult (and sometimes impossible) to enforce, however, when a chemical reaction takes place in a fluctuating environment such as a liquid. High-accuracy approximations to the rate are well known when the solvent forces are treated using stochastic representations, though again, exact no-recrossing surfaces have not been available. To generalize the exact limit of TST to reactive systems driven by noise, we introduce a time-dependent dividing surface that is stochastically moving in phase space, such that it is crossed once and only once by each transition path.

Energy landscapes and properties of biomolecules

Thermodynamic and dynamic properties of biomolecules can be calculated using a coarse-grained approach based upon sampling stationary points of the underlying potential energy surface. The superposition approximation provides an overall partition function as a sum of contributions from the local minima, and hence functions such as internal energy, entropy, free energy and the heat capacity. To obtain rates we must also sample transition states that link the local minima, and the discrete path sampling method provides a systematic means to achieve this goal. A coarse-grained picture is also helpful in locating the global minimum using the basin-hopping approach. Here we can exploit a fictitious dynamics between the basins of attraction of local minima, since the objective is to find the lowest minimum, rather than to reproduce the thermodynamics or dynamics.

Simultaneous computation of free energies and kinetics of rare events

We introduce a method to evaluate simultaneously the reaction rate constant and the free energy profile of a process in a complex environment. The method employs the partial path transition interface sampling technique we recently developed for the calculation of rate constants in diffusive systems. We illustrate the applicability of the technique by studying a simple dimer in a repulsive fluid, and show that the free energy can be obtained at no additional computational cost.

A proficient enzyme: insights on the mechanism of orotidine monophosphate decarboxylase from computer simulations

Decarboxylation of orotidine 5'-monophosphate (Omp) to uridine 5'-monophosphate by orotidine 5'-monophosphate decarboxylase (ODCase) is currently the object of vivid debate. Here, we clarify its enzymatic activity with long time scale classical molecular dynamics and hybrid ab initio Car-Parrinello/molecular mechanics simulations. The lack of structural (experimental) information on the ground state of ODCase/Omp complex is overcome by a careful construction of the model and the analysis of three different strains of the enzyme. We find that the ODCase/substrate complex is characterized by a very stable charged network Omp-Lys-Asp-Lys-Asp, which is incompatible with the previously proposed direct decarboxylation driven by a ground-state destabilization. A direct decarboxylation induced by a transition-state electrostatic stabilization is consistent with our findings. The calculated activation free energy for the direct decarboxylation with the formation of a C6 carboanionic intermediate yields an overall rate enhancement by the enzyme (k(cat)/k(wat) = 3.5 x 10(16)) in agreement with experiments (k(cat)/k(wat) = 1.7 x 10(17)). The decarboxylation is accompanied by the movement of a fully conserved lysine residue toward the developing negative charge at the C6 position.

Adaptive importance sampling Monte Carlo simulation of rare transition events

We develop a general theoretical framework for the recently proposed importance sampling method for enhancing the efficiency of rare-event simulations [W. Cai, M. H. Kalos, M. de Koning, and V. V. Bulatov, Phys. Rev. E 66, 046703 (2002)], and discuss practical aspects of its application. We define the success/fail ensemble of all possible successful and failed transition paths of any duration and demonstrate that in this formulation the rare-event problem can be interpreted as a "hit-or-miss" Monte Carlo quadrature calculation of a path integral. The fact that the integrand contributes significantly only for a very tiny fraction of all possible paths then naturally leads to a "standard" importance sampling approach to Monte Carlo (MC) quadrature and the existence of an optimal importance function. In addition to showing that the approach is general and expected to be applicable beyond the realm of Markovian path simulations, for which the method was originally proposed, the formulation reveals a conceptual analogy with the variational MC (VMC) method. The search for the optimal importance function in the former is analogous to finding the ground-state wave function in the latter. In two model problems we discuss practical aspects of finding a suitable approximation for the optimal importance function. For this purpose we follow the strategy that is typically adopted in VMC calculations: the selection of a trial functional form for the optimal importance function, followed by the optimization of its adjustable parameters. The latter is accomplished by means of an adaptive optimization procedure based on a combination of steepest-descent and genetic algorithms.

Mechanism of Ligand Exchange Studied Using Transition Path Sampling

The mechanism of intermolecular ligand exchange has been studied using transition path sampling (TPS) based molecular dynamics (MD) simulations. Specifically, the exchange of solvent molecules bound to unsaturated Cr(CO)5 in methanol solution has been investigated. The results of the TPS simulations have shown that there are multiple steps in the reaction mechanism. The first involves partial dissociation of the coordinated solvent from the Cr metal center followed by association with a new methanol molecule between the normally void first and second solvent layers. After diffusive motion of the exchanging ligands, the last step involves the originally bound methanol molecule moving into the bath continuum followed by solvation of the Cr metal fragment by the exchanging ligand. It has been found that the reaction center (defined as the organometallic fragment and two exchanging ligands only) and the solvent bath have favorable interactions. This is likely due to the adiabatic nature of the ligand exchange transition. The ability to understand the microscopic molecular dynamics of a chemical process based on a free energy analysis is also discussed.

Tools for channels: moving towards molecular calculations of gating and permeation in ion channel biophysics

Recent X-ray structures of voltage gated potassium channels provide an exciting opportunity to connect molecular structures with measured biological function. Two of the most important connections for these channels are: first, to the molecular basis behind selectivity and the associated free energy profile underlying ionic current flow and, second, to a true molecular understanding of the large-scale conformational transitions that underlie voltage dependent gating. But, existing computational tools need to be further developed to reach these goals. In this contribution to the symposia on sampling methods we outline our dynamic importance sampling method for sampling large-scale conformational transitions as well as our studies with non-equilibrium work events and equilibrium overlap sampling (OS) methods for sampling events related to the calculation of relative free energies.

Kinetic pathways of beta-hairpin (un)folding in explicit solvent

We examine the dynamical (un)folding pathways of the C-terminal beta-hairpin of protein G-B1 at room temperature in explicit solvent, by employing transition path sampling algorithms. The path ensembles contain information on the folding kinetics, including solvent motion. We determine the transition state ensembles for the two main transitions: 1), the hydrophobic collapse; and 2), the backbone hydrogen bond formation. In both cases the transition state ensembles are characterized by a layer (1) or a strip (2) of water molecules in between the two hairpin strands, supporting the hypothesis of the solvent as lubricant in the folding process. The transition state ensembles do not correspond with saddle points in the equilibrium free-energy landscapes. The kinetic pathways are thus not completely determined by the free-energy landscape. This phenomenon can occur if the order parameters obey different timescales. Using the transition interface sampling technique, we calculate the rate constants for (un)folding and find them in reasonable agreement with experiments, thus supporting the validation of using all-atom force fields to study protein folding.

Path integral approach to Brownian motion driven with an ac force

Brownian motion in a periodic potential driven by an ac (oscillatory) force is investigated for the full range of damping constant from the overdamped limit to the underdamped limit. The path (functional) integral approach is advanced to produce formulas for the probability distribution function and for the current of the Brownian particle in response to an ac driving force. The negative friction Langevin dynamics technique is employed to evaluate the dc current for various parameters without invoking the overdamped or the underdamped approximation. The dc current is found to have nonlinear dependence upon the damping constant, the potential parameter, and the ac force magnitude and frequency.

Folding of the GB1 hairpin peptide from discrete path sampling

The discrete path sampling technique is used to calculate folding pathways of the 16-amino acid beta hairpin-forming sequence from residues 41-56 of the B1 domain of protein G. The folding time is obtained using master equation dynamics and kinetic Monte Carlo simulations, and the time evolution of different order parameters and occupation probabilities of groups of minima are calculated and used to characterize intermediates on the folding pathway.

Rate constants for diffusive processes by partial path sampling

We introduce a path sampling method for the computation of rate constants for complex systems with a highly diffusive character. Based on the recently developed transition interface sampling (TIS) algorithm this procedure increases the efficiency by sampling only parts of complete transition trajectories. The algorithm assumes the loss of memory for diffusive progression along the reaction coordinate. We compare the new partial path technique to the TIS method for a simple diatomic system and show that the computational effort of the new method scales linearly, instead of quadratically, with the width of the diffusive barrier. The validity of the memory loss assumption is also discussed.

Transient response of a Brownian particle with general damping

We study the transient response of a Brownian particle with general damping in a system of metastable potential well. The escape rate is evaluated as a function of time after an infinite wall is removed from the potential barrier. It takes a relaxation time for the rate to reach its limit value and this rate relaxation time differs from the relaxation time of the majority of the probability around the bottom of the potential well. The rate relaxation time is found to depend on the temperature as well as the damping constant. It involves the diffusion time and the instanton time, in general agreement with recent studies of the overdamped case by Bier et al. [Phys. Rev. E 59, 6422 (1999)].

From transition paths to transition states and rate coefficients

Transition states are defined as points in configuration space with the highest probability that trajectories passing through them are reactive (i.e., form transition paths between reactants and products). In the high-friction (diffusive) limit of Langevin dynamics, the resulting ensemble of transition states is shown to coincide with the separatrix formed by points of equal commitment (or splitting) probabilities for reaching the product and reactant regions. Transition states according to the new criterion can be identified directly from equilibrium trajectories, or indirectly by calculating probability densities in the equilibrium and transition-path ensembles using umbrella and transition-path sampling, respectively. An algorithm is proposed to calculate rate coefficients from the transition-path and equilibrium ensembles by estimating the frequency of transitions between reactants and products.

Implementation of an adaptive umbrella sampling method for the calculation of multidimensional potential of mean force of chemical reactions in solution

We describe the implementation of an adaptive umbrella sampling method, making use of the weighted histogram analysis method, for computing multidimensional potential of mean force for chemical reaction in solution. The approach is illustrated by investigating the effect of aqueous solution on the free energy surface for the proton transfer reaction of [H(3)N-H-NH(3)](+) using a combined quantum mechanical and molecular mechanical AM1/TIP3P potential.

Transition-path sampling of beta-hairpin folding

We examine the dynamical folding pathways of the C-terminal beta-hairpin of protein G-B1 in explicit solvent at room temperature by means of a transition-path sampling algorithm. In agreement with previous free-energy calculations, the resulting path ensembles reveal a folding mechanism in which the hydrophobic residues collapse first followed by backbone hydrogen-bond formation, starting with the hydrogen bonds inside the hydrophobic core. In addition, the path ensembles contain information on the folding kinetics, including solvent motion. Using the recently developed transition interface sampling technique, we calculate the rate constant for unfolding of the protein fragment and find it to be in reasonable agreement with experiments. The results support the validation of using all-atom force fields to study protein folding.

Exploring potential energy surfaces for chemical reactions: an overview of some practical methods

Potential energy surfaces form a central concept in the application of electronic structure methods to the study of molecular structures, properties, and reactivities. Recent advances in tools for exploring potential energy surfaces are surveyed. Methods for geometry optimization of equilibrium structures, searching for transition states, following reaction paths and ab initio molecular dynamics are discussed. For geometry optimization, topics include methods for large molecules, QM/MM calculations, and simultaneous optimization of the wave function and the geometry. Path optimization methods and dynamics based techniques for transition state searching and reaction path following are outlined. Developments in the calculation of ab initio classical trajectories in the Born-Oppenheimer and Car-Parrinello approaches are described.

Quantum mechanical methods for enzyme kinetics

This review discusses methods for the incorporation of quantum mechanical effects into enzyme kinetics simulations in which the enzyme is an explicit part of the model. We emphasize three aspects: (a) use of quantum mechanical electronic structure methods such as molecular orbital theory and density functional theory, usually in conjunction with molecular mechanics; (b) treating vibrational motions quantum mechanically, either in an instantaneous harmonic approximation, or by path integrals, or by a three-dimensional wave function coupled to classical nuclear motion; (c) incorporation of multidimensional tunneling approximations into reaction rate calculations.

Ab initio molecular dynamics with density functional theory

Recent applications of density functional theory base ab initio molecular dynamics in chemical relevant systems are reviewed. The emphasis is on the dynamical aspect in the study of structures, reaction mechanisms, and electronic properties in both the molecular and condensed phases. Examples were chosen from fluxional molecules, solution reactions, and biological systems to illustrate the broad potential applications and unique information that can be obtained from ab initio molecular dynamics calculations. Recent advances in the development of efficient numerical algorithms for the prediction of spectroscopic properties are highlighted.

Transition path sampling: throwing ropes over rough mountain passes, in the dark

This article reviews the concepts and methods of transition path sampling. These methods allow computational studies of rare events without requiring prior knowledge of mechanisms, reaction coordinates, and transition states. Based upon a statistical mechanics of trajectory space, they provide a perspective with which time dependent phenomena, even for systems driven far from equilibrium, can be examined with the same types of importance sampling tools that in the past have been applied so successfully to static equilibrium properties.