Sampling of Rare Events Using Hidden Restraints

A Robust and Versatile QM/MM Interface for Molecular Dynamics in GROMOS

The integration of quantum mechanics and molecular mechanics (QM/MM) within molecular dynamics simulations is crucial to accurately model complex biochemical systems. Here, we present an enhanced implementation of the QM/MM interface in the GROMOS simulation package, introducing significant improvements in functionality and user control. We present new features, including the link atom scheme, which allows the modeling of QM regions as a part of bigger molecules. Benchmark tests on various systems, including QM water in water, amino acids in water, and tripeptides validate the reliability of the new functionalities. Performance evaluations demonstrate that the updated implementation is efficient, with the primary computational burden attributed to the QM program rather than the QM/MM interface or the MD program itself. The improved QM/MM interface enables more advanced investigations into biomolecular reactivity, enzyme catalysis, and other phenomena requiring detailed quantum mechanical treatment within classical simulations. This work represents a significant advancement in the capabilities of GROMOS, providing enhanced tools to explore complex molecular systems.

Methods for Classical-Mechanical Molecular Simulation in Chemistry: Achievements, Limitations, Perspectives

More than a half century ago it became feasible to simulate, using classical-mechanical equations of motion, the dynamics of molecular systems on a computer. Since then classical-physical molecular simulation has become an integral part of chemical research. It is widely applied in a variety of branches of chemistry and has significantly contributed to the development of chemical knowledge. It offers understanding and interpretation of experimental results, semiquantitative predictions for measurable and nonmeasurable properties of substances, and allows the calculation of properties of molecular systems under conditions that are experimentally inaccessible. Yet, molecular simulation is built on a number of assumptions, approximations, and simplifications which limit its range of applicability and its accuracy. These concern the potential-energy function used, adequate sampling of the vast statistical-mechanical configurational space of a molecular system and the methods used to compute particular properties of chemical systems from statistical-mechanical ensembles. During the past half century various methodological ideas to improve the efficiency and accuracy of classical-physical molecular simulation have been proposed, investigated, evaluated, implemented in general simulation software or were abandoned. The latter because of fundamental flaws or, while being physically sound, computational inefficiency. Some of these methodological ideas are briefly reviewed and the most effective methods are highlighted. Limitations of classical-physical simulation are discussed and perspectives are sketched.

On the use of intra-molecular distance and angle constraints to lengthen the time step in molecular and stochastic dynamics simulations of proteins

Computer simulation of proteins in aqueous solution at the atomic level of resolution is still limited in time span and system size due to limited computing power available and thus employs a variety of time‐saving techniques that trade some accuracy against computational effort. An example of such a time‐saving technique is the application of constraints to particular degrees of freedom when integrating Newton's or Langevin's equations of motion in molecular dynamics (MD) or stochastic dynamics (SD) simulations, respectively. The application of bond‐length constraints is standard practice in protein simulations and allows for a lengthening of the time step by a factor of three. Applying recently proposed algorithms to constrain bond angles or dihedral angles, it is investigated, using the protein trypsin inhibitor as test molecule, whether bond angles and dihedral angles involving hydrogen atoms or even stiff proper (torsional) dihedral angles as well as improper ones (maintaining particular tetrahedral or planar geometries) may be constrained without generating too many artificial side effects. Constraining the relative positions of the hydrogen atoms in the protein allows for a lengthening of the time step by a factor of two. Additionally constraining the improper dihedral angles and the stiff proper (torsional) dihedral angles in the protein does not allow for an increase of the MD or SD time step.

Algorithms to apply dihedral-angle constraints in molecular or stochastic dynamics simulations

Various algorithms to apply dihedral-angle constraints in molecular dynamics or stochastic dynamics simulations of molecular systems are presented, investigated, and tested. They use Cartesian coordinates and determine the Lagrangian multipliers necessary for maintaining the constraints iteratively. The most suitable algorithm to maintain a dihedral-angle constraint is numerically compared to the alternative to use distance constraints to this end. It can easily be used to obtain a potential of mean force along a dihedral-angle coordinate.

Simulation of Reversible Protein-Protein Binding and Calculation of Binding Free Energies Using Perturbed Distance Restraints

Virtually all biological processes depend on the interaction between proteins at some point. The correct prediction of biomolecular binding free-energies has many interesting applications in both basic and applied pharmaceutical research. While recent advances in the field of molecular dynamics (MD) simulations have proven the feasibility of the calculation of protein–protein binding free energies, the large conformational freedom of proteins and complex free energy landscapes of binding processes make such calculations a difficult task. Moreover, convergence and reversibility of resulting free-energy values remain poorly described. In this work, an easy-to-use, yet robust approach for the calculation of standard-state protein–protein binding free energies using perturbed distance restraints is described. In the binding process the conformations of the proteins were restrained, as suggested earlier. Two approaches to avoid end-state problems upon release of the conformational restraints were compared. The method was evaluated by practical application to a small model complex of ubiquitin and the very flexible ubiquitin-binding domain of human DNA polymerase ι (UBM2). All computed free energy differences were closely monitored for convergence, and the calculated binding free energies had a mean unsigned deviation of only 1.4 or 2.5 kJ·mol^–1 from experimental values. Statistical error estimates were in the order of thermal noise. We conclude that the presented method has promising potential for broad applicability to quantitatively describe protein–protein and various other kinds of complex formation.

Multistate λ-local-elevation umbrella-sampling (MS-λ-LEUS): method and application to the complexation of cations by crown ethers

An extension of the λ-local-elevation umbrella-sampling (λ-LEUS) scheme [ Bieler et al. J. Chem. Theory Comput. 2014 , 10 , 3006 ] is proposed to handle the multistate (MS) situation, i.e. the calculation of the relative free energies of multiple physical states based on a single simulation. The key element of the MS-λ-LEUS approach is to use a single coupling variable Λ controlling successive pairwise mutations between the states of interest in a cyclic fashion. The Λ variable is propagated dynamically as an extended-system variable, using a coordinate transformation with plateaus and a memory-based biasing potential as in λ-LEUS. Compared to other available MS schemes (one-step perturbation, enveloping distribution sampling and conventional λ-dynamics) the proposed method presents a number of important advantages, namely: (i) the physical states are visited explicitly and over finite time periods; (ii) the extent of unphysical space required to ensure transitions is kept minimal and, in particular, one-dimensional; (iii) the setup protocol solely requires the topologies of the physical states; and (iv) the method only requires limited modifications in a simulation code capable of handling two-state mutations. As an initial application, the absolute binding free energies of five alkali cations to three crown ethers in three different solvents are calculated. The results are found to reproduce qualitatively the main experimental trends and, in particular, the experimental selectivity of 18C6 for K(+) in water and methanol, which is interpreted in terms of opposing trends along the cation series between the solvation free energy of the cation and the direct electrostatic interactions within the complex.

A Multiscale Treatment of Angeli's Salt Decomposition

Sodium trioxodinitrate's (Na2N2O3, Angeli's salt) unique cardiovascular effects have been associated with its ability to yield HNO upon dissociation under physiological conditions. Due to its potential applications in new therapies for heart failure, the dissociation of Angeli's salt has recently received increased attention. The decomposition mechanism has been previously studied by quantum mechanical methods using a continuum approximation (PCM) for the solvent effects. In this work we use our recently developed interface of the Amber and Gaussian packages via the PUPIL package to study Angeli's salt dissociation in a hybrid QM/MM scheme where the water solvent molecules are treated explicitly with classical mechanics while the solute is treated with full quantum mechanics (UB3LYP/6-31+G(d) and UMP2/6-31+G(d)) level. Multiple steered molecular dynamics was used with the Jarzynski relationship to extract the free energy profile for the process. We obtain 4.8 kcal mol(-1) and 6.4 kcal mol(-1) free energy barriers for the N-N bond breaking for UB3LYP and UMP2, respectively. The geometries and Mulliken charges for reactant, transition state, and products have been characterized through a number of hybrid QM/MM molecular dynamics runs with the N-N distance restrained to representative values of each species. The results highlight the role of individual solvent molecules for the reaction energetics and provide a comparison point against implicit solvation methods.

Flexible Boundaries for Multiresolution Solvation: An Algorithm for Spatial Multiscaling in Molecular Dynamics Simulations

An algorithm is proposed for performing molecular dynamics (MD) simulations of a biomolecular solute represented at atomistic resolution surrounded by a surface layer of atomistic fine-grained (FG) solvent molecules within a bulk represented by coarse-grained (CG) solvent beads. The method, called flexible boundaries for multiresolution solvation (FBMS), is based on: (i) a three-region layering of the solvent around the solute, involving an FG layer surrounded by a mixed FG-CG buffer layer, itself surrounded by a bulk CG region; (ii) a definition of the layer boundary that relies on an effective distance to the solute surface and is thus adapted to the shape of the solute as well as adjusts to its conformational changes. The effective surface distance is defined by inverse-nth power averaging over the distances to all non-hydrogen solute atoms, and the layering is enforced by means of half-harmonic distance restraints, attractive for the FG molecules and repulsive for the CG beads. A restraint-free region at intermediate distances enables the formation of the buffer layer, where the FG and CG solvents can mix freely. The algorithm is tested and validated using the GROMOS force field and the associated FG (SPC) and CG (polarizable CGW) water models. The test systems include pure-water systems, where one FG molecule plays the role of a solute, and a deca-alanine peptide with two widely different solute shapes considered, α-helical and fully extended. In particular, as the peptide unfolds, the number of FG molecules required to fill its close-range solvation layer increases, with the additional molecules being provided by the buffer layer. Further validation involves simulations of four proteins in multiresolution FG/CG mixtures. The resulting structural, energetic, and solvation properties are found to be similar to those observed in corresponding pure FG simulations.

Practical Aspects of Free-Energy Calculations: A Review

Free-energy calculations in the framework of classical molecular dynamics simulations are nowadays used in a wide range of research areas including solvation thermodynamics, molecular recognition, and protein folding. The basic components of a free-energy calculation, that is, a suitable model Hamiltonian, a sampling protocol, and an estimator for the free energy, are independent of the specific application. However, the attention that one has to pay to these components depends considerably on the specific application. Here, we review six different areas of application and discuss the relative importance of the three main components to provide the reader with an organigram and to make nonexperts aware of the many pitfalls present in free energy calculations.

Accelerated adaptive integration method

Conformational changes that occur upon ligand binding may be too slow to observe on the time scales routinely accessible using molecular dynamics simulations. The adaptive integration method (AIM) leverages the notion that when a ligand is either fully coupled or decoupled, according to λ, barrier heights may change, making some conformational transitions more accessible at certain λ values. AIM adaptively changes the value of λ in a single simulation so that conformations sampled at one value of λ seed the conformational space sampled at another λ value. Adapting the value of λ throughout a simulation, however, does not resolve issues in sampling when barriers remain high regardless of the λ value. In this work, we introduce a new method, called Accelerated AIM (AcclAIM), in which the potential energy function is flattened at intermediate values of λ, promoting the exploration of conformational space as the ligand is decoupled from its receptor. We show, with both a simple model system (Bromocyclohexane) and the more complex biomolecule Thrombin, that AcclAIM is a promising approach to overcome high barriers in the calculation of free energies, without the need for any statistical reweighting or additional processors.

Free-energy differences between states with different conformational ensembles

Multiple conformations separated by high-energy barriers represent a challenging problem in free-energy calculations due to the difficulties in achieving adequate sampling. We present an application of thermodynamic integration (TI) in conjunction with the local elevation umbrella sampling (LE/US) method to improve convergence in alchemical free-energy calculations. TI-LE/US was applied to the guanosine triphosphate (GTP) to 8-Br-GTP perturbation, molecules that present high-energy barriers between the anti and syn states and that have inverted preferences for those states. The convergence and reliability of TI-LE/US was assessed by comparing with previous results using the enhanced-sampling one-step perturbation (OSP) method. A linear interpolation of the end-state biasing potentials was sufficient to dramatically improve sampling along the chosen reaction coordinate. Conformational free-energy differences were also computed for the syn and anti states and compared to experimental and theoretical results. Additionally, a coupled OSP with LE/US was carried out, allowing the calculation of conformational and alchemical free energies of GTP and 8-substituted GTP analogs.

Protein-Ligand Binding from Distancefield Distances and Hamiltonian Replica Exchange Simulations

The calculation of protein-ligand binding free energies is an important goal in the field of computational chemistry. Applying path-sampling methods for this purpose involves calculating the associated potential of mean force (PMF) and gives insight into the binding free energy along the binding process. Without a priori knowledge about the binding path, sampling reversible binding can be difficult to achieve. To alleviate this problem, we introduce the distancefield (DF) as a reaction coordinate for such calculations. DF is a grid-based method in which the shortest distance between the binding site and a ligand is determined avoiding routes that pass through the protein. Combining this reaction coordinate with Hamiltonian replica exchange molecular dynamics (HREMD) allows for the reversible binding of the ligand to the protein. A comparison is made between umbrella sampling using regular distance restraints and HREMD with DF restraints to study aspirin binding to the protein phospholipase A2. Although the free energies of binding are similar for both methods, the increased sampling with HREMD has a significant influence on the shape of the PMF. A remarkable agreement between the calculated binding free energies from the PMF and the experimental estimate is obtained.

A reoptimized GROMOS force field for hexopyranose-based carbohydrates accounting for the relative free energies of ring conformers, anomers, epimers, hydroxymethyl rotamers, and glycosidic linkage conformers

This article presents a reoptimization of the GROMOS 53A6 force field for hexopyranose-based carbohydrates (nearly equivalent to 45A4 for pure carbohydrate systems) into a new version 56A(CARBO) (nearly equivalent to 53A6 for non-carbohydrate systems). This reoptimization was found necessary to repair a number of shortcomings of the 53A6 (45A4) parameter set and to extend the scope of the force field to properties that had not been included previously into the parameterization procedure. The new 56A(CARBO) force field is characterized by: (i) the formulation of systematic build-up rules for the automatic generation of force-field topologies over a large class of compounds including (but not restricted to) unfunctionalized polyhexopyranoses with arbritrary connectivities; (ii) the systematic use of enhanced sampling methods for inclusion of experimental thermodynamic data concerning slow or unphysical processes into the parameterization procedure; and (iii) an extensive validation against available experimental data in solution and, to a limited extent, theoretical (quantum-mechanical) data in the gas phase. At present, the 56A(CARBO) force field is restricted to compounds of the elements C, O, and H presenting single bonds only, no oxygen functions other than alcohol, ether, hemiacetal, or acetal, and no cyclic segments other than six-membered rings (separated by at least one intermediate atom). After calibration, this force field is shown to reproduce well the relative free energies of ring conformers, anomers, epimers, hydroxymethyl rotamers, and glycosidic linkage conformers. As a result, the 56A(CARBO) force field should be suitable for: (i) the characterization of the dynamics of pyranose ring conformational transitions (in simulations on the microsecond timescale); (ii) the investigation of systems where alternative ring conformations become significantly populated; (iii) the investigation of anomerization or epimerization in terms of free-energy differences; and (iv) the design of simulation approaches accelerating the anomerization process along an unphysical pathway.

Efficient free energy calculations for compounds with multiple stable conformations separated by high energy barriers

Compounds with high intramolecular energy barriers represent challenging targets for free energy calculations because of the difficulty to obtain sufficient conformational sampling. Existing approaches are therefore computationally very demanding, thus preventing practical applications for such compounds. We present an enhanced sampling-one step perturbation method (ES-OS) to tackle this problem in a highly efficient way. A single molecular dynamics simulation of a judiciously chosen reference state (using two sets of soft-core interactions) is sufficient to determine conformational distributions of chemically similar compounds and the free energy differences between them. The ES-OS method is applied to a set of five biologically relevant 8-substituted GTP analogs having high energy barriers between the anti and the syn conformations of the base with respect to the ribose part. The reliability of ES-OS is verified by comparing the results to Hamiltonian replica exchange simulations of GTP and 8-Br-GTP and the experimentally determined 3J(C4,H1') coupling constant for GMP in water. Additional simulations in vacuum and octanol allow us to calculate differences in the solvation free energies and in lipophilicities (log P). Free energy contributions from individual conformational regions are also calculated, and their relationship with the overall free energy is derived leading to a set of multiconformational free energy formulas. These relationships are of general applicability and can be used in free energy calculations for a more diverse set of compounds.

Conformation, dynamics, solvation and relative stabilities of selected β-hexopyranoses in water: a molecular dynamics study with the gromos 45A4 force field

The present article reports long timescale (200 ns) simulations of four beta-D-hexopyranoses (beta-D-glucose, beta-D-mannose, beta-D-galactose and beta-D-talose) using explicit-solvent (water) molecular dynamics and vacuum stochastic dynamics simulations together with the GROMOS 45A4 force field. Free-energy and solvation free-energy differences between the four compounds are also calculated using thermodynamic integration. Along with previous experimental findings, the present results suggest that the formation of intramolecular hydrogen-bonds in water is an 'opportunistic' consequence of the close proximity of hydrogen-bonding groups, rather than a major conformational driving force promoting this proximity. In particular, the conformational preferences of the hydroxymethyl group in aqueous environment appear to be dominated by 1,3-syn-diaxial repulsion, with gauche and solvation effects being secondary, and intramolecular hydrogen-bonding essentially negligible. The rotational dynamics of the exocyclic hydroxyl groups, which cannot be probed experimentally, is found to be rapid (10-100 ps timescale) and correlated (flip-flop hydrogen-bonds interconverting preferentially through an asynchronous disrotatory pathway). Structured solvent environments are observed between the ring and lactol oxygen atoms, as well as between the 4-OH and hydroxymethyl groups. The calculated stability differences between the four compounds are dominated by intramolecular effects, while the corresponding differences in solvation free energies are small. An inversion of the stereochemistry at either C(2) or C(4) from equatorial to axial is associated with a raise in free energy. Finally, the particularly low hydrophilicity of beta-D-talose appears to be caused by the formation of a high-occurrence hydrogen-bonded bridge between the 1,3-syn-diaxial 2-OH and 4-OH groups. Overall, good agreement is found with available experimental and theoretical data on the structural, dynamical, solvation and energetic properties of these compounds. However, this detailed comparison also reveals some discrepancies, suggesting the need (and providing a solid basis) for further refinement.

Free energy calculations using flexible-constrained, hard-constrained and non-constrained molecular dynamics simulations

A comparison of different treatments of bond-stretching interactions in molecular dynamics simulation is presented. Relative free energies from simulations using rigid bonds maintained with the SHAKE algorithm, using partially rigid bonds maintained with a recently introduced flexible constraints algorithm, and using fully flexible bonds are compared in a multi-configurational thermodynamic integration calculation of changing liquid water into liquid methanol. The formula for the free energy change due to a changing flexible constraint in a flexible constraint simulation is derived. To allow for a more direct comparison between these three methods, three different pairs of models for water and methanol were used: a flexible model (simulated without constraints and with flexible constraints), a rigid model (simulated with standard hard constraints), and an alternative flexible model (simulated with flexible constraints and standard hard constraints) in which the ideal or constrained bond lengths correspond to the average bond lengths obtained from a short simulation of the unconstrained flexible model. The particular treatment of the bonds induces differences of up to 2 % in the liquid densities, whereas (excess) free energy differences of up to 5.7 (4.3) kJ mol(-1) are observed. These values are smaller than the differences observed between the three different pairs of methanol/water models: up to 5 % in density and up to 8.5 kJ mol(-1) in (excess) free energy.

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.

On searching in, sampling of, and dynamically moving through conformational space of biomolecular systems: A review

Methods to search for low-energy conformations, to generate a Boltzmann-weighted ensemble of configurations, or to generate classical-dynamical trajectories for molecular systems in the condensed liquid phase are briefly reviewed with an eye to application to biomolecular systems. After having chosen the degrees of freedom and method to generate molecular configurations, the efficiency of the search or sampling can be enhanced in various ways: (i) efficient calculation of the energy function and forces, (ii) application of a plethora of search enhancement techniques, (iii) use of a biasing potential energy term, and (iv) guiding the sampling using a reaction or transition pathway. The overview of the available methods should help the reader to choose the combination that is most suitable for the biomolecular system, degrees of freedom, interaction function, and molecular or thermodynamic properties of interest.