A machine learning based intramolecular potential for a flexible organic molecule.

Quantum mechanical predictive modelling in chemistry and biology is often hindered by the long time scales and large system sizes required of the computational model. Here, we employ the kernel regression machine learning technique to construct an analytical potential, using the Gaussian Approximation Potential software and framework, that reproduces the quantum mechanical potential energy surface of a small, flexible, drug-like molecule, 3-(benzyloxy)pyridin-2-amine. Challenges linked to the high dimensionality of the configurational space of the molecule are overcome by developing an iterative training protocol and employing a representation that separates short and long range interactions. The analytical model is connected to the MCPRO simulation software, which allows us to perform Monte Carlo simulations of the small molecule bound to two proteins, p38 MAP kinase and leukotriene A4 hydrolase, as well as in water. We demonstrate that our machine learning based intramolecular model is transferable to the condensed phase, and demonstrate that the use of a faithful representation of the quantum mechanical potential energy surface can result in corrections to absolute protein-ligand binding free energies of up to 2 kcal mol-1 in the example studied here.

Ab initio molecular dynamics study of solvated electrons in methanol clusters.

We performed a series of ab initio molecular dynamics simulations to investigate the physical properties of small methanol cluster anions, [(CH3OH)n]-, (n = 8-32). An excess electron was attached to neutral clusters that were prepared to accommodate the electron in interior cavity states or surface bound states. The computed initial binding energies of the electrons to these clusters indicate appealing similarity to the experimentally observed vertical detachment energies. The tendency of the interior state clusters parallels that of the clusters with strong electron binding in the experiments, while the simulated unrelaxed surface state anions are similar to the observed weakly bound species. This assignment is consistent with a previous identification based on hybrid quantum-classical simulations. The time evolution of the cluster anions suggests that interior state electrons slowly move to and relax on the surface, in excess electronic states that appear significantly more stable than the experimentally assigned putative surface states. Based on this result we predict the existence of relaxed surface state isomers of small methanol cluster anions. Due to the kinetic metastability of the experimentally found weakly bound species, we anticipate a serious technical challenge to prepare and identify small methanol cluster anions with relaxed surface states. These more strongly binding surface states are stabilized by dangling hydroxyl hydrogen atoms pointing to the excess electron's charge distribution. In addition, methyl hydrogens also appear to contribute to the stability of these states. During its transition to the surface, the interior excess electron maintains its initial solvent cavity. No signs of non-cavity interior states are observed in the present first principles ab initio molecular dynamics simulations.

Preconditioners for the geometry optimisation and saddle point search of molecular systems.

A class of preconditioners is introduced to enhance geometry optimisation and transition state search of molecular systems. We start from the Hessian of molecular mechanical terms, decompose it and retain only its positive definite part to construct a sparse preconditioner matrix. The construction requires only the computation of the gradient of the corresponding molecular mechanical terms that are already available in popular force field software packages. For molecular crystals, the preconditioner can be combined straightforwardly with the exponential preconditioner recently introduced for periodic systems. The efficiency is demonstrated on several systems using empirical, semiempirical and ab initio potential energy surfaces.

Excess electrons in methanol clusters: Beyond the one-electron picture.

We performed a series of comparative quantum chemical calculations on various size negatively charged methanol clusters, CH3OHn-. The clusters are examined in their optimized geometries (n = 2-4), and in geometries taken from mixed quantum-classical molecular dynamics simulations at finite temperature (n = 2-128). These latter structures model potential electron binding sites in methanol clusters and in bulk methanol. In particular, we compute the vertical detachment energy (VDE) of an excess electron from increasing size methanol cluster anions using quantum chemical computations at various levels of theory including a one-electron pseudopotential model, several density functional theory (DFT) based methods, MP2 and coupled-cluster CCSD(T) calculations. The results suggest that at least four methanol molecules are needed to bind an excess electron on a hydrogen bonded methanol chain in a dipole bound state. Larger methanol clusters are able to form stronger interactions with an excess electron. The two simulated excess electron binding motifs in methanol clusters, interior and surface states, correlate well with distinct, experimentally found VDE tendencies with size. Interior states in a solvent cavity are stabilized significantly stronger than electron states on cluster surfaces. Although we find that all the examined quantum chemistry methods more or less overestimate the strength of the experimental excess electron stabilization, MP2, LC-BLYP, and BHandHLYP methods with diffuse basis sets provide a significantly better estimate of the VDE than traditional DFT methods (BLYP, B3LYP, X3LYP, PBE0). A comparison to the better performing many electron methods indicates that the examined one-electron pseudopotential can be reasonably used in simulations for systems of larger size.

Exploration, Sampling, And Reconstruction of Free Energy Surfaces with Gaussian Process Regression.

Practical free energy reconstruction algorithms involve three separate tasks: biasing, measuring some observable, and finally reconstructing the free energy surface from those measurements. In more than one dimension, adaptive schemes make it possible to explore only relatively low lying regions of the landscape by progressively building up the bias toward the negative of the free energy surface so that free energy barriers are eliminated. Most schemes use the final bias as their best estimate of the free energy surface. We show that large gains in computational efficiency, as measured by the reduction of time to solution, can be obtained by separating the bias used for dynamics from the final free energy reconstruction itself. We find that biasing with metadynamics, measuring a free energy gradient estimator, and reconstructing using Gaussian process regression can give an order of magnitude reduction in computational cost.

A universal preconditioner for simulating condensed phase materials.

We introduce a universal sparse preconditioner that accelerates geometry optimisation and saddle point search tasks that are common in the atomic scale simulation of materials. Our preconditioner is based on the neighbourhood structure and we demonstrate the gain in computational efficiency in a wide range of materials that include metals, insulators, and molecular solids. The simple structure of the preconditioner means that the gains can be realised in practice not only when using expensive electronic structure models but also for fast empirical potentials. Even for relatively small systems of a few hundred atoms, we observe speedups of a factor of two or more, and the gain grows with system size. An open source Python implementation within the Atomic Simulation Environment is available, offering interfaces to a wide range of atomistic codes.

The adaptive buffered force QM/MM method in the CP2K and AMBER software packages.

The implementation and validation of the adaptive buffered force (AdBF) quantum-mechanics/molecular-mechanics (QM/MM) method in two popular packages, CP2K and AMBER are presented. The implementations build on the existing QM/MM functionality in each code, extending it to allow for redefinition of the QM and MM regions during the simulation and reducing QM-MM interface errors by discarding forces near the boundary according to the buffered force-mixing approach. New adaptive thermostats, needed by force-mixing methods, are also implemented. Different variants of the method are benchmarked by simulating the structure of bulk water, water autoprotolysis in the presence of zinc and dimethyl-phosphate hydrolysis using various semiempirical Hamiltonians and density functional theory as the QM model. It is shown that with suitable parameters, based on force convergence tests, the AdBF QM/MM scheme can provide an accurate approximation of the structure in the dynamical QM region matching the corresponding fully QM simulations, as well as reproducing the correct energetics in all cases. Adaptive unbuffered force-mixing and adaptive conventional QM/MM methods also provide reasonable results for some systems, but are more likely to suffer from instabilities and inaccuracies.

The role of reorganization energy in rational enzyme design.

Computational design is becoming an integral component in developing novel enzymatic activities. Catalytic efficiencies of man-made enzymes however are far behind their natural counterparts. The discrepancy between laboratory and naturally evolved enzymes suggests that a major catalytic factor is still missing in the computational process. Reorganization energy, which is the origin of catalytic power of natural enzymes, has not been exploited yet for design. As exemplified in case of KE07 Kemp eliminase, this quantity is optimized by directed evolution. Mutations beneficial for evolution, but without direct impact on catalysis can be identified based on contributions to reorganization energy. We propose to incorporate the reorganization energy in scaffold selection to provide highly evolvable initial designs.

Tests of an adaptive QM/MM calculation on free energy profiles of chemical reactions in solution.

We present reaction free energy calculations using the adaptive buffered force mixing quantum mechanics/molecular mechanics (bf-QM/MM) method. The bf-QM/MM method combines nonadaptive electrostatic embedding QM/MM calculations with extended and reduced QM regions to calculate accurate forces on all atoms, which can be used in free energy calculation methods that require only the forces and not the energy. We calculate the free energy profiles of two reactions in aqueous solution: the nucleophilic substitution reaction of methyl chloride with a chloride anion and the deprotonation reaction of the tyrosine side chain. We validate the bf-QM/MM method against a full QM simulation, and show that it correctly reproduces both geometrical properties and free energy profiles of the QM model, while the electrostatic embedding QM/MM method using a static QM region comprising only the solute is unable to do so. The bf-QM/MM method is not explicitly dependent on the details of the QM and MM methods, so long as it is possible to compute QM forces in a small region and MM forces in the rest of the system, as in a conventional QM/MM calculation. It is simple, with only a few parameters needed to control the QM calculation sizes, and allows (but does not require) a varying and adapting QM region which is necessary for simulating solutions.

Empirical valence bond simulations of the chemical mechanism of ATP to cAMP conversion by anthrax edema factor.

The two-metal catalysis by the adenylyl cyclase domain of the anthrax edema factor toxin was simulated using the empirical valence bond (EVB) quantum mechanical/molecular mechanical approach. These calculations considered the energetics of the nucleophile deprotonation and the formation of a new P-O bond in aqueous solution and in the enzyme-substrate complex present in the crystal structure models of the reactant and product states of the reaction. Our calculations support a reaction pathway that involves metal-assisted transfer of a proton from the nucleophile to the bulk aqueous solution followed by subsequent formation of an unstable pentavalent intermediate that decomposes into cAMP and pyrophosphate (PPi). This pathway involves ligand exchange in the first solvation sphere of the catalytic metal. At 12.9 kcal/mol, the barrier for the last step of the reaction, the cleavage of the P-O bond to PPi, corresponds to the highest point on the free energy profile for this reaction pathway. However, this energy is too close to the value of 11.4 kcal/mol calculated for the barrier of the nucleophilic attack step to reach a definitive conclusion about the rate-limiting step. The calculated reaction mechanism is supported by reasonable agreement between the experimental and calculated catalytic rate constant decrease caused by the mutation of the active site lysine 346 to arginine.

Topologically invariant reaction coordinates for simulating multistate chemical reactions.

Evaluating free energy profiles of chemical reactions in complex environments such as solvents and enzymes requires extensive sampling, which is usually performed by potential of mean force (PMF) techniques. The reliability of the sampling depends not only on the applied PMF method but also the reaction coordinate space within the dynamics is biased. In contrast to simple geometrical collective variables that depend only on the positions of the atomic coordinates of the reactants, the E(gap) reaction coordinate (the energy difference obtained by evaluating a suitable force field using reactant and product state topologies) has the unique property that it is able to take environmental effects into account leading to better convergence, a more faithful description of the transition state ensemble and therefore more accurate free energy profiles. However, E(gap) requires predefined topologies and is therefore inapplicable for multistate reactions, in which the barrier between the chemically equivalent topologies is comparable to the reaction activation barrier, because undesired "side reactions" occur. In this article, we introduce a new energy-based collective variable by generalizing the E(gap) reaction coordinate such that it becomes invariant to equivalent topologies and show that it yields more well behaved free energy profiles than simpler geometrical reaction coordinates.

Quantum-classical simulation of electron localization in negatively charged methanol clusters.

A series of quantum molecular dynamics simulations have been performed to investigate the energetic, structural, dynamic, and spectroscopic properties of methanol cluster anions, [(CH(3)OH)(n)](-), (n = 50-500). Consistent with the inference from photo-electron imaging experiments, we find two main localization modes of the excess electron in equilibrated methanol clusters at ∼200 K. The two different localization patterns have strikingly different physical properties, consistent with experimental observations, and are manifest in comparable cluster sizes to those observed. Smaller clusters (n ≤ 128) tend to localize the electron in very weakly bound, diffuse electronic states on the surface of the cluster, while in larger ones the electron is stabilized in solvent cavities, in compact interior-bound states. The interior states exhibit properties that largely resemble and smoothly extrapolate to those simulated for a solvated electron in bulk methanol. The surface electronic states of methanol cluster anions are significantly more weakly bound than the surface states of the anionic water clusters. The key source of the difference is the lack of stabilizing free hydroxyl groups on a relaxed methanol cluster surface. We also provide a mechanistic picture that illustrates the essential role of the interactions of the excess electron with the hydroxyl groups in the dynamic process of the transition of the electron from surface-bound states to interior-bound states.

Analysis of localization sites for an excess electron in neutral methanol clusters using approximate pseudopotential quantum-mechanical calculations.

We have used a recently developed electron-methanol molecule pseudopotential in approximate quantum mechanical calculations to evaluate and statistically analyze the physical properties of an excess electron in the field of equilibrated neutral methanol clusters ((CH(3)OH)(n), n=50-500). The methanol clusters were generated in classical molecular dynamics simulations at nominal 100 and 200 K temperatures. Topological analysis of the neutral clusters indicates that methyl groups cover the surface of the clusters almost exclusively, while the associated hydroxyl groups point inside. Since the initial neutral clusters are lacking polarity on the surface and compact inside, the excess electron can barely attach to these structures. Nevertheless, most of the investigated cluster configurations do support weakly stabilized cluster anion states. We find that similarly to water clusters, the pre-existing instantaneous dipole moment of the neutral clusters binds the electron. The localizing electrons occupy diffuse, weakly bound surface states that largely engulf the cluster although their centers are located outside the cluster molecular frame. The initial localization of the excess electron is reflected in its larger radius compared to water due to the lack of free OH hydrogens on the cluster surface. The stabilization of the excess electron increases, while the radius decreases monotonically as the clusters grow in size. Stable, interior bound states of the excess electron are not observed to form neither in finite size methanol clusters nor in the equilibrium bulk.

A new electron-methanol molecule pseudopotential and its application for the solvated electron in methanol.

A new electron-methanol molecule pseudopotential is developed and tested in the present paper. The formal development of the potential is based on quantum mechanical calculations on the electron-methanol molecule model in the static exchange approximation. The computational model includes a steep confining potential that keeps the otherwise unbound excess electron in the vicinity of the methanol molecule. Using the Phillips-Kleinman theorem we introduce a smooth pseudowave function of the excess electron with the exact eigenenergy and correct asymptotic behavior. The nonlocal potential energy operator of the model Hamiltonian is then replaced to a local potential that reproduces the ground-state properties of the excess electron satisfactorily. The pseudopotential is then optimized in an analytically simple functional form to fit this approximate local potential in conjunction with the point charges and the geometry of a classical, all-site methanol-methanol interaction potential. Of the adjustable parameters, the parameters for the carbon and the methyl hydrogen atoms are optimized, while those for the oxygen and the hydroxyl hydrogen are taken from a previous electron-water molecule pseudopotential. A polarization term is added to the potential a posteriori. The polarization parameters are chosen to reproduce the experimental position of the optical absorption spectrum of an excess electron in mixed quantum-classical molecular dynamics simulations. The energetic, structural and spectroscopic properties of the solvated electron in a methanol bath are simulated at 300 K and compared with previous solvated electron simulations and available experimental data.

On the divalent metal ion dependence of DNA cleavage by restriction endonucleases of the EcoRI family.

Restriction endonucleases of the PD...D/EXK family need Mg(2+) for DNA cleavage. Whereas Mg(2+) (or Mn(2+)) promotes catalysis, Ca(2+) (without Mg(2+)) only supports DNA binding. The role of Mg(2+) in DNA cleavage by restriction endonucleases has elicited many hypotheses, differing mainly in the number of Mg(2+) involved in catalysis. To address this problem, we measured the Mg(2+) and Mn(2+) concentration dependence of DNA cleavage by BamHI, BglII, Cfr10I, EcoRI, EcoRII (catalytic domain), MboI, NgoMIV, PspGI, and SsoII, which were reported in co-crystal structure analyses to bind one (BglII and EcoRI) or two (BamHI and NgoMIV) Me(2+) per active site. DNA cleavage experiments were carried out at various Mg(2+) and Mn(2+) concentrations at constant ionic strength. All enzymes show a qualitatively similar Mg(2+) and Mn(2+) concentration dependence. In general, the Mg(2+) concentration optimum (between approximately 1 and 10 mM) is higher than the Mn(2+) concentration optimum (between approximately 0.1 and 1 mM). At still higher Mg(2+) or Mn(2+) concentrations, the activities of all enzymes tested are reduced but can be reactivated by Ca(2+). Based on these results, we propose that one Mg(2+) or Mn(2+) is critical for restriction enzyme activation, and binding of a second Me(2+) plays a role in modulating the activity. Steady-state kinetics carried out with EcoRI and BamHI suggest that binding of a second Mg(2+) or Mn(2+) mainly leads to an increase in K(m), such that the inhibitory effect of excess Mg(2+) or Mn(2+) can be overcome by increasing the substrate concentration. Our conclusions are supported by molecular dynamics simulations and are consistent with the structural observations of both one and two Me(2+) binding to these enzymes.

The energy gap as a universal reaction coordinate for the simulation of chemical reactions.

The selection of a proper reaction coordinate is a major bottleneck in simulations of chemical reactions in complex systems. Increasing the number of variables that are used to bias the reaction largely affects the convergence and leads to an unbearable increase in computational price. This problem can be overcome by employing a complex reaction coordinate that depends on many geometrical variables of the system, such as the energy gap (EGAP) in the empirical valence bond (EVB) method. EGAP depends on all of the coordinates of the system, and its robustness has been demonstrated for a variety of enzymatic reactions. In this work, we demonstrate that EGAP, derived from a classical representation, can be used as a reaction coordinate in systems described with any quantum chemistry Hamiltonian. Benefits of using EGAP as a reaction coordinate as compared to a traditional geometrical variable are illustrated in the case of a symmetric nucleophilic substitution reaction in water solution. EGAP is shown to provide a significantly more efficient sampling and allows a better localization of the transition state as compared to a geometrical reaction coordinate.

Probing the two-metal ion mechanism in the restriction endonuclease BamHI.

The choreography of restriction endonuclease catalysis is a long-standing paradigm in molecular biology. Bivalent metal ions are required almost for all PD..D/ExK type enzymes, but the number of cofactors essential for the DNA backbone scission remained ambiguous. On the basis of crystal structures and biochemical data for various restriction enzymes, three models have been developed that assign critical roles for one, two, or three metal ions during the phosphodiester hydrolysis. To resolve this apparent controversy, we investigated the mechanism of BamHI catalysis using quantum mechanical/molecular mechanical simulation techniques and determined the activation barriers of three possible pathways that involve a Glu-113 or a neighboring water molecule as a general base or an external nucleophile that penetrated from bulk solution. The extrinsic mechanism was found to be the most favorable with an activation free energy of 23.4 kcal/mol, in reasonable agreement with the experimental data. On the basis of the effect of the individual metal ions on the activation barrier, metal ion A was concluded to be pivotal for the reaction, while the enzyme lacking metal ion B still has moderate efficiency. Thus, we propose that the catalytic scheme of BamHI does not involve a general base for nucleophile generation and requires one obligatory metal ion for catalysis that stabilizes the attacking nucleophile and coordinates it throughout the nucleophilic attack. Such a model may also explain the variation in the number of metal ions in the crystal structures and thus could serve as a framework for a unified catalytic scheme of type II restriction endonucleases.

Metal-binding sites at the active site of restriction endonuclease BamHI can conform to a one-ion mechanism.

The number of metal ions required for phosphoryl transfer in restriction endonucleases is still an unresolved question in molecular biology. The two Ca(2+) and Mn(2+) ions observed in the pre- and post-reactive complexes of BamHI conform to the classical two-metal ion choreography. We probed the Mg(2+) cofactor positions at the active site of BamHI by molecular dynamics simulations with one and two metal ions present and identified several catalytically relevant sites. These can mark the pathway of a single ion during catalysis, suggesting its critical role, while a regulatory function is proposed for a possible second ion.