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.
I. Kim appreciates Professor Alexei Stuchebrukhov for his introduction to Jarzynski work-fluctuation theorem and to the shifted harmonic oscillator model used in this paper. We acknowledge National Science Foundation (NSF) Award Nos. MCB-0546768 and MCB-1052477 for funding.
I. INTRODUCTION II. METHODS A. Umbrella sampling B. Bennett acceptance ratio and histogram analysis enhanced by US III. NUMERICAL TESTS A. Hummer 1D asymmetric potential B. A strongly (but linearly) shifted harmonic oscillator IV. APPLICATIONS A. The electrostatic contribution to the excess chemical potential of water B. Ionic selectivity of the S2 site in the KcsA potassium channel V. CONCLUSION