The Conveyor Belt Umbrella Sampling (CBUS) Scheme: Principle and Application to the Calculation of the Absolute Binding Free Energies of Alkali Cations to Crown Ethers.

The accurate prediction of absolute protein-ligand binding free energies is one of the grand challenge problems of computational science. Binding free energy measures the strength of binding between a ligand and a protein, and an algorithm that would allow its accurate prediction would be a powerful tool for rational drug design. Here we present the development of a new method that allows for the absolute binding free energy of a protein-ligand complex to be calculated from first principles, using a single simulation. Our method involves the use of a novel reaction coordinate that swaps a ligand bound to a protein with an equivalent volume of bulk water. This water-swap reaction coordinate is built using an identity constraint, which identifies a cluster of water molecules from bulk water that occupies the same volume as the ligand in the protein active site. A dual topology algorithm is then used to swap the ligand from the active site with the identified water cluster from bulk water. The free energy is then calculated using replica exchange thermodynamic integration. This returns the free energy change of simultaneously transferring the ligand to bulk water, as an equivalent volume of bulk water is transferred back to the protein active site. This, directly, is the absolute binding free energy. It should be noted that while this reaction coordinate models the binding process directly, an accurate force field and sufficient sampling are still required to allow for the binding free energy to be predicted correctly. In this paper we present the details and development of this method, and demonstrate how the potential of mean force along the water-swap coordinate can be improved by calibrating the soft-core Coulomb and Lennard-Jones parameters used for the dual topology calculation. The optimal parameters were applied to calculations of protein-ligand binding free energies of a neuraminidase inhibitor (oseltamivir), with these results compared to experiment. These results demonstrate that the water-swap coordinate provides a viable and potentially powerful new route for the prediction of protein-ligand binding free energies.

The standard (absolute) binding free energy of the antibiotic sparsomycin with the 50S bacteria ribosomal subunit is calculated using molecular dynamics (MD) free energy perturbation (FEP) simulations with restraining potentials developed by Wang et al. [Biophys. J. 91:2798-2814 (2006)]. In the simulation protocol, restraining potentials are activated for the orientational and translational movements of the ligand relative to the binding site when it is decoupled from the binding pocket, and then released once the ligand fully interacts with the rest of the system. A reduced system is simulated to decrease the computational cost of the FEP/MD calculations and the effects of the surrounding atoms are incorporated using the generalized solvent boundary potential (GSBP) method. The loss of conformational freedom of the ligand upon binding is characterized using the potential of mean force (PMF) as a function of the root-mean-square deviation (RMSD) relative to the bound conformation. The number of water molecules in the binding pocket is allowed to fluctuate dynamically in response to the ligand during the calculations by combining FEP/MD with grand canonical Monte Carlo (GCMC) simulations. The calculated binding free energy is about -6 kcal/mol, which is in reasonable agreement with the experimental value. The information gleaned from this study provides new insight on the recognition of ribosome by sparsomycin and highlights the challenges in calculations of absolute binding free energies in these systems.

We introduce a novel method called restrain-free energy perturbation-release (R-FEP-R) to estimate conformational free energy changes via an alchemical path, which for some conformational landscapes like those associated with cellular signaling proteins in the kinase family is more direct and readily converged than the corresponding free energy changes along the physical path. The R-FEP-R method was developed from the dual topology free energy perturbation method that is widely applied to estimate the binding free energy difference between two ligands. In R-FEP-R, the free energy change between two conformational basins is calculated by free energy perturbations that remove those atoms involved in the conformational change from their initial conformational basin while simultaneously growing them back according to the final conformational basin. Both the initial and final dual topology states are unphysical, but they are designed in a way such that the unphysical contributions to the initial and final partition functions cancel. Compared with other advanced sampling algorithms such as umbrella sampling and metadynamics, the R-FEP-R method does not require predetermined transition pathways or reaction coordinates that connect the two conformational states. As a first illustration, the R-FEP-R method was applied to calculate the free energy change between conformational basins for alanine dipeptide in solution and for a side chain in the binding pocket of T4 lysozyme. The results obtained by R-FEP-R agree with the benchmarks very well.

The emergence of multidrug-resistant pathogens led to a critical need for new antibiotics. A key property of effective antibiotics against Gram-negative bacteria is their ability to permeate through the bacterial outer membrane via transmembrane porin proteins. Molecular dynamics (MD) simulations are, in principle, capable of modeling antibiotic permeation across outer membrane porins (OMPs). However, owing to sampling problems, it has remained challenging to obtain converged potentials of mean force (PMFs) for antibiotic permeation across OMPs. Here, we investigated the convergence of PMFs along a single collective variable aimed at quantifying the permeation of the antibiotic fosmidomycin across the OprO porin. We compared standard umbrella sampling (US) with three advanced flavors of the US technique: (i) Hamiltonian replica exchange with solute tempering in combination with US, (ii) simulated tempering-enhanced US, and (iii) replica-exchange US. To quantify the PMF convergence and to reveal hysteresis problems, we computed several independent sets of US simulations starting from pulling simulations in the outward and inward permeation directions. We find that replica-exchange US in combination with well-chosen restraints is highly successful for obtaining converged PMFs of fosmidomycin permeation through OprO, reaching PMFs converged to subkilocalorie per mole accuracy.

Membrane channels are an essential part of any life form. They conduct the selective flux across the cell membrane of many important molecules that would otherwise not permeate. Experimental studies on membrane channels have led to the structural and functional characterization of many of them, yet many underlying physico-chemical mechanisms are somewhat out of reach. The aim of this thesis is to gain quantitative understanding on the structural and functional properties of these proteins by means of computational methods, such as Molecular Dynamics (MD) and free energy calculations. One of the most common approaches to study the selectivity and permeation mechanisms of a channel is the calculation of the Potential of Mean Force (PMF) for solute permeation across the pore. Usually, PMFs are calculated via MD simulations, which requires a significant amount of computational power. Hence, we compared the capability of MD with that of 3-Dimensional Reference Interaction-Site Model (3D-RISM), allegedly as accurate as MD but much more computationally efficient, to compute PMFs of solute permeation across Urea Transporter B (UT-B) and Aquaporin 1 (AQP1). We found a remarkable agreement between the PMFs for water permeation calculated from both techniques. However, for the rest of tested solutes, namely ammonia, urea, molecular oxygen, and methanol, we found critical discrepancies between 3D-RISM and with MD, which were found to be independent of the closure relation, the choice of the reaction coordinate, or the fluctuations of the protein. This suggests that, whilst 3D-RISM may provide reasonable approximations on PMFs for the permeation of water, it is not appropriate to study the permeation of uncharged non-water solutes. We further investigated, via a combination of MD simulations and free energy calculations, the structure and function of the fluoride-specific channel Fluc-Bpe. The free energy calculations allowed us to ascertain the specific nature of five isolated electron densities found in the crystal structure of Fluc, four of which were provisionaly assigned to fluoride, and the remaining one to sodium. We conducted two different kinds of binding free energy calculations: i) relative binding free energy differences ∆∆Gbind, and ii) absolute binding free energy ∆Gbind. Notably, the calculation of ∆∆Gbind allowed us to determine, between two putative molecular species, namely water and fluoride, which species was more likely to bind at a certain binding site. The resulting free energies were partly dependent on fluoride-phenylalanine interactions, which we found to be underestimated by ~ 30 kJ mol-1 in current additive force-fields. Thus, the disctimination of one species over the other was only possible because the ∆∆Gbind. values largely deviated from zero. In turn, the calculation of ∆Gbind allowed us to confirm whether a certain species would bind per se to Fluc-Bpe. Besides, short, free MD simulations proved to be key to assess the structural stability of the channel in different conditions, which, together with the free energy calculations, indicated that the four densities assigned to fluoride rather corresponded to ordered water molecules, and that the last electron density corresponded to a structural sodium. We finally evaluated, using MD simulations, the response of Fluc-Bpe to the presence of fluoride ions restrained at the permeating pore. The results suggested that the channel would undergo an opening transition, after which water molecules enter the pore to solvate the ions. Then, we calculated the PMFs for the permeation of water, fluoride and chloride using Umbrella Sampling (US) simulations. The profiles of solute permeation across the open structure indicated that water, fluoride, chloride would efficiently permeate the channel, being in stark contrast with the experimental evidence, which demonstrates that Fluc channels permeate fluoride by a ~ 100-fold ratio over chloride. We suspect that our results might be affected by the inaccurate modelling of ion-protein contacts highlighted before. The proper modelling of ion-protein interactions is extremely important for the establishment of salt-bridges, the structural stability of proteins, or the permeation of ions. Therefore, we conclude that our results regarding the permeation mechanism in Fluc-Bpe mainly reflect the imperfections of current additive force-fields, and that the usage of polarizable force-fields and development of accurate ion-protein interactions may certainly aid future research.

A practical approach that enables one to calculate the standard free energy of binding from a one-dimensional potential of mean force (PMF) is proposed. Umbrella sampling and the weighted histogram analysis method are used to generate a PMF along the reaction coordinate of binding. At each point, a restraint is applied orthogonal to the reaction coordinate to make possible the determination of the volume sampled by the ligand. The free energy of binding from an arbitrary unbound volume to the restrained bound form is calculated from the ratio of the PMF integrated over the bound region to that of the unbound. Adding the free energy changes from the standard-state volume to the unbound volume and from the restrained to the unrestrained bound state gives the standard free energy of binding. Exploration of the best choice of binding paths is also made. This approach is first demonstrated on a model binding system and then tested on the benzamidine-trypsin system for which reasonable agreement with experiment is found. A comparison is made with other methods to obtain the standard free energy of binding from the PMF.

Characterizing protein-protein association quantitatively has been a longstanding challenge for computer simulations. Here, a theoretical framework is put forth that addresses this challenge on the basis of detailed all-atom molecular dynamics simulations with explicit solvent. The proposed methodology relies upon independent potential of mean force (PMF) free-energy calculations carried out sequentially, wherein the biological objects are restrained in the conformation, position and orientation of the bound state, using adequately chosen biasing potentials. These restraints systematically narrow down the configurational entropy available to the system and effectively guarantee that the relevant network of interactions is properly sampled as the two proteins reversibly associate. Decomposition of the binding process into consecutive, well-delineated stages, for both the protein complex and the individual, unbound partners, offers a rigorous definition of the standard state, from which the absolute binding free energy can be determined. The method is applied to the difficult case of the extracellular ribonuclease barnase binding to its intracellular inhibitor barstar. The calculated binding free energy is -21.0 ± 1.4 kcal/mol, which compares well with the experimental value of -19.0 ± 0.2 kcal/mol. The relatively small statistical error reflects the precision and convergence afforded by the PMF-based simulation methodology. In addition to providing an accurate reproduction of the standard binding free energy, the proposed strategy offers a detailed picture of the protein-protein interface, illuminating the thermodynamic forces that underlie reversible association. The application of the present formal framework to barnase:barstar binding provides a foundation for tackling nearly any protein-protein complex.

The recent advances in relative protein-ligand binding free energy calculations have shown the value of alchemical methods in drug discovery. Accurately assessing absolute binding free energies, although highly desired, remains a challenging endeavour, mostly limited to small model cases. Here, we demonstrate accurate first principles based absolute binding free energy estimates for 128 pharmaceutically relevant targets. We use a novel rigorous method to generate protein-ligand ensembles for the ligand in its decoupled state. Not only do the calculations deliver accurate protein-ligand binding affinity estimates, but they also provide detailed physical insight into the structural determinants of binding. We identify subtle rotamer rearrangements between apo and holo states of a protein that are crucial for binding. When compared to relative binding free energy calculations, obtaining absolute binding free energies is considerably more challenging in large part due to the need to explicitly account for the protein in its apo state. In this work we present several approaches to obtain apo state ensembles for accurate absolute ΔG calculations, thus outlining protocols for prospective application of the methods for drug discovery.

Absolute Binding Free Energy Calculations Using Molecular Dynamics Simulations with Restraining Potentials

Recent advances in relative protein-ligand binding free energy calculations have shown the value of alchemical methods in drug discovery. Accurately assessing absolute binding free energies remains a challenging endeavour, mostly limited to small model cases. We demonstrate accurate absolute binding free energy estimates for 128 pharmaceutically relevant ligands across 7 proteins using a highly parallelizable non-equilibrium method. These calculations also provide detailed physical insight into the structural determinants of binding, identifying subtle rotamer rearrangements between protein apo and holo states that are crucial for binding. The challenge behind absolute binding free energy calculations stems in large part from the need to explicitly account for the protein’s apo state. In this work we present several approaches to obtain apo state ensembles, including a novel rigorous method to generate protein-ligand ensembles for the ligand in its decoupled state. Altogether, we present an effective open-source protocol for prospective application in drug discovery.

We present an approach combining alchemical modifications and physical-pathway methods to calculate absolute binding free energies. The employed physical-pathway method is either a stratified umbrella sampling to calculate a potential of mean force or nonequilibrium pulling. We devised two basic approaches: the simultaneous approach (S-approach), where, along the physical unbinding pathway, an alchemical transformation of ligand-protein interactions is installed and deinstalled, and the prior-plus-simultaneous approach (PPS-approach), where, prior to the physical-pathway simulation, an alchemical transformation of ligand-protein interactions is installed in the binding site and deinstalled during the physical-pathway simulation. Using a mutant of T4 lysozyme with a benzene ligand as an example, we show that installation and deinstallation of soft-core interactions concurrent with physical ligand unbinding (S-approach) allow successful potential of mean force calculations and nonequilibrium pulling simulations despite the problems posed by the occluded nature of the lysozyme binding pocket. Good agreement between the potential of the mean-force-based S-approach and double decoupling simulations as well as a remarkable efficiency and accuracy of the nonequilibrium-pulling-based S-approach is found. The latter turned out to be more compute-efficient than the potential of mean force calculation by approximately 70%. Furthermore, we illustrate the merits of reducing ligand-protein interactions prior to potential of mean force calculations using the murine double minute homologue protein MDM2 with a p53-derived peptide ligand (PPS-approach). Here, the problem of breaking strong interactions in the binding pocket is transferred to a prior alchemical transformation that reduces the free-energy barrier between the bound and unbound state in the potential of mean force. Besides, disentangling physical ligand displacement from the deinstallation of ligand-protein interactions was seen to allow a more uniform sampling of distance histograms in the umbrella sampling. In the future, physical ligand unbinding combined with simultaneous alchemical modifications may prove useful in the calculation of protein-protein binding free energies, where sampling problems posed by multiple, possibly sticky interactions and potential steric clashes can thus be reduced.

We report a large thermodynamic effect of solvent buffer on the standard binding free energy for a large hydrophobic ligand and show that a realistic comparison with the experimental binding affinity requires correctly accounting for the solvent reference state and ligand reorganization. We focus on lenacapavir (LEN; MW ≈ 1 kDa), an HIV-1 capsid inhibitor with very low aqueous solubility. Using several absolute binding free energy (ABFE) protocols including double-decoupling method (DDM), potential-of-mean-force (PMF) approaches, and the Alchemical Transfer Method (ATM), we obtained standard binding free energy values in neat water of ΔGbind0 = -26.4 to -30.0 kcal/mol for LEN binding to the HIV-1 capsid CA dimer, much stronger than the SPR-derived affinity measured in 5% DMSO buffer (≈-13.4 kcal/mol). We analyze the discrepancy and identify two dominant contributors to the calculated overbinding. (i) Solvent reference state: A thermodynamic cycle analysis and solvation free energy calculations reveal that the 5% DMSO buffer stabilizes the free ligand by ∼-4 kcal/mol relative to neat water. Additionally, we propose that the enrichment of the hydrophobic cosolvent in the apo binding pocket makes it more energetically costly to displace DMSO upon ligand binding. (ii) Ligand reorganization: incomplete treatment of LEN's internal conformational reorganization in the ABFE protocols leads to overbinding; a DDM variant with explicit ligand reorganization reduces the overestimate by ∼6 kcal/mol. Together, considerations of these effects significantly reduce the discrepancy between the ABFE calculations and experiments. Our results suggest that for large, hydrophobic ligands, quantitative agreement between ABFEs and experiments requires (a) reporting ΔGbind0 in the appropriate assay buffer (not simply in water) and (b) explicit treatment of ligand reorganization.

An important task of biomolecular simulation is the calculation of relative binding free energies upon chemical modification of partner molecules in a biomolecular complex. The potential of mean force (PMF) along a reaction coordinate for association or dissociation of the complex can be used to estimate binding affinities. A free energy perturbation approach, termed umbrella sampling (US) perturbation, has been designed that allows an efficient calculation of the change of the PMF upon modification of a binding partner based on the trajectories obtained for the wild type reference complex. The approach was tested on the interaction of modified water molecules in aqueous solution and applied to in silico alanine scanning of a peptide-protein complex. For the water interaction test case, excellent agreement with an explicit PMF calculation for each modification was obtained as long as no long range electrostatic perturbations were considered. For the alanine scanning, the experimentally determined ranking and binding affinity changes upon alanine substitutions could be reproduced within 0.1-2.0 kcal/mol. In addition, good agreement with explicitly calculated PMFs was obtained mostly within the sampling uncertainty. The combined US and perturbation approach yields, under the condition of sufficiently small system modifications, rigorously derived changes in free energy and is applicable to any PMF calculation.

Molecular dynamics (MD) simulations are widely applied to estimate absolute binding free energies of protein-ligand and protein-protein complexes. A routinely used method for binding free energy calculations with MD is umbrella sampling (US), which calculates the potential of mean force (PMF) along a single reaction coordinate. Surprisingly, in spite of its widespread use, few validation studies have focused on the convergence of the free energy computed along a single path for specific cases, not addressing the reproducibility of such calculations in general. In this work, we therefore investigate the reproducibility and convergence of US along a standard distance-based reaction coordinate for various protein-protein and protein-ligand complexes, following commonly used guidelines for the setup. We show that repeating the complete US workflow can lead to differences of 2-20 kcal/mol in computed binding free energies. We attribute those discrepancies to small differences in the binding pathways. While these differences are unavoidable in the established US protocol, the popularity of the latter could hint at a lack of awareness of such reproducibility problems. To test if the convergence of PMF profiles can be improved if multiple pathways are sampled simultaneously, we performed additional simulations with an adaptive-biasing method, here the accelerated weight histogram (AWH) approach. Indeed, the PMFs obtained from AHW simulations are consistent and reproducible for the systems tested. To the best of our knowledge, our work is the first to attempt a systematic assessment of the pitfalls in one the most widely used protocols for computing binding affinities. We anticipate therefore that our results will provide an incentive for a critical reassessment of the validity of PMFs computed with US, and make a strong case to further benchmark the performance of adaptive-biasing methods for computing binding affinities.

The calculation of binding free energies between highly charged species is a major challenge for free energy simulations. In this study, we applied a combination of molecular dynamics simulations and continuum electrostatics, along with normal-mode analysis, to compute the absolute free energies of binding between cobalt (III) hexammine and two RNA fragments for which NMR structures have been determined. The predicted affinities, using the finite-difference Poisson−Boltzmann method with a solute dielectric constant of 1 to treat solvation, were overall underestimated relative to the experimental values. However, internal consistency in the calculated energies was maintained between the different trajectories, and the structures in the simulations gave excellent agreement with NMR data. Various models for obtaining the electrostatic contributions were analyzed, including the effects of solute dielectric constants and van der Waals radii, linear and nonlinear salt contributions, as well as results from gene...