Minimum Free Energy Path Research Articles

Calculating solution redox free energies with ab initio quantum mechanical/molecular mechanical minimum free energy path method

A quantum mechanical/molecular mechanical minimum free energy path (QM/MM-MFEP) method was developed to calculate the redox free energies of large systems in solution with greatly enhanced efficiency for conformation sampling. The QM/MM-MFEP method describes the thermodynamics of a system on the potential of mean force surface of the solute degrees of freedom. The molecular dynamics (MD) sampling is only carried out with the QM subsystem fixed. It thus avoids "on-the-fly" QM calculations and thus overcomes the high computational cost in the direct QM/MM MD sampling. In the applications to two metal complexes in aqueous solution, the new QM/MM-MFEP method yielded redox free energies in good agreement with those calculated from the direct QM/MM MD method. Two larger biologically important redox molecules, lumichrome and riboflavin, were further investigated to demonstrate the efficiency of the method. The enhanced efficiency and uncompromised accuracy are especially significant for biochemical systems. The QM/MM-MFEP method thus provides an efficient approach to free energy simulation of complex electron transfer reactions.

Mechanism of Cdc25B Phosphatase with the Small Molecule Substrate p-Nitrophenyl Phosphate from QM/MM-MFEP Calculations

Cdc25B is a dual-specificity phosphatase that catalyzes the dephosphorylation of the Cdk2/CycA protein complex. This enzyme is an important regulator of the human cell cycle and has been identified as a potential anticancer target. In general, protein tyrosine phosphatases are thought to bind the dianionic form of the phosphate and employ general acid catalysis via the Asp residue in the highly conserved WPD-loop. However, the Cdc25 phosphatases form a special subfamily based on their distinct differences from other protein tyrosine phosphatases. Although Cdc25B contains the (H/V)CX(5)R catalytic motif present in all other protein tyrosine phosphatases, it lacks an analogous catalytic acid residue. No crystallographic data currently exist for the complex of Cdc25B with Cdk2/CycA, so in addition to its natural protein substrate, experimental and theoretical studies are often carried out with small molecule substrates. In an effort to gain understanding of the dephosphorylation mechanism of Cdc25B with a commonly used small molecule substrate, we have performed simulations of the rate-limiting step of the reaction catalyzed by Cdc25B with the substrate p-nitrophenyl phosphate using the recently developed QM/MM Minimum Free Energy Path method (Hu et al. J. Chem. Phys. 2008, 034105). We have simulated the first step of the reaction with both the monoanionic and the dianionic forms of the substrate, and our calculations favor a mechanism involving the monoanionic form. Thus, Cdc25 may employ a unique dephosphorylation mechanism among protein tyrosine phosphatases, at least in the case of the small molecule substrate p-nitrophenyl phosphate.

Density functional theory of inhomogeneous liquids. III. Liquid-vapor nucleation

The process of nucleation of vapor bubbles from a superheated liquid and of liquid droplets from a supersaturated vapor is investigated using the modified-core van der Waals model density functional theory [J. F. Lutsko, J. Chem. Phys. 128, 184711 (2008)]. A novel approach is developed whereby nucleation is viewed as a transition from a metastable state to a stable state via the minimum free energy path, which is identified using the nudged elastic-band method for exploring free energy surfaces. This allows for the unbiased calculation of the properties of sub- and super-critical clusters, as well as of the critical cluster. For Lennard-Jones fluids, the results compare well to simulation and support the view that at low to moderate supersaturation nucleation proceeds smoothly rather than via spinodal-like instabilities as has been suggested recently. The possibility that the nucleation barrier vanishes at high supersaturation is, however, not ruled out.

Mechanistic Analysis of the Base-Catalyzed HF Elimination from 4-Fluoro-4-(4'-nitrophenyl)butane-2-one Based on Liquid-Phase Kinetic Isotope Effects Calculated by Dynamics Modeling with Multidimensional Tunneling.

The primary and secondary deuterium kinetic isotope effects as well as leaving-group fluorine kinetic isotope effects have been calculated for the base-promoted elimination of hydrogen fluoride from 4-fluoro-4-(4'-nitrophenyl)butane-2-one in 75% aqueous methanol solution. The elimination was studied for both formate and imidazole as the catalytic base; and reactant and transition state structures and vibrational frequencies have been calculated by including the base explicitly and by including the solvent by an implicit solvation model that includes both electrostatics by class IV charges and first-solvation-shell effects by atomic surface tensions. We used the M06-L density functional for all calculations. The optimized stationary points, the geometry changes along the solution-phase minimum free energy path, and the solution-phase free energy profile indicate that the elimination reaction occurs concertedly but asynchronously via an E1cb-like transition state. Reaction rates were calculated by the equilibrium solvation path method, using variational transition state theory with multidimensional tunneling. The primary deuterium kinetic isotope effects are calculated to be large: 1.67 and 5.13 for formate and imidazole, respectively. The corresponding C4-secondary deuterium kinetic isotope effects are 1.044 and 1.044, and the leaving group fluorine kinetic isotope effects are respectively 1.020 and 1.015.

Mechanism of OMP decarboxylation in orotidine 5'-monophosphate decarboxylase.

Despite extensive experimental and theoretical studies, the detailed catalytic mechanism of orotidine 5'-monophosphate decarboxylase (ODCase) remains controversial. In particular simulation studies using high level quantum mechanics have failed to reproduce experimental activation free energy. One common feature of many previous simulations is that there is a water molecule in the vicinity of the leaving CO2 group whose presence was only observed in the inhibitor bound complex of ODCase/BMP. Various roles have even been proposed for this water molecule from the perspective of stabilizing the transition state and/or intermediate state. We hypothesize that this water molecule is not present in the active ODCase/OMP complex. Based on QM/MM minimum free energy path simulations with accurate density functional methods, we show here that in the absence of this water molecule the enzyme functions through a simple direct decarboxylation mechanism. Analysis of the interactions in the active site indicates multiple factors contributing to the catalysis, including the fine-tuned electrostatic environment of the active site and multiple hydrogen-bonding interactions. To understand better the interactions between the enzyme and the inhibitor BMP molecule, simulations were also carried out to determine the binding free energy of this special water molecule in the ODCase/BMP complex. The results indicate that the water molecule in the active site plays a significant role in the binding of BMP by contributing approximately -3 kcal/mol to the binding free energy of the complex. Therefore, the complex of BMP plus a water molecule, instead of the BMP molecule alone, better represents the tight binding transition state analogue of ODCase. Our simulation results support the direct decarboxylation mechanism and highlight the importance of proper recognition of protein bound water molecules in the protein-ligand binding and the enzyme catalysis.

Kinetics of phase transitions in two dimensional Ising models studied with the string method

The kinetics of phase transitions in the two dimensional Ising model under different conditions is studied using the string method. The key idea is to work in collective variables, consisting of block of spins, which allow for a continuous approximation of the collective variable state-space. The string method computes the minimum free energy path (MFEP) in this collective variable space, which is shown to explain the mechanism of the phase transformation (in particular, an approximation of its committor function, its free energy and its transition state). In this paper the theoretical background of the technique as well as its computational aspects are discussed in details. The string method is then used to analyze phase transition in the Ising model with imposed boundary conditions and in a periodic system under an external field of increasing magnitude. In each case, the mechanism of the phase transformation is elucidated.

Quantum mechanics/molecular mechanics minimum free-energy path for accurate reaction energetics in solution and enzymes: Sequential sampling and optimization on the potential of mean force surface

To accurately determine the reaction path and its energetics for enzymatic and solution-phase reactions, we present a sequential sampling and optimization approach that greatly enhances the efficiency of the ab initio quantum mechanics/molecular mechanics minimum free-energy path (QM/MM-MFEP) method. In the QM/MM-MFEP method, the thermodynamics of a complex reaction system is described by the potential of mean force (PMF) surface of the quantum mechanical (QM) subsystem with a small number of degrees of freedom, somewhat like describing a reaction process in the gas phase. The main computational cost of the QM/MM-MFEP method comes from the statistical sampling of conformations of the molecular mechanical (MM) subsystem required for the calculation of the QM PMF and its gradient. In our new sequential sampling and optimization approach, we aim to reduce the amount of MM sampling while still retaining the accuracy of the results by first carrying out MM phase-space sampling and then optimizing the QM subsystem in the fixed-size ensemble of MM conformations. The resulting QM optimized structures are then used to obtain more accurate sampling of the MM subsystem. This process of sequential MM sampling and QM optimization is iterated until convergence. The use of a fixed-size, finite MM conformational ensemble enables the precise evaluation of the QM potential of mean force and its gradient within the ensemble, thus circumventing the challenges associated with statistical averaging and significantly speeding up the convergence of the optimization process. To further improve the accuracy of the QM/MM-MFEP method, the reaction path potential method developed by Lu and Yang [Z. Lu and W. Yang, J. Chem. Phys. 121, 89 (2004)] is employed to describe the QM/MM electrostatic interactions in an approximate yet accurate way with a computational cost that is comparable to classical MM simulations. The new method was successfully applied to two example reaction processes, the classical SN2 reaction of Cl-+CH3Cl in solution and the second proton transfer step of the reaction catalyzed by the enzyme 4-oxalocrotonate tautomerase. The activation free energies calculated with this new sequential sampling and optimization approach to the QM/MM-MFEP method agree well with results from other simulation approaches such as the umbrella sampling technique with direct QM/MM dynamics sampling, demonstrating the accuracy of the iterative QM/MM-MFEP method.

Coarse-Grained Free Energy Functions for Studying Protein Conformational Changes: A Double-Well Network Model

In this work, a double-well network model (DWNM) is presented for generating a coarse-grained free energy function that can be used to study the transition between reference conformational states of a protein molecule. Compared to earlier work that uses a single, multidimensional double-well potential to connect two conformational states, the DWNM uses a set of interconnected double-well potentials for this purpose. The DWNM free energy function has multiple intermediate states and saddle points, and is hence a “rough” free energy landscape. In this implementation of the DWNM, the free energy function is reduced to an elastic-network model representation near the two reference states. The effects of free energy function roughness on the reaction pathways of protein conformational change is demonstrated by applying the DWNM to the conformational changes of two protein systems: the coil-to-helix transition of the DB-loop in G-actin and the open-to-closed transition of adenylate kinase. In both systems, the rough free energy function of the DWNM leads to the identification of distinct minimum free energy paths connecting two conformational states. These results indicate that while the elastic-network model captures the low-frequency vibrational motions of a protein, the roughness in the free energy function introduced by the DWNM can be used to characterize the transition mechanism between protein conformations.

Generalized gradient-augmented harmonic Fourier beads method with multiple atomic and/or center-of-mass positional restraints

We describe a generalization of the gradient-augmented harmonic Fourier beads method for finding minimum free-energy transition path ensembles and similarly minimum potential energy paths to allow positional restraints on the centers of mass of selected atoms. The generalized gradient-augmented harmonic Fourier beads (ggaHFB) method further extends the scope of the HFB methodology to studying molecule transport across various mobile phases such as lipid membranes. Furthermore, the new implementation improves the applicability of the HFB method to studies of ligand binding, protein folding, and enzyme catalysis as well as modeling equilibrium pulling experiments. Like its predecessor, the ggaHFB method provides accurate energy profiles along the specified paths and in certain simple cases avoids the need for path optimization. The utility of the ggaHFB method is demonstrated with an application to the water permeation through a single-wall (5,5) carbon nanotube with a diameter of 6.78 A and length of 16.0 A. We provide a simple rationale as to why water enters the hydrophobic nanotube and why it does so in pulses and in wire assembly.

On-the-fly string method for minimum free energy paths calculation

An improved and simplified version of the string method in collective variables for computing minimum free energy paths is proposed. The string is discretized into a finite number of images and in the new method these images are evolved concurrently with replicas of the original system while keeping their inter-distances equal via reparametrization. There is no need to compute the mean force by time-averaging, nor to project the force perpendicularly to the string. In this Letter, the algorithmic aspects of the on-the-fly string method are presented in detail and the method is tested on the solvated alanine dipeptide molecule.

QM/MM Minimum Free Energy Path: Methodology and Application to Triosephosphate Isomerase.

Structural and energetic changes are two important characteristic properties of a chemical reaction process. In the condensed phase, studying these two properties is very challenging because of the great computational cost associated with the quantum mechanical calculations and phase space sampling. Although the combined quantum mechanics/molecular mechanics (QM/MM) approach significantly reduces the amount of the quantum mechanical calculations and facilitates the simulation of solution phase and enzyme catalyzed reactions, the required quantum mechanical calculations remain quite expensive and extensive sampling can be achieved routinely only with semiempirical quantum mechanical methods. QM/MM simulations with ab initio QM methods, therefore, are often restricted to narrow regions of the potential energy surface such as the reactant, product and transition state, or the minimum energy path. Such ab initio QM/MM calculations have previously been performed with the QM/MM-Free Energy (QM/MM-FE) method of Zhang et al.1 to generate the free energy profile along the reaction coordinate using free energy perturbation calculations at fixed structures of the QM subsystems. Results obtained with the QM/MM-FE method depend on the determination of the minimum energy reaction path, which is based on local conformations of the protein/solvent environment and can be difficult to obtain in practice. To overcome the difficulties associated with the QM/MM-FE method and to further enhance the sampling of the MM environment conformations, we develop here a new method to determine the QM/MM minimum free energy path (QM/MM-MFEP) for chemical reaction processes in solution and in enzymes. Within the QM/MM framework, we express the free energy of the system as a function of the QM conformation, thus leading to a simplified potential of mean force (PMF) description for the thermodynamics of the system. The free energy difference between two QM conformations is evaluated by the QM/MM free energy perturbation method. The free energy gradients with respect to the QM degrees of freedom are calculated from molecular dynamics simulations at given QM conformations. With the free energy and free energy gradients in hand, we further implement chain-of-conformation optimization algorithms in the search for the reaction path on the free energy surface without specifying a reaction coordinate. This method thus efficiently provides a unique minimum free energy path for solution and enzyme reactions, with structural and energetic properties being determined simultaneously. To further incorporate the dynamic contributions of the QM subsystem into the simulations, we develop the reaction path potential of Lu, et al.2 for the minimum free energy path. The combination of the methods developed here presents a comprehensive and accurate treatment for the simulation of reaction processes in solution and in enzymes with ab initio QM/MM methods. The method has been demonstrated on the first step of the reaction of the enzyme triosephosphate isomerase with good agreement with previous studies.

Harmonic Fourier beads method for studying rare events on rugged energy surfaces

We present a robust, distributable method for computing minimum free energy paths of large molecular systems with rugged energy landscapes. The method, which we call harmonic Fourier beads (HFB), exploits the Fourier representation of a path in an appropriate coordinate space and proceeds iteratively by evolving a discrete set of harmonically restrained path points-beads-to generate positions for the next path. The HFB method does not require explicit knowledge of the free energy to locate the path. To compute the free energy profile along the final path we employ an umbrella sampling method in two generalized dimensions. The proposed HFB method is anticipated to aid the study of rare events in biomolecular systems. Its utility is demonstrated with an application to conformational isomerization of the alanine dipeptide in gas phase.

String method in collective variables: Minimum free energy paths and isocommittor surfaces

A computational technique is proposed which combines the string method with a sampling technique to determine minimum free energy paths. The technique only requires to compute the mean force and another conditional expectation locally along the string, and therefore can be applied even if the number of collective variables kept in the free energy calculation is large. This is in contrast with other free energy sampling techniques which aim at mapping the full free energy landscape and whose cost increases exponentially with the number of collective variables kept in the free energy. Provided that the number of collective variables is large enough, the new technique captures the mechanism of transition in that it allows to determine the committor function for the reaction and, in particular, the transition state region. The new technique is illustrated on the example of alanine dipeptide, in which we compute the minimum free energy path for the isomerization transition using either two or four dihedral angles as collective variables. It is shown that the mechanism of transition can be captured using the four dihedral angles, but it cannot be captured using only two of them.

Tracing the minimum-energy path on the free-energy surface

The free-energy profile of a reaction can be estimated in a molecular-dynamics approach by imposing a mechanical constraint along a reaction coordinate (RC). Many recent studies have shown that the temperature can greatly influence the path followed by the reactants. Here, we propose a practical way to construct the minimum-energy path directly on the free-energy surface at a given temperature. First, we follow the blue-moon ensemble method to derive the expression of the free-energy gradient for a given RC. These derivatives are then used to find the actual minimum-energy reaction path at finite temperature, in a way similar to the intrinsic reaction path of Fukui on the potential-energy surface. [K. Fukui, J. Phys. Chem. 74, 4161 (1970)]. Once the path is known, one can calculate the free-energy profile using thermodynamic integration. We also show that the mass-metric correction cancels for many types of constraints, making the procedure easy to use. Finally, the minimum-free-energy path at 300 K for the addition of CCl2 to ethylene is compared with a path based on a simple one-dimensional reaction coordinate. A comparison is also given with the reaction path at 0 K.

Angular distribution functions and specific local structures in liquid water

A method of structural analysis, based on restricted averages of the angular pair correlation function, is proposed and applied to liquid water. Minimum free energy paths for a molecule moving in an average local structure are obtained. Around a central molecule, accepting and donating neighbors form hydrogen bonds having the same energy minima but differing for the availability of low free energy states. Some bounded states have specific bridging functions and connect hydrogen-bond and transition states along negative free energy pathways, thus favoring local diffusion. Bridging and transition states play the role of structural defects and include bifurcated hydrogen-bond and interstitial water geometries. Our methods and results can be relevant to understand, from a structural point of view, the anomalous behavior of stressed and confined water.

On Stochastic Processes Defined by Differential Equations with a Small Parameter

On Stochastic Processes Defined by Differential Equations with a Small Parameter