Accurate Free Energy Calculations Research Articles

Grid-Based Surface Generalized Born Model for Calculation of Electrostatic Binding Free Energies

Fast and accurate calculation of solvation free energies is central to many applications, such as rational drug design. In this study, we present a grid-based molecular surface implementation of "R6" flavor of the generalized Born (GB) implicit solvent model, named GBNSR6. The speed, accuracy relative to numerical Poisson-Boltzmann treatment, and sensitivity to grid surface parameters are tested on a set of 15 small protein-ligand complexes and a set of biomolecules in the range of 268 to 25099 atoms. Our results demonstrate that the proposed model provides a relatively successful compromise between the speed and accuracy of computing polar components of the solvation free energies (ΔGpol) and binding free energies (ΔΔGpol). The model tolerates a relatively coarse grid size h = 0.5 Å, where the grid artifact error in computing ΔΔGpol remains in the range of kBT ∼ 0.6 kcal/mol. The estimated ΔΔGpols are well correlated (r2 = 0.97) with the numerical Poisson-Boltzmann reference, while showing virtually no systematic bias and RMSE = 1.43 kcal/mol. The grid-based GBNSR6 model is available in Amber (AmberTools) package of molecular simulation programs.

Binding free energy predictions of farnesoid X receptor (FXR) agonists using a linear interaction energy (LIE) approach with reliability estimation: application to the D3R Grand Challenge 2

Computational protein binding affinity prediction can play an important role in drug research but performing efficient and accurate binding free energy calculations is still challenging. In the context of phase 2 of the Drug Design Data Resource (D3R) Grand Challenge 2 we used our automated eTOX ALLIES approach to apply the (iterative) linear interaction energy (LIE) method and we evaluated its performance in predicting binding affinities for farnesoid X receptor (FXR) agonists. Efficiency was obtained by our pre-calibrated LIE models and molecular dynamics (MD) simulations at the nanosecond scale, while predictive accuracy was obtained for a small subset of compounds. Using our recently introduced reliability estimation metrics, we could classify predictions with higher confidence by featuring an applicability domain (AD) analysis in combination with protein–ligand interaction profiling. The outcomes of and agreement between our AD and interaction-profile analyses to distinguish and rationalize the performance of our predictions highlighted the relevance of sufficiently exploring protein–ligand interactions during training and it demonstrated the possibility to quantitatively and efficiently evaluate if this is achieved by using simulation data only.

Using force-based adaptive resolution simulations to calculate solvation free energies of amino acid sidechain analogues.

The calculation of free energy differences is a crucial step in the characterization and understanding of the physical properties of biological molecules. In the development of efficient methods to compute these quantities, a promising strategy is that of employing a dual-resolution representation of the solvent, specifically using an accurate model in the proximity of a molecule of interest and a simplified description elsewhere. One such concurrent multi-resolution simulation method is the Adaptive Resolution Scheme (AdResS), in which particles smoothly change their resolution on-the-fly as they move between different subregions. Before using this approach in the context of free energy calculations, however, it is necessary to make sure that the dual-resolution treatment of the solvent does not cause undesired effects on the computed quantities. Here, we show how AdResS can be used to calculate solvation free energies of small polar solutes using Thermodynamic Integration (TI). We discuss how the potential-energy-based TI approach combines with the force-based AdResS methodology, in which no global Hamiltonian is defined. The AdResS free energy values agree with those calculated from fully atomistic simulations to within a fraction of kBT. This is true even for small atomistic regions whose size is on the order of the correlation length, or when the properties of the coarse-grained region are extremely different from those of the atomistic region. These accurate free energy calculations are possible because AdResS allows the sampling of solvation shell configurations which are equivalent to those of fully atomistic simulations. The results of the present work thus demonstrate the viability of the use of adaptive resolution simulation methods to perform free energy calculations and pave the way for large-scale applications where a substantial computational gain can be attained.

The Uhlenbeck-Ford model: Exact virial coefficients and application as a reference system in fluid-phase free-energy calculations.

The Uhlenbeck-Ford (UF) model was originally proposed for the theoretical study of imperfect gases, given that all its virial coefficients can be evaluated exactly, in principle. Here, in addition to computing the previously unknown coefficients B11 through B13, we assess its applicability as a reference system in fluid-phase free-energy calculations using molecular simulation techniques. Our results demonstrate that, although the UF model itself is too soft, appropriately scaled Uhlenbeck-Ford (sUF) models provide robust reference systems that allow accurate fluid-phase free-energy calculations without the need for an intermediate reference model. Indeed, in addition to the accuracy with which their free energies are known and their convenient scaling properties, the fluid is the only thermodynamically stable phase for a wide range of sUF models. This set of favorable properties may potentially put the sUF fluid-phase reference systems on par with the standard role that harmonic and Einstein solids play as reference systems for solid-phase free-energy calculations.

An efficient protocol for obtaining accurate hydration free energies using quantum chemistry and reweighting from molecular dynamics simulations

The non-Boltzmann Bennett (NBB) free energy estimator method is applied to 21 molecules from the blind subset of the SAMPL4 challenge. When NBB is applied with the SMD implicit solvent model, and the OLYP/DZP level of quantum chemistry, highly accurate hydration free energy calculations are obtained with respect to experiment (RMSD=0.89kcal·mol−1). Other quantum chemical methods are also tested, and the effects of solvent model, density functional, basis set are explored in this benchmarking study, providing a framework for improvements in calculating hydration free energies. We provide a practical guide for using the best QM-NBB protocols that are consistently more accurate than either pure QM or pure MM alone. In situations where high accuracy hydration free energy predictions are needed, the QM-NBB method with SMD implicit solvent should be the first choice of quantum chemists.

Extended Adaptive Biasing Force Algorithm. An On-the-Fly Implementation for Accurate Free-Energy Calculations

Proper use of the adaptive biasing force (ABF) algorithm in free-energy calculations needs certain prerequisites to be met, namely, that the Jacobian for the metric transformation and its first derivative be available and the coarse variables be independent and fully decoupled from any holonomic constraint or geometric restraint, thereby limiting singularly the field of application of the approach. The extended ABF (eABF) algorithm circumvents these intrinsic limitations by applying the time-dependent bias onto a fictitious particle coupled to the coarse variable of interest by means of a stiff spring. However, with the current implementation of eABF in the popular molecular dynamics engine NAMD, a trajectory-based post-treatment is necessary to derive the underlying free-energy change. Usually, such a posthoc analysis leads to a decrease in the reliability of the free-energy estimates due to the inevitable loss of information, as well as to a drop in efficiency, which stems from substantial read-write accesses to file systems. We have developed a user-friendly, on-the-fly code for performing eABF simulations within NAMD. In the present contribution, this code is probed in eight illustrative examples. The performance of the algorithm is compared with traditional ABF, on the one hand, and the original eABF implementation combined with a posthoc analysis, on the other hand. Our results indicate that the on-the-fly eABF algorithm (i) supplies the correct free-energy landscape in those critical cases where the coarse variables at play are coupled to either each other or to geometric restraints or holonomic constraints, (ii) greatly improves the reliability of the free-energy change, compared to the outcome of a posthoc analysis, and (iii) represents a negligible additional computational effort compared to regular ABF. Moreover, in the proposed implementation, guidelines for choosing two parameters of the eABF algorithm, namely the stiffness of the spring and the mass of the fictitious particles, are proposed. The present on-the-fly eABF implementation can be viewed as the second generation of the ABF algorithm, expected to be widely utilized in the theoretical investigation of recognition and association phenomena relevant to physics, chemistry, and biology.

Unconstrained enhanced sampling for free energy calculations of biomolecules: a review

Free energy calculations are central to understanding the structure, dynamics and function of biomolecules. Yet insufficient sampling of biomolecular configurations is often regarded as one of the main sources of error. Many enhanced sampling techniques have been developed to address this issue. Notably, enhanced sampling methods based on biasing collective variables (CVs), including the widely used umbrella sampling, adaptive biasing force and metadynamics, have been discussed in a recent excellent review (Abrams and Bussi, Entropy, 2014). Here, we aim to review enhanced sampling methods that do not require predefined system-dependent CVs for biomolecular simulations and as such do not suffer from the hidden energy barrier problem as encountered in the CV-biasing methods. These methods include, but are not limited to, replica exchange/parallel tempering, self-guided molecular/Langevin dynamics, essential energy space random walk and accelerated molecular dynamics. While it is overwhelming to describe all details of each method, we provide a summary of the methods along with the applications and offer our perspectives. We conclude with challenges and prospects of the unconstrained enhanced sampling methods for accurate biomolecular free energy calculations.

Ab-Initio MD Simulations of Redox Reactions of Liquid Electrolytes and SEI Formation

Atomistic understanding of solid electrolyte interphase (SEI) formed at the interface between electrode and liquid electrolyte is still an issue of greatest importance. However, difficulties in experimental in-situ observations and simulations considering the liquid dynamics (fluctuation) and reaction have inhibited elucidation of physics and chemistry of SEIs. We have addressed such unresolved issues by means of accurate large-scale ab-initio molecular dynamics (AIMD) samplings and free energy calculations on supercomputers. The SEI formation is initiated by electrochemical decomposition of the electrolyte molecules such as solvents and additives. Slight change of the electrolyte component often exhibits a large impact on the SEI quality. We first examined the reductive decomposition of the electrolyte molecules at charging (reductive) environment, using a typical electrolyte with ethylene carbonate (EC) solvents and vinylene carbonate (VC) additives. We carried out AIMD free energy calculations for possible decomposition pathways, and found that VC passivation to EC anion radical is a more appropriate mechanism than the conventional scenario with sacrificial reductive decomposition of VC (Fig. (a))1. This novel mechanism can explain the decrease of irreversible capacity and the CO2 evolution with a small amount of VC additives observed experimentally. Next, we extracted stable organic SEI film components (SFCs) from the products of the reductive decompositions of the electrolyte. We then investigated the characteristic of the probable SFC amorphous aggregates on the adhesion to the graphite electrode and the interfacial electronic states. The results showed that the SFC aggregates have characters of “unstable adhesion” to the graphite surface and “high electronic insulation”. With these results, we discussed that a near-shore aggregation mechanism of SFCs is promising for the SEI film formation.2 Very recently, we examined inorganic SFCs as well to elucidate their roles in the SEI formation.3 These new findings with supercomputer simulations give important perspectives for more efficient and reliable batteries. 1. K. Ushirogata, K. Sodeyama, Y. Okuno, Y. Tateyama, J. Am. Chem. Soc. 135, 11967-11974 (2013). 2. K. Ushirogata, K. Sodeyama, Z. Futera, Y. Tateyama, Y. Okuno, J. Electrochem. Soc. 162, A2670-2678 (2015). 3. Y. Okuno, K. Ushirogata, K. Sodeyama, Y. Tateyama, submitted. Figure (a) Proposed role of VC additive in reductive decomposition of the EC-based electrolyte. (b) Schematic picture of the near-shore aggregation mechanism for SEI formation on a graphitic electrode. Figure 1

A Riemannian Framework for Tackling Large-Scale Conformational Changes of Proteins using All-ATOM Molecular Dynamics Simulations

Simulation studies of biomolecular processes often suffer from the curse of dimensionality both in sampling and in analysis. Dimensionality reduction, either using intuitive reaction coordinates or automated algorithms, is thus an inevitable part of advanced biomolecular simulation studies. Most free energy calculation methods and path-finding algorithms, implicitly or explicitly, rely on an effective diffusion model in a low-dimensional space as their underlying theoretical framework. Dimensionality reduction, however, comes with a price and the effective dynamics of the reduced system could be nontrivial. For instance, assuming the diffusivity of the effective dynamics, the diffusion tensor is likely to be position-dependent and anisotropic.We have developed a Riemannian framework for free energy calculation methods and path-finding algorithms in molecular dynamics simulations. We show that a Riemannian diffusion model may describe the dynamics of the reduced system such that a Riemannian metric replaces the position-dependent diffusion tensor. In addition, the Riemannian formulation provides a suitable mathematical framework for generalizing one-dimensional free energy relations to multi-dimensional spaces in a rigorous manner such that the model is invariant under any smooth coordinate transformation.Within the proposed mathematical framework, a novel combination of free energy calculation methods and path-finding algorithms is developed which can be used for accurate calculation of free energies along protein conformational transition pathways using all-atom molecular dynamics simulations. The method is particularly useful for describing large-scale and complex conformational changes of proteins. We note that the simplified Euclidean version of the method was recently employed successfully in reconstructing the transport cycle of a membrane transporter at an atomic level (see M. Moradi, et al, “Atomic-level characterization of transport cycle thermodynamics in the glycerol-3-phosphate:phosphate antiporter.”, Nat. Commun., 6:8393, 2015).

Calculation of molecular free energies in classical potentials

Free energies of molecules can be calculated by quantum chemistry computations or by normal mode classical calculations. However, the first can be computationally impractical for large molecules and the second is based on the assumption of harmonic dynamics. We present a novel, accurate and complete calculation of molecular free energies in standard classical potentials. In this method we transform the molecule by relaxing potential terms which depend on the coordinates of a group of atoms in that molecule and calculate the free energy difference associated with the transformation. Then, since the transformed molecule can be treated as non-interacting systems, the free energy associated with these atoms is analytically or numerically calculated. This two-step calculation can be applied to calculate free energies of molecules or free energy difference between (possibly large) molecules in a general environment. We demonstrate the method in free energy calculations for methanethiol and butane molecules in vacuum and solvent. We suggest the potential application of free energy calculation of chemical reactions in classical molecular simulations.

Hydration Free Energies of Molecular Ions from Theory and Simulation.

We present a theoretical/computational framework for accurate calculation of hydration free energies of ionized molecular species. The method is based on a molecular theory, 3D-RISM, combined with a recently developed pressure correction (PC+). The 3D-RISM/PC+ model can provide ∼3 kcal/mol hydration free energy accuracy for a large variety of ionic compounds, provided that the Galvani potential of water is taken into account. The results are compared with direct atomistic simulations. Several methodological aspects of hydration free energy calculations for charged species are discussed.

Accurate Calculation of Solvation Free Energies in Supercritical Fluids by Fully Atomistic Simulations: Probing the Theory of Solutions in Energy Representation.

Accurate calculation of solvation free energies (SFEs) is a fundamental problem of theoretical chemistry. In this work we perform a careful validation of the theory of solutions in energy representation (ER method) developed by Matubayasi et al. [J. Chem. Phys. 2000, 113, 6070-6081] for SFE calculations in supercritical solvents. This method can be seen as a bridge between the molecular simulations and the classical (not quantum) density functional theory (DFT) formulated in energy representation. We performed extensive calculations of SFEs of organic molecules of different chemical natures in pure supercritical CO2 (sc-CO2) and in sc-CO2 with addition of 6 mol % of ethanol, acetone, and n-hexane as cosolvents. We show that the ER method reproduces SFE data calculated by a method free of theoretical approximations (the Bennett's acceptance ratio) with the mean absolute error of only 0.05 kcal/mol. However, the ER method requires by an order less computational resources. Also, we show that the quality of ER calculations should be carefully monitored since the lack of sampling can result into a considerable bias in predictions. The present calculations reproduce the trends in the cosolvent-induced solubility enhancement factors observed in experimental data. Thus, we think that molecular simulations coupled with the ER method can be used for quick calculations of the effect of variation of temperature, pressure, and cosolvent concentration on SFE and hence solubility of bioactive compounds in supercritical fluids. This should dramatically reduce the burden of experimental work on optimizing solvency of supercritical solvents.

2P021 高精度自由エネルギー計算によるドラッグデザイン(01B. 蛋白質:構造機能相関,ポスター,第52回日本生物物理学会年会(2014年度))

2P021 高精度自由エネルギー計算によるドラッグデザイン(01B. 蛋白質:構造機能相関,ポスター,第52回日本生物物理学会年会(2014年度))

Control of KirBac3.1 Potassium Channel Gating at the Interface between Cytoplasmic Domains

KirBac channels are prokaryotic homologs of mammalian inwardly rectifying potassium (Kir) channels, and recent structures of KirBac3.1 have provided important insights into the structural basis of gating in Kir channels. In this study, we demonstrate that KirBac3.1 channel activity is strongly pH-dependent, and we used x-ray crystallography to determine the structural changes that arise from an activatory mutation (S205L) located in the cytoplasmic domain (CTD). This mutation stabilizes a novel energetically favorable open conformation in which changes at the intersubunit interface in the CTD also alter the electrostatic potential of the inner cytoplasmic cavity. These results provide a structural explanation for the activatory effect of this mutation and provide a greater insight into the role of the CTD in Kir channel gating.

First-principles calculation of entropy for liquid metals

We demonstrate the accurate calculation of entropies and free energies for a variety of liquid metals using an extension of the two-phase thermodynamic (2PT) model based on a decomposition of the velocity autocorrelation function into gas-like (hard sphere) and solid-like (harmonic) subsystems. The hard sphere model for the gas-like component is shown to give systematically high entropies for liquid metals as a direct result of the unphysical Lorentzian high-frequency tail. Using a memory function framework we derive a generally applicable velocity autocorrelation and frequency spectrum for the diffusive component which recovers the low-frequency (long-time) behavior of the hard sphere model while providing for realistic short-time coherence and high-frequency tails to the spectrum. This approach provides a significant increase in the accuracy of the calculated entropies for liquid metals and is compared to ambient pressure data for liquid sodium, aluminum, gallium, tin, and iron. The use of this method for the determination of melt boundaries is demonstrated with a calculation of the high-pressure bcc melt boundary for sodium. With the significantly improved accuracy available with the memory function treatment for softer interatomic potentials, the 2PT model for entropy calculations should find broader application in high energy density science, warm dense matter, planetary science, geophysics, and material science.

Temperature dependent effective potential method for accurate free energy calculations of solids

We have developed a thorough and accurate method of determining anharmonic free energies, the temperature dependent effective potential technique (TDEP). It is based on \emph{ab initio} molecular dynamics followed by a mapping onto a model Hamiltonian that describes the lattice dynamics. The formalism and the numerical aspects of the technique are described in details. A number of practical examples are given, and results are presented, which confirm the usefulness of TDEP within \emph{ab initio} and classical molecular dynamics frameworks. In particular, we examine from first-principles the behavior of force constants upon the dynamical stabilization of body centered phase of Zr, and show that they become more localized. We also calculate phase diagram for $^4$He modeled with the Aziz \emph{et al.} potential and obtain results which are in favorable agreement both with respect to experiment and established techniques.

Mechanism of Membrane Tube Formation Induced by Adhesive Nanocomponents

We report numerical simulations of membrane tubulation driven by large colloidal particles. Using Monte Carlo simulations we study how the process depends on particle size and binding strength, and present accurate free energy calculations to sort out how tube formation compares with the competing budding process. We find that tube formation is a result of the collective behavior of the particles adhering on the surface, and it occurs for binding strengths that are smaller than those required for budding. We also find that long linear aggregates of particles forming on the membrane surface act as nucleation seeds for tubulation by lowering the free energy barrier associated to the process.

Configurational Preferences of Arylamide α-Helix Mimetics via Alchemical Free Energy Calculations of Relative Binding Affinities

We use molecular docking and free energy calculations to estimate the relative free energy of binding of six arylamide compounds designed to inhibit the hDM2-p53 interaction. We show that using docking methods to predict or rank the binding affinity of a series of arylamide inhibitors of the hDM2-p53 interaction is problematic. However, using free energy calculations, we show that we can achieve levels of accuracy that can guide the development of novel arylamide compounds. We perform alchemical free energy calculations using the Desmond molecular dynamics package with the same arylamide inhibitors of the hDM2-p53 system and illustrate the challenges of performing accurate free energy calculations for realistic systems. To our knowledge, these are the first calculations for inhibitors of the hDM2 system that employ a full treatment of statistical mechanics including explicit water representation and full protein flexibility. We show that mutating three functional groups in a single transformation can be more efficient than mutating the groups one by one if proper intermediates are used. We also show that Hamiltonian exchanges can improve the efficiency of the calculation compared to standard alchemical methods, with a novel use of the phase space overlap to monitor sampling extent. We show that, despite sampling limitations, this approach can achieve levels of accuracy sufficient to bias further inhibitor modification toward binding, and identifies antiparallel configurations as stable or more stable than the parallel configurations that are typically considered.

Insights from Free-Energy Calculations: Protein Conformational Equilibrium, Driving Forces, and Ligand-Binding Modes

Accurate free-energy calculations provide mechanistic insights into molecular recognition and conformational equilibrium. In this work, we performed free-energy calculations to study the thermodynamic properties of different states of molecular systems in their equilibrium basin, and obtained accurate absolute binding free-energy calculations for protein-ligand binding using a newly developed M2 algorithm. We used a range of Asp-Phe-Gly (DFG)-in/out p38α mitogen-activated protein kinase inhibitors as our test cases. We also focused on the flexible DFG motif, which is closely connected to kinase activation and inhibitor binding. Our calculations explain the coexistence of DFG-in and DFG-out states of the loop and reveal different components (e.g., configurational entropy and enthalpy) that stabilize the apo p38α conformations. To study novel ligand-binding modes and the key driving forces behind them, we computed the absolute binding free energies of 30 p38α inhibitors, including analogs with unavailable experimental structures. The calculations revealed multiple stable, complex conformations and changes in p38α and inhibitor conformations, as well as balance in several energetic terms and configurational entropy loss. The results provide relevant physics that can aid in designing inhibitors and understanding protein conformational equilibrium. Our approach is fast for use with proteins that contain flexible regions for structure-based drug design.

New Soft-Core Potential Function for Molecular Dynamics Based Alchemical Free Energy Calculations.

The fields of rational drug design and protein engineering benefit from accurate free energy calculations based on molecular dynamics simulations. A thermodynamic integration scheme is often used to calculate changes in the free energy of a system by integrating the change of the system's Hamiltonian with respect to a coupling parameter. These methods exploit nonphysical pathways over thermodynamic cycles involving particle introduction and annihilation. Such alchemical transitions require the modification of the classical nonbonded potential energy terms by applying soft-core potential functions to avoid singularity points. In this work, we propose a novel formulation for a soft-core potential to be applied in nonequilibrium free energy calculations that alleviates singularities, numerical instabilities, and additional minima in the potential energy for all combinations of nonbonded interactions at all intermediate alchemical states. The method was validated by application to (a) the free energy calculations of a closed thermodynamic cycle, (b) the mutation influence on protein thermostability, (c) calculations of small ligand solvation free energies, and (d) the estimation of binding free energies of trypsin inhibitors. The results show that the novel soft-core function provides a robust and accurate general purpose solution to alchemical free energy calculations.