Efficient exploration of reactive potential energy surfaces using Car-Parrinello molecular dynamics

Finite-temperature effects in enzymatic reactions — Insights from QM/MM free-energy simulations

We report potential-energy and free-energy data for three enzymatic reactions: carbon–halogen bond formation in fluorinase, hydrogen abstraction from camphor in cytochrome P450cam, and chorismate-to-prephenate Claisen rearrangement in chorismate mutase. The results were obtained by combined quantum mechanics/molecular mechanics (QM/MM) optimizations and two types of QM/MM free-energy simulations (free-energy perturbation and umbrella sampling) using semi-empirical or density-functional QM methods. Based on these results and our previously published free-energy data on electrophilic substitution in para-hydroxybenzoate hydroxylase, we discuss the importance of finite-temperature effects in the chemical step of enzyme reactions. We find that the entropic contribution to the activation barrier is generally rather small, usually of the order of 5 kJ mol–1 or less, consistent with the notion that enzymes bind and pre-organize the reactants in the active site. A somewhat larger entropic contribution is encountered in the case of chorismate mutase where the pericyclic transition state is intrinsically more rigid than the chorismate reactant (also in the enzyme). The present results suggest that barriers from QM/MM geometry optimization may often be close to free-energy barriers for the chemical step in enzymatic reactions.

Some recent techniques for free energy calculations

A few recent techniques to calculate free energies in the context of molecular dynamics simulations are discussed: temperature-accelerated molecular dynamics, which is a method to explore fast the important regions in the free energy landscape associated with a set of continuous collective variables without having to know where these regions are beforehand; the single sweep method, which is a variational method to interpolate the free energy globally given a set of mean forces (i.e., a set of gradients of the free energy) calculated at specific points, or centers, on the free energy landscape; and a Voronoi-based free energy method for the calculation of the free energy of the Voronoi tessellation associated with a set of centers. We also discuss how this last technique can be used in conjunction with the string method, and how kinetic information such as reaction rates can be calculated by milestoning using the edges of a Voronoi tessellation as milestones.

Conformational free energy surface of alpha-N-acetylneuraminic acid: an interplay between hydrogen bonding and solvation

The conformational free energy surface of alpha-N-acetylneuraminic acid (Neu5Ac, sialic acid) in the space of ring-puckering coordinates was calculated using the metadynamics method. Free energy surfaces in vacuum and with an explicit solvent were calculated in GLYCAM 06 force field. In vacuum three structures are almost equivalently populated, namely, the (2)C(5) chair and the B(3,6)/(2)S(6) and (O)S(3) boat/skew-boat conformations. The B(3,6)/(2)S(6) structure is stabilized by an ionic hydrogen bond between the amide N-H bond and the carboxylic group. However, this structure is unfavorable in a water environment in which the experimentally observed (2)C(5) chair conformation is predicted to be more stable than the other structures. These results indicate that environment significantly influences conformation of Neu5Ac and that Neu5Ac-processing enzymes might modify a conformation of their substrates solely by a changing polarity of the environment. The structure of Neu5Ac bound in influenza neuraminidase ((4)S(2)/B(2,5)) belongs to conformations preferred in a water environment. The free energy penalty of this conformational change was calculated (relative to (2)C(5)) as 10.2 +/- 2.0 and 17.3 +/- 2.0 kJ/mol for (4,O)B/(O)S(3) and (4)S(2), respectively. This result indicates that mimicking of the enzyme-bound conformation is likely to be a viable strategy for the design of neuraminidase inhibitors.

Simulation of structural phase transitions by metadynamics

We describe here in detail the recently introduced methodology for simulation of structural transitions in crystals. The applications of the new scheme are illus trated on various kinds of crystals and the advantages with respect to previous schemes are emphasized. The relevance of the new method for the problem of crystal structure prediction is also discussed.

CPMD: Car-Parrinello molecular dynamics

We present the ab initio molecular dynamics code CPMD. This plane wave/pseudopotential based Kohn-Sham density functional code with a rich set of features is successfully applied to calculate static and dynamic properties for a wide range of molecules and materials. Its flexibility and high performance on many computer platforms make it an optimal tool for the study of liquids, surfaces, crystals as well as biomolecules.

QM/MM methods for biomolecular systems

Combined quantum-mechanics/molecular-mechanics (QM/MM) approaches have become the method of choice for modeling reactions in biomolecular systems. Quantum-mechanical (QM) methods are required for describing chemical reactions and other electronic processes, such as charge transfer or electronic excitation. However, QM methods are restricted to systems of up to a few hundred atoms. However, the size and conformational complexity of biopolymers calls for methods capable of treating up to several 100,000 atoms and allowing for simulations over time scales of tens of nanoseconds. This is achieved by highly efficient, force-field-based molecular mechanics (MM) methods. Thus to model large biomolecules the logical approach is to combine the two techniques and to use a QM method for the chemically active region (e.g., substrates and co-factors in an enzymatic reaction) and an MM treatment for the surroundings (e.g., protein and solvent). The resulting schemes are commonly referred to as combined or hybrid QM/MM methods. They enable the modeling of reactive biomolecular systems at a reasonable computational effort while providing the necessary accuracy.

Single-sweep methods for free energy calculations

A simple, efficient, and accurate method is proposed to map multidimensional free energy landscapes. The method combines the temperature-accelerated molecular dynamics (TAMD) proposed in [L. Maragliano and E. Vanden-Eijnden, Chem. Phys. Lett. 426, 168 (2006)] with a variational reconstruction method using radial-basis functions for the representation of the free energy. TAMD is used to rapidly sweep through the important regions of the free energy landscape and to compute the gradient of the free energy locally at points in these regions. The variational method is then used to reconstruct the free energy globally from the mean force at these points. The algorithmic aspects of the single-sweep method are explained in detail, and the method is tested on simple examples and used to compute the free energy of the solvated alanine dipeptide in two and four dihedral angles.

Adaptive biasing force method for scalar and vector free energy calculations

In free energy calculations based on thermodynamic integration, it is necessary to compute the derivatives of the free energy as a function of one (scalar case) or several (vector case) order parameters. We derive in a compact way a general formulation for evaluating these derivatives as the average of a mean force acting on the order parameters, which involves first derivatives with respect to both Cartesian coordinates and time. This is in contrast with the previously derived formulas, which require first and second derivatives of the order parameter with respect to Cartesian coordinates. As illustrated in a concrete example, the main advantage of this new formulation is the simplicity of its use, especially for complicated order parameters. It is also straightforward to implement in a molecular dynamics code, as can be seen from the pseudocode given at the end. We further discuss how the approach based on time derivatives can be combined with the adaptive biasing force method, an enhanced sampling technique that rapidly yields uniform sampling of the order parameters, and by doing so greatly improves the efficiency of free energy calculations. Using the backbone dihedral angles Phi and Psi in N-acetylalanyl-N'-methylamide as a numerical example, we present a technique to reconstruct the free energy from its derivatives, a calculation that presents some difficulties in the vector case because of the statistical errors affecting the derivatives.

Adaptively biased molecular dynamics for free energy calculations

We present an adaptively biased molecular dynamics (ABMD) method for the computation of the free energy surface of a reaction coordinate using nonequilibrium dynamics. The ABMD method belongs to the general category of umbrella sampling methods with an evolving biasing potential and is inspired by the metadynamics method. The ABMD method has several useful features, including a small number of control parameters and an O(t) numerical cost with molecular dynamics time t. The ABMD method naturally allows for extensions based on multiple walkers and replica exchange, where different replicas can have different temperatures and/or collective variables. This is beneficial not only in terms of the speed and accuracy of a calculation, but also in terms of the amount of useful information that may be obtained from a given simulation. The workings of the ABMD method are illustrated via a study of the folding of the Ace-GGPGGG-Nme peptide in a gaseous and solvated environment.

Mechanism and energy landscape of domain swapping in the B1 domain of protein G

Three-dimensional domain swapping has emerged as a ubiquitous process for homo-oligomer formation in many unrelated proteins, but the molecular mechanism of this process is still poorly understood. Here we present a mechanism for the swapping reaction in the B1 domain of the immunoglobulin G binding protein from group G of Streptococcus (GB1). This is a particularly attractive system for investigating the swapping process, as the swapped dimer formed by the quadruple mutant (L5V/F30V/Y33F/A34F) of GB1 was recently shown to exist in equilibrium with a monomer-like conformation over time scales of minutes. According to our mechanism, swapping in GB1 starts from the C-terminus of the polypeptide chain and progresses by exchanging an increasing portion of the chains until a stable conformational state is reached. This exchange process does not involve unfolding. Rather, the conformational changes of individual monomers and their association are tightly coupled to minimize solvent exposure and maximize the total number of native contacts at all times, thereby closely approximating the minimum energy path of the reaction. Using detailed atomic descriptions, we compute the complete free-energy profiles of the exchange reaction for the GB1 quadruple mutant that forms swapped dimers and for the wild-type protein, which is monomeric. In both GB1 forms, intermediates sample a surprisingly wide range of nearly isoenergetic association modes and hinge conformations, indicating that the exchange reaction is a non-specific process akin to encounter complex formation where the amino acid sequence plays a marginal role. The main role of the mutations in the swapping process is to destabilize the GB1 monomer state, while stabilizing the swapped dimer conformation, with non-native intersubunit interactions, fostered by mutant side chains, contributing significantly to this stabilization. Our findings are rationalized in terms of a generic swapping mechanism that involves the association of activated molecular species, and it is argued that a similar mechanism may apply to swapping in other protein systems.

Ab initio calculations of free-energy reaction barriers

The theoretical description of chemical reactions was until recently limited to a 'static' approach in which important parameters such as the rate constant are deduced from the local topology of the potential energy surface close to minima and saddle points. Such an approach has, however, serious limitations. The growing computational power allows us now to use advanced simulation techniques to determine entropic effects accurately for medium-sized systems at ab initio level. Recently, we have implemented free-energy simulation techniques based on molecular dynamics, in particular on the blue-moon ensemble technique and on metadynamics, in the popular DFT code VASP. In the thermodynamic integration (blue-moon ensemble) technique, the free-energy profile is calculated as the path integral over the restoring forces along a parametrized reaction coordinate. In metadynamics, an image of the free-energy surface is constructed on the fly during the simulation by adding small repulsive Gaussian-shaped hills to the Lagrangian driving the dynamics. The two methods are tested on a simple chemical reaction-the nucleophilic substitution of methyl chloride by a chlorine anion.

The role of Li + , Na + , and K + in the ligand binding inside the human acetylcholinesterase gorge

Alkali cations can affect the catalytic efficiency of enzymes. This is particularly true when dealing with enzymes whose substrate bears a formal positive charge. Computational and biochemical approaches have been combined to shed light on the atomic aspects of the role of Li(+), Na(+), and K(+) on human acetylcholinesterase (hAChE) ligand binding. In this respect, molecular dynamics simulations and our recently developed metadynamics method were applied to study the entrance of the three cations in the gorge of hAChE, and their effect on the dynamical motion of a ligand (tetramethylammonium) from the bulk of the solvent into the deep narrow enzyme gorge. Furthermore, in order to support the theoretical results, K(M) and k(cat) for the acetylcholine hydrolysis in the presence of the three cations were evaluated by using an approach based on the Ellman's method. The combination of computational and biochemical experiments clearly showed that Li(+), Na(+), and K(+) may influence the ligand binding at the hAChE gorge.

First principles molecular dynamics study of catalytic reactions of biological macromolecular systems: toward analyses with QM/MM hybrid molecular simulations

First principles molecular dynamics simulations performed on a fully solvated RNA model structure allowed us to investigate the mechanism for enzymatic cleavage reactions, in vitro, of RNA enzymes (ribozymes). The concerted action of two metal catalysts turns out to be the most efficient way to promote, on the one hand, the proton abstraction from 2(')-OH that triggers the nucleophilic attack and, on the other hand, the cleavage of the P-O(5(')) bond. In fact, the elimination of one of the two metal cations leads to an increase in the activation energy of the reaction. The simulated pathway shows that an OH(-) in the coordination shell of the Mg(2+) close to O(2(')) promotes the initial proton abstraction and prevents its transfer to the ribozyme. This suggests that, in a real ribozyme, the double-metal-ion reaction mechanism is preferred with respect to single-metal-ion mechanisms either in the presence or in absence of the OH(-) anion. Finally, an insight into the importance of hybrid quantum mechanics/molecular mechanics (QM/MM) schemes is discussed in view of the modelling of a realistic system carrying all the features of a true ribozyme.

First-principles molecular dynamics study of proton transfer mechanism in bovine cytochrome c oxidase

Density functional based first-principles molecular dynamics calculations, performed on a model system extracted from the bovine cytochrome c oxidase, have been performed in an attempt to inspect the proton transfer mechanism across a peptide group. Our model system includes the specific Tyr440-Ser441 peptide group involved in a novel proton transfer path and shows that the Y440-S441 enol peptide group [-C(OH) = N-], which is a structural isomer of a keto form [-CO-NH-], is the product of the deprotonation of an imidic acid [-C(OH)-NH-] occurring in the vicinity of the deprotonated aspartic acid residue. For the subsequent enol-to-keto tautomerization, a direct H(+) transfer path in the Y440-S441 peptide group has been identified, in which the transition state takes a distorted four-membered ring structure.

The Conformational Free Energy Landscape of β-d-Glucopyranose. Implications for Substrate Preactivation in β-Glucoside Hydrolases

Using ab initio metadynamics we have computed the conformational free energy landscape of beta-D-glucopyranose as a function of the puckering coordinates. We show that the correspondence between the free energy and the Stoddard's pseudorotational itinerary for the system is rather poor. The number of free energy minima (9) is smaller than the number of ideal structures (13). Moreover, only six minima correspond to a canonical conformation. The structural features, the electronic properties, and the relative stability of the predicted conformers permit the rationalization of the occurrence of distorted sugar conformations in all the available X-ray structures of beta-glucoside hydrolase Michaelis complexes. We show that these enzymes recognize the most stable distorted conformers of the isolated substrate and at the same time the ones better prepared for catalysis in terms of bond elongation/shrinking and charge distribution. This suggests that the factors governing the distortions present in these complexes are largely dictated by the intrinsic properties of a single glucose unit.

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.

Towards a Rational Design of Ruthenium CO2 Hydrogenation Catalysts by Ab Initio Metadynamics

Complete reaction pathways relevant to CO2 hydrogenation by using a homogeneous ruthenium dihydride catalyst ([Ru(dmpe)2H2], dmpe=Me2PCH2CH2PMe2) have been investigated by ab initio metadynamics. This approach has allowed reaction intermediates to be identified and free-energy profiles to be calculated, which provide new insights into the experimentally observed reaction pathway. Our simulations indicate that CO2 insertion, which leads to the formation of formate complexes, proceeds by a concerted insertion mechanism. It is a rapid and direct process with a relatively low activation barrier, which is in agreement with experimental observations. Subsequent H2 insertion into the formate--Ru complex, which leads to the formation of formic acid, instead occurs via an intermediate [Ru(eta2-H2)] complex in which the molecular hydrogen coordinates to the ruthenium center and interacts weakly with the formate group. This step has been identified as the rate-limiting step. The reaction completes by hydrogen transfer from the [Ru(eta2-H2)] complex to the formate oxygen atom, which forms a dihydrogen-bonded Ru--HHO(CHO) complex. The activation energy for the H2 insertion step is lower for the trans isomer than for the cis isomer. A simple measure of the catalytic activity was proposed based on the structure of the transition state of the identified rate-limiting step. From this measure, the relationship between catalysts with different ligands and their experimental catalytic activities can be explained.

Possible Mechanism of Proton Transfer through Peptide Groups in the H-Pathway of the Bovine CytochromecOxidase

The peptide group connecting Tyr440 and Ser441 of the bovine cytochrome c oxidase is involved in a recently proposed proton-transfer path (H-path) where, at variance with other pathways (D- and K-paths), a usual hydrogen-bond network is interrupted, thus making this proton propagation rather unconventional. Our density-functional based molecular dynamics simulations show that, despite this anomaly and provided that a proton can reach a nearby water, a multistep proton-transfer pathway can become a viable pathway for such a reaction: a proton is initially transferred to the carbonyl oxygen of a keto form of the Tyr440-Ser441 peptide group [-CO-NH-], producing an imidic acid [-C(OH)-NH-] as a metastable state; the amide proton of the imidic acid is then transferred, spontaneously to the deprotonated carboxyl group of the Asp51 side chain, leading to the formation of an enol form [-C(OH)=N-] of the Tyr440-Ser441 peptide group. Then a subsequent enol-to-keto tautomerization occurs via a double proton-transfer path realized in the two adjacent Tyr440-Ser441 and Ser441-Asp442 peptide groups. An analysis of this multistep proton-transfer pathway shows that each elementary process occurs through the shortest distance, no permanent conformational changes are induced, thus preserving the X-ray crystal structure, and the reaction path is characterized by a reasonable activation barrier.

Free energy calculation of water addition coupled to reduction of aqueous RuO4-

Free energy calculations were carried out for water addition coupled reduction of aqueous ruthenate, RuO4-+H2O+e--->[RuO3(OH)2]2-, using Car-Parrinello molecular dynamics simulations. The full reaction is divided into the reduction of the tetrahedral monoanion, RuO4-+e--->RuO4(2-), followed by water addition, RuO4(2-)+H2O-->[RuO3(OH)2]2-. The free energy of reduction is computed from the fluctuations of the vertical energy gap using the MnO4-+e--->MnO(4)2- reaction as reference. The free energy for water addition is estimated using constrained molecular dynamics methods. While the description of this complex reaction, in principle, involves multiple reaction coordinates, we found that reversible transformation of the reactant into the product can be achieved by control of a single reaction coordinate consisting of a suitable linear combination of atomic distances. The free energy difference of the full reaction is computed to be -0.62 eV relative to the normal hydrogen electrode. This is in good agreement with the experimental value of -0.59 eV, lending further support to the hypothesis that, contrary to the ruthenate monoanion, the dianion is not tetrahedral but forms a trigonal-bipyramidal dihydroxo complex in aqueous solution. We construct an approximate two-dimensional free energy surface using the coupling parameter for reduction and the mechanical constraint for water addition as variables. Analyzing this surface we find that in the most favorable reaction pathway the reduction reaction precedes water addition. The latter takes place via the protonated complex [RuO3(OH)]- and subsequent transport of the created hydroxide ion to the fifth coordination site of Ru.

Collective Langevin dynamics of conformational motions in proteins

Functionally relevant slow conformational motions of proteins are, at present, in most cases inaccessible to molecular dynamics (MD) simulations. The main reason is that the major part of the computational effort is spend for the accurate description of a huge number of high frequency motions of the protein and the surrounding solvent. The accumulated influence of these fluctuations is crucial for a correct treatment of the conformational dynamics; however, their details can be considered irrelevant for most purposes. To accurately describe long time protein dynamics we here propose a reduced dimension approach, collective Langevin dynamics (CLD), which evolves the dynamics of the system within a small subspace of relevant collective degrees of freedom. The dynamics within the low-dimensional conformational subspace is evolved via a generalized Langevin equation which accounts for memory effects via memory kernels also extracted from short explicit MD simulations. To determine the memory kernel with differing levels of regularization, we propose and evaluate two methods. As a first test, CLD is applied to describe the conformational motion of the peptide neurotensin. A drastic dimension reduction is achieved by considering one single curved conformational coordinate. CLD yielded accurate thermodynamical and dynamical behaviors. In particular, the rate of transitions between two conformational states agreed well with a rate obtained from a 150 ns reference molecular dynamics simulation, despite the fact that the time scale of the transition (approximately 50 ns) was much longer than the 1 ns molecular dynamics simulation from which the memory kernel was extracted.

On the convergence improvement in the metadynamics simulations: a Wang-Landau recursion approach

As a popular tool in exploring free energy landscapes, the metadynamics method has been widely applied to elucidate various chemical or biochemical processes. As deeply discussed by Laio et al. [J. Phys. Chem. B 109, 6714 (2005)], the size of the updating Gaussian function is pivotal to the free energy convergence toward the target free energy surface. For instance, a greater Gaussian height can facilitate the quick visit of a conformation region of interest; however, it may lead to a larger error of the calculated free energy surface. In contrast, a lower Gaussian height can guarantee a better resolution of the calculated free energy surface; however, it will take longer time for such a simulation to navigate through the defined conformational region. In order to reconcile such confliction, the authors present a method by implementing the Wang-Landau recursion scheme in the metadynamics simulations to adaptively update the height of the unit Gaussian function. As demonstrated in their model studies on both a toy system, and a realistic molecular system treated with the hybrid quantum mechanical and molecular mechanical (QMMM) potential, the present approach can quickly result in more decently converged free energy surfaces, compared with the classical metadynamics simulations employing the fixed Gaussian heights.