Permeation and gating in proteins: kinetic Monte Carlo reaction path following

We present a new Monte Carlo technique, kinetic Monte Carlo reaction path following (kMCRPF), for the computer simulation of permeation and large-scale gating transitions in protein channels. It combines ideas from Metropolis Monte Carlo (MMC) and kinetic Monte Carlo (kMC) algorithms, and is particularly suitable when a reaction coordinate is well defined. Evolution of transition proceeds on the reaction coordinate by small jumps (kMC technique) toward the nearest lowest-energy uphill or downhill states, with the jumps thermally activated (constrained MMC). This approach permits navigation among potential minima on an energy surface, finding the minimum-energy paths and determining their associated free-energy profiles. The methodological and algorithmic strategies underlying the kMCRPF method are described. We have tested it using an analytical model and applied it to study permeation through the curvilinear ClC chloride and aquaporin pores and to gating in the gramicidin A channel. These studies of permeation and gating in real proteins provide extensive procedural tests of the method.     

An algorithm to find minimum free-energy paths using umbrella integration

The calculation of free-energy barriers by umbrella sampling and many other methods is hampered by the necessity for an a priori choice of the reaction coordinate along which to sample. We avoid this problem by providing a method to search for saddle points on the free-energy surface in many coordinates. The necessary gradients and Hessians of the free energy are obtained by multidimensional umbrella integration. We construct the minimum free-energy path by following the gradient down to minima on the free-energy surface. The change of free energy along the path is obtained by integrating out all coordinates orthogonal to the path. While we expect the method to be applicable to large systems, we test it on the alanine dipeptide in vacuum. The minima, transition states, and free-energy barriers agree well with those obtained previously with other methods.     

The reaction path intrinsic reaction coordinate method and the Hamilton-Jacobi theory

The definition and location of an intrinsic reaction coordinate path is of crucial importance in many areas of theoretical chemistry. Differential equations used to define the path hitherto are complemented in this study with a variational principle of Fermat type, as Fukui [Int. J. Quantum Chem., Quantum Chem. Symp. 15, 633 (1981)] reported in a more general form some time ago. This definition is more suitable for problems where initial and final points are given. The variational definition can naturally be recast into a Hamilton-Jacobi equation. The character of the variational solution is studied via the Weierstrass necessary and sufficient conditions. The characterization of the local minima character of the intrinsic reaction coordinate is proved. Such result leads to a numerical algorithm to find intrinsic reaction coordinate paths based on the successive minimizations of the Weierstrass E-function evaluated on a guess curve connecting the initial and final points of the desired path.     

One-dimensional free-energy profiles of complex systems: progress variables that preserve the barriers

We show that the balanced minimum-cut procedure introduced in PNAS 2004, 101, 14766 can be reinterpreted as a method for solving the constrained optimization problem of finding the minimum cut among the cuts with a particular value of an additive function of the nodes on either side of the cut. Such an additive function (e.g., the partition function of the reactant region) can be used as a progress coordinate to determine a one-dimensional profile (FEP) of the free-energy surface of the protein-folding reaction as well as other complex reactions. The algorithm is based on the network (obtained from an equilibrium molecular dynamics simulation) that represents the calculated reaction behavior. The resulting FEP gives the exact values of the free energy as a function of the progress coordinate; i.e., at each value of the progress coordinate, the profile is obtained from the surface with the minimal partition function among the surfaces that divide the full free-energy surface between two chosen end points. In many cases, the balanced minimum-cut procedure gives results for only a limited set of points. An approximate method based on p(fold) is shown to provide the profile for a more complete set of values of the progress coordinate. Applications of the approach to model problems and to realistic systems (beta-hairpin of protein G, LJ38 cluster) are presented.     

A recipe for the computation of the free energy barrier and the lowest free energy path of concerted reactions

The recently introduced hills method (Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562) is a powerful tool to compute the multidimensional free energy surface of intrinsically concerted reactions. We have extended this method by focusing our attention on localizing the lowest free energy path that connects the stable reactant and product states. This path represents the most probable reaction mechanism, similar to the zero temperature intrinsic reaction coordinate, but also includes finite temperature effects. The transformation of the multidimensional problem to a one-dimensional reaction coordinate allows for accurate convergence of the free energy profile along the lowest free energy path using standard free energy methods. Here we apply the hills method, our lowest free energy path search algorithm, and umbrella sampling to the prototype S(N)2 reaction. The hills method replaces the in many cases difficult problem of finding a good reaction coordinate with choosing relatively simple collective variables, such as the bond lengths of the broken and formed chemical bonds. The second part of the paper presents a guide to using the hills method, in which we test and fine-tune the method for optimal accuracy and efficiency using the umbrella sampling results as a reference.     

RISM-SCF study of the free-energy profile of the Menshutkin-type reaction NH3+CH3Cl→NH3CH3++Cl− in aqueous solution

The free-energy profile for the Menshutkin-type reaction NH₃ + CH₃Cl → NH₃CH₃ ⁺ + Cl⁻ in aqueous solution is studied using the RISM-SCF method. The effect of electron correlation on the free-energy profile is estimated by the RISM-MP2 method at the HF optimized geometries along the reaction coordinate. Solvation was found to have a large influence on the vibrational frequencies at the reactant, transition state and product; these vibrational frequencies are utilized to calculate the zero-point energy correction of the free-energy profile. The computed barrier height and reaction exothermicity are in reasonable agreement with those of experiment and previous calculations. The change of solvation structure along the reaction path is represented by radial distribution functions between solute-solvent atomic sites. The mechanisms of the reaction are discussed from the view points of solute electronic and solvation structures. Received: 26 June 1998/Accepted: 28 August 1998 / Published online: 2 November 1998     

Two-dimensional free-energy surface on the exchange reaction of alkyl chloride/chloride using the QM/MM-MC method

Two-dimensional free-energy surfaces are calculated for alkyl chloride/chloride exchange/inversion reactions: Cl(-)+RCl (R=Me and t-Bu) surrounded by one hundred H(2)O molecules as a model of solvent. The methodology of free-energy calculation by perturbation theory based on a mixed-Hamiltonian model (QM/MM) combined with Monte Carlo sampling of the solvent configurations was used to obtain the changes in solvation free energy. We devised a special procedure to analyze the two-dimensional free-energy surfaces to gain unique insight into the differences in the reaction mechanisms between the two systems. The inversion reaction path for R=t-Bu on the free-energy surface is found to proceed in an asynchronous way within a concerted framework via the ion-pair region. This is in contrast to the R=Me system that proceeds as a typical S(N)2 reaction.     

An internal coordinate invariant reaction pathway

In this work we show that some properties of a potential energy surface are not independent of the choice of the coordinate frame. So the reaction pathway often described as steepest descent way does not correspond to an invariant curve under coordinate transformations. We propose an internal intrinsic reaction pathway by using some quasi-dynamical considerations (like instantaneous internal acceleration). Our work precises the intrinsic-reaction coordinates of Fukui to any set of 3N-6 internal parameters. Finally, from the equations of motion we deduce the form of the normal reaction coordinates frame anywhere along the postulated reaction pathway.     

Thermodynamics of catalytic formation of dimethyl ether from methanol in acidic zeolites

We present a theoretical study of the formation of the first intermediate, dimethyl ether, in the methanol to gasoline conversion within the framework of an ab initio molecular dynamics approach. The study is performed under conditions that closely resemble the reaction conditions in the zeolite catalyst including the full topology of the framework. The use of the method of thermodynamic integration allows us to extract the free-energy profile along the reaction coordinate. We find that the entropic contribution qualitatively alters the free-energy profile relative to the total energy profile. Different transition states are found from the internal and free energy profiles. The entropy contribution varies significantly along the reaction coordinate and is responsible for stabilizing the products and for lowering the energy barrier. The hugely inhomogeneous variation of the entropy can be understood in terms of elementary processes that take place during the chemical reaction. Our simulations provide new insights into the complex nature of this chemical reaction.     

On coordinate transformations in steepest descent path and stationary point locations

Two aspects of the problems of calculating steepest descent paths and locating stationary points on surfaces E( X ), which are sources of some confusion in the literature, are addressed. These include writing proper expressions for the gradient and Hessian, and their transformation properties relative to coordinate transformations, based on the invariance of the surface E( X ). The appropriate transformation is derived, based on a constrained energy minimization condition, to achieve what we call the Hessian eigenvalue representation. This not only allows decoupling of the variables, but also points to the minimization direction and preserves the eigenvalues of the Hessian. These results allow one to use the steepest descent path and stationary point location algorithms in any coordinate system and obtain invariant results. The validity of these considerations are also confirmed through numerical examples. The stationary condition with constrained kinematic path length is also shown to yield a Hessian eigenvalue representation for the normal modes for small vibrations. Lastly, we have constructed a mathematically consistent definition of mass-weighted Cartesians where the intrinsic reaction path of Fukui is a steepest descent path. © 1992 John Wiley & Sons, Inc.     

Topography of the free-energy landscape probed via mechanical unfolding of proteins

Single-molecule experiments in which proteins are unfolded by applying mechanical stretching forces generally force unfolding to proceed along a reaction coordinate that is different from that in chemical or thermal denaturation. Here we simulate the mechanical unfolding and refolding of a minimalist off-lattice model of the protein ubiquitin to explore in detail the slice of the multidimensional free-energy landscape that is accessible via mechanical pulling experiments. We find that while the free-energy profile along typical "chemical" reaction coordinates may exhibit two minima, corresponding to the native and denatured states, the free energy G(z) is typically a monotonic function of the mechanical coordinate z equal to the protein extension. Application of a stretching force along z tilts the free-energy landscape resulting in a bistable (or multistable) free energy G(z)-fz probed in mechanical unfolding experiments. We construct a two-dimensional free-energy surface as a function of both chemical and mechanical reaction coordinates and examine the coupling between the two. We further study the refolding trajectories after the protein has been prestretched by a large force, as well as the mechanical unfolding trajectories in the presence of a large stretching force. We demonstrate that the stretching forces required to destabilize the native state thermodynamically are larger than those expected on the basis of previous experimental estimates of G(z). This finding is consistent with the recent experimental studies, indicating that proteins may refold even in the presence of a substantial stretching force. Finally, we show that for certain temperatures the free energy of a polyprotein chain consisting of multiple domains is a linear function of the chain extension. We propose that the recently observed "slow phase" in the refolding of proteins under mechanical tension may be viewed as downhill diffusion in such a linear potential.     

Theoretical study on the alkaline hydrolysis of methyl thioacetate in aqueous solution

A base-catalyzed hydrolysis reaction of thiolester has been studied in both gas and solution phases using two ab initio quantum mechanics calculations such as Gaussian09 and CPMD. The free-energy surface along the reaction path is also constructed using a configuration sampling technique, namely, the metadynamics method. While there are two different reaction paths obtained for the potential profile of the base-catalyzed hydrolysis reaction for thiolester in the gas phase, a triple-well reaction path is computed for the reaction in the solution phase by two quantum mechanics calculations. Unlike the S(N)2 mechanism (a concerted mechanism) found for the gas-phase reaction, a nucleophilic attack from the hydroxide ion on the carbonyl carbon to yield a tetrahedral intermediate (a stepwise mechanism) is observed for the solution-phase reaction. Moreover, the energy profiles computed by these two theoretical calculations are found to be very comparable with those determined experimentally.     

气态烯丙醇光异构化反应的理论研究

本文采用自洽场分子轨道UHF/6-31G从头计算法, 辅以能量梯度法研究气相烯丙醇光异构化反应的机理。全部优化了T1态势能面上反应物、过渡态、中间体和产物的几何构型。基于Fukui提出的内禀反应坐标理论(IRC)计算这一体系的反应途径, 并针对各驻点进行MP2/6-31G的相关能校正, 得到该反应在激发态进行并为一经历双自由基中间体的分步反应的结论。支持了实验工作者提出的机理。     

Better understanding of the ring-cleavage process of cyanocyclopropyl anionic derivatives. A theoretical study based on the electron localization function

[reaction: see text] Theoretical calculations at the B3LYP/6-31+G(d), MP2/6-31+G(d), and G3(MP2) levels have been carried out to understand the alternative reaction pathways (the cyclopropyl ring cleavage (RC) and the retrocycloaddition reaction (rCA)) of a constrained tricyanocyclopropyl anionic derivative. The more energetically favorable path is found to be the RC process, a formally "forbidden" rearrangement (Leiviers, M.; Tam, I.; Groves, K.; Leung, D.; Xie, Y.; Breslow, R. Org. Lett. 2003, 5, 19, 3407) yielding an allylic anion system via a concerted transition structure, in agreement with experimental outcomes. rCA is more energetically favorable along a two-stage mechanism, via an intermediate, than a synchronous concerted process. By using isodesmic reactions, we have found that B3LYP presents limitations when it calculates carbon-carbon bond-breaking processes along the present rCA reaction. A detailed analysis of the nature of the topology of the reactive potential energy surface for the RC process points out the presence of a valley-ridge inflection point in the uphill part. An explanation for the low-energy barrier associated with RC is furnished on the analysis of the evolution of the twisting (dis-/conrotatory) motions of cyano substituents in the cyclopropyl ring as well as on the number and type of electron pairs provided by the electron localization function (ELF).     

Quadratic string method for determining the minimum-energy path based on multiobjective optimization

Based on a multiobjective optimization framework, we develop a new quadratic string method for finding the minimum-energy path. In the method, each point on the minimum-energy path is minimized by integration in the descent direction perpendicular to path. Each local integration is done on a quadratic surface approximated by a damped Broyden-Fletcher-Goldfarb-Shanno updated Hessian, allowing the algorithm to take many steps between energy and gradient calls. The integration is performed with an adaptive step-size solver, which is restricted in length to the trust radius of the approximate Hessian. The full algorithm is shown to be capable of practical superlinear convergence, in contrast to the linear convergence of other methods. The method also eliminates the need for predetermining such parameters as step size and spring constants, and is applicable to reactions with multiple barriers. The effectiveness of this method is demonstrated for the Muller-Brown potential, a seven-atom Lennard-Jones cluster, and the enolation of acetaldehyde to vinyl alcohol.     

Obtaining reaction coordinates by likelihood maximization

We present a new approach for calculating reaction coordinates in complex systems. The new method is based on transition path sampling and likelihood maximization. It requires fewer trajectories than a single iteration of existing procedures, and it applies to both low and high friction dynamics. The new method screens a set of candidate collective variables for a good reaction coordinate that depends on a few relevant variables. The Bayesian information criterion determines whether additional variables significantly improve the reaction coordinate. Additionally, we present an advantageous transition path sampling algorithm and an algorithm to generate the most likely transition path in the space of collective variables. The method is demonstrated on two systems: a bistable model potential energy surface and nucleation in the Ising model. For the Ising model of nucleation, we quantify for the first time the role of nuclei surface area in the nucleation reaction coordinate. Surprisingly, increased surface area increases the stability of nuclei in two dimensions but decreases nuclei stability in three dimensions.     

A comparison of methods to compute the potential of mean force.

Most processes occurring in a system are determined by the relative free energy between two or more states because the free energy is a measure of the probability of finding the system in a given state. When the two states of interest are connected by a pathway, usually called reaction coordinate, along which the free-energy profile is determined, this profile or potential of mean force (PMF) will also yield the relative free energy of the two states. Twelve different methods to compute a PMF are reviewed and compared, with regard to their precision, for a system consisting of a pair of methane molecules in aqueous solution. We analyze all combinations of the type of sampling (unbiased, umbrella-biased or constraint-biased), how to compute free energies (from density of states or force averaging) and the type of coordinate system (internal or Cartesian) used for the PMF degree of freedom. The method of choice is constraint-bias simulation combined with force averaging for either an internal or a Cartesian PMF degree of freedom.     

Hydration properties of magnesium and calcium ions from constrained first principles molecular dynamics

We studied the solvation structures of the divalent metal cations Mg(2+) and Ca(2+) in ambient water by applying a Car-Parrinello-based constrained molecular dynamics method. By employing the metal-water oxygen coordination number as a reaction coordinate, we could identify distinct aqua complexes characterized by structural variations of the first coordination shell. In particular, our estimated free-energy profile clearly shows that the global minimum for Mg(2+) is represented by a rather stable sixfold coordination in the octahedral arrangement, in agreement with experiments. Conversely, for Ca(2+) the free-energy curve shows several shallow local minima, suggesting that the hydration structure of Ca(2+) is highly variable. Implications for water exchange reactions are also discussed.     

Exploring chemical reactivity of complex systems with path‐based coordinates: Role of the distance metric

Path‐based reaction coordinates constitute a valuable tool for free‐energy calculations in complex processes. When a reference path is defined by means of collective variables, a nonconstant distance metric that incorporates the nonorthonormality of these variables should be taken into account. In this work, we show that, accounting for the correct metric tensor, these kind of variables can provide iso‐hypersurfaces that coincide with the iso‐committor surfaces and that activation free energies equal the value that would be obtained if the committor function itself were used as reaction coordinate. The advantages of the incorporation of the variable metric tensor are illustrated with the analysis of the enzymatic reaction catalyzed by isochorismate‐pyruvate lyase. Hybrid QM/MM techniques are used to obtain the free energy profile and to analyze reactive trajectories initiated at the transition state. For this example, the committor histogram is peaked at 0.5 only when a variable metric tensor is incorporated in the definition of the path‐based coordinate. © 2014 Wiley Periodicals, Inc.