Bennett Acceptance Ratio Method Research Articles

Prediction of cosolvent effect on solvation free energies and solubilities of organic compounds in supercritical carbon dioxide based on fully atomistic molecular simulations.

The solubility of organic compounds in supercritical fluids can be dramatically affected by addition of a suitable cosolvent (entrainer) at small concentrations. This makes the screening of the best-suited cosolvent an important task for the supercritical technology. The present study aims to improve our fundamental understanding of solvation in supercritical CO2 with cosolvents. We address the following questions: (1) How does the solvation free energy depend on the chemical class of an organic solute and the chemical nature of co-solvents? (2) Which intermolecular interactions determine the effect of a cosolvent on the solubility of organic compounds? We performed extensive calculations of solvation free energies of monofunctional organic molecules at infinite dilution in supercritical media by the Bennett's acceptance ratio method based on fully atomistic molecular dynamics sampling. Sixteen monofunctional organic molecules were solvated in pure sc-CO2 and sc-CO2 with addition of 6 molar % of cosolvents of different chemical nature: ethanol, acetone, and n-hexane. Cosolvent-induced solubility enhancement (CISE) factors were also calculated. It was found that formation of significant number of hydrogen bonds between a solute and cosolvent molecules leads to a profound solubility enhancement. The cosolvent effect is proportional to the number of hydrogen bonds. When polar cosolvents do not form hydrogen bonds with solutes, the CISE correlates with the dipole moment of solute molecules. However, the electrostatic interactions have a small impact on the solubility enhancement compared to hydrogen bonding. Addition of a nonpolar cosolvent, n-hexane, has a very little effect on the solvation Gibbs free energy of studied small organic molecules. The observed trends were discussed in line with available experimental data.

Defining the free-energy landscape of curvature-inducing proteins on membrane bilayers.

Curvature-sensing and curvature-remodeling proteins, such as Amphiphysin, Epsin, and Exo70, are known to reshape cell membranes, and this remodeling event is essential for key biophysical processes such as tubulation, exocytosis, and endocytosis. Curvature-inducing proteins can act as curvature sensors; they aggregate to membrane regions matching their intrinsic curvature; as well as induce curvature in cell membranes to stabilize emergent high curvature, nonspherical, structures such as tubules, discs, and caveolae. A definitive understanding of the interplay between protein recruitment and migration, the evolution of membrane curvature, and membrane morphological transitions is emerging but remains incomplete. Here, within a continuum framework and using the machinery of Monte Carlo simulations, we introduce and compare three free-energy methods to delineate the free-energy landscape of curvature-inducing proteins on bilayer membranes. We demonstrate the utility of the Widom test particle (or field) insertion methodology in computing the excess chemical potentials associated with curvature-inducing proteins on the membrane-in particular, we use this method to track the onset of morphological transitions in the membrane at elevated protein densities. We validate this approach by comparing the results from the Widom method with those of thermodynamic integration and Bennett acceptance ratio methods. Furthermore, the predictions from the Widom method have been tested against analytical calculations of the excess chemical potential at infinite dilution. Our results are useful in precisely quantifying the free-energy landscape, and also in determining the phase boundaries associated with curvature-induction, curvature-sensing, and morphological transitions. This approach can be extended to studies exploring the role of thermal fluctuations and other external (control) variables, such as membrane excess area, in shaping curvature-mediated interactions on bilayer membranes.

Hydration Free Energy from Orthogonal Space Random Walk and Polarizable Force Field.

The orthogonal space random walk (OSRW) method has shown enhanced sampling efficiency in free energy calculations from previous studies. In this study, the implementation of OSRW in accordance with the polarizable AMOEBA force field in TINKER molecular modeling software package is discussed and subsequently applied to the hydration free energy calculation of 20 small organic molecules, among which 15 are positively charged and five are neutral. The calculated hydration free energies of these molecules are compared with the results obtained from the Bennett acceptance ratio method using the same force field, and overall an excellent agreement is obtained. The convergence and the efficiency of the OSRW are also discussed and compared with BAR. Combining enhanced sampling techniques such as OSRW with polarizable force fields is very promising for achieving both accuracy and efficiency in general free energy calculations.

Simulation of rare events in quantum error correction

We consider the problem of calculating the logical error probability for a stabilizer quantum code subject to random Pauli errors. To access the regime of large code distances where logical errors are extremely unlikely we adopt the splitting method widely used in Monte Carlo simulations of rare events and Bennett's acceptance ratio method for estimating the free energy difference between two canonical ensembles. To illustrate the power of these methods in the context of error correction, we calculate the logical error probability $P_L$ for the 2D surface code on a square lattice with a pair of holes for all code distances $d\le 20$ and all error rates $p$ below the fault-tolerance threshold. Our numerical results confirm the expected exponential decay $P_L\sim \exp{[-\alpha(p)d]}$ and provide a simple fitting formula for the decay rate $\alpha(p)$. Both noiseless and noisy syndrome readout circuits are considered.

Nonequilibrium work by charge control in a Josephson junction

We consider a single Josephson junction in the presence of time varying gate charge, and examine the nonequilibrium work done by the charge control in the framework of fluctuation theorems. Assuming first a high quality junction with negligible Ohmic current, we obtain the probability distribution functions of the work and confirm the Crooks relation to give the estimation of the free energy changes ΔF=0. The reliability of ΔF estimated from the Jarzynksi equality is crucially dependent on protocol parameters, while the Bennett's acceptance ratio method yields consistently ΔF=0. We examine the behaviors of the work average and point out its relation to heat and entropy production associated with the circuit control. Finally considering finite tunnel resistance we discuss dissipation effects on the work statistics.

A comparative study of methods to compute the free energy of an ordered assembly by molecular simulation

We present a comparative study of methods to compute the absolute free energy of a crystalline assembly of hard particles by molecular simulation. We consider all combinations of three choices defining the methodology: (1) the reference system: Einstein crystal (EC), interacting harmonic (IH), or r(-12) soft spheres (SS); (2) the integration path: Frenkel-Ladd (FL) or penetrable ramp (PR); and (3) the free-energy method: overlap-sampling free-energy perturbation (OS) or thermodynamic integration (TI). We apply the methods to FCC hard spheres at the melting state. The study shows that, in the best cases, OS and TI are roughly equivalent in efficiency, with a slight advantage to TI. We also examine the multistate Bennett acceptance ratio method, and find that it offers no advantage for this particular application. The PR path shows advantage in general over FL, providing results of the same precision with 2-9 times less computation, depending on the choice of a common reference. The best combination for the FL path is TI+EC, which is how the FL method is usually implemented. For the PR path, the SS system (with either TI or OS) proves to be most effective; it gives equivalent precision to TI+FL+EC with about 6 times less computation (or 12 times less, if discounting the computational effort required to establish the SS reference free energy). Both the SS and IH references show great advantage in capturing finite-size effects, providing a variation in free-energy difference with system size that is about 10 times less than EC. This result further confirms previous work for soft-particle crystals, and suggests that free-energy calculations for a structured assembly be performed using a hybrid method, in which the finite-system free-energy difference is added to the extrapolated (1/N→0) absolute free energy of the reference system, to obtain a result that is nearly independent of system size.

High performance computing for drug development on K computer

Massively parallel computations (MP-CAFEE) ware developed to calculate absolute binding free energies of small molecules bound to a protein by all-atom molecular dynamics. It uses the nonequilibrium work measurement and Bennett acceptance ratio methods to calculate the free energy difference between the bound and unbound states. The FUJI force field was developed in order to assign force field parameters to arbitrary organic molecules in a unified manner including proteins and nucleic acids. Its dihedral parameters agree with the torsion energy profiles calculated by high-level ab initio molecular orbital theory for the model systems of protein backbone. Comparing with various force fields it agrees well with recent observations by vibrational spectroscopy on Ramachandran angle's population of alanine dipeptide in water. MP-CAFEE with FUJI force field has an efficient parallel algorithm and enough accuracy for computer aided drug design.

Comparison of thermodynamic integration and Bennett acceptance ratio for calculating relative protein‐ligand binding free energies

The performances of Bennett's acceptance ratio method and thermodynamic integration (TI) for the calculation of free energy differences in protein simulations are compared. For the latter, the standard trapezoidal rule, Simpson's rule, and Clenshaw-Curtis integration are used as numerical integration methods. We evaluate the influence of the number and definition of intermediate states on the precision, accuracy, and efficiency of the free energy calculations. Our results show that non-equidistantly spaced intermediate states are in some cases beneficial for the TI methods. Using several combinations of softness parameters and the λ power dependence, it is shown that these benefits are strongly dependent on the shape of the integrand. Although TI is more user-friendly due to its simplicity, it was found that Bennett's acceptance ratio method is the more efficient method. It is also the least dependent on the choice of the intermediate states, making it more robust than TI.

Perturbation method to calculate the density of states

Monte Carlo switching moves ("perturbations") are defined between two or more classical Hamiltonians sharing a common ground-state energy. The ratio of the density of states (DOS) of one system to that of another is related to the ensemble averages of the microcanonical acceptance probabilities of switching between these Hamiltonians, analogously to the case of Bennett's acceptance ratio method for the canonical ensemble [C. H. Bennett, J. Comput. Phys. 22, 245 (1976)]. Thus, if the DOS of one of the systems is known, one obtains those of the others and, hence, the partition functions. As a simple test case, the vapor pressure of an anharmonic Einstein crystal is computed, using the harmonic Einstein crystal as the reference system in one dimension; an auxiliary calculation is also performed in three dimensions. As a further example of the algorithm, the energy dependence of the ratio of the DOS of the square-well and hard-sphere tetradecamers is determined, from which the temperature dependence of the constant-volume heat capacity of the square-well system is calculated and compared with canonical Metropolis Monte Carlo estimates. For these cases and reference systems, the perturbation calculations exhibit a higher degree of convergence per Monte Carlo cycle than Wang-Landau (WL) sampling, although for the one-dimensional oscillator the WL sampling is ultimately more efficient for long runs. Last, we calculate the vapor pressure of liquid gold using an empirical Sutton-Chen many-body potential and the ideal gas as the reference state. Although this proves the general applicability of the method, by its inherent perturbation approach the algorithm is suitable for those particular cases where the properties of a related system are well known.

Estimating time-correlation functions by sampling and unbiasing dynamically activated events

Transition path sampling is a rare-event method that estimates state-to-state time-correlation functions in many-body systems from samples of short trajectories. In this framework, it is proposed to bias the importance function using the lowest Jacobian eigenvalue moduli along the dynamical trajectory. A lowest eigenvalue modulus is related to the lowest eigenvalue of the Hessian matrix and is evaluated here using the Lanczos algorithm as in activation-relaxation techniques. This results in favoring the sampling of activated trajectories and enhancing the occurrence of the rare reactive trajectories of interest, those corresponding to transitions between locally stable states. Estimating the time-correlation functions involves unbiasing the sample of simulated trajectories which is done using the multi-state Bennett acceptance ratio (MBAR) method. To assess the performance of our procedure, we compute the time-correlation function associated with the migration of a vacancy in α-iron. The derivative of the estimated time-correlation function yields a migration rate in agreement with the one given by transition state theory. Besides, we show that the information relative to rejected trajectories can be recycled within MBAR, resulting in a substantial speed-up. Unlike original transition path-sampling, our approach does not require computing the reversible work to confine the trajectory endpoints to a reactive state.

Comparison of free-energy estimators and their dependence on dissipated work

The estimate of free energy changes based on Bennett's acceptance ratio method is examined in several limiting cases and compared with other estimates based on the Jarzynski equality and on the Crooks relation. While the absolute amount of the dissipated work, defined as the surplus of the average work over the free energy difference, limits the practical applicability of Jarzynski's and Crooks' methods, the reliability of Bennett's approach is restricted by the difference of the dissipated works in the forward and the backward processes. We illustrate these points by considering a Gaussian chain and a hairpin chain which both are extended during the forward and accordingly compressed during the backward protocols. The reliability of the Crooks relation predominantly depends on the sample size; for the Jarzynski estimator the slowness of the work protocol is crucial, and the Bennett method is shown to give precise estimates irrespective of the pulling speed and sample size as long as the dissipated works are the same for the forward and the backward processes as it is the case for Gaussian work distributions. With an increasing dissipated work difference the Bennett estimator also acquires a bias which increases roughly in proportion to this difference. A substantial simplification of the Bennett estimator is provided by the 1/2 formula which expresses the free energy difference by the algebraic average of the Jarzynski estimates for the forward and the backward processes. It agrees with the Bennett estimate in all cases when the Jarzynski and the Crooks estimates fail to give reliable results.

Minimum Free Energy Path of Ligand-Induced Transition in Adenylate Kinase

Large-scale conformational changes in proteins involve barrier-crossing transitions on the complex free energy surfaces of high-dimensional space. Such rare events cannot be efficiently captured by conventional molecular dynamics simulations. Here we show that, by combining the on-the-fly string method and the multi-state Bennett acceptance ratio (MBAR) method, the free energy profile of a conformational transition pathway in Escherichia coli adenylate kinase can be characterized in a high-dimensional space. The minimum free energy paths of the conformational transitions in adenylate kinase were explored by the on-the-fly string method in 20-dimensional space spanned by the 20 largest-amplitude principal modes, and the free energy and various kinds of average physical quantities along the pathways were successfully evaluated by the MBAR method. The influence of ligand binding on the pathways was characterized in terms of rigid-body motions of the lid-shaped ATP-binding domain (LID) and the AMP-binding (AMPbd) domains. It was found that the LID domain was able to partially close without the ligand, while the closure of the AMPbd domain required the ligand binding. The transition state ensemble of the ligand bound form was identified as those structures characterized by highly specific binding of the ligand to the AMPbd domain, and was validated by unrestrained MD simulations. It was also found that complete closure of the LID domain required the dehydration of solvents around the P-loop. These findings suggest that the interplay of the two different types of domain motion is an essential feature in the conformational transition of the enzyme.

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.

Free-energy molecular simulations of the inclusion complex of Ne with fullerene C60 in water

Binding free energy of Ne in the endohedral fullerene Ne@C60 in water was calculated using closed thermodynamic cycle that involves annihilation in silico of Ne in the free and bound states. The underlying microscopic force field molecular dynamics simulations were performed at 300 K in the constant-volume 27.9 Å periodic box of the TIP3P water molecules. The reference binding free energy of a dummy atom to C60 was calculated from the ideal gas theory equations. The free energies of atom annihilation were computed by means of the free-energy perturbation and Bennett acceptance ratio methods. Introduction of intermediate particle of small size in the mutation scheme in C60 substantially improved free-energy convergence. The predicted binding free energies of Ne and Ne cavity to C60 at the standard state concentration are about − 5 kJ mol− 1 and +22 kJ mol− 1, respectively. The free energy of reversible annihilation of Ne is controlled by the insertion step for simulations both in water and in C60.

Improving the Efficiency of Protein-Ligand Binding Free-Energy Calculations by System Truncation.

We have studied whether the efficiency of alchemical free-energy calculations with the Bennett acceptance ratio method of protein-ligand binding energies can be improved by simulating only part of the protein. To this end, we solvated the full protein in a spherical droplet with a radius of 46 Å, surrounded by a vacuum. Then, we systematically reduced the size of the droplet and at the same time ignored protein residues that were outside the droplet. Radii of 40-15 Å were tested. Ten inhibitors of the blood clotting factor Xa were studied, and the results were compared to an earlier study in which the protein was solvated in a periodic box, showing complete agreement between the two sets of calculations within statistical uncertainty. We then show that the simulated system can be truncated down to 15 Å, without changing the calculated affinities by more than 0.5 kJ/mol on average (maximum difference of 1.4 kJ/mol). Moreover, we show that reducing the number of intermediate states in the calculations from eleven to three gave deviations that, on average, were only 0.5 kJ/mol (maximum of 1.4 kJ/mol). Together, these results show that truncation is an appropriate way to improve the efficiency of free-energy calculations for small mutations that preserve the net charge of the ligand. In fact, each calculation of a relative binding affinity requires only six simulations, each of which takes ∼15 CPU h of computation on a single processor.

Conformational Transition Pathways of Adenylate Kinase Explored by the String Method

Large-scale conformational changes in proteins involve barrier-crossing transitions on the complex free energy surfaces of high-dimensional space. Here we show that, by combining the on-the-fly string method and the multi-state Bennett acceptance ratio (MBAR) method, the free energy profile of a conformational transition pathway in Escherichia coli adenylate kinase can be accurately determined. The minimum free energy paths (MFEPs) of adenylate kinase were searched by the on-the-fly string method in 20-dimensional space spanned by the 20 largest-amplitude principal modes, and the free energy and various kinds of average physical quantities along the MFEPs were successfully evaluated by the MBAR method. The influence of ligand binding on the MFEPs was characterized in terms of rigid-body motions of the LID and AMPbd domains. It was found that the LID domain was able to partially close without the ligand, while the closure of the AMPbd domain required the ligand binding. The transition state ensemble of the ligand bound form was identified as those structures characterized by highly specific binding of the ligand to the AMPbd domain, and was validated by unrestrained MD simulations. It was also found that complete closure of the LID domain required the dehydration of solvents around the P-loop. These findings suggest that the interplay of the two different types of domain motion is an feature in the conformational transition of the enzyme.

Predicting binding affinities of host-guest systems in the SAMPL3 blind challenge: the performance of relative free energy calculations

Relative free energy calculations based on molecular dynamics simulations are combined with available experimental binding free energies to predict unknown binding affinities of acyclic Cucurbituril complexes in the blind SAMPL3 competition. The predictions yield root mean square errors between 2.6 and 3.2 kcal/mol for seven host-guest systems. Those deviations are comparable to results for solvation free energies of small organic molecules. However, the standard deviations found in our simulations range from 0.4 to 2.4 kcal/mol, which indicates the need for better sampling. Three different approaches are compared. Bennett's Acceptance Ratio Method and thermodynamic integration based on the trapezoidal rule with 12 λ-points exhibit a root mean square error of 2.6 kcal/mol, while thermodynamic integration with Simpson's rule and 11 λ-points leads to a root mean square error of 3.2 kcal/mol. In terms of absolute median errors, Bennett's Acceptance Ratio Method performs better than thermodynamic integration with the trapezoidal rule (1.7 vs. 2.9 kcal/mol). Simulations of the deprotonated forms of the guest molecules exhibit a poorer correspondence to experimental results with a root mean square error of 5.2 kcal/mol. In addition, a decrease of the buffer concentration by approximately 20 mM in the simulations raises the root mean square error to 3.8 kcal/mol.

Small Molecule Solvation Free Energy: Enhanced Conformational Sampling Using Expanded Ensemble Molecular Dynamics Simulation

We present an efficient expanded ensemble molecular dynamics method to calculate the solvation free energy (or residual chemical potential) of small molecules with complex topologies. The methodology is validated by computing the solvation free energy of ibuprofen in water, methanol, and ethanol at 300 K and 1 bar and comparing to reference simulation results using Bennett's acceptance ratio method. Difficulties with ibuprofen using conventional molecular dynamics methods stem from an inadequate sampling of the carboxylic acid functional group, which, for the present study, is subject to free energy barriers of rotation of 14-20 kBT. While several advances have been made to overcome such weaknesses, we demonstrate how this shortcoming is easily overcome by using an expanded ensemble methodology to facilitate conformational sampling. Not only does the method enhance conformational sampling but it also boosts the rate of exploration of the configurational phase space and requires only a single simulation to calculate the solvation free energy. Agreement between the expanded ensemble and the reference calculations is good for all three solvents, with the reported uncertainties of the expanded ensemble being comparable to the uncertainties of the reference calculations, while requiring less simulation time; the reduced simulation time demonstrates the improved performance gained from the expanded ensemble method.

Avoiding the van der Waals endpoint problem using serial atomic insertion

In the past analyses of the so-called van der Waals end point problem focused on thermodynamic integration. Here we investigate which of the recommendations, such as the need for soft-core potentials, are still valid when Bennett's acceptance ratio method is used. We show that in combination with Bennett's acceptance ratio method intermediate states characterized by the coupling parameter λ can be replaced by intermediate states in which Lennard-Jones interactions are turned on or off on an "atom by atom" basis. By doing so, there is no necessity to use soft-core potentials. In fact, one can compute free energy differences without dedicated code, making it possible to use any molecular dynamics program to compute alchemical free energy differences. Such an approach, which we illustrate by several examples, makes it possible to exploit the tremendous computational power of the graphics processing unit.

Efficient Solvation Free Energy Calculations of Amino Acid Analogs by Expanded Ensemble Molecular Simulation.

We present an efficient, automated expanded ensemble method to calculate the residual chemical potential or solvation free energy by molecular dynamics simulation. The methodology is validated by computing the residual chemical potential of 13 amino acid analogs in water at 300 K and 1 bar and comparing to reference simulation data. Overall agreement is good, with the methodology of the present study reaching limiting precisions of less than 0.1 kBT in half of the total simulation time of the reference simulation study which utilized Bennett's acceptance ratio method. The apparent difference in the efficiencies is a result of the inherent advantages of the expanded ensemble method, which creates an improved decorrelation of simulation data and improves the sampling of the important regions of the configurational phase space of each subensemble. The present adaptation utilizes histograms of proposed transition energies collected throughout the entire simulation, to make extremely precise calculations of the relative free energy between neighboring subensembles.