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.

Predicting absolute aqueous solubility by applying a machine learning model for an artificially liquid-state as proxy for the solid-state

In this study, we use machine learning algorithms with QM-derived COSMO-RS descriptors, along with Morgan fingerprints, to predict the absolute solubility of drug-like compounds. The QM-derived descriptors account for the molecular properties of the solute, i.e., the solute-solute interactions in an artificial-liquid-state (super-cooled liquid), and the solute-solvent interactions in solution. We employ two main approaches to predict solubility: (i) a hypothetical pathway that involves melting the solute at room temperature T = T¯ ([Formula: see text]) and mixing the artificially liquid solute into the solvent ([Formula: see text]). In this approach [Formula: see text] is predicted using machine learning models, and the [Formula: see text] is obtained from COSMO-RS calculations; (ii) direct solubility prediction using machine learning algorithms. The models were trained on a large number of Bayer in-house compounds for which water solubility data is available at physiological pH of 6.5 and ambient temperature. We also evaluated our models using external datasets from a solubility challenge. Our models present great improvements compared to the absolute solubility prediction with the QSAR model for the artificial liquid state as implemented in the COSMOtherm software, for both in-house and external datasets. We are furthermore able to demonstrate the superiority of QM-derived descriptors compared to cheminformatics descriptors. We finally present low-cost alternative models using fragment-based COSMOquick calculations with only marginal reduction in the quality of predicted solubility.

Efficient Alchemical Intermediate States in Free Energy Calculations Using λ-Enveloping Distribution Sampling

Alchemical free energy calculations generally require intermediate states along a coupling parameter λ to establish sufficient phase space overlap for obtaining converged results. Such intermediate states can also be engineered to lower the energy barriers and, consequently, reduce the required sampling time. The recently introduced λ-enveloping distribution sampling (λ-EDS) scheme combines the properties of the minimum variance pathway and the EDS methods to improve sampling and allow for larger steps along the alchemical pathway compared to conventional approaches. This scheme also eliminates the need for soft-core potentials and retains the behavior of conventional λ-intermediate states as a limiting case. In this study, an automated procedure is developed to select the parameters of λ-EDS for optimal performance. The underlying theory is illustrated based on simulations of simple test systems (bond length changes in harmonic oscillators, mutations of dihedral angles, and charge creation in water), as well as on the calculation of the absolute hydration free energies of 12 small organic molecules.

An Alternative to Conventional λ-Intermediate States in Alchemical Free Energy Calculations: λ-Enveloping Distribution Sampling

Alchemical free energy calculations typically rely on intermediate states to bridge between the relevant phase spaces of the two end states. These intermediate states are usually created by mixing the energies or parameters of the end states according to a coupling parameter λ. The choice of the procedure has a strong impact on the efficiency of the calculation, as it affects both the encountered energy barriers and the phase space overlap between the states. The present work builds on the connection between the minimum variance pathway (MVP) and enveloping distribution sampling (EDS). It is shown that both methods can be regarded as special cases of a common scheme referred to as λ-EDS, which can also reproduce the behavior of conventional λ-intermediate states. A particularly attractive feature of λ-EDS is its ability to emulate the use of soft core potentials (SCP) while avoiding the associated computational overhead when applying efficient free energy estimators such as the multistate Bennett's acceptance ratio (MBAR). The method is illustrated for both relative and absolute free energy calculations considering five benchmark systems. The first two systems (charge inversion and cavity creation in a dipolar solvent) demonstrate the use of λ-EDS as an alternative coupling scheme in the context of thermodynamic integration (TI). The three other systems (change of bond length, change of dihedral angles, and cavity creation in water) investigate the efficiency and optimal choice of parameters in the context of free energy perturbation (FEP) and Bennett's acceptance ratio (BAR). It is shown that λ-EDS allows larger steps along the alchemical pathway than conventional intermediate states.

Toward Automated Free Energy Calculation with Accelerated Enveloping Distribution Sampling (A-EDS)

Free-energy perturbation (FEP) methods are commonly used in drug design to calculate relative binding free energies of different ligands to a common host protein. Alchemical ligand transformations are usually performed in multiple steps which need to be chosen carefully to ensure sufficient phase-space overlap between neighboring states. With one-step or single-step FEP techniques, a single reference state is designed that samples phase-space not only representative of a full transformation but also ideally resembles multiple ligand end states and hence allows for efficient multistate perturbations. Enveloping distribution sampling (EDS) is one example for such a method in which the reference state is created by a mathematical combination of the different ligand end states based on solid statistical mechanics. We have recently proposed a novel approach to EDS which enables efficient barrier crossing between the different end states, termed accelerated EDS (A-EDS). In this work, we further simplify the parametrization of the A-EDS reference state and demonstrate the automated calculation of multiple free-energy differences between different ligands from a single simulation in three different well-described drug design model systems.

Alchemical Free-Energy Calculations by Multiple-Replica λ-Dynamics: The Conveyor Belt Thermodynamic Integration Scheme

A new method is proposed to calculate alchemical free-energy differences based on molecular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at different values of the alchemical coupling parameter λ. The number K is taken to be even, and the replicas are equally spaced on a forward-turn-backward-turn path, akin to a conveyor belt (CB) between the two physical end-states; and as in λ-dynamics (λD), the λ-values associated with the individual systems evolve in time along the simulation. However, they do so in a concerted fashion, determined by the evolution of a single dynamical variable Λ of period 2π controlling the advance of the entire CB. Thus, a change of Λ is always associated with K/2 equispaced replicas moving forward and K/2 equispaced replicas moving backward along λ. As a result, the effective free-energy profile of the replica system along Λ is periodic of period 2 πK^-1, and the magnitude of its variations decreases rapidly upon increasing K, at least as K^-1 in the limit of large K. When a sufficient number of replicas is used, these variations become small, which enables a complete and quasi-homogeneous coverage of the λ-range by the replica system, without application of any biasing potential. If desired, a memory-based biasing potential can still be added to further homogenize the sampling, the preoptimization of which is computationally inexpensive. The final free-energy profile along λ is calculated similarly to TI, by binning of the Hamiltonian λ-derivative as a function of λ considering all replicas simultaneously, followed by quadrature integration. The associated quadrature error can be kept very low owing to the continuous and quasi-homogeneous λ-sampling. The CBTI scheme can be viewed as a continuous/deterministic/dynamical analog of the Hamiltonian replica-exchange/permutation (HRE/HRP) schemes or as a correlated multiple-replica analog of the λD or λ-local elevation umbrella sampling (λ-LEUS) schemes. Compared to TI, it shares the advantage of the latter schemes in terms of enhanced orthogonal sampling, i.e. the availability of variable-λ paths to circumvent conformational barriers present at specific λ-values. Compared to HRE/HRP, it permits a deterministic and continuous sampling of the λ-range, is expected to be less sensitive to possible artifacts of the thermo- and barostating schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt frequency. Compared to λ-LEUS, it eliminates (or drastically reduces) the dead time associated with the preoptimization of a biasing potential. The goal of this article is to provide the mathematical/physical formulation of the proposed CBTI scheme, along with an initial application of the method to the calculation of the hydration free energy of methanol.

Comparison of free-energy methods using a tripeptide-water model system

We investigate the ability of several free-energy calculation methods to combine two alchemical changes. We use Bennett acceptance ratio (BAR), thermodynamic integration (TI), extended TI (X-TI), and enveloping distribution sampling (EDS) to perturb a water molecule, which is restrained to an amino acid that is also being perturbed. In addition to these pairwise methods, we present two two-dimensional approaches, EDS-TI and two-dimensional TI (2D-TI). We compare feasibility, efficiency and usability of these methods in regard to our simple model system, which mimics the displacement of a water molecule in the active site of a protein on residue mutation. The correct treatment of structural water has been shown to greatly aid binding affinity calculations in some cases that remained elusive otherwise. This is of broad interest in, for example, drug design, and we conclude that thus far, the pairwise method BAR and also the newer X-TI remain the most suitable methods to treat this problem as long as few end states are involved. © 2018 Wiley Periodicals, Inc.

Validation and Comparison of Force Fields for Native Cyclodextrins in Aqueous Solution

Molecular dynamics simulations of native α-, β-, and γ-cyclodextrin in aqueous solution have been conducted with the goal to investigate the performance of the CHARMM36 force field, the AMBER-compatible q4md-CD force field, and five variants of the GROMOS force field. The properties analyzed are structural parameters derived from X-ray diffraction and NMR experiments as well as hydrogen bonds and hydration patterns, including hydration free enthalpies. Recent revisions of the torsional-angle parameters for carbohydrate systems within the GROMOS family of force fields lead to a significant improvement of the agreement between simulated and experimental NMR data. Therefore, we recommend using the variant 53A6_GLYC instead of 53A6 and 56A6_{CARBO_R} or 2016H66 instead of 56A6_CARBO to simulate cyclodextrins in solution. The CHARMM36 and q4md-CD force fields show a similar performance as the three recommended GROMOS parameter sets. A significant difference is the more flexible nature of the cyclodextrins modeled with the CHARMM36 and q4md-CD force fields compared to the three recommended GROMOS parameter sets.

Efficient Round-Trip Time Optimization for Replica-Exchange Enveloping Distribution Sampling (RE-EDS)

Replica-exchange enveloping distribution sampling (RE-EDS) allows the efficient estimation of free-energy differences between multiple end-states from a single molecular dynamics (MD) simulation. In EDS, a reference state is sampled, which can be tuned by two types of parameters, i.e., smoothness parameters(s) and energy offsets, such that all end-states are sufficiently sampled. However, the choice of these parameters is not trivial. Replica exchange (RE) or parallel tempering is a widely applied technique to enhance sampling. By combining EDS with the RE technique, the parameter choice problem could be simplified and the challenge shifted toward an optimal distribution of the replicas in the smoothness-parameter space. The choice of a certain replica distribution can alter the sampling efficiency significantly. In this work, global round-trip time optimization (GRTO) algorithms are tested for the use in RE-EDS simulations. In addition, a local round-trip time optimization (LRTO) algorithm is proposed for systems with slowly adapting environments, where a reliable estimate for the round-trip time is challenging to obtain. The optimization algorithms were applied to RE-EDS simulations of a system of nine small-molecule inhibitors of phenylethanolamine N-methyltransferase (PNMT). The energy offsets were determined using our recently proposed parallel energy-offset (PEOE) estimation scheme. While the multistate GRTO algorithm yielded the best replica distribution for the ligands in water, the multistate LRTO algorithm was found to be the method of choice for the ligands in complex with PNMT. With this, the 36 alchemical free-energy differences between the nine ligands were calculated successfully from a single RE-EDS simulation 10 ns in length. Thus, RE-EDS presents an efficient method for the estimation of relative binding free energies.

Replica exchange enveloping distribution sampling (RE-EDS): A robust method to estimate multiple free-energy differences from a single simulation

In molecular dynamics (MD) simulations, free-energy differences are often calculated using free energy perturbation or thermodynamic integration (TI) methods. However, both techniques are only suited to calculate free-energy differences between two end states. Enveloping distribution sampling (EDS) presents an attractive alternative that allows to calculate multiple free-energy differences in a single simulation. In EDS, a reference state is simulated which "envelopes" the end states. The challenge of this methodology is the determination of optimal reference-state parameters to ensure equal sampling of all end states. Currently, the automatic determination of the reference-state parameters for multiple end states is an unsolved issue that limits the application of the methodology. To resolve this, we have generalised the replica-exchange EDS (RE-EDS) approach, introduced by Lee et al. [J. Chem. Theory Comput. 10, 2738 (2014)] for constant-pH MD simulations. By exchanging configurations between replicas with different reference-state parameters, the complexity of the parameter-choice problem can be substantially reduced. A new robust scheme to estimate the reference-state parameters from a short initial RE-EDS simulation with default parameters was developed, which allowed the calculation of 36 free-energy differences between nine small-molecule inhibitors of phenylethanolamine N-methyltransferase from a single simulation. The resulting free-energy differences were in excellent agreement with values obtained previously by TI and two-state EDS simulations.

A Virtual Mixture Approach to the Study of Multistate Equilibrium: Application to Constant pH Simulation in Explicit Water

Chemical and thermodynamic equilibrium of multiple states is a fundamental phenomenon in biology systems and has been the focus of many experimental and computational studies. This work presents a simulation method to directly study the equilibrium of multiple states. This method constructs a virtual mixture of multiple states (VMMS) to sample the conformational space of all chemical states simultaneously. The VMMS system consists of multiple subsystems, one for each state. The subsystem contains a solute and a solvent environment. The solute molecules in all subsystems share the same conformation but have their own solvent environments. Transition between states is implicated by the change of their molar fractions. Simulation of a VMMS system allows efficient calculation of relative free energies of all states, which in turn determine their equilibrium molar fractions. For systems with a large number of state transition sites, an implicit site approximation is introduced to minimize the cost of simulation. A direct application of the VMMS method is for constant pH simulation to study protonation equilibrium. Applying the VMMS method to a heptapeptide of 3 ionizable residues, we calculated the pKas of those residues both with all explicit states and with implicit sites and obtained consistent results. For mouse epidermal growth factor of 9 ionizable groups, our VMMS simulations with implicit sites produced pKas of all 9 ionizable groups and the results agree qualitatively with NMR measurement. This example demonstrates the VMMS method can be applied to systems of a large number of ionizable groups and the computational cost scales linearly with the number of ionizable groups. For one of the most challenging systems in constant pH calculation, SNase Δ+PHS/V66K, our VMMS simulation shows that it is the state-dependent water penetration that causes the large deviation in lysine66’s pKa.

Computer simulation plays an important role to understand molecular systems and has been applied to problems of increasing complexity. Multistate equilibrium is a fundamental concept behind the structure and function of biological systems. Due to the limit in computing resources and lack of good alternative methods, computer simulation has been conducted for systems in a single state, sampling from one state to another to infer equilibrium properties. This sequential approach has been successful in many cases such as protonation equilibrium with implicit solvation model. However, state transition is difficult when explicit solvent is used for more accurate solvation description. Many efforts have been dedicated to overcome this difficulty. Analogous to real multistate systems, we proposed a virtual mixture of multiple states (VMMS) to directly simulate the equilibrium. State transitions are replaced by changes in state molar fractions. Mimicking a test tube environment, all states are simulated in parallel to equilibrate with each other. Application to constant pH simulation in explicit water demonstrates the capability of this method. It is expected that the VMMS method will find more applications in biological problems related to the equilibrium of competing states.

In silico concurrent multisite pH titration in proteins

The concurrent proton binding at multiple sites in macromolecules such as proteins and nucleic acids is an important yet challenging problem in biochemistry. We develop an efficient generalized Hamiltonian approach to attack this issue. Based on the previously developed generalized-ensemble methods, an effective potential energy is constructed which combines the contributions of all (relevant) protonation states of the molecule. The effective potential preserves important phase regions of all states and, thus, allows efficient sampling of these regions in one simulation. The need for intermediate states in alchemical free energy simulations is greatly reduced. Free energy differences between different protonation states can be determined accurately and enable one to construct the grand canonical partition function. Therefore, the complicated concurrent multisite proton titration process of protein molecules can be satisfactorily simulated. Application of this method to the simulation of the pKa of Glu49, Asp50, and C-terminus of bovine pancreatic trypsin inhibitor shows reasonably good agreement with published experimental work. This method provides an unprecedented vivid picture of how different protonation states change their relative population upon pH titration. We believe that the method will be very useful in deciphering the molecular mechanism of pH-dependent biomolecular processes in terms of a detailed atomistic description.

Integrated Hamiltonian sampling: a simple and versatile method for free energy simulations and conformational sampling

Motivated by specific applications and the recent work of Gao and co-workers on integrated tempering sampling (ITS), we have developed a novel sampling approach referred to as integrated Hamiltonian sampling (IHS). IHS is straightforward to implement and complementary to existing methods for free energy simulation and enhanced configurational sampling. The method carries out sampling using an effective Hamiltonian constructed by integrating the Boltzmann distributions of a series of Hamiltonians. By judiciously selecting the weights of the different Hamiltonians, one achieves rapid transitions among the energy landscapes that underlie different Hamiltonians and therefore an efficient sampling of important regions of the conformational space. Along this line, IHS shares similar motivations as the enveloping distribution sampling (EDS) approach of van Gunsteren and co-workers, although the ways that distributions of different Hamiltonians are integrated are rather different in IHS and EDS. Specifically, we report efficient ways for determining the weights using a combination of histogram flattening and weighted histogram analysis approaches, which make it straightforward to include many end-state and intermediate Hamiltonians in IHS so as to enhance its flexibility. Using several relatively simple condensed phase examples, we illustrate the implementation and application of IHS as well as potential developments for the near future. The relation of IHS to several related sampling methods such as Hamiltonian replica exchange molecular dynamics and λ-dynamics is also briefly discussed.

Relative free enthalpies for point mutations in two proteins with highly similar sequences but different folds

Enveloping distribution sampling was used to calculate free-enthalpy changes associated with single amino acid mutations for a pair of proteins, GA95 and GB95, that show 95% sequence identity yet fold into topologically different structures. Of the L → A, I → F, and L → Y mutations at positions 20, 30, and 45, respectively, of the 56-residue sequence, the first and the last contribute the most to the free-enthalpy difference between the native and non-native sequence-structure combinations, in agreement with the experimental findings for this protein pair. The individual free-enthalpy changes are almost sequence-independent in the four-strand/one-helix structure, the stable form of GB95, while in the three-helix bundle structure, the stable form of GA95, an interplay between residues 20 and 45 is observed.

Efficient Combination of Environment Change and Alchemical Perturbation within the Enveloping Distribution Sampling (EDS) Scheme: Twin-System EDS and Application to the Determination of Octanol-Water Partition Coefficients

The methodology of Enveloping Distribution Sampling (EDS) is extended to probe a single-simulation alternative to the thermodynamic cycle that is standardly used for measuring the effect of a modification of a chemical compound, e.g. from a given species to a chemical derivative for a ligand or solute molecule, on the free-enthalpy change associated with a change in environment, e.g. from the unbound state to the bound state for a protein-ligand system or from one solvent to another one for a solute molecule. This alternative approach relies on the coupled simulation of two systems (computational boxes) 1 and 2, and the method is therefore referred to as twin-system EDS. Systems 1 and 2 account for the two choices of environment. The end states of the alchemical perturbation for the twin-system associate the two alternative forms X and Y of the molecule to systems 1 and 2 or 2 and 1, respectively. In this way, the processes of transforming one molecule into the other are carried out simultaneously in opposite directions in the two environments, leading to a change in free enthalpy that is smaller than for the two individual processes and to an energy-difference distribution that is more symmetric. As an illustration, the method is applied to the calculation of octanol-water partition coefficients for C4 to C8 alkanes, 1-hexanol and 1,2-dimethoxyethane. It is shown in particular that the consideration of the residual hydration of octanol leads to calculated partition coefficients that are in better agreement with reported experimental numbers.

Combination of Enveloping Distribution Sampling (EDS) of a Soft-Core Reference-State Hamiltonian with One-Step Perturbation to Predict the Effect of Side Chain Substitution on the Relative Stability of Right- and Left-Helical Folds of β-Peptides

Folding free enthalpies of many not too different polypeptides can be efficiently and accurately predicted with the one-step perturbation (OSP) method using only one or a few molecular dynamics (MD) simulations. In this article, we introduce a combination of enveloping distribution sampling (EDS) and the OSP method (EDS-OSP) and apply it to predict the free enthalpy differences between a right-handed 2.710/12-helix and a left-handed 314-helix for 16 β-peptides with slightly different side-chain substitution patterns. An EDS simulation of a designed soft-core reference-state peptide was carried out in which both helices were sampled. Then, the soft-core atoms were perturbed into physical atoms. Thus, free enthalpy differences between the two helices for the 16 β-peptides can be predicted from only one simulation. The results predicted by EDS-OSP and a previous OSP study are very similar, i.e., the deviations between the results of the 16 peptides are mostly within the order of kBT, and the average absolute deviation is 1.2 kJ mol(-1). Together with the EDS parameter update simulation, about 128 ns of MD simulations needed to be carried out using the EDS-OSP method, while 700 ns of MD simulations were required in the previous OSP study where two separate reference-state simulations and an additional long time MD simulation of one of the 16 β-peptides were carried out. Thus, the computational effort was significantly reduced, i.e., by more than a factor of 5, using the EDS-OSP method. Hence, we consider this method an efficient tool to predict conformational free enthalpy differences from MD simulations.

Free enthalpies of replacing water molecules in protein binding pockets

Water molecules in the binding pocket of a protein and their role in ligand binding have increasingly raised interest in recent years. Displacement of such water molecules by ligand atoms can be either favourable or unfavourable for ligand binding depending on the change in free enthalpy. In this study, we investigate the displacement of water molecules by an apolar probe in the binding pocket of two proteins, cyclin-dependent kinase 2 and tRNA-guanine transglycosylase, using the method of enveloping distribution sampling (EDS) to obtain free enthalpy differences. In both cases, a ligand core is placed inside the respective pocket and the remaining water molecules are converted to apolar probes, both individually and in pairs. The free enthalpy difference between a water molecule and a CH(3) group at the same location in the pocket in comparison to their presence in bulk solution calculated from EDS molecular dynamics simulations corresponds to the binding free enthalpy of CH(3) at this location. From the free enthalpy difference and the enthalpy difference, the entropic contribution of the displacement can be obtained too. The overlay of the resulting occupancy volumes of the water molecules with crystal structures of analogous ligands shows qualitative correlation between experimentally measured inhibition constants and the calculated free enthalpy differences. Thus, such an EDS analysis of the water molecules in the binding pocket may give valuable insight for potency optimization in drug design.

Bennett's acceptance ratio and histogram analysis methods enhanced by umbrella sampling along a reaction coordinate in configurational space

Free energy perturbation, a method for computing the free energy difference between two states, is often combined with non-Boltzmann biased sampling techniques in order to accelerate the convergence of free energy calculations. Here we present a new extension of the Bennett acceptance ratio (BAR) method by combining it with umbrella sampling (US) along a reaction coordinate in configurational space. In this approach, which we call Bennett acceptance ratio with umbrella sampling (BAR-US), the conditional histogram of energy difference (a mapping of the 3N-dimensional configurational space via a reaction coordinate onto 1D energy difference space) is weighted for marginalization with the associated population density along a reaction coordinate computed by US. This procedure produces marginal histograms of energy difference, from forward and backward simulations, with higher overlap in energy difference space, rendering free energy difference estimations using BAR statistically more reliable. In addition to BAR-US, two histogram analysis methods, termed Bennett overlapping histograms with US (BOH-US) and Bennett-Hummer (linear) least square with US (BHLS-US), are employed as consistency and convergence checks for free energy difference estimation by BAR-US. The proposed methods (BAR-US, BOH-US, and BHLS-US) are applied to a 1-dimensional asymmetric model potential, as has been used previously to test free energy calculations from non-equilibrium processes. We then consider the more stringent test of a 1-dimensional strongly (but linearly) shifted harmonic oscillator, which exhibits no overlap between two states when sampled using unbiased Brownian dynamics. We find that the efficiency of the proposed methods is enhanced over the original Bennett's methods (BAR, BOH, and BHLS) through fast uniform sampling of energy difference space via US in configurational space. We apply the proposed methods to the calculation of the electrostatic contribution to the absolute solvation free energy (excess chemical potential) of water. We then address the controversial issue of ion selectivity in the K(+) ion channel, KcsA. We have calculated the relative binding affinity of K(+) over Na(+) within a binding site of the KcsA channel for which different, though adjacent, K(+) and Na(+) configurations exist, ideally suited to these US-enhanced methods. Our studies demonstrate that the significant improvements in free energy calculations obtained using the proposed methods can have serious consequences for elucidating biological mechanisms and for the interpretation of experimental data.

Enhanced Sampling in Free Energy Calculations: Combining SGLD with the Bennett's Acceptance Ratio and Enveloping Distribution Sampling Methods

One of the key requirements for the accurate calculation of free energy differences is proper sampling of conformational space. Especially in biological applications, molecular dynamics simulations are often confronted with rugged energy surfaces and high energy barriers, leading to insufficient sampling and, in turn, poor convergence of the free energy results. In this work, we address this problem by employing enhanced sampling methods. We explore the possibility of using self-guided Langevin dynamics (SGLD) to speed up the exploration process in free energy simulations. To obtain improved free energy differences from such simulations, it is necessary to account for the effects of the bias due to the guiding forces. We demonstrate how this can be accomplished for the Bennett's acceptance ratio (BAR) and the enveloping distribution sampling (EDS) methods. While BAR is considered among the most efficient methods available for free energy calculations, the EDS method developed by Christ and van Gunsteren is a promising development that reduces the computational costs of free energy calculations by simulating a single reference state. To evaluate the accuracy of both approaches in connection with enhanced sampling, EDS was implemented in CHARMM. For testing, we employ benchmark systems with analytical reference results and the mutation of alanine to serine. We find that SGLD with reweighting can provide accurate results for BAR and EDS where conventional molecular dynamics simulations fail. In addition, we compare the performance of EDS with other free energy methods. We briefly discuss the implications of our results and provide practical guidelines for conducting free energy simulations with SGLD.

Comparison of enveloping distribution sampling and thermodynamic integration to calculate binding free energies of phenylethanolamine N-methyltransferase inhibitors

The relative binding free energy between two ligands to a specific protein can be obtained using various computational methods. The more accurate and also computationally more demanding techniques are the so-called free energy methods which use conformational sampling from molecular dynamics or Monte Carlo simulations to generate thermodynamic averages. Two such widely applied methods are the thermodynamic integration (TI) and the recently introduced enveloping distribution sampling (EDS) methods. In both cases relative binding free energies are obtained through the alchemical perturbations of one ligand into another in water and inside the binding pocket of the protein. TI requires many separate simulations and the specification of a pathway along which the system is perturbed from one ligand to another. Using the EDS approach, only a single automatically derived reference state enveloping both end states needs to be sampled. In addition, the choice of an optimal pathway in TI calculations is not trivial and a poor choice may lead to poor convergence along the pathway. Given this, EDS is expected to be a valuable and computationally efficient alternative to TI. In this study, the performances of these two methods are compared using the binding of ten tetrahydroisoquinoline derivatives to phenylethanolamine N-transferase as an example. The ligands involve a diverse set of functional groups leading to a wide range of free energy differences. In addition, two different schemes to determine automatically the EDS reference state parameters and two different topology approaches are compared.

Free energy calculations of protein-ligand interactions

In the calculation of free energies of binding for protein-ligand complexes, we distinguish endpoint methods, methods involving alchemical modifications and methods that physically displace the ligand from the protein. Most methodological advances seem to come from a clever combination of multiple existing methods to enhance the sampling or to utilize specific advantages of various approaches. The coupling parameters common in thermodynamic integration and in Hamiltonian replica exchange are for instance combined to yield replica exchange thermodynamic integration. As new methods mostly aim to improve efficiency or to attain more complete sampling, there are good prospects to understand and tackle the sampling problem better and to shift the focus towards the scoring problem in the context of more robust and accurate force fields.

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.