Acceptance Ratio Method Research Articles

The Variationally Derived Intermediates Method to Calculate Free Energy Differences

Accurate and fast free energy calculations, based on atomistic simulations, are key to a first principles understanding of many-particle systems in biophysics, and enable quantitative descriptions of molecular recognition processes, drug binding, transmembrane transport, or functional conformational motions of proteins or complexes. To this aim, alchemical methods, such as the Free Energy Perturbation or the Bennett Acceptance Ratio method, require ‘morphing’ between the Hamiltonians of the involved molecular end states, typically in terms of a linear interpolation. However, linear transformations and soft-core variants thereof are still a very special case amongst all possible transformations. Here, we present the Variationally Derived Intermediates (VI) method that generalizes this approach and yields, under the assumption of uncorrelated sampling points, the intermediates with optimal sampling accuracy, which differ markedly from the established ones. We present applications to the challenging case of solvation free energies, yielding improved accuracies for decoupling electrostatic interactions, and similar ones for Lennard-Jones interactions, compared to state of the art soft-core transformations.

Computing Absolute Free Energy with Deep Generative Models.

Fast and accurate evaluation of free energy has broad applications from drug design to material engineering. Computing the absolute free energy is of particular interest since it allows the assessment of the relative stability between states without intermediates. Here, we introduce a general framework for calculating the absolute free energy of a state. A key step of the calculation is the definition of a reference state with tractable deep generative models using locally sampled configurations. The absolute free energy of this reference state is zero by design. The free energy for the state of interest can then be determined as the difference from the reference. We applied this approach to both discrete and continuous systems and demonstrated its effectiveness. It was found that the Bennett acceptance ratio method provides more accurate and efficient free energy estimations than approximate expressions based on work. We anticipate the method presented here to be a valuable strategy for computing free energy differences.

Determining Free-Energy Differences Through Variationally Derived Intermediates.

Free-energy calculations based on atomistic Hamiltonians and sampling are key to a first-principles understanding of biomolecular processes, material properties, and macromolecular chemistry. Here, we generalize the free-energy perturbation method and derive nonlinear Hamiltonian transformation sequences yielding free-energy estimates with minimal mean squared error with respect to the exact values. Our variational approach applies to finite sampling and holds for any finite number of intermediate states. We show that our sequences are also optimal for the Bennett acceptance ratio (BAR) method, thereby generalizing BAR to small sampling sizes and non-Gaussian error distributions.

Quantifying Protein-Protein Interactions in Molecular Simulations.

Interactions among proteins, nucleic acids, and other macromolecules are essential for their biological functions and shape the physicochemcial properties of the crowded environments inside living cells. Binding interactions are commonly quantified by dissociation constants Kd, and both binding and nonbinding interactions are quantified by second osmotic virial coefficients B2. As a measure of nonspecific binding and stickiness, B2 is receiving renewed attention in the context of so-called liquid–liquid phase separation in protein and nucleic acid solutions. We show that Kd is fully determined by B2 and the fraction of the dimer observed in molecular simulations of two proteins in a box. We derive two methods to calculate B2. From molecular dynamics or Monte Carlo simulations using implicit solvents, we can determine B2 from insertion and removal energies by applying Bennett’s acceptance ratio (BAR) method or the (binless) weighted histogram analysis method (WHAM). From simulations using implicit or explicit solvents, one can estimate B2 from the probability that the two molecules are within a volume large enough to cover their range of interactions. We validate these methods for coarse-grained Monte Carlo simulations of three weakly binding proteins. Our estimates for Kd and B2 allow us to separate out the contributions of nonbinding interactions to B2. Comparison of calculated and measured values of Kd and B2 can be used to (re-)parameterize and improve molecular force fields by calibrating specific affinities, overall stickiness, and nonbinding interactions. The accuracy and efficiency of Kd and B2 calculations make them well suited for high-throughput studies of large interactomes.

Identifying and estimating bias in overlap-sampling free-energy calculations

ABSTRACT We examine several methods for detecting or quantifying bias in free energy calculations performed using Bennett's acceptance ratio (BAR) method for combining data from molecular simulations of two systems. The metrics for bias detection studied here are one based on the emergence, at the BAR estimate, of a zero-slope inflection of the free energy with respect to its optimisation parameter, and another that makes use of calculation of the dissipated work. Bias estimation methods studied here are two approaches based on consideration of the maximum possible contribution from a configuration that may be rarely sampled, and another that appeared recently in the literature that is derived explicitly using the neglected-tail ansatz. All approaches are tested as a function of sampling amount in application to a model multiharmonic system. The dissipated-work method provides the more reliable indicator of the presence of bias, and the maximal sample techniques yield the tightest bounds on the bias, among the methods tested.

An Insight into Solubility of H2S in Choline Based Ionic Liquids from Molecular Dynamics Simulation

Solvation free energy of six choline based ILs were computed under molecular dynamics (MD) simulation using Bennet Acceptance Ratio (BAR) method. H2S removal from the natural gas can be accomplished by absorption process using a solvent. Degree elimination of gases in liquid can be estimated by excess chemical potential (μex), Henry’s law constant (kH), enthalpy (hex), and entropy (sex). The Optimised Potential Liquids Simulation (OPLS) force fields used for ILs were validated against experimental density in order to get valid results of solvation free energy. The percentage error of the predicted density of six choline based ILs were below than 2.0% which is a good agreement with the previous study. The highest value of kH was 109.23 atm gained by choline thiocyanate [Chl][tcy] and the lowest value of kH was 23.25 atm obtained by choline acetate [Chl][act]. There are many features that influence the solubility of H2S for example length of alkyl chain, free volume effect, molar density, type of anions and hydrogen bonding interaction. In order to understand the behaviour of H2S molecules in ILs system, the radial distribution function (RDF) of H2S-ILs and RDF for specific atom were analysed. Basically, the shorter the alkyl chain, increase in molar density of ILs, reduce the free volume thus provide less space for H2S to occupy and decrease insolubility.

Building a Macro-mixing Dual-basin Go Model using the Multistate Bennett Acceptance Ratio

Molecular dynamics (MD) simulations of biomolecules are widely used to investigate conformational dynamics and structure-function relationship. All-atom (AA) models provide the most accurate description of the underlying dynamics. However at present, even with the fastest computers, the time-scale obtainable with an atomistic simulation is up to a millisecond, whereas biologically relevant motions occur at the time scale of milliseconds to seconds. To overcome this, coarse-grained (CG) modeling can be utilized. The use of CG models reduces the computational time by several orders of magnitude, allowing access to biologically relevant time-scales. In this work, we use the dual-basin Gō-model, which is a structural-based CG model, for simulating conformational transitions between two known structures of a protein. The dual-basin model is formed by mixing two single-basin potentials and includes system-dependent parameters which determine the barrier height and the energetic offsets between the two basins. The determination of parameters however, is usually not straightforward and can be time consuming. Here, we have developed an efficient scheme to determine the mixing parameters using the Multistate Bennett Acceptance Ratio (MBAR) method after short simulations with a set of parameters. In the scheme, MBAR allows us to predict observables at various unsimulated conditions, which are useful to improve the mixing parameters in the next round of iterative simulations. The number of iterations that are necessary for obtaining the converged mixing parameters is significantly reduced in the scheme. We applied the scheme to several proteins, for showing the effectiveness in parameter determination. After obtaining the converged parameters, the proteins show frequent conformational transitions between open and closed states, providing the theoretical basis to investigate structure-dynamics-function relationships of the proteins.

Insight into structural properties of polyethylene glycol monolaurate in water and alcohols from molecular dynamics studies.

By means of molecular dynamics (MD) simulations, we explored the structural properties of polyethylene glycol monolaurate (PEGML) in water and in various aliphatic alcohols (methanol, ethanol, 2-propanol, 2-butanol, tert-butanol, and 1-pentanol). The PEGML and the alcohols were simulated using the optimized potentials for liquid simulations, all-atom (OPLS-AA) force field and water using the extended simple point charge (SPC/E) model. From the isothermal-isobaric (NPT, constant number of particles, constant pressure, and constant temperature) ensemble, we extracted the densities from the simulations and compared them with those from experimental results in order to confirm the validity of the selected force fields. The densities from MD simulations are in good agreement with the experimental values. To gain more insight into the nature of interactions between the PEGML and the solvent molecules, we analyzed the hydrogen-bonds, the electrostatic (Coulomb) interactions, and the van der Waals (Lennard-Jones) interaction energies extracted from MD simulations. The results were further strengthened by computing the solvation free energy by employing the free energy perturbation (FEP) approach. In this method, the free energy difference was computed by using the Bennet Acceptance Ratio (BAR) method. Moreover, the radial distribution functions were analyzed in order to gain more understanding of the solution behavior at the molecular level.

Algorithms for estimating the partition function of restricted Boltzmann machines

Accurate estimates of the normalization constants (partition functions) of energy-based probabilistic models (Markov random fields) are highly important, for example, for assessing the performance of models, monitoring training progress, and conducting likelihood ratio tests. Several algorithms for estimating the partition function (in relation to a reference distribution) have been introduced, including Annealed Importance Sampling (AIS) and Bennett's Acceptance Ratio method (BAR). However, their conceptual similarities and differences have not been worked out so far and systematic comparisons of their behavior in practice have been missing. We devise a unifying theoretical framework for these algorithms, which comprises existing variants and suggests new approaches. It is based on a generalized form of Crooks' equality linking the expectation over a distribution of samples generated by a transition operator to the expectation over the distribution induced by the reversed operator. The framework covers different ways of generating samples, such as parallel tempering and path sampling. An empirical comparison revealed the differences between the methods when estimating the partition function of restricted Boltzmann machines and Ising models. In our experiments, BAR using parallel tempering worked well with a small number of bridging distributions, while path sampling based AIS performed best when many bridging distributions were available. Because BAR gave the overall best results, we favor it over AIS. Furthermore, the experiments showed the importance of choosing a proper reference distribution.

Replica-Exchange Umbrella Sampling Combined with Gaussian Accelerated Molecular Dynamics for Free-Energy Calculation of Biomolecules.

We have developed an enhanced conformational sampling method combining replica-exchange umbrella sampling (REUS) with Gaussian accelerated molecular dynamics (GaMD). REUS enhances the sampling along predefined reaction coordinates, while GaMD accelerates the conformational dynamics by adding a boost potential to the system energy. The method, which we call GaREUS (Gaussian accelerated replica-exchange umbrella sampling), enhances the sampling more efficiently than REUS or GaMD, while the computational resource for GaREUS is the same as that required for REUS. The two-step reweighting procedure using the multistate Bennett acceptance ratio method and the cumulant expansion for the exponential average is applied to the simulation trajectories for obtaining the unbiased free-energy landscapes. We apply GaREUS to the calculations of free-energy landscapes for three different cases: conformational equilibria of N-glycan, folding of chignolin, and conformational change of adenyl kinase. We show that GaREUS speeds up the convergences of free-energy calculations using the same amount of computational resources as REUS. The free-energy landscapes reweighted from the trajectories of GaREUS agree with previously reported ones. GaREUS is applicable to free-energy calculations of various biomolecular dynamics and functions with reasonable computational costs.

Development of a Robust Indirect Approach for MM → QM Free Energy Calculations That Combines Force-Matched Reference Potential and Bennett's Acceptance Ratio Methods.

We use the PBE0/6-31G* density functional method to perform ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations under periodic boundary conditions with rigorous electrostatics using the ambient potential composite Ewald method in order to test the convergence of MM → QM/MM free energy corrections for the prediction of 17 small-molecule solvation free energies and eight ligand binding free energies to T4 lysozyme. The "indirect" thermodynamic cycle for calculating free energies is used to explore whether a series of reference potentials improve the statistical quality of the predictions. Specifically, we construct a series of reference potentials that optimize a molecular mechanical (MM) force field's parameters to reproduce the ab initio QM/MM forces from a QM/MM simulation. The optimizations form a systematic progression of successively expanded parameters that include bond, angle, dihedral, and charge parameters. For each reference potential, we calculate benchmark quality reference values for the MM → QM/MM correction by performing the mixed MM and QM/MM Hamiltonians at 11 intermediate states, each for 200 ps. We then compare forward and reverse application of Zwanzig's relation, thermodynamic integration (TI), and Bennett's acceptance ratio (BAR) methods as a function of reference potential, simulation time, and the number of simulated intermediate states. We find that Zwanzig's equation is inadequate unless a large number of intermediate states are explicitly simulated. The TI and BAR mean signed errors are very small even when only the end-state simulations are considered, and the standard deviations of the TI and BAR errors are decreased by choosing a reference potential that optimizes the bond and angle parameters. We find a robust approach for the data sets of fairly rigid molecules considered here is to use bond + angle reference potential together with the end-state-only BAR analysis. This requires QM/MM simulations to be performed in order to generate reference data to parametrize the bond + angle reference potential, and then this same simulation serves a dual purpose as the full QM/MM end state. The convergence of the results with respect to time suggests that computational resources may be used more efficiently by running multiple simulations for no more than 50 ps, rather than running one long simulation.

In silico solvation free energy and thermodynamics properties of H2S in cholinium-based amino acid ionic liquids

Solvation free energy of H2S in five cholinium-based amino acid ILs was computed using Bennet Acceptance Ratio (BAR) method under molecular dynamics (MD) simulation. From the predicted free energy, degree of removal of acid gases in liquids system can be measured and quantified by looking at excess chemical potential (μex) and Henry's law constant (kH). In order to obtain accurate value of predicted free energy, the optimised potential for liquid simulation (OPLS) force field used was validated against experimental density. The predicted density of five cholinium-based amino acid ILs showed a great agreement with the experimental results showing percentage error of <3.0%. The lowest value of Henry's law constant (kH) of H2S was obtained by cholinium phenylalanate [Chl][PHY] (11.20 atm) while the highest value was 44.12 atm gained by cholinium serinate [Chl][SRI]. To support our understanding, the radial distribution function (RDF) of pure ILs and RDF for a specific atom of H2S in these ILs was evaluated. As observed, there were other factors that can affect the solubility of H2S in these ILs such as the length of alkyl chain of anion, interaction energy between cation and anion as well as the hydrogen bonding interaction between polar sites of ILs with H2S molecules. Significantly, the longer hydrocarbon chain of anion, the weaker the interaction energy between the two ionic species. This in turn increases the free volume between the two ions and consequently promotes an excellent solubility of H2S in ILs. Furthermore, strong hydrogen bond interaction between the H2S and ILs molecules also contributes to high solubility of H2S in this ionic environment.

QresFEP: An Automated Protocol for Free Energy Calculations of Protein Mutations in Q.

Predicting the effect of single-point mutations on protein stability or protein-ligand binding is a major challenge in computational biology. Free energy calculations constitute the most rigorous approach to this problem, though the estimation of converged values for amino acid mutations remains challenging. To overcome this limitation, we developed tailored protocols to calculate free energy shifts associated with single-point mutations. We herein describe the QresFEP protocol, which includes an extension of our recent protocols to cover all amino acids mutations, based on the latest versions of the OPLS-AA force field. QresFEP is implemented in an application programming interface framework and the graphic interface QGui, for the molecular dynamics software Q. The complete protocol is benchmarked in several model systems, optimizing a number of sampling parameters and the implementation of Zwanzig's exponential formula and Bennet's acceptance ratio methods. QresFEP shows an excellent performance on estimating the hydration free energies of amino acid side-chain mimics, including their charged analogues. We also examined its performance on a protein-ligand binding problem of pharmaceutical relevance, the antagonism of neuropeptide Y1 G protein-coupled receptor. Here, the calculations show very good agreement with the experimental effect of 16 mutations on the binding of antagonists BIBP3226, in line with our recent applications in this field. Finally, the characterization of 43 mutations of T4-lysozyme reveals the capacity of our protocol to assess variations of the thermal stability of proteins, achieving a similar performance to alternative free energy perturbation (FEP) approaches. In summary, QresFEP is a robust, versatile, and user-friendly computational FEP protocol to examine biochemical effects of single-point mutations with high accuracy.

Q-Canonical Monte Carlo Sampling for Modeling the Linkage between Charge Regulation and Conformational Equilibria of Peptides.

The overall charge content and the patterning of charged residues have a profound impact on the conformational ensembles adopted by intrinsically disordered proteins. These parameters can be altered by charge regulation, which refers to the effects of post-translational modifications, pH-dependent changes to charge, and conformational fluctuations that modify the pKa values of ionizable residues. Although atomistic simulations have played a prominent role in uncovering the major sequence-ensemble relationships of IDPs, most simulations assume fixed charge states for ionizable residues. This may lead to erroneous estimates for conformational equilibria if they are linked to charge regulation. Here, we report the development of a new method we term q-canonical Monte Carlo sampling for modeling the linkage between charge regulation and conformational equilibria. The method, which is designed to be interoperable with the ABSINTH implicit solvation model, operates as follows: For a protein sequence with n ionizable residues, we start with all 2n charge microstates and use a criterion based on model compound pKa values to prune down to a subset of thermodynamically relevant charge microstates. This subset is then grouped into mesostates, where all microstates that belong to a mesostate have the same net charge. Conformational distributions, drawn from a canonical ensemble, are generated for each of the charge microstates that make up a mesostate using a method we designate as proton walk sampling. This method combines Metropolis Monte Carlo sampling in conformational space with an auxiliary Markov process that enables interconversions between charge microstates along a mesostate. Proton walk sampling helps identify the most likely charge microstate per mesostate. We then use thermodynamic integration aided by the multistate Bennett acceptance ratio method to estimate the free energies for converting between mesostates. These free energies are then combined with the per-microstate weights along each mesostate to estimate standard state free energies and pH-dependent free energies for all thermodynamically relevant charge microstates. The results provide quantitative estimates of the probabilities and preferred conformations associated with every thermodynamically accessible charge microstate. We showcase the application of q-canonical sampling using two model systems. The results establish the soundness of the method and the importance of charge regulation in systems characterized by conformational heterogeneity.

Free energy of formation of SPC/E-water and TIP4P-water using BAR and TI in MD and MC

Accurate free energy difference calculation based on the molecular dynamics (MD) and Monte-Carlo (MC) simulations is still a challenge in many different areas of science. In this paper, the free energy difference of the SPC/E-water and the TIP4P-water between two subsequent systems will be calculated using the thermodynamic integration (TI) and the Bennett acceptance ratio (BAR) methods. We compare the accuracy and efficiency of the two methods by calculating the free energy difference at three different temperatures, namely, T = 220, 240, 260, 280 and 300 K. The simulated results show that the free energy difference of small clusters scale nicely to the inverse of TCT−1, with the corresponding right critical temperature, TC, of SPC/E-water of 625 K and of TIP4P-water of 588 K.

Calculating free energies from the vibrational density of states function: Validation and critical assessment.

We explore and show the usefulness of the density of states function for computing vibrational free energies and free energy differences between small systems. Therefore, we compare this density of states integration method (DSI) to more established schemes such as Bennett's Acceptance Ratio method (BAR), the Normal Mode Analysis (NMA), and the Quasiharmonic Analysis (QHA). The strengths and shortcomings of all methods are highlighted with three numerical examples. Furthermore, the free energy of the ionization of ammonia and the mutation from serine to cysteine are computed using extensive ab initio molecular dynamics simulations. We conclude that DSI improves upon the other frequency-based methods (NMA and QHA) regarding the treatment of anharmonicity and yielding results comparable to BAR in all cases without the need for alchemical transformations. Low-frequency modes lead to larger errors indicating that long simulation times might be required for larger systems. In addition, we introduce the use of DSI for the localization of the vibrational free energy to specific atoms or residues, leading to insights into the underlying process, a unique feature that is only offered by this method.

All-atom molecular dynamics study on the non-solvent induced phase separation: Thermodynamics of adding water to poly(vinylidene fluoride)/N-methyl-2-pyrrolidone solution.

The change in the thermodynamics when adding water in poly(vinylidene fluoride) (PVDF)/N-methyl-2-pyrrolidone (NMP) solution is studied from all atom molecular dynamics (MD) simulations. This is done by estimating the free energy of mixing of PVDF/NMP solution with increasing volume fraction of water (ϕw) using an appropriately chosen thermodynamic cycle and the Bennett acceptance ratio method. The MD calculations predict the thermodynamic phase separation point of water/NMP/PVDF to be at ϕw = 0.08, in close agreement with the experimental cloud point measurement (ϕw = 0.05). Examining the enthalpic and entropic components of the free energy of mixing reveals that at low concentrations of water, the enthalpy term has the most significant contribution to the miscibility of the ternary system, whereas at higher concentrations of water, the entropy term dominates. Finally, the free energy of mixing was compared with the Flory-Huggins (FH) free energy of mixing by computing the concentration-dependent interaction parameters from MD simulations. The FH model inadequately predicted the miscibility of the PVDF solution, mainly due to its negligence of the excess entropy of mixing.

Hydroxycinnamic acids in supercritical carbon dioxide. The dependence of cosolvent-induced solubility enhancement on the selective solvation

Simulation of hydroxycinnamic acids (caffeic, p-coumaric, ferulic, synapic) in pure and modified by polar cosolvents (methanol, acetone) supercritical carbon dioxide has been carried out using classical molecular dynamics. Dependences of solvation free energies of the acids on the pressure and cosolvent mole fraction have been calculated using Bennett’s acceptance ratio method. When a polar cosolvent is added, the solvation free energy decreases sharply, and with a subsequent increase in its concentration, it changes slightly. The cosolvent-induced effect decreases with increasing pressure. The structural parameters characterizing the selective solvation of the acids have been calculated. A correlation between the cosolvent-induced solubility enhancement and the average number of solute – cosolvent hydrogen bonds has been found.

Determining Free Energy Differences Through Non-Linear Morphing

The Bennett Acceptance Ratio (BAR) method is widely used to estimate the free energy difference between two given states of a many-body system, such as a macromolecule, where Molecular Dynamics or Monte Carlo simulations are typically used for sampling. Because insufficient phase space overlap of the two states can lead to large biases in the results, the calculation is commonly split up into several smaller steps using intermediate states defined along a so-called morphing path — typically a linear interpolation between the start and the end states' Hamiltonians. Here, using BAR, we develop an optimization scheme to find the path that yields the most accurate free energy values from the much more general class of paths that includes also non-linear interpolations. For a broad range of one-dimensional test systems, we show that the free energy sampling error can be decreased by 25% to 35% compared to linear interpolations, with a theoretical upper limit of 50%. Further, correlations arising from using the same sample in one state to evaluate differences to more than one other state can be taken into account. Finally, we demonstrate that the optimization is still valid for a large number of intermediate states with very sparse sampling for each, enabling further optimization by trading off sampling versus number of intermediate states.

Introducing a Novel Multi-Level Method for Simulating the pH Dependence of Charge State Fluctuations and Conformational Ensembles of Intrinsically Disordered Proteins

Simulations of intrinsically disordered proteins (IDPs) use fixed charge models where the charge states of ionizable residues are set using model compound pKa values. Recent experiments have questioned the validity of fixed charge models. For a 100-residue IDP without sequence biases there can be roughly 103 or more unique possible charge states. Accurate and fully ergodic simulations require converged conformational sampling of each of the charge states. To facilitate this we have developed a multi-level Monte Carlo and free energy based method that leverages the ABSINTH implicit solvation model to enable efficient simulations of the conformational ensembles of all possible charge states for a given sequence. This is achieved by grouping charge states based on their net charge and adding auxiliary Markov chains that enable protonation/deprotonation of ionizable residues. At a given pH, the populations of thermodynamically accessible charge states and conformations are calculated using the Bennett Acceptance Ratio method and data from simulations obtained for all possible charge states. Key features of the ABSINTH model allow reductions in the number of replicas needed for free energy calculations. We also estimate shifts in free energies upon protonation / deprotonation relative to model compounds. Using these estimates, we quantify charge state fluctuations - known as charge regulation - to estimate pKa shifts within IDPs and extract the pH dependencies of the populations of charge states and conformational ensembles. In peptides excised from IDP sequences we find that local sequence contexts cause shifts in pKa values that can be larger than one pH unit. The method, with the incorporation of heuristics to enable increased efficiencies, is being deployed for simulations of charge states and conformational ensembles of a spectrum of IDPs.